SymmetryUtils.prejava

#include "macros.h"

import com.donhatchsw.util.Arrays;
import com.donhatchsw.util.VecMath;

final class SymmetryUtils {

    private SymmetryUtils(){ throw new AssertionError(); } // non-instantiatable util class

    // Note, it turns out this is called with only leftRightSymmetry=false, currently...
    // i.e. it always returns a quad, not a triangle.
    public static double[][] getFundamentalRegionVerts(int p, int q, boolean leftRightSymmetry)
    {
        CHECK(false);
        int verboseLevel = 0;
        if (verboseLevel >= 1) System.out.println("        in getFundamentalRegionVerts(p="+p+", q="+q+", leftRightSymmetry="+leftRightSymmetry+")");
        CHECK_GE(p, 2);
        CHECK_GE(q, 2);
        double answer[][];
        if ((p-2)*(q-2) < 4)
        {
            // Verts of schwarz triangle:
            //     V2 = +z axis, angle pi/p
            //     V0 = between +z and +y axis, angle pi/q
            //     V1 = somewhere in +x +y +z, angle pi/2
            // By https://en.wikipedia.org/wiki/Solution_of_triangles#Three_angles_given_.28spherical_AAA.29,
            // the corresponding spherical side lengths are:
            //     S01 = acos((cos(pi/p) + cos(pi/q)cos(pi/2)) / (sin(pi/q)sin(pi/2)))
            //         = acos((cos(pi/p) / sin(pi/q))
            //     S12 = acos((cos(pi/q) + cos(pi/2)cos(pi/p)) / (sin(pi/2)sin(pi/p)))
            //         = acos((cos(pi/q) / sin(pi/p))
            //         = (pi-dihedral{p,q})/2
            //     S20 = acos((cos(pi/2) + cos(pi/p)cos(pi/q)) / (sin(pi/p)sin(pi/q))
            //         = acos(cos(pi/p)cos(pi/q) / (sin(pi/p)sin(pi/q))
            //         = acos(cot(pi/p)cot(pi/q))
            double cosSqrdS12 = ExactTrig.cosSquaredPiTimes(1,q) / ExactTrig.sinSquaredPiTimes(1,p);
            double cosS12 = Math.sqrt(cosSqrdS12);
            double sinS12 = Math.sqrt(1.-cosSqrdS12);
            double cosSqrdS20 = ExactTrig.cotSquaredPiTimes(1,p) * ExactTrig.cotSquaredPiTimes(1,q);  // increases with p,q
            double cosS20 = Math.sqrt(cosSqrdS20);
            double sinS20 = Math.sqrt(1.-cosSqrdS20);

            double V2[] = {0,0,1};
            // V1 is: +z axis, rotated +z->+y by s12, rotated +y->+x by pi/p
            // = {0,cosS12,sinS12} rotated +y->+x by pi/p
            double V1[] = VecMath.vxm(new double[]{0,sinS12,cosS12},
                                      VecMath.makeRowRotMat(3, 1, 0, Math.PI/p));
            // V0 is: +z axis, rotated +z->+y by S20
            double V0[] = {0,sinS20,cosS20};


            if (leftRightSymmetry)
            {
                // fundamental region is a triangle
                answer = new double[][] {V2, V1, V0};
            }
            else
            {
                // fundamental region is a quad
                answer = new double[][] {V2, V1, V0, new double[]{-V1[0],V1[1],V1[2]}};
            }
        }
        else
        {
          // WORK IN PROGRESS
          // Hmm, I'm interpreting [2] as h here.  That's probably different from the above, and probably not going to work.
          double cosSquaredPiOverQ = ExactTrig.cosSquaredPiTimes(1,q);
          double cosSquaredPiOverP = ExactTrig.cosSquaredPiTimes(1,p);
          double[] V2 = {0,0,0};
          double[] V1 = {cosSquaredPiOverP*.5, Math.sqrt(cosSquaredPiOverQ*cosSquaredPiOverP)*.5,Double.NaN};
          double[] V0 = {0,.5,Double.NaN};
          V0[2] = -(SQR(V0[0])+SQR(V0[1]))/2.;
          V1[2] = -(SQR(V1[0])+SQR(V1[1]))/2.;
          answer = new double[][] {V2, V1, V0, new double[]{-V1[0],V1[1],V1[2]}};
        }

        if (verboseLevel >= 1) System.out.println("        out getFundamentalRegionVerts(p="+p+", q="+q+", leftRightSymmetry="+leftRightSymmetry+"), returning "+Arrays.toStringCompact(answer));
        return answer;
    } // getFundamentalRegionVerts

