[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 42e04b6 002/113: Recent changes in master merg
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 42e04b6 002/113: Recent changes in master merged |
|
Date: |
Fri, 16 Apr 2021 10:33:29 -0400 (EDT) |
branch: master
commit 42e04b65c3615a05a02b9fd44a4832e43e1c097c
Merge: 3df6490 a84aa87
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>
Recent changes in master merged
Several bugs were fixed and one feature was added on master since this
branch was created. So with this commit, we are merging them into this
branch so it better corresponds to master.
---
bin/arithmetic/args.h | 14 +++
bin/arithmetic/main.h | 5 +-
bin/arithmetic/operands.c | 6 +-
bin/arithmetic/ui.c | 25 ++--
bin/arithmetic/ui.h | 16 +++
bin/cosmiccal/ui.c | 16 +--
bin/crop/onecrop.c | 2 +-
bin/crop/ui.c | 3 +-
bin/fits/fits.c | 2 +-
bin/mkcatalog/ui.c | 8 +-
doc/gnuastro.texi | 26 ++++-
doc/release-checklist.txt | 22 +++-
lib/arithmetic-binary.c | 32 ++---
lib/arithmetic-onlyint.c | 31 ++---
lib/arithmetic.c | 103 +++++++++-------
lib/data.c | 10 +-
lib/fits.c | 253 ++++++++++++++++++++++++++++++++--------
lib/gnuastro/fits.h | 3 +-
lib/wcs.c | 56 ++++++---
tests/arithmetic/or.sh | 2 +-
tests/mkcatalog/aperturephot.sh | 2 +-
21 files changed, 469 insertions(+), 168 deletions(-)
diff --git a/bin/arithmetic/args.h b/bin/arithmetic/args.h
index f6f73a8..2bf0e25 100644
--- a/bin/arithmetic/args.h
+++ b/bin/arithmetic/args.h
@@ -28,6 +28,20 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Definition of program-specific options. */
struct argp_option program_options[] =
{
+ {
+ "globalhdu",
+ UI_KEY_GLOBALHDU,
+ "STR",
+ 0,
+ "Use this HDU for all inputs, ignore `--hdu'.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->globalhdu,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+
{0}
};
diff --git a/bin/arithmetic/main.h b/bin/arithmetic/main.h
index 78ac2fb..f55451c 100644
--- a/bin/arithmetic/main.h
+++ b/bin/arithmetic/main.h
@@ -68,11 +68,12 @@ struct arithmeticparams
struct gal_options_common_params cp; /* Common parameters. */
/* Input: */
- gal_list_str_t *hdus; /* List of all given HDU strings. */
- gal_list_str_t *tokens; /* List of all arithmetic tokens. */
+ gal_list_str_t *hdus; /* List of all given HDU strings. */
+ gal_list_str_t *tokens; /* List of all arithmetic tokens. */
size_t addcounter; /* The number of FITS images added. */
size_t popcounter; /* The number of FITS images popped. */
gal_data_t refdata; /* Container for information of the data. */
+ char *globalhdu; /* Single HDU for all inputs. */
/* Operating mode: */
diff --git a/bin/arithmetic/operands.c b/bin/arithmetic/operands.c
index 3d6b77a..c069f34 100644
--- a/bin/arithmetic/operands.c
+++ b/bin/arithmetic/operands.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro-internal/checkset.h>
#include "main.h"
@@ -79,7 +80,10 @@ operands_add(struct arithmeticparams *p, char *filename,
gal_data_t *data)
if(filename != NULL && gal_fits_name_is_fits(filename))
{
/* Set the HDU for this filename. */
- newnode->hdu=gal_list_str_pop(&p->hdus);
+ if(p->globalhdu)
+ gal_checkset_allocate_copy(p->globalhdu, &newnode->hdu);
+ else
+ newnode->hdu=gal_list_str_pop(&p->hdus);
/* Increment the FITS counter. */
++p->addcounter;
diff --git a/bin/arithmetic/ui.c b/bin/arithmetic/ui.c
index ce995e6..e3cbac3 100644
--- a/bin/arithmetic/ui.c
+++ b/bin/arithmetic/ui.c
@@ -272,14 +272,21 @@ ui_check_options_and_arguments(struct arithmeticparams *p)
token->v[0]='-';
}
- /* Count the number of HDU values and check if its not less than the
- number of input FITS images. */
- for(hdu=p->hdus; hdu!=NULL; hdu=hdu->next) ++numhdus;
- if(numhdus<numfits)
- error(EXIT_FAILURE, 0, "not enough HDUs. There are %zu input FITS "
- "files, but only %zu HDUs. You can use the `--hdu' (`-h') option "
- "to specify the number or name of a HDU for each FITS file",
- numfits, numhdus);
+ /* Count the number of HDU values (if globalhdu isn't given) and check if
+ its not less than the number of input FITS images. */
+ if(p->globalhdu)
+ { if(p->hdus) { gal_list_str_free(p->hdus, 1); p->hdus=NULL; }; }
+ else
+ {
+ for(hdu=p->hdus; hdu!=NULL; hdu=hdu->next) ++numhdus;
+ if(numhdus<numfits)
+ error(EXIT_FAILURE, 0, "not enough HDUs. There are %zu input FITS "
+ "files, so the `--hdu' (`-h') option must be called %zu "
+ "times (once for each FITS file). If the HDU value is the "
+ "same for all the files, you may use `--globalhdu' (`-g') to "
+ "specify a single HDU to be used for any number of input "
+ "files", numfits, numfits);
+ }
}
@@ -382,7 +389,9 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
arithmeticparams *p)
void
freeandreport(struct arithmeticparams *p, struct timeval *t1)
{
+ /* Free the simple strings. */
free(p->cp.output);
+ if(p->globalhdu) free(p->globalhdu);
/* If there are any remaining HDUs in the hdus linked list, then
free them. */
diff --git a/bin/arithmetic/ui.h b/bin/arithmetic/ui.h
index 7ef7468..63ecc17 100644
--- a/bin/arithmetic/ui.h
+++ b/bin/arithmetic/ui.h
@@ -23,6 +23,22 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef UI_H
#define UI_H
+
+/* Available letters for short options:
+
+ a b c d e f i j k l m n p r s t u v w x y z
+ A B C E G H J L O Q R W X Y
+*/
+enum option_keys_enum
+{
+ /* With short-option version. */
+ UI_KEY_GLOBALHDU = 'g',
+
+ /* Only with long version (start with a value 1000, the rest will be set
+ automatically). */
+};
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct arithmeticparams *p);
diff --git a/bin/cosmiccal/ui.c b/bin/cosmiccal/ui.c
index c2864de..a7d8d98 100644
--- a/bin/cosmiccal/ui.c
+++ b/bin/cosmiccal/ui.c
@@ -195,16 +195,18 @@ parse_opt(int key, char *arg, struct argp_state *state)
static void
ui_read_check_only_options(struct cosmiccalparams *p)
{
- /* Check if the density fractions add up to 1. */
- if( (p->olambda + p->omatter + p->oradiation) != 1.0f )
- error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, but %f. "
+ double sum=p->olambda + p->omatter + p->oradiation;
+
+ /* Check if the density fractions add up to 1 (within floating point
+ error). */
+ if( sum > (1+1e-8) || sum < (1-1e-8) )
+ error(EXIT_FAILURE, 0, "sum of fractional densities is not 1, but %.8f. "
"The cosmological constant (`olambda'), matter (`omatter') "
- "and radiation (`oradiation') densities are given as %f, %f, %f",
- p->olambda + p->omatter + p->oradiation, p->olambda, p->omatter,
- p->oradiation);
+ "and radiation (`oradiation') densities are given as %.8f, %.8f, "
+ "%.8f", sum, p->olambda, p->omatter, p->oradiation);
/* The curvature fractional density: */
- p->ocurv=1-(p->olambda+p->omatter+p->oradiation);
+ p->ocurv=1-sum;
/* Convert H0 from km/sec/Mpc to 1/sec: */
p->H0s=p->H0/1000/GSL_CONST_MKSA_PARSEC;
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index 46ae330..64c9ada 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -838,7 +838,7 @@ onecrop_center_filled(struct onecropparams *crp)
if(checkcenter==0) return GAL_BLANK_UINT8;
/* Get the final size of the output image. */
- gal_fits_img_info(ofp, &type, &ndim, &dsize);
+ gal_fits_img_info(ofp, &type, &ndim, &dsize, NULL, NULL);
naxes[0]=dsize[1];
naxes[1]=dsize[0];
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index e766bfa..d60ea37 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -831,7 +831,8 @@ ui_preparations(struct cropparams *p)
img=&p->imgs[--input_counter];
img->name=gal_list_str_pop(&p->inputs);
tmpfits=gal_fits_hdu_open_format(img->name, p->cp.hdu, 0);
- gal_fits_img_info(tmpfits, &p->type, &img->ndim, &img->dsize);
+ gal_fits_img_info(tmpfits, &p->type, &img->ndim, &img->dsize,
+ NULL, NULL);
img->wcs=gal_wcs_read_fitsptr(tmpfits, p->hstartwcs, p->hendwcs,
&img->nwcs);
if(img->wcs)
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index 4045b5c..3af2ce5 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -138,7 +138,7 @@ fits_print_extension_info(struct fitsparams *p)
switch(hdutype)
{
case IMAGE_HDU:
- gal_fits_img_info(fptr, &type, &ndim, &dsize);
+ gal_fits_img_info(fptr, &type, &ndim, &dsize, NULL, NULL);
tstr=gal_type_name(type , 1);
break;
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 7096feb..87f95a1 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -500,6 +500,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
static void
ui_preparations_read_keywords(struct mkcatalogparams *p)
{
+ char *msg;
gal_data_t *tmp;
gal_data_t *keys=gal_data_array_calloc(2);
char *stdfile=p->stdfile ? p->stdfile : p->inputname;
@@ -575,8 +576,11 @@ ui_preparations_read_keywords(struct mkcatalogparams *p)
gal_fits_key_read(clumpsfile, p->clumpshdu, keys, 0, 0);
if(keys[0].status) p->clumpsn=NAN;
if(keys[1].status)
- error(EXIT_FAILURE, 0, "couldn't find NCLUMPS in the header of "
- "%s (hdu: %s).", p->clumpsfile, p->clumpshdu);
+ {
+ asprintf(&msg, "couldn't find/read NUMLABS in the header of "
+ "%s (hdu: %s), see error above", clumpsfile, p->clumpshdu);
+ gal_fits_io_error(keys[1].status, msg);
+ }
}
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 34c83b6..ae2b02a 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -4251,6 +4251,9 @@ can accept specific formats, for example
@ref{ConvertType} also accepts
@file{.fits}: The standard file name ending of a FITS image.
@item
+@file{.fit}: Alternative (3 character) FITS suffix.
+
+@item
@file{.fits.Z}: A FITS image compressed with @command{compress}.
@item
@@ -8896,6 +8899,13 @@ extra HDUs in the configuration files and by default
several HDUs with
a value of @option{0} are kept in the system-wide configuration file
when you install Gnuastro.
+@item -g INT/STR
+@itemx --globalhdu INT/STR
+Use the value to this option as the HDU of all input FITS files. This
+option is very convenient when you have many input files and the dataset of
+interest is in the same HDU of all the files. When this option is called,
+any values given to the @option{--hdu} option (explained above) are ignored
+and will not be used.
@end table
Arithmetic accepts two kinds of input: images and numbers. Images are
@@ -19033,8 +19043,9 @@ default string after the CFITSIO error.
@deftypefun int gal_fits_name_is_fits (char @code{*name})
If the @code{name} is an acceptable CFITSIO FITS filename return @code{1}
(one), otherwise return @code{0} (zero). The currently acceptable FITS
-suffixes are @file{.fits}, @file{.fits.gz}, @file{.fits.Z}, @file{.imh},
-@file{.fits.fz}. IMH is the IRAF format which is acceptable to CFITSIO.
+suffixes are @file{.fits}, @file{.fit}, @file{.fits.gz}, @file{.fits.Z},
+@file{.imh}, @file{.fits.fz}. IMH is the IRAF format which is acceptable to
+CFITSIO.
@end deftypefun
@deftypefun int gal_fits_suffix_is_fits (char @code{*suffix})
@@ -19293,7 +19304,7 @@ filename and HDU as input instead of an already opened
CFITSIO
@deftypefun void gal_fits_key_list_add (gal_fits_list_key_t @code{**list},
uint8_t @code{type}, char @code{*keyname}, int @code{kfree}, void
@code{*value}, int @code{vfree}, char @code{*comment}, int @code{cfree}, char
@code{*unit})
-Add on keyword to the top of list of header keywords that need to be
+Add a keyword to the top of list of header keywords that need to be
written into 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 (hence
the presence of the arguments ending in @code{free}). If the space for the
@@ -19359,10 +19370,13 @@ formats that is stored in FITS files. Only one image
may be stored in each
FITS HDU/extension. The functions described here can be used to get the
information of, read, or write images in FITS files.
-@deftypefun void gal_fits_img_info (fitsfile @code{*fptr}, int @code{*type},
size_t @code{*ndim}, size_t @code{**dsize})
+@deftypefun void gal_fits_img_info (fitsfile @code{*fptr}, int @code{*type},
size_t @code{*ndim}, size_t @code{**dsize}, char @code{**name}, char
@code{**unit})
Read the type (see @ref{Library data types}), number of dimensions, and
-size of the array along each dimension of the CFITSIO @code{fitsfile} into
-the @code{type}, @code{ndim}, and @code{dsize} pointers respectively.
+size along each dimension of the CFITSIO @code{fitsfile} into the
+@code{type}, @code{ndim}, and @code{dsize} pointers respectively. If
+@code{name} and @code{unit} are not @code{NULL} (point to a @code{char *}),
+then if the image has a name and units, the respective string will be put
+in these pointers.
@end deftypefun
@deftypefun {gal_data_t *} gal_fits_img_read (char @code{*filename}, char
@code{*hdu}, size_t @code{minmapsize})
diff --git a/doc/release-checklist.txt b/doc/release-checklist.txt
index 1e6b130..2a45564 100644
--- a/doc/release-checklist.txt
+++ b/doc/release-checklist.txt
@@ -130,6 +130,12 @@ Steps necessary to Package Gnuastro for Debian.
$ git clone git://anonscm.debian.org/debian-astro/packages/gnuastro.git
+ - Pull any possible changes that already exist:
+
+ $ cd gnuastro
+ $ git pull
+ $ cd ..
+
- A `gnuastro' directory will be built, for the moment don't go in it,
we'll keep all temporary files in this parent directory. Clean any
existing ones (if this process was done earlier):
@@ -185,10 +191,18 @@ Steps necessary to Package Gnuastro for Debian.
$ tar xf ../gnuastro_$ver.orig.tar.gz --strip-components=1
- - Update the ChangeLog (similar to previous entries):
+ - Update the ChangeLog (similar to previous entries).
+
+ IMPORTANT: if this is not a final release (just for making sure Gnuastro
+ builds on all systems), change `unstable' to `experimental'.
$ emacs debian/changelog
+ - Check the current Debian policy version at
+ https://www.debian.org/doc/debian-policy/ and update debian/control:
+
+ $ emacs debian/control
+
- Update your version of `pbuilder':
$ sudo pbuilder update
@@ -209,4 +223,8 @@ Steps necessary to Package Gnuastro for Debian.
$ git status # Add anything that remains.
$ git commit -m "Gnuastro $ver"
- - Push all the changes to the repository and announce on Debian Astro.
+ - Push all the changes to the repository:
+
+ $ git push --all origin
+
+ - Inform Debian Astro.
diff --git a/lib/arithmetic-binary.c b/lib/arithmetic-binary.c
index db2e2f2..02d7ce1 100644
--- a/lib/arithmetic-binary.c
+++ b/lib/arithmetic-binary.c
@@ -359,6 +359,7 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t
*lo,
l=gal_arithmetic_convert_to_compiled_type(lo, flags);
r=gal_arithmetic_convert_to_compiled_type(ro, flags);
+
/* Set the output type. For the comparison operators, the output type is
either 0 or 1, so we will set the output type to `unsigned char' for
efficient memory and CPU usage. Since the number of operators without
@@ -367,6 +368,7 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t
*lo,
operatrs are given, it will be chosen based on the input types.*/
otype=gal_arithmetic_binary_out_type(operator, l, r);
+
/* Set the output sizes. */
minmapsize = ( l->minmapsize < r->minmapsize
? l->minmapsize : r->minmapsize );
@@ -415,19 +417,6 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t
*lo,
}
- /* The type of the output dataset (`o->type') was chosen from `l' and `r'
- (copies of the orignal operands but in a compiled type, not
- necessarily the original `lo' and `ro' data structures). So we need to
- to get the final type based on the original operands and check if the
- final output needs changing. */
- if( o->type != final_otype )
- {
- tmp_o=gal_data_copy_to_new_type(o, final_otype);
- gal_data_free(o);
- o=tmp_o;
- }
-
-
/* Clean up. Note that if the input arrays can be freed, and any of right
or left arrays needed conversion, `BINARY_CONVERT_TO_COMPILED_TYPE'
has already freed the input arrays, so only `r' and `l' need
@@ -447,6 +436,23 @@ arithmetic_binary(int operator, uint8_t flags, gal_data_t
*lo,
if(r!=ro) gal_data_free(r);
}
+
+ /* The type of the output dataset (`o->type') was chosen from `l' and `r'
+ (copies of the orignal operands but in a compiled type, not
+ necessarily the original `lo' and `ro' data structures). So we need to
+ to get the final type based on the original operands and check if the
+ final output needs changing.
+
+ IMPORTANT: This has to be done after (possibly) freeing the left and
+ right operands because this step can change the `o' pointer which
+ they may depend on (when working in-place). */
+ if( o->type != final_otype )
+ {
+ tmp_o=gal_data_copy_to_new_type(o, final_otype);
+ gal_data_free(o);
+ o=tmp_o;
+ }
+
/* Return */
return o;
}
diff --git a/lib/arithmetic-onlyint.c b/lib/arithmetic-onlyint.c
index 564c849..9e26382 100644
--- a/lib/arithmetic-onlyint.c
+++ b/lib/arithmetic-onlyint.c
@@ -368,19 +368,6 @@ arithmetic_onlyint_binary(int operator, unsigned char
flags,
}
- /* The type of the output dataset (`o->type') was chosen from `l' and `r'
- (copies of the orignal operands but in a compiled type, not
- necessarily the original `lo' and `ro' data structures). So we need to
- to get the final type based on the original operands and check if the
- final output needs changing. */
- if( o->type != final_otype )
- {
- tmp_o=gal_data_copy_to_new_type(o, final_otype);
- gal_data_free(o);
- o=tmp_o;
- }
-
-
/* Clean up. Note that if the input arrays can be freed, and any of right
or left arrays needed conversion, `BINOIN_CONVERT_TO_COMPILED_TYPE'
has already freed the input arrays, so only `r' and `l' need
@@ -400,6 +387,24 @@ arithmetic_onlyint_binary(int operator, unsigned char
flags,
if(r!=ro) gal_data_free(r);
}
+
+ /* The type of the output dataset (`o->type') was chosen from `l' and `r'
+ (copies of the orignal operands but in a compiled type, not
+ necessarily the original `lo' and `ro' data structures). So we need to
+ to get the final type based on the original operands and check if the
+ final output needs changing.
+
+ IMPORTANT: This has to be done after (possibly) freeing the left and
+ right operands because this step can change the `o' pointer which
+ they may depend on (when working in-place). */
+ if( o->type != final_otype )
+ {
+ tmp_o=gal_data_copy_to_new_type(o, final_otype);
+ gal_data_free(o);
+ o=tmp_o;
+ }
+
+
/* Return */
return o;
}
diff --git a/lib/arithmetic.c b/lib/arithmetic.c
index 7fe313b..6e20539 100644
--- a/lib/arithmetic.c
+++ b/lib/arithmetic.c
@@ -673,18 +673,23 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/***********************************************************************/
#define MULTIOPERAND_MIN(TYPE) { \
TYPE p, max; \
+ size_t n, j=0; \
gal_type_max(list->type, &max); \
do /* Loop over each pixel */ \
{ \
+ n=0; \
p=max; \
for(i=0;i<dnum;++i) /* Loop over each array. */ \
- { /* Only for integer types, does b==b. */ \
- if(hasblank[i] && b==b) \
- { if( *a[i] != b ) p = *a[i] < p ? *a[i] : p; \
- else p = *a[i] < p ? *a[i] : p; } \
- ++a[i]; \
+ { /* Only for integer types, b==b. */ \
+ if( hasblank[i] && b==b) \
+ { \
+ if( a[i][j] != b ) \
+ { p = a[i][j] < p ? a[i][j] : p; ++n; } \
+ } \
+ else { p = a[i][j] < p ? a[i][j] : p; ++n; } \
} \
- *o++=p; \
+ *o++ = n ? p : b; /* No usable elements: set to blank. */ \
+ ++j; \
} \
while(o<of); \
}
@@ -695,18 +700,23 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_MAX(TYPE) { \
TYPE p, min; \
+ size_t n, j=0; \
gal_type_min(list->type, &min); \
do /* Loop over each pixel */ \
{ \
+ n=0; \
p=min; \
for(i=0;i<dnum;++i) /* Loop over each array. */ \
- { /* Only for integer types, does b==b. */ \
- if(hasblank[i] && b==b) \
- { if( *a[i] != b ) p = *a[i] > p ? *a[i] : p; \
- else p = *a[i] > p ? *a[i] : p; } \
- ++a[i]; \
+ { /* Only for integer types, b==b. */ \
+ if( hasblank[i] && b==b) \
+ { \
+ if( a[i][j] != b ) \
+ { p = a[i][j] > p ? a[i][j] : p; ++n; } \
+ } \
+ else { p = a[i][j] > p ? a[i][j] : p; ++n; } \
} \
- *o++=p; \
+ *o++ = n ? p : b; /* No usable elements: set to blank. */ \
+ ++j; \
} \
while(o<of); \
}
@@ -716,7 +726,8 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_NUM { \
- int n, use; \
+ int use; \
+ size_t n, j=0; \
do /* Loop over each pixel */ \
{ \
n=0; \
@@ -725,14 +736,15 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/* Only integers and non-NaN floats: v==v is 1. */ \
if(hasblank[i]) \
use = ( b==b \
- ? ( *a[i]!=b ? 1 : 0 ) /* Integer */ \
- : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
+ : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
else use=1; \
\
- /* a[i] must be incremented to next pixel in any case. */ \
- if(use) ++n; else ++a[i]; \
+ /* Increment counter if necessary. */ \
+ if(use) ++n; \
} \
*o++ = n; \
+ ++j; \
} \
while(o<of); \
}
@@ -742,8 +754,9 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_SUM { \
+ int use; \
double sum; \
- int n, use; \
+ size_t n, j=0; \
do /* Loop over each pixel */ \
{ \
n=0; \
@@ -753,14 +766,15 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/* Only integers and non-NaN floats: v==v is 1. */ \
if(hasblank[i]) \
use = ( b==b \
- ? ( *a[i]!=b ? 1 : 0 ) /* Integer */ \
- : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
+ : ( a[i][j]==*a[i] ? 1 : 0 ) ); /* Float */ \
else use=1; \
\
- /* a[i] must be incremented to next pixel in any case. */ \
- if(use) { sum += *a[i]++; ++n; } else ++a[i]; \
+ /* Use in sum if necessary. */ \
+ if(use) { sum += a[i][j]; ++n; } \
} \
*o++ = n ? sum : b; \
+ ++j; \
} \
while(o<of); \
}
@@ -770,8 +784,9 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_MEAN { \
+ int use; \
double sum; \
- int n, use; \
+ size_t n, j=0; \
do /* Loop over each pixel */ \
{ \
n=0; \
@@ -781,14 +796,15 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/* Only integers and non-NaN floats: v==v is 1. */ \
if(hasblank[i]) \
use = ( b==b \
- ? ( *a[i]!=b ? 1 : 0 ) /* Integer */ \
- : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
+ : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
else use=1; \
\
- /* a[i] must be incremented to next pixel in any case. */ \
- if(use) { sum += *a[i]++; ++n; } else ++a[i]; \
+ /* Calculate the mean if necessary. */ \
+ if(use) { sum += a[i][j]; ++n; } \
} \
*o++ = n ? sum/n : b; \
+ ++j; \
} \
while(o<of); \
}
@@ -798,7 +814,8 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_STD { \
- int n, use; \
+ int use; \
+ size_t n, j=0; \
double sum, sum2; \
do /* Loop over each pixel */ \
{ \
@@ -809,20 +826,20 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/* Only integers and non-NaN floats: v==v is 1. */ \
if(hasblank[i]) \
use = ( b==b \
- ? ( *a[i]!=b ? 1 : 0 ) /* Integer */ \
- : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
+ : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
else use=1; \
\
- /* a[i] must be incremented to next pixel in any case. */ \
+ /* Calculate the necessary parameters if necessary. */ \
if(use) \
{ \
- sum2 = *a[i] * *a[i]; \
- sum += *a[i]++; \
+ sum2 += a[i][j] * a[i][j]; \
+ sum += a[i][j]; \
++n; \
} \
- else ++a[i]; \
} \
*o++ = n ? sqrt( (sum2-sum*sum/n)/n ) : b; \
+ ++j; \
} \
while(o<of); \
}
@@ -832,12 +849,14 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_MEDIAN(TYPE, QSORT_F) { \
- int n, use; \
+ int use; \
+ size_t n, j=0; \
TYPE *pixs=gal_data_malloc_array(list->type, dnum, __func__, "pixs"); \
\
/* Loop over each pixel */ \
do \
{ \
+ /* Initialize. */ \
n=0; \
\
/* Loop over each array. */ \
@@ -846,12 +865,12 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
/* Only integers and non-NaN floats: v==v is 1. */ \
if(hasblank[i]) \
use = ( b==b \
- ? ( *a[i]!=b ? 1 : 0 ) /* Integer */ \
- : ( *a[i]==*a[i] ? 1 : 0 ) ); /* Float */ \
+ ? ( a[i][j]!=b ? 1 : 0 ) /* Integer */ \
+ : ( a[i][j]==a[i][j] ? 1 : 0 ) ); /* Float */ \
else use=1; \
\
- /* a[i] must be incremented to next pixel in any case. */ \
- if(use) pixs[n++]=*a[i]++; else ++a[i]; \
+ /* Put the value into the array of values. */ \
+ if(use) pixs[n++]=a[i][j]; \
} \
\
/* Sort all the values for this pixel and return the median. */ \
@@ -862,6 +881,7 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
} \
else \
*o++=b; \
+ ++j; \
} \
while(o<of); \
\
@@ -874,7 +894,8 @@ arithmetic_where(unsigned char flags, gal_data_t *out,
gal_data_t *cond,
#define MULTIOPERAND_TYPE_SET(TYPE, QSORT_F) { \
- TYPE b, **a, *o=out->array, *of=out->array+out->size; \
+ TYPE b, **a, *o=out->array, *of=o+out->size; \
+ size_t i=0; /* Different from the `i' in the main function. */ \
\
/* Allocate space to keep the pointers to the arrays of each. */ \
/* Input data structure. The operators will increment these */ \
@@ -1517,7 +1538,7 @@ gal_arithmetic(int operator, unsigned char flags, ...)
break;
- /* When operator is not recognized. */
+ /* When the operator is not recognized. */
default:
error(EXIT_FAILURE, 0, "%s: the argument \"%d\" could not be "
"interpretted as an operator", __func__, operator);
diff --git a/lib/data.c b/lib/data.c
index 0b1589d..68178b2 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -409,11 +409,11 @@ gal_data_free_contents(gal_data_t *data)
"`gal_data_free_contents' was a NULL pointer", __func__);
/* Free all the possible allocations. */
- if(data->name) { free(data->name); data->name=NULL; }
- if(data->unit) { free(data->unit); data->unit=NULL; }
- if(data->dsize) { free(data->dsize); data->dsize=NULL; }
- if(data->wcs) { wcsfree(data->wcs); data->wcs=NULL; }
- if(data->comment) { free(data->comment); data->comment=NULL; }
+ if(data->name) { free(data->name); data->name = NULL; }
+ if(data->unit) { free(data->unit); data->unit = NULL; }
+ if(data->dsize) { free(data->dsize); data->dsize = NULL; }
+ if(data->wcs) { wcsfree(data->wcs); data->wcs = NULL; }
+ if(data->comment) { free(data->comment); data->comment = NULL; }
/* If the data type is string, then each element in the array is actually
a pointer to the array of characters, so free them before freeing the
diff --git a/lib/fits.c b/lib/fits.c
index 0f42120..abb30cd 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -102,7 +102,8 @@ gal_fits_name_is_fits(char *name)
{
size_t len;
len=strlen(name);
- if ( ( len>=4 && strcmp(&name[len-4], "fits") == 0 )
+ if ( ( len>=3 && strcmp(&name[len-3], "fit") == 0 )
+ || ( len>=4 && strcmp(&name[len-4], "fits") == 0 )
|| ( len>=7 && strcmp(&name[len-7], "fits.gz") == 0 )
|| ( len>=6 && strcmp(&name[len-6], "fits.Z") == 0 )
|| ( len>=3 && strcmp(&name[len-3], "imh") == 0 )
@@ -123,10 +124,11 @@ gal_fits_name_is_fits(char *name)
int
gal_fits_suffix_is_fits(char *suffix)
{
- if (strcmp(suffix, "fits") == 0 || strcmp(suffix, ".fits") == 0
+ if (strcmp(suffix, "fit") == 0 || strcmp(suffix, ".fit") == 0
+ || strcmp(suffix, "fits") == 0 || strcmp(suffix, ".fits") == 0
|| strcmp(suffix, "fits.gz") == 0 || strcmp(suffix, ".fits.gz") == 0
- || strcmp(suffix, "fits.Z") == 0 || strcmp(suffix, ".fits.Z") == 0
- || strcmp(suffix, "imh") == 0 || strcmp(suffix, ".imh") == 0
+ || strcmp(suffix, "fits.Z") == 0 || strcmp(suffix, ".fits.Z") == 0
+ || strcmp(suffix, "imh") == 0 || strcmp(suffix, ".imh") == 0
|| strcmp(suffix, "fits.fz") == 0 || strcmp(suffix, ".fits.fz") == 0)
return 1;
else
@@ -308,17 +310,22 @@ gal_fits_type_to_datatype(uint8_t type)
else if( sizeof(int) == w ) return TINT;
break;
+ /* On 32-bit systems, the length of `int' and `long' are both
+ 32-bits. But CFITSIO's LONG type is preferred because it is designed
+ to be 32-bit. Its `INT' type is not clearly defined and caused
+ problems when reading keywords.*/
case GAL_TYPE_UINT32:
w=4;
- if ( sizeof(int) == w ) return TUINT;
- else if( sizeof(long) == w ) return TULONG;
+ if ( sizeof(long) == w ) return TULONG;
+ else if( sizeof(int) == w ) return TUINT;
else if( sizeof(short) == w ) return TUSHORT;
break;
+ /* Similar to UINT32 above. */
case GAL_TYPE_INT32:
w=4;
- if ( sizeof(int) == w ) return TINT;
- else if( sizeof(long) == w ) return TLONG;
+ if ( sizeof(long) == w ) return TLONG;
+ else if( sizeof(int) == w ) return TINT;
else if( sizeof(short) == w ) return TSHORT;
break;
@@ -459,6 +466,55 @@ gal_fits_datatype_to_type(int datatype, int
is_table_column)
+/* When there is a BZERO or TZERO and BSCALE or TSCALE keywords, then the
+ type that must be used to store the actual values of the pixels may be
+ different from the type from BITPIX. This function does the necessary
+ corrections.*/
+void
+fits_type_correct(int *type, double bscale, double bzero)
+{
+ int tofloat=1;
+
+ /* Work based on type, for the default conversions defined by CFITSIO,
+ make the proper correction, otherwise set the type to float. */
+ if(bscale==1.0f)
+ switch(*type)
+ {
+ case GAL_TYPE_UINT8:
+ if(bzero == -128.0f) { *type = GAL_TYPE_INT8; tofloat=0; }
+ break;
+
+ case GAL_TYPE_INT16:
+ /* Defined by the FITS standard: */
+ if(bzero == 32768) { *type = GAL_TYPE_UINT16; tofloat=0; }
+ break;
+
+ case GAL_TYPE_INT32:
+ /* Defined by the FITS standard: */
+ if(bzero == 2147483648) { *type = GAL_TYPE_UINT32; tofloat=0; }
+ break;
+
+ case GAL_TYPE_INT64:
+ /* Defined by the FITS standard: */
+ if(bzero == 9223372036854775808.0f)
+ {*type = GAL_TYPE_UINT64; tofloat=0;}
+ break;
+
+ /* For the other types (when `BSCALE=1.0f'), currently no correction is
+ necessary, maybe later we can check if the scales are integers and
+ set the integer output type to the smallest type that can allow the
+ scaled values. */
+ default: tofloat=0;
+ }
+
+ /* If the type must be a float, then do the conversion. */
+ if(tofloat) *type=GAL_TYPE_FLOAT32;
+}
+
+
+
+
+
@@ -1226,14 +1282,19 @@ gal_fits_key_write_version(fitsfile *fptr,
gal_fits_list_key_t *headers,
/* Note that the FITS standard defines any array as an `image',
irrespective of how many dimensions it has. This function will return
- the Gnuastro-type and also the number of dimensions and size along each
- dimension of the image. Note that `*dsize' will be allocated here, so it
- must not point to any already allocated space. */
+ the Gnuastro-type, the number of dimensions and size along each
+ dimension of the image along with its name and units if necessary (not
+ NULL). Note that `*dsize' will be allocated here, so it must not point
+ to any already allocated space. */
void
-gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize)
+gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize,
+ char **name, char **unit)
{
- size_t i;
+ char **str;
+ size_t i, dsize_key=1;
int bitpix, status=0, naxis;
+ double bzero=NAN, bscale=NAN;
+ gal_data_t *key, *keysll=NULL;
long naxes[GAL_FITS_MAX_NDIM];
/* Get the BITPIX, number of dimensions and size of each dimension. */
@@ -1245,11 +1306,57 @@ gal_fits_img_info(fitsfile *fptr, int *type, size_t
*ndim, size_t **dsize)
/* Convert bitpix to Gnuastro's known types. */
*type=gal_fits_bitpix_to_type(bitpix);
+ /* Define the names of the possibly existing important keywords about the
+ dataset. We are defining these in the opposite order to be read by
+ CFITSIO. The way Gnuastro writes the FITS keywords, the output will
+ first have `BZERO', then `BSCALE', then `EXTNAME', then, `BUNIT'.*/
+ gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_STRING, 1, &dsize_key,
+ NULL, 0, -1, "BUNIT", NULL, NULL);
+ gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_STRING, 1, &dsize_key,
+ NULL, 0, -1, "EXTNAME", NULL, NULL);
+ gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_FLOAT64, 1, &dsize_key,
+ NULL, 0, -1, "BSCALE", NULL, NULL);
+ gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_FLOAT64, 1, &dsize_key,
+ NULL, 0, -1, "BZERO", NULL, NULL);
+ gal_fits_key_read_from_ptr(fptr, keysll, 0, 0);
+
+
+ /* Read the special keywords. */
+ i=1;
+ for(key=keysll;key!=NULL;key=key->next)
+ {
+ /* Recall that the order is the opposite (this is a last-in-first-out
+ list. */
+ if(key->status==0)
+ {
+ switch(i)
+ {
+ case 4: if(unit) {str = key->array; *unit = *str; *str=NULL;} break;
+ case 3: if(name) {str = key->array; *name = *str; *str=NULL;} break;
+ case 2: bscale = *(double *)(key->array); break;
+ case 1: bzero = *(double *)(key->array); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. For some reason, there are more "
+ "keywords requested ", __func__, PACKAGE_BUGREPORT);
+ }
+ }
+ ++i;
+ }
+
+ if( !isnan(bscale) || !isnan(bzero) )
+ fits_type_correct(type, bscale, bzero);
+
+
/* Allocate the array to keep the dimension size and fill it in, note
that its order is the opposite of naxes. */
*dsize=gal_data_malloc_array(GAL_TYPE_INT64, *ndim, __func__, "dsize");
for(i=0; i<*ndim; ++i)
(*dsize)[i]=naxes[*ndim-1-i];
+
+
+ /* Clean up. */
+ gal_list_data_free(keysll);
}
@@ -1261,13 +1368,12 @@ gal_data_t *
gal_fits_img_read(char *filename, char *hdu, size_t minmapsize)
{
void *blank;
- int anyblank;
long *fpixel;
fitsfile *fptr;
- int status=0, type;
- gal_data_t *img, *keysll=NULL;
- char **str, *name=NULL, *unit=NULL;
- size_t i, ndim, *dsize, dsize_key=1;
+ gal_data_t *img;
+ size_t i, ndim, *dsize;
+ char *name=NULL, *unit=NULL;
+ int status=0, type, anyblank;
/* Check HDU for realistic conditions: */
@@ -1275,7 +1381,7 @@ gal_fits_img_read(char *filename, char *hdu, size_t
minmapsize)
/* Get the info and allocate the data structure. */
- gal_fits_img_info(fptr, &type, &ndim, &dsize);
+ gal_fits_img_info(fptr, &type, &ndim, &dsize, &name, &unit);
/* Check if there is any dimensions (the first header can sometimes have
@@ -1290,30 +1396,20 @@ gal_fits_img_read(char *filename, char *hdu, size_t
minmapsize)
hdu);
- /* Set the fpixel array (first pixel in all dimensions): */
+ /* Set the fpixel array (first pixel in all dimensions). Note that the
+ `long' type will not be larger than 64-bits, so, we'll just assume it
+ is 64-bits for space allocation. On 32-bit systems, this won't be a
+ problem, the space will be written/read as 32-bit `long' any way,
+ we'll just have a few empty bytes that will be freed anyway at the end
+ of this function. */
fpixel=gal_data_malloc_array(GAL_TYPE_INT64, ndim, __func__, "fpixel");
for(i=0;i<ndim;++i) fpixel[i]=1;
- /* Read the possibly existing useful keywords. Note that the values are
- in allocated strings in the keys[i] data structures. Note that we need
- the linked list of keys to keep the `name' and `unit' pointers. We can
- free the linked list after `gal_data_alloc' has read/copied the
- values.*/
- gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_STRING, 1, &dsize_key,
- NULL, 0, -1, "EXTNAME", NULL, NULL);
- gal_list_data_add_alloc(&keysll, NULL, GAL_TYPE_STRING, 1, &dsize_key,
- NULL, 0, -1, "BUNIT", NULL, NULL);
- gal_fits_key_read_from_ptr(fptr, keysll, 0, 0);
- if(keysll->status==0) {str=keysll->array; unit=*str; }
- if(keysll->next->status==0) {str=keysll->next->array; name=*str; }
-
-
/* Allocate the space for the array and for the blank values. */
img=gal_data_alloc(NULL, type, ndim, dsize, NULL, 0, minmapsize,
name, unit, NULL);
blank=gal_blank_alloc_write(type);
- gal_list_data_free(keysll);
free(dsize);
@@ -1426,32 +1522,97 @@ fitsfile *
gal_fits_img_write_to_ptr(gal_data_t *input, char *filename)
{
void *blank;
- char *wcsstr;
+ int64_t *i64;
fitsfile *fptr;
+ uint64_t *u64, *u64f;
long fpixel=1, *naxes;
size_t i, ndim=input->ndim;
- gal_data_t *towrite, *block=gal_tile_block(input);
- int nkeyrec, status=0, datatype=gal_fits_type_to_datatype(block->type);
+ int nkeyrec, hasblank, status=0, datatype=0;
+ char *wcsstr, *u64key;
+ gal_data_t *i64data, *towrite, *block=gal_tile_block(input);
/* If the input is a tile (isn't a contiguous region of memory), then
copy it into a contiguous region. */
towrite = input==block ? input : gal_data_copy(input);
+ hasblank=gal_blank_present(towrite, 0);
/* Allocate the naxis area. */
naxes=gal_data_malloc_array( ( sizeof(long)==8
? GAL_TYPE_INT64
: GAL_TYPE_INT32 ), ndim, __func__, "naxes");
+
/* Open the file for writing */
fptr=gal_fits_open_to_write(filename);
+
/* Fill the `naxes' array (in opposite order, and `long' type): */
for(i=0;i<ndim;++i) naxes[ndim-1-i]=towrite->dsize[i];
- /* Create the FITS file. */
- fits_create_img(fptr, gal_fits_type_to_bitpix(towrite->type),
- ndim, naxes, &status);
- gal_fits_io_error(status, NULL);
+
+ /* Create the FITS file. Unfortunately CFITSIO doesn't have a macro for
+ UINT64, TLONGLONG is only for (signed) INT64. So if the dataset has
+ that type, we'll have to convert it to `INT64' and in the mean-time
+ shift its zero, we will then have to write the BZERO and BSCALE
+ keywords accordingly. */
+ if(block->type==GAL_TYPE_UINT64)
+ {
+ /* Allocate the necessary space. */
+ i64data=gal_data_alloc(NULL, GAL_TYPE_INT64, ndim, towrite->dsize,
+ NULL, 0, block->minmapsize, NULL, NULL, NULL);
+
+ /* Copy the values while making the conversion. */
+ i64=i64data->array;
+ u64f=(u64=towrite->array)+towrite->size;
+ if(hasblank)
+ {
+ do *i64++ = ( *u64==GAL_BLANK_UINT64
+ ? GAL_BLANK_INT64
+ : (*u64 + INT64_MIN) );
+ while(++u64<u64f);
+ }
+ else
+ do *i64++ = (*u64 + INT64_MIN); while(++u64<u64f);
+
+ /* We can now use CFITSIO's signed-int64 type macros. */
+ datatype=TLONGLONG;
+ fits_create_img(fptr, LONGLONG_IMG, ndim, naxes, &status);
+ gal_fits_io_error(status, NULL);
+
+
+ /* Write the image into the file. */
+ fits_write_img(fptr, datatype, fpixel, i64data->size, i64data->array,
+ &status);
+ gal_fits_io_error(status, NULL);
+
+
+ /* We need to write the BZERO and BSCALE keywords manually. VERY
+ IMPORTANT: this has to be done after writing the array. We cannot
+ write this huge integer as a variable, so we'll simply write the
+ full record/card. It is just important that the string be larger
+ than 80 characters, CFITSIO will trim the rest of the string. */
+ u64key="BZERO = 9223372036854775808 / Offset of data
";
+ fits_write_record(fptr, u64key, &status);
+ u64key="BSCALE = 1 / Default scaling factor
";
+ fits_write_record(fptr, u64key, &status);
+ gal_fits_io_error(status, NULL);
+ }
+ else
+ {
+ /* Set the datatype */
+ datatype=gal_fits_type_to_datatype(block->type);
+
+ /* Create the FITS file. */
+ fits_create_img(fptr, gal_fits_type_to_bitpix(towrite->type),
+ ndim, naxes, &status);
+ gal_fits_io_error(status, NULL);
+
+ /* Write the image into the file. */
+ fits_write_img(fptr, datatype, fpixel, towrite->size, towrite->array,
+ &status);
+ gal_fits_io_error(status, NULL);
+ }
+
/* Remove the two comment lines put by CFITSIO. Note that in some cases,
it might not exist. When this happens, the status value will be
@@ -1461,13 +1622,10 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char
*filename)
fits_delete_key(fptr, "COMMENT", &status);
status=0;
- /* Write the image into the file. */
- fits_write_img(fptr, datatype, fpixel, towrite->size, towrite->array,
- &status);
/* If we have blank pixels, we need to define a BLANK keyword when we are
dealing with integer types. */
- if(gal_blank_present(towrite, 0))
+ if(hasblank)
switch(towrite->type)
{
case GAL_TYPE_FLOAT32:
@@ -1484,14 +1642,17 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char
*filename)
free(blank);
}
+
/* Write the extension name to the header. */
if(towrite->name)
fits_write_key(fptr, TSTRING, "EXTNAME", towrite->name, "", &status);
+
/* Write the units to the header. */
if(towrite->unit)
fits_write_key(fptr, TSTRING, "BUNIT", towrite->unit, "", &status);
+
/* Write comments if they exist. */
if(towrite->comment)
fits_write_comment(fptr, towrite->comment, &status);
diff --git a/lib/gnuastro/fits.h b/lib/gnuastro/fits.h
index 05299b3..b7b67c8 100644
--- a/lib/gnuastro/fits.h
+++ b/lib/gnuastro/fits.h
@@ -202,7 +202,8 @@ gal_fits_key_write_version(fitsfile *fptr,
gal_fits_list_key_t *headers,
****************** Array functions *****************
*************************************************************/
void
-gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize);
+gal_fits_img_info(fitsfile *fptr, int *type, size_t *ndim, size_t **dsize,
+ char **name, char **unit);
gal_data_t *
gal_fits_img_read(char *filename, char *hdu, size_t minmapsize);
diff --git a/lib/wcs.c b/lib/wcs.c
index bb321bf..d27282c 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -70,8 +70,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs, size_t
hendwcs,
int *nwcs)
{
/* Declaratins: */
- struct wcsprm *wcs;
int nkeys=0, status=0;
+ struct wcsprm *wcs=NULL;
char *fullheader, *to, *from;
int relax = WCSHDR_all; /* Macro: use all informal WCS extensions. */
int ctrl = 0; /* Don't report why a keyword wasn't used. */
@@ -111,7 +111,8 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs,
size_t hendwcs,
/*******************************************************/
}
- /* WCSlib function */
+
+ /* WCSlib function to parse the FITS headers. */
status=wcspih(fullheader, nkeys, relax, ctrl, &nreject, nwcs, &wcs);
if(status)
{
@@ -125,22 +126,45 @@ gal_wcs_read_fitsptr(fitsfile *fptr, size_t hstartwcs,
size_t hendwcs,
gal_fits_io_error(status, "problem in fitsarrayvv.c for freeing "
"the memory used to keep all the headers");
+
/* Set the internal structure: */
- status=wcsset(wcs);
- if(status)
+ if(wcs)
{
- fprintf(stderr, "\n##################\n"
- "WCSLIB Warning: wcsset ERROR %d: %s.\n"
- "##################\n",
- status, wcs_errmsg[status]);
- wcs=NULL; *nwcs=0;
+ /* CTYPE is a mandatory WCS keyword, so if it hasn't been given (its
+ '\0'), then the headers didn't have a WCS structure. However,
+ WCSLIB still fills in the basic information (for example the
+ dimensionality of the dataset). */
+ if(wcs->ctype[0][0]=='\0')
+ {
+ wcsfree(wcs);
+ wcs=NULL;
+ *nwcs=0;
+ }
+ else
+ {
+ /* Set the WCS structure. */
+ status=wcsset(wcs);
+ if(status)
+ {
+ fprintf(stderr, "\n##################\n"
+ "WCSLIB Warning: wcsset ERROR %d: %s.\n"
+ "##################\n",
+ status, wcs_errmsg[status]);
+ wcsfree(wcs);
+ wcs=NULL;
+ *nwcs=0;
+ }
+ else
+ /* A correctly useful WCS is present. When no PC matrix
+ elements were present in the header, the default PC matrix
+ (a unity matrix) is used. In this case WCSLIB doesn't set
+ `altlin' (and gives it a value of 0). In Gnuastro, later on,
+ we might need to know the type of the matrix used, so in
+ such a case, we will set `altlin' to 1. */
+ if(wcs->altlin==0) wcs->altlin=1;
+ }
}
- /* Initialize the wcsprm struct
- if ((status = wcsset(*wcs)))
- error(EXIT_FAILURE, 0, "%s: wcsset ERROR %d: %s.\n", __func__,
- status, wcs_errmsg[status]);
- */
/* Return the WCS structure. */
return wcs;
@@ -294,13 +318,13 @@ gal_wcs_warp_matrix(struct wcsprm *wcs)
__func__, size*sizeof *out);
/* Fill in the array. */
- if(wcs->altlin & 0x1) /* Has a PCi_j array. */
+ if(wcs->altlin & 0x1) /* Has a PCi_j array. */
{
for(i=0;i<wcs->naxis;++i)
for(j=0;j<wcs->naxis;++j)
out[i*wcs->naxis+j] = wcs->cdelt[i] * wcs->pc[i*wcs->naxis+j];
}
- else if(wcs->altlin & 0x2) /* Has CDi_j array */
+ else if(wcs->altlin & 0x2) /* Has CDi_j array. */
{
for(i=0;i<size;++i)
out[i]=wcs->cd[i];
diff --git a/tests/arithmetic/or.sh b/tests/arithmetic/or.sh
index eb08426..263dcbf 100755
--- a/tests/arithmetic/or.sh
+++ b/tests/arithmetic/or.sh
@@ -49,4 +49,4 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# Actual test script
# ==================
-$execname $img 6 eq $img 3 eq or -h1 -h1 --output=or.fits
+$execname $img 6 eq $img 3 eq or -h2 -h2 --output=or.fits
diff --git a/tests/mkcatalog/aperturephot.sh b/tests/mkcatalog/aperturephot.sh
index 5004ae8..130be3f 100755
--- a/tests/mkcatalog/aperturephot.sh
+++ b/tests/mkcatalog/aperturephot.sh
@@ -53,4 +53,4 @@ if [ ! -f $objimg ]; then echo "$objimg does not exist";
exit 77; fi
# ==================
$execname $img --objectsfile=$objimg --objectshdu=1 \
--output=aperturephot.txt \
- --ids --x --y --ra --dec --magnitude --sn
+ --objid --x --y --ra --dec --magnitude --sn
- [gnuastro-commits] master updated (89db01d -> dd4d43e), Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 17d6f65 004/113: Merged recent corrections in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 42e04b6 002/113: Recent changes in master merged,
Mohammad Akhlaghi <=
- [gnuastro-commits] master f25b55e 005/113: Merged with recent changes in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master a11e71c 003/113: Merged with recent changes in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 3df6490 001/113: Crop works in 3D, except for --polygon, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master d14384a 007/113: Minor corrections in MakeProfiles manual, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 5abb468 008/113: Other minor corrections in the documentation, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 3e76e1f 010/113: MakeProfiles --kernel builds 3D kernels also, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 81f3f65 023/113: More --coordcol options acceptable in Crop, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e47f8db 024/113: Merged recent work in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8fa5ff1 026/113: Minor edit in book (part added in last commit), Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 40f0a56 013/113: Minor corrections to MakeProfiles continued, Mohammad Akhlaghi, 2021/04/16