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