#include <fstream>
#include <iostream>
#include <deque>
#include <set>
#include <string>
#include <iterator>
#include <utility>
#include <algorithm>
#include <limits>
using namespace std;
#include <opencv2/core/core.hpp>
#include <opencv2/imgproc/imgproc.hpp>
#include <opencv2/highgui/highgui.hpp>
using namespace cv;
const Point2d origin(0, 0);
const Scalar White = CV_RGB (255, 255, 255);
const Scalar Black = CV_RGB (0, 0, 0);
const Scalar Red = CV_RGB (255, 0, 0);
const Scalar Lime = CV_RGB (0, 255, 0);
const Scalar Green = CV_RGB (0, 128, 0);
const Scalar Blue = CV_RGB (0, 0, 255);
const Scalar Yellow = CV_RGB (255, 255, 0);
const Scalar Orange = CV_RGB (255, 128, 0);
const Scalar Magenta = CV_RGB (255, 0, 255);
const Scalar Cyan = CV_RGB (0, 255, 255);
const Scalar Olive = CV_RGB (128,128,0);
const Scalar Gray = CV_RGB (128,128,128);
const Scalar Silver = CV_RGB (192,192,192);
const Scalar Maroon = CV_RGB (128,0,0);
const Scalar Purple = CV_RGB (128,0,128);
const Scalar Teal = CV_RGB (0,128,128);
const Scalar Navy = CV_RGB (0,0,128);
const Scalar Brown = CV_RGB (165,42,42);
const Scalar Coral = CV_RGB (255,127,80);
const Scalar Salmon = CV_RGB (250,128,114);
const Scalar Khaki = CV_RGB (240,230,140);
const Scalar Indigo = CV_RGB (75,0,130);
const Scalar Plum = CV_RGB  (221,160,221);
const Scalar orchid = CV_RGB (218,112,214);
const Scalar beige = CV_RGB (245,245,220);
const Scalar peru = CV_RGB (205,133,63);
const Scalar sienna = CV_RGB (160,82,45);
const Scalar lavender = CV_RGB (230,230,250);
const Scalar mocassin = CV_RGB (255,228,181);
const Scalar honeydew = CV_RGB (240,255,240);
const Scalar ivory = CV_RGB (255,255,240);
const Scalar azure = CV_RGB (240,255,255);
const Scalar crimson = CV_RGB (220,20,60);
const Scalar gold = CV_RGB (255,215,0);
const Scalar sky_blue = CV_RGB (135,206,235);
const Scalar aqua_marine = CV_RGB (127,255,212);
const Scalar bisque = CV_RGB (255,228,196);
const Scalar peach_puff = CV_RGB (255,218,185);
const Scalar corn_silk = CV_RGB (255,248,220);
const Scalar wheat = CV_RGB (245,222,179);
const Scalar violet = CV_RGB (238,130,238);
const Scalar lawn_green  = CV_RGB (124,252,0);
vector<CvScalar> Colors = { Lime, Orange, Gray, Cyan, Yellow, Red, Green, Olive, Magenta, Teal, Silver, Coral, Salmon, Khaki, Plum, orchid, beige, lavender, mocassin, honeydew, ivory, azure, crimson, gold, sky_blue, aqua_marine, bisque, peach_puff, corn_silk, wheat, violet, lawn_green};
// dark colors: Blue, Black, Brown, Maroon, Navy, Purple, Indigo, peru, sienna,
#include "boost/graph/adjacency_list.hpp"
#include "boost/graph/topological_sort.hpp"
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/dijkstra_shortest_paths.hpp>
#include <boost/graph/graphviz.hpp>
using namespace boost;
#include <CGAL/enum.h>
#include <CGAL/Delaunay_triangulation_2.h>
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Triangulation_vertex_base_with_info_2.h>
#include <CGAL/Triangulation_face_base_with_info_2.h>
using namespace CGAL;

class Node
{
public:
    int index; // index of triangle in FaceHandles
    int uplink;
    double birth;
    int height;
    int bar; // index of the region containing the triangle of this node
    vector<int> live; // indices of live triangles from Tree of this node
    //int weight;
    
    Node ()
    {
        index = 0;  // index in FaceHandles: 0 means external boundary
        uplink = 0; // index of the parent node in FaceHandles
        birth = 0; //
        height = 1; // no nodes below the new node
        bar = -1;
        live.clear();
    }
};

typedef Exact_predicates_inexact_constructions_kernel K;
typedef Triangulation_vertex_base_with_info_2<int, K> VB;
typedef Triangulation_face_base_with_info_2<Node, K> FB;
typedef Triangulation_data_structure_2<VB,FB> TDS;
typedef Delaunay_triangulation_2<K, TDS> DT;
typedef K::Point_2 P2;
typedef DT::Edge DE;
typedef DT::Face DF;
typedef DT::Vertex_handle VH;
typedef DT::Face_handle FH;
typedef DT::Vertex_iterator VI;
typedef DT::Edge_iterator EI;
typedef DT::Finite_vertices_iterator FVI;
typedef DT::Finite_edges_iterator FEI;
typedef DT::All_faces_iterator AFI;
typedef DT::Finite_faces_iterator FFI;

pair<VH,VH> Vertices (DE const& e)
{
    // Edges are a triplet; (face, i, j). Get the i-th (either 0, 1, or 2) vertex of a face using vertex(i)
    VH v1 = e.first -> vertex ( (e.second + 1) % 3 );
    VH v2 = e.first -> vertex ( (e.second + 2) % 3 );
    return make_pair (v1, v2);
}

pair<FH,FH> Faces (DE const& e)
{
    //if you have an edge e, then e.first and e.first.neighbor( e.second ) are the two incident facets.
    FH f1 = e.first;
    FH f2 = e.first->neighbor( e.second );
    return make_pair (f1, f2);
}//

double SquaredLength (DE const& e)
{
    VH v1 = Vertices(e).first;
    VH v2 = Vertices(e).second;
    return squared_distance( v1->point(), v2->point() );
}

class Edge
{
public:
    DE edge;
    VH end1, end2;
    int face1, face2;
    double length;
    bool boundary;
    
    Edge () {}
    Edge (DE const& e)
    {
        edge = e;
        pair<VH,VH> ends = Vertices(e);
        end1 = ends.first;
        end2 = ends.second;
        pair<FH,FH> faces = Faces(e);
        FH face = faces.first;
        face1 = face->info().index;
        face = faces.second;
        face2 = face->info().index;
        length = sqrt( SquaredLength(e) );
        boundary = false;
        //UpdateLength();
    }
    void UpdateLength() { length = sqrt( squared_distance( end1->point(), end1->point() ) ); }
};

bool DecreasingLengths (Edge const& e1, Edge const& e2)
{
    return e2.length < e1.length;
}

bool IncreasingLengths (Edge const& e1, Edge const& e2)
{
    return e2.length > e1.length;
}

void WriteImage (String const& name, Mat const& image)
{
    try{ imwrite( name, image ); }
    catch (runtime_error& ex) { fprintf(stderr, "Exception converting to PNG: %s\n", ex.what()); }
}

void LoadCloud (int SizeCloud, String const& FileName, vector<P2>& PointCloud)
{
    ifstream input( FileName );
    string line_data;
    char comma;
    while ( getline(input, line_data) )
    {
        stringstream line_stream (line_data);
        double x,y;
        if (line_stream >> x >> comma >> y) PointCloud.push_back( P2(x,y) );
    }
}

void load_clouds( String const& InputFolder, vector<String>const& file_names, String const& add_ext, vector< vector<P2>>& clouds )
{
    clouds.resize( file_names.size() );
    for (size_t f=0; f<file_names.size(); f++)
    {
        cout<<"\nOpen file"<<InputFolder + file_names[f] + add_ext;
        ifstream input( InputFolder + file_names[f] + add_ext, ifstream::in );
        string line_data;
        getline( input, line_data ); // ignore first line
        cout<<line_data<<endl;
        while ( getline( input, line_data ) )
        {
            //cout<<line_data<<endl;
            stringstream line_stream( line_data );
            double x,y;
            if (line_stream >> x >> y) clouds[f].push_back( P2(x,y) );
        }
        input.close();
    }
}

Point2i ScaleShift(P2 p, double scale, Point2i shift)
{
    return Point2i ( int( scale * p.x() + shift.x ), int( scale * p.y() + shift.y ) );
}

