[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master fb3660f 037/113: MakeCatalog works in 3D
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master fb3660f 037/113: MakeCatalog works in 3D |
|
Date: |
Fri, 16 Apr 2021 10:33:39 -0400 (EDT) |
branch: master
commit fb3660fccaa196e29e896432a638f60aa0784d5d
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>
MakeCatalog works in 3D
The positional measurements of MakeProfile were updated to allow its
operation of 3D datasets. The WCS conversion libraries were also updated to
allow varying number of dimensions.
---
bin/crop/wcsmode.c | 44 ++--
bin/mkcatalog/args.h | 112 ++++++++++
bin/mkcatalog/columns.c | 518 ++++++++++++++++++++++++++++++++-----------
bin/mkcatalog/main.h | 19 +-
bin/mkcatalog/mkcatalog.c | 115 ++++++----
bin/mkcatalog/ui.c | 37 +---
bin/mkcatalog/ui.h | 10 +-
bin/mkcatalog/upperlimit.c | 15 +-
bin/mkprof/astmkprof-3d.conf | 2 +-
bin/mkprof/ui.c | 52 +++--
doc/gnuastro.texi | 146 +++++++-----
lib/fits.c | 2 +-
lib/gnuastro/wcs.h | 10 +-
lib/wcs.c | 264 +++++++++++++++-------
tests/Makefile.am | 4 +-
tests/mkcatalog/simple-3d.sh | 53 +++++
16 files changed, 1003 insertions(+), 400 deletions(-)
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index e438dda..f884620 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -438,41 +438,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/args.h b/bin/mkcatalog/args.h
index 30db17b..36a7ad5 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -446,6 +446,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "z",
+ UI_KEY_Z,
+ 0,
+ 0,
+ "Flux weighted center in third FITS axis.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"geox",
UI_KEY_GEOX,
0,
@@ -474,6 +488,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "geoz",
+ UI_KEY_GEOZ,
+ 0,
+ 0,
+ "Geometric center in third FITS axis.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"clumpsx",
UI_KEY_CLUMPSX,
0,
@@ -502,6 +530,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "clumpsz",
+ UI_KEY_CLUMPSZ,
+ 0,
+ 0,
+ "Flux.wht center of all clumps in obj. (Z).",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"clumpsgeox",
UI_KEY_CLUMPSGEOX,
0,
@@ -530,6 +572,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "clumpsgeoz",
+ UI_KEY_CLUMPSGEOZ,
+ 0,
+ 0,
+ "Geometric center of all clumps in obj. (Z).",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"ra",
UI_KEY_RA,
0,
@@ -586,6 +642,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "w3",
+ UI_KEY_W3,
+ 0,
+ 0,
+ "Flux weighted center in third WCS axis.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"geow1",
UI_KEY_GEOW1,
0,
@@ -614,6 +684,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "geow3",
+ UI_KEY_GEOW2,
+ 0,
+ 0,
+ "Geometric center in third WCS axis.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"clumpsw1",
UI_KEY_CLUMPSW1,
0,
@@ -642,6 +726,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "clumpsw3",
+ UI_KEY_CLUMPSW3,
+ 0,
+ 0,
+ "Flux.wht center of all clumps in 3rd WCS.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"clumpsgeow1",
UI_KEY_CLUMPSGEOW1,
0,
@@ -670,6 +768,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "clumpsgeow3",
+ UI_KEY_CLUMPSGEOW3,
+ 0,
+ 0,
+ "Geometric center of all clumps in 3rd WCS.",
+ ARGS_GROUP_COLUMNS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"brightness",
UI_KEY_BRIGHTNESS,
0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index a3e08f4..bb2d9a1 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,27 +122,49 @@ 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");
+}
+
+
+
+
@@ -267,6 +251,70 @@ columns_wcs_preparation(struct mkcatalogparams *p)
+static void
+columns_sanity_check(struct mkcatalogparams *p)
+{
+ gal_list_i32_t *colcode;
+
+ /* If there is any columns that need WCS, the input image needs to have a
+ WCS in its headers. So before anything, we need to check if a WCS is
+ present or not. This can't be done after the initial setting of column
+ properties because the WCS-related columns use information that is
+ based on it (for units and names). */
+ columns_wcs_preparation(p);
+
+ /* Check for dimension-specific columns. */
+ switch(p->input->ndim)
+ {
+ case 2:
+ for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+ switch(colcode->v)
+ {
+ case UI_KEY_Z:
+ case UI_KEY_GEOZ:
+ case UI_KEY_CLUMPSZ:
+ case UI_KEY_CLUMPSGEOZ:
+ case UI_KEY_W3:
+ case UI_KEY_GEOW3:
+ case UI_KEY_CLUMPSW3:
+ case UI_KEY_CLUMPSGEOW3:
+ error(EXIT_FAILURE, 0, "%s (hdu %s) is a 2D dataset, so columns "
+ "relating to a third dimension cannot be requested",
+ p->inputname, p->cp.hdu);
+ }
+ break;
+
+ case 3:
+ for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+ switch(colcode->v)
+ {
+ case UI_KEY_SEMIMAJOR:
+ case UI_KEY_SEMIMINOR:
+ case UI_KEY_AXISRATIO:
+ case UI_KEY_POSITIONANGLE:
+ case UI_KEY_GEOSEMIMAJOR:
+ case UI_KEY_GEOSEMIMINOR:
+ case UI_KEY_GEOAXISRATIO:
+ case UI_KEY_GEOPOSITIONANGLE:
+ error(EXIT_FAILURE, 0, "columns requiring second moment "
+ "calculations (like semi-major, semi-minor, axis ratio "
+ "and position angle) are not yet supported for 3D "
+ "inputs");
+ }
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. MakeCatalog should not reach this point "
+ "when the input has %zu dimensions", __func__,
+ PACKAGE_BUGREPORT, p->input->ndim);
+ }
+}
+
+
+
+
+
/* Set the necessary parameters for each output column and allocate the
space necessary to keep the values. */
void
@@ -278,12 +326,8 @@ columns_define_alloc(struct mkcatalogparams *p)
char *name=NULL, *unit=NULL, *ocomment=NULL, *ccomment=NULL;
uint8_t otype=GAL_TYPE_INVALID, ctype=GAL_TYPE_INVALID, *oiflag, *ciflag;
- /* If there is any columns that need WCS, the input image needs to have a
- WCS in its headers. So before anything, we need to check if a WCS is
- present or not. This can't be done after the initial setting of column
- properties because the WCS-related columns use information that is
- based on it (for units and names). */
- columns_wcs_preparation(p);
+ /* Do a sanity check on the columns given the input dataset. */
+ columns_sanity_check(p);
/* Allocate the array for which intermediate parameters are
necessary. The basic issue is that higher-level calculations require a
@@ -399,7 +443,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 +457,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;
@@ -425,9 +469,23 @@ columns_define_alloc(struct mkcatalogparams *p)
ciflag[ CCOL_VY ] = 1;
break;
+ case UI_KEY_Z:
+ name = "Z";
+ unit = "pixel";
+ ocomment = "Flux weighted center (FITS axis 3).";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
+ disp_width = 10;
+ disp_precision = 3;
+ oiflag[ OCOL_VZ ] = 1;
+ ciflag[ CCOL_VZ ] = 1;
+ break;
+
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 +499,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;
@@ -453,9 +511,23 @@ columns_define_alloc(struct mkcatalogparams *p)
ciflag[ CCOL_GY ] = 1;
break;
+ case UI_KEY_GEOZ:
+ name = "GEO_Z";
+ unit = "pixel";
+ ocomment = "Geometric center (FITS axis 3).";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
+ disp_width = 10;
+ disp_precision = 3;
+ oiflag[ OCOL_GZ ] = 1;
+ ciflag[ CCOL_GZ ] = 1;
+ break;
+
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 +540,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;
@@ -479,9 +551,22 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_C_VY ] = 1;
break;
+ case UI_KEY_CLUMPSZ:
+ name = "CLUMPS_Z";
+ unit = "pixel";
+ ocomment = "Flux weighted center of clumps (FITS axis 3).";
+ ccomment = NULL;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_INVALID;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
+ disp_width = 10;
+ disp_precision = 3;
+ oiflag[ OCOL_C_VZ ] = 1;
+ break;
+
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 +579,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;
@@ -505,40 +590,80 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_C_GY ] = 1;
break;
+ case UI_KEY_CLUMPSGEOZ:
+ name = "CLUMPS_GEO_Z";
+ unit = "pixel";
+ ocomment = "Geometric center of clumps (FITS axis 3).";
+ ccomment = NULL;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_INVALID;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_FLOAT;
+ disp_width = 10;
+ disp_precision = 3;
+ oiflag[ OCOL_C_GZ ] = 1;
+ break;
+
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[ OCOL_VZ ] = 1;
+ oiflag[ CCOL_VX ] = 1;
+ oiflag[ CCOL_VY ] = 1;
+ oiflag[ CCOL_VZ ] = 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[ OCOL_VZ ] = 1;
+ oiflag[ CCOL_VX ] = 1;
+ oiflag[ CCOL_VY ] = 1;
+ oiflag[ CCOL_VZ ] = 1;
+ columns_alloc_radec(p);
+ break;
+
+ case UI_KEY_W3:
+ name = p->ctype[2];
+ unit = p->input->wcs->cunit[2];
+ ocomment = "Flux weighted center (WCS axis 3).";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT64;
+ ctype = GAL_TYPE_FLOAT64;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 13;
+ disp_precision = 7;
+ oiflag[ OCOL_VX ] = 1;
+ oiflag[ OCOL_VY ] = 1;
+ oiflag[ OCOL_VZ ] = 1;
+ oiflag[ CCOL_VX ] = 1;
+ oiflag[ CCOL_VY ] = 1;
+ oiflag[ CCOL_VZ ] = 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;
@@ -547,15 +672,17 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_GX ] = 1;
oiflag[ OCOL_GY ] = 1;
+ oiflag[ OCOL_GZ ] = 1;
ciflag[ CCOL_GX ] = 1;
ciflag[ CCOL_GY ] = 1;
+ ciflag[ CCOL_GZ ] = 1;
columns_alloc_georadec(p);
break;
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;
@@ -564,15 +691,36 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_GX ] = 1;
oiflag[ OCOL_GY ] = 1;
+ oiflag[ OCOL_GZ ] = 1;
+ ciflag[ CCOL_GX ] = 1;
+ ciflag[ CCOL_GY ] = 1;
+ ciflag[ CCOL_GZ ] = 1;
+ columns_alloc_georadec(p);
+ break;
+
+ case UI_KEY_GEOW3:
+ name = gal_checkset_malloc_cat("GEO_", p->ctype[2]);
+ unit = p->input->wcs->cunit[2];
+ ocomment = "Geometric center (WCS axis 3).";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT64;
+ ctype = GAL_TYPE_FLOAT64;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 13;
+ disp_precision = 7;
+ oiflag[ OCOL_GX ] = 1;
+ oiflag[ OCOL_GY ] = 1;
+ oiflag[ OCOL_GZ ] = 1;
ciflag[ CCOL_GX ] = 1;
ciflag[ CCOL_GY ] = 1;
+ ciflag[ CCOL_GZ ] = 1;
columns_alloc_georadec(p);
break;
case UI_KEY_CLUMPSW1:
name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[0]);
unit = p->input->wcs->cunit[0];
- ocomment = "Flux.wht center of all clumps in 1st WCS axis.";
+ ocomment = "Flux.wht center of all clumps (WCS axis 1).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -581,13 +729,14 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_C_VX ] = 1;
oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_VZ ] = 1;
columns_alloc_clumpsradec(p);
break;
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 2).";
ccomment = NULL;
otype = GAL_TYPE_FLOAT64;
ctype = GAL_TYPE_INVALID;
@@ -596,13 +745,30 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_C_VX ] = 1;
oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_VZ ] = 1;
+ columns_alloc_clumpsradec(p);
+ break;
+
+ case UI_KEY_CLUMPSW3:
+ name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[2]);
+ unit = p->input->wcs->cunit[2];
+ ocomment = "Flux.wht center of all clumps (WCS axis 3).";
+ ccomment = NULL;
+ otype = GAL_TYPE_FLOAT64;
+ ctype = GAL_TYPE_INVALID;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 15;
+ disp_precision = 7;
+ oiflag[ OCOL_C_VX ] = 1;
+ oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_VZ ] = 1;
columns_alloc_clumpsradec(p);
break;
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;
@@ -611,13 +777,14 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_C_GX ] = 1;
oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
columns_alloc_clumpsgeoradec(p);
break;
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;
@@ -626,6 +793,23 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_precision = 7;
oiflag[ OCOL_C_GX ] = 1;
oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
+ columns_alloc_clumpsgeoradec(p);
+ break;
+
+ case UI_KEY_CLUMPSGEOW3:
+ name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[2]);
+ unit = p->input->wcs->cunit[2];
+ ocomment = "Geometric center of all clumps (WCS axis 3).";
+ ccomment = NULL;
+ otype = GAL_TYPE_FLOAT64;
+ ctype = GAL_TYPE_INVALID;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 13;
+ disp_precision = 7;
+ oiflag[ OCOL_C_GX ] = 1;
+ oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
columns_alloc_clumpsgeoradec(p);
break;
@@ -1254,6 +1438,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)
@@ -1295,6 +1484,11 @@ columns_fill(struct mkcatalog_passparams *pp)
OCOL_VY, OCOL_GY);
break;
+ case UI_KEY_Z:
+ ((float *)colarr)[oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL,
+ OCOL_VZ, OCOL_GZ);
+ break;
+
case UI_KEY_GEOX:
((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
break;
@@ -1303,6 +1497,10 @@ columns_fill(struct mkcatalog_passparams *pp)
((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
break;
+ case UI_KEY_GEOZ:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_GZ], oi[OCOL_NUMALL] );
+ break;
+
case UI_KEY_CLUMPSX:
((float *)colarr)[oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
OCOL_C_VX, OCOL_C_GX);
@@ -1313,6 +1511,11 @@ columns_fill(struct mkcatalog_passparams *pp)
OCOL_C_VY, OCOL_C_GY);
break;
+ case UI_KEY_CLUMPSZ:
+ ((float *)colarr)[oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
+ OCOL_C_VZ, OCOL_C_GZ);
+ break;
+
case UI_KEY_CLUMPSGEOX:
((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_C_GX],
oi[OCOL_C_NUMALL] );
@@ -1323,32 +1526,51 @@ columns_fill(struct mkcatalog_passparams *pp)
oi[OCOL_C_NUMALL] );
break;
+ case UI_KEY_CLUMPSGEOZ:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_C_GZ],
+ oi[OCOL_C_NUMALL] );
+ break;
+
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);
+ case UI_KEY_W3:
+ 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);
+ if(p->input->ndim==3)
+ vo[2][oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMALL, OCOL_VZ,
+ OCOL_GZ);
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] );
+ case UI_KEY_GEOW3:
+ go[0][oind] = MKC_RATIO( oi[OCOL_GX], oi[OCOL_NUMALL] );
+ go[1][oind] = MKC_RATIO( oi[OCOL_GY], oi[OCOL_NUMALL] );
+ if(p->input->ndim==3)
+ go[2][oind] = MKC_RATIO( oi[OCOL_GZ], 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);
+ case UI_KEY_CLUMPSW3:
+ 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);
+ if(p->input->ndim==3)
+ vcc[2][oind] = POS_V_G(oi, OCOL_C_SUMWHT, OCOL_C_NUMALL,
+ OCOL_C_VZ, OCOL_C_GZ);
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] );
+ case UI_KEY_CLUMPSGEOW3:
+ 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] );
+ if(p->input->ndim==3)
+ gcc[2][oind] = MKC_RATIO( oi[OCOL_C_GZ], oi[OCOL_C_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
@@ -1478,6 +1700,11 @@ columns_fill(struct mkcatalog_passparams *pp)
CCOL_VY, CCOL_GY);
break;
+ case UI_KEY_Z:
+ ((float *)colarr)[cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
+ CCOL_VZ, CCOL_GZ);
+ break;
+
case UI_KEY_GEOX:
((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GX],
ci[CCOL_NUMALL] );
@@ -1488,18 +1715,30 @@ columns_fill(struct mkcatalog_passparams *pp)
ci[CCOL_NUMALL] );
break;
+ case UI_KEY_GEOZ:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GZ],
+ ci[CCOL_NUMALL] );
+ break;
+
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);
+ case UI_KEY_W3:
+ 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);
+ if(p->input->ndim==3)
+ vc[2][cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL, CCOL_VZ,
+ CCOL_GZ);
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] );
+ case UI_KEY_GEOW3:
+ gc[0][cind] = MKC_RATIO( ci[CCOL_GX], ci[CCOL_NUMALL] );
+ gc[1][cind] = MKC_RATIO( ci[CCOL_GY], ci[CCOL_NUMALL] );
+ if(p->input->ndim==3)
+ gc[2][cind] = MKC_RATIO( ci[CCOL_GZ], ci[CCOL_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
@@ -1517,7 +1756,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 +1846,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..73db98f 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -71,8 +71,10 @@ enum objectcols
OCOL_SUM, /* Sum of (value-sky) in object. */
OCOL_VX, /* Sum of (value-sky) * x. */
OCOL_VY, /* Sum of (value-sky) * y. */
+ OCOL_VZ, /* Sum of (value-sky) * z. */
OCOL_SX, /* Shift along X axis. */
OCOL_SY, /* Shift along Y axis. */
+ OCOL_SZ, /* Shift along Z axis. */
OCOL_VXX, /* Sum of (value-sky) * x * x. */
OCOL_VYY, /* Sum of (value-sky) * y * y. */
OCOL_VXY, /* Sum of (value-sky) * x * y. */
@@ -82,6 +84,7 @@ enum objectcols
OCOL_NUMWHT, /* Number of positive pixels used for wht. */
OCOL_GX, /* Geometric center of object in X. */
OCOL_GY, /* Geometric center of object in Y. */
+ OCOL_GZ, /* Geometric center of object in Z. */
OCOL_GXX, /* Second order geometric variable: X*X. */
OCOL_GYY, /* Second order geometric variable: Y*Y. */
OCOL_GXY, /* Second order geometric variable: X*Y. */
@@ -91,8 +94,10 @@ enum objectcols
OCOL_C_SUM, /* Brightness in object clumps. */
OCOL_C_VX, /* Sum of (value-sky)*x on clumps. */
OCOL_C_VY, /* Sum of (value-sky)*y on obj. clumps. */
+ OCOL_C_VZ, /* Sum of (value-sky)*z on obj. clumps. */
OCOL_C_GX, /* Geometric center of clumps in object X. */
OCOL_C_GY, /* Geometric center of clumps in object Y. */
+ OCOL_C_GZ, /* Geometric center of clumps in object Z. */
OCOL_C_SUMWHT, /* Sum of positive image pixels for wht. */
OCOL_C_NUMWHT, /* Num of positive image pixels for wht. */
@@ -105,6 +110,7 @@ enum clumpcols
CCOL_NUM, /* Area of this clump. */
CCOL_VX, /* Sum of (value-sky) * x. */
CCOL_VY, /* Sum of (value-sky) * y. */
+ CCOL_VZ, /* Sum of (value-sky) * z. */
CCOL_VXX, /* Sum of flux*x*x of this clump. */
CCOL_VYY, /* Sum of flux*y*y of this clump. */
CCOL_VXY, /* Sum of flux*x*y of this clump. */
@@ -117,6 +123,7 @@ enum clumpcols
CCOL_NUMWHT, /* Num of positive image pixels for wht. */
CCOL_GX, /* Geometric center of clump in X. */
CCOL_GY, /* Geometric center of clump in Y. */
+ CCOL_GZ, /* Geometric center of clump in Y. */
CCOL_GXX, /* Second order geometric moment. */
CCOL_GYY, /* Second order geometric moment. */
CCOL_GXY, /* Second order geometric moment. */
@@ -188,12 +195,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..8ce056d 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -171,15 +171,19 @@ mkcatalog_first_pass(struct mkcatalog_passparams *pp)
/* Do the general geometric (independent of pixel value)
calculations. */
oi[ OCOL_NUMALL ]++;
- oi[ OCOL_GX ] += c[1];
- oi[ OCOL_GY ] += c[0];
+ oi[ OCOL_GX ] += c[ndim-1];
+ oi[ OCOL_GY ] += c[ndim-2];
+ if(ndim==3)
+ oi[ OCOL_GZ ] += c[ndim-3];
oi[ OCOL_GXX ] += sc[1] * sc[1];
oi[ OCOL_GYY ] += sc[0] * sc[0];
oi[ OCOL_GXY ] += sc[1] * sc[0];
if(p->clumps && *C>0)
{
- oi[ OCOL_C_GX ] += c[1];
- oi[ OCOL_C_GY ] += c[0];
+ oi[ OCOL_C_GX ] += c[ndim-1];
+ oi[ OCOL_C_GY ] += c[ndim-2];
+ if(ndim==3)
+ oi[ OCOL_C_GZ ] += c[ndim-3];
}
@@ -217,8 +221,10 @@ mkcatalog_first_pass(struct mkcatalog_passparams *pp)
{
oi[ OCOL_NUMWHT ]++;
oi[ OCOL_SUMWHT ] += ss;
- oi[ OCOL_VX ] += ss * c[1];
- oi[ OCOL_VY ] += ss * c[0];
+ oi[ OCOL_VX ] += ss * c[ndim-1];
+ oi[ OCOL_VY ] += ss * c[ndim-2];
+ if(ndim==3)
+ oi[ OCOL_VZ ] += ss * c[ndim-3];
oi[ OCOL_VXX ] += ss * sc[1] * sc[1];
oi[ OCOL_VYY ] += ss * sc[0] * sc[0];
oi[ OCOL_VXY ] += ss * sc[1] * sc[0];
@@ -226,8 +232,10 @@ mkcatalog_first_pass(struct mkcatalog_passparams *pp)
{
oi[ OCOL_C_NUMWHT ]++;
oi[ OCOL_C_SUMWHT ] += ss;
- oi[ OCOL_C_VX ] += ss * c[1];
- oi[ OCOL_C_VY ] += ss * c[0];
+ oi[ OCOL_C_VX ] += ss * c[ndim-1];
+ oi[ OCOL_C_VY ] += ss * c[ndim-2];
+ if(ndim==3)
+ oi[ OCOL_C_VZ ] += ss * c[ndim-3];
}
}
}
@@ -308,8 +316,10 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
/* Geometric measurements (independent of pixel value). */
ci[ CCOL_NUMALL ]++;
- ci[ CCOL_GX ] += c[1];
- ci[ CCOL_GY ] += c[0];
+ ci[ CCOL_GX ] += c[ndim-1];
+ ci[ CCOL_GY ] += c[ndim-2];
+ if(ndim==3)
+ ci[ CCOL_GZ ] += c[ndim-3];
ci[ CCOL_GXX ] += sc[1] * sc[1];
ci[ CCOL_GYY ] += sc[0] * sc[0];
ci[ CCOL_GXY ] += sc[1] * sc[0];
@@ -328,8 +338,10 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
{
ci[ CCOL_NUMWHT ]++;
ci[ CCOL_SUMWHT ] += ss;
- ci[ CCOL_VX ] += ss * c[1];
- ci[ CCOL_VY ] += ss * c[0];
+ ci[ CCOL_VX ] += ss * c[ndim-1];
+ ci[ CCOL_VY ] += ss * c[ndim-2];
+ if(ndim==3)
+ ci[ CCOL_VZ ] += ss * c[ndim-3];
ci[ CCOL_VXX ] += ss * sc[1] * sc[1];
ci[ CCOL_VYY ] += ss * sc[0] * sc[0];
ci[ CCOL_VXY ] += ss * sc[1] * sc[0];
@@ -559,38 +571,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 +613,24 @@ 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_W3: c=p->wcs_vo->next->next; break;
+ case UI_KEY_GEOW1: c=p->wcs_go; break;
+ case UI_KEY_GEOW2: c=p->wcs_go->next; break;
+ case UI_KEY_GEOW3: c=p->wcs_go->next->next; break;
+ case UI_KEY_CLUMPSW1: c=p->wcs_vcc; break;
+ case UI_KEY_CLUMPSW2: c=p->wcs_vcc->next; break;
+ case UI_KEY_CLUMPSW3: c=p->wcs_vcc->next->next; break;
+ case UI_KEY_CLUMPSGEOW1: c=p->wcs_gcc; break;
+ case UI_KEY_CLUMPSGEOW2: c=p->wcs_gcc->next; break;
+ case UI_KEY_CLUMPSGEOW3: c=p->wcs_gcc->next->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 +645,18 @@ 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_W3: c=p->wcs_vc->next->next; break;
+ case UI_KEY_GEOW1: c=p->wcs_gc; break;
+ case UI_KEY_GEOW2: c=p->wcs_gc->next; break;
+ case UI_KEY_GEOW3: c=p->wcs_gc->next->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 +901,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..5c68f4b 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -383,7 +383,7 @@ ui_preparations_read_inputs(struct mkcatalogparams *p)
/* Currently MakeCatalog is only implemented for 2D images. */
- if(p->input->ndim!=2)
+ if(p->input->ndim!=2 && p->input->ndim!=3)
error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, MakeCatalog "
"currently only supports 2D inputs", p->inputname, p->cp.hdu,
p->input->ndim);
@@ -688,8 +688,10 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
/* For a check.
for(i=0;i<p->numobjects;++i)
- printf("%zu: (%zu, %zu) --> (%zu, %zu)\n", i+1, minmax[i*width],
- minmax[i*width+1], minmax[i*width+2], minmax[i*width+3]);
+ printf("%zu: (%zu, %zu, %zu) --> (%zu, %zu, %zu)\n", i+1, minmax[i*width],
+ minmax[i*width+1], minmax[i*width+2], minmax[i*width+3],
+ minmax[i*width+4], minmax[i*width+5]);
+ exit(0);
*/
/* Make the tiles. */
@@ -1010,28 +1012,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/ui.h b/bin/mkcatalog/ui.h
index 5173056..a3cd08b 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -29,7 +29,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Available letters for short options:
- f g k l u v w x y z
+ f g k l u v w x y
H J L W X Y
*/
enum option_keys_enum
@@ -49,6 +49,7 @@ enum option_keys_enum
UI_KEY_AREA = 'a',
UI_KEY_X = 'x',
UI_KEY_Y = 'y',
+ UI_KEY_Z = 'z',
UI_KEY_RA = 'r',
UI_KEY_DEC = 'd',
UI_KEY_BRIGHTNESS = 'b',
@@ -82,18 +83,25 @@ enum option_keys_enum
UI_KEY_WEIGHTAREA,
UI_KEY_GEOX,
UI_KEY_GEOY,
+ UI_KEY_GEOZ,
UI_KEY_CLUMPSX,
UI_KEY_CLUMPSY,
+ UI_KEY_CLUMPSZ,
UI_KEY_CLUMPSGEOX,
UI_KEY_CLUMPSGEOY,
+ UI_KEY_CLUMPSGEOZ,
UI_KEY_W1,
UI_KEY_W2,
+ UI_KEY_W3,
UI_KEY_GEOW1,
UI_KEY_GEOW2,
+ UI_KEY_GEOW3,
UI_KEY_CLUMPSW1,
UI_KEY_CLUMPSW2,
+ UI_KEY_CLUMPSW3,
UI_KEY_CLUMPSGEOW1,
UI_KEY_CLUMPSGEOW2,
+ UI_KEY_CLUMPSGEOW3,
UI_KEY_CLUMPSBRIGHTNESS,
UI_KEY_NORIVERBRIGHTNESS,
UI_KEY_CLUMPSMAGNITUDE,
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index 0068b8f..1e932a6 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -179,14 +179,18 @@ 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] = gsl_rng_uniform_int(pp->rng, dsize[d]-tile->dsize[d]-1);
+ rcoord[d] = ( (dsize[d]-tile->dsize[d])
+ ? gsl_rng_uniform_int(pp->rng,
+ dsize[d]-tile->dsize[d]-1)
+ : 0 );
/* Set the tile's new starting pointer. */
tile->array = gal_data_ptr_increment(p->input->array,
@@ -207,9 +211,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/astmkprof-3d.conf b/bin/mkprof/astmkprof-3d.conf
index 0bf19e5..2df9040 100644
--- a/bin/mkprof/astmkprof-3d.conf
+++ b/bin/mkprof/astmkprof-3d.conf
@@ -50,7 +50,7 @@
# WCS:
crpix 1,1,1
- crval 1,1,1
+ crval 1,1,1e-10
cdelt 0.2/3600,0.2/3600,1.25e-10
pc -1,0,0,0,1,0,0,0,1
cunit deg,deg,m
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 78cc20e..06dba85 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -1580,9 +1580,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
@@ -1590,30 +1591,47 @@ ui_finalize_coordinates(struct mkprofparams *p)
need to change them into actual image coordinates. */
if(p->mode==MKPROF_MODE_WCS)
{
- /* `gal_wcs_world_to_img' API needs to be changed to allow any number
- of dimensions. */
- if(ndim!=2)
- error(EXIT_FAILURE, 0, "%s: conversion from WCS coordinates is not "
- "yet implemented for %zu dimensions", __func__, ndim);
+ /* 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 = ndim==2 ? p->y : p->z; break;
+ case 1: arr = ndim==2 ? p->x : p->y; break;
+ case 2: arr = p->x; break;
+ default:
+ error(EXIT_FAILURE, 0, "conversion from WCS to image "
+ "coordinates is not 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);
- /* 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);
/* 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 2868fec..96c3bc3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -13667,12 +13667,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
@@ -13703,7 +13702,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}
@@ -13716,21 +13715,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
+@mymath{\overline{x^2}} 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
+@mymath{\overline{x^2}} 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
+@mymath{\sim10000}. 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
@@ -14306,22 +14308,31 @@ reported in this column (see @option{--geox}). You
can use
The flux weighted center of all objects and clumps along the second FITS
axis (vertical when viewed in SAO ds9). See @option{--x}.
+@item -z
+@itemx --z
+The flux weighted center of all objects and clumps along the third FITS
+axis. See @option{--x}.
+
@item --geox
-The geometric center of all objects and clumps along the first FITS
-axis axis. The geometric center is the average pixel positions
-irrespective of their pixel values.
+The geometric center of all objects and clumps along the first FITS axis
+axis. The geometric center is the average pixel positions irrespective of
+their pixel values.
@item --geoy
-The geometric center of all objects and clumps along the second FITS
-axis axis, see @option{--geox}.
+The geometric center of all objects and clumps along the second FITS axis
+axis, see @option{--geox}.
@item --clumpsx
-[Objects] The flux weighted center of all the clumps in this
-object along the first FITS axis. See @option{--x}.
+[Objects] The flux weighted center of all the clumps in this object along
+the first FITS axis. See @option{--x}.
@item --clumpsy
-[Objects] The flux weighted center of all the clumps in this
-object along the second FITS axis. See @option{--x}.
+[Objects] The flux weighted center of all the clumps in this object along
+the second FITS axis. See @option{--x}.
+
+@item --clumpsz
+[Objects] The flux weighted center of all the clumps in this object along
+the third FITS axis. See @option{--x}.
@item --clumpsgeox
[Objects] The geometric center of all the clumps in this object along
@@ -14331,6 +14342,10 @@ the first FITS axis. See @option{--geox}.
[Objects] The geometric center of all the clumps in this object along
the second FITS axis. See @option{--geox}.
+@item --clumpsgeoz
+[Objects] The geometric center of all the clumps in this object along
+the third FITS axis. See @option{--geoz}.
+
@item -r
@itemx --ra
Flux weighted right ascension of all objects or clumps, see
@@ -14359,6 +14374,11 @@ Flux weighted second WCS axis of all objects or
clumps, see
@option{--x}. The second WCS axis is commonly used as declination in
images.
+@item --w3
+Flux weighted third WCS axis of all objects or clumps, see
+@option{--x}. The third WCS axis is commonly used as wavelength in integral
+field unit data cubes.
+
@item --geow1
Geometric center in first WCS axis of all objects or clumps, see
@option{--geox}. The first WCS axis is commonly used as right ascension in
@@ -14369,26 +14389,41 @@ Geometric center in second WCS axis of all objects or
clumps, see
@option{--geox}. The second WCS axis is commonly used as declination in
images.
+@item --geow3
+Geometric center in third WCS axis of all objects or clumps, see
+@option{--geox}. The third WCS axis is commonly used as wavelength in
+integral field unit data cubes.
+
@item --clumpsw1
[Objects] Flux weighted center in first WCS axis of all clumps in this
object, see @option{--x}. The first WCS axis is commonly used as right
ascension in images.
@item --clumpsw2
-[Objects] Flux weighted declination of all clumps in this object, see
-@option{--x}. The second WCS axis is commonly used as declination in
-images.
+[Objects] Flux weighted center in third WCS axis of all clumps in this
+object, see @option{--x}. The second WCS axis is commonly used as
+declination in images.
+
+@item --clumpsw3
+[Objects] Flux weighted center in third WCS axis of all clumps in this
+object, see @option{--x}. The third WCS axis is commonly used as wavelength
+in integral field unit data cubes.
@item --clumpsgeow1
-[Objects] Geometric center right ascension of all clumps in this object,
+[Objects] Geometric center in first WCS axis of all clumps in this object,
see @option{--geox}. The first WCS axis is commonly used as right ascension
in images.
@item --clumpsgeow2
-[Objects] Geometric center declination of all clumps in this object, see
-@option{--geox}. The second WCS axis is commonly used as declination in
+[Objects] Geometric center in second WCS axis of all clumps in this object,
+see @option{--geox}. The second WCS axis is commonly used as declination in
images.
+@item --clumpsgeow3
+[Objects] Geometric center in third WCS axis of all clumps in this object,
+see @option{--geox}. The third WCS axis is commonly used as wavelength in
+integral field unit data cubes.
+
@item -b
@itemx --brightness
The brightness (sum of all pixel values), see @ref{Flux Brightness and
@@ -20452,25 +20487,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
-@deftypefun 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
-@code{*x==NULL} 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.
+@deftypefun {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,
+see@ref{Linked 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
+@code{coords} (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
-@deftypefun 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
-@code{gal_wcs_world_to_img} for more.
+@deftypefun {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..7c974f1 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -620,59 +620,166 @@ 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, %g, %g), [stat: %d]\n",
+ world[i*3], world[i*3+1 ], world[i*3+2],
+ pixcrd[i*3], pixcrd[i*3+1], pixcrd[i*3+2], 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 +787,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 +820,27 @@ 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, %g, %g), [stat: %d]\n",
+ pixcrd[i*3], pixcrd[i*3+1], pixcrd[i*3+2],
+ world[i*3], world[i*3+1], world[i*3+2], 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 +848,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;
}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 213f94e..2285ee2 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -118,9 +118,11 @@ if COND_FITS
fits/copyhdu.sh: fits/write.sh.log mkprof/mosaic2.sh.log
endif
if COND_MKCATALOG
- MAYBE_MKCATALOG_TESTS = mkcatalog/simple.sh mkcatalog/aperturephot.sh
+ MAYBE_MKCATALOG_TESTS = mkcatalog/simple.sh mkcatalog/simple-3d.sh \
+ mkcatalog/aperturephot.sh
mkcatalog/simple.sh: noisechisel/noisechisel.sh.log
+ mkcatalog/simple-3d.sh: noisechisel/noisechisel-3d.sh.log
mkcatalog/aperturephot.sh: noisechisel/noisechisel.sh.log \
mkprof/clearcanvas.sh.log
endif
diff --git a/tests/mkcatalog/simple-3d.sh b/tests/mkcatalog/simple-3d.sh
new file mode 100755
index 0000000..e89b205
--- /dev/null
+++ b/tests/mkcatalog/simple-3d.sh
@@ -0,0 +1,53 @@
+# Make a simple catalog for NoiseChisel's output.
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <akhlaghi@gnu.org>
+# Contributing author(s):
+#
+# Copying and distribution of this file, with or without modification,
+# are permitted in any medium without royalty provided the copyright
+# notice and this notice are preserved. This file is offered as-is,
+# without any warranty.
+
+
+
+
+
+# Preliminaries
+# =============
+#
+# Set the variables (The executable is in the build tree). Do the
+# basic checks to see if the executable is made or if the defaults
+# file exists (basicchecks.sh is in the source tree).
+prog=mkcatalog
+execname=../bin/$prog/ast$prog
+img=3d-cat_noised_labeled.fits
+
+
+
+
+
+# Skip?
+# =====
+#
+# If the dependencies of the test don't exist, then skip it. There are two
+# types of dependencies:
+#
+# - The executable was not made (for example due to a configure option),
+#
+# - The input data was not made (for example the test that created the
+# data file failed).
+if [ ! -f $execname ]; then echo "$execname not created."; exit 77; fi
+if [ ! -f $img ]; then echo "$img does not exist."; exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+$execname $img --x --y --z --w1 --w2 --w3 --area --brightness --sn \
+ --upperlimit --tableformat=txt
- [gnuastro-commits] master 81f3f65 023/113: More --coordcol options acceptable in Crop, (continued)
- [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
- [gnuastro-commits] master 8372486 020/113: NoiseChisel's detection complete in 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6cc3d25 027/113: No 3D projections in function to inspect NoiseChisel outputs, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e990860 029/113: Merged recent updates in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6bdc5d6 030/113: Oversegmentation connectivity one less for 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master be3bfd8 033/113: NoiseChisel configuration file in 3D updated, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 3b6c15d 036/113: Merged with recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master fb3660f 037/113: MakeCatalog works in 3D,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 73fdf0c 006/113: MakeProfiles builds 3D ellipsoids, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 03d73fe 011/113: Some minor corrections in code and comments, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 87ab805 014/113: Identifiers for integer constants in BZERO comparisons, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 1483201 009/113: Noised cube created in make check, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master bdeaba9 012/113: Corrected name of 3D catalog for tarball inclusion, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6b53c5e 015/113: Bug fixes in master: commits 1ff1c25 and 540af65, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master adeb1c6 016/113: NoiseChisel's convolution step in 3D completed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master ace4160 017/113: New --mcolisbrightness option for MakeProfiles, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master f06acc5 018/113: Merged recent work from master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 61aec7b 019/113: Convolve name and dimensionality checks, Mohammad Akhlaghi, 2021/04/16