[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9e553f0 022/125: All old arithmetic operators
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 9e553f0 022/125: All old arithmetic operators are now implemented |
|
Date: |
Sun, 23 Apr 2017 22:36:29 -0400 (EDT) |
branch: master
commit 9e553f0e05c8982b4432cd2e7c850b5b781c1a78
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
All old arithmetic operators are now implemented
The final remaining operators in Arithmetic were those that worked on any
number of operands (like `min', or `median'). In the old versions, these
operators would pop every operand. But now, they read one operand first
(which must be a number and will tell the operator how many operands to
pop). Then it will pop the given number of operands and define a linked
list of data structures which will be fed into `gal_data_arithmetic' for
the actual operation.
The Arithmetic chapter of the book was also updated to include all the
updates. With this commit, the first draft of Arithmetic using the new root
data structure has been implemented and we can go onto the rest of the
programs/librarires.
---
bin/arithmetic/arithmetic.c | 160 +++++++--
bin/arithmetic/operands.c | 46 ++-
bin/arithmetic/operators.c | 825 --------------------------------------------
bin/arithmetic/operators.h | 96 ------
doc/gnuastro.texi | 303 +++++++++-------
lib/data-arithmetic-other.c | 374 ++++++++++++++++++--
lib/data-arithmetic-other.h | 4 +
lib/data.c | 117 +++++++
lib/gnuastro/data.h | 23 ++
lib/gnuastro/qsort.h | 52 +++
lib/qsort.c | 105 ++++++
11 files changed, 997 insertions(+), 1108 deletions(-)
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 4389893..5b75361 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -54,6 +54,98 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
+/***************************************************************/
+/************* Internal functions *************/
+/***************************************************************/
+size_t
+set_number_of_operands(struct imgarithparams *p, gal_data_t *data,
+ char *token_string)
+{
+ unsigned char *uc;
+ char *c;
+ unsigned short *us;
+ short *s;
+ unsigned int *ui;
+ int *i;
+ unsigned long *ul;
+ long *l;
+ LONGLONG *L;
+
+ /* Check if its a number. */
+ if(data->size>1)
+ error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" "
+ "operator must be a number, not an array", token_string);
+
+ /* Check its type and return the value. */
+ switch(data->type)
+ {
+
+ /* For the integer types, if they are unsigned, then just pass their
+ value, if they are signed, you have to make sure they are zero or
+ positive. */
+ case GAL_DATA_TYPE_UCHAR:
+ uc=data->array; if(*uc>0) return *uc;
+ break;
+ case GAL_DATA_TYPE_CHAR:
+ c=data->array; if(*c>0) return *c;
+ break;
+ case GAL_DATA_TYPE_USHORT:
+ us=data->array; if(*us>0) return *us;
+ break;
+ case GAL_DATA_TYPE_SHORT:
+ s=data->array; if(*s>0) return *s;
+ break;
+ case GAL_DATA_TYPE_UINT:
+ ui=data->array; if(*ui>0) return *ui;
+ break;
+ case GAL_DATA_TYPE_INT:
+ i=data->array; if(*i>0) return *i;
+ break;
+ case GAL_DATA_TYPE_ULONG:
+ ul=data->array; if(*ul>0) return *ul;
+ break;
+ case GAL_DATA_TYPE_LONG:
+ l=data->array; if(*l>0) return *l;
+ break;
+ case GAL_DATA_TYPE_LONGLONG:
+ L=data->array; if(*L>0) return *L;
+ break;
+
+ /* Floating point numbers are not acceptable in this context. */
+ case GAL_DATA_TYPE_FLOAT:
+ case GAL_DATA_TYPE_DOUBLE:
+ error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" "
+ "operator must be an integer type", token_string);
+
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`set_number_of_operands", data->type);
+ }
+
+ /* If control reaches here, then the number must have been a negative
+ value, so print an error. */
+ error(EXIT_FAILURE, 0, "the first popped operand to the \"%s\" operator "
+ "cannot be zero or a negative number", token_string);
+ return 0;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/***************************************************************/
@@ -68,6 +160,7 @@ void
reversepolish(struct imgarithparams *p)
{
int op=0, nop=0;
+ unsigned int numop, i;
struct gal_linkedlist_stll *token;
gal_data_t *d1=NULL, *d2=NULL, *d3=NULL;
unsigned char flags = ( GAL_DATA_ARITH_INPLACE | GAL_DATA_ARITH_FREE
@@ -117,7 +210,7 @@ reversepolish(struct imgarithparams *p)
else if (!strcmp(token->v, "log10"))
{ op=GAL_DATA_OPERATOR_LOG10; nop=1; }
- /* Statistical operators. */
+ /* Statistical/higher-level operators. */
else if (!strcmp(token->v, "minvalue"))
{ op=GAL_DATA_OPERATOR_MINVAL; nop=1; }
else if (!strcmp(token->v, "maxvalue"))
@@ -126,6 +219,8 @@ reversepolish(struct imgarithparams *p)
{ op=GAL_DATA_OPERATOR_MIN; nop=-1; }
else if (!strcmp(token->v, "max"))
{ op=GAL_DATA_OPERATOR_MAX; nop=-1; }
+ else if (!strcmp(token->v, "sum"))
+ { op=GAL_DATA_OPERATOR_SUM; nop=-1; }
else if (!strcmp(token->v, "average"))
{ op=GAL_DATA_OPERATOR_AVERAGE; nop=-1; }
else if (!strcmp(token->v, "median"))
@@ -197,35 +292,50 @@ reversepolish(struct imgarithparams *p)
token->v);
-
/* Pop the necessary number of operators. Note that the operators
are poped from a linked list (which is last-in-first-out). So
for the operators which need a specific order, the first poped
operand is actally the last (right most, in in-fix notation)
input operand.*/
- if(nop==1)
- d1=pop_operand(p, token->v);
- else if(nop==2)
- {
- d2=pop_operand(p, token->v);
- d1=pop_operand(p, token->v);
- }
- else if(nop==3)
- {
- d3=pop_operand(p, token->v);
- d2=pop_operand(p, token->v);
- d1=pop_operand(p, token->v);
- }
- else if(nop==-1)
- /* It might be possible to add a `next' pointer to the data
- structure so it can act like a linked list too. In that
- case, t */
- error(EXIT_FAILURE, 0, "Operators with an unknown number of "
- "operands are not yet implemented.");
- else
- error(EXIT_FAILURE, 0, "No operators need %d operands", nop);
-
- /* Do the arithemtic operation. */
+ switch(nop)
+ {
+ case 1:
+ d1=pop_operand(p, token->v);
+ break;
+
+ case 2:
+ d2=pop_operand(p, token->v);
+ d1=pop_operand(p, token->v);
+ break;
+
+ case 3:
+ d3=pop_operand(p, token->v);
+ d2=pop_operand(p, token->v);
+ d1=pop_operand(p, token->v);
+ break;
+
+ case -1:
+ /* This case is when the number of operands is itself an
+ operand. So the first popped operand must be an integer
+ number, we will use that to construct a linked list of any
+ number of operands within the single `d1' pointer. */
+ d1=NULL;
+ numop=set_number_of_operands(p, pop_operand(p, token->v),
+ token->v);
+ for(i=0;i<numop;++i)
+ gal_data_add_to_ll(&d1, pop_operand(p, token->v));
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "No operators need %d operands", nop);
+ }
+
+
+ /* Run the arithmetic operation. Note that `gal_data_arithmetic'
+ is a variable argument function (like printf). So the number
+ of arguments it uses depend on the operator. So when the
+ operator doesn't need three operands, the extra arguments will
+ be ignored. */
add_operand(p, NULL, gal_data_arithmetic(op, flags, d1, d2, d3));
}
}
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 1153941..d18754e 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -60,28 +60,36 @@ add_operand(struct imgarithparams *p, char *filename,
gal_data_t *data)
{
struct operand *newnode;
- /* Allocate space for the new operand. */
- errno=0;
- newnode=malloc(sizeof *newnode);
- if(newnode==NULL)
- error(EXIT_FAILURE, errno, "imgarith.c: Making new element in operand");
-
- /* Fill in the values. */
- newnode->data=data;
- newnode->filename=filename;
-
- if(filename != NULL && gal_fits_name_is_fits(filename))
+ /* Some operators might not actually return any dataset (for example the
+ multi-operand datasets when the number of operands is 0, which can
+ happen in a script), in such cases, we shouldn't add anything to the
+ stack. */
+ if(data)
{
- /* Set the HDU for this filename. */
- gal_linkedlist_pop_from_stll(&p->hdus, &newnode->hdu);
+ /* Allocate space for the new operand. */
+ errno=0;
+ newnode=malloc(sizeof *newnode);
+ if(newnode==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for newnode in"
+ "add_operand", sizeof *newnode);
+
+ /* Fill in the values. */
+ newnode->data=data;
+ newnode->filename=filename;
+
+ if(filename != NULL && gal_fits_name_is_fits(filename))
+ {
+ /* Set the HDU for this filename. */
+ gal_linkedlist_pop_from_stll(&p->hdus, &newnode->hdu);
- /* Increment the FITS counter. */
- ++p->addcounter;
- }
+ /* Increment the FITS counter. */
+ ++p->addcounter;
+ }
- /* Make the link to the previous list. */
- newnode->next=p->operands;
- p->operands=newnode;
+ /* Make the link to the previous list. */
+ newnode->next=p->operands;
+ p->operands=newnode;
+ }
}
diff --git a/bin/arithmetic/operators.c b/bin/arithmetic/operators.c
deleted file mode 100644
index 04349d0..0000000
--- a/bin/arithmetic/operators.c
+++ /dev/null
@@ -1,825 +0,0 @@
-/*********************************************************************
-Arithmetic - Do arithmetic operations on images.
-Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#include <config.h>
-
-#include <math.h>
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <string.h>
-#include <stdlib.h>
-
-#include <gnuastro/data.h>
-#include <gnuastro/array.h>
-#include <gnuastro/statistics.h>
-
-#include "main.h"
-
-#include "operands.h"
-
-
-
-
-
-
-
-
-
-
-void
-sum(struct imgarithparams *p)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- char *operator="+";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in the first
- input. */
- gal_array_dsum_arrays(f->array, s->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
-
- /* Clean up. */
- gal_data_free(s);
- }
- else if(f) /* Only the first is an array. */
- {
- gal_array_dsum_const(f->array, f->size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else if(s) /* Only the second is an array. */
- {
- gal_array_dsum_const(s->array, s->size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, fnum+snum, NOOPTDATA);
-}
-
-
-
-
-
-void
-subtract(struct imgarithparams *p)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- char *operator="-";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
-
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_dsubtract_arrays(s->array, f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s->array);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else if(f) /* Only the first is an array. */
- {
- gal_array_dconst_subtract(f->array, f->size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f->array);
- }
- else if(s) /* Only the second is an array. */
- {
- gal_array_dsubtract_const(s->array, s->size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, snum-fnum, NOOPTDATA);
-}
-
-
-
-
-
-void
-multiply(struct imgarithparams *p)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- char *operator="*";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in farr. */
- gal_array_dmultip_arrays(f->array, s->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
-
- /* Clean up. */
- gal_data_free(s);
- }
- else if(f) /* Only the first is an array. */
- {
- gal_array_dmultip_const(f->array, f->size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else if(s) /* Only the first is an array. */
- {
- gal_array_dmultip_const(s->array, s->size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, fnum*snum, NOOPTDATA);
-}
-
-
-
-
-
-void
-divide(struct imgarithparams *p)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- char *operator="/";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_ddivide_arrays(s->array, f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else if(f) /* Only the first is an array. */
- {
- gal_array_dconst_divide(f->array, f->size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else if(s) /* Only the first is an array. */
- {
- gal_array_ddivide_const(s->array, s->size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, snum/fnum, NOOPTDATA);
-}
-
-
-
-
-
-void
-topower(struct imgarithparams *p, char *op)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- char *operator = op?op:"pow";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
-
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_dpower_arrays(s->array, f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else if(f) /* Only the first is an array. */
- {
- gal_array_dconst_power(f->array, f->size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else if(s) /* Only the first is an array. */
- {
- gal_array_dpower_const(s->array, s->size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, pow(snum, fnum), NOOPTDATA);
-}
-
-
-
-
-
-void
-alloppixs(struct imgarithparams *p, char *operator)
-{
- int i, j; /* Integer to allow negative checks. */
- double num; /* Temporary number holder. */
- int firstisdata=-1; /* ==-1: unset. ==1: array. ==0: number. */
- double *allpixels=NULL; /* Array for all values in one pixel. */
- double **allarrays=NULL; /* Pointer to all the arrays. */
- gal_data_t **alldata=NULL; /* Array for pointers to input datasets. */
- size_t numop=num_operands(p);
- double (*thisfunction)(double *, size_t)=NULL;
-
-
- /* First set the appropriate function to call. */
- if(!strcmp(operator, "min"))
- thisfunction = &gal_statistics_double_min_return;
- else if(!strcmp(operator, "max"))
- thisfunction = &gal_statistics_double_max_return;
- else if(!strcmp(operator, "median"))
- thisfunction = &gal_statistics_median_double_in_place;
- else if(!strcmp(operator, "average"))
- thisfunction = &gal_statistics_double_average;
- else
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we "
- "can address the problem. The value of `operator' in "
- "alloppixs (%s) is not recognized",
- PACKAGE_BUGREPORT, operator);
-
-
- /* Allocate the array of pointers to all input arrays and also the
- array to temporarily keep all values for each pixel */
- errno=0;
- alldata=malloc(numop*sizeof *alldata);
- if(alldata==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for alldata in alloppixs",
- numop*sizeof *alldata);
- errno=0;
- allarrays=malloc(numop*sizeof *allarrays);
- if(allarrays==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for allarrays in alloppixs",
- numop*sizeof *allarrays);
- errno=0;
- allpixels=malloc(numop*sizeof *allpixels);
- if(allpixels==NULL)
- error(EXIT_FAILURE, errno, "%zu bytes for allpixels in alloppixs",
- numop*sizeof *allpixels);
-
-
- /* Prepare all the inputs, note that since it is a linked list, the
- operands pop from the last to first. Here order is not important,
- but in other cases that it might be (and also for debugging),
- here the inputs are put in order so the indexs start from the
- last to first. */
- for(i=numop-1;i>=0;--i)
- {
- /* Pop out the operand. */
- pop_operand(p, &num, &alldata[i], operator);
-
- /* Do the appropriate action if it is an array or a number. */
- if(alldata[i])
- switch(firstisdata)
- {
- case -1:
- firstisdata=1;
- allarrays[i]=alldata[i]->array;
- break;
- case 1:
- allarrays[i]=alldata[i]->array;
- break;
- case 0:
- error(EXIT_FAILURE, 0, "for the %s operator, all operands "
- "must be either an array or number", operator);
- break;
- default:
- error(EXIT_FAILURE, 0, "a Bug! Please contact us at %s so we "
- "can address the problem. The value of firstisdata (%d) "
- "in the alloppixs function is not recognized",
- PACKAGE_BUGREPORT, firstisdata);
- }
- else
- switch(firstisdata)
- {
- case -1:
- firstisdata=0;
- allpixels[i]=num;
- break;
- case 0:
- allpixels[i]=num;
- break;
- case 1:
- error(EXIT_FAILURE, 0, "for the %s operator, all operands "
- "must be either an array or number", operator);
- break;
- default:
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we "
- "can address the problem. The value of firstisdata (%d) "
- "in the alloppixs function is not recognized",
- PACKAGE_BUGREPORT, firstisdata);
- }
- }
-
-
- /* Do the operation and report the result: */
- if(firstisdata==1)
- {
- /* For each pixel, find the value and replace it with the first
- operand. */
- for(i=0;i<alldata[0]->size;++i)
- {
- /* Go over all the inputs and put the appropriate pixel
- values in the allpixels array. */
- for(j=0;j<numop;++j)
- allpixels[j]=allarrays[j][i];
-
- /* Do the appropriate action. */
- allarrays[0][i]=(*thisfunction)(allpixels, numop);
- }
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, alldata[0]);
-
- /* Free all the extra operands */
- for(i=1;i<numop;++i) gal_data_free(alldata[i]);
- }
- else
- add_operand(p, NOOPTFILENAME, (*thisfunction)(allpixels, numop),
- NOOPTDATA);
-
-
- /* Clean up. Reminder: the allocated arrays were freed above as part of
- `gal_data_free'. */
- free(alldata);
- free(allarrays);
- free(allpixels);
-}
-
-
-
-
-
-void
-takesqrt(struct imgarithparams *p)
-{
- char *operator="sqrt";
-
- /* Add a 0.5 number to the operand stack */
- add_operand(p, NOOPTFILENAME, 0.5f, NOOPTDATA);
-
- /* Call the power operator. */
- topower(p, operator);
-}
-
-
-
-
-
-void
-takelog(struct imgarithparams *p)
-{
- double fnum;
- gal_data_t *f;
- char *operator="log";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_dlog_array(f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, log(fnum), NOOPTDATA);
-}
-
-
-
-
-
-void
-takelog10(struct imgarithparams *p)
-{
- double fnum;
- gal_data_t *f;
- char *operator="log10";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_dlog10_array(f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, log10(fnum), NOOPTDATA);
-}
-
-
-
-
-
-void
-takeabs(struct imgarithparams *p)
-{
- double fnum;
- gal_data_t *f;
- char *operator="abs";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_array_dabs_array(f->array, f->size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fabs(fnum), NOOPTDATA);
-}
-
-
-
-
-
-void
-findmin(struct imgarithparams *p)
-{
- gal_data_t *f;
- double min, fnum;
- char *operator="min";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_statistics_double_min(f->array, f->size, &min);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, min, NOOPTDATA);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fnum, NOOPTDATA);
-}
-
-
-
-
-
-void
-findmax(struct imgarithparams *p)
-{
- gal_data_t *f;
- double max, fnum;
- char *operator="max";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- gal_statistics_double_max(f->array, f->size, &max);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, max, NOOPTDATA);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fnum, NOOPTDATA);
-}
-
-
-
-
-
-int
-lessthan(double left, double right)
-{ return left<right; }
-
-int
-lessequal(double left, double right)
-{ return left<=right; }
-
-int
-greaterthan(double left, double right)
-{ return left>right; }
-
-int
-greaterequal(double left, double right)
-{ return left>=right; }
-
-int
-equal(double left, double right)
-{ return left==right; }
-
-int
-notequal(double left, double right)
-{ return left!=right; }
-
-
-
-
-
-void
-conditionals(struct imgarithparams *p, char *operator)
-{
- double fnum, snum; /* First or second number. */
- gal_data_t *f, *s; /* First or second dataset. */
- double *fp, *sp, *ff, *ss;
- int (*thisfunction)(double, double)=NULL;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Set the function to use. */
- if(!strcmp(operator, "lt")) thisfunction = &lessthan;
- else if(!strcmp(operator, "le")) thisfunction = &lessequal;
- else if(!strcmp(operator, "gt")) thisfunction = &greaterthan;
- else if(!strcmp(operator, "ge")) thisfunction = &greaterequal;
- else if(!strcmp(operator, "eq")) thisfunction = &equal;
- else if(!strcmp(operator, "neq")) thisfunction = ¬equal;
- else
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we "
- "can address the problem. The value of `operator' in "
- "conditionals (%s) is not recognized",
- PACKAGE_BUGREPORT, operator);
-
- /* Do the operation: */
- if(f && s) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- fp=f->array;
- ss=(sp=s->array)+s->size;
- do *sp = thisfunction(*sp, *fp++); while(++sp<ss);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
-
- /* Clean up. */
- gal_data_free(f);
- }
- else if(f) /* Only the first is an array. */
- {
- ff=(fp=f->array)+f->size;
- do *fp = thisfunction(snum, *fp); while(++fp<ff);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else if(s) /* Only the first is an array. */
- {
- ss=(sp=s->array)+s->size;
- do *sp = thisfunction(*sp, fnum); while(++sp<ss);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, s);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, thisfunction(snum, fnum), NOOPTDATA);
-}
-
-
-
-
-
-void
-andor(struct imgarithparams *p, char *operator)
-{
- double fnum, snum;
- gal_data_t *f, *s;
- double *fp, *sp, *ff;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
- pop_operand(p, &snum, &s, operator);
-
- /* Do a small sanity check: */
- if(strcmp(operator, "and") && strcmp(operator, "or") )
- error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we "
- "can address the problem. The value of `operator' in "
- "`andor' (%s) is not recognized", PACKAGE_BUGREPORT, operator);
-
- /* Do the operation: */
- if(f && s)
- {
- /* Fill the first array with the result. IMPORTANT: It is important
- that the second array pointer is the first checked array, since it
- is also incremented. In the `or' operation, if `*f' is successful,
- `*s++' is never called and so not incremented. */
- sp = s->array;
- ff = (fp=f->array) + f->size;
- if(!strcmp(operator, "and"))
- do *fp = *sp++ && *fp; while(++fp<ff);
- else
- do *fp = *sp++ || *fp; while(++fp<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
-
- /* Clean up */
- gal_data_free(s);
- }
- else if(f==NULL || s==NULL)
- error(EXIT_FAILURE, 0, "The `and' and `or' operators need two operators "
- "of the same type: either both images or both numbers.");
- else
- add_operand(p, NOOPTFILENAME,
- !strcmp(operator, "and") ? fnum && snum : fnum || snum,
- NOOPTDATA);
-}
-
-
-
-
-void
-notfunc(struct imgarithparams *p)
-{
- double fnum;
- gal_data_t *f;
- double *fp, *ff;
- char *operator="not";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f)
- {
- /* Fill the array with the output values. */
- ff = (fp=f->array) + f->size;
- do *fp = !(*fp); while(++fp<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else
- add_operand(p, NOOPTFILENAME, !fnum, NOOPTDATA);
-}
-
-
-
-
-
-/* In order to not conflict with the internal C `is...' functions, and in
- particular the `isblank' function, we are calling this function
- opisblank for operator-isblank. */
-void
-opisblank(struct imgarithparams *p)
-{
- gal_data_t *f;
- char *operator="isblank";
- double *fp, *ff, fnum;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &f, operator);
-
- /* Do the operation: */
- if(f) /* Operand is array. */
- {
- /* Do the operation, note that the output is stored in the first
- input. Also note that since the linked list is
- first-in-first-out, the second operand should be put first
- here. */
- ff=(fp=f->array)+f->size;
- do *fp = isnan(*fp); while(++fp<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, f);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, isnan(fnum), NOOPTDATA);
-}
-
-
-
-
-
-/* Replace the pixels in the second popped element with the first. While
- choosing the pixels that are selected from the third The third popped
- element. The third is considered to be an array that can only be filled
- with 0 or 1. */
-void
-where(struct imgarithparams *p)
-{
- gal_data_t *n, *c, *i; /* New, condition, input datasets. */
- double nnum, cnum, inum; /* First, second, or third number. */
- double *np, *cp, *ip, *ii;
-
- /* ORDER IS VERY IMPORTANT HERE. Pop out the number of operands needed. */
- pop_operand(p, &nnum, &n, "where"); /* New value. */
- pop_operand(p, &cnum, &c, "where"); /* Condition. */
- pop_operand(p, &inum, &i, "where"); /* Input. */
-
- /* Do the operation: */
- if(i && c) /* Both are arrays. */
- {
- cp=c->array;
- ii=(ip=i->array)+i->size;
- if(n) /* `new' is an array, not number. */
- {
- np=n->array;
- /* Note that we need to increment "fp" after the check and
- replace. If the increment is inside the conditional replace,
- then when t==0, the increment won't work.*/
- do {*ip = *cp++ ? *np : *ip; ++np;} while(++ip<ii);
- }
- else /* `new' is a number, not array. */
- do *ip = *cp++ ? nnum : *ip; while(++ip<ii);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, i);
-
- /* Clean up. */
- gal_data_free(c);
- gal_data_free(n);
- }
- else if ( i==NULL && c==NULL && n==NULL )
- add_operand(p, NOOPTFILENAME, cnum ? nnum : inum, NOOPTDATA);
- else
- error(EXIT_FAILURE, 0, "the first and second arguments (second and "
- "third popped elements) to `where' have to be arrays, or all have "
- "to be numbers.");
-}
diff --git a/bin/arithmetic/operators.h b/bin/arithmetic/operators.h
deleted file mode 100644
index 6ce946c..0000000
--- a/bin/arithmetic/operators.h
+++ /dev/null
@@ -1,96 +0,0 @@
-/*********************************************************************
-Arithmetic - Do arithmetic operations on images.
-Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
-
-Original author:
- Mohammad Akhlaghi <address@hidden>
-Contributing author(s):
-Copyright (C) 2015, Free Software Foundation, Inc.
-
-Gnuastro is free software: you can redistribute it and/or modify it
-under the terms of the GNU General Public License as published by the
-Free Software Foundation, either version 3 of the License, or (at your
-option) any later version.
-
-Gnuastro is distributed in the hope that it will be useful, but
-WITHOUT ANY WARRANTY; without even the implied warranty of
-MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
-General Public License for more details.
-
-You should have received a copy of the GNU General Public License
-along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
-**********************************************************************/
-#ifndef OPERATORS_H
-#define OPERATORS_H
-
-
-void
-sum(struct imgarithparams *p);
-
-void
-subtract(struct imgarithparams *p);
-
-void
-multiply(struct imgarithparams *p);
-
-void
-divide(struct imgarithparams *p);
-
-void
-topower(struct imgarithparams *p, char *op);
-
-void
-alloppixs(struct imgarithparams *p, char *operator);
-
-void
-takesqrt(struct imgarithparams *p);
-
-void
-takelog(struct imgarithparams *p);
-
-void
-takelog10(struct imgarithparams *p);
-
-void
-takeabs(struct imgarithparams *p);
-
-void
-findmin(struct imgarithparams *p);
-
-void
-findmax(struct imgarithparams *p);
-
-int
-lessthan(double left, double right);
-
-int
-lessequal(double left, double right);
-
-int
-greaterthan(double left, double right);
-
-int
-greaterequal(double left, double right);
-
-int
-equal(double left, double right);
-
-int
-notequal(double left, double right);
-
-void
-conditionals(struct imgarithparams *p, char *operator);
-
-void
-andor(struct imgarithparams *p, char *operator);
-
-void
-notfunc(struct imgarithparams *p);
-
-void
-opisblank(struct imgarithparams *p);
-
-void
-where(struct imgarithparams *p);
-
-#endif
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 785f822..4e97e3c 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -7044,16 +7044,18 @@ verbose mode on the command-line for each row.
@node Arithmetic, Convolve, ImageCrop, Image manipulation
@section Arithmetic
-It is commonly necessary to do arithmetic operations on the
-astronomical data. For example in the reduction of raw data it is
-necessary to subtract the Sky value (@ref{Sky value}) from each image
-image. Later (once the images as warped into a single grid using
-ImageWarp for example, see @ref{ImageWarp}), the images can be co-added
-or the output pixel grid is the average of the pixels of the
-individual input images. Arithmetic currently uses the reverse
-polish or postfix notation, see @ref{Reverse polish notation}, for
-more information on how to run Arithmetic, please see
address@hidden astarithmetic}.
+It is commonly necessary to do operations on some or all of the elements of
+a dataset independently (pixels in an image). For example, in the reduction
+of raw data it is necessary to subtract the Sky value (@ref{Sky value})
+from each image image. Later (once the images as warped into a single grid
+using ImageWarp for example, see @ref{ImageWarp}), the images are co-added
+(the output pixel grid is the average of the pixels of the individual input
+images). Arithmetic is Gnuastro's program for such operations on your
+datasets directly from the command-line. It currently uses the reverse
+polish or postfix notation, see @ref{Reverse polish notation} and will work
+on the native data types of the input images/data to reduce CPU and RAM
+resources, see @ref{Data types}. For more information on how to run
+Arithmetic, please see @ref{Invoking astarithmetic}.
@menu
@@ -7068,19 +7070,19 @@ more information on how to run Arithmetic, please see
@cindex Postfix notation
@cindex Reverse Polish Notation
-The most common notation for arithmetic operations is the infix
address@hidden@url{https://en.wikipedia.org/wiki/Infix_notation}}
-where the operator goes between the two operands, for example
address@hidden While the infix notation is the preferred way in most
-programming languages, currently Arithmetic does not use it since it
-will require parenthesis which can complicate the implementation of
-the code. In the near future we do plan to adopt this
address@hidden@url{https://savannah.gnu.org/task/index.php?13867}},
-but for the time being (due to time constraints on the developers),
-Arithmetic uses the postfix or reverse polish
address@hidden@url{https://en.wikipedia.org/wiki/Reverse_Polish_notation}}. The
-referenced Wikipedia article provides some excellent explanation on
-this notation but here we will give a short summary for
+The most common notation for arithmetic operations is the
address@hidden://en.wikipedia.org/wiki/Infix_notation, infix notation} where
+the operator goes between the two operands, for example @mymath{4+5}. While
+the infix notation is the preferred way in most programming languages,
+currently Arithmetic does not use it since it will require parenthesis
+which can complicate the implementation of the code. In the near future we
+do plan to adopt this
address@hidden@url{https://savannah.gnu.org/task/index.php?13867}}, but
+for the time being (due to time constraints on the developers), Arithmetic
+uses the postfix or
address@hidden://en.wikipedia.org/wiki/Reverse_Polish_notation, reverse polish
+notation}. The Wikipedia article provides some excellent explanation on
+this notation but here we will give a short summary here for
self-sufficiency.
In the postfix notation, the operator is placed after the operands, as we
@@ -7100,10 +7102,9 @@ operations that are done are:
@item
@command{6} is an operand, so it is pushed to the top of the stack.
@item
address@hidden is a binary operator, so pull the top two elements of the
-stack and perform addition on them (the order is @command{5+6} in the
-example above). The result is @command{11}, push it on top of the
-stack.
address@hidden is a binary operator, so pull the top two elements of the stack
+and perform addition on them (the order is @mymath{5+6} in the example
+above). The result is @command{11}, push it on top of the stack.
@item
@command{2} is an operand so push it onto the top of the stack.
@item
@@ -7112,12 +7113,13 @@ the stack (top-most is @command{2}, then @command{11})
and divide the
second one by the first.
@end enumerate
-In Gnuastro's Arithmetic, the operands can be FITS images or numbers. As
+In the Arithmetic program, the operands can be FITS images or numbers. As
you can see, very complicated procedures can be created without the need
-for parenthesis. Even functions which take an arbitrary number of arguments
-can be defined in this notation. For example the Postscript
address@hidden the EPS and PDF part of @ref{Recognized file types}
-for a little more on the Postscript language.} (the programming language in
+for parenthesis or worrying about precedence. Even functions which take an
+arbitrary number of arguments can be defined in this notation. This is a
+very powerfull notation and is used in languages like Postscript
address@hidden the EPS and PDF part of @ref{Recognized file types} for a
+little more on the Postscript language.} (the programming language in
Postscript and compiled into PDF files) uses this notation.
@@ -7132,11 +7134,11 @@ Postscript and compiled into PDF files) uses this
notation.
At the lowest level, the computer stores everything in terms of @code{1} or
@code{0}. For example, each program in Gnuastro, or each astronomical image
you take with the telescope is actually a string of millions of these zeros
-and ones. The space required to keep either of these two values is the
-smallest unit of storage, and is known as a @emph{bit}. Understanding and
-manipulating this string of bits is extremely hard for most people,
-therefore, we define packages with a certain number of these bits along
-with a standard on how to interpret these bits as @emph{type}s.
+and ones. The space required to keep a zero or one is the smallest unit of
+storage, and is known as a @emph{bit}. Understanding and manipulating this
+string of bits is extremely hard for most people, therefore, we define
+packages of these bits along with a standard on how to interpret the bits
+in each package as @emph{type}s.
The most basic standard for reading the bits is integer numbers. Since we
use CFITSIO, we define the @code{char}, @code{short}, or @code{long}
@@ -7154,7 +7156,7 @@ doesn't involve negative numbers (for example labeled
images, where pixels
can only have zero or positive integer values), it is best to use the
@code{unsigned} types.
-The second standard of converting a given number of bits to numbers is the
+Another standard of converting a given number of bits to numbers is the
floating point standard, this standard can approximately store any real
number with a given precision. There are two floating point types in
CFITSIO: @code{float} and @code{double}, for single and double precision
@@ -7191,7 +7193,8 @@ ordered on the command-line. The operands to all
operators can be a data
array (for example a FITS image) or a number, the output will be an array
or number according to the inputs. For example a number multiplied by an
array will produce an array. The conditional operators will return pixel,
-or numerical values of 0 (false) or 1 (true).
+or numerical values of 0 (false) or 1 (true) and stored in an
address@hidden char} data type (see @ref{Data types}).
@table @command
@@ -7228,61 +7231,85 @@ Absolute value of first operand, so address@hidden
abs}'' is
equivalent to @mymath{|4|}.
@item pow
-First operand to the power of the second, so address@hidden 5
-pow}'' is equivalent to @mymath{4^5}.
+First operand to the power of the second, so address@hidden 5f pow}'' is
+equivalent to @mymath{4.3^{5}}. Currently @code{pow} will only work on
+single or double precision floating point numbers or images. To be sure
+that a number is read as a floating point (even if it doesn't have any
+non-zero decimals) put an @code{f} after it.
@item sqrt
-The square root of the first operand, so address@hidden sqrt}'' is
-equivalent to @mymath{\sqrt{4}}.
+The square root of the first operand, so address@hidden sqrt}'' is equivalent
+to @mymath{\sqrt{5}}. The output type is determined from the input, so the
+output of this example will be @command{2} (since @command{5} doesn't have
+any non-zero decimal digits). If you want @command{2.23607}, run
address@hidden sqrt} instead, the @command{f} will ensure that a number will
+be read as a floating point number, even if it doesn't have decimal
+digits. If the input image has an integer type, you should explicitly
+convert the image to floating point, for example @command{a.fits float
+sqrt}, see the type conversion operators below.
@item log
-Natural logarithm of first operand, so address@hidden log}'' is
-equivalent to @mymath{\ln(4)}.
+Natural logarithm of first operand, so address@hidden log}'' is equivalent to
address@hidden(4)}. The output type is determined from the input, see the
+explanation under @command{sqrt} for more.
@item log10
-Base-10 logarithm of first operand, so address@hidden log10}'' is
-equivalent to @mymath{\log(4)}.
+Base-10 logarithm of first operand, so address@hidden log10}'' is equivalent
+to @mymath{\log(4)}. The output type is determined from the input, see the
+explanation under @command{sqrt} for more.
@item minvalue
Minimum value in the top operand on the stack, so address@hidden
-min}'' will push the the minimum pixel value in this image onto the
-stack. Therefore this operator is mainly intended for data (for
-example images), if the top operand is a number, this operator just
-returns it without any change. So note that when this operator acts on
-a single image, the output will no longer be an image, but a number.
+minvalue}'' will push the the minimum pixel value in this image onto the
+stack. Therefore this operator is mainly intended for data (for example
+images), if the top operand is a number, this operator just returns it
+without any change. So note that when this operator acts on a single image,
+the output will no longer be an image, but a number.
@item maxvalue
Maximum value of first operand, similar to @command{minvalue}.
@cindex NaN
@item min
-If there are several operands (images), each pixel of the output of
-this operator will be set to the minimum value of all the operands
-(images) in that pixel. If the operands are all numbers, then their
-minimum value will be put in the stack. So address@hidden b.fits
-c.fits min}'' will produce an image with the same size as the inputs
-but each output pixel is set to the minimum respective pixel value of
-the three input images. Important notes:
address@hidden
+The first popped operand to this operator must be a positive integer number
+which specifies how many further operands should be popped from the
+stack. The given number of operands must have the same type and size. Each
+pixel of the output of this operator will be set to the minimum value of
+the given number of operands (images) in that pixel.
address@hidden
-All the operands must be either images or numbers, they cannot be
-mixed.
+For example the following command will produce an image with the same size
+and type as the inputs but each output pixel is set to the minimum
+respective pixel value of the three input images.
+
address@hidden
+$ astarithmetic a.fits b.fits c.fits 3 min
address@hidden example
+
+Important notes:
address@hidden
@item
NaN/blank pixels will be ignored, see @ref{Blank
pixels}.
@item
-All available operands on the stack will be pulled for this operator
-to generate its output. So after calling this operator, there will
-only be one operand on the stack, see @ref{Reverse polish notation}.
+The output will have the same type as the inputs. This is natural for the
address@hidden and @command{max} operators, but for other similar operators
+(for example @command{sum}, or @command{average}) the per-pixel operations
+will be done in double precision floating point and then stored back in the
+input type. Therefore, if the input was an integer, C's internal type
+conversion will be used.
+
@end itemize
@item max
Similar to @command{min}, but the pixels of the output will contain
the maximum of the respective pixels in all operands in the stack.
address@hidden sum
+Similar to @command{min}, but the pixels of the output will contain the sum
+of the respective pixels in all input operands.
+
@item average
Similar to @command{min}, but the pixels of the output will contain
the average of the respective pixels in all operands in the stack.
@@ -7299,7 +7326,8 @@ return a value of 0. If both operands are images, then
all the pixels will
be compared with their counterparts in the other image. If only one operand
is an image, then all the pixels will be compared with the the single value
(number) of the other operand. Finally if both are numbers, then the output
-is also just one number (0 or 1).
+is also just one number (0 or 1). When the output is not a single number,
+it will be stored as an @code{unsigned char} type.
@item le
Less or equal: similar to @code{lt} (`less than' operator), but returning 1
@@ -7349,32 +7377,38 @@ blank pixels. See the ``Blank pixels'' box below for
more on Blank pixels
in Arithmetic.
@item where
-Change the input (pixel) value `where' a certain condition holds. The
-conditional operators above can be used to define the condition. Three
-operands are required for @command{where}. The input format is demonstrated
-in this simplified example:
+Change the input (pixel) value @emph{where}/if a certain condition
+holds. The conditional operators above can be used to define the
+condition. Three operands are required for @command{where}. The input
+format is demonstrated in this simplified example:
@example
-$ astarithmetic in.fits condition.fits new.fits where
+$ astarithmetic modify.fits cond.fits if-true.fits where
@end example
-The third and second popped operands (@file{in.fits} and
address@hidden above, see @ref{Reverse polish notation}) have to be
+The third and second popped operands (@file{modify.fits} and
address@hidden above, see @ref{Reverse polish notation}) have to be
images or all three operands have to be numbers. In the former case (when
-the first two are images), the third operand can be an image or a
-number. The value of any pixel in @file{in.fits} that corresponds to a
-non-zero pixel of @file{condition.fits} will be changed to the value of the
-same pixel in @file{new.fits} (or a fixed number if the first popped
-operand is a number). Commonly you won't be dealing with an actual FITS
-file of a conditional image. You will probably define the condition in the
-same run based on some other reference image and use the conditional and
-logical operators above to make a true/false (or one/zero) image for you
-internally. For example the case below:
+the first two are images), the third popped operand can be an image or a
+number. The value of any pixel in @file{out.fits} that corresponds to a
+non-zero pixel of @file{cond.fits} will be changed to the value of the same
+pixel in @file{if-true.fits} (or a fixed number). Note than the
address@hidden has to have @code{unsigned char} type.
+
+Commonly you won't be dealing with an actual FITS file of a conditional
+image. You will probably define the condition in the same run based on some
+other reference image and use the conditional and logical operators above
+to make a true/false (or one/zero) image for you internally. For example
+the case below:
@example
$ astarithmetic in.fits reference.fits 100 gt new.fits where
@end example
+In the example above, any of the @file{in.fits} pixels that has a value in
address@hidden greater than @command{100}, will be replaced with the
+corresponding pixel in @file{new.fits}.
+
@item bitand
Bitwise AND operator: only bits with values of 1 in both popped operands
will get the value of 1, the rest will be set to 0. For example (assuming
@@ -7463,13 +7497,13 @@ data types.
@cartouche
@noindent
address@hidden pixels in Arithmetic:} Currently all Arithmetic operations
-are internally done in double precision floating points. So all blank
-pixels in the image (see @ref{Blank pixels}) will be stored as IEEE NaN
-values. One particular aspect of NaN values is that by definition they will
-fail on @emph{any} comparison. Hence both equal and not-equal operators
-will fail when both their operands are NaN! The only way to select blank
-pixels is through the @command{isblank} operator explained above.
address@hidden pixels in Arithmetic:} Blank pixels in the image (see
address@hidden pixels}) will be stored based on the data type. When the input
+is floating point type, blank values are NaN. One aspect of NaN values is
+that by definition they will fail on @emph{any} comparison. Hence both
+equal and not-equal operators will fail when both their operands are NaN!
+The only way to guarantee that select blank pixels is through the
address@hidden operator explained above.
One way you can exploit this property of the NaN value to your advantage is
when you want a fully zero-valued image (even over the blank pixels) based
@@ -7490,9 +7524,10 @@ floating point number in Gnuastro isn't case-sensitive.
@node Invoking astarithmetic, , Arithmetic operators, Arithmetic
@subsection Invoking Arithmetic
-Arithmetic will do pixel to pixel arithmetic operations on the
-individual pixels of input data and/or numbers. All input data files
-must have the same dimensions. The general template is:
+Arithmetic will do pixel to pixel arithmetic operations on the individual
+pixels of input data and/or numbers. Any single element input (number or
+single pixel FITS image file) will be read as a number, the rest of the
+inputs must have the same dimensions. The general template is:
@example
$ astarithmetic [OPTION...] ASTRdata1 [ASTRdata2] OPERATOR ...
@@ -7515,35 +7550,61 @@ $ astarithmetic img1.fits img2.fits img3.fits median
\
Arithmetic accepts two kinds of input: images and numbers. Images are
considered to be any of the inputs that is a file name of a recognized type
-(see @ref{Arguments}). Numbers will be read into the smallest type that can
-store them, so @command{-2} will be read as a @code{char} type (which is
-signed on most systems and can thus keep negative values), @command{2500}
-will be read as an @code{unsigned short} (all positive numbers will be read
-as unsigned), while @code{3.14159265358979323846} will be read as a
address@hidden and @code{3.14} will be read as a @code{float}.
-
-All the operators will work on the native type of their operands, so the C
-programming language's internal type conversion will be used if the operand
-types of of multi-operand operators are different. Some operators can only
-work on integral types (for example bitwise operators) while others only
-work on floating point types, (for example @command{log}, or
address@hidden). In such cases if the operand type is different an error
-will be printed. Arithmetic also comes with internal type conversion
-operators which you can use to convert the data into the appropriate
-format, see @ref{Arithmetic operators}.
-
-If the output is an image, and the @option{--output} option is not
-given, automatic output will use the name of the first FITS image
-encountered to generate an output file name, see @ref{Automatic
-output}. Also, output WCS information will be taken from the first
-input image encountered. If the output is a single number, that number
-will be printed in the standard output. See @ref{Reverse polish
-notation} for the notation used to mix operands and operators on the
-command-line.
+(see @ref{Arguments}) and has more than one element. Numbers will be read
+into the smallest type that can store them, so @command{-2} will be read as
+a @code{char} type (which is signed on most systems and can thus keep
+negative values), @command{2500} will be read as an @code{unsigned short}
+(all positive numbers will be read as unsigned), while
address@hidden will be read as a @code{double} and
address@hidden will be read as a @code{float}. To force a number to be read as
+float, add a @code{f} after it, so @command{5f} will be added to the stack
+as @code{float} (see @ref{Reverse polish notation}).
+
+Unless otherwise stated (in @ref{Arithmetic operators}), the operators can
+deal with multiple datatypes. For example in address@hidden b.fits +}'',
+the image types can be @code{long} and @code{float}. In such cases, C's
+internal type conversion will be used, and the output will have the larger
+type of the two inputs. Unsigned integer types have smaller ranking than
+their signed counterparts and floating points types have higher ranking
+than the integer types. So the internal C type conversions done in the
+example above are equivalent to this piece of C:
address@hidden Mask image
-All the operators will work on the native types of the input data. Unless
-otherwise stated for each operator. To ignore certain pixels (see
address@hidden
+size_t i;
+long a[100];
+float b[100], out;
+for(i=0;i<100;++i) out[i]=a[i]+b[i];
address@hidden example
+
address@hidden
+Relying on the default C type conversion significantly speeds up the
+processing and also requires less RAM (when using very large
+images). However this great advantage comes at the cost of preparing for
+all the combinations of types while building/compiling Gnuastro, see the
+list of @option{--enable-bin-op-XXXX} options in @ref{Gnuastro configure
+options}. With the full list of types, compilation can take up to an hour
+or more. When a type isn't enabled for binary operators, the input data
+will be internally converted to the smallest larger type that was
+enabled. This can slow down your processing and consume more RAM, so if you
+often deal with data of a specific types, it is much better to make the
+one-time investment at compilation time and reap the benefits each time you
+run Gnuastro.
+
+Some operators can only work on integer types (of any length, for example
+bitwise operators) while others only work on floating point types,
+(currently only the @code{pow} operator). In such cases, if the operand
+type(s) are different an error will be printed and internal conversion
+won't occur. Arithmetic also comes with internal type conversion operators
+which you can use to convert the data into the appropriate format, see
address@hidden operators}.
+
+If the output is an image, and the @option{--output} option is not given,
+automatic output will use the name of the first FITS image encountered to
+generate an output file name, see @ref{Automatic output}. Also, output WCS
+information will be taken from the first input image encountered. If the
+output is a single number, that number will be printed in the standard
+output. See @ref{Reverse polish notation} for the notation used to mix
+operands and operators on the command-line. To ignore certain pixels (see
@ref{Blank pixels}), you can specify a mask image, see @ref{Mask image}. In
that case, when reading the first FITS image, all the masked pixels will be
set to a blank value depending on the type of the image.
diff --git a/lib/data-arithmetic-other.c b/lib/data-arithmetic-other.c
index 6b392b9..8c7d9bb 100644
--- a/lib/data-arithmetic-other.c
+++ b/lib/data-arithmetic-other.c
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <gnuastro/data.h>
+#include <gnuastro/qsort.h>
@@ -306,17 +307,21 @@ check_float_input(gal_data_t *in, int operator, char
*numstr)
/*************** Unary functions **************/
/***********************************************************************/
-#define UNIFUNC_MINVALUE(ISINT) { \
- int hasblank= (ISINT && gal_data_has_blank(in)) ? 1 : 0; \
- if(hasblank) \
+/* Note that for floating point types *b!=*b (by definition of NaN). And in
+ such cases, even if there are blank values, the smaller and larger
+ condition checked will fail, therefore the final result will be what we
+ want (to ignore the blank values). */
+#define UNIFUNC_MINVALUE { \
+ int blankeq = (*b==*b && gal_data_has_blank(in)) ? 1 : 0; \
+ if(blankeq) \
do if(*ia!=*b) *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf); \
else \
do *oa = *ia < *oa ? *ia : *oa; while(++ia<iaf); \
}
-#define UNIFUNC_MAXVALUE(ISINT) { \
- int hasblank= (ISINT && gal_data_has_blank(in)) ? 1 : 0; \
- if(hasblank) \
+#define UNIFUNC_MAXVALUE { \
+ int blankeq = (*b==*b && gal_data_has_blank(in)) ? 1 : 0; \
+ if(blankeq) \
do if(*ia!=*b) *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf); \
else \
do *oa = *ia > *oa ? *ia : *oa; while(++ia<iaf); \
@@ -331,22 +336,22 @@ check_float_input(gal_data_t *in, int operator, char
*numstr)
do *oa++ = OP(*ia++); while(ia<iaf); \
}
-#define UNIFUNC_RUN_FUNCTION_ON_ARRAY(IT, ISINT){ \
+#define UNIFUNC_RUN_FUNCTION_ON_ARRAY(IT){ \
+ IT *b = gal_data_alloc_blank(in->type); \
IT *ia=in->array, *oa=o->array, *iaf=ia + in->size; \
- IT *b = ISINT ? gal_data_alloc_blank(in->type) : NULL; \
switch(operator) \
{ \
case GAL_DATA_OPERATOR_MINVAL: \
- UNIFUNC_MINVALUE(ISINT); \
+ UNIFUNC_MINVALUE; \
break; \
case GAL_DATA_OPERATOR_MAXVAL: \
- UNIFUNC_MAXVALUE(ISINT); \
+ UNIFUNC_MAXVALUE; \
break; \
default: \
error(EXIT_FAILURE, 0, "the operator code %d is not " \
"recognized in UNIFUNC_RUN_FUNCTION_ON_ARRAY", operator); \
} \
- if(ISINT) free(b); \
+ free(b); \
}
@@ -402,37 +407,37 @@ check_float_input(gal_data_t *in, int operator, char
*numstr)
switch(in->type) \
{ \
case GAL_DATA_TYPE_UCHAR: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned char, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned char) \
break; \
case GAL_DATA_TYPE_CHAR: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(char, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(char) \
break; \
case GAL_DATA_TYPE_USHORT: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned short, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned short) \
break; \
case GAL_DATA_TYPE_SHORT: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(short, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(short) \
break; \
case GAL_DATA_TYPE_UINT: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned int, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned int) \
break; \
case GAL_DATA_TYPE_INT: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(int, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(int) \
break; \
case GAL_DATA_TYPE_ULONG: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned long, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(unsigned long) \
break; \
case GAL_DATA_TYPE_LONG: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(long, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(long) \
break; \
case GAL_DATA_TYPE_LONGLONG: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(LONGLONG, 1) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(LONGLONG) \
break; \
case GAL_DATA_TYPE_FLOAT: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(float, 0) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(float) \
break; \
case GAL_DATA_TYPE_DOUBLE: \
- UNIFUNC_RUN_FUNCTION_ON_ARRAY(double, 0) \
+ UNIFUNC_RUN_FUNCTION_ON_ARRAY(double) \
break; \
default: \
error(EXIT_FAILURE, 0, "type code %d not recognized in " \
@@ -833,3 +838,328 @@ data_arithmetic_where(unsigned char flags, gal_data_t
*out,
gal_data_free(iftrue);
}
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
+/*************** Multiple operand operators **************/
+/***********************************************************************/
+#define MULTIOPERAND_MIN(TYPE) { \
+ TYPE p, max; \
+ gal_data_type_max(list->type, &max); \
+ do /* Loop over each pixel */ \
+ { \
+ p=max; \
+ for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ { /* Only for integer types, does *b==*b. */ \
+ if(hasblank[i] && *b==*b) \
+ { if( *a[i] != *b ) p = *a[i] < p ? *a[i] : p; \
+ else p = *a[i] < p ? *a[i] : p; } \
+ ++a[i]; \
+ } \
+ *o++=p; \
+ } \
+ while(o<of); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_MAX(TYPE) { \
+ TYPE p, min; \
+ gal_data_type_min(list->type, &min); \
+ do /* Loop over each pixel */ \
+ { \
+ p=min; \
+ for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ { /* Only for integer types, does *b==*b. */ \
+ if(hasblank[i] && *b==*b) \
+ { if( *a[i] != *b ) p = *a[i] > p ? *a[i] : p; \
+ else p = *a[i] > p ? *a[i] : p; } \
+ ++a[i]; \
+ } \
+ *o++=p; \
+ } \
+ while(o<of); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_SUM { \
+ double sum; \
+ int n, use; \
+ do /* Loop over each pixel */ \
+ { \
+ n=0; \
+ use=1; \
+ sum=0.0f; \
+ for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ { \
+ /* Only integers and non-NaN floats: v==v is 1. */ \
+ if(hasblank[i]) \
+ use = ( *b==*b \
+ ? ( *a[i]!=*b ? 1 : 0 ) /* Integer */ \
+ : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ \
+ /* a[i] must be incremented to next pixel in any case. */ \
+ if(use) { sum += *a[i]++; ++n; } else ++a[i]; \
+ } \
+ *o++ = n ? sum : *b; \
+ } \
+ while(o<of); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_AVERAGE { \
+ double sum; \
+ int n, use; \
+ do /* Loop over each pixel */ \
+ { \
+ n=0; \
+ use=1; \
+ sum=0.0f; \
+ for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ { \
+ /* Only integers and non-NaN floats: v==v is 1. */ \
+ if(hasblank[i]) \
+ use = ( *b==*b \
+ ? ( *a[i]!=*b ? 1 : 0 ) /* Integer */ \
+ : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ \
+ /* a[i] must be incremented to next pixel in any case. */ \
+ if(use) { sum += *a[i]++; ++n; } else ++a[i]; \
+ } \
+ *o++ = n ? sum/n : *b; \
+ } \
+ while(o<of); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) { \
+ TYPE *pixs; \
+ int n, use; \
+ \
+ errno=0; \
+ pixs=malloc(dnum*sizeof *pixs); \
+ if(pixs==NULL) \
+ error(EXIT_FAILURE, 0, "%zu bytes in MULTIOPERAND_MEDIAN", \
+ dnum*sizeof *pixs); \
+ \
+ do /* Loop over each pixel */ \
+ { \
+ n=0; \
+ use=1; \
+ for(i=0;i<dnum;++i) /* Loop over each array. */ \
+ { \
+ /* Only integers and non-NaN floats: v==v is 1. */ \
+ if(hasblank[i]) \
+ use = ( *b==*b \
+ ? ( *a[i]!=*b ? 1 : 0 ) /* Integer */ \
+ : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ \
+ /* a[i] must be incremented to next pixel in any case. */ \
+ if(use) pixs[n++]=*a[i]++; else ++a[i]; \
+ } \
+ \
+ /* Sort all the values for this pixel and return the median. */ \
+ if(n) \
+ { \
+ qsort(pixs, n, sizeof *pixs, QSORT_F); \
+ *o++ = n%2 ? pixs[n/2] : (pixs[n/2] + pixs[n/2-1])/2 ; \
+ } \
+ else \
+ *o++=*b; \
+ } \
+ while(o<of); \
+ }
+
+
+
+
+
+#define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) { \
+ TYPE *o=out->array, *of=out->array+out->size; \
+ TYPE **a, *b=gal_data_alloc_blank(list->type); \
+ \
+ /* Allocate space to keep the pointers to the arrays of each. */ \
+ /* Input data structure. The operators will increment these */ \
+ /* pointers while parsing them. */ \
+ errno=0; \
+ a=malloc(dnum*sizeof *a); \
+ if(a==NULL) \
+ error(EXIT_FAILURE, 0, "%zu bytes for `arrays' in " \
+ "MULTIOPERAND_TYPE_SET", dnum*sizeof *a); \
+ \
+ /* Fill in the array pointers. */ \
+ for(tmp=list;tmp!=NULL;tmp=tmp->next) \
+ a[i++]=tmp->array; \
+ \
+ /* Do the operation. */ \
+ switch(operator) \
+ { \
+ case GAL_DATA_OPERATOR_MIN: \
+ MULTIOPERAND_MIN(TYPE); \
+ break; \
+ \
+ case GAL_DATA_OPERATOR_MAX: \
+ MULTIOPERAND_MAX(TYPE); \
+ break; \
+ \
+ case GAL_DATA_OPERATOR_SUM: \
+ MULTIOPERAND_SUM; \
+ break; \
+ \
+ case GAL_DATA_OPERATOR_AVERAGE: \
+ MULTIOPERAND_AVERAGE; \
+ break; \
+ \
+ case GAL_DATA_OPERATOR_MEDIAN: \
+ MULTIOPERAND_MEDIAN(TYPE, QSORT_F); \
+ break; \
+ \
+ default: \
+ error(EXIT_FAILURE, 0, "the operator code %d not recognized " \
+ "in MULTIOPERAND_TYPE_SET", operator); \
+ } \
+ \
+ /* Clean up. */ \
+ free(b); \
+ free(a); \
+ }
+
+
+
+
+
+/* The single operator in this function is assumed to be a linked list. The
+ number of operators is determined from the fact that the last node in
+ the linked list must have a NULL pointer as its `next' element.*/
+gal_data_t *
+data_arithmetic_multioperand(int operator, unsigned char flags,
+ gal_data_t *list)
+{
+ int *hasblank;
+ size_t i=0, dnum=1;
+ gal_data_t *out, *tmp, *ttmp;
+
+
+ /* For generality, `list' can be a NULL pointer, in that case, this
+ function will return a NULL pointer and avoid further processing. */
+ if(list==NULL) return NULL;
+
+
+ /* Do a simple sanity check, comparing the operand on top of the list to
+ the rest of the operands within the list. */
+ for(tmp=list->next;tmp!=NULL;tmp=tmp->next)
+ {
+ /* Increment the number of structures. */
+ ++dnum;
+
+ /* Check the types. */
+ if(tmp->type!=list->type)
+ error(EXIT_FAILURE, 0, "the types of all operands to the %s "
+ "operator must be same",
+ gal_data_operator_string(operator));
+
+ /* Check the sizes. */
+ if( gal_data_dsize_is_different(list, tmp) )
+ error(EXIT_FAILURE, 0, "the sizes of all operands to the %s "
+ "operator must be same",
+ gal_data_operator_string(operator));
+ }
+
+
+ /* Set the output data structure. */
+ if(flags & GAL_DATA_ARITH_INPLACE)
+ out = list; /* The top element in the list. */
+ else
+ out = gal_data_alloc(NULL, list->type, list->ndim, list->dsize,
+ list->wcs, 0, list->minmapsize);
+
+
+ /* hasblank is used to see if a blank value should be checked or not. */
+ errno=0;
+ hasblank=malloc(dnum*sizeof *hasblank);
+ if(hasblank==NULL)
+ error(EXIT_FAILURE, 0, "%zu bytes for hasblank in "
+ "`data_arithmetic_multioperand", dnum*sizeof *hasblank);
+
+
+ /* Fill in the hasblank array. */
+ for(tmp=list;tmp!=NULL;tmp=tmp->next)
+ hasblank[i++]=gal_data_has_blank(tmp);
+
+
+ /* Start the operation. */
+ switch(list->type)
+ {
+ case GAL_DATA_TYPE_UCHAR:
+ MULTIOPERAND_TYPE_SET(unsigned char, gal_qsort_uchar_increasing);
+ case GAL_DATA_TYPE_CHAR:
+ MULTIOPERAND_TYPE_SET(char, gal_qsort_char_increasing);
+ case GAL_DATA_TYPE_USHORT:
+ MULTIOPERAND_TYPE_SET(unsigned short, gal_qsort_ushort_increasing);
+ case GAL_DATA_TYPE_SHORT:
+ MULTIOPERAND_TYPE_SET(short, gal_qsort_short_increasing);
+ case GAL_DATA_TYPE_UINT:
+ MULTIOPERAND_TYPE_SET(unsigned int, gal_qsort_uint_increasing);
+ case GAL_DATA_TYPE_INT:
+ MULTIOPERAND_TYPE_SET(int, gal_qsort_int_increasing);
+ case GAL_DATA_TYPE_ULONG:
+ MULTIOPERAND_TYPE_SET(unsigned long, gal_qsort_ulong_increasing);
+ case GAL_DATA_TYPE_LONG:
+ MULTIOPERAND_TYPE_SET(long, gal_qsort_long_increasing);
+ case GAL_DATA_TYPE_LONGLONG:
+ MULTIOPERAND_TYPE_SET(LONGLONG, gal_qsort_longlong_increasing);
+ case GAL_DATA_TYPE_FLOAT:
+ MULTIOPERAND_TYPE_SET(float, gal_qsort_float_increasing);
+ case GAL_DATA_TYPE_DOUBLE:
+ MULTIOPERAND_TYPE_SET(double, gal_qsort_double_increasing);
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`data_arithmetic_multioperand'", list->type);
+ }
+
+
+ /* Clean up and return. Note that the operation might have been done in
+ place. In that case, the top most list element was used. So we need to
+ check before freeing each data structure. */
+ if(flags & GAL_DATA_ARITH_FREE)
+ {
+ tmp=list;
+ while(tmp!=NULL)
+ {
+ ttmp=tmp->next;
+ if(tmp!=out) gal_data_free(tmp);
+ tmp=ttmp;
+ }
+ }
+ free(hasblank);
+ return out;
+}
diff --git a/lib/data-arithmetic-other.h b/lib/data-arithmetic-other.h
index f233d43..dec1451 100644
--- a/lib/data-arithmetic-other.h
+++ b/lib/data-arithmetic-other.h
@@ -45,4 +45,8 @@ void
data_arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond, gal_data_t *iftrue);
+gal_data_t *
+data_arithmetic_multioperand(int operator, unsigned char flags,
+ gal_data_t *list);
+
#endif
diff --git a/lib/data.c b/lib/data.c
index 77af0b4..f862bca 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -490,6 +490,111 @@ gal_data_free(gal_data_t *data)
+
+
+/*********************************************************************/
+/************* Data structure as a linked list ******************/
+/*********************************************************************/
+/* Add a new data structure to the top of an existing linked list of data
+ structures. Note that if the new node is its self a list, all its nodes
+ will be added to the list. */
+void
+gal_data_add_to_ll(gal_data_t **list, gal_data_t *newnode)
+{
+ gal_data_t *tmp=newnode, *toadd;
+
+ /* Check if newnode is itself a list or not. */
+ if(newnode->next)
+ {
+ /* Go onto the last node in newnode's existing list. */
+ while(tmp->next) tmp=tmp->next;
+
+ /* Set the last node as the node to add to the list. */
+ toadd=tmp;
+ }
+ else
+ /* Its not a list, so just set it to `toadd'. */
+ toadd=newnode;
+
+
+ /* Set the next element of toadd and update what list points to.*/
+ toadd->next=*list;
+ *list=toadd;
+}
+
+
+
+
+
+gal_data_t *
+gal_data_pop_from_ll(gal_data_t **list)
+{
+ struct gal_data_t *out;
+ out=*list;
+ *list=out->next;
+ return out;
+}
+
+
+
+
+
+size_t
+gal_data_num_in_ll(gal_data_t *list)
+{
+ size_t num=0;
+ while(list!=NULL)
+ {
+ ++num;
+ list=list->next;
+ }
+ return num;
+}
+
+
+
+
+
+gal_data_t **
+gal_data_ll_to_array_of_ptrs(gal_data_t *list, size_t *num)
+{
+ size_t i=0;
+ gal_data_t **out;
+
+ /* Find the number of nodes in the list. */
+ *num=gal_data_num_in_ll(list);
+
+ /* Allocate space for the array. */
+ errno=0;
+ out=malloc(*num * sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, 0, "%zu bytes for the output pointer in "
+ "`gal_data_ll_to_array_of_ptrs'", *num * sizeof *out);
+
+ /* Fill in the array with pointers to each data-structure. Note that we
+ don't need the list pointer any more, so we can just increment it.*/
+ while(list!=NULL) { out[i++]=list; list=list->next; }
+
+ /* Return the allocated array. */
+ return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/*************************************************************
************** Blank data ***************
*************************************************************/
@@ -1449,6 +1554,9 @@ gal_data_string_to_number(char *string)
/*************************************************************
************** Type minimum and maximums ***************
*************************************************************/
+/* Put the minimum (or maximum for the `gal_data_type_max') value for the
+ type in the space (that must already be allocated before the call to
+ this function) pointed to by in. */
void
gal_data_type_min(int type, void *in)
{
@@ -1655,6 +1763,15 @@ gal_data_arithmetic(int operator, unsigned char flags,
...)
out=data_arithmetic_abs(flags, d1);
break;
+ case GAL_DATA_OPERATOR_MIN:
+ case GAL_DATA_OPERATOR_MAX:
+ case GAL_DATA_OPERATOR_SUM:
+ case GAL_DATA_OPERATOR_AVERAGE:
+ case GAL_DATA_OPERATOR_MEDIAN:
+ d1 = va_arg(va, gal_data_t *);
+ out=data_arithmetic_multioperand(operator, flags, d1);
+ break;
+
/* Binary function operators. */
case GAL_DATA_OPERATOR_POW:
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 74e5939..138e28f 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -167,6 +167,7 @@ enum gal_data_operators
GAL_DATA_OPERATOR_MAXVAL, /* Maximum value of array. */
GAL_DATA_OPERATOR_MIN, /* Minimum per pixel of multiple arrays. */
GAL_DATA_OPERATOR_MAX, /* Maximum per pixel of multiple arrays. */
+ GAL_DATA_OPERATOR_SUM, /* Sum per pixel of multiple arrays. */
GAL_DATA_OPERATOR_AVERAGE, /* Average per pixel of multiple arrays. */
GAL_DATA_OPERATOR_MEDIAN, /* Median per pixel of multiple arrays. */
@@ -241,6 +242,9 @@ gal_data_malloc_array(int type, size_t size);
void *
gal_data_calloc_array(int type, size_t size);
+void *
+gal_data_alloc_number(int type, void *number);
+
gal_data_t *
gal_data_alloc(void *array, int type, size_t ndim, long *dsize,
struct wcsprm *wcs, int clear, size_t minmapsize);
@@ -252,6 +256,25 @@ gal_data_free(gal_data_t *data);
+/*********************************************************************/
+/************* Data structure as a linked list ******************/
+/*********************************************************************/
+void
+gal_data_add_to_ll(gal_data_t **list, gal_data_t *newnode);
+
+gal_data_t *
+gal_data_pop_from_ll(struct gal_data_t **list);
+
+size_t
+gal_data_num_in_ll(struct gal_data_t *list);
+
+gal_data_t **
+gal_data_ll_to_array_of_ptrs(gal_data_t *list, size_t *num);
+
+
+
+
+
/*************************************************************
************** Blank data ***************
diff --git a/lib/gnuastro/qsort.h b/lib/gnuastro/qsort.h
index 2d71b43..cd59910 100644
--- a/lib/gnuastro/qsort.h
+++ b/lib/gnuastro/qsort.h
@@ -55,6 +55,40 @@ extern float *gal_qsort_index_arr;
int
gal_qsort_index_float_decreasing(const void * a, const void * b);
+
+
+
+
+int
+gal_qsort_uchar_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_uchar_increasing(const void * a, const void * b);
+
+int
+gal_qsort_char_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_char_increasing(const void * a, const void * b);
+
+int
+gal_qsort_ushort_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_ushort_increasing(const void * a, const void * b);
+
+int
+gal_qsort_short_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_short_increasing(const void * a, const void * b);
+
+int
+gal_qsort_uint_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_uint_increasing(const void * a, const void * b);
+
int
gal_qsort_int_decreasing(const void * a, const void * b);
@@ -62,6 +96,24 @@ int
gal_qsort_int_increasing(const void * a, const void * b);
int
+gal_qsort_ulong_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_ulong_increasing(const void * a, const void * b);
+
+int
+gal_qsort_long_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_long_increasing(const void * a, const void * b);
+
+int
+gal_qsort_longlong_decreasing(const void * a, const void * b);
+
+int
+gal_qsort_longlong_increasing(const void * a, const void * b);
+
+int
gal_qsort_float_decreasing(const void * a, const void * b);
int
diff --git a/lib/qsort.c b/lib/qsort.c
index 655cadc..889cd94 100644
--- a/lib/qsort.c
+++ b/lib/qsort.c
@@ -24,6 +24,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
+#include <fitsio.h>
+
#include <gnuastro/qsort.h>
/* Initialize the array for sorting indexs to NULL. */
@@ -40,6 +42,72 @@ gal_qsort_index_float_decreasing(const void * a, const void
* b)
+
+
+
+
+
+
+int
+gal_qsort_uchar_decreasing(const void * a, const void * b)
+{
+ return ( *(unsigned char*)b - *(unsigned char*)a );
+}
+
+int
+gal_qsort_uchar_increasing(const void * a, const void * b)
+{
+ return ( *(unsigned char*)a - *(unsigned char*)b );
+}
+
+int
+gal_qsort_char_decreasing(const void * a, const void * b)
+{
+ return ( *(char*)b - *(char*)a );
+}
+
+int
+gal_qsort_char_increasing(const void * a, const void * b)
+{
+ return ( *(char*)a - *(char*)b );
+}
+
+int
+gal_qsort_ushort_decreasing(const void * a, const void * b)
+{
+ return ( *(unsigned short*)b - *(unsigned short*)a );
+}
+
+int
+gal_qsort_ushort_increasing(const void * a, const void * b)
+{
+ return ( *(unsigned short*)a - *(unsigned short*)b );
+}
+
+int
+gal_qsort_short_decreasing(const void * a, const void * b)
+{
+ return ( *(short*)b - *(short*)a );
+}
+
+int
+gal_qsort_short_increasing(const void * a, const void * b)
+{
+ return ( *(short*)a - *(short*)b );
+}
+
+int
+gal_qsort_uint_decreasing(const void * a, const void * b)
+{
+ return ( *(unsigned int*)b - *(unsigned int*)a );
+}
+
+int
+gal_qsort_uint_increasing(const void * a, const void * b)
+{
+ return ( *(unsigned int*)a - *(unsigned int*)b );
+}
+
int
gal_qsort_int_decreasing(const void * a, const void * b)
{
@@ -53,6 +121,43 @@ gal_qsort_int_increasing(const void * a, const void * b)
}
int
+gal_qsort_ulong_decreasing(const void * a, const void * b)
+{
+ return ( *(unsigned long*)b - *(unsigned long*)a );
+}
+
+int
+gal_qsort_ulong_increasing(const void * a, const void * b)
+{
+ return ( *(unsigned long*)a - *(unsigned long*)b );
+}
+
+
+int
+gal_qsort_long_decreasing(const void * a, const void * b)
+{
+ return ( *(long*)b - *(long*)a );
+}
+
+int
+gal_qsort_long_increasing(const void * a, const void * b)
+{
+ return ( *(long*)a - *(long*)b );
+}
+
+int
+gal_qsort_longlong_decreasing(const void * a, const void * b)
+{
+ return ( *(LONGLONG*)b - *(LONGLONG*)a );
+}
+
+int
+gal_qsort_longlong_increasing(const void * a, const void * b)
+{
+ return ( *(LONGLONG*)a - *(LONGLONG*)b );
+}
+
+int
gal_qsort_float_decreasing(const void * a, const void * b)
{
float ta=*(float*)a;
- [gnuastro-commits] master 23a541a 056/125: Single correction of string keyword values in FITS tables, (continued)
- [gnuastro-commits] master 23a541a 056/125: Single correction of string keyword values in FITS tables, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4c28d13 042/125: Problem in reading blank FITS ASCII table fixed, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 08927b8 044/125: New Table formats section in manual, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master eab25b1 059/125: Option descriptions also printed with values, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 909fa0d 050/125: Table info printing in libraries, updates to Table program, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 36bcedd 066/125: Fixed automatic output checking in Arithmetic, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master d2ed7ea 060/125: Options setdirconf and setusrconf implemented, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 08147ce 062/125: All mandatory options not given are listed with error, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master a9092e9 036/125: The table library can read ASCII inputs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 77b7910 057/125: Program specific global variables for in options library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 9e553f0 022/125: All old arithmetic operators are now implemented,
Mohammad Akhlaghi <=
- [gnuastro-commits] master b470ee4 023/125: More efficient macro implementation for binary operators, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 6b51397 077/125: Added -j8 to make examples in book, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master b5c9a64 063/125: Science and its tools and ImageCrop corrections in book, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master edd5612 053/125: Single function to read all configuration files, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master acad34c 049/125: Sanity checks for reading txt tables, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 5b61ae3 032/125: Arithmetic functions/macros out of data.h, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 232e059 041/125: First draft of FITS table writing function complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 56f8c18 071/125: Cosmetic changes in option printing style, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 07a043c 026/125: Starting to use the new data structure to Table, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master f47a4c2 072/125: MakeProfiles now reads the columns, Mohammad Akhlaghi, 2017/04/23