1 /*********************************************************************
2 dimension -- Functions for multi-dimensional operations.
3 This is part of GNU Astronomy Utilities (Gnuastro) package.
4
5 Original author:
6 Mohammad Akhlaghi <mohammad@akhlaghi.org>
7 Contributing author(s):
8 Copyright (C) 2017-2021, Free Software Foundation, Inc.
9
10 Gnuastro is free software: you can redistribute it and/or modify it
11 under the terms of the GNU General Public License as published by the
12 Free Software Foundation, either version 3 of the License, or (at your
13 option) any later version.
14
15 Gnuastro is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of
17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 General Public License for more details.
19
20 You should have received a copy of the GNU General Public License
21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
22 **********************************************************************/
23 #include <config.h>
24
25 #include <math.h>
26 #include <stdio.h>
27 #include <errno.h>
28 #include <error.h>
29 #include <stdlib.h>
30
31 #include <gnuastro/wcs.h>
32 #include <gnuastro/pointer.h>
33 #include <gnuastro/dimension.h>
34
35
36
37
38
39 /************************************************************************/
40 /******************** Info **********************/
41 /************************************************************************/
42 size_t
gal_dimension_total_size(size_t ndim,size_t * dsize)43 gal_dimension_total_size(size_t ndim, size_t *dsize)
44 {
45 size_t i, num=1;
46 for(i=0;i<ndim;++i) num *= dsize[i];
47 return num;
48 }
49
50
51
52
53
54 int
gal_dimension_is_different(gal_data_t * first,gal_data_t * second)55 gal_dimension_is_different(gal_data_t *first, gal_data_t *second)
56 {
57 size_t i;
58
59 /* First make sure that the dimensionality is the same. */
60 if(first->ndim!=second->ndim)
61 return 1;
62
63 /* If the sizes are not zero, check if each dimension also has the same
64 length. */
65 if(first->size==0 && first->size==second->size)
66 return 0;
67 else
68 for(i=0;i<first->ndim;++i)
69 if( first->dsize[i] != second->dsize[i] )
70 return 1;
71
72 /* If it got to here, we know the dimensions have the same length. */
73 return 0;
74 }
75
76
77
78
79
80 /* Calculate the values necessary to increment/decrement along each
81 dimension of a dataset with size 'dsize'. */
82 size_t *
gal_dimension_increment(size_t ndim,size_t * dsize)83 gal_dimension_increment(size_t ndim, size_t *dsize)
84 {
85 int i;
86 size_t *out=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__, "out");
87
88 /* Along the fastest dimension, it is 1. */
89 out[ndim-1]=1;
90
91 /* For the rest of the dimensions, it is the multiple of the faster
92 dimension's length and the value for the previous dimension. */
93 if(ndim>1)
94 for(i=ndim-2;i>=0;--i)
95 out[i]=dsize[i+1]*out[i+1];
96
97 /* Return the allocated array. */
98 return out;
99 }
100
101
102
103
104
105 size_t
gal_dimension_num_neighbors(size_t ndim)106 gal_dimension_num_neighbors(size_t ndim)
107 {
108 if(ndim)
109 return pow(3, ndim)-1;
110 else
111 error(EXIT_FAILURE, 0, "%s: ndim cannot be zero", __func__);
112 return 0;
113 }
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132 /************************************************************************/
133 /******************** Coordinates **********************/
134 /************************************************************************/
135 void
gal_dimension_add_coords(size_t * c1,size_t * c2,size_t * out,size_t ndim)136 gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim)
137 {
138 size_t *end=c1+ndim;
139 do *out++ = *c1++ + *c2++; while(c1<end);
140 }
141
142
143
144
145
146 /* Return the index of an element from its coordinates. The index is the
147 position in the contiguous array (assuming it is a 1D arrray). */
148 size_t
gal_dimension_coord_to_index(size_t ndim,size_t * dsize,size_t * coord)149 gal_dimension_coord_to_index(size_t ndim, size_t *dsize, size_t *coord)
150 {
151 size_t i, d, ind=0, in_all_faster_dim;
152
153 switch(ndim)
154 {
155 case 0:
156 error(EXIT_FAILURE, 0, "%s: doesn't accept 0 dimensional arrays",
157 __func__);
158
159 case 1:
160 ind=coord[0];
161 break;
162
163 case 2:
164 ind=coord[0]*dsize[1]+coord[1];
165 break;
166
167 default:
168 for(d=0;d<ndim;++d)
169 {
170 /* First, find the number of elements in all dimensions faster
171 than this one. */
172 in_all_faster_dim=1;
173 for(i=d+1;i<ndim;++i)
174 in_all_faster_dim *= dsize[i];
175
176 /* Multiply it by the coordinate value of this dimension and add
177 to the index. */
178 ind += coord[d] * in_all_faster_dim;
179 }
180 }
181
182 /* Return the derived index. */
183 return ind;
184 }
185
186
187
188
189
190 /* You know the index ('ind') of a point/tile in an n-dimensional ('ndim')
191 array which has 'dsize[i]' elements along dimension 'i'. You want to
192 know the coordinates of that point along each dimension. The output is
193 not actually returned, it must be allocated ('ndim' elements) before
194 calling this function. This function will just fill it. The reason for
195 this is that this function will often be called with a loop and a single
196 allocated space would be enough for the whole loop. */
197 void
gal_dimension_index_to_coord(size_t index,size_t ndim,size_t * dsize,size_t * coord)198 gal_dimension_index_to_coord(size_t index, size_t ndim, size_t *dsize,
199 size_t *coord)
200 {
201 size_t d, *dinc;
202
203 switch(ndim)
204 {
205 case 0:
206 error(EXIT_FAILURE, 0, "%s: a 0-dimensional dataset is not defined",
207 __func__);
208
209 /* One dimensional dataset. */
210 case 1:
211 coord[0] = index;
212 break;
213
214 /* 2D dataset. */
215 case 2:
216 coord[0] = index / dsize[1];
217 coord[1] = index % dsize[1];
218 break;
219
220 /* Higher dimensional datasets. */
221 default:
222 /* Set the incrementation values for each dimension. */
223 dinc=gal_dimension_increment(ndim, dsize);
224
225 /* We start with the slowest dimension (first in the C standard) and
226 continue until (but not including) the fastest dimension. This is
227 because except for the fastest (coniguous) dimension, the other
228 coordinates can be found by division. */
229 for(d=0;d<ndim;++d)
230 {
231 /* Set the coordinate value for this dimension. */
232 coord[d] = index / dinc[d];
233
234 /* Replace the index with its remainder with the number of
235 elements in all faster dimensions. */
236 index %= dinc[d];
237 }
238
239 /* Clean up. */
240 free(dinc);
241 }
242 }
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263 /************************************************************************/
264 /******************** Distances **********************/
265 /************************************************************************/
266 float
gal_dimension_dist_manhattan(size_t * a,size_t * b,size_t ndim)267 gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t ndim)
268 {
269 size_t i, out=0;
270 for(i=0;i<ndim;++i) out += (a[i] > b[i]) ? (a[i]-b[i]) : (b[i]-a[i]);
271 return out;
272 }
273
274
275
276
277
278 float
gal_dimension_dist_radial(size_t * a,size_t * b,size_t ndim)279 gal_dimension_dist_radial(size_t *a, size_t *b, size_t ndim)
280 {
281 size_t i, out=0;
282 for(i=0;i<ndim;++i) out += (a[i]-b[i])*(a[i]-b[i]);
283 return sqrt(out);
284 }
285
286
287
288 /* Elliptical distance of a point from a given center. */
289 #define DEGREESTORADIANS M_PI/180.0
290 float
gal_dimension_dist_elliptical(double * center,double * pa_deg,double * q,size_t ndim,double * point)291 gal_dimension_dist_elliptical(double *center, double *pa_deg, double *q,
292 size_t ndim, double *point)
293 {
294 double Xr, Yr, Zr; /* Rotated x, y, z. */
295 double q1=q[0], q2;
296 double c1=cos( pa_deg[0] * DEGREESTORADIANS ), c2, c3;
297 double s1= sin( pa_deg[0] * DEGREESTORADIANS ), s2, s3;
298 double x=center[0]-point[0], y=center[1]-point[1], z;
299
300 /* Find the distance depending on the dimension. */
301 switch(ndim)
302 {
303 case 2:
304 /* The parenthesis aren't necessary, but help in readability and
305 avoiding human induced bugs. */
306 Xr = x * ( c1 ) + y * ( s1 );
307 Yr = x * ( -1.0f*s1 ) + y * ( c1 );
308 return sqrt( Xr*Xr + Yr*Yr/q1/q1 );
309 break;
310
311 case 3:
312 /* Define the necessary parameters. */
313 q2=q[1];
314 z=center[2]-point[2];
315 c2 = cos( pa_deg[1] * DEGREESTORADIANS );
316 s2 = sin( pa_deg[1] * DEGREESTORADIANS );
317 c3 = cos( pa_deg[2] * DEGREESTORADIANS );
318 s3 = sin( pa_deg[2] * DEGREESTORADIANS );
319
320 /* Calculate the distance. */
321 Xr = x*( c3*c1 - s3*c2*s1 ) + y*( c3*s1 + s3*c2*c1) + z*( s3*s2 );
322 Yr = x*( -1*s3*c1 - c3*c2*s1 ) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
323 Zr = x*( s1*s2 ) + y*(-1*s2*c1 ) + z*( c2 );
324 return sqrt( Xr*Xr + Yr*Yr/q1/q1 + Zr*Zr/q2/q2 );
325 break;
326
327 default:
328 error(EXIT_FAILURE, 0, "%s: currently only 2 and 3 dimensional "
329 "distances are supported, you have asked for an %zu-dimensional "
330 "dataset", __func__, ndim);
331
332 }
333
334 /* Control should never reach here. */
335 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address the "
336 "problem. Control should not reach the end of this function",
337 __func__, PACKAGE_BUGREPORT);
338 return NAN;
339 }
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359 /************************************************************************/
360 /******************** Collapsing a dimension **********************/
361 /************************************************************************/
362 enum dimension_collapse_operation
363 {
364 DIMENSION_COLLAPSE_INVALID, /* ==0 by C standard. */
365
366 DIMENSION_COLLAPSE_SUM,
367 DIMENSION_COLLAPSE_MAX,
368 DIMENSION_COLLAPSE_MIN,
369 DIMENSION_COLLAPSE_MEAN,
370 DIMENSION_COLLAPSE_NUMBER,
371 };
372
373
374
375
376
377 static gal_data_t *
dimension_collapse_sanity_check(gal_data_t * in,gal_data_t * weight,size_t c_dim,int hasblank,size_t * cnum,double ** warr)378 dimension_collapse_sanity_check(gal_data_t *in, gal_data_t *weight,
379 size_t c_dim, int hasblank, size_t *cnum,
380 double **warr)
381 {
382 gal_data_t *wht=NULL;
383
384 /* The requested dimension to collapse cannot be larger than the input's
385 number of dimensions. */
386 if( c_dim > (in->ndim-1) )
387 error(EXIT_FAILURE, 0, "%s: the input has %zu dimension(s), but you have "
388 "asked to collapse dimension %zu", __func__, in->ndim, c_dim);
389
390 /* If there is no blank value, there is no point in calculating the
391 number of points in each collapsed dataset (when necessary). In that
392 case, 'cnum!=0'. */
393 if(hasblank==0)
394 *cnum=in->dsize[c_dim];
395
396 /* Weight sanity checks */
397 if(weight)
398 {
399 if( weight->ndim!=1 )
400 error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu dimensions, "
401 "it must be one-dimensional", __func__, weight->ndim);
402 if( in->dsize[c_dim]!=weight->size )
403 error(EXIT_FAILURE, 0, "%s: the weight dataset has %zu elements, "
404 "but the input dataset has %zu elements in dimension %zu",
405 __func__, weight->size, in->dsize[c_dim], c_dim);
406 wht = ( weight->type == GAL_TYPE_FLOAT64
407 ? weight
408 : gal_data_copy_to_new_type(weight, GAL_TYPE_FLOAT64) );
409 *warr = wht->array;
410 }
411
412 /* Return the weight data structure. */
413 return wht;
414 }
415
416
417
418
419 /* Set the collapsed output sizes. */
420 static void
dimension_collapse_sizes(gal_data_t * in,size_t c_dim,size_t * outndim,size_t * outdsize)421 dimension_collapse_sizes(gal_data_t *in, size_t c_dim, size_t *outndim,
422 size_t *outdsize)
423 {
424 size_t i, a=0;
425
426 if(in->ndim==1)
427 *outndim=outdsize[0]=1;
428 else
429 {
430 *outndim=in->ndim-1;
431 for(i=0;i<in->ndim;++i)
432 if(i!=c_dim) outdsize[a++]=in->dsize[i];
433 }
434 }
435
436
437
438
439
440 /* Depending on the operator, write the result into the output. */
441 #define COLLAPSE_WRITE(OIND,IIND) { \
442 /* Sum */ \
443 if(farr) \
444 farr[ OIND ] += (warr ? warr[w] : 1) * inarr[ IIND ]; \
445 \
446 /* Number */ \
447 if(iarr) \
448 { \
449 if(num->type==GAL_TYPE_UINT8) iarr[ OIND ] = 1; \
450 else ++iarr[ OIND ]; \
451 } \
452 \
453 /* Sum of weights. */ \
454 if(wsumarr) wsumarr[ OIND ] += warr[w]; \
455 \
456 /* Minimum or maximum. */ \
457 if(mmarr) \
458 mmarr[OIND] = ( max1_min0 \
459 ? (mmarr[OIND]>inarr[IIND]?mmarr[OIND]:inarr[IIND]) \
460 : (mmarr[OIND]<inarr[IIND]?mmarr[OIND]:inarr[IIND]) ); \
461 }
462
463
464 /* Deal properly with blanks. */
465 #define COLLAPSE_CHECKBLANK(OIND,IIND) { \
466 if(hasblank) \
467 { \
468 if(B==B) /* An integer type: blank can be checked with '=='. */ \
469 { \
470 if( inarr[IIND] != B ) COLLAPSE_WRITE(OIND,IIND); \
471 } \
472 else /* A floating point type where NAN != NAN. */ \
473 { \
474 if( inarr[IIND] == inarr[IIND] ) COLLAPSE_WRITE(OIND,IIND); \
475 } \
476 } \
477 else COLLAPSE_WRITE(OIND,IIND); \
478 }
479
480
481
482 #define COLLAPSE_DIM(IT) { \
483 IT m, B, *inarr=in->array, *mmarr=minmax?minmax->array:NULL; \
484 if(hasblank) gal_blank_write(&B, in->type); \
485 \
486 /* Initialize the array for minimum or maximum. */ \
487 if(mmarr) \
488 { \
489 if(max1_min0) gal_type_min(in->type, &m); \
490 else gal_type_max(in->type, &m); \
491 for(i=0;i<minmax->size;++i) mmarr[i]=m; \
492 } \
493 \
494 /* Collapse the dataset. */ \
495 switch(in->ndim) \
496 { \
497 /* 1D input dataset. */ \
498 case 1: \
499 for(i=0;i<in->dsize[0];++i) \
500 { \
501 if(weight) w=i; \
502 COLLAPSE_CHECKBLANK(0,i); \
503 } \
504 break; \
505 \
506 /* 2D input dataset. */ \
507 case 2: \
508 for(i=0;i<in->dsize[0];++i) \
509 for(j=0;j<in->dsize[1];++j) \
510 { \
511 /* In a more easy to understand format: \
512 dim==0 --> a=j; \
513 dim==1 --> a=i; */ \
514 a = c_dim==0 ? j : i; \
515 if(weight) w = c_dim == 0 ? i : j; \
516 COLLAPSE_CHECKBLANK(a, i*in->dsize[1] + j); \
517 } \
518 break; \
519 \
520 /* 3D input dataset. */ \
521 case 3: \
522 slice=in->dsize[1]*in->dsize[2]; \
523 for(i=0;i<in->dsize[0];++i) \
524 for(j=0;j<in->dsize[1];++j) \
525 for(k=0;k<in->dsize[2];++k) \
526 { \
527 /* In a more easy to understand format: \
528 dim==0 --> a=j; b=k; \
529 dim==1 --> a=i; b=k; \
530 dim==2 --> a=i; b=j; */ \
531 a = c_dim==0 ? j : i; \
532 b = c_dim==2 ? j : k; \
533 if(weight) w = c_dim==0 ? i : (c_dim==1 ? j : k); \
534 COLLAPSE_CHECKBLANK(a*outdsize[1]+b, \
535 i*slice + j*in->dsize[2] + k); \
536 } \
537 break; \
538 \
539 /* Input dataset's dimensionality not yet supported. */ \
540 default: \
541 error(EXIT_FAILURE, 0, "%s: %zu-dimensional datasets not yet " \
542 "supported, please contact us at %s to add this feature", \
543 __func__, in->ndim, PACKAGE_BUGREPORT); \
544 } \
545 \
546 /* For minimum or maximum, elements with no input must be blank. */ \
547 if(mmarr && iarr) \
548 for(i=0;i<minmax->size;++i) if(iarr[i]==0) mmarr[i]=B; \
549 }
550
551
552
553
554
555 gal_data_t *
gal_dimension_collapse_sum(gal_data_t * in,size_t c_dim,gal_data_t * weight)556 gal_dimension_collapse_sum(gal_data_t *in, size_t c_dim, gal_data_t *weight)
557 {
558 int max1_min0=0;
559 double *wsumarr=NULL;
560 uint8_t *ii, *iarr=NULL;
561 size_t a, b, i, j, k, w=-1, cnum=0;
562 size_t outdsize[10], slice, outndim;
563 int hasblank=gal_blank_present(in, 0);
564 double *dd, *df, *warr=NULL, *farr=NULL;
565 gal_data_t *sum=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
566
567 /* Basic sanity checks. */
568 wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
569 &cnum, &warr);
570
571 /* Set the size of the collapsed output. */
572 dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
573
574 /* Allocate the sum (output) dataset. */
575 sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, outndim, outdsize, in->wcs,
576 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
577
578 /* The number dataset (when there are blank values).*/
579 if(hasblank)
580 num=gal_data_alloc(NULL, GAL_TYPE_INT8, outndim, outdsize, NULL,
581 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
582
583 /* Set the array pointers. */
584 if(sum) farr=sum->array;
585 if(num) iarr=num->array;
586
587 /* Parse the dataset. */
588 switch(in->type)
589 {
590 case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
591 case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
592 case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
593 case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
594 case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
595 case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
596 case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
597 case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
598 case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
599 case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
600 default:
601 error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
602 __func__, in->type);
603 }
604
605 /* If 'num' is zero on any element, set its sum to NaN. */
606 if(num)
607 {
608 ii = num->array;
609 df = (dd=sum->array) + sum->size;
610 do if(*ii++==0) *dd=NAN; while(++dd<df);
611 }
612
613 /* Remove the respective dimension in the WCS structure also (if any
614 exists). Note that 'sum->ndim' has already been changed. So we'll use
615 'in->wcs'. */
616 gal_wcs_remove_dimension(sum->wcs, in->ndim-c_dim);
617
618 /* Clean up and return. */
619 if(wht!=weight) gal_data_free(wht);
620 if(num) gal_data_free(num);
621 return sum;
622 }
623
624
625
626
627
628 gal_data_t *
gal_dimension_collapse_mean(gal_data_t * in,size_t c_dim,gal_data_t * weight)629 gal_dimension_collapse_mean(gal_data_t *in, size_t c_dim,
630 gal_data_t *weight)
631 {
632 int max1_min0=0;
633 double wsum=NAN;
634 double *wsumarr=NULL;
635 int32_t *ii, *iarr=NULL;
636 size_t a, b, i, j, k, w=-1, cnum=0;
637 size_t outdsize[10], slice, outndim;
638 int hasblank=gal_blank_present(in, 0);
639 double *dd, *dw, *df, *warr=NULL, *farr=NULL;
640 gal_data_t *sum=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
641
642
643 /* Basic sanity checks. */
644 wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
645 &cnum, &warr);
646
647 /* Set the size of the collapsed output. */
648 dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
649
650 /* The sum array. */
651 sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, outndim, outdsize, in->wcs,
652 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
653
654 /* If a weighted mean is requested. */
655 if( weight )
656 {
657 /* There are blank values, so we'll need to keep the sums of the
658 weights for each collapsed dimension */
659 if( hasblank )
660 wsumarr=gal_pointer_allocate(GAL_TYPE_FLOAT64, sum->size, 1,
661 __func__, "wsumarr");
662
663 /* There aren't any blank values, so one summation over the
664 weights is enough to calculate the weighted mean. */
665 else
666 {
667 wsum=0.0f;
668 df=(dd=weight->array)+weight->size;
669 do wsum += *dd++; while(dd<df);
670 }
671 }
672 /* No weight is given, so we'll need the number of elements. */
673 else if( hasblank )
674 num=gal_data_alloc(NULL, GAL_TYPE_INT32, outndim, outdsize, NULL,
675 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
676
677 /* Set the array pointers. */
678 if(sum) farr=sum->array;
679 if(num) iarr=num->array;
680
681 /* Parse the dataset. */
682 switch(in->type)
683 {
684 case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
685 case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
686 case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
687 case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
688 case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
689 case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
690 case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
691 case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
692 case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
693 case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
694 default:
695 error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
696 __func__, in->type);
697 }
698
699 /* If 'num' is zero on any element, set its sum to NaN. */
700 if(num)
701 {
702 ii = num->array;
703 df = (dd=sum->array) + sum->size;
704 do if(*ii++==0) *dd=NAN; while(++dd<df);
705 }
706
707 /* Divide the sum by the number. */
708 df = (dd=sum->array) + sum->size;
709 if(weight)
710 {
711 if(hasblank) { dw=wsumarr; do *dd /= *dw++; while(++dd<df); }
712 else do *dd /= wsum; while(++dd<df);
713 }
714 else
715 if(num) { ii = num->array; do *dd /= *ii++; while(++dd<df); }
716 else do *dd /= cnum; while(++dd<df);
717
718 /* Correct the WCS, clean up and return. */
719 gal_wcs_remove_dimension(sum->wcs, in->ndim-c_dim);
720 if(wht!=weight) gal_data_free(wht);
721 if(wsumarr) free(wsumarr);
722 gal_data_free(num);
723 return sum;
724 }
725
726
727
728
729
730 gal_data_t *
gal_dimension_collapse_number(gal_data_t * in,size_t c_dim)731 gal_dimension_collapse_number(gal_data_t *in, size_t c_dim)
732 {
733 int max1_min0=0;
734 double *wsumarr=NULL;
735 double *warr=NULL, *farr=NULL;
736 int32_t *ii, *iif, *iarr=NULL;
737 size_t a, b, i, j, k, w, cnum=0;
738 size_t outdsize[10], slice, outndim;
739 int hasblank=gal_blank_present(in, 0);
740 gal_data_t *weight=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
741
742 /* Basic sanity checks. */
743 wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
744 &cnum, &warr);
745
746 /* Set the size of the collapsed output. */
747 dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
748
749 /* The number dataset (when there are blank values).*/
750 num=gal_data_alloc(NULL, GAL_TYPE_INT32, outndim, outdsize, in->wcs,
751 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
752
753 /* Set the array pointers. */
754 iarr=num->array;
755
756 /* Parse the input dataset (if necessary). */
757 if(hasblank)
758 switch(in->type)
759 {
760 case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
761 case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
762 case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
763 case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
764 case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
765 case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
766 case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
767 case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
768 case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
769 case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
770 default:
771 error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
772 __func__, in->type);
773 }
774 else
775 {
776 iif=(ii=num->array)+num->size;
777 do *ii++ = cnum; while(ii<iif);
778 }
779
780 /* Remove the respective dimension in the WCS structure also (if any
781 exists). Note that 'sum->ndim' has already been changed. So we'll use
782 'in->wcs'. */
783 gal_wcs_remove_dimension(num->wcs, in->ndim-c_dim);
784
785 /* Return. */
786 if(wht!=weight) gal_data_free(wht);
787 return num;
788 }
789
790
791
792
793
794 gal_data_t *
gal_dimension_collapse_minmax(gal_data_t * in,size_t c_dim,int max1_min0)795 gal_dimension_collapse_minmax(gal_data_t *in, size_t c_dim, int max1_min0)
796 {
797 uint8_t *iarr=NULL;
798 double *wsumarr=NULL;
799 double *warr=NULL, *farr=NULL;
800 size_t a, b, i, j, k, w, cnum=0;
801 size_t outdsize[10], slice, outndim;
802 int hasblank=gal_blank_present(in, 0);
803 gal_data_t *weight=NULL, *wht=NULL, *num=NULL, *minmax=NULL;
804
805 /* Basic sanity checks. */
806 wht=dimension_collapse_sanity_check(in, weight, c_dim, hasblank,
807 &cnum, &warr);
808
809 /* Set the size of the collapsed output. */
810 dimension_collapse_sizes(in, c_dim, &outndim, outdsize);
811
812 /* Allocate the necessary datasets. If there are blank pixels, we'll need
813 to count how many elements whent into the calculation so we can set
814 them to blank. */
815 minmax=gal_data_alloc(NULL, in->type, outndim, outdsize, in->wcs,
816 0, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
817 if(hasblank)
818 {
819 num=gal_data_alloc(NULL, GAL_TYPE_UINT8, outndim, outdsize, in->wcs,
820 1, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
821 iarr=num->array;
822 }
823
824 /* Parse the input dataset (if necessary). */
825 switch(in->type)
826 {
827 case GAL_TYPE_UINT8: COLLAPSE_DIM( uint8_t ); break;
828 case GAL_TYPE_INT8: COLLAPSE_DIM( int8_t ); break;
829 case GAL_TYPE_UINT16: COLLAPSE_DIM( uint16_t ); break;
830 case GAL_TYPE_INT16: COLLAPSE_DIM( int16_t ); break;
831 case GAL_TYPE_UINT32: COLLAPSE_DIM( uint32_t ); break;
832 case GAL_TYPE_INT32: COLLAPSE_DIM( int32_t ); break;
833 case GAL_TYPE_UINT64: COLLAPSE_DIM( uint64_t ); break;
834 case GAL_TYPE_INT64: COLLAPSE_DIM( int64_t ); break;
835 case GAL_TYPE_FLOAT32: COLLAPSE_DIM( float ); break;
836 case GAL_TYPE_FLOAT64: COLLAPSE_DIM( double ); break;
837 default:
838 error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
839 __func__, in->type);
840 }
841
842 /* Remove the respective dimension in the WCS structure also (if any
843 exists). Note that 'sum->ndim' has already been changed. So we'll use
844 'in->wcs'. */
845 gal_wcs_remove_dimension(minmax->wcs, in->ndim-c_dim);
846
847 /* Clean up and return. */
848 if(wht!=weight) gal_data_free(wht);
849 if(num) gal_data_free(num);
850 return minmax;
851 }
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872 /************************************************************************/
873 /******************** Other **********************/
874 /************************************************************************/
875 size_t
gal_dimension_remove_extra(size_t ndim,size_t * dsize,struct wcsprm * wcs)876 gal_dimension_remove_extra(size_t ndim, size_t *dsize, struct wcsprm *wcs)
877 {
878 size_t i, j;
879
880 for(i=0;i<ndim;++i)
881 if(dsize[i]==1)
882 {
883 /* Correct the WCS. */
884 if(wcs) gal_wcs_remove_dimension(wcs, ndim-i);
885
886 /* Shift all subsequent dimensions to replace this one. */
887 for(j=i;j<ndim-1;++j) dsize[j]=dsize[j+1];
888
889 /* Decrement the 'i' and the total number of dimension. */
890 --i;
891 --ndim;
892 }
893
894 /* Return the number of dimensions. */
895 return ndim;
896 }
897