1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 * Copyright (C) 2006 Claus-Justus Heine
4 *
5 * This file is part of Geomview.
6 *
7 * Geomview is free software; you can redistribute it and/or modify it
8 * under the terms of the GNU Lesser General Public License as published
9 * by the Free Software Foundation; either version 2, or (at your option)
10 * any later version.
11 *
12 * Geomview is distributed in the hope that it will be useful, but
13 * WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
15 * Lesser General Public License for more details.
16 *
17 * You should have received a copy of the GNU Lesser General Public
18 * License along with Geomview; see the file COPYING. If not, write
19 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
20 * USA, or visit http://www.gnu.org.
21 */
22
23 /* Authors: Charlie Gunn, Pat Hanrahan, Stuart Levy, Tamara Munzner,
24 * Mark Phillips, Claus-Justus Heine */
25
26 #ifndef _GV_HPOINTN_H_
27 #define _GV_HPOINTN_H_
28
29 #if HAVE_CONFIG_H
30 # include "config.h"
31 #endif
32
33 #if 0
34 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
35 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
36 #endif
37
38 #include <stdlib.h>
39 #include <math.h>
40 #include <ooglutil.h>
41 #include "geomtypes.h"
42 #include "freelist.h"
43
44 extern DEF_FREELIST(HPointN);
45
46 static inline HPointN *HPtNCreate(int dim, const HPtNCoord *vec);
47 static inline void HPtNDelete(HPointN *pt);
48 static inline HPointN *HPtNCopy(const HPointN *pt1, HPointN *pt2);
49 static inline HPointN *Pt4ToHPtN(const HPoint3 *v4, HPointN *vN);
50 static inline HPointN *Pt3ToHPtN(Point3 *v3, HPointN *vN);
51 static inline HPointN *HPt3ToHPtN(const HPoint3 *v4, int *perm, HPointN *vN);
52 static inline HPoint3 *HPtNToHPt3(const HPointN *from, int *axes, HPoint3 *hp3);
53 static inline HPointN *HPtNPad(const HPointN *pt1, int dim2, HPointN *pt2);
54 static inline HPointN *HPtNAdd(const HPointN *pt1, const HPointN *pt2,
55 HPointN *sum);
56 static inline int HPtNSpace(const HPointN *pt);
57 static inline HPointN *HPtNSetSpace(HPointN *pt, int space);
58 static inline HPointN *HPtNComb(HPtNCoord s1, const HPointN *pt1,
59 HPtNCoord s2, const HPointN *pt2,
60 HPointN *sum);
61 static inline HPtNCoord HPtNDot(const HPointN *p1, const HPointN *p2);
62 static inline HPtNCoord HPtNDehomogenize(const HPointN *from, HPointN *to);
63 static inline HPointN *HPtNTransform3(Transform3 T, int *perm,
64 const HPointN *from, HPointN *to);
65 static inline HPointN *HPtNTransform(const TransformN *T,
66 const HPointN *from, HPointN *to);
67 static inline HPointN *HPtNTTransform(const TransformN *T,
68 const HPointN *from, HPointN *to);
69 static inline HPoint3 *HPtNTransProj(const TransformN *T,
70 const HPointN *from,
71 HPoint3 *result);
72 static inline HPoint3 *HPtNTransformComponents(const TransformN *T,
73 const HPointN *from,
74 int *perm,
75 HPoint3 *results);
76 static inline HPointN *HPt3NTransform(const TransformN *T,
77 const HPoint3 *from, HPointN *to);
78 static inline HPointN *Pt4NTransform(const TransformN *T,
79 const HPoint3 *from, HPointN *to);
80 static inline HPt3Coord HPtNNTransPt3(TransformN *TN, int *axes,
81 const HPointN *ptN, Point3 *result);
82 static inline HPt3Coord HPt3NTransPt3(TransformN *TN, int *axes,
83 const HPoint3 *hpt4, int v4d,
84 Point3 *result);
85 static inline void HPt3NTransHPt3(TransformN *TN, int *axes,
86 const HPoint3 *hpt4, int v4d,
87 HPoint3 *result);
88 static inline void HPtNMinMax(HPointN *min, HPointN *max, HPointN *other);
89
90 #include "transform3.h"
91
92 #ifndef max
93 # define _gv_hptn_max_
94 # define max(a,b) ((a) > (b) ? (a) : (b))
95 #endif
96 #ifndef min
97 # define _gv_hptn_min_
98 # define min(a,b) ((a) < (b) ? (a) : (b))
99 #endif
100
HPtNFreeListPrune(void)101 static inline void HPtNFreeListPrune(void)
102 {
103 FreeListNode *old;
104 HPointN *oldpt;
105 size_t size = 0;
106
107 while (HPointNFreeList) {
108 old = HPointNFreeList;
109 HPointNFreeList = old->next;
110 oldpt = (HPointN *)old;
111 if (oldpt->size && oldpt->v) {
112 OOGLFree(oldpt->v);
113 size += oldpt->size * sizeof(HPtNCoord);
114 }
115 OOGLFree(old);
116 size += sizeof(HPointN);
117 }
118 OOGLWarn("Freed %ld bytes.\n", size);
119 }
120
121 static inline HPointN *
HPtNCreate(int dim,const HPtNCoord * vec)122 HPtNCreate(int dim, const HPtNCoord *vec)
123 {
124 HPointN *pt;
125
126 FREELIST_NEW(HPointN, pt);
127
128 if(dim <= 0) {
129 dim = 1;
130 }
131 pt->dim = dim;
132 pt->flags = 0; /* for now */
133 if (pt->size < dim) {
134 pt->v = OOGLRenewNE(HPtNCoord, pt->v, dim, "new HPointN data");
135 pt->size = dim;
136 }
137 if(vec == NULL) {
138 memset(pt->v+1, 0, (dim-1)*sizeof(HPtNCoord));
139 pt->v[0] = 1.0;
140 } else {
141 memcpy(pt->v, vec, dim*sizeof(HPtNCoord));
142 }
143 return pt;
144 }
145
146 static inline void
HPtNDelete(HPointN * pt)147 HPtNDelete(HPointN *pt)
148 {
149 if(pt) {
150 /* if(pt->v) OOGLFree(pt->v); */
151 FREELIST_FREE(HPointN, pt);
152 }
153 }
154
155 static inline HPointN *
HPtNCopy(const HPointN * pt1,HPointN * pt2)156 HPtNCopy(const HPointN *pt1, HPointN *pt2)
157 {
158 if(pt2 == NULL) {
159 pt2 = HPtNCreate(pt1->dim, pt1->v);
160 } else {
161 if(pt2->dim != pt1->dim) {
162 pt2->v = OOGLRenewNE(HPtNCoord, pt2->v, pt1->dim, "renew HPointN");
163 pt2->dim = pt1->dim;
164 }
165 memcpy(pt2->v, pt1->v, pt1->dim*sizeof(HPtNCoord));
166 }
167 return pt2;
168 }
169
170 /* Convert a HPoint3 into a HPointN while interpreting the HPoint3 as
171 * a 4-point. This means that we do NOT perfom dehomogenization here.
172 */
173 static inline HPointN *
Pt4ToHPtN(const HPoint3 * v4,HPointN * vN)174 Pt4ToHPtN(const HPoint3 *v4, HPointN *vN)
175 {
176 int i;
177
178 if (!vN) {
179 vN = HPtNCreate(5, NULL);
180 } else if (vN->dim < 5) {
181 vN->v = OOGLRenewNE(HPtNCoord, vN->v, 5, "renew HPointN");
182 vN->dim = 5;
183 }
184 vN->v[0] = 1.0;
185 for (i = 0; i < 4; ++i) {
186 vN->v[i+1] = ((HPt3Coord *)v4)[i];
187 }
188 for (++i; i < vN->dim; i++) {
189 vN->v[i] = 0.0;
190 }
191 return vN;
192 }
193
194 /* Convert a Point3 into a HPointN */
195 static inline HPointN *
Pt3ToHPtN(Point3 * v3,HPointN * vN)196 Pt3ToHPtN(Point3 *v3, HPointN *vN)
197 {
198 int i;
199
200 if (!vN) {
201 vN = HPtNCreate(4, NULL);
202 } else if (vN->dim < 4) {
203 vN->v = OOGLRenewNE(HPtNCoord, vN->v, 4, "renew HPointN");
204 vN->dim = 5;
205 }
206 vN->v[0] = 1.0;
207 for (i = 0; i < 3; ++i) {
208 vN->v[i+1] = ((HPt3Coord *)v3)[i];
209 }
210 for (++i; i < vN->dim; i++) {
211 vN->v[i] = 0.0;
212 }
213
214 return vN;
215 }
216
217 /* Convert a HPoint3 into a HPointN while interpreting the HPoint3 as
218 * a 3-point. The homogeneous component of v4 is moved to vN->v[0].
219 */
220 static inline HPointN *
HPt3ToHPtN(const HPoint3 * v4,int * perm,HPointN * vN)221 HPt3ToHPtN(const HPoint3 *v4, int *perm, HPointN *vN)
222 {
223 const int d3 = 4;
224 int i;
225 int perm_dim;
226
227 if (!perm) {
228 perm_dim = 4;
229 } else {
230 perm_dim = perm[0];
231 for (i = 0; i < d3; i++) {
232 perm_dim = max(perm[i], perm_dim);
233 }
234 ++perm_dim;
235 }
236 if (!vN) {
237 vN = HPtNCreate(perm_dim, NULL);
238 } else if (vN->dim < perm_dim) {
239 vN->v = OOGLRenewNE(HPtNCoord, vN->v, perm_dim, "renew HPointN");
240 vN->dim = perm_dim;
241 }
242 if (!perm) {
243 vN->v[0] = v4->w;
244 vN->v[1] = v4->x;
245 vN->v[2] = v4->y;
246 vN->v[3] = v4->z;
247 for (i = 4; i < vN->dim; i++) {
248 vN->v[i] = 0.0;
249 }
250 } else {
251 memset(vN->v, 0, vN->dim*sizeof(HPtNCoord));
252 for (i = 0; i < d3; i++) {
253 vN->v[perm[i]] = ((HPt3Coord *)v4)[i];
254 }
255 }
256 return vN;
257 }
258
259 /* Transform a 4-point to a 3-point according to the mapping defined
260 * in "axes".
261 */
262 static inline HPoint3 *
HPtNToHPt3(const HPointN * from,int * axes,HPoint3 * hp3)263 HPtNToHPt3(const HPointN *from, int *axes, HPoint3 *hp3)
264 {
265 HPt3Coord *to = (HPt3Coord *)hp3;
266 int i, dim = from->dim;
267
268 if (!axes) {
269 hp3->w = from->v[0];
270 hp3->x = from->v[1];
271 hp3->y = from->v[2];
272 hp3->z = from->v[3];
273 } else {
274 for (i = 0; i < 4; i++) {
275 if (axes[i] > dim-1) {
276 to[i] = 0.0;
277 } else {
278 to[i] = from->v[axes[i]];
279 }
280 }
281 }
282 return hp3;
283 }
284
285 /* Copy pt1 to pt2 and make sure that pt2 has dimension dim2. If
286 * pt1->dim < dim2, then pt1 is implicitly padded with zeros.
287 */
288 static inline HPointN *
HPtNPad(const HPointN * pt1,int dim2,HPointN * pt2)289 HPtNPad(const HPointN *pt1, int dim2, HPointN *pt2)
290 {
291 int dim1 = pt1->dim;
292
293 if( dim2 < 1 ) { /* Uhh? */
294 return(NULL);
295 }
296
297 if(pt1 == NULL) {
298 pt2 = HPtNCreate(dim2, NULL);
299 return pt2;
300 }
301
302 if(pt1 != pt2) {
303 if(pt2 == NULL) {
304 pt2 = HPtNCreate(dim2,NULL);
305 } else if (pt2->dim != dim2) {
306 pt2->v = OOGLRenewNE(HPtNCoord, pt2->v, dim2, "renew HPointN");
307 pt2->dim = dim2;
308 }
309 if (dim1 <= dim2) {
310 memcpy(pt2->v, pt1->v, dim1* sizeof(HPtNCoord));
311 memset(pt2->v+dim1, 0, (dim2-dim1)*sizeof(HPtNCoord));
312 } else {
313 memcpy(pt2->v, pt1->v, dim2*sizeof(HPtNCoord));
314 }
315 } else {
316 /* now that the homogeneous component is at 0 this is a simple
317 * padding operation.
318 */
319 if (pt2->dim != dim2) {
320 pt2->v = OOGLRenewNE(HPtNCoord, pt2->v, dim2, "renew HPointN");
321 }
322 if (dim2 > pt2->dim) {
323 memset(pt2->v+pt2->dim, 0, (dim2-pt2->dim)*sizeof(HPtNCoord));
324 }
325 }
326 return pt2;
327 }
328
329 /* Add two HPointN's */
330 static inline HPointN *
HPtNAdd(const HPointN * pt1,const HPointN * pt2,HPointN * sum)331 HPtNAdd(const HPointN *pt1, const HPointN *pt2, HPointN *sum)
332 {
333 int dim1 = pt1->dim, dim2 = pt2->dim;
334 HPtNCoord c1 = pt1->v[0], c2 = pt2->v[0];
335 int i;
336
337 if (dim1 == dim2) {
338 if (sum == NULL) {
339 sum = HPtNCreate(dim1, NULL);
340 } else if(sum->dim < dim1) {
341 sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
342 sum->dim = dim1;
343 }
344 sum->v[0] = c1*c2;
345 for (i = 1; i < dim1; i++) {
346 sum->v[i] = c2*pt1->v[i] + c1*pt2->v[i];
347 }
348 } else {
349 /* make sure pt1 is the larger one */
350 if (dim1 < dim2) {
351 const HPointN *swap = pt1;
352
353 pt1 = pt2; pt2 = swap; dim2 = dim1; dim1 = pt1->dim;
354 }
355
356 if (sum == NULL) {
357 sum = HPtNCreate(dim1,NULL);
358 } else if(sum->dim < dim1) {
359 sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
360 sum->dim = dim1;
361 }
362 sum->v[0] = c1*c2;
363 for (i = 1; i < dim2; i++) {
364 sum->v[i] = c2*pt1->v[i] + c1*pt2->v[i];
365 }
366 for (; i < dim1; i++) {
367 sum->v[i] = c2*pt1->v[i];
368 }
369 }
370 return sum;
371 }
372
373 /* Space */
374 static inline int
HPtNSpace(const HPointN * pt)375 HPtNSpace(const HPointN *pt)
376 {
377 (void)pt;
378 return TM_EUCLIDEAN;
379 }
380
381 static inline HPointN *
HPtNSetSpace(HPointN * pt,int space)382 HPtNSetSpace(HPointN *pt, int space)
383 {
384 if (space != TM_EUCLIDEAN) {
385 OOGLError(1, "Non-Euclidean space not support in higher dimensions.\n");
386 return NULL;
387 }
388 return pt;
389 }
390
391 /* Linear combination */
392 static inline HPointN *
HPtNComb(HPtNCoord s1,const HPointN * pt1,HPtNCoord s2,const HPointN * pt2,HPointN * sum)393 HPtNComb(HPtNCoord s1, const HPointN *pt1,
394 HPtNCoord s2, const HPointN *pt2,
395 HPointN *sum)
396 {
397 int dim1 = pt1->dim, dim2 = pt2->dim;
398 HPtNCoord c1 = pt1->v[0], c2 = pt2->v[0];
399 int i;
400
401 if (dim1 == dim2) {
402 if (sum == NULL) {
403 sum = HPtNCreate(dim1, NULL);
404 } else if(sum->dim < dim1) {
405 sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
406 sum->dim = dim1;
407 }
408 sum->v[0] = c1*c2;
409 for (i = 1; i < dim1; i++) {
410 sum->v[i] = c2*s1*pt1->v[i] + c1*s2*pt2->v[i];
411 }
412 } else {
413 /* make sure pt1 is the larger one */
414 if (dim1 < dim2) {
415 const HPointN *swap = pt1;
416 pt1 = pt2; pt2 = swap; dim2 = dim1; dim1 = pt1->dim;
417 }
418
419 if (sum == NULL) {
420 sum = HPtNCreate(dim1,NULL);
421 } else if(sum->dim < dim1) {
422 sum->v = OOGLRenewNE(HPtNCoord, sum->v, dim1, "renew HPointN");
423 sum->dim = dim1;
424 }
425 sum->v[0] = c1*c2;
426 for (i = 1; i < dim2; i++) {
427 sum->v[i] = c2*s1*pt1->v[i] + c1*s2*pt2->v[i];
428 }
429 for (; i < dim1; i++) {
430 sum->v[i] = c2*s1*pt1->v[i];
431 }
432 }
433 return sum;
434 }
435
436 /* Dot product of two vectors */
437 static inline HPtNCoord
HPtNDot(const HPointN * p1,const HPointN * p2)438 HPtNDot(const HPointN *p1, const HPointN *p2)
439 {
440 HPtNCoord result;
441 int i;
442 int dim = p1->dim;
443
444 if (p2->dim < dim) {
445 dim = p2->dim;
446 }
447
448 result = 0;
449 for(i = 1; i< dim; i++)
450 result += p1->v[i] * p2->v[i];
451
452 return result / (p1->v[0] * p2->v[0]);
453 }
454
455 /* Dehomogenize */
456 static inline HPtNCoord
HPtNDehomogenize(const HPointN * from,HPointN * to)457 HPtNDehomogenize(const HPointN *from, HPointN *to)
458 {
459 int dim = from->dim;
460 HPtNCoord c = from->v[0], inv = 1.0 / c;
461 int i;
462
463 if (c == 1.0 || c == 0.0) {
464 if (from != to) {
465 HPtNCopy(from, to);
466 }
467 return (HPtNCoord)0.0;
468 }
469
470 if (to == NULL) {
471 to = HPtNCreate(dim, NULL);
472 } else if (to->dim != dim) {
473 to->v = OOGLRenewNE(HPtNCoord, to->v, dim, "renew HPointN");
474 to->dim = dim;
475 }
476
477 for( i=1; i < dim; i++)
478 to->v[i] = from->v[i] * inv;
479 to->v[0] = 1.0;
480
481 return c;
482 }
483
484 /* Transform an HPointN according to a 3d transform acting only on the
485 * sub-space defined by "axes". The standard axes should look like
486 * {dx,dy,dz,0} because the homogeneous divisor of HPoint3's is
487 * located at index 3.
488 */
489 static inline HPointN *
HPtNTransform3(Transform3 T,int * perm,const HPointN * from,HPointN * to)490 HPtNTransform3(Transform3 T, int *perm, const HPointN *from, HPointN *to)
491 {
492 const int d3 = 4;
493 int i;
494 HPoint3 from3;
495 HPt3Coord *from3data = HPoint3Data(&from3);
496 int perm_dim;
497
498 if (!perm) {
499 perm_dim = 4;
500 } else {
501 perm_dim = perm[0];
502 for (i = 0; i < d3; i++) {
503 perm_dim = max(perm[i], perm_dim);
504 }
505 ++perm_dim;
506 }
507 HPtNToHPt3(from, perm, &from3);
508
509 HPt3Transform(T, &from3, &from3);
510 if (from->dim < perm_dim) {
511 to = HPtNPad(from, perm_dim, to);
512 } else {
513 to = HPtNCopy(from, to);
514 }
515 if (perm) {
516 for (i = 0; i < 4; i++) {
517 to->v[perm[i]] = from3data[i];
518 }
519 } else {
520 to->v[0] = from3data[3];
521 to->v[1] = from3data[0];
522 to->v[2] = from3data[1];
523 to->v[3] = from3data[2];
524 }
525 return to;
526 }
527
528 /* Apply a TransformN to an HPointN.
529 *
530 * If from->dim < T->idim, then from is implicitly padded with zeros,
531 * if from->dim > T->idim, then T is interpreted as identity on the
532 * sub-space defined by T->idim < idx < T->odim, the remaining input
533 * components are mapped to zero.
534 */
535 static inline HPointN *
HPtNTransform(const TransformN * T,const HPointN * from,HPointN * to)536 HPtNTransform(const TransformN *T, const HPointN *from, HPointN *to)
537 {
538 int idim, odim, dim = from->dim;
539 int i, j;
540 VARARRAY(vspace, HPtNCoord, dim);
541 HPtNCoord *v;
542
543 if (!T) {
544 return HPtNCopy(from, to);
545 }
546
547 idim = T->idim;
548 odim = T->odim;
549
550 if (from == to) {
551 v = vspace;
552 for (i = 0; i < dim; i++) {
553 v[i] = from->v[i];
554 }
555 } else {
556 v = from->v;
557 }
558
559 if(to == NULL) {
560 to = HPtNCreate(odim, NULL);
561 } else if (to->dim != odim) {
562 to->v = OOGLRenewNE(HPtNCoord, to->v, odim, "renew HPointN");
563 to->dim = odim;
564 }
565
566 if (idim == dim) {
567 /* the easy case */
568 for (i = 0; i < odim; i++) {
569 to->v[i] = 0;
570 for (j = 0; j < idim; j++) {
571 to->v[i] += v[j] * T->a[j*odim+i];
572 }
573 }
574 } else if (idim > dim) {
575 /* pad with zeroes, the homogeneous component sits at index zero
576 and is automatically handled correctly. */
577 for(i = 0; i < odim; i++) {
578 to->v[i] = 0;
579 for (j = 0; j < dim; j++) {
580 to->v[i] += v[j] * T->a[j*odim+i];
581 }
582 }
583 } else { /* obviously the case idim < dim */
584 /* implicitly pad the matrix, i.e. T acts as unity on the missing
585 * dimens+ions.
586 */
587 for (i = 0; i < odim; i++) {
588 to->v[i] = 0;
589 for (j = 0; j < idim; j++) {
590 to->v[i] += v[j] * T->a[j*odim+i];
591 }
592 if (i >= idim && i < dim) {
593 to->v[i] += v[i];
594 }
595 }
596 }
597
598 return to;
599 }
600
601 /* Apply the transpose of a TransformN to an HPointN.
602 *
603 * If from->dim < T->odim, then form is implicitly padded with zeros,
604 * if from->dim > T->odim, then T is interpreted as identity on the
605 * sub-space defined by T->odim < idx < T->idim, the remaining input
606 * components are mapped to zero.
607 */
608 static inline HPointN *
HPtNTTransform(const TransformN * T,const HPointN * from,HPointN * to)609 HPtNTTransform(const TransformN *T, const HPointN *from, HPointN *to)
610 {
611 int idim, odim, dim = from->dim;
612 int i, j;
613 VARARRAY(vspace, HPtNCoord, dim);
614 HPtNCoord *v;
615
616 if (!T) {
617 return HPtNCopy(from, to);
618 }
619
620 idim = T->idim;
621 odim = T->odim;
622
623 if (from == to) {
624 v = vspace;
625 for (i = 0; i < dim; i++) {
626 v[i] = from->v[i];
627 }
628 } else {
629 v = from->v;
630 }
631
632 if(to == NULL) {
633 to = HPtNCreate(odim, NULL);
634 } else if (to->dim != idim) {
635 to->v = OOGLRenewNE(HPtNCoord, to->v, idim, "renew HPointN");
636 to->dim = idim;
637 }
638
639 if (odim == dim) {
640 /* the easy case */
641 for (i = 0; i < idim; i++) {
642 to->v[i] = 0;
643 for (j = 0; j < odim; j++) {
644 to->v[i] += v[j] * T->a[i*odim+j];
645 }
646 }
647 } else if (odim > dim) {
648 /* pad with zeroes, the homogeneous component sits at index zero
649 and is automatically handled correctly. */
650 for(i = 0; i < idim; i++) {
651 to->v[i] = 0;
652 for (j = 0; j < dim; j++) {
653 to->v[i] += v[j] * T->a[i*odim+j];
654 }
655 }
656 } else { /* obviously the case odim < dim */
657 /* implicitly pad the matrix, i.e. T acts as unity on the missing
658 * dimens+ions.
659 */
660 for (i = 0; i < idim; i++) {
661 to->v[i] = 0;
662 for (j = 0; j < odim; j++) {
663 to->v[i] += v[j] * T->a[i*odim+j];
664 }
665 if (i >= odim && i < dim) {
666 to->v[i] += v[i];
667 }
668 }
669 }
670
671 return to;
672 }
673
674 /* A short and efficient routine. This routine assumes that T is the
675 * result of a call to TmNProject() and is therefor non-NULL. "from"
676 * may have more or less dimensions than T->idim, it is padded with
677 * zeroes or truncated as necessary. from must have dimension 1 at
678 * least.
679 *
680 * T already contains all necessary permutations and the projection to
681 * map "from" to result "properly".
682 */
683 static inline HPoint3 *
HPtNTransProj(const TransformN * T,const HPointN * from,HPoint3 * result)684 HPtNTransProj(const TransformN *T,
685 const HPointN *from,
686 HPoint3 *result)
687 {
688 HPtNCoord *v = from->v;
689 int idim;
690 const int odim = 4; /* must be 4 */
691 int i;
692
693 idim = T->idim > from->dim ? from->dim : T->idim;
694
695 result->x = v[0] * T->a[0*odim+0];
696 result->y = v[0] * T->a[0*odim+1];
697 result->z = v[0] * T->a[0*odim+2];
698 result->w = v[0] * T->a[0*odim+3];
699 for (i = 1; i < idim; i++) {
700 result->x += v[i] * T->a[i*odim+0];
701 result->y += v[i] * T->a[i*odim+1];
702 result->z += v[i] * T->a[i*odim+2];
703 result->w += v[i] * T->a[i*odim+3];
704 }
705 return result;
706 }
707
708 /* Transform p by T and then project to the sub-space defined by
709 * "perm".
710 */
711 static inline HPoint3 *
HPtNTransformComponents(const TransformN * T,const HPointN * from,int * perm,HPoint3 * results)712 HPtNTransformComponents(const TransformN *T,
713 const HPointN *from,
714 int *perm,
715 HPoint3 *results)
716 {
717 int idim = T->idim, odim = T->odim, dim = from->dim;
718 int i, j, k;
719 HPt3Coord *res = (HPt3Coord *)results;
720
721 if (idim == dim) {
722 /* the easy case */
723 for (k = 0; k < 4; k++) {
724 i = perm[k];
725 if (i > odim)
726 continue;
727 res[k] = 0;
728 for (j = 0; j < idim; j++) {
729 res[k] += from->v[j] * T->a[j*odim+i];
730 }
731 }
732 } else if (idim > dim) {
733 /* pad with zeroes, the homogeneous component sits at index zero
734 and is automatically handled correctly. */
735 for (k = 0; k < 4; k++) {
736 i = perm[k];
737 if (i > odim)
738 continue;
739 res[k] = 0;
740 for (j = 0; j < dim; j++) {
741 res[k] += from->v[j] * T->a[j*odim+i];
742 }
743 }
744 } else { /* obviously the case idim < dim */
745 /* implicitly pad the matrix, i.e. T acts as unity on the missing
746 * dimens+ions.
747 */
748 for (k = 0; k < 4; k++) {
749 i = perm[k];
750 if (i > odim)
751 continue;
752 res[k] = 0;
753 for (j = 0; j < idim; j++) {
754 res[k] += from->v[j] * T->a[j*odim+i];
755 }
756 if (i >= idim && i < dim) {
757 res[k] += from->v[i];
758 }
759 }
760 }
761
762 return results;
763 }
764
765 static inline HPointN *
HPt3NTransform(const TransformN * T,const HPoint3 * from,HPointN * to)766 HPt3NTransform(const TransformN *T, const HPoint3 *from, HPointN *to)
767 {
768 int idim, odim;
769 int i, j;
770 HPtNCoord *v;
771
772 if (!T) {
773 return HPt3ToHPtN(from, NULL, to);
774 }
775
776 v = (HPtNCoord *)from;
777
778 idim = T->idim;
779 odim = T->odim;
780
781 if(to == NULL) {
782 to = HPtNCreate(odim, NULL);
783 } else if (to->dim != odim) {
784 to->v = OOGLRenewNE(HPtNCoord, to->v, odim, "renew HPointN");
785 to->dim = odim;
786 }
787
788 if (idim == 4) {
789 /* the easy case */
790 for (i = 0; i < odim; i++) {
791 to->v[i] = 0;
792 for (j = 0; j < idim; j++) {
793 to->v[i] += v[(j+3)%4] * T->a[j*odim+i];
794 }
795 }
796 } else if (idim > 4) {
797 /* pad with zeroes, the homogeneous component sits at index zero
798 and is automatically handled correctly. */
799 for(i = 0; i < odim; i++) {
800 to->v[i] = 0;
801 for (j = 0; j < 4; j++) {
802 to->v[i] += v[(j+3)%4] * T->a[j*odim+i];
803 }
804 }
805 } else { /* obviously the case idim < dim */
806 /* implicitly pad the matrix, i.e. T acts as unity on the missing
807 * dimens+ions.
808 */
809 for (i = 0; i < odim; i++) {
810 to->v[i] = 0;
811 for (j = 0; j < idim; j++) {
812 to->v[i] += v[(j+3) % 4] * T->a[j*odim+i];
813 }
814 if (i >= idim && i < 4) {
815 to->v[i] += v[i];
816 }
817 }
818 }
819
820 return to;
821 }
822
823 static inline HPointN *
Pt4NTransform(const TransformN * T,const HPoint3 * from,HPointN * to)824 Pt4NTransform(const TransformN *T, const HPoint3 *from, HPointN *to)
825 {
826 int idim, odim;
827 int i, j;
828 HPtNCoord *v;
829
830 if (!T) {
831 return Pt4ToHPtN(from, to);
832 }
833
834 v = (HPtNCoord *)from;
835
836 idim = T->idim;
837 odim = T->odim;
838
839 if(to == NULL) {
840 to = HPtNCreate(odim, NULL);
841 } else if (to->dim != odim) {
842 to->v = OOGLRenewNE(HPtNCoord, to->v, odim, "renew HPointN");
843 to->dim = odim;
844 }
845
846 if (idim == 5) {
847 /* the easy case */
848 for (i = 0; i < odim; i++) {
849 to->v[i] = T->a[i];
850 for (j = 1; j < idim; j++) {
851 to->v[i] += v[j-1] * T->a[j*odim+i];
852 }
853 }
854 } else if (idim > 5) {
855 /* pad with zeroes, the homogeneous component sits at index zero
856 and is automatically handled correctly. */
857 for(i = 0; i < odim; i++) {
858 to->v[i] = T->a[i];
859 for (j = 1; j < 5; j++) {
860 to->v[i] += v[j-1] * T->a[j*odim+i];
861 }
862 }
863 } else { /* obviously the case idim < dim */
864 /* implicitly pad the matrix, i.e. T acts as unity on the missing
865 * dimens+ions.
866 */
867 for (i = 0; i < odim; i++) {
868 to->v[i] = T->a[i];
869 for (j = 0; j < idim; j++) {
870 to->v[i] += v[j-1] * T->a[j*odim+i];
871 }
872 if (i >= idim && i < 5) {
873 to->v[i] += v[i-1];
874 }
875 }
876 }
877
878 return to;
879 }
880
881 static inline HPt3Coord
HPtNNTransPt3(TransformN * TN,int * axes,const HPointN * ptN,Point3 * result)882 HPtNNTransPt3(TransformN *TN, int *axes, const HPointN *ptN, Point3 *result)
883 {
884 HPoint3 tmp;
885
886 HPtNTransformComponents(TN, ptN, axes, &tmp);
887
888 result->x = tmp.x / tmp.w;
889 result->y = tmp.y / tmp.w;
890 result->z = tmp.z / tmp.w;
891
892 return tmp.w;
893 }
894
895 static inline HPt3Coord
HPt3NTransPt3(TransformN * TN,int * axes,const HPoint3 * hpt4,int v4d,Point3 * result)896 HPt3NTransPt3(TransformN *TN, int *axes,
897 const HPoint3 *hpt4, int v4d, Point3 *result)
898 {
899 HPointN *tmp;
900 HPt3Coord retval;
901
902 /* axes[3] should be 0 */
903
904 if (v4d) {
905 tmp = Pt4NTransform(TN, hpt4, NULL);
906 } else {
907 tmp = HPt3NTransform(TN, hpt4, NULL);
908 }
909 result->x = tmp->v[axes[0]] / tmp->v[axes[3]];
910 result->y = tmp->v[axes[1]] / tmp->v[axes[3]];
911 result->z = tmp->v[axes[2]] / tmp->v[axes[3]];
912
913 retval = tmp->v[axes[3]];
914
915 HPtNDelete(tmp);
916
917 return retval;
918 }
919
920 static inline void
HPt3NTransHPt3(TransformN * TN,int * axes,const HPoint3 * hpt4,int v4d,HPoint3 * result)921 HPt3NTransHPt3(TransformN *TN, int *axes,
922 const HPoint3 *hpt4, int v4d, HPoint3 *result)
923 {
924 HPointN *tmp;
925
926 if (v4d) {
927 tmp = Pt4NTransform(TN, hpt4, NULL);
928 } else {
929 tmp = HPt3NTransform(TN, hpt4, NULL);
930 }
931
932 result->w = tmp->v[axes[3]];
933 result->x = tmp->v[axes[0]];
934 result->y = tmp->v[axes[1]];
935 result->z = tmp->v[axes[2]];
936
937 HPtNDelete(tmp);
938 }
939
940 /* Utility function for bounding box computations. We assume that min
941 * and max are dehomogenized (third part of the #if clause below), and
942 * that they are large enough (min->dim >= other->dim <= max->dim)
943 */
HPtNMinMax(HPointN * min,HPointN * max,HPointN * other)944 static inline void HPtNMinMax(HPointN *min, HPointN *max, HPointN *other)
945 {
946 #if 0
947 int i;
948
949 for (i = 1; i < other->dim; i++) {
950 if (min->v[i] > other->v[i]) {
951 min->v[i] = other->v[i];
952 } else if (max->v[i] < other->v[i]) {
953 max->v[i] = other->v[i];
954 }
955 }
956 #elif 0
957 int i;
958 HPtNCoord c = min->v[0], C = max->v[0], oc = other->v[0];
959
960 for (i = 1; i < other->dim; i++) {
961 if (oc * min->v[i] > c * other->v[i]) {
962 min->v[i] = other->v[i] * c / oc ;
963 } else if (oc * max->v[i] < C * other->v[i]) {
964 max->v[i] = other->v[i] * C / oc;
965 }
966 }
967 #else
968 int i;
969 HPtNCoord oc = other->v[0];
970
971 for (i = 1; i < other->dim; i++) {
972 if (oc * min->v[i] > other->v[i]) {
973 min->v[i] = other->v[i] / oc ;
974 } else if (oc * max->v[i] < other->v[i]) {
975 max->v[i] = other->v[i] / oc;
976 }
977 }
978 #endif
979 }
980
981 #ifdef _gv_hptn_max_
982 # undef max
983 #endif
984 #ifdef _gv_hptn_min_
985 # undef min
986 #endif
987
988 #endif /* _GV_HPOINTN_H_ */
989
990 /*
991 * Local Variables: ***
992 * c-basic-offset: 2 ***
993 * End: ***
994 */
995