/*
    Hyperbolic Tessellations
    Copyright (C) 2007 Dmitry Brant
    http://dmitrybrant.com

    Based largely on code from
    http://www.plunk.org/~hatch/HyperbolicApplet/
*/
//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
//---------------------------------------------------------------------------

#include "hyperbolicdraw.h"

void DrawLine2fPoincare(MyGraphics* mg,
                   complex& x, complex& y,
                   double tolerance,
                   TessellationSettings* settings)
{

    if (settings->doRandomJitter){
        double jit = .02;
        x[0] += (2*(double(rand()%32)/32)-1.) * jit;
        x[1] += (2*(double(rand()%32)/32)-1.) * jit;
        y[0] += (2*(double(rand()%32)/32)-1.) * jit;
        y[1] += (2*(double(rand()%32)/32)-1.) * jit;
        double len0Sqrd = x[0]*x[0]+x[1]*x[1];
        double len1Sqrd = y[0]*y[0]+y[1]*y[1];
        if (len0Sqrd > .99){
            double scale = .99 / sqrt(len0Sqrd);
            x *= scale;
        }
        if (len1Sqrd > .99){
            double scale = .99 / sqrt(len1Sqrd);
            y *= scale;
        }
    }


    int savedColor;
    if(settings->fadeOutColors){
        savedColor = mg->getColor();
        complex temp(.4,0.);
        double dist = sqrt(x[0]*x[0] + x[1]*x[1]);
        dist = cos(dist * PI/2.);
        if(dist < 0.)
            dist = 0.;
        int r=(savedColor & 0xFF),g=(savedColor>>8)&0xFF,b=(savedColor>>16)&0xFF;
        r = int(double(r)*dist);
        g = int(double(g)*dist);
        b = int(double(b)*dist);
        mg->setColor((b<<16)+(g<<8)+r);
    }


    if (settings->hyperbolicModel == KLEIN_DISK){

        complex x_(x), y_(y);
        p2k(x_); p2k(y_);
        mg->drawLine(x_[0],x_[1],y_[0],y_[1]);

    }else if (settings->hyperbolicModel == POINCARE_DISK){

        if (!settings->doCurveLines){

            mg->drawLine(x[0],x[1],y[0],y[1]);

        }else{

            if (settings->doAlternateDraw){
                int maxIntermediatePoints = 9;
                Vector<complex> intermediatePoints(maxIntermediatePoints+2);
                intermediatePoints[0] = x;
                intermediatePoints[maxIntermediatePoints+1] = y;
                int totalPoints = 2 + subdivideArc(intermediatePoints[0], intermediatePoints[maxIntermediatePoints+1], tolerance, maxIntermediatePoints, 1, intermediatePoints);
                intermediatePoints[totalPoints-1] = y;
                for(int i=0; i<totalPoints; i++){
                    mg->drawPoint(intermediatePoints[i][0], intermediatePoints[i][1],
                                 i==0||i==totalPoints-1 ? 3 : 1);
                }

            }else{

                if (((y[0]-x[0])*(y[0]-x[0])) + ((y[1]-x[1])*(y[1]-x[1])) > 1e-9) // XXX a bit ad-hoc, but don't need to draw tiny stuff and math is probably unstable
                {
                    // Need a third point on the circle,
                    // so we can compute its curvature and stuff.
                    complex mid;
                    hlerp2_poincare(x, y, .5, mid);

                    int sign; // -1 if CW, 1 if CCW
                    {
                        complex seg0 = mid - x;
                        complex seg1 = y - mid;
                        sign = vxv2(seg0,seg1) < 0. ? -1 : 1;
                    }

                    double curvature, start, end, focusAngleRadians;
                    {
                        // From http://mathworld.wolfram.com/Circumradius.html
                        // We want the inverse of the circumradius
                        // (so this should work even if the triangle is degenerate,
                        // as long as all three points are distinct,
                        // which they will be since p0 != p1 with some tolerance
                        double a = dist(x,mid);
                        double b = dist(mid,y);
                        double c = dist(x,y);
                        double s = (a+b+c)/2;

                        double temp = s*(a+b-s)*(a+c-s)*(b+c-s);
                        if (temp < 0.)
                            temp = 0.; // XXX ARGH! it galls me to have to do stuff like this, since it means we're using a bad formula
                        curvature = 4*sqrt(temp)/(a*b*c);

                        start = sign * - a * asin_over_x(a*.5*curvature);
                        end   = sign *   b * asin_over_x(b*.5*curvature);

                        //assert(ABS(end-start) < 2.01); // sanity check

                        double dirAtMid = atan2(mid[1]-x[1],mid[0]-x[0]);
                        dirAtMid += sign * asin(a*.5*curvature); // angle the arc makes with the segment from p0 to mid
                        focusAngleRadians = dirAtMid - sign * (PI/2);
                    }

                    mg->smartDrawArc(mid[0], mid[1], focusAngleRadians, curvature, start, end);
                }
            }
        }
    }
    else
    {

    }

    if(settings->fadeOutColors){
        mg->setColor(savedColor);
    }

} // DrawLine2fPoincare
//---------------------------------------------------------------------------





