[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9959e45 004/125: Arithmetic using new data str
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 9959e45 004/125: Arithmetic using new data structure |
|
Date: |
Sun, 23 Apr 2017 22:36:24 -0400 (EDT) |
branch: master
commit 9959e45d7cf4d3788b40f36c1661d60f27f1b0a9
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Arithmetic using new data structure
Many functions were changed/added in `fits.h' and `data.h', so when
this branch is complete, go through the full header and update
them. Since the job is ongoing now and further changes are possible
as the data structure is implemented in all the utilities, other
changes might be made soon, so the manual is not being updated with
this commit.
The new data structure is now fully implemented in Arithmetic. Also,
there is now the possibility that data can be stored on the hdd/ssd
instead of kept in RAM. Several changes were made that are listed
below:
- The operands and operators functions were moved to two separate
`.c' and `.h' files in Arithmetic's source: `operands.[ch]' and
`operators.[ch]'.
- The man-pages were built unconditionally, even when Gnuastro was
configured to only build a limited number of utilities. Now, the
man-page will only be created if the program was built.
- A new set of internal library functions for changing data type
arrays was defined in `lib/changetype.[ch]'. This was done to keep
the `lib/data.c' file clean from these repetative functions and
currently its not intended for installation, there is a high-level
`gal_data_copy_to_new_type' for use by library users.
Changed option values:
- `--type' in Arithmetic and MakeProfiles now uses `uchar' instead of
`byte' for the `unsigned char' type array (equivalent to CFITSIO's
`BYTE_IMG').
New user-visible macros:
- GAL_CONFIG_HAVE_WCSLIB_VERSION
New user-visible functions:
- gal_linkedlist_free_stll
---
bin/arithmetic/Makefile.am | 2 +-
bin/arithmetic/args.h | 8 +-
bin/arithmetic/arithmetic.c | 1040 +----------------------------
bin/arithmetic/arithmetic.h | 4 +-
bin/arithmetic/main.c | 7 +-
bin/arithmetic/main.h | 35 +-
bin/arithmetic/operands.c | 177 +++++
bin/arithmetic/{ui.h => operands.h} | 19 +-
bin/arithmetic/operators.c | 826 +++++++++++++++++++++++
bin/arithmetic/operators.h | 96 +++
bin/arithmetic/ui.c | 45 +-
bin/arithmetic/ui.h | 7 +-
bin/mkprof/args.h | 2 +-
configure.ac | 9 +-
doc/Makefile.am | 60 +-
lib/Makefile.am | 14 +-
lib/changetype.c | 363 ++++++++++
lib/changetype.h | 60 ++
lib/checkset.c | 43 +-
lib/checkset.h | 5 +-
lib/config.h.in | 1 +
lib/data.c | 633 +++++++++++++++++-
lib/fits.c | 1244 +++++++++++------------------------
lib/gnuastro/data.h | 59 +-
lib/gnuastro/fits.h | 199 +++---
lib/gnuastro/linkedlist.h | 3 +
lib/gnuastro/threads.h | 15 +-
lib/linkedlist.c | 18 +
lib/mesh.c | 47 +-
29 files changed, 2964 insertions(+), 2077 deletions(-)
diff --git a/bin/arithmetic/Makefile.am b/bin/arithmetic/Makefile.am
index 97a004a..cf714b6 100644
--- a/bin/arithmetic/Makefile.am
+++ b/bin/arithmetic/Makefile.am
@@ -29,7 +29,7 @@ AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib
bin_PROGRAMS = astarithmetic
astarithmetic_SOURCES = main.c main.h cite.h ui.c ui.h args.h \
-arithmetic.c arithmetic.h
+arithmetic.c arithmetic.h operands.c operands.h operators.c operators.h
astarithmetic_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la \
-lgnuastro
diff --git a/bin/arithmetic/args.h b/bin/arithmetic/args.h
index 135732e..2be11b8 100644
--- a/bin/arithmetic/args.h
+++ b/bin/arithmetic/args.h
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -137,7 +137,7 @@ static struct argp_option options[] =
'T',
"STR",
0,
- "byte, short, long, longlong, float, double.",
+ "uchar, short, long, longlong, float, double.",
2
},
@@ -207,7 +207,7 @@ parse_opt(int key, char *arg, struct argp_state *state)
/* Output */
case 'T':
- gal_checkset_known_types(arg, &p->type, NULL, 0);
+ gal_checkset_known_types(arg, &p->outtype, NULL, 0);
p->up.typeset=1;
break;
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index cff9ee4..efb10aa 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -36,981 +36,10 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <checkset.h>
#include "main.h"
-#include "arithmetic.h" /* needs main.h. */
-
-
-
-
-/***************************************************************/
-/************* Operand linked list functions *************/
-/***************************************************************/
-size_t
-num_operands(struct imgarithparams *p)
-{
- size_t counter=0;
- struct operand *tmp=NULL;
- for(tmp=p->operands;tmp!=NULL;tmp=tmp->next)
- ++counter;
- return counter;
-}
-
-
-
-
-
-void
-add_operand(struct imgarithparams *p, char *filename, double number,
- double *array)
-{
- 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->array=array;
- newnode->number=number;
- 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;
- }
-
- /* Make the link to the previous list. */
- newnode->next=p->operands;
- p->operands=newnode;
-}
-
-
-
-
-
-void
-pop_operand(struct imgarithparams *p, double *number, double **array,
- char *operator)
-{
- int bitpix;
- size_t s0, s1;
- struct uiparams *up=&p->up;
- struct operand *operands=p->operands;
- char *maskname, *mhdu, *filename, *hdu;
-
- /* If the operand linked list has finished, then give an error and
- exit. */
- if(operands==NULL)
- error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" operator",
- operator);
-
-
- /* Set the array output. If filename is present then read the file
- and fill in the array, if not then just set the array. */
- if(strlen(operands->filename))
- {
- hdu=operands->hdu;
- filename=operands->filename;
-
- /* In case this is the first image that is read, then read the
- WCS information and set the mask name so masked pixels can be
- set to NaN. For the other images, the mask can be completely
- ignored. */
- if(p->popcounter) /* This is not the first FITS file. */
- {
- maskname=NULL;
- mhdu=NULL;
- }
- else
- {
- mhdu=up->mhdu;
- maskname=up->maskname;
- gal_fits_read_wcs(filename, hdu, 0, 0, &p->nwcs, &p->wcs);
- }
- gal_fits_file_to_double(filename, maskname, hdu, mhdu,
- array, &bitpix, &p->anyblank, &s0, &s1);
-
- /* If the output size was not set yet, then set it. Otherwise,
- make sure the size of this image is the same as the previous
- images. */
- if(p->s0==0 && p->s1==0)
- {
- p->s0=s0;
- p->s1=s1;
- }
- else
- {
- if(p->s0!=s0 || p->s1!=s1)
- error(EXIT_FAILURE, 0, "%s (hdu=%s): has size of %zu x %zu. "
- "However, previous images had a size of %zu x %zu. All "
- "the images must be the same size in order for "
- "ImageArithmetic to work", filename, hdu, s0, s1,
- p->s0, p->s1);
- }
-
- /* Free the HDU string: */
- free(hdu);
-
- /* Add to the number of popped FITS images: */
- ++p->popcounter;
-
- /* Report the read image if desired: */
- if(p->cp.verb) printf("%s is read.\n", filename);
- }
- else
- *array=operands->array;
-
- /* Set the number: */
- *number=operands->number;
-
- /* Remove this node from the queue. */
- p->operands=operands->next;
- free(operands);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/***************************************************************/
-/************* Operators *************/
-/***************************************************************/
-void
-sum(struct imgarithparams *p)
-{
- size_t size;
- char *operator="+";
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr && sarr) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in the first
- input. */
- gal_array_dsum_arrays(farr, sarr, size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
-
- /* Clean up. */
- free(sarr);
- }
- else if(farr) /* Only the first is an array. */
- {
- gal_array_dsum_const(farr, size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- gal_array_dsum_const(sarr, size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, fnum+snum, NOOPTARRAY);
-}
-
-
-
-
-
-void
-subtract(struct imgarithparams *p)
-{
- size_t size;
- char *operator="-";
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr && sarr) /* 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(sarr, farr, size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
-
- /* Clean up. */
- free(farr);
- }
- else if(farr) /* Only the first is an array. */
- {
- gal_array_dconst_subtract(farr, size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- gal_array_dsubtract_const(sarr, size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, snum-fnum, NOOPTARRAY);
-}
-
-
-
-
-
-void
-multiply(struct imgarithparams *p)
-{
- size_t size;
- char *operator="*";
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr && sarr) /* Both are arrays. */
- {
- /* Do the operation, note that the output is stored in farr. */
- gal_array_dmultip_arrays(farr, sarr, size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
-
- /* Clean up. */
- free(sarr);
- }
- else if(farr) /* Only the first is an array. */
- {
- gal_array_dmultip_const(farr, size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- gal_array_dmultip_const(sarr, size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, fnum*snum, NOOPTARRAY);
-}
-
-
-
-
-
-void
-divide(struct imgarithparams *p)
-{
- size_t size;
- char *operator="/";
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr && sarr) /* 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(sarr, farr, size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
-
- /* Clean up. */
- free(farr);
- }
- else if(farr) /* Only the first is an array. */
- {
- gal_array_dconst_divide(farr, size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- gal_array_ddivide_const(sarr, size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, snum/fnum, NOOPTARRAY);
-}
-
-
-
-
-
-void
-topower(struct imgarithparams *p, char *op)
-{
- size_t size;
- char *operator = op?op:"pow";
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr && sarr) /* 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(sarr, farr, size);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
-
-
- /* Clean up. */
- free(farr);
- }
- else if(farr) /* Only the first is an array. */
- {
- gal_array_dconst_power(farr, size, snum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- gal_array_dpower_const(sarr, size, fnum);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, pow(snum, fnum), NOOPTARRAY);
-}
-
-
-
-
-
-void
-alloppixs(struct imgarithparams *p, char *operator)
-{
- int i, j; /* Integer to allow negative checks. */
- double num; /* Temporary number holder. */
- double *arr; /* Temporary array holder. */
- int firstisarray=-1; /* ==-1: unset. ==1: array. ==0: number. */
- double *allpixels=NULL; /* Array for all values in one pixel. */
- double **allarrays=NULL; /* Array for pointers to input arrays. */
- size_t size, 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;
- 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, &arr, operator);
-
- /* Do the appropriate action if it is an array or a number. */
- if(arr)
- switch(firstisarray)
- {
- case -1:
- firstisarray=1;
- allarrays[i]=arr;
- break;
- case 1:
- allarrays[i]=arr;
- 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 firstisarray (%d) "
- "in the alloppixs function is not recognized",
- PACKAGE_BUGREPORT, firstisarray);
- }
- else
- switch(firstisarray)
- {
- case -1:
- firstisarray=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 firstisarray (%d) "
- "in the alloppixs function is not recognized",
- PACKAGE_BUGREPORT, firstisarray);
- }
- }
-
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
-
- /* Do the operation and report the result: */
- if(arr)
- {
- /* Find the value and replace it with the first operand. */
- for(i=0;i<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, allarrays[0]);
-
- /* Free all the extra operands */
- for(i=1;i<numop;++i) free(allarrays[i]);
- }
- else
- add_operand(p, NOOPTFILENAME, (*thisfunction)(allpixels, numop),
- NOOPTARRAY);
-
- /* Clean up: */
- 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, NOOPTARRAY);
-
- /* Call the power operator. */
- topower(p, operator);
-}
-
-
-
-
-
-void
-takelog(struct imgarithparams *p)
-{
- char *operator="log";
- double fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr) /* 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(farr, p->s0*p->s1);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, log(fnum), NOOPTARRAY);
-}
-
-
-
-
-
-void
-takelog10(struct imgarithparams *p)
-{
- char *operator="log10";
- double fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr) /* 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(farr, p->s0*p->s1);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, log10(fnum), NOOPTARRAY);
-}
-
-
-
-
-
-void
-takeabs(struct imgarithparams *p)
-{
- char *operator="abs";
- double fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr) /* 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(farr, p->s0*p->s1);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fabs(fnum), NOOPTARRAY);
-}
-
-
-
-
-
-void
-findmin(struct imgarithparams *p)
-{
- char *operator="min";
- double min, fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr) /* 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(farr, p->s0*p->s1, &min);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, min, NOOPTARRAY);
-
- /* Clean up. */
- free(farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fnum, NOOPTARRAY);
-}
-
-
-
-
-
-void
-findmax(struct imgarithparams *p)
-{
- char *operator="max";
- double max, fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr) /* 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(farr, p->s0*p->s1, &max);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, max, NOOPTARRAY);
-
- /* Clean up. */
- free(farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, fnum, NOOPTARRAY);
-}
-
-
-
-
-
-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)
-{
- size_t size;
- double fnum, snum; /* First or second number. */
- double *farr, *sarr; /* First or second array. */
- double *f, *s, *ff, *ss;
- int (*thisfunction)(double, double)=NULL;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- 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(farr && sarr) /* 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. */
- f=farr;
- ss=(s=sarr)+size;
- do *s = thisfunction(*s, *f++); while(++s<ss);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
-
- /* Clean up. */
- free(farr);
- }
- else if(farr) /* Only the first is an array. */
- {
- ff=(f=farr)+size;
- do *f = thisfunction(snum, *f); while(++f<ff);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else if(sarr) /* Only the first is an array. */
- {
- ss=(s=sarr)+size;
- do *s = thisfunction(*s, fnum); while(++s<ss);
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, sarr);
- }
- else /* Both are numbers. */
- add_operand(p, NOOPTFILENAME, thisfunction(snum, fnum), NOOPTARRAY);
-}
-
-
-
-
-
-void
-andor(struct imgarithparams *p, char *operator)
-{
- double *f, *s, *ff;
- double fnum, snum, *farr, *sarr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
- pop_operand(p, &snum, &sarr, 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(farr && sarr)
- {
- /* 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. */
- s = sarr;
- ff = (f=farr) + p->s0*p->s1;
- if(!strcmp(operator, "and"))
- do *f = *s++ && *f; while(++f<ff);
- else
- do *f = *s++ || *f; while(++f<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
-
- /* Clean up */
- free(sarr);
- }
- else if(farr==NULL || sarr==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,
- NOOPTARRAY);
-}
-
-
-
-
-void
-notfunc(struct imgarithparams *p)
-{
- double *f, *ff;
- double fnum, *farr;
- char *operator="not";
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Do the operation: */
- if(farr)
- {
- /* Fill the array with the output values. */
- ff = (f=farr) + p->s0*p->s1;
- do *f = !(*f); while(++f<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else
- add_operand(p, NOOPTFILENAME, !fnum, NOOPTARRAY);
-}
-
-
-
-
-
-/* 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)
-{
- size_t size;
- char *operator="isblank";
- double *f, *ff, fnum, *farr;
-
- /* Pop out the number of operands needed. */
- pop_operand(p, &fnum, &farr, operator);
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(farr) /* 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=(f=farr)+size;
- do *f = isnan(*f); while(++f<ff);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, farr);
- }
- else /* Operand is a number. */
- add_operand(p, NOOPTFILENAME, isnan(fnum), NOOPTARRAY);
-}
-
-
-
-
-
-/* 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)
-{
- size_t size;
- double *n, *c, *i, *ii;
- double nnum, cnum, inum; /* First, second, or third number. */
- double *narr, *carr, *iarr; /* First, second, or third array. */
-
- /* ORDER IS VERY IMPORTANT HERE. Pop out the number of operands needed. */
- pop_operand(p, &nnum, &narr, "where"); /* New value. */
- pop_operand(p, &cnum, &carr, "where"); /* Condition. */
- pop_operand(p, &inum, &iarr, "where"); /* Input. */
-
- /* Set the total number of pixels, note that we can't do this in the
- definition of the variable because p->s0 and p->s1 will be set in
- pop_operand for the first image. */
- size=p->s0*p->s1;
-
- /* Do the operation: */
- if(iarr && carr) /* Both are arrays. */
- {
- c=carr;
- ii=(i=iarr)+size;
- if(narr) /* `new' is an array, not number. */
- {
- n=narr;
- /* Note that we need to increment "f" after the check and
- replace. If the increment is inside the conditional replace,
- then when t==0, the increment won't work.*/
- do {*i = *c++ ? *n : *i; ++n;} while(++i<ii);
- }
- else /* `new' is a number, not array. */
- do *i = *c++ ? nnum : *i; while(++i<ii);
-
- /* Push the output onto the stack. */
- add_operand(p, NOOPTFILENAME, NOOPTNUMBER, iarr);
-
- /* Clean up. */
- free(carr);
- free(narr);
- }
- else if ( iarr==NULL && carr==NULL && narr==NULL )
- add_operand(p, NOOPTFILENAME, cnum ? nnum : inum, NOOPTARRAY);
- 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.");
-}
-
+#include "operands.h"
+#include "operators.h"
+#include "arithmetic.h"
@@ -1038,13 +67,11 @@ where(struct imgarithparams *p)
void
reversepolish(struct imgarithparams *p)
{
- void *array;
double number;
- char *tokeepvalue;
+ gal_data_t *ddata;
struct gal_linkedlist_stll *token;
/* Prepare the processing: */
- p->s0=p->s1=0;
p->operands=NULL;
p->addcounter=p->popcounter=0;
@@ -1055,9 +82,9 @@ reversepolish(struct imgarithparams *p)
linked list. Otherwise, pull out two members and do the
specified operation on them. */
if(gal_fits_name_is_fits(token->v))
- add_operand(p, token->v, NOOPTNUMBER, NOOPTARRAY);
+ add_operand(p, token->v, NOOPTNUMBER, NOOPTDATA);
else if(gal_checkset_str_is_double(token->v, &number))
- add_operand(p, NOOPTFILENAME, number, NOOPTARRAY);
+ add_operand(p, NOOPTFILENAME, number, NOOPTDATA);
else
{
if (!strcmp(token->v, "+")) sum(p);
@@ -1092,48 +119,46 @@ reversepolish(struct imgarithparams *p)
}
}
- /* If there is more than one node in the operands stack, then the
- user has given too many operands and there is an error. */
+ /* If there is more than one node in the operands stack then the user has
+ given too many operands which is an error. */
if(p->operands->next!=NULL)
- error(EXIT_FAILURE, 0, "there are too many operands for the operators "
- "in the given expression");
-
+ error(EXIT_FAILURE, 0, "too many operands");
/* If the remaining operand is an array then save the array as a
FITS image, if not, simply print the floating point number. */
- if(p->operands->array)
+ if(p->operands->data)
{
/* Internally, all arrays were double type. But the user could set
the the output type using the type option. So if the user has
asked for anything other than a double, we will have to convert
the arrays.*/
- if(p->type==DOUBLE_IMG)
- gal_fits_array_to_file(p->cp.output, "astimgarith",
- DOUBLE_IMG, p->operands->array,
- p->s0, p->s1, p->anyblank,
- p->wcs, NULL, SPACK_STRING);
+ if(p->outtype==GAL_DATA_TYPE_DOUBLE)
+ ddata=p->operands->data;
else
{
- gal_fits_change_type(p->operands->array, DOUBLE_IMG, p->s0*p->s1,
- p->anyblank, &array, p->type);
- gal_fits_array_to_file(p->cp.output, "astimgarith", p->type,
- array, p->s0, p->s1, p->anyblank, p->wcs,
- NULL, SPACK_STRING);
- free(array);
+ ddata=gal_data_copy_to_new_type(p->operands->data, p->outtype);
+ gal_data_free(p->operands->data);
}
- free(p->operands->array);
+
+ /* Copy the WCS structure into the final dataset and write it into
+ the output file. */
+ ddata->wcs=p->refdata.wcs;
+ gal_fits_write_img(ddata, p->cp.output, "Arithmetic", NULL,
+ SPACK_STRING);
+
+ /* Clean up: */
+ gal_data_free(ddata);
+ free(p->refdata.dsize);
+ wcsfree(p->refdata.wcs);
}
else
printf("%g\n", p->operands->number);
-
- /* If there are any remaining HDUs in the hdus linked list, then
- free them. */
- while(p->hdus!=NULL)
- {
- gal_linkedlist_pop_from_stll(&p->hdus, &tokeepvalue);
- free(tokeepvalue);
- }
+ /* Clean up. Note that the tokens were taken from the command-line
+ arguments, so the string within each token linked list must not be
+ freed. */
+ gal_linkedlist_free_stll(p->tokens, 0);
+ free(p->operands);
}
@@ -1160,5 +185,6 @@ reversepolish(struct imgarithparams *p)
void
imgarith(struct imgarithparams *p)
{
+ /* Parse the arguments */
reversepolish(p);
}
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 64f44a0..2b72b44 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
diff --git a/bin/arithmetic/main.c b/bin/arithmetic/main.c
index ecb4618..ad723a2 100644
--- a/bin/arithmetic/main.c
+++ b/bin/arithmetic/main.c
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -48,6 +48,9 @@ main (int argc, char *argv[])
/* Run MakeProfiles */
imgarith(&p);
+ /* Free any allocated space */
+ freeandreport(&p, &t1);
+
/* Return successfully.*/
return 0;
}
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index f33a2b0..3aa5000 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -41,7 +41,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Constants: */
#define NEGDASHREPLACE 11 /* A vertical tab (ASCII=11) for negative dash */
-#define NOOPTARRAY NULL
+#define NOOPTDATA NULL
#define NOOPTNUMBER NAN
#define NOOPTFILENAME ""
@@ -50,20 +50,19 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
1. A file name which is not read yet. When inactive, NOOPTFILENAME.
2. A number which is not used yet. When inactive, NOOPTNUMBER.
- 3. An array (operator output). When inactive, NOOPTARRAY.
+ 3. A dataset (operator output). When inactive, NOOPTDATA.
- In every node of the operand linked list, only one of these should
- be active. The other two should be inactive. Otherwise it will be a
- bug and will cause problems. All the operands operate on this
- premise.
+ In every node of the operand linked list, only one of these should be
+ active. The other two should be inactive. Otherwise it will be a bug and
+ will cause problems. All the operands operate on this premise.
*/
struct operand
{
- char *filename;
- char *hdu;
- double number;
- double *array;
- struct operand *next;
+ char *filename; /* !=NULL if the operand is a filename. */
+ char *hdu; /* !=NULL if the operand is a filename. */
+ double number; /* Value used if not a dataset. */
+ gal_data_t *data; /* !=NULL if the operand is a dataset. */
+ struct operand *next; /* Pointer to next operand. */
};
@@ -88,19 +87,13 @@ struct imgarithparams
/* Input: */
struct gal_linkedlist_stll *hdus; /* List of all given HDU strings. */
struct gal_linkedlist_stll *tokens; /* List of all arithmetic tokens. */
- float *array; /* Main array to keep results. */
- float *tmp; /* Secondary array for temporary reading. */
size_t numfits; /* Total number of input FITS images. */
size_t addcounter; /* The number of FITS images added. */
size_t popcounter; /* The number of FITS images popped. */
- size_t s0; /* Length of image along first C axis. */
- size_t s1; /* Length of image along second C axis. */
- int nwcs; /* The number of WCS coordinates. */
- struct wcsprm *wcs; /* The WCS structure. */
- int anyblank; /* If there are blank pixels in the image. */
+ gal_data_t refdata; /* Container for information of the data. */
/* Output: */
- int type; /* User's desired output bixpix value. */
+ int outtype; /* User's desired output bixpix value. */
/* Operating mode: */
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
new file mode 100644
index 0000000..369917c
--- /dev/null
+++ b/bin/arithmetic/operands.c
@@ -0,0 +1,177 @@
+/*********************************************************************
+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/fits.h>
+#include <gnuastro/array.h>
+
+#include "main.h"
+
+#include "operands.h"
+
+
+
+
+
+
+size_t
+num_operands(struct imgarithparams *p)
+{
+ size_t counter=0;
+ struct operand *tmp=NULL;
+ for(tmp=p->operands;tmp!=NULL;tmp=tmp->next)
+ ++counter;
+ return counter;
+}
+
+
+
+
+
+void
+add_operand(struct imgarithparams *p, char *filename, double number,
+ 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->number=number;
+ 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;
+ }
+
+ /* Make the link to the previous list. */
+ newnode->next=p->operands;
+ p->operands=newnode;
+}
+
+
+
+
+
+void
+pop_operand(struct imgarithparams *p, double *number, gal_data_t **data,
+ char *operator)
+{
+ size_t i;
+ struct uiparams *up=&p->up;
+ struct operand *operands=p->operands;
+ char *maskname, *mhdu, *filename, *hdu;
+
+ /* If the operand linked list has finished, then give an error and
+ exit. */
+ if(operands==NULL)
+ error(EXIT_FAILURE, 0, "not enough operands for the \"%s\" operator",
+ operator);
+
+
+ /* Set the dataset. If filename is present then read the file
+ and fill in the array, if not then just set the array. */
+ if(strlen(operands->filename))
+ {
+ hdu=operands->hdu;
+ filename=operands->filename;
+
+ /* In case this is the first image that is read, then read the
+ WCS information and set the mask name so masked pixels can be
+ set to NaN. For the other images, the mask can be completely
+ ignored. */
+ mhdu = p->popcounter ? NULL : up->mhdu;
+ maskname = p->popcounter ? NULL : up->maskname;
+ if(p->popcounter==0)
+ gal_fits_read_wcs(filename, hdu, 0, 0, &p->refdata.nwcs,
+ &p->refdata.wcs);
+
+ /* Read the input image as a double type array. */
+ *data=gal_fits_read_to_type(filename, maskname, hdu, mhdu,
+ GAL_DATA_TYPE_DOUBLE);
+
+ /* When the reference data structure's dimensionality is non-zero, it
+ means that this is not the first image read. So, write its basic
+ information into the reference data structure for future
+ checks. */
+ if(p->refdata.ndim)
+ {
+ if(gal_data_dsize_is_different(&p->refdata, *data))
+ error(EXIT_FAILURE, 0, "%s (hdu=%s): has a different size "
+ "compared to previous images. All the images must be "
+ "the same size in order for Arithmetic to work",
+ filename, hdu);
+ }
+ else
+ {
+ /* Set the dimensionality. */
+ p->refdata.ndim=(*data)->ndim;
+
+ /* Allocate the dsize array. */
+ errno=0;
+ p->refdata.dsize=malloc(p->refdata.ndim * sizeof *p->refdata.dsize);
+ if(p->refdata.dsize==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for p->refdata.dsize in "
+ "`operands.c'", p->refdata.ndim * sizeof *p->refdata.dsize);
+
+ /* Write the values into it. */
+ for(i=0;i<p->refdata.ndim;++i)
+ p->refdata.dsize[i]=(*data)->dsize[i];
+ }
+
+ /* Free the HDU string: */
+ free(hdu);
+
+ /* Add to the number of popped FITS images: */
+ ++p->popcounter;
+
+ /* Report the read image if desired: */
+ if(p->cp.verb) printf("%s is read.\n", filename);
+ }
+ else
+ *data=operands->data;
+
+ /* Set the number: */
+ *number=operands->number;
+
+ /* Remove this node from the queue. */
+ p->operands=operands->next;
+ free(operands);
+}
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/operands.h
similarity index 68%
copy from bin/arithmetic/ui.h
copy to bin/arithmetic/operands.h
index 7bf7564..4f3798d 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/operands.h
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -20,10 +20,19 @@ 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 IMCROPUI_H
-#define IMCROPUI_H
+#ifndef OPERANDS_H
+#define OPERANDS_H
+
+size_t
+num_operands(struct imgarithparams *p);
+
+void
+add_operand(struct imgarithparams *p, char *filename, double number,
+ gal_data_t *data);
void
-setparams(int argc, char *argv[], struct imgarithparams *p);
+pop_operand(struct imgarithparams *p, double *number, gal_data_t **data,
+ char *operator);
+
#endif
diff --git a/bin/arithmetic/operators.c b/bin/arithmetic/operators.c
new file mode 100644
index 0000000..dfd8d4f
--- /dev/null
+++ b/bin/arithmetic/operators.c
@@ -0,0 +1,826 @@
+/*********************************************************************
+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"
+#include "operators.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
new file mode 100644
index 0000000..6ce946c
--- /dev/null
+++ b/bin/arithmetic/operators.h
@@ -0,0 +1,96 @@
+/*********************************************************************
+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/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index 8caf5cf..556b910 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -125,7 +125,7 @@ readconfig(char *filename, struct imgarithparams *p)
else if(strcmp(name, "type")==0)
{
if(p->up.typeset) continue;
- gal_checkset_known_types(value, &p->type, filename, lineno);
+ gal_checkset_known_types(value, &p->outtype, filename, lineno);
p->up.typeset=1;
}
@@ -174,7 +174,7 @@ printvalues(FILE *fp, struct imgarithparams *p)
if(cp->outputset)
fprintf(fp, CONF_SHOWFMT"%s\n", "output", cp->output);
if(up->typeset)
- gal_configfiles_print_type(fp, p->type);
+ gal_configfiles_print_type(fp, p->outtype);
/* For the operating mode, first put the macro to print the common
@@ -394,3 +394,40 @@ setparams(int argc, char *argv[], struct imgarithparams *p)
/* Do a sanity check. */
sanitycheck(p);
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**************************************************************/
+/************ Free allocated, report *************/
+/**************************************************************/
+void
+freeandreport(struct imgarithparams *p, struct timeval *t1)
+{
+ free(p->cp.output);
+
+ /* If there are any remaining HDUs in the hdus linked list, then
+ free them. */
+ if(p->hdus)
+ gal_linkedlist_free_stll(p->hdus, 1);
+
+ /* Report the duration of the job */
+ if(p->cp.verb)
+ gal_timing_report(t1, SPACK_NAME" finished in", 0);
+}
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index 7bf7564..f327640 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -1,6 +1,6 @@
/*********************************************************************
-ImageArithmetic - Do arithmetic operations on images.
-ImageArithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
+Arithmetic - Do arithmetic operations on images.
+Arithmetic is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
Mohammad Akhlaghi <address@hidden>
@@ -26,4 +26,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
void
setparams(int argc, char *argv[], struct imgarithparams *p);
+void
+freeandreport(struct imgarithparams *p, struct timeval *t1);
+
#endif
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index 9cb3bc4..3c7e69b 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -151,7 +151,7 @@ static struct argp_option options[] =
'T',
"STR",
0,
- "byte, short, long, longlong, float, double.",
+ "uchar, short, long, longlong, float, double.",
2
},
diff --git a/configure.ac b/configure.ac
index 9b6416d..a1bbc2e 100644
--- a/configure.ac
+++ b/configure.ac
@@ -152,9 +152,12 @@ AC_SEARCH_LIBS([wcspih], [wcs], [],
# These are secondary tests for more fine-grained control in libraries that
# have already been checked. We don't need to add them to the LIBS
# variable, so we are using AC_CHECK_LIB for these tests.
-AC_CHECK_LIB([wcs], [wcslib_version],
- [AC_DEFINE([HAVE_WCSLIB_VERSION], [1], [Has wcslib_version])],
- [], [-lcfitsio -lm])
+AC_CHECK_LIB([wcs], [wcslib_version], [has_wcslib_version=1],
+ [has_wcslib_version=0], [-lcfitsio -lm])
+AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_WCSLIB_VERSION], [$has_wcslib_version],
+ [WCSLIB comes with wcslib_version])
+AC_SUBST(HAVE_WCSLIB_VERSION, [$has_wcslib_version])
+
AC_CHECK_LIB([pthread], [pthread_barrier_wait], [has_pthread_barrier=1],
[has_pthread_barrier=0])
AC_DEFINE_UNQUOTED([GAL_CONFIG_HAVE_PTHREAD_BARRIER], [$has_pthread_barrier],
diff --git a/doc/Makefile.am b/doc/Makefile.am
index 7b29d2a..b314c29 100644
--- a/doc/Makefile.am
+++ b/doc/Makefile.am
@@ -76,11 +76,61 @@ dist_infognuastro_DATA =
$(top_srcdir)/doc/gnuastro-figures/*
-## Man pages: to keep things managable, sort the utilities alphabetically.
-dist_man_MANS = man/astarithmetic.1 man/astconvertt.1 man/astconvolve.1 \
-man/astcosmiccal.1 man/astheader.1 man/astimgcrop.1 man/astimgstat.1 \
-man/astimgwarp.1 man/astmkcatalog.1 man/astmknoise.1 man/astmkprof.1 \
-man/astnoisechisel.1 man/astsubtractsky.1 man/asttable.1
+## Man pages to build
+## ==================
+##
+## Installing all the programs is not mandatory in Gnuastro. If a program
+## is not built, we shouldn't build its man-page either.
+if COND_ARITHMETIC
+ MAYBE_ARITHMETIC_MAN = man/astarithmetic.1
+endif
+if COND_CONVERTT
+ MAYBE_CONVERTT_MAN = man/astconvertt.1
+endif
+if COND_CONVOLVE
+ MAYBE_CONVOLVE_MAN = man/astconvolve.1
+endif
+if COND_COSMICCAL
+ MAYBE_COSMICCAL_MAN = man/astcosmiccal.1
+endif
+if COND_HEADER
+ MAYBE_HEADER_MAN = man/astheader.1
+endif
+if COND_IMGCROP
+ MAYBE_IMGCROP_MAN = man/astimgcrop.1
+endif
+if COND_IMGSTAT
+ MAYBE_IMGSTAT_MAN = man/astimgstat.1
+endif
+if COND_IMGWARP
+ MAYBE_IMGWARP_MAN = man/astimgwarp.1
+endif
+if COND_MKCATALOG
+ MAYBE_MKCATALOG_MAN = man/astmkcatalog.1
+endif
+if COND_MKNOISE
+ MAYBE_MKNOISE_MAN = man/astmknoise.1
+endif
+if COND_MKPROF
+ MAYBE_MKPROF_MAN = man/astmkprof.1
+endif
+if COND_NOISECHISEL
+ MAYBE_NOISECHISEL_MAN = man/astnoisechisel.1
+endif
+if COND_SUBTRACTSKY
+ MAYBE_SUBTRACTSKY_MAN = man/astsubtractsky.1
+endif
+if COND_TABLE
+ MAYBE_TABLE_MAN = man/asttable.1
+endif
+#if COND_TEMPLATE
+# MAYBE_TEMPLATE_MAN = man/astTEMPLATE.1
+#endif
+dist_man_MANS = $(MAYBE_ARITHMETIC_MAN) $(MAYBE_CONVERTT_MAN) \
+ $(MAYBE_CONVOLVE_MAN) $(MAYBE_COSMICCAL_MAN) $(MAYBE_HEADER_MAN) \
+ $(MAYBE_IMGCROP_MAN) $(MAYBE_IMGSTAT_MAN) $(MAYBE_IMGWARP_MAN) \
+ $(MAYBE_MKCATALOG_MAN) $(MAYBE_MKNOISE_MAN) $(MAYBE_MKPROF_MAN) \
+ $(MAYBE_NOISECHISEL_MAN) $(MAYBE_SUBTRACTSKY_MAN) $(MAYBE_TABLE_MAN)
## See if help2man is present or not. When help2man doesn't exist, we don't
diff --git a/lib/Makefile.am b/lib/Makefile.am
index 136a449..e1704f5 100644
--- a/lib/Makefile.am
+++ b/lib/Makefile.am
@@ -31,9 +31,10 @@ libgnuastro_la_LDFLAGS = -version-info $(GAL_LT_VERSION)
# Specify the library .c files
-libgnuastro_la_SOURCES = array.c box.c checkset.c configfiles.c data.c \
- fits.c git.c linkedlist.c mesh.c mode.c polygon.c qsort.c \
- spatialconvolve.c statistics.c threads.c timing.c txtarray.c wcs.c
+libgnuastro_la_SOURCES = array.c box.c changetype.c checkset.c \
+ configfiles.c data.c fits.c git.c linkedlist.c mesh.c mode.c polygon.c \
+ qsort.c spatialconvolve.c statistics.c threads.c timing.c txtarray.c \
+ wcs.c
# Installed headers, note that we are not blindly including all `.h' files
@@ -53,9 +54,9 @@ pkginclude_HEADERS = gnuastro/config.h $(headersdir)/array.h
\
# and if they are not explicitly mentioned somewhere in the Makefile, they
# will not distributed, so we need to explicitly tell Automake to
# distribute them here.
-EXTRA_DIST = gnuastro.pc.in config.h.in checkset.h commonargs.h \
- commonparams.h configfiles.h fixedstringmacros.h mode.h neighbors.h \
- timing.h $(headersdir)/README
+EXTRA_DIST = gnuastro.pc.in changetype.h config.h.in checkset.h \
+ commonargs.h commonparams.h configfiles.h fixedstringmacros.h mode.h \
+ neighbors.h timing.h $(headersdir)/README
# Definitions for Gnuastro's the pkg-config file (inspired from GSL's
@@ -72,6 +73,7 @@ gnuastro/config.h: Makefile config.h.in
$(MKDIR_P) gnuastro
$(SED) -e 's|@address@hidden|$(VERSION)|g' \
-e 's|@address@hidden|$(HAVE_LIBGIT2)|g' \
+ -e 's|@address@hidden|$(HAVE_WCSLIB_VERSION)|g' \
-e 's|@address@hidden|$(HAVE_PTHREAD_BARRIER)|g' \
$(srcdir)/config.h.in >> address@hidden
chmod a-w address@hidden
diff --git a/lib/changetype.c b/lib/changetype.c
new file mode 100644
index 0000000..21f3220
--- /dev/null
+++ b/lib/changetype.c
@@ -0,0 +1,363 @@
+/*********************************************************************
+changetype -- Changing the array of one data type to another.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2016, 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 <errno.h>
+#include <error.h>
+
+#include <gnuastro/data.h>
+
+#include <changetype.h>
+
+
+
+
+
+
+
+
+
+
+#define BASIC_DEFINITIONS \
+ unsigned char *iuc=in->array, *iiuc=in->array; \
+ char *ic=in->array, *iic=in->array; \
+ unsigned short *ius=in->array, *iius=in->array; \
+ short *is=in->array, *iis=in->array; \
+ unsigned int *iui=in->array, *iiui=in->array; \
+ int *ii=in->array, *iii=in->array; \
+ unsigned long *iul=in->array, *iiul=in->array; \
+ long *il=in->array, *iil=in->array; \
+ LONGLONG *iL=in->array, *iiL=in->array; \
+ float *iif=in->array; \
+ double *id=in->array;
+
+
+
+
+
+#define IN_INTEGER \
+ case GAL_DATA_TYPE_UCHAR: \
+ do *o=*iuc++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iiuc++ == GAL_DATA_BLANK_UCHAR ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_CHAR: \
+ do *o=*ic++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iic++ == GAL_DATA_BLANK_CHAR ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_USHORT: \
+ do *o=*ius++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iius++ == GAL_DATA_BLANK_USHORT ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_SHORT: \
+ do *o=*is++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iis++ == GAL_DATA_BLANK_SHORT ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_UINT: \
+ do *o=*iui++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iiui++ == GAL_DATA_BLANK_UINT ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_INT: \
+ do *o=*ii++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iii++ == GAL_DATA_BLANK_INT ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_ULONG: \
+ do *o=*iul++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iiul++ == GAL_DATA_BLANK_ULONG ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_LONG: \
+ do *o=*il++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iil++ == GAL_DATA_BLANK_LONG ? blank : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_LONGLONG: \
+ do *o=*iL++; while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = *iiL++==GAL_DATA_BLANK_LONGLONG ? blank : *ob; \
+ while(++ob<of); } \
+ return;
+
+
+
+
+
+#define OUT_INT_IN_FLOAT \
+ case GAL_DATA_TYPE_FLOAT: \
+ do *o=roundf(*iif++); while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = isnan(*iiif++) ? GAL_DATA_BLANK_UCHAR : *ob; \
+ while(++ob<of); } \
+ return; \
+ \
+ case GAL_DATA_TYPE_DOUBLE: \
+ do *o=round(*id++); while(++o<of); \
+ if(in->anyblank && in->type!=out->type) \
+ { do *ob = isnan(*iid++) ? GAL_DATA_BLANK_UCHAR : *ob; \
+ while(++ob<of); } \
+ return;
+
+
+
+
+
+#define OUT_FLOAT_IN_FLOAT \
+ case GAL_DATA_TYPE_FLOAT: \
+ do *o=*iif++; while(++o<of); \
+ return; \
+ \
+ case GAL_DATA_TYPE_DOUBLE: \
+ do *o=*id++; while(++o<of); \
+ return;
+
+
+
+
+
+#define DEFAULT_PRINT_ERROR \
+ case GAL_DATA_TYPE_STRING: \
+ error(EXIT_FAILURE, 0, "type conversion can't be done on " \
+ "string arrays."); \
+ \
+ default: \
+ error(EXIT_FAILURE, 0, "type %d was not recognized in " \
+ "`change_type_out_is_%s' (data.c)", in->type, tforerr);
+
+
+
+
+
+#define CHANGE_OUT_IS_INTEGER \
+ BASIC_DEFINITIONS; \
+ float *iiif=in->array; \
+ double *iid=in->array; \
+ \
+ switch(in->type) \
+ { \
+ IN_INTEGER; \
+ \
+ OUT_INT_IN_FLOAT; \
+ \
+ DEFAULT_PRINT_ERROR; \
+ }
+
+
+
+
+
+#define CHANGE_OUT_IS_FLOATING \
+ BASIC_DEFINITIONS; \
+ \
+ switch(in->type) \
+ { \
+ IN_INTEGER; \
+ \
+ OUT_FLOAT_IN_FLOAT; \
+ \
+ DEFAULT_PRINT_ERROR; \
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+void
+gal_changetype_out_is_uchar(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="uchar";
+ unsigned char *o = out->array, *ob = out->array;
+ unsigned char blank=GAL_DATA_BLANK_UCHAR, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_char(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="char";
+ char *o = out->array, *ob = out->array;
+ char blank=GAL_DATA_BLANK_CHAR, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+void
+gal_changetype_out_is_ushort(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="ushort";
+ unsigned short *o = out->array, *ob = out->array;
+ unsigned short blank=GAL_DATA_BLANK_USHORT, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_short(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="short";
+ short *o = out->array, *ob = out->array;
+ short blank=GAL_DATA_BLANK_SHORT, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+void
+gal_changetype_out_is_uint(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="uint";
+ unsigned int *o = out->array, *ob = out->array;
+ unsigned int blank=GAL_DATA_BLANK_UINT, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_int(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="int";
+ int *o = out->array, *ob = out->array;
+ int blank=GAL_DATA_BLANK_INT, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_ulong(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="ulong";
+ unsigned long *o = out->array, *ob = out->array;
+ unsigned long blank=GAL_DATA_BLANK_LONG, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_long(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="long";
+ long *o = out->array, *ob = out->array;
+ long blank=GAL_DATA_BLANK_LONG, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_longlong(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="longlong";
+ LONGLONG *o = out->array, *ob = out->array;
+ LONGLONG blank=GAL_DATA_BLANK_LONGLONG, *of = o + out->size;
+
+ CHANGE_OUT_IS_INTEGER;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_float(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="float";
+ float *o = out->array, *ob = out->array;
+ float blank=GAL_DATA_BLANK_FLOAT, *of = o + out->size;
+
+ CHANGE_OUT_IS_FLOATING;
+}
+
+
+
+
+
+void
+gal_changetype_out_is_double(gal_data_t *in, gal_data_t *out)
+{
+ char *tforerr="double";
+ double *o = out->array, *ob = out->array;
+ double blank=GAL_DATA_BLANK_DOUBLE, *of = o + out->size;
+
+ CHANGE_OUT_IS_FLOATING;
+}
diff --git a/lib/changetype.h b/lib/changetype.h
new file mode 100644
index 0000000..8385fed
--- /dev/null
+++ b/lib/changetype.h
@@ -0,0 +1,60 @@
+/*********************************************************************
+changetype -- Changing the array of one data type to another.
+This is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2016, 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 __GAL_CHANGETYPE_H__
+#define __GAL_CHANGETYPE_H__
+
+void
+gal_changetype_out_is_uchar(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_char(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_ushort(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_short(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_uint(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_int(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_ulong(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_long(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_longlong(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_float(gal_data_t *in, gal_data_t *out);
+
+void
+gal_changetype_out_is_double(gal_data_t *in, gal_data_t *out);
+
+
+#endif
diff --git a/lib/checkset.c b/lib/checkset.c
index cacacee..de075e8 100644
--- a/lib/checkset.c
+++ b/lib/checkset.c
@@ -25,7 +25,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <errno.h>
#include <error.h>
-#include <dirent.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
@@ -33,6 +32,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <fitsio.h>
+#include <gnuastro/data.h>
+
#include "checkset.h"
@@ -508,27 +509,29 @@ gal_checkset_any_double(char *optarg, double *var, char
*lo, char so,
/* Check if the value to the `--type' option is recognized, if so set the
integer value. */
void
-gal_checkset_known_types(char *optarg, int *bitpix, char *filename,
+gal_checkset_known_types(char *optarg, int *type, char *filename,
size_t lineno)
{
/* First check if the value is one of the accepted types. */
- if (strcmp(optarg, "byte")==0) *bitpix=BYTE_IMG;
- else if(strcmp(optarg, "short")==0) *bitpix=SHORT_IMG;
- else if(strcmp(optarg, "long")==0) *bitpix=LONG_IMG;
- else if(strcmp(optarg, "longlong")==0) *bitpix=LONGLONG_IMG;
- else if(strcmp(optarg, "float")==0) *bitpix=FLOAT_IMG;
- else if(strcmp(optarg, "double")==0) *bitpix=DOUBLE_IMG;
+ if (strcmp(optarg, "uchar")==0) *type = GAL_DATA_TYPE_UCHAR;
+ else if(strcmp(optarg, "short")==0) *type = GAL_DATA_TYPE_SHORT;
+ else if(strcmp(optarg, "long")==0) *type = GAL_DATA_TYPE_LONG;
+ else if(strcmp(optarg, "longlong")==0) *type = GAL_DATA_TYPE_LONGLONG;
+ else if(strcmp(optarg, "float")==0) *type = GAL_DATA_TYPE_FLOAT;
+ else if(strcmp(optarg, "double")==0) *type = GAL_DATA_TYPE_DOUBLE;
else
{
if(filename)
error_at_line(EXIT_FAILURE, 0, filename, lineno, "given value of "
"the `type' option (`%s') is not recognized. It must "
- "be `byte', `short', `long', `longlong', `float', or "
- "`double'.", optarg);
+ "be `uchar', `short', `long', `longlong', `float', or "
+ "`double'. The FITS standard only defines these types "
+ "for image arrays", optarg);
else
error(EXIT_FAILURE, 0, "given value of the `--type' (`-T') option "
"(`%s') is not recognized. It must be `byte', `short', `long' "
- "`longlong', `float', or `double'.", optarg);
+ "`longlong', `float', or `double'. The FITS standard only "
+ "defines these types for image arrays", optarg);
}
}
@@ -996,3 +999,21 @@ gal_checkset_check_dir_write_add_slash(char **dirname)
free(*dirname);
*dirname=tmpname;
}
+
+
+
+
+
+/* If the given directory exists, then nothing is done, if it doesn't, it
+ will be created. */
+void
+gal_checkset_mkdir(char *dirname)
+{
+ struct stat st={0};
+ if( stat(dirname, &st) == -1 )
+ {
+ errno=0;
+ if( mkdir(dirname, 0700) == -1 )
+ error(EXIT_FAILURE, errno, "making %s", dirname);
+ }
+}
diff --git a/lib/checkset.h b/lib/checkset.h
index 7954a7e..3930e16 100644
--- a/lib/checkset.h
+++ b/lib/checkset.h
@@ -181,7 +181,7 @@ gal_checkset_any_double(char *optarg, double *var, char
*lo, char so,
/********** Check fixed strings ************/
/**************************************************************/
void
-gal_checkset_known_types(char *optarg, int *bitpix, char *filename,
+gal_checkset_known_types(char *optarg, int *type, char *filename,
size_t lineno);
@@ -233,6 +233,9 @@ gal_checkset_not_dir_part(char *input);
void
gal_checkset_check_dir_write_add_slash(char **dirname);
+void
+gal_checkset_mkdir(char *dirname);
+
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/config.h.in b/lib/config.h.in
index f378244..e742ee6 100644
--- a/lib/config.h.in
+++ b/lib/config.h.in
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Macros. */
#define GAL_CONFIG_VERSION "@VERSION@"
#define GAL_CONFIG_HAVE_LIBGIT2 @HAVE_LIBGIT2@
+#define GAL_CONFIG_HAVE_WCSLIB_VERSION @HAVE_WCSLIB_VERSION@
#define GAL_CONFIG_HAVE_PTHREAD_BARRIER @HAVE_PTHREAD_BARRIER@
diff --git a/lib/data.c b/lib/data.c
index 8d3291a..68a1230 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -24,12 +24,17 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <errno.h>
#include <error.h>
+#include <fcntl.h>
+#include <unistd.h>
#include <stdlib.h>
#include <string.h>
+#include <dirent.h>
+#include <sys/mman.h>
#include <gnuastro/data.h>
-
+#include <checkset.h>
+#include <changetype.h>
@@ -39,8 +44,30 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/*********************************************************************/
-/************* Type size and allocation *******************/
+/************* Size and allocation *******************/
/*********************************************************************/
+int
+gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second)
+{
+ size_t i;
+
+ /* First make sure that the dimensionality is the same. */
+ if(first->ndim!=second->ndim)
+ return 1;
+
+ /* Check if the sizes along all dimensions are the same: */
+ for(i=0;i<first->ndim;++i)
+ if( first->dsize[i] != second->dsize[i] )
+ return 1;
+
+ /* If it got to here, we know the dimensions have the same length. */
+ return 0;
+}
+
+
+
+
+
size_t
gal_data_sizeof(int type)
{
@@ -48,8 +75,8 @@ gal_data_sizeof(int type)
switch(type)
{
case GAL_DATA_TYPE_BIT:
- error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support TBIT "
- "datatype, please get in touch with us to implement it.");
+ error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
+ "types, please get in touch with us to implement it.");
/* The parenthesis after sizeof is not a function, it is actually a
type cast, so we have put a space between size of and the
@@ -83,7 +110,7 @@ gal_data_sizeof(int type)
return sizeof (long);
case GAL_DATA_TYPE_LONGLONG:
- return sizeof (long long);
+ return sizeof (LONGLONG);
case GAL_DATA_TYPE_FLOAT:
if( sizeof (float) != 4 )
@@ -111,8 +138,8 @@ gal_data_sizeof(int type)
}
error(EXIT_FAILURE, 0, "Control has reached the end of `gal_data_sizeof' "
- "This is a bug! Please contact us at %s so we can find the cause of "
- "the problem.", PACKAGE_BUGREPORT);
+ "This is a bug! Please contact us at %s so we can find the cause "
+ "of the problem.", PACKAGE_BUGREPORT);
return -1;
}
@@ -125,14 +152,32 @@ gal_data_sizeof(int type)
bytes each element needs will be determined internaly by this function
using the datatype argument, so you don't have to worry about it. */
void *
-gal_data_alloc(int type, size_t size)
+gal_data_malloc_array(int type, size_t size)
{
void *array;
errno=0;
array=malloc( size * gal_data_sizeof(type) );
if(array==NULL)
- error(EXIT_FAILURE, errno, "array of %zu bytes in gal_data_alloc",
+ error(EXIT_FAILURE, errno, "array of %zu bytes in gal_data_malloc_array",
+ size * gal_data_sizeof(type));
+
+ return array;
+}
+
+
+
+
+
+void *
+gal_data_calloc_array(int type, size_t size)
+{
+ void *array;
+
+ errno=0;
+ array=calloc( size, gal_data_sizeof(type) );
+ if(array==NULL)
+ error(EXIT_FAILURE, errno, "array of %zu bytes in gal_data_calloc_array",
size * gal_data_sizeof(type));
return array;
@@ -142,6 +187,178 @@ gal_data_alloc(int type, size_t size)
+void
+gal_data_mmap(gal_data_t *data)
+{
+ int filedes;
+ char *filename;
+ unsigned char uc=0;
+ size_t bsize=data->size*gal_data_sizeof(data->type);
+
+ /* Check if the .gnuastro folder exists, write the file there. If it
+ doesn't exist, then make the .gnuastro directory.*/
+ gal_checkset_mkdir(".gnuastro");
+
+ /* Set the filename */
+ gal_checkset_allocate_copy("./.gnuastro/mmap_XXXXXX", &filename);
+
+ /* Create a zero-sized file and keep its descriptor. */
+ errno=0;
+ /*filedes=open(filename, O_RDWR | O_CREAT | O_EXCL | O_TRUNC );*/
+ filedes=mkstemp(filename);
+ if(filedes==-1)
+ error(EXIT_FAILURE, errno, "%s couldn't be created", filename);
+
+
+ /* Make enough space to keep the array data. */
+ errno=0;
+ if( lseek(filedes, bsize, SEEK_SET) == -1 )
+ error(EXIT_FAILURE, errno, "%s: unable to change file position by "
+ "%zu bytes", filename, bsize);
+
+
+ /* Write to the newly set file position so the space is allocated. */
+ if( write(filedes, &uc, 1) == -1)
+ error(EXIT_FAILURE, errno, "%s: unable to write one byte at the "
+ "%zu-th position", filename, bsize);
+
+
+ /* Map the memory. */
+ data->array=mmap(NULL, bsize, PROT_READ | PROT_WRITE, MAP_SHARED,
+ filedes, 0);
+
+ /* Close the file. */
+ if( close(filedes) == -1 )
+ error(EXIT_FAILURE, errno, "%s couldn't be closed", filename);
+
+ /* Set the mmaped flag to 1 and keep the filename. */
+ data->mmapped=1;
+ data->mmapname=filename;
+}
+
+
+
+
+
+gal_data_t *
+gal_data_alloc(int type, size_t ndim, long *dsize, int clear, int map)
+{
+ size_t i;
+ gal_data_t *out;
+
+ /* Allocate the space for the actual structure. */
+ errno=0;
+ out=malloc(sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for gal_data_t in gal_data_alloc",
+ sizeof *out);
+
+
+ /* Set the basic information we know so far */
+ out->ndim=ndim;
+ out->type=type;
+
+ /* Initialize the other values */
+ out->anyblank=0;
+ out->wcs=NULL;
+
+
+ /* Allocate space for the dsize array: */
+ errno=0;
+ out->dsize=malloc(ndim*sizeof *out->dsize);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for dsize in gal_data_alloc",
+ ndim*sizeof *out->dsize);
+
+
+ /* Fill in the dsize values: */
+ out->size=1;
+ for(i=0;i<ndim;++i)
+ {
+ /* Do a small sanity check. */
+ if(dsize[i]==0)
+ error(EXIT_FAILURE, 0, "the size of a dimension cannot be zero. "
+ "dsize[%zu] in `gal_data_alloc' has a value of 0", i);
+
+ /* Write this dimension's size, also correct the total number of
+ elements. */
+ out->size *= ( out->dsize[i] = dsize[i] );
+ }
+
+ /* Allocate space for the array, clear it if necessary: */
+ if(map)
+ gal_data_mmap(out);
+ else
+ {
+ /* Allocate the space for the array. */
+ if(clear)
+ out->array = gal_data_calloc_array(out->type, out->size);
+ else
+ out->array = gal_data_malloc_array(out->type, out->size);
+
+ /* Set the values. */
+ out->mmapped=0;
+ out->mmapname=NULL;
+ }
+
+ /* Return the final structure. */
+ return out;
+}
+
+
+
+
+
+void
+gal_data_free(gal_data_t *data)
+{
+ /* If there is a WCS structure, then free it. */
+ if(data->wcs)
+ wcsfree(data->wcs);
+
+ /* Free all the allocated space and finally the data structure itself. */
+ free(data->dsize);
+
+ if(data->mmapped)
+ {
+ /* Delete the file keeping the array. */
+ remove(data->mmapname);
+
+ /* If there is nothing else in the .gnuastro directory, then delete
+ the .gnuastro directory too. */
+
+ /* Free the file name space */
+ free(data->mmapname);
+ }
+ else
+ free(data->array);
+
+ free(data);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
+ ************** Blank data ***************
+ *************************************************************/
+
+/* Allocate space for one blank value of the given type and put the value
+ in it. */
void *
gal_data_alloc_blank(int type)
{
@@ -156,14 +373,14 @@ gal_data_alloc_blank(int type)
int *i;
unsigned long *ul;
long *l;
- long long *L;
+ LONGLONG *L;
float *f;
double *d;
gsl_complex_float *cx;
gsl_complex *dcx;
/* Allocate the space for the blank value: */
- allocated=gal_data_alloc(1, type);
+ allocated=gal_data_malloc_array(type, 1);
/* Put the blank value into it. */
errno=0;
@@ -254,3 +471,397 @@ gal_data_alloc_blank(int type)
return NULL;
}
+
+
+
+
+
+/* Any non-zero pixel in the `mask' array is set as blank in the `in'
+ array. */
+void
+gal_data_apply_mask(gal_data_t *in, gal_data_t *mask)
+{
+ float *mpt, *m, *mf;
+ gal_data_t *converted;
+ int hasblank=0;
+
+ unsigned char *uc = in->array, *ucf = uc + in->size;
+ char *c = in->array, *cf = c + in->size;
+ char **str = in->array, **strf = str + in->size;
+ unsigned short *us = in->array, *usf = us + in->size;
+ short *s = in->array, *sf = s + in->size;
+ unsigned int *ui = in->array, *uif = ui + in->size;
+ int *ii = in->array, *iif = ii + in->size;
+ unsigned long *ul = in->array, *ulf = ul + in->size;
+ long *l = in->array, *lf = l + in->size;
+ LONGLONG *L = in->array, *Lf = L + in->size;
+ float *f = in->array, *ff = f + in->size;
+ double *d = in->array, *df = d + in->size;
+ gsl_complex_float *cx = in->array, *cxf = cx + in->size;
+ gsl_complex *dcx = in->array, *dcxf = dcx + in->size;
+
+
+ /* Make sure the mask and input image have the same dimentionality: */
+ if(in->ndim!=mask->ndim)
+ error(EXIT_FAILURE, 0, "the `in' and `mask' data structures given "
+ "to `gal_data_apply_mask' do not have the same dimensionality: "
+ "%zu and %zu respectively", in->ndim, mask->ndim);
+
+
+ /* Make sure the mask and input image have the same size along each
+ dimension: */
+ if(gal_data_dsize_is_different(in, mask))
+ error(EXIT_FAILURE, 0, "the `in' and `mask' data structures given "
+ "to `gal_data_apply_mask' do not have the same size along each "
+ "dimension");
+
+
+ /* First convert the mask to floating point. Note that although the
+ standard definition of a mask is for it to be an integer, in some
+ situations, the users might want to specify a floating point image as
+ a mask. Such a mask might have values between 0 and 1 (for example
+ they have made a mock profiles and want to mask all pixels covered by
+ a profile). If we simply convert the mask image to an integer, all
+ pixels with values between zero and one will be set to 0. So we need
+ to internally convert the mask image to float to preserve this. */
+ if(mask->type==GAL_DATA_TYPE_FLOAT)
+ converted=mask;
+ else
+ converted=gal_data_copy_to_new_type(mask, GAL_DATA_TYPE_FLOAT);
+ mpt=converted->array;
+
+
+ /* Go over all the pixels and apply the mask. But first check if there
+ actually are any blank values, it might be the case that there isn't
+ any and in that case we can ignore the checks: */
+ mf=(m=mpt)+in->size; do if(*m++ != 0.0f) hasblank=1; while(m<mf);
+ if(hasblank)
+ {
+ in->anyblank=1;
+ switch(in->type)
+ {
+ case GAL_DATA_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support blank "
+ "values for `GAL_DATA_TYPE_BIT', please get in touch with "
+ "us to see how we can implement it.");
+
+ case GAL_DATA_TYPE_UCHAR:
+ do *uc = *mpt++ == 0.0f ? *uc : GAL_DATA_BLANK_UCHAR;
while(++uc<ucf);
+ break;
+
+ /* CFITSIO says "int for keywords, char for table columns". Here we
+ are only assuming table columns. So in practice this also applies
+ to TSBYTE.*/
+ case GAL_DATA_TYPE_CHAR: case GAL_DATA_TYPE_LOGICAL:
+ do *c = *mpt++ == 0.0f ? *c : GAL_DATA_BLANK_CHAR; while(++c<cf);
+ break;
+
+ case GAL_DATA_TYPE_STRING:
+ do *str = *mpt++ == 0.0f ? *str : GAL_DATA_BLANK_STRING;
+ while(++str<strf);
+ break;
+
+ case GAL_DATA_TYPE_USHORT:
+ do *us = *mpt++ == 0.0f ? *us : GAL_DATA_BLANK_USHORT;
+ while(++us<usf);
+ break;
+
+ case GAL_DATA_TYPE_SHORT:
+ do *s = *mpt++ == 0.0f ? *s : GAL_DATA_BLANK_SHORT; while(++s<sf);
+ break;
+
+ case GAL_DATA_TYPE_UINT:
+ do *ui = *mpt++ == 0.0f ? *ui : GAL_DATA_BLANK_UINT; while(++ui<uif);
+ break;
+
+ case GAL_DATA_TYPE_INT:
+ do *ii = *mpt++ == 0.0f ? *ii : GAL_DATA_BLANK_INT; while(++ii<iif);
+ break;
+
+ case GAL_DATA_TYPE_ULONG:
+ do *ul = *mpt++ == 0.0f ? *ul : GAL_DATA_BLANK_ULONG;
while(++ul<ulf);
+ break;
+
+ case GAL_DATA_TYPE_LONG:
+ do *l = *mpt++ == 0.0f ? *l : GAL_DATA_BLANK_LONG; while(++l<lf);
+ break;
+
+ case GAL_DATA_TYPE_LONGLONG:
+ do *L = *mpt++ == 0.0f ? *L : GAL_DATA_BLANK_LONGLONG; while(++L<Lf);
+ break;
+
+ case GAL_DATA_TYPE_FLOAT:
+ do *f = *mpt++ == 0.0f ? *f : GAL_DATA_BLANK_FLOAT; while(++f<ff);
+ break;
+
+ case GAL_DATA_TYPE_DOUBLE:
+ do *d = *mpt++ == 0.0f ? *d : GAL_DATA_BLANK_DOUBLE; while(++d<df);
+ break;
+
+ case GAL_DATA_TYPE_COMPLEX:
+ do
+ if(*mpt++ == 0.0f)
+ GSL_SET_COMPLEX(cx, GAL_DATA_BLANK_FLOAT, GAL_DATA_BLANK_FLOAT);
+ while(++cx<cxf);
+ break;
+
+ case GAL_DATA_TYPE_DCOMPLEX:
+ do
+ if(*mpt++ == 0.0f)
+ GSL_SET_COMPLEX(dcx, GAL_DATA_BLANK_DOUBLE,
+ GAL_DATA_BLANK_DOUBLE);
+ while(++dcx<dcxf);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "type value of %d not recognized in "
+ "`gal_data_alloc_blank'", in->type);
+ }
+ }
+
+ /* Free the converted mask data if it was allocated: */
+ if(converted!=mask)
+ free(converted);
+}
+
+
+
+
+
+/* Change all blank values in `data' to the value pointed to by `value'. */
+void
+gal_data_blank_to_value(gal_data_t *data, void *value)
+{
+ /* 'value' will only be read from one of these based on the
+ datatype. Which the caller assigned. If there is any problem, it is
+ their responsability, not this function's.*/
+ void *A=data->array;
+ size_t S=data->size;
+ unsigned char *uc = A, *ucf = A+S, ucv = *(unsigned char *) value;
+ char *c = A, *cf = A+S, cv = *(char *) value;
+ char **str = A, **strf = A+S, *strv = *(char **) value;
+ unsigned short *us = A, *usf = A+S, usv = *(unsigned short *) value;
+ short *s = A, *sf = A+S, sv = *(short *) value;
+ unsigned int *ui = A, *uif = A+S, uiv = *(unsigned int *) value;
+ int *in = A, *inf = A+S, inv = *(int *) value;
+ unsigned long *ul = A, *ulf = A+S, ulv = *(unsigned long *) value;
+ long *l = A, *lf = A+S, lv = *(int32_t *) value;
+ LONGLONG *L = A, *Lf = A+S, Lv = *(int64_t *) value;
+ float *f = A, *ff = A+S, fv = *(float *) value;
+ double *d = A, *df = A+S, dv = *(double *) value;
+ gsl_complex_float *cx = A, *cxf = A+S, cxv = *(gsl_complex_float *)
value;
+ gsl_complex *dcx = A, *dcxf = A+S, dcxv = *(gsl_complex *) value;
+
+ switch(data->type)
+ {
+ case GAL_DATA_TYPE_BIT:
+ error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support bit "
+ "datatype, please get in touch with us to implement it.");
+
+ case GAL_DATA_TYPE_UCHAR:
+ do if(*uc==GAL_DATA_BLANK_UCHAR) *uc++=ucv; while(uc<ucf);
+ break;
+
+
+ case GAL_DATA_TYPE_CHAR: case GAL_DATA_TYPE_LOGICAL:
+ do if(*c==GAL_DATA_BLANK_CHAR) *c++=cv; while(c<cf);
+ break;
+
+
+ case GAL_DATA_TYPE_STRING:
+ do if(*str==GAL_DATA_BLANK_STRING) *str++=strv; while(str<strf);
+ break;
+
+
+ case GAL_DATA_TYPE_USHORT:
+ do if(*us==GAL_DATA_BLANK_USHORT) *us++=usv; while(us<usf);
+ break;
+
+
+ case GAL_DATA_TYPE_SHORT:
+ do if(*s==GAL_DATA_BLANK_SHORT) *s++=sv; while(s<sf);
+ break;
+
+
+ case GAL_DATA_TYPE_UINT:
+ do if(*ui==GAL_DATA_BLANK_UINT) *ui++=uiv; while(ui<uif);
+ break;
+
+
+ case GAL_DATA_TYPE_INT:
+ do if(*in==GAL_DATA_BLANK_INT) *in++=inv; while(in<inf);
+ break;
+
+
+ case GAL_DATA_TYPE_ULONG:
+ do if(*ul==GAL_DATA_BLANK_ULONG) *ul++=ulv; while(ul<ulf);
+ break;
+
+
+ case GAL_DATA_TYPE_LONG:
+ do if(*l==GAL_DATA_BLANK_LONG) *l++=lv; while(l<lf);
+ break;
+
+
+ case GAL_DATA_TYPE_LONGLONG:
+ do if(*L==GAL_DATA_BLANK_LONGLONG) *L++=Lv; while(L<Lf);
+ break;
+
+
+ /* Note that a NaN value is not equal to another NaN value, so we
+ can't use the easy check for cases were the blank value is
+ NaN. Also note that `isnan' is actually a macro, so it works for
+ both float and double types.*/
+ case GAL_DATA_TYPE_FLOAT:
+ if(isnan(GAL_DATA_BLANK_FLOAT))
+ do if(isnan(*f)) *f++=fv; while(f<ff);
+ else
+ do if(*f==GAL_DATA_BLANK_FLOAT) *f++=fv; while(f<ff);
+ break;
+
+
+ case GAL_DATA_TYPE_DOUBLE:
+ if(isnan(GAL_DATA_BLANK_DOUBLE))
+ do if(isnan(*d)) *d++=dv; while(d<df);
+ else
+ do if(*d==GAL_DATA_BLANK_FLOAT) *d++=dv; while(d<df);
+ break;
+
+
+ case GAL_DATA_TYPE_COMPLEX:
+ if(isnan(GAL_DATA_BLANK_FLOAT))
+ do
+ if(isnan(GSL_COMPLEX_P_REAL(cx))
+ && isnan(GSL_COMPLEX_P_IMAG(cx)) )
+ GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
+ GSL_COMPLEX_P_IMAG(&cxv));
+ while(++cx<cxf);
+ else
+ do
+ if( GSL_COMPLEX_P_REAL(cx) == GAL_DATA_BLANK_FLOAT
+ && GSL_COMPLEX_P_IMAG(cx) == GAL_DATA_BLANK_FLOAT)
+ GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
+ GSL_COMPLEX_P_IMAG(&cxv));
+ while(++cx<cxf);
+ break;
+
+
+ case GAL_DATA_TYPE_DCOMPLEX:
+ if(isnan(GAL_DATA_BLANK_DOUBLE))
+ do
+ if(isnan(GSL_COMPLEX_P_REAL(dcx))
+ && isnan(GSL_COMPLEX_P_IMAG(dcx)) )
+ GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
+ GSL_COMPLEX_P_IMAG(&dcxv));
+ while(++dcx<dcxf);
+ else
+ do
+ if( GSL_COMPLEX_P_REAL(dcx) == GAL_DATA_BLANK_FLOAT
+ && GSL_COMPLEX_P_IMAG(dcx) == GAL_DATA_BLANK_FLOAT)
+ GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
+ GSL_COMPLEX_P_IMAG(&dcxv));
+ while(++dcx<dcxf);
+ break;
+
+
+ default:
+ error(EXIT_FAILURE, 0, "a bug! type value (%d) not recognized "
+ "in `gal_data_blank_to_value'", data->type);
+ }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
+ ************** Change type ***************
+ *************************************************************/
+int
+gal_data_out_type(gal_data_t *first, gal_data_t *second)
+{
+ return first->type > second->type ? first->type : second->type;
+}
+
+
+
+
+
+gal_data_t *
+gal_data_copy_to_new_type(gal_data_t *in, int newtype)
+{
+ gal_data_t *out;
+
+ /* Allocate space for the output type */
+ out=gal_data_alloc(newtype, in->ndim, in->dsize, 0, in->mmapped);
+
+ /* Fill in the output data array while doing the conversion */
+ switch(newtype)
+ {
+ case GAL_DATA_TYPE_UCHAR:
+ gal_changetype_out_is_uchar(in, out);
+ break;
+
+ case GAL_DATA_TYPE_CHAR:
+ gal_changetype_out_is_char(in, out);
+ break;
+
+ case GAL_DATA_TYPE_USHORT:
+ gal_changetype_out_is_ushort(in, out);
+ break;
+
+ case GAL_DATA_TYPE_SHORT:
+ gal_changetype_out_is_short(in, out);
+ break;
+
+ case GAL_DATA_TYPE_UINT:
+ gal_changetype_out_is_uint(in, out);
+ break;
+
+ case GAL_DATA_TYPE_INT:
+ gal_changetype_out_is_int(in, out);
+ break;
+
+ case GAL_DATA_TYPE_ULONG:
+ gal_changetype_out_is_ulong(in, out);
+ break;
+
+ case GAL_DATA_TYPE_LONG:
+ gal_changetype_out_is_long(in, out);
+ break;
+
+ case GAL_DATA_TYPE_LONGLONG:
+ gal_changetype_out_is_longlong(in, out);
+ break;
+
+ case GAL_DATA_TYPE_FLOAT:
+ gal_changetype_out_is_float(in, out);
+ break;
+
+ case GAL_DATA_TYPE_DOUBLE:
+ gal_changetype_out_is_double(in, out);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "type %d not recognized in "
+ "gal_data_copy_to_new_type", newtype);
+ }
+
+ /* Return the created array */
+ return out;
+}
diff --git a/lib/fits.c b/lib/fits.c
index e4b7f31..852876b 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -37,6 +37,14 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "fixedstringmacros.h"
+
+
+
+
+
+
+
+
/*************************************************************
************** Reporting errors: ***************
*************************************************************/
@@ -74,7 +82,7 @@ gal_fits_io_error(int status, char *message)
/*************************************************************
- ************** Acceptable FITS names ***************
+ ************** FITS names ***************
*************************************************************/
/* IMPORTANT NOTE: if other compression suffixes are add to this function,
@@ -121,6 +129,49 @@ gal_fits_suffix_is_fits(char *suffix)
+/* We have the name of the input file. But in most cases, the files
+ that should be used (for example a mask image) are other extensions
+ in the same file. So the user only has to give the HDU. The job of
+ this function is to determine which is the case and set othername
+ to the appropriate value. */
+void
+gal_fits_file_or_ext_name(char *inputname, char *inhdu, int othernameset,
+ char **othername, char *ohdu, int ohduset,
+ char *type)
+{
+ if(othernameset)
+ {
+ /* In some cases, for example a mask image, both the name and
+ HDU are optional. So just to be safe, we will check this all
+ the time. */
+ if(ohduset==0)
+ error(EXIT_FAILURE, 0, "a %s image was specified (%s). However, "
+ "no HDU is given for it. Please add a HDU. If you regularly "
+ "use the same HDU as %s, you may consider adding it to "
+ "the configuration file. For more information, please see "
+ "the `Configuration files' section of the %s manual by "
+ "running ` info gnuastro ' on the command-line", type,
+ *othername, type, PACKAGE_NAME);
+ if(strcmp(inputname, *othername)==0)
+ {
+ if(strcmp(ohdu, inhdu)==0)
+ error(EXIT_FAILURE, 0, "the specified %s name and "
+ "input image name (%s) are the same while the input "
+ "image hdu name and mask hdu are also identical (%s)",
+ type, inputname, inhdu);
+ }
+ }
+ else if(ohduset && strcmp(ohdu, inhdu))
+ *othername=inputname;
+ else
+ *othername=NULL;
+}
+
+
+
+
+
+
@@ -162,8 +213,42 @@ gal_fits_bitpix_to_type(int bitpix)
case DOUBLE_IMG:
return GAL_DATA_TYPE_DOUBLE;
default:
- error(EXIT_FAILURE, 0, "bitpix value of %d not recognized",
- bitpix);
+ error(EXIT_FAILURE, 0, "bitpix value of %d not recognized in "
+ "gal_fits_bitpix_to_type", bitpix);
+ }
+ return 0;
+}
+
+
+
+
+
+int
+gal_fits_type_to_bitpix(int type)
+{
+ switch(type)
+ {
+ case GAL_DATA_TYPE_UCHAR:
+ return BYTE_IMG;
+ case GAL_DATA_TYPE_CHAR:
+ return SBYTE_IMG;
+ case GAL_DATA_TYPE_USHORT:
+ return USHORT_IMG;
+ case GAL_DATA_TYPE_SHORT:
+ return SHORT_IMG;
+ case GAL_DATA_TYPE_ULONG:
+ return ULONG_IMG;
+ case GAL_DATA_TYPE_LONG:
+ return LONG_IMG;
+ case GAL_DATA_TYPE_LONGLONG:
+ return LONGLONG_IMG;
+ case GAL_DATA_TYPE_FLOAT:
+ return FLOAT_IMG;
+ case GAL_DATA_TYPE_DOUBLE:
+ return DOUBLE_IMG;
+ default:
+ error(EXIT_FAILURE, 0, "type value of %d not recognized in "
+ "gal_fits_type_to_bitpix", type);
}
return 0;
}
@@ -352,474 +437,6 @@ gal_fits_datatype_to_type(int datatype)
-void
-gal_fits_img_bitpix_size(fitsfile *fptr, int *bitpix, long *naxes)
-{
- int status=0, maxdim=10, naxis;
-
- if( fits_get_img_param(fptr, maxdim, bitpix, &naxis, naxes, &status) )
- gal_fits_io_error(status, NULL);
-
- if(naxis!=2)
- error(EXIT_FAILURE, 0, "currently only a 2 dimensional image array "
- "is supported. Your array is %d dimension(s). %s", naxis,
- naxis ? "Please contact us to add this feature." : "This might "
- "Be due to the fact that the data in images with multiple "
- "extensions are sometimes put on the second extension. If this "
- "is the case, try changing the hdu (maybe to --hdu=1)");
-}
-
-
-
-
-
-void
-gal_fits_blank_to_value(void *array, int datatype, size_t size, void *value)
-{
- /* 'value' will only be read from one of these based on the
- datatype. Which the caller assigned. If there is any problem, it is
- their responsability, not this function's.*/
- unsigned char *b, *bf, bv = *(uint8_t *) value;
- char *c, *cf, cv = *(char *) value;
- char **str, **strf, *strv = *(char **) value;
- short *s, *sf, sv = *(int16_t *) value;
- long *l, *lf, lv = *(int32_t *) value;
- LONGLONG *L, *Lf, Lv = *(int64_t *) value;
- float *f, *ff, fv = *(float *) value;
- double *d, *df, dv = *(double *) value;
- gsl_complex_float *cx, *cxf, cxv = *(gsl_complex_float *) value;
- gsl_complex *dcx, *dcxf, dcxv = *(gsl_complex *) value;
- int *in, *inf, inv = *(int *) value;
- unsigned int *ui, *uif, uiv = *(unsigned int *) value;
- unsigned short *us, *usf, usv = *(unsigned short *) value;
- unsigned long *ul, *ulf, ulv = *(unsigned long *) value;
-
- switch(datatype)
- {
- case TBIT:
- error(EXIT_FAILURE, 0, "Currently Gnuastro doesn't support TBIT "
- "datatype, please get in touch with us to implement it.");
-
- case TBYTE:
- bf=(b=array)+size;
- do if(*b==GAL_DATA_BLANK_UCHAR) *b++=bv; while(b<bf);
- break;
-
-
- case TLOGICAL: case TSBYTE:
- cf=(c=array)+size;
- do if(*c==GAL_DATA_BLANK_LOGICAL) *c++=cv; while(c<cf);
- break;
-
-
- case TSTRING:
- strf=(str=array)+size;
- do if(*str==GAL_DATA_BLANK_STRING) *str++=strv; while(str<strf);
- break;
-
-
- case TSHORT:
- sf=(s=array)+size;
- do if(*s==GAL_DATA_BLANK_SHORT) *s++=sv; while(s<sf);
- break;
-
-
- case TLONG:
- lf=(l=array)+size;
- do if(*l==GAL_DATA_BLANK_LONG) *l++=lv; while(l<lf);
- break;
-
-
- case TLONGLONG:
- Lf=(L=array)+size;
- do if(*L==GAL_DATA_BLANK_LONGLONG) *L++=Lv; while(L<Lf);
- break;
-
-
- /* Note that a NaN value is not equal to another NaN value, so we
- can't use the easy check for cases were the blank value is
- NaN. Also note that `isnan' is actually a macro, so it works for
- both float and double types.*/
- case TFLOAT:
- ff=(f=array)+size;
- if(isnan(GAL_DATA_BLANK_FLOAT))
- do if(isnan(*f)) *f++=fv; while(f<ff);
- else
- do if(*f==GAL_DATA_BLANK_FLOAT) *f++=fv; while(f<ff);
- break;
-
-
- case TDOUBLE:
- df=(d=array)+size;
- if(isnan(GAL_DATA_BLANK_DOUBLE))
- do if(isnan(*d)) *d++=dv; while(d<df);
- else
- do if(*d==GAL_DATA_BLANK_FLOAT) *d++=dv; while(d<df);
- break;
-
-
- case TCOMPLEX:
- cxf=(cx=array)+size;
- if(isnan(GAL_DATA_BLANK_FLOAT))
- do
- if(isnan(GSL_COMPLEX_P_REAL(cx))
- && isnan(GSL_COMPLEX_P_IMAG(cx)) )
- GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
- GSL_COMPLEX_P_IMAG(&cxv));
- while(++cx<cxf);
- else
- do
- if( GSL_COMPLEX_P_REAL(cx) == GAL_DATA_BLANK_FLOAT
- && GSL_COMPLEX_P_IMAG(cx) == GAL_DATA_BLANK_FLOAT)
- GSL_SET_COMPLEX(cx, GSL_COMPLEX_P_REAL(&cxv),
- GSL_COMPLEX_P_IMAG(&cxv));
- while(++cx<cxf);
- break;
-
-
- case TDBLCOMPLEX:
- dcxf=(dcx=array)+size;
- if(isnan(GAL_DATA_BLANK_DOUBLE))
- do
- if(isnan(GSL_COMPLEX_P_REAL(dcx))
- && isnan(GSL_COMPLEX_P_IMAG(dcx)) )
- GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
- GSL_COMPLEX_P_IMAG(&dcxv));
- while(++dcx<dcxf);
- else
- do
- if( GSL_COMPLEX_P_REAL(dcx) == GAL_DATA_BLANK_FLOAT
- && GSL_COMPLEX_P_IMAG(dcx) == GAL_DATA_BLANK_FLOAT)
- GSL_SET_COMPLEX(dcx, GSL_COMPLEX_P_REAL(&dcxv),
- GSL_COMPLEX_P_IMAG(&dcxv));
- while(++dcx<dcxf);
- break;
-
-
- case TINT:
- inf=(in=array)+size;
- do if(*in==GAL_DATA_BLANK_INT) *in++=inv; while(in<inf);
- break;
-
-
- case TUINT:
- uif=(ui=array)+size;
- do if(*ui==GAL_DATA_BLANK_UINT) *ui++=uiv; while(ui<uif);
- break;
-
-
- case TUSHORT:
- usf=(us=array)+size;
- do if(*us==GAL_DATA_BLANK_USHORT) *us++=usv; while(us<usf);
- break;
-
-
- case TULONG:
- ulf=(ul=array)+size;
- do if(*ul==GAL_DATA_BLANK_ULONG) *ul++=ulv; while(ul<ulf);
- break;
-
- default:
- error(EXIT_FAILURE, 0, "a bug! datatype value of %d not recognized. "
- "This should not happen here (blanktovalue in fitsarrayvv.c). "
- "Please contact us at %s to see how this happened", datatype,
- PACKAGE_BUGREPORT);
- }
-}
-
-
-
-
-
-void
-gal_fits_change_type(void *in, int inbitpix, size_t size, int anyblank,
- void **out, int outbitpix)
-{
- size_t i=0;
- unsigned char *b, *bf, *ib=in, *iib=in;
- short *s, *sf, *is=in, *iis=in;
- long *l, *lf, *il=in, *iil=in;
- LONGLONG *L, *Lf, *iL=in, *iiL=in;
- float *f, *ff, *iif=in, *iiif=in;
- double *d, *df, *id=in, *iid=in;
-
- /* Allocate space for the output and start filling it. */
- *out=gal_data_alloc(gal_fits_bitpix_to_type(outbitpix), size);
- switch(outbitpix)
- {
- case BYTE_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- bf=(b=*out)+size; do *b=*ib++; while(++b<bf); return;
- case SHORT_IMG:
- bf=(b=*out)+size; do *b=*is++; while(++b<bf);
- if(anyblank)
- {b=*out; do {b[i]=(iis[i]==GAL_DATA_BLANK_SHORT)
- ?GAL_DATA_BLANK_UCHAR:b[i];}
- while(++i!=size);}
- return;
- case LONG_IMG:
- bf=(b=*out)+size; do *b=*il++; while(++b<bf);
- if(anyblank)
- {b=*out; do {b[i]=(iil[i]==GAL_DATA_BLANK_LONG)
- ?GAL_DATA_BLANK_UCHAR:b[i];}
- while(++i!=size);}
- return;
- case LONGLONG_IMG:
- bf=(b=*out)+size; do *b=*iL++; while(++b<bf);
- if(anyblank)
- {b=*out; do {b[i]=(iiL[i]==GAL_DATA_BLANK_LONGLONG)
- ?GAL_DATA_BLANK_UCHAR:b[i];}
- while(++i!=size);}
- return;
- case FLOAT_IMG:
- bf=(b=*out)+size; do *b=roundf(*iif++); while(++b<bf);
- if(anyblank)
- {b=*out; do {b[i]=isnan(iiif[i])?GAL_DATA_BLANK_UCHAR:b[i];}
- while(++i!=size);}
- return;
- case DOUBLE_IMG:
- bf=(b=*out)+size; do *b=round(*id++); while(++b<bf);
- if(anyblank)
- {b=*out; do {b[i]=isnan(iid[i])?GAL_DATA_BLANK_UCHAR:b[i];}
- while(++i!=size);}
- return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
- case SHORT_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- sf=(s=*out)+size; do *s=*ib++; while(++s<sf);
- if(anyblank)
- {s=*out; do {s[i]=(iib[i]==GAL_DATA_BLANK_UCHAR)
- ?GAL_DATA_BLANK_SHORT:s[i];}
- while(++i!=size);}
- return;
- case SHORT_IMG:
- sf=(s=*out)+size; do *s=*is++; while(++s<sf); return;
- case LONG_IMG:
- sf=(s=*out)+size; do *s=*il++; while(++s<sf);
- if(anyblank)
- {s=*out; do {s[i]=(iil[i]==GAL_DATA_BLANK_LONG)
- ?GAL_DATA_BLANK_SHORT:s[i];}
- while(++i!=size);}
- return;
- case LONGLONG_IMG:
- sf=(s=*out)+size; do *s=*iL++; while(++s<sf);
- if(anyblank)
- {s=*out; do {s[i]=(iiL[i]==GAL_DATA_BLANK_LONGLONG)
- ?GAL_DATA_BLANK_SHORT:s[i];}
- while(++i!=size);}
- return;
- case FLOAT_IMG:
- sf=(s=*out)+size; do *s=roundf(*iif++); while(++s<sf);
- if(anyblank)
- {s=*out; do {s[i]=isnan(iiif[i])?GAL_DATA_BLANK_SHORT:s[i];}
- while(++i!=size);}
- return;
- case DOUBLE_IMG:
- sf=(s=*out)+size; do *s=round(*id++); while(++s<sf);
- if(anyblank)
- {s=*out; do {s[i]=isnan(iid[i])?GAL_DATA_BLANK_SHORT:s[i];}
- while(++i!=size);}
- return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
- case LONG_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- lf=(l=*out)+size; do *l=*ib++; while(++l<lf);
- if(anyblank)
- {l=*out; do {l[i]=(iib[i]==GAL_DATA_BLANK_UCHAR)
- ?GAL_DATA_BLANK_LONG:l[i];}
- while(++i!=size);}
- return;
- case SHORT_IMG:
- lf=(l=*out)+size; do *l=*is++; while(++l<lf);
- if(anyblank)
- {l=*out; do {l[i]=(iis[i]==GAL_DATA_BLANK_SHORT)
- ?GAL_DATA_BLANK_LONG:l[i];}
- while(++i!=size);}
- return;
- case LONG_IMG:
- lf=(l=*out)+size; do *l=*il++; while(++l<lf); return;
- case LONGLONG_IMG:
- lf=(l=*out)+size; do *l=*iL++; while(++l<lf);
- if(anyblank)
- {l=*out; do {l[i]=(iiL[i]==GAL_DATA_BLANK_LONGLONG)
- ?GAL_DATA_BLANK_LONG:l[i];}
- while(++i!=size);}
- return;
- case FLOAT_IMG:
- lf=(l=*out)+size; do *l=roundf(*iif++); while(++l<lf);
- if(anyblank)
- {l=*out; do {l[i]=isnan(iiif[i])?GAL_DATA_BLANK_LONG:l[i];}
- while(++i!=size);}
- return;
- case DOUBLE_IMG:
- lf=(l=*out)+size; do *l=round(*id++); while(++l<lf);
- if(anyblank)
- {l=*out; do {l[i]=isnan(iid[i])?GAL_DATA_BLANK_LONG:l[i];}
- while(++i!=size);}
- return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
- case LONGLONG_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- Lf=(L=*out)+size; do *L=*ib++; while(++L<Lf);
- if(anyblank)
- {L=*out; do {L[i]=(iib[i]==GAL_DATA_BLANK_UCHAR)
- ?GAL_DATA_BLANK_LONGLONG:L[i];}
- while(++i!=size);}
- return;
- case SHORT_IMG:
- Lf=(L=*out)+size; do *L=*is++; while(++L<Lf);
- if(anyblank)
- {L=*out; do {L[i]=(iis[i]==GAL_DATA_BLANK_SHORT)
- ?GAL_DATA_BLANK_LONGLONG:L[i];}
- while(++i!=size);}
- return;
- case LONG_IMG:
- Lf=(L=*out)+size; do *L=*il++; while(++L<Lf);
- if(anyblank)
- {L=*out; do {L[i]=(iil[i]==GAL_DATA_BLANK_LONG)
- ?GAL_DATA_BLANK_LONGLONG:L[i];}
- while(++i!=size);}
- return;
- case LONGLONG_IMG:
- Lf=(L=*out)+size; do *L=*iL++; while(++L<Lf); return;
- case FLOAT_IMG:
- Lf=(L=*out)+size; do *L=roundf(*iif++); while(++L<Lf);
- if(anyblank)
- {L=*out; do {L[i]=isnan(iiif[i])?GAL_DATA_BLANK_LONGLONG:L[i];}
- while(++i!=size);}
- return;
- case DOUBLE_IMG:
- Lf=(L=*out)+size; do *L=round(*id++); while(++L<Lf);
- if(anyblank)
- {L=*out; do {L[i]=isnan(iid[i])?GAL_DATA_BLANK_LONGLONG:L[i];}
- while(++i!=size);}
- return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
- case FLOAT_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- ff=(f=*out)+size; do *f=*ib++; while(++f<ff);
- if(anyblank)
- {f=*out; do {f[i]=iib[i]==GAL_DATA_BLANK_UCHAR
- ?GAL_DATA_BLANK_FLOAT:f[i];}
- while(++i!=size);}
- return;
- case SHORT_IMG:
- ff=(f=*out)+size; do *f=*is++; while(++f<ff);
- if(anyblank)
- {f=*out; do {f[i]=iis[i]==GAL_DATA_BLANK_SHORT
- ?GAL_DATA_BLANK_FLOAT:f[i];}
- while(++i!=size);}
- return;
- case LONG_IMG:
- ff=(f=*out)+size; do *f=*il++; while(++f<ff);
- if(anyblank)
- {f=*out; do {f[i]=iil[i]==GAL_DATA_BLANK_LONG
- ?GAL_DATA_BLANK_FLOAT:f[i];}
- while(++i!=size);}
- return;
- case LONGLONG_IMG:
- ff=(f=*out)+size; do *f=*iL++; while(++f<ff);
- if(anyblank)
- {f=*out; do {f[i]=iiL[i]==GAL_DATA_BLANK_LONGLONG
- ?GAL_DATA_BLANK_FLOAT:f[i];}
- while(++i!=size);}
- return;
- case FLOAT_IMG:
- ff=(f=*out)+size; do *f=*iif++; while(++f<ff); return;
- case DOUBLE_IMG:
- ff=(f=*out)+size; do *f=*id++; while(++f<ff); return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
- case DOUBLE_IMG:
- switch(inbitpix)
- {
- case BYTE_IMG:
- df=(d=*out)+size; do *d=*ib++; while(++d<df);
- if(anyblank)
- {d=*out; do {d[i]=iib[i]==GAL_DATA_BLANK_UCHAR
- ?GAL_DATA_BLANK_FLOAT:d[i];}
- while(++i!=size);}
- return;
- case SHORT_IMG:
- df=(d=*out)+size; do *d=*is++; while(++d<df);
- if(anyblank)
- {d=*out; do {d[i]=iis[i]==GAL_DATA_BLANK_SHORT
- ?GAL_DATA_BLANK_FLOAT:d[i];}
- while(++i!=size);}
- return;
- case LONG_IMG:
- df=(d=*out)+size; do *d=*il++; while(++d<df);
- if(anyblank)
- {d=*out; do {d[i]=iil[i]==GAL_DATA_BLANK_LONG
- ?GAL_DATA_BLANK_FLOAT:d[i];}
- while(++i!=size);}
- return;
- case LONGLONG_IMG:
- df=(d=*out)+size; do *d=*iL++; while(++d<df);
- if(anyblank)
- {d=*out; do {d[i]=iiL[i]==GAL_DATA_BLANK_LONGLONG
- ?GAL_DATA_BLANK_FLOAT:d[i];}
- while(++i!=size);}
- return;
- case FLOAT_IMG:
- df=(d=*out)+size; do *d=*iif++; while(++d<df); return;
- case DOUBLE_IMG:
- df=(d=*out)+size; do *d=*id++; while(++d<df); return;
- default:
- error(EXIT_FAILURE, 0, "a bug! In gal_fits_change_type "
- "(fitsarrayvv.c). BITPIX=%d of input not recognized. "
- "Please contact us so we can fix it", inbitpix);
- }
- break;
-
-
- default:
- error(EXIT_FAILURE, 0, "a bug! Output Bitpix value of %d is not "
- "recognized. This should not happen here "
- "(gal_fits_change_type in fitsarrayvv.c). Please "
- "contact us to see how this happened", outbitpix);
- }
-}
-
@@ -836,7 +453,7 @@ gal_fits_change_type(void *in, int inbitpix, size_t size,
int anyblank,
/*************************************************************
- ************** Number of extensions: ***************
+ ************** Get information ***************
*************************************************************/
void
gal_fits_num_hdus(char *filename, int *numhdu)
@@ -861,6 +478,40 @@ gal_fits_num_hdus(char *filename, int *numhdu)
+/* Note that the FITS standard defines any array as an `image',
+ irrespective of how many dimensions it has. */
+void
+gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, long **dsize)
+{
+ size_t i;
+ int bitpix, status=0, naxis;
+ long naxes[GAL_DATA_MAXDIM];
+
+ /* Get the BITPIX, number of dimensions and size of each dimension. */
+ if( fits_get_img_param(fptr, GAL_DATA_MAXDIM, &bitpix, &naxis,
+ naxes, &status) )
+ gal_fits_io_error(status, NULL);
+ *ndim=naxis;
+
+ /* Convert bitpix to Gnuastro's known types. */
+ *type=gal_fits_bitpix_to_type(bitpix);
+
+ /* Allocate space for the size along each dimension. */
+ errno=0;
+ *dsize=malloc( *ndim * sizeof **dsize );
+ if(*dsize==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for dsize in gal_fits_img_info",
+ *ndim * sizeof **dsize);
+
+ /* Put the size of each dimention into the output array. */
+ for(i=0;i<*ndim;++i)
+ (*dsize)[i]=naxes[i];
+}
+
+
+
+
+
@@ -877,7 +528,7 @@ gal_fits_num_hdus(char *filename, int *numhdu)
/**************************************************************/
-/********** Check FITS image HDUs: ************/
+/********** HDU ************/
/**************************************************************/
static char *
hdutypestring(int hdutype)
@@ -926,11 +577,10 @@ gal_fits_read_hdu(char *filename, char *hdu, unsigned
char img0_tab1,
gal_fits_io_error(status, "reading this FITS file");
fptr=*outfptr;
- /* Check the Type of the given HDU: */
+ /* Check the type of the given HDU: */
if (fits_get_hdu_type(fptr, &hdutype, &status) )
gal_fits_io_error(status, NULL);
-
/* Check if the type of the HDU is the expected type. We could have
written these as && conditions, but this is easier to read, it makes
no meaningful difference to the compiler. */
@@ -955,6 +605,24 @@ gal_fits_read_hdu(char *filename, char *hdu, unsigned char
img0_tab1,
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**************************************************************/
+/********** Header keywords ************/
+/**************************************************************/
/* Read keywords from a FITS file. The gal_fits_key pointer is an array of
gal_fits_key structures, which keep the basic information for each
keyword that is to be read and also stores the value in the appropriate
@@ -966,8 +634,7 @@ gal_fits_read_hdu(char *filename, char *hdu, unsigned char
img0_tab1,
`gal_fits_key' structure (to be `FLEN_VALUE' characters, `FLEN_VALUE' is
defined by CFITSIO). So if the value is necessary where `gal_fits_key'
is no longer available, then you have to allocate space dynamically and
- copy the string there.
-*/
+ copy the string there. */
void
gal_fits_read_keywords(char *filename, char *hdu, struct gal_fits_key *keys,
size_t num)
@@ -997,37 +664,37 @@ gal_fits_read_keywords(char *filename, char *hdu, struct
gal_fits_key *keys,
keys[i].status=0;
/* Set the value-pointer based on the required type. */
- switch(keys[i].datatype)
+ switch(keys[i].type)
{
- case TSTRING:
- valueptr=keys[i].str;
- break;
- case TBYTE:
+ case GAL_DATA_TYPE_UCHAR:
valueptr=&keys[i].u;
break;
- case TSHORT:
+ case GAL_DATA_TYPE_STRING:
+ valueptr=keys[i].str;
+ break;
+ case GAL_DATA_TYPE_SHORT:
valueptr=&keys[i].s;
break;
- case TLONG:
+ case GAL_DATA_TYPE_LONG:
valueptr=&keys[i].l;
break;
- case TLONGLONG:
+ case GAL_DATA_TYPE_LONGLONG:
valueptr=&keys[i].L;
break;
- case TFLOAT:
+ case GAL_DATA_TYPE_FLOAT:
valueptr=&keys[i].f;
break;
- case TDOUBLE:
+ case GAL_DATA_TYPE_DOUBLE:
valueptr=&keys[i].d;
break;
default:
error(EXIT_FAILURE, 0, "the value of keys[%zu].datatype (=%d) "
- "is not recognized", i, keys[i].datatype);
+ "is not recognized", i, keys[i].type);
}
/* Read the keyword and place its value in the poitner. */
- fits_read_key(fptr, keys[i].datatype, keys[i].keyname,
- valueptr, NULL, &keys[i].status);
+ fits_read_key(fptr, gal_fits_type_to_datatype(keys[i].type),
+ keys[i].keyname, valueptr, NULL, &keys[i].status);
/* In some cases, the caller might be fine with some kinds of errors,
@@ -1058,30 +725,12 @@ gal_fits_read_keywords(char *filename, char *hdu, struct
gal_fits_key *keys,
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/*************************************************************
- ****************** Header ******************
- *************************************************************/
/* Add on keyword to the list of header keywords that need to be added
to a FITS file. In the end, the keywords will have to be freed, so
it is important to know before hand if they were allocated or
not. If not, they don't need to be freed. */
void
-gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int datatype,
+gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int type,
char *keyname, int kfree, void *value, int vfree,
char *comment, int cfree, char *unit)
{
@@ -1093,7 +742,7 @@ gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int
datatype,
if(newnode==NULL)
error(EXIT_FAILURE, errno,
"linkedlist: new element in gal_fits_key_ll");
- newnode->datatype=datatype;
+ newnode->type=type;
newnode->keyname=keyname;
newnode->value=value;
newnode->comment=comment;
@@ -1111,7 +760,7 @@ gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int
datatype,
void
-gal_fits_add_to_key_ll_end(struct gal_fits_key_ll **list, int datatype,
+gal_fits_add_to_key_ll_end(struct gal_fits_key_ll **list, int type,
char *keyname, int kfree, void *value, int vfree,
char *comment, int cfree, char *unit)
{
@@ -1123,7 +772,7 @@ gal_fits_add_to_key_ll_end(struct gal_fits_key_ll **list,
int datatype,
if(newnode==NULL)
error(EXIT_FAILURE, errno,
"linkedlist: new element in gal_fits_key_ll");
- newnode->datatype=datatype;
+ newnode->type=type;
newnode->keyname=keyname;
newnode->value=value;
newnode->comment=comment;
@@ -1271,14 +920,14 @@ gal_fits_update_keys(fitsfile *fptr, struct
gal_fits_key_ll **keylist)
/* Write the information: */
if(tmp->value)
{
- if( fits_update_key(fptr, tmp->datatype, tmp->keyname,
- tmp->value, tmp->comment, &status) )
+ if( fits_update_key(fptr, gal_fits_type_to_datatype(tmp->type),
+ tmp->keyname, tmp->value, tmp->comment,
+ &status) )
gal_fits_io_error(status, NULL);
}
else
{
- if(fits_update_key_null(fptr, tmp->keyname, tmp->comment,
- &status))
+ if(fits_update_key_null(fptr, tmp->keyname, tmp->comment, &status))
gal_fits_io_error(status, NULL);
}
if(tmp->unit && fits_write_key_unit(fptr, tmp->keyname,
@@ -1318,7 +967,7 @@ gal_fits_write_keys_version(fitsfile *fptr, struct
gal_fits_key_ll *headers,
defined. Sometime in the future were everyone has moved to more
recent versions of WCSLIB, we can remove this macro and its check
in configure.ac.*/
-#ifdef HAVE_WCSLIB_VERSION
+#ifdef GAL_CONFIG_HAVE_WCSLIB_VERSION
int wcslibvers[3];
char wcslibversion[20];
const char *wcslibversion_const;
@@ -1359,7 +1008,7 @@ gal_fits_write_keys_version(fitsfile *fptr, struct
gal_fits_key_ll *headers,
"CFITSIO version.", &status);
/* Write the WCSLIB version */
-#ifdef HAVE_WCSLIB_VERSION
+#ifdef GAL_CONFIG_HAVE_WCSLIB_VERSION
wcslibversion_const=wcslib_version(wcslibvers);
strcpy(wcslibversion, wcslibversion_const);
fits_update_key(fptr, TSTRING, "WCSLIB", wcslibversion,
@@ -1400,8 +1049,11 @@ gal_fits_write_keys_version(fitsfile *fptr, struct
gal_fits_key_ll *headers,
+
+
+
/*************************************************************
- *********** FITS to array functions: ***********
+ *********** Read WCS from FITS pointer ***********
*************************************************************/
/* Read the WCS information from the header. Unfortunately, WCS lib is
not thread safe, so it needs a mutex. In case you are not using
@@ -1524,174 +1176,305 @@ gal_fits_read_wcs(char *filename, char *hdu, size_t
hstartwcs,
+
+
+
+
+
+
+
+
+
+
+
+/*************************************************************
+ *********** Array functions ***********
+ *************************************************************/
+
/* Read a FITS image into an array corresponding to fitstype and also
save the size of the array.
If the image has any null pixels, their number is returned by this
function. The value that is placed for those pixels is defined by
the macros in fitsarrayvv.h and depends on the type of the data.*/
-int
-gal_fits_hdu_to_array(char *filename, char *hdu, int *bitpix,
- void **array, size_t *s0, size_t *s1)
+gal_data_t *
+gal_fits_read_img_hdu(char *filename, char *hdu)
{
- void *bitblank;
+ void *blank;
+ size_t i, ndim;
fitsfile *fptr;
- long naxes[2], fpixel[]={1,1};
- int status=0, anyblank=0, type;
+ gal_data_t *data;
+ int status=0, type;
+ long *fpixel, *dsize;
+
/* Check HDU for realistic conditions: */
gal_fits_read_hdu(filename, hdu, 0, &fptr);
- /* Get the bitpix and size of the image: */
- gal_fits_img_bitpix_size(fptr, bitpix, naxes);
- *s0=naxes[1];
- *s1=naxes[0];
- /* Allocate space for the array. */
- type=gal_fits_bitpix_to_type(*bitpix);
- bitblank=gal_data_alloc_blank(type);
- *array=gal_data_alloc(type, *s0 * *s1);
+ /* Get the info and allocate the data structure. */
+ gal_fits_img_info(fptr, &type, &ndim, &dsize);
+
+
+ /* Check if there is any dimensions (the first header can sometimes have
+ no images). */
+ if(ndim==0)
+ error(EXIT_FAILURE, 0, "%s (hdu: %s) has 0 dimensions! The most common "
+ "cause for this is a wrongly specified HDU: in some FITS images, "
+ "the first HDU doesn't have any data. So probably reading the "
+ "second HDU (with `--hdu=1' or `-h1') will solve the problem. Note "
+ "that currently HDU counting starts from 0." , filename, hdu);
+
+
+ /* Set the fpixel array (first pixel in all dimensions): */
+ errno=0;
+ fpixel=malloc(ndim*sizeof *fpixel);
+ if(fpixel==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for fpixel in gal_fits_read_img_hdu",
+ ndim*sizeof *fpixel);
+ for(i=0;i<ndim;++i) fpixel[i]=1;
+
+
+ /* Allocate the space for the array and for the blank values. */
+ printf("alloc for %s\n", filename);
+ data=gal_data_alloc(type, (long)ndim, dsize, 0, 1);
+ blank=gal_data_alloc_blank(type);
+ free(dsize);
+
/* Read the image into the allocated array: */
fits_read_pix(fptr, gal_fits_type_to_datatype(type), fpixel,
- *s0 * *s1, bitblank, *array, &anyblank, &status);
+ data->size, blank, data->array, &data->anyblank, &status);
if(status) gal_fits_io_error(status, NULL);
- free(bitblank);
+ free(fpixel);
+ free(blank);
- /* Close the FITS file: */
+
+ /* Close the FITS file, and return the data pointer */
fits_close_file(fptr, &status);
gal_fits_io_error(status, NULL);
-
- /* Return the number of blank pixels: */
- return anyblank;
+ return data;
}
+/* The user has specified an input file and a mask file. In the
+ processing, all masked pixels will be converted to NaN pixels in
+ the input image so we only have to deal with one array. Also since
+ all processing is done on floating point arrays, the input is
+ converted to floating point, irrespective of its input type. The
+ original input bitpix will be stored so if you want to, you can
+ return it back to the input type if you please. */
+gal_data_t *
+gal_fits_read_to_type(char *inputname, char *maskname, char *inhdu,
+ char *mhdu, int type)
+{
+ gal_data_t *in, *mask, *converted;
+
+ /* Read the specified input image HDU. */
+ in=gal_fits_read_img_hdu(inputname, inhdu);
+
+ /* If the input had another type, convert it to float. */
+ if(in->type!=type)
+ {
+ converted=gal_data_copy_to_new_type(in, type);
+ gal_data_free(in);
+ in=converted;
+ }
+
+ /* If a mask was specified, read it as a float image, then set all
+ the corresponding pixels of the input image to NaN. */
+ if(maskname)
+ {
+ /* Read the mask HDU. */
+ mask=gal_fits_read_img_hdu(maskname, mhdu);
+
+ /* Apply the mask on the input. */
+ gal_data_apply_mask(in, mask);
+ /* Free the mask space */
+ gal_data_free(mask);
+ }
+ return in;
+}
+gal_data_t *
+gal_fits_read_float_kernel(char *inputname, char *inhdu, float **outkernel,
+ size_t *ins0, size_t *ins1)
+{
+ size_t i;
+ int check=0;
+ double sum=0;
+ gal_data_t *kernel;
+ float *f, *fp, tmp;
+
+ /* Read the image as a float */
+ kernel=gal_fits_read_to_type(inputname, NULL, inhdu, NULL,
+ GAL_DATA_TYPE_FLOAT);
+
+ /* Check if the size along each dimension of the kernel is an odd
+ number. If they are all an odd number, then the for each dimension,
+ check will be incremented once. */
+ for(i=0;i<kernel->ndim;++i)
+ check=kernel->dsize[i]%2;
+ if(check!=kernel->ndim)
+ error(EXIT_FAILURE, 0, "the kernel image has to have an odd number "
+ "of pixels in all dimensions (there has to be one element/pixel "
+ "in the center). At least one of the dimensions of %s (hdu: %s) "
+ "doesn't have an odd number of pixels", inputname, inhdu);
+ /* If there are any NaN pixels, set them to zero and normalize it. There
+ are no blank pixels any more.*/
+ fp=(f=kernel->array)+kernel->size;
+ do
+ {
+ if(isnan(*f)) *f=0.0f;
+ else sum+=*f;
+ }
+ while(++f<fp);
+ kernel->anyblank=0;
+ f=kernel->array; do *f++ *= 1/sum; while(f<fp);
+ /* Flip the kernel about the center (necessary for non-symmetric
+ kernels). */
+ f=kernel->array;
+ for(i=0;i<kernel->size/2;++i)
+ {
+ tmp=f[i];
+ f[i]=f[ kernel->size - i - 1 ];
+ f[ kernel->size - i - 1 ]=tmp;
+ }
+ /* Return the kernel*/
+ return kernel;
+}
-/*************************************************************
- ****************** Array to FITS ******************
- *************************************************************/
-void
-gal_fits_array_to_file(char *filename, char *extname, int bitpix, void *array,
- size_t s0, size_t s1, int anyblank, struct wcsprm *wcs,
- struct gal_fits_key_ll *headers, char *spack_string)
+/* This function will write all the data array information (including its
+ WCS information) into a FITS file, but will not close it. Instead it
+ will pass along the FITS pointer for further modification. */
+fitsfile *
+gal_fits_write_img_fitsptr(gal_data_t *data, char *filename, char *extname)
{
- int nkeyrec;
void *blank;
+ long fpixel=1;
fitsfile *fptr;
char *wcsheader;
- int status=0, type, datatype;
- long fpixel=1, naxis=2, nelements, naxes[]={s1,s0};
-
- type=gal_fits_bitpix_to_type(bitpix);
- datatype=gal_fits_type_to_datatype(type);
- nelements=naxes[0]*naxes[1];
+ int nkeyrec, status=0, datatype=gal_fits_type_to_datatype(data->type);
+ /* Check if the file already exists. If it does, we want to add the array
+ as a new extension. */
if(access(filename,F_OK) != -1 )
fits_open_file(&fptr,filename, READWRITE, &status);
else
fits_create_file(&fptr, filename, &status);
- fits_create_img(fptr, bitpix, naxis, naxes, &status);
- fits_write_img(fptr, datatype, fpixel, nelements, array, &status);
-
- if(anyblank)
- if(bitpix==BYTE_IMG || bitpix==SHORT_IMG
- || bitpix==LONG_IMG || bitpix==LONGLONG_IMG)
+ /* Create the FITS file and put the image into it. */
+ fits_create_img(fptr, gal_fits_type_to_bitpix(data->type),
+ data->ndim, data->dsize, &status);
+ fits_write_img(fptr, datatype, fpixel, data->size, data->array, &status);
+
+ /* If we have blank pixels, we need to define a BLANK keyword when we are
+ dealing with integer types. */
+ if(data->anyblank)
+ if( data->type==GAL_DATA_TYPE_UCHAR || data->type==GAL_DATA_TYPE_CHAR
+ || data->type==GAL_DATA_TYPE_USHORT || data->type==GAL_DATA_TYPE_SHORT
+ || data->type==GAL_DATA_TYPE_UINT || data->type==GAL_DATA_TYPE_INT
+ || data->type==GAL_DATA_TYPE_ULONG || data->type==GAL_DATA_TYPE_LONG
+ || data->type==GAL_DATA_TYPE_LONGLONG)
{
- blank=gal_data_alloc_blank(type);
+ blank=gal_data_alloc_blank(data->type);
if(fits_write_key(fptr, datatype, "BLANK", blank,
"Pixels with no data.", &status) )
gal_fits_io_error(status, "adding the BLANK keyword");
free(blank);
}
+ /* Write the extension name to the header. */
fits_write_key(fptr, TSTRING, "EXTNAME", extname, "", &status);
gal_fits_io_error(status, NULL);
- if(wcs)
+ /* If a WCS structure is present, write it in */
+ if(data->wcs)
{
/* Convert the WCS information to text. */
- status=wcshdo(WCSHDO_safe, wcs, &nkeyrec, &wcsheader);
+ status=wcshdo(WCSHDO_safe, data->wcs, &nkeyrec, &wcsheader);
if(status)
error(EXIT_FAILURE, 0, "wcshdo ERROR %d: %s", status,
wcs_errmsg[status]);
gal_fits_add_wcs_to_header(fptr, wcsheader, nkeyrec);
}
- gal_fits_write_keys_version(fptr, headers, spack_string);
-
- fits_close_file(fptr, &status);
+ /* Report any errors if we had any */
gal_fits_io_error(status, NULL);
+ return fptr;
}
-/* `atof' stands for "array to file". This is essentially the same as
- `gal_fits_array_to_file' except that the WCS structure's CRPIX values
- have changed. */
void
-gal_fits_atof_correct_wcs(char *filename, char *hdu, int bitpix, void *array,
- size_t s0, size_t s1, char *wcsheader, int
wcsnkeyrec,
- double *crpix, char *spack_string)
+gal_fits_write_img(gal_data_t *data, char *filename, char *extname,
+ struct gal_fits_key_ll *headers, char *spack_string)
{
+ int status=0;
fitsfile *fptr;
- int status=0, type, datatype;
- long fpixel=1, naxis=2, nelements, naxes[]={s1,s0};
- type=gal_fits_bitpix_to_type(bitpix);
- datatype=gal_fits_type_to_datatype(type);
- nelements=naxes[0]*naxes[1];
+ /* Write the data array into a FITS file and keep it open: */
+ fptr=gal_fits_write_img_fitsptr(data, filename, extname);
- if(access(filename,F_OK) != -1 )
- fits_open_file(&fptr,filename, READWRITE, &status);
- else
- fits_create_file(&fptr, filename, &status);
-
- fits_create_img(fptr, bitpix, naxis, naxes, &status);
- fits_write_img(fptr, datatype, fpixel, nelements, array, &status);
+ /* Write all the headers and the version information. */
+ gal_fits_write_keys_version(fptr, headers, spack_string);
- fits_write_key(fptr, TSTRING, "EXTNAME", hdu, "", &status);
+ /* Close the FITS file. */
+ fits_close_file(fptr, &status);
gal_fits_io_error(status, NULL);
+}
+
+
- fits_delete_key(fptr, "COMMENT", &status);
- fits_delete_key(fptr, "COMMENT", &status);
- gal_fits_io_error(status, NULL);
- if(wcsheader)
+
+void
+gal_fits_write_img_update_crpix(gal_data_t *data, char *filename,
+ char *extname,
+ struct gal_fits_key_ll *headers,
+ double *crpix, char *spack_string)
+{
+ int status=0;
+ fitsfile *fptr;
+
+ /* Write the data array into a FITS file and keep it open: */
+ fptr=gal_fits_write_img_fitsptr(data, filename, extname);
+
+ /* Update the CRPIX keywords. Note that we don't want to change the
+ values in the WCS information of gal_data_t. Because, it often happens
+ that things are done in parallel, so we don't want to touch the
+ original version, we just want to change the copied version. */
+ if(crpix)
{
- gal_fits_add_wcs_to_header(fptr, wcsheader, wcsnkeyrec);
- if(crpix)
- {
- fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix[0],
- NULL, &status);
- fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix[1],
- NULL, &status);
- gal_fits_io_error(status, NULL);
- }
+ fits_update_key(fptr, TDOUBLE, "CRPIX1", &crpix[0],
+ NULL, &status);
+ fits_update_key(fptr, TDOUBLE, "CRPIX2", &crpix[1],
+ NULL, &status);
+ gal_fits_io_error(status, NULL);
}
+ /* Write all the headers and the version information. */
gal_fits_write_keys_version(fptr, NULL, spack_string);
+ /* Close the file and return. */
fits_close_file(fptr, &status);
gal_fits_io_error(status, NULL);
}
@@ -1774,270 +1557,3 @@ gal_fits_table_type(fitsfile *fptr)
"the end of the function! This must not happen", PACKAGE_BUGREPORT);
return -1;
}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-/**************************************************************/
-/********** Check prepare file ************/
-/**************************************************************/
-/* We have the name of the input file. But in most cases, the files
- that should be used (for example a mask image) are other extensions
- in the same file. So the user only has to give the HDU. The job of
- this function is to determine which is the case and set othername
- to the appropriate value. */
-void
-gal_fits_file_or_ext_name(char *inputname, char *inhdu, int othernameset,
- char **othername, char *ohdu, int ohduset,
- char *type)
-{
- if(othernameset)
- {
- /* In some cases, for example a mask image, both the name and
- HDU are optional. So just to be safe, we will check this all
- the time. */
- if(ohduset==0)
- error(EXIT_FAILURE, 0, "a %s image was specified (%s). However, "
- "no HDU is given for it. Please add a HDU. If you regularly "
- "use the same HDU as %s, you may consider adding it to "
- "the configuration file. For more information, please see "
- "the `Configuration files' section of the %s manual by "
- "running ` info gnuastro ' on the command-line", type,
- *othername, type, PACKAGE_NAME);
- if(strcmp(inputname, *othername)==0)
- {
- if(strcmp(ohdu, inhdu)==0)
- error(EXIT_FAILURE, 0, "the specified %s name and "
- "input image name (%s) are the same while the input "
- "image hdu name and mask hdu are also identical (%s)",
- type, inputname, inhdu);
- }
- }
- else if(ohduset && strcmp(ohdu, inhdu))
- *othername=inputname;
- else
- *othername=NULL;
-}
-
-
-
-
-
-/* The user has specified an input file and a mask file. In the
- processing, all masked pixels will be converted to NaN pixels in
- the input image so we only have to deal with one array. Also since
- all processing is done on floating point arrays, the input is
- converted to floating point, irrespective of its input type. The
- original input bitpix will be stored so if you want to, you can
- return it back to the input type if you please. */
-void
-gal_fits_file_to_float(char *inputname, char *maskname, char *inhdu,
- char *mhdu, float **img, int *inbitpix,
- int *anyblank, size_t *ins0, size_t *ins1)
-{
- void *array;
- int maskbitpix;
- float *mask, *f, *ff, *fp;
- size_t maskanyblank, s0, s1;
-
- /* Read the input array and convert it to float. */
- *anyblank=gal_fits_hdu_to_array(inputname, inhdu, inbitpix,
- &array, ins0, ins1);
- if(*inbitpix==FLOAT_IMG)
- *img=array;
- else
- {
- gal_fits_change_type(array, *inbitpix, *ins0 * *ins1, *anyblank,
- (void **)img, FLOAT_IMG);
- free(array);
- }
-
- /* If a mask was specified, read it as a float image, then set all
- the corresponding pixels of the input image to NaN. */
- if(maskname)
- {
- maskanyblank=gal_fits_hdu_to_array(maskname, mhdu, &maskbitpix,
- &array, &s0, &s1);
-
- if(maskbitpix==FLOAT_IMG || maskbitpix==DOUBLE_IMG)
- fprintf(stderr, "WARNING: the mask image (%s, hdu: %s) has a %s "
- "precision floating point data type (BITPIX=%d). The mask "
- "image is usually an integer type. Therefore this might "
- "be due to a mistake in the inputs and the results might "
- "not be what you intended. However, the program will not "
- "abort and continue working only with zero valued pixels "
- "in the given masked image.", maskname, mhdu,
- maskbitpix==FLOAT_IMG ? "single" : "double", maskbitpix);
-
- if(s0!=*ins0 || s1!=*ins1)
- error(EXIT_FAILURE, 0, "the input image %s (hdu: %s) has size: "
- "%zu x %zu. The mask image %s (hdu: %s) has size %zu x %zu. "
- "The two images have to have the same size", inputname,
- inhdu, *ins1, *ins0, maskname, mhdu, s1, s0);
-
- if(maskbitpix==FLOAT_IMG)
- mask=array;
- else
- {
- gal_fits_change_type(array, maskbitpix, *ins0 * *ins1,
- maskanyblank, (void **)(&mask), FLOAT_IMG);
- free(array);
- }
-
- ff=mask;
- fp=(f=*img)+s0*s1;
- do if(*ff++!=0.0f) {*f=NAN; ++(*anyblank);} while(++f<fp);
- free(mask);
- }
-}
-
-
-
-
-/* Similar to filetofloat, but for double type */
-void
-gal_fits_file_to_double(char *inputname, char *maskname, char *inhdu,
- char *mhdu, double **img, int *inbitpix,
- int *anyblank, size_t *ins0, size_t *ins1)
-{
- void *array;
- int maskbitpix;
- double *mask, *f, *ff, *fp;
- size_t maskanyblank, s0, s1;
-
- /* Read the input array and convert it to double. */
- *anyblank=gal_fits_hdu_to_array(inputname, inhdu, inbitpix,
- &array, ins0, ins1);
- if(*inbitpix==DOUBLE_IMG)
- *img=array;
- else
- {
- gal_fits_change_type(array, *inbitpix, *ins0 * *ins1, *anyblank,
- (void **)img, DOUBLE_IMG);
- free(array);
- }
-
- /* If a mask was specified, read it as a double image, then set all
- the corresponding pixels of the input image to NaN. */
- if(maskname)
- {
- maskanyblank=gal_fits_hdu_to_array(maskname, mhdu, &maskbitpix,
- &array, &s0, &s1);
-
- if(maskbitpix==FLOAT_IMG || maskbitpix==DOUBLE_IMG)
- fprintf(stderr, "WARNING: the mask image (%s, hdu: %s) has a %s "
- "precision floating point data type (BITPIX=%d). The mask "
- "image is usually an integer type. Therefore this might "
- "be due to a mistake in the inputs and the results might "
- "not be what you intended. However, the program will not "
- "abort and continue working only with zero valued pixels in "
- "the given masked image.", maskname, mhdu,
- maskbitpix==FLOAT_IMG ? "single" : "double", maskbitpix);
-
- if(s0!=*ins0 || s1!=*ins1)
- error(EXIT_FAILURE, 0, "the input image %s (hdu: %s) has size: "
- "%zu x %zu. The mask image %s (hdu: %s) has size %zu x %zu. "
- "The two images have to have the same size", inputname,
- inhdu, *ins1, *ins0, maskname, mhdu, s1, s0);
-
- if(maskbitpix==DOUBLE_IMG)
- mask=array;
- else
- {
- gal_fits_change_type(array, maskbitpix, *ins0 * *ins1,
- maskanyblank, (void **)(&mask),
- DOUBLE_IMG);
- free(array);
- }
-
- ff=mask;
- fp=(f=*img)+s0*s1;
- do if(*ff++!=0.0f) {*f=NAN; ++(*anyblank);} while(++f<fp);
- free(mask);
- }
-}
-
-
-
-
-
-void
-gal_fits_file_to_long(char *inputname, char *inhdu, long **img,
- int *inbitpix, int *anyblank, size_t *ins0,
- size_t *ins1)
-{
- void *array;
-
- /* Read the input array and convert it to float. */
- *anyblank=gal_fits_hdu_to_array(inputname, inhdu, inbitpix,
- &array, ins0, ins1);
- if(*inbitpix==LONG_IMG)
- *img=array;
- else
- {
- gal_fits_change_type(array, *inbitpix, *ins0 * *ins1, *anyblank,
- (void **)img, LONG_IMG);
- free(array);
- }
-}
-
-
-
-
-
-void
-gal_fits_prep_float_kernel(char *inputname, char *inhdu, float **outkernel,
- size_t *ins0, size_t *ins1)
-{
- size_t i, size;
- double sum=0.0f;
- int bitpix, anyblank;
- float *f, *fp, *kernel, tmp;
-
- /* Read the kernel as a float array: */
- gal_fits_file_to_float(inputname, NULL, inhdu, NULL, outkernel, &bitpix,
- &anyblank, ins0, ins1);
- size = *ins0 * *ins1;
- kernel=*outkernel;
-
- /* Simple sanity check: */
- if(*ins0%2==0 || *ins1%2==0)
- error(EXIT_FAILURE, 0, "the kernel image has to have an odd number "
- "of pixels on both sides (there has to be on pixel in the "
- "center). %s (hdu: %s) is %zu by %zu", inputname, inhdu,
- *ins1, *ins0);
-
- /* If there are any NaN pixels, set them to zero and normalize it.*/
- fp=(f=kernel)+size;
- do
- {
- if(isnan(*f)) *f=0.0f;
- else sum+=*f;
- }
- while(++f<fp);
- f=kernel; do *f++ *= 1/sum; while(f<fp);
-
- /* Flip the kernel: */
- for(i=0;i<size/2;++i)
- {
- tmp=kernel[i];
- kernel[i]=kernel[size-i-1];
- kernel[size-i-1]=tmp;
- }
-}
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 32deba8..ea1efa8 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -29,6 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <limits.h>
#include <stdint.h>
+#include <fitsio.h> /* Only for the LONGLONG type */
#include <wcslib/wcs.h>
#include <gsl/gsl_complex.h>
@@ -56,7 +57,7 @@ __BEGIN_C_DECLS /* From C++ preparations */
/* Macros: */
/* The maximum dimensionality of datasets. */
-#define GAL_DATA_MAXDIM 30
+#define GAL_DATA_MAXDIM 999
/* Blank values: Note that for the unsigned types or small types (like
char), the maximum value is considered as a blank value, since the
@@ -108,34 +109,80 @@ enum gal_data_alltypes
/* Main data structure
- If mmaped==0, it is assumed that the data is allocated (using malloc).
+ If mmaped==0, it is assumed that the data is allocated (using
+ malloc). The `dsize' array is in the `long' type because CFITSIO uses
+ the long type and this will make it easier to call CFITSIO functions.
*/
typedef struct
{
void *array; /* Array keeping data elements. */
int type; /* Type of data (from `gal_data_alltypes'). */
size_t ndim; /* Number of dimensions in the array. */
- size_t size[GAL_DATA_MAXDIM]; /* Length along each dimension. */
+ long *dsize; /* Size of array along each dimension. */
+ size_t size; /* Total number of data-elements. */
int mmapped; /* ==1: not in physical RAM, it is mmap'd. */
char *mmapname; /* File name of the mmap. */
int anyblank; /* ==1: has blank values. */
+ int nwcs; /* for WCSLIB: no. coord. representations. */
struct wcsprm *wcs; /* WCS information for this dataset. */
-} gal_data;
+} gal_data_t;
-/* Functions */
+/*********************************************************************/
+/************* Size and allocation *******************/
+/*********************************************************************/
+int
+gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second);
+
size_t
gal_data_sizeof(int type);
void *
-gal_data_alloc(int type, size_t size);
+gal_data_malloc_array(int type, size_t size);
+
+void *
+gal_data_calloc_array(int type, size_t size);
+
+gal_data_t *
+gal_data_alloc(int type, size_t ndim, long *dsize, int clear, int map);
+
+void
+gal_data_free(gal_data_t *data);
+
+
+
+
+
+/*************************************************************
+ ************** Blank data ***************
+ *************************************************************/
void *
gal_data_alloc_blank(int type);
+void
+gal_data_apply_mask(gal_data_t *in, gal_data_t *mask);
+
+void
+gal_data_blank_to_value(gal_data_t *data, void *value);
+
+
+
+
+
+/*************************************************************
+ ************** Convert types ***************
+ *************************************************************/
+int
+gal_data_out_type(gal_data_t *first, gal_data_t *second);
+
+gal_data_t *
+gal_data_copy_to_new_type(gal_data_t *in, int newtype);
+
+
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 1cd8531..6d99114 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -23,6 +23,17 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef __GAL_FITS_H__
#define __GAL_FITS_H__
+/* When we are within Gnuastro's building process, `IN_GNUASTRO_BUILD' is
+ defined. In the build process, installation information (in particular
+ `GAL_CONFIG_HAVE_WCSLIB_VERION' that we need in `fits.c') is kept in
+ `config.h'. When building a user's programs, this information is kept in
+ `gnuastro/config.h'. Note that all `.c' files must start with the
+ inclusion of `config.h' and that `gnuastro/config.h' is only created at
+ installation time (not present during the building of Gnuastro).*/
+#ifndef IN_GNUASTRO_BUILD
+#include <gnuastro/config.h>
+#endif
+
/* Include other headers if necessary here. Note that other header files
must be included before the C++ preparations below */
#include <math.h>
@@ -56,32 +67,13 @@ __BEGIN_C_DECLS /* From C++ preparations */
-/*************************************************************
- ****************** Basic *******************
- *************************************************************/
-void
-gal_fits_io_error(int status, char *message);
-
-int
-gal_fits_name_is_fits(char *name);
-
-int
-gal_fits_suffix_is_fits(char *suffix);
-
-
-
-
-
-/*************************************************************
- ****************** Header ******************
- *************************************************************/
/* To create a linked list of headers. */
struct gal_fits_key_ll
{
int kfree; /* ==1, free keyword name. */
int vfree; /* ==1, free keyword value. */
int cfree; /* ==1, free comment. */
- int datatype; /* Keyword value datatype. */
+ int type; /* Keyword value type. */
char *keyname; /* Keyword Name. */
void *value; /* Keyword value. */
char *comment; /* Keyword comment. */
@@ -97,7 +89,7 @@ struct gal_fits_key
{
int status; /* CFITSIO status. */
char *keyname; /* Name of keyword. */
- int datatype; /* Type of keyword value. */
+ int type; /* Type of keyword value. */
char str[FLEN_VALUE]; /* String value. */
unsigned char u; /* Byte value. */
short s; /* Short integer value. */
@@ -109,37 +101,29 @@ struct gal_fits_key
-
-
+/*************************************************************
+ ************** Reporting errors: ***************
+ *************************************************************/
void
-gal_fits_read_keywords(char *filename, char *hdu, struct gal_fits_key *out,
- size_t num);
+gal_fits_io_error(int status, char *message);
-void
-gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int datatype,
- char *keyname, int kfree, void *value, int vfree,
- char *comment, int cfree, char *unit);
-void
-gal_fits_add_to_key_ll_end(struct gal_fits_key_ll **list, int datatype,
- char *keyname, int kfree, void *value, int vfree,
- char *comment, int cfree, char *unit);
-
-void
-gal_fits_file_name_in_keywords(char *keynamebase, char *filename,
- struct gal_fits_key_ll **list);
-void
-gal_fits_add_wcs_to_header(fitsfile *fptr, char *wcsheader, int nkeyrec);
-void
-gal_fits_update_keys(fitsfile *fptr, struct gal_fits_key_ll **keylist);
-void
-gal_fits_write_keys_version(fitsfile *fptr, struct gal_fits_key_ll *headers,
- char *spack_string);
+/*************************************************************
+ ************** FITS names ***************
+ *************************************************************/
+int
+gal_fits_name_is_fits(char *name);
+int
+gal_fits_suffix_is_fits(char *suffix);
+void
+gal_fits_file_or_ext_name(char *inputname, char *inhdu, int othernameset,
+ char **othername, char *ohdu, int ohduset,
+ char *type);
@@ -150,6 +134,9 @@ int
gal_fits_bitpix_to_type(int bitpix);
int
+gal_fits_type_to_bitpix(int type);
+
+int
gal_fits_tform_to_type(char tform);
int
@@ -163,100 +150,120 @@ gal_fits_datatype_to_type(int type);
/*************************************************************
- ****************** Read/Write *****************
+ ************** Get information ***************
*************************************************************/
void
-gal_fits_convert_blank(void *array, int bitpix, size_t size, void *value);
+gal_fits_num_hdus(char *filename, int *numhdu);
void
-gal_fits_blank_to_value(void *array, int datatype, size_t size, void *value);
+gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, long **dsize);
-void
-gal_fits_img_bitpix_size(fitsfile *fptr, int *bitpix, long *naxis);
+
+
+
+
+/**************************************************************/
+/********** HDU ************/
+/**************************************************************/
void
gal_fits_read_hdu(char *filename, char *hdu, unsigned char img0_tab1,
fitsfile **outfptr);
+
+
+
+
+/**************************************************************/
+/********** Header keywords ************/
+/**************************************************************/
void
-gal_fits_change_type(void *in, int inbitpix, size_t size, int anyblank,
- void **out, int outbitpix);
+gal_fits_read_keywords(char *filename, char *hdu, struct gal_fits_key *out,
+ size_t num);
void
-gal_fits_num_hdus(char *filename, int *numhdu);
+gal_fits_add_to_key_ll(struct gal_fits_key_ll **list, int datatype,
+ char *keyname, int kfree, void *value, int vfree,
+ char *comment, int cfree, char *unit);
void
-gal_fits_read_wcs_from_pointer(fitsfile *fptr, int *nwcs,
- struct wcsprm **wcs,
- size_t hstart, size_t hend);
+gal_fits_add_to_key_ll_end(struct gal_fits_key_ll **list, int datatype,
+ char *keyname, int kfree, void *value, int vfree,
+ char *comment, int cfree, char *unit);
void
-gal_fits_read_wcs(char *filename, char *hdu, size_t hstartwcs,
- size_t hendwcs, int *nwcs, struct wcsprm **wcs);
+gal_fits_file_name_in_keywords(char *keynamebase, char *filename,
+ struct gal_fits_key_ll **list);
-int
-gal_fits_hdu_to_array(char *filename, char *hdu, int *bitpix,
- void **array, size_t *s0, size_t *s1);
+void
+gal_fits_add_wcs_to_header(fitsfile *fptr, char *wcsheader, int nkeyrec);
void
-gal_fits_array_to_file(char *filename, char *hdu, int bitpix,
- void *array, size_t s0, size_t s1, int anyblank,
- struct wcsprm *wcs, struct gal_fits_key_ll *headers,
- char *spack_string);
+gal_fits_update_keys(fitsfile *fptr, struct gal_fits_key_ll **keylist);
void
-gal_fits_atof_correct_wcs(char *filename, char *hdu, int bitpix,
- void *array, size_t s0, size_t s1,
- char *wcsheader, int wcsnkeyrec,
- double *crpix, char *spack_string);
+gal_fits_write_keys_version(fitsfile *fptr, struct gal_fits_key_ll *headers,
+ char *spack_string);
-/**************************************************************/
-/********** Table ************/
-/**************************************************************/
+/*************************************************************
+ *********** Read WCS from FITS pointer ***********
+ *************************************************************/
void
-gal_fits_table_size(fitsfile *fitsptr, size_t *nrows, size_t *ncols);
+gal_fits_read_wcs_from_pointer(fitsfile *fptr, int *nwcs, struct wcsprm **wcs,
+ size_t hstartwcs, size_t hendwcs);
-int
-gal_fits_table_type(fitsfile *fptr);
+void
+gal_fits_read_wcs(char *filename, char *hdu, size_t hstartwcs,
+ size_t hendwcs, int *nwcs, struct wcsprm **wcs);
-/**************************************************************/
-/********** Check prepare file ************/
-/**************************************************************/
-void
-gal_fits_file_or_ext_name(char *inputname, char *inhdu, int othernameset,
- char **othername, char *ohdu, int ohduset,
- char *type);
+/*************************************************************
+ ****************** Array functions *****************
+ *************************************************************/
+gal_data_t *
+gal_fits_read_img_hdu(char *filename, char *hdu);
-void
-gal_fits_set_mask_name(char *inputname, char **maskname, char *inhdu,
- char *mhdu);
+gal_data_t *
+gal_fits_read_to_type(char *inputname, char *maskname, char *inhdu,
+ char *mhdu, int type);
-void
-gal_fits_file_to_double(char *inputname, char *maskname, char *inhdu,
- char *mhdu, double **img, int *inbitpix,
- int *anyblank, size_t *ins0, size_t *ins1);
+gal_data_t *
+gal_fits_read_float_kernel(char *inputname, char *inhdu, float **outkernel,
+ size_t *ins0, size_t *ins1);
+
+fitsfile *
+gal_fits_write_img_fitsptr(gal_data_t *data, char *filename, char *extname);
void
-gal_fits_file_to_float(char *inputname, char *maskname, char *inhdu,
- char *mhdu, float **img, int *inbitpix,
- int *anyblank, size_t *ins0, size_t *ins1);
+gal_fits_write_img(gal_data_t *data, char *filename, char *extname,
+ struct gal_fits_key_ll *headers, char *spack_string);
void
-gal_fits_file_to_long(char *inputname, char *inhdu, long **img,
- int *inbitpix, int *anyblank, size_t *ins0,
- size_t *ins1);
+gal_fits_write_img_update_crpix(gal_data_t *data, char *filename,
+ char *extname,
+ struct gal_fits_key_ll *headers,
+ double *crpix, char *spack_string);
+
+
+
+
+/**************************************************************/
+/********** Table ************/
+/**************************************************************/
void
-gal_fits_prep_float_kernel(char *inputname, char *inhdu, float **kernel,
- size_t *ins0, size_t *ins1);
+gal_fits_table_size(fitsfile *fitsptr, size_t *nrows, size_t *ncols);
+
+int
+gal_fits_table_type(fitsfile *fptr);
+
diff --git a/lib/gnuastro/linkedlist.h b/lib/gnuastro/linkedlist.h
index 34155d1..87f4952 100644
--- a/lib/gnuastro/linkedlist.h
+++ b/lib/gnuastro/linkedlist.h
@@ -138,6 +138,9 @@ gal_linkedlist_print_stll(struct gal_linkedlist_stll *list);
size_t
gal_linkedlist_num_in_stll(struct gal_linkedlist_stll *list);
+void
+gal_linkedlist_free_stll(struct gal_linkedlist_stll *list, int freevalue);
+
/******************* size_t: */
diff --git a/lib/gnuastro/threads.h b/lib/gnuastro/threads.h
index 1e8a15b..48d3b1e 100644
--- a/lib/gnuastro/threads.h
+++ b/lib/gnuastro/threads.h
@@ -27,14 +27,13 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
must be included before the C++ preparations below */
#include <pthread.h>
-/* When this header is included within Gnuastro's building process,
- `IN_GNUASTRO_BUILD' is defined. In the build process, installation
- information (in particular `GAL_CONFIG_HAVE_PTHREAD_BARRIER' that we
- need below) is kept in `config.h'. When building a user's programs, this
- information is kept in `gnuastro/config.h'. Note that all `.c' files
- must start with the inclusion of `config.h' and that
- `gnuastro/config.h' is only created at installation time (not present
- during the building of Gnuastro).*/
+/* When we are within Gnuastro's building process, `IN_GNUASTRO_BUILD' is
+ defined. In the build process, installation information (in particular
+ `GAL_CONFIG_HAVE_PTHREAD_BARRIER' that we need below) is kept in
+ `config.h'. When building a user's programs, this information is kept in
+ `gnuastro/config.h'. Note that all `.c' files must start with the
+ inclusion of `config.h' and that `gnuastro/config.h' is only created at
+ installation time (not present during the building of Gnuastro).*/
#ifndef IN_GNUASTRO_BUILD
#include <gnuastro/config.h>
#endif
diff --git a/lib/linkedlist.c b/lib/linkedlist.c
index b09df81..8ad9ee0 100644
--- a/lib/linkedlist.c
+++ b/lib/linkedlist.c
@@ -381,6 +381,24 @@ gal_linkedlist_num_in_stll(struct gal_linkedlist_stll
*list)
+void
+gal_linkedlist_free_stll(struct gal_linkedlist_stll *list, int freevalue)
+{
+ struct gal_linkedlist_stll *tmp;
+ while(list!=NULL)
+ {
+ tmp=list->next;
+ if(freevalue)
+ free(list->v);
+ free(list);
+ list=tmp;
+ }
+}
+
+
+
+
+
diff --git a/lib/mesh.c b/lib/mesh.c
index 9f56bc4..1f89a6f 100644
--- a/lib/mesh.c
+++ b/lib/mesh.c
@@ -453,8 +453,18 @@ gal_mesh_value_file(struct gal_mesh_params *mp, char
*filename,
char *extname1, char *extname2, struct wcsprm *wcs,
char *spack_string)
{
+ long dsize[2];
+ gal_data_t data;
float *tmp1=NULL, *tmp2=NULL;
+ /* Set the basic shared properties of the dataset. Note that the float
+ image is written as a float FITS file, so anyblank is irrelevant, so
+ is `size'.*/
+ data.ndim=2;
+ data.anyblank=0;
+ data.dsize=dsize;
+ data.type=GAL_DATA_TYPE_FLOAT;
+
if(mp->meshbasedcheck)
{
/* We want one pixel per mesh. If the last operation was on
@@ -463,29 +473,32 @@ gal_mesh_value_file(struct gal_mesh_params *mp, char
*filename,
used for this job. In cgarray the meshs are ordered
differently. */
if(mp->garray1==mp->cgarray1) gal_mesh_full_garray(mp, 0);
- gal_fits_array_to_file(filename, extname1, FLOAT_IMG,
- mp->fgarray1, mp->gs0*mp->nch2,
- mp->gs1*mp->nch1, 0, wcs, NULL,
- spack_string);
+ data.wcs=NULL; /* This is not the original image size to have same WCS */
+ data.array=mp->fgarray1;
+ data.dsize[0]=mp->gs0*mp->nch2;
+ data.dsize[1]=mp->gs1*mp->nch1;
+ gal_fits_write_img(&data, filename, extname1, NULL, spack_string);
if(mp->ngarrays==2)
- /* Note that gal_mesh_full_garray will correct both the meshs if
- there are two.*/
- gal_fits_array_to_file(filename, extname2, FLOAT_IMG,
- mp->fgarray2, mp->gs0*mp->nch2,
- mp->gs1*mp->nch1, 0, wcs, NULL,
- spack_string);
-
+ {
+ /* Note that gal_mesh_full_garray will correct both the meshs if
+ there are two.*/
+ data.array=mp->fgarray2;
+ gal_fits_write_img(&data, filename, extname2, NULL, spack_string);
+ }
}
else
{
gal_mesh_check_garray(mp, &tmp1, &tmp2);
- gal_fits_array_to_file(filename, extname1, FLOAT_IMG, tmp1,
- mp->s0, mp->s1, 0, wcs, NULL,
- spack_string);
+ data.wcs=wcs;
+ data.array=tmp1;
+ data.dsize[0]=mp->s0;
+ data.dsize[1]=mp->s1;
+ gal_fits_write_img(&data, filename, extname1, NULL, spack_string);
if(mp->ngarrays==2)
- gal_fits_array_to_file(filename, extname2, FLOAT_IMG, tmp2,
- mp->s0, mp->s1, 0, wcs, NULL,
- spack_string);
+ {
+ data.array=tmp2;
+ gal_fits_write_img(&data, filename, extname2, NULL, spack_string);
+ }
free(tmp1);
free(tmp2);
}
- [gnuastro-commits] master e353932 109/125: NaN magnitude will only blank the profile's area, (continued)
- [gnuastro-commits] master e353932 109/125: NaN magnitude will only blank the profile's area, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 96ab804 098/125: Working on tiles implemented in tile checking, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c0dd6db 102/125: The noedgecorrection option in spatial domain convolution, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dfdc34 069/125: Option values directly set in library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db4c16d 076/125: Pulled in all the recent work in master, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 64f8db5 080/125: ConvertType working except for text inputs and outputs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4d853b4 112/125: Old neighbors.h header removed, used new in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 9959e45 004/125: Arithmetic using new data structure,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 0a2fd0e 089/125: New implementations for histogram and CFP, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 213d1c3 106/125: gal_data_copy_* functions can now work on tiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4f6a8c4 095/125: Fits prints HDU info of FITS file with no options, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 63043e5 105/125: Statistical functions can accept tiled input, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 4839a11 104/125: Preparations for tile interpolation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master d6aa57f 081/125: Text library can write images also, besides tables, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 67cc39d 107/125: Structure to keep tessellation information, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master b03f7f9 113/125: New inplementation of interpolation complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 5fa37eb 096/125: HDU modification now fully implemented in Fits program, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 38f401b 121/125: Clump S/N threshold is now calculated on Sky clumps, Mohammad Akhlaghi, 2017/04/23