    // 3d symmetry groups.
    // Rotation of order p around z axis,
    // and rotation of order q around some point between the +z and +y axis.
    // Note, we don't return a minimal set;
    // rather, we add in mathematically redundant generators
    // for all powers of p and q, since the math for generating them
    // is a bit more accurate than we'd get by products of the minimal set of generators.
    //
    // Instead of returning a flat list, we return a list of sub-lists;
    // each of the sub-lists are assumed to be powers of a single generator,
    // so there's no need to try more than one thing out of a given sub-list consecutively.
    //
    private static double[][][/*4*/][/*4*/] computeSymmetryGroupGenerators3d(int p, int q, boolean leftRightSymmetry, boolean sphereCentralSymmetry)
    {
        double genGroups[][][][] = new double[0][][][];

        // add reflection first, so it will get favored when generating,
        // since multiplying by it doesn't add any roundoff error.
        if (leftRightSymmetry)
            genGroups = (double[][][][])Arrays.append(genGroups, new double[][][]{{{-1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}}});

        if (sphereCentralSymmetry)
            genGroups = (double[][][][])Arrays.append(genGroups, new double[][][]{{{-1,0,0,0},{0,-1,0,0},{0,0,-1,0},{0,0,0,1}}});



        if (p > 1)
        {
            genGroups = (double[][][][])Arrays.append(genGroups, new double[p-1][][]);
            double subGens[][][] = genGroups[genGroups.length-1];

            for (int i = 1; i < p; ++i) // skip 0
            {
                double c = ExactTrig.cosPiTimes(2*i,p); // cos(2*pi * i/p)
                double s = ExactTrig.sinPiTimes(2*i,p); // sin(2*pi * i/p)
                double gen[][] = {
                    { c,s,0,0},
                    {-s,c,0,0},
                    { 0,0,1,0},
                    { 0,0,0,1},
                };
                subGens[i-1] = gen;
            }
        }
        if (q > 1)
        {
            genGroups = (double[][][][])Arrays.append(genGroups, new double[q-1][][]);
            double subGens[][][] = genGroups[genGroups.length-1];

            // The rotation of order p is centered on the z axis;
            // where do we want to center the rotation of order q?
            // Figure out the cosine and sine of the tilt,
            // that is, distance from vert to face center of spherical {p,q}.
            double cosSqrdTilt = MIN(p,q)==1 ? 1. : // angle = 0; just produce identity matrix so no effect   XXX think about this again
                                 ExactTrig.cotSquaredPiTimes(1,p)*ExactTrig.cotSquaredPiTimes(1,q);
            double cosTilt = Math.sqrt(cosSqrdTilt);
            double sinTilt = Math.sqrt(1.-cosSqrdTilt);
            // tilt from +z towards +y
            double tiltMat[][] = {
                {1,      0,       0,0},
                {0,cosTilt,-sinTilt,0},
                {0,sinTilt, cosTilt,0},
                {0,      0,       0,1},
            };
            for (int i = 1; i < q; ++i) // skip 0
            {
                double c = ExactTrig.cosPiTimes(2*i,q); // cos(2*pi * i/q)
                double s = ExactTrig.sinPiTimes(2*i,q); // sin(2*pi * i/q)
                double gen[][] = {
                    { c,s,0,0},
                    {-s,c,0,0},
                    { 0,0,1,0},
                    { 0,0,0,1},
                };
                gen = VecMath.mxmxm(VecMath.transpose(tiltMat), gen, tiltMat);
                subGens[i-1] = gen;
            }
        }

        return genGroups;
    } // computeSymmetryGroupGenerators3d