void DrawTiling(MyGraphics* mg,
           Vector<Isometry>& isometries,
           Vector<complex>& localVerts,
           int nSegments,
           Vector< Complex<int> >& segments,
           int** segmentColors, // can be null  [nsegments][3]
           bool floodFill,
           double tolerance,
           TessellationSettings* settings)
{

    int nIsometries = isometries.length;
    int nLocalVerts = localVerts.length;


    complex label0Position, labelPosition;
    if (settings->showVertexLabels){
        if (nIsometries >= 2){
            complex v0(isometries[0].P[0], isometries[0].P[1]);
            complex v1(isometries[1].P[0], isometries[1].P[1]);
            complex temp;
            double hDist = .25 * hlerp2_poincare(v0, v1, 0., temp);
            double eDist = h2eNorm(hDist);
            if ((eDist - 0.) > 1e-3){
                label0Position[0] = eDist;
                label0Position[1] = .3*eDist;
            }else{
                label0Position[0] = .08;
                label0Position[1] = .03; // arbitrary
            }
        }
        else{
            label0Position[0] = .08;
            label0Position[1] = .03; // arbitrary
        }
    }

    Vector<complex> verts(nLocalVerts);
    int i, j;

    for(i=0; i<nIsometries; i++)
    {
        for(j=0; j<nLocalVerts; j++)
            isometries[i].apply(localVerts[j], verts[j]);


        for(j=0; j<nSegments; j++){

            if (settings->doRemoveDups){
                // implement this!
            }

            DrawLine2fPoincare(mg,
                              verts[segments[j][0]],
                              verts[segments[j][1]],
                              tolerance,
                              settings);

            if(settings->showSegmentLabels){
                int savedColor = mg->getColor();
                mg->setColor(0xFF);
                complex v0(verts[segments[j][0]]);
                complex v1(verts[segments[j][1]]);
                hlerp2_poincare(v0, v1, .5, labelPosition);
                if(settings->hyperbolicModel == KLEIN_DISK) p2k(labelPosition);
                String label = String(i) + ":" + String(j);
                if (isometries[i].R == -1)
                    label = "["+label+"]";
                mg->drawStringCentered(label, labelPosition[0], labelPosition[1]);
                mg->setColor(savedColor);
            }

        }

        if(floodFill){
            if(!mg->isTrivial(verts[0], verts[1], verts[2])){
                complex floodFillPoint(0.,0.);
                if(settings->doCurveLines && (settings->hyperbolicModel != KLEIN_DISK)){
                    complex temp;
                    hlerp2_poincare(verts[0], verts[1], 0.5, temp);
                    hlerp2_poincare(verts[0], verts[2], 0.5, floodFillPoint);
                    hlerp2_poincare(temp, floodFillPoint, 0.5, floodFillPoint);
                }else{
                    floodFillPoint[0] += verts[0][0]; floodFillPoint[1] += verts[0][1];
                    floodFillPoint[0] += verts[1][0]; floodFillPoint[1] += verts[1][1];
                    floodFillPoint[0] += verts[2][0]; floodFillPoint[1] += verts[2][1];
                    floodFillPoint[0] /= 3.; floodFillPoint[1] /= 3.;
                }
                if(settings->hyperbolicModel == KLEIN_DISK) p2k(floodFillPoint);
                mg->floodFill(floodFillPoint[0], floodFillPoint[1]);

                int savedColor = mg->getColor();
                mg->setColor(0xFFFFFF);
                mg->drawPoint(floodFillPoint[0], floodFillPoint[1], 1);
                mg->setColor(savedColor);
            }
        }


        if (settings->showVertexLabels){
            int savedColor = mg->getColor();
            mg->setColor(0xFFFFFF);
            isometries[i].apply(label0Position, labelPosition);

            if(settings->hyperbolicModel == KLEIN_DISK){
                p2k(labelPosition);
            }

            String label = String(i);
            if (isometries[i].R == -1)
                label = "["+label+"]";
            mg->drawStringCentered(label, labelPosition[0], labelPosition[1]);
            mg->setColor(savedColor);
        }

    }
} // DrawTiling



