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 <stdio.h>
32 #include <stdlib.h>
33 #include <math.h>
34 #include "vec4.h"
35 #include "enum.h"
36 #include "transform.h"
37 #include "dgflag.h"
38
39
40 extern int debug, stringent, metric;
41 static float epsilon = .005;
42
43 void traverse_list();
44 int is_same();
45
46
47 struct node {
48 Transform t;
49 struct node *l,*r, *p; /* left, right, parent */
50 struct node *same; /* if it's "close" to this one */
51 int num;
52 int type; /* process but don't print? */
53 float norm;
54 char *word;
55 } ;
56
57 static struct node *head;
58
59 double
getnorm(int metric,Transform t)60 getnorm(int metric, Transform t)
61 {
62 switch (metric) {
63 case DG_EUCLIDEAN:
64 return(sqrt(t[3][0]*t[3][0] + t[3][1]*t[3][1] + t[3][2]*t[3][2]));
65 break;
66 case DG_HYPERBOLIC:
67 if (ABS(t[3][3]) < 1.0) return(0.0);
68 return ( acosh( ABS( t[3][3])) );
69 break;
70 case DG_SPHERICAL:
71 {
72 float sum = 0;
73 int i,j;
74 for (i=0; i<4; ++i)
75 for (j=0; j<4; ++j)
76 sum += fabs(t[i][j] - (i==j));
77 return(sum);
78 }
79 break;
80 }
81 return (double)0;
82 }
83
84 static
85 struct node *
alloc_node()86 alloc_node()
87 {
88 struct node *n;
89 n = (struct node *) malloc(sizeof(struct node) );
90 if (n == NULL) exit(printf("Unable to allocate: alloc_node\n"));
91 n->p = n->l = n->r = n->same = NULL;
92 n->num = 1;
93 return(n);
94 }
95
96 int
insert_or_match_mat(mat,mode)97 insert_or_match_mat(mat, mode)
98 Transform mat;
99 int mode; /* insert or match? */
100 {
101 struct node *n = NULL, *p;
102 float d;
103 struct node tnode;
104 double getnorm();
105
106 if (debug == 4) traverse_list(head);
107
108 if (mode & INSERT ) n = alloc_node();
109 else if (mode & MATCH) n = &tnode;
110
111 else n->type = 0;
112
113 n->norm = getnorm(metric, mat);
114
115 TmCopy(mat, n->t);
116
117 if (head == NULL)
118 {
119 if (mode & MATCH) return(0); /* no match */
120 else if (mode & INSERT) {
121 head = n;
122 return(1); /* successful insert */
123 }
124 }
125 /* ...else... */
126 for (p = head; p != NULL; ) {
127 d = fabs(p->norm - n->norm);
128 if (ABS(d) < epsilon) { /* insert! */
129 if (mode & INSERT) {
130 /* and stick it at the end of the cluster */
131 p->num++; /* count this one */
132 for ( ; p->same != NULL; p = p->same) ;
133 /* back up one */
134 p->same = n;
135 n->p = p; /* mark the parent */
136 return(1);
137 }
138 else { /* look here for a match */
139 for ( ; p != NULL; p = p->same)
140 if (is_same(p->t, n->t))
141 return( 1 ); /* matched! */
142 return(0);
143 }
144 }
145 else if (d > 0) { /* go right */
146 if (p->r == NULL) {
147 if (mode & INSERT) {
148 p->r = n;
149 n->p = p;
150 return(1);
151 }
152 else { /* no match */
153 return(0);
154 }
155 }
156 else p = p->r; /* continue search */
157 }
158 else if (d < 0) { /* go left */
159 if (p->l == NULL) {
160 if (mode & INSERT) {
161 p->l = n;
162 n->p = p;
163 return(1);
164 }
165 else { /* no match */
166 return(0);
167 }
168 }
169 else p = p->l; /* continue search */
170 }
171 }
172 return 0;
173 }
174
175 int
is_same(t0,t1)176 is_same(t0, t1)
177 Transform t0, t1;
178 {
179 int i, j, same = 1;
180
181 if (stringent) {
182 float factor, fepsilon;
183 Transform tt1, tt2;
184 Tm3Invert(t0, tt1);
185 Tm3Concat(t1,tt1,tt2);
186 /* is tt2 the identity, or a multiple of the identity? */
187 factor = tt2[0][0];
188 fepsilon = ABS(factor * epsilon);
189 for (i=0; i<4; ++i)
190 for (j=0; j<4; ++j)
191 {
192 /* check against identity matrix */
193 if ( ABS( tt2[i][j] - factor * ((i == j) ? 1 : 0) ) > fepsilon )
194 {
195 same = 0;
196 goto OUT;
197 }
198 }
199 }
200 else {
201 for (i=0; i<4; ++i)
202 {
203 for (j=0; j<4; ++j)
204 if (ABS(t0[i][j] - t1[i][j]) > epsilon)
205 {
206 same = 0;
207 goto OUT;
208 }
209 }
210 }
211 OUT:
212 return(same);
213 }
214
215 int
is_new(t)216 is_new(t)
217 Transform t;
218 {
219 if ( insert_or_match_mat(t, MATCH) )
220 {
221 return(0);
222 }
223 return(DG_CONSTRAINT_NEW);
224 }
225
226
227 void
traverse_list(n)228 traverse_list(n)
229 struct node *n;
230 {
231 if (n == NULL) return;
232 traverse_list(n->l);
233 fprintf(stderr,"%10f\t%d\n",n->t[3][3], n->num);
234 traverse_list(n->r);
235 }
236
237 static void
_delete_list(n)238 _delete_list(n)
239 struct node *n;
240 {
241 struct node *nt, *ot;
242 if (n == NULL) return;
243 _delete_list(n->l);
244 _delete_list(n->r);
245 /* erase the cluster of equal values at this node */
246 /* first go to the end of the cluster ... */
247 for (ot = n ; ot->same != NULL; ot = ot->same) ;
248 /* ... then climb back up, erasing as you go */
249 if (ot != n) for ( ; ot->p != n; ) {
250 nt = ot->p;
251 free(ot);
252 ot = nt;
253 }
254 /* erase the evidence of this node in its parent */
255 if (n->p) {
256 if (n->p->l == n) n->p->l = NULL;
257 if (n->p->r == n) n->p->r = NULL;
258 }
259 free(n);
260 }
261
262 void
delete_list()263 delete_list()
264 {
265 _delete_list(head);
266 head = NULL;
267 }
268
269