1 /***********************************************************************/
2 /* Open Visualization Data Explorer */
3 /* (C) Copyright IBM Corp. 1989,1999 */
4 /* ALL RIGHTS RESERVED */
5 /* This code licensed under the */
6 /* "IBM PUBLIC LICENSE - Open Visualization Data Explorer" */
7 /***********************************************************************/
8
9 #include <dxconfig.h>
10
11
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <string.h>
16 #include <dx/dx.h>
17 #include "tetrasClass.h"
18
19
20 #define SEARCH_LIMIT 4
21
22 #undef DIAG
23
24 #define FALSE 0
25 #define TRUE 1
26
27 typedef struct
28 {
29 ArrayHandle points;
30 TetNeighbors *neighbors;
31 Tetrahedron *tetras;
32 int *visited;
33 int vCount;
34 Point *point;
35 Point *p0, *p1, *p2, *p3;
36 BaryCoord *bary;
37 float fuzz;
38 int fuzzFlag;
39 int limit;
40 Point *mm;
41 } TetraArgs;
42
43 static int _dxfbarycentric_coords(TetraArgs *);
44 static void _dxfCleanup(TetrasInterpolator);
45 static int _dxfInitializeTask(Pointer);
46 static int _dxfInitialize(TetrasInterpolator);
47 static int _dxfbaryExit(int, int, TetraArgs *);
48 static int _dxfTetraSearch(TetraArgs *, int);
49 static int _dxfBBoxTest(TetraArgs *, int);
50
51 #define DIAGNOSTIC(str) \
52 DXMessage("tetrahedra interpolator failure: %s", (str))
53
54 /*
55 * The following are the returned by _dxfbarycentric_coords indicating
56 * whether the sample point lay inside the tetrahedra, within FUZZ of
57 * the tetrahedra, well outside the tetrahedra, or that the tetrahedra
58 * was of zero volume.
59 */
60
61 #define INSIDE 0
62 #define POS_OUT 1
63 #define NEG_OUT -1
64 #define POS_FUZZ 2
65 #define NEG_FUZZ -2
66 #define ZERO_VOLUME 3
67
68 #define WithinFuzz(s) ((s) == POS_FUZZ || (s) == NEG_FUZZ)
69 #define NegSide(s) ((s) == NEG_FUZZ || (s) == NEG_OUT)
70 #define PosSide(s) ((s) == POS_FUZZ || (s) == POS_OUT)
71
72 int
_dxfRecognizeTetras(Field field)73 _dxfRecognizeTetras(Field field)
74 {
75 Array array;
76 Type type;
77 Category cat;
78 int rank, shape[32];
79 int status;
80
81 CHECK(field, CLASS_FIELD);
82
83 status = OK;
84
85 ELT_TYPECHECK(field, "tetrahedra");
86
87 array = (Array)DXGetComponentValue(field, "positions");
88 if (! array)
89 {
90 DIAGNOSTIC("missing positions component");
91 status = ERROR;
92 }
93 else
94 {
95 if (! DXGetArrayInfo(array, NULL, NULL, &cat, &rank, shape))
96 {
97 DIAGNOSTIC("error accessing positions info");
98 status = ERROR;
99 }
100
101 if (cat != CATEGORY_REAL)
102 {
103 DIAGNOSTIC("only CATEGORY_REAL positions supported");
104 status = ERROR;
105 }
106 }
107
108 array = (Array)DXGetComponentValue(field, "connections");
109 if (! array)
110 {
111 DIAGNOSTIC("missing connections component");
112 status = ERROR;
113 }
114 else
115 {
116 if (! DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape))
117 {
118 DIAGNOSTIC("error accessing connections info");
119 status = ERROR;
120 }
121
122 if (type != TYPE_INT)
123 {
124 DIAGNOSTIC("only TYPE_INT connections supported");
125 status = ERROR;
126 }
127
128 if (cat != CATEGORY_REAL)
129 {
130 DIAGNOSTIC("only CATEGORY_REAL connections supported");
131 status = ERROR;
132 }
133 }
134
135 return status;
136 }
137
138 TetrasInterpolator
_dxfNewTetrasInterpolator(Field field,enum interp_init initType,double fuzz,Matrix * m)139 _dxfNewTetrasInterpolator(Field field,
140 enum interp_init initType, double fuzz, Matrix *m)
141 {
142 return (TetrasInterpolator)_dxf_NewTetrasInterpolator(field,
143 initType, fuzz, m, &_dxdtetrasinterpolator_class);
144 }
145
146 TetrasInterpolator
_dxf_NewTetrasInterpolator(Field field,enum interp_init initType,float fuzz,Matrix * m,struct tetrasinterpolator_class * class)147 _dxf_NewTetrasInterpolator(Field field,
148 enum interp_init initType, float fuzz, Matrix *m,
149 struct tetrasinterpolator_class *class)
150 {
151 TetrasInterpolator ti;
152 float *mm, *MM;
153
154 ti = (TetrasInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
155 (struct fieldinterpolator_class *)class);
156
157 if (! ti)
158 return NULL;
159
160 mm = ((Interpolator)ti)->min;
161 MM = ((Interpolator)ti)->max;
162
163 if (((MM[0] - mm[0]) * (MM[1] - mm[1]) * (MM[2] - mm[2])) == 0.0)
164 {
165 DXDelete((Object)ti);
166 return NULL;
167 }
168
169 ti->points = NULL;
170 ti->pointsArray = NULL;
171 ti->data = NULL;
172 ti->dataArray = NULL;
173 ti->mmArray = NULL;
174 ti->tetras = NULL;
175 ti->tetrasArray = NULL;
176 ti->neighbors = NULL;
177 ti->neighborsArray = NULL;
178 ti->visited = NULL;
179
180 ti->vCount = 1;
181
182 if (initType == INTERP_INIT_PARALLEL)
183 {
184 if (! DXAddTask(_dxfInitializeTask, (Pointer)&ti, sizeof(ti), 1.0))
185 {
186 DXDelete((Object)ti);
187 return NULL;
188 }
189 }
190 else if (initType == INTERP_INIT_IMMEDIATE)
191 {
192 if (! _dxfInitialize(ti))
193 {
194 DXDelete((Object)ti);
195 return NULL;
196 }
197 }
198
199 return ti;
200 }
201
202 Object
_dxfTetrasInterpolator_Copy(TetrasInterpolator old,enum _dxd_copy copy)203 _dxfTetrasInterpolator_Copy(TetrasInterpolator old, enum _dxd_copy copy)
204 {
205 TetrasInterpolator new;
206
207 new = (TetrasInterpolator)
208 _dxf_NewObject((struct object_class *)&_dxdtetrasinterpolator_class);
209
210 if (!(_dxf_CopyTetrasInterpolator(new, old, copy)))
211 {
212 DXDelete((Object)new);
213 return NULL;
214 }
215 else
216 return (Object)new;
217 }
218
219 TetrasInterpolator
_dxf_CopyTetrasInterpolator(TetrasInterpolator new,TetrasInterpolator old,enum _dxd_copy copy)220 _dxf_CopyTetrasInterpolator(TetrasInterpolator new,
221 TetrasInterpolator old, enum _dxd_copy copy)
222 {
223 if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
224 (FieldInterpolator)old, copy))
225 return NULL;
226
227 if (! old->fieldInterpolator.initialized)
228 {
229 new->points = NULL;
230 new->tetras = NULL;
231 new->data = NULL;
232 new->neighbors = NULL;
233 new->pointsArray = NULL;
234 new->neighborsArray = NULL;
235 new->tetrasArray = NULL;
236 new->dataArray = NULL;
237 new->mmArray = NULL;
238 new->nPoints = old->nPoints;
239 new->nTetras = old->nTetras;
240 new->nElements = old->nElements;
241 new->visited = NULL;
242 }
243 else
244 {
245 new->pointsArray = (Array)DXReference((Object)old->pointsArray);
246 new->neighborsArray = (Array)DXReference((Object)old->neighborsArray);
247 new->tetrasArray = (Array)DXReference((Object)old->tetrasArray);
248 new->dataArray = (Array)DXReference((Object)old->dataArray);
249 new->mmArray = (Array)DXReference((Object)old->mmArray);
250
251 new->nPoints = old->nPoints;
252 new->nTetras = old->nTetras;
253 new->nElements = old->nElements;
254
255 new->hint = -1;
256
257 new->points = DXCreateArrayHandle(new->pointsArray);
258 if (! new->points)
259 return NULL;
260
261 new->data = DXCreateArrayHandle(new->dataArray);
262 if (! new->data)
263 return NULL;
264
265 if (new->fieldInterpolator.localized)
266 {
267 new->tetras =
268 (Tetrahedron *)DXGetArrayDataLocal(new->tetrasArray);
269 new->neighbors =
270 (TetNeighbors *)DXGetArrayDataLocal(new->neighborsArray);
271 new->mm =(Point *)DXGetArrayDataLocal(new->mmArray);
272 }
273 else
274 {
275 new->tetras =(Tetrahedron *)DXGetArrayData(new->tetrasArray);
276 new->neighbors =
277 (TetNeighbors *)DXGetArrayData(new->neighborsArray);
278 new->mm =(Point *)DXGetArrayData(new->mmArray);
279 }
280
281 if (old->gridFlag)
282 {
283 new->gridFlag = 1;
284 new->grid = _dxfCopySearchGrid(old->grid);
285 if (! new->grid)
286 return NULL;
287 }
288 else
289 {
290 new->gridFlag = 0;
291 }
292
293 }
294
295 new->vCount = 1;
296 new->searchLimit = old->searchLimit;
297
298 new->visited = (int *)DXAllocate(new->nTetras*sizeof(int));
299 if (! new->visited)
300 return NULL;
301
302 memset(new->visited, 0, new->nTetras*sizeof(int));
303
304 if (DXGetError())
305 return NULL;
306
307 return (TetrasInterpolator)new;
308 }
309
310 static int
_dxfInitializeTask(Pointer p)311 _dxfInitializeTask(Pointer p)
312 {
313 return _dxfInitialize(*(TetrasInterpolator *)p);
314 }
315
316 static int
_dxfInitialize(TetrasInterpolator ti)317 _dxfInitialize(TetrasInterpolator ti)
318 {
319 int i, j;
320 Field field;
321 float len, vol;
322 int *t;
323 Point *mm;
324 float fuzz = 0;
325
326 ti->hint = -1;
327
328 field = (Field)((Interpolator)ti)->dataObject;
329
330 /*
331 * De-reference data
332 */
333 ti->dataArray = (Array)DXGetComponentValue(field, "data");
334 if (!ti->dataArray)
335 {
336 DXSetError(ERROR_MISSING_DATA, "#10240", "data");
337 return 0;
338 }
339 DXReference((Object)ti->dataArray);
340
341 ti->data = DXCreateArrayHandle(ti->dataArray);
342 if (! ti->data)
343 return 0;
344
345 DXGetArrayInfo(ti->dataArray, NULL, &((Interpolator)ti)->type,
346 &((Interpolator)ti)->category,
347 &((Interpolator)ti)->rank, ((Interpolator)ti)->shape);
348
349 /*
350 * Get the grid.
351 */
352 ti->pointsArray = (Array)DXGetComponentValue(field, "positions");
353 if (!ti->pointsArray)
354 {
355 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
356 return 0;
357 }
358 DXReference((Object)ti->pointsArray);
359
360 ti->points = DXCreateArrayHandle(ti->pointsArray);
361 if (! ti->points)
362 return ERROR;
363
364 /*
365 * Get the tetrahedra.
366 */
367 ti->tetrasArray = (Array)DXGetComponentValue(field, "connections");
368 if (!ti->tetrasArray)
369 {
370 DXSetError(ERROR_MISSING_DATA, "#10240", "connections");
371 return 0;
372 }
373 DXReference((Object)ti->tetrasArray);
374
375 DXGetArrayInfo(ti->tetrasArray, &ti->nTetras, NULL, NULL, NULL, NULL);
376
377 /*
378 * Don't worry about maintaining shape of input; just determine how
379 * many values to interpolate.
380 */
381 ti->nElements = DXGetItemSize(ti->dataArray) /
382 DXTypeSize(((Interpolator)ti)->type);
383
384 /*
385 * Get the neighbors array
386 */
387 ti->neighborsArray = DXNeighbors((Field)((Interpolator)ti)->dataObject);
388 if (!ti->neighborsArray)
389 return 0;
390
391 DXReference((Object)ti->neighborsArray);
392
393 if (ti->fieldInterpolator.localized)
394 ti->neighbors = (TetNeighbors *)DXGetArrayDataLocal(ti->neighborsArray);
395 else
396 ti->neighbors = (TetNeighbors *)DXGetArrayData(ti->neighborsArray);
397
398 if (ti->fieldInterpolator.localized)
399 ti->tetras = (Tetrahedron *)DXGetArrayDataLocal(ti->tetrasArray);
400 else
401 ti->tetras = (Tetrahedron *)DXGetArrayData(ti->tetrasArray);
402
403 ti->mmArray = DXNewArray(TYPE_FLOAT, CATEGORY_REAL, 1, 3);
404 if (! ti->mmArray)
405 return 0;
406
407 DXReference((Object)ti->mmArray);
408
409 if (! DXAddArrayData(ti->mmArray, 0, 2*ti->nTetras, NULL))
410 return 0;
411
412 ti->mm = (Point *)DXGetArrayData(ti->mmArray);
413
414 ti->grid = _dxfMakeSearchGrid(((Interpolator)ti)->max,
415 ((Interpolator)ti)->min, ti->nTetras, 3);
416 if (! ti->grid)
417 return 0;
418
419 if (ti->nTetras)
420 {
421 vol = 1.0;
422 for (i = 0; i < 3; i++)
423 vol *= ((Interpolator)ti)->max[i] - ((Interpolator)ti)->min[i];
424
425 vol /= ti->nTetras;
426
427 len = pow((double)vol, 0.3333333);
428
429 fuzz = ti->fieldInterpolator.fuzz * len;
430 }
431
432 t = (int *)ti->tetras;
433
434 mm = ti->mm;
435 for (i = 0; i < ti->nTetras; i++, t += 4)
436 {
437 float *v;
438 register float x,X,y,Y,z,Z;
439 float *pt[4];
440 Point pBuf[4];
441
442 pt[0] = v = (float *)DXGetArrayEntry(ti->points, t[0],
443 (Pointer)(pBuf));
444
445 x = X = *v++;
446 y = Y = *v++;
447 z = Z = *v++;
448
449 for (j = 1; j < 4; j++)
450 {
451 register float tt;
452
453 pt[j] = v = (float *)DXGetArrayEntry(ti->points, t[j],
454 (Pointer)(pBuf + j));
455
456 tt = *v++;
457 if (x > tt) x = tt;
458 if (X < tt) X = tt;
459
460 tt = *v++;
461 if (y > tt) y = tt;
462 if (Y < tt) Y = tt;
463
464 tt = *v;
465 if (z > tt) z = tt;
466 if (Z < tt) Z = tt;
467 }
468
469 mm->x = x - fuzz;
470 mm->y = y - fuzz;
471 mm->z = z - fuzz;
472 mm++;
473
474 mm->x = X + fuzz;
475 mm->y = Y + fuzz;
476 mm->z = Z + fuzz;
477 mm++;
478
479 if (! _dxfAddItemToSearchGrid(ti->grid, pt, 4, i))
480 break;
481 }
482
483 if (i != ti->nTetras)
484 return 0;
485
486 ti->gridFlag = 1;
487
488 ti->visited = (int *)DXAllocate(ti->nTetras*sizeof(int));
489 if (! ti->visited)
490 return ERROR;
491
492 memset(ti->visited, 0, ti->nTetras*sizeof(int));
493
494 ti->fieldInterpolator.initialized = 1;
495
496 return OK;
497 }
498
499 Error
_dxfTetrasInterpolator_Delete(TetrasInterpolator ti)500 _dxfTetrasInterpolator_Delete(TetrasInterpolator ti)
501 {
502 _dxfCleanup(ti);
503 _dxfFieldInterpolator_Delete((FieldInterpolator) ti);
504
505 return OK;
506 }
507
508 #ifdef DIAG
509 static int nT;
510 static int nB;
511 #endif
512
513 static int
_dxfBBoxTest(TetraArgs * args,int i)514 _dxfBBoxTest(TetraArgs *args, int i)
515 {
516 float *mm, *p;
517
518 mm = (float *)(args->mm + 2*i);
519 p = (float *) args->point;
520
521 if (args->fuzzFlag)
522 {
523 float f;
524
525 f = args->fuzz * (mm[3] - mm[0]);
526 if (*p < mm[0]-f || *p > mm[3]+f) goto out;
527 p++; mm++;
528
529 f = args->fuzz * (mm[3] - mm[0]);
530 if (*p < mm[0]-f || *p > mm[3]+f) goto out;
531 p++; mm++;
532
533 f = args->fuzz * (mm[3] - mm[0]);
534 if (*p < mm[0]-f || *p > mm[3]+f) goto out;
535 }
536 else
537 {
538 if (*p < mm[0] || *p > mm[3]) goto out;
539 p++; mm++;
540
541 if (*p < mm[0] || *p > mm[3]) goto out;
542 p++; mm++;
543
544 if (*p < mm[0] || *p > mm[3]) goto out;
545 }
546
547 return 1;
548
549 out:
550
551 #ifdef DIAG
552 nB ++;
553 #endif
554
555 return 0;
556 }
557
558 int
_dxfTetrasInterpolator_PrimitiveInterpolate(TetrasInterpolator ti,int * n,float ** points,Pointer * values,int fuzzFlag)559 _dxfTetrasInterpolator_PrimitiveInterpolate(TetrasInterpolator ti, int *n,
560 float **points, Pointer *values, int fuzzFlag)
561 {
562 BaryCoord bary;
563 int i;
564 int found;
565 Pointer v;
566 Point *p;
567 float xMin, xMax;
568 float yMin, yMax;
569 float zMin, zMax;
570 float fuzz;
571 TetraArgs args;
572 int dep;
573 int itemSize;
574 ubyte *dbuf;
575 InvalidComponentHandle icH = ((FieldInterpolator)ti)->invCon;
576 Matrix *xform;
577
578 if (!((FieldInterpolator)ti)->initialized)
579 {
580 if (! _dxfInitialize(ti))
581 {
582 _dxfCleanup(ti);
583 return 0;
584 }
585
586 ((FieldInterpolator)ti)->initialized = 1;
587 }
588
589 xMin = ((Interpolator)ti)->min[0];
590 yMin = ((Interpolator)ti)->min[1];
591 zMin = ((Interpolator)ti)->min[2];
592
593 xMax = ((Interpolator)ti)->max[0];
594 yMax = ((Interpolator)ti)->max[1];
595 zMax = ((Interpolator)ti)->max[2];
596
597 v = (Pointer)*values;
598 p = (Point *)*points;
599
600 fuzz = ((FieldInterpolator)ti)->fuzz;
601
602 args.fuzz = ti->fieldInterpolator.fuzz;
603 args.bary = &bary;
604 args.points = ti->points;
605 args.tetras = ti->tetras;
606 args.neighbors = ti->neighbors;
607 args.visited = ti->visited;
608 args.mm = ti->mm;
609
610 dep = ti->fieldInterpolator.data_dependency;
611
612 itemSize = ti->nElements*DXTypeSize(((Interpolator)ti)->type);
613
614 dbuf = (ubyte *)DXAllocate(itemSize*4);
615 if (! dbuf)
616 return 0;
617
618 if (((FieldInterpolator)ti)->xflag)
619 xform = &(((FieldInterpolator)ti)->xform);
620 else
621 xform = NULL;
622
623 /*
624 * For each point in the input, attempt to interpolate the point.
625 * When a point cannot be interpolated, quit.
626 */
627 while(*n != 0)
628 {
629 Point xpt;
630 Point *pPtr;
631
632 if (xform)
633 {
634 xpt.x = p->x*xform->A[0][0] +
635 p->y*xform->A[1][0] +
636 p->z*xform->A[2][0] +
637 xform->b[0];
638
639 xpt.y = p->x*xform->A[0][1] +
640 p->y*xform->A[1][1] +
641 p->z*xform->A[2][1] +
642 xform->b[1];
643
644 xpt.z = p->x*xform->A[0][2] +
645 p->y*xform->A[1][2] +
646 p->z*xform->A[2][2] +
647 xform->b[2];
648
649 pPtr = &xpt;
650 }
651 else
652 pPtr = p;
653
654
655 if (fuzzFlag == FUZZ_ON)
656 {
657 if ((pPtr->x > xMax+fuzz ||
658 pPtr->x < xMin-fuzz) ||
659 (pPtr->y > yMax+fuzz ||
660 pPtr->y < yMin-fuzz) ||
661 (pPtr->z > zMax+fuzz ||
662 pPtr->z < zMin-fuzz))
663 break;
664 }
665 else
666 {
667 if ((pPtr->x > xMax ||
668 pPtr->x < xMin) ||
669 (pPtr->y > yMax ||
670 pPtr->y < yMin) ||
671 (pPtr->z > zMax ||
672 pPtr->z < zMin))
673 break;
674 }
675
676 #ifdef DIAG
677 nT = 0;
678 nB = 0;
679 #endif
680 found = -1;
681
682 args.fuzzFlag = fuzzFlag;
683 args.point = pPtr;
684 args.vCount = ti->vCount ++;
685 args.limit = SEARCH_LIMIT;
686
687 /*
688 * Check the hint first, if one exists...
689 */
690 if (ti->hint != -1 && _dxfBBoxTest(&args, ti->hint))
691 {
692 found = _dxfTetraSearch(&args, ti->hint);
693 }
694
695 /*
696 * If it wasn't found, try the search grid.
697 */
698 if (found == -1)
699 {
700 if (ti->gridFlag)
701 {
702 int primNum;
703
704 _dxfInitGetSearchGrid(ti->grid, (float *)pPtr);
705 while((primNum = _dxfGetNextSearchGrid(ti->grid)) != -1
706 && found == -1)
707 if (ti->visited[primNum] != ti->vCount &&
708 _dxfBBoxTest(&args, primNum))
709 found = _dxfTetraSearch(&args, primNum);
710 }
711 }
712
713 if (found == -1 || (icH && DXIsElementInvalid(icH, found)))
714 break;
715
716 ti->hint = found;
717
718 #define INTERPOLATE(type, round) \
719 { \
720 type *d0, *d1, *d2, *d3, *r; \
721 Tetrahedron *tet; \
722 \
723 tet = ti->tetras + found; \
724 \
725 r = (type *)v; \
726 \
727 d0 = (type *)DXGetArrayEntry(ti->data, \
728 tet->p, (Pointer)(dbuf)); \
729 \
730 d1 = (type *)DXGetArrayEntry(ti->data, \
731 tet->q, (Pointer)(dbuf+itemSize)); \
732 \
733 d2 = (type *)DXGetArrayEntry(ti->data, \
734 tet->r, (Pointer)(dbuf+2*itemSize)); \
735 \
736 d3 = (type *)DXGetArrayEntry(ti->data, \
737 tet->s, (Pointer)(dbuf+3*itemSize)); \
738 \
739 for (i = 0; i < ti->nElements; i++) \
740 *r++ = (*d0++ * bary.p) + (*d1++ * bary.q) + \
741 (*d2++ * bary.r) + (*d3++ * bary.s) + \
742 round; \
743 \
744 v = (Pointer)r; \
745 }
746
747
748
749 if (((FieldInterpolator)ti)->cstData)
750 {
751 memcpy(v, ((FieldInterpolator)ti)->cstData, itemSize);
752 v = (Pointer)(((char *)v) + itemSize);
753 }
754 else if (dep == DATA_POSITIONS_DEPENDENT)
755 {
756 Type dataType;
757
758 if ((dataType = ((Interpolator)ti)->type) == TYPE_FLOAT)
759 {
760 INTERPOLATE(float, 0.0);
761 }
762 else if (dataType == TYPE_DOUBLE)
763 {
764 INTERPOLATE(double, 0.0);
765 }
766 else if (dataType == TYPE_INT)
767 {
768 INTERPOLATE(int, 0.5);
769 }
770 else if (dataType == TYPE_SHORT)
771 {
772 INTERPOLATE(short, 0.5);
773 }
774 else if (dataType == TYPE_USHORT)
775 {
776 INTERPOLATE(ushort, 0.5);
777 }
778 else if (dataType == TYPE_UINT)
779 {
780 INTERPOLATE(uint, 0.5);
781 }
782 else if (dataType == TYPE_BYTE)
783 {
784 INTERPOLATE(byte, 0.5);
785 }
786 else if (dataType == TYPE_UBYTE)
787 {
788 INTERPOLATE(ubyte, 0.5);
789 }
790 else
791 {
792 INTERPOLATE(unsigned char, 0.5);
793 }
794 }
795 else
796 {
797 memcpy(v, DXGetArrayEntry(ti->data,
798 found, (Pointer)dbuf), itemSize);
799 v = (Pointer)(((char *)v) + itemSize);
800 }
801
802
803 /*
804 * Only use fuzz on first point
805 */
806 fuzzFlag = FUZZ_OFF;
807
808 p ++;
809 *n -= 1;
810 }
811
812 DXFree((Pointer)dbuf);
813
814 *values = (Pointer)v;
815 *points = (float *)p;
816
817 return OK;
818 }
819
820 static void
_dxfCleanup(TetrasInterpolator ti)821 _dxfCleanup(TetrasInterpolator ti)
822 {
823 /*
824 * Note: only free up the array data if its allocated
825 * in local memory
826 */
827
828 if (((FieldInterpolator)ti)->localized)
829 {
830 if (ti->tetras)
831 {
832 DXFreeArrayDataLocal(ti->tetrasArray, (Pointer)ti->tetras);
833 ti->tetras = NULL;
834 }
835
836 if (ti->neighbors)
837 {
838 DXFreeArrayDataLocal(ti->neighborsArray, (Pointer)ti->neighbors);
839 ti->neighbors = NULL;
840 }
841 }
842
843 if (ti->points)
844 {
845 DXFreeArrayHandle(ti->points);
846 ti->points = NULL;
847 }
848
849 if (ti->data)
850 {
851 DXFreeArrayHandle(ti->data);
852 ti->data = NULL;
853 }
854
855 if (ti->pointsArray)
856 {
857 DXDelete((Object)ti->pointsArray);
858 ti->pointsArray = NULL;
859 }
860
861 if (ti->dataArray)
862 {
863 DXDelete((Object)ti->dataArray);
864 ti->dataArray = NULL;
865 }
866
867 if (ti->tetrasArray)
868 {
869 DXDelete((Object)ti->tetrasArray);
870 ti->tetrasArray = NULL;
871 }
872
873 if (ti->neighborsArray)
874 {
875 DXDelete((Object)ti->neighborsArray);
876 ti->neighborsArray = NULL;
877 }
878
879 if (ti->gridFlag)
880 _dxfFreeSearchGrid(ti->grid);
881
882 if (ti->visited)
883 {
884 DXFree((Pointer)ti->visited);
885 ti->visited = NULL;
886 }
887
888 if (ti->mmArray)
889 {
890 DXDelete((Object)ti->mmArray);
891 ti->mm = NULL;
892 }
893
894 }
895
896 /*
897 * This procedure determines the barycentric coordinates of the point pt with
898 * respect to the points p0, p1, p2, and p3. It does so by finding five
899 * volumes: vol, vol0, vol1, vol2, vol3. The code to determine these
900 * volumes could be replaced by the following 5 lines of code:
901 * vol0 = tetra_volume (pt, p1, p2, p3);
902 * vol1 = tetra_volume (p0, pt, p2, p3);
903 * vol2 = tetra_volume (p0, p1, pt, p3);
904 * vol3 = tetra_volume (p0, p1, p2, pt);
905 * vol = vol0 + vol1 + vol2 + vol3;
906 * I.e., vol is the scaled oriented volume of the original tetra, and voli
907 * is the scaled oriented volume of the tetra formed by substituting the point
908 * pt for the point pi in the ordering p0, p1, p2, p3. If any of the latter
909 * four volumes have orientation opposite to the orientation of the original
910 * tetra, then the points pt and pi are on opposite sides of a face of the
911 * tetrahedron, and therefore pt is outside of the tetrahedron. Since the
912 * above volumes are calculated with VERY similar determinants, this
913 * procedure takes all the common multiplies and subtracts from all those
914 * determinants, and calculates them each only once. Using 5 separate volume
915 * procedure calls, this would take 36 multiplies and 56 add/subtracts. Using
916 * this procedure, this takes 32 multiplies and 39 add/subtracts, which should
917 * be a better than 10% savings on one of the most costly procedures (because
918 * of the number of times it gets called) in the whole streamline algorithm.
919 */
920
_dxfbarycentric_coords(TetraArgs * args)921 static int _dxfbarycentric_coords (TetraArgs *args)
922 {
923 float vol0, vol1, vol2, vol3, vol;
924 float xt, x0, x1, x2, x3;
925 float yt, y0, y1, y2, y3;
926 float zt, z0, z1, z2, z3;
927 float x01, x02, x03, x0t, xt1, xt2, xt3;
928 float y01, y02, y03, y0t, yt1, yt2, yt3;
929 float z01, z02, z03, z0t, zt1, zt2, zt3;
930 float yz012, yz013, yz023, yzt12, yzt13;
931 float yzt23, yz0t2, yz0t3, yz01t, yz02t;
932 Point *pt, *p0, *p1, *p2, *p3;
933 BaryCoord *b;
934 float fuzz;
935
936 pt = args->point;
937 p0 = args->p0;
938 p1 = args->p1;
939 p2 = args->p2;
940 p3 = args->p3;
941 b = args->bary;
942
943 xt = pt->x; yt = pt->y; zt = pt->z;
944 x0 = p0->x; y0 = p0->y; z0 = p0->z;
945 x1 = p1->x; y1 = p1->y; z1 = p1->z;
946 x2 = p2->x; y2 = p2->y; z2 = p2->z;
947 x3 = p3->x; y3 = p3->y; z3 = p3->z;
948
949 x01 = x1 - x0; y01 = y1 - y0; z01 = z1 - z0;
950 x02 = x2 - x0; y02 = y2 - y0; z02 = z2 - z0;
951 x03 = x3 - x0; y03 = y3 - y0; z03 = z3 - z0;
952 x0t = xt - x0; y0t = yt - y0; z0t = zt - z0;
953 xt1 = x1 - xt; yt1 = y1 - yt; zt1 = z1 - zt;
954 xt2 = x2 - xt; yt2 = y2 - yt; zt2 = z2 - zt;
955 xt3 = x3 - xt; yt3 = y3 - yt; zt3 = z3 - zt;
956
957 yz012 = y01 * z02 - y02 * z01;
958 yz013 = y01 * z03 - y03 * z01;
959 yz023 = y02 * z03 - y03 * z02;
960 yzt12 = yt1 * zt2 - yt2 * zt1;
961 yzt13 = yt1 * zt3 - yt3 * zt1;
962 yzt23 = yt2 * zt3 - yt3 * zt2;
963 yz0t2 = y0t * z02 - y02 * z0t;
964 yz0t3 = y0t * z03 - y03 * z0t;
965 yz01t = y01 * z0t - y0t * z01;
966 yz02t = y02 * z0t - y0t * z02;
967
968 vol0 = xt1 * yzt23 - xt2 * yzt13 + xt3 * yzt12;
969 vol1 = x0t * yz023 - x02 * yz0t3 + x03 * yz0t2;
970 vol2 = x01 * yz0t3 - x0t * yz013 + x03 * yz01t;
971 vol3 = x01 * yz02t - x02 * yz01t + x0t * yz012;
972
973 vol = vol0 + vol1 + vol2 + vol3;
974
975 if (vol == 0.0)
976 return ZERO_VOLUME;
977
978 b->p = vol0;
979 b->q = vol1;
980 b->r = vol2;
981 b->s = vol3;
982
983 #if 1
984 if (((b->p >= 0.0)&&(b->q >= 0.0)&&(b->r >= 0.0)&&(b->s >= 0.0)) ||
985 ((b->p <= 0.0)&&(b->q <= 0.0)&&(b->r <= 0.0)&&(b->s <= 0.0)))
986 #else
987 if (((b->p >= 0.0)&&(b->q >= 0.0)&&(b->r >= 0.0)&&(b->s >= 0.0)) ||
988 ((b->p <= 0.0)&&(b->q <= 0.0)&&(b->r <= 0.0)&&(b->s <= 0.0)) ||
989 ((b->p == 0.0)||(b->q == 0.0)||(b->r == 0.0)||(b->s == 0.0)))
990 #endif
991 {
992 vol = 1.0 / vol;
993 b->p *= vol;
994 b->q *= vol;
995 b->r *= vol;
996 b->s *= vol;
997 return INSIDE;
998 }
999
1000 fuzz = args->fuzz * vol;
1001
1002 if (((b->p >= -fuzz)&&(b->q >= -fuzz)&&(b->r >= -fuzz)&&(b->s >= -fuzz)) ||
1003 ((b->p <= fuzz)&&(b->q <= fuzz)&&(b->r <= fuzz)&&(b->s <= fuzz)))
1004 {
1005 vol = 1.0 / vol;
1006 b->p *= vol;
1007 b->q *= vol;
1008 b->r *= vol;
1009 b->s *= vol;
1010 return POS_FUZZ;
1011 }
1012
1013 if (vol > 0)
1014 return POS_OUT;
1015 else
1016 return NEG_OUT;
1017 }
1018
1019
1020 static int
_dxfbaryExit(int side,int index,TetraArgs * args)1021 _dxfbaryExit(int side, int index, TetraArgs *args)
1022 {
1023 float best; /* best choice so far */
1024 int f;
1025 int *nbrs;
1026 float *be;
1027 int *v, vc;
1028 int n;
1029 float d;
1030
1031 nbrs = (int *)(args->neighbors + index);
1032 be = (float *)args->bary;
1033 v = args->visited;
1034 vc = args->vCount;
1035
1036 f = -1;
1037
1038 if (NegSide(side))
1039 {
1040 best = -DXD_MAX_FLOAT;
1041
1042 n = *nbrs++; d = *be++;
1043 if (n >= 0 && v[n] != vc && d > 0.0 && d > best)
1044 {
1045 best = d; f = n;
1046 }
1047
1048 n = *nbrs++; d = *be++;
1049 if (n >= 0 && v[n] != vc && d > 0.0 && d > best)
1050 {
1051 best = d; f = n;
1052 }
1053
1054 n = *nbrs++; d = *be++;
1055 if (n >= 0 && v[n] != vc && d > 0.0 && d > best)
1056 {
1057 best = d; f = n;
1058 }
1059
1060 n = *nbrs; d = *be;
1061 if (n >= 0 && v[n] != vc && d > 0.0 && d > best)
1062 {
1063 best = d; f = n;
1064 }
1065 }
1066 else
1067 {
1068 best = DXD_MAX_FLOAT;
1069
1070 n = *nbrs++; d = *be++;
1071 if (n >= 0 && v[n] != vc && d < 0.0 && d < best)
1072 {
1073 best = d; f = n;
1074 }
1075
1076 n = *nbrs++; d = *be++;
1077 if (n >= 0 && v[n] != vc && d < 0.0 && d < best)
1078 {
1079 best = d; f = n;
1080 }
1081
1082 n = *nbrs++; d = *be++;
1083 if (n >= 0 && v[n] != vc && d < 0.0 && d < best)
1084 {
1085 best = d; f = n;
1086 }
1087
1088 n = *nbrs++; d = *be++;
1089 if (n >= 0 && v[n] != vc && d < 0.0 && d < best)
1090 {
1091 best = d; f = n;
1092 }
1093 }
1094
1095 return f;
1096 }
1097
1098 static int
_dxfTetraSearch(TetraArgs * args,int index)1099 _dxfTetraSearch(TetraArgs *args, int index)
1100 {
1101 int nbr;
1102 int side;
1103 Tetrahedron *tet;
1104 int knt;
1105 int *vPtr;
1106 ArrayHandle points;
1107 Tetrahedron *tetras;
1108 int *visited, vCount;
1109 int limit;
1110 int fuzzFlag;
1111 Point ptBuf[4];
1112
1113 points = args->points;
1114 tetras = args->tetras;
1115 visited = args->visited;
1116 vCount = args->vCount;
1117 limit = args->limit;
1118 fuzzFlag = args->fuzzFlag;
1119
1120 nbr = 0;
1121
1122 for (knt = 0; index != -1 && knt < limit; knt++, index = nbr)
1123 {
1124 if (visited)
1125 {
1126 vPtr = visited + index;
1127
1128 if (*vPtr == vCount)
1129 break;
1130 else
1131 *vPtr = vCount;
1132 }
1133
1134 tet = tetras + index;
1135
1136 args->p0 = (Point *)
1137 DXGetArrayEntry(points, tet->p, (Pointer)(ptBuf+0));
1138 args->p1 = (Point *)
1139 DXGetArrayEntry(points, tet->q, (Pointer)(ptBuf+1));
1140 args->p2 = (Point *)
1141 DXGetArrayEntry(points, tet->r, (Pointer)(ptBuf+2));
1142 args->p3 = (Point *)
1143 DXGetArrayEntry(points, tet->s, (Pointer)(ptBuf+3));
1144
1145 #ifdef DIAG
1146 nT ++;
1147 #endif
1148
1149 side = _dxfbarycentric_coords(args);
1150
1151 /*
1152 * If the point was contained in the tetrahedra, return its index.
1153 * Otherwise, if its within fuzz of the tetrahedra and the tetrahedra
1154 * has a neighbor in that direction, continue. If its within fuzz of
1155 * the primitive but the tetra has no neighbor in that direction,
1156 * return the primitive's index. If the point is not within fuzz of
1157 * the primitive, continue the neighbor search.
1158 *
1159 * But first, if the tetra is of zero volume, neighbor walking will
1160 * not work.
1161 */
1162 if (side == ZERO_VOLUME)
1163 break;
1164
1165 if (side == INSIDE)
1166 return index;
1167
1168 /*
1169 * If the fuzz flag is on and we are within fuzz and there's
1170 * no neighbor, this one will do.
1171 */
1172 if (fuzzFlag == FUZZ_ON && WithinFuzz(side))
1173 return index;
1174
1175 /*
1176 * Otherwise, if there's no neighbor, we failed. If there is,
1177 * we keep going.
1178 */
1179
1180 if (knt < (limit-1))
1181 {
1182 nbr = _dxfbaryExit(side, index, args);
1183 if (nbr < 0)
1184 return -1;
1185 }
1186
1187 if (! _dxfBBoxTest(args, nbr))
1188 return -1;
1189 }
1190
1191 return -1;
1192 }
1193
1194 Interpolator
_dxfTetrasInterpolator_LocalizeInterpolator(TetrasInterpolator ti)1195 _dxfTetrasInterpolator_LocalizeInterpolator(TetrasInterpolator ti)
1196 {
1197 if (ti->fieldInterpolator.localized)
1198 return (Interpolator)ti;
1199
1200 ti->fieldInterpolator.localized = 1;
1201
1202 if (ti->fieldInterpolator.initialized)
1203 {
1204 ti->data = DXGetArrayDataLocal(ti->dataArray);
1205 ti->tetras = (Tetrahedron *)DXGetArrayDataLocal(ti->tetrasArray);
1206 ti->neighbors = (TetNeighbors *)DXGetArrayDataLocal(ti->neighborsArray);
1207 }
1208
1209 if (DXGetError())
1210 return NULL;
1211 else
1212 return (Interpolator)ti;
1213 }
1214