void DrawOmnitruncatedTiling(MyGraphics* mg,
                        TessellationSettings* settings,
                        Isometry& F0,
                        int p, int q,
                        Vector<double>& wythoffCoeffs,
                        Isometry& return_smallestIsometry,
                        PickableSchwarzPolygon& return_pickableSchwarz)
{
    double tolerance = 1e-4;
    double halfEdgeLength = calcRegularTilingHalfEdgeLength(p, q);
    double R = h2eNorm(2*halfEdgeLength);
    double r = h2eNorm(halfEdgeLength);
    // find generators for the group.  B is rotation around the vertex,
    // A is rotation around the face center.
    Isometry A(complex(cos(PI-2*PI/q), sin(PI-2*PI/q)), complex(R,0), 1);
    Isometry B(complex(cos(2*PI/q), sin(2*PI/q)), zero, 1);

    Vector<Isometry> gens(2);
    gens[0] = A;
    gens[1] = B;

    double tooSmall = 1e-4; // XXX what does this mean?
    double tooFar = 1 - tooSmall;

    //double x0 = 0., y0 = 0., x1 = R, y1 = 0.;
    Vector<Isometry> list(settings->maxIsometries);
    int nIsometries;

    mg->beginDraw();
    try{

        // woops, enumerateForUniformTiling expects generators
        // to be different...
        // XXX clean this up
        {
            Isometry C; // rotation around center of edge 0
            {
                C = Isometry::pureTranslation(R,0) * Isometry::pureRotation(PI);
            }

            gens.Resize(q);
            for(int iGen=0; iGen<q; iGen++){
                Isometry rotate = Isometry::pow(B, iGen);
                Isometry temp = C * rotate.inverse();
                gens[iGen] = rotate * temp;
            }
        }

        Vector<int> pp(1);
        pp[0] = p;
        Vector<int> backEdgeInds(1);
        backEdgeInds[0] = ~0;

        Vector<int> CompositionLengths(settings->maxIsometries);

        nIsometries = UnCachedEnumerateIsometry2GroupForUniformTiling(
                            F0,
                            gens.length,
                            gens,
                            pp,
                            q,
                            backEdgeInds, // rotation everywhere
                            settings->maxIsometries,
                            list,
                            CompositionLengths, // don't care about composition lengths
                            return_smallestIsometry,
                            settings->maxLevels,
                            tooFar,
                            tooSmall);
        // now we have only 1 isometry per vertex.  we want q per vertex.

        {
            int newNIsometries = nIsometries * q;
            Vector<Isometry> newIsometries(newNIsometries);
            Vector<Isometry> rotates(q);
            for(int iArm=0; iArm<q; iArm++)
                rotates[iArm] = Isometry::pow(B, iArm);
            for(int iIsometry=0; iIsometry<nIsometries; iIsometry++){
                for(int iArm=0; iArm<q; iArm++)
                    newIsometries[iIsometry*q + iArm] = list[iIsometry] * rotates[iArm];
            }
            nIsometries = newNIsometries;
            list = newIsometries;
        }

        Vector<complex> elementCenters(3);     // vert, edge, face
        elementCenters[0][0] = 0.;
        elementCenters[0][1] = 0.;
        elementCenters[1][0] = r;
        elementCenters[1][1] = 0.;
        calcFaceCenter(p, q, elementCenters[2]);

        return_pickableSchwarz.setVertices(elementCenters, F0);


        if(settings->drawGrid){
            mg->setColor(0x8000);
            int subDiv = 20;
            for(int i=0; i<subDiv+1; i++){
                mg->drawLine(LERP(-1., 1., double(i)/double(subDiv)), -1.,
                            LERP(-1., 1., double(i)/double(subDiv)), 1.);
                mg->drawLine(-1., LERP(-1., 1., double(i)/double(subDiv)),
                            1., LERP(-1., 1., double(i)/double(subDiv)));
            }
        }


        if (settings->drawDual){
            mg->setColor(settings->dualColor);

            Vector<complex> P(3);
            P[0][0] = 0.; P[0][1] = 0.;
            P[1][0] = r;  P[1][1] = 0.;
            calcFaceCenter(p, q, P[2]);

            Vector<Complex <int> > segments(4);
            int nSegments = 0;
            if (wythoffCoeffs[0] != 0.){
                segments[nSegments][0] = 1;
                segments[nSegments][1] = 2;
                nSegments++;
            }
            if (wythoffCoeffs[1] != 0.)
            {
                segments[nSegments][0] = 0;
                segments[nSegments][1] = 2;
                nSegments++;
            }
            if (wythoffCoeffs[2] != 0.)
            {
                segments[nSegments][0] = 0;
                segments[nSegments][1] = 1;
                nSegments++;
            }

            DrawTiling(mg, list, P, nSegments, segments, NULL, /*nSegments > 2 ? true :*/ false, tolerance, settings);
        }

        if (settings->drawSnub){
            mg->setColor(settings->snubColor);

            Vector<complex> P(6);
            try{
            if (!hBary2(P[0], wythoffCoeffs, elementCenters))
                return; // XXX print warning?
            }catch(...){ }

            {
                B.apply(P[0], P[1]);
                B.inverse().apply(P[0], P[2]);
                A.apply(P[1], P[3]);
                A.apply(P[0], P[4]);
                A.inverse().apply(P[0], P[5]);
            }
            if (settings->snubParity)
            {
                for(int i=0; i<P.length; i++)
                    P[i][1] *= -1;
            }

            Vector< Complex<int> > segments(4);
            segments[0][0] = 0; segments[0][1] = 1;
            segments[1][0] = 0; segments[1][1] = 2;
            segments[2][0] = 0; segments[2][1] = 3;
            segments[3][0] = 0; segments[3][1] = 4;
            int nSegments = segments.length;
            DrawTiling(mg, list, P, nSegments, segments, NULL, false, tolerance, settings);
        }

        if (settings->drawPrimal){
            mg->setColor(settings->primalColor);

            Vector<complex> P(7);
            complex temp, P0;
            try{
            if (!hBary2(P0, wythoffCoeffs, elementCenters))
                return; // XXX print warning?
            }catch(...){ }
            P[0] = P0;
            for(int i=0; i<3; i++){
                calcClosestPointOnLine2(temp, P0, elementCenters[i], elementCenters[(i+1)%3]);
                P[i+1] = temp;

                hlerp2_poincare(P0, temp, 2., temp);
                P[i+1+3] = temp;
            }

            Vector< Complex<int> > segments(6);
            int nSegments = 0;

            for(int i=0; i<3; i++)
                if (wythoffCoeffs[(i+2)%3] != 0.)
                {
                    segments[nSegments][0] = 0;
                    segments[nSegments][1] = i+1;
                    nSegments++;
                    if (i==2 || wythoffCoeffs[2*i] != 0.) // XXX should try to make this intelligible
                    {
                        segments[nSegments][0] = i+1;
                        segments[nSegments][1] = i+1 + 3;
                        nSegments++;
                    }
                }

            DrawTiling(mg, list, P, nSegments, segments, NULL, false, tolerance, settings);
        }

        if (settings->drawTriangle1){
            mg->setColor(0xFFFFFF);

            Vector<complex> P(5);
            P[0][0] = 0.; P[0][1] = 0.;
            P[1][0] = r;  P[1][1] = 0.;
            calcFaceCenter(p, q, P[2]);
            hlerp2_poincare(P[0], P[2], 0.5, P[3]);
            A.apply(P[3], P[4]);


            Vector<Complex <int> > segments(10);
            int nSegments = 0;

            segments[nSegments][0] = 0;
            segments[nSegments++][1] = 3;

            segments[nSegments][0] = 0;
            segments[nSegments++][1] = 1;

            segments[nSegments][0] = 1;
            segments[nSegments++][1] = 3;

            segments[nSegments][0] = 3;
            segments[nSegments++][1] = 2;

            segments[nSegments][0] = 3;
            segments[nSegments++][1] = 4;

            segments[nSegments][0] = 1;
            segments[nSegments++][1] = 4;

            DrawTiling(mg, list, P, nSegments, segments, NULL, false, tolerance, settings);
        }


        if (settings->drawTriangle2){
            mg->setColor(0xFFFFFF);

            Vector<complex> P(15);
            complex fc;

            P[0][0] = 0.; P[0][1] = 0.;
            calcFaceCenter(p, q, fc);
            P[4] = fc;
            hlerp2_poincare(P[0], fc, 0.25, P[1]);
            hlerp2_poincare(P[0], fc, 0.5, P[2]);
            hlerp2_poincare(P[0], fc, 0.75, P[3]);

            A.apply(P[0], P[8]);
            hlerp2_poincare(P[8], fc, 0.25, P[7]);
            hlerp2_poincare(P[8], fc, 0.5, P[6]);
            hlerp2_poincare(P[8], fc, 0.75, P[5]);

            P[10][0] = r; P[10][1] = 0.;
            hlerp2_poincare(P[0], P[10], 0.5, P[11]);
            hlerp2_poincare(P[10], P[8], 0.5, P[9]);

            hlerp2_poincare(P[2], P[6], 0.5, P[14]);
            hlerp2_poincare(P[2], P[10], 0.5, P[12]);
            hlerp2_poincare(P[10], P[6], 0.5, P[13]);

            Vector<Complex <int> > segments(32);
            int nSegments = 0;

            segments[nSegments][0] = 3; segments[nSegments++][1] = 4;
            segments[nSegments][0] = 3; segments[nSegments++][1] = 5;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 3;
            segments[nSegments][0] = 3; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 5;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 6;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 2;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 6;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 7;
            segments[nSegments][0] = 0; segments[nSegments++][1] = 1;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 11;
            segments[nSegments][0] = 11; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 10;
            segments[nSegments][0] = 10; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 9;
            segments[nSegments][0] = 9; segments[nSegments++][1] = 7;
            segments[nSegments][0] = 0; segments[nSegments++][1] = 11;
            segments[nSegments][0] = 11; segments[nSegments++][1] = 10;

            DrawTiling(mg, list, P, nSegments, segments, NULL, false, tolerance, settings);
        }


        if (settings->drawExperimental){
            mg->setColor(0xFFFFFF);

            /*
            Vector<complex> P(12);
            complex fc;

            P[8][0] = r; P[8][1] = 0.;
            calcFaceCenter(p, q, fc);
            P[4] = fc;
            A.apply(P[8], P[0]);

            hlerp2_poincare(P[0], fc, 0.25, P[1]);
            hlerp2_poincare(P[0], fc, 0.5, P[2]);
            hlerp2_poincare(P[0], fc, 0.75, P[3]);

            hlerp2_poincare(P[8], fc, 0.25, P[7]);
            hlerp2_poincare(P[8], fc, 0.5, P[6]);
            hlerp2_poincare(P[8], fc, 0.75, P[5]);

            hlerp2_poincare(P[0], P[8], 0.25, P[11]);
            hlerp2_poincare(P[0], P[8], 0.5, P[10]);
            hlerp2_poincare(P[0], P[8], 0.75, P[9]);

            Vector<Complex <int> > segments(11);
            int nSegments = 0;

            segments[nSegments][0] = 0;
            segments[nSegments++][1] = 4;

            segments[nSegments][0] = 0;
            segments[nSegments++][1] = 10;

            segments[nSegments][0] = 11;
            segments[nSegments++][1] = 5;

            segments[nSegments][0] = 10;
            segments[nSegments++][1] = 6;

            segments[nSegments][0] = 9;
            segments[nSegments++][1] = 7;

            segments[nSegments][0] = 1;
            segments[nSegments++][1] = 11;

            segments[nSegments][0] = 2;
            segments[nSegments++][1] = 10;

            segments[nSegments][0] = 3;
            segments[nSegments++][1] = 9;

            segments[nSegments][0] = 3;
            segments[nSegments++][1] = 5;

            segments[nSegments][0] = 2;
            segments[nSegments++][1] = 6;

            segments[nSegments][0] = 1;
            segments[nSegments++][1] = 7;

            DrawTiling(mg, list, P, nSegments, segments, NULL, tolerance, settings);
            */

            /*
            Vector<complex> P(19);
            complex fc;

            P[8][0] = r; P[8][1] = 0.;
            calcFaceCenter(p, q, fc);
            P[4] = fc;
            A.apply(P[8], P[0]);

            hlerp2_poincare(P[0], fc, 0.25, P[1]);
            hlerp2_poincare(P[0], fc, 0.5, P[2]);
            hlerp2_poincare(P[0], fc, 0.75, P[3]);

            hlerp2_poincare(P[8], fc, 0.25, P[7]);
            hlerp2_poincare(P[8], fc, 0.5, P[6]);
            hlerp2_poincare(P[8], fc, 0.75, P[5]);

            hlerp2_poincare(P[0], P[8], 0.25, P[11]);
            hlerp2_poincare(P[0], P[8], 0.5, P[10]);
            hlerp2_poincare(P[0], P[8], 0.75, P[9]);

            hlerp2_poincare(P[2], P[6], 0.5, P[14]);
            hlerp2_poincare(P[2], P[10], 0.5, P[12]);
            hlerp2_poincare(P[10], P[6], 0.5, P[13]);


            //hlerp2_poincare(P[12], P[11], 2., P[15]);
            //hlerp2_poincare(P[1], P[11], 2., P[16]);
            //hlerp2_poincare(P[11], P[16], 2., P[18]);
            //hlerp2_poincare(P[0], P[15], 2., P[17]);

            Isometry J = Isometry::reflectionAcrossLine(P[0],P[8]);
            J.apply(P[1], P[15]);
            J.apply(P[2], P[17]);
            J.apply(P[12], P[16]);
            J.apply(P[14], P[18]);


            Vector<Complex <int> > segments(48);
            int nSegments = 0;

            segments[nSegments][0] = 3; segments[nSegments++][1] = 4;
            segments[nSegments][0] = 3; segments[nSegments++][1] = 5;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 3;
            segments[nSegments][0] = 3; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 5;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 6;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 2;
            segments[nSegments][0] = 2; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 14;
            segments[nSegments][0] = 14; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 6;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 7;
            segments[nSegments][0] = 0; segments[nSegments++][1] = 1;
            segments[nSegments][0] = 1; segments[nSegments++][1] = 11;
            segments[nSegments][0] = 11; segments[nSegments++][1] = 12;
            segments[nSegments][0] = 12; segments[nSegments++][1] = 10;
            segments[nSegments][0] = 10; segments[nSegments++][1] = 13;
            segments[nSegments][0] = 13; segments[nSegments++][1] = 9;
            segments[nSegments][0] = 9; segments[nSegments++][1] = 7;
            segments[nSegments][0] = 0; segments[nSegments++][1] = 11;
            segments[nSegments][0] = 11; segments[nSegments++][1] = 10;
            segments[nSegments][0] = 0; segments[nSegments++][1] = 15;
            segments[nSegments][0] = 15; segments[nSegments++][1] = 11;
            segments[nSegments][0] = 11; segments[nSegments++][1] = 16;
            segments[nSegments][0] = 16; segments[nSegments++][1] = 10;
            segments[nSegments][0] = 15; segments[nSegments++][1] = 16;
            segments[nSegments][0] = 15; segments[nSegments++][1] = 17;
            segments[nSegments][0] = 17; segments[nSegments++][1] = 16;
            segments[nSegments][0] = 16; segments[nSegments++][1] = 18;

            DrawTiling(mg, list, P, nSegments, segments, NULL, tolerance, settings);
            */

        }



        if(settings->drawSchwarzPolygon){
            mg->setColor(0xFFFFFF);
            for(int i=0; i<return_pickableSchwarz.SchwarzPolygon.length; i++)
                DrawLine2fPoincare(mg,
                                   return_pickableSchwarz.SchwarzPolygon[i],
                                   return_pickableSchwarz.SchwarzPolygon[(i+1)%return_pickableSchwarz.SchwarzPolygon.length],
                                   tolerance,
                                   settings);

            complex floodFillPoint(0.,0.);
            if(settings->doCurveLines){
                hlerp2_poincare(return_pickableSchwarz.SchwarzPolygon[0], return_pickableSchwarz.SchwarzPolygon[1], 0.5, floodFillPoint);
                hlerp2_poincare(floodFillPoint, return_pickableSchwarz.SchwarzPolygon[2], 0.5, floodFillPoint);
            }else{
                floodFillPoint[0] += return_pickableSchwarz.SchwarzPolygon[0][0]; floodFillPoint[1] += return_pickableSchwarz.SchwarzPolygon[0][1];
                floodFillPoint[0] += return_pickableSchwarz.SchwarzPolygon[1][0]; floodFillPoint[1] += return_pickableSchwarz.SchwarzPolygon[1][1];
                floodFillPoint[0] += return_pickableSchwarz.SchwarzPolygon[2][0]; floodFillPoint[1] += return_pickableSchwarz.SchwarzPolygon[2][1];
                floodFillPoint[0] /= 3.; floodFillPoint[1] /= 3.;
            }
            mg->floodFill(floodFillPoint[0], floodFillPoint[1]);

        }

        if(settings->drawBoundary){
            mg->setColor(settings->boundaryColor);
            mg->drawArc(-1.,-1.,2.,2.,0.,2*PI);
        }

    }__finally{
        mg->endDraw();
    }
} // DrawOmniTruncatedTiling
//------------------------------------------------------------------------