    // 3d symmetry groups that hold origin fixed.
    // rotation of order p around z axis,
    // and rotation of order q around some point between the +z and +y axis.
    public static double[][/*4*/][/*4*/] computeSymmetryGroup3d(int p, int q, boolean leftRightSymmetry, boolean sphereCentralSymmetry, boolean q346meanRepeatRegardlessOfP)
    {
        int verboseLevel = 0;
        if (verboseLevel >= 1) System.out.println("        in computeSymmetryGroup3d(p="+p+", q="+q+", leftRightSymmetry="+leftRightSymmetry+", sphereCentralSymmetry="+sphereCentralSymmetry+", q346meanRepeatRegardlessOfP="+q346meanRepeatRegardlessOfP+")");

        int repeatQ = 1;
        if ((q346meanRepeatRegardlessOfP && (q==3||q==4||q==6)) || (p-2)*(q-2)==4) {
           // in this case, we do *not* handle q in the initial part, we handle it in the repeat logic at the end of the function.
           repeatQ = q;
           q = 1;
        }

        double gens[][][][] = computeSymmetryGroupGenerators3d(p, q, leftRightSymmetry, sphereCentralSymmetry);

        if (false)
            System.out.println("gens = "+Arrays.toStringNonCompact(gens,
                                                                   "", // indentString
                                                                   "    "));  // indentIncr

        int nExpected = (MIN(p,q) == 1 ? MAX(p,q) :
                         MIN(p,q) == 2 ? 2*MAX(p,q) :
                         MIN(p,q)==3 && MAX(p,q)==3 ? 12 :
                         MIN(p,q)==3 && MAX(p,q)==4 ? 24 :
                         MIN(p,q)==3 && MAX(p,q)==5 ? 60 : -1);
        CHECK_NE(nExpected, -1);
        if (leftRightSymmetry)
            nExpected *= 2;
        if (sphereCentralSymmetry)
            nExpected *= 2;
        // In some cases, leftRightSymmetry is redundant with sphereCentralSymmetry.
        // I don't know of any rhyme or reason to this, it just is.
        if (leftRightSymmetry && sphereCentralSymmetry)
            if ((MIN(p,q)==2 && MAX(p,q)%2==0) || MIN(p,q)==3&&MAX(p,q)>=4)
            {
                nExpected /= 2;
            }

        double group[][][] = new double[nExpected][4][4];
        int lastSubgroupIndex[] = new int[group.length];

        int n = 0;
        VecMath.identitymat(group[n]);
        lastSubgroupIndex[n] = -1; // so nothing will match it
        n++;

        double scratch[][] = new double[4][4];
        FORI (i, n) // while n is growing
        {
            FORI (iSubgroup, gens.length)
            {
                if (iSubgroup == lastSubgroupIndex[i])
                    continue; // no need to look at two in a row from the same subgroup
                double subgroupGens[][][] = gens[iSubgroup];
                FORI (iGen, subgroupGens.length)
                {
                    VecMath.mxm(scratch, subgroupGens[iGen], group[i]); // or other order? not sure it matters
                    int j;
                    FOR (j, n)
                        if (VecMath.equals(scratch, group[j], 1e-3))
                            break;
                    if (j == n) // if didn't find it
                    {
                        VecMath.copymat(group[n], scratch);
                        lastSubgroupIndex[n] = iSubgroup;
                        n++;
                    }
                }
            }
        }
        CHECK_EQ(n, group.length);

        if (repeatQ != 1)
        {
            //  THINK ABOUT THIS:
            //  What's a homogeneous matrix that pans in x and keeps points the same offset from the parabola?
            //  So, homogeneous:
            //     -2,4,1 -> -1,1,1
            //     -1,1,1 ->  0,0,1
            //      0,0,1 ->  1,1,1
            //      1,1,1 ->  2,4,1
            //      2,4,1 ->  3,9,1
            //      3,9,1 ->  4,16,1
            //
            //  So the matrix must be:
            //      ? ? ?
            //      ? ? ?
            //      1 1 1
            //
            //  Oh let's see, here's a clue: x^2+x is an offset paraboloid, right?
            //  So, if we have x^2 and x, can we express x+1, (x+1)^2?  
            //  Well, (x+1)^2 = x^2 + 2x + 1
            //  So we want x,x^2,1 -> x+1,x^2+2x+1,1
            //             a,b,1 -> a+1,2a+b+1,1
            //             1,0,1 -> 2,3,1
            //             0,1,1 -> 1,2,1
            //             0,0,1 -> 1,1,1
            //
            //             1,0,0 -> 1,2,0
            //             0,1,0 -> 0,1,0
            //             0,0,1 -> 1,1,1

            //  So, is that the answer?
            //  Yes, it is!!
            //  And, by 2 instead of 1:
            //    {{1,2,0},{0,1,0},{1,1,1}}^2
            //      1 4 0
            //      0 1 0
            //      2 4 1
            //    {{1,2,0},{0,1,0},{1,1,1}}^3
            //      1 6 0
            //      0 1 0
            //      3 9 1
            //  So, to move by dx:
            //      1 2*dx 0
            //      0  1   0
            //     dx dx^2 1
            //  In one higher dimension:
            //      1   0     2*dx   0
            //      0   1     2*dy   0
            //      0   0      1     0
            //     dx  dy  dx^2+dy^2 1
            // Ok, but that was the paraboloid z=x^2+y^2, which has curvature 2 at the origin.
            // We want the one that has curvature 1 at the origin.
            // So that's:  steepen * m * flatten
            // I think that's:
            //      1   0     dx         0
            //      0   1     dy         0
            //      0   0      1         0
            //     dx  dy  (dx^2+dy^2)/2 1
            // but let's make sure.
            //    {{1,0,0},{0,2,0},{0,0,1}} . {{1,2,0},{0,1,0},{1,1,1}} . {{1,0,0},{0,1/2,0},{0,0,1}}
            //        1  1  0
            //        0  1  0
            //        1 1/2 1
            // and 2 units is:
            //        1 2 0
            //        0 1 0
            //        2 2 1
            // so for general dx it's:
            //        1  dx    0
            //        0  1     0
            //       dx dx^2/2 1
            // oh and actually the curvature has to be negative:
            //        1   -dx   0
            //        0    1    0
            //       dx -dx^2/2 1


            if (verboseLevel >= 1) System.out.println("          doing repeat stuff");
            if (verboseLevel >= 1) System.out.println("              group was "+Arrays.toStringNonCompact(group, "", "    "));
            // WORK IN PROGRESS
            // Experiment with planar repeat.

            double[][][] translations;

            CHECK(repeatQ == 3 || repeatQ == 4 || repeatQ == 6);
            if (repeatQ == 3 || repeatQ == 6)
            {
                int r = repeatQ==6 ? 3 : 2;
                // Put the "q" vertex at 0,.5
                double tileWidth = repeatQ==6 ? 1/2. : Math.sqrt(3.)/2.;
                translations = new double[HEXED(r+1)][4][4];
                int nTranslations = 0;
                for (int iy = -r; iy <= r; ++iy)
                for (int ix = (iy<0 ? -r-iy : -r); ix <= (iy<0 ? r : r-iy); ++ix) {
                    double dx = ix*tileWidth + iy*tileWidth/2.;
                    double dy = iy*tileWidth * Math.sqrt(3)/2.;
                    if (repeatQ == 6) {
                      double temp;
                      SWAP(dx, dy, temp);
                      dy *= -1;
                      if (MOD(iy-ix, 3) == 1) continue;
                    }
                    translations[nTranslations++] = new double[][] {
                        {1,0,-dx,0},
                        {0,1,-dy,0},
                        {0,0,1,0},
                        {dx,dy,-(dx*dx+dy*dy)/2.,1},
                    };

                    // OH!  It's wrong when repeatQ is 6; in this case,
                    // it doesn't suffice to do pure rotations; in this case
                    // we must also turn half of them upside down.
                    if (repeatQ == 6)
                    {
                        if (MOD(iy-ix, 3) == 2)
                        {
                            double[][] flip = {
                              {-1,0,0,0},
                              {0,-1,0,0},
                              {0,0,1,0},
                              {0,0,0,1},
                            };
                            translations[nTranslations-1] = VecMath.mxm(flip, translations[nTranslations-1]);
                        }
                    }
                }
                if (repeatQ == 6)
                {
                    // I don't know nor care the exact formula for how many points were omitted.
                    translations = (double[][][])Arrays.subarray(translations, 0, nTranslations);
                }
                CHECK_EQ(nTranslations, translations.length);
            }
            else
            {
                int r = 2;  // radius, in units where a unit is from one point to the next
                double tileWidth = 1.;
                int rx = r;
                int ry = r;
                translations = new double[(2*rx+1)*(2*ry+1)][4][4];
                int nTranslations = 0;
                for (int ix = -rx; ix <= rx; ++ix)
                for (int iy = -ry; iy <= ry; ++iy)
                {
                    //double dx = ix*tileWidth;
                    //double dy = iy*tileWidth;
                    double dx = (ix*tileWidth + iy*tileWidth)*.5;
                    double dy = (iy*tileWidth - ix*tileWidth)*.5;

                    translations[nTranslations++] = new double[][] {
                        {1,0,-dx,0},
                        {0,1,-dy,0},
                        {0,0,1,0},
                        {dx,dy,-(dx*dx+dy*dy)/2.,1},
                    };
                }
                CHECK_EQ(nTranslations, translations.length);
            }
            // Put the middle one (identity transform) first; that's what the caller will expect.
            Arrays.swap(translations,0, translations,(translations.length-1)/2);

            double bigGroup[][][] = new double[translations.length * group.length][4][4];
            FORI (iTranslation, translations.length)
            {
               FORI (iGroup, group.length)
               {
                  VecMath.mxm(bigGroup[iTranslation*group.length+iGroup], group[iGroup], translations[iTranslation]);
               }
            }
            group = bigGroup;

            if (verboseLevel >= 1) System.out.println("          done repeat stuff");
        }
        if (verboseLevel >= 1) System.out.println("        out computeSymmetryGroup3d(p="+p+", q="+q+", leftRightSymmetry="+leftRightSymmetry+", sphereCentralSymmetry="+sphereCentralSymmetry+", q346meanRepeatRegardlessOfP="+q346meanRepeatRegardlessOfP+"), returning "+Arrays.toStringNonCompact(group, /*indentString=*/"", /*indentIncr=*/"    "));
        return group;
    } // computeSymmetryGroup3d