Point2i ScaleShift(Point2d p, double scale, Point2i shift)
{
    return Point2i ( int( scale * p.x + shift.x ), int( scale * p.y + shift.y ) );
}

void draw_cloud( vector<P2>const& cloud, int MaxIndex, double scale, Point2d shift, Mat& image )
{
    for (size_t i=0; i<std::min((int)cloud.size(),MaxIndex+1); i++)
        circle( image, ScaleShift( cloud[i], scale, shift), 2, Black, -1);
}

void DrawVertices ( DT const& Del, double scale, Point2d shift, Mat& image )
{
    for( FVI it = Del.finite_vertices_begin(); it != Del.finite_vertices_end(); ++it)
    {
        Point2i p = ScaleShift( it->point(), scale, shift);
        circle( image, p,  2, Red, -1);
        String s = "("+to_string( (int)it->point().x() )+","+to_string( (int)it->point().y() )+")";
        putText(image, s, p, FONT_HERSHEY_PLAIN, 1.0, Red);
    }
}

void DrawEdge ( P2 v1, P2 v2, double scale, Point2d shift, CvScalar c, Mat& image )
{
    line( image, ScaleShift( v1, scale, shift), ScaleShift( v2, scale, shift), c, 1);
}

void DrawEdges ( DT const& Del, double scale, Point2d shift, Mat& image )
{
    for( FEI it = Del.finite_edges_begin(); it != Del.finite_edges_end(); ++it)
        DrawEdge ( Vertices(*it).first->point(), Vertices(*it).second->point(), scale, shift, Blue, image);
}

String Round (double v)
{
    char b[9];
    sprintf( b, "%.1f", v );
    return string(b);
}

void DrawFaces ( DT const& Del, double scale, Point2d shift, Mat& image )
{
    for(FFI iF = Del.finite_faces_begin(); iF != Del.finite_faces_end(); ++iF)
    {
        P2 p0 = iF->vertex(0)->point();
        P2 p1 = iF->vertex(1)->point();
        P2 p2 = iF->vertex(2)->point();
        P2 p ( (p0.x()+p1.x()+p2.x())/3, (p0.y()+p1.y()+p2.y())/3 );
        String s = to_string( iF->info().index ) + "," + Round( iF->info().birth );
        putText(image, s, ScaleShift(p , scale, shift), FONT_HERSHEY_PLAIN, 1.0, Black);
        //P2 p01 ( (p0.x()+p1.x())/2, (p0.y()+p1.y())/2 );
        //putText(image, to_string(iF->info().index), ScaleShift(p , scale, shift), FONT_HERSHEY_PLAIN, 1.0, Black);
    }
}

void DrawTriangulation (DT const& Del, double scale, Point2d shift, Mat& image)
{
    //DrawFaces( Del, scale, shift, image );
    DrawEdges( Del, scale, shift, image );
    //DrawVertices( Del, scale, shift, image );
}

bool FindCloudBox ( vector<P2>const& PointCloud, pair<Point2d, Point2d>& CloudBox ) {
    CloudBox.first = Point2d (1e+64, 1e+64);
    CloudBox.second = Point2d (-1e+64, -1e+64);
    for (size_t i=0; i<PointCloud.size(); i++)
    {
        if (PointCloud[i].x() < CloudBox.first.x) CloudBox.first.x = PointCloud[i].x();
        if (PointCloud[i].x() > CloudBox.second.x) CloudBox.second.x = PointCloud[i].x();
        if (PointCloud[i].y() < CloudBox.first.y) CloudBox.first.y = PointCloud[i].y();
        if (PointCloud[i].y() > CloudBox.second.y) CloudBox.second.y = PointCloud[i].y();
    }
    return true;
}

bool FindCloudBox ( vector< pair<int, P2> >const& PointCloud, pair<Point2d, Point2d>& CloudBox ) {
    CloudBox.first = Point2d (1e+64, 1e+64);
    CloudBox.second = Point2d (-1e+64, -1e+64);
    for (size_t i=0; i<PointCloud.size(); i++)
    {
        if (PointCloud[i].second.x() < CloudBox.first.x) CloudBox.first.x = PointCloud[i].second.x();
        if (PointCloud[i].second.x() > CloudBox.second.x) CloudBox.second.x = PointCloud[i].second.x();
        if (PointCloud[i].second.y() < CloudBox.first.y) CloudBox.first.y = PointCloud[i].second.y();
        if (PointCloud[i].second.y() > CloudBox.second.y) CloudBox.second.y = PointCloud[i].second.y();
    }
    return true;
}

int FindRoot (int node, vector<FH>const& FaceHandles)
{
    int next = FaceHandles[node]->info().uplink;
    if (next == node) return node;
    return FindRoot( next, FaceHandles );
}

bool is_acute (K::Triangle_2 const& Triangle, double& radius)
{
    P2 p = circumcenter( Triangle.vertex(0), Triangle.vertex(1), Triangle.vertex(2) );
    //radius = distance ( p, Triangle.vertex(0) );
    radius = sqrt( squared_distance( p, Triangle.vertex(0) ) );
    if ( Triangle.bounded_side(p) == ON_BOUNDED_SIDE ) return true;
    return false;
}

class Life
{
public:
    int index;
    double birth;
    double death;
    double span;
    int edge; // index of critical edge in DelEdges
    Life (int i, double b, double d, int e) { index = i; birth = b; death = d; span = d-b; edge = e; }
};

class Region
{
public:
    int index; // index of this region in Map
    int edge; // index of the critical edge in DelEdges that created this region
    double birth;
    double death;
    double span;
    vector<int> core; // indices of triangles that were last alive before this region died
    int heir; // index of the root triangle in the older region that absorbed this region
    int supr; // index of the region (in Map) defined by the triangle heir, computed at the end
    Region (int i, double b, double d, vector<int>const& c, int h, int s, int e)
    { index = i; birth = b; death = d; span = d-b; core = c; heir = h; supr = s; edge = e; }
};

bool DecreasingSpans (Region const& r1, Region const& r2) { return r1.span >= r2.span; }

void DrawHopes (vector<Edge>const& DelEdges, vector<int>const& NegEdges, vector<int>const& PosEdges,
                double scale, Point2d shift, Mat& iHopes)
{
    for (size_t i=0; i<NegEdges.size(); i++)
        DrawEdge( DelEdges[NegEdges[i]].end1->point(), DelEdges[NegEdges[i]].end2->point(),scale,shift,Black,iHopes);
    for (size_t i=0; i<PosEdges.size(); i++)
        DrawEdge( DelEdges[PosEdges[i]].end1->point(), DelEdges[PosEdges[i]].end2->point(), scale, shift, Red,iHopes);
}

void DrawMap (vector<FH>const& FaceHandles, vector<Region>const& Map, double scale, Point2d shift, Mat& image)
{
    int ind = -1;
    for (int j=0; j<Map.size(); j++)
    {
        if ( Map[j].core.size() == 0) continue; else ind++;
        Scalar c = Colors[ ind % Colors.size() ];
        if ( Map[j].index == Map.size()-1 ) c = White; // external region
        for (size_t i=0; i<Map[j].core.size(); i++)
        {
            if (Map[j].core[i]==0) continue; // external triangle doesn't exist
            FH iF = FaceHandles[ Map[j].core[i] ];
            Point Poly[3];
            Poly[0] = ScaleShift( iF->vertex(0)->point(), scale, shift );
            Poly[1] = ScaleShift( iF->vertex(1)->point(), scale, shift );
            Poly[2] = ScaleShift( iF->vertex(2)->point(), scale, shift );
            //cout<<"\nR"<<j<<"T"<<Map[j].core[i]<<"("<<iF->vertex(0)->point()<<")("<<iF->vertex(1)->point()<<")("<<iF->vertex(2)->point()<<")";
            fillConvexPoly( image, Poly, 3, c );
            //P2 p ( (p0.x()+p1.x()+p2.x())/3, (p0.y()+p1.y()+p2.y())/3 );
            //String s = "R" + to_string(j) + "," + to_string( iF->info().index );
            //putText(image, s, ScaleShift(p , scale, shift), FONT_HERSHEY_PLAIN, 1.0, Black);
        }
    }
}

void FindBoundary (int nRegions, vector<Region>const& Map, vector<Edge>& DelEdges)
{
    vector<int> regions( nRegions );
    for (int i=0; i<Map.size(); i++)
        for (size_t j=0; j<Map[i].core.size(); j++)
            regions[ Map[i].core[j] ] = i; // this triangle belongs to i-th region in Map
    for (size_t e=0; e<DelEdges.size(); e++)
        if ( regions[ DelEdges[e].face1 ] != regions[ DelEdges[e].face2 ] )
            DelEdges[e].boundary = true;
}

