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 <string.h>
15 #include <dx/dx.h>
16 #include "math.h"
17 #include "cubesRRClass.h"
18
19 static Error _dxfInitialize(CubesRRInterpolator);
20 static Error _dxfInitializeTask(Pointer);
21 static void _dxfCleanup(CubesRRInterpolator);
22
23 #define DIAGNOSTIC(str) \
24 DXMessage("regular cubes interpolator failure: %s", (str))
25
26 int
_dxfRecognizeCubesRR(Field field)27 _dxfRecognizeCubesRR(Field field)
28 {
29 Array array;
30 int nDim;
31 Type type;
32 Category cat;
33 int rank, shape[32];
34
35 CHECK(field, CLASS_FIELD);
36
37 ELT_TYPECHECK(field, "cubes");
38
39 array = (Array)DXGetComponentValue(field, "positions");
40 if (!array)
41 return ERROR;
42 else
43 {
44 if (!DXQueryGridPositions(array, &nDim, NULL, NULL, NULL))
45 return ERROR;
46
47 if (! DXGetArrayInfo(array, NULL, &type, &cat, &rank, shape))
48 return ERROR;
49
50 if (cat != CATEGORY_REAL)
51 return ERROR;
52
53 if (nDim != 3)
54 return ERROR;
55 }
56
57 array = (Array)DXGetComponentValue(field, "connections");
58 if (!array)
59 return ERROR;
60 else
61 {
62 if (!DXQueryGridConnections(array, &nDim, NULL))
63 return ERROR;
64
65 if (nDim != 3)
66 return ERROR;
67 }
68
69 return OK;
70 }
71
72 CubesRRInterpolator
_dxfNewCubesRRInterpolator(Field field,enum interp_init initType,double fuzz,Matrix * m)73 _dxfNewCubesRRInterpolator(Field field,
74 enum interp_init initType, double fuzz, Matrix *m)
75 {
76 return (CubesRRInterpolator)_dxf_NewCubesRRInterpolator(field,
77 initType, fuzz, m, &_dxdcubesrrinterpolator_class);
78 }
79
80
81 CubesRRInterpolator
_dxf_NewCubesRRInterpolator(Field field,enum interp_init initType,float fuzz,Matrix * m,struct cubesrrinterpolator_class * class)82 _dxf_NewCubesRRInterpolator(Field field,
83 enum interp_init initType, float fuzz, Matrix *m,
84 struct cubesrrinterpolator_class *class)
85 {
86 CubesRRInterpolator ci;
87
88 ci = (CubesRRInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
89 (struct fieldinterpolator_class *)class);
90 if (! ci)
91 return NULL;
92
93 ci->pointsArray = NULL;
94 ci->dataArray = NULL;
95 ci->data = NULL;
96
97 /*
98 * The grid should be regular
99 */
100 {
101 float localOrigin[3], displacement[3];
102 Array connections, positions;
103 Matrix mIn;
104
105 positions = (Array)DXGetComponentValue(field, "positions");
106 connections = (Array)DXGetComponentValue(field, "connections");
107
108 if (!positions || !connections)
109 {
110 DXDelete((Object)(ci));
111 return NULL;
112 }
113
114 DXQueryGridPositions(positions, NULL, NULL, localOrigin, (float *)mIn.A);
115 DXGetMeshOffsets((MeshArray)connections, ci->meshOffsets);
116
117 displacement[0] = ci->meshOffsets[0]*mIn.A[0][0] +
118 ci->meshOffsets[1]*mIn.A[1][0] +
119 ci->meshOffsets[2]*mIn.A[2][0];
120 displacement[1] = ci->meshOffsets[0]*mIn.A[0][1] +
121 ci->meshOffsets[1]*mIn.A[1][1] +
122 ci->meshOffsets[2]*mIn.A[2][1];
123 displacement[2] = ci->meshOffsets[0]*mIn.A[0][2] +
124 ci->meshOffsets[1]*mIn.A[1][2] +
125 ci->meshOffsets[2]*mIn.A[2][2];
126
127 mIn.b[0] = localOrigin[0] - displacement[0];
128 mIn.b[1] = localOrigin[1] - displacement[1];
129 mIn.b[2] = localOrigin[2] - displacement[2];
130
131 ((Interpolator)ci)->matrix = DXInvert(mIn);
132 }
133
134 if (initType == INTERP_INIT_PARALLEL)
135 {
136 if (! DXAddTask(_dxfInitializeTask, (Pointer)&ci, sizeof(ci), 1.0))
137 {
138 DXDelete((Object)ci);
139 return NULL;
140 }
141 }
142 else if (initType == INTERP_INIT_IMMEDIATE)
143 {
144 if (! _dxfInitialize(ci))
145 {
146 DXDelete((Object)ci);
147 return NULL;
148 }
149 }
150
151 return ci;
152 }
153
154 static Error
_dxfInitializeTask(Pointer p)155 _dxfInitializeTask(Pointer p)
156 {
157 return _dxfInitialize(*(CubesRRInterpolator *) p);
158 }
159
160 static Error
_dxfInitialize(CubesRRInterpolator ci)161 _dxfInitialize(CubesRRInterpolator ci)
162 {
163 Field field;
164 int nDim;
165 int i;
166 Type dataType;
167 Category dataCategory;
168
169 field = (Field)((Interpolator)ci)->dataObject;
170
171 ci->fieldInterpolator.initialized = 1;
172
173 /*
174 * De-reference data
175 */
176 ci->dataArray = (Array)DXGetComponentValue(field, "data");
177 if (!ci->dataArray)
178 {
179 DXSetError(ERROR_MISSING_DATA, "#10240", "data");
180 return ERROR;
181 }
182 DXReference((Object)ci->dataArray);
183
184 DXGetArrayInfo(ci->dataArray, NULL, &((Interpolator)ci)->type,
185 &((Interpolator)ci)->category, &((Interpolator)ci)->rank,
186 ((Interpolator)ci)->shape);
187
188 ci->data = DXCreateArrayHandle(ci->dataArray);
189 if (! ci->data)
190 return ERROR;
191
192 /*
193 * Get the grid.
194 */
195 ci->pointsArray = (Array)DXGetComponentValue(field, "positions");
196 if (!ci->pointsArray)
197 {
198 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
199 return ERROR;
200 }
201 DXReference((Object)ci->pointsArray);
202
203 /*
204 * get info about data values
205 */
206 DXGetArrayInfo(ci->dataArray, NULL, &dataType, &dataCategory, NULL, NULL);
207
208 /*
209 * Don't worry about maintaining shape of input; just determine how
210 * many values to interpolate.
211 */
212 ci->nElements = DXGetItemSize(ci->dataArray) / DXTypeSize(dataType);
213
214 DXQueryGridPositions(ci->pointsArray, &nDim, ci->counts, NULL, NULL);
215
216 /*
217 * Figure out how big to increment pointers along each dimension.
218 * Use either number of positions or number of primitives along the
219 * axis, depending on whether the data is positions or connections
220 * dependant.
221 */
222
223 if (ci->fieldInterpolator.data_dependency == DATA_POSITIONS_DEPENDENT)
224 {
225 ci->size[2] = 1;
226 for (i = 1; i >= 0; i--)
227 ci->size[i] = ci->counts[i+1] * ci->size[i+1];
228 }
229 else
230 {
231 ci->size[2] = 1;
232 for (i = 1; i >= 0; i--)
233 ci->size[i] = (ci->counts[i+1] - 1) * ci->size[i+1];
234 }
235
236 ci->eltStrides[2] = 1;
237 ci->eltStrides[1] = ci->counts[2] - 1;
238 ci->eltStrides[0] = (ci->counts[1] - 1) * ci->eltStrides[1];
239
240 return OK;
241 }
242
243 Error
_dxfCubesRRInterpolator_Delete(CubesRRInterpolator ci)244 _dxfCubesRRInterpolator_Delete(CubesRRInterpolator ci)
245 {
246 _dxfCleanup(ci);
247 return _dxfFieldInterpolator_Delete((FieldInterpolator) ci);
248 }
249
250 static void
_dxfCleanup(CubesRRInterpolator ci)251 _dxfCleanup(CubesRRInterpolator ci)
252 {
253 if (ci->data)
254 {
255 DXFreeArrayHandle(ci->data);
256 ci->data = NULL;
257 }
258
259 if (ci->pointsArray)
260 {
261 DXDelete((Object)ci->pointsArray);
262 ci->pointsArray = NULL;
263 }
264
265 if (ci->dataArray)
266 {
267 DXDelete((Object)ci->dataArray);
268 ci->dataArray = NULL;
269 }
270 }
271
272 #define INTERPOLATE(type, round) \
273 { \
274 type *ptr000, *ptr001, *ptr010, *ptr011; \
275 type *ptr100, *ptr101, *ptr110, *ptr111; \
276 type *r = (type *)v; \
277 \
278 ptr000 = (type *)DXGetArrayEntry(ci->data, baseIndex, \
279 (Pointer)(dbuf)); \
280 \
281 ptr100 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz0, \
282 (Pointer)(dbuf+itemSize)); \
283 \
284 ptr010 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz1, \
285 (Pointer)(dbuf+2*itemSize)); \
286 \
287 ptr110 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz1+sz0, \
288 (Pointer)(dbuf+3*itemSize)); \
289 \
290 ptr001 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz2, \
291 (Pointer)(dbuf+4*itemSize)); \
292 \
293 ptr101 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz2+sz0, \
294 (Pointer)(dbuf+5*itemSize)); \
295 \
296 ptr011 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz2+sz1, \
297 (Pointer)(dbuf+6*itemSize)); \
298 \
299 ptr111 = (type *)DXGetArrayEntry(ci->data, baseIndex+sz2+sz1+sz0, \
300 (Pointer)(dbuf+7*itemSize)); \
301 \
302 for (i = 0; i < ci->nElements; i++) \
303 *r++ = weight000*(*ptr000++) + weight001*(*ptr001++) + \
304 weight010*(*ptr010++) + weight011*(*ptr011++) + \
305 weight100*(*ptr100++) + weight101*(*ptr101++) + \
306 weight110*(*ptr110++) + weight111*(*ptr111++) + round; \
307 \
308 v = (Pointer)r; \
309 }
310
311 int
_dxfCubesRRInterpolator_PrimitiveInterpolate(CubesRRInterpolator ci,int * n,float ** points,Pointer * values,int fuzzFlag)312 _dxfCubesRRInterpolator_PrimitiveInterpolate(CubesRRInterpolator ci,
313 int *n, float **points, Pointer *values, int fuzzFlag)
314 {
315 Point pt;
316 int ix, iy, iz;
317 float dx, dy, dz;
318 int sz0, sz1, sz2;
319 int ixMax, iyMax, izMax;
320 float weight000, weight001, weight010, weight011;
321 float weight100, weight101, weight110, weight111;
322 float A, B, C, D;
323 int baseIndex, i;
324 Point *p;
325 Pointer v;
326 float *mA, *mb;
327 float fuzz;
328 int dep;
329 int itemSize;
330 int typeSize;
331 int *mo;
332 int esz0, esz1;
333 InvalidComponentHandle icH = ((FieldInterpolator)ci)->invCon;
334 ubyte *dbuf = NULL;
335 Matrix *xform = NULL;
336
337 if (! ci->fieldInterpolator.initialized)
338 {
339 if (! _dxfInitialize(ci))
340 {
341 _dxfCleanup(ci);
342 return 0;
343 }
344
345 ci->fieldInterpolator.initialized = 1;
346 }
347
348 /*
349 * De-reference indexing info from interpolator
350 */
351 sz0 = ci->size[0];
352 sz1 = ci->size[1];
353 sz2 = ci->size[2];
354 mA = (float *)((Interpolator)ci)->matrix.A;
355 mb = (float *)((Interpolator)ci)->matrix.b;
356 mo = ci->meshOffsets;
357 esz0 = ci->eltStrides[0];
358 esz1 = ci->eltStrides[1];
359
360 dep = ci->fieldInterpolator.data_dependency;
361 typeSize = DXTypeSize(((Interpolator)ci)->type);
362 itemSize = ci->nElements * typeSize;
363
364 ixMax = ci->counts[0] - 1;
365 iyMax = ci->counts[1] - 1;
366 izMax = ci->counts[2] - 1;
367
368 /*
369 * De-reference bounding box
370 */
371 /* xMax = ((Interpolator)ci)->max[0]; */
372 /* xMin = ((Interpolator)ci)->min[0]; */
373 /* yMax = ((Interpolator)ci)->max[1]; */
374 /* yMin = ((Interpolator)ci)->min[1]; */
375 /* zMax = ((Interpolator)ci)->max[2]; */
376 /* zMin = ((Interpolator)ci)->min[2]; */
377
378 p = (Point *)*points;
379 v = (Pointer)*values;
380
381 fuzz = ((FieldInterpolator)ci)->fuzz;
382
383 dbuf = (ubyte *)DXAllocate(8*itemSize);
384 if (! dbuf)
385 return 0;
386
387 if (((FieldInterpolator)ci)->xflag)
388 xform = &(((FieldInterpolator)ci)->xform);
389
390 while (*n > 0)
391 {
392
393 if (((FieldInterpolator)ci)->xflag)
394 {
395 float x, y, z;
396
397 x = p->x*xform->A[0][0] +
398 p->y*xform->A[1][0] +
399 p->z*xform->A[2][0] +
400 xform->b[0];
401
402 y = p->x*xform->A[0][1] +
403 p->y*xform->A[1][1] +
404 p->z*xform->A[2][1] +
405 xform->b[1];
406
407 z = p->x*xform->A[0][2] +
408 p->y*xform->A[1][2] +
409 p->z*xform->A[2][2] +
410 xform->b[2];
411
412 /*
413 * This matrix transforms the point in (XYZ) space into a point
414 * in index space.
415 */
416 pt.x = (x*mA[0] + y*mA[3] + z*mA[6] + mb[0]);
417 pt.y = (x*mA[1] + y*mA[4] + z*mA[7] + mb[1]);
418 pt.z = (x*mA[2] + y*mA[5] + z*mA[8] + mb[2]);
419 }
420 else
421 {
422 /*
423 * This matrix transforms the point in (XYZ) space into a point
424 * in index space.
425 */
426 pt.x = (p->x*mA[0] + p->y*mA[3] + p->z*mA[6] + mb[0]);
427 pt.y = (p->x*mA[1] + p->y*mA[4] + p->z*mA[7] + mb[1]);
428 pt.z = (p->x*mA[2] + p->y*mA[5] + p->z*mA[8] + mb[2]);
429 }
430
431 pt.x -= mo[0];
432 pt.y -= mo[1];
433 pt.z -= mo[2];
434
435 if (fuzzFlag == FUZZ_ON)
436 {
437 if ((pt.x < -fuzz ||
438 pt.x > ixMax+fuzz) ||
439 (pt.y < -fuzz ||
440 pt.y > iyMax+fuzz) ||
441 (pt.z < -fuzz ||
442 pt.z > izMax+fuzz)) break;
443 }
444 else
445 {
446 if ((pt.x < 0 ||
447 pt.x > ixMax) ||
448 (pt.y < 0 ||
449 pt.y > iyMax) ||
450 (pt.z < 0 ||
451 pt.z > izMax)) break;
452 }
453
454 if (((FieldInterpolator)ci)->cstData)
455 {
456 memcpy(v, ((FieldInterpolator)ci)->cstData, itemSize);
457 v = (Pointer)(((char *)v) + itemSize);
458 }
459 else if (dep == DATA_POSITIONS_DEPENDENT)
460 {
461 Type dataType;
462
463 ix = pt.x;
464 dx = pt.x - ix;
465
466 if (fuzzFlag != FUZZ_ON && (ix>ixMax || (ix==ixMax && dx>0.0) || ix<0))
467 goto not_found;
468
469 if (ix >= ixMax)
470 {
471 ix = ixMax - 1;
472 dx = 1.0;
473 }
474 else if (ix < 0 || dx < 0.0)
475 {
476 ix = 0;
477 dx = 0.0;
478 }
479
480 iy = pt.y;
481 dy = pt.y - iy;
482
483 if (fuzzFlag != FUZZ_ON && (iy>iyMax || (iy==iyMax && dy>0.0) || iy<0))
484 goto not_found;
485
486 if (iy >= iyMax)
487 {
488 iy = iyMax - 1;
489 dy = 1.0;
490 }
491 else if (iy < 0 || dy < 0.0)
492 {
493 iy = 0;
494 dy = 0.0;
495 }
496
497 iz = pt.z;
498 dz = pt.z - iz;
499
500 if (fuzzFlag != FUZZ_ON && (iz>izMax || (iz==izMax && dz>0.0) || iz<0))
501 goto not_found;
502
503 if (iz >= izMax)
504 {
505 iz = izMax - 1;
506 dz = 1.0;
507 }
508 else if (iz < 0 || dz < 0.0)
509 {
510 iz = 0;
511 dz = 0.0;
512 }
513
514 if (icH)
515 {
516 baseIndex = ix*esz0 + iy*esz1 + iz;
517 if (DXIsElementInvalid(icH, baseIndex))
518 goto not_found;
519 }
520
521 baseIndex = (ix*sz0 + iy*sz1 + iz*sz2);
522
523 A = dx * dy;
524 B = dx * dz;
525 C = dy * dz;
526 D = dz * A;
527
528 weight111 = D;
529 weight011 = C - D;
530 weight101 = B - D;
531 weight001 = dz - B - C + D;
532 weight110 = A - D;
533 weight010 = dy - C - A + D;
534 weight100 = dx - B - A + D;
535 weight000 = 1.0 - dz - dy - dx + A + B + C - D;
536
537 if ((dataType = ((Interpolator)ci)->type) == TYPE_FLOAT)
538 {
539 INTERPOLATE(float, 0.0);
540 }
541 else if (dataType == TYPE_DOUBLE)
542 {
543 INTERPOLATE(double, 0.0);
544 }
545 else if (dataType == TYPE_INT)
546 {
547 INTERPOLATE(int, 0.5);
548 }
549 else if (dataType == TYPE_SHORT)
550 {
551 INTERPOLATE(short, 0.5);
552 }
553 else if (dataType == TYPE_USHORT)
554 {
555 INTERPOLATE(ushort, 0.5);
556 }
557 else if (dataType == TYPE_UINT)
558 {
559 INTERPOLATE(uint, 0.5);
560 }
561 else if (dataType == TYPE_BYTE)
562 {
563 INTERPOLATE(byte, 0.5);
564 }
565 else if (dataType == TYPE_UBYTE)
566 {
567 INTERPOLATE(ubyte, 0.5);
568 }
569 else
570 {
571 INTERPOLATE(unsigned char, 0.5);
572 }
573 }
574 else
575 {
576 ix = pt.x;
577 if (ix > ixMax || ix < 0)
578 goto not_found;
579
580 if (ix >= ixMax)
581 ix = ixMax - 1;
582 else if (ix < 0)
583 ix = 0;
584
585 iy = pt.y;
586 if (iy > iyMax || iy < 0)
587 goto not_found;
588
589 if (iy >= iyMax)
590 iy = iyMax - 1;
591 else if (iy < 0)
592 iy = 0;
593
594 iz = pt.z;
595 if (iz > izMax || iz < 0)
596 goto not_found;
597
598 if (iz >= izMax)
599 iz = izMax - 1;
600 else if (iz < 0)
601 iz = 0;
602
603 if (icH)
604 {
605 baseIndex = ix*esz0 + iy*esz1 + iz;
606 if (DXIsElementInvalid(icH, baseIndex))
607 goto not_found;
608 }
609
610 baseIndex = (ix*sz0 + iy*sz1 + iz*sz2);
611
612 memcpy(v, DXGetArrayEntry(ci->data,
613 baseIndex, (Pointer)dbuf), itemSize);
614 v = (Pointer)((char *)v + itemSize);
615 }
616
617 /*
618 * Only use fuzz on first primitive
619 */
620 fuzzFlag = FUZZ_OFF;
621
622 p ++;
623 (*n)--;
624 }
625
626 not_found:
627
628 DXFree((Pointer)dbuf);
629
630 *points = (float *)p;
631 *values = (Pointer)v;
632
633 return OK;
634 }
635
636 Object
_dxfCubesRRInterpolator_Copy(CubesRRInterpolator old,enum _dxd_copy copy)637 _dxfCubesRRInterpolator_Copy(CubesRRInterpolator old, enum _dxd_copy copy)
638 {
639 CubesRRInterpolator new;
640
641 new = (CubesRRInterpolator)
642 _dxf_NewObject((struct object_class *)&_dxdcubesrrinterpolator_class);
643
644 if (!(_dxf_CopyCubesRRInterpolator(new, old, copy)))
645 {
646 DXDelete((Object)new);
647 return NULL;
648 }
649 else
650 return (Object)new;
651 }
652
653 CubesRRInterpolator
_dxf_CopyCubesRRInterpolator(CubesRRInterpolator new,CubesRRInterpolator old,enum _dxd_copy copy)654 _dxf_CopyCubesRRInterpolator(CubesRRInterpolator new,
655 CubesRRInterpolator old, enum _dxd_copy copy)
656 {
657
658 if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
659 (FieldInterpolator)old, copy))
660 return NULL;
661
662 new->nElements = old->nElements;
663
664 memcpy((char *)new->size, (char *)old->size, sizeof(old->size));
665 memcpy((char *)new->counts, (char *)old->counts, sizeof(old->counts));
666 memcpy((char *)new->eltStrides, (char *)old->eltStrides, sizeof(old->eltStrides));
667 memcpy((char *)new->meshOffsets, (char *)old->meshOffsets, sizeof(old->meshOffsets));
668
669 new->fuzz = old->fuzz;
670
671 if (new->fieldInterpolator.initialized)
672 {
673 new->pointsArray = (Array)DXReference((Object)old->pointsArray);
674 new->dataArray = (Array)DXReference((Object)old->dataArray);
675 new->data = DXCreateArrayHandle(old->dataArray);
676 }
677 else
678 {
679 new->pointsArray = NULL;
680 new->dataArray = NULL;
681
682 new->data = NULL;
683 }
684
685 if (DXGetError())
686 return NULL;
687
688 return new;
689 }
690
691 Interpolator
_dxfCubesRRInterpolator_LocalizeInterpolator(CubesRRInterpolator ci)692 _dxfCubesRRInterpolator_LocalizeInterpolator(CubesRRInterpolator ci)
693 {
694 ci->fieldInterpolator.localized = 1;
695 return (Interpolator)ci;
696 }
697