[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 808c95d: Fits: default size with --wcsdistorti
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 808c95d: Fits: default size with --wcsdistortion=SIP and no input data |
Date: |
Tue, 23 Jun 2020 20:53:19 -0400 (EDT) |
branch: master
commit 808c95dc56baf023928eeab3edf8bc6e3f572de0
Author: Sachin Kumar Singh <sachinkumarsingh092@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Fits: default size with --wcsdistortion=SIP and no input data
Until now, when the Fits program was called with '--wcsdistortion=SIP'
(from a TPV-based WCS) and the input HDU didn't have any data, it would
crash, complaining that a default size is necessary.
With this commit, a set of steps have been added in the function called by
'--wcsdistortion' to allow a default size when the input HDU has no
size.
In the process, a bug was fixed in reading the keywords by
`wcsdistortion_get_tpvparams`. The problem was fixed by reading directly
from 'disseq' parameter in wcs, which also made its reading style similar
to that of `wcsdistortion_get_sipparams`.
---
bin/fits/keywords.c | 38 +++++++++++++++++++++-----
doc/gnuastro.texi | 4 ++-
lib/gnuastro/wcs.h | 3 ++-
lib/wcs.c | 6 ++++-
lib/wcsdistortion.c | 78 ++++++++++++++++++++++++++++++-----------------------
5 files changed, 87 insertions(+), 42 deletions(-)
diff --git a/bin/fits/keywords.c b/bin/fits/keywords.c
index a4e47f3..c941539 100644
--- a/bin/fits/keywords.c
+++ b/bin/fits/keywords.c
@@ -442,9 +442,10 @@ keywords_distortion_wcs(struct fitsparams *p)
{
int nwcs;
size_t ndim, *insize;
- gal_data_t *data=NULL;
char *suffix, *output;
+ gal_data_t *data=NULL;
struct wcsprm *inwcs, *outwcs;
+ size_t *dsize, defaultsize[2]={2000,2000};
/* If the extension has any data, read it, otherwise just make an empty
array. */
@@ -461,11 +462,36 @@ keywords_distortion_wcs(struct fitsparams *p)
p->cp.quietmmap);
}
- /* Read the input WCS structure and convert it to the desired output
- distortion. */
+ /* Read the input's WCS and make sure one exists. */
inwcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
- outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid,
- data?data->dsize:NULL);
+ if(inwcs==NULL)
+ error(EXIT_FAILURE, 0, "%s (hdu %s): doesn't have any WCS structure "
+ "for converting its distortion",
+ p->filename, p->cp.hdu);
+
+ /* In case there is no dataset and the conversion is between TPV to SIP,
+ we need to set a default size and use that for the conversion, but we
+ also need to warn the user. */
+ if(data==NULL)
+ {
+ if( !p->cp.quiet
+ && gal_wcs_distortion_identify(inwcs)==GAL_WCS_DISTORTION_TPV
+ && p->distortionid==GAL_WCS_DISTORTION_SIP )
+ error(0, 0, "no data associated with WCS for distortion "
+ "converstion.\n\n"
+ "The requested conversion can't be done analytically, so a "
+ "solution has to be found by fitting the parameters over a "
+ "grid of pixels. We will use a default grid of %zux%zu pixels "
+ "and will proceed with the conversion. But it would be more "
+ "accurate if it is the size of the image that this WCS is "
+ "associated with",
+ defaultsize[1], defaultsize[0]);
+ dsize=defaultsize;
+ }
+ else dsize=data->dsize;
+
+ /* Do the conversion. */
+ outwcs=gal_wcs_distortion_convert(inwcs, p->distortionid, dsize);
/* Set the output filename. */
if(p->cp.output)
@@ -491,7 +517,7 @@ keywords_distortion_wcs(struct fitsparams *p)
gal_data_free(data);
}
else
- gal_wcs_write(outwcs, output, NULL, PROGRAM_NAME);
+ gal_wcs_write(outwcs, output, p->wcsdistortion, NULL, PROGRAM_NAME);
/* Clean up. */
wcsfree(inwcs);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 4b22347..0c17cf6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -22192,9 +22192,11 @@ number of coordinate representations found into the
space that @code{nwcs}
points to. Please see @code{gal_wcs_read_fitsptr} for more.
@end deftypefun
-@deftypefun gal_wcs_write (struct wcsprm @code{*wcs}, char @code{*filename},
gal_fits_list_key_t @code{*headers}, char @code{*program_string})
+@deftypefun gal_wcs_write (struct wcsprm @code{*wcs}, char @code{*filename},
char @code{*extname}, gal_fits_list_key_t @code{*headers}, char
@code{*program_string})
Write the given WCS structure into the second extension of an empty FITS
header.
The first/primary extension will be empty like the default format of all
Gnuastro outputs.
+When @code{extname!=NULL} it will be used as the FITS extension name.
+Any set of extra headers can also be written through the @code{headers} list
and if @code{program_string!=NULL} it will be used in a commented keyword title
just above the written version information.
@end deftypefun
@deftypefun {struct wcsprm *} gal_wcs_copy (struct wcsprm @code{*wcs})
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 5821e3c..b76eabd 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -93,7 +93,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
*************************************************************/
void
gal_wcs_write(struct wcsprm *wcs, char *filename,
- gal_fits_list_key_t *headers, char *program_string);
+ char *extname, gal_fits_list_key_t *headers,
+ char *program_string);
diff --git a/lib/wcs.c b/lib/wcs.c
index a19d00d..1170b01 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -276,7 +276,8 @@ gal_wcs_read(char *filename, char *hdu, size_t hstartwcs,
*************************************************************/
void
gal_wcs_write(struct wcsprm *wcs, char *filename,
- gal_fits_list_key_t *headers, char *program_string)
+ char *extname, gal_fits_list_key_t *headers,
+ char *program_string)
{
char *wcsstr;
size_t ndim=0;
@@ -306,6 +307,9 @@ gal_wcs_write(struct wcsprm *wcs, char *filename,
fits_delete_key(fptr, "COMMENT", &status);
status=0;
+ if(extname)
+ fits_write_key(fptr, TSTRING, "EXTNAME", extname, "", &status);
+
/* Decompose the 'PCi_j' matrix and 'CDELTi' vector. */
gal_wcs_decompose_pc_cdelt(wcs);
diff --git a/lib/wcsdistortion.c b/lib/wcsdistortion.c
index 56a3b5a..31212f7 100644
--- a/lib/wcsdistortion.c
+++ b/lib/wcsdistortion.c
@@ -74,59 +74,72 @@ static void
wcsdistortion_get_tpvparams(struct wcsprm *wcs, double cd[2][2],
double *pv1, double *pv2)
{
- size_t pv_m=0;
- size_t i, j, k,index=0;
+ const char *cp;
+ double *temp_cd=NULL;
+ struct dpkey *keyp=NULL;
+ size_t i, j, k=0, index=0;
+ struct disprm *disseq=NULL;
/* Make sure a WCS structure is actually given. */
if(wcs==NULL)
error(EXIT_FAILURE, 0, "%s: input WCS structure is NULL", __func__);
- for(i=0,k=0; i<2; ++i)
+ /* For easy reading. */
+ disseq=wcs->lin.disseq;
+ keyp=disseq->dp;
+
+
+ /* Fill the 2-times allocated CD array (cd[][]). Note that the required
+ CD matix is extracted using the `gal_wcs_wrap_matrix` as a single
+ allocated array (temp_cd[]), that is then used to fill cd[][]. */
+ temp_cd=gal_wcs_warp_matrix(wcs);
+ for(i=0; i<2; ++i)
for(j=0; j<2; ++j)
{
/* If a value is present store it, else store 0.*/
- if((wcs->cd)[i] != 0) cd[i][j]=(wcs->cd)[k++];
- else cd[i][j]=0;
+ if(temp_cd[k] != 0) cd[i][j]=temp_cd[k++];
+ else cd[i][j]=0;
/* For a check:
- printf("CD%ld_%ld\t%.10lf\n", i+1, j+1, cd[i][j]);
+ printf("CD%ld_%ld\t%.10E\n", i+1, j+1, cd[i][j]);
*/
}
- for(j=0; j < wcs->npvmax; ++j)
+ for(i=0; i<disseq->ndp; ++i, ++keyp)
{
- if(wcs->pv[j].i == 1)
+ /* For axis1. */
+ if (keyp->j == 1)
{
- /*pv_m is used to check the index in the header.*/
- pv_m=wcs->pv[j].m;
-
- /* `index` is the index of the pv* array.*/
- index = pv_m;
+ /* Ignore any missing keyvalues. */
+ if ( keyp->field == NULL ) continue;
+ cp = strchr(keyp->field, '.') + 1;
+ if (strncmp(cp, "TPV.", 4) != 0) continue;
+ cp += 4;
- if( wcs->pv[pv_m].value != 0 && index == pv_m )
- pv1[index]=wcs->pv[j].value;
- else
- pv1[index]=0;
+ sscanf(cp, "%ld", &index);
+ pv1[index]=disseq->dp[i].value.f;
/* For a check
- printf("PV1_%d\t%.10f\n", pv_m, pv1[k]);
+ printf("PV1_%ld\t%.10f\n", index, pv1[index]);
*/
}
- else if(wcs->pv[j].i == 2)
+
+ /* For axis2. */
+ else if (keyp->j == 2)
{
- /*pv_m is used to check the index in the header.*/
- pv_m=wcs->pv[j].m;
+ /* Ignore any missing keyvalues. */
+ if ( keyp->field == NULL ) continue;
+
+ cp = strchr(keyp->field, '.') + 1;
+ if (strncmp(cp, "TPV.", 4) != 0) continue;
+ cp += 4;
- /* `index` is the index of the pv* array.*/
- index = pv_m;
+ sscanf(cp, "%ld", &index);
- if( wcs->pv[pv_m].value != 0 && index == pv_m )
- pv2[index]=wcs->pv[j].value;
- else
- pv2[index]=0;
+ pv2[index]=disseq->dp[i].value.f;
/* For a check
- printf("PV2_%d\t%.10f\n", pv_m, pv2[k]);
+ printf("PV2_%ld\t%.10f\n", index, pv2[index]);
*/
}
else
@@ -513,7 +526,7 @@ wcsdistortion_intermidate_tpveq(double cd[2][2], double
*pv1, double *pv2,
/* For a check:
{
- size i, j;
+ size_t i, j;
for(i=0; i<=4; ++i)
for(j=0;j<=4;++j)
{
@@ -1063,14 +1076,13 @@ wcsdistortion_calcsip(size_t axis, size_t m, size_t n,
double tpvu[8][8],
double sip_coeff;
if( axis == 1) sip_coeff=tpvu[m][n];
else if(axis == 2) sip_coeff=tpvv[m][n];
- else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
+ else error(EXIT_FAILURE, 0, "%s: axis %zu does not exists! ",
+ __func__, axis);
if( (axis == 1) && (m == 1) && (n == 0) ) sip_coeff = sip_coeff - 1.0;
else if( (axis == 2) && (m == 0) && (n == 1) ) sip_coeff = sip_coeff - 1.0;
- else error(EXIT_FAILURE, 0, "%s: axis does not exists! ", __func__);
return sip_coeff;
-
}
@@ -1432,7 +1444,7 @@ wcsdistortion_add_sipkeywords(struct wcsprm *wcs, size_t
*fitsize,
int size = wcs->naxis;
size_t a_order=0, b_order=0;
size_t ap_order=0, bp_order=0;
- size_t m, n, num=0, numkey=100;
+ size_t m, n, num=0, numkey=200;
double ap_coeff[5][5], bp_coeff[5][5];
char *fullheader, fmt[50], sipkey[8], keyaxis[9], pcaxis[10];
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 808c95d: Fits: default size with --wcsdistortion=SIP and no input data,
Mohammad Akhlaghi <=