void DrawBoundary (vector<Edge>const& DelEdges, double scale, Point2d shift, Mat& image)
{
    for (size_t e=0; e<DelEdges.size(); e++) if ( DelEdges[e].boundary )
    {
        P2 p1 = DelEdges[e].end1->point();
        P2 p2 = DelEdges[e].end2->point();
        line( image, ScaleShift( p1, scale, shift ), ScaleShift( p2, scale, shift ), Black, 2 );
    }
}

void DrawDiagram (vector<Life>const& Persistence, int iAbove, double CritLength, double MaxRadius, double scale, Point2i shift, Mat& image)
{
    size_t nDots = Persistence.size();
    double d = MaxRadius;
    Point Poly[4];
    double s = Persistence[iAbove].span, c = 0.5*CritLength, y = 0;
    if (iAbove>0) y = Persistence[iAbove-1].span;
    Poly[0] = ScaleShift( Point2d( 0, d - s ), scale, shift);
    Poly[1] = ScaleShift( Point2d( 0, d - y ), scale, shift);
    Poly[2] = ScaleShift( Point2d( d - y, 0 ), scale, shift);
    Poly[3] = ScaleShift( Point2d( d - s, 0 ), scale, shift);
    fillConvexPoly( image, Poly, 4, Yellow);
    /* Prepare VerSubDiagram
    int iRight = -1; double xRight = 1e+8;
    for (int i=iAbove; i<nDots; i++)
        if (Persistence[i].birth > c)
            if (Persistence[i].birth < xRight)
            {
                xRight = Persistence[i].birth;
                iRight = i;
            }
    if ( iRight < 0 ) // Draw triangular VerSubDiagram
    {
        Poly[0] = ScaleShift( Point2d( c, d-s-c), scale,shift);
        Poly[1] = ScaleShift( Point2d( c, 0 ), scale, shift);
        Poly[2] = ScaleShift( Point2d( d-s, 0 ), scale, shift);
        fillConvexPoly( image, Poly, 3, Lime);
    }
    else
    {   // Draw VerSubDiagram as a trapezium
        double r = Persistence[iRight].birth;
        Poly[0] = ScaleShift( Point2d( c, d-s-c), scale, shift);
        Poly[1] = ScaleShift( Point2d( c, 0 ), scale, shift);
        Poly[2] = ScaleShift( Point2d( r, 0 ), scale, shift);
        Poly[3] = ScaleShift( Point2d( r, d-s-r), scale, shift);
        fillConvexPoly( image, Poly, 4, Lime);
    }*/
    for (size_t i=0; i<nDots; i++) // Draw red dots
    {
        Point2d p ( Persistence[i].birth, d - Persistence[i].death );
        int size = 4;
        if (Persistence[i].birth > c) size = 2;
        if (i<iAbove) size = 2;
        circle( image, ScaleShift( p, scale, shift), size, Red, -1);
    }
    line( image, ScaleShift( Point2d(0,d), scale, shift), ScaleShift( Point2d(0,0), scale, shift), Black, 1);
    line( image, ScaleShift( Point2d(0,d), scale, shift), ScaleShift( Point2d(d,d), scale, shift), Black, 1);
    line( image, ScaleShift( Point2d(0,d), scale, shift), ScaleShift( Point2d(d,0), scale, shift), Black, 1);
}

bool IncreasingSpans (Life const& l1, Life const& l2) { return l1.span < l2.span; }

bool IncreasingBirths (Life const& l1, Life const& l2) { return l1.birth < l2.birth; }

class IndexValue
{
public:
    int index;
    double value;
    IndexValue () { value = 0; }
    IndexValue (int i, double v) { index=i; value=v; }
};

bool DecreasingValues (IndexValue const& p1, IndexValue const& p2) { return p1.value > p2.value; }

bool IncreasingValues (IndexValue const& p1, IndexValue const& p2) { return p1.value < p2.value; }

void BuildStaircase ( vector<Life>const& Persistence, int ind, vector<IndexValue>& Staircase)
{
    size_t nValues = Persistence.size();
    vector<IndexValue> Endpoints ( 2*(nValues-ind) );
    for (size_t i=ind; i<nValues; i++)
    {
        Endpoints[2*(i-ind)+0] = IndexValue( 1, Persistence[i].birth );
        Endpoints[2*(i-ind)+1] = IndexValue( -1, Persistence[i].death );
    }
    sort( Endpoints.begin(), Endpoints.end(), IncreasingValues);
    //cout<<endl; for (size_t i=2*ind; i<2*nValues; i++) cout<<" "<<Endpoints[i-2*ind].value;
    Staircase = Endpoints;
    for (size_t i=1; i<2*(nValues-ind); i++)
        Staircase[i].index += Staircase[i-1].index;
    //for (size_t i=0; i<2*(nValues-ind); i++) cout<<" s["<<Staircase[i].value<<"]="<<Staircase[i].index;
}

void DrawStaircase (vector<IndexValue>const& Staircase, Point2i shift, int SizeImage, Mat& image )
{
    int MaxHoles = 0;
    int nPoints = (int)Staircase.size();
    for (size_t i=0; i<nPoints; i++)
        if (MaxHoles < Staircase[i].index) MaxHoles = Staircase[i].index;
    double MaxValue = Staircase[nPoints-1].value;
    double xScale = 1.0 * SizeImage / MaxValue;
    double yScale = 1.0 * SizeImage / MaxHoles;
    //cout<<"\n MaxHoles="<<MaxHoles<<" MaxValue="<<MaxValue<<" xScale="<<xScale<<" yScale="<<yScale<<endl;
    line( image, Point2i(0,SizeImage)+shift, Point2i(0,0)+shift, Green, 1);
    line( image, Point2i(0,SizeImage)+shift, Point2i(SizeImage,SizeImage)+shift, Green, 1);
    int x0 = int( xScale * Staircase[0].value );
    int y0 = int( yScale * (MaxHoles - Staircase[0].index) );
    line( image, Point2i(0,SizeImage)+shift, Point2i(x0,SizeImage)+shift, Red, 2);
    line( image, Point2i(x0,SizeImage)+shift, Point2i(x0,y0)+shift, Red, 1);
    for (size_t i=0; i+1<nPoints; i++)
    {
        int x1 = int( xScale * Staircase[i].value );
        int x2 = int( xScale * Staircase[i+1].value );
        int y1 = int( yScale * (MaxHoles - Staircase[i].index) );
        int y2 = int( yScale * (MaxHoles - Staircase[i+1].index) );
        line( image, Point2i(x1,y1)+shift, Point2i(x2,y1)+shift, Red, 2);
        line( image, Point2i(x2,y1)+shift, Point2i(x2,y2)+shift, Red, 1);
    }
}

void CannyDetector( Mat const& InputImage, int LowThreshold, double ratio, Mat& CannyEdges)
{
    int SizeKernel = 3;
    Mat EdgesImage;
    blur( InputImage, EdgesImage, Size(3,3) );
    Canny( EdgesImage, CannyEdges, LowThreshold, ratio * LowThreshold, SizeKernel );
}

void LoadCloud (Mat const& ImageCloud, vector<P2>& PointCloud)
{
    for (int j=0; j<ImageCloud.rows; j++)
        for (int i=0; i<ImageCloud.cols; i++)
            if ( ImageCloud.at<uchar>(j, i) != 0 )
                PointCloud.push_back( P2(i,j) );
}

void BinaryToCloud (Mat const& image, double noise, double ratio, vector<P2>& PointCloud)
{
    //cout<<" sizes: rows="<<image.rows<<" cols="<<image.cols;
    double r = ratio * RAND_MAX;
    double factor = noise * std::min( image.rows, image.cols );
    for (int j=0; j<image.rows; j++)
        for (int i=0; i<image.cols; i++)
            if ( int( image.at<uchar>(j, i) ) == 255 || rand() > r ) continue;
            else
            {
                //if (j==20) cout<<" c("<<i<<","<<j<<")="<<int( image.at<uchar>(j, i) );
                double x = i + factor * rand() / RAND_MAX;
                double y = j + factor * rand() / RAND_MAX;
                PointCloud.push_back( P2(x,y) );
            }
}

