1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 * Copyright (C) 2006-20007 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 #if HAVE_CONFIG_H
24 # include "config.h"
25 #endif
26
27 #if 0
28 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
29 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
30 #endif
31
32 #include "mgP.h"
33 #include "hpointn.h"
34 #include "../common/drawer.h"
35
36 #include <stdlib.h>
37 #ifndef alloca
38 #include <alloca.h>
39 #endif
40
41 #if 0
42 static void compute_mg_factor(NDstuff *nds);
43 static inline void fast_map_ND_point(NDstuff *nds, HPointN *from, HPoint3 *to);
44 #endif
45
46 /*
47 * Arguments:
48 * mginfo - pointer to an "ndstuff" structure (from drawer.h), containing:
49 * T - transformation from object to N-D camera space
50 * axes - array of 4 axis indices, giving the subspace seen by this view
51 * Tc - (dimension) x (n-color-axes) transform from object space to
52 * color axis space. Use NULL if not projecting colors.
53 * Given vectors in camera space, could get Tc with:
54 * Tc = TmNConcat(T, {matrix of projection vectors}, NULL)
55 * ncm - number of colormaps
56 * cm - colormap array
57 * p - the N-D HPointN to be transformed, in object space
58 * np - resulting 3-D HPoint3 point, in camera space
59 * c - resulting color; only changed if ncm > 0
60 *
61 * Results: returns 1 if color was set, else 0.
62 */
63
64 static int
map_ND_point(mgNDctx * mgNDctx,HPointN * p,HPoint3 * np,ColorA * c)65 map_ND_point(mgNDctx *mgNDctx, HPointN *p, HPoint3 *np, ColorA *c)
66 {
67 NDstuff *nds = (NDstuff *)mgNDctx;
68 #define iT nds->T
69 #define iTc nds->Tc
70 #define incm nds->ncm
71 #define icm nds->cm
72 #define ihc nds->hc
73 ColorA ci;
74 float t;
75 int i;
76
77 #if 0
78 /* cH: here is still room for optimization: provided that iT has
79 * maximal rank (i.e. rank 4) it would be possible to factor out an
80 * ordinary Transform and handle that one over to the MG layer. The
81 * remaining N-Transform would act as identity on a 4x4 sub-space
82 * which would reduce the number of floating point operations to
83 * perform the projection. The MG layer uses an object-camera
84 * transform anyway, so this would speed up things.
85 *
86 * Well. The code at the end of this files hacks this, but: the
87 * resulting 4x4 factor is -- of course -- not orthogonal. Therefore
88 * I stopped working at this optimization. The code for the
89 * factorization is at the end of this file. cH.
90 *
91 */
92 fast_map_ND_point(nds, p, np);
93 #else
94 HPtNTransProj(iT, p, np);
95 HPt3Dehomogenize(np, np);
96 #endif
97
98 if(incm > 0 && !(_mgc->astk->ap.flag & APF_KEEPCOLOR) && c != NULL) {
99 ci.r = ci.g = ci.b = ci.a = 0;
100
101 HPtNTransform(iTc, p, ihc);
102 for(i = 0; i < incm; i++) {
103 cent *ce = VVEC(icm[i].cents, cent);
104 float v = ihc->v[i];
105 if(!(v > ce->v)) {
106 ci.r += ce->c.r;
107 ci.g += ce->c.g;
108 ci.b += ce->c.b;
109 ci.a += ce->c.a;
110 } else {
111 do ce++; while(v > ce->v);
112 if((ce-1)->interp) {
113 t = (v - (ce-1)->v) / (ce->v - (ce-1)->v);
114 } else {
115 t = 1;
116 }
117 ci.r += t*ce->c.r + (1-t)*(ce-1)->c.r;
118 ci.g += t*ce->c.g + (1-t)*(ce-1)->c.g;
119 ci.b += t*ce->c.b + (1-t)*(ce-1)->c.b;
120 ci.a += t*ce->c.a + (1-t)*(ce-1)->c.a;
121 }
122 }
123 /* XXX Clamp colors to range 0..1 here?? */
124 if(ci.r < 0) ci.r = 0; else if(ci.r > 1) ci.r = 1;
125 if(ci.g < 0) ci.g = 0; else if(ci.g > 1) ci.g = 1;
126 if(ci.b < 0) ci.b = 0; else if(ci.b > 1) ci.b = 1;
127 if(ci.a < 0) ci.a = 0; else if(ci.a > 1) ci.a = 1;
128 *c = ci;
129 return 1;
130 }
131 return 0;
132
133 #undef iT
134 #undef iTc
135 #undef iaxes
136 #undef icm
137 #undef incm
138 #undef ihc
139 }
140
141 /* push/pop could be done much more efficiently, but ND is inefficient
142 * anyway, so what.
143 */
144 typedef struct savedCTX
145 {
146 TransformN *T;
147 TransformN *Tc;
148 TransformN *rest;
149 Transform MGFactor;
150 int perm[1];
151 } SavedCTX;
152
saveCTX(mgNDctx * NDctx)153 static void *saveCTX(mgNDctx *NDctx)
154 {
155 NDstuff *nds = (NDstuff *)NDctx;
156 SavedCTX *savedCTX =
157 (SavedCTX *)OOGLNewNE(char,
158 sizeof(SavedCTX)+(nds->T->idim-1)*sizeof(int),
159 "SavedCTX");
160
161 savedCTX->T = nds->T;
162 nds->T = TmNCopy(nds->T, NULL);
163 if (nds->Tc) {
164 savedCTX->Tc = nds->Tc;
165 nds->Tc = TmNCopy(nds->Tc, NULL);
166 } else {
167 savedCTX->Tc = NULL;
168 }
169 #if 0
170 savedCTX->rest = nds->rest;
171 nds->rest = TmNCopy(nds->rest, NULL);
172 TmCopy(nds->MGFactor, savedCTX->MGFactor);
173 memcpy(savedCTX->perm, nds->perm, nds->T->idim * sizeof(int));
174
175 mgpushtransform();
176 #endif
177
178 return savedCTX;
179 }
180
pushTN(mgNDctx * NDctx,TransformN * TN)181 static void pushTN(mgNDctx *NDctx, TransformN *TN)
182 {
183 NDstuff *nds = (NDstuff *)NDctx;
184
185 TmNConcat(TN, nds->T, nds->T);
186
187 if (nds->Tc) {
188 TmNConcat(TN, nds->Tc, nds->Tc);
189 }
190
191 #if 0
192 compute_mg_factor(nds);
193 mgtransform(nds->MGFactor);
194 #endif
195 }
196
pushT(mgNDctx * NDctx,Transform T)197 static void pushT(mgNDctx *NDctx, Transform T)
198 {
199 NDstuff *nds = (NDstuff *)NDctx;
200
201 TmNApplyT3TN(T, NULL, nds->T);
202
203 if (nds->Tc) {
204 TmNApplyT3TN(T, NULL, nds->Tc);
205 }
206
207 #if 0
208 compute_mg_factor(nds);
209 mgtransform(nds->MGFactor);
210 #endif
211 }
212
restoreCTX(mgNDctx * NDctx,void * vsavedCTX)213 static void restoreCTX(mgNDctx *NDctx, void *vsavedCTX)
214 {
215 NDstuff *nds = (NDstuff *)NDctx;
216 SavedCTX *savedCTX = (SavedCTX *)vsavedCTX;
217
218 #if 0
219 mgpoptransform();
220 #endif
221
222 TmNDelete(nds->T);
223 nds->T = savedCTX->T;
224 if (nds->Tc) {
225 TmNDelete(nds->Tc);
226 }
227 nds->Tc = savedCTX->Tc;
228 #if 0
229 TmNDelete(nds->rest);
230 nds->rest = savedCTX->rest;
231 TmCopy(savedCTX->MGFactor, nds->MGFactor);
232 memcpy(nds->perm, savedCTX->perm, nds->T->idim * sizeof(int));
233 #endif
234
235 OOGLFree(savedCTX);
236 }
237
238 static mgNDctx NDctx_proto = {
239 map_ND_point,
240 saveCTX,
241 pushTN,
242 pushT,
243 restoreCTX,
244 NULL,
245 };
246
drawer_init_ndstuff(DView * dv,TransformN * W2C,TransformN * W2U)247 NDstuff *drawer_init_ndstuff(DView *dv, TransformN *W2C, TransformN *W2U)
248 {
249 NDstuff *nds = OOGLNewE(NDstuff, "new NDstuff");
250 int dim = W2C->idim;
251
252 #if 0
253 mgpushtransform();
254 #endif
255
256 nds->mgNDctx = NDctx_proto;
257
258 /* Initialize the color transform */
259 nds->ncm = dv->nNDcmap;
260 nds->cm = dv->NDcmap;
261 nds->W2c = NULL;
262 nds->Tc = NULL;
263 nds->hc = NULL;
264
265 if(dv->nNDcmap > 0) {
266 HPointN *caxis = NULL;
267 int i, j;
268
269 /* Build array of N-D-to-color projection vectors: it becomes a
270 * matrix, multiplied by N-D row vector on the left, yielding an
271 * array of dv->nNDcmap projections used for coloring. So it has
272 * (dimension) rows, (dv->nNDcmap) columns.
273 */
274 nds->W2c = TmNCreate(dim, dv->nNDcmap, NULL);
275 nds->hc = HPtNCreate(dim, NULL);
276
277 for(i = 0; i < dv->nNDcmap; i++) {
278 cmap *cm = &dv->NDcmap[i];
279 int cdim = cm->axis->dim;
280 int mindim = (dim < cdim) ? dim : cdim;
281 HPointN *our_caxis = cm->axis;
282 TransformN *TxC;
283
284 if(cm->coords != UNIVERSE) {
285 TxC = drawer_get_ND_transform(cm->coords, UNIVERSE);
286 if(TxC) {
287 our_caxis = caxis = HPtNTransform(TxC, cm->axis, caxis);
288 TmNDelete(TxC);
289 }
290 }
291 for(j = 0; j < mindim; j++)
292 nds->W2c->a[j*dv->nNDcmap + i] = our_caxis->v[j];
293 }
294 TmNConcat(W2U, nds->W2c, nds->W2c);
295 HPtNDelete(caxis);
296 }
297
298 /* initialize the W2C geometry transform. Note that nds->W2C is a
299 * projection matrix, i.e. it includes the mapping defined in
300 * dv->NDperm. This spares a lot of unnecessary floating point
301 * operations.
302 */
303 nds->T = NULL;
304
305 #if 0
306 nds->rest = NULL;
307 nds->perm = OOGLNewNE(int, dim, "permutation for fast ND mapping");
308 #endif
309
310 nds->W2C = TmNProject(W2C, dv->NDPerm, NULL);
311
312 return nds;
313 }
314
drawer_destroy_ndstuff(NDstuff * nds)315 void drawer_destroy_ndstuff(NDstuff *nds)
316 {
317 TmNDelete(nds->Tc);
318 TmNDelete(nds->W2c);
319
320 TmNDelete(nds->T);
321 TmNDelete(nds->W2C);
322
323 #if 0
324 TmNDelete(nds->rest);
325 OOGLFree(nds->perm);
326 #endif
327
328 HPtNDelete(nds->hc);
329
330 /* ... and finally do not forget to destroy nds. Gnah. */
331 OOGLFree(nds);
332
333 #if 0
334 mgpoptransform();
335 #endif
336 }
337
338 /* install a new object transformation, this is simply a matter of
339 * matrix multiplication; the resulting O2C transfrom will contain the
340 * proper sub-space projection.
341 */
drawer_transform_ndstuff(NDstuff * nds,TransformN * T)342 void drawer_transform_ndstuff(NDstuff *nds, TransformN *T)
343 {
344 nds->T = TmNConcat(T, nds->W2C, nds->T);
345 if (nds->W2c) {
346 nds->Tc = TmNConcat(T, nds->W2c, nds->Tc);
347 }
348
349 #if 0
350 compute_mg_factor(nds);
351 mgtransform(nds->MGFactor);
352 #endif
353 }
354
355 #if 0
356
357 /* The stuff below is incomplete and does not work in its current
358 * shape. See the comment at the start of map_ND_point().
359 */
360
361 /* Factor out a 4x4 matrix such that
362 *
363 * / I \
364 * (nds->T) = P | --- | A
365 * \ B /
366 *
367 * with a permutation matrix P, an (Nd x 4) matrix B and a 4x4 matrix
368 * A. The goal is to hand over A to the MG-layer (and thus possibly to
369 * the hardware). This should speed up the projection process
370 * considerably, because we need 16 FLOPS (matrix multiplication) less
371 * than before. A is added to the object->screen transform which is
372 * applied anyway.
373 */
374 static inline void swap_rows(int r1, int r2, int n, HPtNCoord *Ta, int *perm)
375 {
376 int swapi, col;
377 HPtNCoord swapf;
378
379 if (r1 == r2) {
380 return;
381 }
382
383 /* record permutation */
384 swapi = perm[r1];
385 perm[r1] = perm[r2];
386 perm[r2] = swapi;
387
388 for (col = 0; col < 4; col++) {
389 swapf = Ta[r1*4 + col];
390 Ta[r1*4 + col] = Ta[r2*4 + col];
391 Ta[r2*4 + col] = swapf;
392 }
393 }
394
395 static inline void
396 swap_cols(int rc1, int rc2, int n, HPtNCoord *Ta, Transform A)
397 {
398 HPtNCoord swapf;
399 int row, col;
400
401 if (rc1 == rc2) {
402 return;
403 }
404
405 /* swap columns in T */
406 for (row = 0; row < n; row++) {
407 swapf = Ta[row*4 + rc1];
408 Ta[row*4 + rc1] = Ta[row*4 + rc2];
409 Ta[row*4 + rc2] = swapf;
410 }
411
412 /* swap rows in A */
413 for (col = 0; col < 4; col++) {
414 swapf = A[rc1][col];
415 A[rc1][col] = A[rc2][col];
416 A[rc2][col] = swapf;
417 }
418 }
419
420 #define DEBUG_MG_FACTOR 1
421
422 static void compute_mg_factor(NDstuff *nds)
423 {
424 int idim = nds->T->idim;
425 int i, row, col;
426 TransformN *T;
427 TransformPtr A = nds->MGFactor;
428 int perm[idim];
429 HPtNCoord max, pivot;
430 int maxrow, maxcol;
431
432 T = nds->rest = TmNCopy(nds->T, nds->rest);
433
434 for (i = 0; i < idim; i++) {
435 perm[i] = i;
436 }
437
438 TmCopy(TM_IDENTITY, A);
439
440 for (i = 0; i < 4; i++) {
441 /* find pivot element */
442 max = 0.0;
443 maxrow = maxcol = i;
444 for (row = i; row < idim; row++) {
445 for (col = i; col < 4; col++) {
446 if (fabs(T->a[row*4+col]) > max) {
447 max = fabs(T->a[row*4+col]);
448 maxrow = row;
449 maxcol = col;
450 }
451 }
452 }
453 /* exchange row i and maxrow */
454 swap_rows(i, maxrow, idim, T->a, perm);
455
456 /* exchange column i and maxcol in T and row i and maxcol in A */
457 swap_cols(i, maxcol, idim, T->a, A);
458
459 /* now perform the Gauss' elimination process; note that we have
460 * to perform the operation on the columns of T as T operates from
461 * the right.
462 *
463 * We perform the inverse operation on the rows of A.
464 */
465
466 pivot = T->a[i*4+i];
467
468 /* Clear i-th row of T, record inverse operation in A */
469 for (col = i+1; col < 4; col++) {
470 HPtNCoord factor;
471
472 factor = T->a[i*4+col] / pivot;
473
474 #if DEBUG_MG_FACTOR
475 /* debugging: also perform the operation on the upper part of T */
476 for (row = 0; row < i+1; row++) {
477 T->a[row*4+col] -= T->a[row*4+i] * factor;
478 }
479 #endif
480
481 for (row = i+1; row < idim; row++) {
482 T->a[row*4+col] -= T->a[row*4+i] * factor;
483 }
484 for (row = 0; row < 4; row++) {
485 A[i][row] += A[col][row] * factor;
486 }
487 }
488 }
489
490 /* At this point T has the form
491 *
492 * / L \
493 * (nds->T) = P^{-1} | --- | A'
494 * \ B' /
495 *
496 * with a lower triangular matrix L, proceed now by converting L
497 * to the 4x4 unity matrix.
498 *
499 * We need to apply the operations only on B' and A'
500 */
501 for (col = 4; --col >= 0;) {
502 HPtNCoord factor;
503 int col2;
504
505 for (col2 = 0; col2 < col; col2++) {
506 /* clear the col-th row */
507 factor = T->a[col*4+col2] / T->a[col*4+col];
508 #if DEBUG_MG_FACTOR
509 /* debugging: also perform the operation on the upper part of T */
510 for (row = 0; row < 4; row++) {
511 T->a[row*4+col2] -= T->a[row*4+col] * factor;
512 }
513 #endif
514 for (row = 4; row < idim; row++) {
515 T->a[row*4+col2] -= T->a[row*4+col] * factor;
516 }
517 for (row = 0; row < 4; row++) {
518 A[col][row] += A[col2][row] * factor;
519 }
520 }
521
522 /* finally scale the diagonal such that we really have an
523 * identity matrix
524 */
525 factor = T->a[col*4+col];
526 for (row = 0; row < 4; row++) {
527 A[col][row] *= factor;
528 }
529 factor = 1.0 / T->a[col*4+col];
530 #if DEBUG_MG_FACTOR
531 /* debugging: also perform the operation on the upper part of T */
532 for (row = 0; row < 4; row++) {
533 T->a[row*4+col] *= factor;
534 }
535 #endif
536 for (row = 4; row < idim; row++) {
537 T->a[row*4+col] *= factor;
538 }
539 }
540
541 for (row = 0; row < idim; row++) {
542 nds->perm[row] = perm[row];
543 }
544
545 for (col = 0; col < 4; col++) {
546 HPt3Coord swap;
547 swap = A[0][col];
548 A[0][col] = A[1][col];
549 A[1][col] = A[2][col];
550 A[2][col] = A[3][col];
551 A[3][col] = swap;
552 }
553
554 for (row = 4; row < idim; row++) {
555 HPt3Coord swap;
556
557 swap = T->a[row*4+0];
558 T->a[row*4+0] = T->a[row*4+1];
559 T->a[row*4+1] = T->a[row*4+2];
560 T->a[row*4+2] = T->a[row*4+3];
561 T->a[row*4+3] = swap;
562 }
563
564 #if DEBUG_MG_FACTOR
565 for (row = 0; row < 4; row++) {
566 HPt3Coord swap;
567
568 swap = T->a[row*4+0];
569 T->a[row*4+0] = T->a[row*4+1];
570 T->a[row*4+1] = T->a[row*4+2];
571 T->a[row*4+2] = T->a[row*4+3];
572 T->a[row*4+3] = swap;
573 }
574 fprintf(stderr, "Permuation: ");
575 for (row = 0; row < idim; row++) {
576 fprintf(stderr, "%d ", perm[row]);
577 }
578 fprintf(stderr, "\nRest:\n");
579 TmNPrint(stderr, T);
580 fprintf(stderr, "MG Factor:\n");
581 TmPrint(stderr, A);
582 fprintf(stderr, "Origimal matrix:\n");
583 TmNPrint(stderr, nds->T);
584
585 fprintf(stderr, "Product of factors:\n");
586 {
587 TransformN *prod, *prod2;
588 int i;
589
590 prod = TmNCreate(idim, 4, NULL);
591 for (row = 0; row < idim; row++) {
592 for (col = 0; col < 4; col++) {
593 for (i = 0; i < 4; i++) {
594 prod->a[row*4+col] += T->a[row*4+i]*A[i][col];
595 }
596 }
597 }
598
599 prod2 = TmNCreate(idim, 4, NULL);
600 for (row = 0; row < idim; row++) {
601 for (col = 0; col < 4; col++) {
602 prod2->a[perm[row]*4+col] = prod->a[row*4+col];
603 }
604 }
605 /*TmNPrint(stderr, prod);*/
606 TmNPrint(stderr, prod2);
607 TmNDelete(prod);
608 TmNDelete(prod2);
609 }
610 #endif
611 }
612
613 /* Use the factorization computed above to map "from" to "to"
614 * efficiently
615 */
616 static inline void fast_map_ND_point(NDstuff *nds, HPointN *from, HPoint3 *to)
617 {
618 HPtNCoord *v = from->v;
619 int idim;
620 const int odim = 4;
621 int *perm = nds->perm;
622 int i;
623
624 idim = nds->T->idim > from->dim ? from->dim : nds->T->idim;
625
626 /* The following 4 lines used to consume 16 multiplications */
627 to->x = v[perm[1]];
628 to->y = v[perm[2]];
629 to->z = v[perm[3]];
630 to->w = v[perm[0]];
631 for (i = 4; i < idim; i++) {
632 to->x += v[perm[i]] * nds->rest->a[i*odim+0];
633 to->y += v[perm[i]] * nds->rest->a[i*odim+1];
634 to->z += v[perm[i]] * nds->rest->a[i*odim+2];
635 to->w += v[perm[i]] * nds->rest->a[i*odim+3];
636 }
637
638 /*HPt3Transform(nds->MGFactor, to, to);*/
639
640 #if 1
641 {
642 HPoint3 tp[2];
643 static HPt3Coord max;
644 HPtNTransProj(nds->T, from, tp);
645
646 HPt3Transform(nds->MGFactor, to, tp+1);
647 /**(tp+1) = *to;*/
648
649 if (HPt3Distance(tp, tp+1) > max) {
650 max = HPt3Distance(tp, tp+1);
651 fprintf(stderr, "max dist: %e\n", max);
652 }
653 }
654 #endif
655
656 }
657 #endif
658
659 /*
660 * Local Variables: ***
661 * c-basic-offset: 2 ***
662 * End: ***
663 */
664