    public static double[][/*4*/][/*4*/] getSymmetryGroup(int p, int q, boolean leftRight, boolean sphereCentral,  boolean wrapAroundSphereFlag, double wrapSphereCurvature, boolean centerSphereFlag, boolean q346meanRepeatRegardlessOfP)
    {
        int verboseLevel = 0;
        if (verboseLevel >= 1) OUT("    in getSymmetryGroup(p="+p+", q="+q+", leftRight="+leftRight+", sphereCentral="+sphereCentral+", wrapAroundSphereFlag="+wrapAroundSphereFlag+", wrapSphereCurvature="+wrapSphereCurvature+", centerSphereFlag="+centerSphereFlag+", q346meanRepeatRegardlessOfP="+q346meanRepeatRegardlessOfP+")");

        double group[][][] = computeSymmetryGroup3d(p, q, leftRight, sphereCentral, q346meanRepeatRegardlessOfP);

        if (!((q346meanRepeatRegardlessOfP && (q==3||q==4||q==6)) || (p-2)*(q-2)==4)) {  // if repeat, then we did the logic in the plane/paraboloid, so blow off any sphere stuff.
          if (!centerSphereFlag) // XXX TODO: wait a minute, this test can't be right when not wrapped, can it?
          {
              // un-center the symmetry group
              double wrapSphereRadius = 1./wrapSphereCurvature;
              double originToSphereCenter[][] = {
                  {1,0,0,0},
                  {0,1,0,0},
                  {0,0,1,0},
                  {0,0,-wrapSphereRadius,1},
              };
              double sphereCenterToOrigin[][] = {
                  {1,0,0,0},
                  {0,1,0,0},
                  {0,0,1,0},
                  {0,0,wrapSphereRadius,1},
              };
              FORI (iGroup, group.length)
              {
                  group[iGroup] = VecMath.mxmxm(sphereCenterToOrigin,
                                                group[iGroup],
                                                originToSphereCenter);
              }
          }
          if (!wrapAroundSphereFlag)
          {
              // unwrap the symmetry group
              // TODO: this assumes the symmetry group was wrapped to begin with... not true if it's all translational, is it??
              double wrapMat[][] = GeomUtils.getWrapAroundSphereMatrix(wrapSphereCurvature, centerSphereFlag);
              double unwrapMat[][] = GeomUtils.getUnwrapAroundSphereMatrix(wrapSphereCurvature, centerSphereFlag);
              if (verboseLevel >= 1) OUT("      wrapMat = "+Arrays.toStringCompact(wrapMat));
              if (verboseLevel >= 1) OUT("      unwrapMat = "+Arrays.toStringCompact(wrapMat));
              FORI (iGroup, group.length)
                  group[iGroup] = VecMath.mxmxm(wrapMat,
                                                group[iGroup],
                                                unwrapMat);
          }
        }
        if (verboseLevel >= 1) OUT("    out getSymmetryGroup(p="+p+", q="+q+", leftRight="+leftRight+", sphereCentral="+sphereCentral+", wrapAroundSphereFlag="+wrapAroundSphereFlag+", wrapSphereCurvature="+wrapSphereCurvature+", centerSphereFlag="+centerSphereFlag+", q346meanRepeatRegardlessOfP="+q346meanRepeatRegardlessOfP+"), returning "+Arrays.toStringNonCompact(group, /*indentString=*/"", /*indentIncr=*/"    "));
        return group;
    } // getTheSymmetryGroup

