[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 171e3c4: WCS conversion functions using gal_da
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 171e3c4: WCS conversion functions using gal_data_t |
|
Date: |
Wed, 25 Oct 2017 20:26:59 -0400 (EDT) |
branch: master
commit 171e3c4f8d652aa685504d217996bbe9ef77262c
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
WCS conversion functions using gal_data_t
Until now, the two WCS conversion functions (`gal_wcs_world_to_img' and
`gal_wcs_img_to_world') would explicitly assume RA and Dec and work based
on input arrays (so for example it was also necessary to give the number of
elements and etc). In other words, they weren't using the advantages of
`gal_data_t'. With this commit, they have been modified to accept
`gal_data_t' as input for the input coordinates.
The calls to these functions in Crop, MakeCatalog and MakeProfiles have
also been modified to work with the new API.
---
NEWS | 7 ++
bin/crop/wcsmode.c | 44 ++++----
bin/mkcatalog/columns.c | 259 ++++++++++++++++++++++----------------------
bin/mkcatalog/main.h | 12 +--
bin/mkcatalog/mkcatalog.c | 73 ++++++-------
bin/mkcatalog/ui.c | 29 ++---
bin/mkcatalog/upperlimit.c | 10 +-
bin/mkprof/ui.c | 47 +++++---
doc/gnuastro.texi | 83 +++++++-------
lib/fits.c | 2 +-
lib/gnuastro/wcs.h | 10 +-
lib/wcs.c | 262 ++++++++++++++++++++++++++++++++-------------
12 files changed, 484 insertions(+), 354 deletions(-)
diff --git a/NEWS b/NEWS
index 48266a9..d297e62 100644
--- a/NEWS
+++ b/NEWS
@@ -91,6 +91,13 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
`gal_wcs_pixel_area_arcsec2' will return NaN (instead of aborting) when
input is unreasonable (not two dimensions or not in units of degrees).
+ `gal_wcs_world_to_img' and `gal_wcs_img_to_world': Until now, these two
+ WCS conversion functions would explicitly assume RA and Dec and work
+ based on input arrays (so for example it was also necessary to give the
+ number of elements and etc). They now accept `gal_data_t' as input for
+ the input coordinates, thus their API has been greatly simplified and
+ their functionality increased.
+
** Bug fixes
ConvertType crash when changing values (bug #52010).
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index e443db4..4065ab7 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -319,41 +319,43 @@ wcsmode_crop_corners(struct onecropparams *crp)
void
fillcrpipolygon(struct onecropparams *crp)
{
- size_t i;
- double *x, *y, *ra, *dec;
struct cropparams *p=crp->p;
+ gal_data_t *tmp, *coords=NULL;
+ size_t i, d, ndim=p->imgs->ndim;
+
+ /* Allocate the necessary arrays for each column. */
+ for(d=0;d<ndim;++d)
+ gal_list_data_add_alloc(&coords, NULL, GAL_TYPE_FLOAT64, 1, &p->nvertices,
+ NULL, 0, -1, NULL, NULL, NULL);
- /* Allocate the necessary arrays. */
- x=gal_data_calloc_array(GAL_TYPE_FLOAT64, p->nvertices, __func__, "x");
- y=gal_data_calloc_array(GAL_TYPE_FLOAT64, p->nvertices, __func__, "y");
- ra=gal_data_calloc_array(GAL_TYPE_FLOAT64, p->nvertices, __func__, "ra");
- dec=gal_data_calloc_array(GAL_TYPE_FLOAT64, p->nvertices, __func__, "dec");
- crp->ipolygon=gal_data_malloc_array(GAL_TYPE_FLOAT64, 2*p->nvertices,
- __func__, "crp->ipolygon");
- /* Fill in the RA and Dec columns. */
+ /* Fill in the world coordinate columns. */
for(i=0;i<p->nvertices;++i)
{
- ra[i]=p->wpolygon[i*2];
- dec[i]=p->wpolygon[i*2+1];
+ d=0;
+ for(tmp=coords;tmp!=NULL;tmp=tmp->next)
+ ((double *)(tmp->array))[i] = p->wpolygon[ i * ndim + d++ ];
}
+
/* Convert them to image coordinates. */
- gal_wcs_world_to_img(p->imgs[crp->in_ind].wcs, ra, dec, &x, &y,
- p->nvertices);
+ gal_wcs_world_to_img(coords, p->imgs[crp->in_ind].wcs, 1);
- /* Put them in the image polygon vertice array. */
+
+ /* Allocate the image polygon array, and put the image polygon vertice
+ values into it. */
+ crp->ipolygon=gal_data_malloc_array(GAL_TYPE_FLOAT64, ndim*p->nvertices,
+ __func__, "crp->ipolygon");
for(i=0;i<p->nvertices;++i)
{
- crp->ipolygon[i*2] = x[i];
- crp->ipolygon[i*2+1] = y[i];
+ d=0;
+ for(tmp=coords;tmp!=NULL;tmp=tmp->next)
+ crp->ipolygon[ i * ndim + d++ ] = ((double *)(tmp->array))[i];
}
+
/* Clean up. */
- free(x);
- free(y);
- free(ra);
- free(dec);
+ gal_list_data_free(coords);
}
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index a3e08f4..576fcd0 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -55,36 +55,21 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
static void
columns_alloc_radec(struct mkcatalogparams *p)
{
- if(p->rd_vo==NULL)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_vo = malloc(p->input->ndim * sizeof *p->rd_vo);
- if(p->rd_vo==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for p->rd_vo",
- __func__, p->input->ndim * sizeof *p->rd_vo );
-
- /* Space for each dimension. */
- p->rd_vo[0] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_vo[0]");
- p->rd_vo[1] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_vo[1]");
- if(p->clumps)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_vc = malloc(p->input->ndim * sizeof *p->rd_vc);
- if(p->rd_vc==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for p->rd_vo",
- __func__, p->input->ndim * sizeof *p->rd_vc );
-
- /* Space for each dimension. */
- p->rd_vc[0]=gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numclumps,
- __func__, "p->rd_vc[0]");
- p->rd_vc[1]=gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numclumps,
- __func__, "p->rd_vc[1]");
- }
- }
+ size_t i;
+
+ /* For objects. */
+ if(p->wcs_vo==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_vo, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numobjects, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
+
+ /* For clumps */
+ if(p->clumps && p->wcs_vc==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_vc, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numclumps, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
}
@@ -95,36 +80,21 @@ columns_alloc_radec(struct mkcatalogparams *p)
static void
columns_alloc_georadec(struct mkcatalogparams *p)
{
- if(p->rd_go==NULL)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_go = malloc(p->input->ndim * sizeof *p->rd_go);
- if(p->rd_go==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `p->rd_go'",
- __func__, p->input->ndim * sizeof *p->rd_go );
-
- /* Space for each dimension. */
- p->rd_go[0] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_go[0]");
- p->rd_go[1] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_go[1]");
- if(p->clumps)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_gc = malloc(p->input->ndim * sizeof *p->rd_gc);
- if(p->rd_gc==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `p->rd_go'",
- __func__, p->input->ndim * sizeof *p->rd_gc );
-
- /* Space for each dimension. */
- p->rd_gc[0]=gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numclumps,
- __func__, "p->rd_gc[0]");
- p->rd_gc[1]=gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numclumps,
- __func__, "p->rd_gc[1]");
- }
- }
+ size_t i;
+
+ /* For objects. */
+ if(p->wcs_go==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_go, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numobjects, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
+
+ /* For clumps */
+ if(p->clumps && p->wcs_gc==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_gc, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numclumps, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
}
@@ -135,21 +105,13 @@ columns_alloc_georadec(struct mkcatalogparams *p)
static void
columns_alloc_clumpsradec(struct mkcatalogparams *p)
{
- if(p->rd_vcc==NULL)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_vcc = malloc(p->input->ndim * sizeof *p->rd_vcc);
- if(p->rd_vcc==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `p->rd_vcc'",
- __func__, p->input->ndim * sizeof *p->rd_vcc );
-
- /* Space for each dimension. */
- p->rd_vcc[0] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_vcc[0]");
- p->rd_vcc[1] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_vcc[1]");
- }
+ size_t i;
+
+ if(p->wcs_vcc==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_vcc, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numobjects, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
}
@@ -160,21 +122,44 @@ columns_alloc_clumpsradec(struct mkcatalogparams *p)
static void
columns_alloc_clumpsgeoradec(struct mkcatalogparams *p)
{
- if(p->rd_gcc==NULL)
- {
- /* Allocate the space for all dimensions. */
- errno=0;
- p->rd_gcc = malloc(p->input->ndim * sizeof *p->rd_gcc);
- if(p->rd_gcc==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `p->rd_gcc'",
- __func__, p->input->ndim * sizeof *p->rd_gcc );
-
- /* Space for each dimension. */
- p->rd_gcc[0] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_gcc[0]");
- p->rd_gcc[1] = gal_data_malloc_array(GAL_TYPE_FLOAT64, p->numobjects,
- __func__, "p->rd_gcc[1]");
- }
+ size_t i;
+
+ if(p->wcs_gcc==NULL)
+ for(i=0;i<p->input->ndim;++i)
+ gal_list_data_add_alloc(&p->wcs_gcc, NULL, GAL_TYPE_FLOAT64, 1,
+ &p->numobjects, NULL, 0, p->cp.minmapsize,
+ NULL, NULL, NULL);
+}
+
+
+
+
+
+/* Set pointers to fascilitate filling in the values. */
+#define SET_WCS_PREPARE(ARR, LIST, ARRNAME) { \
+ d=0; \
+ errno=0; \
+ (ARR)=malloc(p->input->ndim * sizeof (ARR) ); \
+ if( (ARR)==NULL ) \
+ error(EXIT_FAILURE, 0, "%s: %zu bytes for %s", __func__, \
+ p->input->ndim * sizeof (ARR), ARRNAME); \
+ for(tmp=(LIST);tmp!=NULL;tmp=tmp->next) (ARR)[d++]=tmp->array; \
+ }
+
+static void
+columns_set_wcs_pointers(struct mkcatalogparams *p, double ***vo,
+ double ***vc, double ***go, double ***gc,
+ double ***vcc, double ***gcc)
+{
+ size_t d;
+ gal_data_t *tmp;
+
+ if(p->wcs_vo) SET_WCS_PREPARE(*vo, p->wcs_vo, "vo" );
+ if(p->wcs_vc) SET_WCS_PREPARE(*vc, p->wcs_vc, "vc" );
+ if(p->wcs_go) SET_WCS_PREPARE(*go, p->wcs_go, "go" );
+ if(p->wcs_gc) SET_WCS_PREPARE(*gc, p->wcs_gc, "gc" );
+ if(p->wcs_vcc) SET_WCS_PREPARE(*vcc, p->wcs_vcc, "vcc");
+ if(p->wcs_gcc) SET_WCS_PREPARE(*gcc, p->wcs_gcc, "gcc");
}
@@ -399,7 +384,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_X:
name = "X";
- unit = "position";
+ unit = "pixel";
ocomment = "Flux weighted center (FITS axis 1).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT32;
@@ -413,7 +398,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_Y:
name = "Y";
- unit = "position";
+ unit = "pixel";
ocomment = "Flux weighted center (FITS axis 2).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT32;
@@ -427,7 +412,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_GEOX:
name = "GEO_X";
- unit = "position";
+ unit = "pixel";
ocomment = "Geometric center (FITS axis 1).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT32;
@@ -441,7 +426,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_GEOY:
name = "GEO_Y";
- unit = "position";
+ unit = "pixel";
ocomment = "Geometric center (FITS axis 2).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT32;
@@ -455,7 +440,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSX:
name = "CLUMPS_X";
- unit = "position";
+ unit = "pixel";
ocomment = "Flux weighted center of clumps (FITS axis 1).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT32;
@@ -468,7 +453,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSY:
name = "CLUMPS_Y";
- unit = "position";
+ unit = "pixel";
ocomment = "Flux weighted center of clumps (FITS axis 2).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT32;
@@ -481,7 +466,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSGEOX:
name = "CLUMPS_GEO_X";
- unit = "position";
+ unit = "pixel";
ocomment = "Geometric center of clumps (FITS axis 1).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT32;
@@ -494,7 +479,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSGEOY:
name = "CLUMPS_GEO_Y";
- unit = "position";
+ unit = "pixel";
ocomment = "Geometric center of clumps (FITS axis 2).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT32;
@@ -508,37 +493,41 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_W1:
name = p->ctype[0];
unit = p->input->wcs->cunit[0];
- ocomment = "Flux weighted center in first WCS axis.";
+ ocomment = "Flux weighted center (WCS axis 1).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
disp_width = 13;
disp_precision = 7;
- oiflag[ OCOL_C_VX ] = 1;
- oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_VX ] = 1;
+ oiflag[ OCOL_VY ] = 1;
+ oiflag[ CCOL_VX ] = 1;
+ oiflag[ CCOL_VY ] = 1;
columns_alloc_radec(p);
break;
case UI_KEY_W2:
name = p->ctype[1];
unit = p->input->wcs->cunit[1];
- ocomment = "Flux weighted center in second WCS axis.";
+ ocomment = "Flux weighted center (WCS axis 2).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
disp_width = 13;
disp_precision = 7;
- oiflag[ OCOL_C_VX ] = 1;
- oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_VX ] = 1;
+ oiflag[ OCOL_VY ] = 1;
+ oiflag[ CCOL_VX ] = 1;
+ oiflag[ CCOL_VY ] = 1;
columns_alloc_radec(p);
break;
case UI_KEY_GEOW1:
name = gal_checkset_malloc_cat("GEO_", p->ctype[0]);
unit = p->input->wcs->cunit[0];
- ocomment = "Geometric center in first WCS axis.";
+ ocomment = "Geometric center (WCS axis 1).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -555,7 +544,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_GEOW2:
name = gal_checkset_malloc_cat("GEO_", p->ctype[1]);
unit = p->input->wcs->cunit[1];
- ocomment = "Geometric center in second WCS axis.";
+ ocomment = "Geometric center (WCS axis 2).";
ccomment = ocomment;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_FLOAT64;
@@ -587,7 +576,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSW2:
name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[1]);
unit = p->input->wcs->cunit[1];
- ocomment = "Flux.wht center of all clumps in 2nd WCS axis.";
+ ocomment = "Flux.wht center of all clumps (WCS axis 1).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -602,7 +591,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSGEOW1:
name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[0]);
unit = p->input->wcs->cunit[0];
- ocomment = "Geometric center of all clumps in 1st WCS axis.";
+ ocomment = "Geometric center of all clumps (WCS axis 1).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -617,7 +606,7 @@ columns_define_alloc(struct mkcatalogparams *p)
case UI_KEY_CLUMPSGEOW2:
name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[1]);
unit = p->input->wcs->cunit[1];
- ocomment = "Geometric center of all clumps in 2nd WCS axis.";
+ ocomment = "Geometric center of all clumps (WCS axis 2).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -1254,6 +1243,11 @@ columns_fill(struct mkcatalog_passparams *pp)
double *ci, *oi=pp->oi;
size_t sr=pp->clumpstartindex, cind, coind;
size_t oind=pp->object-1; /* IDs start from 1, indexs from 0. */
+ double **vo=NULL, **vc=NULL, **go=NULL, **gc=NULL, **vcc=NULL, **gcc=NULL;
+
+ /* If a WCS column is requested (check will be done inside the function),
+ then set the pointers. */
+ columns_set_wcs_pointers(p, &vo, &vc, &go, &gc, &vcc, &gcc);
/* Go over all the object columns and fill in the information. */
for(column=p->objectcols; column!=NULL; column=column->next)
@@ -1325,30 +1319,30 @@ columns_fill(struct mkcatalog_passparams *pp)
case UI_KEY_W1:
case UI_KEY_W2:
- p->rd_vo[0][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL,
- OCOL_VX, OCOL_GX);
- p->rd_vo[1][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL,
- OCOL_VY, OCOL_GY);
+ vo[0][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL, OCOL_VX,
+ OCOL_GX);
+ vo[1][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL, OCOL_VY,
+ OCOL_GY);
break;
case UI_KEY_GEOW1:
case UI_KEY_GEOW2:
- p->rd_go[0][oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
- p->rd_go[1][oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
+ go[0][oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
+ go[1][oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
break;
case UI_KEY_CLUMPSW1:
case UI_KEY_CLUMPSW2:
- p->rd_vcc[0][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
- OCOL_C_VX, OCOL_C_GX);
- p->rd_vcc[1][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
- OCOL_C_VY, OCOL_C_GY);
+ vcc[0][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL, OCOL_C_VX,
+ OCOL_C_GX);
+ vcc[1][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL, OCOL_C_VY,
+ OCOL_C_GY);
break;
case UI_KEY_CLUMPSGEOW1:
case UI_KEY_CLUMPSGEOW2:
- p->rd_gcc[0][oind] = MKC_RATIO( oi[OCOL_C_GX], oi[OCOL_C_NUMALL] );
- p->rd_gcc[1][oind] = MKC_RATIO( oi[OCOL_C_GY], oi[OCOL_C_NUMALL] );
+ gcc[0][oind] = MKC_RATIO( oi[OCOL_C_GX], oi[OCOL_C_NUMALL] );
+ gcc[1][oind] = MKC_RATIO( oi[OCOL_C_GY], oi[OCOL_C_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
@@ -1490,16 +1484,16 @@ columns_fill(struct mkcatalog_passparams *pp)
case UI_KEY_W1:
case UI_KEY_W2:
- p->rd_vc[0][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
- CCOL_VX, CCOL_GX);
- p->rd_vc[1][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
- CCOL_VY, CCOL_GY);
+ vc[0][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL, CCOL_VX,
+ CCOL_GX);
+ vc[1][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL, CCOL_VY,
+ CCOL_GY);
break;
case UI_KEY_GEOW1:
case UI_KEY_GEOW2:
- p->rd_gc[0][cind] = MKC_RATIO( ci[CCOL_GX], ci[CCOL_NUMALL] );
- p->rd_gc[1][cind] = MKC_RATIO( ci[CCOL_GY], ci[CCOL_NUMALL] );
+ gc[0][cind] = MKC_RATIO( ci[CCOL_GX], ci[CCOL_NUMALL] );
+ gc[1][cind] = MKC_RATIO( ci[CCOL_GY], ci[CCOL_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
@@ -1517,7 +1511,8 @@ columns_fill(struct mkcatalog_passparams *pp)
break;
case UI_KEY_NORIVERBRIGHTNESS:
- ((float *)colarr)[cind] = ci[ CCOL_NUM ]>0.0f ? ci[ CCOL_SUM ]:NAN;
+ ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
+ ? ci[ CCOL_SUM ] : NAN );
break;
case UI_KEY_MAGNITUDE: /* Similar: brightness for clumps */
@@ -1606,4 +1601,12 @@ columns_fill(struct mkcatalog_passparams *pp)
key);
}
}
+
+ /* Clean up. */
+ if(vo) free(vo);
+ if(vc) free(vc);
+ if(go) free(go);
+ if(gc) free(gc);
+ if(vcc) free(vcc);
+ if(gcc) free(gcc);
}
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 39db6f6..178b0f5 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -188,12 +188,12 @@ struct mkcatalogparams
uint64_t seed; /* Random number generator seed. */
const char *rngname; /* Name of random number generator. */
- double **rd_vo; /* Object RA-Dec flux weighted X, Y. */
- double **rd_vc; /* Clump RA-Dec flux weighted X, Y. */
- double **rd_go; /* Object RA-Dec geometric X,Y. */
- double **rd_gc; /* Clump RA-Dec geometric X, Y. */
- double **rd_vcc; /* All clumps RA-Dec flx. wht. X, Y. */
- double **rd_gcc; /* All clumps RA-Dec geometric X, Y. */
+ gal_data_t *wcs_vo; /* Object RA-Dec flux weighted X, Y. */
+ gal_data_t *wcs_vc; /* Clump RA-Dec flux weighted X, Y. */
+ gal_data_t *wcs_go; /* Object RA-Dec geometric X,Y. */
+ gal_data_t *wcs_gc; /* Clump RA-Dec geometric X, Y. */
+ gal_data_t *wcs_vcc; /* All clumps RA-Dec flx. wht. X, Y. */
+ gal_data_t *wcs_gcc; /* All clumps RA-Dec geometric X, Y. */
char **ctype; /* Type of WCS axis. */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index e41278d..fb167ca 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -559,38 +559,35 @@ mkcatalog_single_object(void *in_prm)
static void
mkcatalog_wcs_conversion(struct mkcatalogparams *p)
{
- double *c, *d, *df;
+ gal_data_t *c;
gal_data_t *column;
/* Flux weighted center positions for clumps and objects. */
- if(p->rd_vo)
+ if(p->wcs_vo)
{
- gal_wcs_img_to_world(p->input->wcs, p->rd_vo[0], p->rd_vo[1],
- &p->rd_vo[0], &p->rd_vo[1], p->numobjects);
- if(p->rd_vc)
- gal_wcs_img_to_world(p->input->wcs, p->rd_vc[0], p->rd_vc[1],
- &p->rd_vc[0], &p->rd_vc[1], p->numclumps);
+ gal_wcs_img_to_world(p->wcs_vo, p->input->wcs, 1);
+ if(p->wcs_vc)
+ gal_wcs_img_to_world(p->wcs_vc, p->input->wcs, 1);
}
+
/* Geometric center positions for clumps and objects. */
- if(p->rd_go)
+ if(p->wcs_go)
{
- gal_wcs_img_to_world(p->input->wcs, p->rd_go[0], p->rd_go[1],
- &p->rd_go[0], &p->rd_go[1], p->numobjects);
- if(p->rd_gc)
- gal_wcs_img_to_world(p->input->wcs, p->rd_gc[0], p->rd_gc[1],
- &p->rd_gc[0], &p->rd_gc[1], p->numclumps);
+ gal_wcs_img_to_world(p->wcs_go, p->input->wcs, 1);
+ if(p->wcs_gc)
+ gal_wcs_img_to_world(p->wcs_gc, p->input->wcs, 1);
}
+
/* All clumps flux weighted center. */
- if(p->rd_vcc)
- gal_wcs_img_to_world(p->input->wcs, p->rd_vcc[0], p->rd_vcc[1],
- &p->rd_vcc[0], &p->rd_vcc[1], p->numobjects);
+ if(p->wcs_vcc)
+ gal_wcs_img_to_world(p->wcs_vcc, p->input->wcs, 1);
+
/* All clumps geometric center. */
- if(p->rd_gcc)
- gal_wcs_img_to_world(p->input->wcs, p->rd_gcc[0], p->rd_gcc[1],
- &p->rd_gcc[0], &p->rd_gcc[1], p->numobjects);
+ if(p->wcs_gcc)
+ gal_wcs_img_to_world(p->wcs_gcc, p->input->wcs, 1);
/* Go over all the object columns and fill in the values. */
@@ -604,18 +601,20 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
probably columns that don't need any correction. */
switch(column->status)
{
- case UI_KEY_W1: c=p->rd_vo[0]; break;
- case UI_KEY_W2: c=p->rd_vo[1]; break;
- case UI_KEY_GEOW1: c=p->rd_go[0]; break;
- case UI_KEY_GEOW2: c=p->rd_go[1]; break;
- case UI_KEY_CLUMPSW1: c=p->rd_vcc[0]; break;
- case UI_KEY_CLUMPSW2: c=p->rd_vcc[1]; break;
- case UI_KEY_CLUMPSGEOW1: c=p->rd_gcc[0]; break;
- case UI_KEY_CLUMPSGEOW2: c=p->rd_gcc[1]; break;
+ case UI_KEY_W1: c=p->wcs_vo; break;
+ case UI_KEY_W2: c=p->wcs_vo->next; break;
+ case UI_KEY_GEOW1: c=p->wcs_go; break;
+ case UI_KEY_GEOW2: c=p->wcs_go->next; break;
+ case UI_KEY_CLUMPSW1: c=p->wcs_vcc; break;
+ case UI_KEY_CLUMPSW2: c=p->wcs_vcc->next; break;
+ case UI_KEY_CLUMPSGEOW1: c=p->wcs_gcc; break;
+ case UI_KEY_CLUMPSGEOW2: c=p->wcs_gcc->next; break;
}
- /* Copy the elements. */
- if(c) { df=(d=column->array)+column->size; do *d=*c++; while(++d<df); }
+ /* Copy the elements into the output column. */
+ if(c)
+ memcpy(column->array, c->array,
+ column->size*gal_type_sizeof(c->type));
}
@@ -630,14 +629,16 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
probably columns that don't need any correction. */
switch(column->status)
{
- case UI_KEY_W1: c=p->rd_vc[0]; break;
- case UI_KEY_W2: c=p->rd_vc[1]; break;
- case UI_KEY_GEOW1: c=p->rd_gc[0]; break;
- case UI_KEY_GEOW2: c=p->rd_gc[1]; break;
+ case UI_KEY_W1: c=p->wcs_vc; break;
+ case UI_KEY_W2: c=p->wcs_vc->next; break;
+ case UI_KEY_GEOW1: c=p->wcs_gc; break;
+ case UI_KEY_GEOW2: c=p->wcs_gc->next; break;
}
- /* Copy the elements. */
- if(c) { df=(d=column->array)+column->size; do *d=*c++; while(++d<df); }
+ /* Copy the elements into the output column. */
+ if(c)
+ memcpy(column->array, c->array,
+ column->size*gal_type_sizeof(c->type));
}
}
@@ -882,7 +883,7 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
- /* OBJECT CATALOG
+ /* CLUMPS CATALOG
============== */
if(p->clumps)
{
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 9b86483..f196824 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -1010,28 +1010,13 @@ ui_free_report(struct mkcatalogparams *p, struct
timeval *t1)
{
size_t d;
- /* The temporary arrays for RA and Dec. */
- if(p->rd_vo || p->rd_vc || p->rd_go || p->rd_gc || p->rd_vcc || p->rd_gcc)
- {
- /* First free the separate dimensions in each value. */
- for(d=0;d<p->input->ndim;++d)
- {
- if(p->rd_vo) free(p->rd_vo[d]);
- if(p->rd_vc) free(p->rd_vc[d]);
- if(p->rd_go) free(p->rd_go[d]);
- if(p->rd_gc) free(p->rd_gc[d]);
- if(p->rd_vcc) free(p->rd_vcc[d]);
- if(p->rd_gcc) free(p->rd_gcc[d]);
- }
-
- /* Then free the container arrays. */
- if(p->rd_vo) free(p->rd_vo);
- if(p->rd_vc) free(p->rd_vc);
- if(p->rd_go) free(p->rd_go);
- if(p->rd_gc) free(p->rd_gc);
- if(p->rd_vcc) free(p->rd_vcc);
- if(p->rd_gcc) free(p->rd_gcc);
- }
+ /* The temporary arrays for WCS coordinates. */
+ if(p->wcs_vo ) gal_list_data_free(p->wcs_vo);
+ if(p->wcs_vc ) gal_list_data_free(p->wcs_vc);
+ if(p->wcs_go ) gal_list_data_free(p->wcs_go);
+ if(p->wcs_gc ) gal_list_data_free(p->wcs_gc);
+ if(p->wcs_vcc) gal_list_data_free(p->wcs_vcc);
+ if(p->wcs_gcc) gal_list_data_free(p->wcs_gcc);
/* Free the types of the WCS coordinates (for catalog meta-data). */
if(p->ctype)
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index 44d4184..136d096 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -179,12 +179,13 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp,
gal_data_t *tile,
st_oc = clumplab ? (int32_t *)(p->clumps->array) + se_inc[0] : NULL;
- /* Continue measuring magnitudes randomly until we get the desired
- number. */
+ /* Continue measuring randomly until we get the desired total number. */
while(tcounter<maxcount && counter<p->upnum)
{
/* Get the random coordinates, note that `gsl_rng_uniform_int'
- returns an inclusive value. */
+ returns an inclusive value. It may happen that the labeled region
+ extends the full range of a dimension. In that case, the only
+ possible want the random starting point would be 0. */
for(d=0;d<ndim;++d)
rcoord[d] = ( (dsize[d]-tile->dsize[d])
? gsl_rng_uniform_int(pp->rng, dsize[d]-tile->dsize[d]-1)
@@ -209,9 +210,6 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp,
gal_data_t *tile,
st_sky = (float *)(p->sky->array) + se_inc[0];
if(p->upmask) st_m = (uint8_t *)(p->upmask->array) + se_inc[0];
-
- /* Starting pointers for the original tile.*/
-
/* Parse over this object/clump. */
while( se_inc[0] + increment <= se_inc[1] )
{
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index fb423ce..16fe676 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -1273,9 +1273,10 @@ ui_prepare_canvas(struct mkprofparams *p)
static void
ui_finalize_coordinates(struct mkprofparams *p)
{
- size_t i, ndim=p->ndim;
- double *x=NULL, *y=NULL;
+ void *arr=NULL;
+ size_t i=0, ndim=p->ndim;
uint8_t os=p->oversample;
+ gal_data_t *tmp, *coords=NULL;
double *cdelt=p->wcs->cdelt, *crpix=p->wcs->crpix;
/* When the user specified RA and Dec columns, the respective values
@@ -1283,24 +1284,46 @@ ui_finalize_coordinates(struct mkprofparams *p)
need to change them into actual image coordinates. */
if(p->mode==MKPROF_MODE_WCS)
{
- /* Note that we read the RA and Dec columns into the `p->x' and `p->y'
- arrays temporarily before. Here, we will convert them, free the old
- ones and replace them with the proper X and Y values. */
- gal_wcs_world_to_img(p->wcs, p->x, p->y, &x, &y, p->num);
+ /* Make list of coordinates for input of `gal_wcs_world_to_img'. */
+ for(i=0;i<ndim;++i)
+ {
+ /* Set the array pointer. Note that we read the WCS columns into
+ the `p->x', `p->y' and `p->z' arrays temporarily before. Here, we
+ will convert them to image coordinates in place. */
+ switch(i)
+ {
+ /* Note that the linked list gets filled in a first-in-last-out
+ order, so the last column added should be the first WCS
+ dimension. */
+ case 0: arr = p->y; break;
+ case 1: arr = p->x; break;
+ default:
+ error(EXIT_FAILURE, 0, "conversion from WCS to image "
+ "coordinates is not currently supported for "
+ "%zu-dimensional datasets", ndim);
+ }
+
+ /* Allocate the list of coordinates. */
+ gal_list_data_add_alloc(&coords, arr, GAL_TYPE_FLOAT64, 1, &p->num,
+ NULL, 0, -1, NULL, NULL, NULL);
+ }
+
+ /* Convert the world coordinates to image coordinates (inplace). */
+ gal_wcs_world_to_img(coords, p->wcs, 1);
+
/* If any conversions created a WCSLIB error, both the outputs will be
set to NaN. */
for(i=0;i<p->num;++i)
- if( isnan(x[i]) )
+ if( isnan(p->x[i]) )
error(EXIT_FAILURE, 0, "catalog row %zu: WCSLIB could not convert "
"(%f, %f) coordinates into image coordinates", i, p->x[i],
p->y[i]);
- /* Free the RA and Dec arrays and put in the new image values. */
- free(p->x);
- free(p->y);
- p->x=x;
- p->y=y;
+ /* We want the actual arrays of each `coords' column. So, first we'll
+ set all the array elements to NULL, then free it. */
+ for(tmp=coords;tmp!=NULL;tmp=tmp->next) tmp->array=NULL;
+ gal_list_data_free(coords);
}
/* Correct the WCS scale. Note that when the WCS is read from a
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index e3e8395..ce5e533 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -13537,12 +13537,11 @@ image separately. Not one value for the whole image.
The shape or morphology of a target is one of the most commonly desired
parameters of a target. Here, we will review the derivation of the most
-basic/simple morphological parameters are estimated: the elliptical
-parameters for a set of labeled pixels. The elliptical parameters are: the
-(semi-)major axis, the (semi-)minor axis and the position angle along with
-the central position of the profile. The derivations below follow the
-SExtractor manual derivations with some added explanations for easier
-reading.
+basic/simple morphological parameters: the elliptical parameters for a set
+of labeled pixels. The elliptical parameters are: the (semi-)major axis,
+the (semi-)minor axis and the position angle along with the central
+position of the profile. The derivations below follow the SExtractor manual
+derivations with some added explanations for easier reading.
@cindex Moments
Let's begin with one dimension for simplicity: Assume we have a set of
@@ -13573,7 +13572,7 @@ the distribution:
@noindent
The last step was done from the definition of @mymath{\overline{x}}. Hence,
the square root of @mymath{\overline{x^2}} is the spatial standard
-deviation (along the one-dimensional) of this particular brightness
+deviation (along the one-dimension) of this particular brightness
distribution (@mymath{B_i}). Crudely (or qualitatively), you can think of
its square root as the distance (from @mymath{\overline{x}}) which contains
a specific amount of the flux (depending on the @mymath{B_i}
@@ -13586,21 +13585,24 @@ other words, it quantifies how ``sharp'' the object's
image is.
@cindex Floating point error
Before continuing to two dimensions and the derivation of the elliptical
parameters, let's pause for an important implementation technicality. You
-can ignore this paragraph if you don't want to implement these
-concepts. The basic definition (first fraction for @mymath{\overline{x^2}})
-can be used without any major problem. However, using this fraction
-requires two runs over the data: one run to find @mymath{\overline{x}} and
-one run to find @mymath{\overline{x^2}}, this can be slow. However, using
-the last fraction above, we can estimate both the first and second moments
-in one run (since the @mymath{-\overline{x}^2} term can easily be added
-later). The logarithmic nature of floating point number digitization
-creates a complication in this approach: suppose the object is located
-between pixels 10000 and 10020. Hence the target's pixels are only
-distributed over 20 pixels (with a standard deviation @mymath{<20}), while
-the mean has a value of @mymath{\sim10000}. The @mymath{\sum_iB_i^2x_i^2}
-will go to very very large values while the individual pixel differences
-will be much smaller, this will lower the accuracy of our calculation due
-to the limited accuracy of floating point operations. The variance only
+can ignore this paragraph and the next two if you don't want to implement
+these concepts. The basic definition (first definition of
address@hidden above) can be used without any major
+problem. However, using this fraction requires two runs over the data: one
+run to find @mymath{\overline{x}} and another run to find
address@hidden from @mymath{\overline{x}}, this can be slow. The
+advantage of the last fraction above, is that we can estimate both the
+first and second moments in one run (since the @mymath{-\overline{x}^2}
+term can easily be added later).
+
+The logarithmic nature of floating point number digitization creates a
+complication however: suppose the object is located between pixels 10000
+and 10020. Hence the target's pixels are only distributed over 20 pixels
+(with a standard deviation @mymath{<20}), while the mean has a value of
address@hidden The @mymath{\sum_iB_i^2x_i^2} will go to very very
+large values while the individual pixel differences will be orders of
+magnitude smaller. This will lower the accuracy of our calculation due to
+the limited accuracy of floating point operations. The variance only
depends on the distance of each point from the mean, so we can shift all
position by a constant/arbitrary @mymath{K} which is much closer to the
mean: @mymath{\overline{x-K}=\overline{x}-K}. Hence we can calculate the
@@ -20186,25 +20188,28 @@ structure is not two dimensional and the units
(@code{CUNIT} keywords) are
not @code{deg} (for degrees), then this function will return a NaN.
@end deftypefun
address@hidden void gal_wcs_world_to_img (struct wcsprm @code{*wcs}, double
@code{*ra}, double @code{*dec}, double @code{**x}, double @code{**y}, size_t
@code{size})
-Convert the arrays of input world coordinates (@code{ra} and @code{dec})
-into arrays of image coordinates (@code{x} and @code{y}). Each is assumed
-to be a separate one-dimensional array of @code{size} elements. If
address@hidden or @code{*y==NULL}, then space will be allocated for them
-by this function, otherwise, it is assumed that they already contain the
-space necessary to write the values.
-
-If you don't need the input values after this conversion any more, you can
-pass the pointers of the @code{ra} and @code{dec} arrays and the outputs
-will be written into them. This can help to avoid extra allocations and
-freeing.
address@hidden {gal_data_t *} gal_wcs_world_to_img (gal_data_t @code{*coords},
struct wcsprm @code{*wcs}, int @code{inplace})
+Convert the linked list of world coordinates in @code{coords} to a linked
+list of image coordinates given the input WCS structure. @code{coords} must
+be a linked list of data structures of float64 (`double') type,
address@hidden lists} and @ref{List of gal_data_t}. The top (first
+popped/read) node of the linked list must be the first WCS coordinate (RA
+in an image usually) and etc. Similarly, the top node of the output will be
+the first image coordinate (in the FITS standard).
+
+If @code{inplace} is zero, then the output will be a newly allocated list
+and the input list will be untouched. However, if @code{inplace} is
+non-zero, the output values will be written into the input's already
+allocated array and the returned pointer will be the same pointer to
address@hidden (in other words, you can ignore the returned value). Note
+that in the latter case, only the values will be changed, things like units
+or name (if present) will be untouched.
@end deftypefun
address@hidden void gal_wcs_img_to_world (struct wcsprm @code{*wcs}, double
@code{*x}, double @code{*y}, double @code{**ra}, double @code{**dec}, size_t
@code{size})
-Convert the arrays of input image coordinates (@code{x} and @code{y}) into
-arrays of world coordinates (@code{ra} and @code{dec}). Each is assumed to
-be a separate one-dimensional array of @code{size} elements. See
address@hidden for more.
address@hidden {gal_data_t *} gal_wcs_img_to_world (gal_data_t @code{*coords},
struct wcsprm @code{*wcs}, int @code{inplace})
+Convert the linked list of image coordinates in @code{coords} to a linked
+list of world coordinates given the input WCS structure. See the
+description of @code{gal_wcs_world_to_img} for more details.
@end deftypefun
diff --git a/lib/fits.c b/lib/fits.c
index 052a5df..33b4f9b 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1670,9 +1670,9 @@ gal_fits_img_write_to_ptr(gal_data_t *input, char
*filename)
fitsfile *fptr;
uint64_t *u64, *u64f;
long fpixel=1, *naxes;
+ char *wcsstr, *u64key;
size_t i, ndim=input->ndim;
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
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index 78c6dec..553b8ae 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -98,13 +98,11 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
/**************************************************************/
/********** Conversion ************/
/**************************************************************/
-void
-gal_wcs_world_to_img(struct wcsprm *wcs, double *ra, double *dec,
- double **x, double **y, size_t size);
+gal_data_t *
+gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm *wcs, int inplace);
-void
-gal_wcs_img_to_world(struct wcsprm *wcs, double *x, double *y,
- double **ra, double **dec, size_t size);
+gal_data_t *
+gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm *wcs, int inplace);
diff --git a/lib/wcs.c b/lib/wcs.c
index 0649e57..38f5734 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -620,59 +620,165 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
/**************************************************************/
/********** Array conversion ************/
/**************************************************************/
-/* Convert an array of world coordinates to image coordinates. Note that in
- Gnuastro, each column is treated independently, so the inputs are
- separate. If `*x==NULL', or `*y==NULL', then space will be allocated for
- them, otherwise, it is assumed that space has already been
- allocated. Note that they must each be a 1 dimensional array.
-
- You can do the conversion in place: just pass the same array as you give
- to RA and Dec to X and Y. */
-void
-gal_wcs_world_to_img(struct wcsprm *wcs, double *ra, double *dec,
- double **x, double **y, size_t size)
+/* Some sanity checks for the WCS conversion functions. */
+static void
+wcs_convert_sanity_check_alloc(gal_data_t *coords, struct wcsprm *wcs,
+ const char *func, int **stat, double **phi,
+ double **theta, double **world,
+ double **pixcrd, double **imgcrd)
{
- size_t i;
- int status, *stat, ncoord=size, nelem=2;
- double *phi, *theta, *world, *pixcrd, *imgcrd;
+ gal_data_t *tmp;
+ size_t ndim=0, firstsize=0, size=coords->size;
+
+ for(tmp=coords; tmp!=NULL; tmp=tmp->next)
+ {
+ /* Count how many coordinates are given. */
+ ++ndim;
+
+ /* Check the type of the input. */
+ if(tmp->type!=GAL_TYPE_FLOAT64)
+ error(EXIT_FAILURE, 0, "%s: input coordinates must have `float64' "
+ "type", func);
+
+ /* Make sure it has a single dimension. */
+ if(tmp->ndim!=1)
+ error(EXIT_FAILURE, 0, "%s: input coordinates for each dimension "
+ "must each be one dimensional. Coordinate dataset %zu of the "
+ "inputs has %zu dimensions", func, ndim, tmp->ndim);
+
+ /* See if all inputs have the same size. */
+ if(ndim==1) firstsize=tmp->size;
+ else
+ if(firstsize!=tmp->size)
+ error(EXIT_FAILURE, 0, "%s: all input coordinates must have the "
+ "same number of elements. Coordinate dataset %zu has %zu "
+ "elements while the first coordinate has %zu", func, ndim,
+ tmp->size, firstsize);
+ }
+
+ /* See if the number of coordinates given corresponds to the dimensions
+ of the WCS structure. */
+ if(ndim!=wcs->naxis)
+ error(EXIT_FAILURE, 0, "%s: the number of input coordinates (%zu) does "
+ "not match the dimensions of the input WCS structure (%d)", func,
+ ndim, wcs->naxis);
/* Allocate all the necessary arrays. */
- phi = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "phi");
- stat = gal_data_calloc_array( GAL_TYPE_INT32, size, __func__, "stat");
- theta = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "theta");
- world = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "world");
- imgcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "imgcrd");
- pixcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "pixcrd");
-
- /* Write the values into the allocated contiguous array. */
- for(i=0;i<size;++i) { world[i*2]=ra[i]; world[i*2+1]=dec[i]; }
-
- /* Use WCSLIB's `wcss2p'. */
+ *phi = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "phi");
+ *stat = gal_data_calloc_array( GAL_TYPE_INT32, size, __func__, "stat");
+ *theta = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "theta");
+ *world = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
+ "world");
+ *imgcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
+ "imgcrd");
+ *pixcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, ndim*size, __func__,
+ "pixcrd");
+}
+
+
+
+
+
+/* In Gnuastro, each column (coordinate for WCS conversion) is treated as a
+ separate array in a `gal_data_t' that are linked through a linked
+ list. But in WCSLIB, the input is a single array (with multiple
+ columns). This function will convert between the two. */
+static void
+wcs_convert_list_to_array(gal_data_t *list, double *array, int *stat,
+ size_t ndim, int listtoarray)
+{
+ size_t i, d=0;
+ gal_data_t *tmp;
+
+ for(tmp=list; tmp!=NULL; tmp=tmp->next)
+ {
+ /* Put all this coordinate's values into the single array that is
+ input into or output from WCSLIB. */
+ for(i=0;i<list->size;++i)
+ {
+ if(listtoarray)
+ array[i*ndim+d] = ((double *)(tmp->array))[i];
+ else
+ ((double *)(tmp->array))[i] = stat[i] ? NAN : array[i*ndim+d];
+ }
+
+ /* Increment the dimension. */
+ ++d;
+ }
+}
+
+
+
+
+
+/* Prepare the output of the WCS conversion functions. */
+static gal_data_t *
+wcs_convert_prepare_out(gal_data_t *coords, struct wcsprm *wcs, int inplace)
+{
+ size_t i;
+ gal_data_t *out=NULL;
+ if(inplace)
+ out=coords;
+ else
+ for(i=0;i<wcs->naxis;++i)
+ gal_list_data_add_alloc(&out, NULL, GAL_TYPE_FLOAT64, 1,
+ &coords->size, NULL, 0, coords->minmapsize,
+ wcs->ctype[i], wcs->cunit[i], NULL);
+ return out;
+}
+
+
+
+
+
+/* Convert world coordinates to image coordinates given the input WCS
+ structure. The input must be a linked list of data structures of float64
+ (`double') type. The top element of the linked list must be the first
+ coordinate and etc. If `inplace' is non-zero, then the output will be
+ written into the input's allocated space. */
+gal_data_t *
+gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm *wcs, int inplace)
+{
+ gal_data_t *out;
+ int status, *stat=NULL, ncoord=coords->size, nelem=wcs->naxis;
+ double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
+
+ /* Some sanity checks. */
+ wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi, &theta,
+ &world, &pixcrd, &imgcrd);
+
+
+ /* Write the values from the input list of separate columns into a single
+ array (WCSLIB input). */
+ wcs_convert_list_to_array(coords, world, stat, wcs->naxis, 1);
+
+
+ /* Use WCSLIB's wcsp2s for the conversion. */
status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd, pixcrd, stat);
if(status)
error(EXIT_FAILURE, 0, "%s: wcss2p ERROR %d: %s", __func__, status,
wcs_errmsg[status]);
- /* For a sanity check:
- printf("\n\ngal_wcs_world_to_img sanity check:\n");
- for(i=0;i<size;++i)
- printf("world (%f, %f) --> pix (%f, %f), [stat: %d]\n",
- world[i*2], world[i*2+1], pixcrd[i*2], pixcrd[i*2+1], stat[i]);
+
+ /* For a sanity check.
+ {
+ size_t i;
+ printf("\n\n%s sanity check:\n", __func__);
+ for(i=0;i<coords->size;++i)
+ printf("(%g, %g) --> (%g, %g), [stat: %d]\n", world[i*2], world[i*2+1],
+ pixcrd[i*2], pixcrd[i*2+1], stat[i]);
+ }
*/
+
/* Allocate the output arrays if they were not already allocated. */
- if(*x==NULL) *x=gal_data_malloc_array(GAL_TYPE_FLOAT64, size, __func__,"x");
- if(*y==NULL) *y=gal_data_malloc_array(GAL_TYPE_FLOAT64, size, __func__,"y");
+ out=wcs_convert_prepare_out(coords, wcs, inplace);
+
+
+ /* Write the output from a single array (WCSLIB output) into the output
+ list of this function. */
+ wcs_convert_list_to_array(out, pixcrd, stat, wcs->naxis, 0);
- /* Put the values into the output arrays. */
- for(i=0;i<size;++i)
- {
- (*x)[i] = stat[i] ? NAN : pixcrd[i*2];
- (*y)[i] = stat[i] ? NAN : pixcrd[i*2+1];
- }
/* Clean up. */
free(phi);
@@ -680,35 +786,32 @@ gal_wcs_world_to_img(struct wcsprm *wcs, double *ra,
double *dec,
free(theta);
free(world);
free(pixcrd);
+
+ /* Return the output list of coordinates. */
+ return out;
}
-/* Similar to `gal_wcs_world_to_img' but converts image coordinates into
- world coordinates. */
-void
-gal_wcs_img_to_world(struct wcsprm *wcs, double *x, double *y,
- double **ra, double **dec, size_t size)
+/* Similar to `gal_wcs_world_to_img'. */
+gal_data_t *
+gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm *wcs, int inplace)
{
- size_t i;
- int status, *stat, ncoord=size, nelem=2;
- double *phi, *theta, *world, *pixcrd, *imgcrd;
+ gal_data_t *out;
+ int status, *stat=NULL, ncoord=coords->size, nelem=wcs->naxis;
+ double *phi=NULL, *theta=NULL, *world=NULL, *pixcrd=NULL, *imgcrd=NULL;
+
+ /* Some sanity checks. */
+ wcs_convert_sanity_check_alloc(coords, wcs, __func__, &stat, &phi, &theta,
+ &world, &pixcrd, &imgcrd);
+
+
+ /* Write the values from the input list of separate columns into a single
+ array (WCSLIB input). */
+ wcs_convert_list_to_array(coords, pixcrd, stat, wcs->naxis, 1);
- /* Allocate all the necessary arrays. */
- phi = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "phi");
- stat = gal_data_calloc_array( GAL_TYPE_INT32, size, __func__, "stat");
- theta = gal_data_malloc_array( GAL_TYPE_FLOAT64, size, __func__, "theta");
- world = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "world");
- imgcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "imgcrd");
- pixcrd = gal_data_malloc_array( GAL_TYPE_FLOAT64, 2*size, __func__,
- "pixcrd");
-
- /* Write the values into the allocated contiguous array. */
- for(i=0;i<size;++i) { pixcrd[i*2]=x[i]; pixcrd[i*2+1]=y[i]; }
/* Use WCSLIB's wcsp2s for the conversion. */
status=wcsp2s(wcs, ncoord, nelem, pixcrd, imgcrd, phi, theta, world, stat);
@@ -716,25 +819,26 @@ gal_wcs_img_to_world(struct wcsprm *wcs, double *x,
double *y,
error(EXIT_FAILURE, 0, "%s: wcsp2s ERROR %d: %s", __func__, status,
wcs_errmsg[status]);
- /* For a sanity check:
- printf("\n\ngal_wcs_img_to_world sanity check:\n");
- for(i=0;i<size;++i)
- printf("img (%f, %f) --> world (%f, %f), [stat: %d]\n",
- pixcrd[i*2], pixcrd[i*2+1], world[i*2], world[i*2+1], stat[i]);
+
+ /* For a sanity check.
+ {
+ size_t i;
+ printf("\n\n%s sanity check:\n", __func__);
+ for(i=0;i<coords->size;++i)
+ printf("(%g, %g) --> (%g, %g), [stat: %d]\n", pixcrd[i*2], pixcrd[i*2+1],
+ world[i*2], world[i*2+1], stat[i]);
+ }
*/
+
/* Allocate the output arrays if they were not already allocated. */
- if(*ra==NULL)
- *ra = gal_data_malloc_array(GAL_TYPE_FLOAT64, size, __func__, "ra");
- if(*dec==NULL)
- *dec = gal_data_malloc_array(GAL_TYPE_FLOAT64, size, __func__, "dec");
+ out=wcs_convert_prepare_out(coords, wcs, inplace);
+
+
+ /* Write the output from a single array (WCSLIB output) into the output
+ list of this function. */
+ wcs_convert_list_to_array(out, world, stat, wcs->naxis, 0);
- /* Put the values into the output arrays. */
- for(i=0;i<size;++i)
- {
- (*ra)[i] = stat[i] ? NAN : world[i*2];
- (*dec)[i] = stat[i] ? NAN : world[i*2+1];
- }
/* Clean up. */
free(phi);
@@ -742,4 +846,8 @@ gal_wcs_img_to_world(struct wcsprm *wcs, double *x, double
*y,
free(theta);
free(world);
free(pixcrd);
+
+
+ /* Return the output list of coordinates. */
+ return out;
}
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 171e3c4: WCS conversion functions using gal_data_t,
Mohammad Akhlaghi <=