void DrawUniformTiling(MyGraphics* mg,
                        TessellationSettings* settings,
                        Isometry& F0,
                        Vector<int>& pp, int q, Vector<int>& backEdgeInds,
                        Vector<double>& wythoffCoeffs,
                        Isometry& return_smallestIsometry,
                        PickableSchwarzPolygon& return_pickableSchwarz)
{
    double tolerance = 1e-3;
    double halfEdgeLength = calcUniformTilingHalfEdgeLength(pp,q,1);
    double R = h2eNorm(2*halfEdgeLength);
    double r = h2eNorm(halfEdgeLength);

    Vector<double> angles(pp.length * q);
    {
        double angle = 0.;
        for(int i=0; i<pp.length*q; i++){
            angles[i] = angle;
            angle += 2*polygonHalfAngle(pp[i%pp.length], halfEdgeLength);
        }
    }

    Vector<complex> neighbors(pp.length * q);
    Vector<complex> edgeCenters(pp.length * q);
    {
        for(int i=0; i<pp.length*q; i++){
            double angle = angles[i];
            neighbors[i][0] = R*cos(angle);
            neighbors[i][1] = R*sin(angle);
            edgeCenters[i][0] = r*cos(angle);
            edgeCenters[i][1] = r*sin(angle);
        }
    }

    Vector<complex> faceCenters(pp.length*q);
    {
        for(int i=0; i<pp.length*q; i++){
            complex P2;
            calcFaceCenterFromHalfEdgeLength(pp[i%pp.length], halfEdgeLength, P2);
            faceCenters[i] = P2;
            // ... and rotate by angles[i]
            Isometry::pureRotation(angles[i]).apply(faceCenters[i], faceCenters[i]);
        }
    }

    Vector<Isometry> gens(pp.length * q);
    {
        for(int i=0; i<pp.length*q; i++){
            bool isRotation = false; // it's a reflection by default
            int j = backEdgeInds[i%pp.length];
            if (j < 0){
                j = ~j;
                isRotation = true;
            }
            j += i/pp.length*pp.length;

            // rotate edge j to x axis...
            Isometry G = Isometry::pureRotation(-angles[j]);
            if (isRotation){
                // rotate x axis to -x axis (180 degrees)
                G = Isometry(complex(-1.,0.), complex(0.,0.), 1) * G;
            }else{
                // reflect x axis to -x axis,
                // i.e. reflect across y axis.
                // This is equivalent to rotation by 180 degrees
                // followed by reflection across x axis.
                G = Isometry(complex(-1.,0.), complex(0.,0.), -1) * G;
            }
            // rotate x axis to edge i...
            G = Isometry::pureRotation(angles[i]) * G;
            // move origin directly to neighbors[i].
            G = Isometry::pureTranslation(neighbors[i][0], neighbors[i][1]) * G;

            gens[i] = G;
        }
    }

//#ifdef NOTYET // XXX
//    double tooFar = globals.dragging ? 1 - 1e-2 : !globals.do_quality ? 1 - 5e-3 : 1 - 1e-4;
//    double tooSmall = globals.dragging ? 1e-2 : !globals.do_quality ? 5e-3 : 1e-4;
//#endif // NOTYET
    double tooSmall = 1e-2; // XXX what does this mean?
    double tooFar = 1 - tooSmall;

    //double x0 = 0, y0 = 0, x1 = R, y1 = 0;
    Vector<Isometry> list(settings->maxIsometries);

    Vector<int> compositionLengths(settings->maxIsometries);

    //if (drawSnub && useNewShellMethod)
    //    compositionLengths = new int[maxIsometries];

    int nIsometries;

    nIsometries = UnCachedEnumerateIsometry2GroupForUniformTiling(
                            F0,
                            gens.length,
                            gens,
                            pp,
                            q,
                            backEdgeInds,
                            settings->maxIsometries,
                            list,
                            compositionLengths,
                            return_smallestIsometry,
                            settings->maxLevels,
                            tooFar,
                            tooSmall);




    mg->beginDraw();
    try{

        /*
        //if (drawGreenGridForDebugging)
        {
            mg->setColor(0x8000);
            int subDiv = 20;
            for(int i=0; i<subDiv+1; i++){
                mg->drawLine(LERP(-1., 1., double(i)/double(subDiv)), -1.,
                            LERP(-1., 1., double(i)/double(subDiv)), 1.);
                mg->drawLine(-1., LERP(-1., 1., double(i)/double(subDiv)),
                            1., LERP(-1., 1., double(i)/double(subDiv)));
            }
        }
        */
        

        if(settings->drawDual){
            mg->setColor(settings->dualColor);

            // Could just connect each face center to the next,
            // But it's better to bisect, because that way
            // things meet up correctly if we happen to be drawing
            // line approximations.
            // meet up d
            Vector<complex> P(2*edgeCenters.length);
            Vector< Complex<int> > segments(2*edgeCenters.length);
            for(int i=0; i<edgeCenters.length; i++){
                P[2*i+0] = edgeCenters[i];
                P[2*i+1] = faceCenters[i];
            }
            for(int i=0; i<2*edgeCenters.length; i++){
                segments[i][0] = i;
                segments[i][1] = (i+1)%segments.length;
            }

            DrawTiling(mg, list, P, segments.length, segments, NULL, false,  tolerance, settings);
        }

        return_pickableSchwarz.setVertices(faceCenters, F0);

        /*
        boolean drawSchwarzPolygon = false; // set to true to debug
        if (drawSchwarzPolygon)
        {
            mg.setColor(java.awt.Color.red);
            Complex sp[] = return_pickableSchwarz.SchwarzPolygon;
            int i;
            FOR (i, sp.length)
                DrawLine2fPoincare(mg,
                                   sp[i].x,
                                   sp[i].y,
                                   sp[(i+1)%sp.length].x,
                                   sp[(i+1)%sp.length].y,
                                   tolerance,
                                   hyperbolicModel,
                                   doCurveLines,
                                   doAntiAlias,
                                   doRandomJitter);
        }
        */


        if (settings->drawSnub){
            mg->setColor(settings->snubColor);

            Vector<complex> P(neighbors.length);
            Vector< Complex<int> > segments(neighbors.length);

            //if (wythoffCoeffs == NULL) // this will be the case the first time around, meaning put the vertex in the center (it was hard to figure out what the wythoff coeffs for that were).  It will also be the case if not allowed to drag the vertex at all.
            //{
            //    P[0] = zero;
            //    for(int i=0; i<neighbors.length; i++){
            //        P[i] = neighbors[i];
            //    }
            //}
            //else
            {
                complex P2;
                if (!hBary2(P2, wythoffCoeffs, faceCenters))
                    return; // XXX print warning?

                //mg.drawPoint(P2[0], P2[1], 3);

                for(int i=0; i<neighbors.length; i++){
                    calcClosestPointOnLine2(P[i], P2, faceCenters[(i-1+neighbors.length)%neighbors.length], faceCenters[i]);
                    hlerp2_poincare(P2, P[i], 2., P[i]);
                }
            }

            for(int i=0; i<neighbors.length; i++){
                segments[i][0] = i;
                segments[i][1] = (i+1)%neighbors.length;
            }

            // For snub, we only want half of the isometries, namely
            // the ones an even distance from the seed.
            // Usually this can be deduced by the sign (i.e. reflection or
            // not) of the isometry, but we want to make it work
            // even when there are some non-reflection isometries
            // (e.g. the weird (6,4,4,4) tiling), so we actually measure
            // distance, i.e. composition length.
            Vector<Isometry> snubbedList(nIsometries);
            int nSnubbedIsometries = 0;
            {
                int desiredParity = 1 - settings->snubParity; // we draw *around* the opposite-parity vertices, so that we retain snubParity vertices.
                for(int i=0; i<nIsometries; i++){

                    //if (compositionLengths != NULL && compositionLengths[i] % 2 == desiredParity)
                    if (compositionLengths[i] % 2 == desiredParity)
                        snubbedList[nSnubbedIsometries++] = list[i];
                }
            }

            DrawTiling(mg, snubbedList, P, segments.length, segments, NULL, false, tolerance, settings);
        }

        if (settings->drawPrimal){
            mg->setColor(settings->primalColor);

            Vector<complex> P(neighbors.length + 1);
            Vector< Complex<int> > segments(neighbors.length);

            //if (wythoffCoeffs == NULL) // this will be the case the first time around, meaning put the vertex in the center (it was hard to figure out what the wythoff coeffs for that were).  It will also be the case if not allowed to drag the vertex at all.
            //{
            //    P[0] = zero;
            //    for(int i=0; i<neighbors.length; i++){
            //        //P[i+1] = neighbors[i];
            //        P[i+1] = edgeCenters[i];
            //    }
            //}
            //else
            {
                complex P2;
                if (!hBary2(P2, wythoffCoeffs, faceCenters))
                    return; // XXX print warning?

                //mg.drawPoint(P2[0], P2[1], 3);

                P[0] = P2;
                for(int i=0; i<neighbors.length; i++){
                    calcClosestPointOnLine2(P[i+1], P[0], faceCenters[(i-1+neighbors.length)%neighbors.length], faceCenters[i]);
                }
            }

            for(int i=0; i<neighbors.length; i++){
                segments[i][0] = 0;
                segments[i][1] = i+1;
            }

            DrawTiling(mg, list, P, segments.length, segments, NULL, false, tolerance, settings);
        }

        if (settings->drawBoundary){
            //
            // Draw circle...
            //
            mg->setColor(settings->boundaryColor);
            mg->drawArc(-1.,-1.,2.,2.,0.,2*PI);
        }

    }__finally{
        mg->endDraw();
    }
} // DrawUniformTiling