void Print (vector<P2>const& PointCloud)
{
    cout<<"Point Cloud: ";
    for (int j=0; j<PointCloud.size(); j++) cout<<" ("<<PointCloud[j].x()<<","<<PointCloud[j].y()<<")";
}

void Print (vector<Life>const& P)
{
    for (int j=0; j<P.size(); j++) cout<<"\nP["<<j<<"]: ind="<<P[j].index<<" b="<<Round(P[j].birth)<<" d="<<Round(P[j].death)<<" span="<<Round(P[j].span);
}

void Print (vector<Region>const& Map)
{
    for (int j=0; j<Map.size(); j++)
    {
        cout<<"\nMap["<<j<<"]: ind="<<Map[j].index<<" b="<<Round(Map[j].birth)<<" d="<<Round(Map[j].death)<<" span="<<Round(Map[j].span)<<" h="<<Map[j].heir<<" s="<<Map[j].supr<<" core size="<<Map[j].core.size();
        for (int i=0; i<Map[j].core.size(); i++) cout<<" "<<Map[j].core[i];
    }
}

class Graph
{
public:
    vector<P2> vertices;
    vector< pair<int,int> > edges; // pairs of indices of vertices
    vector<double> edgelengths;
    vector< vector<int> > vedges; // indices of edges around a vertex
    Graph() {}
    Graph( vector<Edge>const& DelEdges, vector<int>const& EdgeIndices, int nVertices )
    {
        vertices.resize (nVertices);
        vedges.resize (nVertices);
        edges.resize( EdgeIndices.size() );
        edgelengths.resize( EdgeIndices.size() );
        for (int i=0; i<EdgeIndices.size(); i++)
        {
            Edge e = DelEdges[ EdgeIndices[i] ];
            edges[i].first = e.end1->info();
            edges[i].second = e.end2->info();
            edgelengths[i] = e.length;
            vertices[ e.end1->info() ] = e.end1->point();
            vertices[ e.end2->info() ] = e.end2->point();
            vedges[ e.end1->info() ].push_back( i );
            vedges[ e.end2->info() ].push_back( i );
        }
    }
    
    Graph (String type, int size)
    {
        if (type=="star")
        {
            if (size<2) size = 2;
            vertices.resize( size+1 );
            vertices[0] = P2 (0,0);
            edges.resize( size );
            edgelengths.resize( size );
            double angle = 2.0 * M_PI / size;
            for(int i=1; i<=size; i++)
            {
                edges[i-1] = make_pair( 0, i );
                edgelengths [i-1] = 1;
                double x = cos( i * angle );
                double y = sin( i * angle );
                vertices[i] = P2 (x,y);
            }
        }
        else if (type=="wheel")
        {
            if (size<3) size = 3;
            vertices.resize( size+1 );
            vertices[0] = P2 (0,0);
            edges.resize( 2*size );
            edgelengths.resize( 2*size );
            double angle = 2.0 * M_PI / size;
            double side = 2 * sin ( angle/2 );
            //cout<<" side = "<<side;
            for(int i=1; i<=size; i++)
            {
                edges[i-1] = make_pair( 0, i );
                edges[ size+i-1 ] = make_pair( i, i+1 );
                if (i==size) edges[ 2*size-1 ].second = 1;
                edgelengths [i-1] = 1;
                edgelengths [size+i-1] = side;
                double x = cos( i * angle );
                double y = sin( i * angle );
                vertices[i] = P2 (x,y);
            }
        }
        else if (type=="lattice")
        {
            if (size<1) size = 1;
            double side = 1.0 / (size-1);
            edgelengths.assign( 2*size*(size-1), side );
            vertices.resize( size * size );
            for(int j=0;j<size;j++)
                for(int i=0;i<size;i++)
                    vertices[ j * size + i ] = P2( i*side, j*side );
            for(int j=0; j<size;j++)
                for(int i=1; i<size;i++)
                    edges.push_back( make_pair( j*size+i, j*size+i-1) ); // from (i,j) to (i-1,j)
            for(int j=1; j<size; j++)
                for(int i=0; i<size; i++)
                    edges.push_back( make_pair( j*size+i, (j-1)*size+i) ); // from (i,j) to (i,j+1)
        }
        else if (type=="split-lattice")
        {
            if (size<1) size = 1;
            vertices.clear();
            edges.clear();
            // only main corners
            for (int i=0; i<=size; i++)
                for(int j=0; j<=size; j++)
                    vertices.push_back( P2(4*i,4*j) );
            // only 2 vertices on each horizontal edge
            for (int i=0;i<size;i++)
                for (int j=0;j<=size;j++)
                {
                    edges.push_back( make_pair ( (size+1)*i+j, vertices.size() ) );
                    vertices.push_back( P2( 4*i+2, 4*j ) );
                    edges.push_back( make_pair( (size+1)*(i+1)+j, vertices.size() ) );
                    vertices.push_back( P2( 4*i+3, 4*j ) );
                }
       	    // only 2 vertices on each vertical edge
            for(int i=0;i<=size;i++)
                for(int j=0;j<size;j++)
                {
                    edges.push_back( make_pair( (size+1)*i+j, vertices.size() ) );
                    vertices.push_back( P2( 4*i, 4*j+2 ) );
                    edges.push_back( make_pair( (size+1)*i+j+1, vertices.size() ) );
                    vertices.push_back( P2( 4*i, 4*j+3 ) );
                }
            // random edges in each small square
            for( int i=0;i<size;i++)
                for( int j=0;j<size;j++)
                {
                    double r = (double)rand() / RAND_MAX;
                    int s = (int)vertices.size();
                    if (r<0.3)
                    {
                        vertices.push_back( P2( 4*i+2, 4*j+1 ) );
                        vertices.push_back( P2( 4*i+2, 4*j+3));
                        edges.push_back( make_pair( s, s+1 ) );
                    }
                    else if (r<0.6)
                    {
                        vertices.push_back( P2( 4*i+1, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+3, 4*j+2 ) );
                        edges.push_back( make_pair( s, s+1 ) );
                    }
                    else if (r<0.8)
                    {
                        vertices.push_back( P2( 4*i+2, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+3, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+2, 4*j+3 ) );
                        edges.push_back( make_pair( s, s+1 ) );
                        edges.push_back( make_pair( s, s+2 ) );
                    }
                    else
                    {
                        vertices.push_back( P2( 4*i+2, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+1, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+2, 4*j+1 ) );
                        vertices.push_back( P2( 4*i+3, 4*j+2 ) );
                        vertices.push_back( P2( 4*i+2, 4*j+3 ) );
                        edges.push_back( make_pair( s, s+1 ) );
                        edges.push_back( make_pair( s, s+2 ) );
                        edges.push_back( make_pair( s, s+3 ) );
                        edges.push_back( make_pair( s, s+4 ) );
                    }
            }
            edgelengths.resize( edges.size() );
            for (size_t e=0; e<edges.size(); e++)
                edgelengths[e] = sqrt( squared_distance( vertices[edges[e].first], vertices[edges[e].second] ) );
        }
    }
    
    void LoadCloud (int SizeCloud, double NoiseLevel, vector<P2>& PointCloud)
    {
        if ( edges.size() == 0)
        {/*
            for(int i=0;i<numPoints;i++)
            {
                int v=(int)(Math.random()*graph.vertex.size());
                double x=graph.vertex.get(v).coord.get(0);
                double y=graph.vertex.get(v).coord.get(1);
                double randx=x+(2*Math.random()-1)*noise;
                double randy=y+(2*Math.random()-1)*noise;
                point.add(new Point(randx,randy));
            }*/
        }
        else
        {
            //double GraphLength = accumulate( edgelengths.begin(), edgelengths.end(), 0);
            double GraphLength = 0;
            for(int i=0; i<edgelengths.size(); i++) GraphLength += edgelengths[i];
            //cout<<" GraphLength="<<GraphLength;
            for(int i=0; i<SizeCloud; i++)
            {
                double length = GraphLength * rand() / RAND_MAX;
                if (length==0) length = GraphLength;
                size_t e=0;
                for (; e<edges.size(); e++)
                    if ( length > edgelengths[e] ) length -= edgelengths[e]; else break;
                //cout<<" l="<<length<<" e"<<e;
                double coef =  length / edgelengths[e];
                int v1 = edges[e].first;
                int v2 = edges[e].second;
                double randx = ( 2.0 * rand() / RAND_MAX - 1) * NoiseLevel;
                double randy = ( 2.0 * rand() / RAND_MAX - 1) * NoiseLevel;
                double x = coef * vertices[v1].x() + (1-coef) * vertices[v2].x() + randx;
                double y = coef * vertices[v1].y() + (1-coef) * vertices[v2].y() + randy;
                PointCloud.push_back( P2(x,y) );
            }
        }
    }

