/* Hyperbolic Tessellations Copyright (C) 2007 Dmitry Brant http://dmitrybrant.com Based largely on code from http://www.plunk.org/~hatch/HyperbolicApplet/ */ //--------------------------------------------------------------------------- #ifndef hyperbolicutilsH #define hyperbolicutilsH #define PI 3.141592653589793238 #define POINCARE_DISK 0 #define KLEIN_DISK 1 #define UPPER_HALF_PLANE 2 #define ANTIKLEIN_DISK 3 #include "mycomplex.h" #include "isometry.h" #include "mathfuncs.h" class LerpFunction; class Hlerp2_poincare; class Hlerp2_uhp; class Hlerp2_polar; class Hlerp2_simple_inversion; class Hlerp2_euclidean; double HYPOTSQRD(double a, double b); // hyperbolic to euclidean norm (distance from 0,0) in Poincare disk. double h2eNorm(double hNorm); // euclidean to hyperbolic norm (distance from 0,0) in Poincare disk double e2hNorm(double eNorm); //-------------------------------------------------- // Given a point p and nVerts verts, // calculate the barycentric coordinates void unHBary2(Vector& result, complex& p, Vector& verts); bool hBary2(complex& result, Vector& b, // barycentric weights Vector& p); // This isn't really barycentric anything, // but it just happens sometimes we know the distances // to all the verts. bool hBary2FromDistances(complex& result, Vector& dists, // barycentric weights Vector& p); // Calculate the intersection of p0p1 with q0q1. bool calcLineLineIntersection2(complex& result, complex& p0, complex& p1, complex& q0, complex& q1); // Calculate the point on the given hyperbolic line that is closest to the origin. // The function return value is the hyperbolic distance from the origin to that point. double calcSmallestPointOnLine2(complex& result, complex& p0, complex& p1); double calcClosestPointOnLine2(complex& result, complex& p, complex& l0, complex& l1); // Calculate the point on the boundary of the disk // obtained by following the geodesic from p0 towards p1. // p0 and p1 must be distinct. void _calcEndOfGeodesic2(complex& result, complex& p0, complex& p1); // Calculate a line originating from the vertex p0 // in the direction of the other vertices p[0..nPoints-1], // whose perpendicular distance to the face opposite p[i] // has sinh proportional to b[i]. // This is equivalent to saying that close to p0, // the distances to the respective faces are proportional to the b[i]'s. void _calcPointOnWeightedSector2(complex& result, complex& p0, Vector& p, Vector& b); double calcContentWithOrigin(complex& v); // calculate unsigned k-dimensional content // spanned by k n-dimensional vectors // together with the origin. // If n == k then this is the absolute value of the determinant. double calcContentWithOrigin(int n, int k, Vector& _v); // calculate content of k-1-dimensional simplex with given k n-dimensional // vertices. (actually (k-1)! times it, but we only care about proportions) double calcSimplexContent(int n, int k, Vector& p); bool pointIsToRightSideOfLinePoincare(complex p, complex l0, complex l1, double eps); void calcFaceCenter(double p, int q, complex& result); void calcFaceCenterFromHalfEdgeLength(double p, double halfEdgeLength, complex& result); double polygonHalfAngle(double p, double halfEdgeLength); double polygonCircumRadius(double p, double halfEdgeLength); double polygonInRadius(double p, double halfEdgeLength); double calcUniformTilingHalfEdgeLength(Vector& p, int q); double calcUniformTilingHalfEdgeLength(Vector& p, double q, double density); double calcUniformTilingHalfEdgeLength(Vector& p, double q, double density); // Calculate the edge length for the hyperbolic // tiling with schlafli symbol {p,q}, // i.e. q p-gons surrounding each vertex. double calcRegularTilingHalfEdgeLength(int p, int q); // // Label a hyperbolic triangle's sides A,B,C // and the opposite angles alpha,beta,gamma. // Given alpha,C,beta, find A. // double solveTriangleAFromAlphaCBeta(double alpha, double C, double beta); // // Label a hyperbolic triangle's sides A,B,C // and the opposite angles alpha,beta,gamma. // Given A,gamma,B, find C. // double solveTriangleCFromAGammaB(double A, double gamma, double B); // // Label a hyperbolic triangle's sides A,B,C // and the opposite angles alpha,beta,gamma. // Given A,B,C, find alpha. // From http://home.hetnet.nl/~vanadovv/Gonio.html double solveTriangleAlphaFromABC(double A, double B, double C); double hangle(complex& p0, complex& p1, complex& p2); // Very special-case and ad-hoc... bool findBaryCoordsThatMakeSnubUniform(Vector& result, Vector& verts); int UnCachedEnumerateIsometry2GroupForUniformTiling( Isometry& F0, int nGens, // equal to valence of a vertex, including rotationalMultiplicity Vector& gens, Vector& pp, // face list, not expanded by rotationalMultiplicity int rotationalMultiplicity, Vector& backEdgeInds, // not expanded by rotationalMultiplicity int maxIsometries, Vector& return_isometries, Vector return_isometry_composition_lengths, Isometry& return_smallestIsometry, int maxLevels, double tooFar, double tooSmall); int UnCachedEnumerateIsometryGroup(Isometry F0, int nGens, Isometry gens[], int rotationalMultiplicity, int maxIsometries, Isometry return_isometries[], Isometry& return_smallestIsometry, int maxLevels, double tooFar, double tooSmall); int subdivideArc(complex& p0, complex& p1, double tolerance, int maxIntermediatePoints, int startI, Vector& result_intermediatePoints); double hlerp2_poincare(complex& p0, complex& p1, double t, complex& result); // the version with aspect double hlerp2_poincare(complex& p0, complex& p1, double t, double aspect, complex& result); double hdist2_poincare(complex p0, complex p1); double hlerp2_uhp(complex& p0, complex& p1, double t, double aspect, complex& result); double hlerp2_simple_inversion(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result); // distance to circle centered at (0,-1). // very accurate for points near the origin. double distToCircleCenteredAtMinus1(double x, double y); // so we can get insanely accurate near the origin! use the above function. double hlerp2_polar(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result); bool hspline2_polar(double t0, complex& p0, double t1, complex& p1, double t2, complex& p2, double t3, complex& p3, double t, double aspect, bool conformal, LerpFunction& lerp, complex& result); // circular-inversion isometry between Poincare disk and upper half plane void canonicalCircularInvert(double x, double y, complex& result); bool hlerp3_polar(complex& p0, complex& p1, double t, double aspect, complex& result); // compose with a reflection of the disk about the x axis, for a more intuitive conversion... void d2uhp(double x, double y, complex& result); void uhp2d(double x, double y, complex& result); class LerpFunction{ public: // Subclass must override this. // returns distance, or -1. on failure/confusion virtual double apply( complex& p0, complex& p1, double t, double aspect, // > 1 means higher arcs, might be meaningless for some kinds of lerp bool conformal, // only meaningful for some kinds of lerp; means do some fudging to make it conformal if possible complex& result); // if null, double dist(complex& p0, complex& p1, double aspect, bool conformal){ return apply(p0, p1, 0.5, aspect, conformal, complex()); } }; //--------------------------------------------------------------------------- class Hlerp2_poincare : public LerpFunction{ public: double apply(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result) { return hlerp2_poincare(p0,p1,t,aspect,result); } }; //--------------------------------------------------------------------------- class Hlerp2_uhp : public LerpFunction{ public: double apply(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result) { return hlerp2_uhp(p0,p1,t,aspect,result); } }; //--------------------------------------------------------------------------- class Hlerp2_polar : public LerpFunction{ public: double apply(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result) { return hlerp2_polar(p0,p1,t,aspect,conformal,result); } }; //--------------------------------------------------------------------------- class Hlerp2_simple_inversion : public LerpFunction{ public: double apply(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result) { return hlerp2_simple_inversion(p0,p1,t,aspect,conformal,result); } }; //--------------------------------------------------------------------------- class Hlerp2_euclidean : public LerpFunction{ public: double apply(complex& p0, complex& p1, double t, double aspect, bool conformal, complex& result) { result = p0*(1.0-t) + p1*t; return sqrt(distsqrd(p0,p1)); } }; //--------------------------------------------------------------------------- class Node{ public: Node(int valence, int sign, int start){ this->valence = valence; neighbors = new int[valence]; backEdgeInds = new int[valence]; for(int i=0; isign = sign; this->start = start; } ~Node(){ delete [] neighbors; delete [] backEdgeInds; } int* neighbors; // index into an array of Nodes int* backEdgeInds; // usually the global ones, but can differ if we reach an edge from two different directions and they don't match int sign; // with respect to original seed node int start; // where to start enumerating neighbors for nice ordering int valence; Node(const Node& from){ sign = from.sign; start = from.start; valence = from.valence; neighbors = new int[valence]; backEdgeInds = new int[valence]; for(int i=0; i