/* * The mesh looks like this, where the centre triangle, #25, is always * filled. Vertices, edges, and triangles each have their own * numbering sequence, so there are three versions of this diagram, * one for each. (Trying to put them all on the same diagram proved * too visually confusing.) * * 0________ 1________ 2________ 3________ 4 * /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ * / \ / \ / \ / \ / \ * / \ / \ / \ / \ / \ * 5________ 6________ 7________ 8________ 9________10 * /\ /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ / \ * / \ / \ / \ / \ / \ / \ * / \ / \ / \ / \ / \ / \ * 11________12________13________14________15________16________17 * \ /\ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / \ / * \ / \ / \ / \ / \ / \ / * \ / \ / \ / \ / \ / \ / * 18________19________20________21________22________23 * \ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / * \ / \ / \ / \ / \ / * \ / \ / \ / \ / \ / * 24________25________26________27________28 * \ /\ /\ /\ / * \ / \ / \ / \ / * \ / \ / \ / \ / * \ / \ / \ / \ / * 29________30________31________32 * \ /\ /\ / * \ / \ / \ / * \ / \ / \ / * \ / \ / \ / * 33________34________35 * * _____0_________1_________2_________3____ * /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ * 4 5 6 7 8 9 10 11 12 13 * / \ / \ / \ / \ / \ * /___14___\/___15___\/___16___\/___17___\/___18___\ * /\ /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ / \ * 19 20 21 22 23 24 25 26 27 28 29 30 * / \ / \ / \ / \ / \ / \ * /___31___\/___32___\/___33___\/___34___\/___35___\/___36___\ * \ /\ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / \ / * 37 38 39 40 41 42 43 44 45 46 47 48 * \ / \ / \ / \ / \ / \ / * \/___49___\/___50___\/___51___\/___52___\/___53___\/ * \ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / * 54 55 56 57 58 59 60 61 62 63 * \ / \ / \ / \ / \ / * \/___64___\/___65___\/___66___\/___67___\/ * \ /\ /\ /\ / * \ / \ / \ / \ / * 68 69 70 71 72 73 74 75 * \ / \ / \ / \ / * \/___76___\/___77___\/___78___\/ * \ /\ /\ / * \ / \ / \ / * 79 80 81 82 83 84 * \ / \ / \ / * \/___85___\/___86___\/ * * ________________________________________ * /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ * / 0 \ 1 / 2 \ 3 / 4 \ 5 / 6 \ 7 / 8 \ * / \ / \ / \ / \ / \ * /________\/________\/________\/________\/________\ * /\ /\ /\ /\ /\ /\ * / \ / \ / \ / \ / \ / \ * / 9 \ 10 / 11 \ 12 / 13 \ 14 / 15 \ 16 / 17 \ 18 / 19 \ * / \ / \ / \ / \ / \ / \ * /________\/________\/________\/________\/________\/________\ * \ /\ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / \ / * \ 20 / 21 \ 22 / 23 \ 24 / 25 \ 26 / 27 \ 28 / 29 \ 30 / * \ / \ / \ / \ / \ / \ / * \/________\/________\/________\/________\/________\/ * \ /\ /\ /\ /\ / * \ / \ / \ / \ / \ / * \ 31 / 32 \ 33 / 34 \ 35 / 36 \ 37 / 38 \ 39 / * \ / \ / \ / \ / \ / * \/________\/________\/________\/________\/ * \ /\ /\ /\ / * \ / \ / \ / \ / * \ 40 / 41 \ 42 / 43 \ 44 / 45 \ 46 / * \ / \ / \ / \ / * \/________\/________\/________\/ * \ /\ /\ / * \ / \ / \ / * \ 47 / 48 \ 49 / 50 \ 51 / * \ / \ / \ / * \/________\/________\/ * * X and Y coordinates of points are in a coordinate system whose * origin is at the centre of gravity of triangle #25, whose X axis is * to the right in the above diagram, and in which the side of a * triangle has length 1. * * A shape is represented by a 52-bit unsigned integer, one bit per * triangle, listing the triangles making it up. Shapes that are * oriented differently, reflected, or translated lead to different * shapes at this level; redundancies due to translation, rotation, * and reflection are removed after generating the list. Triangle N * is represented by bit 1ULL << N. */ #include #include #include #include #define N_PTS 36 #define N_LNS 87 #define N_TRIS 52 #define MAX_SIZE 6 #define MAX_SHAPES 1000 #define SQRT_3_4 .86602540378443864676372317 /* sqrt(3/4) */ typedef struct pt PT; typedef struct ln LN; typedef struct tri TRI; typedef unsigned long long int SHAPE; struct pt { double x; double y; } ; struct ln { int pts[2]; } ; struct tri { int pts[3]; int lns[3]; int ngs[3]; double cx; double cy; } ; static PT pts[N_PTS] #define I(x,y) { (x)*.5, ((y)-(1.0/3))*SQRT_3_4 } = { I(-4, 3), I(-2, 3), I( 0, 3), I( 2, 3), I( 4, 3), /* 0 1 2 3 4 */ I(-5, 2), I(-3, 2), I(-1, 2), I( 1, 2), I( 3, 2), I( 5, 2), /* 5 6 7 8 9 10 */ I(-6, 1), I(-4, 1), I(-2, 1), I( 0, 1), I( 2, 1), I( 4, 1), I( 6, 1), /* 11 12 13 14 15 16 17 */ I(-5, 0), I(-3, 0), I(-1, 0), I( 1, 0), I( 3, 0), I( 5, 0), /* 18 19 20 21 22 23 */ I(-4,-1), I(-2,-1), I( 0,-1), I( 2,-1), I( 4,-1), /* 24 25 26 27 28 */ I(-3,-2), I(-1,-2), I( 1,-2), I( 3,-2), /* 29 30 31 32 */ I(-2,-3), I( 0,-3), I( 2,-3), /* 33 34 35 */ }; #undef I static LN lns[N_LNS] = { {{ 0, 1}}, {{ 1, 2}}, {{ 2, 3}}, {{ 3, 4}}, /* 0.. 3 */ {{ 0, 5}},{{ 0, 6}},{{ 1, 6}},{{ 1, 7}},{{ 2, 7}},{{ 2, 8}},{{ 3, 8}},{{ 3, 9}},{{ 4, 9}},{{ 4,10}}, /* 4..13 */ {{ 5, 6}}, {{ 6, 7}}, {{ 7, 8}}, {{ 8, 9}}, {{ 9,10}}, /* 14..18 */ {{ 5,11}},{{ 5,12}},{{ 6,12}},{{ 6,13}},{{ 7,13}},{{ 7,14}},{{ 8,14}},{{ 8,15}},{{ 9,15}},{{ 9,16}},{{10,16}},{{10,17}}, /* 19..30 */ {{11,12}}, {{12,13}}, {{13,14}}, {{14,15}}, {{15,16}}, {{16,17}}, /* 31..36 */ {{11,18}},{{12,18}},{{12,19}},{{13,19}},{{13,20}},{{14,20}},{{14,21}},{{15,21}},{{15,22}},{{16,22}},{{16,23}},{{17,23}}, /* 37..48 */ {{18,19}}, {{19,20}}, {{20,21}}, {{21,22}}, {{22,23}}, /* 49..53 */ {{18,24}},{{19,24}},{{19,25}},{{20,25}},{{20,26}},{{21,26}},{{21,27}},{{22,27}},{{22,28}},{{23,28}}, /* 54..63 */ {{24,25}}, {{25,26}}, {{26,27}}, {{27,28}}, /* 64..67 */ {{24,29}},{{25,29}},{{25,30}},{{26,30}},{{26,31}},{{27,31}},{{27,32}},{{28,32}}, /* 68..75 */ {{29,30}}, {{30,31}}, {{31,32}}, /* 76..78 */ {{29,33}},{{30,33}},{{30,34}},{{31,34}},{{31,35}},{{32,35}}, /* 79..84 */ {{33,34}}, {{34,35}}, /* 85..86 */ }; static TRI tris[N_TRIS] = { {{ 0, 5, 6}}, {{ 0, 1, 6}}, {{ 1, 6, 7}}, {{ 1, 2, 7}}, {{ 2, 7, 8}}, {{ 2, 3, 8}}, {{ 3, 8, 9}}, {{ 3, 4, 9}}, {{ 4, 9,10}}, /* 0.. 8 */ {{ 5,11,12}}, {{ 5, 6,12}}, {{ 6,12,13}}, {{ 6, 7,13}}, {{ 7,13,14}}, {{ 7, 8,14}}, {{ 8,14,15}}, {{ 8, 9,15}}, {{ 9,15,16}}, {{ 9,10,16}}, {{10,16,17}}, /* 9..19 */ {{11,12,18}}, {{12,18,19}}, {{12,13,19}}, {{13,19,20}}, {{13,14,20}}, {{14,20,21}}, {{14,15,21}}, {{15,21,22}}, {{15,16,22}}, {{16,22,23}}, {{16,17,23}}, /* 20..30 */ {{18,19,24}}, {{19,24,25}}, {{19,20,25}}, {{20,25,26}}, {{20,21,26}}, {{21,26,27}}, {{21,22,27}}, {{22,27,28}}, {{22,23,28}}, /* 31..39 */ {{24,25,29}}, {{25,29,30}}, {{25,26,30}}, {{26,30,31}}, {{26,27,31}}, {{27,31,32}}, {{27,28,32}}, /* 40..46 */ {{29,30,33}}, {{30,33,34}}, {{30,31,34}}, {{31,34,35}}, {{31,32,35}}, }; /* 47..51 */ #define MID_TRI 25 static SHAPE shapes[MAX_SHAPES]; static char shapedup[MAX_SHAPES]; static int nshapes; #define THASHSIZE 66 static int thash[THASHSIZE]; static int find_ln(int p1, int p2) { int i; for (i=N_LNS-1;i>=0;i--) { if ( ((lns[i].pts[0] == p1) && (lns[i].pts[1] == p2)) || ((lns[i].pts[1] == p1) && (lns[i].pts[0] == p2)) ) return(i); } abort(); } static int find_ng(int p1, int p2, int exc) { int i; for (i=N_TRIS-1;i>=0;i--) { if (i == exc) continue; if ( ((tris[i].pts[0] == p1) && (tris[i].pts[1] == p2)) || ((tris[i].pts[1] == p1) && (tris[i].pts[0] == p2)) || ((tris[i].pts[1] == p1) && (tris[i].pts[2] == p2)) || ((tris[i].pts[2] == p1) && (tris[i].pts[1] == p2)) || ((tris[i].pts[2] == p1) && (tris[i].pts[0] == p2)) || ((tris[i].pts[0] == p1) && (tris[i].pts[2] == p2)) ) return(i); } return(-1); } static int hash_tri(double cx, double cy) { return((11*(3+(int)rint(cy/SQRT_3_4)))+5+(int)rint(cx/.5)); } static void gen_mesh(void) { int i; double d; int errs; int h; errs = 0; for (i=N_LNS-1;i>=0;i--) { d = hypot( pts[lns[i].pts[0]].x - pts[lns[i].pts[1]].x, pts[lns[i].pts[0]].y - pts[lns[i].pts[1]].y ); if ((d < .99) || (d > 1.01)) { printf("ln %d: %d-%d: d=%g\n",i,lns[i].pts[0],lns[i].pts[1],d); errs = 1; } } for (i=THASHSIZE-1;i>=0;i--) thash[i] = -1; for (i=N_TRIS-1;i>=0;i--) { tris[i].lns[0] = find_ln(tris[i].pts[0],tris[i].pts[1]); tris[i].lns[1] = find_ln(tris[i].pts[0],tris[i].pts[2]); tris[i].lns[2] = find_ln(tris[i].pts[1],tris[i].pts[2]); tris[i].ngs[0] = find_ng(tris[i].pts[1],tris[i].pts[2],i); tris[i].ngs[1] = find_ng(tris[i].pts[0],tris[i].pts[2],i); tris[i].ngs[2] = find_ng(tris[i].pts[0],tris[i].pts[1],i); tris[i].cx = (pts[tris[i].pts[0]].x + pts[tris[i].pts[1]].x + pts[tris[i].pts[2]].x) / 3; tris[i].cy = (pts[tris[i].pts[0]].y + pts[tris[i].pts[1]].y + pts[tris[i].pts[2]].y) / 3; h = hash_tri(tris[i].cx,tris[i].cy); if ((h < 0) || (h >= THASHSIZE) || (thash[h] >= 0)) abort(); thash[h] = i; } if (errs) exit(1); } static void print_shape(SHAPE s) { int i; const char *pre; printf("<"); pre = ""; for (i=0;i"); } static void gen_shapes(void) { int size; int prevmin; int prevmax; int curmin; int curmax; int sx; SHAPE s; int tx; int nx; int nt; SHAPE ns; int i; shapes[0] = 1ULL << MID_TRI; nshapes = 1; curmin = 0; curmax = 0; for (size=2;size<=MAX_SIZE;size++) { prevmin = curmin; prevmax = curmax; /* printf("size=%d (prev min=%d max=%d)\n",size,prevmin,prevmax); */ curmin = MAX_SHAPES; curmax = -1; for (sx=prevmin;sx<=prevmax;sx++) { s = shapes[sx]; /* printf("shape "); print_shape(s); printf(" from %d\n",sx); */ for (tx=N_TRIS-1;tx>=0;tx--) { if (s & (1ULL << tx)) { for <"ngh"> (nx=3-1;nx>=0;nx--) { nt = tris[tx].ngs[nx]; if (nt < 0) continue; if (s & (1ULL << nt)) continue; ns = s | (1ULL << nt); /* printf("tx=%d nx=%d nt=%d ns=",tx,nx,nt); */ /* print_shape(ns);*/ for (i=curmin;i<=curmax;i++) { if (shapes[i] == ns) { /* printf(" duplicate of %d\n",i); */ continue <"ngh">; } } /* printf(" saving as %d\n",nshapes); */ if (nshapes >= MAX_SHAPES) abort(); shapes[nshapes] = ns; if (nshapes < curmin) curmin = nshapes; if (nshapes > curmax) curmax = nshapes; nshapes ++; } } } } /* printf("size loop done, min=%d max=%d\n",curmin,curmax); */ } } /* * Compute the transformation of the plane that takes fp1 to tp1, fp2 * to tp2, and fp3 to tp3. Transform s with this transformation; if * the result is not equal to s, call (*fn) with the new SHAPE. * * To do this, we solve for M and b in * * [tp1.x] [ Mxx Mxy ] [ fp1.x ] [b.x] * [ ] = [ ] [ ] + [ ] * [tp1.y] [ Myx Myy ] [ fp1.y ] [b.y] * * [tp2.x] [ Mxx Mxy ] [ fp2.x ] [b.x] * [ ] = [ ] [ ] + [ ] * [tp2.y] [ Myx Myy ] [ fp2.y ] [b.y] * * [tp3.x] [ Mxx Mxy ] [ fp3.x ] [b.x] * [ ] = [ ] [ ] + [ ] * [tp3.y] [ Myx Myy ] [ fp3.y ] [b.y] * * This gives six equations in six unknowns: * * [1 0 0 0 0 0] [tp1.x] [fp1.x fp1.y 0 0 1 0 ] [Mxx] * [0 1 0 0 0 0] [tp1.y] [ 0 0 fp1.x fp1.y 0 1 ] [Mxy] * [0 0 1 0 0 0] [tp2.x] = [fp2.x fp2.y 0 0 1 0 ] [Myx] * [0 0 0 1 0 0] [tp2.y] [ 0 0 fp2.x fp2.y 0 1 ] [Myy] * [0 0 0 0 1 0] [tp3.x] [fp3.x fp3.y 0 0 1 0 ] [b.x] * [0 0 0 0 0 1] [tp3.y] [ 0 0 fp3.x fp3.y 0 1 ] [b.y] * * We solve this with traditional Gauss-Jordan elimination. Because * of the way our arguments are generated, it is a can't-happen for * the matrix on the RHS above to be anywhere near ill-determined. * * Once we have this M and b, we identify the transforms of all the * triangles in the shape, by transforming the centre of the source * triangle according to M and b, then using hash_tri on it. If the * result is outside thash, if the thash bucket indicated is empty, or * if the centre of the triangle indicated by the thash[] bucket is * more than a small distance from the transformed centre * (mathematically zero, but numerical errors occur), then the * transformed triangle is outside the mesh; this is a can't-happen. * Otherwise, this gives us a list of triangles making up the * transformed piece. It is another can't-happen for any particular * triangle to occur more than once in this list. We construct the * corresponding SHAPE; if it is different from s, we pass it to * (*fn)(). */ static void maybe_call_transformed( void (*fn)(SHAPE), SHAPE s, int fp1, int fp2, int fp3, int tp1, int tp2, int tp3 ) { double eq[12][6]; int x; int y; int my; double p; int i; double t; double Mb[6]; double tpv[6]; int tx; double fcx; double fcy; double tcx; double tcy; int ttx; int hb; SHAPE s2; static void dump_eq(const char *tag __attribute__((__unused__))) { /* int xi; int yi; printf("%s\n",tag); for (yi=0;yi<6;yi++) { for (xi=0;xi<12;xi++) { printf(" %10.5f",eq[xi][yi]); } printf("\n"); } */ } for (x=12-1;x>=0;x--) for (y=6-1;y>=0;y--) eq[x][y] = 0; eq[0][0] = pts[fp1].x; eq[2][1] = pts[fp1].x; eq[1][0] = pts[fp1].y; eq[3][1] = pts[fp1].y; eq[0][2] = pts[fp2].x; eq[2][3] = pts[fp2].x; eq[1][2] = pts[fp2].y; eq[3][3] = pts[fp2].y; eq[0][4] = pts[fp3].x; eq[2][5] = pts[fp3].x; eq[1][4] = pts[fp3].y; eq[3][5] = pts[fp3].y; eq[4][0] = 1; eq[5][1] = 1; eq[4][2] = 1; eq[5][3] = 1; eq[4][4] = 1; eq[5][5] = 1; for (x=6-1;x>=0;x--) eq[x+6][x] = 1; dump_eq("initial"); for (x=6-1;x>=0;x--) { my = x; for (y=x-1;y>=0;y--) if (fabs(eq[x][y]) > fabs(eq[x][my])) my = y; p = eq[x][my]; /* printf("x=%d my=%d p=%g\n",x,my,p); */ if (fabs(p) < .0001) abort(); if (my != x) { for (i=12-1;i>=0;i--) { t = eq[i][my]; eq[i][my] = eq[i][x]; eq[i][x] = t; } dump_eq("after row swap"); } for (i=12-1;i>=0;i--) eq[i][x] /= p; dump_eq("after dividing"); for (y=6-1;y>=0;y--) { if (y != x) { t = eq[x][y]; for (i=12-1;i>=0;i--) eq[i][y] -= t * eq[i][x]; } } dump_eq("after subtracting"); } tpv[0] = pts[tp1].x; tpv[1] = pts[tp1].y; tpv[2] = pts[tp2].x; tpv[3] = pts[tp2].y; tpv[4] = pts[tp3].x; tpv[5] = pts[tp3].y; for (y=6-0;y>=0;y--) { t = 0; for (x=6-1;x>=0;x--) t += tpv[x] * eq[x+6][y]; Mb[y] = t; } /* printf("Mb = %g %g %g %g %g %g\n",Mb[0],Mb[1],Mb[2],Mb[3],Mb[4],Mb[5]); */ s2 = 0; for (tx=N_TRIS-1;tx>=0;tx--) { if (s & (1ULL << tx)) { fcx = tris[tx].cx; fcy = tris[tx].cy; tcx = (Mb[0] * fcx) + (Mb[1] * fcy) + Mb[4]; tcy = (Mb[2] * fcx) + (Mb[3] * fcy) + Mb[5]; hb = hash_tri(tcx,tcy); if ((hb < 0) || (hb >= THASHSIZE) || (thash[hb] < 0)) abort(); ttx = thash[hb]; if ((fabs(tris[ttx].cx-tcx) > .01) || (fabs(tris[ttx].cy-tcy) > .01)) abort(); if (s2 & (1ULL << ttx)) abort(); s2 |= 1ULL << ttx; } } if (s2 != s) (*fn)(s2); } static void map_varieties(SHAPE s, void (*fn)(SHAPE)) { int t; for (t=N_TRIS-1;t>=0;t--) { if (s & (1ULL << t)) { maybe_call_transformed( fn, s, tris[t].pts[0], tris[t].pts[1], tris[t].pts[2], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); maybe_call_transformed( fn, s, tris[t].pts[0], tris[t].pts[2], tris[t].pts[1], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); maybe_call_transformed( fn, s, tris[t].pts[1], tris[t].pts[2], tris[t].pts[0], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); maybe_call_transformed( fn, s, tris[t].pts[1], tris[t].pts[0], tris[t].pts[2], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); maybe_call_transformed( fn, s, tris[t].pts[2], tris[t].pts[0], tris[t].pts[1], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); maybe_call_transformed( fn, s, tris[t].pts[2], tris[t].pts[1], tris[t].pts[0], tris[MID_TRI].pts[0], tris[MID_TRI].pts[1], tris[MID_TRI].pts[2] ); } } } static void filter_shapes(void) { int i; int j; SHAPE s; void drop(SHAPE s2) { int k; for (k=i+1;k