    // does at most each vert generator times each xform, omitting dups.
    //  CBB: should use a fuzzy point hash table if a lot of verts and/or xforms.
    public static double[][] generateVertices(double vertGenerators[][],
                                              double xforms[][][])
    {
        int verboseLevel = 1;  // 1 is a good number
        if (verboseLevel >= 1) OUT("    in generateVertices");
        double tol = 1e-9;
        double verts[][] = new double[vertGenerators.length*xforms.length][4];
        int nVerts = 0;
        double vert4in[] = {0,0,0,1};
        double vert4out[] = new double[4];
        double vert3out[] = new double[3];

        FORI (ixform, xforms.length)
        FORI (iVertGenerator, vertGenerators.length)
        {
            CHECK_EQ(vertGenerators[iVertGenerator].length, 3);
            // TODO: really need homogeneous stuff
            VecMath.copyvec(3, vert4in, vertGenerators[iVertGenerator]);
            VecMath.vxm(vert4out,
                        vert4in,
                        xforms[ixform]);
            VecMath.vxs(3, vert3out, vert4out, 1./vert4out[3]);

            int jVert;
            FOR (jVert, nVerts)
                if (VecMath.equals(3, vert3out, verts[jVert], tol))
                    break;
            if (jVert == nVerts)  // didn't find it
            {
                VecMath.copyvec(3, verts[nVerts], vert3out);
                verts[nVerts][3] = 1.;
                nVerts++;
            }
        }
        verts = (double[][])Arrays.subarray(verts, 0, nVerts);
        if (verboseLevel >= 1) OUT("      verts = "+Arrays.toStringNonCompact(verts, "        ", "    "));

        if (verboseLevel >= 1) OUT("    out generateVertices, returning "+verts.length+"/"+(vertGenerators.length*xforms.length)+" verts");
        return verts;
    } // generateVertices

} // class SymmetryUtils