    void Print()
    {
        for (size_t i=0; i<vedges.size();i++)
        {
            cout<<"\n v"<<i<<"("<<vertices[i]<<") has e";
            for (size_t j=0; j<vedges[i].size();j++) cout<<" "<<vedges[i][j];
        }
        for (size_t i=0; i<edges.size();i++) cout<<"\n e"<<i<<"-v"<<edges[i].first<<"("<<vertices[edges[i].first]<<")-v"<<edges[i].second<<"("<<vertices[edges[i].second]<<") l="<<edgelengths[i];
    }

    void Draw( double scale, Point2i shift, Mat& image )
    {
        for (size_t i=0; i<edges.size(); i++)
        {
            int v1 = edges[i].first, v2 = edges[i].second;
            P2 p1 = vertices[v1], p2 = vertices[v2];
            DrawEdge( p1, p2, scale, shift, Blue, image);
            //circle( image, ScaleShift(p1, scale, shift), 2, Red, -1);
            //circle( image, ScaleShift(p2, scale, shift), 2, Red, -1);
            //Point2d p ( 0.5 * ( p1.x()+p2.x() ), 0.5 * ( p1.y()+p2.y() ) );
            //putText(image, to_string(i), ScaleShift(p, scale, shift), FONT_HERSHEY_PLAIN, 1.0, Black);
            //putText(image, to_string(v1), ScaleShift(p1, scale, shift), FONT_HERSHEY_PLAIN, 1.0, Red);
            //putText(image, to_string(v2), ScaleShift(p2, scale, shift), FONT_HERSHEY_PLAIN, 1.0, Red);
        }
    }
    //
    void Draw( double scale, Point2i shift, CvScalar c, Mat& image )
    {
        for (size_t i=0; i<edges.size(); i++)
        {
            int v1 = edges[i].first, v2 = edges[i].second;
            P2 p1 = vertices[v1], p2 = vertices[v2];
            DrawEdge ( p1, p2, scale, shift, c, image);
        }
        for (size_t i=0; i<vertices.size(); i++)
            circle( image, ScaleShift(vertices[i], scale, shift), 2, c, -1);
    }
    
    bool CollapseShortestEdge (double CritLength)
    {
        if (vertices.size()<=10) return false; // too small graph
        int e=0;
        double MinLength = 1e+8;
        for (int i=1; i<edges.size(); i++)
        {
            double l = squared_distance( vertices[edges[i].first], vertices[edges[i].second] );
            if ( l < MinLength ) { MinLength = l; e = i; }
        }
        if ( MinLength > CritLength) return false;
        int v1 = edges[e].first; // move to mid point
        int v2 = edges[e].second; // remove
        double x = 0.5 * (vertices[v1].x() + vertices[v2].x() );
        double y = 0.5 * (vertices[v1].y() + vertices[v2].y() );
        P2 v(x,y); // mid point is the new vertex
        vertices[v1] = v;
        int ev1=0, ev2=0;
        for (; ev1<vedges[v1].size(); ev1++)
            if (vedges[v1][ev1] == e) break;
        if (ev1>vedges[v1].size()) cout<<" ev1="<<ev1<<">=vedges[v1].size()"<<vedges[v1].size();
        vedges[v1].erase( vedges[v1].begin()+ev1 ); // erase ref to the removed edge
        for (int j=0; j<vedges[v2].size(); j++)
            if (vedges[v2][j] == e) ev2=j; // old edge from v2 to v1
            else
            {
                int ind = vedges[v2][j];
                vedges[v1].push_back (ind); // new edge from v2 added to v1
                if (edges[ind].first == v2) edges[ind].first = v1;
                else edges[ind].second = v1; // new endpoints instead of v2
            }
        vedges[v2].erase( vedges[v2].begin()+ev2 ); // erase ref to the removed edge
        edges.erase( edges.begin()+e );
        vedges.erase( vedges.begin()+v2 ); // no more edges from old v2
        for (size_t i=0; i<vedges.size();i++)
            for (size_t j=0; j<vedges[i].size();j++)
                if (vedges[i][j]>e) vedges[i][j]--;
                else if (vedges[i][j]==e) cout<<" vedges["<<i<<"]["<<j<<"]="<<e;
        vertices.erase( vertices.begin()+v2 );
        for (size_t i=0; i<edges.size();i++)
        {
            if (edges[i].first > v2) edges[i].first--;
            else if (edges[i].first == v2) cout<<" edges["<<i<<"].first="<<v2;
            if (edges[i].second > v2) edges[i].second--;
            else if (edges[i].second == v2) cout<<" edges["<<i<<"].second="<<v2;
        }
        return true;
    }
};

