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 * $Header: /src/master/dx/src/exec/dxmods/_refinereg.c,v 1.5 2000/08/24 20:04:19 davidt Exp $
14 */
15
16 #include <stdio.h>
17 #include <string.h>
18 #include "math.h"
19 #include <dx/dx.h>
20 #include "_refine.h"
21
22 #define MAXDIM 16
23
24 static Array RefineDepPos(Array, int, int, int *, int *);
25 static Array RefineDepPosIrreg(Array, int, int, int *, int *);
26 static Array RefineDepCon(Array, int, int, int *, int *);
27 static Array RefineDepConIrreg(Array, int, int, int *, int *);
28 static Array RefineRegReferences(Array, int, int, int *, int *, int);
29
30 #if 0
31 static Array RefineDepConBoolean(Array, int, int, int *, int *);
32 #endif
33
34 #define REF_POSITIONS 1
35 #define REF_CONNECTIONS 2
36
37 Field
_dxfRefineReg(Field field,int n)38 _dxfRefineReg(Field field, int n)
39 {
40 Array inArray, outArray;
41 int nDim;
42 int inCounts[MAXDIM], outCounts[MAXDIM];
43 int i;
44 char *name;
45 Object attr;
46 char *str;
47 float origin[MAXDIM];
48 float deltas[MAXDIM*MAXDIM];
49 int counts[MAXDIM*MAXDIM];
50 int regDim, totItems;
51 int meshOffsets[MAXDIM];
52 int irreg_positions;
53
54 /*
55 * Create refined regular positions and connections. In the subsequent
56 * loop we will traverse the list of component arrays refining all
57 * non-grid arrays in place. Since we cannot refine grid arrays in place,
58 * we need to do so before we enter the later loop. We now that
59 * connections are grid or we wouldn't have gotten here; we need to test
60 * positions.
61 */
62 inArray = (Array)DXGetComponentValue(field, "connections");
63 if (! inArray)
64 goto error;
65
66 if (DXGetArrayClass(inArray) != CLASS_MESHARRAY)
67 goto error;
68
69 DXQueryGridConnections(inArray, &nDim, inCounts);
70
71 totItems = 1;
72 for (i = 0; i < nDim; i++)
73 {
74 int k;
75
76 outCounts[i] = ((inCounts[i] - 1) * (1 << n)) + 1;
77
78 k = totItems * outCounts[i];
79 if (k / totItems != outCounts[i])
80 {
81 int j;
82 for (j = 1; j < 32; j++)
83 {
84 outCounts[i] = ((inCounts[i] - 1) * (1 << n)) + 1;
85 k = totItems * outCounts[i];
86 if (k / totItems != outCounts[i])
87 break;
88 }
89
90 DXSetError(ERROR_BAD_PARAMETER, "#10040", "level", 0, (n-1));
91 goto error;
92 }
93
94 totItems = k;
95 }
96
97 if (NULL == (outArray = DXMakeGridConnectionsV(nDim, outCounts)))
98 goto error;
99
100 if (DXGetMeshOffsets((MeshArray)inArray, meshOffsets))
101 {
102 for (i = 0; i < nDim; i++)
103 meshOffsets[i] = meshOffsets[i] << n;
104
105 DXSetMeshOffsets((MeshArray)outArray, meshOffsets);
106 }
107
108 DXSetComponentValue(field, "connections", (Object)outArray);
109
110 inArray = (Array)DXGetComponentValue(field, "positions");
111 if (! inArray)
112 goto error;
113
114 if (DXQueryGridPositions(inArray, ®Dim, counts, origin, deltas))
115 {
116 for (i = 0; i < regDim*regDim; i++)
117 deltas[i] /= (1 << n);
118
119 totItems = 1;
120 for (i = 0; i < regDim; i++)
121 if (counts[i])
122 {
123 int k;
124
125 counts[i] = ((counts[i] - 1) * (1 << n)) + 1;
126
127 k = totItems * counts[i];
128 if (k / totItems != counts[i])
129 {
130 int j;
131 for (j = 1; j < 32; j++)
132 {
133 outCounts[i] = ((inCounts[i] - 1) * (1 << n)) + 1;
134 k = totItems * outCounts[i];
135 if (k / totItems != outCounts[i])
136 break;
137 }
138
139 DXSetError(ERROR_BAD_PARAMETER, "#10040", "level", 0, (n-1));
140 goto error;
141 }
142
143 totItems = k;
144 }
145
146 outArray = DXMakeGridPositionsV(regDim, counts, origin, deltas);
147 if (NULL == outArray)
148 goto error;
149
150 DXSetComponentValue(field, "positions", (Object)outArray);
151 irreg_positions = 0;
152 }
153 else
154 irreg_positions = 1;
155
156 i = 0;
157 while(NULL !=
158 (inArray=(Array)DXGetEnumeratedComponentValue(field, i++, &name)))
159 {
160 Array out = inArray;
161
162 /*
163 * Skip connections and regular positions
164 */
165 if (!strcmp(name, "connections") ||
166 !strcmp(name, "box") ||
167 (!strcmp(name, "positions") && !irreg_positions))
168 continue;
169
170 attr = DXGetComponentAttribute(field, name, "dep");
171 if (attr)
172 {
173 if (! DXGetComponentAttribute(field, name, "ref"))
174 {
175 str = DXGetString((String)attr);
176
177 /*
178 * _dxfRefine it in place if it is dep on positions
179 * or connections and is not positions (unless
180 * positions are irregular) or connections
181 */
182
183 if (!strcmp(str, "positions"))
184 {
185 if (strcmp(name, "positions") || irreg_positions)
186 {
187 out = RefineDepPos(inArray, n, nDim,
188 inCounts, outCounts);
189
190 if (! out)
191 goto error;
192 }
193 }
194 else if (strcmp(name, "connections")
195 && !strcmp(str, "connections"))
196 {
197 out = RefineDepCon(inArray, n, nDim, inCounts, outCounts);
198
199 if (! out)
200 goto error;
201 }
202 }
203 }
204 else if (NULL != (attr = DXGetComponentAttribute(field, name, "ref")))
205 {
206 str = DXGetString((String)attr);
207
208 if (!strcmp(str, "positions"))
209 {
210 out = RefineRegReferences(inArray, nDim, n,
211 inCounts, outCounts, REF_POSITIONS);
212
213 if (! out)
214 goto error;
215 }
216 else if (strcmp(name, "connections")
217 && !strcmp(str, "connections"))
218 {
219 out = RefineRegReferences(inArray, nDim, n,
220 inCounts, outCounts, REF_CONNECTIONS);
221
222 if (! out)
223 goto error;
224 }
225 }
226
227 if (! DXSetComponentValue(field, name, (Object)out))
228 {
229 DXDelete((Object)out);
230 goto error;
231 }
232 }
233
234 DXDeleteComponent(field, "data statistics");
235
236 if (! DXEndField(field))
237 goto error;
238
239 return field;
240
241 error:
242
243 return NULL;
244 }
245
246 static Array
RefineDepCon(Array inArray,int n,int nDim,int * inCounts,int * outCounts)247 RefineDepCon(Array inArray, int n, int nDim, int *inCounts, int *outCounts)
248 {
249 Array outArray;
250
251 outArray = NULL;
252
253 /*
254 * If its a constant array, then so's the result
255 */
256 if (DXQueryConstantArray(inArray, NULL, NULL))
257 {
258 int num, i, rank, shape[32];
259 Type type;
260 Category cat;
261
262 DXGetArrayInfo(inArray, NULL, &type, &cat, &rank, shape);
263
264 num = 1;
265 for (i = 0; i < nDim; i++)
266 num *= outCounts[i];
267
268 outArray = (Array)DXNewConstantArrayV(num,
269 DXGetConstantArrayData(inArray), type, cat, rank, shape);
270
271 if (! outArray)
272 return NULL;
273 }
274 /*
275 * If its a varying regular array and
276 * nDim == 1, we can create a regular output.
277 */
278 else if (DXGetArrayClass(inArray) == CLASS_REGULARARRAY && nDim == 1)
279 {
280 int num, i, size, aDim;
281 Type type;
282 Pointer o, d;
283
284 o = d = NULL;
285
286 size = DXGetItemSize(inArray);
287
288 o = DXAllocate(size);
289 d = DXAllocate(size);
290 if (! o || ! d)
291 goto block_error;
292
293 DXGetRegularArrayInfo((RegularArray)inArray, NULL, o, d);
294 DXGetArrayInfo(inArray, NULL, &type, NULL, NULL, &aDim);
295
296
297 #define REDUCE_DELTAS(type) \
298 { \
299 type *ptr = (type *)d; \
300 for (i = 0; i < aDim; i++) \
301 ptr[i] = ((*inCounts - 1)) / ((*outCounts) - 1); \
302 }
303
304 if (*inCounts > 1)
305 switch(type) {
306 case TYPE_DOUBLE: REDUCE_DELTAS(double); break;
307 case TYPE_FLOAT: REDUCE_DELTAS(float); break;
308 case TYPE_INT: REDUCE_DELTAS(int); break;
309 case TYPE_UINT: REDUCE_DELTAS(uint); break;
310 case TYPE_SHORT: REDUCE_DELTAS(short); break;
311 case TYPE_USHORT: REDUCE_DELTAS(ushort); break;
312 case TYPE_BYTE: REDUCE_DELTAS(byte); break;
313 case TYPE_UBYTE: REDUCE_DELTAS(ubyte); break;
314 default: break;
315 }
316
317 #undef REDUCE_DELTAS
318
319 num = 1;
320 for (i = 0; i < nDim; i++)
321 num *= outCounts[i];
322
323 outArray = (Array)DXNewRegularArray(type, aDim, num, o, d);
324
325 block_error:
326 DXFree((Pointer)d);
327 DXFree((Pointer)o);
328 }
329 else
330 {
331 outArray = RefineDepConIrreg(inArray, n, nDim, inCounts, outCounts);
332 }
333
334 return outArray;
335 }
336
337 static Array
RefineDepConIrreg(Array inArray,int n,int nDim,int * inCounts,int * outCounts)338 RefineDepConIrreg(Array inArray, int n, int nDim, int *inCounts, int *outCounts)
339 {
340 Array outArray;
341 Type t; Category c; int r, s[32];
342 int i, j, k, nInItems, nOutItems, itemSize, count;
343 int oCC[32], skip[32];
344 int iknt[32], oknt[32];
345 char *ibase[32], *obase[32];
346 char *inBase, *outBase, *oPtr, *inBuffer;
347
348 inBuffer = NULL;
349 outArray = NULL;
350
351 /*
352 * Get stats on input array and create an output array of the same type
353 */
354 DXGetArrayInfo(inArray, &nInItems, &t, &c, &r, s);
355 outArray = DXNewArrayV(t, c, r, s);
356 if (! outArray)
357 goto error;
358
359 /*
360 * Count is the number of output samples along a single axis
361 * produced from each input sample.
362 */
363 count = 1 << n;
364
365 /*
366 * Number of output samples is the number of input samples times the
367 * number of output samples derived from each. This, in turn, is the
368 * count along one axis raised to the power of the dimensionality of
369 * the space.
370 */
371 nOutItems = nInItems * pow((double)count, (double)nDim);
372
373 /*
374 * Get number of output values along each axis. Since we are dep on
375 * connections, this is the outCounts - 1, since the outCounts count the
376 * number of positions along the axis. While we are at it, reverse the
377 * indices.
378 */
379 for (i = 0; i < nDim; i++)
380 oCC[(nDim-1) - i] = outCounts[i] - 1;
381
382 itemSize = DXGetItemSize(inArray);
383
384 /*
385 * Get the input data.
386 */
387 inBuffer = (char *)DXGetArrayDataLocal(inArray);
388 if (! inBuffer)
389 {
390 DXResetError();
391 inBuffer = (char *)DXGetArrayData(inArray);
392 }
393
394 /*
395 * DXAllocate the output array buffer
396 */
397 if (! DXAddArrayData(outArray, 0, nOutItems, NULL))
398 goto error;
399
400 /*
401 * The outermost loop will index through the input samples. Within
402 * this loop we will follow a base that represents the first sample
403 * location in the output in which to copy the input sample. Given
404 * this base we enter a nDim-dimensional loop which stuffs the input
405 * sample into each appropriate output slot. To do so we need to set
406 * up some skip sizes. Here we do the index reversal: remember, Z
407 * varies fastest. So the Z skip is the item size, the Y skip is the
408 * item size times the Z count (eg count[nDim-1]) and so on.
409 */
410 skip[0] = itemSize;
411 for (i = 1; i < nDim; i++)
412 skip[i] = skip[i-1] * oCC[i-1];
413
414 /*
415 * And off we go.
416 */
417 inBase = inBuffer;
418 outBase = (char * )DXGetArrayData(outArray);
419
420 for (i = 0; i < nDim; i++)
421 {
422 obase[i] = outBase;
423 oknt[i] = 0;
424 }
425
426 for (i = 0; i < nInItems; i++)
427 {
428 for (j = 0; j < nDim; j++)
429 {
430 ibase[j] = outBase;
431 iknt[j] = 0;
432 }
433
434 while(1)
435 {
436 oPtr = ibase[0];
437 for (j = 0; j < count; j++)
438 {
439 memcpy(oPtr, inBase, itemSize);
440 oPtr += itemSize;
441 }
442
443 for (j = 1; j < nDim; j++)
444 {
445 iknt[j] ++;
446 if (iknt[j] < count)
447 break;
448 }
449
450 if (j >= nDim)
451 break;
452
453 ibase[j] += skip[j];
454
455 for (k = j-1; k >= 0; --k)
456 {
457 ibase[k] = ibase[j];
458 iknt[k] = 0;
459 }
460 }
461
462 inBase += itemSize;
463
464 for (j = 0; j < nDim; j++)
465 {
466 oknt[j] += count;
467 if (oknt[j] < oCC[j])
468 break;
469 }
470
471 if (j >= nDim)
472 break;
473
474 obase[j] += count * skip[j];
475
476 for (k = j-1; k >= 0; k--)
477 {
478 obase[k] = obase[j];
479 oknt[k] = 0;
480 }
481
482 outBase = obase[0];
483 }
484
485 if (inBuffer != (char *)DXGetArrayData(inArray))
486 DXFreeArrayDataLocal(inArray, (Pointer)inBuffer);
487
488 return outArray;
489
490 error:
491
492 DXDelete((Object)outArray);
493 DXFree((Pointer)inBuffer);
494
495 return NULL;
496 }
497
498 #if 0
499 static Array
500 RefineDepConBoolean(Array inArray, int n,
501 int nDim, int *inCounts, int *outCounts)
502 {
503 Array outArray;
504 Type t; Category c; int r, s[32];
505 int i, j, k, nInItems, nOutItems, itemSize, zSize, count;
506 int oCC[32], skip[32];
507 int iknt[32], oknt[32];
508 char *ibase[32], *obase[32];
509 char *inBase, *outBase, *oPtr, *inBuffer;
510
511 inBuffer = NULL;
512 outArray = NULL;
513
514 /*
515 * Get stats on input array and create an output array of the same type
516 */
517 DXGetArrayInfo(inArray, &nInItems, &t, &c, &r, s);
518
519 if ((t != TYPE_BYTE && t != TYPE_UBYTE) || c != CATEGORY_REAL)
520 {
521 DXSetError(ERROR_DATA_INVALID,
522 "invalid data must be byte or ubyte reals");
523 goto error;
524 }
525
526 outArray = DXNewArrayV(t, c, r, s);
527 if (! outArray)
528 goto error;
529
530 /*
531 * Count is the number of output samples along a single axis
532 * produced from each input sample.
533 */
534 count = 1 << n;
535
536 /*
537 * Number of output samples is the number of input samples times the
538 * number of output samples derived from each. This, in turn, is the
539 * count along one axis raised to the power of the dimensionality of
540 * the space.
541 */
542 nOutItems = nInItems * pow((double)count, (double)nDim);
543
544 /*
545 * Get number of output values along each axis. Since we are dep on
546 * connections, this is the outCounts - 1, since the outCounts count the
547 * number of positions along the axis. While we are at it, reverse the
548 * indices.
549 */
550 for (i = 0; i < nDim; i++)
551 oCC[(nDim-1) - i] = outCounts[i] - 1;
552
553 itemSize = DXGetItemSize(inArray);
554
555 if (itemSize != 1)
556 {
557 DXSetError(ERROR_DATA_INVALID,
558 "invalid data must be scalar or 1-vector");
559 goto error;
560 }
561
562 /*
563 * Get the input data.
564 */
565 inBuffer = (char *)DXGetArrayDataLocal(inArray);
566 if (! inBuffer)
567 {
568 DXResetError();
569 inBuffer = (char *)DXGetArrayData(inArray);
570 }
571
572 /*
573 * DXAllocate the output array buffer
574 */
575 if (! DXAddArrayData(outArray, 0, nOutItems, NULL))
576 goto error;
577
578 /*
579 * The outermost loop will index through the input samples. Within
580 * this loop we will follow a base that represents the first sample
581 * location in the output in which to copy the input sample. Given
582 * this base we enter a nDim-dimensional loop which stuffs the input
583 * sample into each appropriate output slot. To do so we need to set
584 * up some skip sizes. Here we do the index reversal: remember, Z
585 * varies fastest. So the Z skip is the item size, the Y skip is the
586 * item size times the Z count (eg count[nDim-1]) and so on.
587 */
588 skip[0] = itemSize;
589 for (i = 1; i < nDim; i++)
590 skip[i] = skip[i-1] * oCC[i-1];
591
592 /*
593 * And off we go.
594 */
595 inBase = inBuffer;
596 outBase = (char * )DXGetArrayData(outArray);
597
598 for (i = 0; i < nDim; i++)
599 {
600 obase[i] = outBase;
601 oknt[i] = 0;
602 }
603
604 for (i = 0; i < nInItems; i++)
605 {
606 for (j = 0; j < nDim; j++)
607 {
608 ibase[j] = outBase;
609 iknt[j] = 0;
610 }
611
612 while(1)
613 {
614 oPtr = ibase[0];
615 for (j = 0; j < count; j++)
616 *oPtr++ = *inBase;
617
618 for (j = 1; j < nDim; j++)
619 {
620 iknt[j] ++;
621 if (iknt[j] < count)
622 break;
623 }
624
625 if (j >= nDim)
626 break;
627
628 ibase[j] += skip[j];
629
630 for (k = j-1; k >= 0; --k)
631 {
632 ibase[k] = ibase[j];
633 iknt[k] = 0;
634 }
635 }
636
637 inBase += itemSize;
638
639 for (j = 0; j < nDim; j++)
640 {
641 oknt[j] += count;
642 if (oknt[j] < oCC[j])
643 break;
644 }
645
646 if (j >= nDim)
647 break;
648
649 obase[j] += count * skip[j];
650
651 for (k = j-1; k >= 0; k--)
652 {
653 obase[k] = obase[j];
654 oknt[k] = 0;
655 }
656
657 outBase = obase[0];
658 }
659
660 if (inBuffer != (char *)DXGetArrayData(inArray))
661 DXFreeArrayDataLocal(inArray, (Pointer)inBuffer);
662
663 return outArray;
664
665 error:
666
667 DXDelete((Object)outArray);
668 DXFree((Pointer)inBuffer);
669
670 return NULL;
671 }
672 #endif
673
674
675 static Array
RefineDepPos(Array inArray,int n,int nDim,int * inCounts,int * outCounts)676 RefineDepPos(Array inArray, int n, int nDim, int *inCounts, int *outCounts)
677 {
678 Array outArray;
679
680 outArray = NULL;
681
682 /*
683 * If its a constant array, then so's the result
684 */
685 if (DXQueryConstantArray(inArray, NULL, NULL))
686 {
687 int num, i, rank, shape[32];
688 Type type;
689 Category cat;
690
691 DXGetArrayInfo(inArray, NULL, &type, &cat, &rank, shape);
692
693 num = 1;
694 for (i = 0; i < nDim; i++)
695 num *= outCounts[i];
696
697 outArray = (Array)DXNewConstantArrayV(num,
698 DXGetConstantArrayData(inArray), type, cat, rank, shape);
699
700 if (! outArray)
701 return NULL;
702 }
703 /*
704 * If its a varying regular array and
705 * nDim == 1, we can create a regular output.
706 */
707 else if (DXGetArrayClass(inArray) == CLASS_REGULARARRAY && nDim == 1)
708 {
709 int num, i, size, aDim;
710 Type type;
711 Pointer o, d;
712
713 o = d = NULL;
714
715 size = DXGetItemSize(inArray);
716
717 o = DXAllocate(size);
718 d = DXAllocate(size);
719 if (! o || ! d)
720 goto block_error;
721
722 DXGetRegularArrayInfo((RegularArray)inArray, NULL, o, d);
723 DXGetArrayInfo(inArray, NULL, &type, NULL, NULL, &aDim);
724
725
726 #define REDUCE_DELTAS(type) \
727 { \
728 type *dptr = (type *)d; \
729 type *optr = (type *)o; \
730 for (i = 0; i < aDim; i++) \
731 dptr[i] = ((*inCounts - 1)) / ((*outCounts) - 1); \
732 optr[i] -= (dptr[i]/2.0); \
733 }
734
735 if (*inCounts > 1)
736 switch(type) {
737 case TYPE_DOUBLE: REDUCE_DELTAS(double); break;
738 case TYPE_FLOAT: REDUCE_DELTAS(float); break;
739 case TYPE_INT: REDUCE_DELTAS(int); break;
740 case TYPE_UINT: REDUCE_DELTAS(uint); break;
741 case TYPE_SHORT: REDUCE_DELTAS(short); break;
742 case TYPE_USHORT: REDUCE_DELTAS(ushort); break;
743 case TYPE_BYTE: REDUCE_DELTAS(byte); break;
744 case TYPE_UBYTE: REDUCE_DELTAS(ubyte); break;
745 default: break;
746 }
747
748 #undef REDUCE_DELTAS
749
750 num = 1;
751 for (i = 0; i < nDim; i++)
752 num *= outCounts[i];
753
754 outArray = (Array)DXNewRegularArray(type, aDim, num, o, d);
755
756 block_error:
757 DXFree((Pointer)d);
758 DXFree((Pointer)o);
759 }
760 else
761 {
762 outArray = RefineDepPosIrreg(inArray, n, nDim, inCounts, outCounts);
763 }
764
765 return outArray;
766 }
767
768 static Array
RefineDepPosIrreg(Array inArray,int n,int nDim,int * inCounts,int * outCounts)769 RefineDepPosIrreg(Array inArray, int n, int nDim, int *inCounts, int *outCounts)
770 {
771 Array outArray;
772 Type type;
773 Category cat;
774 int nInItems, nOutItems, r, s[MAXDIM];
775 byte *dIn;
776 float *dOut;
777 float *src;
778 float *dst;
779 int i, j, k, skip[MAXDIM], counts[MAXDIM];
780 int level;
781 int olim, oskp;
782 int nBytes, nValues;
783 int oddCount, o;
784 float div;
785 float *buffer;
786 struct
787 {
788 int inc;
789 int start;
790 int limit;
791 int skip;
792 float *ptr;
793 } outerLoop[MAXDIM], innerLoop[MAXDIM], *loop0, *loop1;
794 static float divisors[] = {1.0/1.0, 1.0/2.0,
795 1.0/4.0, 1.0/8.0,
796 1.0/16.0, 1.0/32.0};
797
798 buffer = NULL;
799 outArray = NULL;
800
801 DXGetArrayInfo(inArray, &nInItems, &type, &cat, &r, s);
802 if (cat != CATEGORY_REAL)
803 {
804 DXSetError(ERROR_NOT_IMPLEMENTED, "#11150", "dependent components");
805 goto error;
806 }
807
808 outArray = DXNewArrayV(type, cat, r, s);
809
810 nBytes = DXGetItemSize(inArray);
811 nValues = nBytes / DXTypeSize(type);
812
813 nOutItems = 1;
814 for (i = 0; i < nDim; i++)
815 nOutItems *= outCounts[i];
816
817 /*
818 * We are going to hit each output sample once for each input sample
819 * that affects it. Create a local buffer (if possible) for the
820 * output values. Leave the input in global memory.
821 */
822 buffer = (float *)DXAllocateLocal(nOutItems * nValues * sizeof(float));
823 if (! buffer)
824 {
825 DXResetError();
826 buffer = (float *)DXAllocate(nOutItems * nValues * sizeof(float));
827 }
828 if (! buffer)
829 goto error;
830
831 /*
832 * We will be accumulating in the buffer, so we must zero it first.
833 */
834 memset(buffer, 0, nOutItems * nValues * sizeof(float));
835
836 dIn = DXGetArrayData(inArray);
837 dOut = buffer;
838
839 for (i = 0; i < nDim; i++)
840 counts[i] = inCounts[(nDim-1)-i];
841
842 skip[0] = nValues * (1 << n);
843 for (i = 1; i < nDim; i++)
844 skip[i] = outCounts[nDim-i] * skip[i-1];
845
846 /*
847 * First we copy the input data into the appropriate slots
848 * of the output data array
849 */
850 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
851 {
852 loop0->ptr = dOut;
853 loop0->skip = skip[i];
854 loop0->inc = 0;
855 loop0->limit = counts[i];
856 }
857
858 olim = outerLoop[0].limit;
859 oskp = skip[0];
860
861 for (;;)
862 {
863 dst = outerLoop[0].ptr;
864
865 #define COPYDATA(type) \
866 { \
867 type *tPtr; \
868 register int j; \
869 \
870 tPtr = (type *)dIn; \
871 for (i = 0; i < olim; i++) \
872 { \
873 for (j = 0; j < nValues; j++) \
874 *dst++ += *tPtr++; \
875 \
876 dst += oskp - nValues; \
877 dIn = (Pointer)((char *)dIn + nBytes); \
878 } \
879 }
880
881 switch(type)
882 {
883 case TYPE_DOUBLE: COPYDATA(double); break;
884 case TYPE_FLOAT: COPYDATA(float); break;
885 case TYPE_UINT: COPYDATA(uint); break;
886 case TYPE_INT: COPYDATA(int); break;
887 case TYPE_USHORT: COPYDATA(ushort); break;
888 case TYPE_SHORT: COPYDATA(short); break;
889 case TYPE_UBYTE: COPYDATA(ubyte); break;
890 case TYPE_BYTE: COPYDATA(byte); break;
891 default: break;
892 }
893
894 #undef COPYDATA
895
896 if (nDim == 0) break;
897
898 for (i = 1, loop0 = outerLoop + 1; i < nDim; i++, loop0++)
899 {
900 loop0->inc ++;
901 loop0->ptr += loop0->skip;
902 if (loop0->inc < loop0->limit)
903 break;
904 }
905
906 if (i == nDim)
907 break;
908
909 for (j = i-1, loop1 = outerLoop + (i-1); j >= 0; j--, loop1--)
910 {
911 loop1->ptr = loop0->ptr;
912 loop1->inc = 0;
913 }
914 }
915
916 /*
917 * Now for each level, we insert midpoints. Assume at this point that
918 * the outer loop is set up to index each input point.
919 */
920 for (level = 0; level < n; level++)
921 {
922
923 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
924 {
925 loop0->ptr = dOut;
926 loop0->inc = 0;
927
928 innerLoop[i].skip = loop0->skip >> 1;
929 }
930
931 olim = outerLoop[0].limit;
932 oskp = outerLoop[0].skip;
933
934 /*
935 * Outer loop: for each input sample....
936 */
937 for (;;)
938 {
939 /*
940 * Inner loop: accrue input samples onto 3x...x3 neighbors
941 */
942 for (i = 0, src = outerLoop[0].ptr; i < olim; i++, src += oskp)
943 {
944 outerLoop[0].inc = i;
945
946 dst = src;
947
948 /*
949 * Set up inner loop structure.
950 * Boundary conditions vary the inner loop starts and ends
951 */
952 loop0 = outerLoop;
953 loop1 = innerLoop;
954 for (j = 0; j < nDim; j++, loop0++, loop1++)
955 {
956
957 if (loop0->inc == 0)
958 loop1->start = 0;
959 else
960 {
961 loop1->start = -1;
962 dst -= loop1->skip;
963 }
964
965 if (loop0->inc == loop0->limit-1)
966 loop1->limit = 1;
967 else
968 loop1->limit = 2;
969 }
970
971 for (j = 0, loop0 = innerLoop; j < nDim; j++, loop0++)
972 {
973 loop0->inc = loop0->start;
974 loop0->ptr = dst;
975 }
976
977 /*
978 * Inner loop: add src values to neighbors.
979 */
980 for (;;)
981 {
982 dst = innerLoop[0].ptr;
983
984 for (j = innerLoop[0].inc; j < innerLoop[0].limit; j++)
985 {
986 if (src != dst) /* neighbors only, not source */
987 for (k = 0; k < nValues; k++)
988 dst[k] += src[k];
989
990 dst += innerLoop[0].skip;
991 }
992
993 if (nDim == 0)
994 break;
995
996 for (j = 1, loop0 = innerLoop + 1; j < nDim; j++, loop0++)
997 {
998 loop0->inc ++;
999 loop0->ptr += loop0->skip;
1000 if (loop0->inc < loop0->limit)
1001 break;
1002 }
1003
1004 if (j == nDim)
1005 break;
1006
1007 loop1 = innerLoop + (j-1);
1008 for (k = j-1; k >= 0; k--, loop1--)
1009 {
1010 loop1->ptr = loop0->ptr;
1011 loop1->inc = loop1->start;
1012 }
1013 }
1014 }
1015
1016 /*
1017 * Update outer loop counters
1018 */
1019 if (nDim == 0) break;
1020
1021 for (i = 1, loop0 = outerLoop + 1; i < nDim; i++, loop0++)
1022 {
1023 loop0->inc ++;
1024 loop0->ptr += loop0->skip;
1025 if (loop0->inc < loop0->limit)
1026 break;
1027 }
1028
1029 if (i == nDim)
1030 break;
1031
1032 for (j = i-1, loop1 = outerLoop + (i-1); j >= 0; j--, loop1--)
1033 {
1034 loop1->ptr = loop0->ptr;
1035 loop1->inc = 0;
1036 }
1037 }
1038
1039 /*
1040 * At this point we have accumulated all midpoint values. Need
1041 * to divide each by number of input samples that accrued there.
1042 * This is related to whether the midpoint is the midpoint of an edge,
1043 * a face, a cube ... Which is determined by how many of its indices
1044 * are odd: 0, is an input point which does not require dividing;
1045 * 1 for an edge midpoint to be divided by 2; 2 for the midpoint of a
1046 * face to be divided by 4; 3, a cube, by 8 and so forth.
1047 *
1048 * Reset outer loop for division pass. We now reset this loop
1049 * structure to index ech valid value in the output array. This
1050 * agrees with the state at the start of the levels loop, which ends
1051 * when the division for this level is done.
1052 */
1053 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
1054 {
1055 loop0->skip = loop0->skip >> 1;
1056 loop0->limit = (loop0->limit << 1) - 1;
1057 loop0->ptr = dOut;
1058 loop0->inc = 0;
1059 }
1060
1061 olim = outerLoop[0].limit;
1062 oskp = outerLoop[0].skip;
1063
1064 oddCount = 0;
1065 for (;;)
1066 {
1067 dst = outerLoop[0].ptr;
1068
1069 for (i = 0; i < olim; i++)
1070 {
1071 o = oddCount + (i & 0x01);
1072
1073 if (o)
1074 {
1075 div = divisors[o];
1076 for (j = 0; j < nValues; j++)
1077 dst[j] *= div;
1078 }
1079
1080 dst += oskp;
1081 }
1082
1083 if (nDim == 0) break;
1084
1085 for (i = 1, loop0 = outerLoop + 1; i < nDim; i++, loop0++)
1086 {
1087 loop0->inc ++;
1088 loop0->ptr += loop0->skip;
1089 if (loop0->inc < loop0->limit)
1090 break;
1091 }
1092
1093 if (i == nDim)
1094 break;
1095
1096 for (j = i-1, loop1 = outerLoop + (i-1); j >= 0; j--, loop1--)
1097 {
1098 loop1->ptr = loop0->ptr;
1099 loop1->inc = 0;
1100 }
1101
1102 oddCount = 0;
1103 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
1104 oddCount += (loop0->inc & 0x01);
1105 }
1106 }
1107
1108 /*
1109 * Results are now complete in buffer. Stick them into the array and
1110 * return. Do a type conversion if necessary.
1111 */
1112
1113 #define COPYOUT_CONVERSION(type) \
1114 { \
1115 register float *limit; \
1116 register float *srcPtr; \
1117 register type *dstPtr; \
1118 \
1119 if (! DXAddArrayData(outArray, 0, nOutItems, NULL)) \
1120 goto error; \
1121 \
1122 dstPtr = (type *)DXGetArrayData(outArray); \
1123 srcPtr = buffer; \
1124 \
1125 limit = srcPtr + nOutItems*nValues; \
1126 \
1127 while(srcPtr < limit) \
1128 *dstPtr++ = (type)*srcPtr++; \
1129 }
1130
1131
1132 switch(type)
1133 {
1134 case TYPE_DOUBLE: COPYOUT_CONVERSION(double); break;
1135 case TYPE_FLOAT:
1136 if (! DXAddArrayData(outArray, 0, nOutItems, (Pointer)buffer))
1137 goto error;
1138 break;
1139 case TYPE_UBYTE: COPYOUT_CONVERSION(ubyte); break;
1140 case TYPE_BYTE: COPYOUT_CONVERSION(byte); break;
1141 case TYPE_USHORT: COPYOUT_CONVERSION(ushort); break;
1142 case TYPE_SHORT: COPYOUT_CONVERSION(short); break;
1143 case TYPE_INT: COPYOUT_CONVERSION(int); break;
1144 case TYPE_UINT: COPYOUT_CONVERSION(uint); break;
1145 default: break;
1146 }
1147
1148 #undef COPYOUT_CONVERSION
1149
1150
1151 DXFree((Pointer)buffer);
1152 return outArray;
1153
1154 error:
1155
1156 DXFree((Pointer)buffer);
1157 return NULL;
1158 }
1159
1160 #if 0
1161 static Array
1162 RefineDepPosIrregBoolean(Array inArray, int n,
1163 int nDim, int *inCounts, int *outCounts)
1164 {
1165 Array outArray;
1166 Type type;
1167 Category cat;
1168 int nInItems, nOutItems, r, s[MAXDIM];
1169 byte *dIn;
1170 byte *dOut;
1171 byte *src;
1172 byte *dst;
1173 int i, j, k, skip[MAXDIM], counts[MAXDIM];
1174 int level;
1175 int olim, oskp;
1176 int nBytes, nValues;
1177 int oddCount, o;
1178 byte *buffer;
1179 struct
1180 {
1181 int inc;
1182 int start;
1183 int limit;
1184 int skip;
1185 byte *ptr;
1186 } outerLoop[MAXDIM], innerLoop[MAXDIM], *loop0, *loop1;
1187
1188 buffer = NULL;
1189 outArray = NULL;
1190
1191 DXGetArrayInfo(inArray, &nInItems, &type, &cat, &r, s);
1192 if (cat != CATEGORY_REAL)
1193 {
1194 DXSetError(ERROR_NOT_IMPLEMENTED, "#11150", "dependent components");
1195 goto error;
1196 }
1197
1198 if (type != TYPE_BYTE && type != TYPE_UBYTE)
1199 {
1200 DXSetError(ERROR_DATA_INVALID, "invalid flags must be byte or ubyte");
1201 goto error;
1202 }
1203
1204 outArray = DXNewArrayV(type, cat, r, s);
1205
1206 nBytes = 1;
1207 nValues = 1;
1208
1209 nOutItems = 1;
1210 for (i = 0; i < nDim; i++)
1211 nOutItems *= outCounts[i];
1212
1213 /*
1214 * We are going to hit each output sample once for each input sample
1215 * that affects it. Create a local buffer (if possible) for the
1216 * output values. Leave the input in global memory.
1217 */
1218 buffer = (byte *)DXAllocateLocal(nOutItems * nValues * sizeof(byte));
1219 if (! buffer)
1220 {
1221 DXResetError();
1222 buffer = (byte *)DXAllocate(nOutItems * nValues * sizeof(byte));
1223 }
1224 if (! buffer)
1225 goto error;
1226
1227 /*
1228 * We will be accumulating in the buffer, so we must zero it first.
1229 */
1230 memset(buffer, 0, nOutItems * nValues * sizeof(byte));
1231
1232 dIn = DXGetArrayData(inArray);
1233 dOut = buffer;
1234
1235 for (i = 0; i < nDim; i++)
1236 counts[i] = inCounts[(nDim-1)-i];
1237
1238 skip[0] = nValues * (1 << n);
1239 for (i = 1; i < nDim; i++)
1240 skip[i] = outCounts[nDim-i] * skip[i-1];
1241
1242 /*
1243 * First we copy the input data into the appropriate slots
1244 * of the output data array
1245 */
1246 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
1247 {
1248 loop0->ptr = dOut;
1249 loop0->skip = skip[i];
1250 loop0->inc = 0;
1251 loop0->limit = counts[i];
1252 }
1253
1254 olim = outerLoop[0].limit;
1255 oskp = skip[0];
1256
1257 for (;;)
1258 {
1259 dst = outerLoop[0].ptr;
1260 src = dIn;
1261
1262 for (i = 0; i < olim; i++)
1263 {
1264 *dst = *src;
1265 src++;
1266 dst += oskp;
1267 }
1268
1269 if (nDim == 0) break;
1270
1271 for (i = 1, loop0 = outerLoop + 1; i < nDim; i++, loop0++)
1272 {
1273 loop0->inc ++;
1274 loop0->ptr += loop0->skip;
1275 if (loop0->inc < loop0->limit)
1276 break;
1277 }
1278
1279 if (i == nDim)
1280 break;
1281
1282 for (j = i-1, loop1 = outerLoop + (i-1); j >= 0; j--, loop1--)
1283 {
1284 loop1->ptr = loop0->ptr;
1285 loop1->inc = 0;
1286 }
1287 }
1288
1289 /*
1290 * Now for each level, we insert midpoints. Assume at this point that
1291 * the outer loop is set up to index each input point.
1292 */
1293 for (level = 0; level < n; level++)
1294 {
1295
1296 for (i = 0, loop0 = outerLoop; i < nDim; i++, loop0++)
1297 {
1298 loop0->ptr = dOut;
1299 loop0->inc = 0;
1300
1301 innerLoop[i].skip = loop0->skip >> 1;
1302 }
1303
1304 olim = outerLoop[0].limit;
1305 oskp = outerLoop[0].skip;
1306
1307 /*
1308 * Outer loop: for each input sample....
1309 */
1310 for (;;)
1311 {
1312 /*
1313 * Inner loop: accrue input samples onto 3x...x3 neighbors
1314 */
1315 for (i = 0, src = outerLoop[0].ptr; i < olim; i++, src += oskp)
1316 {
1317 outerLoop[0].inc = i;
1318
1319 dst = src;
1320
1321 /*
1322 * Set up inner loop structure.
1323 * Boundary conditions vary the inner loop starts and ends
1324 */
1325 loop0 = outerLoop;
1326 loop1 = innerLoop;
1327 for (j = 0; j < nDim; j++, loop0++, loop1++)
1328 {
1329
1330 if (loop0->inc == 0)
1331 loop1->start = 0;
1332 else
1333 {
1334 loop1->start = -1;
1335 dst -= loop1->skip;
1336 }
1337
1338 if (loop0->inc == loop0->limit-1)
1339 loop1->limit = 1;
1340 else
1341 loop1->limit = 2;
1342 }
1343
1344 for (j = 0, loop0 = innerLoop; j < nDim; j++, loop0++)
1345 {
1346 loop0->inc = loop0->start;
1347 loop0->ptr = dst;
1348 }
1349
1350 /*
1351 * Inner loop: add src values to neighbors.
1352 */
1353 for (;;)
1354 {
1355 dst = innerLoop[0].ptr;
1356
1357 for (j = innerLoop[0].inc; j < innerLoop[0].limit; j++)
1358 {
1359 if (src != dst)
1360 if (*src)
1361 *dst = 1;
1362
1363 dst += innerLoop[0].skip;
1364 }
1365
1366 if (nDim == 0)
1367 break;
1368
1369 for (j = 1, loop0 = innerLoop + 1; j < nDim; j++, loop0++)
1370 {
1371 loop0->inc ++;
1372 loop0->ptr += loop0->skip;
1373 if (loop0->inc < loop0->limit)
1374 break;
1375 }
1376
1377 if (j == nDim)
1378 break;
1379
1380 loop1 = innerLoop + (j-1);
1381 for (k = j-1; k >= 0; k--, loop1--)
1382 {
1383 loop1->ptr = loop0->ptr;
1384 loop1->inc = loop1->start;
1385 }
1386 }
1387 }
1388
1389 /*
1390 * Update outer loop counters
1391 */
1392 if (nDim == 0) break;
1393
1394 for (i = 1, loop0 = outerLoop + 1; i < nDim; i++, loop0++)
1395 {
1396 loop0->inc ++;
1397 loop0->ptr += loop0->skip;
1398 if (loop0->inc < loop0->limit)
1399 break;
1400 }
1401
1402 if (i == nDim)
1403 break;
1404
1405 for (j = i-1, loop1 = outerLoop + (i-1); j >= 0; j--, loop1--)
1406 {
1407 loop1->ptr = loop0->ptr;
1408 loop1->inc = 0;
1409 }
1410 }
1411 }
1412
1413 if (! DXAddArrayData(outArray, 0, nOutItems, (Pointer)buffer))
1414 goto error;
1415
1416 DXFree((Pointer)buffer);
1417 return outArray;
1418
1419 error:
1420
1421 DXFree((Pointer)buffer);
1422 return NULL;
1423 }
1424 #endif
1425
1426 #define INDICES(ref, indices, strides, ndim) \
1427 { \
1428 int i, j = (ref); \
1429 for (i = 0; i < (ndim)-1; i++) { \
1430 (indices)[i] = j / (strides)[i]; \
1431 j = j % (strides)[i]; \
1432 } \
1433 (indices)[i] = j; \
1434 }
1435
1436 #define REFERENCE(indices, ref, strides, ndim) \
1437 { \
1438 int i, j = 0; \
1439 for (i = 0; i < (ndim)-1; i++) \
1440 j += (indices)[i]*(strides)[i]; \
1441 (ref) = j + (indices)[(ndim)-1]; \
1442 }
1443
1444
1445 static Array
RefineRegReferences(Array in,int nDim,int levels,int * inKnts,int * outKnts,int ref)1446 RefineRegReferences(Array in, int nDim, int levels,
1447 int *inKnts, int *outKnts, int ref)
1448 {
1449 Array out = NULL;
1450 Type type; Category cat;
1451 int i, nRefs, rank, shape[32], *ptr, *refs;
1452 int indices[32], max[32], min[32], inS[32], outS[32];
1453 SegListSegment *segment;
1454 SegList *segList = DXNewSegList(sizeof(int), 0, 0);
1455 if (! segList)
1456 goto error;
1457
1458 DXGetArrayInfo(in, &nRefs, &type, &cat, &rank, shape);
1459 if (type != TYPE_INT && cat != CATEGORY_REAL)
1460 {
1461 DXSetError(ERROR_DATA_INVALID, "ref arrays must be real integers");
1462 goto error;
1463 }
1464
1465 nRefs *= DXGetItemSize(in) / sizeof(int);
1466
1467 refs = (int *)DXGetArrayData(in);
1468 if (! refs)
1469 goto error;
1470
1471 DXGetArrayInfo(in, &nRefs, NULL, NULL, NULL, NULL);
1472
1473 if (ref == REF_POSITIONS)
1474 {
1475 inS[nDim-1] = outS[nDim-1] = 1;
1476 for (i = nDim-2; i >= 0; i--)
1477 {
1478 inS[i] = inKnts[i+1]* inS[i+1];
1479 outS[i] = outKnts[i+1]*outS[i+1];
1480 }
1481 }
1482 else
1483 {
1484 inS[nDim-1] = outS[nDim-1] = 1;
1485 for (i = nDim-2; i >= 0; i--)
1486 {
1487 inS[i] = (inKnts[i+1]-1)* inS[i+1];
1488 outS[i] = (outKnts[i+1]-1)*outS[i+1];
1489 }
1490 }
1491
1492 for (i = 0; i < nRefs; i++)
1493 {
1494 int overlap, nOutRefs, j, r;
1495 SegListSegment *seg;
1496
1497 if ((r = refs[i]) != -1)
1498 {
1499 INDICES(r, indices, inS, nDim);
1500
1501 for (j = 0; j < nDim; j++)
1502 indices[j] *= (1 << levels);
1503
1504 if (ref == REF_POSITIONS)
1505 {
1506 overlap = (1 << levels) - 1;
1507
1508 nOutRefs = 1;
1509 for (j = 0; j < nDim; j++)
1510 {
1511 max[j] = indices[j] + overlap;
1512 if (max[j] >= outKnts[j])
1513 max[j] = outKnts[j] - 1;
1514
1515 min[j] = indices[j] - overlap;
1516 if (min[j] < 0)
1517 min[j] = 0;
1518
1519 indices[j] = min[j];
1520
1521 nOutRefs *= (max[j] - min[j]) + 1;
1522 }
1523 }
1524 else
1525 {
1526 overlap = (1 << levels) - 1;
1527
1528 nOutRefs = 1;
1529 for (j = 0; j < nDim; j++)
1530 {
1531 min[j] = indices[j];
1532 max[j] = indices[j] + overlap;
1533
1534 nOutRefs *= (max[j] - min[j]) + 1;
1535 }
1536 }
1537
1538 seg = DXNewSegListSegment(segList, nOutRefs);
1539 if (! seg)
1540 goto error;
1541
1542 ptr = (int *)DXGetSegListSegmentPointer(seg);
1543
1544 for ( ;; )
1545 {
1546 REFERENCE(indices, *ptr, outS, nDim);
1547 ptr ++;
1548
1549 for (j = 0; j < nDim; j++)
1550 {
1551 indices[j] ++;
1552 if (indices[j] > max[j])
1553 indices[j] = min[j];
1554 else
1555 break;
1556 }
1557
1558 if (j == nDim)
1559 break;
1560 }
1561 }
1562 }
1563
1564 out = DXNewArrayV(type, cat, rank, shape);
1565 if (! out)
1566 goto error;
1567
1568 if (! DXAddArrayData(out, 0, DXGetSegListItemCount(segList), NULL))
1569 goto error;
1570
1571 ptr = (int *)DXGetArrayData(out);
1572
1573 DXInitGetNextSegListSegment(segList);
1574 while(NULL != (segment = DXGetNextSegListSegment(segList)))
1575 {
1576 int knt = DXGetSegListSegmentItemCount(segment);
1577 memcpy((char *)ptr, DXGetSegListSegmentPointer(segment),
1578 knt*sizeof(int));
1579 ptr += knt;
1580 }
1581
1582 DXDeleteSegList(segList);
1583
1584 return out;
1585
1586 error:
1587 DXDeleteSegList(segList);
1588 DXDelete((Object)out);
1589
1590 return NULL;
1591 }
1592