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 "cubesIIClass.h"
18
19 #define FALSE 0
20 #define TRUE 1
21
22 #define SEARCH_LIMIT 10
23
24 typedef struct
25 {
26 ArrayHandle pHandle;
27 int *nbrs;
28 ArrayHandle cHandle;
29 int *visited;
30 int vCount;
31 Point *point;
32 Point *p0, *p1, *p2, *p3;
33 int tetra[4];
34 float *bary;
35 float fuzz;
36 int fuzzFlag;
37 int limit;
38 int *gridCounts;
39 int Cyz, Cz;
40 CubesIIInterpolator ci;
41 } CubesIIArgs;
42
43 #define NBR_0 0x01
44 #define NBR_1 0x02
45 #define NBR_2 0x04
46 #define NBR_3 0x08
47 #define NBR_4 0x10
48 #define NBR_5 0x20
49
50 static Error _dxfInitializeTask(Pointer);
51 static Error _dxfInitialize(CubesIIInterpolator);
52 static void _dxfCleanup(CubesIIInterpolator);
53 static int CubesIISearch(CubesIIArgs *, int);
54 static int IsInCube(CubesIIArgs *, int);
55 static int IsInTetra(CubesIIArgs *);
56 static int GetRegularNeighbor(CubesIIArgs *, int, int);
57
58 #define DIAGNOSTIC(str) \
59 DXMessage("cubes interpolator failure: %s", (str))
60
61 /*
62 * Return values from IsInTetra, IsInCube (plus 0-5)
63 */
64 #define INSIDE -2
65 #define OUTSIDE -1
66
67 int
_dxfRecognizeCubesII(Field field)68 _dxfRecognizeCubesII(Field field)
69 {
70 Array array;
71 Type type;
72 Category cat;
73 int rank, shape[32];
74 int status;
75
76 CHECK(field, CLASS_FIELD);
77
78 ELT_TYPECHECK(field, "cubes");
79
80 status = OK;
81
82 array = (Array)DXGetComponentValue(field, "positions");
83 if (!array)
84 {
85 DIAGNOSTIC("missing positions component");
86 status = ERROR;
87 }
88 else
89 {
90 if (! DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape))
91 {
92 DIAGNOSTIC("error accessing positions info");
93 status = ERROR;
94 }
95
96 if (cat != CATEGORY_REAL)
97 {
98 DIAGNOSTIC("only CATEGORY_REAL positions supported");
99 status = ERROR;
100 }
101 }
102
103 array = (Array)DXGetComponentValue(field, "connections");
104 if (!array)
105 {
106 DIAGNOSTIC("missing connections component");
107 status = ERROR;
108 }
109 else
110 {
111 if (! DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape))
112 {
113 DIAGNOSTIC("error accessing connections info");
114 status = ERROR;
115 }
116
117 if (type != TYPE_INT)
118 {
119 DIAGNOSTIC("only TYPE_INT connections supported");
120 status = ERROR;
121 }
122
123 if (cat != CATEGORY_REAL)
124 {
125 DIAGNOSTIC("only CATEGORY_REAL connections supported");
126 status = ERROR;
127 }
128 }
129
130 return status;
131 }
132
133 CubesIIInterpolator
_dxfNewCubesIIInterpolator(Field field,enum interp_init initType,double fuzz,Matrix * m)134 _dxfNewCubesIIInterpolator(Field field,
135 enum interp_init initType, double fuzz, Matrix *m)
136 {
137 return (CubesIIInterpolator)_dxf_NewCubesIIInterpolator(field,
138 initType, fuzz, m, &_dxdcubesiiinterpolator_class);
139 }
140
141 CubesIIInterpolator
_dxf_NewCubesIIInterpolator(Field field,enum interp_init initType,float fuzz,Matrix * m,struct cubesiiinterpolator_class * class)142 _dxf_NewCubesIIInterpolator(Field field,
143 enum interp_init initType, float fuzz, Matrix *m,
144 struct cubesiiinterpolator_class *class)
145 {
146 CubesIIInterpolator ci;
147
148 if (DXEmptyField(field))
149 return NULL;
150
151 ci = (CubesIIInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
152 (struct fieldinterpolator_class *)class);
153 if (! ci)
154 return NULL;
155
156 ci->pHandle = NULL;
157 ci->pointsArray = NULL;
158 ci->dHandle = NULL;
159 ci->dataArray = NULL;
160 ci->cHandle = NULL;
161 ci->cubesArray = NULL;
162 ci->visited = NULL;
163 ci->grid = NULL;
164 ci->vCount = 1;
165
166 if (initType == INTERP_INIT_PARALLEL)
167 {
168 if (! DXAddTask(_dxfInitializeTask, (Pointer)&ci, sizeof(ci), 1.0))
169 {
170 DXDelete((Object)ci);
171 return NULL;
172 }
173 }
174 else if (initType == INTERP_INIT_IMMEDIATE)
175 {
176 if (! _dxfInitialize(ci))
177 {
178 DXDelete((Object)ci);
179 return NULL;
180 }
181 }
182
183 ci->hint = -1;
184
185 return ci;
186 }
187
188 static Error
_dxfInitializeTask(Pointer p)189 _dxfInitializeTask(Pointer p)
190 {
191 int status;
192
193 status = _dxfInitialize(*(CubesIIInterpolator *)p);
194 return status;
195 }
196
197 #define ISOBUF 32
198
199 static Error
_dxfInitialize(CubesIIInterpolator ci)200 _dxfInitialize(CubesIIInterpolator ci)
201 {
202 Field field;
203 Type dataType;
204 Category dataCategory;
205 register int i;
206
207 field = (Field) ((Interpolator)ci)->dataObject;
208
209 ci->fieldInterpolator.initialized = 1;
210
211 /*
212 * De-reference data
213 */
214 ci->dataArray = (Array)DXGetComponentValue(field, "data");
215 if (!ci->dataArray)
216 {
217 DXSetError(ERROR_MISSING_DATA, "#10240", "data");
218 return ERROR;
219 }
220
221 DXGetArrayInfo(ci->dataArray, NULL, &((Interpolator)ci)->type,
222 &((Interpolator)ci)->category, &((Interpolator)ci)->rank,
223 ((Interpolator)ci)->shape);
224
225 ci->dHandle = DXCreateArrayHandle(ci->dataArray);
226 if (! ci->dHandle)
227 return ERROR;
228
229 /*
230 * Get the grid.
231 */
232 ci->pointsArray = (Array)DXGetComponentValue(field, "positions");
233 if (!ci->pointsArray)
234 {
235 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
236 return ERROR;
237 }
238
239 ci->pHandle = DXCreateArrayHandle(ci->pointsArray);
240 if (! ci->pHandle)
241 return ERROR;
242
243 /*
244 * Get the cubes.
245 */
246 ci->cubesArray = (Array)DXGetComponentValue(field, "connections");
247 if (!ci->cubesArray)
248 {
249 DXSetError(ERROR_MISSING_DATA, "#10240", "connections");
250 return ERROR;
251 }
252
253 if (! DXGetArrayInfo(ci->cubesArray, &ci->nCubes, NULL, NULL, NULL, NULL))
254 return ERROR;
255
256 if (NULL == (ci->cHandle = DXCreateArrayHandle(ci->cubesArray)))
257 return ERROR;
258
259 /*
260 * If cubes are a mesh array (eg. regular connections), get the mesh info.
261 * Otherwise, get the neighbors, generating one if one exists. If NULL
262 * is returned, well, OK.
263 */
264 if (DXQueryGridConnections(ci->cubesArray, NULL, ci->gridCounts))
265 {
266 ci->nbrsArray = NULL;
267 ci->gridCounts[0] -= 1;
268 ci->gridCounts[1] -= 1;
269 ci->gridCounts[2] -= 1;
270 ci->Cyz = ci->gridCounts[1]*ci->gridCounts[2];
271 ci->Cz = ci->gridCounts[2];
272 }
273 else
274 ci->nbrsArray = (Array)DXNeighbors(field);
275
276 /*
277 * get info about data values
278 */
279 DXGetArrayInfo(ci->dataArray, NULL, &dataType, &dataCategory, NULL, NULL);
280
281 /*
282 * Don't worry about maintaining shape of input; just determine how
283 * many values to interpolate.
284 */
285 ci->nElements = DXGetItemSize(ci->dataArray) / DXTypeSize(dataType);
286
287 if (ci->fieldInterpolator.localized)
288 {
289 if (ci->nbrsArray)
290 ci->nbrs = (int *)DXGetArrayDataLocal(ci->nbrsArray);
291 else
292 ci->nbrs = NULL;
293 }
294 else
295 {
296 if (ci->nbrsArray)
297 ci->nbrs = (int *)DXGetArrayData(ci->nbrsArray);
298 else
299 ci->nbrs = NULL;
300 }
301
302 if (DXGetError())
303 return ERROR;
304
305 /*
306 * Create the search grid
307 */
308 ci->grid = _dxfMakeSearchGrid(((Interpolator)ci)->max,
309 ((Interpolator)ci)->min, ci->nCubes, 3);
310
311 if (! ci->grid)
312 return ERROR;
313
314 for (i = 0; i < ci->nCubes; i++)
315 {
316 float *points[8], pbuf[24];
317 register int j;
318 int cBuf[8], *c;
319
320 c = (int *)DXGetArrayEntry(ci->cHandle, i, (Pointer)cBuf);
321
322 for (j = 0; j < 8; j++)
323 points[j] = (float *)
324 DXGetArrayEntry(ci->pHandle, c[j], (Pointer)(pbuf+j*3));
325
326 if (! _dxfAddItemToSearchGrid(ci->grid, points, 8, i))
327 break;
328 }
329
330 ci->gridFlag = 1;
331
332 ci->visited = (int *)DXAllocate(ci->nCubes * sizeof(int));
333 if (! ci->visited)
334 return ERROR;
335
336 memset(ci->visited, 0, ci->nCubes*sizeof(int));
337
338 return OK;
339 }
340
341 Error
_dxfCubesIIInterpolator_Delete(CubesIIInterpolator ci)342 _dxfCubesIIInterpolator_Delete(CubesIIInterpolator ci)
343 {
344 _dxfCleanup(ci);
345 _dxfFieldInterpolator_Delete((FieldInterpolator) ci);
346
347 return OK;
348 }
349
350 int
_dxfCubesIIInterpolator_PrimitiveInterpolate(CubesIIInterpolator ci,int * n,float ** points,Pointer * values,int fuzzFlag)351 _dxfCubesIIInterpolator_PrimitiveInterpolate(CubesIIInterpolator ci,
352 int *n, float **points, Pointer *values, int fuzzFlag)
353 {
354 int i;
355 int found;
356 int cubeI;
357 Point *p;
358 Pointer v;
359 float fuzz;
360 float bary[4];
361 CubesIIArgs args;
362 int *visited;
363 int dep;
364 float xMin, xMax;
365 float yMin, yMax;
366 float zMin, zMax;
367 int itemSize;
368 ubyte *dbuf = NULL;
369 InvalidComponentHandle icH = ((FieldInterpolator)ci)->invCon;
370 Matrix *xform;
371
372
373 if (! ci->fieldInterpolator.initialized)
374 {
375 if (! _dxfInitialize(ci))
376 {
377 _dxfCleanup(ci);
378 return 0;
379 }
380
381 ci->fieldInterpolator.initialized = 1;
382 }
383
384 xMin = ((Interpolator)ci)->min[0];
385 yMin = ((Interpolator)ci)->min[1];
386 zMin = ((Interpolator)ci)->min[2];
387
388 xMax = ((Interpolator)ci)->max[0];
389 yMax = ((Interpolator)ci)->max[1];
390 zMax = ((Interpolator)ci)->max[2];
391
392 dep = ci->fieldInterpolator.data_dependency;
393
394 itemSize = ci->nElements*DXTypeSize(((Interpolator)ci)->type);
395
396 args.pHandle = ci->pHandle;
397 args.cHandle = ci->cHandle;
398 args.nbrs = ci->nbrs;
399 args.fuzz = fuzz = ci->fieldInterpolator.fuzz;
400 args.visited = ci->visited;
401 args.bary = bary;
402 args.limit = SEARCH_LIMIT;
403
404 args.gridCounts = ci->gridCounts;
405 args.Cyz = ci->Cyz;
406 args.Cz = ci->Cz;
407
408 args.ci = ci;
409
410 visited = ci->visited;
411
412 dbuf = (ubyte *)DXAllocate(4*itemSize);
413 if (! dbuf)
414 return 0;
415
416 if (((FieldInterpolator)ci)->xflag)
417 xform = &((FieldInterpolator)ci)->xform;
418 else
419 xform = NULL;
420
421 p = (Point *)*points;
422 v = (Pointer)*values;
423
424
425 while (*n)
426 {
427 Point xpt;
428 Point *pPtr;
429
430 if (xform)
431 {
432 xpt.x = p->x*xform->A[0][0] +
433 p->y*xform->A[1][0] +
434 p->z*xform->A[2][0] +
435 xform->b[0];
436
437 xpt.y = p->x*xform->A[0][1] +
438 p->y*xform->A[1][1] +
439 p->z*xform->A[2][1] +
440 xform->b[1];
441
442 xpt.z = p->x*xform->A[0][2] +
443 p->y*xform->A[1][2] +
444 p->z*xform->A[2][2] +
445 xform->b[2];
446
447 pPtr = &xpt;
448 }
449 else
450 pPtr = p;
451
452
453 if (fuzzFlag == FUZZ_ON)
454 {
455 if ((pPtr->x > xMax+fuzz ||
456 pPtr->x < xMin-fuzz) ||
457 (pPtr->y > yMax+fuzz ||
458 pPtr->y < yMin-fuzz) ||
459 (pPtr->z > zMax+fuzz ||
460 pPtr->z < zMin-fuzz))
461 break;
462 }
463 else
464 {
465 if ((pPtr->x > xMax ||
466 pPtr->x < xMin) ||
467 (pPtr->y > yMax ||
468 pPtr->y < yMin) ||
469 (pPtr->z > zMax ||
470 pPtr->z < zMin))
471 break;
472 }
473
474 found = OUTSIDE;
475
476 args.fuzzFlag = fuzzFlag;
477 args.point = pPtr;
478 args.vCount = ci->vCount++;
479
480 /*
481 * Check the hint first, if one exists...
482 */
483
484 cubeI = ci->hint;
485
486 /*
487 fuzz = fuzzFlag * fuzz;
488 */
489
490 if (cubeI >= 0)
491 {
492 found = CubesIISearch(&args, cubeI);
493 }
494
495 if (found == OUTSIDE)
496 {
497 /*
498 * If the sample was not found in the hint, look in the bin
499 *
500 * For each cube in the bin...
501 */
502 if (ci->gridFlag)
503 {
504 _dxfInitGetSearchGrid(ci->grid, (float *)p);
505 while((cubeI = _dxfGetNextSearchGrid(ci->grid)) != -1
506 && found==OUTSIDE)
507 {
508 if (visited[cubeI] != args.vCount)
509 {
510 found = CubesIISearch(&args, cubeI);
511 }
512 }
513 }
514 else
515 {
516 for (cubeI = 0; cubeI < ci->nCubes && found==OUTSIDE; cubeI++)
517 if (visited[cubeI] != args.vCount)
518 {
519 found = CubesIISearch(&args, cubeI);
520 }
521 }
522
523 }
524
525 if (found == OUTSIDE || (icH && DXIsElementInvalid(icH, found)))
526 break;
527
528 ci->hint = found;
529
530 #define INTERPOLATE(type, round) \
531 { \
532 type *d0, *d1, *d2, *d3, *r; \
533 \
534 r = (type *)v; \
535 \
536 /* \
537 * Get pointers to first data elements for each corner \
538 */ \
539 d0 = (type *)DXGetArrayEntry(ci->dHandle, \
540 args.tetra[0], (Pointer)(dbuf)); \
541 \
542 d1 = (type *)DXGetArrayEntry(ci->dHandle, \
543 args.tetra[1], (Pointer)(dbuf+itemSize)); \
544 \
545 d2 = (type *)DXGetArrayEntry(ci->dHandle, \
546 args.tetra[2], (Pointer)(dbuf+2*itemSize)); \
547 \
548 d3 = (type *)DXGetArrayEntry(ci->dHandle, \
549 args.tetra[3], (Pointer)(dbuf+3*itemSize)); \
550 \
551 for (i = 0; i < ci->nElements; i++) \
552 { \
553 *r++ = round + (bary[0] * (*d0++)) + (bary[1] * (*d1++)) + \
554 (bary[2] * (*d2++)) + (bary[3] * (*d3++)); \
555 } \
556 \
557 v = (Pointer)r; \
558 }
559
560 if (((FieldInterpolator)ci)->cstData)
561 {
562 memcpy(v, ((FieldInterpolator)ci)->cstData, itemSize);
563 v = (Pointer)(((char *)v) + itemSize);
564 }
565 else if (dep == DATA_POSITIONS_DEPENDENT)
566 {
567 Type dataType;
568
569 if ((dataType = ((Interpolator)ci)->type) == TYPE_FLOAT)
570 {
571 INTERPOLATE(float, 0.0);
572 }
573 else if (dataType == TYPE_DOUBLE)
574 {
575 INTERPOLATE(double, 0.0);
576 }
577 else if (dataType == TYPE_INT)
578 {
579 INTERPOLATE(int, 0.5);
580 }
581 else if (dataType == TYPE_SHORT)
582 {
583 INTERPOLATE(short, 0.5);
584 }
585 else if (dataType == TYPE_USHORT)
586 {
587 INTERPOLATE(ushort, 0.5);
588 }
589 else if (dataType == TYPE_UINT)
590 {
591 INTERPOLATE(uint, 0.5);
592 }
593 else if (dataType == TYPE_BYTE)
594 {
595 INTERPOLATE(byte, 0.5);
596 }
597 else if (dataType == TYPE_UBYTE)
598 {
599 INTERPOLATE(ubyte, 0.5);
600 }
601 else
602 {
603 INTERPOLATE(unsigned char, 0.5);
604 }
605 }
606 else
607 {
608 memcpy(v, DXGetArrayEntry(ci->dHandle,
609 found, (Pointer)dbuf), itemSize);
610 v = (Pointer)(((char *)v) + itemSize);
611 }
612
613 (*n) --;
614 p++;
615
616 /*
617 * Only use fuzz the first time
618 */
619 fuzzFlag = FUZZ_OFF;
620 }
621
622 if (dbuf)
623 DXFree((Pointer)dbuf);
624
625 *points = (float *)p;
626 *values = (Pointer)v;
627
628 return OK;
629 }
630
631 static void
_dxfCleanup(CubesIIInterpolator ci)632 _dxfCleanup(CubesIIInterpolator ci)
633 {
634 if (ci->fieldInterpolator.localized)
635 {
636 if (ci->nbrs)
637 DXFreeArrayDataLocal(ci->nbrsArray, (Pointer)ci->nbrs);
638 }
639
640 if (ci->dHandle)
641 {
642 DXFreeArrayHandle(ci->dHandle);
643 ci->dHandle = NULL;
644 }
645
646 if (ci->pHandle)
647 {
648 DXFreeArrayHandle(ci->pHandle);
649 ci->pHandle = NULL;
650 }
651
652 if (ci->cHandle)
653 {
654 DXFreeArrayHandle(ci->cHandle);
655 ci->cHandle = NULL;
656 }
657
658 DXFree((Pointer)ci->visited);
659
660 _dxfFreeSearchGrid(ci->grid);
661 }
662
663 Object
_dxfCubesIIInterpolator_Copy(CubesIIInterpolator old,enum _dxd_copy copy)664 _dxfCubesIIInterpolator_Copy(CubesIIInterpolator old, enum _dxd_copy copy)
665 {
666 CubesIIInterpolator new;
667
668 new = (CubesIIInterpolator)
669 _dxf_NewObject((struct object_class *)&_dxdcubesiiinterpolator_class);
670
671 if (!(_dxf_CopyCubesIIInterpolator(new, old, copy)))
672 {
673 DXDelete((Object)new);
674 return NULL;
675 }
676 else
677 return (Object)new;
678 }
679
680 CubesIIInterpolator
_dxf_CopyCubesIIInterpolator(CubesIIInterpolator new,CubesIIInterpolator old,enum _dxd_copy copy)681 _dxf_CopyCubesIIInterpolator(CubesIIInterpolator new,
682 CubesIIInterpolator old, enum _dxd_copy copy)
683 {
684
685 if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
686 (FieldInterpolator)old, copy))
687 return NULL;
688
689 new->nPoints = old->nPoints;
690 new->nCubes = old->nCubes;
691 new->nElements = old->nElements;
692 new->hint = old->hint;
693 new->visited = NULL;
694
695 if (new->fieldInterpolator.initialized)
696 {
697 new->pointsArray = old->pointsArray;
698 new->cubesArray = old->cubesArray;
699 new->dataArray = old->dataArray;
700 new->nbrsArray = old->nbrsArray;
701
702 new->cHandle = DXCreateArrayHandle(old->cubesArray);
703 new->pHandle = DXCreateArrayHandle(old->pointsArray);
704 new->dHandle = DXCreateArrayHandle(old->dataArray);
705 if (! new->cHandle || ! new->pHandle || ! new->dHandle)
706 return NULL;
707
708 if (new->fieldInterpolator.localized)
709 {
710 if (old->nbrsArray)
711 new->nbrs = (int *)DXGetArrayDataLocal(new->nbrsArray);
712 }
713 else
714 {
715 if (old->nbrsArray)
716 new->nbrs = (int *)DXGetArrayData(new->nbrsArray);
717 }
718
719 if (old->gridFlag)
720 {
721 new->gridFlag = 1;
722 new->grid = _dxfCopySearchGrid(old->grid);
723 if (! new->grid)
724 return NULL;
725 }
726 else
727 new->gridFlag = 0;
728 }
729 else
730 {
731 new->pointsArray = NULL;
732 new->dataArray = NULL;
733 new->cubesArray = NULL;
734 new->nbrsArray = NULL;
735
736 new->cHandle = NULL;
737 new->dHandle = NULL;
738 new->pHandle = NULL;
739
740 new->nbrs = NULL;
741 new->visited = NULL;
742
743 new->gridFlag = 0;
744 }
745
746 if (! old->nbrsArray)
747 {
748 memcpy(new->gridCounts, old->gridCounts, sizeof(old->gridCounts));
749 new->Cyz = old->Cyz;
750 new->Cz = old->Cz;
751 }
752
753 new->vCount = 1;
754
755 new->visited = (int *)DXAllocate(new->nCubes*sizeof(int));
756 if (! new->visited)
757 return NULL;
758
759 memset(new->visited, 0, new->nCubes*sizeof(int));
760
761 if (DXGetError())
762 return NULL;
763
764 return new;
765 }
766
767 Interpolator
_dxfCubesIIInterpolator_LocalizeInterpolator(CubesIIInterpolator ci)768 _dxfCubesIIInterpolator_LocalizeInterpolator(CubesIIInterpolator ci)
769 {
770 ci->fieldInterpolator.localized = 1;
771
772 if (ci->fieldInterpolator.initialized)
773 {
774 if (ci->nbrs)
775 ci->nbrs = (int *)DXGetArrayDataLocal(ci->nbrsArray);
776 }
777
778 if (DXGetError())
779 return NULL;
780 else
781 return (Interpolator)ci;
782 }
783
784 #define TET0 0x01
785 #define TET1 0x02
786 #define TET2 0x04
787 #define TET3 0x08
788 #define TET4 0x10
789 #define TET5 0x20
790
791 static int
IsInCube(CubesIIArgs * args,int index)792 IsInCube(CubesIIArgs *args, int index)
793 {
794 int *cube, cubeBuf[8];
795 char been_here;
796 int *tet;
797
798 #ifdef TEST
799 nCubeSearches ++;
800 #endif
801
802 cube = DXGetArrayEntry(args->cHandle, index, (Pointer)cubeBuf);
803 tet = args->tetra;
804
805 been_here = 0;
806
807 goto tetra1;
808
809 tetra0:
810 if (been_here & TET0) goto done;
811
812 been_here |= TET0;
813
814 tet[0] = cube[0]; tet[1] = cube[3]; tet[2] = cube[1]; tet[3] = cube[5];
815 switch(IsInTetra(args)) {
816 case INSIDE: return INSIDE;
817 case 0: return 5;
818 case 1: return 2;
819 case 2: goto tetra3;
820 case 3: return 0;
821 }
822
823 tetra1:
824 if (been_here & TET1) goto done;
825
826 been_here |= TET1;
827
828 tet[0] = cube[5]; tet[1] = cube[2]; tet[2] = cube[3]; tet[3] = cube[7];
829 switch(IsInTetra(args)) {
830 case INSIDE: return INSIDE;
831 case 0: return 3;
832 case 1: return 5;
833 case 2: goto tetra2;
834 case 3: goto tetra3;
835 }
836
837 tetra2:
838 if (been_here & TET2) goto done;
839
840 been_here |= TET2;
841
842 tet[0] = cube[5]; tet[1] = cube[4]; tet[2] = cube[2]; tet[3] = cube[7];
843 switch(IsInTetra(args)) {
844 case INSIDE: return INSIDE;
845 case 0: goto tetra5;
846 case 1: goto tetra1;
847 case 2: return 1;
848 case 3: goto tetra4;
849 }
850
851 tetra3:
852 if (been_here & TET3) goto done;
853
854 been_here |= TET3;
855
856 tet[0] = cube[5]; tet[1] = cube[3]; tet[2] = cube[2]; tet[3] = cube[0];
857 switch(IsInTetra(args)) {
858 case INSIDE: return INSIDE;
859 case 0: return 0;
860 case 1: goto tetra4;
861 case 2: goto tetra0;
862 case 3: goto tetra1;
863 }
864
865 tetra4:
866 if (been_here & TET4) goto done;
867
868 been_here |= TET4;
869
870 tet[0] = cube[5]; tet[1] = cube[2]; tet[2] = cube[4]; tet[3] = cube[0];
871 switch(IsInTetra(args)) {
872 case INSIDE: return INSIDE;
873 case 0: return 4;
874 case 1: return 2;
875 case 2: goto tetra3;
876 case 3: goto tetra2;
877 }
878
879 tetra5:
880 if (been_here & TET5) goto done;
881
882 been_here |= TET5;
883
884 tet[0] = cube[4]; tet[1] = cube[2]; tet[2] = cube[7]; tet[3] = cube[6];
885 switch(IsInTetra(args)) {
886 case INSIDE: return INSIDE;
887 case 0: return 3;
888 case 1: return 1;
889 case 2: return 4;
890 case 3: goto tetra2;
891 }
892
893 return OUTSIDE;
894
895 done:
896 return OUTSIDE;
897 }
898
899 static int
IsInTetra(CubesIIArgs * args)900 IsInTetra(CubesIIArgs *args)
901 {
902 Point *pt;
903 Vector *vert, vbuf;
904 float *bary;
905 float v0, v1, v2, v3, v;
906 float xt, x0, x1, x2, x3;
907 float yt, y0, y1, y2, y3;
908 float zt, z0, z1, z2, z3;
909 float x01, x02, x03, x0t, xt1, xt2, xt3;
910 float y01, y02, y03, y0t, yt1, yt2, yt3;
911 float z01, z02, z03, z0t, zt1, zt2, zt3;
912 float yz012, yz013, yz023, yzt12, yzt13;
913 float yzt23, yz0t2, yz0t3, yz01t, yz02t;
914 int p, q, r, s;
915 float opp, fuzz;
916 int k;
917
918 pt = args->point;
919 bary = args->bary;
920 fuzz = args->fuzz;
921
922 p = args->tetra[0];
923 q = args->tetra[1];
924 r = args->tetra[2];
925 s = args->tetra[3];
926
927 xt = pt->x; yt = pt->y; zt = pt->z;
928
929 vert = (Vector *)DXGetArrayEntry(args->pHandle, p, (Pointer)&vbuf);
930 x0 = vert->x; y0 = vert->y; z0 = vert->z;
931
932 vert = (Vector *)DXGetArrayEntry(args->pHandle, q, (Pointer)&vbuf);
933 x1 = vert->x; y1 = vert->y; z1 = vert->z;
934
935 vert = (Vector *)DXGetArrayEntry(args->pHandle, r, (Pointer)&vbuf);
936 x2 = vert->x; y2 = vert->y; z2 = vert->z;
937
938 vert = (Vector *)DXGetArrayEntry(args->pHandle, s, (Pointer)&vbuf);
939 x3 = vert->x; y3 = vert->y; z3 = vert->z;
940
941 x01 = x1 - x0; y01 = y1 - y0; z01 = z1 - z0;
942 x02 = x2 - x0; y02 = y2 - y0; z02 = z2 - z0;
943 x03 = x3 - x0; y03 = y3 - y0; z03 = z3 - z0;
944 x0t = xt - x0; y0t = yt - y0; z0t = zt - z0;
945 xt1 = x1 - xt; yt1 = y1 - yt; zt1 = z1 - zt;
946 xt2 = x2 - xt; yt2 = y2 - yt; zt2 = z2 - zt;
947 xt3 = x3 - xt; yt3 = y3 - yt; zt3 = z3 - zt;
948
949 yz012 = y01 * z02 - y02 * z01;
950 yz013 = y01 * z03 - y03 * z01;
951 yz023 = y02 * z03 - y03 * z02;
952 yzt12 = yt1 * zt2 - yt2 * zt1;
953 yzt13 = yt1 * zt3 - yt3 * zt1;
954 yzt23 = yt2 * zt3 - yt3 * zt2;
955 yz0t2 = y0t * z02 - y02 * z0t;
956 yz0t3 = y0t * z03 - y03 * z0t;
957 yz01t = y01 * z0t - y0t * z01;
958 yz02t = y02 * z0t - y0t * z02;
959
960 v0 = xt1 * yzt23 - xt2 * yzt13 + xt3 * yzt12;
961 v1 = x0t * yz023 - x02 * yz0t3 + x03 * yz0t2;
962 v2 = x01 * yz0t3 - x0t * yz013 + x03 * yz01t;
963 v3 = x01 * yz02t - x02 * yz01t + x0t * yz012;
964
965 v = x01 * yz023 - x02 * yz013 + x03 * yz012;
966
967 fuzz *= v;
968 if (fuzz < 0.0) fuzz = -fuzz;
969
970 if (v == 0.0)
971 {
972 bary[0] = 0.0;
973 bary[1] = 0.0;
974 bary[2] = 0.0;
975 bary[3] = 0.0;
976 return OUTSIDE;
977 }
978
979 if (((v0 >= -fuzz) && (v1 >= -fuzz) && (v2 >= -fuzz) && (v3 >= -fuzz))
980 || ((v0 <= fuzz) && (v1 <= fuzz) && (v2 <= fuzz) && (v3 <= fuzz)))
981 {
982 v = 1.0 / v;
983 bary[0] = v0 * v;
984 bary[1] = v1 * v;
985 bary[2] = v2 * v;
986 bary[3] = 1.0 - (bary[0] + bary[1] + bary[2]);
987 return INSIDE;
988 }
989
990 if (v > 0)
991 {
992 k = 0; opp = v0;
993 if (v1 < opp) {opp = v1; k = 1;}
994 if (v2 < opp) {opp = v2; k = 2;}
995 if (v3 < opp) { k = 3;}
996 }
997 else
998 {
999 k = 0; opp = v0;
1000 if (v1 > opp) {opp = v1; k = 1;}
1001 if (v2 > opp) {opp = v2; k = 2;}
1002 if (v3 > opp) { k = 3;}
1003 }
1004
1005 return k;
1006 }
1007
1008 static int
CubesIISearch(CubesIIArgs * args,int index)1009 CubesIISearch(CubesIIArgs *args, int index)
1010 {
1011 int *nbrs;
1012 int *visited, vCount;
1013 int limit;
1014 int k, new;
1015 int *vPtr;
1016
1017 visited = args->visited;
1018 vCount = args->vCount;
1019 limit = args->limit;
1020 nbrs = args->nbrs;
1021
1022 for (k = 0; index != OUTSIDE && k < limit; k++)
1023 {
1024 if (visited)
1025 {
1026 vPtr = visited + index;
1027
1028 if (*vPtr == vCount)
1029 break;
1030 else
1031 *vPtr = vCount;
1032 }
1033
1034 new = IsInCube(args, index);
1035
1036 /*
1037 * Results of the IsInCube routine are either INSIDE, OUTSIDE,
1038 * or the exit side for the neighbors walk.
1039 */
1040 switch(new)
1041 {
1042 case INSIDE: return index;
1043 case OUTSIDE: index = OUTSIDE; break;
1044 default: if (nbrs)
1045 index = nbrs[(index*6) + new];
1046 else
1047 index = GetRegularNeighbor(args, index, new);
1048 }
1049 }
1050
1051 return OUTSIDE;
1052 }
1053
1054 static int
GetRegularNeighbor(CubesIIArgs * args,int index,int exit)1055 GetRegularNeighbor(CubesIIArgs *args, int index, int exit)
1056 {
1057 int i;
1058
1059 switch(exit) {
1060 case 0:
1061 i = index / args->Cyz;
1062 if (i > 0)
1063 return index - args->Cyz;
1064 else
1065 return -1;
1066
1067 case 1:
1068 i = index / args->Cyz;
1069 if (i < (args->gridCounts[0] - 1))
1070 return index + args->Cyz;
1071 else
1072 return -1;
1073
1074 case 2:
1075 i = (index % args->Cyz) / args->Cz;
1076 if (i > 0)
1077 return index - args->Cz;
1078 else
1079 return -1;
1080
1081 case 3:
1082 i = (index % args->Cyz) / args->Cz;
1083 if (i < (args->gridCounts[1] - 1))
1084 return index + args->Cz;
1085 else
1086 return -1;
1087
1088 case 4:
1089 i = index % args->gridCounts[2];
1090 if (i > 0)
1091 return index - 1;
1092 else
1093 return -1;
1094
1095 case 5:
1096 i = index % args->gridCounts[2];
1097 if (i < (args->gridCounts[2] - 1))
1098 return index + 1;
1099 else
1100 return -1;
1101
1102 default:
1103 DXMessage("illegal exit");
1104 return -1;
1105 }
1106 }
1107