[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master ea821eb 1/2: BZERO and BSCALE used to determin
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master ea821eb 1/2: BZERO and BSCALE used to determine type of input |
|
Date: |
Sun, 25 Jun 2017 08:54:57 -0400 (EDT) |
branch: master
commit ea821ebc8d2d2b2bed1f1d0b9ac4b49c8aaa7e12
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
BZERO and BSCALE used to determine type of input
According to the FITS standard, for images, BITPIX doesn't have any code
for `int8_t', `uint16_t', `uint32_t' and `uint64_t' types. So the FITS
standard allows for the `BZERO' and `BSCALE' keywords to be used to create
these types and CFITSIO does the proper conversion. But until now, we
hadn't taken these two keywords into account.
So for example, when `BITPIX=16' (FITS code for `int16'), `BZERO=32768' and
`BSCALE=1', the dataset is actually `uint16'. However, until this commit a
signed 16-bit image would be created to store the values. Hence, if the
image contained large values, we would get a CFITSIO numerical overflow
error (CFITSIO error code 412).
With this commit, the function that reads the image information (including
its numerical data type: `gal_fits_img_info') now also looks for `BZERO'
and `BSCALE' and if they are given and correspond to the standard type
conversions described above, the type of the image is corrected
appropriately. Since it involved reading keywords, to speed things up, it
also reads the name and units if given and so it takes two extra arguments.
In the process of playing with the `BZERO' and `BSCALE', Gnuastro can now
also write `uint64_t' datasets in FITS.
This issue was reported by David Valls-Gaboud.
Improvements:
- Programs can write unsigned 64-bit integer FITS files.
Changes:
- `gal_fits_img_info' now also returns the name and units of the dataset
(if they aren't NULL). So it takes two extra arguments.
---
bin/crop/onecrop.c | 2 +-
bin/crop/ui.c | 3 +-
bin/fits/fits.c | 2 +-
doc/gnuastro.texi | 11 ++-
lib/fits.c | 228 +++++++++++++++++++++++++++++++++++++++++++---------
lib/gnuastro/fits.h | 3 +-
6 files changed, 203 insertions(+), 46 deletions(-)
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index d9f824e..a4634f9 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -822,7 +822,7 @@ iscenterfilled(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 9b9e5ac..9b8cf8e 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -676,7 +676,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/doc/gnuastro.texi b/doc/gnuastro.texi
index c583877..3ce6c99 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -19328,7 +19328,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
@@ -19394,10 +19394,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.
address@hidden void gal_fits_img_info (fitsfile @code{*fptr}, int @code{*type},
size_t @code{*ndim}, size_t @code{**dsize})
address@hidden 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
address@hidden, @code{ndim}, and @code{dsize} pointers respectively. If
address@hidden 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/lib/fits.c b/lib/fits.c
index 0f42120..d5f555d 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -459,6 +459,53 @@ 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;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized", __func__,
+ *type);
+ }
+
+ /* If the type must be a float, then do the conversion. */
+ if(tofloat) *type=GAL_TYPE_FLOAT32;
+}
+
+
+
+
+
@@ -1226,14 +1273,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 +1297,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;} break;
+ case 3: if(name) {str = key->array; *name = *str;} 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 +1359,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 +1372,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 +1387,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 +1513,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 +1613,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 +1633,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);