1 /*********************************************************************
2 Arithmetic operations on data structures
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) 2016-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 <errno.h>
26 #include <error.h>
27 #include <stdio.h>
28 #include <string.h>
29 #include <stdlib.h>
30 #include <stdarg.h>
31
32 #include <gsl/gsl_rng.h>
33 #include <gsl/gsl_randist.h>
34
35 #include <gnuastro/box.h>
36 #include <gnuastro/list.h>
37 #include <gnuastro/blank.h>
38 #include <gnuastro/units.h>
39 #include <gnuastro/qsort.h>
40 #include <gnuastro/pointer.h>
41 #include <gnuastro/threads.h>
42 #include <gnuastro/dimension.h>
43 #include <gnuastro/statistics.h>
44 #include <gnuastro/arithmetic.h>
45
46 #include <gnuastro-internal/checkset.h>
47 #include <gnuastro-internal/arithmetic-internal.h>
48
49 /* Headers for each binary operator. Since they heavily involve macros,
50 their compilation can be very large if they are in a single function and
51 file. So there is a separate C source and header file for each of these
52 functions.*/
53 #include <gnuastro-internal/arithmetic-lt.h>
54 #include <gnuastro-internal/arithmetic-le.h>
55 #include <gnuastro-internal/arithmetic-gt.h>
56 #include <gnuastro-internal/arithmetic-ge.h>
57 #include <gnuastro-internal/arithmetic-eq.h>
58 #include <gnuastro-internal/arithmetic-ne.h>
59 #include <gnuastro-internal/arithmetic-or.h>
60 #include <gnuastro-internal/arithmetic-and.h>
61 #include <gnuastro-internal/arithmetic-plus.h>
62 #include <gnuastro-internal/arithmetic-minus.h>
63 #include <gnuastro-internal/arithmetic-bitor.h>
64 #include <gnuastro-internal/arithmetic-bitand.h>
65 #include <gnuastro-internal/arithmetic-bitxor.h>
66 #include <gnuastro-internal/arithmetic-bitlsh.h>
67 #include <gnuastro-internal/arithmetic-bitrsh.h>
68 #include <gnuastro-internal/arithmetic-modulo.h>
69 #include <gnuastro-internal/arithmetic-divide.h>
70 #include <gnuastro-internal/arithmetic-multiply.h>
71
72
73
74
75
76
77
78
79
80
81 /***********************************************************************/
82 /*************** Internal checks **************/
83 /***********************************************************************/
84 /* Some functions are only for a floating point operand, so if the input
85 isn't floating point, inform the user to change the type explicitly,
86 doing it implicitly/internally puts too much responsability on the
87 program.
88 static void
89 arithmetic_check_float_input(gal_data_t *in, int operator, char *numstr)
90 {
91 switch(in->type)
92 {
93 case GAL_TYPE_FLOAT32:
94 case GAL_TYPE_FLOAT64:
95 break;
96 default:
97 error(EXIT_FAILURE, 0, "the %s operator can only accept single or "
98 "double precision floating point numbers as its operand. The "
99 "%s operand has type %s. You can use the 'float' or 'double' "
100 "operators before this operator to explicitly convert to the "
101 "desired precision floating point type. If the operand was "
102 "originally a typed number (string of characters), add an 'f' "
103 "after it so it is directly read into the proper precision "
104 "floating point number (based on the number of non-zero "
105 "decimals it has)", gal_arithmetic_operator_string(operator),
106 numstr, gal_type_name(in->type, 1));
107 }
108 }
109 */
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129 /***********************************************************************/
130 /*************** Unary functions/operators **************/
131 /***********************************************************************/
132 /* Change input data structure type. */
133 static gal_data_t *
arithmetic_change_type(gal_data_t * data,int operator,int flags)134 arithmetic_change_type(gal_data_t *data, int operator, int flags)
135 {
136 int type=-1;
137 gal_data_t *out;
138
139 /* Set the output type. */
140 switch(operator)
141 {
142 case GAL_ARITHMETIC_OP_TO_UINT8: type=GAL_TYPE_UINT8; break;
143 case GAL_ARITHMETIC_OP_TO_INT8: type=GAL_TYPE_INT8; break;
144 case GAL_ARITHMETIC_OP_TO_UINT16: type=GAL_TYPE_UINT16; break;
145 case GAL_ARITHMETIC_OP_TO_INT16: type=GAL_TYPE_INT16; break;
146 case GAL_ARITHMETIC_OP_TO_UINT32: type=GAL_TYPE_UINT32; break;
147 case GAL_ARITHMETIC_OP_TO_INT32: type=GAL_TYPE_INT32; break;
148 case GAL_ARITHMETIC_OP_TO_UINT64: type=GAL_TYPE_UINT64; break;
149 case GAL_ARITHMETIC_OP_TO_INT64: type=GAL_TYPE_INT64; break;
150 case GAL_ARITHMETIC_OP_TO_FLOAT32: type=GAL_TYPE_FLOAT32; break;
151 case GAL_ARITHMETIC_OP_TO_FLOAT64: type=GAL_TYPE_FLOAT64; break;
152 default:
153 error(EXIT_FAILURE, 0, "%s: operator value of %d not recognized",
154 __func__, operator);
155 }
156
157 /* Copy to the new type. */
158 out=gal_data_copy_to_new_type(data, type);
159
160 /* Delete the input structure if the user asked for it. */
161 if(flags & GAL_ARITHMETIC_FLAG_FREE)
162 gal_data_free(data);
163
164 /* Return */
165 return out;
166 }
167
168
169
170
171
172 /* Return an array of value 1 for any zero valued element and zero for any
173 non-zero valued element. */
174 #define TYPE_CASE_FOR_NOT(CTYPE) { \
175 CTYPE *a=data->array, *af=a+data->size; \
176 do *o++ = !(*a); while(++a<af); \
177 }
178
179 static gal_data_t *
arithmetic_not(gal_data_t * data,int flags)180 arithmetic_not(gal_data_t *data, int flags)
181 {
182 uint8_t *o;
183 gal_data_t *out;
184
185 /* Allocate the output array. */
186 out=gal_data_alloc(NULL, GAL_TYPE_UINT8, data->ndim, data->dsize,
187 data->wcs, 0, data->minmapsize, data->quietmmap,
188 data->name, data->unit, data->comment);
189 o=out->array;
190
191
192 /* Go over the pixels and set the output values. */
193 switch(data->type)
194 {
195 case GAL_TYPE_UINT8: TYPE_CASE_FOR_NOT( uint8_t ); break;
196 case GAL_TYPE_INT8: TYPE_CASE_FOR_NOT( int8_t ); break;
197 case GAL_TYPE_UINT16: TYPE_CASE_FOR_NOT( uint16_t ); break;
198 case GAL_TYPE_INT16: TYPE_CASE_FOR_NOT( int16_t ); break;
199 case GAL_TYPE_UINT32: TYPE_CASE_FOR_NOT( uint32_t ); break;
200 case GAL_TYPE_INT32: TYPE_CASE_FOR_NOT( int32_t ); break;
201 case GAL_TYPE_UINT64: TYPE_CASE_FOR_NOT( uint64_t ); break;
202 case GAL_TYPE_INT64: TYPE_CASE_FOR_NOT( int64_t ); break;
203 case GAL_TYPE_FLOAT32: TYPE_CASE_FOR_NOT( float ); break;
204 case GAL_TYPE_FLOAT64: TYPE_CASE_FOR_NOT( double ); break;
205
206 case GAL_TYPE_BIT:
207 error(EXIT_FAILURE, 0, "%s: bit datatypes are not yet supported, "
208 "please get in touch with us to implement it.", __func__);
209
210 default:
211 error(EXIT_FAILURE, 0, "%s: type value (%d) not recognized",
212 __func__, data->type);
213 }
214
215 /* Delete the input structure if the user asked for it. */
216 if(flags & GAL_ARITHMETIC_FLAG_FREE)
217 gal_data_free(data);
218
219 /* Return */
220 return out;
221 }
222
223
224
225
226
227 /* Bitwise not operator. */
228 static gal_data_t *
arithmetic_bitwise_not(int flags,gal_data_t * in)229 arithmetic_bitwise_not(int flags, gal_data_t *in)
230 {
231 gal_data_t *o;
232 uint8_t *iu8 = in->array, *iu8f = iu8 + in->size, *ou8;
233 int8_t *ii8 = in->array, *ii8f = ii8 + in->size, *oi8;
234 uint16_t *iu16 = in->array, *iu16f = iu16 + in->size, *ou16;
235 int16_t *ii16 = in->array, *ii16f = ii16 + in->size, *oi16;
236 uint32_t *iu32 = in->array, *iu32f = iu32 + in->size, *ou32;
237 int32_t *ii32 = in->array, *ii32f = ii32 + in->size, *oi32;
238 uint64_t *iu64 = in->array, *iu64f = iu64 + in->size, *ou64;
239 int64_t *ii64 = in->array, *ii64f = ii64 + in->size, *oi64;
240
241 /* Check the type */
242 switch(in->type)
243 {
244 case GAL_TYPE_FLOAT32:
245 case GAL_TYPE_FLOAT64:
246 error(EXIT_FAILURE, 0, "%s: bitwise not (one's complement) "
247 "operator can only work on integer types", __func__);
248 }
249
250 /* If we want inplace output, set the output pointer to the input
251 pointer, for every pixel, the operation will be independent. */
252 if(flags & GAL_ARITHMETIC_FLAG_INPLACE)
253 o = in;
254 else
255 o = gal_data_alloc(NULL, in->type, in->ndim, in->dsize, in->wcs,
256 0, in->minmapsize, in->quietmmap, NULL, NULL, NULL);
257
258 /* Start setting the types. */
259 switch(in->type)
260 {
261 case GAL_TYPE_UINT8:
262 ou8=o->array; do *ou8++ = ~(*iu8++); while(iu8<iu8f); break;
263
264 case GAL_TYPE_INT8:
265 oi8=o->array; do *oi8++ = ~(*ii8++); while(ii8<ii8f); break;
266
267 case GAL_TYPE_UINT16:
268 ou16=o->array; do *ou16++ = ~(*iu16++); while(iu16<iu16f); break;
269
270 case GAL_TYPE_INT16:
271 oi16=o->array; do *oi16++ = ~(*ii16++); while(ii16<ii16f); break;
272
273 case GAL_TYPE_UINT32:
274 ou32=o->array; do *ou32++ = ~(*iu32++); while(iu32<iu32f); break;
275
276 case GAL_TYPE_INT32:
277 oi32=o->array; do *oi32++ = ~(*ii32++); while(ii32<ii32f); break;
278
279 case GAL_TYPE_UINT64:
280 ou64=o->array; do *ou64++ = ~(*iu64++); while(iu64<iu64f); break;
281
282 case GAL_TYPE_INT64:
283 oi64=o->array; do *oi64++ = ~(*ii64++); while(ii64<ii64f); break;
284
285 default:
286 error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
287 __func__, in->type);
288 }
289
290
291 /* Clean up (if necessary). */
292 if( (flags & GAL_ARITHMETIC_FLAG_FREE) && o!=in)
293 gal_data_free(in);
294
295 /* Return */
296 return o;
297 }
298
299
300
301
302
303 /* We don't want to use the standard function for unary functions in the
304 case of the absolute operator. This is because there are multiple
305 versions of this function in the C library for different types, which
306 can greatly improve speed. */
307 #define ARITHMETIC_ABS_SGN(CTYPE, FUNC) { \
308 CTYPE *o=out->array, *a=in->array, *af=a+in->size; \
309 do *o++ = FUNC(*a); while(++a<af); \
310 }
311 static gal_data_t *
arithmetic_abs(int flags,gal_data_t * in)312 arithmetic_abs(int flags, gal_data_t *in)
313 {
314 gal_data_t *out;
315
316 /* Set the output array. */
317 if(flags & GAL_ARITHMETIC_FLAG_INPLACE)
318 out=in;
319 else
320 out = gal_data_alloc(NULL, in->type, in->ndim, in->dsize,
321 in->wcs, 0, in->minmapsize, in->quietmmap,
322 in->name, in->unit, in->comment);
323
324 /* Put the absolute value depending on the type. */
325 switch(in->type)
326 {
327 /* Unsigned types are already positive, so if the input is not to be
328 freed (the output must be a separate array), just copy the whole
329 array. */
330 case GAL_TYPE_UINT8:
331 case GAL_TYPE_UINT16:
332 case GAL_TYPE_UINT32:
333 case GAL_TYPE_UINT64:
334 if(out!=in) gal_data_copy_to_allocated(in, out);
335 break;
336
337 /* For the signed types, we actually have to go over the data and
338 calculate the absolute value. There are unique functions for
339 different types, so we will be using them.*/
340 case GAL_TYPE_INT8: ARITHMETIC_ABS_SGN( int8_t, abs ); break;
341 case GAL_TYPE_INT16: ARITHMETIC_ABS_SGN( int16_t, abs ); break;
342 case GAL_TYPE_INT32: ARITHMETIC_ABS_SGN( int32_t, labs ); break;
343 case GAL_TYPE_INT64: ARITHMETIC_ABS_SGN( int64_t, llabs ); break;
344 case GAL_TYPE_FLOAT32: ARITHMETIC_ABS_SGN( float, fabsf ); break;
345 case GAL_TYPE_FLOAT64: ARITHMETIC_ABS_SGN( double, fabs ); break;
346 default:
347 error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
348 __func__, in->type);
349 }
350
351 /* Clean up and return */
352 if( (flags & GAL_ARITHMETIC_FLAG_FREE) && out!=in)
353 gal_data_free(in);
354 return out;
355 }
356
357
358
359
360
361 /* Wrapper functions for RA/Dec strings. */
362 static char *
arithmetic_units_degree_to_ra(double decimal)363 arithmetic_units_degree_to_ra(double decimal)
364 { return gal_units_degree_to_ra(decimal, 0); }
365
366 static char *
arithmetic_units_degree_to_dec(double decimal)367 arithmetic_units_degree_to_dec(double decimal)
368 { return gal_units_degree_to_dec(decimal, 0); }
369
370 #define UNIFUNC_RUN_FUNCTION_ON_ELEMENT(OT, IT, OP, BEFORE, AFTER){ \
371 OT *oa=o->array; \
372 IT *ia=in->array, *iaf=ia + in->size; \
373 do *oa++ = OP( *ia++ BEFORE ) AFTER; while(ia<iaf); \
374 }
375
376 #define UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(IT, OP, BEFORE, AFTER) \
377 switch(o->type) \
378 { \
379 case GAL_TYPE_UINT8: \
380 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint8_t, IT, OP, BEFORE, AFTER) \
381 break; \
382 case GAL_TYPE_INT8: \
383 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int8_t, IT, OP, BEFORE, AFTER) \
384 break; \
385 case GAL_TYPE_UINT16: \
386 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint16_t, IT, OP, BEFORE, AFTER) \
387 break; \
388 case GAL_TYPE_INT16: \
389 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int16_t, IT, OP, BEFORE, AFTER) \
390 break; \
391 case GAL_TYPE_UINT32: \
392 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint32_t, IT, OP, BEFORE, AFTER) \
393 break; \
394 case GAL_TYPE_INT32: \
395 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int32_t, IT, OP, BEFORE, AFTER) \
396 break; \
397 case GAL_TYPE_UINT64: \
398 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(uint64_t, IT, OP, BEFORE, AFTER) \
399 break; \
400 case GAL_TYPE_INT64: \
401 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(int64_t, IT, OP, BEFORE, AFTER) \
402 break; \
403 case GAL_TYPE_FLOAT32: \
404 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(float, IT, OP, BEFORE, AFTER) \
405 break; \
406 case GAL_TYPE_FLOAT64: \
407 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(double, IT, OP, BEFORE, AFTER) \
408 break; \
409 default: \
410 error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
411 "UNIARY_FUNCTION_ON_ELEMENT", in->type); \
412 }
413
414 #define UNIARY_FUNCTION_ON_ELEMENT_OUTPUT_STRING(OP) \
415 switch(in->type) \
416 { \
417 case GAL_TYPE_UINT8: \
418 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, uint8_t, OP, +0, +0) \
419 break; \
420 case GAL_TYPE_INT8: \
421 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, int8_t, OP, +0, +0) \
422 break; \
423 case GAL_TYPE_UINT16: \
424 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, uint16_t, OP, +0, +0) \
425 break; \
426 case GAL_TYPE_INT16: \
427 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, int16_t, OP, +0, +0) \
428 break; \
429 case GAL_TYPE_UINT32: \
430 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, uint32_t, OP, +0, +0) \
431 break; \
432 case GAL_TYPE_INT32: \
433 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, int32_t, OP, +0, +0) \
434 break; \
435 case GAL_TYPE_UINT64: \
436 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, uint64_t, OP, +0, +0) \
437 break; \
438 case GAL_TYPE_INT64: \
439 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, int64_t, OP, +0, +0) \
440 break; \
441 case GAL_TYPE_FLOAT32: \
442 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, float, OP, +0, +0) \
443 break; \
444 case GAL_TYPE_FLOAT64: \
445 UNIFUNC_RUN_FUNCTION_ON_ELEMENT(char *, double, OP, +0, +0) \
446 break; \
447 default: \
448 error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
449 "UNIARY_FUNCTION_ON_ELEMENT_OUTPUT_STRING", in->type); \
450 }
451
452 #define UNIARY_FUNCTION_ON_ELEMENT(OP, BEFORE, AFTER) \
453 switch(in->type) \
454 { \
455 case GAL_TYPE_UINT8: \
456 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint8_t, OP, BEFORE, AFTER) \
457 break; \
458 case GAL_TYPE_INT8: \
459 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int8_t, OP, BEFORE, AFTER) \
460 break; \
461 case GAL_TYPE_UINT16: \
462 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint16_t, OP, BEFORE, AFTER) \
463 break; \
464 case GAL_TYPE_INT16: \
465 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int16_t, OP, BEFORE, AFTER) \
466 break; \
467 case GAL_TYPE_UINT32: \
468 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint32_t, OP, BEFORE, AFTER) \
469 break; \
470 case GAL_TYPE_INT32: \
471 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int32_t, OP, BEFORE, AFTER) \
472 break; \
473 case GAL_TYPE_UINT64: \
474 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(uint64_t, OP, BEFORE, AFTER) \
475 break; \
476 case GAL_TYPE_INT64: \
477 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(int64_t, OP, BEFORE, AFTER) \
478 break; \
479 case GAL_TYPE_FLOAT32: \
480 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(float, OP, BEFORE, AFTER) \
481 break; \
482 case GAL_TYPE_FLOAT64: \
483 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_INSET(double, OP, BEFORE, AFTER) \
484 break; \
485 default: \
486 error(EXIT_FAILURE, 0, "%s: type code %d not recognized", \
487 "UNIARY_FUNCTION_ON_ELEMENT", in->type); \
488 }
489
490
491 #define UNIFUNC_RUN_FUNCTION_ON_ELEMENT_STRING(OT, OP){ \
492 OT *oa=o->array; \
493 char **ia=in->array, **iaf=ia + in->size; \
494 do *oa++ = OP(*ia++); while(ia<iaf); \
495 }
496
497 static gal_data_t *
arithmetic_function_unary(int operator,int flags,gal_data_t * in)498 arithmetic_function_unary(int operator, int flags, gal_data_t *in)
499 {
500 uint8_t otype;
501 int inplace=0;
502 gal_data_t *o;
503 double pi=3.14159265358979323846264338327;
504
505 /* See if the operation should be done in place. The output of these
506 operators is defined in the floating point space. So even if the input
507 is integer type and user requested inplace opereation, if its not a
508 floating point type, it will not be in-place. */
509 if( (flags & GAL_ARITHMETIC_FLAG_INPLACE)
510 && ( in->type==GAL_TYPE_FLOAT32 || in->type==GAL_TYPE_FLOAT64 )
511 && ( operator != GAL_ARITHMETIC_OP_RA_TO_DEGREE
512 && operator != GAL_ARITHMETIC_OP_DEC_TO_DEGREE
513 && operator != GAL_ARITHMETIC_OP_DEGREE_TO_RA
514 && operator != GAL_ARITHMETIC_OP_DEGREE_TO_DEC ) )
515 inplace=1;
516
517 /* Set the output pointer. */
518 if(inplace)
519 {
520 o = in;
521 otype=in->type;
522 }
523 else
524 {
525 /* Check for operators which have fixed output types */
526 if( operator == GAL_ARITHMETIC_OP_RA_TO_DEGREE
527 || operator == GAL_ARITHMETIC_OP_DEC_TO_DEGREE )
528 otype = GAL_TYPE_FLOAT64;
529 else if( operator == GAL_ARITHMETIC_OP_DEGREE_TO_RA
530 || operator == GAL_ARITHMETIC_OP_DEGREE_TO_DEC )
531 otype = GAL_TYPE_STRING;
532 else
533 otype = ( in->type==GAL_TYPE_FLOAT64
534 ? GAL_TYPE_FLOAT64
535 : GAL_TYPE_FLOAT32 );
536
537 /* Set the final output type. */
538 o = gal_data_alloc(NULL, otype, in->ndim, in->dsize, in->wcs,
539 0, in->minmapsize, in->quietmmap,
540 NULL, NULL, NULL);
541 }
542
543 /* Start setting the operator and operands. */
544 switch(operator)
545 {
546 case GAL_ARITHMETIC_OP_SQRT:
547 UNIARY_FUNCTION_ON_ELEMENT( sqrt, +0, +0); break;
548 case GAL_ARITHMETIC_OP_LOG:
549 UNIARY_FUNCTION_ON_ELEMENT( log, +0, +0); break;
550 case GAL_ARITHMETIC_OP_LOG10:
551 UNIARY_FUNCTION_ON_ELEMENT( log10, +0, +0); break;
552 case GAL_ARITHMETIC_OP_SIN:
553 UNIARY_FUNCTION_ON_ELEMENT( sin, *pi/180.0f, +0); break;
554 case GAL_ARITHMETIC_OP_COS:
555 UNIARY_FUNCTION_ON_ELEMENT( cos, *pi/180.0f, +0); break;
556 case GAL_ARITHMETIC_OP_TAN:
557 UNIARY_FUNCTION_ON_ELEMENT( tan, *pi/180.0f, +0); break;
558 case GAL_ARITHMETIC_OP_ASIN:
559 UNIARY_FUNCTION_ON_ELEMENT( asin, +0, *180.0f/pi); break;
560 case GAL_ARITHMETIC_OP_ACOS:
561 UNIARY_FUNCTION_ON_ELEMENT( acos, +0, *180.0f/pi); break;
562 case GAL_ARITHMETIC_OP_ATAN:
563 UNIARY_FUNCTION_ON_ELEMENT( atan, +0, *180.0f/pi); break;
564 case GAL_ARITHMETIC_OP_SINH:
565 UNIARY_FUNCTION_ON_ELEMENT( sinh, +0, +0); break;
566 case GAL_ARITHMETIC_OP_COSH:
567 UNIARY_FUNCTION_ON_ELEMENT( cosh, +0, +0); break;
568 case GAL_ARITHMETIC_OP_TANH:
569 UNIARY_FUNCTION_ON_ELEMENT( tanh, +0, +0); break;
570 case GAL_ARITHMETIC_OP_ASINH:
571 UNIARY_FUNCTION_ON_ELEMENT( asinh, +0, +0); break;
572 case GAL_ARITHMETIC_OP_ACOSH:
573 UNIARY_FUNCTION_ON_ELEMENT( acosh, +0, +0); break;
574 case GAL_ARITHMETIC_OP_ATANH:
575 UNIARY_FUNCTION_ON_ELEMENT( atanh, +0, +0); break;
576 case GAL_ARITHMETIC_OP_AU_TO_PC:
577 UNIARY_FUNCTION_ON_ELEMENT( gal_units_au_to_pc, +0, +0); break;
578 case GAL_ARITHMETIC_OP_PC_TO_AU:
579 UNIARY_FUNCTION_ON_ELEMENT( gal_units_pc_to_au, +0, +0); break;
580 case GAL_ARITHMETIC_OP_LY_TO_PC:
581 UNIARY_FUNCTION_ON_ELEMENT( gal_units_ly_to_pc, +0, +0); break;
582 case GAL_ARITHMETIC_OP_PC_TO_LY:
583 UNIARY_FUNCTION_ON_ELEMENT( gal_units_pc_to_ly, +0, +0); break;
584 case GAL_ARITHMETIC_OP_LY_TO_AU:
585 UNIARY_FUNCTION_ON_ELEMENT( gal_units_ly_to_au, +0, +0); break;
586 case GAL_ARITHMETIC_OP_AU_TO_LY:
587 UNIARY_FUNCTION_ON_ELEMENT( gal_units_au_to_ly, +0, +0); break;
588 case GAL_ARITHMETIC_OP_RA_TO_DEGREE:
589 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_STRING(double, gal_units_ra_to_degree);
590 break;
591 case GAL_ARITHMETIC_OP_DEC_TO_DEGREE:
592 UNIFUNC_RUN_FUNCTION_ON_ELEMENT_STRING(double, gal_units_dec_to_degree);
593 break;
594 case GAL_ARITHMETIC_OP_DEGREE_TO_RA:
595 UNIARY_FUNCTION_ON_ELEMENT_OUTPUT_STRING(arithmetic_units_degree_to_ra);
596 break;
597 case GAL_ARITHMETIC_OP_DEGREE_TO_DEC:
598 UNIARY_FUNCTION_ON_ELEMENT_OUTPUT_STRING(arithmetic_units_degree_to_dec);
599 break;
600 default:
601 error(EXIT_FAILURE, 0, "%s: operator code %d not recognized",
602 __func__, operator);
603 }
604
605 /* Clean up. Note that if the input arrays can be freed, and any of right
606 or left arrays needed conversion, 'UNIFUNC_CONVERT_TO_COMPILED_TYPE'
607 has already freed the input arrays, and we only have 'r' and 'l'
608 allocated in any case. Alternatively, when the inputs shouldn't be
609 freed, the only allocated spaces are the 'r' and 'l' arrays if their
610 types weren't compiled for binary operations, we can tell this from
611 the pointers: if they are different from the original pointers, they
612 were allocated. */
613 if( (flags & GAL_ARITHMETIC_FLAG_FREE) && o!=in)
614 gal_data_free(in);
615
616 /* Return */
617 return o;
618 }
619
620
621
622
623
624 /* Call functions in the 'gnuastro/statistics' library. */
625 static gal_data_t *
arithmetic_from_statistics(int operator,int flags,gal_data_t * input)626 arithmetic_from_statistics(int operator, int flags, gal_data_t *input)
627 {
628 gal_data_t *out=NULL;
629 int ip=(flags & GAL_ARITHMETIC_FLAG_INPLACE) || (flags & GAL_ARITHMETIC_FLAG_FREE);
630
631 switch(operator)
632 {
633 case GAL_ARITHMETIC_OP_MINVAL: out=gal_statistics_minimum(input);break;
634 case GAL_ARITHMETIC_OP_MAXVAL: out=gal_statistics_maximum(input);break;
635 case GAL_ARITHMETIC_OP_NUMBERVAL:out=gal_statistics_number(input); break;
636 case GAL_ARITHMETIC_OP_SUMVAL: out=gal_statistics_sum(input); break;
637 case GAL_ARITHMETIC_OP_MEANVAL: out=gal_statistics_mean(input); break;
638 case GAL_ARITHMETIC_OP_STDVAL: out=gal_statistics_std(input); break;
639 case GAL_ARITHMETIC_OP_MEDIANVAL:
640 out=gal_statistics_median(input, ip); break;
641 default:
642 error(EXIT_FAILURE, 0, "%s: operator code %d not recognized",
643 __func__, operator);
644 }
645
646 /* If the input is to be freed, then do so and return the output. */
647 if( flags & GAL_ARITHMETIC_FLAG_FREE ) gal_data_free(input);
648 return out;
649 }
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670 /***********************************************************************/
671 /*************** Adding noise **************/
672 /***********************************************************************/
673
674 /* The size operator. Reports the size along a given dimension. */
675 static gal_data_t *
arithmetic_mknoise(int operator,int flags,gal_data_t * in,gal_data_t * arg)676 arithmetic_mknoise(int operator, int flags, gal_data_t *in, gal_data_t *arg)
677 {
678 gsl_rng *rng;
679 const char *rng_name;
680 double *d, *df, arg_v;
681 unsigned long rng_seed;
682 gal_data_t *out, *targ;
683
684 /* Column counter in case '--envseed' is given and we have multiple
685 columns. */
686 static unsigned long colcounter=0;
687
688 /* Sanity checks. */
689 if(arg->size!=1)
690 error(EXIT_FAILURE, 0, "the first popped operand to the '%s' "
691 "operator should be a single number (specifying the fixed "
692 "sigma, or background level for Poisson noise), but it "
693 "has %zu elements, in %zu dimension(s)",
694 gal_arithmetic_operator_string(operator), arg->size,
695 arg->ndim);
696 if(in->type==GAL_TYPE_STRING)
697 error(EXIT_FAILURE, 0, "the input dataset to the '%s' operator "
698 "should have a numerical data type (integer or float), but "
699 "it has a string type",
700 gal_arithmetic_operator_string(operator));
701
702 /* Convert the input and argument into 'double' (and immediately free it
703 if it is no longer necessary). */
704 if(in->type==GAL_TYPE_FLOAT64) out=in;
705 else
706 {
707 out=gal_data_copy_to_new_type(in, GAL_TYPE_FLOAT64);
708 if(flags & GAL_ARITHMETIC_FLAG_FREE)
709 { gal_data_free(in); in=NULL; }
710 }
711 targ=gal_data_copy_to_new_type(arg, GAL_TYPE_FLOAT64);
712 arg_v=((double *)(targ->array))[0];
713 gal_data_free(targ);
714
715 /* Make sure the noise identifier is positive. */
716 if(arg_v<0)
717 error(EXIT_FAILURE, 0, "the noise identifier (sigma for "
718 "'mknoise-sigma', background for 'mknoise-poisson', or "
719 "range for 'mknoise-uniform') must be positive (it is %g)",
720 arg_v);
721
722 /* Setup the random number generator. For 'envseed', we want to pass a
723 boolean value: either 0 or 1. However, when we say 'flags &
724 GAL_ARITHMETIC_ENVSEED', the returned value is the integer positioning
725 of the envseed bit (for example if its on the fourth bit, the value
726 will be 8). This can cause problems if it is on the 8th bit (or any
727 multiple of 8). So to avoid issues with the bit-positioning of the
728 'ENVSEED', we will return the conditional to see if the result of the
729 bit-wise '&' is larger than 0 or not (binary). */
730 rng=gal_checkset_gsl_rng( (flags & GAL_ARITHMETIC_FLAG_ENVSEED)>0,
731 &rng_name, &rng_seed);
732
733 /* If '--envseed' was called, we need to add the column counter to the
734 requested seed. */
735 if(flags & GAL_ARITHMETIC_FLAG_ENVSEED)
736 {
737 rng_seed += colcounter++;
738 gsl_rng_set(rng, rng_seed);
739 }
740
741 /* Print the basic RNG information if requested. */
742 if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
743 {
744 printf(" - Parameters used for '%s':\n",
745 gal_arithmetic_operator_string(operator));
746 printf(" - Random number generator name: %s\n", rng_name);
747 printf(" - Random number generator seed: %lu\n", rng_seed);
748 }
749
750 /* Add the noise. */
751 df=(d=out->array)+out->size;
752 switch(operator)
753 {
754 case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
755 do *d += gsl_ran_gaussian(rng, arg_v); while(++d<df);
756 break;
757 case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
758 do
759 *d += arg_v + gsl_ran_gaussian(rng, sqrt( arg_v + *d ));
760 while(++d<df);
761 break;
762 case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
763 do
764 *d += ( (gsl_rng_uniform(rng)*arg_v) - (arg_v/2) );
765 while(++d<df);
766 break;
767 default:
768 error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
769 "fix the problem. The operator code %d isn't recognized "
770 "in this function", __func__, PACKAGE_BUGREPORT, operator);
771 }
772
773 /* Clean up and return */
774 if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(arg);
775 return out;
776 }
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795 /***********************************************************************/
796 /*************** Metadata **************/
797 /***********************************************************************/
798
799 /* The size operator. Reports the size along a given dimension. */
800 static gal_data_t *
arithmetic_size(int operator,int flags,gal_data_t * in,gal_data_t * arg)801 arithmetic_size(int operator, int flags, gal_data_t *in, gal_data_t *arg)
802 {
803 size_t one=1, arg_val;
804 gal_data_t *usearg=NULL, *out=NULL;
805
806 /* Sanity checks on argument (dimension number): it should be an integer,
807 and have a size of 1. */
808 if(arg->type==GAL_TYPE_FLOAT32 || arg->type==GAL_TYPE_FLOAT64)
809 error(EXIT_FAILURE, 0, "%s: size operator's dimension argument"
810 "must have an integer type", __func__);
811 if(arg->size!=1)
812 error(EXIT_FAILURE, 0, "%s: size operator's dimension argument"
813 "must be a single number, but it has %zu elements", __func__,
814 arg->size);
815
816
817 /* Convert 'arg' to 'size_t' and read it. Note that we can only free the
818 'arg' array (while changing its type), when the freeing flag has been
819 set. */
820 if(flags & GAL_ARITHMETIC_FLAG_FREE)
821 {
822 arg=gal_data_copy_to_new_type_free(arg, GAL_TYPE_SIZE_T);
823 arg_val=*(size_t *)(arg->array);
824 gal_data_free(arg);
825 }
826 else
827 {
828 usearg=gal_data_copy_to_new_type(arg, GAL_TYPE_SIZE_T);
829 arg_val=*(size_t *)(usearg->array);
830 gal_data_free(usearg);
831 }
832
833
834 /* Sanity checks on the value of the given argument.*/
835 if(arg_val>in->ndim)
836 error(EXIT_FAILURE, 0, "%s: size operator's dimension argument "
837 "(given %zu) cannot be larger than the dimensions of the "
838 "given input (%zu)", __func__, arg_val, in->ndim);
839 if(arg_val==0)
840 error(EXIT_FAILURE, 0, "%s: size operator's dimension argument "
841 "(given %zu) cannot be zero: dimensions are counted from 1",
842 __func__, arg_val);
843
844
845 /* Allocate the output array and write the desired dimension. Note that
846 'dsize' is in the C order, while the output must be in FITS/Fortran
847 order. Also that C order starts from 0, while the FITS order starts
848 from 1. */
849 out=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &one, NULL, 0,
850 in->minmapsize, 0, NULL, NULL, NULL);
851 *(size_t *)(out->array)=in->dsize[in->ndim-arg_val];
852
853
854 /* Clean up and return */
855 if(flags & GAL_ARITHMETIC_FLAG_FREE)
856 gal_data_free(in);
857 return out;
858 }
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879 /***********************************************************************/
880 /*************** Where **************/
881 /***********************************************************************/
882 /* When the 'iftrue' dataset only has one element and the element is blank,
883 then it will be replaced with the blank value of the type of the output
884 data. */
885 #define DO_WHERE_OPERATION(ITT, OT) { \
886 ITT *it=iftrue->array; \
887 OT b, *o=out->array, *of=o+out->size; \
888 gal_blank_write(&b, out->type); \
889 if(iftrue->size==1) \
890 { \
891 if( gal_blank_is(it, iftrue->type) ) \
892 { \
893 do{*o = (chb && *c==cb) ? *o : (*c ? b : *o); ++c; } \
894 while(++o<of); \
895 } \
896 else \
897 do {*o = (chb && *c==cb) ? *o : (*c ? *it : *o); ++c; } \
898 while(++o<of); \
899 } \
900 else \
901 do {*o = (chb && *c==cb) ? *o : (*c ? *it : *o); ++it; ++c;} \
902 while(++o<of); \
903 }
904
905
906
907
908
909 #define WHERE_OUT_SET(OT) \
910 switch(iftrue->type) \
911 { \
912 case GAL_TYPE_UINT8: DO_WHERE_OPERATION( uint8_t, OT); break; \
913 case GAL_TYPE_INT8: DO_WHERE_OPERATION( int8_t, OT); break; \
914 case GAL_TYPE_UINT16: DO_WHERE_OPERATION( uint16_t, OT); break; \
915 case GAL_TYPE_INT16: DO_WHERE_OPERATION( int16_t, OT); break; \
916 case GAL_TYPE_UINT32: DO_WHERE_OPERATION( uint32_t, OT); break; \
917 case GAL_TYPE_INT32: DO_WHERE_OPERATION( int32_t, OT); break; \
918 case GAL_TYPE_UINT64: DO_WHERE_OPERATION( uint64_t, OT); break; \
919 case GAL_TYPE_INT64: DO_WHERE_OPERATION( int64_t, OT); break; \
920 case GAL_TYPE_FLOAT32: DO_WHERE_OPERATION( float, OT); break; \
921 case GAL_TYPE_FLOAT64: DO_WHERE_OPERATION( double, OT); break; \
922 default: \
923 error(EXIT_FAILURE, 0, "%s: type code %d not recognized for the " \
924 "'iftrue' dataset", "WHERE_OUT_SET", iftrue->type); \
925 }
926
927
928
929
930
931 static void
arithmetic_where(int flags,gal_data_t * out,gal_data_t * cond,gal_data_t * iftrue)932 arithmetic_where(int flags, gal_data_t *out, gal_data_t *cond,
933 gal_data_t *iftrue)
934 {
935 int chb; /* Read as: "Condition-Has-Blank" */
936 unsigned char *c=cond->array, cb=GAL_BLANK_UINT8;
937
938 /* The condition operator has to be unsigned char. */
939 if(cond->type!=GAL_TYPE_UINT8)
940 error(EXIT_FAILURE, 0, "%s: the condition operand must be an "
941 "'uint8' type, but the given condition operand has a "
942 "'%s' type", __func__, gal_type_name(cond->type, 1));
943
944 /* The dimension and sizes of the out and condition data sets must be the
945 same. */
946 if( gal_dimension_is_different(out, cond) )
947 error(EXIT_FAILURE, 0, "%s: the output and condition datasets "
948 "must have the same size", __func__);
949
950 /* See if the condition array has blank values. */
951 chb=gal_blank_present(cond, 0);
952
953 /* Do the operation. */
954 switch(out->type)
955 {
956 case GAL_TYPE_UINT8: WHERE_OUT_SET( uint8_t ); break;
957 case GAL_TYPE_INT8: WHERE_OUT_SET( int8_t ); break;
958 case GAL_TYPE_UINT16: WHERE_OUT_SET( uint16_t ); break;
959 case GAL_TYPE_INT16: WHERE_OUT_SET( int16_t ); break;
960 case GAL_TYPE_UINT32: WHERE_OUT_SET( uint32_t ); break;
961 case GAL_TYPE_INT32: WHERE_OUT_SET( int32_t ); break;
962 case GAL_TYPE_UINT64: WHERE_OUT_SET( uint64_t ); break;
963 case GAL_TYPE_INT64: WHERE_OUT_SET( int64_t ); break;
964 case GAL_TYPE_FLOAT32: WHERE_OUT_SET( float ); break;
965 case GAL_TYPE_FLOAT64: WHERE_OUT_SET( double ); break;
966 default:
967 error(EXIT_FAILURE, 0, "%s: type code %d not recognized for the 'out'",
968 __func__, out->type);
969 }
970
971 /* Clean up if necessary. */
972 if(flags & GAL_ARITHMETIC_FLAG_FREE)
973 {
974 gal_data_free(cond);
975 gal_data_free(iftrue);
976 }
977 }
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998 /***********************************************************************/
999 /*************** Multiple operand operators **************/
1000 /***********************************************************************/
1001 struct multioperandparams
1002 {
1003 gal_data_t *list; /* List of input datasets. */
1004 gal_data_t *out; /* Output dataset. */
1005 size_t dnum; /* Number of input dataset. */
1006 int operator; /* Operator to use. */
1007 uint8_t *hasblank; /* Array of 0s or 1s for each input. */
1008 float p1; /* Sigma-cliping parameter 1. */
1009 float p2; /* Sigma-cliping parameter 2. */
1010 };
1011
1012
1013
1014
1015
1016 #define MULTIOPERAND_MIN(TYPE) { \
1017 size_t n, j=0; \
1018 TYPE t, max, *o=p->out->array; \
1019 gal_type_max(p->list->type, &max); \
1020 \
1021 /* Go over all the pixels assigned to this thread. */ \
1022 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1023 { \
1024 /* Initialize, 'j' is desired pixel's index. */ \
1025 n=0; \
1026 t=max; \
1027 j=tprm->indexs[tind]; \
1028 \
1029 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1030 { /* Only for integer types, b==b. */ \
1031 if( p->hasblank[i] && b==b) \
1032 { \
1033 if( a[i][j] != b ) \
1034 { t = a[i][j] < t ? a[i][j] : t; ++n; } \
1035 } \
1036 else { t = a[i][j] < t ? a[i][j] : t; ++n; } \
1037 } \
1038 o[j] = n ? t : b; /* No usable elements: set to blank. */ \
1039 } \
1040 }
1041
1042
1043
1044
1045
1046 #define MULTIOPERAND_MAX(TYPE) { \
1047 size_t n, j=0; \
1048 TYPE t, min, *o=p->out->array; \
1049 gal_type_min(p->list->type, &min); \
1050 \
1051 /* Go over all the pixels assigned to this thread. */ \
1052 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1053 { \
1054 /* Initialize, 'j' is desired pixel's index. */ \
1055 n=0; \
1056 t=min; \
1057 j=tprm->indexs[tind]; \
1058 \
1059 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1060 { /* Only for integer types, b==b. */ \
1061 if( p->hasblank[i] && b==b) \
1062 { \
1063 if( a[i][j] != b ) \
1064 { t = a[i][j] > t ? a[i][j] : t; ++n; } \
1065 } \
1066 else { t = a[i][j] > t ? a[i][j] : t; ++n; } \
1067 } \
1068 o[j] = n ? t : b; /* No usable elements: set to blank. */ \
1069 } \
1070 }
1071
1072
1073
1074
1075
1076 #define MULTIOPERAND_NUM { \
1077 int use; \
1078 size_t n, j; \
1079 uint32_t *o=p->out->array; \
1080 \
1081 /* Go over all the pixels assigned to this thread. */ \
1082 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1083 { \
1084 /* Initialize, 'j' is desired pixel's index. */ \
1085 n=0; \
1086 j=tprm->indexs[tind]; \
1087 \
1088 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1089 { \
1090 /* Only integers and non-NaN floats: v==v is 1. */ \
1091 if(p->hasblank[i]) \
1092 use = ( b==b \
1093 ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
1094 : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
1095 else use=1; \
1096 \
1097 /* Increment counter if necessary. */ \
1098 if(use) ++n; \
1099 } \
1100 o[j] = n; \
1101 } \
1102 }
1103
1104
1105
1106
1107
1108 #define MULTIOPERAND_SUM { \
1109 int use; \
1110 double sum; \
1111 size_t n, j; \
1112 float *o=p->out->array; \
1113 \
1114 /* Go over all the pixels assigned to this thread. */ \
1115 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1116 { \
1117 /* Initialize, 'j' is desired pixel's index. */ \
1118 n=0; \
1119 sum=0.0f; \
1120 j=tprm->indexs[tind]; \
1121 \
1122 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1123 { \
1124 /* Only integers and non-NaN floats: v==v is 1. */ \
1125 if(p->hasblank[i]) \
1126 use = ( b==b \
1127 ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
1128 : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
1129 else use=1; \
1130 \
1131 /* Use in sum if necessary. */ \
1132 if(use) { sum += a[i][j]; ++n; } \
1133 } \
1134 o[j] = n ? sum : b; \
1135 } \
1136 }
1137
1138
1139
1140
1141
1142 #define MULTIOPERAND_MEAN { \
1143 int use; \
1144 double sum; \
1145 size_t n, j; \
1146 float *o=p->out->array; \
1147 \
1148 /* Go over all the pixels assigned to this thread. */ \
1149 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1150 { \
1151 /* Initialize, 'j' is desired pixel's index. */ \
1152 n=0; \
1153 sum=0.0f; \
1154 j=tprm->indexs[tind]; \
1155 \
1156 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1157 { \
1158 /* Only integers and non-NaN floats: v==v is 1. */ \
1159 if(p->hasblank[i]) \
1160 use = ( b==b \
1161 ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
1162 : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
1163 else use=1; \
1164 \
1165 /* Calculate the mean if necessary. */ \
1166 if(use) { sum += a[i][j]; ++n; } \
1167 } \
1168 o[j] = n ? sum/n : b; \
1169 } \
1170 }
1171
1172
1173
1174
1175
1176 #define MULTIOPERAND_STD { \
1177 int use; \
1178 size_t n, j; \
1179 double sum, sum2; \
1180 float *o=p->out->array; \
1181 \
1182 /* Go over all the pixels assigned to this thread. */ \
1183 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1184 { \
1185 /* Initialize, 'j' is desired pixel's index. */ \
1186 n=0; \
1187 sum=sum2=0.0f; \
1188 j=tprm->indexs[tind]; \
1189 \
1190 for(i=0;i<p->dnum;++i) /* Loop over each array. */ \
1191 { \
1192 /* Only integers and non-NaN floats: v==v is 1. */ \
1193 if(p->hasblank[i]) \
1194 use = ( b==b \
1195 ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
1196 : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
1197 else use=1; \
1198 \
1199 /* Calculate the necessary parameters if necessary. */ \
1200 if(use) \
1201 { \
1202 sum2 += a[i][j] * a[i][j]; \
1203 sum += a[i][j]; \
1204 ++n; \
1205 } \
1206 } \
1207 o[j] = n ? sqrt( (sum2-sum*sum/n)/n ) : b; \
1208 } \
1209 }
1210
1211
1212
1213
1214
1215 #define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) { \
1216 int use; \
1217 size_t n, j; \
1218 float *o=p->out->array; \
1219 TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
1220 __func__, "pixs"); \
1221 \
1222 /* Go over all the pixels assigned to this thread. */ \
1223 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1224 { \
1225 /* Initialize, 'j' is desired pixel's index. */ \
1226 n=0; \
1227 j=tprm->indexs[tind]; \
1228 \
1229 /* Loop over each array: 'i' is input dataset's index. */ \
1230 for(i=0;i<p->dnum;++i) \
1231 { \
1232 /* Only integers and non-NaN floats: v==v is 1. */ \
1233 if(p->hasblank[i]) \
1234 use = ( b==b \
1235 ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
1236 : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
1237 else use=1; \
1238 \
1239 /* Put the value into the array of values. */ \
1240 if(use) pixs[n++]=a[i][j]; \
1241 } \
1242 \
1243 /* Sort all the values for this pixel and return the median. */ \
1244 if(n) \
1245 { \
1246 qsort(pixs, n, sizeof *pixs, QSORT_F); \
1247 o[j] = n%2 ? pixs[n/2] : (pixs[n/2] + pixs[n/2-1])/2 ; \
1248 } \
1249 else \
1250 o[j]=b; \
1251 } \
1252 \
1253 /* Clean up. */ \
1254 free(pixs); \
1255 }
1256
1257
1258
1259
1260
1261 #define MULTIOPERAND_QUANTILE(TYPE) { \
1262 size_t n, j; \
1263 gal_data_t *quantile; \
1264 TYPE *o=p->out->array; \
1265 TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
1266 __func__, "pixs"); \
1267 gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
1268 NULL, 0, -1, 1, NULL, NULL, NULL); \
1269 \
1270 /* Go over all the pixels assigned to this thread. */ \
1271 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1272 { \
1273 /* Initialize, 'j' is desired pixel's index. */ \
1274 n=0; \
1275 j=tprm->indexs[tind]; \
1276 \
1277 /* Read the necessay values from each input. */ \
1278 for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j]; \
1279 \
1280 /* If there are any elements, measure the */ \
1281 if(n) \
1282 { \
1283 /* Calculate the quantile and put it in the output. */ \
1284 quantile=gal_statistics_quantile(cont, p->p1, 1); \
1285 memcpy(&o[j], quantile->array, \
1286 gal_type_sizeof(p->list->type)); \
1287 gal_data_free(quantile); \
1288 \
1289 /* Since we are doing sigma-clipping in place, the size, */ \
1290 /* and flags need to be reset. */ \
1291 cont->flag=0; \
1292 cont->size=cont->dsize[0]=p->dnum; \
1293 } \
1294 else \
1295 o[j]=b; \
1296 } \
1297 \
1298 /* Clean up (note that 'pixs' is inside of 'cont'). */ \
1299 gal_data_free(cont); \
1300 }
1301
1302
1303
1304
1305
1306 #define MULTIOPERAND_SIGCLIP(TYPE) { \
1307 size_t n, j; \
1308 gal_data_t *sclip; \
1309 uint32_t *N=p->out->array; \
1310 float *sarr, *o=p->out->array; \
1311 TYPE *pixs=gal_pointer_allocate(p->list->type, p->dnum, 0, \
1312 __func__, "pixs"); \
1313 gal_data_t *cont=gal_data_alloc(pixs, p->list->type, 1, &p->dnum, \
1314 NULL, 0, -1, 1, NULL, NULL, NULL); \
1315 \
1316 /* Go over all the pixels assigned to this thread. */ \
1317 for(tind=0; tprm->indexs[tind] != GAL_BLANK_SIZE_T; ++tind) \
1318 { \
1319 /* Initialize, 'j' is desired pixel's index. */ \
1320 n=0; \
1321 j=tprm->indexs[tind]; \
1322 \
1323 /* Read the necessay values from each input. */ \
1324 for(i=0;i<p->dnum;++i) pixs[n++]=a[i][j]; \
1325 \
1326 /* If there are any elements, measure the */ \
1327 if(n) \
1328 { \
1329 /* Calculate the sigma-clip and write it in. */ \
1330 sclip=gal_statistics_sigma_clip(cont, p->p1, p->p2, 1, 1); \
1331 sarr=sclip->array; \
1332 switch(p->operator) \
1333 { \
1334 case GAL_ARITHMETIC_OP_SIGCLIP_STD: o[j]=sarr[3]; break;\
1335 case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: o[j]=sarr[2]; break;\
1336 case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: o[j]=sarr[1]; break;\
1337 case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: N[j]=sarr[0]; break;\
1338 default: \
1339 error(EXIT_FAILURE, 0, "%s: a bug! the code %d is not " \
1340 "valid for sigma-clipping results", __func__, \
1341 p->operator); \
1342 } \
1343 gal_data_free(sclip); \
1344 \
1345 /* Since we are doing sigma-clipping in place, the size, */ \
1346 /* and flags need to be reset. */ \
1347 cont->flag=0; \
1348 cont->size=cont->dsize[0]=p->dnum; \
1349 } \
1350 else \
1351 o[j]=b; \
1352 } \
1353 \
1354 /* Clean up (note that 'pixs' is inside of 'cont'). */ \
1355 gal_data_free(cont); \
1356 }
1357
1358
1359
1360
1361
1362 #define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) { \
1363 TYPE b, **a; \
1364 gal_data_t *tmp; \
1365 size_t i=0, tind; \
1366 /* Allocate space to keep the pointers to the arrays of each. */ \
1367 /* Input data structure. The operators will increment these */ \
1368 /* pointers while parsing them. */ \
1369 errno=0; \
1370 a=malloc(p->dnum*sizeof *a); \
1371 if(a==NULL) \
1372 error(EXIT_FAILURE, 0, "%s: %zu bytes for 'a'", \
1373 "MULTIOPERAND_TYPE_SET", p->dnum*sizeof *a); \
1374 \
1375 /* Fill in the array pointers and the blank value for this type. */ \
1376 gal_blank_write(&b, p->list->type); \
1377 for(tmp=p->list;tmp!=NULL;tmp=tmp->next) \
1378 a[i++]=tmp->array; \
1379 \
1380 /* Do the operation. */ \
1381 switch(p->operator) \
1382 { \
1383 case GAL_ARITHMETIC_OP_MIN: \
1384 MULTIOPERAND_MIN(TYPE); \
1385 break; \
1386 \
1387 case GAL_ARITHMETIC_OP_MAX: \
1388 MULTIOPERAND_MAX(TYPE); \
1389 break; \
1390 \
1391 case GAL_ARITHMETIC_OP_NUMBER: \
1392 MULTIOPERAND_NUM; \
1393 break; \
1394 \
1395 case GAL_ARITHMETIC_OP_SUM: \
1396 MULTIOPERAND_SUM; \
1397 break; \
1398 \
1399 case GAL_ARITHMETIC_OP_MEAN: \
1400 MULTIOPERAND_MEAN; \
1401 break; \
1402 \
1403 case GAL_ARITHMETIC_OP_STD: \
1404 MULTIOPERAND_STD; \
1405 break; \
1406 \
1407 case GAL_ARITHMETIC_OP_MEDIAN: \
1408 MULTIOPERAND_MEDIAN(TYPE, QSORT_F); \
1409 break; \
1410 \
1411 case GAL_ARITHMETIC_OP_QUANTILE: \
1412 MULTIOPERAND_QUANTILE(TYPE); \
1413 break; \
1414 \
1415 case GAL_ARITHMETIC_OP_SIGCLIP_STD: \
1416 case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: \
1417 case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: \
1418 case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: \
1419 MULTIOPERAND_SIGCLIP(TYPE); \
1420 break; \
1421 \
1422 default: \
1423 error(EXIT_FAILURE, 0, "%s: operator code %d not recognized", \
1424 "MULTIOPERAND_TYPE_SET", p->operator); \
1425 } \
1426 \
1427 /* Clean up. */ \
1428 free(a); \
1429 }
1430
1431
1432
1433
1434
1435 /* Worker function on each thread. */
1436 void *
multioperand_on_thread(void * in_prm)1437 multioperand_on_thread(void *in_prm)
1438 {
1439 /* Low-level definitions to be done first. */
1440 struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
1441 struct multioperandparams *p=(struct multioperandparams *)tprm->params;
1442
1443 /* Do the operation on each thread. */
1444 switch(p->list->type)
1445 {
1446 case GAL_TYPE_UINT8:
1447 MULTIOPERAND_TYPE_SET(uint8_t, gal_qsort_uint8_i);
1448 break;
1449 case GAL_TYPE_INT8:
1450 MULTIOPERAND_TYPE_SET(int8_t, gal_qsort_int8_i);
1451 break;
1452 case GAL_TYPE_UINT16:
1453 MULTIOPERAND_TYPE_SET(uint16_t, gal_qsort_uint16_i);
1454 break;
1455 case GAL_TYPE_INT16:
1456 MULTIOPERAND_TYPE_SET(int16_t, gal_qsort_int16_i);
1457 break;
1458 case GAL_TYPE_UINT32:
1459 MULTIOPERAND_TYPE_SET(uint32_t, gal_qsort_uint32_i);
1460 break;
1461 case GAL_TYPE_INT32:
1462 MULTIOPERAND_TYPE_SET(int32_t, gal_qsort_int32_i);
1463 break;
1464 case GAL_TYPE_UINT64:
1465 MULTIOPERAND_TYPE_SET(uint64_t, gal_qsort_uint64_i);
1466 break;
1467 case GAL_TYPE_INT64:
1468 MULTIOPERAND_TYPE_SET(int64_t, gal_qsort_int64_i);
1469 break;
1470 case GAL_TYPE_FLOAT32:
1471 MULTIOPERAND_TYPE_SET(float, gal_qsort_float32_i);
1472 break;
1473 case GAL_TYPE_FLOAT64:
1474 MULTIOPERAND_TYPE_SET(double, gal_qsort_float64_i);
1475 break;
1476 default:
1477 error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
1478 __func__, p->list->type);
1479 }
1480
1481 /* Wait for all the other threads to finish, then return. */
1482 if(tprm->b) pthread_barrier_wait(tprm->b);
1483 return NULL;
1484 }
1485
1486
1487
1488
1489
1490 /* The single operator in this function is assumed to be a linked list. The
1491 number of operators is determined from the fact that the last node in
1492 the linked list must have a NULL pointer as its 'next' element. */
1493 static gal_data_t *
arithmetic_multioperand(int operator,int flags,gal_data_t * list,gal_data_t * params,size_t numthreads)1494 arithmetic_multioperand(int operator, int flags, gal_data_t *list,
1495 gal_data_t *params, size_t numthreads)
1496 {
1497 size_t i=0, dnum=1;
1498 float p1=NAN, p2=NAN;
1499 struct multioperandparams p;
1500 gal_data_t *out, *tmp, *ttmp;
1501 uint8_t *hasblank, otype=GAL_TYPE_INVALID;
1502
1503
1504 /* For generality, 'list' can be a NULL pointer, in that case, this
1505 function will return a NULL pointer and avoid further processing. */
1506 if(list==NULL) return NULL;
1507
1508
1509 /* If any parameters are given, prepare them. */
1510 for(tmp=params; tmp!=NULL; tmp=tmp->next)
1511 {
1512 /* Basic sanity checks. */
1513 if(tmp->size>1)
1514 error(EXIT_FAILURE, 0, "%s: parameters must be a single number",
1515 __func__);
1516 if(tmp->type!=GAL_TYPE_FLOAT32)
1517 error(EXIT_FAILURE, 0, "%s: parameters must be float32 type",
1518 __func__);
1519
1520 /* Write them */
1521 if(isnan(p1)) p1=((float *)(tmp->array))[0];
1522 else p2=((float *)(tmp->array))[0];
1523
1524 /* Operator specific, parameter sanity checks. */
1525 switch(operator)
1526 {
1527 case GAL_ARITHMETIC_OP_QUANTILE:
1528 if(p1<0 || p1>1)
1529 error(EXIT_FAILURE, 0, "%s: the parameter given to the "
1530 "'quantile' operator must be between (and including) "
1531 "0 and 1. The given value is: %g", __func__, p1);
1532 break;
1533 }
1534 }
1535
1536
1537 /* Do a simple sanity check, comparing the operand on top of the list to
1538 the rest of the operands within the list. */
1539 for(tmp=list->next;tmp!=NULL;tmp=tmp->next)
1540 {
1541 /* Increment the number of structures. */
1542 ++dnum;
1543
1544 /* Check the types. */
1545 if(tmp->type!=list->type)
1546 error(EXIT_FAILURE, 0, "%s: the types of all operands to the '%s' "
1547 "operator must be same", __func__,
1548 gal_arithmetic_operator_string(operator));
1549
1550 /* Check the sizes. */
1551 if( gal_dimension_is_different(list, tmp) )
1552 error(EXIT_FAILURE, 0, "%s: the sizes of all operands to the '%s' "
1553 "operator must be same", __func__,
1554 gal_arithmetic_operator_string(operator));
1555 }
1556
1557
1558 /* Set the output dataset type */
1559 switch(operator)
1560 {
1561 case GAL_ARITHMETIC_OP_MIN: otype=list->type; break;
1562 case GAL_ARITHMETIC_OP_MAX: otype=list->type; break;
1563 case GAL_ARITHMETIC_OP_NUMBER: otype=GAL_TYPE_UINT32; break;
1564 case GAL_ARITHMETIC_OP_SUM: otype=GAL_TYPE_FLOAT32; break;
1565 case GAL_ARITHMETIC_OP_MEAN: otype=GAL_TYPE_FLOAT32; break;
1566 case GAL_ARITHMETIC_OP_STD: otype=GAL_TYPE_FLOAT32; break;
1567 case GAL_ARITHMETIC_OP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
1568 case GAL_ARITHMETIC_OP_QUANTILE: otype=list->type; break;
1569 case GAL_ARITHMETIC_OP_SIGCLIP_STD: otype=GAL_TYPE_FLOAT32; break;
1570 case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: otype=GAL_TYPE_FLOAT32; break;
1571 case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: otype=GAL_TYPE_FLOAT32; break;
1572 case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: otype=GAL_TYPE_UINT32; break;
1573 default:
1574 error(EXIT_FAILURE, 0, "%s: operator code %d isn't recognized",
1575 __func__, operator);
1576 }
1577
1578
1579 /* Set the output data structure. */
1580 if( (flags & GAL_ARITHMETIC_FLAG_INPLACE) && otype==list->type)
1581 out = list; /* The top element in the list. */
1582 else
1583 out = gal_data_alloc(NULL, otype, list->ndim, list->dsize,
1584 list->wcs, 0, list->minmapsize, list->quietmmap,
1585 NULL, NULL, NULL);
1586
1587
1588 /* hasblank is used to see if a blank value should be checked for each
1589 list element or not. */
1590 hasblank=gal_pointer_allocate(GAL_TYPE_UINT8, dnum, 0, __func__,
1591 "hasblank");
1592 for(tmp=list;tmp!=NULL;tmp=tmp->next)
1593 hasblank[i++]=gal_blank_present(tmp, 0);
1594
1595
1596 /* Set the parameters necessary for multithreaded operation and spin them
1597 off to do apply the operator. */
1598 p.p1=p1;
1599 p.p2=p2;
1600 p.out=out;
1601 p.list=list;
1602 p.dnum=dnum;
1603 p.operator=operator;
1604 p.hasblank=hasblank;
1605 gal_threads_spin_off(multioperand_on_thread, &p, out->size, numthreads,
1606 list->minmapsize, list->quietmmap);
1607
1608
1609 /* Clean up and return. Note that the operation might have been done in
1610 place. In that case, a list element was used. So we need to check
1611 before freeing each data structure. If we are on the designated output
1612 dataset, we should set its 'next' pointer to NULL so it isn't treated
1613 as a list any more by future functions. */
1614 if(flags & GAL_ARITHMETIC_FLAG_FREE)
1615 {
1616 tmp=list;
1617 while(tmp!=NULL)
1618 {
1619 ttmp=tmp->next;
1620 if(tmp==out) tmp->next=NULL; else gal_data_free(tmp);
1621 tmp=ttmp;
1622 }
1623 if(params) gal_list_data_free(params);
1624 }
1625 free(hasblank);
1626 return out;
1627 }
1628
1629
1630
1631
1632
1633
1634
1635
1636
1637
1638
1639
1640
1641
1642
1643
1644
1645
1646
1647
1648 /**********************************************************************/
1649 /**************** Binary operators *****************/
1650 /**********************************************************************/
1651 /* See if we should check for blanks. When both types are floats, blanks
1652 don't need to be checked (the floating point standard will do the job
1653 for us). It is also not necessary to check blanks in bitwise operators,
1654 but bitwise operators have their own macro
1655 ('BINARY_OP_INCR_OT_RT_LT_SET') which doesn' use 'checkblanks'.*/
1656 int
gal_arithmetic_binary_checkblank(gal_data_t * l,gal_data_t * r)1657 gal_arithmetic_binary_checkblank(gal_data_t *l, gal_data_t *r)
1658 {
1659 return ((((l->type!=GAL_TYPE_FLOAT32 && l->type!=GAL_TYPE_FLOAT64)
1660 || (r->type!=GAL_TYPE_FLOAT32 && r->type!=GAL_TYPE_FLOAT64))
1661 && (gal_blank_present(l, 1) || gal_blank_present(r, 1)))
1662 ? 1 : 0 );
1663 }
1664
1665
1666
1667
1668
1669 static int
arithmetic_binary_out_type(int operator,gal_data_t * l,gal_data_t * r)1670 arithmetic_binary_out_type(int operator, gal_data_t *l, gal_data_t *r)
1671 {
1672 switch(operator)
1673 {
1674 case GAL_ARITHMETIC_OP_LT:
1675 case GAL_ARITHMETIC_OP_LE:
1676 case GAL_ARITHMETIC_OP_GT:
1677 case GAL_ARITHMETIC_OP_GE:
1678 case GAL_ARITHMETIC_OP_EQ:
1679 case GAL_ARITHMETIC_OP_NE:
1680 case GAL_ARITHMETIC_OP_AND:
1681 case GAL_ARITHMETIC_OP_OR:
1682 return GAL_TYPE_UINT8;
1683
1684 default:
1685 return gal_type_out(l->type, r->type);
1686 }
1687 return -1;
1688 }
1689
1690
1691
1692
1693
1694 static gal_data_t *
arithmetic_binary(int operator,int flags,gal_data_t * l,gal_data_t * r)1695 arithmetic_binary(int operator, int flags, gal_data_t *l, gal_data_t *r)
1696 {
1697 /* Read the variable arguments. 'lo' and 'ro' keep the original data, in
1698 case their type isn't built (based on configure options are configure
1699 time). */
1700 int32_t otype;
1701 gal_data_t *o=NULL;
1702 size_t out_size, minmapsize;
1703 int quietmmap=l->quietmmap && r->quietmmap;
1704
1705
1706 /* Simple sanity check on the input sizes */
1707 if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (l->size==1 || r->size==1))
1708 && gal_dimension_is_different(l, r) )
1709 error(EXIT_FAILURE, 0, "%s: the non-number inputs to '%s' don't "
1710 "have the same dimension/size", __func__,
1711 gal_arithmetic_operator_string(operator));
1712
1713
1714 /* Set the output type. For the comparison operators, the output type is
1715 either 0 or 1, so we will set the output type to 'unsigned char' for
1716 efficient memory and CPU usage. Since the number of operators without
1717 a fixed output type (like the conditionals) is less, by 'default' we
1718 will set the output type to 'unsigned char', and if any of the other
1719 operatrs are given, it will be chosen based on the input types.*/
1720 otype=arithmetic_binary_out_type(operator, l, r);
1721
1722
1723 /* Set the output sizes. */
1724 minmapsize = ( l->minmapsize < r->minmapsize
1725 ? l->minmapsize : r->minmapsize );
1726 out_size = l->size > r->size ? l->size : r->size;
1727
1728
1729 /* If we want inplace output, set the output pointer to one input. Note
1730 that the output type can be different from both inputs. */
1731 if(flags & GAL_ARITHMETIC_FLAG_INPLACE)
1732 {
1733 if (l->type==otype && out_size==l->size) o = l;
1734 else if(r->type==otype && out_size==r->size) o = r;
1735 }
1736
1737
1738 /* If the output pointer was not set above for any of the possible
1739 reasons, allocate it. For 'mmapsize', note that since its 'size_t', it
1740 will always be positive. The '-1' that is recommended to give when you
1741 want the value in RAM is actually the largest possible memory
1742 location. So we just have to choose the smaller minmapsize of the two
1743 to decide if the output array should be in RAM or not. */
1744 if(o==NULL)
1745 o = gal_data_alloc(NULL, otype,
1746 l->size>1 ? l->ndim : r->ndim,
1747 l->size>1 ? l->dsize : r->dsize,
1748 l->size>1 ? l->wcs : r->wcs,
1749 0, minmapsize, quietmmap, NULL, NULL, NULL );
1750
1751
1752 /* Call the proper function for the operator. Since they heavily involve
1753 macros, their compilation can be very large if they are in a single
1754 function and file. So there is a separate C source and header file for
1755 each of these functions. */
1756 switch(operator)
1757 {
1758 case GAL_ARITHMETIC_OP_PLUS: arithmetic_plus(l, r, o); break;
1759 case GAL_ARITHMETIC_OP_MINUS: arithmetic_minus(l, r, o); break;
1760 case GAL_ARITHMETIC_OP_MULTIPLY: arithmetic_multiply(l, r, o); break;
1761 case GAL_ARITHMETIC_OP_DIVIDE: arithmetic_divide(l, r, o); break;
1762 case GAL_ARITHMETIC_OP_LT: arithmetic_lt(l, r, o); break;
1763 case GAL_ARITHMETIC_OP_LE: arithmetic_le(l, r, o); break;
1764 case GAL_ARITHMETIC_OP_GT: arithmetic_gt(l, r, o); break;
1765 case GAL_ARITHMETIC_OP_GE: arithmetic_ge(l, r, o); break;
1766 case GAL_ARITHMETIC_OP_EQ: arithmetic_eq(l, r, o); break;
1767 case GAL_ARITHMETIC_OP_NE: arithmetic_ne(l, r, o); break;
1768 case GAL_ARITHMETIC_OP_AND: arithmetic_and(l, r, o); break;
1769 case GAL_ARITHMETIC_OP_OR: arithmetic_or(l, r, o); break;
1770 case GAL_ARITHMETIC_OP_BITAND: arithmetic_bitand(l, r, o); break;
1771 case GAL_ARITHMETIC_OP_BITOR: arithmetic_bitor(l, r, o); break;
1772 case GAL_ARITHMETIC_OP_BITXOR: arithmetic_bitxor(l, r, o); break;
1773 case GAL_ARITHMETIC_OP_BITLSH: arithmetic_bitlsh(l, r, o); break;
1774 case GAL_ARITHMETIC_OP_BITRSH: arithmetic_bitrsh(l, r, o); break;
1775 case GAL_ARITHMETIC_OP_MODULO: arithmetic_modulo(l, r, o); break;
1776 default:
1777 error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to address "
1778 "the problem. %d is not a valid operator code", __func__,
1779 PACKAGE_BUGREPORT, operator);
1780 }
1781
1782
1783 /* Clean up if necessary. Note that if the operation was requested to be
1784 in place, then the output might be one of the inputs. */
1785 if(flags & GAL_ARITHMETIC_FLAG_FREE)
1786 {
1787 if (o==l) gal_data_free(r);
1788 else if(o==r) gal_data_free(l);
1789 else { gal_data_free(l); gal_data_free(r); }
1790 }
1791
1792
1793 /* Return */
1794 return o;
1795 }
1796
1797
1798
1799
1800
1801 #define BINFUNC_RUN_FUNCTION(OT, RT, LT, OP, AFTER){ \
1802 LT *la=l->array; \
1803 RT *ra=r->array; \
1804 OT *oa=o->array, *of=oa + o->size; \
1805 if(l->size==r->size) do *oa = OP(*la++, *ra++) AFTER; while(++oa<of); \
1806 else if(l->size==1) do *oa = OP(*la, *ra++) AFTER; while(++oa<of); \
1807 else do *oa = OP(*la++, *ra ) AFTER; while(++oa<of); \
1808 }
1809
1810
1811 #define BINFUNC_F_OPERATOR_LEFT_RIGHT_SET(RT, LT, OP, AFTER) \
1812 switch(o->type) \
1813 { \
1814 case GAL_TYPE_FLOAT32: \
1815 BINFUNC_RUN_FUNCTION(float, RT, LT, OP, AFTER); \
1816 break; \
1817 case GAL_TYPE_FLOAT64: \
1818 BINFUNC_RUN_FUNCTION(double, RT, LT, OP, AFTER); \
1819 break; \
1820 default: \
1821 error(EXIT_FAILURE, 0, "%s: type %d not recognized for o->type ", \
1822 "BINFUNC_F_OPERATOR_LEFT_RIGHT_SET", o->type); \
1823 }
1824
1825
1826 #define BINFUNC_F_OPERATOR_LEFT_SET(LT, OP, AFTER) \
1827 switch(r->type) \
1828 { \
1829 case GAL_TYPE_FLOAT32: \
1830 BINFUNC_F_OPERATOR_LEFT_RIGHT_SET(float, LT, OP, AFTER); \
1831 break; \
1832 case GAL_TYPE_FLOAT64: \
1833 BINFUNC_F_OPERATOR_LEFT_RIGHT_SET(double, LT, OP, AFTER); \
1834 break; \
1835 default: \
1836 error(EXIT_FAILURE, 0, "%s: type %d not recognized for r->type", \
1837 "BINFUNC_F_OPERATOR_LEFT_SET", r->type); \
1838 }
1839
1840
1841 #define BINFUNC_F_OPERATOR_SET(OP, AFTER) \
1842 switch(l->type) \
1843 { \
1844 case GAL_TYPE_FLOAT32: \
1845 BINFUNC_F_OPERATOR_LEFT_SET(float, OP, AFTER); \
1846 break; \
1847 case GAL_TYPE_FLOAT64: \
1848 BINFUNC_F_OPERATOR_LEFT_SET(double, OP, AFTER); \
1849 break; \
1850 default: \
1851 error(EXIT_FAILURE, 0, "%s: type %d not recognized for l->type", \
1852 "BINFUNC_F_OPERATOR_SET", l->type); \
1853 }
1854
1855
1856 static gal_data_t *
arithmetic_function_binary_flt(int operator,int flags,gal_data_t * il,gal_data_t * ir)1857 arithmetic_function_binary_flt(int operator, int flags, gal_data_t *il,
1858 gal_data_t *ir)
1859 {
1860 int final_otype;
1861 size_t out_size, minmapsize;
1862 gal_data_t *l, *r, *o=NULL;
1863 double pi=3.14159265358979323846264338327;
1864 int quietmmap=il->quietmmap && ir->quietmmap;
1865
1866 /* Simple sanity check on the input sizes */
1867 if( !( (flags & GAL_ARITHMETIC_FLAG_NUMOK) && (il->size==1 || ir->size==1))
1868 && gal_dimension_is_different(il, ir) )
1869 error(EXIT_FAILURE, 0, "%s: the input datasets don't have the same "
1870 "dimension/size", __func__);
1871
1872
1873 /* Convert the values to double precision floating point if they are
1874 integer. */
1875 l = ( (il->type==GAL_TYPE_FLOAT32 || il->type==GAL_TYPE_FLOAT64)
1876 ? il : gal_data_copy_to_new_type(il, GAL_TYPE_FLOAT64) );
1877 r = ( (ir->type==GAL_TYPE_FLOAT32 || ir->type==GAL_TYPE_FLOAT64)
1878 ? ir : gal_data_copy_to_new_type(ir, GAL_TYPE_FLOAT64) );
1879
1880
1881 /* Set the output type. */
1882 final_otype = gal_type_out(l->type, r->type);
1883
1884
1885 /* Set the output sizes. */
1886 minmapsize = ( l->minmapsize < r->minmapsize
1887 ? l->minmapsize
1888 : r->minmapsize );
1889 out_size = l->size > r->size ? l->size : r->size;
1890
1891
1892 /* If we want inplace output, set the output pointer to one input. Note
1893 that the output type can be different from both inputs. */
1894 if(flags & GAL_ARITHMETIC_FLAG_INPLACE)
1895 {
1896 if (l->type==final_otype && out_size==l->size) o = l;
1897 else if(r->type==final_otype && out_size==r->size) o = r;
1898 }
1899
1900
1901 /* If the output pointer was not set for any reason, allocate it. For
1902 'mmapsize', note that since its 'size_t', it will always be
1903 Positive. The '-1' that is recommended to give when you want the value
1904 in RAM is actually the largest possible memory location. So we just
1905 have to choose the smaller minmapsize of the two to decide if the
1906 output array should be in RAM or not. */
1907 if(o==NULL)
1908 o = gal_data_alloc(NULL, final_otype,
1909 l->size>1 ? l->ndim : r->ndim,
1910 l->size>1 ? l->dsize : r->dsize,
1911 l->size>1 ? l->wcs : r->wcs, 0, minmapsize,
1912 quietmmap, NULL, NULL, NULL);
1913
1914
1915 /* Start setting the operator and operands. */
1916 switch(operator)
1917 {
1918 case GAL_ARITHMETIC_OP_POW:
1919 BINFUNC_F_OPERATOR_SET( pow, +0 ); break;
1920 case GAL_ARITHMETIC_OP_ATAN2:
1921 BINFUNC_F_OPERATOR_SET( atan2, *180.0f/pi ); break;
1922 case GAL_ARITHMETIC_OP_COUNTS_TO_MAG:
1923 BINFUNC_F_OPERATOR_SET( gal_units_counts_to_mag, +0 ); break;
1924 case GAL_ARITHMETIC_OP_MAG_TO_COUNTS:
1925 BINFUNC_F_OPERATOR_SET( gal_units_mag_to_counts, +0 ); break;
1926 case GAL_ARITHMETIC_OP_COUNTS_TO_JY:
1927 BINFUNC_F_OPERATOR_SET( gal_units_counts_to_jy, +0 ); break;
1928 default:
1929 error(EXIT_FAILURE, 0, "%s: operator code %d not recognized",
1930 __func__, operator);
1931 }
1932
1933
1934 /* Clean up. Note that if the input arrays can be freed, and any of right
1935 or left arrays needed conversion, 'BINFUNC_CONVERT_TO_COMPILED_TYPE'
1936 has already freed the input arrays, and we only have 'r' and 'l'
1937 allocated in any case. Alternatively, when the inputs shouldn't be
1938 freed, the only allocated spaces are the 'r' and 'l' arrays if their
1939 types weren't compiled for binary operations, we can tell this from
1940 the pointers: if they are different from the original pointers, they
1941 were allocated. */
1942 if(flags & GAL_ARITHMETIC_FLAG_FREE)
1943 {
1944 /* Clean the main used (temporarily allocated) datasets. */
1945 if (o==l) gal_data_free(r);
1946 else if(o==r) gal_data_free(l);
1947 else { gal_data_free(l); gal_data_free(r); }
1948
1949 /* Clean the raw inputs, if they weren't equal to the datasets. */
1950 if (o==il) { if(ir!=r) gal_data_free(ir); }
1951 else if(o==ir) { if(il!=l) gal_data_free(il); }
1952 else { if(il!=l) gal_data_free(il);
1953 if(ir!=r) gal_data_free(ir); }
1954 }
1955 else
1956 {
1957 /* Input datasets should be kept, but we don't want the temporary
1958 datasets, so if they were allocated (they don't equal the input
1959 pointers, free them). */
1960 if (l!=il) gal_data_free(l);
1961 if (r!=ir) gal_data_free(r);
1962 }
1963
1964 /* Return */
1965 return o;
1966 }
1967
1968
1969
1970
1971
1972 gal_data_t *
arithmetic_box_around_ellipse(gal_data_t * d1,gal_data_t * d2,gal_data_t * d3,int operator,int flags)1973 arithmetic_box_around_ellipse(gal_data_t *d1, gal_data_t *d2,
1974 gal_data_t *d3, int operator, int flags)
1975 {
1976 size_t i;
1977 double *a, *b, *pa, out[2];
1978 gal_data_t *a_data, *b_data, *pa_data;
1979
1980 /* Basic sanity check. */
1981 if( d1->size != d2->size || d1->size != d3->size )
1982 error(EXIT_FAILURE, 0, "%s: the three operands to this "
1983 "function don't have the same number of elements",
1984 __func__);
1985
1986 /* Convert the two inputs into double. Note that if the user doesn't want
1987 to free the inputs, we should make a copy of 'a_data' and 'b_data'
1988 because the output will also be written in them. */
1989 a_data=( ( d1->type==GAL_TYPE_FLOAT64
1990 || flags & GAL_ARITHMETIC_FLAG_FREE )
1991 ? d1
1992 : gal_data_copy_to_new_type(d1, GAL_TYPE_FLOAT64) );
1993 b_data=( ( d2->type==GAL_TYPE_FLOAT64
1994 || flags & GAL_ARITHMETIC_FLAG_FREE )
1995 ? d2
1996 : gal_data_copy_to_new_type(d2, GAL_TYPE_FLOAT64) );
1997 pa_data=( d3->type==GAL_TYPE_FLOAT64
1998 ? d3
1999 : gal_data_copy_to_new_type(d3, GAL_TYPE_FLOAT64) );
2000
2001 /* Set the double pointers and do the calculation for each. */
2002 a=a_data->array;
2003 b=b_data->array;
2004 pa=pa_data->array;
2005 for(i=0;i<a_data->size;++i)
2006 {
2007 /* If the minor axis has a larger value, print a warning and set the
2008 value to NaN. */
2009 if(b[i]>a[i])
2010 {
2011 if( (flags & GAL_ARITHMETIC_FLAG_QUIET)==0 )
2012 error(EXIT_SUCCESS, 0, "%s: the minor axis (%g) is larger "
2013 "than the major axis (%g), output for this element "
2014 "will be NaN", __func__, b[i], a[i]);
2015 out[0]=out[1]=NAN;
2016 }
2017 else gal_box_bound_ellipse_extent(a[i], b[i], pa[i], out);
2018
2019 /* Write the output in the inputs (we don't need the inputs any
2020 more). */
2021 a[i]=out[0];
2022 b[i]=out[1];
2023 }
2024
2025 /* Clean up. */
2026 gal_data_free(pa_data);
2027 if(flags & GAL_ARITHMETIC_FLAG_FREE)
2028 {
2029 if(a_data !=d1) gal_data_free(d1);
2030 if(b_data !=d2) gal_data_free(d2);
2031 if(pa_data!=d3) gal_data_free(d3);
2032 }
2033
2034 /* Make the output as a list and return it. */
2035 b_data->next=a_data;
2036 return b_data;
2037 }
2038
2039
2040
2041
2042
2043 /* Make a new dataset. */
2044 gal_data_t *
arithmetic_makenew(gal_data_t * sizes)2045 arithmetic_makenew(gal_data_t *sizes)
2046 {
2047 gal_data_t *out, *tmp, *ttmp;
2048 int quietmmap=sizes->quietmmap;
2049 size_t minmapsize=sizes->minmapsize;
2050 size_t i, *dsize, ndim=gal_list_data_number(sizes);
2051
2052 /* Make sure all the elements are a single, integer number. */
2053 for(tmp=sizes; tmp!=NULL; tmp=tmp->next)
2054 {
2055 if(tmp->size!=1)
2056 error(EXIT_FAILURE, 0, "%s: operands given to 'makenew' operator "
2057 "should only be a single number. However, at least one of "
2058 "the input operands has %zu elements", __func__, tmp->size);
2059
2060 if( tmp->type==GAL_TYPE_FLOAT32 || tmp->type==GAL_TYPE_FLOAT64)
2061 error(EXIT_FAILURE, 0, "%s: operands given to 'makenew' operator "
2062 "should have integer types. However, at least one of "
2063 "the input operands is floating point", __func__);
2064 }
2065
2066 /* Fill the 'dsize' array based on the given values. */
2067 i=ndim-1;
2068 tmp=sizes;
2069 dsize=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 1, __func__, "dsize");
2070 while(tmp!=NULL)
2071 {
2072 /* Set the next pointer and conver this one to size_t. */
2073 ttmp=tmp->next;
2074 tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_SIZE_T);
2075
2076 /* Write the dimension's length into 'dsize'. */
2077 dsize[i--] = ((size_t *)(tmp->array))[0];
2078
2079 /* Free 'tmp' and re-set it to the next element. */
2080 free(tmp);
2081 tmp=ttmp;
2082 }
2083
2084 /* allocate the necessary dataset. */
2085 out=gal_data_alloc(NULL, GAL_TYPE_UINT8, ndim, dsize, NULL, 1, minmapsize,
2086 quietmmap, "EMPTY", "NOT-SET",
2087 "Empty dataset created by arithmetic.");
2088
2089 /* Clean up and return. */
2090 free(dsize);
2091 return out;
2092 }
2093
2094
2095
2096
2097
2098
2099
2100
2101
2102
2103
2104
2105
2106
2107
2108
2109
2110
2111
2112 /**********************************************************************/
2113 /**************** High-level functions *****************/
2114 /**********************************************************************/
2115 /* Order is the same as in the manual. */
2116 int
gal_arithmetic_set_operator(char * string,size_t * num_operands)2117 gal_arithmetic_set_operator(char *string, size_t *num_operands)
2118 {
2119 int op;
2120
2121 /* Simple arithmetic operators. */
2122 if (!strcmp(string, "+" ))
2123 { op=GAL_ARITHMETIC_OP_PLUS; *num_operands=2; }
2124 else if (!strcmp(string, "-" ))
2125 { op=GAL_ARITHMETIC_OP_MINUS; *num_operands=2; }
2126 else if (!strcmp(string, "x" ))
2127 { op=GAL_ARITHMETIC_OP_MULTIPLY; *num_operands=2; }
2128 else if (!strcmp(string, "/" ))
2129 { op=GAL_ARITHMETIC_OP_DIVIDE; *num_operands=2; }
2130 else if (!strcmp(string, "%" ))
2131 { op=GAL_ARITHMETIC_OP_MODULO; *num_operands=2; }
2132
2133 /* Mathematical Operators. */
2134 else if (!strcmp(string, "abs"))
2135 { op=GAL_ARITHMETIC_OP_ABS; *num_operands=1; }
2136 else if (!strcmp(string, "pow"))
2137 { op=GAL_ARITHMETIC_OP_POW; *num_operands=2; }
2138 else if (!strcmp(string, "sqrt"))
2139 { op=GAL_ARITHMETIC_OP_SQRT; *num_operands=1; }
2140 else if (!strcmp(string, "log"))
2141 { op=GAL_ARITHMETIC_OP_LOG; *num_operands=1; }
2142 else if (!strcmp(string, "log10"))
2143 { op=GAL_ARITHMETIC_OP_LOG10; *num_operands=1; }
2144
2145 /* Trigonometric functions. */
2146 else if( !strcmp(string, "sin"))
2147 { op=GAL_ARITHMETIC_OP_SIN; *num_operands=1; }
2148 else if( !strcmp(string, "cos"))
2149 { op=GAL_ARITHMETIC_OP_COS; *num_operands=1; }
2150 else if( !strcmp(string, "tan"))
2151 { op=GAL_ARITHMETIC_OP_TAN; *num_operands=1; }
2152 else if( !strcmp(string, "asin"))
2153 { op=GAL_ARITHMETIC_OP_ASIN; *num_operands=1; }
2154 else if( !strcmp(string, "acos"))
2155 { op=GAL_ARITHMETIC_OP_ACOS; *num_operands=1; }
2156 else if( !strcmp(string, "atan"))
2157 { op=GAL_ARITHMETIC_OP_ATAN; *num_operands=1; }
2158 else if( !strcmp(string, "atan2"))
2159 { op=GAL_ARITHMETIC_OP_ATAN2; *num_operands=2; }
2160 else if( !strcmp(string, "sinh"))
2161 { op=GAL_ARITHMETIC_OP_SINH; *num_operands=1; }
2162 else if( !strcmp(string, "cosh"))
2163 { op=GAL_ARITHMETIC_OP_COSH; *num_operands=1; }
2164 else if( !strcmp(string, "tanh"))
2165 { op=GAL_ARITHMETIC_OP_TANH; *num_operands=1; }
2166 else if( !strcmp(string, "asinh"))
2167 { op=GAL_ARITHMETIC_OP_ASINH; *num_operands=1; }
2168 else if( !strcmp(string, "acosh"))
2169 { op=GAL_ARITHMETIC_OP_ACOSH; *num_operands=1; }
2170 else if( !strcmp(string, "atanh"))
2171 { op=GAL_ARITHMETIC_OP_ATANH; *num_operands=1; }
2172
2173 /* Units conversion functions */
2174 else if (!strcmp(string, "ra-to-degree"))
2175 { op=GAL_ARITHMETIC_OP_RA_TO_DEGREE; *num_operands=1; }
2176 else if (!strcmp(string, "dec-to-degree"))
2177 { op=GAL_ARITHMETIC_OP_DEC_TO_DEGREE; *num_operands=1; }
2178 else if (!strcmp(string, "degree-to-ra"))
2179 { op=GAL_ARITHMETIC_OP_DEGREE_TO_RA; *num_operands=1; }
2180 else if (!strcmp(string, "degree-to-dec"))
2181 { op=GAL_ARITHMETIC_OP_DEGREE_TO_DEC; *num_operands=1; }
2182 else if (!strcmp(string, "counts-to-mag"))
2183 { op=GAL_ARITHMETIC_OP_COUNTS_TO_MAG; *num_operands=2; }
2184 else if (!strcmp(string, "mag-to-counts"))
2185 { op=GAL_ARITHMETIC_OP_MAG_TO_COUNTS; *num_operands=2; }
2186 else if (!strcmp(string, "counts-to-jy"))
2187 { op=GAL_ARITHMETIC_OP_COUNTS_TO_JY; *num_operands=2; }
2188 else if( !strcmp(string, "au-to-pc"))
2189 { op=GAL_ARITHMETIC_OP_AU_TO_PC; *num_operands=1; }
2190 else if( !strcmp(string, "pc-to-au"))
2191 { op=GAL_ARITHMETIC_OP_PC_TO_AU; *num_operands=1; }
2192 else if( !strcmp(string, "ly-to-pc"))
2193 { op=GAL_ARITHMETIC_OP_LY_TO_PC; *num_operands=1; }
2194 else if( !strcmp(string, "pc-to-ly"))
2195 { op=GAL_ARITHMETIC_OP_PC_TO_LY; *num_operands=1; }
2196 else if( !strcmp(string, "ly-to-au"))
2197 { op=GAL_ARITHMETIC_OP_LY_TO_AU; *num_operands=1; }
2198 else if( !strcmp(string, "au-to-ly"))
2199 { op=GAL_ARITHMETIC_OP_AU_TO_LY; *num_operands=1; }
2200
2201 /* Statistical/higher-level operators. */
2202 else if (!strcmp(string, "minvalue"))
2203 { op=GAL_ARITHMETIC_OP_MINVAL; *num_operands=1; }
2204 else if (!strcmp(string, "maxvalue"))
2205 { op=GAL_ARITHMETIC_OP_MAXVAL; *num_operands=1; }
2206 else if (!strcmp(string, "numbervalue"))
2207 { op=GAL_ARITHMETIC_OP_NUMBERVAL; *num_operands=1; }
2208 else if (!strcmp(string, "sumvalue"))
2209 { op=GAL_ARITHMETIC_OP_SUMVAL; *num_operands=1; }
2210 else if (!strcmp(string, "meanvalue"))
2211 { op=GAL_ARITHMETIC_OP_MEANVAL; *num_operands=1; }
2212 else if (!strcmp(string, "stdvalue"))
2213 { op=GAL_ARITHMETIC_OP_STDVAL; *num_operands=1; }
2214 else if (!strcmp(string, "medianvalue"))
2215 { op=GAL_ARITHMETIC_OP_MEDIANVAL; *num_operands=1; }
2216 else if (!strcmp(string, "min"))
2217 { op=GAL_ARITHMETIC_OP_MIN; *num_operands=-1; }
2218 else if (!strcmp(string, "max"))
2219 { op=GAL_ARITHMETIC_OP_MAX; *num_operands=-1; }
2220 else if (!strcmp(string, "number"))
2221 { op=GAL_ARITHMETIC_OP_NUMBER; *num_operands=-1; }
2222 else if (!strcmp(string, "sum"))
2223 { op=GAL_ARITHMETIC_OP_SUM; *num_operands=-1; }
2224 else if (!strcmp(string, "mean"))
2225 { op=GAL_ARITHMETIC_OP_MEAN; *num_operands=-1; }
2226 else if (!strcmp(string, "std"))
2227 { op=GAL_ARITHMETIC_OP_STD; *num_operands=-1; }
2228 else if (!strcmp(string, "median"))
2229 { op=GAL_ARITHMETIC_OP_MEDIAN; *num_operands=-1; }
2230 else if (!strcmp(string, "quantile"))
2231 { op=GAL_ARITHMETIC_OP_QUANTILE; *num_operands=-1; }
2232 else if (!strcmp(string, "sigclip-number"))
2233 { op=GAL_ARITHMETIC_OP_SIGCLIP_NUMBER; *num_operands=-1; }
2234 else if (!strcmp(string, "sigclip-mean"))
2235 { op=GAL_ARITHMETIC_OP_SIGCLIP_MEAN; *num_operands=-1; }
2236 else if (!strcmp(string, "sigclip-median"))
2237 { op=GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN; *num_operands=-1; }
2238 else if (!strcmp(string, "sigclip-std"))
2239 { op=GAL_ARITHMETIC_OP_SIGCLIP_STD; *num_operands=-1; }
2240
2241 /* Adding noise operators. */
2242 else if (!strcmp(string, "mknoise-sigma"))
2243 { op=GAL_ARITHMETIC_OP_MKNOISE_SIGMA; *num_operands=2; }
2244 else if (!strcmp(string, "mknoise-poisson"))
2245 { op=GAL_ARITHMETIC_OP_MKNOISE_POISSON; *num_operands=2; }
2246 else if (!strcmp(string, "mknoise-uniform"))
2247 { op=GAL_ARITHMETIC_OP_MKNOISE_UNIFORM; *num_operands=2; }
2248
2249 /* The size operator */
2250 else if (!strcmp(string, "size"))
2251 { op=GAL_ARITHMETIC_OP_SIZE; *num_operands=2; }
2252
2253 /* Conditional operators. */
2254 else if (!strcmp(string, "lt" ))
2255 { op=GAL_ARITHMETIC_OP_LT; *num_operands=2; }
2256 else if (!strcmp(string, "le"))
2257 { op=GAL_ARITHMETIC_OP_LE; *num_operands=2; }
2258 else if (!strcmp(string, "gt" ))
2259 { op=GAL_ARITHMETIC_OP_GT; *num_operands=2; }
2260 else if (!strcmp(string, "ge"))
2261 { op=GAL_ARITHMETIC_OP_GE; *num_operands=2; }
2262 else if (!strcmp(string, "eq"))
2263 { op=GAL_ARITHMETIC_OP_EQ; *num_operands=2; }
2264 else if (!strcmp(string, "ne"))
2265 { op=GAL_ARITHMETIC_OP_NE; *num_operands=2; }
2266 else if (!strcmp(string, "and"))
2267 { op=GAL_ARITHMETIC_OP_AND; *num_operands=2; }
2268 else if (!strcmp(string, "or"))
2269 { op=GAL_ARITHMETIC_OP_OR; *num_operands=2; }
2270 else if (!strcmp(string, "not"))
2271 { op=GAL_ARITHMETIC_OP_NOT; *num_operands=1; }
2272 else if (!strcmp(string, "isblank"))
2273 { op=GAL_ARITHMETIC_OP_ISBLANK; *num_operands=1; }
2274 else if (!strcmp(string, "where"))
2275 { op=GAL_ARITHMETIC_OP_WHERE; *num_operands=3; }
2276
2277 /* Bitwise operators. */
2278 else if (!strcmp(string, "bitand"))
2279 { op=GAL_ARITHMETIC_OP_BITAND; *num_operands=2; }
2280 else if (!strcmp(string, "bitor"))
2281 { op=GAL_ARITHMETIC_OP_BITOR; *num_operands=2; }
2282 else if (!strcmp(string, "bitxor"))
2283 { op=GAL_ARITHMETIC_OP_BITXOR; *num_operands=2; }
2284 else if (!strcmp(string, "lshift"))
2285 { op=GAL_ARITHMETIC_OP_BITLSH; *num_operands=2; }
2286 else if (!strcmp(string, "rshift"))
2287 { op=GAL_ARITHMETIC_OP_BITRSH; *num_operands=2; }
2288 else if (!strcmp(string, "bitnot"))
2289 { op=GAL_ARITHMETIC_OP_BITNOT; *num_operands=1; }
2290
2291 /* Type conversion. */
2292 else if (!strcmp(string, "uint8"))
2293 { op=GAL_ARITHMETIC_OP_TO_UINT8; *num_operands=1; }
2294 else if (!strcmp(string, "int8"))
2295 { op=GAL_ARITHMETIC_OP_TO_INT8; *num_operands=1; }
2296 else if (!strcmp(string, "uint16"))
2297 { op=GAL_ARITHMETIC_OP_TO_UINT16; *num_operands=1; }
2298 else if (!strcmp(string, "int16"))
2299 { op=GAL_ARITHMETIC_OP_TO_INT16; *num_operands=1; }
2300 else if (!strcmp(string, "uint32"))
2301 { op=GAL_ARITHMETIC_OP_TO_UINT32; *num_operands=1; }
2302 else if (!strcmp(string, "int32"))
2303 { op=GAL_ARITHMETIC_OP_TO_INT32; *num_operands=1; }
2304 else if (!strcmp(string, "uint64"))
2305 { op=GAL_ARITHMETIC_OP_TO_UINT64; *num_operands=1; }
2306 else if (!strcmp(string, "int64"))
2307 { op=GAL_ARITHMETIC_OP_TO_INT64; *num_operands=1; }
2308 else if (!strcmp(string, "float32"))
2309 { op=GAL_ARITHMETIC_OP_TO_FLOAT32; *num_operands=1; }
2310 else if (!strcmp(string, "float64"))
2311 { op=GAL_ARITHMETIC_OP_TO_FLOAT64; *num_operands=1; }
2312
2313 /* Surrounding box. */
2314 else if (!strcmp(string, "box-around-ellipse"))
2315 { op=GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE;*num_operands=3; }
2316
2317 /* New dataset. */
2318 else if (!strcmp(string, "makenew"))
2319 { op=GAL_ARITHMETIC_OP_MAKENEW; *num_operands=-1; }
2320
2321 /* Operator not defined. */
2322 else
2323 { op=GAL_ARITHMETIC_OP_INVALID; *num_operands=GAL_BLANK_INT; }
2324 return op;
2325 }
2326
2327
2328
2329
2330
2331 char *
gal_arithmetic_operator_string(int operator)2332 gal_arithmetic_operator_string(int operator)
2333 {
2334 switch(operator)
2335 {
2336 case GAL_ARITHMETIC_OP_PLUS: return "+";
2337 case GAL_ARITHMETIC_OP_MINUS: return "-";
2338 case GAL_ARITHMETIC_OP_MULTIPLY: return "x";
2339 case GAL_ARITHMETIC_OP_DIVIDE: return "/";
2340 case GAL_ARITHMETIC_OP_MODULO: return "%";
2341
2342 case GAL_ARITHMETIC_OP_LT: return "lt";
2343 case GAL_ARITHMETIC_OP_LE: return "le";
2344 case GAL_ARITHMETIC_OP_GT: return "gt";
2345 case GAL_ARITHMETIC_OP_GE: return "ge";
2346 case GAL_ARITHMETIC_OP_EQ: return "eq";
2347 case GAL_ARITHMETIC_OP_NE: return "ne";
2348 case GAL_ARITHMETIC_OP_AND: return "and";
2349 case GAL_ARITHMETIC_OP_OR: return "or";
2350 case GAL_ARITHMETIC_OP_NOT: return "not";
2351 case GAL_ARITHMETIC_OP_ISBLANK: return "isblank";
2352 case GAL_ARITHMETIC_OP_WHERE: return "where";
2353
2354 case GAL_ARITHMETIC_OP_BITAND: return "bitand";
2355 case GAL_ARITHMETIC_OP_BITOR: return "bitor";
2356 case GAL_ARITHMETIC_OP_BITXOR: return "bitxor";
2357 case GAL_ARITHMETIC_OP_BITLSH: return "lshift";
2358 case GAL_ARITHMETIC_OP_BITRSH: return "rshift";
2359 case GAL_ARITHMETIC_OP_BITNOT: return "bitnot";
2360
2361 case GAL_ARITHMETIC_OP_ABS: return "abs";
2362 case GAL_ARITHMETIC_OP_POW: return "pow";
2363 case GAL_ARITHMETIC_OP_SQRT: return "sqrt";
2364 case GAL_ARITHMETIC_OP_LOG: return "log";
2365 case GAL_ARITHMETIC_OP_LOG10: return "log10";
2366
2367 case GAL_ARITHMETIC_OP_SIN: return "sin";
2368 case GAL_ARITHMETIC_OP_COS: return "cos";
2369 case GAL_ARITHMETIC_OP_TAN: return "tan";
2370 case GAL_ARITHMETIC_OP_ASIN: return "asin";
2371 case GAL_ARITHMETIC_OP_ACOS: return "acos";
2372 case GAL_ARITHMETIC_OP_ATAN: return "atan";
2373 case GAL_ARITHMETIC_OP_SINH: return "sinh";
2374 case GAL_ARITHMETIC_OP_COSH: return "cosh";
2375 case GAL_ARITHMETIC_OP_TANH: return "tanh";
2376 case GAL_ARITHMETIC_OP_ASINH: return "asinh";
2377 case GAL_ARITHMETIC_OP_ACOSH: return "acosh";
2378 case GAL_ARITHMETIC_OP_ATANH: return "atanh";
2379 case GAL_ARITHMETIC_OP_ATAN2: return "atan2";
2380
2381 case GAL_ARITHMETIC_OP_RA_TO_DEGREE: return "ra-to-degree";
2382 case GAL_ARITHMETIC_OP_DEC_TO_DEGREE: return "dec-to-degree";
2383 case GAL_ARITHMETIC_OP_DEGREE_TO_RA: return "degree-to-ra";
2384 case GAL_ARITHMETIC_OP_DEGREE_TO_DEC: return "degree-to-dec";
2385 case GAL_ARITHMETIC_OP_COUNTS_TO_MAG: return "counts-to-mag";
2386 case GAL_ARITHMETIC_OP_MAG_TO_COUNTS: return "mag-to-counts";
2387 case GAL_ARITHMETIC_OP_COUNTS_TO_JY: return "counts-to-jy";
2388 case GAL_ARITHMETIC_OP_AU_TO_PC: return "au-to-pc";
2389 case GAL_ARITHMETIC_OP_PC_TO_AU: return "pc-to-au";
2390 case GAL_ARITHMETIC_OP_LY_TO_PC: return "ly-to-pc";
2391 case GAL_ARITHMETIC_OP_PC_TO_LY: return "pc-to-ly";
2392 case GAL_ARITHMETIC_OP_LY_TO_AU: return "ly-to-au";
2393 case GAL_ARITHMETIC_OP_AU_TO_LY: return "au-to-ly";
2394
2395 case GAL_ARITHMETIC_OP_MINVAL: return "minvalue";
2396 case GAL_ARITHMETIC_OP_MAXVAL: return "maxvalue";
2397 case GAL_ARITHMETIC_OP_NUMBERVAL: return "numbervalue";
2398 case GAL_ARITHMETIC_OP_SUMVAL: return "sumvalue";
2399 case GAL_ARITHMETIC_OP_MEANVAL: return "meanvalue";
2400 case GAL_ARITHMETIC_OP_STDVAL: return "stdvalue";
2401 case GAL_ARITHMETIC_OP_MEDIANVAL: return "medianvalue";
2402
2403 case GAL_ARITHMETIC_OP_MIN: return "min";
2404 case GAL_ARITHMETIC_OP_MAX: return "max";
2405 case GAL_ARITHMETIC_OP_NUMBER: return "number";
2406 case GAL_ARITHMETIC_OP_SUM: return "sum";
2407 case GAL_ARITHMETIC_OP_MEAN: return "mean";
2408 case GAL_ARITHMETIC_OP_STD: return "std";
2409 case GAL_ARITHMETIC_OP_MEDIAN: return "median";
2410 case GAL_ARITHMETIC_OP_QUANTILE: return "quantile";
2411 case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER: return "sigclip-number";
2412 case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN: return "sigclip-median";
2413 case GAL_ARITHMETIC_OP_SIGCLIP_MEAN: return "sigclip-mean";
2414 case GAL_ARITHMETIC_OP_SIGCLIP_STD: return "sigclip-number";
2415
2416 case GAL_ARITHMETIC_OP_MKNOISE_SIGMA: return "mknoise-sigma";
2417 case GAL_ARITHMETIC_OP_MKNOISE_POISSON: return "mknoise-poisson";
2418 case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM: return "mknoise-uniform";
2419
2420 case GAL_ARITHMETIC_OP_SIZE: return "size";
2421
2422 case GAL_ARITHMETIC_OP_TO_UINT8: return "uchar";
2423 case GAL_ARITHMETIC_OP_TO_INT8: return "char";
2424 case GAL_ARITHMETIC_OP_TO_UINT16: return "ushort";
2425 case GAL_ARITHMETIC_OP_TO_INT16: return "short";
2426 case GAL_ARITHMETIC_OP_TO_UINT32: return "uint";
2427 case GAL_ARITHMETIC_OP_TO_INT32: return "int";
2428 case GAL_ARITHMETIC_OP_TO_UINT64: return "ulong";
2429 case GAL_ARITHMETIC_OP_TO_INT64: return "long";
2430 case GAL_ARITHMETIC_OP_TO_FLOAT32: return "float32";
2431 case GAL_ARITHMETIC_OP_TO_FLOAT64: return "float64";
2432
2433 case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE: return "box-around-ellipse";
2434
2435 case GAL_ARITHMETIC_OP_MAKENEW: return "makenew";
2436
2437 default: return NULL;
2438 }
2439 return NULL;
2440 }
2441
2442
2443
2444
2445
2446 gal_data_t *
gal_arithmetic(int operator,size_t numthreads,int flags,...)2447 gal_arithmetic(int operator, size_t numthreads, int flags, ...)
2448 {
2449 va_list va;
2450 gal_data_t *d1, *d2, *d3, *out=NULL;
2451
2452 /* Prepare the variable arguments (starting after the flags argument). */
2453 va_start(va, flags);
2454
2455 /* Depending on the operator, do the job: */
2456 switch(operator)
2457 {
2458 /* Binary operators with any data type. */
2459 case GAL_ARITHMETIC_OP_PLUS:
2460 case GAL_ARITHMETIC_OP_MINUS:
2461 case GAL_ARITHMETIC_OP_MULTIPLY:
2462 case GAL_ARITHMETIC_OP_DIVIDE:
2463 case GAL_ARITHMETIC_OP_LT:
2464 case GAL_ARITHMETIC_OP_LE:
2465 case GAL_ARITHMETIC_OP_GT:
2466 case GAL_ARITHMETIC_OP_GE:
2467 case GAL_ARITHMETIC_OP_EQ:
2468 case GAL_ARITHMETIC_OP_NE:
2469 case GAL_ARITHMETIC_OP_AND:
2470 case GAL_ARITHMETIC_OP_OR:
2471 d1 = va_arg(va, gal_data_t *);
2472 d2 = va_arg(va, gal_data_t *);
2473 out=arithmetic_binary(operator, flags, d1, d2);
2474 break;
2475
2476 case GAL_ARITHMETIC_OP_NOT:
2477 d1 = va_arg(va, gal_data_t *);
2478 out=arithmetic_not(d1, flags);
2479 break;
2480
2481 case GAL_ARITHMETIC_OP_ISBLANK:
2482 d1 = va_arg(va, gal_data_t *);
2483 out = gal_blank_flag(d1);
2484 if(flags & GAL_ARITHMETIC_FLAG_FREE) gal_data_free(d1);
2485 break;
2486
2487 case GAL_ARITHMETIC_OP_WHERE:
2488 d1 = va_arg(va, gal_data_t *); /* To modify value/array. */
2489 d2 = va_arg(va, gal_data_t *); /* Condition (unsigned char). */
2490 d3 = va_arg(va, gal_data_t *); /* If true value/array. */
2491 arithmetic_where(flags, d1, d2, d3);
2492 out=d1;
2493 break;
2494
2495 /* Unary function operators. */
2496 case GAL_ARITHMETIC_OP_SQRT:
2497 case GAL_ARITHMETIC_OP_LOG:
2498 case GAL_ARITHMETIC_OP_LOG10:
2499 case GAL_ARITHMETIC_OP_SIN:
2500 case GAL_ARITHMETIC_OP_COS:
2501 case GAL_ARITHMETIC_OP_TAN:
2502 case GAL_ARITHMETIC_OP_ASIN:
2503 case GAL_ARITHMETIC_OP_ACOS:
2504 case GAL_ARITHMETIC_OP_ATAN:
2505 case GAL_ARITHMETIC_OP_SINH:
2506 case GAL_ARITHMETIC_OP_COSH:
2507 case GAL_ARITHMETIC_OP_TANH:
2508 case GAL_ARITHMETIC_OP_ASINH:
2509 case GAL_ARITHMETIC_OP_ACOSH:
2510 case GAL_ARITHMETIC_OP_ATANH:
2511 case GAL_ARITHMETIC_OP_AU_TO_PC:
2512 case GAL_ARITHMETIC_OP_PC_TO_AU:
2513 case GAL_ARITHMETIC_OP_LY_TO_PC:
2514 case GAL_ARITHMETIC_OP_PC_TO_LY:
2515 case GAL_ARITHMETIC_OP_LY_TO_AU:
2516 case GAL_ARITHMETIC_OP_AU_TO_LY:
2517 case GAL_ARITHMETIC_OP_RA_TO_DEGREE:
2518 case GAL_ARITHMETIC_OP_DEC_TO_DEGREE:
2519 case GAL_ARITHMETIC_OP_DEGREE_TO_RA:
2520 case GAL_ARITHMETIC_OP_DEGREE_TO_DEC:
2521 d1 = va_arg(va, gal_data_t *);
2522 out=arithmetic_function_unary(operator, flags, d1);
2523 break;
2524
2525 /* Binary function operators. */
2526 case GAL_ARITHMETIC_OP_POW:
2527 case GAL_ARITHMETIC_OP_ATAN2:
2528 case GAL_ARITHMETIC_OP_COUNTS_TO_MAG:
2529 case GAL_ARITHMETIC_OP_MAG_TO_COUNTS:
2530 case GAL_ARITHMETIC_OP_COUNTS_TO_JY:
2531 d1 = va_arg(va, gal_data_t *);
2532 d2 = va_arg(va, gal_data_t *);
2533 out=arithmetic_function_binary_flt(operator, flags, d1, d2);
2534 break;
2535
2536 /* Statistical operators that return one value. */
2537 case GAL_ARITHMETIC_OP_MINVAL:
2538 case GAL_ARITHMETIC_OP_MAXVAL:
2539 case GAL_ARITHMETIC_OP_NUMBERVAL:
2540 case GAL_ARITHMETIC_OP_SUMVAL:
2541 case GAL_ARITHMETIC_OP_MEANVAL:
2542 case GAL_ARITHMETIC_OP_STDVAL:
2543 case GAL_ARITHMETIC_OP_MEDIANVAL:
2544 d1 = va_arg(va, gal_data_t *);
2545 out=arithmetic_from_statistics(operator, flags, d1);
2546 break;
2547
2548 /* Absolute operator. */
2549 case GAL_ARITHMETIC_OP_ABS:
2550 d1 = va_arg(va, gal_data_t *);
2551 out=arithmetic_abs(flags, d1);
2552 break;
2553
2554 /* Multi-operand operators */
2555 case GAL_ARITHMETIC_OP_MIN:
2556 case GAL_ARITHMETIC_OP_MAX:
2557 case GAL_ARITHMETIC_OP_NUMBER:
2558 case GAL_ARITHMETIC_OP_SUM:
2559 case GAL_ARITHMETIC_OP_MEAN:
2560 case GAL_ARITHMETIC_OP_STD:
2561 case GAL_ARITHMETIC_OP_MEDIAN:
2562 case GAL_ARITHMETIC_OP_QUANTILE:
2563 case GAL_ARITHMETIC_OP_SIGCLIP_STD:
2564 case GAL_ARITHMETIC_OP_SIGCLIP_MEAN:
2565 case GAL_ARITHMETIC_OP_SIGCLIP_MEDIAN:
2566 case GAL_ARITHMETIC_OP_SIGCLIP_NUMBER:
2567 d1 = va_arg(va, gal_data_t *);
2568 d2 = va_arg(va, gal_data_t *);
2569 out=arithmetic_multioperand(operator, flags, d1, d2, numthreads);
2570 break;
2571
2572 /* Binary operators that only work on integer types. */
2573 case GAL_ARITHMETIC_OP_BITAND:
2574 case GAL_ARITHMETIC_OP_BITOR:
2575 case GAL_ARITHMETIC_OP_BITXOR:
2576 case GAL_ARITHMETIC_OP_BITLSH:
2577 case GAL_ARITHMETIC_OP_BITRSH:
2578 case GAL_ARITHMETIC_OP_MODULO:
2579 d1 = va_arg(va, gal_data_t *);
2580 d2 = va_arg(va, gal_data_t *);
2581 out=arithmetic_binary(operator, flags, d1, d2);
2582 break;
2583
2584 case GAL_ARITHMETIC_OP_BITNOT:
2585 d1 = va_arg(va, gal_data_t *);
2586 out=arithmetic_bitwise_not(flags, d1);
2587 break;
2588
2589 /* Adding noise. */
2590 case GAL_ARITHMETIC_OP_MKNOISE_SIGMA:
2591 case GAL_ARITHMETIC_OP_MKNOISE_POISSON:
2592 case GAL_ARITHMETIC_OP_MKNOISE_UNIFORM:
2593 d1 = va_arg(va, gal_data_t *);
2594 d2 = va_arg(va, gal_data_t *);
2595 out=arithmetic_mknoise(operator, flags, d1, d2);
2596 break;
2597
2598 /* Size operator */
2599 case GAL_ARITHMETIC_OP_SIZE:
2600 d1 = va_arg(va, gal_data_t *);
2601 d2 = va_arg(va, gal_data_t *);
2602 out=arithmetic_size(operator, flags, d1, d2);
2603 break;
2604
2605 /* Conversion operators. */
2606 case GAL_ARITHMETIC_OP_TO_UINT8:
2607 case GAL_ARITHMETIC_OP_TO_INT8:
2608 case GAL_ARITHMETIC_OP_TO_UINT16:
2609 case GAL_ARITHMETIC_OP_TO_INT16:
2610 case GAL_ARITHMETIC_OP_TO_UINT32:
2611 case GAL_ARITHMETIC_OP_TO_INT32:
2612 case GAL_ARITHMETIC_OP_TO_UINT64:
2613 case GAL_ARITHMETIC_OP_TO_INT64:
2614 case GAL_ARITHMETIC_OP_TO_FLOAT32:
2615 case GAL_ARITHMETIC_OP_TO_FLOAT64:
2616 d1 = va_arg(va, gal_data_t *);
2617 out=arithmetic_change_type(d1, operator, flags);
2618 break;
2619
2620 /* Calculate the width and height of a box surrounding an ellipse with
2621 a certain major axis, minor axis and position angle. */
2622 case GAL_ARITHMETIC_OP_BOX_AROUND_ELLIPSE:
2623 d1 = va_arg(va, gal_data_t *);
2624 d2 = va_arg(va, gal_data_t *);
2625 d3 = va_arg(va, gal_data_t *);
2626 out=arithmetic_box_around_ellipse(d1, d2, d3, operator, flags);
2627 break;
2628
2629 /* Build dataset from scratch. */
2630 case GAL_ARITHMETIC_OP_MAKENEW:
2631 d1 = va_arg(va, gal_data_t *);
2632 out=arithmetic_makenew(d1);
2633 break;
2634
2635 /* When the operator is not recognized. */
2636 default:
2637 error(EXIT_FAILURE, 0, "%s: the argument '%d' could not be "
2638 "interpretted as an operator", __func__, operator);
2639 }
2640
2641 /* End the variable argument structure and return. */
2642 va_end(va);
2643 return out;
2644 }
2645