#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 CvScalar White = CV_RGB (255, 255, 255);
const CvScalar Black = CV_RGB (0, 0, 0);
const CvScalar Red = CV_RGB (255, 0, 0);
const CvScalar Lime = CV_RGB (0, 255, 0);
const CvScalar Green = CV_RGB (0, 128, 0);
const CvScalar Blue = CV_RGB (0, 0, 255);
const CvScalar Yellow = CV_RGB (255, 255, 0);
const CvScalar Orange = CV_RGB (255, 128, 0);
const CvScalar Magenta = CV_RGB (255, 0, 255);
const CvScalar Cyan = CV_RGB (0, 255, 255);
const CvScalar Olive = CV_RGB (128,128,0);
const CvScalar Gray = CV_RGB (128,128,128);
const CvScalar Silver = CV_RGB (192,192,192);
const CvScalar Maroon = CV_RGB (128,0,0);
const CvScalar Purple = CV_RGB (128,0,128);
const CvScalar Teal = CV_RGB (0,128,128);
const CvScalar Navy = CV_RGB (0,0,128);
vector<CvScalar> Colors = { Cyan, Yellow, Lime, Orange, Red, Blue, Green, Gray, Olive, Cyan, Magenta, Silver, Maroon, Purple, Teal, Navy, Black };
#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;
    
    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) );
        //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 (String type, int SizeCloud, String const& FileName, vector<P2>& PointCloud)
{
    if ( type=="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) );
        return;
    }
    if ( type=="file")
    {
        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) );
        }
    }
}

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 DrawPointCloud ( 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), 1, 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, "%.4f", 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 r2.span <= r1.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 )
{
    for (int j=0; j<Map.size(); j++)
    {
        CvScalar c = Colors[ j % Colors.size() ];
        if ( Map[j].index == Map.size()-1 ) c = White; // external region
        for (size_t i=0; i<Map[j].core.size(); i++)
        {
            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 );
            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 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, 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);
    int iRight = -1;
    double c = 0.5*CritLength, 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 = 3;
        if (Persistence[i].birth > c) size = 2;
        if (i<iAbove) size = 1;
        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 LoadPointCloud (Mat const& ImageCloud, vector<P2>& PointCloud)
{
    size_t cols = ImageCloud.cols;
    size_t rows = ImageCloud.rows;
    for (int j=0; j<rows; j++)
        for (int i=0; i<cols; i++)
            if ( ImageCloud.at<uchar>(j, i) != 0 )
                PointCloud.push_back( P2(i,j) );
}

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<Region>const& Map)
{
    for (int j=0; j<Map.size(); j++)
    {
        cout<<" Map["<<j<<"]: b="<<Round(Map[j].birth)<<" d="<<Round(Map[j].death)<<" h="<<Map[j].heir<<" s="<<Map[j].supr<<" core:";
        for (int i=0; i<Map[j].core.size(); i++) cout<<" "<<Map[j].core[i]; cout<<endl;
    }
}

class Graph
{
public:
    vector<P2> vertices;
    vector< pair<int,int> > edges;
    vector< vector<int> > vedges; // indices of edges around a vertex
    Graph() {}
    Graph( vector<Edge>const& DelEdges, vector<int>const& EdgeIndices, int nVertices )
    {
        edges.resize( EdgeIndices.size() );
        vertices.resize (nVertices);
        vedges.resize (nVertices);
        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();
            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 );
        }
    }
    
    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]<<")";
    }

    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);
        }
    }
    //
    
    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()
{
    int SizeImage = 532;
    int BorderImage = 16;
    Point2d ShiftImage ( BorderImage, BorderImage );
    String InputFolder = "/Users/kurlin/C++/Datasets/IconicImages/";
    String OutputFolder = "/Users/kurlin/C++/Output/AlphaComplexes2d/";
    String InputFile = "figure-eight"; // "cloud-ears";
    String ImageFile = "cameraman300"; // "cat490"; "peppers"; "mandrill";
    String FileExt = "txt"; // "jpg"; "png";
    Mat InputImage, CannyEdges;
    String type = "manual"; // "file"; "image";
    int CannyLowThreshold = 150;
    double CannyRatio = 3;
    String sCanny = OutputFolder;
    int DiagGap = 1;
    int VerGap = 1;
    String sGaps = "d"+to_string(DiagGap)+"v"+to_string(VerGap);
    
    // Load PointCloud
    int SizeCloud = 10;
    vector<P2> PointCloud;
    if (type == "image")
    {
        InputImage = imread( InputFolder + ImageFile + "." + FileExt);
        if( !InputImage.data ) { cout<<"\nFile not found:"<<InputFolder + ImageFile; return -1; }
        CannyDetector( InputImage, CannyLowThreshold, CannyRatio, CannyEdges);
        sCanny += ImageFile + "T"+to_string(CannyLowThreshold)+"R"+to_string((int)(CannyRatio*10));
        WriteImage( sCanny + "Canny.png", CannyEdges );
        LoadPointCloud ( CannyEdges, PointCloud );
    }
    else
    {
        InputFolder = "/Users/kurlin/C++/Datasets/AlphaComplexes2d/";
        LoadCloud ( type, SizeCloud, InputFolder + InputFile + ".txt", PointCloud );
        sCanny += InputFile;
    }
    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 );
    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);
    Mat iCloud = Image.clone();
    int MaxIndex = (int)PointCloud.size();
    DrawPointCloud ( PointCloud, MaxIndex, Scale, Shift, iCloud );
    WriteImage( sCloud + "Cloud.png", iCloud );

    // 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<<"\nMaxRadius="<<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;
    //cout<<"\n Acute triangles: ";
    
    // 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);
        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++;
            if ( node_u != root_u ) cout<<"\n E2: "<<" 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++;
            if ( node_v != root_v ) cout<<"\n E2: "<<" 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++;
            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 (Root_v->bar < 0) heir = root_v; // root_u still alive and will be the heir
                else supr = Root_v->bar; // bar index of the already dead region containing T_u
                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 new 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 (Root_u->bar < 0) heir = root_u; // root_u still alive and will be the heir
                else supr = Root_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<<"\n 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<<" MST 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, -1, -1) ); // last heir irrelevant
    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"; }
    
    // Find MaxDeath and ScaleDiagram
    size_t nDots = Persistence.size();
    double ScaleDiagram = (SizeImage-2.0*BorderImage) / MaxRadius;
    cout<<"\nPersistence diagram has "<<nDots<<" dots"; //, 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];
    
    // Create the diagonal subdiagram above dgap
    vector<Life> DiagSub;
    DiagSub.insert( DiagSub.begin(), Persistence.begin()+IndicesAbove[DiagGap-1], Persistence.end());
    int nDiagSub = (int)DiagSub.size(); //(int)nDots-AboveLowestOfWidest[DiagGap-1];
    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;
    int nSimEdges = (nVertices-1) + (int)PositiveEdges.size();
    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();
    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 );
    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( CritLength ) ) nSimEdges--;
    cout<<"\nSimHoPeS 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( 20*CritLength ) ) nSimEdges--;
    cout<<"\nSimHoPeS2 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 is 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;

    /* to be completed by late June 2015: Draw the Map of initial regions
    Mat iMap = Image.clone();
    DrawMap ( FaceHandles, Map, Scale, Shift, iMap );
    WriteImage( sCloud + "Regions" + to_string(nDiagSub)+".png", iMap );
     */
    return 0;
}