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 "linesRI1DClass.h"
18
19 #define FALSE 0
20 #define TRUE 1
21
22 static Error _dxfInitialize(LinesRI1DInterpolator);
23 static Error _dxfInitializeTask(Pointer);
24 static int linear_coords(float, float, float, LinearCoord *, float);
25 static int _dxfCleanup(LinesRI1DInterpolator);
26
27 int
_dxfRecognizeLinesRI1D(Field field)28 _dxfRecognizeLinesRI1D(Field field)
29 {
30 Array array;
31 int nDim;
32 int rank, shape[32];
33 Type type;
34 Category cat;
35 int nPts;
36
37 CHECK(field, CLASS_FIELD);
38
39 ELT_TYPECHECK(field, "lines");
40
41 if (DXGetComponentValue(field, "invalid connections") ||
42 DXGetComponentValue(field, "invalid positions"))
43 return 0;
44
45 array = (Array)DXGetComponentValue(field, "connections");
46 if (!array)
47 {
48 DXSetError(ERROR_MISSING_DATA, "#10240", "connections");
49 return 0;
50 }
51
52 if (! DXQueryGridConnections(array, &nDim, &nPts))
53 return 0;
54
55 if (nDim != 1)
56 {
57 DXSetError(ERROR_DATA_INVALID, "#11004", "lines", 2);
58 return 0;
59 }
60
61 array = (Array)DXGetComponentValue(field, "positions");
62 if (!array)
63 {
64 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
65 return 0;
66 }
67
68 if (DXQueryGridPositions(array, &rank, NULL, NULL, NULL))
69 {
70 if (rank != 1)
71 {
72 DXSetError(ERROR_DATA_INVALID, "#11003",
73 "lines", 1);
74 return 0;
75 }
76 }
77 else
78 {
79 float *pts = (float *)DXGetArrayData(array);
80 int i, n;
81
82 DXGetArrayInfo(array, &n, &type, &cat, &rank, shape);
83
84 if (rank != 0 && (rank != 1 || shape[0] != 1))
85 {
86 DXSetError(ERROR_DATA_INVALID, "#11003", "lines", 1);
87 return 0;
88 }
89
90 if (nPts > n)
91 {
92 DXSetError(ERROR_DATA_INVALID, "not enough positions");
93 return 0;
94 }
95
96 if (nPts > 2)
97 {
98 for (i = 1; i < nPts; i++)
99 if (pts[i] != pts[i-1])
100 break;
101
102 if (i != nPts)
103 {
104 if (pts[i-1] > pts[i])
105 {
106 for (i = 2; i < nPts; i++)
107 if (pts[i-1] < pts[i])
108 {
109 DXSetError(ERROR_DATA_INVALID, "non-monotonic map");
110 return 0;
111 }
112 }
113 else
114 {
115 for (i = 2; i < nPts; i++)
116 if (pts[i-1] > pts[i])
117 {
118 DXSetError(ERROR_DATA_INVALID, "non-monotonic map");
119 return 0;
120 }
121 }
122 }
123 }
124 }
125
126 return 1;
127 }
128
129 LinesRI1DInterpolator
_dxfNewLinesRI1DInterpolator(Field field,enum interp_init initType,double fuzz,Matrix * m)130 _dxfNewLinesRI1DInterpolator(Field field,
131 enum interp_init initType, double fuzz, Matrix *m)
132 {
133 return (LinesRI1DInterpolator)_dxf_NewLinesRI1DInterpolator(field,
134 initType, fuzz, m, &_dxdlinesri1dinterpolator_class);
135 }
136
137 LinesRI1DInterpolator
_dxf_NewLinesRI1DInterpolator(Field field,enum interp_init initType,float fuzz,Matrix * m,struct linesri1dinterpolator_class * class)138 _dxf_NewLinesRI1DInterpolator(Field field,
139 enum interp_init initType, float fuzz, Matrix *m,
140 struct linesri1dinterpolator_class *class)
141 {
142 LinesRI1DInterpolator li;
143
144 li = (LinesRI1DInterpolator)_dxf_NewFieldInterpolator(field, fuzz, m,
145 (struct fieldinterpolator_class *)class);
146
147 if (! li)
148 return NULL;
149
150 li->pointsArray = NULL;
151 li->dataArray = NULL;
152 li->dHandle = NULL;
153 li->linesArray = NULL;
154
155 li->hint = 0;
156
157 if (initType == INTERP_INIT_PARALLEL)
158 {
159 if (! DXAddTask(_dxfInitializeTask, (Pointer)&li, sizeof(li), 1.0))
160 {
161 DXDelete((Object)li);
162 return NULL;
163 }
164 }
165 else if (initType == INTERP_INIT_IMMEDIATE)
166 {
167 if (! _dxfInitialize(li))
168 {
169 DXDelete((Object)li);
170 return NULL;
171 }
172 }
173
174 return li;
175 }
176
177 static Error
_dxfInitializeTask(Pointer p)178 _dxfInitializeTask(Pointer p)
179 {
180 return _dxfInitialize(*(LinesRI1DInterpolator *)p);
181 }
182
183 static Error
_dxfInitialize(LinesRI1DInterpolator li)184 _dxfInitialize(LinesRI1DInterpolator li)
185 {
186 Field field;
187 int seg, i, i0, i1, dir;
188 float p0, p1;
189 float len;
190
191 li->fieldInterpolator.initialized = 1;
192
193 field = (Field)((Interpolator)li)->dataObject;
194
195 /*
196 * De-reference data
197 */
198 li->dataArray = (Array)DXGetComponentValue(field, "data");
199 if (!li->dataArray)
200 {
201 DXSetError(ERROR_MISSING_DATA, "#10240", "data");
202 return ERROR;
203 }
204 DXReference((Object)li->dataArray);
205
206 DXGetArrayInfo(li->dataArray, NULL, &((Interpolator)li)->type,
207 &((Interpolator)li)->category,
208 &((Interpolator)li)->rank, ((Interpolator)li)->shape);
209
210 li->dHandle = DXCreateArrayHandle(li->dataArray);
211 if (! li->dHandle)
212 return ERROR;
213
214 /*
215 * Get the grid.
216 */
217 li->pointsArray = (Array)DXGetComponentValue(field, "positions");
218 if (!li->pointsArray)
219 {
220 DXSetError(ERROR_MISSING_DATA, "#10240", "positions");
221 return ERROR;
222 }
223 DXReference((Object)li->pointsArray);
224
225 li->pHandle = DXCreateArrayHandle(li->pointsArray);
226 if (! li->pHandle)
227 return ERROR;
228
229 /*
230 * Get the lines.
231 */
232 li->linesArray = (Array)DXGetComponentValue(field, "connections");
233 if (!li->linesArray)
234 {
235 DXSetError(ERROR_MISSING_DATA, "#10240", "connections");
236 return ERROR;
237 }
238
239 if (DXQueryGridConnections(li->linesArray, NULL, &li->nLines))
240 {
241 li->linesArray = NULL;
242
243 /*
244 * The above call returns the number of points in the sequence, not
245 * the number of segments.
246 */
247 li->nLines -= 1;
248 }
249 else
250 {
251 Type type;
252 Category cat;
253 int rank, shape[32];
254
255 if (! DXGetArrayInfo(li->linesArray, &li->nLines,
256 &type, &cat, &rank, shape))
257 return ERROR;
258
259 if (type != TYPE_INT || cat != CATEGORY_REAL ||
260 rank != 1 || shape[0] != 2)
261 return ERROR;
262
263 DXReference((Object)li->linesArray);
264 }
265
266 /*
267 * Don't worry about maintaining shape of input; just determine how
268 * many values to interpolate.
269 */
270 li->nElements = DXGetItemSize(li->dataArray) /
271 DXTypeSize(((Interpolator)li)->type);
272
273 if (li->fieldInterpolator.localized)
274 {
275 if (li->linesArray)
276 li->lines = (int *)DXGetArrayDataLocal(li->linesArray);
277 else
278 li->lines = NULL;
279 }
280 else
281 {
282 if (li->linesArray)
283 li->lines = (int *)DXGetArrayDataLocal(li->linesArray);
284 else
285 li->lines = NULL;
286 }
287
288 if (DXGetError())
289 return ERROR;
290
291 /*
292 * Check to verify that the positions are monotonic. If so,
293 * remember which way they go.
294 */
295 dir = 0;
296 for (seg = 0; seg < li->nLines; seg++)
297 {
298 i = seg << 1;
299
300 if (li->lines)
301 {
302 i0 = li->lines[i];
303 i1 = li->lines[i+1];
304 }
305 else
306 {
307 i0 = seg;
308 i1 = seg + 1;
309 }
310
311 p0 = *(float *)DXGetArrayEntry(li->pHandle, i0, (Pointer)&p0);
312 p1 = *(float *)DXGetArrayEntry(li->pHandle, i1, (Pointer)&p1);
313
314 if (p0 < p1)
315 {
316 if (dir == -1)
317 {
318 DXSetError(ERROR_DATA_INVALID, "#11851");
319 return ERROR;
320 }
321 else
322 dir = 1;
323 }
324 else if (p0 > p1)
325 {
326 if (dir == 1)
327 {
328 DXSetError(ERROR_DATA_INVALID, "#11851");
329 return ERROR;
330 }
331 else
332 dir = -1;
333 }
334 }
335
336 li->direction = dir;
337
338 /*
339 * Figure fuzz as a proportion of average primitive length
340 */
341 len = (((Interpolator)li)->max[0] - ((Interpolator)li)->min[0]);
342 if (li->nLines)
343 len /= li->nLines;
344 else
345 len = 0.0;
346
347 li->fieldInterpolator.fuzz *= len;
348
349 li->fieldInterpolator.initialized = 1;
350
351 return OK;
352 }
353
354 Error
_dxfLinesRI1DInterpolator_Delete(LinesRI1DInterpolator li)355 _dxfLinesRI1DInterpolator_Delete(LinesRI1DInterpolator li)
356 {
357 _dxfCleanup(li);
358 return _dxfFieldInterpolator_Delete((FieldInterpolator) li);
359 }
360
361 int
_dxfLinesRI1DInterpolator_PrimitiveInterpolate(LinesRI1DInterpolator li,int * n,float ** points,Pointer * values,int fuzzFlag)362 _dxfLinesRI1DInterpolator_PrimitiveInterpolate(LinesRI1DInterpolator li,
363 int *n, float **points, Pointer *values, int fuzzFlag)
364 {
365 float p0, p1;
366 LinearCoord linesCoord;
367 int seg;
368 int i;
369 int found;
370 int outside;
371 Pointer v;
372 float *p;
373 int dir;
374 int i0 = 0, i1 = 0;
375 float fuzz;
376 int dep;
377 int itemSize;
378 ubyte *dbuf;
379 InvalidComponentHandle icH = li->fieldInterpolator.invCon;
380 Matrix *xform;
381
382
383 if (! li->fieldInterpolator.initialized)
384 {
385 if (! _dxfInitialize(li))
386 {
387 _dxfCleanup(li);
388 return 0;
389 }
390
391 li->fieldInterpolator.initialized = 1;
392 }
393
394 v = *values;
395 p = (float *)*points;
396
397 fuzz = li->fieldInterpolator.fuzz;
398 dep = li->fieldInterpolator.data_dependency;
399 itemSize = li->nElements*DXTypeSize(((Interpolator)li)->type);
400
401 dbuf = (ubyte *)DXAllocate(2*itemSize);
402 if (! dbuf)
403 return 0;
404
405 if (((FieldInterpolator)li)->xflag)
406 xform = &((FieldInterpolator)li)->xform;
407 else
408 xform = NULL;
409
410
411 /*
412 * For each point in the input, attempt to interpolate the point.
413 * When a point cannot be interpolated, quit.
414 */
415 while(*n != 0)
416 {
417 float xpt;
418 float *pPtr;
419
420 if (xform)
421 {
422 xpt = (*p)*xform->A[0][0] + xform->b[0];
423
424 pPtr = &xpt;
425 }
426 else
427 pPtr = p;
428
429 seg = li->hint;
430
431 /*
432 * If no hint, start in the middle.
433 */
434 if (seg < 0)
435 seg = li->nLines >> 1;
436
437 found = FALSE;
438 outside = FALSE;
439 dir = 1;
440
441 while (!found && !outside)
442 {
443 i = seg << 1;
444
445 if (li->lines)
446 {
447 i0 = li->lines[i];
448 i1 = li->lines[i+1];
449 }
450 else
451 {
452 i0 = seg;
453 i1 = seg + 1;
454 }
455
456 p0 = *(float *)DXGetArrayEntry(li->pHandle, i0, (Pointer)&p0);
457 p1 = *(float *)DXGetArrayEntry(li->pHandle, i1, (Pointer)&p1);
458
459 if (li->direction == -1)
460 {
461 float t = p0;
462 p0 = p1;
463 p1 = t;
464 }
465
466 if (p0 == p1)
467 {
468 if (p0 == *pPtr)
469 {
470 linesCoord.p = 1.0;
471 linesCoord.q = 0.0;
472 found = 1;
473 break;
474 }
475 else if (*pPtr > p0)
476 dir = li->direction;
477 else
478 dir = -li->direction;
479 }
480 else
481 {
482 /*
483 * linear_coords will return 0 for inside (inclusive),
484 * +-1 for outside the fuzz zone, +-2 for inside the fuzz
485 * zone.
486 */
487 if (! (dir=linear_coords(*pPtr, p0, p1, &linesCoord, fuzz)))
488 {
489 found = 1;
490 break;
491 }
492
493 if (li->direction == -1)
494 dir = -dir;
495
496 /*
497 * If we are inside the fuzz zone and in the corresponding
498 * extreme of the set of segments, we are done.
499 */
500 if (fuzzFlag == FUZZ_ON &&
501 ((dir == -2 && seg == 0) ||
502 (dir == 2 && seg == (li->nLines-1))))
503 {
504 found = 1;
505 break;
506 }
507 }
508
509 if (dir < 0)
510 seg --;
511 else
512 seg ++;
513
514 if (seg < 0 || seg >= li->nLines) {
515 if (dir == 2 || dir == -2)
516 {
517 seg -= dir;
518 outside = FALSE;
519 found = 1;
520 }
521 else
522 outside = TRUE;
523 }
524 }
525
526 if (outside == TRUE || (icH && DXIsElementInvalid(icH, seg)))
527 break;
528
529 li->hint = seg;
530
531 #define INTERPOLATE(type, round) \
532 { \
533 type *d0, *d1, *r; \
534 \
535 r = (type *)v; \
536 \
537 d0 = (type *)DXGetArrayEntry(li->dHandle, i0, (Pointer)dbuf); \
538 d1 = (type *)DXGetArrayEntry(li->dHandle, i1, \
539 (Pointer)(dbuf+itemSize)); \
540 for (i = 0; i < li->nElements; i++) \
541 *r++ = (*d0++ * linesCoord.p) + (*d1++ * linesCoord.q) + round; \
542 \
543 v = (Pointer)r; \
544 }
545
546 if (((FieldInterpolator)li)->cstData)
547 {
548 memcpy(v, ((FieldInterpolator)li)->cstData, itemSize);
549 v = (Pointer)(((char *)v) + itemSize);
550 }
551 else if (dep == DATA_POSITIONS_DEPENDENT)
552 {
553 Type dataType;
554
555 if ((dataType = ((Interpolator)li)->type) == TYPE_FLOAT)
556 {
557 INTERPOLATE(float, 0.0);
558 }
559 else if (dataType == TYPE_DOUBLE)
560 {
561 INTERPOLATE(double, 0.0);
562 }
563 else if (dataType == TYPE_INT)
564 {
565 INTERPOLATE(int, 0.5);
566 }
567 else if (dataType == TYPE_SHORT)
568 {
569 INTERPOLATE(short, 0.5);
570 }
571 else if (dataType == TYPE_USHORT)
572 {
573 INTERPOLATE(ushort, 0.5);
574 }
575 else if (dataType == TYPE_UINT)
576 {
577 INTERPOLATE(uint, 0.5);
578 }
579 else if (dataType == TYPE_BYTE)
580 {
581 INTERPOLATE(byte, 0.5);
582 }
583 else if (dataType == TYPE_UBYTE)
584 {
585 INTERPOLATE(ubyte, 0.5);
586 }
587 else
588 {
589 INTERPOLATE(unsigned char, 0.5);
590 }
591 }
592 else
593 {
594 memcpy(v, DXGetArrayEntry(li->dHandle, seg, (Pointer)dbuf),
595 itemSize);
596 v = (Pointer)(((char *)v) + itemSize);
597 }
598
599 /*
600 * Only use fuzz first time around
601 */
602 fuzzFlag = FUZZ_OFF;
603
604 p ++;
605 *n -= 1;
606 }
607
608 DXFree(dbuf);
609
610 *values = (Pointer)v;
611 *points = (float *)p;
612
613 return OK;
614 }
615
616 static int
_dxfCleanup(LinesRI1DInterpolator li)617 _dxfCleanup(LinesRI1DInterpolator li)
618 {
619 if (li->pHandle)
620 {
621 DXFreeArrayHandle(li->pHandle);
622 li->pHandle = NULL;
623 }
624
625 if (li->dHandle)
626 {
627 DXFreeArrayHandle(li->dHandle);
628 li->dHandle = NULL;
629 }
630
631 if (li->dataArray)
632 {
633 DXDelete((Object)li->dataArray);
634 li->dataArray = NULL;
635 }
636
637 if (li->pointsArray)
638 {
639 DXDelete((Object)li->pointsArray);
640 li->pointsArray = NULL;
641 }
642
643 if (li->lines)
644 {
645 DXFreeArrayDataLocal(li->linesArray, (Pointer)li->lines);
646 li->lines = NULL;
647 }
648
649 if (li->linesArray)
650 {
651 DXDelete((Object)li->linesArray);
652 li->linesArray = NULL;
653 }
654
655 return OK;
656 }
657
658 static int
linear_coords(float pt,float p0,float p1,LinearCoord * b,float fuzz)659 linear_coords (float pt, float p0, float p1, LinearCoord *b, float fuzz)
660 {
661 int code;
662
663 b->p = 0.0;
664 b->q = 0.0;
665
666 if (p0-fuzz > pt)
667 code = -1;
668 else if (p1+fuzz < pt)
669 code = 1;
670 else
671 {
672 if (p0 == p1)
673 {
674 b->p = 1.0;
675 b->q = 0.0;
676 }
677 else
678 {
679 b->q = (pt - p0) / (p1 - p0);
680 b->p = 1.0 - b->q;
681 }
682
683 if (p0 > pt)
684 code = -2;
685 else if (p1 <= pt)
686 code = 2;
687 else
688 code = 0;
689 }
690
691 return code;
692 }
693
694
695 Interpolator
_dxfLinesRI1DInterpolator_LocalizeInterpolator(LinesRI1DInterpolator li)696 _dxfLinesRI1DInterpolator_LocalizeInterpolator(LinesRI1DInterpolator li)
697 {
698 if (li->fieldInterpolator.localized)
699 return (Interpolator)li;
700
701 li->fieldInterpolator.localized = 1;
702
703 if (li->fieldInterpolator.initialized)
704 {
705 li->lines = (int *)DXGetArrayDataLocal(li->linesArray);
706 }
707 return (Interpolator)li;
708 }
709
710 Object
_dxfLinesRI1DInterpolator_Copy(LinesRI1DInterpolator old,enum _dxd_copy copy)711 _dxfLinesRI1DInterpolator_Copy(LinesRI1DInterpolator old, enum _dxd_copy copy)
712 {
713 LinesRI1DInterpolator new;
714
715 new = (LinesRI1DInterpolator)
716 _dxf_NewObject((struct object_class *)&_dxdlinesri1dinterpolator_class);
717
718 if (!(_dxf_CopyLinesRI1DInterpolator(new, old, copy)))
719 {
720 DXDelete((Object)new);
721 return NULL;
722 }
723 else
724 return (Object)new;
725 }
726
727 LinesRI1DInterpolator
_dxf_CopyLinesRI1DInterpolator(LinesRI1DInterpolator new,LinesRI1DInterpolator old,enum _dxd_copy copy)728 _dxf_CopyLinesRI1DInterpolator(LinesRI1DInterpolator new,
729 LinesRI1DInterpolator old, enum _dxd_copy copy)
730 {
731 if (! _dxf_CopyFieldInterpolator((FieldInterpolator)new,
732 (FieldInterpolator)old, copy))
733 return NULL;
734
735 new->nPoints = old->nPoints;
736 new->nLines = old->nLines;
737 new->nElements = old->nElements;
738 new->hint = old->hint;
739 new->direction = old->direction;
740
741 if (new->fieldInterpolator.initialized)
742 {
743 new->pointsArray = (Array)DXReference((Object)old->pointsArray);
744 new->linesArray = (Array)DXReference((Object)old->linesArray);
745 new->dataArray = (Array)DXReference((Object)old->dataArray);
746
747 if (new->fieldInterpolator.localized)
748 {
749 new->lines = (int *)DXGetArrayDataLocal(new->linesArray);
750 }
751 else
752 {
753 new->lines = (int *)DXGetArrayData(new->linesArray);
754 }
755
756 new->pHandle = DXCreateArrayHandle(new->pointsArray);
757 if (! new->pHandle)
758 return NULL;
759
760 new->dHandle = DXCreateArrayHandle(old->dataArray);
761 if (! new->dHandle)
762 return NULL;
763 }
764 else
765 {
766 new->pointsArray = NULL;
767 new->linesArray = NULL;
768 new->dataArray = NULL;
769 new->lines = NULL;
770 new->dHandle = NULL;
771 new->pHandle = NULL;
772 }
773
774 if (DXGetError())
775 return NULL;
776 else
777 return new;
778 }
779
780