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 */
10
11 #include <dxconfig.h>
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <string.h>
16 #include <dx/dx.h>
17 #include "_connectgrids.h"
18
19 struct arg_radius {
20 Object ino;
21 Field base;
22 float radius;
23 float exponent;
24 Array missing;
25 };
26
27 struct arg_nearest {
28 Object ino;
29 Field base;
30 int number;
31 float *radius;
32 float radiusvalue;
33 float exponent;
34 Array missing;
35 };
36
37 struct arg_scatter {
38 Object ino;
39 Field base;
40 Array missing;
41 };
42
43
44 typedef struct matrix2D {
45 float A[2][2];
46 float b[2];
47 } Matrix2D;
48
49 typedef struct vector2D {
50 float x, y;
51 } Vector2D;
52
53 #define PI 3.14156
54 #define E .0001 /* fuzz */
55 #define FLOOR(x) ((x)+E>0? (int)((x)+E) : (int)((x)+E)-1) /* fuzzy floor */
56 #define CEIL(x) ((x)-E>0? (int)((x)-E)+1 : (int)((x)-E)) /* fuzzy ceiling */
57 #define ABS(x) ((x)<0? -(x) : (x))
58 #define RND(x) (x<0? x-0.5 : x+0.5)
59 #define POW(x,y) (pow((double)x, (double)y))
60
61 #define INSERT_NEAREST_BUFFER(TYPE) \
62 { \
63 TYPE *holder=(TYPE *)bufdata; \
64 float distancesq=0; \
65 switch(shape_base[0]) { \
66 case (1): \
67 distancesq= \
68 (gridposition1D-inputposition_ptr[l])* \
69 (gridposition1D-inputposition_ptr[l]); \
70 break; \
71 case (2): \
72 distancesq = (gridposition2D.x- \
73 inputposition_ptr[2*l])* \
74 (gridposition2D.x- \
75 inputposition_ptr[2*l]) + \
76 (gridposition2D.y- \
77 inputposition_ptr[2*l+1])* \
78 (gridposition2D.y- \
79 inputposition_ptr[2*l+1]); \
80 break; \
81 case (3): \
82 distancesq = (gridposition.x- \
83 inputposition_ptr[3*l])* \
84 (gridposition.x- \
85 inputposition_ptr[3*l]) + \
86 (gridposition.y- \
87 inputposition_ptr[3*l+1])* \
88 (gridposition.y- \
89 inputposition_ptr[3*l+1]) + \
90 (gridposition.z- \
91 inputposition_ptr[3*l+2])* \
92 (gridposition.z- \
93 inputposition_ptr[3*l+2]); \
94 break; \
95 } \
96 /* first try to insert into the already present buffer */ \
97 for (m=0; m<numvalid; m++) { \
98 if (distancesq < bufdistance[m]) { \
99 if (numvalid < numnearest) numvalid++; \
100 for (n=numvalid-1; n>m; n--) { \
101 bufdistance[n] = bufdistance[n-1]; \
102 memcpy(&holder[shape[0]*n], &holder[shape[0]*(n-1)], \
103 DXGetItemSize(newcomponent)); \
104 } \
105 bufdistance[m]=distancesq; \
106 memcpy(&holder[shape[0]*m], tmpdata, \
107 DXGetItemSize(newcomponent)); \
108 goto filled_buf##TYPE; \
109 } \
110 } \
111 if (numvalid == 0) { \
112 bufdistance[0] = distancesq; \
113 memcpy(&holder[0], tmpdata, DXGetItemSize(newcomponent)); \
114 numvalid++; \
115 } \
116 filled_buf##TYPE: \
117 tmpdata += datastep; \
118 continue; \
119 } \
120
121
122 #define ACCUMULATE_DATA(TYPE) \
123 { \
124 TYPE *holder; \
125 float dis; \
126 float radiussq = radiusvalue*radiusvalue; \
127 int dc; \
128 for (ii=0, holder=(TYPE *)bufdata; ii<numvalid; ii++) { \
129 dis = sqrt(bufdistance[ii]); \
130 if (bufdistance[ii] != 0.0) { \
131 if (radius) { \
132 if (bufdistance[ii] <= radiussq) { \
133 distanceweight += 1.0/POW(dis, exponent); \
134 for (dc=0; dc<shape[0]; dc++) { \
135 data[dc] += \
136 (float)(holder[dc])/POW(dis, exponent); \
137 } \
138 } \
139 } \
140 else { \
141 distanceweight += 1.0/(POW(dis,exponent)); \
142 for (dc=0; dc<shape[0]; dc++) { \
143 data[dc] += (float)holder[dc]/POW(dis, exponent); \
144 } \
145 } \
146 } \
147 else { \
148 distanceweight = 1.0; \
149 for (dc=0; dc<shape[0]; dc++) { \
150 data[dc] = (float)holder[dc]; \
151 } \
152 break; \
153 } \
154 holder+=shape[0]; \
155 } \
156 } \
157
158 #define CHECK_AND_ALLOCATE_MISSING(missing,size,type,type_enum,shape,ptr) \
159 { \
160 ptr = (type *)DXAllocateLocal(size); \
161 if (!ptr) \
162 goto error; \
163 if(missing){ \
164 if (!DXExtractParameter((Object)missing, type_enum, \
165 shape, 1, (Pointer)ptr)) { \
166 DXSetError(ERROR_DATA_INVALID, \
167 "missing must match all position-dependent components; does not match %s", name); \
168 goto error; \
169 } \
170 } \
171 else{ \
172 DXWarning("Grid missing values assigned to 0 and marked invalid"); \
173 for(i=0;i<shape;i++) ptr[i]=0; \
174 } \
175 } \
176
177 static Error ConnectNearest(Group, Field, int, float *, float, float, Array);
178 static Error FillBuffersIrreg(int, int, ubyte *, float *, float *,
179 float *, InvalidComponentHandle,
180 InvalidComponentHandle, int *, int, float);
181 static Error FillDataArray(Type, Array, Array, ubyte *, float *,
182 int, int, int, float, Array, char *);
183 static Error ConnectRadiusField(Pointer);
184 static Error ConnectNearestField(Pointer);
185 static Error ConnectScatterField(Pointer);
186 static Error ConnectScatter(Field, Field, Array);
187 static Error ConnectRadius(Field, Field, float, float, Array);
188 static Error ConnectRadiusRegular(Field, Field, float, float, Array);
189 static Error ConnectRadiusIrregular(Field, Field, float, float, Array);
190
Vec2D(double x,double y)191 static Vector2D Vec2D(double x, double y)
192 {
193 Vector2D v;
194 v.x = x; v.y = y;
195 return v;
196 }
197
Apply2D(Vector2D vec,Matrix2D mat)198 static Vector2D Apply2D(Vector2D vec, Matrix2D mat)
199 {
200 Vector2D result;
201
202 result.x = vec.x*mat.A[0][0] + vec.y*mat.A[1][0] + mat.b[0];
203 result.y = vec.x*mat.A[0][1] + vec.y*mat.A[1][1] + mat.b[1];
204 return result;
205 }
206
Invert2D(Matrix2D matrix2d)207 static Matrix2D Invert2D(Matrix2D matrix2d)
208 {
209 float a,b,c,d,notsing;
210 Matrix2D invert;
211 Vector2D vec2d, invertb;
212
213 a = matrix2d.A[0][0];
214 b = matrix2d.A[0][1];
215 c = matrix2d.A[1][0];
216 d = matrix2d.A[1][1];
217
218 /* check singularity */
219 notsing = a*d-b*c;
220 if (notsing) {
221 invert.A[0][0] = (1.0/notsing)*(d);
222 invert.A[0][1] = (1.0/notsing)*(-b);
223 invert.A[1][0] = (1.0/notsing)*(-c);
224 invert.A[1][1] = (1.0/notsing)*(a);
225 invert.b[0] = 0;
226 invert.b[1] = 0;
227 }
228
229 invertb.x = matrix2d.b[0];
230 invertb.y = matrix2d.b[1];
231 vec2d = Apply2D(invertb, invert);
232 invert.b[0] = -vec2d.x;
233 invert.b[1] = -vec2d.y;
234 return invert;
235 }
236
237
_dxfConnectNearestObject(Object ino,Object base,int numnearest,float * radius,float exponent,Array missing)238 Error _dxfConnectNearestObject(Object ino, Object base,
239 int numnearest,
240 float *radius, float exponent,
241 Array missing)
242 {
243 /* recurse on base to the field level */
244 struct arg_nearest arg;
245 Object subbase;
246 int i;
247
248 switch (DXGetObjectClass(base)) {
249 case CLASS_FIELD:
250 /* create the task group */
251 arg.ino = ino;
252 arg.base = (Field)base;
253 arg.number = numnearest;
254 arg.radius = radius;
255 arg.exponent = exponent;
256 if (radius)
257 arg.radiusvalue = *radius;
258 else
259 arg.radiusvalue = 0.0;
260 arg.missing = missing;
261 if (!DXAddTask(ConnectNearestField,(Pointer)&arg, sizeof(arg),0.0))
262 goto error;
263 break;
264 case CLASS_GROUP:
265 for (i=0; (subbase=DXGetEnumeratedMember((Group)base, i, NULL)); i++) {
266 if (!(_dxfConnectNearestObject((Object)ino, (Object)subbase,
267 numnearest, radius, exponent, missing)))
268 goto error;
269 }
270 break;
271 default:
272 DXSetError(ERROR_DATA_INVALID,"base must be a field or a group");
273 goto error;
274 }
275 return OK;
276
277 error:
278 return ERROR;
279
280 }
281
282
_dxfConnectRadiusObject(Object ino,Object base,float radius,float exponent,Array missing)283 Error _dxfConnectRadiusObject(Object ino, Object base,
284 float radius, float exponent,
285 Array missing)
286 /* recurse on the base to the field level */
287 {
288 struct arg_radius arg;
289 Object subbase;
290 int i;
291
292
293 switch (DXGetObjectClass(base)) {
294 case CLASS_FIELD:
295 /* create the task group */
296 arg.ino = ino;
297 arg.base = (Field)base;
298 arg.radius = radius;
299 arg.exponent = exponent;
300 arg.missing = missing;
301 if (!DXAddTask(ConnectRadiusField,(Pointer)&arg, sizeof(arg),0.0))
302 goto error;
303 break;
304 case CLASS_GROUP:
305 for (i=0; (subbase=DXGetEnumeratedMember((Group)base, i, NULL)); i++) {
306 if (!(_dxfConnectRadiusObject((Object)ino, (Object)subbase,
307 radius, exponent, missing)))
308 goto error;
309 }
310 break;
311 default:
312 DXSetError(ERROR_DATA_INVALID,"base must be a field or a group");
313 goto error;
314 }
315 return OK;
316
317 error:
318 return ERROR;
319
320 }
321
ConnectNearestField(Pointer p)322 static Error ConnectNearestField(Pointer p)
323 {
324 struct arg_nearest *arg = (struct arg_nearest *)p;
325 Object ino, subino, newino=NULL;
326 int i, numnearest;
327 Field base;
328 float *radius;
329 float radiusvalue;
330 float exponent;
331 Array missing;
332
333 ino = arg->ino;
334 newino = DXApplyTransform(ino,NULL);
335
336 switch (DXGetObjectClass(newino)) {
337
338 case CLASS_FIELD:
339
340 base = arg->base;
341 radius = arg->radius;
342 radiusvalue = arg->radiusvalue;
343 missing = arg->missing;
344 numnearest = arg->number;
345 exponent = arg->exponent;
346 if (!ConnectNearest((Group)newino, base, numnearest, radius,
347 radiusvalue, exponent, missing))
348 goto error;
349 break;
350
351 case CLASS_GROUP:
352
353 switch (DXGetGroupClass((Group)newino)) {
354 case CLASS_COMPOSITEFIELD:
355 /* for now, can't handle a composite field */
356 DXSetError(ERROR_NOT_IMPLEMENTED,"input cannot be a composite field");
357 goto error;
358 base = arg->base;
359 radius = arg->radius;
360 radiusvalue = arg->radiusvalue;
361 missing = arg->missing;
362 numnearest = arg->number;
363 exponent = arg->exponent;
364 if (!ConnectNearest((Group)newino, base, numnearest, radius,
365 radiusvalue, exponent, missing))
366 goto error;
367 break;
368 default:
369 DXSetError(ERROR_NOT_IMPLEMENTED,"input must be a single field");
370 goto error;
371 /* recurse on the members of the group */
372 for (i=0; (subino=DXGetEnumeratedMember((Group)newino, i, NULL)); i++) {
373 base = arg->base;
374 radius = arg->radius;
375 radiusvalue = arg->radiusvalue;
376 missing = arg->missing;
377 numnearest = arg->number;
378 exponent = arg->exponent;
379 if (!ConnectNearest((Group)subino, base, numnearest, radius,
380 radiusvalue, exponent, missing))
381 goto error;
382 }
383 break;
384 }
385 break;
386 default:
387 DXSetError(ERROR_DATA_INVALID,"input must be a field or a group");
388 goto error;
389 }
390 DXDelete((Object)newino);
391 return OK;
392
393 error:
394 DXDelete((Object)newino);
395 return ERROR;
396
397 }
398
399
ConnectRadiusField(Pointer p)400 static Error ConnectRadiusField(Pointer p)
401 {
402 struct arg_radius *arg = (struct arg_radius *)p;
403 Object ino, newino=NULL;
404 Field base;
405 float radius;
406 float exponent;
407 Array missing;
408
409 /* this one recurses on object ino */
410
411 ino = arg->ino;
412 newino = DXApplyTransform(ino,NULL);
413
414 switch (DXGetObjectClass(newino)) {
415
416 case CLASS_FIELD:
417 base = arg->base;
418 radius = arg->radius;
419 missing = arg->missing;
420 exponent = arg->exponent;
421 if (!ConnectRadius((Field)newino, base, radius, exponent, missing))
422 goto error;
423 break;
424
425 case CLASS_GROUP:
426 DXSetError(ERROR_NOT_IMPLEMENTED,"input cannot be a group");
427 goto error;
428 break;
429
430 default:
431 DXSetError(ERROR_DATA_INVALID,"input must be a field");
432 goto error;
433 }
434 DXDelete((Object)newino);
435 return OK;
436
437 error:
438 DXDelete((Object)newino);
439 return ERROR;
440
441 }
442
ConnectNearest(Group ino,Field base,int numnearest,float * radius,float radiusvalue,float exponent,Array missing)443 static Error ConnectNearest(Group ino, Field base, int numnearest,
444 float *radius, float radiusvalue,
445 float exponent, Array missing)
446 {
447 int componentsdone, anyinvalid=0;
448 Matrix matrix;
449 Matrix2D matrix2D;
450 Array positions_base, positions_ino, oldcomponent=NULL, newcomponent;
451 InvalidComponentHandle out_invalidhandle=NULL, in_invalidhandle=NULL;
452 int counts[8], rank_base, shape_base[8], shape_ino[8], rank_ino;
453 int rank, shape[8], i, j, k, l, m, n, numvalid, pos_count=0;
454 int numitemsbase, compcount, numitemsino, numitems;
455 float distanceweight, distance=0, *data=NULL, *dp, *scratch=NULL;
456 ArrayHandle handle=NULL;
457 Type type;
458 Point gridposition = {0, 0, 0};
459 Vector2D gridposition2D = {0, 0};
460 float gridposition1D=0, origin=0, delta=0;
461 float *bufdistance=NULL;
462 char *bufdata=NULL;
463 char *name;
464 Category category;
465 float *inputposition_ptr;
466 int data_count, invalid_count=0, ii, datastep, dc, regular;
467 char *old_data_ptr=NULL, *tmpdata=NULL;
468 float *new_data_ptr_f=NULL, *missingval_f=NULL;
469 int *new_data_ptr_i=NULL, *missingval_i=NULL;
470 short *new_data_ptr_s=NULL, *missingval_s=NULL;
471 double *new_data_ptr_d=NULL, *missingval_d=NULL;
472 byte *new_data_ptr_b=NULL, *missingval_b=NULL;
473 ubyte *new_data_ptr_ub=NULL, *missingval_ub=NULL;
474 ushort *new_data_ptr_us=NULL, *missingval_us=NULL;
475 uint *new_data_ptr_ui=NULL, *missingval_ui=NULL;
476
477 if (DXEmptyField((Field)ino))
478 return OK;
479
480 positions_base = (Array)DXGetComponentValue((Field)base, "positions");
481 if (!positions_base) {
482 DXSetError(ERROR_MISSING_DATA,"base field has no positions");
483 goto error;
484 }
485 positions_ino = (Array)DXGetComponentValue((Field)ino, "positions");
486 if (!positions_ino) {
487 DXSetError(ERROR_MISSING_DATA,"input field has no positions");
488 goto error;
489 }
490
491 DXGetArrayInfo(positions_base, &numitemsbase, &type, &category,
492 &rank_base,
493 shape_base);
494
495 regular = 0;
496 if ((DXQueryGridPositions((Array)positions_base,NULL,counts,
497 (float *)&matrix.b,(float *)&matrix.A) )) {
498 regular = 1;
499 if (shape_base[0]==2) {
500 matrix2D.A[0][0] = matrix.A[0][0];
501 matrix2D.A[1][0] = matrix.A[0][1];
502 matrix2D.A[0][1] = matrix.A[0][2];
503 matrix2D.A[1][1] = matrix.A[1][0];
504 matrix2D.b[0] = matrix.b[0];
505 matrix2D.b[1] = matrix.b[1];
506 }
507 else if (shape_base[0]==1) {
508 origin = matrix.b[0];
509 delta = matrix.A[0][0];
510 }
511 }
512 else {
513 /* irregular positions */
514 if (!(handle = DXCreateArrayHandle(positions_base)))
515 goto error;
516 scratch = DXAllocateLocal(DXGetItemSize(positions_base));
517 if (!scratch) goto error;
518 }
519
520 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
521 DXSetError(ERROR_DATA_INVALID,
522 "positions component must be type float, category real");
523 goto error;
524 }
525
526
527
528 out_invalidhandle = DXCreateInvalidComponentHandle((Object)base, NULL,
529 "positions");
530 if (!out_invalidhandle)
531 goto error;
532 DXSetAllValid(out_invalidhandle);
533 in_invalidhandle = DXCreateInvalidComponentHandle((Object)ino, NULL,
534 "positions");
535 if (!in_invalidhandle)
536 return ERROR;
537
538
539 DXGetArrayInfo(positions_ino, &numitemsino, &type, &category, &rank_ino,
540 shape_ino);
541
542
543 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
544 DXSetError(ERROR_DATA_INVALID,
545 "positions component must be type float, category real");
546 goto error;
547 }
548 if (rank_base == 0)
549 shape_base[0] = 1;
550 if (rank_base > 1) {
551 DXSetError(ERROR_DATA_INVALID,
552 "rank of base positions cannot be greater than 1");
553 goto error;
554 }
555 if (rank_ino == 0)
556 shape_ino[0] = 1;
557 if (rank_ino > 1) {
558 DXSetError(ERROR_DATA_INVALID,
559 "rank of input positions cannot be greater than 1");
560 goto error;
561 }
562
563 if ((shape_base[0]<0)||(shape_base[0]>3)) {
564 DXSetError(ERROR_DATA_INVALID,
565 "only 1D, 2D, and 3D positions for base supported");
566 goto error;
567 }
568
569 if (shape_base[0] != shape_ino[0]) {
570 DXSetError(ERROR_DATA_INVALID,
571 "dimensionality of base (%dD) does not match dimensionality of input (%dD)",
572 shape_base[0], shape_ino[0]);
573 goto error;
574 }
575
576
577 if ((shape_ino[0]<0)||(shape_ino[0]>3)) {
578 DXSetError(ERROR_DATA_INVALID,
579 "only 1D, 2D, and 3D positions for input supported");
580 goto error;
581 }
582 /* first set the extra counts to one to simplify the code for handling
583 1, 2, and 3 dimensions */
584 if (regular) {
585 if (shape_base[0]==1) {
586 counts[1]=1;
587 counts[2]=1;
588 }
589 else if (shape_base[0]==2) {
590 counts[2]=1;
591 }
592 }
593 else {
594 counts[0]=numitemsbase;
595 counts[1]=1;
596 counts[2]=1;
597 }
598
599 /* need to get the components we are going to copy into the new field */
600 /* copied wholesale from gda */
601 compcount = 0;
602 componentsdone = 0;
603 while (NULL !=
604 (oldcomponent=(Array)DXGetEnumeratedComponentValue((Field)ino,
605 compcount,
606 &name))) {
607 Object attr;
608
609 data_count=0;
610 invalid_count=0;
611
612
613 if (DXGetObjectClass((Object)oldcomponent) != CLASS_ARRAY)
614 goto component_done;
615
616 if (!strcmp(name,"positions") ||
617 !strcmp(name,"connections") ||
618 !strcmp(name,"invalid positions"))
619 goto component_done;
620
621 /* if the component refs anything, leave it */
622 if (DXGetComponentAttribute((Field)ino,name,"ref"))
623 goto component_done;
624
625 attr = DXGetComponentAttribute((Field)ino,name,"der");
626 if (attr)
627 goto component_done;
628
629 attr = DXGetComponentAttribute((Field)ino,name,"dep");
630 if (!attr) {
631 /* does it der? if so, throw it away */
632 attr = DXGetComponentAttribute((Field)ino,name,"der");
633 if (attr)
634 goto component_done;
635 /* it doesn't ref, doesn't dep and doesn't der.
636 Copy it to the output and quit */
637 if (!DXSetComponentValue((Field)base,name,(Object)oldcomponent))
638 goto error;
639 goto component_done;
640 }
641
642 if (DXGetObjectClass(attr) != CLASS_STRING) {
643 DXSetError(ERROR_DATA_INVALID, "dependency attribute not type string");
644 goto error;
645 }
646
647 if (strcmp(DXGetString((String)attr),"positions")){
648 DXWarning("Component '%s' is not dependent on positions!",name);
649 DXWarning("Regrid will remove '%s' and output the base grid 'data' component if it exists",name);
650 goto component_done;
651 }
652
653
654 /* if the component happens to be "invalid positions" ignore it */
655 if (!strcmp(name,"invalid positions"))
656 goto component_done;
657
658
659 componentsdone++;
660
661 DXGetArrayInfo((Array)oldcomponent,&numitems, &type,
662 &category, &rank, shape);
663 if (rank==0) shape[0]=1;
664
665
666
667 /* check that the missing component, if present, matches the component
668 we're working on */
669 if (missing)
670 switch (type) {
671 case (TYPE_FLOAT):
672 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
673 float, type, shape[0], missingval_f);
674 break;
675 case (TYPE_INT):
676 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
677 int, type, shape[0], missingval_i);
678 break;
679 case (TYPE_DOUBLE):
680 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
681 double, type, shape[0], missingval_d);
682 break;
683 case (TYPE_SHORT):
684 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
685 short, type, shape[0], missingval_s);
686 break;
687 case (TYPE_BYTE):
688 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
689 byte, type, shape[0], missingval_b);
690 break;
691 case (TYPE_UINT):
692 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
693 uint, type, shape[0], missingval_ui);
694 break;
695 case (TYPE_USHORT):
696 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
697 ushort, type, shape[0], missingval_us);
698 break;
699 case (TYPE_UBYTE):
700 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
701 ubyte, type, shape[0], missingval_ub);
702 break;
703 default:
704 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
705 goto error;
706 }
707
708
709 if (rank == 0)
710 shape[0] = 1;
711 if (rank > 1) {
712 DXSetError(ERROR_NOT_IMPLEMENTED,"rank larger than 1 not implemented");
713 goto error;
714 }
715 if (numitems != numitemsino) {
716 DXSetError(ERROR_BAD_PARAMETER,
717 "number of items in %s component not equal to number of positions",
718 name);
719 goto error;
720 }
721 newcomponent = DXNewArrayV(type, category, rank, shape);
722 if (!DXAddArrayData(newcomponent, 0, numitemsbase, NULL))
723 goto error;
724 old_data_ptr = (char *)DXGetArrayData(oldcomponent);
725
726 switch (type) {
727 case (TYPE_FLOAT):
728 new_data_ptr_f = (float *)DXGetArrayData(newcomponent);
729 break;
730 case (TYPE_INT):
731 new_data_ptr_i = (int *)DXGetArrayData(newcomponent);
732 break;
733 case (TYPE_SHORT):
734 new_data_ptr_s = (short *)DXGetArrayData(newcomponent);
735 break;
736 case (TYPE_BYTE):
737 new_data_ptr_b = (byte *)DXGetArrayData(newcomponent);
738 break;
739 case (TYPE_DOUBLE):
740 new_data_ptr_d = (double *)DXGetArrayData(newcomponent);
741 break;
742 case (TYPE_UINT):
743 new_data_ptr_ui = (uint *)DXGetArrayData(newcomponent);
744 break;
745 case (TYPE_UBYTE):
746 new_data_ptr_ub = (ubyte *)DXGetArrayData(newcomponent);
747 break;
748 case (TYPE_USHORT):
749 new_data_ptr_us = (ushort *)DXGetArrayData(newcomponent);
750 break;
751 default:
752 DXSetError(ERROR_DATA_INVALID,"unrecognized data type");
753 goto error;
754 }
755
756
757
758 /* need to allocate a buffer of length n to hold the nearest neighbors*/
759 bufdata =
760 (char *)DXAllocateLocal(numnearest*DXGetItemSize(newcomponent)*shape[0]);
761 if (!bufdata) goto error;
762 bufdistance = (float *)DXAllocateLocal(numnearest*sizeof(float));
763 if (!bufdistance) goto error;
764 data = (float *)DXAllocateLocalZero(sizeof(float)*shape[0]);
765 if (!data) goto error;
766 datastep = DXGetItemSize(newcomponent);
767
768 /* loop over all the grid points */
769 inputposition_ptr = (float *)DXGetArrayData(positions_ino);
770 for (i=0; i<counts[0]; i++) {
771 for (j=0; j<counts[1]; j++) {
772 for (k=0; k<counts[2]; k++) {
773 if (regular) {
774 switch (shape_base[0]) {
775 case (3):
776 gridposition = DXApply(DXVec(i,j,k), matrix);
777 break;
778 case (2):
779 gridposition2D = Apply2D(Vec2D(i,j), matrix2D);
780 break;
781 case (1):
782 gridposition1D = origin + i*delta;
783 }
784 }
785 else {
786 switch(shape_base[0]) {
787 case (3):
788 dp = DXGetArrayEntry(handle, pos_count, scratch);
789 gridposition = DXPt(dp[0], dp[1], dp[2]);
790 break;
791 case (2):
792 dp = DXGetArrayEntry(handle, pos_count, scratch);
793 gridposition2D.x = dp[0];
794 gridposition2D.y = dp[1];
795 break;
796 case (1):
797 dp = DXGetArrayEntry(handle, pos_count, scratch);
798 gridposition1D = dp[0];
799 break;
800 }
801 pos_count++;
802 }
803 for (ii=0; ii<numnearest; ii++) {
804 bufdistance[ii] = DXD_MAX_FLOAT;
805 }
806 /* numvalid is how many valid data points are in the buffer
807 for this particular grid point */
808 numvalid = 0;
809 for (l=0, tmpdata=old_data_ptr; l<numitemsino; l++) {
810 /* don't do anything if the input is invalid */
811 if (DXIsElementValid(in_invalidhandle, l)) {
812 switch (type) {
813 case (TYPE_FLOAT):
814 INSERT_NEAREST_BUFFER(float);
815 break;
816 case (TYPE_INT):
817 INSERT_NEAREST_BUFFER(int);
818 break;
819 case (TYPE_SHORT):
820 INSERT_NEAREST_BUFFER(short);
821 break;
822 case (TYPE_DOUBLE):
823 INSERT_NEAREST_BUFFER(double);
824 break;
825 case (TYPE_BYTE):
826 INSERT_NEAREST_BUFFER(byte);
827 break;
828 case (TYPE_UBYTE):
829 INSERT_NEAREST_BUFFER(ubyte);
830 break;
831 case (TYPE_USHORT):
832 INSERT_NEAREST_BUFFER(ushort);
833 break;
834 case (TYPE_UINT):
835 INSERT_NEAREST_BUFFER(uint);
836 break;
837 default:
838 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
839 goto error;
840 }
841 }
842 }
843 /* have now looped over all the scattered points */
844 distanceweight = 0;
845 memset(data, 0, shape[0]*sizeof(float));
846
847
848 switch (type) {
849 case (TYPE_FLOAT):
850 ACCUMULATE_DATA(float);
851 break;
852 case (TYPE_INT):
853 ACCUMULATE_DATA(int);
854 break;
855 case (TYPE_BYTE):
856 ACCUMULATE_DATA(byte);
857 break;
858 case (TYPE_DOUBLE):
859 ACCUMULATE_DATA(double);
860 break;
861 case (TYPE_SHORT):
862 ACCUMULATE_DATA(short);
863 break;
864 case (TYPE_UINT):
865 ACCUMULATE_DATA(uint);
866 break;
867 case (TYPE_USHORT):
868 ACCUMULATE_DATA(ushort);
869 break;
870 case (TYPE_UBYTE):
871 ACCUMULATE_DATA(ubyte);
872 break;
873 default:
874 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
875 goto error;
876 }
877
878 /* now stick the data value into the new data component */
879 switch (type) {
880 case (TYPE_FLOAT):
881 if (distanceweight != 0.) {
882 for (dc=0; dc<shape[0]; dc++) {
883 new_data_ptr_f[data_count++] =
884 (float)(data[dc]/distanceweight);
885 }
886 }
887 else {
888 if (missing) {
889 for (dc=0; dc<shape[0]; dc++) {
890 new_data_ptr_f[data_count++] = missingval_f[dc];
891 }
892 }
893 else {
894 /* use the invalid positions array and set the data
895 value to zero */
896 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
897 goto error;
898 for (dc=0; dc<shape[0]; dc++) {
899 new_data_ptr_f[data_count++] = 0;
900 }
901 }
902 }
903 break;
904 case (TYPE_INT):
905 if (distanceweight != 0) {
906 for (dc=0; dc<shape[0]; dc++) {
907 new_data_ptr_i[data_count++] =
908 (int)(RND(data[dc]/distanceweight));
909 }
910 }
911 else {
912 if (missing) {
913 for (dc=0; dc<shape[0]; dc++) {
914 new_data_ptr_i[data_count++] = missingval_i[dc];
915 }
916 }
917 else {
918 /* use the invalid positions array and set the data
919 value to zero */
920 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
921 goto error;
922 for (dc=0; dc<shape[0]; dc++) {
923 new_data_ptr_i[data_count++] = 0;
924 }
925 }
926 }
927 break;
928 case (TYPE_SHORT):
929 if (distanceweight != 0) {
930 for (dc=0; dc<shape[0]; dc++) {
931 new_data_ptr_s[data_count++] =
932 (short)(RND(data[dc]/distanceweight));
933 }
934 }
935 else {
936 if (missing) {
937 for (dc=0; dc<shape[0]; dc++) {
938 new_data_ptr_s[data_count++] = missingval_s[dc];
939 }
940 }
941 else {
942 /* use the invalid positions array and set the data
943 value to zero */
944 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
945 goto error;
946 for (dc=0; dc<shape[0]; dc++) {
947 new_data_ptr_s[data_count++] = 0;
948 }
949 }
950 }
951 break;
952 case (TYPE_BYTE):
953 if (distanceweight != 0) {
954 for (dc=0; dc<shape[0]; dc++) {
955 new_data_ptr_b[data_count++] =
956 (byte)(RND(data[dc]/distanceweight));
957 }
958 }
959 else {
960 if (missing) {
961 for (dc=0; dc<shape[0]; dc++) {
962 new_data_ptr_b[data_count++] = missingval_b[dc];
963 }
964 }
965 else {
966 /* use the invalid positions array and set the data
967 value to zero */
968 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
969 goto error;
970 for (dc=0; dc<shape[0]; dc++) {
971 new_data_ptr_b[data_count++] = 0;
972 }
973 }
974 }
975 break;
976 case (TYPE_DOUBLE):
977 if (distanceweight != 0) {
978 for (dc=0; dc<shape[0]; dc++) {
979 new_data_ptr_d[data_count++] =
980 (double)(data[dc]/distanceweight);
981 }
982 }
983 else {
984 if (missing) {
985 for (dc=0; dc<shape[0]; dc++) {
986 new_data_ptr_d[data_count++] = missingval_d[dc];
987 }
988 }
989 else {
990 /* use the invalid positions array and set the data
991 value to zero */
992 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
993 goto error;
994 for (dc=0; dc<shape[0]; dc++) {
995 new_data_ptr_d[data_count++] = 0;
996 }
997 }
998 }
999 break;
1000 case (TYPE_UBYTE):
1001 if (distanceweight != 0) {
1002 for (dc=0; dc<shape[0]; dc++) {
1003 new_data_ptr_ub[data_count++] =
1004 (ubyte)(RND(data[dc]/distanceweight));
1005 }
1006 }
1007 else {
1008 if (missing) {
1009 for (dc=0; dc<shape[0]; dc++) {
1010 new_data_ptr_ub[data_count++] = missingval_ub[dc];
1011 }
1012 }
1013 else {
1014 /* use the invalid positions array and set the data
1015 value to zero */
1016 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
1017 goto error;
1018 for (dc=0; dc<shape[0]; dc++) {
1019 new_data_ptr_ub[data_count++] = 0;
1020 }
1021 }
1022 }
1023 break;
1024 case (TYPE_USHORT):
1025 if (distanceweight != 0) {
1026 for (dc=0; dc<shape[0]; dc++) {
1027 new_data_ptr_us[data_count++] =
1028 (ushort)(RND(data[dc]/distanceweight));
1029 }
1030 }
1031 else {
1032 if (missing) {
1033 for (dc=0; dc<shape[0]; dc++) {
1034 new_data_ptr_us[data_count++] = missingval_us[dc];
1035 }
1036 }
1037 else {
1038 /* use the invalid positions array and set the data
1039 value to zero */
1040 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
1041 goto error;
1042 for (dc=0; dc<shape[0]; dc++) {
1043 new_data_ptr_us[data_count++] = 0;
1044 }
1045 }
1046 }
1047 break;
1048 case (TYPE_UINT):
1049 if (distanceweight != 0) {
1050 for (dc=0; dc<shape[0]; dc++) {
1051 new_data_ptr_ui[data_count++] =
1052 (uint)(RND(data[dc]/distanceweight));
1053 }
1054 }
1055 else {
1056 if (missing) {
1057 for (dc=0; dc<shape[0]; dc++) {
1058 new_data_ptr_ui[data_count++] = missingval_ui[dc];
1059 }
1060 }
1061 else {
1062 /* use the invalid positions array and set the data
1063 value to zero */
1064 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
1065 goto error;
1066 for (dc=0; dc<shape[0]; dc++) {
1067 new_data_ptr_ui[data_count++] = 0;
1068 }
1069 }
1070 }
1071 break;
1072 default:
1073 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
1074 goto error;
1075 }
1076 /* increment the invalid positions array pointer */
1077 invalid_count++;
1078 }
1079 }
1080 }
1081
1082
1083 if (!DXChangedComponentValues((Field)base,name))
1084 goto error;
1085 if (!DXSetComponentValue((Field)base,name,(Object)newcomponent))
1086 goto error;
1087 /* XXX and radius also (otherwise inv pos not nec) */
1088 /* should inv pos be done each time? (for each component) */
1089 if ((!missing)&&(componentsdone==1)) {
1090 if (!DXSaveInvalidComponent(base, out_invalidhandle))
1091 goto error;
1092
1093 }
1094 newcomponent = NULL;
1095 if (!DXSetComponentAttribute((Field)base, name,
1096 "dep",
1097 (Object)DXNewString("positions")))
1098 goto error;
1099
1100 component_done:
1101 compcount++;
1102 DXFree((Pointer)missingval_f);
1103 missingval_f=NULL;
1104 DXFree((Pointer)missingval_s);
1105 missingval_s=NULL;
1106 DXFree((Pointer)missingval_i);
1107 missingval_i=NULL;
1108 DXFree((Pointer)missingval_d);
1109 missingval_d=NULL;
1110 DXFree((Pointer)missingval_b);
1111 missingval_b=NULL;
1112 DXFree((Pointer)missingval_ui);
1113 missingval_ui=NULL;
1114 DXFree((Pointer)missingval_ub);
1115 missingval_ub=NULL;
1116 DXFree((Pointer)missingval_us);
1117 missingval_us=NULL;
1118 DXFree((Pointer)bufdistance);
1119 bufdistance=NULL;
1120 DXFree((Pointer)bufdata);
1121 bufdata=NULL;
1122 DXFree((Pointer)data);
1123 data=NULL;
1124 }
1125 if (componentsdone==0) {
1126 if (missing) {
1127 DXWarning("no position dependent components found; missing value given cannot be used");
1128 }
1129 else {
1130 /* need to allocate a buffer of length n to hold the nearest neighbors*/
1131 bufdistance = (float *)DXAllocateLocal(numnearest*sizeof(float));
1132 if (!bufdistance) goto error;
1133
1134 /* loop over all the grid points */
1135 inputposition_ptr = (float *)DXGetArrayData(positions_ino);
1136 for (i=0; i<counts[0]; i++) {
1137 for (j=0; j<counts[1]; j++) {
1138 for (k=0; k<counts[2]; k++) {
1139 if (regular) {
1140 switch (shape_base[0]) {
1141 case (3):
1142 gridposition = DXApply(DXVec(i,j,k), matrix);
1143 break;
1144 case (2):
1145 gridposition2D = Apply2D(Vec2D(i,j), matrix2D);
1146 break;
1147 case (1):
1148 gridposition1D = origin + i*delta;
1149 break;
1150 }
1151 }
1152 else {
1153 switch(shape_base[0]) {
1154 case (3):
1155 dp = DXGetArrayEntry(handle, pos_count, scratch);
1156 gridposition = DXPt(dp[0], dp[1], dp[2]);
1157 break;
1158 case (2):
1159 dp = DXGetArrayEntry(handle, pos_count, scratch);
1160 gridposition2D.x = dp[0];
1161 gridposition2D.y = dp[1];
1162 break;
1163 case (1):
1164 dp = DXGetArrayEntry(handle, pos_count, scratch);
1165 gridposition1D = dp[0];
1166 break;
1167 }
1168 pos_count++;
1169 }
1170 for (ii=0; ii<numnearest; ii++) {
1171 bufdistance[ii] = DXD_MAX_FLOAT;
1172 }
1173 /* numvalid is how many valid data points are in the buffer
1174 for this particular grid point */
1175 numvalid = 0;
1176 for (l=0, tmpdata=old_data_ptr; l<numitemsino; l++) {
1177 /* don't do anything if the input is invalid */
1178 if (DXIsElementValid(in_invalidhandle, l)) {
1179 switch(shape_base[0]) {
1180 case (1):
1181 distance= ABS(gridposition1D-inputposition_ptr[l]);
1182 break;
1183 case (2):
1184 distance = sqrt(POW(gridposition2D.x-
1185 inputposition_ptr[2*l],2) +
1186 POW(gridposition2D.y-
1187 inputposition_ptr[2*l+1],2));
1188 break;
1189 case (3):
1190 distance = DXLength(DXSub(gridposition,
1191 DXPt(inputposition_ptr[3*l],
1192 inputposition_ptr[3*l+1],
1193 inputposition_ptr[3*l+2])));
1194 break;
1195 }
1196 /* first try to insert into the already present buffer */
1197 for (m=0; m<numvalid; m++) {
1198 if (distance < bufdistance[m]) {
1199 if (numvalid < numnearest) numvalid++;
1200 for (n=numvalid-1; n>m; n--) {
1201 bufdistance[n] = bufdistance[n-1];
1202 }
1203 bufdistance[m]=distance;
1204 goto filled_buf;
1205 }
1206 }
1207 if (numvalid == 0) {
1208 bufdistance[0] = distance;
1209 numvalid++;
1210 }
1211 filled_buf:
1212 tmpdata += datastep;
1213 continue;
1214 }
1215 }
1216
1217 /* have now looped over all the scattered points */
1218
1219 if (radius && (bufdistance[0]>radiusvalue)) {
1220 anyinvalid=1;
1221 if (!DXSetElementInvalid(out_invalidhandle, invalid_count))
1222 goto error;
1223 }
1224 /* increment the invalid positions array pointer */
1225 invalid_count++;
1226 }
1227 }
1228 }
1229 /* put the invalid component in the field */
1230 if (anyinvalid) {
1231 if (!DXSaveInvalidComponent(base, out_invalidhandle))
1232 goto error;
1233
1234 }
1235 DXFree((Pointer)bufdistance);
1236 }
1237 }
1238
1239
1240 DXFreeInvalidComponentHandle(out_invalidhandle);
1241 DXFreeInvalidComponentHandle(in_invalidhandle);
1242 DXFree((Pointer)scratch);
1243 DXFreeArrayHandle(handle);
1244 if (!DXEndField(base)) goto error;
1245 return OK;
1246
1247
1248 error:
1249 DXFreeInvalidComponentHandle(out_invalidhandle);
1250 DXFreeInvalidComponentHandle(in_invalidhandle);
1251 DXFree((Pointer)scratch);
1252 DXFreeArrayHandle(handle);
1253 DXFree((Pointer)bufdistance);
1254 DXFree((Pointer)bufdata);
1255 DXFree((Pointer)data);
1256 DXFree((Pointer)missingval_f);
1257 DXFree((Pointer)missingval_s);
1258 DXFree((Pointer)missingval_i);
1259 DXFree((Pointer)missingval_d);
1260 DXFree((Pointer)missingval_b);
1261 DXFree((Pointer)missingval_ui);
1262 DXFree((Pointer)missingval_ub);
1263 DXFree((Pointer)missingval_us);
1264 return ERROR;
1265 }
1266
1267
1268
ConnectRadius(Field ino,Field base,float radius,float exponent,Array missing)1269 static Error ConnectRadius(Field ino, Field base, float radius,
1270 float exponent, Array missing)
1271 {
1272 Array positions_base;
1273
1274 if (DXEmptyField((Field)ino))
1275 return OK;
1276
1277 positions_base = (Array)DXGetComponentValue((Field)base,"positions");
1278 if (!positions_base) {
1279 DXSetError(ERROR_MISSING_DATA,"missing positions component in base");
1280 goto error;
1281 }
1282
1283 if (DXQueryGridPositions((Array)positions_base,NULL,NULL,NULL,NULL)) {
1284 if (!ConnectRadiusRegular(ino,base,radius,exponent, missing))
1285 goto error;
1286 }
1287 else {
1288 if (!ConnectRadiusIrregular(ino,base,radius,exponent, missing))
1289 goto error;
1290 }
1291 if (!DXEndField(base)) goto error;
1292 return OK;
1293 error:
1294 return ERROR;
1295 }
1296
ConnectRadiusIrregular(Field ino,Field base,float radius,float exponent,Array missing)1297 static Error ConnectRadiusIrregular(Field ino, Field base, float radius,
1298 float exponent, Array missing)
1299 {
1300 int componentsdone, valid;
1301 Array positions_base, positions_ino;
1302 Array oldcomponent=NULL, newcomponent=NULL;
1303 InvalidComponentHandle out_invalidhandle=NULL, in_invalidhandle=NULL;
1304 int numitemsbase, numitemsino;
1305 int rank_base, rank_ino, shape_base[30], shape_ino[30], i, j;
1306 int anyinvalid=1, compcount, rank, shape[30];
1307 int numitems;
1308 Type type;
1309 ubyte *flagbuf=NULL;
1310 float *distancebuf=NULL;
1311 Category category;
1312 float *pos_ino_ptr, *pos_base_ptr;
1313 char *name;
1314
1315 /* all points within radius of a given grid position will be averaged
1316 * (weighted) to give a data point to that position. If there are are
1317 * none within that range, the missing value is used. If missing is
1318 * NULL, then an invalid position indication is placed there and cull
1319 * will be used at the end */
1320
1321 in_invalidhandle = DXCreateInvalidComponentHandle((Object)ino, NULL,
1322 "positions");
1323 if (!in_invalidhandle)
1324 goto error;
1325 if (!missing) {
1326 out_invalidhandle = DXCreateInvalidComponentHandle((Object)base,
1327 NULL, "positions");
1328 if (!out_invalidhandle)
1329 goto error;
1330 anyinvalid = 0;
1331 }
1332
1333 positions_base = (Array)DXGetComponentValue((Field)base, "positions");
1334 if (!positions_base) {
1335 DXSetError(ERROR_MISSING_DATA,"base field has no positions");
1336 goto error;
1337 }
1338 positions_ino = (Array)DXGetComponentValue((Field)ino, "positions");
1339 if (!positions_ino) {
1340 DXSetError(ERROR_MISSING_DATA,"input field has no positions");
1341 goto error;
1342 }
1343 DXGetArrayInfo(positions_base, &numitemsbase, &type, &category, &rank_base,
1344 shape_base);
1345 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
1346 DXSetError(ERROR_DATA_INVALID,
1347 "only real, floating point positions for base accepted");
1348 goto error;
1349 }
1350 if ((rank_base > 1 ) || (shape_base[0] > 3)) {
1351 DXSetError(ERROR_DATA_INVALID,
1352 "only 1D, 2D, or 3D positions for base accepted");
1353 goto error;
1354 }
1355 DXGetArrayInfo(positions_ino, &numitemsino, &type, &category, &rank_ino,
1356 shape_ino);
1357 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
1358 DXSetError(ERROR_DATA_INVALID,
1359 "only real, floating point positions for input accepted");
1360 goto error;
1361 }
1362 if ((rank_ino > 1 ) || (shape_ino[0] > 3)) {
1363 DXSetError(ERROR_DATA_INVALID,
1364 "only 1D, 2D, or 3D positions for input accepted");
1365 goto error;
1366 }
1367
1368 if (shape_base[0] != shape_ino[0]) {
1369 DXSetError(ERROR_DATA_INVALID,
1370 "base positions are %dD; input positions are %dD",
1371 shape_base[0], shape_ino[0]);
1372 goto error;
1373 }
1374
1375
1376 /* regular arrays? XXX */
1377 pos_ino_ptr = (float *)DXGetArrayData(positions_ino);
1378
1379 /* make a buffer to hold the distances and flags */
1380 flagbuf = (ubyte *)DXAllocateLocalZero((numitemsbase*numitemsino*sizeof(ubyte)));
1381 if (!flagbuf) goto error;
1382 distancebuf =
1383 (float *)DXAllocateLocalZero((numitemsbase*numitemsino*sizeof(float)));
1384 if (!distancebuf) goto error;
1385
1386 if (out_invalidhandle) {
1387 DXSetAllValid(out_invalidhandle);
1388 }
1389
1390 /* call irregular method */
1391 pos_base_ptr = (float *)DXGetArrayData(positions_base);
1392 if (!(FillBuffersIrreg(numitemsbase, numitemsino, flagbuf, distancebuf,
1393 pos_ino_ptr, pos_base_ptr, out_invalidhandle,
1394 in_invalidhandle,
1395 &anyinvalid, shape_base[0], radius)))
1396 goto error;
1397
1398
1399 /* need to get the components we are going to copy into the new field */
1400 /* copied wholesale from gda */
1401 compcount = 0;
1402 componentsdone=0;
1403 while (NULL !=
1404 (oldcomponent=(Array)DXGetEnumeratedComponentValue((Field)ino,
1405 compcount,
1406 &name))) {
1407 Object attr;
1408
1409 if (DXGetObjectClass((Object)oldcomponent) != CLASS_ARRAY)
1410 goto component_done;
1411
1412 if (!strcmp(name,"positions") ||
1413 !strcmp(name,"connections") ||
1414 !strcmp(name,"invalid positions"))
1415 goto component_done;
1416
1417 /* if the component refs anything, leave it */
1418 if (DXGetComponentAttribute((Field)ino,name,"ref"))
1419 goto component_done;
1420
1421 attr = DXGetComponentAttribute((Field)ino,name,"der");
1422 if (attr)
1423 goto component_done;
1424
1425 attr = DXGetComponentAttribute((Field)ino,name,"dep");
1426 if (!attr) {
1427 /* does it der? if so, throw it away */
1428 attr = DXGetComponentAttribute((Field)ino,name,"der");
1429 if (attr)
1430 goto component_done;
1431 /* it doesn't ref, doesn't dep and doesn't der.
1432 Copy it to the output and quit */
1433 if (!DXSetComponentValue((Field)base,name,(Object)oldcomponent))
1434 goto error;
1435 goto component_done;
1436 }
1437
1438
1439 if (DXGetObjectClass(attr) != CLASS_STRING) {
1440 DXSetError(ERROR_DATA_INVALID, "dependency attribute not type string");
1441 goto error;
1442 }
1443
1444 if (strcmp(DXGetString((String)attr),"positions")){
1445 DXWarning("Component '%s' is not dependent on positions!",name);
1446 DXWarning("Regrid will remove '%s' and output the base grid 'data' component if it exists",name);
1447 goto component_done;
1448 }
1449
1450
1451 /* if the component happens to be "invalid positions" ignore it */
1452 if (!strcmp(name,"invalid positions"))
1453 goto component_done;
1454
1455
1456 componentsdone++;
1457
1458 DXGetArrayInfo((Array)oldcomponent, &numitems, &type, &category,
1459 &rank, shape);
1460
1461
1462
1463
1464 if (rank == 0)
1465 shape[0] = 1;
1466 if (rank > 1) {
1467 DXSetError(ERROR_NOT_IMPLEMENTED,
1468 "rank larger than 1 not implemented");
1469 goto error;
1470 }
1471 if (numitems != numitemsino) {
1472 DXSetError(ERROR_BAD_PARAMETER,
1473 "number of items in %s component not equal to number of positions", name);
1474 goto error;
1475 }
1476 newcomponent = DXNewArrayV(type, category, rank, shape);
1477
1478
1479 if (!DXAddArrayData(newcomponent, 0, numitemsbase, NULL))
1480 goto error;
1481
1482
1483 /* here's where I use my preset buffers to fill the new array */
1484 /* XXX here do something about missing */
1485 if (!FillDataArray(type, oldcomponent, newcomponent,
1486 flagbuf, distancebuf,
1487 numitemsbase, numitemsino,
1488 shape[0], exponent, missing, name))
1489 goto error;
1490
1491 if (!DXChangedComponentValues((Field)base,name))
1492 goto error;
1493 if (!DXSetComponentValue((Field)base,name,(Object)newcomponent))
1494 goto error;
1495 newcomponent = NULL;
1496 if (!DXSetComponentAttribute((Field)base, name,
1497 "dep",
1498 (Object)DXNewString("positions")))
1499 goto error;
1500 component_done:
1501 compcount++;
1502 }
1503
1504
1505 if (componentsdone==0) {
1506 if (missing) {
1507 DXWarning("no position-dependent components found; missing value provided cannot be used");
1508 }
1509 else {
1510 for (i=0; i<numitemsbase; i++) {
1511 valid = 0;
1512 for (j=0; j<numitemsino; j++) {
1513 if (flagbuf[i*numitemsino+j]) {
1514 /* there's at least one valid */
1515 valid = 1;
1516 break;
1517 }
1518 }
1519 if (!valid) {
1520 DXSetElementInvalid(out_invalidhandle, i);
1521 anyinvalid=1;
1522 }
1523 }
1524 }
1525 }
1526 else {
1527 if (!missing) {
1528 for (i=0; i<numitemsbase; i++) {
1529 valid = 0;
1530 for (j=0; j<numitemsino; j++) {
1531 if (flagbuf[i*numitemsino+j]) {
1532 /* there's at least one valid */
1533 valid = 1;
1534 break;
1535 }
1536 }
1537 if (!valid) {
1538 DXSetElementInvalid(out_invalidhandle, i);
1539 anyinvalid=1;
1540 }
1541 }
1542 }
1543 }
1544
1545 if (anyinvalid && out_invalidhandle) {
1546 if (!DXSaveInvalidComponent(base, out_invalidhandle))
1547 goto error;
1548 }
1549
1550 DXFree((Pointer)flagbuf);
1551 DXFree((Pointer)distancebuf);
1552 DXFreeInvalidComponentHandle(out_invalidhandle);
1553 DXFreeInvalidComponentHandle(in_invalidhandle);
1554 return OK;
1555
1556 error:
1557 DXDelete((Object)newcomponent);
1558 DXFree((Pointer)flagbuf);
1559 DXFree((Pointer)distancebuf);
1560 DXFreeInvalidComponentHandle(out_invalidhandle);
1561 DXFreeInvalidComponentHandle(in_invalidhandle);
1562 return ERROR;
1563
1564 }
1565
ConnectRadiusRegular(Field ino,Field base,float radius,float exponent,Array missing)1566 static Error ConnectRadiusRegular(Field ino, Field base, float radius,
1567 float exponent, Array missing)
1568 {
1569 int componentsdone;
1570 Array positions_base, positions_ino;
1571 Array oldcomponent, newcomponent=NULL;
1572 InvalidComponentHandle in_invalidhandle=NULL, out_invalidhandle=NULL;
1573 int numitemsbase, numitemsino, cntr;
1574 int rank_base, rank_ino, shape_base[30], shape_ino[30], i, j, k, l, c;
1575 int compcount, rank, shape[30];
1576 float tmpval;
1577 int numitems, counts[8], knt;
1578 Type type;
1579 float xpos, ypos, zpos, xmin, ymin, zmin, xmax, ymax, zmax, pos;
1580 float min, max;
1581 Vector minvec, maxvec, posvec;
1582 Vector2D minvec2D, maxvec2D, posvec2D;
1583 Matrix matrix, invmatrix;
1584 Matrix2D matrix2D, invmatrix2D;
1585 Category category;
1586 float *pos_ino_ptr, distance, *weights_ptr=NULL;
1587 float *data_ptr = NULL;
1588 float rsquared;
1589 byte *exact_hit;
1590 Point box[8];
1591 char *name;
1592 float *data_in_f=NULL, *data_f=NULL, *missingval_f=NULL;
1593 byte *data_in_b=NULL, *data_b=NULL, *missingval_b=NULL;
1594 int *data_in_i=NULL, *data_i=NULL, *missingval_i=NULL;
1595 short *data_in_s=NULL, *data_s=NULL, *missingval_s=NULL;
1596 double *data_in_d=NULL, *data_d=NULL, *missingval_d=NULL;
1597 ubyte *data_in_ub=NULL, *data_ub=NULL, *missingval_ub=NULL;
1598 ushort *data_in_us=NULL, *data_us=NULL, *missingval_us=NULL;
1599 uint *data_in_ui=NULL, *data_ui=NULL, *missingval_ui=NULL;
1600
1601
1602
1603
1604 positions_base = (Array)DXGetComponentValue((Field)base, "positions");
1605 if (!positions_base) {
1606 DXSetError(ERROR_MISSING_DATA,"base field has no positions");
1607 goto error;
1608 }
1609 if (!(DXQueryGridPositions((Array)positions_base,NULL,counts,
1610 (float *)&matrix.b,(float *)&matrix.A) ))
1611 goto error;
1612 /* all points within radius of a given grid position will be averaged
1613 * (weighted) to give a data point to that position. If there are are
1614 * none within that range, the missing value is used. If missing is
1615 * NULL, then an invalid position indication is placed there */
1616 DXGetArrayInfo(positions_base, &numitemsbase, &type, &category,
1617 &rank_base, shape_base);
1618 out_invalidhandle=DXCreateInvalidComponentHandle((Object)base,
1619 NULL,"positions");
1620 if (!out_invalidhandle) goto error;
1621
1622
1623 DXSetAllInvalid(out_invalidhandle);
1624
1625 /* if radius==-1, that's infinity. I really ought to handle it separately */
1626 /* XXX */
1627 if (radius == -1) {
1628 if (!DXBoundingBox((Object)ino, box))
1629 goto error;
1630 radius = 2*DXLength(DXSub(box[0], box[7]));
1631 }
1632
1633 positions_ino = (Array)DXGetComponentValue((Field)ino, "positions");
1634 if (!positions_ino) {
1635 DXSetError(ERROR_MISSING_DATA,"input field has no positions");
1636 goto error;
1637 }
1638 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
1639 DXSetError(ERROR_DATA_INVALID,
1640 "only real, floating point positions for base accepted");
1641 goto error;
1642 }
1643 if ((rank_base > 1 ) || (shape_base[0] > 3)) {
1644 DXSetError(ERROR_DATA_INVALID,
1645 "only 1D, 2D, or 3D positions for base accepted");
1646 goto error;
1647 }
1648 DXGetArrayInfo(positions_ino, &numitemsino, &type, &category,
1649 &rank_ino, shape_ino);
1650 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
1651 DXSetError(ERROR_DATA_INVALID,
1652 "only real, floating point positions for input accepted");
1653 goto error;
1654 }
1655 if ((rank_ino > 1 ) || (shape_ino[0] > 3)) {
1656 DXSetError(ERROR_DATA_INVALID,
1657 "only 1D, 2D, or 3D positions for input accepted");
1658 goto error;
1659 }
1660
1661 if (shape_base[0] != shape_ino[0]) {
1662 DXSetError(ERROR_DATA_INVALID,
1663 "base positions are %dD; input positions are %dD",
1664 shape_base[0], shape_ino[0]);
1665 goto error;
1666 }
1667
1668
1669 /* regular arrays? XXX */
1670 pos_ino_ptr = (float *)DXGetArrayData(positions_ino);
1671
1672
1673
1674 /* need to get the components we are going to copy into the new field */
1675 /* copied wholesale from gda */
1676 compcount = 0;
1677 componentsdone=0;
1678 while (NULL != (oldcomponent=(Array)DXGetEnumeratedComponentValue((Field)ino,
1679 compcount,
1680 &name))) {
1681 Object attr;
1682
1683 if (DXGetObjectClass((Object)oldcomponent) != CLASS_ARRAY)
1684 goto component_done;
1685
1686 if (!strcmp(name,"positions") ||
1687 !strcmp(name,"connections") ||
1688 !strcmp(name,"invalid positions"))
1689 goto component_done;
1690
1691 /* if the component refs anything, leave it */
1692 if (DXGetComponentAttribute((Field)ino,name,"ref"))
1693 goto component_done;
1694
1695 attr = DXGetComponentAttribute((Field)ino,name,"der");
1696 if (attr)
1697 goto component_done;
1698
1699 attr = DXGetComponentAttribute((Field)ino,name,"dep");
1700 if (!attr) {
1701 /* does it der? if so, throw it away */
1702 attr = DXGetComponentAttribute((Field)ino,name,"der");
1703 if (attr)
1704 goto component_done;
1705 /* it doesn't ref, doesn't dep and doesn't der.
1706 Copy it to the output and quit */
1707 if (!DXSetComponentValue((Field)base,name,(Object)oldcomponent))
1708 goto error;
1709 goto component_done;
1710 }
1711
1712
1713
1714 if (DXGetObjectClass(attr) != CLASS_STRING) {
1715 DXSetError(ERROR_DATA_INVALID,
1716 "dependency attribute not type string");
1717 goto error;
1718 }
1719
1720 if (strcmp(DXGetString((String)attr),"positions")){
1721 DXWarning("Component '%s' is not dependent on positions!",name);
1722 DXWarning("Regrid will remove '%s' and output the base grid 'data' component if it exists",name);
1723 goto component_done;
1724 }
1725
1726
1727 /* if the component happens to be "invalid positions" ignore it */
1728 if (!strcmp(name,"invalid positions"))
1729 goto component_done;
1730
1731 componentsdone++;
1732
1733 DXGetArrayInfo((Array)oldcomponent,&numitems, &type, &category,
1734 &rank, shape);
1735 if (rank==0) shape[0]=1;
1736
1737 /* check that the missing component, if present, matches the component
1738 we're working on */
1739 if (missing)
1740 switch (type) {
1741 case (TYPE_FLOAT):
1742 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
1743 float, type, shape[0], missingval_f);
1744 break;
1745 case (TYPE_INT):
1746 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), int,
1747 type, shape[0], missingval_i);
1748 break;
1749 case (TYPE_BYTE):
1750 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), byte,
1751 type, shape[0], missingval_b);
1752 break;
1753 case (TYPE_DOUBLE):
1754 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
1755 double, type, shape[0], missingval_d);
1756 break;
1757 case (TYPE_SHORT):
1758 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), short,
1759 type, shape[0], missingval_s);
1760 break;
1761 case (TYPE_UBYTE):
1762 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), ubyte,
1763 type, shape[0], missingval_ub);
1764 break;
1765 case (TYPE_USHORT):
1766 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
1767 ushort, type, shape[0], missingval_us);
1768 break;
1769 case (TYPE_UINT):
1770 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), uint,
1771 type, shape[0], missingval_ui);
1772 break;
1773 default:
1774 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
1775 goto error;
1776 }
1777
1778
1779
1780
1781
1782 if (rank == 0)
1783 shape[0] = 1;
1784 if (rank > 1) {
1785 DXSetError(ERROR_NOT_IMPLEMENTED,"rank larger than 1 not implemented");
1786 goto error;
1787 }
1788 if (numitems != numitemsino) {
1789 DXSetError(ERROR_BAD_PARAMETER,
1790 "number of items in %s component not equal to number of positions", name);
1791 goto error;
1792 }
1793
1794 weights_ptr = (float *)DXAllocateLocalZero(numitemsbase*sizeof(float));
1795 /* I need a floating point buffer to hold the weighted sum data */
1796 data_ptr =
1797 (float *)DXAllocateLocalZero(numitemsbase*sizeof(float)*shape[0]);
1798 if (!weights_ptr) goto error;
1799 if (!data_ptr) goto error;
1800
1801 /* new component is for the eventual output data */
1802 newcomponent = DXNewArrayV(type, category, rank, shape);
1803 if (!DXAddArrayData(newcomponent, 0, numitemsbase, NULL))
1804 goto error;
1805
1806 /* allocate the exact hit buf */
1807 exact_hit = (byte *)DXAllocateLocalZero(numitemsbase*sizeof(byte));
1808
1809
1810 switch (type) {
1811 case (TYPE_FLOAT):
1812 data_in_f = (float *)DXGetArrayData(oldcomponent);
1813 data_f = (float *)DXGetArrayData(newcomponent);
1814 break;
1815 case (TYPE_INT):
1816 data_in_i = (int *)DXGetArrayData(oldcomponent);
1817 data_i = (int *)DXGetArrayData(newcomponent);
1818 break;
1819 case (TYPE_DOUBLE):
1820 data_in_d = (double *)DXGetArrayData(oldcomponent);
1821 data_d = (double *)DXGetArrayData(newcomponent);
1822 break;
1823 case (TYPE_SHORT):
1824 data_in_s = (short *)DXGetArrayData(oldcomponent);
1825 data_s = (short *)DXGetArrayData(newcomponent);
1826 break;
1827 case (TYPE_BYTE):
1828 data_in_b = (byte *)DXGetArrayData(oldcomponent);
1829 data_b = (byte *)DXGetArrayData(newcomponent);
1830 break;
1831 case (TYPE_UBYTE):
1832 data_in_ub = (ubyte *)DXGetArrayData(oldcomponent);
1833 data_ub = (ubyte *)DXGetArrayData(newcomponent);
1834 break;
1835 case (TYPE_USHORT):
1836 data_in_us = (ushort *)DXGetArrayData(oldcomponent);
1837 data_us = (ushort *)DXGetArrayData(newcomponent);
1838 break;
1839 case (TYPE_UINT):
1840 data_in_ui = (uint *)DXGetArrayData(oldcomponent);
1841 data_ui = (uint *)DXGetArrayData(newcomponent);
1842 break;
1843 default:
1844 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
1845 goto error;
1846 }
1847
1848 #define INCREMENT_DATA_1D(type, data_in) \
1849 { \
1850 for (i=0; i<numitemsino; i++) { \
1851 xpos = pos_ino_ptr[i]; \
1852 xmin = xpos-radius; \
1853 xmax = xpos+radius; \
1854 min = (xmin/matrix.A[0][0])-(matrix.b[0]/matrix.A[0][0]); \
1855 max = (xmax/matrix.A[0][0])-(matrix.b[0]/matrix.A[0][0]); \
1856 if (min > max) { \
1857 tmpval = min; \
1858 min = max; \
1859 max = tmpval; \
1860 } \
1861 /* minvec and maxvec are the indices corresponding to the \
1862 bounding box around the sample point */ \
1863 for (j=CEIL(min); j <= FLOOR(max) ; j++) \
1864 if (j>=0 && j<counts[0]) { \
1865 /* here are the only ones which are in the bounding \
1866 * box */ \
1867 pos = j*matrix.A[0][0]+matrix.b[0]; \
1868 distance = ABS(pos-xpos); \
1869 if (distance <= radius) { \
1870 DXSetElementValid(out_invalidhandle, j); \
1871 if ((distance > 0.0)&&(!exact_hit[j])) { \
1872 weights_ptr[j] += 1.0/POW(distance,exponent); \
1873 for (c = 0; c<shape[0]; c++) \
1874 data_ptr[j*shape[0]+c] += \
1875 (float)data_in[i*shape[0]+c]/POW(distance,exponent); \
1876 } \
1877 else if (!exact_hit[j]){ \
1878 /* we've hit right on; no more averaging here */ \
1879 for (c = 0; c<shape[0]; c++) \
1880 data_ptr[j*shape[0]+c] = \
1881 (float)data_in[i*shape[0]+c]; \
1882 weights_ptr[j]=1; \
1883 exact_hit[j]=1; \
1884 } \
1885 } \
1886 } \
1887 } \
1888 } \
1889
1890 #define INCREMENT_DATA_2D(type, data_in) \
1891 { \
1892 Matrix2D tmpmat; \
1893 matrix2D.A[0][0] = matrix.A[0][0]; \
1894 matrix2D.A[1][0] = matrix.A[0][1]; \
1895 matrix2D.A[0][1] = matrix.A[0][2]; \
1896 matrix2D.A[1][1] = matrix.A[1][0]; \
1897 matrix2D.b[0] = matrix.b[0]; \
1898 matrix2D.b[1] = matrix.b[1]; \
1899 invmatrix2D = Invert2D(matrix2D); \
1900 /* loop over all scattered points */ \
1901 rsquared = radius* radius; \
1902 for (i=0; i<numitemsino; i++) { \
1903 xpos = pos_ino_ptr[2*i]; \
1904 ypos = pos_ino_ptr[2*i+1]; \
1905 xmin = xpos-radius; \
1906 xmax = xpos+radius; \
1907 ymin = ypos-radius; \
1908 ymax = ypos+radius; \
1909 tmpmat = invmatrix2D; \
1910 minvec2D = Apply2D(Vec2D(xmin,ymin), tmpmat); \
1911 maxvec2D = Apply2D(Vec2D(xmax,ymax), tmpmat); \
1912 if (minvec2D.x > maxvec2D.x) { \
1913 tmpval = minvec2D.x; \
1914 minvec2D.x = maxvec2D.x; \
1915 maxvec2D.x = tmpval; \
1916 } \
1917 if (minvec2D.y > maxvec2D.y) { \
1918 tmpval = minvec2D.y; \
1919 minvec2D.y = maxvec2D.y; \
1920 maxvec2D.y = tmpval; \
1921 } \
1922 /* minvec and maxvec are the indices corresponding to \
1923 the bounding box around the sample point */ \
1924 tmpmat = matrix2D; \
1925 for (j=CEIL(minvec2D.x); j <= FLOOR(maxvec2D.x) ; j++) \
1926 if (j>=0 && j<counts[0]) \
1927 for (k=CEIL(minvec2D.y); k <= FLOOR(maxvec2D.y) ; k++) \
1928 if (k>=0 && k<counts[1]) { \
1929 /* here are the only ones which are in the \
1930 bounding box */ \
1931 posvec2D = Apply2D(Vec2D(j,k),tmpmat); \
1932 distance = (posvec2D.x-xpos)*(posvec2D.x-xpos) + \
1933 (posvec2D.y-ypos)*(posvec2D.y-ypos); \
1934 if (distance <= rsquared) { \
1935 distance = sqrt(distance); \
1936 cntr = j*counts[1]+k; \
1937 DXSetElementValid(out_invalidhandle, cntr); \
1938 if ((distance > 0.0)&&(!exact_hit[cntr])) { \
1939 weights_ptr[cntr] += 1.0/POW(distance,exponent); \
1940 for (c=0; c<shape[0]; c++) \
1941 data_ptr[j*counts[1]*shape[0] + k*shape[0] +c] \
1942 += (float)data_in[i*shape[0]+c]/POW(distance,exponent); \
1943 } \
1944 else if (!exact_hit[cntr]) { \
1945 for (c=0; c<shape[0]; c++) \
1946 data_ptr[j*counts[1]*shape[0] + k*shape[0]] \
1947 = (float)data_in[i*shape[0]+c]; \
1948 weights_ptr[cntr]=1; \
1949 exact_hit[cntr]=1; \
1950 } \
1951 } \
1952 } \
1953 } \
1954 } \
1955
1956 #define INCREMENT_DATA_3D(type, data_in) \
1957 { \
1958 Matrix tmpmat; \
1959 invmatrix = DXInvert(matrix); \
1960 /* loop over all scattered points */ \
1961 for (knt=0, i=0; i<numitemsino; i++) { \
1962 xpos = pos_ino_ptr[knt]; \
1963 ypos = pos_ino_ptr[knt+1]; \
1964 zpos = pos_ino_ptr[knt+2]; \
1965 knt+=3; \
1966 xmin = xpos-radius; \
1967 xmax = xpos+radius; \
1968 ymin = ypos-radius; \
1969 ymax = ypos+radius; \
1970 zmin = zpos-radius; \
1971 zmax = zpos+radius; \
1972 tmpmat=invmatrix; \
1973 minvec = DXApply(DXVec(xmin,ymin,zmin), tmpmat); \
1974 maxvec = DXApply(DXVec(xmax,ymax,zmax), tmpmat); \
1975 if (minvec.x > maxvec.x) { \
1976 tmpval = minvec.x; \
1977 minvec.x = maxvec.x; \
1978 maxvec.x = tmpval; \
1979 } \
1980 if (minvec.y > maxvec.y) { \
1981 tmpval = minvec.y; \
1982 minvec.y = maxvec.y; \
1983 maxvec.y = tmpval; \
1984 } \
1985 if (minvec.z > maxvec.z) { \
1986 tmpval = minvec.z; \
1987 minvec.z = maxvec.z; \
1988 maxvec.z = tmpval; \
1989 } \
1990 /* minvec and maxvec are the indices corresponding to the */ \
1991 /* bounding box around the sample point */ \
1992 tmpmat=matrix; \
1993 rsquared = radius * radius; \
1994 for (j=CEIL(minvec.x); j <= FLOOR(maxvec.x) ; j++) \
1995 if (j>=0 && j<counts[0]) \
1996 for (k=CEIL(minvec.y); k <= FLOOR(maxvec.y) ; k++) \
1997 if (k>=0 && k<counts[1]) \
1998 for (l=CEIL(minvec.z); l <= FLOOR(maxvec.z) ; l++) \
1999 if (l>=0 && l<counts[2]) { \
2000 /* here are the only ones which are in the \
2001 bounding box */ \
2002 posvec = DXApply(DXVec(j,k,l),tmpmat); \
2003 distance = (posvec.x-xpos)*(posvec.x-xpos) + \
2004 (posvec.y-ypos)*(posvec.y-ypos) + \
2005 (posvec.z-zpos)*(posvec.z-zpos); \
2006 if (distance <= rsquared) { \
2007 distance = sqrt(distance); \
2008 cntr = j*counts[2]*counts[1] + \
2009 k*counts[2] + l; \
2010 DXSetElementValid(out_invalidhandle, cntr); \
2011 if ((distance > 0.0)&&(!exact_hit[cntr])) { \
2012 weights_ptr[cntr] += 1.0/POW(distance,exponent); \
2013 for (c=0; c<shape[0]; c++) \
2014 data_ptr[j*counts[2]*counts[1]*shape[0] + \
2015 k*counts[2]*shape[0] + l*shape[0]+c] \
2016 += (float)data_in[i*shape[0]+c]/POW(distance,exponent); \
2017 } \
2018 else if (!exact_hit[cntr]){ \
2019 for (c=0; c<shape[0]; c++) \
2020 data_ptr[j*counts[2]*counts[1]*shape[0] + \
2021 k*counts[2]*shape[0] + l*shape[0] +c] \
2022 = (float)data_in[i*shape[0]+c]; \
2023 weights_ptr[cntr]=1; \
2024 exact_hit[cntr] = 1; \
2025 } \
2026 } \
2027 } \
2028 } \
2029 } \
2030
2031
2032 if (shape_base[0] == 1) {
2033 switch (type) {
2034 case(TYPE_FLOAT):
2035 INCREMENT_DATA_1D(float, data_in_f);
2036 break;
2037 case(TYPE_INT):
2038 INCREMENT_DATA_1D(int, data_in_i);
2039 break;
2040 case(TYPE_BYTE):
2041 INCREMENT_DATA_1D(byte, data_in_b);
2042 break;
2043 case(TYPE_SHORT):
2044 INCREMENT_DATA_1D(short, data_in_s);
2045 break;
2046 case(TYPE_DOUBLE):
2047 INCREMENT_DATA_1D(double, data_in_d);
2048 break;
2049 case(TYPE_UBYTE):
2050 INCREMENT_DATA_1D(ubyte, data_in_ub);
2051 break;
2052 case(TYPE_UINT):
2053 INCREMENT_DATA_1D(uint, data_in_ui);
2054 break;
2055 case(TYPE_USHORT):
2056 INCREMENT_DATA_1D(ushort, data_in_us);
2057 break;
2058 default:
2059 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2060 goto error;
2061 }
2062 }
2063 else if (shape_base[0] == 2) {
2064 switch (type) {
2065 case(TYPE_FLOAT):
2066 INCREMENT_DATA_2D(float, data_in_f);
2067 break;
2068 case(TYPE_INT):
2069 INCREMENT_DATA_2D(int, data_in_i);
2070 break;
2071 case(TYPE_BYTE):
2072 INCREMENT_DATA_2D(byte, data_in_b);
2073 break;
2074 case(TYPE_SHORT):
2075 INCREMENT_DATA_2D(short, data_in_s);
2076 break;
2077 case(TYPE_DOUBLE):
2078 INCREMENT_DATA_2D(double, data_in_d);
2079 break;
2080 case(TYPE_UBYTE):
2081 INCREMENT_DATA_2D(ubyte, data_in_ub);
2082 break;
2083 case(TYPE_UINT):
2084 INCREMENT_DATA_2D(uint, data_in_ui);
2085 break;
2086 case(TYPE_USHORT):
2087 INCREMENT_DATA_2D(ushort, data_in_us);
2088 break;
2089 default:
2090 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2091 goto error;
2092 }
2093 }
2094 else if (shape_base[0] ==3) {
2095 switch (type) {
2096 case(TYPE_FLOAT):
2097 INCREMENT_DATA_3D(float, data_in_f);
2098 break;
2099 case(TYPE_INT):
2100 INCREMENT_DATA_3D(int, data_in_i);
2101 break;
2102 case(TYPE_BYTE):
2103 INCREMENT_DATA_3D(byte, data_in_b);
2104 break;
2105 case(TYPE_SHORT):
2106 INCREMENT_DATA_3D(short, data_in_s);
2107 break;
2108 case(TYPE_DOUBLE):
2109 INCREMENT_DATA_3D(double, data_in_d);
2110 break;
2111 case(TYPE_UBYTE):
2112 INCREMENT_DATA_3D(ubyte, data_in_ub);
2113 break;
2114 case(TYPE_UINT):
2115 INCREMENT_DATA_3D(uint, data_in_ui);
2116 break;
2117 case(TYPE_USHORT):
2118 INCREMENT_DATA_3D(ushort, data_in_us);
2119 break;
2120 default:
2121 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2122 goto error;
2123 }
2124 }
2125 else {
2126 DXSetError(ERROR_DATA_INVALID,
2127 "only 1D, 2D, and 3D input positions accepted");
2128 goto error;
2129 }
2130
2131
2132 /* now fill the actual data array (int or float or whatever)
2133 with the computed sums */
2134 switch (type) {
2135 case (TYPE_FLOAT):
2136 for (i=0; i<numitemsbase; i++) {
2137 if (DXIsElementValid(out_invalidhandle, i))
2138 for (c=0; c<shape[0]; c++)
2139 data_f[i*shape[0]+c] = data_ptr[i*shape[0]+c]/weights_ptr[i];
2140 else {
2141 /* XXX here do something about missing */
2142 if (missing)
2143 for (c=0; c<shape[0]; c++)
2144 data_f[i*shape[0]+c] = missingval_f[c];
2145 else
2146 for (c=0; c<shape[0]; c++)
2147 data_f[i*shape[0]+c] = 0;
2148 }
2149 }
2150 break;
2151 case (TYPE_INT):
2152 for (i=0; i<numitemsbase; i++) {
2153 if (DXIsElementValid(out_invalidhandle, i))
2154 for (c=0; c<shape[0]; c++)
2155 data_i[i*shape[0]+c] =
2156 (int)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2157 else {
2158 if (missing)
2159 for (c=0; c<shape[0]; c++)
2160 data_i[i*shape[0]+c] = (int)missingval_i[c];
2161 else
2162 for (c=0; c<shape[0]; c++)
2163 data_i[i*shape[0]+c] = (int)0;
2164 }
2165 }
2166 break;
2167 case (TYPE_SHORT):
2168 for (i=0; i<numitemsbase; i++) {
2169 if (DXIsElementValid(out_invalidhandle, i))
2170 for (c=0; c<shape[0]; c++)
2171 data_s[i*shape[0]+c]
2172 = (short)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2173 else {
2174 if (missing)
2175 for (c=0; c<shape[0]; c++)
2176 data_s[i*shape[0]+c] =(short)missingval_s[c];
2177 else
2178 for (c=0; c<shape[0]; c++)
2179 data_s[i*shape[0]+c] =(short)0;
2180 }
2181 }
2182 break;
2183 case (TYPE_BYTE):
2184 for (i=0; i<numitemsbase; i++) {
2185 if (DXIsElementValid(out_invalidhandle, i))
2186 for (c=0; c<shape[0]; c++)
2187 data_b[i*shape[0]+c] =
2188 (byte)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2189 else {
2190 if (missing)
2191 for (c=0; c<shape[0]; c++)
2192 data_b[i*shape[0]+c] = (byte)missingval_b[c];
2193 else
2194 for (c=0; c<shape[0]; c++)
2195 data_b[i*shape[0]+c] = (byte)0;
2196 }
2197 }
2198 break;
2199 case (TYPE_DOUBLE):
2200 for (i=0; i<numitemsbase; i++) {
2201 if (DXIsElementValid(out_invalidhandle, i))
2202 for (c=0; c<shape[0]; c++)
2203 data_d[i*shape[0]+c] =
2204 (double)(data_ptr[i*shape[0]+c]/weights_ptr[i]);
2205 else {
2206 if (missing)
2207 for (c=0; c<shape[0]; c++)
2208 data_d[i*shape[0]+c] = (double)missingval_d[c];
2209 else
2210 for (c=0; c<shape[0]; c++)
2211 data_d[i*shape[0]+c] = (double)0;
2212 }
2213 }
2214 break;
2215 case (TYPE_UINT):
2216 for (i=0; i<numitemsbase; i++) {
2217 if (DXIsElementValid(out_invalidhandle, i))
2218 for (c=0; c<shape[0]; c++)
2219 data_ui[i*shape[0]+c] =
2220 (uint)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2221 else {
2222 if (missing)
2223 for (c=0; c<shape[0]; c++)
2224 data_ui[i*shape[0]+c] = (uint)missingval_ui[c];
2225 else
2226 for (c=0; c<shape[0]; c++)
2227 data_ui[i*shape[0]+c] = (uint)0;
2228 }
2229 }
2230 break;
2231 case (TYPE_USHORT):
2232 for (i=0; i<numitemsbase; i++) {
2233 if (DXIsElementValid(out_invalidhandle, i))
2234 for (c=0; c<shape[0]; c++)
2235 data_us[i*shape[0]+c] =
2236 (ushort)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2237 else {
2238 if (missing)
2239 for (c=0; c<shape[0]; c++)
2240 data_us[i*shape[0]+c] = (ushort)missingval_us[c];
2241 else
2242 for (c=0; c<shape[0]; c++)
2243 data_us[i*shape[0]+c] = (ushort)0;
2244 }
2245 }
2246 break;
2247 case (TYPE_UBYTE):
2248 for (i=0; i<numitemsbase; i++) {
2249 if (DXIsElementValid(out_invalidhandle, i))
2250 for (c=0; c<shape[0]; c++)
2251 data_ub[i*shape[0]+c] =
2252 (ubyte)(RND(data_ptr[i*shape[0]+c]/weights_ptr[i]));
2253 else {
2254 if (missing)
2255 for (c=0; c<shape[0]; c++)
2256 data_ub[i*shape[0]+c] = (ubyte)missingval_ub[c];
2257 else
2258 for (c=0; c<shape[0]; c++)
2259 data_ub[i*shape[0]+c] = (ubyte)0;
2260 }
2261 }
2262 break;
2263 default:
2264 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2265 goto error;
2266 }
2267
2268 if (!DXSetComponentValue((Field)base,name,(Object)newcomponent))
2269 goto error;
2270 if (!DXSetComponentAttribute((Field)base, name,
2271 "dep",
2272 (Object)DXNewString("positions")))
2273 goto error;
2274 newcomponent = NULL;
2275 if (!DXChangedComponentValues((Field)base,name))
2276 goto error;
2277 component_done:
2278 compcount++;
2279 DXFree((Pointer)missingval_f);
2280 missingval_f=NULL;
2281 DXFree((Pointer)missingval_b);
2282 missingval_b=NULL;
2283 DXFree((Pointer)missingval_i);
2284 missingval_i=NULL;
2285 DXFree((Pointer)missingval_s);
2286 missingval_s=NULL;
2287 DXFree((Pointer)missingval_d);
2288 missingval_d=NULL;
2289 DXFree((Pointer)missingval_ub);
2290 missingval_ub=NULL;
2291 DXFree((Pointer)missingval_us);
2292 missingval_us=NULL;
2293 DXFree((Pointer)missingval_ui);
2294 missingval_ui=NULL;
2295 }
2296
2297
2298 if (componentsdone==0) {
2299 if (missing) {
2300 DXWarning("no position-dependent components found; cannot use missing value provided");
2301 }
2302 else {
2303 switch (shape_base[0]) {
2304 case (1):
2305 for (i=0; i<numitemsino; i++) {
2306 xpos = pos_ino_ptr[i];
2307 xmin = xpos-radius;
2308 xmax = xpos+radius;
2309 min = (xmin/matrix.A[0][0])-(matrix.b[0]/matrix.A[0][0]);
2310 max = (xmax/matrix.A[0][0])-(matrix.b[0]/matrix.A[0][0]);
2311 /* minvec and maxvec are the indices corresponding to the
2312 bounding box around the sample point */
2313 for (j=CEIL(min); j <= FLOOR(max) ; j++)
2314 if (j>=0 && j<counts[0]) {
2315 /* here are the only ones which are in the bounding
2316 * box */
2317 pos = j*matrix.A[0][0]+matrix.b[0];
2318 distance = pos-xpos;
2319 if (distance <= radius) {
2320 DXSetElementValid(out_invalidhandle, j);
2321 }
2322 }
2323 }
2324 break;
2325 case (2): {
2326 Matrix2D tmpmat;
2327 matrix2D.A[0][0] = matrix.A[0][0];
2328 matrix2D.A[1][0] = matrix.A[0][1];
2329 matrix2D.A[0][1] = matrix.A[0][2];
2330 matrix2D.A[1][1] = matrix.A[1][0];
2331 matrix2D.b[0] = matrix.b[0];
2332 matrix2D.b[1] = matrix.b[1];
2333 invmatrix2D = Invert2D(matrix2D);
2334 /* loop over all scattered points */
2335 for (i=0; i<numitemsino; i++) {
2336 xpos = pos_ino_ptr[2*i];
2337 ypos = pos_ino_ptr[2*i+1];
2338 xmin = xpos-radius;
2339 xmax = xpos+radius;
2340 ymin = ypos-radius;
2341 ymax = ypos+radius;
2342 tmpmat = invmatrix2D;
2343 minvec2D = Apply2D(Vec2D(xmin,ymin), tmpmat);
2344 maxvec2D = Apply2D(Vec2D(xmax,ymax), tmpmat);
2345 /* minvec and maxvec are the indices corresponding to
2346 the bounding box around the sample point */
2347 tmpmat = matrix2D;
2348 for (j=CEIL(minvec2D.x); j <= FLOOR(maxvec2D.x) ; j++)
2349 if (j>=0 && j<counts[0])
2350 for (k=CEIL(minvec2D.y); k <= FLOOR(maxvec2D.y) ; k++)
2351 if (k>=0 && k<counts[1]) {
2352 /* here are the only ones which are in the
2353 bounding box */
2354 posvec2D = Apply2D(Vec2D(j,k),tmpmat);
2355 distance = POW(posvec2D.x-xpos, 2) +
2356 POW(posvec2D.y-ypos, 2);
2357 distance = sqrt(distance);
2358 if (distance <= radius) {
2359 cntr = j*counts[1]+k;
2360 DXSetElementValid(out_invalidhandle, cntr);
2361 }
2362 }
2363 }
2364 break;
2365 }
2366 case (3): {
2367 Matrix tmpmat;
2368 invmatrix = DXInvert(matrix);
2369 /* loop over all scattered points */
2370 for (i=0; i<numitemsino; i++) {
2371 xpos = pos_ino_ptr[3*i];
2372 ypos = pos_ino_ptr[3*i+1];
2373 zpos = pos_ino_ptr[3*i+2];
2374 xmin = xpos-radius;
2375 xmax = xpos+radius;
2376 ymin = ypos-radius;
2377 ymax = ypos+radius;
2378 zmin = zpos-radius;
2379 zmax = zpos+radius;
2380 tmpmat=invmatrix;
2381 minvec = DXApply(DXVec(xmin,ymin,zmin), tmpmat);
2382 maxvec = DXApply(DXVec(xmax,ymax,zmax), tmpmat);
2383 /* minvec and maxvec are the indices corresponding to the */
2384 /* bounding box around the sample point */
2385 tmpmat=matrix;
2386 for (j=CEIL(minvec.x); j <= FLOOR(maxvec.x) ; j++)
2387 if (j>=0 && j<counts[0])
2388 for (k=CEIL(minvec.y); k <= FLOOR(maxvec.y) ; k++)
2389 if (k>=0 && k<counts[1])
2390 for (l=CEIL(minvec.z); l <= FLOOR(maxvec.z) ; l++)
2391 if (l>=0 && l<counts[2]) {
2392 /* here are the only ones which are in the
2393 bounding box */
2394 posvec = DXApply(DXVec(j,k,l),tmpmat);
2395 distance = POW(posvec.x-xpos, 2) +
2396 POW(posvec.y-ypos, 2) +
2397 POW(posvec.z-zpos, 2);
2398 distance = sqrt(distance);
2399 if (distance <= radius) {
2400 cntr = j*counts[2]*counts[1] +
2401 k*counts[2] + l;
2402 DXSetElementValid(out_invalidhandle, cntr);
2403 }
2404 }
2405 }
2406 break;
2407 }
2408 default:
2409 DXSetError(ERROR_DATA_INVALID,"grid must be 1-, 2-, or 3-dimensional");
2410 goto error;
2411 break;
2412 }
2413 }
2414 }
2415
2416
2417
2418 if (!missing) {
2419 if (!DXSaveInvalidComponent(base,out_invalidhandle))
2420 goto error;
2421 }
2422
2423 DXFreeInvalidComponentHandle(out_invalidhandle);
2424 DXFreeInvalidComponentHandle(in_invalidhandle);
2425 DXFree((Pointer)weights_ptr);
2426 return OK;
2427
2428 error:
2429 DXFreeInvalidComponentHandle(out_invalidhandle);
2430 DXFreeInvalidComponentHandle(in_invalidhandle);
2431 DXDelete((Object)newcomponent);
2432 DXFree((Pointer)weights_ptr);
2433 DXFree((Pointer)missingval_f);
2434 DXFree((Pointer)missingval_b);
2435 DXFree((Pointer)missingval_i);
2436 DXFree((Pointer)missingval_s);
2437 DXFree((Pointer)missingval_d);
2438 DXFree((Pointer)missingval_ub);
2439 DXFree((Pointer)missingval_us);
2440 DXFree((Pointer)missingval_ui);
2441 return ERROR;
2442
2443 }
2444
2445
FillDataArray(Type type,Array oldcomponent,Array newcomponent,ubyte * flagbuf,float * distancebuf,int numitemsbase,int numitemsino,int shape,float exponent,Array missing,char * name)2446 static Error FillDataArray(Type type, Array oldcomponent,
2447 Array newcomponent,
2448 ubyte *flagbuf, float *distancebuf,
2449 int numitemsbase, int numitemsino, int shape,
2450 float exponent, Array missing, char *name)
2451 {
2452 float weight, distance;
2453 int i, j, k;
2454 float *buffer=NULL;
2455 float *d_in_f=NULL, *d_f=NULL, *missingval_f=NULL;
2456 int *d_in_i=NULL, *d_i=NULL, *missingval_i=NULL;
2457 short *d_in_s=NULL, *d_s=NULL, *missingval_s=NULL;
2458 double *d_in_d=NULL, *d_d=NULL, *missingval_d=NULL;
2459 byte *d_in_b=NULL, *d_b=NULL, *missingval_b=NULL;
2460 uint *d_in_ui=NULL, *d_ui=NULL, *missingval_ui=NULL;
2461 ushort *d_in_us=NULL, *d_us=NULL, *missingval_us=NULL;
2462 ubyte *d_in_ub=NULL, *d_ub=NULL, *missingval_ub=NULL;
2463
2464
2465 buffer = (float *)DXAllocateLocal(shape*sizeof(float));
2466 if (!buffer) goto error;
2467
2468 /* check that the missing component, if present, matches the component
2469 we're working on */
2470 if (missing)
2471 switch (type) {
2472 case (TYPE_FLOAT):
2473 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), float,
2474 type, shape, missingval_f);
2475 break;
2476 case (TYPE_INT):
2477 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), int,
2478 type, shape, missingval_i);
2479 break;
2480 case (TYPE_BYTE):
2481 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), byte,
2482 type, shape, missingval_b);
2483 break;
2484 case (TYPE_DOUBLE):
2485 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), double,
2486 type, shape, missingval_d);
2487 break;
2488 case (TYPE_SHORT):
2489 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), short,
2490 type, shape, missingval_s);
2491 break;
2492 case (TYPE_UBYTE):
2493 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), ubyte,
2494 type, shape, missingval_ub);
2495 break;
2496 case (TYPE_USHORT):
2497 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), ushort,
2498 type, shape, missingval_us);
2499 break;
2500 case (TYPE_UINT):
2501 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent), uint,
2502 type, shape, missingval_ui);
2503 break;
2504 default:
2505 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2506 goto error;
2507 }
2508
2509
2510 switch (type) {
2511 case (TYPE_FLOAT):
2512 d_in_f = (float *)DXGetArrayData(oldcomponent);
2513 d_f = (float *)DXGetArrayData(newcomponent);
2514 break;
2515 case (TYPE_INT):
2516 d_in_i = (int *)DXGetArrayData(oldcomponent);
2517 d_i = (int *)DXGetArrayData(newcomponent);
2518 break;
2519 case (TYPE_SHORT):
2520 d_in_s = (short *)DXGetArrayData(oldcomponent);
2521 d_s = (short *)DXGetArrayData(newcomponent);
2522 break;
2523 case (TYPE_DOUBLE):
2524 d_in_d = (double *)DXGetArrayData(oldcomponent);
2525 d_d = (double *)DXGetArrayData(newcomponent);
2526 break;
2527 case (TYPE_BYTE):
2528 d_in_b = (byte *)DXGetArrayData(oldcomponent);
2529 d_b = (byte *)DXGetArrayData(newcomponent);
2530 break;
2531 case (TYPE_UBYTE):
2532 d_in_ub = (ubyte *)DXGetArrayData(oldcomponent);
2533 d_ub = (ubyte *)DXGetArrayData(newcomponent);
2534 break;
2535 case (TYPE_USHORT):
2536 d_in_us = (ushort *)DXGetArrayData(oldcomponent);
2537 d_us = (ushort *)DXGetArrayData(newcomponent);
2538 break;
2539 case (TYPE_UINT):
2540 d_in_ui = (uint *)DXGetArrayData(oldcomponent);
2541 d_ui = (uint *)DXGetArrayData(newcomponent);
2542 break;
2543 default:
2544 DXSetError(ERROR_DATA_INVALID,"unrecognized data type");
2545 goto error;
2546 }
2547
2548
2549 for (i=0; i<numitemsbase; i++) {
2550 weight = 0;
2551 for (k=0; k<shape; k++) {
2552 buffer[k] = 0.0;
2553 }
2554 for (j=0; j<numitemsino; j++) {
2555 if (flagbuf[i*numitemsino+j]) {
2556 distance = distancebuf[i*numitemsino+j];
2557 if (distance != 0) {
2558 for (k=0; k<shape; k++) {
2559 switch (type) {
2560 case (TYPE_FLOAT):
2561 buffer[k] += (float)d_in_f[j*shape+k]/POW(distance,exponent);
2562 break;
2563 case (TYPE_INT):
2564 buffer[k] += (float)d_in_i[j*shape+k]/POW(distance,exponent);
2565 break;
2566 case (TYPE_SHORT):
2567 buffer[k] += (float)d_in_s[j*shape+k]/POW(distance,exponent);
2568 break;
2569 case (TYPE_BYTE):
2570 buffer[k] += (float)d_in_b[j*shape+k]/POW(distance,exponent);
2571 break;
2572 case (TYPE_DOUBLE):
2573 buffer[k] += (float)d_in_d[j*shape+k]/POW(distance,exponent);
2574 break;
2575 case (TYPE_UBYTE):
2576 buffer[k] += (float)d_in_ub[j*shape+k]/POW(distance,exponent);
2577 break;
2578 case (TYPE_USHORT):
2579 buffer[k] += (float)d_in_us[j*shape+k]/POW(distance,exponent);
2580 break;
2581 case (TYPE_UINT):
2582 buffer[k] += (float)d_in_ui[j*shape+k]/POW(distance,exponent);
2583 break;
2584 default:
2585 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2586 goto error;
2587 }
2588 }
2589 weight += 1/POW(distance,exponent);
2590 }
2591 else {
2592 /* if distance == 0, then simply use the data value */
2593 for (k=0; k<shape; k++) {
2594 switch (type) {
2595 case (TYPE_FLOAT):
2596 buffer[k] = (float)d_in_f[j];
2597 break;
2598 case (TYPE_INT):
2599 buffer[k] = (float)d_in_i[j];
2600 break;
2601 case (TYPE_BYTE):
2602 buffer[k] = (float)d_in_b[j];
2603 break;
2604 case (TYPE_DOUBLE):
2605 buffer[k] = (float)d_in_d[j];
2606 break;
2607 case (TYPE_SHORT):
2608 buffer[k] = (float)d_in_s[j];
2609 break;
2610 case (TYPE_UINT):
2611 buffer[k] = (float)d_in_ui[j];
2612 break;
2613 case (TYPE_USHORT):
2614 buffer[k] = (float)d_in_us[j];
2615 break;
2616 case (TYPE_UBYTE):
2617 buffer[k] = (float)d_in_us[j];
2618 break;
2619 default:
2620 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2621 goto error;
2622 }
2623 }
2624 weight = 1;
2625 continue;
2626 }
2627 }
2628 }
2629 if (weight > 0) {
2630 for (k=0; k<shape; k++) {
2631 switch (type) {
2632 case (TYPE_FLOAT):
2633 d_f[shape*i+k] = (float)(buffer[k]/weight);
2634 break;
2635 case (TYPE_DOUBLE):
2636 d_d[shape*i+k] = (double)(buffer[k]/weight);
2637 break;
2638 case (TYPE_INT):
2639 d_i[shape*i+k] = (int)(RND(buffer[k]/weight));
2640 break;
2641 case (TYPE_SHORT):
2642 d_s[shape*i+k] = (short)(RND(buffer[k]/weight));
2643 break;
2644 case (TYPE_BYTE):
2645 d_b[shape*i+k] = (byte)(RND(buffer[k]/weight));
2646 break;
2647 case (TYPE_UBYTE):
2648 d_ub[shape*i+k] = (ubyte)(RND(buffer[k]/weight));
2649 break;
2650 case (TYPE_UINT):
2651 d_ui[shape*i+k] = (uint)(RND(buffer[k]/weight));
2652 break;
2653 case (TYPE_USHORT):
2654 d_us[shape*i+k] = (ushort)(RND(buffer[k]/weight));
2655 break;
2656 default:
2657 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2658 goto error;
2659 }
2660 }
2661 }
2662 else {
2663 if (missing)
2664 for (k=0; k<shape; k++) {
2665 switch (type) {
2666 case (TYPE_FLOAT):
2667 d_f[i*shape+k] = (float)missingval_f[k];
2668 break;
2669 case (TYPE_INT):
2670 d_i[i*shape+k] = (int)missingval_i[k];
2671 break;
2672 case (TYPE_BYTE):
2673 d_b[i*shape+k] = (byte)missingval_b[k];
2674 break;
2675 case (TYPE_DOUBLE):
2676 d_d[i*shape+k] = (double)missingval_d[k];
2677 break;
2678 case (TYPE_SHORT):
2679 d_s[i*shape+k] = (short)missingval_s[k];
2680 break;
2681 case (TYPE_UBYTE):
2682 d_ub[i*shape+k] = (ubyte)missingval_ub[k];
2683 break;
2684 case (TYPE_UINT):
2685 d_ui[i*shape+k] = (uint)missingval_ui[k];
2686 break;
2687 case (TYPE_USHORT):
2688 d_us[i*shape+k] = (ushort)missingval_us[k];
2689 break;
2690 default:
2691 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2692 goto error;
2693 }
2694 }
2695 else
2696 for (k=0; k<shape; k++) {
2697 switch (type) {
2698 case (TYPE_FLOAT):
2699 d_f[i*shape+k] = 0.0;
2700 break;
2701 case (TYPE_INT):
2702 d_i[i*shape+k] = 0.0;
2703 break;
2704 case (TYPE_BYTE):
2705 d_b[i*shape+k] = 0.0;
2706 break;
2707 case (TYPE_DOUBLE):
2708 d_d[i*shape+k] = 0.0;
2709 break;
2710 case (TYPE_SHORT):
2711 d_s[i*shape+k] = 0.0;
2712 break;
2713 case (TYPE_UBYTE):
2714 d_ub[i*shape+k] = 0.0;
2715 break;
2716 case (TYPE_UINT):
2717 d_ui[i*shape+k] = 0.0;
2718 break;
2719 case (TYPE_USHORT):
2720 d_us[i*shape+k] = 0.0;
2721 break;
2722 default:
2723 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
2724 goto error;
2725 }
2726 }
2727 }
2728 }
2729 DXFree((Pointer)missingval_f);
2730 missingval_f=NULL;
2731 DXFree((Pointer)missingval_i);
2732 missingval_i=NULL;
2733 DXFree((Pointer)missingval_d);
2734 missingval_d=NULL;
2735 DXFree((Pointer)missingval_s);
2736 missingval_s=NULL;
2737 DXFree((Pointer)missingval_b);
2738 missingval_b=NULL;
2739 DXFree((Pointer)missingval_ub);
2740 missingval_ub=NULL;
2741 DXFree((Pointer)missingval_ui);
2742 missingval_ui=NULL;
2743 DXFree((Pointer)missingval_us);
2744 missingval_us=NULL;
2745 return OK;
2746 error:
2747 DXFree ((Pointer)buffer);
2748 DXFree((Pointer)missingval_f);
2749 DXFree((Pointer)missingval_i);
2750 DXFree((Pointer)missingval_d);
2751 DXFree((Pointer)missingval_s);
2752 DXFree((Pointer)missingval_b);
2753 DXFree((Pointer)missingval_ub);
2754 DXFree((Pointer)missingval_ui);
2755 DXFree((Pointer)missingval_us);
2756 return ERROR;
2757 }
2758
2759
2760
2761
2762
2763
FillBuffersIrreg(int numitemsbase,int numitemsino,ubyte * flagbuf,float * distancebuf,float * pos_ino_ptr,float * pos_base_ptr,InvalidComponentHandle out_invalidhandle,InvalidComponentHandle in_invalidhandle,int * anyinvalid,int shape_base,float radius)2764 static Error FillBuffersIrreg(int numitemsbase, int numitemsino,
2765 ubyte *flagbuf,
2766 float *distancebuf, float *pos_ino_ptr,
2767 float *pos_base_ptr,
2768 InvalidComponentHandle out_invalidhandle,
2769 InvalidComponentHandle in_invalidhandle,
2770 int *anyinvalid, int shape_base, float radius)
2771 {
2772 int i, j, anyvalid;
2773 float distance=0;
2774
2775 /* there's a special check in here for radius == -1 (infinity) */
2776
2777
2778 /* fill the buffers */
2779 if (shape_base == 1) {
2780 for (i=0; i<numitemsbase; i++) {
2781 anyvalid = 0;
2782 for (j=0; j<numitemsino;j++) {
2783 if (DXIsElementValid(in_invalidhandle, j)) {
2784 if (radius != -1) {
2785 distance = ABS(pos_ino_ptr[j]-pos_base_ptr[i]);
2786 if (distance <= radius) {
2787 anyvalid = 1;
2788 flagbuf[i*numitemsino + j]=1;
2789 distancebuf[i*numitemsino + j]=distance;
2790 }
2791 }
2792 else {
2793 anyvalid = 1;
2794 flagbuf[i*numitemsino + j]=1;
2795 distancebuf[i*numitemsino + j]=distance;
2796 }
2797 }
2798 }
2799 if (!anyvalid) {
2800 if (out_invalidhandle)
2801 DXSetElementInvalid(out_invalidhandle, i);
2802 *anyinvalid = 1;
2803 }
2804 }
2805 }
2806
2807 else if (shape_base == 2) {
2808 for (i=0; i<numitemsbase; i++) {
2809 anyvalid = 0;
2810 for (j=0; j<numitemsino;j++) {
2811 if (DXIsElementValid(in_invalidhandle, j)) {
2812 if (radius != -1) {
2813 if (((pos_ino_ptr[2*j]-pos_base_ptr[2*i])<radius)&&
2814 ((pos_ino_ptr[2*j+1]-pos_base_ptr[2*i+1])<radius)) {
2815 distance = POW(pos_ino_ptr[2*j]-pos_base_ptr[2*i], 2) +
2816 POW(pos_ino_ptr[2*j+1]-pos_base_ptr[2*i+1], 2);
2817 distance = sqrt(distance);
2818 if (distance <= radius) {
2819 anyvalid = 1;
2820 flagbuf[i*numitemsino + j]=1;
2821 distancebuf[i*numitemsino + j]=distance;
2822 }
2823 }
2824 }
2825 else {
2826 anyvalid = 1;
2827 flagbuf[i*numitemsino + j]=1;
2828 distancebuf[i*numitemsino + j]=distance;
2829 }
2830 }
2831 }
2832 if (!anyvalid) {
2833 if (out_invalidhandle)
2834 DXSetElementInvalid(out_invalidhandle, i);
2835 *anyinvalid = 1;
2836 }
2837 }
2838 }
2839
2840
2841 else if (shape_base==3) {
2842 for (i=0; i<numitemsbase; i++) {
2843 anyvalid = 0;
2844 for (j=0; j<numitemsino;j++) {
2845 if (DXIsElementValid(in_invalidhandle, j)) {
2846 if (radius != -1) {
2847 if (((pos_ino_ptr[3*j]-pos_base_ptr[3*i])<radius)&&
2848 ((pos_ino_ptr[3*j+1]-pos_base_ptr[3*i+1])<radius) &&
2849 ((pos_ino_ptr[3*j+2]-pos_base_ptr[3*i+2])<radius)) {
2850 distance = POW(pos_ino_ptr[3*j]-pos_base_ptr[3*i], 2) +
2851 POW(pos_ino_ptr[3*j+1]-pos_base_ptr[3*i+1], 2) +
2852 POW(pos_ino_ptr[3*j+2]-pos_base_ptr[3*i+2], 2);
2853 distance = sqrt(distance);
2854 if (distance <= radius) {
2855 anyvalid = 1;
2856 flagbuf[i*numitemsino + j]=1;
2857 distancebuf[i*numitemsino + j]=distance;
2858 }
2859 }
2860 }
2861 else {
2862 anyvalid = 1;
2863 flagbuf[i*numitemsino + j]=1;
2864 distancebuf[i*numitemsino + j]=distance;
2865 }
2866 }
2867 }
2868 if (!anyvalid) {
2869 if (out_invalidhandle)
2870 DXSetElementInvalid(out_invalidhandle, 1);
2871 *anyinvalid = 1;
2872 }
2873 }
2874 }
2875 else {
2876 DXSetError(ERROR_BAD_PARAMETER,"only 1D 2D 3D positions accepted");
2877 goto error;
2878 }
2879 return OK;
2880 error:
2881 return ERROR;
2882 }
2883
2884
2885 /*Added by Jeff Braun to put scattered data on closest grid point*/
2886
_dxfConnectScatterObject(Object ino,Object base,Array missing)2887 Error _dxfConnectScatterObject(Object ino, Object base, Array missing)
2888 {
2889 struct arg_scatter arg;
2890 Object subbase;
2891 int i;
2892
2893 switch (DXGetObjectClass(base)) {
2894 case CLASS_FIELD:
2895 /* create the task group */
2896 arg.ino = ino;
2897 arg.base = (Field)base;
2898 arg.missing = missing;
2899 if (!DXAddTask(ConnectScatterField,(Pointer)&arg, sizeof(arg),0.0))
2900 goto error;
2901 break;
2902 case CLASS_GROUP:
2903 for (i=0; (subbase=DXGetEnumeratedMember((Group)base, i, NULL)); i++) {
2904 if (!(_dxfConnectScatterObject((Object)ino, (Object)subbase, missing)))
2905 goto error;
2906 }
2907 break;
2908 default:
2909 DXSetError(ERROR_DATA_INVALID,"base must be a field or a group");
2910 goto error;
2911 }
2912 return OK;
2913
2914 error:
2915 return ERROR;
2916
2917 }
2918
ConnectScatterField(Pointer p)2919 static Error ConnectScatterField(Pointer p)
2920 {
2921 struct arg_scatter *arg = (struct arg_scatter *)p;
2922 Object ino, newino=NULL;
2923 Field base;
2924 Array missing;
2925
2926 ino = arg->ino;
2927 newino = DXApplyTransform(ino,NULL);
2928
2929 switch (DXGetObjectClass(newino)) {
2930
2931 case CLASS_FIELD:
2932 base = arg->base;
2933 missing = arg->missing;
2934 if (!ConnectScatter((Field)newino, base, missing))
2935 goto error;
2936 break;
2937
2938 case CLASS_GROUP:
2939 DXSetError(ERROR_NOT_IMPLEMENTED,"input cannot be a group");
2940 goto error;
2941 break;
2942
2943 default:
2944 DXSetError(ERROR_DATA_INVALID,"input must be a field");
2945 goto error;
2946 }
2947 DXDelete((Object)newino);
2948 return OK;
2949
2950 error:
2951 DXDelete((Object)newino);
2952 return ERROR;
2953
2954 }
2955
ConnectScatter(Field ino,Field base,Array missing)2956 static Error ConnectScatter(Field ino, Field base, Array missing)
2957 {
2958 Matrix matrix, inverse;
2959 Matrix2D matrix2D, inverse2D;
2960 Type type;
2961 Category category;
2962 Array positions_base, positions_ino, oldcomponent=NULL, newcomponent;
2963 InvalidComponentHandle out_invalidhandle=NULL, in_invalidhandle=NULL;
2964
2965 int rank, shape[8];
2966 int counts[8], rank_base, shape_base[8], shape_ino[8], rank_ino;
2967 int nitemsbase, nitemsino, nitems, compcount, componentsdone;
2968 float *ino_positions;
2969 char *name;
2970
2971 int i,j,k, ijk[3], index, *count=NULL;
2972 float xyz[3], sum, eps=.000001;
2973
2974
2975 float *old_data_f=NULL, *new_data_ptr_f=NULL, *missingval_f=NULL;
2976 int *old_data_i=NULL, *new_data_ptr_i=NULL, *missingval_i=NULL;
2977 short *old_data_s=NULL, *new_data_ptr_s=NULL, *missingval_s=NULL;
2978 double *old_data_d=NULL, *new_data_ptr_d=NULL, *missingval_d=NULL;
2979 byte *old_data_b=NULL, *new_data_ptr_b=NULL, *missingval_b=NULL;
2980 ubyte *old_data_ub=NULL, *new_data_ptr_ub=NULL, *missingval_ub=NULL;
2981 ushort *old_data_us=NULL, *new_data_ptr_us=NULL, *missingval_us=NULL;
2982 uint *old_data_ui=NULL, *new_data_ptr_ui=NULL, *missingval_ui=NULL;
2983
2984 if (DXEmptyField((Field)ino))
2985 return OK;
2986
2987 positions_base = (Array)DXGetComponentValue((Field)base, "positions");
2988 if (!positions_base) {
2989 DXSetError(ERROR_MISSING_DATA,"base field has no positions");
2990 goto error;
2991 }
2992 positions_ino = (Array)DXGetComponentValue((Field)ino, "positions");
2993 if (!positions_ino) {
2994 DXSetError(ERROR_MISSING_DATA,"input field has no positions");
2995 goto error;
2996 }
2997
2998 /*Grid Position Data*/
2999 DXGetArrayInfo(positions_base, &nitemsbase, &type, &category,
3000 &rank_base, shape_base);
3001 if ((DXQueryGridPositions((Array)positions_base,NULL,counts,
3002 (float *)&matrix.b,(float *)&matrix.A) )) {
3003 if (shape_base[0]==2) {
3004 matrix2D.A[0][0] = matrix.A[0][0];
3005 matrix2D.A[1][0] = matrix.A[0][1];
3006 matrix2D.A[0][1] = matrix.A[0][2];
3007 matrix2D.A[1][1] = matrix.A[1][0];
3008 matrix2D.b[0] = matrix.b[0];
3009 matrix2D.b[1] = matrix.b[1];
3010 inverse2D=Invert2D(matrix2D);
3011 inverse.A[0][0]=inverse2D.A[0][0];
3012 inverse.A[0][1]=inverse2D.A[0][1];
3013 inverse.A[1][0]=inverse2D.A[1][0];
3014 inverse.A[1][1]=inverse2D.A[1][1];
3015 }
3016 else if (shape_base[0]==1) {
3017 /*origin = matrix.b[0];*/
3018 /*delta = matrix.A[0][0];*/
3019 DXSetError(ERROR_DATA_INVALID, "Only 2D and 3D grids supported");
3020 goto error;
3021 }
3022 }
3023 else {
3024 /* irregular positions */
3025 DXSetError(ERROR_DATA_INVALID, "irregular grids not currently supported for radius=0");
3026 goto error;
3027 }
3028
3029 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
3030 DXSetError(ERROR_DATA_INVALID,
3031 "positions component must be type float, category real");
3032 goto error;
3033 }
3034
3035 /*Scatter Position Data*/
3036 DXGetArrayInfo(positions_ino, &nitemsino, &type, &category, &rank_ino,
3037 shape_ino);
3038
3039 if ((type != TYPE_FLOAT)||(category != CATEGORY_REAL)) {
3040 DXSetError(ERROR_DATA_INVALID,
3041 "positions component must be type float, category real");
3042 goto error;
3043 }
3044
3045 if (rank_base == 0)
3046 shape_base[0] = 1;
3047 if (rank_base > 1) {
3048 DXSetError(ERROR_DATA_INVALID,
3049 "rank of base positions cannot be greater than 1");
3050 goto error;
3051 }
3052 if (rank_ino == 0)
3053 shape_ino[0] = 1;
3054 if (rank_ino > 1) {
3055 DXSetError(ERROR_DATA_INVALID,
3056 "rank of input positions cannot be greater than 1");
3057 goto error;
3058 }
3059
3060 if ((shape_base[0]<0)||(shape_base[0]>3)) {
3061 DXSetError(ERROR_DATA_INVALID,
3062 "only 1D, 2D, and 3D positions for base supported");
3063 goto error;
3064 }
3065
3066 if (shape_base[0] != shape_ino[0]) {
3067 DXSetError(ERROR_DATA_INVALID,
3068 "dimensionality of base (%dD) does not match dimensionality of input (%dD)",
3069 shape_base[0], shape_ino[0]);
3070 goto error;
3071 }
3072
3073 /* first set the extra counts to one to simplify the code for handling
3074 1, 2, and 3 dimensions */
3075
3076 if (shape_base[0]==1) {
3077 counts[1]=1;
3078 counts[2]=1;
3079 }
3080 else if (shape_base[0]==2) {
3081 counts[2]=1;
3082 }
3083
3084 /*Will mark as invalid if 'missing' not specified*/
3085
3086 if (!missing){
3087 out_invalidhandle = DXCreateInvalidComponentHandle((Object)base, NULL,
3088 "positions");
3089 if (!out_invalidhandle)
3090 goto error;
3091 DXSetAllInvalid(out_invalidhandle);
3092 }
3093 in_invalidhandle = DXCreateInvalidComponentHandle((Object)ino, NULL,
3094 "positions");
3095 if (!in_invalidhandle)
3096 return ERROR;
3097
3098 /* need to get the components we are going to copy into the new field */
3099 /* copied portions from ConnectNearest -JAB */
3100 compcount = 0;
3101 componentsdone = 0;
3102 while (NULL !=
3103 (oldcomponent=(Array)DXGetEnumeratedComponentValue(ino, compcount,
3104 &name))) {
3105 Object attr;
3106
3107 if (DXGetObjectClass((Object)oldcomponent) != CLASS_ARRAY)
3108 goto component_done;
3109
3110 if (!strcmp(name,"positions") ||
3111 !strcmp(name,"connections") ||
3112 !strcmp(name,"invalid positions"))
3113 goto component_done;
3114
3115 /* if the component refs anything, leave it */
3116 if (DXGetComponentAttribute((Field)ino,name,"ref"))
3117 goto component_done;
3118
3119 attr = DXGetComponentAttribute((Field)ino,name,"der");
3120 if (attr)
3121 goto component_done;
3122
3123 attr = DXGetComponentAttribute((Field)ino,name,"dep");
3124 if (!attr) {
3125 /* does it der? if so, throw it away */
3126 attr = DXGetComponentAttribute((Field)ino,name,"der");
3127 if (attr)
3128 goto component_done;
3129 /* it doesn't ref, doesn't dep and doesn't der.
3130 Copy it to the output and quit */
3131 if (!DXSetComponentValue((Field)base,name,(Object)oldcomponent))
3132 goto error;
3133 goto component_done;
3134 }
3135
3136 if (DXGetObjectClass(attr) != CLASS_STRING) {
3137 DXSetError(ERROR_DATA_INVALID, "dependency attribute not type string");
3138 goto error;
3139 }
3140
3141 if (strcmp(DXGetString((String)attr),"positions")){
3142 DXWarning("Component '%s' is not dependent on positions!",name);
3143 DXWarning("Regrid will remove '%s' and output the base grid 'data' component if it exists",name);
3144 goto component_done;
3145 }
3146
3147
3148 /* if the component happens to be "invalid positions" ignore it */
3149 if (!strcmp(name,"invalid positions"))
3150 goto component_done;
3151
3152
3153 componentsdone++;
3154
3155 DXGetArrayInfo((Array)oldcomponent,&nitems, &type,
3156 &category, &rank, shape);
3157 if (rank==0) shape[0]=1;
3158
3159 /* check that the missing component, if present, matches the component
3160 we're working on. If missing not present, 0 is used as default */
3161 switch (type) {
3162 case (TYPE_FLOAT):
3163 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3164 float, type, shape[0], missingval_f);
3165 break;
3166 case (TYPE_INT):
3167 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3168 int, type, shape[0], missingval_i);
3169 break;
3170 case (TYPE_DOUBLE):
3171 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3172 double, type, shape[0], missingval_d);
3173 break;
3174 case (TYPE_SHORT):
3175 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3176 short, type, shape[0], missingval_s);
3177 break;
3178 case (TYPE_BYTE):
3179 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3180 byte, type, shape[0], missingval_b);
3181 break;
3182 case (TYPE_UINT):
3183 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3184 uint, type, shape[0], missingval_ui);
3185 break;
3186 case (TYPE_USHORT):
3187 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3188 ushort, type, shape[0], missingval_us);
3189 break;
3190 case (TYPE_UBYTE):
3191 CHECK_AND_ALLOCATE_MISSING(missing, DXGetItemSize(oldcomponent),
3192 ubyte, type, shape[0], missingval_ub);
3193 break;
3194 default:
3195 DXSetError(ERROR_DATA_INVALID,"unsupported data type");
3196 goto error;
3197 }
3198
3199 if (rank > 1) {
3200 DXSetError(ERROR_NOT_IMPLEMENTED,"rank larger than 1 not implemented");
3201 goto error;
3202 }
3203 if (nitems != nitemsino) {
3204 DXSetError(ERROR_BAD_PARAMETER,
3205 "number of items in %s component not equal to number of positions",
3206 name);
3207 goto error;
3208 }
3209
3210
3211 /*Intialize newcomponent with missing value*/
3212 newcomponent = DXNewArrayV(type, category, rank, shape);
3213 if (!DXAddArrayData(newcomponent, 0, nitemsbase, NULL))
3214 goto error;
3215
3216 switch (type) {
3217 case (TYPE_FLOAT):
3218 old_data_f = (float *)DXGetArrayData(oldcomponent);
3219 new_data_ptr_f = (float *)DXGetArrayData(newcomponent);
3220 for(i=0;i<nitemsbase;i++){
3221 for(j=0;j<shape[0];j++){
3222 new_data_ptr_f[i*shape[0]+j] = missingval_f[j];
3223 }
3224 }
3225 break;
3226 case (TYPE_INT):
3227 old_data_i = (int *)DXGetArrayData(oldcomponent);
3228 new_data_ptr_i = (int *)DXGetArrayData(newcomponent);
3229 for(i=0;i<nitemsbase;i++){
3230 for(j=0;j<shape[0];j++){
3231 new_data_ptr_i[i*shape[0]+j] = missingval_i[j];
3232 }
3233 }
3234 break;
3235 case (TYPE_SHORT):
3236 old_data_s = (short *)DXGetArrayData(oldcomponent);
3237 new_data_ptr_s = (short *)DXGetArrayData(newcomponent);
3238 for(i=0;i<nitemsbase;i++){
3239 for(j=0;j<shape[0];j++){
3240 new_data_ptr_s[i*shape[0]+j] = missingval_s[j];
3241 }
3242 }
3243 break;
3244 case (TYPE_BYTE):
3245 old_data_b = (byte *)DXGetArrayData(oldcomponent);
3246 new_data_ptr_b = (byte *)DXGetArrayData(newcomponent);
3247 for(i=0;i<nitemsbase;i++){
3248 for(j=0;j<shape[0];j++){
3249 new_data_ptr_b[i*shape[0]+j] = missingval_b[j];
3250 }
3251 }
3252 break;
3253 case (TYPE_DOUBLE):
3254 old_data_d = (double *)DXGetArrayData(oldcomponent);
3255 new_data_ptr_d = (double *)DXGetArrayData(newcomponent);
3256 for(i=0;i<nitemsbase;i++){
3257 for(j=0;j<shape[0];j++){
3258 new_data_ptr_d[i*shape[0]+j] = missingval_d[j];
3259 }
3260 }
3261 break;
3262 case (TYPE_UINT):
3263 old_data_ui = (uint *)DXGetArrayData(oldcomponent);
3264 new_data_ptr_ui = (uint *)DXGetArrayData(newcomponent);
3265 for(i=0;i<nitemsbase;i++){
3266 for(j=0;j<shape[0];j++){
3267 new_data_ptr_ui[i*shape[0]+j] = missingval_ui[j];
3268 }
3269 }
3270 break;
3271 case (TYPE_UBYTE):
3272 old_data_ub = (ubyte *)DXGetArrayData(oldcomponent);
3273 new_data_ptr_ub = (ubyte *)DXGetArrayData(newcomponent);
3274 for(i=0;i<nitemsbase;i++){
3275 for(j=0;j<shape[0];j++){
3276 new_data_ptr_ub[i*shape[0]+j] = missingval_ub[j];
3277 }
3278 }
3279 break;
3280 case (TYPE_USHORT):
3281 old_data_us = (ushort *)DXGetArrayData(oldcomponent);
3282 new_data_ptr_us = (ushort *)DXGetArrayData(newcomponent);
3283 for(i=0;i<nitemsbase;i++){
3284 for(j=0;j<shape[0];j++){
3285 new_data_ptr_us[i*shape[0]+j] = missingval_us[j];
3286 }
3287 }
3288 break;
3289 default:
3290 DXSetError(ERROR_DATA_INVALID,"unrecognized data type");
3291 goto error;
3292 }
3293
3294
3295 /*Need to allocate count and compute inverse of grid matrix*/
3296 count = (int *) DXAllocateLocalZero(sizeof(int)*nitemsbase);
3297 if (! count)
3298 goto error;
3299 if(shape_base[0]==3){
3300 inverse = DXInvert(matrix);
3301 }
3302
3303 ino_positions = (float *)DXGetArrayData(positions_ino);
3304
3305 /*Put Scattered data into grid data array*/
3306 for(i=0;i<nitemsino;i++){ /*loops through scatter data */
3307 if(DXIsElementValid(in_invalidhandle,i)){ /*seems to work fine without this*/
3308 index=0;
3309 for(j=0;j<shape_base[0];j++){
3310 xyz[j]=ino_positions[shape_base[0]*i+j]-matrix.b[j];
3311 }
3312 for(j=0;j<shape_base[0];j++){
3313 sum=0;
3314 for(k=0;k<shape_base[0];k++){
3315 sum=sum+inverse.A[k][j]*xyz[k];
3316 }
3317 ijk[j]=rint(sum + eps);
3318 if(ijk[j]<0 || ijk[j]>=counts[j]){
3319 index=-1;
3320 }
3321 }
3322
3323 /* Finds data array index number for grid position*/
3324 if(index !=-1){
3325 for(j=0;j<shape_base[0];j++){
3326 for(k=shape_base[0]-1;k>j;k--){
3327 ijk[j]=ijk[j]*counts[k];
3328 }
3329 index=index+ijk[j];
3330 }
3331 }
3332
3333
3334 if(index >= 0 && index <nitemsbase){
3335 if(!missing){
3336 if(!DXSetElementValid(out_invalidhandle, index)) goto error;
3337 }
3338 switch (type){
3339 case(TYPE_FLOAT):
3340 for(j=0;j<shape[0];j++){
3341 if(count[index] == 0)
3342 new_data_ptr_f[shape[0]*index+j]=old_data_f[shape[0]*i+j];
3343 else
3344 new_data_ptr_f[shape[0]*index+j]=
3345 (new_data_ptr_f[shape[0]*index+j]*count[index]
3346 +old_data_f[shape[0]*i+j])/(count[index]+1);
3347 }
3348 break;
3349 case (TYPE_INT):
3350 for(j=0;j<shape[0];j++){
3351 if(count[index] == 0)
3352 new_data_ptr_i[shape[0]*index+j]=old_data_i[shape[0]*i+j];
3353 else
3354 new_data_ptr_i[shape[0]*index+j]=
3355 (new_data_ptr_i[shape[0]*index+j]*count[index]
3356 +old_data_i[shape[0]*i+j])/(count[index]+1);
3357 }
3358 break;
3359 case (TYPE_SHORT):
3360 for(j=0;j<shape[0];j++){
3361 if(count[index] == 0)
3362 new_data_ptr_s[shape[0]*index+j]=old_data_s[shape[0]*i+j];
3363 else
3364 new_data_ptr_s[shape[0]*index+j]=
3365 (new_data_ptr_s[shape[0]*index+j]*count[index]
3366 +old_data_s[shape[0]*i+j])/(count[index]+1);
3367 }
3368 break;
3369 case (TYPE_DOUBLE):
3370 for(j=0;j<shape[0];j++){
3371 if(count[index] == 0)
3372 new_data_ptr_d[shape[0]*index+j]=old_data_d[shape[0]*i+j];
3373 else
3374 new_data_ptr_d[shape[0]*index+j]=
3375 (new_data_ptr_d[shape[0]*index+j]*count[index]
3376 +old_data_d[shape[0]*i+j])/(count[index]+1);
3377 }
3378 break;
3379 case (TYPE_BYTE):
3380 for(j=0;j<shape[0];j++){
3381 if(count[index] == 0)
3382 new_data_ptr_b[shape[0]*index+j]=old_data_b[shape[0]*i+j];
3383 else
3384 new_data_ptr_b[shape[0]*index+j]=
3385 (new_data_ptr_b[shape[0]*index+j]*count[index]
3386 +old_data_b[shape[0]*i+j])/(count[index]+1);
3387 }
3388 break;
3389 case (TYPE_UINT):
3390 for(j=0;j<shape[0];j++){
3391 if(count[index] == 0)
3392 new_data_ptr_ui[shape[0]*index+j]=old_data_ui[shape[0]*i+j];
3393 else
3394 new_data_ptr_ui[shape[0]*index+j]=
3395 (new_data_ptr_ui[shape[0]*index+j]*count[index]
3396 +old_data_ui[shape[0]*i+j])/(count[index]+1);
3397 }
3398 break;
3399 case (TYPE_USHORT):
3400 for(j=0;j<shape[0];j++){
3401 if(count[index] == 0)
3402 new_data_ptr_us[shape[0]*index+j]=old_data_us[shape[0]*i+j];
3403 else
3404 new_data_ptr_us[shape[0]*index+j]=
3405 (new_data_ptr_us[shape[0]*index+j]*count[index]
3406 +old_data_us[shape[0]*i+j])/(count[index]+1);
3407 }
3408 break;
3409 case (TYPE_UBYTE):
3410 for(j=0;j<shape[0];j++){
3411 if(count[index] == 0)
3412 new_data_ptr_ub[shape[0]*index+j]=old_data_ub[shape[0]*i+j];
3413 else
3414 new_data_ptr_ub[shape[0]*index+j]=
3415 (new_data_ptr_ub[shape[0]*index+j]*count[index]
3416 +old_data_ub[shape[0]*i+j])/(count[index]+1);
3417 }
3418 break;
3419 default:
3420 DXSetError(ERROR_DATA_INVALID,"unrecognized data type");
3421 goto error;
3422 }
3423 count[index]++;
3424 }
3425 }
3426 }
3427 if (!DXChangedComponentValues(base,name))
3428 goto error;
3429 if (!DXSetComponentValue(base,name,(Object)newcomponent))
3430 goto error;
3431 newcomponent = NULL;
3432 if (!DXSetComponentAttribute((Field)base, name, "dep",
3433 (Object)DXNewString("positions")))
3434 goto error;
3435
3436 component_done:
3437 compcount++;
3438 DXFree((Pointer)missingval_f);
3439 missingval_f=NULL;
3440 DXFree((Pointer)missingval_s);
3441 missingval_s=NULL;
3442 DXFree((Pointer)missingval_i);
3443 missingval_i=NULL;
3444 DXFree((Pointer)missingval_d);
3445 missingval_d=NULL;
3446 DXFree((Pointer)missingval_b);
3447 missingval_b=NULL;
3448 DXFree((Pointer)missingval_ui);
3449 missingval_ui=NULL;
3450 DXFree((Pointer)missingval_ub);
3451 missingval_ub=NULL;
3452 DXFree((Pointer)missingval_us);
3453 missingval_us=NULL;
3454 DXFree((Pointer)count);
3455 count=NULL;
3456 }
3457
3458 if(!missing){
3459 for(i=0;i<nitemsbase;i++){
3460 if (DXIsElementInvalid(out_invalidhandle,i)) {
3461 i=nitemsbase;
3462 if(!DXSaveInvalidComponent(base, out_invalidhandle)) goto error;
3463 }
3464 }
3465 }
3466 DXFreeInvalidComponentHandle(out_invalidhandle);
3467 DXFreeInvalidComponentHandle(in_invalidhandle);
3468 if(!DXEndField(base)) goto error;
3469 return OK;
3470
3471 error:
3472 DXFreeInvalidComponentHandle(out_invalidhandle);
3473 DXFreeInvalidComponentHandle(in_invalidhandle);
3474 DXFree((Pointer)missingval_f);
3475 DXFree((Pointer)missingval_s);
3476 DXFree((Pointer)missingval_i);
3477 DXFree((Pointer)missingval_d);
3478 DXFree((Pointer)missingval_b);
3479 DXFree((Pointer)missingval_ui);
3480 DXFree((Pointer)missingval_ub);
3481 DXFree((Pointer)missingval_us);
3482 DXFree((Pointer)count);
3483 return ERROR;
3484
3485 }
3486