1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #if 0
27 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
28 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
29 #endif
30
31 #include "geom.h"
32 #include "polylistP.h"
33 #include "discgrpP.h"
34 #include "point.h"
35 #include "winged_edge.h"
36 #include "transform.h"
37 #include "math.h"
38
39 static WEpolyhedron *wepoly1, **wepoly2;
40 HPoint3 DGorigin = {0,0,0,1};
41 HPoint3 DGrandom = {.1,.2,.3,.4};
42
43 extern void do_weeks_code();
44
45 static int
is_id(Transform t)46 is_id(Transform t)
47 {
48 int i,j;
49
50 for (i=0; i<4; ++i)
51 for (j=0; j<4; ++j)
52 if (fabs(t[i][j] - (i == j)) > .0005) return(0);
53 return(1);
54 }
55
56 /* make sure that the center point attached to the discrete group
57 isn't fixed by some generator; if it is; redefine the center
58 point to be the center of gravity of the orbit by the generators
59 */
60 void
DiscGrpCheckCPoint(DiscGrp * dg)61 DiscGrpCheckCPoint(DiscGrp *dg)
62 {
63 int i, cnt, fixed;
64 HPoint3 tmp, out;
65 float d;
66
67 if (dg->gens == NULL) {
68 return;
69 }
70
71 /* go through the generators, checking for fixed-pointed-ness */
72 for (i = 0, fixed = 0 ; i < dg->gens->num_el; ++i )
73 {
74 HPt3Transform(dg->gens->el_list[i].tform, &dg->cpoint, &tmp);
75 d = HPt3SpaceDistance(&dg->cpoint, &tmp, dg->attributes & DG_METRIC_BITS);
76 if (fabs(d) < .0005) {
77 fixed = 1;
78 break;
79 }
80 }
81
82 /* no fixed points */
83 if (fixed == 0) return;
84
85 /* clean out the special bit */
86 for (i = 0 ; i < dg->gens->num_el; ++i )
87 dg->gens->el_list[i].attributes &= ~DG_TMP;
88
89 out.x = out.y = out.z = out.w = 0.0;
90 /* don't average but one of each generator, inverse pair */
91 for (cnt = 0, i = 0 ; i < dg->gens->num_el; ++i ) {
92 if (!(dg->gens->el_list[i].attributes & DG_TMP)) {
93 HPt3Transform(dg->gens->el_list[i].tform, &DGrandom, &tmp);
94 HPt3Add(&tmp, &out, &out);
95 dg->gens->el_list[i].inverse->attributes |= DG_TMP;
96 cnt++;
97 }
98 }
99 HPt3Dehomogenize(&out, &out);
100 /* return it or set cpoint?? */
101 dg->cpoint = out;
102 }
103
104 void
DiscGrpAddInverses(DiscGrp * discgrp)105 DiscGrpAddInverses(DiscGrp *discgrp)
106 {
107 int i, j, found = 0;
108 Transform product;
109 DiscGrpElList *newgens;
110
111 /* remove all identity matrices */
112 for (j=0, i=0; i<discgrp->gens->num_el; ++i) {
113 if ( !is_id(discgrp->gens->el_list[i].tform) ) {
114 /* ought to have a DiscGrpElCopy() */
115 discgrp->gens->el_list[j] = discgrp->gens->el_list[i];
116 TmCopy(discgrp->gens->el_list[i].tform,
117 discgrp->gens->el_list[j].tform);
118 j++;
119 }
120 }
121 /* store the new count */
122 discgrp->gens->num_el = j;
123
124 for (i=0; i<discgrp->gens->num_el; ++i) {
125 if (discgrp->gens->el_list[i].inverse == NULL) {
126 /* look for inverse among the existing generators */
127 for (j=i; j<discgrp->gens->num_el; ++j) {
128 TmConcat(discgrp->gens->el_list[i].tform, discgrp->gens->el_list[j].tform, product);
129 if ( is_id(product) ) {
130 discgrp->gens->el_list[i].inverse = &discgrp->gens->el_list[j];
131 discgrp->gens->el_list[j].inverse = &discgrp->gens->el_list[i];
132 found++;
133 }
134 }
135 }
136 else found++;
137 }
138
139 newgens = OOGLNew(DiscGrpElList);
140 newgens->num_el = 2 * discgrp->gens->num_el - found;
141 newgens->el_list = OOGLNewN(DiscGrpEl, newgens->num_el );
142
143 memcpy(newgens->el_list, discgrp->gens->el_list,
144 sizeof(DiscGrpEl) * discgrp->gens->num_el);
145
146 /* now go through looking for group elements without inverses */
147 {
148 char c;
149 j = discgrp->gens->num_el;
150 for (i=0; i<discgrp->gens->num_el; ++i) {
151 if (newgens->el_list[i].inverse == NULL) {
152 newgens->el_list[j+i] = newgens->el_list[i];
153 /* make the symbol of the inverse be the 'other case' */
154 c = newgens->el_list[i].word[0];
155 if (c < 'a') newgens->el_list[j+i].word[0] = c + 32;
156 else newgens->el_list[j+i].word[0] = c - 32;
157 TmInvert( newgens->el_list[i].tform, newgens->el_list[j+i].tform);
158 newgens->el_list[j+i].inverse = &newgens->el_list[i];
159 newgens->el_list[i].inverse = &newgens->el_list[j+i];
160 }
161 else j--;
162 }
163 }
164
165 DiscGrpElListDelete(discgrp->gens);
166 discgrp->gens = newgens;
167 }
168
169 void
DiscGrpSetupDirdom(DiscGrp * discgrp)170 DiscGrpSetupDirdom(DiscGrp *discgrp)
171 {
172 WEpolyhedron *dd;
173
174 if ( discgrp->nhbr_list ) {
175 OOGLFree(discgrp->nhbr_list->el_list);
176 OOGLFree(discgrp->nhbr_list);
177 }
178
179 /* worry about fixed points */
180 DiscGrpCheckCPoint(discgrp);
181 dd = DiscGrpMakeDirdom(discgrp, &discgrp->cpoint, 0);
182 discgrp->nhbr_list = DiscGrpExtractNhbrs(dd);
183 }
184
185 /* find the group element whose 'center point' is closest to the point poi */
186 DiscGrpEl *
DiscGrpClosestGroupEl(DiscGrp * discgrp,HPoint3 * poi)187 DiscGrpClosestGroupEl(DiscGrp *discgrp, HPoint3 *poi)
188 {
189 int count, i, closeri;
190 int metric;
191 float min = 0, d;
192 HPoint3 pt0, pt1;
193 DiscGrpEl *closer = NULL, *closest = OOGLNew(DiscGrpEl);
194 Transform cinv;
195
196 TmIdentity(closest->tform);
197 closest->attributes = 0;
198
199 if (!discgrp->nhbr_list) DiscGrpSetupDirdom(discgrp);
200
201 metric = discgrp->attributes & (DG_METRIC_BITS);
202
203 /* iterate until we're in the fundamental domain */
204 count = 0;
205 closeri = -1;
206 pt0 = *poi;
207 while (count < 1000 && closeri != 0) {
208 for (i = 0; i<discgrp->nhbr_list->num_el; ++i) {
209 HPt3Transform(discgrp->nhbr_list->el_list[i].tform, &discgrp->cpoint, &pt1);
210 d = HPt3SpaceDistance(&pt0, &pt1, metric);
211 if (i==0) {
212 min = d;
213 closer = &discgrp->nhbr_list->el_list[i];
214 closeri = i;
215 }
216 else if (d < min) {
217 min = d;
218 closer = &discgrp->nhbr_list->el_list[i];
219 closeri = i;
220 }
221 }
222 count++;
223 if (closeri) {
224 TmConcat(closer->tform, closest->tform, closest->tform);
225 /* move the point of interest by the inverse of the closest nhbr
226 and iterate */
227 TmInvert(closest->tform, cinv);
228 HPt3Transform(cinv, poi, &pt0);
229 }
230 /*
231 if (needstuneup(closest->tform))
232 tuneup(closest->tform, metric);
233 */
234 }
235 if (is_id(closest->tform)) closest->attributes |= DGEL_IS_IDENTITY;
236 return (closest);
237 }
238 static ColorA white = {1,1,1,1};
239
240 /* return a list of group elements corresponding to the faces of the
241 dirichlet domain */
242 DiscGrpElList *
DiscGrpExtractNhbrs(WEpolyhedron * wepoly)243 DiscGrpExtractNhbrs( WEpolyhedron *wepoly )
244 {
245 int i,j,k;
246 WEface *fptr;
247 DiscGrpElList *mylist;
248 ColorA GetCmapEntry();
249
250 if (!wepoly) return(NULL);
251
252 /* should use realloc() here to take care of large groups...*/
253 mylist = OOGLNew(DiscGrpElList);
254 mylist->el_list = OOGLNewN(DiscGrpEl, wepoly->num_faces + 1);
255 mylist->num_el = wepoly->num_faces + 1;
256
257 /* include the identity matrix */
258 TmIdentity( mylist->el_list[0].tform);
259 mylist->el_list[0].color = white;
260 mylist->el_list[0].attributes = DGEL_IS_IDENTITY;
261
262 /* read the matrices corresponding to the faces of dirichlet domain */
263 for (fptr = wepoly->face_list, k = 1;
264 k<=wepoly->num_faces && fptr != NULL;
265 k++, fptr = fptr->next)
266 {
267 for (i=0; i<4; ++i)
268 for (j=0; j<4; ++j)
269 /* the group elements stored in fptr are transposed! */
270 mylist->el_list[k].tform[j][i] = fptr->group_element[i][j];
271 mylist->el_list[k].color = GetCmapEntry(fptr->fill_tone);
272 }
273 if (mylist->num_el != k)
274 OOGLError(1,"Incorrect number of nhbrs.\n");;
275
276 return(mylist);
277 }
278
279 /* attempt to create a scaled copy of fundamental domain for spherical
280 groups by taking weighted sum of vertices with the distinguished point */
281 static void
DiscGrpScalePolyList(DiscGrp * dg,PolyList * dirdom,HPoint3 * pt0,float scale)282 DiscGrpScalePolyList(DiscGrp *dg, PolyList *dirdom, HPoint3 *pt0, float scale)
283 {
284 int i, metric;
285 HPoint3 tmp1, tmp2, tpt0, tpt1;
286 HPt3Copy(pt0, &tpt0);
287 metric = dg->attributes & DG_METRIC_BITS;
288 if (metric != DG_EUCLIDEAN) {
289 HPt3SpaceNormalize(&tpt0, metric);
290 HPt3Scale( 1.0 - scale, &tpt0, &tmp2);
291 for (i=0; i<dirdom->n_verts; ++i) {
292 HPt3Copy(&dirdom->vl[i].pt, &tpt1);
293 HPt3SpaceNormalize(&tpt1, metric);
294 HPt3SpaceNormalize(&tpt1, metric);
295 HPt3Scale( scale, &tpt1, &tmp1);
296 HPt3Add(&tmp1, &tmp2, &dirdom->vl[i].pt);
297 }
298 }
299 else {
300 Transform T, TT, ITT, tmp;
301 static HPoint3 average = {0,0,0,0};
302 /* compute average */
303 for (i=0; i<dirdom->n_verts; ++i)
304 HPt3Add(&average, &dirdom->vl[i].pt, &average);
305 HPt3Dehomogenize(&average, &average);
306
307 TmTranslate(TT, average.x, average.y, average.z );
308 TmInvert(TT, ITT);
309 TmScale(T, scale, scale, scale);
310 TmConcat(ITT, T, tmp);
311 TmConcat(tmp, TT, tmp);
312 for (i=0; i<dirdom->n_verts; ++i)
313 HPt3Transform(tmp, &dirdom->vl[i].pt, &dirdom->vl[i].pt);
314 }
315 }
316
317 Geom *small_dd, *large_dd;
318 Geom *
DiscGrpDirDom(DiscGrp * dg)319 DiscGrpDirDom(DiscGrp *dg)
320 {
321 Geom *oogldirdom;
322 WEpolyhedron *dd;
323 extern Geom *WEPolyhedronToPolyList();
324 Geom *mylist, *smlist;
325
326 if (dg->flag & DG_DDBEAM) {
327 WEpolyhedron *poly = DiscGrpMakeDirdom( dg, &dg->cpoint, 0);
328 Geom *beams;
329 beams = WEPolyhedronToBeams(poly, dg->scale);
330 return(beams);
331 }
332 else {
333 float scale;
334 /* first a full-size, wireframe version of dd */
335 dd = DiscGrpMakeDirdom(dg, &dg->cpoint, 0);
336 if (dd) {
337 oogldirdom = WEPolyhedronToPolyList(dd);
338 scale = 1.0;
339 DiscGrpScalePolyList(dg, (PolyList *)oogldirdom, &dg->cpoint, scale);
340 large_dd = oogldirdom;
341 large_dd->ap = ApCreate(AP_DO, APF_EDGEDRAW, AP_DONT, APF_FACEDRAW, AP_END);
342 }
343 else return((Geom *) NULL);
344 /* next a scaled version with cusps cut off */
345 dd = DiscGrpMakeDirdom(dg, &dg->cpoint, 1);
346 if (dd) {
347 oogldirdom = WEPolyhedronToPolyList(dd);
348 DiscGrpScalePolyList(dg, (PolyList *)oogldirdom, &dg->cpoint, dg->scale);
349 small_dd = oogldirdom;
350 small_dd->ap = ApCreate(AP_DONT, APF_EDGEDRAW, AP_DO, APF_FACEDRAW, AP_END);
351 }
352 else return((Geom *) NULL);
353 smlist = GeomCreate("list", CR_GEOM, small_dd, CR_END);
354 mylist = GeomCreate("list", CR_GEOM, large_dd, CR_CDR, smlist, CR_END);
355 return(mylist);
356 }
357 }
358
359 WEpolyhedron *
DiscGrpMakeDirdom(DiscGrp * gamma,HPoint3 * poi,int slice)360 DiscGrpMakeDirdom(DiscGrp *gamma, HPoint3 *poi, int slice)
361 {
362 int i, j, k;
363 proj_matrix *gen_list;
364 point origin;
365 int metric, transp;
366
367 transp = gamma->attributes & DG_TRANSPOSED;
368 /* transform from floating point to double, essentially */
369 gen_list = OOGLNewNE(proj_matrix, gamma->gens->num_el, "DiscGrp gens");
370 /* jeff week's code basically uses transposed matrices, so
371 if the transposed flag is set, we do nothing, otherwise we
372 have to transpose! */
373 for (i=0; i<gamma->gens->num_el; ++i)
374 for (j=0; j<4; ++j)
375 for (k=0; k<4; ++k)
376 {
377 if (transp) gen_list[i][j][k] = gamma->gens->el_list[i].tform[j][k];
378 else gen_list[i][k][j] = gamma->gens->el_list[i].tform[j][k];
379 }
380 origin[0] = poi->x;
381 origin[1] = poi->y;
382 origin[2] = poi->z;
383 origin[3] = poi->w;
384
385 wepoly2 = &wepoly1;
386 metric = (gamma->attributes & DG_METRIC_BITS);
387 do_weeks_code(wepoly2, origin, gen_list, gamma->gens->num_el, metric,slice);
388
389 OOGLFree(gen_list);
390
391 /* turn off the 'new dirdom' bit */
392 gamma->flag &= ~DG_NEWDIRDOM;
393 return(*wepoly2);
394 }
395