int main ()
{
    srand( (int)time(0) );
    int SizeImage = 532;
    int BorderImage = 16;
    Point2d ShiftImage( BorderImage, BorderImage );
    String InputFolder = "/Users/kurlin/C++/Datasets/BSD500/";
    //String InputFolder = "/Users/kurlin/C++/Datasets/CloudFiles/";
    //String InputFolder = "/Users/kurlin/C++/Datasets/TestImages/";
    //String InputFolder = "/Users/kurlin/C++/Datasets/Glasgow/";
    String OutputFolder = "/Users/kurlin/C++/AlphaComplexes2d/Output/";
    //String InputType = "image";
    //String InputType = "binary";
    //String InputType = "file";
    String InputType = "graph"; //"Glasgow";
    //String GraphType = "wheel";
    String GraphType = "lattice";
    //String GraphType = "split-lattice";
    //String InputFile = "figure-eight";
    String InputFile = "cloud-ears";
    //String ImageFile = "cameraman300"; // "cat490"; "peppers"; "mandrill";
    String ImageFile = "372019"; //"302003"; //"15088"; //"86016"; //"302008"; //"cameraman300";
    //String ImageFile = "butterfly-wings"; //"concentric-arcs"; //"broken-circle"; //"power-button"; "human-face"; "tooth-contour";
    //String FileExt = "txt";
    String FileExt = "jpg";
    //String FileExt = "png";
    Mat InputImage, CannyEdges;
    Graph GivenGraph;
    int CannyLowThreshold = 75;
    double CannyRatio = 3;
    String sCanny = OutputFolder;
    int DiagGap = 1;
    int VerGap = 1;
    String sGaps = "d"+to_string(DiagGap)+"v"+to_string(VerGap);
    int SizeGraph = 3;
    double NoiseLevel = 0.05, PointRatio = 0.2;
    String sNoise = "n" + to_string( (int)(NoiseLevel*100) );
    String sGraph = OutputFolder + GraphType + to_string(SizeGraph);
    int SizeCloud = 1000;
    vector<String> file_names = {
        "sc99.phi35.Er45.0Rr6.0",
        "sc98.phi35.Er37.5Rr5.0",
        "sc38.phi35.Er3.75Rr0.5",
        "sc29.phi35.Er9.375Rr1.25",
        "sc28.phi35.Er13.125Rr1.75",
        "sc11.phi35.Er22.5Rr3.0",
        "sc2.phi35.Er60.0Rr8.0",
        "sc1.phi35.Er7.5Rr1.0",
        "s30.phi35.Er11.5Rr1.5",
        "s29.phi35.Er15.0Rr2.0",
        "s28.phi35.Er30.0Rr4.0" };
    String add_ext = "_t_10000000.c";
    
    // Load PointCloud
    vector<P2> PointCloud;
    vector< vector<P2> > clouds;
    if (InputType == "image")
    {
        InputImage = imread( InputFolder + ImageFile + "." + FileExt);
        if( !InputImage.data ) { cout<<"\nFile not found:"<<InputFolder + ImageFile+"."+FileExt; return -1; }
        CannyDetector( InputImage, CannyLowThreshold, CannyRatio, CannyEdges);
        sCanny += ImageFile + "T"+to_string(CannyLowThreshold)+"R"+to_string( (int)(CannyRatio*10) );
        WriteImage( sCanny + "Canny.png", CannyEdges );
        LoadCloud( CannyEdges, PointCloud );
    }
    else if (InputType == "binary")
    {
        InputImage = imread( InputFolder + ImageFile + "." + FileExt);
        if( !InputImage.data ) { cout<<"\nFile not found:"<<InputFolder + ImageFile+"."+FileExt; return -1; }
        //cout<<" sizes: rows="<<InputImage.rows<<" cols="<<InputImage.cols;
        cvtColor( InputImage, InputImage, CV_RGB2GRAY );
        BinaryToCloud( InputImage, NoiseLevel, PointRatio, PointCloud );
        sCanny += ImageFile;
    }
    else if (InputType == "graph")
    {
        GivenGraph = Graph (GraphType, SizeGraph );
        GivenGraph.LoadCloud( SizeCloud, NoiseLevel, PointCloud );
    }
    else if (InputType == "file")
    {
        InputFolder = "/Users/kurlin/C++/Datasets/CloudFiles/";
        LoadCloud( SizeCloud, InputFolder + InputFile + ".txt", PointCloud );
        sCanny += InputFile;
    }
    else if (InputType == "Glasgow")
    {
        cout<<"Load "<<file_names.size()<<" clouds";
        load_clouds( InputFolder, file_names, add_ext, clouds );
        PointCloud = clouds[0];
        sCanny += file_names[0];
    }
    else
    { // InputType=="manual"
        PointCloud.push_back( P2(-1,4) );
        PointCloud.push_back( P2(2,4) );
        PointCloud.push_back( P2(-3,2) );
        PointCloud.push_back( P2(3,2) );
        PointCloud.push_back( P2(-2,0) );
        PointCloud.push_back( P2(2,0) );
        PointCloud.push_back( P2(-3,-2) );
        PointCloud.push_back( P2(3,-2) );
        PointCloud.push_back( P2(-2,-4) );
        PointCloud.push_back( P2(1,-4) );
    }
    SizeCloud = (int)PointCloud.size();
    if (SizeCloud<9) { cout<<"Tiny cloud of size "<<PointCloud.size(); return 1; }
    String sCloud = sCanny + "C" + to_string(SizeCloud);
    pair<Point2d, Point2d> CloudBox;
    FindCloudBox( PointCloud, CloudBox );
    //cout<<"\nCloudBox: first = "<<CloudBox.first<<" second = "<<CloudBox.second;
    Point2d SizesCloud( CloudBox.second.x - CloudBox.first.x, CloudBox.second.y - CloudBox.first.y );
    double MaxSizeCloud = std::max( SizesCloud.x, SizesCloud.y );
    double Scale = (SizeImage-2*BorderImage) / MaxSizeCloud;
    Point2d Shift = - Point2d( Scale * CloudBox.first.x, Scale * CloudBox.first.y );
    Shift += ShiftImage;
    Mat Image( Point2i( SizeImage, SizeImage ), CV_8UC3, White );
    if (InputType=="graph")
    {
        //GivenGraph.Print();
        Mat iGraph = Image.clone();
        GivenGraph.Draw(Scale, Shift, Green, iGraph);
        WriteImage( sGraph+".png", iGraph );
        sCloud = sGraph + sNoise + "C" + to_string(SizeCloud);
    }
    Mat iCloud = Image.clone();
    int MaxIndex = (int)PointCloud.size();
    draw_cloud( PointCloud, MaxIndex, Scale, Shift, iCloud );
    WriteImage( sCloud + ".png", iCloud );
    //return 4;

    // Build Delaunay triangulation Del
    DT Del;
    Del.insert( PointCloud.begin(), PointCloud.end() );
    int nVertices = (int)Del.number_of_vertices(); // finite
    size_t nFaces = Del.number_of_faces(); // finite
    int nEdges = int(nVertices + nFaces - 1);
    cout<<"Triangulation:"<<" #vertices = "<<nVertices<<" #edges = "<<nEdges<<" #faces = "<<nFaces;
    
    //Update info at vertices and faces: parent=itself, birth = circumradius for acute triangles
    int ind = 0, MaxInd = 1, nAcute = 0;
    for (VI iV = Del.vertices_begin(); iV != Del.vertices_end(); ++iV, ind++) iV->info() = ind;
    vector<FH> FaceHandles (nFaces+1);
    double radius = 0, MaxRadius = 0;
    ind = 0;
    for(FFI iF = Del.finite_faces_begin(); iF != Del.finite_faces_end(); ++iF)
    {
        ind++;
        FaceHandles[ind] = iF;
        FaceHandles[ind]->info().index = ind; // set index
        FaceHandles[ind]->info().uplink = ind; // parent = itself
        FaceHandles[ind]->info().live.resize(1);
        FaceHandles[ind]->info().live[0] = ind; // the only live triangle is the current
        FaceHandles[ind]->info().birth = 0;
        K::Triangle_2 Triangle (iF->vertex(0)->point(), iF->vertex(1)->point(), iF->vertex(2)->point() );
        if ( !is_acute( Triangle, radius ) ) continue; // non-acute triangle
        FaceHandles[ind]->info().birth = radius;
        if (radius > MaxRadius) { MaxRadius = radius; MaxInd = ind; }
        nAcute++;
    }
    cout<<" MaxRadius="<<MaxRadius; //<<" MaxInd="<<MaxInd; //<<" LargestTriangle: "<<FaceHandles[MaxInd]->vertex(0)->point().x();
    FaceHandles[0] = Del.infinite_face(); // default handle isn't used for external region
    FaceHandles[0]->info().birth = MaxRadius+1e-2; // external region is born at alpha = MaxRadius
    FaceHandles[0]->info().index = 0;
    FaceHandles[0]->info().uplink = 0;
    FaceHandles[0]->info().live.resize(1);
    FaceHandles[0]->info().live[0] = 0;
    
    // Build vector DelEdges and draw a triangulation
    vector<Edge> DelEdges;
    for(FEI iE = Del.finite_edges_begin(); iE != Del.finite_edges_end(); ++iE)
        DelEdges.push_back( Edge (*iE) );
    Mat iDelaunay = Image.clone();
    DrawTriangulation( Del, Scale, Shift, iDelaunay );
    WriteImage( sCloud + "DelT.png", iDelaunay );

    // Cycle over edges to find all critical values
    sort( DelEdges.begin(), DelEdges.end(), DecreasingLengths );
    double alpha=0;
    vector<Life> Persistence;
    vector<Region> Map;
    vector<int> NegativeEdges; NegativeEdges.clear();
    vector<int> PositiveEdges; PositiveEdges.clear();
    int iEdge=0, Links=0, nCases1=0, nCases2=0, nCases3=0, nCases4=0;
    for (; iEdge<DelEdges.size(); iEdge++)
    {
        Edge* CurEdge = &DelEdges[iEdge];
        alpha = 0.5 * CurEdge->length;
        int node_u = CurEdge->face1;
        int node_v = CurEdge->face2;
        int root_u = FindRoot (node_u, FaceHandles);
        int root_v = FindRoot (node_v, FaceHandles);
        //cout<<"\nE="<<iEdge<<" u"<<node_u<<" v"<<node_v<<" ru"<<root_u<<" rv"<<root_v;
        if (root_u == root_v) { nCases1++; NegativeEdges.push_back(iEdge); continue; } // Case 1
        Node* Node_u = &FaceHandles[node_u]->info();
        Node* Node_v = &FaceHandles[node_v]->info();
        Node_u->uplink = root_u; // shorten path to root
        Node_v->uplink = root_v; // shorten path to root
        Node* Root_u = &FaceHandles[root_u]->info();
        Node* Root_v = &FaceHandles[root_v]->info();
        double alpha_u = Root_u->birth; //FaceHandles[node_u]->info().birth;
        double alpha_v = Root_v->birth;
        if ( (alpha_u==0) && (alpha_v>0) ) // Case 2: node u joins the region of node v
        {
            nCases2++; //cout<<" C2: ";
            if ( node_u != root_u ) cout<<"\n Error2: "<<" u"<<node_u<<" v"<<node_v<<" ru"<<root_u<<" rv"<<root_v;
            Node_u->uplink = root_v; // root_v becomes parent of root_u
            Node_u->birth = alpha_v;
            if ( Root_v->height == 1 ) Root_v->height = 2; // node_u added to the tree at root_v
            int r = Node_v->bar; // index of the region containing the triangle T_v
            if ( r >= 0) // the region is already dead
            {
                Node_u->bar = r; // node_u joins the dead region of node_v
                Map[r].core.push_back(node_u);
            }
            else Root_v->live.push_back(node_u); // node_u becomes a live triangle at root_v
            Node_u->live.clear();
        }
        else if ( (alpha_u>0) && (alpha_v==0) ) // symmetric Case 2: node_v joins the region of node_u
        {
            nCases2++; //cout<<" C2': ";
            if ( node_v != root_v ) cout<<"\n Error2: "<<" u"<<node_u<<" v"<<node_v<<" ru"<<root_u<<" rv"<<root_v;
            Node_v->uplink = root_u;
            Node_v->birth = alpha_u;
            if ( Root_u->height == 1 ) Root_u->height = 2; // node_v added to the tree at root_u
            int r = Node_u->bar; // index the region containing the triangle T_u
            if ( r >= 0) // the region is already dead
            {
                Node_v->bar = r; // node_v joins the dead region of node_u
                Map[r].core.push_back(node_v);
            }
            else Root_u->live.push_back(node_v); // node_v becomes a live triangle at root_u
            Node_v->live.clear();
        }
        else if ( (alpha_u==0) && (alpha_v==0) ) // Case 3: two right-angled triangles at node_u, node_v
        {
            nCases3++; //cout<<" C3: ";
            if ( node_u != root_u || node_v != root_v )
                cout<<"\n E3: "<<" u"<<node_u<<" v"<<node_v<<" ru"<<root_u<<" rv"<<root_v;
            Node_v->uplink = node_u; // root_u = node_u becomes parent of root_v = node_v
            if (Node_u->height == 1) Node_u->height = 2; // node_v is now attached to its parent node_u
            Node_u->birth = alpha;
            Node_v->birth = alpha;
            Node_u->live.push_back(node_v);
            Node_v->live.clear();
        }
        else if ( (alpha_u>0) && (alpha_v>0) ) // Case 4: merge the regions of node_u, node_v
        {
            nCases4++; //cout<<" C4: ";
            double death = std::min(alpha_u,alpha_v);
            int heir = -1, supr = -1, ind = (int)Map.size(); // index of the newly added region
            Persistence.push_back( Life( ind, alpha, death, iEdge ) ); // younger class dies and is recorded
            // Create a new region (in Map) corresponding to a dot in the persistence diagram
            if (Root_u->birth <= Root_v->birth) // u is younger than v
            {   // set bar index for all live triangles at root_u
                for (int i=0; i<Root_u->live.size(); i++) FaceHandles[ Root_u->live[i] ]->info().bar = ind;
                if (Node_v->bar < 0) heir = node_v; // node_v is still alive and will be the heir
                else supr = Node_v->bar; // bar index of the already dead region containing T_v
                Map.push_back( Region( ind, alpha, death, Root_u->live, heir, supr, iEdge ) );
                if ( Root_u->height > Root_v->height )
                {   // u is younger, but higher than v
                    Root_v->uplink = root_u;
                    Root_u->live = Root_v->live; // triangles from Root_v->live are saved to the higher Root_u->live
                    Root_u->birth = Root_v->birth; // the new root now has the earliest birth
                }
                else
                {   // u is younger and lower than v
                    Root_u->uplink = root_v;
                    if ( Root_u->height == Root_v->height ) Root_v->height ++;
                }
            }
            else
            {   // v is younger (than u) and dies
                for (int i=0; i<Root_v->live.size(); i++) FaceHandles[ Root_v->live[i] ]->info().bar = ind;
                if (Node_u->bar < 0) heir = node_u; // root_u still alive and will be the heir
                else supr = Node_u->bar; // bar index of the already dead region containing T_u
                Map.push_back( Region( ind, alpha, death, Root_v->live, heir, supr, iEdge ) );
                if ( Root_u->height >= Root_v->height )
                {   // v is younger and lower than v
                    Root_v->uplink = root_u;
                    if ( Root_u->height == Root_v->height ) Root_u->height ++;
                }
                else
                {   // v is younger, but higher than u
                    Root_u->uplink = root_v;
                    Root_v->live = Root_u->live;
                    Root_v->birth = Root_u->birth;
                }
            }
        }
        Links++;
        if ( Links >= nFaces ) break;
    }
    cout<<"\nCases1="<<nCases1<<" nCases2="<<nCases2<<" nCases3="<<nCases3<<" nCases4="<<nCases4<<" C2+C3+C4="<<nCases2+nCases3+nCases4;// =nFaces
    for (int i=iEdge+1; i<DelEdges.size();i++) NegativeEdges.push_back( i );
    cout<<"\nMST has "<<NegativeEdges.size()<<" edges";
    if ( NegativeEdges.size() != nVertices-1 ) { cout<<"\nError with NegativeEdges.\n"; return -1; }
    
    // Add the final unbounded region to Map
    ind = (int)Map.size();
    int root = FindRoot (0, FaceHandles); //cout<<" root="<<root;
    Node* N = &FaceHandles[root]->info();
    for (int i=0; i < N->live.size(); i++) FaceHandles[ N->live[i] ]->info().bar = ind;
    Map.push_back( Region( ind, alpha, N->birth, N->live, -1, ind, -1 ) ); // last heir irrelevant
    
    //Output intermediate results
    cout<<" Acute="<<nAcute<<" nRegions="<<Map.size();
    int SumCores = 0;
    for (int j=0; j<Map.size()-1; j++) SumCores += Map[j].core.size();
    int ExtCore = (int)Map[ Map.size()-1 ].core.size();
    if (SumCores+ExtCore != nFaces+1 ) { cout<<"\nError with cores of initial regions.\n"; }
    int nDots = (int)Persistence.size();
    cout<<"\nPersistence diagram has nDots = "<<nDots; //, MaxDeath="<<MaxDeath;
    
    // Find k-th widest diagonal gap in Persistence, k = DiagGap
    sort( Persistence.begin(), Persistence.end(), IncreasingSpans );
    vector<IndexValue> DiagonalGaps (nDots);
    DiagonalGaps[0] = IndexValue( 0, Persistence[0].span ); //smallest gap
    for(int i=1; i<nDots; i++)
        DiagonalGaps[i] = IndexValue( i, Persistence[i].span - Persistence[i-1].span );
    sort( DiagonalGaps.begin(), DiagonalGaps.end(), DecreasingValues );
    vector<int> IndicesAbove( nDots );
    IndicesAbove[0] = DiagonalGaps[0].index; // index of dots above the 1st widest diagonal gap
    for(size_t i=1; i<nDots; i++)
        IndicesAbove[i] = std::min( DiagonalGaps[i].index, IndicesAbove[i-1] );
    for(size_t i=0; i<DiagGap; i++) cout<<" A["<<i+1<<"]="<<IndicesAbove[i];
    int nDiagSub = nDots - IndicesAbove[DiagGap-1];
    cout<<"\nThe cloud has nDiagSub = "<<nDiagSub<<" persistent holes";
    
    // Create the diagonal subdiagram above dgap
    vector<Life> DiagSub;
    DiagSub.insert( DiagSub.begin(), Persistence.begin()+IndicesAbove[DiagGap-1], Persistence.end());
    cout<<"\nDiagonal subdiagram above gap"<<DiagGap<<" has "<<nDiagSub<<" dots ";
    if (nDiagSub==1 || VerGap == 1) // take the 1st widest vertical gap to the right of all points above dgap
        for (size_t i = IndicesAbove[DiagGap-1]; i<nDots; i++)
            PositiveEdges.push_back( Persistence[i].edge );
    else
    {   // Find l-th widest vertical gap in DiagSub, l = VerGap
        sort( DiagSub.begin(), DiagSub.end(), IncreasingBirths );
        vector<IndexValue> VerticalGaps (nDiagSub-1);
        for(int i=0; i<nDiagSub-1; i++)
            VerticalGaps[i] = IndexValue ( i, DiagSub[i+1].birth - DiagSub[i].birth );
        sort( VerticalGaps.begin(), VerticalGaps.end(), DecreasingValues );
        vector<int> IndicesOnLeft(nDiagSub);
        IndicesOnLeft[0] = VerticalGaps[0].index; // index of dot on the left bound of vgap
        for(size_t i=1; i<nDiagSub-1; i++)
            IndicesOnLeft[i] = std::min( VerticalGaps[i].index, IndicesOnLeft[i-1] );
        for(size_t i=0; i<nDiagSub-1; i++) cout<<" L["<<i+2<<"]="<<IndicesOnLeft[i];
        cout<<"\nVertical subdiagram to the left of vgap"<<VerGap<<" has "<<IndicesOnLeft[VerGap-2]+1<<" dots";
        for (size_t i = 0; i<=IndicesOnLeft[VerGap-2]; i++)
            PositiveEdges.push_back( DiagSub[i].edge );
    }
    
    // Find the doubled scale 2vs_{kl} equal to the length of a longest positive edge
    double CritLength = 0;
    for (size_t i=0; i<PositiveEdges.size(); i++)
        if ( DelEdges[ PositiveEdges[i] ].length > CritLength )
            CritLength = DelEdges[ PositiveEdges[i] ].length;
    //cout<<"\nCritLength = "<<Round(CritLength);

    // Draw diagram
    Mat iDiagram = Image.clone();
    double ScaleDiagram = (SizeImage-2.0*BorderImage) / MaxRadius;
    DrawDiagram( Persistence, IndicesAbove[DiagGap-1], CritLength, MaxRadius, ScaleDiagram, ShiftImage, iDiagram );
    WriteImage( sCloud + "P" + sGaps + ".png", iDiagram );
    
    /* Draw HoPeS
    Mat iHopes = Image.clone();
    DrawHopes (DelEdges, NegativeEdges, PositiveEdges, Scale, Shift, iHopes);
    WriteImage( sCloud + "HoPeS"+sGaps+".png", iHopes );
    int nSimEdges = (nVertices-1) + (int)PositiveEdges.size();
    cout<<"\nHoPeS has "<<nSimEdges<<" edges";
    
    // Build HoPeS
    vector<int> EdgeIndices;
    EdgeIndices.insert ( EdgeIndices.begin(), NegativeEdges.begin(), NegativeEdges.end() );
    Graph MST (DelEdges, EdgeIndices, nVertices);
    Mat iMST = Image.clone();
    MST.Draw (Scale, Shift, iMST);
    WriteImage( sCloud + "MST.png", iMST );
    EdgeIndices.insert ( EdgeIndices.begin(), PositiveEdges.begin(), PositiveEdges.end() );
    Graph Hopes (DelEdges, EdgeIndices, nVertices);
    Graph SimHopes = Hopes;

    // Simplify HoPeS
    while ( SimHopes.CollapseShortestEdge( 0.5 * CritLength ) ) nSimEdges--;
    cout<<" SimHoPeS has "<<nSimEdges<<" edges";
    Mat iSimHopes = Image.clone();
    SimHopes.Draw (Scale, Shift, iSimHopes);
    WriteImage( sCloud + "HoPeS"+sGaps+"sim"+to_string(nSimEdges)+".png", iSimHopes );
    
    while ( SimHopes.CollapseShortestEdge( CritLength ) ) nSimEdges--;
    cout<<" SimHoPeS2 has "<<nSimEdges<<" edges";
    Mat iSimHopes2 = Image.clone();
    SimHopes.Draw (Scale, Shift, iSimHopes2);
    WriteImage( sCloud + "HoPeS"+sGaps+"sim"+to_string(nSimEdges)+".png", iSimHopes2 );
   
    // Build Staircase
    vector<IndexValue> Staircase;
    BuildStaircase ( Persistence, IndicesAbove[DiagGap-1], Staircase );
    Mat iStaircase = Image.clone();
    DrawStaircase ( Staircase, ShiftImage, SizeImage-2*BorderImage, iStaircase );
    WriteImage( sCloud + "Staircase.png", iStaircase );
  
    // Compute probabilities for the number of holes
    int MaxHoles = 0;
    for (size_t i=0; i<Staircase.size(); i++)
        if ( Staircase[i].index > MaxHoles ) MaxHoles = Staircase[i].index;
    vector<IndexValue> Probabilities (MaxHoles+1);
    for (int i=0; i<=MaxHoles; i++)
    {
        Probabilities[i].index = i;
        Probabilities[i].value = 0;
    }
    for (size_t i=0; i<Staircase.size()-1; i++)
    {
        Probabilities[ Staircase[i].index ].value += Staircase[i+1].value - Staircase[i].value;
        //cout<<" P["<<Staircase[i].index<<"]+="<<Staircase[i+1].value - Staircase[i].value;
    }
    cout<<"\nLet the scale be random and uniform over [ "
        <<Staircase[0].value<<", "<<Staircase[ Staircase.size()-1 ].value<<" ].\nThen";
    for (size_t i=0; i<=MaxHoles; i++)
        Probabilities[i].value /= Staircase[ Staircase.size()-1 ].value - Staircase[0].value;
    sort( Probabilities.begin(), Probabilities.end(), DecreasingValues );
    for (size_t i=0; i<=std::min(10, MaxHoles); i++)
        cout<<" Prob("<<Probabilities[i].index<<" holes)="<<Probabilities[i].value;
     */
    // Draw the Map of initial regions
    Mat iMap = Image.clone();
    DrawMap( FaceHandles, Map, Scale, Shift, iMap );
    WriteImage( sCloud + "InitRegions" + to_string( Map.size()-1 )+".png", iMap );
    
    // Sort Map by life spans and build correspondce OldIndex -> NewIndex
    //sort( Map.begin(), Map.end(), DecreasingSpans );
    //for (int f=0; f<=nFaces; f++) cout<<" FH["<<f<<"].bar="<<FaceHandles[f]->info().bar; Print( Persistence );
    // Set supr (indices of superior components in Map)
    for (int k=0; k<=nDots; k++) if (Map[k].supr < 0)
    {
        if ( Map[k].heir >= 0 )
            Map[k].supr = FaceHandles[ Map[k].heir ]->info().bar; // old index of superior region
        else cout<<"\n*error1 with heir and supr: k="<<k;
        //cout<<"\n supr<0 k="<<k<<" now: supr="<<Map[k].supr;
        //if ( Map[ Map[k].supr ].span <= Map[k].span ) cout<<"\n! superior less persistent: k="<<k<<" supr="<<Map[k].supr;
    }
    //    Map[k].supr = NewIndices[ FaceHandles[ Map[k].heir ]->info().bar ];
    //cout<<"\nMap after setting supr"; Print( Map );
    int SumRegions = 0;
    for (int k=nDots-nDiagSub; k<=nDots; k++)
    {
        int old_ind = nDots;
        if (k<nDots) old_ind = Persistence[k].index;
        int s = (int)Map[ old_ind ].core.size();
        //cout<<"\nPersistent region old_ind="<<old_ind<<" includes "<<s<<" initial regions";
        SumRegions += s;
    }
    // Update supr
    vector<int> NewIndices( nDots+1 );
    for (int k=0; k<nDots; k++) NewIndices[ Persistence[k].index ] = k;
    NewIndices[ nDots ] = nDots;
    for (int k=0; k<nDots-nDiagSub; k++)
    {
        int old_ind = Persistence[k].index;
        int s = Map[old_ind].supr;
        while( NewIndices[s] < nDots-nDiagSub && s != Map[s].supr )
        {
            s = Map[s].supr; //cout<<"?";
        }
        Map[old_ind].supr = s;
    }
    //cout<<"\nMap after updating supr"; Print( Map );
    for (int k=0; k<nDots-nDiagSub; k++) // go over only low persistent regions ordered by increasing spans
    {
        int old_ind = Persistence[k].index; // old index of a low persistent initial region
        int superior = Map[old_ind].supr;
        if ( superior<0 ) cout<<"\n*error2 with heir and supr";
        //cout<<"\n k="<<k<<" old_ind="<<old_ind<<" superior="<<superior;
        Map[superior].core.insert( Map[superior].core.end(), Map[old_ind].core.begin(), Map[old_ind].core.end() );
        Map[old_ind].core.clear();
    }
    //Print(Map);
    SumRegions = 0;
    for (int k=nDots-nDiagSub; k<=nDots; k++)
    {
        int old_ind = nDots;
        if (k<nDots) old_ind = Persistence[k].index;
        int s = (int)Map[ old_ind ].core.size();
        //cout<<"\nFinal region old_ind="<<old_ind<<" includes "<<s<<" initial regions";
        SumRegions += s;
    }
    cout<<"\nTotal sum of initial regions ater merging = "<<SumRegions;
    DrawMap( FaceHandles, Map, Scale, Shift, iMap );
    FindBoundary( (int)nFaces+1, Map, DelEdges );
    DrawBoundary( DelEdges, Scale, Shift, iMap );
    WriteImage( sCloud + "FinalRegions" + to_string(nDiagSub)+".png", iMap );
    return 0;
}