[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 63d062a: 3D features in Convolve, Match, MakeC
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 63d062a: 3D features in Convolve, Match, MakeCatalog, MakeProfiles and library |
|
Date: |
Wed, 30 Oct 2019 14:58:24 -0400 (EDT) |
branch: master
commit 63d062ac3194574acbb5d68f025240ffe968d086
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
3D features in Convolve, Match, MakeCatalog, MakeProfiles and library
After some time of working on a development branch on 3D datasets, with
this commit, I am importing the finished 3D work in the programs listed
above and the library into the `master' branch. NoiseChisel and Segment are
still under development, so I haven't brought them in yet.
---
NEWS | 37 +++
bin/convolve/ui.c | 16 +-
bin/match/args.h | 2 +-
bin/match/ui.c | 140 ++++++++++-
bin/mkcatalog/args.h | 181 ++++++++++++++
bin/mkcatalog/columns.c | 425 ++++++++++++++++++++++++++++++---
bin/mkcatalog/main.h | 17 +-
bin/mkcatalog/mkcatalog.c | 135 ++++++++---
bin/mkcatalog/mkcatalog.h | 1 +
bin/mkcatalog/parse.c | 508 +++++++++++++++++++++++++++++++++++++---
bin/mkcatalog/ui.c | 145 +++++++++++-
bin/mkcatalog/ui.h | 13 +
bin/mkcatalog/upperlimit.c | 85 +++++--
bin/mkprof/Makefile.am | 2 +-
bin/mkprof/args.h | 43 +++-
bin/mkprof/astmkprof-3d.conf | 57 +++++
bin/mkprof/main.h | 17 +-
bin/mkprof/mkprof.c | 69 +++++-
bin/mkprof/mkprof.h | 20 +-
bin/mkprof/oneprofile.c | 102 ++++++--
bin/mkprof/ui.c | 431 +++++++++++++++++++++++++++++-----
bin/mkprof/ui.h | 3 +
bin/warp/ui.c | 2 +-
doc/gnuastro.texi | 321 +++++++++++++++++++++++--
lib/binary.c | 49 +++-
lib/box.c | 176 +++++++++++++-
lib/data.c | 1 -
lib/gnuastro-internal/options.h | 1 +
lib/gnuastro/box.h | 7 +
lib/label.c | 6 +-
lib/match.c | 208 ++++++++++++----
lib/options.c | 1 -
lib/wcs.c | 10 +-
tests/Makefile.am | 26 +-
tests/mkcatalog/simple-3d.sh | 53 +++++
tests/mknoise/addnoise-3d.sh | 59 +++++
tests/mkprof/3d-cat.sh | 48 ++++
tests/mkprof/3d-cat.txt | 5 +
tests/mkprof/3d-kernel.sh | 47 ++++
39 files changed, 3126 insertions(+), 343 deletions(-)
diff --git a/NEWS b/NEWS
index 7cb6478..53cce0e 100644
--- a/NEWS
+++ b/NEWS
@@ -26,11 +26,46 @@ See the end of the file for license conditions.
the default colormap of Python's Matplotlib, and is available in many
other plotting tools like LaTeX's PGFPlots.
+ Convolve:
+ - Spatial domain convolution now possible on 3D data cubes (with a 3D
+ kernel).
+
CosmicCalculator:
--lineatz: return the observed wavelength of a line if it was emitted at
the redshift given to CosmicCalculator. You can either use known line
names, or directly give a number as any emitted line's wavelength.
+ MakeCatalog:
+ - Catalogs from 3D inputs now available, with the following new options,
+ see book for more.
+ --spectrum: label's spectrum (across the third dimension).
+ --z: Flux-weighted position in 3rd dimension.
+ --geoz: Geometric center in third FITS axis.
+ --minz: Minimum third FITS axis position.
+ --maxz: Maximum third FITS axis position.
+ --clumpsz: Flux weighted center of all clumps in object in 3rd dim.
+ --clumpsgeoz: Geometric center of all clumps in object in 3rd dim.
+ --w3: Flux weighted center in third WCS axis.
+ --geow3: Geometric center in third WCS axis.
+ --clumpsw3: Flux wheighted center of all clumps in 3rd dim.
+ --clumpsgeow3: Geometric center of all clumps in 3rd dim.
+ --areaxy: Projected area in first two dimentions.
+ --geoareaxy: Projected geoarea in first two dimensions.
+
+ MakeProfiles:
+ - Can produce mock ellipsoids in a datacube (using X-Z-X Euler angles
+ for 3D orientation), the following options have been added, see the
+ book for more details.
+ --p2col: Second Euler angle (X-Z-X order).
+ --p3col: Third Euler angle (X-Z-X order).
+ --q2col: Axis ratio (major/dim3 in 3D).
+ - The `--kernel' option can build 3D kernels, see the description of
+ this option in the book for examples and details on how to run it.
+
+ Match:
+ - Matching of catalogs now possible using 3 coordinates (on catalogs
+ generated from 3D data cubes), see book for more.
+
Statistics:
--contour: compute a contour plot which can be directly fed into the
PGFPlots package of LaTeX for plotting the contours. Support for more
@@ -50,6 +85,8 @@ See the end of the file for license conditions.
Library:
- gal_binary_connected_indexs: store indexs of connected components.
+ - gal_box_bound_ellipsoid_extent: Extent of 3D ellipsoid.
+ - gal_box_bound_ellipsoid: Bounding box for a 3D ellipsoid.
** Removed features
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 6722c69..c6c4f2e 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -503,16 +503,24 @@ ui_preparations(struct convolveparams *p)
ui_read_input(p);
- /* Currently Convolve only works on 1D and 2D datasets. */
- if(p->input->ndim>2)
+ /* Currently Convolve only works on 1D, 2D and 3D datasets. */
+ if(p->input->ndim>3)
error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions. Currently "
- "Convolve only operates on 1D (table column) and 2D (image) "
- "datasets", p->filename, cp->hdu, p->input->ndim);
+ "Convolve only operates on 1D (table column, spectrum), 2D "
+ "(image), and 3D (data cube) datasets", p->filename, cp->hdu,
+ p->input->ndim);
/* Domain-specific checks. */
if(p->domain==CONVOLVE_DOMAIN_FREQUENCY)
{
+ /* Check the dimensionality. */
+ if(p->input->ndim!=2)
+ error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions. Frequency "
+ "domain convolution currently only operates on 2D images",
+ p->filename, cp->hdu, p->input->ndim);
+
+ /* Blank values. */
if( gal_blank_present(p->input, 1) )
fprintf(stderr, "\n----------------------------------------\n"
"######## %s WARNING ########\n"
diff --git a/bin/match/args.h b/bin/match/args.h
index 657cf73..5dd2424 100644
--- a/bin/match/args.h
+++ b/bin/match/args.h
@@ -145,7 +145,7 @@ struct argp_option program_options[] =
{
"aperture",
UI_KEY_APERTURE,
- "FLT[,FLT[,FLT]]",
+ "FLT[,...]",
0,
"Acceptable aperture for matching.",
UI_GROUP_CATALOGMATCH,
diff --git a/bin/match/ui.c b/bin/match/ui.c
index d9eb99f..dacb0d3 100644
--- a/bin/match/ui.c
+++ b/bin/match/ui.c
@@ -388,7 +388,7 @@ ui_read_columns_aperture_2d(struct matchparams *p)
{
size_t apersize=3;
gal_data_t *newaper=NULL;
- double *naper, *oaper=p->aperture->array;
+ double *naper=NULL, *oaper=p->aperture->array;
/* A general sanity check: the first two elements of aperture cannot be
zero or negative. */
@@ -400,7 +400,7 @@ ui_read_columns_aperture_2d(struct matchparams *p)
"zero or negative");
/* Will be needed in more than one case. */
- if(p->aperture->size!=3)
+ if(p->aperture->size!=apersize)
{
newaper=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &apersize, NULL,
0, -1, 1, NULL, NULL, NULL);
@@ -447,6 +447,136 @@ ui_read_columns_aperture_2d(struct matchparams *p)
+/* The 3D aperture must finally have these elements.
+
+ p->aperture[0] = Major axis length.
+ p->aperture[1] = Second axis ratio.
+ p->aperture[2] = Third axis ratio.
+ p->aperture[3] = First ZXZ Euler angle rotation.
+ p->aperture[4] = Second ZXZ Euler angle rotation.
+ p->aperture[5] = Third ZXZ Euler angle rotation. */
+static void
+ui_read_columns_aperture_3d(struct matchparams *p)
+{
+ size_t apersize=6;
+ gal_data_t *newaper=NULL;
+ double *naper=NULL, *oaper=p->aperture->array;
+
+ /* A general sanity check: the first two elements of aperture cannot be
+ zero or negative. */
+ if( oaper[0]<=0 )
+ error(EXIT_FAILURE, 0, "the first value of `--aperture' cannot be "
+ "zero or negative");
+ if( p->aperture->size>2 && (oaper[1]<=0 || oaper[2]<=0) )
+ error(EXIT_FAILURE, 0, "the second and third values of `--aperture' "
+ "cannot be zero or negative");
+
+ /* Will be needed in more than one case. */
+ if(p->aperture->size!=apersize)
+ {
+ newaper=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &apersize, NULL,
+ 0, -1, 1, NULL, NULL, NULL);
+ naper=newaper->array;
+ }
+
+ /* Different based on */
+ switch(p->aperture->size)
+ {
+ case 1:
+ naper[0] = oaper[0];
+ naper[1] = naper[2] = 1;
+ naper[3] = naper[4] = naper[5] = 0;
+ break;
+
+ case 3:
+ /* Major axis is along the first dimension, no rotation necessary. */
+ if(oaper[0]>=oaper[1] && oaper[0]>=oaper[2])
+ {
+ naper[0] = oaper[0];
+ naper[1] = oaper[1] / oaper[0];
+ naper[2] = oaper[2] / oaper[0];
+ naper[3] = naper[4] = naper[5] = 0;
+ }
+
+ /* Major axis is along the second dimension. So we want `X' to be in
+ the direction of `y'. Therefore, just the first Eurler ZXZ
+ rotation is necessary by 90 degrees. Here is how the rotated
+ coordinates (X,Y,Z) look like (after one rotation about Z).
+
+ |z (before) Z| /Y
+ | | /
+ | ==> |/
+ .-------- y .------X
+ /
+ x /
+
+ You see how the major axis (X) now lies in the second original
+ direction (y). The length along `x' is now along `Y' and the third
+ hasn't changed. Note that we are talking about the semi-axis
+ lengths (a scalar, not a vector), so +- is irrelevant. */
+ else if(oaper[1]>=oaper[0] && oaper[1]>=oaper[2])
+ {
+ naper[0] = oaper[1];
+ naper[1] = oaper[0] / oaper[1];
+ naper[2] = oaper[2] / oaper[1];
+ naper[3] = 90;
+ naper[4] = naper[5] = 0;
+ }
+
+ /* The major axis is along the third dimension. So we want `X' to
+ point in the direction of `z'. To get to that point, we need 90
+ degree rotations in all three Euler ZXZ rotations.
+
+ |z (before) z'| /y' |y'' |X
+ | | / | |
+ | ==> |/ ==> | ==> |
+ .------ y .------x' .------x' Y------.
+ / / /
+ x / z''/ Z/
+
+ We thus see that the length along the second original direction
+ (y) doesn't change stays in the second direction (Y). But the
+ original first axis direction (x) is now represented by (Z). Note
+ that we are talking about the semi-axis lengths (a scalar, not a
+ vector), so +- is irrelevant. */
+ else
+ {
+ naper[0] = oaper[2];
+ naper[1] = oaper[1] / oaper[2];
+ naper[2] = oaper[0] / oaper[2];
+ naper[3] = naper[4] = naper[5] = 90;
+ }
+
+ break;
+
+ case 6:
+ if(oaper[1]>1 || oaper[2]>1)
+ error(EXIT_FAILURE, 0, "atleast one of the second or third values "
+ "to `--aperture' is larger than one. When size numbers are "
+ "given to this option (in 3D), the second and third are the "
+ "axis ratios (which must always be less than 1).");
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%zu values given to `--aperture'. In 3D, this "
+ "option can only take 1, 3, or 6 values. See the description of "
+ "this option in the book for more with this command:\n\n"
+ " $ info astmatch\n", p->aperture->size);
+ }
+
+ /* If a new aperture was defined, then replace it with the exitsting
+ one. */
+ if(newaper)
+ {
+ gal_data_free(p->aperture);
+ p->aperture=newaper;
+ }
+}
+
+
+
+
+
static size_t
ui_set_columns_sanity_check_read_aperture(struct matchparams *p)
{
@@ -489,10 +619,8 @@ ui_set_columns_sanity_check_read_aperture(struct
matchparams *p)
p->aperture->size);
break;
- case 2:
- ui_read_columns_aperture_2d(p);
- break;
-
+ case 2: ui_read_columns_aperture_2d(p); break;
+ case 3: ui_read_columns_aperture_3d(p); break;
default:
error(EXIT_FAILURE, 0, "%zu dimensional matches are not currently "
"supported (maximum is 2 dimensions). The number of "
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index a45ed91..1a2ad60 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -231,6 +231,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
+ {
+ "spectrum",
+ UI_KEY_SPECTRUM,
+ 0,
+ 0,
+ "Object spectrum for cube (3D) datasets.",
+ GAL_OPTIONS_GROUP_OUTPUT,
+ &p->spectrum,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
@@ -455,6 +468,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "z",
+ UI_KEY_Z,
+ 0,
+ 0,
+ "Flux weighted center in third FITS axis.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 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,
@@ -483,6 +510,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "geoz",
+ UI_KEY_GEOZ,
+ 0,
+ 0,
+ "Geometric center in third FITS axis.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"minx",
UI_KEY_MINX,
0,
@@ -539,6 +580,34 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "minz",
+ UI_KEY_MINZ,
+ 0,
+ 0,
+ "Minimum third FITS axis position.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "maxz",
+ UI_KEY_MAXZ,
+ 0,
+ 0,
+ "Maximum third FITS axis position.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 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,
@@ -567,6 +636,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).",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 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,
@@ -594,6 +677,20 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_column_codes_ll
},
+ {
+ "clumpsgeoz",
+ UI_KEY_CLUMPSGEOZ,
+ 0,
+ 0,
+ "Geometric center of all clumps in obj. (Z).",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
@@ -662,6 +759,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "w3",
+ UI_KEY_W3,
+ 0,
+ 0,
+ "Flux weighted center in third WCS axis.",
+ UI_GROUP_COLUMNS_POSITION_WCS,
+ 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,
@@ -690,6 +801,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "geow3",
+ UI_KEY_GEOW2,
+ 0,
+ 0,
+ "Geometric center in third WCS axis.",
+ UI_GROUP_COLUMNS_POSITION_WCS,
+ 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,
@@ -718,6 +843,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.",
+ UI_GROUP_COLUMNS_POSITION_WCS,
+ 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,
@@ -745,6 +884,20 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
ui_column_codes_ll
},
+ {
+ "clumpsgeow3",
+ UI_KEY_CLUMPSGEOW3,
+ 0,
+ 0,
+ "Geometric center of all clumps in 3rd WCS.",
+ UI_GROUP_COLUMNS_POSITION_WCS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
@@ -1075,6 +1228,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "areaxy",
+ UI_KEY_AREAXY,
+ 0,
+ 0,
+ "Projected area in first two dimentions.",
+ UI_GROUP_COLUMNS_MORPHOLOGY,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"clumpsarea",
UI_KEY_CLUMPSAREA,
0,
@@ -1117,6 +1284,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "geoareaxy",
+ UI_KEY_GEOAREAXY,
+ 0,
+ 0,
+ "Projected geoarea in first two dimensions.",
+ UI_GROUP_COLUMNS_MORPHOLOGY,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"semimajor",
UI_KEY_SEMIMAJOR,
0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index ac4640d..6fac3a7 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -182,7 +182,6 @@ columns_set_wcs_pointers(struct mkcatalogparams *p, double
***vo,
-
/******************************************************************/
/********** Column definition/allocation ***************/
/******************************************************************/
@@ -254,23 +253,88 @@ 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->objects->ndim)
+ {
+ case 2:
+ for(colcode=p->columnids; colcode!=NULL; colcode=colcode->next)
+ switch(colcode->v)
+ {
+ case UI_KEY_AREAXY:
+ case UI_KEY_GEOAREAXY:
+ case UI_KEY_Z:
+ case UI_KEY_MINZ:
+ case UI_KEY_MAXZ:
+ 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->objectsfile, 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->objects->ndim);
+ }
+}
+
+
+
+
+
/* Set the necessary parameters for each output column and allocate the
space necessary to keep the values. */
void
columns_define_alloc(struct mkcatalogparams *p)
{
gal_list_i32_t *colcode;
+ size_t ndim=p->objects->ndim;
gal_list_str_t *strtmp, *noclumpimg=NULL;
int disp_fmt=0, disp_width=0, disp_precision=0;
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
@@ -356,6 +420,19 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_NUM ] = ciflag[ CCOL_NUM ] = 1;
break;
+ case UI_KEY_AREAXY:
+ name = "AREAXY";
+ unit = "counter";
+ ocomment = "Projected valued pixels in first two dimensions.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_INT32;
+ ctype = GAL_TYPE_INT32;
+ disp_fmt = 0;
+ disp_width = 6;
+ disp_precision = 0;
+ oiflag[ OCOL_NUMXY ] = ciflag[ CCOL_NUMXY ] = 1;
+ break;
+
case UI_KEY_CLUMPSAREA:
name = "AREA_CLUMPS";
unit = "counter";
@@ -395,6 +472,19 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
break;
+ case UI_KEY_GEOAREAXY:
+ name = "AREA_FULL_XY";
+ unit = "counter";
+ ocomment = "Project number in first two dimensions.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_INT32;
+ ctype = GAL_TYPE_INT32;
+ disp_fmt = 0;
+ disp_width = 6;
+ disp_precision = 0;
+ oiflag[ OCOL_NUMALLXY ] = ciflag[ CCOL_NUMALLXY ] = 1;
+ break;
+
case UI_KEY_X:
name = "X";
unit = "pixel";
@@ -427,6 +517,20 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_NUMWHT ] = ciflag[ CCOL_NUMWHT ] = 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 = "pixel";
@@ -455,6 +559,20 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 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 = "pixel";
@@ -487,6 +605,19 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_C_NUMWHT ] = 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 = "pixel";
@@ -515,6 +646,18 @@ columns_define_alloc(struct mkcatalogparams *p)
oiflag[ OCOL_C_NUMALL ] = 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;
+
case UI_KEY_MINX:
name = "MIN_X";
unit = "pixel";
@@ -567,6 +710,32 @@ columns_define_alloc(struct mkcatalogparams *p)
ciflag[ CCOL_MAXY ] = 1;
break;
+ case UI_KEY_MINZ:
+ name = "MIN_Z";
+ unit = "pixel";
+ ocomment = "Minimum Z axis pixel position.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_UINT32;
+ ctype = GAL_TYPE_UINT32;
+ disp_fmt = 0;
+ disp_width = 10;
+ disp_precision = 0;
+ ciflag[ CCOL_MINZ ] = 1;
+ break;
+
+ case UI_KEY_MAXZ:
+ name = "MAX_Z";
+ unit = "pixel";
+ ocomment = "Maximum Z axis pixel position.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_UINT32;
+ ctype = GAL_TYPE_UINT32;
+ disp_fmt = 0;
+ disp_width = 10;
+ disp_precision = 0;
+ ciflag[ CCOL_MAXZ ] = 1;
+ break;
+
case UI_KEY_W1:
name = p->ctype[0];
unit = p->objects->wcs->cunit[0];
@@ -578,12 +747,17 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_radec(p);
- oiflag[ OCOL_VX ] = ciflag[ CCOL_VX ] = 1;
- oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
- oiflag[ OCOL_VY ] = ciflag[ CCOL_VY ] = 1;
- oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
- oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
- oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ oiflag[ OCOL_VX ] = ciflag[ CCOL_VX ] = 1;
+ oiflag[ OCOL_VY ] = ciflag[ CCOL_VY ] = 1;
+ oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
+ oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
+ oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ if(ndim==3)
+ {
+ oiflag[ OCOL_VZ ] = ciflag[ CCOL_VZ ] = 1;
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
+ }
break;
case UI_KEY_W2:
@@ -597,10 +771,36 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_radec(p);
+ oiflag[ OCOL_VX ] = ciflag[ CCOL_VX ] = 1;
+ oiflag[ OCOL_VY ] = ciflag[ CCOL_VY ] = 1;
+ oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
+ oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
+ oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ if(ndim==3)
+ {
+ oiflag[ OCOL_VZ ] = ciflag[ CCOL_VZ ] = 1;
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
+ }
+ break;
+
+ case UI_KEY_W3:
+ name = p->ctype[2];
+ unit = p->objects->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_FLOAT;
+ disp_width = 13;
+ disp_precision = 7;
+ columns_alloc_radec(p);
oiflag[ OCOL_VX ] = ciflag[ CCOL_VX ] = 1;
- oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
oiflag[ OCOL_VY ] = ciflag[ CCOL_VY ] = 1;
+ oiflag[ OCOL_VZ ] = ciflag[ CCOL_VZ ] = 1;
+ oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
oiflag[ OCOL_SUMWHT ] = ciflag[ CCOL_SUMWHT ] = 1;
oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
break;
@@ -616,9 +816,11 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_georadec(p);
- oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
- oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
- oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
+ oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ if(ndim==3)
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
break;
case UI_KEY_GEOW2:
@@ -632,8 +834,27 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_georadec(p);
+ oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
+ oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
+ if(ndim==3)
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
+ break;
+
+ case UI_KEY_GEOW3:
+ name = gal_checkset_malloc_cat("GEO_", p->ctype[2]);
+ unit = p->objects->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_FLOAT;
+ disp_width = 13;
+ disp_precision = 7;
+ columns_alloc_georadec(p);
oiflag[ OCOL_GX ] = ciflag[ CCOL_GX ] = 1;
oiflag[ OCOL_GY ] = ciflag[ CCOL_GY ] = 1;
+ oiflag[ OCOL_GZ ] = ciflag[ CCOL_GZ ] = 1;
oiflag[ OCOL_NUMALL ] = ciflag[ CCOL_NUMALL ] = 1;
break;
@@ -648,12 +869,17 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_clumpsradec(p);
- oiflag[ OCOL_C_VX ] = 1;
- oiflag[ OCOL_C_VY ] = 1;
- oiflag[ OCOL_C_GX ] = 1;
- oiflag[ OCOL_C_GY ] = 1;
- oiflag[ OCOL_C_SUMWHT ] = 1;
- oiflag[ OCOL_C_NUMALL ] = 1;
+ oiflag[ OCOL_C_VX ] = 1;
+ oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_GX ] = 1;
+ oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_SUMWHT ] = 1;
+ oiflag[ OCOL_C_NUMALL ] = 1;
+ if(ndim==3)
+ {
+ oiflag[ OCOL_C_VZ ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
+ }
break;
case UI_KEY_CLUMPSW2:
@@ -667,12 +893,41 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 15;
disp_precision = 7;
columns_alloc_clumpsradec(p);
+ oiflag[ OCOL_C_VX ] = 1;
+ oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_GX ] = 1;
+ oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_SUMWHT ] = 1;
+ oiflag[ OCOL_C_NUMALL ] = 1;
+ if(ndim==3)
+ {
+ oiflag[ OCOL_C_VZ ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
+ }
+ break;
+
+ case UI_KEY_CLUMPSW3:
+ name = gal_checkset_malloc_cat("CLUMPS_", p->ctype[2]);
+ unit = p->objects->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_FLOAT;
+ disp_width = 15;
+ disp_precision = 7;
+ columns_alloc_clumpsradec(p);
oiflag[ OCOL_C_VX ] = 1;
oiflag[ OCOL_C_VY ] = 1;
+ oiflag[ OCOL_C_VZ ] = 1;
oiflag[ OCOL_C_GX ] = 1;
oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
oiflag[ OCOL_C_SUMWHT ] = 1;
oiflag[ OCOL_C_NUMALL ] = 1;
+ if(ndim==3)
+ {
+ }
break;
case UI_KEY_CLUMPSGEOW1:
@@ -686,9 +941,11 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_clumpsgeoradec(p);
- oiflag[ OCOL_C_GX ] = 1;
- oiflag[ OCOL_C_GY ] = 1;
- oiflag[ OCOL_C_NUMALL ] = 1;
+ oiflag[ OCOL_C_GX ] = 1;
+ oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_NUMALL ] = 1;
+ if(ndim==3)
+ oiflag[ OCOL_C_GZ ] = 1;
break;
case UI_KEY_CLUMPSGEOW2:
@@ -702,8 +959,27 @@ columns_define_alloc(struct mkcatalogparams *p)
disp_width = 13;
disp_precision = 7;
columns_alloc_clumpsgeoradec(p);
+ oiflag[ OCOL_C_GX ] = 1;
+ oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_NUMALL ] = 1;
+ if(ndim==3)
+ oiflag[ OCOL_C_GZ ] = 1;
+ break;
+
+ case UI_KEY_CLUMPSGEOW3:
+ name = gal_checkset_malloc_cat("CLUMPS_GEO", p->ctype[2]);
+ unit = p->objects->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_FLOAT;
+ disp_width = 13;
+ disp_precision = 7;
+ columns_alloc_clumpsgeoradec(p);
oiflag[ OCOL_C_GX ] = 1;
oiflag[ OCOL_C_GY ] = 1;
+ oiflag[ OCOL_C_GZ ] = 1;
oiflag[ OCOL_C_NUMALL ] = 1;
break;
@@ -1447,6 +1723,7 @@ columns_clump_brightness(double *ci)
static uint32_t
columns_xy_extrema(struct mkcatalog_passparams *pp, size_t *coord, int key)
{
+ size_t ndim=pp->tile->ndim;
gal_data_t *tile=pp->tile, *block=tile->block;
/* We only want to do the coordinate estimation once: in `columns_fill',
@@ -1460,13 +1737,15 @@ columns_xy_extrema(struct mkcatalog_passparams *pp,
size_t *coord, int key)
block->ndim, block->dsize, coord);
/* Return the proper value: note that `coord' is in C standard: starting
- from the slowest dimension and counting from zero.*/
+ from the slowest dimension and counting from zero. */
switch(key)
{
- case UI_KEY_MINX: return coord[1] + 1; break;
- case UI_KEY_MAXX: return coord[1] + tile->dsize[1]; break;
- case UI_KEY_MINY: return coord[0] + 1; break;
- case UI_KEY_MAXY: return coord[0] + tile->dsize[0]; break;
+ case UI_KEY_MINX: return coord[ndim-1] + 1; break;
+ case UI_KEY_MAXX: return coord[ndim-1] + tile->dsize[ndim-1]; break;
+ case UI_KEY_MINY: return coord[ndim-2] + 1; break;
+ case UI_KEY_MAXY: return coord[ndim-2] + tile->dsize[ndim-2]; break;
+ case UI_KEY_MINZ: return coord[ndim-3] + 1; break;
+ case UI_KEY_MAXZ: return coord[ndim-3] + tile->dsize[ndim-3]; break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
"problem. The value %d is not a recognized value", __func__,
@@ -1538,7 +1817,7 @@ columns_fill(struct mkcatalog_passparams *pp)
void *colarr;
gal_data_t *column;
double *ci, *oi=pp->oi;
- size_t coord[2]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
+ size_t coord[3]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
size_t sr=pp->clumpstartindex, cind, coind;
size_t oind=pp->object-1; /* IDs start from 1, indexs from 0. */
@@ -1575,6 +1854,10 @@ columns_fill(struct mkcatalog_passparams *pp)
((int32_t *)colarr)[oind] = oi[OCOL_NUM];
break;
+ case UI_KEY_AREAXY:
+ ((int32_t *)colarr)[oind] = oi[OCOL_NUMXY];
+ break;
+
case UI_KEY_CLUMPSAREA:
((int32_t *)colarr)[oind] = oi[OCOL_C_NUM];
break;
@@ -1587,6 +1870,10 @@ columns_fill(struct mkcatalog_passparams *pp)
((int32_t *)colarr)[oind] = oi[OCOL_NUMALL];
break;
+ case UI_KEY_GEOAREAXY:
+ ((int32_t *)colarr)[oind] = oi[OCOL_NUMALLXY];
+ break;
+
case UI_KEY_X:
((float *)colarr)[oind] = POS_V_G(oi, OCOL_SUMWHT, OCOL_NUMWHT,
OCOL_VX, OCOL_GX);
@@ -1597,6 +1884,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;
@@ -1605,6 +1897,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_NUMWHT,
OCOL_C_VX, OCOL_C_GX);
@@ -1615,6 +1911,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] );
@@ -1625,39 +1926,59 @@ 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] );
+
case UI_KEY_MINX:
case UI_KEY_MAXX:
case UI_KEY_MINY:
case UI_KEY_MAXY:
+ case UI_KEY_MINZ:
+ case UI_KEY_MAXZ:
((uint32_t *)colarr)[oind]=columns_xy_extrema(pp, coord, key);
break;
case UI_KEY_W1:
case UI_KEY_W2:
+ 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->objects->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:
+ 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->objects->ndim==3)
+ go[2][oind] = MKC_RATIO( oi[OCOL_GZ], oi[OCOL_NUMALL] );
break;
case UI_KEY_CLUMPSW1:
case UI_KEY_CLUMPSW2:
+ 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->objects->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:
+ 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->objects->ndim==3)
+ gcc[2][oind] = MKC_RATIO( oi[OCOL_C_GZ], oi[OCOL_C_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
@@ -1817,6 +2138,10 @@ columns_fill(struct mkcatalog_passparams *pp)
((int32_t *)colarr)[cind]=ci[CCOL_NUM];
break;
+ case UI_KEY_AREAXY:
+ ((int32_t *)colarr)[cind]=ci[CCOL_NUMXY];
+ break;
+
case UI_KEY_WEIGHTAREA:
((int32_t *)colarr)[cind]=ci[CCOL_NUMWHT];
break;
@@ -1825,6 +2150,10 @@ columns_fill(struct mkcatalog_passparams *pp)
((int32_t *)colarr)[cind]=ci[CCOL_NUMALL];
break;
+ case UI_KEY_GEOAREAXY:
+ ((int32_t *)colarr)[cind]=ci[CCOL_NUMALLXY];
+ break;
+
case UI_KEY_X:
((float *)colarr)[cind] = POS_V_G(ci, CCOL_SUMWHT, CCOL_NUMALL,
CCOL_VX, CCOL_GX);
@@ -1835,6 +2164,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] );
@@ -1845,34 +2179,53 @@ columns_fill(struct mkcatalog_passparams *pp)
ci[CCOL_NUMALL] );
break;
- case UI_KEY_MINX:
+ case UI_KEY_GEOZ:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GZ],
+ ci[CCOL_NUMALL] );
+
+ case UI_KEY_MINX:
((uint32_t *)colarr)[cind] = ci[CCOL_MINX];
break;
- case UI_KEY_MAXX:
+ case UI_KEY_MAXX:
((uint32_t *)colarr)[cind] = ci[CCOL_MAXX];
break;
- case UI_KEY_MINY:
+ case UI_KEY_MINY:
((uint32_t *)colarr)[cind] = ci[CCOL_MINY];
break;
- case UI_KEY_MAXY:
+ case UI_KEY_MAXY:
((uint32_t *)colarr)[cind] = ci[CCOL_MAXY];
break;
+ case UI_KEY_MINZ:
+ ((uint32_t *)colarr)[cind] = ci[CCOL_MINZ];
+ break;
+
+ case UI_KEY_MAXZ:
+ ((uint32_t *)colarr)[cind] = ci[CCOL_MAXZ];
+ break;
+
case UI_KEY_W1:
case UI_KEY_W2:
+ 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->objects->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:
+ 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->objects->ndim==3)
+ gc[2][cind] = MKC_RATIO( ci[CCOL_GZ], ci[CCOL_NUMALL] );
break;
case UI_KEY_BRIGHTNESS:
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 0b7459b..00eba86 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -71,12 +71,15 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
enum objectcols
{
OCOL_NUMALL, /* Number of all pixels with this label. */
- OCOL_NUM, /* Number of pixels with a value. */
+ OCOL_NUMALLXY, /* Number of all pixels in first two dims. */
+ OCOL_NUM, /* Number of values used in this object. */
+ OCOL_NUMXY, /* Number of values in the first two dims. */
OCOL_SUM, /* Sum of (value-sky) in object. */
OCOL_SUM_VAR, /* Varience of sum (for brightness error). */
OCOL_MEDIAN, /* Median of value in object. */
OCOL_VX, /* Sum of (value-sky) * x. */
OCOL_VY, /* Sum of (value-sky) * y. */
+ OCOL_VZ, /* Sum of (value-sky) * z. */
OCOL_VXX, /* Sum of (value-sky) * x * x. */
OCOL_VYY, /* Sum of (value-sky) * y * y. */
OCOL_VXY, /* Sum of (value-sky) * x * y. */
@@ -86,6 +89,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. */
@@ -98,8 +102,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. */
@@ -109,7 +115,9 @@ enum objectcols
enum clumpcols
{
CCOL_NUMALL, /* Number of pixels in clump. */
+ CCOL_NUMALLXY, /* Number of pixels in first two dims. */
CCOL_NUM, /* Number of values used in clump. */
+ CCOL_NUMXY, /* Number of values only in first two dims. */
CCOL_SUM, /* River subtracted brightness. */
CCOL_SUM_VAR, /* Variance of sum (for brightness error). */
CCOL_MEDIAN, /* Median of values in clump. */
@@ -118,6 +126,7 @@ enum clumpcols
CCOL_RIV_SUM_VAR, /* Variance of sum (for error measurements). */
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. */
@@ -127,6 +136,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. */
@@ -134,6 +144,8 @@ enum clumpcols
CCOL_MAXX, /* Maximum X value of clump. */
CCOL_MINY, /* Minimum Y value of clump. */
CCOL_MAXY, /* Maximum Y value of clump. */
+ CCOL_MINZ, /* Minimum Z value of clump. */
+ CCOL_MAXZ, /* Maximum Z value of clump. */
CCOL_UPPERLIMIT_B, /* Upper limit brightness. */
CCOL_UPPERLIMIT_S, /* Upper limit one-sigma value. */
CCOL_UPPERLIMIT_Q, /* Quantile of object in random distribution.*/
@@ -171,6 +183,7 @@ struct mkcatalogparams
uint8_t subtractsky; /* ==1: subtract the Sky from values. */
float sfmagnsigma; /* Surface brightness multiple of sigma.*/
float sfmagarea; /* Surface brightness area (arcsec^2). */
+ uint8_t spectrum; /* Object spectrum for 3D datasets. */
char *upmaskfile; /* Name of upper limit mask file. */
char *upmaskhdu; /* HDU of upper limit mask file. */
@@ -213,6 +226,8 @@ struct mkcatalogparams
uint8_t uprangewarning; /* A warning must be printed. */
size_t *hostobjid_c; /* To sort the clumps table by Obj.ID. */
size_t *numclumps_c; /* To sort the clumps table by Obj.ID. */
+ gal_data_t *specsliceinfo; /* Slice information for spectra. */
+ gal_data_t *spectra; /* Array of datasets containing spectra.*/
char *usedvaluesfile; /* Ptr to final name used for values. */
char *usedclumpsfile; /* Ptr to final name used for clumps. */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index fab5d33..5d12634 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -42,6 +42,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/permutation.h>
#include <gnuastro-internal/timing.h>
+#include <gnuastro-internal/checkset.h>
#include "main.h"
#include "mkcatalog.h"
@@ -147,9 +148,10 @@ mkcatalog_single_object(void *in_prm)
{
/* For easy reading. Note that the object IDs start from one while
the array positions start from 0. */
- pp.ci = NULL;
- pp.object = tprm->indexs[i] + 1;
- pp.tile = &p->tiles[ tprm->indexs[i] ];
+ pp.ci = NULL;
+ pp.object = tprm->indexs[i] + 1;
+ pp.tile = &p->tiles[ tprm->indexs[i] ];
+ pp.spectrum = &p->spectra[ tprm->indexs[i] ];
/* Initialize the parameters for this object/tile. */
parse_initialize(&pp);
@@ -273,12 +275,16 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
{
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 into the output column. */
@@ -301,8 +307,10 @@ mkcatalog_wcs_conversion(struct mkcatalogparams *p)
{
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 into the output column. */
@@ -606,37 +614,84 @@ sort_clumps_by_objid(struct mkcatalogparams *p)
static void
mkcatalog_write_outputs(struct mkcatalogparams *p)
{
- /*char *str;*/
+ size_t i, scounter;
+ char str[200], *fname;
gal_list_str_t *comments;
+ int outisfits=gal_fits_name_is_fits(p->objectsout);
-
- /* OBJECT catalog */
- comments=mkcatalog_outputs_same_start(p, 0, "Detection");
-
- /* Reverse the comments list (so it is printed in the same order here),
- write the objects catalog and free the comments. */
- gal_list_str_reverse(&comments);
- gal_table_write(p->objectcols, comments, p->cp.tableformat,
- p->objectsout, "OBJECTS", 0);
- gal_list_str_free(comments, 1);
-
-
- /* CLUMPS catalog */
- if(p->clumps)
+ /* If a catalog is to be generated. */
+ if(p->objectcols)
{
- /* Make the comments. */
- comments=mkcatalog_outputs_same_start(p, 1, "Clumps");
+ /* OBJECT catalog */
+ comments=mkcatalog_outputs_same_start(p, 0, "Detection");
- /* Reverse the comments and write the catalog. */
+ /* Reverse the comments list (so it is printed in the same order
+ here), write the objects catalog and free the comments. */
gal_list_str_reverse(&comments);
- gal_table_write(p->clumpcols, comments, p->cp.tableformat,
- p->clumpsout, "CLUMPS", 0);
+ gal_table_write(p->objectcols, comments, p->cp.tableformat,
+ p->objectsout, "OBJECTS", 0);
gal_list_str_free(comments, 1);
+
+
+ /* CLUMPS catalog */
+ if(p->clumps)
+ {
+ /* Make the comments. */
+ comments=mkcatalog_outputs_same_start(p, 1, "Clumps");
+
+ /* Write objects catalog
+ ---------------------
+
+ Reverse the comments list (so it is printed in the same order
+ here), write the objects catalog and free the comments. */
+ gal_list_str_reverse(&comments);
+ gal_table_write(p->clumpcols, comments, p->cp.tableformat,
+ p->clumpsout, "CLUMPS", 0);
+ gal_list_str_free(comments, 1);
+ }
}
+ /* Spectra. */
+ if(p->spectra)
+ {
+ /* Inform the user (Writing many FITS extensions can take
+ long). */
+ if(p->objectcols && outisfits)
+ printf(" - Catalog(s) complete, writing spectra.\n");
+
+ /* Start counting and writing the files. Note that due to some
+ conditions (for example in debugging), a `p->spectra[i]' may not
+ actually contain any data. So we'll also count the number of
+ spectra that are written. */
+ scounter=0;
+ for(i=0;i<p->numobjects;++i)
+ if(p->spectra[i].ndim>0)
+ {
+ /* Increment the written spectra-counter. */
+ ++scounter;
+
+ /* Write the spectra based on the requested format. */
+ if(outisfits)
+ {
+ /* Write the table. */
+ sprintf(str, "SPECTRUM_%zu", i+1);
+ gal_table_write(&p->spectra[i], NULL, GAL_TABLE_FORMAT_BFITS,
+ p->objectsout, str, 0);
+ }
+ else
+ {
+ sprintf(str, "-spec-%zu.txt", i+1);
+ fname=gal_checkset_automatic_output(&p->cp, p->objectsout,
+ str);
+ gal_table_write(&p->spectra[i], NULL, GAL_TABLE_FORMAT_TXT,
+ fname, NULL, 0);
+ free(fname);
+ }
+ }
+ }
/* Configuration information. */
- if(gal_fits_name_is_fits(p->objectsout))
+ if(outisfits)
{
gal_fits_key_write_filename("input", p->objectsfile, &p->cp.okeys, 1);
gal_fits_key_write_config(&p->cp.okeys, "MakeCatalog configuration",
@@ -647,14 +702,34 @@ mkcatalog_write_outputs(struct mkcatalogparams *p)
/* Inform the user */
if(!p->cp.quiet)
{
- if(p->clumpsout && strcmp(p->clumpsout,p->objectsout))
+ if(p->objectcols)
{
- printf(" - Output objects catalog: %s\n", p->objectsout);
- if(p->clumps)
- printf(" - Output clumps catalog: %s\n", p->clumpsout);
+ if(p->clumpsout && strcmp(p->clumpsout,p->objectsout))
+ {
+ printf(" - Output objects catalog: %s\n", p->objectsout);
+ if(p->clumps)
+ printf(" - Output clumps catalog: %s\n", p->clumpsout);
+ }
+ else
+ printf(" - Catalog written to %s\n", p->objectsout);
+
+ }
+
+ if(p->spectra)
+ {
+ if(outisfits)
+ {
+ if(p->objectcols)
+ printf(" - Spectra in %zu extensions named `SPECTRUM_NN'.\n",
+ p->numobjects);
+ else
+ printf(" - Output: %s (Spectra in %zu extensions named "
+ "`SPECTRUM_NN').\n)", p->objectsout, p->numobjects);
+ }
+ else
+ printf(" - Spectra in %zu files with `-spec-NN.txt' suffix.\n",
+ p->numobjects);
}
- else
- printf(" - Catalog written to %s\n", p->objectsout);
}
}
diff --git a/bin/mkcatalog/mkcatalog.h b/bin/mkcatalog/mkcatalog.h
index a8de93b..730274b 100644
--- a/bin/mkcatalog/mkcatalog.h
+++ b/bin/mkcatalog/mkcatalog.h
@@ -41,6 +41,7 @@ struct mkcatalog_passparams
gsl_rng *rng; /* Random number generator. */
size_t clumpstartindex; /* Clump starting row in final catalog. */
gal_data_t *up_vals; /* Container for upper-limit values. */
+ gal_data_t *spectrum; /* Spectrum of each object. */
};
void
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index 879a7d1..1229de9 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -106,6 +106,331 @@ parse_initialize(struct mkcatalog_passparams *pp)
+static size_t *
+parse_spectrum_pepare(struct mkcatalog_passparams *pp, size_t *start_end_inc,
+ int32_t **st_o, float **st_v, float **st_std)
+{
+ size_t *tsize;
+ gal_data_t *spectile;
+ struct mkcatalogparams *p=pp->p;
+ size_t coord[3], minmax[6], numslices=p->objects->dsize[0];
+ gal_data_t *area, *sum, *esum, *proj, *eproj, *oarea, *osum, *oesum;
+
+ /* Get the coordinates of the spectral tile's starting element, then make
+ the tile. */
+ gal_dimension_index_to_coord(gal_pointer_num_between(p->objects->array,
+ pp->tile->array,
+ p->objects->type),
+ p->objects->ndim, p->objects->dsize, coord);
+ minmax[0]=0; /* Changed to first slice.*/
+ minmax[1]=coord[1];
+ minmax[2]=coord[2];
+ minmax[3]=p->objects->dsize[0]-1; /* Changed to last slice. */
+ minmax[4]=coord[1]+pp->tile->dsize[1]-1;
+ minmax[5]=coord[2]+pp->tile->dsize[2]-1;
+ spectile=gal_tile_series_from_minmax(p->objects, minmax, 1);
+
+ /* Find the starting (and ending) pointers on each of the datasets. */
+ *st_o = gal_tile_start_end_ind_inclusive(spectile, p->objects,
+ start_end_inc);
+ *st_v = (float *)(p->values->array) + start_end_inc[0];
+ *st_std = ( p->std
+ ? ( p->std->size==p->objects->size
+ ? (float *)(p->std->array) + start_end_inc[0]
+ : NULL )
+ : NULL );
+
+ /* Allocate the columns. */
+ area = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "AREA",
+ "counter", "Area of object in a slice");
+ sum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "SUM",
+ p->values->unit, "Sum of values with this label.");
+ esum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "SUM_ERR",
+ p->values->unit, "Error in SUM column.");
+ proj = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "SUM_PROJECTED",
+ p->values->unit, "Sum of full projected 2D area on "
+ "a slice.");
+ eproj = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap,
"SUM_PROJECTED_ERR",
+ p->values->unit, "Error in SUM_PROJECTED column.");
+ oarea = gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "AREA_OTHER",
+ "counter", "Area covered by other labels in a
slice.");
+ osum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "SUM_OTHER",
+ p->values->unit, "Sum of values in other labels on "
+ "a slice.");
+ oesum = gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 1,
+ p->cp.minmapsize, p->cp.quietmmap, "SUM_OTHER_ERR",
+ p->values->unit, "Error in SUM_OTHER column.");
+
+ /* Fill up the contents of the first element (note that the first
+ `gal_data_t' is actually in an array, so the skeleton is already
+ allocated, we just have to allocate its contents. */
+ gal_data_initialize(pp->spectrum, NULL, p->specsliceinfo->type, 1,
+ &numslices, NULL, 0, p->cp.minmapsize,
+ p->cp.quietmmap, NULL, NULL, NULL);
+ gal_data_copy_to_allocated(p->specsliceinfo, pp->spectrum);
+ pp->spectrum->next=gal_data_copy(p->specsliceinfo->next);
+
+
+ /* Add all the other columns in the final spectrum table. */
+ pp->spectrum->next->next = area;
+ area->next = sum;
+ area->next->next = esum;
+ area->next->next->next = proj;
+ area->next->next->next->next = eproj;
+ area->next->next->next->next->next = oarea;
+ area->next->next->next->next->next->next = osum;
+ area->next->next->next->next->next->next->next = oesum;
+
+ /* Clean up and return. */
+ tsize=spectile->dsize;
+ spectile->dsize=NULL;
+ gal_data_free(spectile);
+ return tsize;
+}
+
+
+
+
+
+/* Since spectra will be a large table for many objects, it is very
+ important to not consume too much space in for columns that don't need
+ it. This function will check the integer columns and if they are smaller
+ than the maximum values of smaller types, */
+static void
+parse_spectrum_uint32_to_best_type(gal_data_t **input)
+{
+ gal_data_t *tmp=gal_statistics_maximum(*input);
+
+ /* If maximum is smaller than UINT8_MAX, convert it to uint8_t */
+ if( *(uint32_t *)(tmp->array) < UINT8_MAX )
+ *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT8);
+
+ /* otherwise, if it is smaller than UINT16_MAX, convert it to uint16_t */
+ else if( *(uint32_t *)(tmp->array) < UINT16_MAX )
+ *input = gal_data_copy_to_new_type_free(*input, GAL_TYPE_UINT16);
+
+ /* Clean up. */
+ gal_data_free(tmp);
+}
+
+
+
+
+
+static void
+parse_spectrum_end(struct mkcatalog_passparams *pp, gal_data_t *xybin)
+{
+ size_t i;
+ double *searr, *pearr, *osearr;
+ struct mkcatalogparams *p=pp->p;
+
+ /* The datasets and their pointers. */
+ gal_data_t *area = pp->spectrum->next->next;
+ gal_data_t *sum = area->next;
+ gal_data_t *esum = area->next->next;
+ gal_data_t *proj = area->next->next->next;
+ gal_data_t *eproj = area->next->next->next->next;
+ gal_data_t *oarea = area->next->next->next->next->next;
+ gal_data_t *osum = area->next->next->next->next->next->next;
+ gal_data_t *oesum = area->next->next->next->next->next->next->next;
+
+ /* Apply corrections to the columns that need it. */
+ searr = esum->array;
+ pearr = eproj->array;
+ osearr = oesum->array;
+ for(i=0; i<p->objects->dsize[0]; ++i)
+ {
+ searr[i] = sqrt( searr[i] );
+ pearr[i] = sqrt( pearr[i] );
+ osearr[i] = sqrt( osearr[i] );
+ }
+
+ /* Convert the `double' type columns to `float'. The extra precision of
+ `double' was necessary when we were summing values in each slice. But
+ afterwards, it is not necessary at all (the measurement error is much
+ larger than a double-precision floating point number (15
+ decimals). But the extra space gained (double) is very useful in not
+ wasting too much memory and hard-disk space or online transfer time.*/
+ sum = gal_data_copy_to_new_type_free(sum, GAL_TYPE_FLOAT32);
+ esum = gal_data_copy_to_new_type_free(esum, GAL_TYPE_FLOAT32);
+ proj = gal_data_copy_to_new_type_free(proj, GAL_TYPE_FLOAT32);
+ eproj = gal_data_copy_to_new_type_free(eproj, GAL_TYPE_FLOAT32);
+ osum = gal_data_copy_to_new_type_free(osum, GAL_TYPE_FLOAT32);
+ oesum = gal_data_copy_to_new_type_free(oesum, GAL_TYPE_FLOAT32);
+
+ /* For the two area columns, find their maximum value and convert the
+ dataset to the smallest type that can hold them. */
+ parse_spectrum_uint32_to_best_type(&area);
+ parse_spectrum_uint32_to_best_type(&oarea);
+
+ /* List the datasets and write them into the pointer for this object
+ (exact copy of the statement in `parse_spectrum_pepare'). */
+ pp->spectrum->next->next = area;
+ area->next = sum;
+ area->next->next = esum;
+ area->next->next->next = proj;
+ area->next->next->next->next = eproj;
+ area->next->next->next->next->next = oarea;
+ area->next->next->next->next->next->next = osum;
+ area->next->next->next->next->next->next->next = oesum;
+}
+
+
+
+
+/* Each spectrum is a multi-column table (note that the slice counter and
+ wavelength are written in the end):
+
+ Column 3: Number of object pixels.
+ Column 4: Sum of object pixel values.
+ Column 5: Error in Column 2.
+ Column 6: Sum over all 2D projection over whole specturm.
+ Column 7: Error in Column 4.
+ Column 8: Area of other labels in this slice.
+ Column 9: Flux by other objects in projected area.
+ Column 10: Error in Column 9.
+ */
+static void
+parse_spectrum(struct mkcatalog_passparams *pp, gal_data_t *xybin)
+{
+ struct mkcatalogparams *p=pp->p;
+
+ gal_data_t *area;
+ float *st_v, *st_std;
+ uint32_t *narr, *oarr;
+ size_t nproj=0, *tsize, start_end_inc[2];
+ uint8_t *xybinarr = xybin ? xybin->array : NULL;
+ int32_t *O, *OO, *st_o, *objarr=p->objects->array;
+ size_t tid, *dsize=p->objects->dsize, num_increment=1;
+ double var, *sarr, *searr, *parr, *pearr, *osarr, *osearr;
+ size_t increment=0, pind=0, sind=0, ndim=p->objects->ndim, c[3];
+ float st, sval, *V=NULL, *ST=NULL, *std=p->std?p->std->array:NULL;
+
+ /* Prepare the columns to write in. */
+ tsize = parse_spectrum_pepare(pp, start_end_inc, &st_o, &st_v, &st_std);
+ area = pp->spectrum->next->next;
+ narr = area->array;
+ sarr = area->next->array;
+ searr = area->next->next->array;
+ parr = area->next->next->next->array;
+ pearr = area->next->next->next->next->array;
+ oarr = area->next->next->next->next->next->array;
+ osarr = area->next->next->next->next->next->next->array;
+ osearr = area->next->next->next->next->next->next->next->array;
+
+ /* If tile-id isn't necessary, set `tid' to a blank value. */
+ tid = (p->std && p->std->size>1 && st_std == NULL) ? 0 : GAL_BLANK_SIZE_T;
+
+ /* Parse each contiguous patch of memory covered by this object. */
+ while( start_end_inc[0] + increment <= start_end_inc[1] )
+ {
+ /* Set the contiguous range to parse. The pixel-to-pixel counting
+ along the fastest dimension will be done over the `O' pointer. */
+ if( p->values ) V = st_v + increment;
+ if( p->std && st_std ) ST = st_std + increment;
+ OO = ( O = st_o + increment ) + pp->tile->dsize[ndim-1];
+
+ /* Parse the tile. */
+ do
+ {
+ /* Only continue if this voxel is useful: it isn't NaN, or its
+ covered by the projected area or object's label. */
+ if( !isnan(*V) && (xybin && xybinarr[pind]==2) )
+ {
+ /* Get the error in measuing this pixel's flux. */
+ if(p->std)
+ {
+ /* If the standard deviation is given on a tile
+ structure, estimate the tile ID. */
+ if(tid != GAL_BLANK_SIZE_T)
+ {
+ gal_dimension_index_to_coord(O-objarr, ndim, dsize, c);
+ tid=gal_tile_full_id_from_coord(&p->cp.tl, c);
+ }
+
+ /* Get the error associated with this voxel. Note that if
+ we are given a variance dataset already, there is no
+ need to use `st*st', we can directly use `sval'. */
+ sval = st_std ? *ST : (p->std->size>1?std[tid]:std[0]);
+ st = p->variance ? sqrt(sval) : sval;
+ var = (p->variance ? sval : st*st) + fabs(*V);
+ }
+ else var = NAN;
+
+
+ /* Projected spectra: see if we have a value of `2' in the
+ `xybin' array (showing that there is atleast one non-blank
+ element there over the whole spectrum. */
+ ++nproj;
+ parr [ sind ] += *V;
+ pearr[ sind ] += var;
+
+ /* Calculate the number of labeled/detected pixels that
+ don't belong to this object. */
+ if(*O>0)
+ {
+ if(*O==pp->object)
+ {
+ ++narr[ sind ];
+ sarr [ sind ] += *V;
+ searr [ sind ] += var;
+ }
+ else
+ {
+ ++oarr [ sind ];
+ osarr [ sind ] += *V;
+ osearr [ sind ] += var;
+ }
+ }
+ }
+
+ /* Increment the pointers. */
+ if( xybin ) ++pind;
+ if( p->values ) ++V;
+ if( p->std && st_std ) ++ST;
+ }
+ while(++O<OO);
+
+ /* Increment to the next contiguous region of this tile. */
+ increment += ( gal_tile_block_increment(p->objects, tsize,
+ num_increment++, NULL) );
+
+ /* Increment the slice number, `sind', and reset the projection (2D)
+ index `pind' if we have just finished parsing a slice. */
+ if( (num_increment-1)%pp->tile->dsize[1]==0 )
+ {
+ /* If there was no measurement, set NaN for the values and their
+ errors (zero is meaningful). */
+ if( nproj ==0 ) parr[sind] = pearr[sind] = NAN;
+ if( narr[sind] ==0 ) sarr[sind] = searr[sind] = NAN;
+ if( oarr[sind] ==0 ) osarr[sind] = osearr[sind] = NAN;
+
+ nproj=pind=0;
+ ++sind;
+ }
+ }
+
+ /* Finalize the spectrum generation and clean up. */
+ parse_spectrum_end(pp, xybin);
+ free(tsize);
+
+ /* For a check.
+ gal_table_write(pp->spectrum, NULL, GAL_TABLE_FORMAT_BFITS,
+ "spectrum.fits", "SPECTRUM", 0);
+ */
+}
+
+
+
+
+
void
parse_objects(struct mkcatalog_passparams *pp)
{
@@ -114,9 +439,11 @@ parse_objects(struct mkcatalog_passparams *pp)
size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
double *oi=pp->oi;
+ gal_data_t *xybin=NULL;
size_t *tsize=pp->tile->dsize;
- size_t d, increment=0, num_increment=1;
+ uint8_t *xybinarr=NULL, *u, *uf;
float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
+ size_t d, pind=0, increment=0, num_increment=1;
int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
@@ -136,10 +463,13 @@ parse_objects(struct mkcatalog_passparams *pp)
/* Coordinate-related columns. */
( oif[ OCOL_GX ]
|| oif[ OCOL_GY ]
+ || oif[ OCOL_GZ ]
|| oif[ OCOL_VX ]
|| oif[ OCOL_VY ]
+ || oif[ OCOL_VZ ]
|| oif[ OCOL_C_GX ]
|| oif[ OCOL_C_GY ]
+ || oif[ OCOL_C_GZ ]
|| sc
/* When the sky and its STD are tiles, we'll also need
the coordinate to find which tile a pixel belongs
@@ -149,6 +479,19 @@ parse_objects(struct mkcatalog_passparams *pp)
: NULL );
+ /* If an XY projection area is necessary, we'll need to allocate an array
+ to keep the projected space. */
+ if( p->spectrum
+ || oif[ OCOL_NUMALLXY ]
+ || oif[ OCOL_NUMXY ] )
+ {
+ xybin=gal_data_alloc(NULL, GAL_TYPE_UINT8, 2, &tsize[1], NULL,
+ 1, p->cp.minmapsize, p->cp.quietmmap,
+ NULL, NULL, NULL);
+ xybinarr=xybin->array;
+ }
+
+
/* Parse each contiguous patch of memory covered by this object. */
while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
{
@@ -174,7 +517,8 @@ parse_objects(struct mkcatalog_passparams *pp)
/* Add to the area of this object. */
- if(oif[ OCOL_NUMALL ]) oi[ OCOL_NUMALL ]++;
+ if(xybin) xybinarr[ pind ]=1;
+ if(oif[ OCOL_NUMALL ]) oi[ OCOL_NUMALL ]++;
/* Geometric coordinate measurements. */
@@ -189,8 +533,9 @@ parse_objects(struct mkcatalog_passparams *pp)
/* Do the general geometric (independent of pixel value)
calculations. */
- if(oif[ OCOL_GX ]) oi[ OCOL_GX ] += c[1]+1;
- if(oif[ OCOL_GY ]) oi[ OCOL_GY ] += c[0]+1;
+ if(oif[ OCOL_GX ]) oi[ OCOL_GX ] += c[ ndim-1 ]+1;
+ if(oif[ OCOL_GY ]) oi[ OCOL_GY ] += c[ ndim-2 ]+1;
+ if(oif[ OCOL_GZ ]) oi[ OCOL_GZ ] += c[ ndim-3 ]+1;
if(pp->shift)
{
/* Calculate the shifted coordinates for second order
@@ -210,8 +555,9 @@ parse_objects(struct mkcatalog_passparams *pp)
if(p->clumps && *C>0)
{
if(oif[ OCOL_C_NUMALL ]) oi[ OCOL_C_NUMALL ]++;
- if(oif[ OCOL_C_GX ]) oi[ OCOL_C_GX ] += c[1]+1;
- if(oif[ OCOL_C_GY ]) oi[ OCOL_C_GY ] += c[0]+1;
+ if(oif[ OCOL_C_GX ]) oi[ OCOL_C_GX ] += c[ ndim-1 ]+1;
+ if(oif[ OCOL_C_GY ]) oi[ OCOL_C_GY ] += c[ ndim-2 ]+1;
+ if(oif[ OCOL_C_GZ ]) oi[ OCOL_C_GZ ] += c[ ndim-3 ]+1;
}
}
@@ -220,6 +566,7 @@ parse_objects(struct mkcatalog_passparams *pp)
if( p->values && !( p->hasblank && isnan(*V) ) )
{
/* General flux summations. */
+ if(xybin) xybinarr[ pind ]=2;
if(oif[ OCOL_NUM ]) oi[ OCOL_NUM ]++;
if(oif[ OCOL_SUM ]) oi[ OCOL_SUM ] += *V;
@@ -236,8 +583,9 @@ parse_objects(struct mkcatalog_passparams *pp)
{
if(oif[ OCOL_NUMWHT ]) oi[ OCOL_NUMWHT ]++;
if(oif[ OCOL_SUMWHT ]) oi[ OCOL_SUMWHT ] += *V;
- if(oif[ OCOL_VX ]) oi[ OCOL_VX ] += *V*(c[1]+1);
- if(oif[ OCOL_VY ]) oi[ OCOL_VY ] += *V*(c[0]+1);
+ if(oif[ OCOL_VX ]) oi[ OCOL_VX ] += *V*(c[ ndim-1 ]+1);
+ if(oif[ OCOL_VY ]) oi[ OCOL_VY ] += *V*(c[ ndim-2 ]+1);
+ if(oif[ OCOL_VZ ]) oi[ OCOL_VZ ] += *V*(c[ ndim-3 ]+1);
if(pp->shift)
{
oi[ OCOL_VXX ] += *V * sc[1] * sc[1];
@@ -249,9 +597,11 @@ parse_objects(struct mkcatalog_passparams *pp)
if(oif[ OCOL_C_NUMWHT ]) oi[ OCOL_C_NUMWHT ]++;
if(oif[ OCOL_C_SUMWHT ]) oi[ OCOL_C_SUMWHT ] += *V;
if(oif[ OCOL_C_VX ])
- oi[ OCOL_C_VX ] += *V * (c[1]+1);
+ oi[ OCOL_C_VX ] += *V * (c[ ndim-1 ]+1);
if(oif[ OCOL_C_VY ])
- oi[ OCOL_C_VY ] += *V * (c[0]+1);
+ oi[ OCOL_C_VY ] += *V * (c[ ndim-2 ]+1);
+ if(oif[ OCOL_C_VZ ])
+ oi[ OCOL_C_VZ ] += *V * (c[ ndim-3 ]+1);
}
}
}
@@ -290,6 +640,7 @@ parse_objects(struct mkcatalog_passparams *pp)
}
/* Increment the other pointers. */
+ if( xybin ) ++pind;
if( p->values ) ++V;
if( p->clumps ) ++C;
if( p->sky && pp->st_sky ) ++SK;
@@ -300,11 +651,43 @@ parse_objects(struct mkcatalog_passparams *pp)
/* Increment to the next contiguous region of this tile. */
increment += ( gal_tile_block_increment(p->objects, tsize,
num_increment++, NULL) );
+
+ /* If a 2D projection is requested, see if we should initialize (set
+ to zero) the projection-index (`pind') not. */
+ if(xybin && (num_increment-1)%tsize[1]==0 )
+ pind=0;
}
+ /* Write the projected area columns. */
+ if(xybin)
+ {
+ /* Any non-zero pixel must be set for NUMALLXY. */
+ uf=(u=xybin->array)+xybin->size;
+ do
+ if(*u)
+ {
+ if(oif[ OCOL_NUMALLXY ] ) oi[ OCOL_NUMALLXY ]++;
+ if(oif[ OCOL_NUMXY ] && *u==2 ) oi[ OCOL_NUMXY ]++;
+ }
+ while(++u<uf);
+
+ /* For a check on the projected 2D areas.
+ if(xybin && pp->object==2)
+ {
+ gal_fits_img_write(xybin, "xybin.fits", NULL, NULL);
+ exit(0);
+ }
+ */
+ }
+
+ /* Generate the Spectrum. */
+ if(p->spectrum)
+ parse_spectrum(pp, xybin);
+
/* Clean up. */
- if(c) free(c);
- if(sc) free(sc);
+ if(c) free(c);
+ if(sc) free(sc);
+ if(xybin) gal_data_free(xybin);
}
@@ -332,12 +715,13 @@ parse_clumps(struct mkcatalog_passparams *pp)
size_t ndim=p->objects->ndim, *dsize=p->objects->dsize;
double *ci, *cir;
- uint8_t *cif=p->ciflag;
+ gal_data_t *xybin=NULL;
size_t *tsize=pp->tile->dsize;
int32_t *O, *OO, *C=NULL, nlab;
+ uint8_t *u, *uf, *cif=p->ciflag;
float var, sval, *V=NULL, *SK=NULL, *ST=NULL;
- size_t i, ii, d, increment=0, num_increment=1;
size_t nngb=gal_dimension_num_neighbors(ndim);
+ size_t i, ii, d, pind=0, increment=0, num_increment=1;
int32_t *objects=p->objects->array, *clumps=p->clumps->array;
float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
@@ -353,14 +737,18 @@ parse_clumps(struct mkcatalog_passparams *pp)
: NULL );
/* If any coordinate columns are requested. */
- size_t *c = ( ( cif[ CCOL_GX ]
- || cif[ CCOL_GY ]
- || cif[ CCOL_VX ]
- || cif[ CCOL_VY ]
+ size_t *c = ( ( cif[ CCOL_GX ]
+ || cif[ CCOL_GY ]
+ || cif[ CCOL_GZ ]
+ || cif[ CCOL_VX ]
+ || cif[ CCOL_VY ]
+ || cif[ CCOL_VZ ]
|| cif[ CCOL_MINX ]
|| cif[ CCOL_MAXX ]
|| cif[ CCOL_MINY ]
|| cif[ CCOL_MAXY ]
+ || cif[ CCOL_MINZ ]
+ || cif[ CCOL_MAXZ ]
|| sc
|| tid==GAL_BLANK_SIZE_T )
? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
@@ -376,6 +764,18 @@ parse_clumps(struct mkcatalog_passparams *pp)
: NULL );
size_t *dinc = ngblabs ? gal_dimension_increment(ndim, dsize) : NULL;
+ /* If an XY projection area is requested, we'll need to allocate an array
+ to keep the projected space.*/
+ if( cif[ CCOL_NUMALLXY ]
+ || cif[ CCOL_NUMXY ] )
+ {
+ xybin=gal_data_array_calloc(pp->clumpsinobj);
+ for(i=0;i<pp->clumpsinobj;++i)
+ gal_data_initialize(&xybin[i], NULL, GAL_TYPE_UINT8, 2, &tsize[1],
+ NULL, 1, p->cp.minmapsize, p->cp.quietmmap,
+ NULL, NULL, NULL);
+ }
+
/* Parse each contiguous patch of memory covered by this object. */
while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
@@ -405,8 +805,11 @@ parse_clumps(struct mkcatalog_passparams *pp)
/* Add to the area of this object. */
if( cif[ CCOL_NUMALL ]
|| cif[ CCOL_MINX ] || cif[ CCOL_MAXX ]
- || cif[ CCOL_MINY ] || cif[ CCOL_MAXY ] )
+ || cif[ CCOL_MINY ] || cif[ CCOL_MAXY ]
+ || cif[ CCOL_MINZ ] || cif[ CCOL_MAXZ ] )
ci[ CCOL_NUMALL ]++;
+ if(cif[ CCOL_NUMALLXY ])
+ ((uint8_t *)(xybin[*C-1].array))[ pind ] = 1;
/* Raw-position related measurements. */
if(c)
@@ -415,10 +818,18 @@ parse_clumps(struct mkcatalog_passparams *pp)
gal_dimension_index_to_coord(O-objects, ndim, dsize, c);
/* Position extrema measurements. */
- if(cif[ CCOL_MINX ]) ci[CCOL_MINX]=CMIN(CCOL_MINX, 1);
- if(cif[ CCOL_MAXX ]) ci[CCOL_MAXX]=CMAX(CCOL_MAXX, 1);
- if(cif[ CCOL_MINY ]) ci[CCOL_MINY]=CMIN(CCOL_MINY, 0);
- if(cif[ CCOL_MAXY ]) ci[CCOL_MAXY]=CMAX(CCOL_MAXY, 0);
+ if(cif[ CCOL_MINX ])
+ ci[CCOL_MINX]=CMIN(CCOL_MINX, ndim-1);
+ if(cif[ CCOL_MAXX ])
+ ci[CCOL_MAXX]=CMAX(CCOL_MAXX, ndim-1);
+ if(cif[ CCOL_MINY ])
+ ci[CCOL_MINY]=CMIN(CCOL_MINY, ndim-2);
+ if(cif[ CCOL_MAXY ])
+ ci[CCOL_MAXY]=CMAX(CCOL_MAXY, ndim-2);
+ if(cif[ CCOL_MINZ ])
+ ci[CCOL_MINZ]=CMIN(CCOL_MINZ, ndim-3);
+ if(cif[ CCOL_MAXZ ])
+ ci[CCOL_MAXZ]=CMAX(CCOL_MAXZ, ndim-3);
/* If we need tile-ID, get the tile ID now. */
if(tid!=GAL_BLANK_SIZE_T)
@@ -426,8 +837,9 @@ parse_clumps(struct mkcatalog_passparams *pp)
/* General geometric (independent of pixel value)
calculations. */
- if(cif[ CCOL_GX ]) ci[ CCOL_GX ] += c[1]+1;
- if(cif[ CCOL_GY ]) ci[ CCOL_GY ] += c[0]+1;
+ if(cif[ CCOL_GX ]) ci[ CCOL_GX ] += c[ ndim-1 ]+1;
+ if(cif[ CCOL_GY ]) ci[ CCOL_GY ] += c[ ndim-2 ]+1;
+ if(cif[ CCOL_GZ ]) ci[ CCOL_GZ ] += c[ ndim-3 ]+1;
if(pp->shift)
{
/* Shifted coordinates for second order moments,
@@ -446,14 +858,20 @@ parse_clumps(struct mkcatalog_passparams *pp)
if( p->values && !( p->hasblank && isnan(*V) ) )
{
/* Fill in the necessary information. */
- if(cif[ CCOL_NUM ]) ci[ CCOL_NUM ]++;
- if(cif[ CCOL_SUM ]) ci[ CCOL_SUM ] += *V;
+ if(cif[ CCOL_NUM ]) ci[ CCOL_NUM ]++;
+ if(cif[ CCOL_SUM ]) ci[ CCOL_SUM ] += *V;
+ if(cif[ CCOL_NUMXY ])
+ ((uint8_t *)(xybin[*C-1].array))[ pind ] = 2;
if( *V > 0.0f )
{
if(cif[ CCOL_NUMWHT ]) ci[ CCOL_NUMWHT ]++;
if(cif[ CCOL_SUMWHT ]) ci[ CCOL_SUMWHT ] += *V;
- if(cif[ CCOL_VX ]) ci[ CCOL_VX ]+=*V*(c[1]+1);
- if(cif[ CCOL_VY ]) ci[ CCOL_VY ]+=*V*(c[0]+1);
+ if(cif[ CCOL_VX ])
+ ci[ CCOL_VX ] += *V * (c[ ndim-1 ]+1);
+ if(cif[ CCOL_VY ])
+ ci[ CCOL_VY ] += *V * (c[ ndim-2 ]+1);
+ if(cif[ CCOL_VZ ])
+ ci[ CCOL_VZ ] += *V * (c[ ndim-3 ]+1);
if(pp->shift)
{
ci[ CCOL_VXX ] += *V * sc[1] * sc[1];
@@ -556,6 +974,7 @@ parse_clumps(struct mkcatalog_passparams *pp)
/* Increment the other pointers. */
++C;
+ if( xybin ) ++pind;
if( p->values ) ++V;
if( p->sky && pp->st_sky ) ++SK;
if( p->std && pp->st_std ) ++ST;
@@ -565,13 +984,44 @@ parse_clumps(struct mkcatalog_passparams *pp)
/* Increment to the next contiguous region of this tile. */
increment += ( gal_tile_block_increment(p->objects, tsize,
num_increment++, NULL) );
+
+ /* If a 2D projection is requested, see if we should initialize (set
+ to zero) the projection-index (`pind') not. */
+ if(xybin && (num_increment-1) % tsize[1]==0 )
+ pind=0;
}
+
+ /* Write the projected area columns. */
+ if(xybin)
+ for(i=0;i<pp->clumpsinobj;++i)
+ {
+ /* Pointer to make things easier. */
+ ci=&pp->ci[ i * CCOL_NUMCOLS ];
+
+ /* Any non-zero pixel must be set for NUMALLXY. */
+ uf=(u=xybin[i].array)+xybin[i].size;
+ do
+ if(*u)
+ {
+ if(cif[ CCOL_NUMALLXY ] ) ci[ CCOL_NUMALLXY ]++;
+ if(cif[ CCOL_NUMXY ] && *u==2 ) ci[ CCOL_NUMXY ]++;
+ }
+ while(++u<uf);
+
+ /* For a check on the projected 2D areas. */
+ if(xybin && pp->object==2)
+ gal_fits_img_write(&xybin[i], "xybin.fits", NULL, NULL);
+
+ }
+
+
/* Clean up. */
if(c) free(c);
if(sc) free(sc);
if(dinc) free(dinc);
if(ngblabs) free(ngblabs);
+ if(xybin) gal_data_array_free(xybin, pp->clumpsinobj, 1);
}
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 8ade64c..fcdae2a 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -634,11 +634,17 @@ ui_read_labels(struct mkcatalogparams *p)
p->objects=gal_data_copy_to_new_type_free(p->objects, GAL_TYPE_INT32);
- /* Currently MakeCatalog is only implemented for 2D images. */
- if(p->objects->ndim!=2)
+ /* Currently MakeCatalog is only implemented for 2D images or 3D cubes. */
+ if(p->objects->ndim!=2 && p->objects->ndim!=3)
error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, MakeCatalog "
- "currently only supports 2D inputs", p->objectsfile, p->cp.hdu,
- p->objects->ndim);
+ "currently only supports 2D or 3D datasets", p->objectsfile,
+ p->cp.hdu, p->objects->ndim);
+
+ /* Make sure the `--spectrum' option is not given on a 2D image. */
+ if(p->spectrum && p->objects->ndim!=3)
+ error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions, but `--spectrum' "
+ "is currently only defined on 3D datasets", p->objectsfile,
+ p->cp.hdu, p->objects->ndim);
/* See if the total number of objects is given in the header keywords. */
keys[0].name="NUMLABS";
@@ -764,6 +770,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
size_t i;
if(p->upperlimit) *values=1;
+ if(p->spectrum) *values=*std=1;
/* Go over all the object columns. Note that the objects and clumps (if
the `--clumpcat' option is given) inputs are mandatory and it is not
@@ -773,12 +780,15 @@ ui_necessary_inputs(struct mkcatalogparams *p, int
*values, int *sky,
switch(i)
{
case OCOL_NUMALL: /* Only object labels. */ break;
+ case OCOL_NUMALLXY: /* Only object labels. */ break;
case OCOL_NUM: *values = 1; break;
+ case OCOL_NUMXY: *values = 1; break;
case OCOL_SUM: *values = 1; break;
case OCOL_SUM_VAR: *values = *std = 1; break;
case OCOL_MEDIAN: *values = 1; break;
case OCOL_VX: *values = 1; break;
case OCOL_VY: *values = 1; break;
+ case OCOL_VZ: *values = 1; break;
case OCOL_VXX: *values = 1; break;
case OCOL_VYY: *values = 1; break;
case OCOL_VXY: *values = 1; break;
@@ -788,6 +798,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
case OCOL_NUMWHT: *values = 1; break;
case OCOL_GX: /* Only object labels. */ break;
case OCOL_GY: /* Only object labels. */ break;
+ case OCOL_GZ: /* Only object labels. */ break;
case OCOL_GXX: /* Only object labels. */ break;
case OCOL_GYY: /* Only object labels. */ break;
case OCOL_GXY: /* Only object labels. */ break;
@@ -800,8 +811,10 @@ ui_necessary_inputs(struct mkcatalogparams *p, int
*values, int *sky,
case OCOL_C_SUM: *values = 1; break;
case OCOL_C_VX: *values = 1; break;
case OCOL_C_VY: *values = 1; break;
+ case OCOL_C_VZ: *values = 1; break;
case OCOL_C_GX: /* Only clump labels. */ break;
case OCOL_C_GY: /* Only clump labels. */ break;
+ case OCOL_C_GZ: /* Only clump labels. */ break;
case OCOL_C_SUMWHT: *values = 1; break;
case OCOL_C_NUMWHT: *values = 1; break;
default:
@@ -818,7 +831,9 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
switch(i)
{
case CCOL_NUMALL: /* Only clump labels. */ break;
+ case CCOL_NUMALLXY: /* Only clump labels. */ break;
case CCOL_NUM: *values = 1; break;
+ case CCOL_NUMXY: *values = 1; break;
case CCOL_SUM: *values = 1; break;
case CCOL_SUM_VAR: *values = *std = 1; break;
case CCOL_MEDIAN: *values = 1; break;
@@ -827,6 +842,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
case CCOL_RIV_SUM_VAR: *values = *std = 1; break;
case CCOL_VX: *values = 1; break;
case CCOL_VY: *values = 1; break;
+ case CCOL_VZ: *values = 1; break;
case CCOL_VXX: *values = 1; break;
case CCOL_VYY: *values = 1; break;
case CCOL_VXY: *values = 1; break;
@@ -836,6 +852,7 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
case CCOL_NUMWHT: *values = 1; break;
case CCOL_GX: /* Only clump labels. */ break;
case CCOL_GY: /* Only clump labels. */ break;
+ case CCOL_GZ: /* Only clump labels. */ break;
case CCOL_GXX: /* Only clump labels. */ break;
case CCOL_GYY: /* Only clump labels. */ break;
case CCOL_GXY: /* Only clump labels. */ break;
@@ -843,6 +860,8 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
case CCOL_MAXX: /* Only clump labels. */ break;
case CCOL_MINY: /* Only clump labels. */ break;
case CCOL_MAXY: /* Only clump labels. */ break;
+ case CCOL_MINZ: /* Only clump labels. */ break;
+ case CCOL_MAXZ: /* Only clump labels. */ break;
case CCOL_UPPERLIMIT_B: *values = 1; break;
case CCOL_UPPERLIMIT_S: *values = 1; break;
case CCOL_UPPERLIMIT_Q: *values = 1; break;
@@ -1395,6 +1414,91 @@ ui_one_tile_per_object(struct mkcatalogparams *p)
+/* When a spectrum is requested, the slice information (slice number and
+ slice WCS) is common to all different spectra. So instead of calculating
+ it every time, we'll just make it once here, then copy it for every
+ object.
+
+ The Slice information is going to be written in every spectrum. So we
+ don't want it to take too much space. Therefore, only when the number of
+ slices is less than 65535 (2^16-1), will we actually use a 32-bit
+ integer type for the slice number column.
+*/
+static void
+ui_preparations_spectrum_wcs(struct mkcatalogparams *p)
+{
+ double *xarr, *yarr, *zarr;
+ gal_data_t *x, *y, *z, *coords;
+ size_t i, numslices=p->objects->dsize[0];
+ size_t slicenumtype=numslices>=65535 ? GAL_TYPE_UINT32 : GAL_TYPE_UINT16;
+
+ /* A small sanity check. */
+ if(p->objects->ndim!=3)
+ error(EXIT_FAILURE, 0, "%s (hdu %s) is a %zuD dataset, but `--spectrum' "
+ "is currently only defined on 3D datasets", p->objectsfile,
+ p->cp.hdu, p->objects->ndim);
+
+ /* Allocate space for the slice number as well as the X and Y positions
+ for WCS conversion. Note that the `z' axis is going to be converted to
+ WCS later, so we'll just give it the basic information now.*/
+ x=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, NULL, NULL, NULL);
+ y=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, NULL, NULL, NULL);
+ z=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numslices, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, p->ctype[2],
+ p->objects->wcs->cunit[2], "Slice WCS coordinates.");
+
+ /* Write values into the 3 coordinates. */
+ xarr=x->array; yarr=y->array; zarr=z->array;
+ for(i=0;i<numslices;++i) { zarr[i]=i+1; xarr[i]=yarr[i]=1; }
+
+
+ /* Convert the coordinates to WCS. We are doing this inplace to avoid too
+ much memory/speed consumption. */
+ coords=x;
+ coords->next=y;
+ coords->next->next=z;
+ gal_wcs_img_to_world(coords, p->objects->wcs, 1);
+
+ /* For a check.
+ for(i=0;i<numslices;++i)
+ printf("%g, %g, %g\n", xarr[i], yarr[i], zarr[i]);
+ exit(0);
+ */
+
+ /* Allocate the slice counter array (we are doing it again because we
+ want it to be in integer type now). */
+ p->specsliceinfo=gal_data_alloc(NULL, slicenumtype, 1, &numslices, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, "SLICE",
+ "counter",
+ "Slice number in cube (counting from 1).");
+ if(p->specsliceinfo->type==GAL_TYPE_UINT16)
+ for(i=0;i<numslices;++i) ((uint16_t *)(p->specsliceinfo->array))[i]=i+1;
+ else
+ for(i=0;i<numslices;++i) ((uint32_t *)(p->specsliceinfo->array))[i]=i+1;
+
+ /* Set the slice WCS column information. Note that `z' is now the WCS
+ coordinate value of the third dimension, and to avoid wasting extra
+ space (this column is repeated one very object's spectrum), we'll
+ convert it to a 32-bit floating point dataset. */
+ p->specsliceinfo->next=gal_data_copy_to_new_type(z, GAL_TYPE_FLOAT32);
+
+ /* For a final check.
+ gal_table_write(p->specsliceinfo, NULL, GAL_TABLE_FORMAT_BFITS,
+ "specsliceinfo.fits", "test-debug",0);
+ */
+
+ /* Clean up. */
+ gal_data_free(x);
+ gal_data_free(y);
+ gal_data_free(z);
+}
+
+
+
+
+
/* Sanity checks and preparations for the upper-limit magnitude. */
static void
ui_preparations_upperlimit(struct mkcatalogparams *p)
@@ -1451,9 +1555,9 @@ void
ui_preparations(struct mkcatalogparams *p)
{
/* If no columns are requested, then inform the user. */
- if(p->columnids==NULL)
- error(EXIT_FAILURE, 0, "no columns requested, please run again with "
- "`--help' for the full list of columns you can ask for");
+ if(p->columnids==NULL && p->spectrum==0)
+ error(EXIT_FAILURE, 0, "no measurements requested! Please run again "
+ "with `--help' for the possible list of measurements");
/* Set the actual filenames to use. */
@@ -1485,6 +1589,14 @@ ui_preparations(struct mkcatalogparams *p)
ui_one_tile_per_object(p);
+ /* If a spectrum is requested, generate the two WCS columns. */
+ if(p->spectrum)
+ {
+ ui_preparations_spectrum_wcs(p);
+ p->spectra=gal_data_array_calloc(p->numobjects);
+ }
+
+
/* Allocate the reference random number generator and seed values. It
will be cloned once for every thread. If the user hasn't called
`envseed', then we want it to be different for every run, so we need
@@ -1667,7 +1779,7 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
mkcatalogparams *p)
void
ui_free_report(struct mkcatalogparams *p, struct timeval *t1)
{
- size_t d;
+ size_t d, i;
/* The temporary arrays for WCS coordinates. */
if(p->wcs_vo ) gal_list_data_free(p->wcs_vo);
@@ -1713,9 +1825,26 @@ ui_free_report(struct mkcatalogparams *p, struct timeval
*t1)
gal_data_free(p->upmask);
gal_data_free(p->clumps);
gal_data_free(p->objects);
+ gal_list_data_free(p->specsliceinfo);
if(p->upcheckout) free(p->upcheckout);
gal_data_array_free(p->tiles, p->numobjects, 0);
+ /* Clean up the spectra. */
+ if(p->spectra)
+ {
+ /* Note that each element of the array is the first node in a list of
+ datasets. So we can't free the first one with
+ `gal_list_data_free', we'll delete all the nodes after it in the
+ loop. */
+ for(i=0;i<p->numobjects;++i)
+ {
+ gal_list_data_free( p->spectra[i].next );
+ p->spectra[i].next=NULL;
+ gal_data_free_contents(&p->spectra[i]);
+ }
+ gal_data_array_free(p->spectra, p->numobjects, 0);
+ }
+
/* If the Sky or its STD image were given in tiles, then we defined a
tile structure to deal with them. The initialization of the tile
structure is checked with its `ndim' element. */
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 4560e65..f4d2f05 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -66,6 +66,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',
@@ -90,6 +91,7 @@ enum option_keys_enum
UI_KEY_SUBTRACTSKY,
UI_KEY_SFMAGNSIGMA,
UI_KEY_SFMAGAREA,
+ UI_KEY_SPECTRUM,
UI_KEY_UPMASKFILE,
UI_KEY_UPMASKHDU,
UI_KEY_UPNUM,
@@ -101,27 +103,38 @@ enum option_keys_enum
UI_KEY_OBJID, /* Catalog columns. */
UI_KEY_IDINHOSTOBJ,
+ UI_KEY_AREAXY,
UI_KEY_CLUMPSAREA,
UI_KEY_WEIGHTAREA,
UI_KEY_GEOAREA,
+ UI_KEY_GEOAREAXY,
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_MINX,
UI_KEY_MAXX,
UI_KEY_MINY,
UI_KEY_MAXY,
+ UI_KEY_MINZ,
+ UI_KEY_MAXZ,
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_BRIGHTNESSERR,
UI_KEY_CLUMPSBRIGHTNESS,
UI_KEY_BRIGHTNESSNORIVER,
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index c0576e7..7e3c5ad 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -295,9 +295,25 @@ upperlimit_write_comments(struct mkcatalogparams *p,
if(p->uprange)
{
- if( asprintf(&str, "Range of random samples about target: %zu, %zu",
- p->uprange[1], p->uprange[0])<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ switch(p->objects->ndim)
+ {
+ case 2:
+ if( asprintf(&str, "Range of random samples about target: "
+ "%zu, %zu", p->uprange[1], p->uprange[0])<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ break;
+ case 3:
+ if( asprintf(&str, "Range of random samples about target: %zu, "
+ "%zu, %zu", p->uprange[2], p->uprange[1],
+ p->uprange[0])<0 )
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "address the problem. The value %zu is not recognized for "
+ "`p->input->ndim'", __func__, PACKAGE_BUGREPORT,
+ p->objects->ndim);
+ }
gal_list_str_add(comments, str, 0);
}
@@ -347,19 +363,21 @@ upperlimit_write_comments(struct mkcatalogparams *p,
/* Write the values into a table for the user */
static void
upperlimit_write_check(struct mkcatalogparams *p, gal_list_sizet_t *check_x,
- gal_list_sizet_t *check_y, gal_list_f32_t *check_s)
+ gal_list_sizet_t *check_y, gal_list_sizet_t *check_z,
+ gal_list_f32_t *check_s)
{
float *sarr;
- gal_data_t *x, *y, *s;
char *tmp=NULL, *tmp2=NULL;
gal_list_str_t *comments=NULL;
- size_t *xarr, *yarr, tnum, num;
+ size_t *xarr, *yarr, *zarr=NULL, tnum, ttnum, num;
+ gal_data_t *x=NULL, *y=NULL, *z=NULL, *s=NULL; /* To avoid warnings. */
/* Convert the lists to an array. */
xarr=gal_list_sizet_to_array(check_x, 1, &num);
yarr=gal_list_sizet_to_array(check_y, 1, &tnum);
- if(tnum!=num)
+ if(check_z) zarr=gal_list_sizet_to_array(check_z, 1, &ttnum);
+ if(tnum!=num || (check_z && ttnum!=num) )
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
"problem. For some reason the size of the input lists don't "
"match (%zu, %zu)", __func__, PACKAGE_BUGREPORT, tnum, num);
@@ -377,16 +395,26 @@ upperlimit_write_check(struct mkcatalogparams *p,
gal_list_sizet_t *check_x,
y=gal_data_alloc(yarr, GAL_TYPE_SIZE_T, 1, &num, NULL, 0, p->cp.minmapsize,
p->cp.quietmmap, "RANDOM_Y", "pixel",
"Y-axis position of random footprint's first pixel.");
+ if(check_z)
+ z=gal_data_alloc(zarr, GAL_TYPE_SIZE_T, 1, &num, NULL, 0,
+ p->cp.minmapsize, p->cp.quietmmap, "RANDOM_Z", "pixel",
+ "Z-axis position of random footprint's first pixel.");
s=gal_data_alloc(sarr, GAL_TYPE_FLOAT32, 1, &num, NULL, 0, p->cp.minmapsize,
p->cp.quietmmap, "RANDOM_SUM",
p->values->unit ? p->values->unit : "input-units",
"Sum of pixel values over random footprint.");
- /* Convert the unsigned 64-bit values to 32-bit because the FITS table
- format doesn't recognize 64-bit integers.*/
- x=gal_data_copy_to_new_type_free(x, GAL_TYPE_UINT32);
- y=gal_data_copy_to_new_type_free(y, GAL_TYPE_UINT32);
+ /* If `size_t' isn't 32-bit on this system, then convert the unsigned
+ 64-bit values to 32-bit because the FITS table format doesn't
+ recognize 64-bit integers.*/
+ if( GAL_TYPE_SIZE_T != GAL_TYPE_UINT32 )
+ {
+ x=gal_data_copy_to_new_type_free( x, GAL_TYPE_UINT32);
+ y=gal_data_copy_to_new_type_free( y, GAL_TYPE_UINT32);
+ if(check_z)
+ z=gal_data_copy_to_new_type_free( z, GAL_TYPE_UINT32);
+ }
/* Write exactly what object/clump this table is for. */
@@ -415,7 +443,8 @@ upperlimit_write_check(struct mkcatalogparams *p,
gal_list_sizet_t *check_x,
/* Define a list from the containers and write them into a table. */
x->next=y;
- y->next=s;
+ if(check_z) { y->next=z; z->next=s; }
+ else { y->next=s; }
gal_list_str_reverse(&comments);
gal_table_write(x, comments, p->cp.tableformat, p->upcheckout,
"UPPERLIMIT_CHECK", 0);
@@ -428,6 +457,7 @@ upperlimit_write_check(struct mkcatalogparams *p,
gal_list_sizet_t *check_x,
gal_data_free(x);
gal_data_free(y);
gal_data_free(s);
+ if(check_z) gal_data_free(z);
}
@@ -549,10 +579,10 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp,
gal_data_t *tile,
struct gal_list_f32_t *check_s=NULL;
size_t d, counter=0, se_inc[2], nfailed=0;
float *V, *st_v, *uparr=pp->up_vals->array;
- size_t min[2], max[2], increment, num_increment;
- struct gal_list_sizet_t *check_x=NULL, *check_y=NULL;
+ size_t min[3], max[3], increment, num_increment;
int32_t *O, *OO, *oO, *st_o, *st_oo, *st_oc, *oC=NULL;
size_t maxfails = p->upnum * MKCATALOG_UPPERLIMIT_MAXFAILS_MULTIP;
+ struct gal_list_sizet_t *check_x=NULL, *check_y=NULL, *check_z=NULL;
size_t *rcoord=gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
"rcoord");
@@ -668,25 +698,42 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp,
gal_data_t *tile,
else ++nfailed;
- /* If a check is necessary, write in the values. */
+ /* If a check is necessary, write in the values (in FITS
+ coordinates). */
if(writecheck)
{
- gal_list_sizet_add(&check_x, rcoord[1]+1);
- gal_list_sizet_add(&check_y, rcoord[0]+1);
+ switch(ndim)
+ {
+ case 2:
+ gal_list_sizet_add(&check_x, rcoord[1]+1);
+ gal_list_sizet_add(&check_y, rcoord[0]+1);
+ break;
+
+ case 3:
+ gal_list_sizet_add(&check_x, rcoord[2]+1);
+ gal_list_sizet_add(&check_y, rcoord[1]+1);
+ gal_list_sizet_add(&check_z, rcoord[0]+1);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+ "to fix the problem. `ndim' value of %zu is not "
+ "recognized", __func__, PACKAGE_BUGREPORT, ndim);
+ }
gal_list_f32_add(&check_s, continueparse ? sum : NAN);
}
}
/* If a check is necessary, then write the values. */
if(writecheck)
- upperlimit_write_check(p, check_x, check_y, check_s);
+ upperlimit_write_check(p, check_x, check_y, check_z, check_s);
/* Do the measurement on the random distribution. */
upperlimit_measure(pp, clumplab, counter==p->upnum);
/* Reset the tile's array pointer, clean up and return. */
- tile->array=tarray;
free(rcoord);
+ tile->array=tarray;
gal_list_f32_free(check_s);
gal_list_sizet_free(check_x);
gal_list_sizet_free(check_y);
diff --git a/bin/mkprof/Makefile.am b/bin/mkprof/Makefile.am
index a8f23c5..884a270 100644
--- a/bin/mkprof/Makefile.am
+++ b/bin/mkprof/Makefile.am
@@ -44,4 +44,4 @@ EXTRA_DIST = main.h authors-cite.h args.h ui.h mkprof.h
oneprofile.h \
## The configuration file (distribute and install).
## NOTE: the man page is created in doc/Makefile.am
-dist_sysconf_DATA = astmkprof.conf
+dist_sysconf_DATA = astmkprof.conf astmkprof-3d.conf
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index b1852e5..be1039c 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -403,7 +403,7 @@ struct argp_option program_options[] =
UI_KEY_PCOL,
"STR/INT",
0,
- "Position angle.",
+ "Position angle (First X-Z-X Euler angle in 3D).",
UI_GROUP_CATALOG,
&p->pcol,
GAL_TYPE_STRING,
@@ -412,11 +412,37 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "p2col",
+ UI_KEY_P2COL,
+ "STR/INT",
+ 0,
+ "Second Euler angle (X-Z-X order).",
+ UI_GROUP_CATALOG,
+ &p->p2col,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "p3col",
+ UI_KEY_P2COL,
+ "STR/INT",
+ 0,
+ "Third Euler angle (X-Z-X order).",
+ UI_GROUP_CATALOG,
+ &p->p3col,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"qcol",
UI_KEY_QCOL,
"STR/INT",
0,
- "Axis ratio.",
+ "Axis ratio (major/dim2 in 3D).",
UI_GROUP_CATALOG,
&p->qcol,
GAL_TYPE_STRING,
@@ -425,6 +451,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "q2col",
+ UI_KEY_Q2COL,
+ "STR/INT",
+ 0,
+ "Axis ratio (major/dim3 in 3D).",
+ UI_GROUP_CATALOG,
+ &p->q2col,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"mcol",
UI_KEY_MCOL,
"STR/INT",
diff --git a/bin/mkprof/astmkprof-3d.conf b/bin/mkprof/astmkprof-3d.conf
new file mode 100644
index 0000000..9ccb1d0
--- /dev/null
+++ b/bin/mkprof/astmkprof-3d.conf
@@ -0,0 +1,57 @@
+# Default parameters (System) for MakeProfiles.
+# MakeProfiles is part of GNU Astronomy Utitlies.
+#
+# Use the long option name of each parameter followed by a value. The name
+# and value should be separated by atleast one white-space character (for
+# example ` '[space], or tab). Lines starting with `#' are ignored.
+#
+# For more information, please run these commands:
+#
+# $ astmkprof --help # Full list of options, short doc.
+# $ astmkprof -P # Print all options and used values.
+# $ info astmkprof # All options and input/output.
+# $ info gnuastro "Configuration files" # How to use configuration files.
+#
+# 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.
+
+#input
+ backhdu 1
+
+# Output:
+ mergedsize 100,100,200
+ oversample 3
+ circumwidth 2
+ type float32
+
+# Profiles:
+ tunitinp 0
+ numrandom 10000
+ tolerance 0.01
+ zeropoint 0.00
+
+# Catalog:
+ mode img
+ ccol 2
+ ccol 3
+ ccol 4
+ fcol 5
+ rcol 6
+ ncol 7
+ pcol 8
+ p2col 9
+ p3col 10
+ qcol 11
+ q2col 12
+ mcol 13
+ tcol 14
+
+# WCS:
+ crpix 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
+ ctype RA---TAN,DEC--TAN,AWAV
\ No newline at end of file
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index d8dfb49..832b1c6 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -137,12 +137,15 @@ struct mkprofparams
char *fcol; /* Column specifying profile function. */
char *rcol; /* Effective radius of profile. */
char *ncol; /* Sersic index column of profile. */
- char *pcol; /* Position angle column of profile. */
- char *qcol; /* Axis ratio column of profile. */
+ char *pcol; /* First Euler angle (X-Z-X order). */
+ char *p2col; /* Second Euler angle (X-Z-X order). */
+ char *p3col; /* Third Euler angle (X-Z-X order). */
+ char *qcol; /* Axis ratio1 (major/2nd dim. radius). */
+ char *q2col; /* Axis ratio2 (major/3rd dim. radius). */
char *mcol; /* Magnitude column. */
char *tcol; /* Truncation of the profiles. */
uint8_t mforflatpix; /* mcol is flat pixel value (f is 4 or 5). */
- uint8_t mcolisbrightness; /* mcol is total brightness not magnitude. */
+ uint8_t mcolisbrightness; /* mcol is total brightness, not magnitude. */
gal_data_t *crpix; /* CRPIX FITS header keywords. */
gal_data_t *crval; /* CRVAL FITS header keywords. */
gal_data_t *cdelt; /* For CDELTi FITS header keywords. */
@@ -161,11 +164,15 @@ struct mkprofparams
size_t num; /* The number of profiles. */
double *x; /* X axis position of profile center. */
double *y; /* Y axis position of profile center. */
+ double *z; /* Z axis position of profile center. */
uint8_t *f; /* Profile function code. */
float *r; /* Radius of profile. */
float *n; /* Index of profile. */
- float *p; /* Position angle of profile. */
- float *q; /* Axis ratio of profile. */
+ float *p1; /* First Euler angle (X-Z-X order). */
+ float *p2; /* Second Euler angle (X-Z-X order). */
+ float *p3; /* Third Euler angle (X-Z-X order). */
+ float *q1; /* Ratio of radius to second axis. */
+ float *q2; /* Ratio of radius to third axis. */
float *m; /* Magnitude of profile. */
float *t; /* Truncation distance. */
gsl_rng *rng; /* Main instance of random number generator.*/
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 7612f77..2d30e99 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -178,6 +178,10 @@ saveindividual(struct mkonthread *mkp)
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT64, "YCENTER", 0,
&p->y[id], 0, "Center of profile in catalog "
"(FITS axis 2)", 0, NULL);
+ if(ndim==3)
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT64, "ZCENTER", 0,
+ &p->z[id], 0, "Center of profile in catalog "
+ "(FITS axis 3)", 0, NULL);
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "RADIUS", 0,
&p->r[id], 0, "Radial parameter in catalog",
0, NULL);
@@ -185,12 +189,33 @@ saveindividual(struct mkonthread *mkp)
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PINDEX", 0,
&p->r[id], 0, "Index (Sersic or Moffat) of profile"
"in catalog", 0, NULL);
- gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA_DEG", 0,
- &p->p[id], 0, "Position angle of profile in catalog",
- 0, "deg");
- gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO", 0,
- &p->q[id], 0, "Axis ratio of profile in catalog",
- 0, NULL);
+ if(ndim==2)
+ {
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA_DEG", 0,
+ &p->p1[id], 0, "Position angle of profile in "
+ "catalog", 0, "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO", 0,
+ &p->q1[id], 0, "Axis ratio of profile in catalog",
+ 0, NULL);
+ }
+ else
+ {
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA1_DEG", 0,
+ &p->p1[id], 0, "First X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA2_DEG", 0,
+ &p->p2[id], 0, "Second X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA3_DEG", 0,
+ &p->p3[id], 0, "Third X-Z-X Euler angle in 3D", 0,
+ "deg");
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO1", 0,
+ &p->q1[id], 0, "Axis ratio along second dim",
+ 0, NULL);
+ gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO2", 0,
+ &p->q2[id], 0, "Axis ratio along third dim",
+ 0, NULL);
+ }
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "MAGNITUDE", 0,
&p->m[id], 0, "Magnitude of profile in catalog",
0, NULL);
@@ -288,7 +313,7 @@ mkprof_build_single(struct mkonthread *mkp, long *fpixel_i,
long *lpixel_i,
void *ptr;
int needs_crop=0;
size_t i, ind, fits_i, ndim=p->ndim;
- size_t start_indiv[2], start_mrg[2], dsize[2], os=p->oversample;
+ size_t start_indiv[3], start_mrg[3], dsize[3], os=p->oversample;
/* Make a copy of the main random number generator to use for this
profile (in this thread). */
@@ -455,11 +480,10 @@ mkprof_build(void *inparam)
{
struct mkonthread *mkp=(struct mkonthread *)inparam;
struct mkprofparams *p=mkp->p;
-
- double center[2];
size_t i, id, ndim=p->ndim;
struct builtqueue *ibq, *fbq=NULL;
- long fpixel_i[2], lpixel_i[2], fpixel_o[2], lpixel_o[2];
+ double center[3], semiaxes[3], euler_deg[3];
+ long fpixel_i[3], lpixel_i[3], fpixel_o[3], lpixel_o[3];
/* Make each profile that was specified for this thread. */
@@ -485,8 +509,28 @@ mkprof_build(void *inparam)
if( p->f[id] == PROFILE_POINT )
mkp->width[0]=mkp->width[1]=1;
else
- gal_box_bound_ellipse(mkp->truncr, mkp->q[0]*mkp->truncr,
- p->p[id], mkp->width);
+ switch(ndim)
+ {
+ case 2:
+ gal_box_bound_ellipse(mkp->truncr, mkp->q[0]*mkp->truncr,
+ p->p1[id], mkp->width);
+ break;
+
+ case 3:
+ euler_deg[0] = p->p1[id];
+ euler_deg[1] = p->p2[id];
+ euler_deg[2] = p->p3[id];
+ semiaxes[0] = mkp->truncr;
+ semiaxes[1] = mkp->truncr * mkp->q[0];
+ semiaxes[2] = mkp->truncr * mkp->q[1];
+ gal_box_bound_ellipsoid(semiaxes, euler_deg, mkp->width);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! please contact us at %s to "
+ "address the issue. %zu is not recognized for `ndim'",
+ __func__, PACKAGE_BUGREPORT, ndim);
+ }
/* Get the overlapping pixels using the starting points (NOT
@@ -495,6 +539,7 @@ mkprof_build(void *inparam)
{
center[0]=p->x[id];
center[1]=p->y[id];
+ if(ndim==3) center[2]=p->z[id];
gal_box_border_from_center(center, ndim, mkp->width, fpixel_i,
lpixel_i);
memcpy(mkp->fpixel_i, fpixel_i, ndim*sizeof *fpixel_i);
diff --git a/bin/mkprof/mkprof.h b/bin/mkprof/mkprof.h
index e10e5bd..8ddfe01 100644
--- a/bin/mkprof/mkprof.h
+++ b/bin/mkprof/mkprof.h
@@ -31,22 +31,22 @@ struct mkonthread
{
/* General parameters: */
double r; /* Elliptical radius at this point. */
- double coord[2]; /* Pixel coordinate. */
- double lower[2]; /* Coordinates of lower pixel position. */
- double higher[2]; /* Coordinates of higher pixel position. */
- double c[1]; /* Cosine of position angle(s). */
- double s[1]; /* Sine of position angle(s). */
- double q[1]; /* Axis ratio(s). */
- double center[2]; /* Center (in FITS) in oversampled image.*/
+ double coord[3]; /* Pixel coordinate. */
+ double lower[3]; /* Coordinates of lower pixel position. */
+ double higher[3]; /* Coordinates of higher pixel position. */
+ double c[3]; /* Cosine of position angle(s). */
+ double s[3]; /* Sine of position angle(s). */
+ double q[2]; /* Axis ratio(s). */
+ double center[3]; /* Center (in FITS) in oversampled image.*/
double (*profile)(struct mkonthread *); /* Function to use. */
double truncr; /* Truncation radius in pixels. */
double intruncr; /* Inner truncation radius in pixels. */
- long width[2]; /* Enclosing box in FITS axes, not C. */
+ long width[3]; /* Enclosing box in FITS axes, not C. */
float peakflux; /* Flux at profile peak. */
float brightness; /* The brightness of the profile. */
- uint8_t func; /* Radial function of the profile. */
+ uint8_t func; /* Radial function code of the profile. */
long *onaxes; /* Sides of the unover-sampled image. */
- long fpixel_i[2]; /* fpixel_i before running overlap. */
+ long fpixel_i[3]; /* fpixel_i before running overlap. */
int correction; /* ==1: correct the pixels afterwards. */
unsigned long rng_seed; /* Seed used to generate this profile. */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index d01ffa2..0aacb3c 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -69,7 +69,7 @@ oneprofile_center_oversampled(struct mkonthread *mkp)
for(i=0;i<p->ndim;++i)
{
- dim = i==0 ? p->x : p->y;
+ dim = i==0 ? p->x : (i==1 ? p->y : p->z);
pixfrac = modf(fabs(dim[id]), &intpart);
val = ( os*(mkp->width[i]/2 + pixfrac)
+ (pixfrac<0.5f ? os/2 : -1*os/2-1) );
@@ -84,7 +84,7 @@ oneprofile_center_oversampled(struct mkonthread *mkp)
static void
oneprofile_set_coord(struct mkonthread *mkp, size_t index)
{
- size_t i, coord_c[2];
+ size_t i, coord_c[3];
uint8_t os=mkp->p->oversample;
size_t ndim=mkp->ibq->image->ndim, *dsize=mkp->ibq->image->dsize;
@@ -103,18 +103,40 @@ oneprofile_set_coord(struct mkonthread *mkp, size_t index)
/* Convert cartesian coordinates to the rotated elliptical radius. See the
- "Defining an ellipse" section of the book for the full derivation. */
+ "Defining an ellipse and ellipsoid" section of the book for the full
+ derivation. */
static void
oneprofile_r_el(struct mkonthread *mkp)
{
- double Xr, Yr; /* Rotated x and y. */
- double q=mkp->q[0];
- double c=mkp->c[0], s=mkp->s[0];
- double x=mkp->coord[0], y=mkp->coord[1];
-
- Xr = x * ( c ) + y * ( s );
- Yr = x * ( -1.0f*s ) + y * ( c );
- mkp->r = sqrt( Xr*Xr + Yr*Yr/q/q );
+ double Xr, Yr, Zr; /* Rotated x, y, z. */
+ double q1=mkp->q[0], q2=mkp->q[1];
+ double c1=mkp->c[0], s1=mkp->s[0];
+ double c2=mkp->c[1], s2=mkp->s[1];
+ double c3=mkp->c[2], s3=mkp->s[2];
+ double x=mkp->coord[0], y=mkp->coord[1], z=mkp->coord[2];
+
+ switch(mkp->p->ndim)
+ {
+ case 2:
+ /* The parenthesis aren't necessary, but help in readability and
+ avoiding human induced bugs. */
+ Xr = x * ( c1 ) + y * ( s1 );
+ Yr = x * ( -1.0f*s1 ) + y * ( c1 );
+ mkp->r = sqrt( Xr*Xr + Yr*Yr/q1/q1 );
+ break;
+
+ case 3:
+ Xr = x*( c3*c1 - s3*c2*s1 ) + y*( c3*s1 + s3*c2*c1) + z*( s3*s2 );
+ Yr = x*( -1*s3*c1 - c3*c2*s1 ) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
+ Zr = x*( s1*s2 ) + y*(-1*s2*c1 ) + z*( c2 );
+ mkp->r = sqrt( Xr*Xr + Yr*Yr/q1/q1 + Zr*Zr/q2/q2 );
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The value %zu is not recognized for "
+ "`mkp->p->ndim'", __func__, PACKAGE_BUGREPORT, mkp->p->ndim);
+ }
}
@@ -128,7 +150,8 @@ oneprofile_r_el(struct mkonthread *mkp)
static float
oneprofile_r_circle(size_t index, struct mkonthread *mkp)
{
- size_t i, c[2];
+
+ size_t i, c[3];
double d, sum=0.0f;
size_t ndim=mkp->ibq->image->ndim, *dsize=mkp->ibq->image->dsize;
@@ -173,9 +196,9 @@ float
oneprofile_randompoints(struct mkonthread *mkp)
{
double r_before=mkp->r;
- double range[2], sum=0.0f;
- double coord_before[2]={mkp->coord[0], mkp->coord[1]};
+ double range[3], sum=0.0f;
size_t i, j, numrandom=mkp->p->numrandom, ndim=mkp->p->ndim;
+ double coord_before[3]={mkp->coord[0], mkp->coord[1], mkp->coord[2]};
/* Set the range in each dimension. */
for(i=0;i<ndim;++i)
@@ -199,6 +222,7 @@ oneprofile_randompoints(struct mkonthread *mkp)
mkp->r=r_before;
mkp->coord[0]=coord_before[0];
mkp->coord[1]=coord_before[1];
+ mkp->coord[2]=coord_before[2];
return sum/numrandom;
}
@@ -309,7 +333,7 @@ oneprofile_center_pix_index(struct mkonthread *mkp)
{
double pixfrac, intpart;
size_t *dsize=mkp->ibq->image->dsize;
- size_t i, coord[2], ndim=mkp->p->ndim;
+ size_t i, coord[3], ndim=mkp->p->ndim;
/* Find the coordinates of the center point. Note `mkp->center' is in
FITS coordinates, while coord must be in C coordinates (to be used in
@@ -520,7 +544,7 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
double sigma;
int tp=p->tunitinp;
- size_t id=mkp->ibq->id;
+ size_t id=mkp->ibq->id, ndim=p->ndim;
/* Fill the most basic profile agnostic parameters. */
mkp->brightness = ( p->mcolisbrightness
@@ -530,14 +554,40 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
mkp->func = mkp->ibq->func = p->f[id];
- /* Shifts were already multiplied with oversample. Just note that
- p->x and p->y are in the FITS ordering, while p->shift is in C
- ordering. */
- mkp->q[0] = p->q[id];
- p->x[id] += p->shift[1]/p->oversample;
- p->y[id] += p->shift[0]/p->oversample;
- mkp->c[0] = cos( p->p[id] * DEGREESTORADIANS );
- mkp->s[0] = sin( p->p[id] * DEGREESTORADIANS );
+ /* Fill in the dimension-dependent parameters. */
+ switch(ndim)
+ {
+ case 2:
+ /* Shifts were already multiplied with oversample. Just note that
+ p->x and p->y are in the FITS ordering, while p->shift is in C
+ ordering. */
+ mkp->q[0] = p->q1[id];
+ p->x[id] += p->shift[1]/p->oversample;
+ p->y[id] += p->shift[0]/p->oversample;
+ mkp->c[0] = cos( p->p1[id] * DEGREESTORADIANS );
+ mkp->s[0] = sin( p->p1[id] * DEGREESTORADIANS );
+ break;
+
+ case 3:
+ /* See comments for 2D. */
+ mkp->q[0] = p->q1[id];
+ mkp->q[1] = p->q2[id];
+ p->x[id] += p->shift[2]/p->oversample;
+ p->y[id] += p->shift[1]/p->oversample;
+ p->z[id] += p->shift[0]/p->oversample;
+ mkp->c[0] = cos( p->p1[id] * DEGREESTORADIANS );
+ mkp->s[0] = sin( p->p1[id] * DEGREESTORADIANS );
+ mkp->c[1] = cos( p->p2[id] * DEGREESTORADIANS );
+ mkp->s[1] = sin( p->p2[id] * DEGREESTORADIANS );
+ mkp->c[2] = cos( p->p3[id] * DEGREESTORADIANS );
+ mkp->s[2] = sin( p->p3[id] * DEGREESTORADIANS );
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "address the problem. The value `%zu' is not recognized for "
+ "`ndim'", __func__, PACKAGE_BUGREPORT, ndim);
+ }
/* Fill the profile-dependent parameters. */
@@ -566,6 +616,8 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
mkp->brightness = 1.0f; /* When the PSF is a separate image, */
p->x[id] = 0.0f; /* it should be centered and have a */
p->y[id] = 0.0f; /* total brightness of 1.0f. */
+ if(ndim==3)
+ p->z[id] = 0.0f;
}
break;
@@ -582,6 +634,8 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
mkp->brightness = 1.0f; /* Same as the explanations for */
p->x[id] = 0.0f; /* The Moffat profile. */
p->y[id] = 0.0f;
+ if(ndim==3)
+ p->z[id] = 0.0f;
}
break;
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 47dbbad..ce000ab 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -191,7 +191,6 @@ ui_initialize_options(struct mkprofparams *p,
p->zeropoint = NAN;
p->cp.type = GAL_TYPE_FLOAT32;
-
/* Modify the common options for this program. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
{
@@ -284,7 +283,7 @@ ui_parse_kernel(struct argp_option *option, char *arg,
double *darray;
gal_data_t *kernel;
size_t i, nc, need=0;
- char *c, *profile, *tailptr;
+ char *c, *dstr, *profile, *tailptr;
char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
/* We want to print the stored values. */
@@ -297,7 +296,15 @@ ui_parse_kernel(struct argp_option *option, char *arg,
/* First write the profile function code into the output string. */
nc=0;
profile=ui_profile_name_write(kernel->status);
- nc += sprintf(sstr+nc, "%s,", profile);
+ switch(kernel->flag)
+ {
+ case 2: nc += sprintf(sstr+nc, "%s,", profile); break;
+ case 3: nc += sprintf(sstr+nc, "%s-3d,", profile); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. %u is not a recognized kernel "
+ "dimensionality", __func__, PACKAGE_BUGREPORT, kernel->flag);
+ }
/* Write the values into a string. */
for(i=0;i<kernel->size;++i)
@@ -343,6 +350,32 @@ ui_parse_kernel(struct argp_option *option, char *arg,
"positive", i+1, darray[i], arg);
+ /* See if a 2D kernel is requested or a 3D kernel and keep the value
+ in `kernel->flag'. If no dimensionality is defined, then by
+ default, we'll assume it is 2D.*/
+ c=profile;while(*c!='\0' && *c!='-') ++c;
+ if(*c=='\0')
+ kernel->flag=2;
+ else
+ {
+ *c='\0';
+ dstr=c+1;
+ if( (dstr[1]!='d' && dstr[1]!='D') || dstr[2]!='\0')
+ error(EXIT_FAILURE, 0, "bad formatting in `--kernel' "
+ "dimensionality. The dimensionality suffix must be either "
+ "2d, 3d (not case sensitive). You have given `%s'", dstr);
+ switch(dstr[0])
+ {
+ case '2': kernel->flag=2; break;
+ case '3': kernel->flag=3; break;
+ default:
+ error(EXIT_FAILURE, 0, "only 2 or 3 dimensional kernels can "
+ "currently be built, you have asked for a %c "
+ "dimensional kernel", dstr[0]);
+ }
+ }
+
+
/* Write the profile type code into `kernel->status'. If it starts
with a digit, then the user might have given the code of the
profile directly. In that case, parse the number. Otherwise,
@@ -367,13 +400,14 @@ ui_parse_kernel(struct argp_option *option, char *arg,
/* Make sure the number of parameters conforms with the profile. */
switch(kernel->status)
{
- case PROFILE_SERSIC: need = 3; break;
- case PROFILE_MOFFAT: need = 3; break;
- case PROFILE_GAUSSIAN: need = 2; break;
- case PROFILE_POINT: need = 0; break;
- case PROFILE_FLAT: need = 1; break;
- case PROFILE_CIRCUMFERENCE: need = 1; break;
- case PROFILE_DISTANCE: need = 1; break;
+
+ case PROFILE_SERSIC: need = kernel->flag==2 ? 3 : 4; break;
+ case PROFILE_MOFFAT: need = kernel->flag==2 ? 3 : 4; break;
+ case PROFILE_GAUSSIAN: need = kernel->flag==2 ? 2 : 3; break;
+ case PROFILE_POINT: need = 0; break;
+ case PROFILE_FLAT: need = kernel->flag==2 ? 1 : 2; break;
+ case PROFILE_CIRCUMFERENCE: need = kernel->flag==2 ? 1 : 2; break;
+ case PROFILE_DISTANCE: need = kernel->flag==2 ? 1 : 2; break;
default:
error_at_line(EXIT_FAILURE, 0, filename, lineno, "%s: a bug! "
"Please contact us at %s to correct the issue. "
@@ -385,9 +419,9 @@ ui_parse_kernel(struct argp_option *option, char *arg,
/* Make sure the number of parameters given are the same number that
are needed. */
if( kernel->size != need )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "as a kernel, "
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "as a %uD kernel, "
"a `%s' profile needs %zu parameters, but %zu "
- "parameter%s given to `--kernel'",
+ "parameter%s given to `--kernel'", kernel->flag,
ui_profile_name_write(kernel->status), need,
kernel->size, kernel->size>1?"s are":" is");
@@ -487,7 +521,6 @@ ui_read_check_only_options(struct mkprofparams *p)
"(i.e., when the contents of `--mcol' must be interpretted as "
"a magnitude, not brightness).");
-
/* Make sure no zero value is given for the `--mergedsize' option (only
when it is necessary). */
if(p->dsize && p->backname==NULL)
@@ -591,7 +624,7 @@ ui_check_options_and_arguments(struct mkprofparams *p)
/*************** Preparations *******************/
/**************************************************************/
static void
-ui_read_cols(struct mkprofparams *p)
+ui_read_cols_2d(struct mkprofparams *p)
{
int checkblank;
size_t i, counter=0;
@@ -724,21 +757,21 @@ ui_read_cols(struct mkprofparams *p)
case 6:
colname="position angle (`pcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
- p->p=corrtype->array;
+ p->p1=corrtype->array;
break;
case 7:
colname="axis ratio (`qcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
- p->q=corrtype->array;
+ p->q1=corrtype->array;
/* Check if there is no negative or >1.0f axis ratio. */
for(i=0;i<p->num;++i)
- if( p->f[i]!=PROFILE_POINT && (p->q[i]<=0.0f || p->q[i]>1.0f) )
+ if( p->f[i]!=PROFILE_POINT && (p->q1[i]<=0.0f || p->q1[i]>1.0f) )
error(EXIT_FAILURE, 0, "%s: row %zu, the axis ratio value %g "
"is not acceptable for a `%s' profile. It has to be >0 "
- "and <=1", p->catname, i+1, p->q[i],
+ "and <=1", p->catname, i+1, p->q1[i],
ui_profile_name_write(p->f[i]));
break;
@@ -819,13 +852,249 @@ ui_read_cols(struct mkprofparams *p)
+/* Read the columns for a 3D profile. */
+static void
+ui_read_cols_3d(struct mkprofparams *p)
+{
+ int checkblank;
+ size_t i, counter=0;
+ char *colname=NULL, **strarr;
+ gal_data_t *cols, *tmp, *corrtype=NULL;
+ gal_list_str_t *lines, *ccol, *colstrs=NULL;
+
+ /* The 3D-specific columns are not mandatory in `args.h', so we need to
+ check here if they are given or not before starting to read them. */
+ if(p->p2col==NULL || p->p3col==NULL || p->q2col==NULL)
+ error(EXIT_FAILURE, 0, "at least one of `--p2col', `--p3col', or "
+ "`--q2col' have not been identified. When building a 3D profile, "
+ "these three columns are also mandatory");
+
+ /* The coordinate columns are a linked list of strings. */
+ ccol=p->ccol;
+ for(i=0; i<p->ndim; ++i)
+ {
+ gal_list_str_add(&colstrs, ccol->v, 0);
+ ccol=ccol->next;
+ }
+
+ /* Put the columns a specific order to read. */
+ gal_list_str_add(&colstrs, p->fcol, 0);
+ gal_list_str_add(&colstrs, p->rcol, 0);
+ gal_list_str_add(&colstrs, p->ncol, 0);
+ gal_list_str_add(&colstrs, p->pcol, 0);
+ gal_list_str_add(&colstrs, p->p2col, 0);
+ gal_list_str_add(&colstrs, p->p3col, 0);
+ gal_list_str_add(&colstrs, p->qcol, 0);
+ gal_list_str_add(&colstrs, p->q2col, 0);
+ gal_list_str_add(&colstrs, p->mcol, 0);
+ gal_list_str_add(&colstrs, p->tcol, 0);
+
+ /* Reverse the order to make the column orders correspond to how we added
+ them here and avoid possible bugs. */
+ gal_list_str_reverse(&colstrs);
+
+ /* Read the desired columns from the file. */
+ lines=gal_options_check_stdin(p->catname, p->cp.stdintimeout, "input");
+ cols=gal_table_read(p->catname, p->cp.hdu, lines, colstrs, p->cp.searchin,
+ p->cp.ignorecase, p->cp.minmapsize, p->cp.quietmmap,
+ NULL);
+ gal_list_str_free(lines, 1);
+
+ /* Set the number of objects. */
+ p->num=cols->size;
+
+ /* Put each column's data in the respective internal array. */
+ while(cols!=NULL)
+ {
+ /* Pop out the top column. */
+ tmp=gal_list_data_pop(&cols);
+
+ /* By default check if the column has blank values, but it can be
+ turned off for some columns. */
+ checkblank=1;
+
+ /* See which column we are currently reading. */
+ switch(++counter)
+ {
+ case 1:
+ case 2:
+ case 3:
+ colname = ( counter==1
+ ? "first coordinate column (`--coordcol')"
+ : ( counter==2
+ ? "second coordinate column (`--coordcol')"
+ : "third coordinate column (`--coordcol')" ) );
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT64);
+ switch(counter)
+ {
+ case 1: p->x=corrtype->array; break;
+ case 2: p->y=corrtype->array; break;
+ case 3: p->z=corrtype->array; break;
+ }
+ break;
+
+ case 4:
+ if(tmp->type==GAL_TYPE_STRING)
+ {
+ p->f=gal_pointer_allocate(GAL_TYPE_UINT8, p->num, 0,
+ __func__, "p->f");
+ strarr=tmp->array;
+ for(i=0;i<p->num;++i)
+ p->f[i]=ui_profile_name_read(strarr[i], i+1);
+ gal_data_free(tmp);
+ corrtype=NULL;
+ }
+ else
+ {
+ /* Read the user's profile codes. */
+ colname="profile function code (`fcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_UINT8);
+ p->f=corrtype->array;
+
+ /* Check if they are in the correct range. */
+ for(i=0;i<p->num;++i)
+ if(p->f[i]<=PROFILE_INVALID || p->f[i]>=PROFILE_MAXIMUM_CODE)
+ error(EXIT_FAILURE, 0, "%s: row %zu, the function "
+ "code is %u. It should be >%d and <%d. Please run "
+ "again with `--help' and check the acceptable "
+ "codes.\n\nAlternatively, you can use alphabetic "
+ "strings to specify the profile functions, see the "
+ "explanations under `fcol' from the command "
+ "below (press the `SPACE' key to go down, and the "
+ "`q' to return back to the command-line):\n\n"
+ " $ info %s\n", p->catname, i+1, p->f[i],
+ PROFILE_INVALID, PROFILE_MAXIMUM_CODE, PROGRAM_EXEC);
+ }
+ break;
+
+ case 5:
+ colname="radius (`rcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->r=corrtype->array;
+
+ /* Check if there is no negative or zero-radius profile. */
+ for(i=0;i<p->num;++i)
+ if(p->f[i]!=PROFILE_POINT && p->r[i]<=0.0f)
+ error(EXIT_FAILURE, 0, "%s: row %zu, the radius value %g is "
+ "not acceptable for a `%s' profile. It has to be larger "
+ "than 0", p->catname, i+1, p->r[i],
+ ui_profile_name_write(p->f[i]));
+ break;
+
+ case 6:
+ colname="index (`ncol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->n=corrtype->array;
+ break;
+
+ case 7:
+ colname="first euler angle (`pcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->p1=corrtype->array;
+ break;
+
+ case 8:
+ colname="second euler angle (`p2col')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->p2=corrtype->array;
+ break;
+
+ case 9:
+ colname="third euler angle (`p3col')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->p3=corrtype->array;
+ break;
+
+ case 10:
+ colname="axis ratio 1 (`qcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->q1=corrtype->array;
+
+ /* Check if there is no negative or >1.0f axis ratio. */
+ for(i=0;i<p->num;++i)
+ if( p->f[i]!=PROFILE_POINT && (p->q1[i]<=0.0f || p->q1[i]>1.0f) )
+ error(EXIT_FAILURE, 0, "%s: row %zu, the first axis ratio "
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be >0 and <=1", p->catname, i+1, p->q1[i],
+ ui_profile_name_write(p->f[i]));
+ break;
+
+ case 11:
+ colname="axis ratio 2 (`q2col')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->q2=corrtype->array;
+
+ /* Check if there is no negative or >1.0f axis ratio. */
+ for(i=0;i<p->num;++i)
+ if( p->f[i]!=PROFILE_POINT && (p->q2[i]<=0.0f || p->q2[i]>1.0f) )
+ error(EXIT_FAILURE, 0, "%s: row %zu, the second axis ratio "
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be >0 and <=1", p->catname, i+1, p->q2[i],
+ ui_profile_name_write(p->f[i]));
+ break;
+
+ case 12:
+ colname="magnitude (`mcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->m=corrtype->array;
+ checkblank=0; /* Magnitude can be NaN: to mask regions. */
+ break;
+
+ case 13:
+ colname="truncation (`tcol')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->t=corrtype->array;
+
+ /* Check if there is no negative or zero truncation radius. */
+ for(i=0;i<p->num;++i)
+ if(p->f[i]!=PROFILE_POINT && p->t[i]<=0.0f)
+ error(EXIT_FAILURE, 0, "%s: row %zu, the truncation radius "
+ "value %g is not acceptable for a `%s' profile. It has "
+ "to be larger than 0", p->catname, i+1, p->t[i],
+ ui_profile_name_write(p->f[i]));
+ break;
+
+ /* If the index isn't recognized, then it is larger, showing that
+ there was more than one match for the given criteria */
+ default:
+ gal_tableintern_error_col_selection(p->catname, p->cp.hdu, "too "
+ "many columns were selected "
+ "by the given values to the "
+ "options ending in `col'.");
+ }
+
+ /* Sanity check and clean up. Note that it might happen that the
+ input structure is already freed. In that case, `corrtype' will be
+ NULL. */
+ if(corrtype)
+ {
+ /* Make sure there are no blank values in this column. */
+ if( checkblank && gal_blank_present(corrtype, 1) )
+ error(EXIT_FAILURE, 0, "%s column has blank values. "
+ "Input columns cannot contain blank values", colname);
+
+ /* Free the unnecessary sturcture information. The correct-type
+ (`corrtype') data structure's array is necessary for later
+ steps, so its pointer has been copied in the main program's
+ structure. Hence, we should set the structure's pointer to
+ NULL so the important data isn't freed.*/
+ corrtype->array=NULL;
+ gal_data_free(corrtype);
+ }
+ }
+}
+
+
+
+
+
/* It is possible to define the internal catalog through a catalog or the
`--kernel' option. This function will do the job. */
static void
ui_prepare_columns(struct mkprofparams *p)
{
double *karr;
- float r, n, t;
+ float r, n, t, q2;
/* If the kernel option was called, then we need to build a series of
single element columns to create an internal catalog. */
@@ -835,15 +1104,22 @@ ui_prepare_columns(struct mkprofparams *p)
p->num=1;
/* Allocate the necessary columns. */
- p->x = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->x");
- p->y = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->y");
- p->f = gal_pointer_allocate(GAL_TYPE_UINT8, 1, 1, __func__, "p->f");
- p->r = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->r");
- p->n = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->n");
- p->p = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->p");
- p->q = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->q");
- p->m = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->m");
- p->t = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->t");
+ p->x = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->x");
+ p->y = gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->y");
+ p->f = gal_pointer_allocate(GAL_TYPE_UINT8, 1, 1, __func__, "p->f");
+ p->r = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->r");
+ p->n = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->n");
+ p->p1 = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->p1");
+ p->q1 = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->q1");
+ p->m = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->m");
+ p->t = gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__, "p->t");
+ if(p->ndim==3)
+ {
+ p->z =gal_pointer_allocate(GAL_TYPE_FLOAT64, 1, 1, __func__, "p->z");
+ p->p2=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__,
"p->p2");
+ p->p3=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__,
"p->p3");
+ p->q2=gal_pointer_allocate(GAL_TYPE_FLOAT32, 1, 1, __func__,
"p->q2");
+ }
/* For profiles that need a different number of input values. Note
that when a profile doesn't need a value, it will be ignored. */
@@ -852,20 +1128,46 @@ ui_prepare_columns(struct mkprofparams *p)
{
r = karr[0];
n = p->kernel->size==2 ? 0.0f : karr[1];
- t = p->kernel->size==1 ? 1.0f : karr[ p->kernel->size - 1 ];
+ t = ( p->ndim==2
+ ? p->kernel->size==1 ? 1.0f : karr[ p->kernel->size - 1 ]
+ : p->kernel->size==1 ? 1.0f : karr[ p->kernel->size - 2 ] );
}
else r=n=t=0.0f;
/* Fill the allocated spaces. */
- p->x[0] = 0.0f;
- p->y[0] = 0.0f;
- p->f[0] = p->kernel->status;
- p->r[0] = r;
- p->n[0] = n;
- p->p[0] = 0.0f;
- p->q[0] = 1.0f;
- p->m[0] = 0.0f;
- p->t[0] = t;
+ p->x[0] = 0.0f;
+ p->y[0] = 0.0f;
+ p->f[0] = p->kernel->status;
+ p->r[0] = r;
+ p->n[0] = n;
+ p->p1[0] = 0.0f;
+ p->q1[0] = 1.0f;
+ p->m[0] = 0.0f;
+ p->t[0] = t;
+ if(p->ndim==3)
+ {
+ /* Parameters for any case. */
+ p->z[0] = 0.0f;
+ q2 = p->kernel->size ? karr[ p->kernel->size - 1 ] : 0.0f;
+
+ /* 3rd-dim axis ratio > 1: Set the major axis in the direction of
+ the 3rd dimension (90 degree rotation for all three
+ rotations). Also set the two axis ratios to the inverse of the
+ requested value. */
+ if(q2>1.0)
+ {
+ p->q1[0] = p->q2[0] = 1/q2;
+ p->p1[0] = p->p2[0] = p->p3[0] = 90.0;
+ }
+
+ /* 3rd-dim axis ratio <=1: No extra rotation is necessary and
+ `q2'can simply be put in the respective column. */
+ else
+ {
+ p->q2[0] = q2;
+ p->p2[0] = p->p3[0] = 0.0;
+ }
+ }
}
else
{
@@ -878,8 +1180,16 @@ ui_prepare_columns(struct mkprofparams *p)
"`--coordcol') given but output has %zu dimensions",
gal_list_str_number(p->ccol), p->ndim);
- /* Call the column-reading function. */
- ui_read_cols(p);
+ /* Call the respective function. */
+ switch(p->ndim)
+ {
+ case 2: ui_read_cols_2d(p); break;
+ case 3: ui_read_cols_3d(p); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "resolve the issue. %zu not recognized for `p->ndim'",
+ __func__, PACKAGE_BUGREPORT, p->ndim);
+ }
}
}
@@ -1024,10 +1334,10 @@ static void
ui_prepare_canvas(struct mkprofparams *p)
{
float *f, *ff;
- double truncr;
- long width[2]={1,1};
+ long width[3]={1,1,1};
size_t tndim, *tdsize;
int status=0, setshift=0;
+ double truncr, semiaxes[3], euler_deg[3];
size_t i, nshift=0, *dsize=NULL, ndim_counter;
/* If a background image is specified, then use that as the output
@@ -1103,8 +1413,19 @@ ui_prepare_canvas(struct mkprofparams *p)
half of it for the shift. */
setshift=1;
truncr = p->tunitinp ? p->t[i] : p->t[i] * p->r[i]/2;
- gal_box_bound_ellipse(truncr, p->q[i]*truncr, p->p[i],
- width);
+ if(p->ndim==2)
+ gal_box_bound_ellipse(truncr, p->q1[i]*truncr,
+ p->p1[i], width);
+ else
+ {
+ euler_deg[0] = p->p1[i];
+ euler_deg[1] = p->p2[i];
+ euler_deg[2] = p->p3[i];
+ semiaxes[0] = truncr;
+ semiaxes[1] = truncr * p->q1[i];
+ semiaxes[2] = truncr * p->q2[i];
+ gal_box_bound_ellipsoid( semiaxes, euler_deg, width);
+ }
}
/* Either set the shifts to zero or to the values set from
@@ -1118,6 +1439,7 @@ ui_prepare_canvas(struct mkprofparams *p)
{
p->shift[0] = (width[0]/2)*p->oversample;
p->shift[1] = (width[1]/2)*p->oversample;
+ if(p->ndim==3) p->shift[2] = (width[2]/2)*p->oversample;
}
}
}
@@ -1209,12 +1531,13 @@ ui_finalize_coordinates(struct mkprofparams *p)
/* 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;
+ 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 currently supported for "
- "%zu-dimensional datasets", ndim);
+ "coordinates is not supported for %zu-dimensional "
+ "datasets", ndim);
}
/* Allocate the list of coordinates. */
@@ -1322,7 +1645,7 @@ ui_read_ndim(struct mkprofparams *p)
if(p->kernel)
{
/* The kernel's dimensionality is fixed. */
- p->ndim=2;
+ p->ndim=p->kernel->flag;
/* Make sure the kernel and background are not given together. */
if(p->backname)
@@ -1362,9 +1685,9 @@ ui_read_ndim(struct mkprofparams *p)
}
/* Make sure the dimensionality is supported. */
- if(p->ndim!=2)
+ if(p->ndim!=2 && p->ndim!=3)
error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions. Currently "
- "only 2 dimensional outputs can be produced",
+ "only 2 or 3 dimensional outputs can be produced",
p->backname, p->backhdu, p->ndim);
}
else
@@ -1375,9 +1698,9 @@ ui_read_ndim(struct mkprofparams *p)
p->ndim=ndim_counter;
/* Make sure the dimensionality is supported. */
- if(p->ndim!=2)
+ if(p->ndim!=2 && p->ndim!=3)
error(EXIT_FAILURE, 0, "%zu values given to `--mergedsize'. "
- "Currently only 2 dimensional outputs can be produced",
+ "Currently only 2 or 3 dimensional outputs can be produced",
p->ndim);
}
}
@@ -1403,10 +1726,10 @@ ui_preparations(struct mkprofparams *p)
if(p->kernel)
{
/* Set the necessary constants. */
- p->ndim=2;
p->nomerged=1;
p->psfinimg=0;
p->individual=1;
+ p->ndim=p->kernel->flag;
/* Set the shift array. */
p->shift=gal_pointer_allocate(GAL_TYPE_SIZE_T, p->ndim, 1,
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index 429fbf3..4b0898c 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -81,7 +81,10 @@ enum option_keys_enum
UI_KEY_RCOL,
UI_KEY_NCOL,
UI_KEY_PCOL,
+ UI_KEY_P2COL,
+ UI_KEY_P3COL,
UI_KEY_QCOL,
+ UI_KEY_Q2COL,
UI_KEY_MCOL,
UI_KEY_TCOL,
UI_KEY_CRPIX,
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index da77151..2a95287 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -766,7 +766,7 @@ ui_matrix_finalize(struct warpparams *p)
error(EXIT_FAILURE, 0, "the determinant of the given matrix "
"is zero");
- /* Note yet implemented: Check if the transformation is spatially
+ /* Not yet implemented: Check if the transformation is spatially
invariant, in other words, if it differs between differet regions of
the output. If it doesn't we can use this information for a more
efficient processing. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 78c8b22..852cac6 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -535,7 +535,7 @@ MakeProfiles
Modeling basics
-* Defining an ellipse:: An ellipse on a pixelated grid.
+* Defining an ellipse and ellipsoid:: Definition of these important shapes.
* PSF:: Radial profiles for the PSF.
* Stars:: Making mock star profiles.
* Galaxies:: Radial profiles for galaxies.
@@ -9577,8 +9577,6 @@ When the crop was not defined by its center (see
@ref{Crop modes}), or @option{-
-
-
@node Arithmetic, Convolve, Crop, Data manipulation
@section Arithmetic
@@ -9876,11 +9874,11 @@ So depending on the type of the input and its nature,
it is recommended to use o
If any WCS is present, the returned dataset will also lack the respective
dimension in its WCS matrix.
Therefore, when the WCS is important for later processing, be sure that the
input is aligned with the respective axes: all non-diagonal elements in the WCS
matrix are zero.
-@cindex IFU
@cindex Data cubes
@cindex 3D data-cubes
@cindex Cubes (3D data)
@cindex Narrow-band image
+@cindex IFU: Integral Field Unit
@cindex Integral field unit (IFU)
One common application of this operator is the creation of pseudo broad-band
or narrow-band 2D images from 3D data cubes.
For example integral field unit (IFU) data products that have two spatial
dimensions (first two FITS dimensions) and one spectral dimension (third FITS
dimension).
@@ -10878,7 +10876,8 @@
C(\omega)=H(\omega)\int_{-\infty}^{\infty}f(\tau)e^{-i{\omega}\tau}d\tau=H(\omeg
@noindent
where @mymath{F(\omega)} is the Fourier transform of @mymath{f(l)}.
So multiplying the Fourier transform of two functions individually, we get the
Fourier transform of their convolution.
-The convolution theorem also proves a relation between the convolutions in the
frequency space. Let's define:
+The convolution theorem also proves a relation between the convolutions in the
frequency space.
+Let's define:
@dispmath{D(\omega){\equiv}F(\omega){\ast}H(\omega)}
@@ -14024,7 +14023,6 @@ The resulting @file{.fits.gz} file can then be fed into
any of Gnuastro's progra
-
@node MakeCatalog, Match, Segment, Data analysis
@section MakeCatalog
@@ -14740,6 +14738,11 @@ The final group of options particular to MakeCatalog
are those that specify whic
The current measurements in MakeCatalog are those which only produce one final
value for each label (for example its total brightness: a single number).
All the different label's measurements can be written as one column in a final
table/catalog that contains other columns for other similar single-number
measurements.
+In this case, all the different label's measurements can be written as one
column in a final table/catalog that contains other columns for other similar
single-number measurements.
+The majority of this section is devoted to MakeCatalog's single-valued
measurements.
+However, MakeCatalog can also do measurements that produce more than one value
for each label.
+Currently the only such measurement is generation of spectra from 3D cubes
with the @option{--spectrum} option and it is discussed in the end of this
section.
+
Command-line options are used to identify which measurements you want in the
final catalog(s) and in what order.
If any of the options below is called on the command line or in any of the
configuration files, it will be included as a column in the output catalog.
The order of the columns is in the same order as the options were seen by
MakeCatalog (see @ref{Configuration file precedence}).
@@ -14776,6 +14779,11 @@ You can use @option{--weightarea} to see which was
used.
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.
@@ -14783,6 +14791,9 @@ The geometric center is the average pixel positions
irrespective of their pixel
@item --geoy
The geometric center of all objects and clumps along the second FITS axis
axis, see @option{--geox}.
+@item --geoz
+The geometric center of all objects and clumps along the third FITS axis axis,
see @option{--geox}.
+
@item --minx
The minimum position of all objects and clumps along the first FITS axis.
@@ -14795,6 +14806,12 @@ The minimum position of all objects and clumps along
the second FITS axis.
@item --maxy
The maximum position of all objects and clumps along the second FITS axis.
+@item --minz
+The minimum position of all objects and clumps along the third FITS axis.
+
+@item --maxz
+The maximum position of all objects and clumps along the third FITS axis.
+
@item --clumpsx
[Objects] The flux weighted center of all the clumps in this object along the
first FITS axis.
See @option{--x}.
@@ -14803,6 +14820,10 @@ 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 the
first FITS axis.
See @option{--geox}.
@@ -14811,6 +14832,10 @@ 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 @option{--x}.
@@ -14833,6 +14858,11 @@ The first WCS axis is commonly used as right ascension
in images.
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 images.
@@ -14841,6 +14871,11 @@ The first WCS axis is commonly used as right ascension
in images.
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.
@@ -14849,6 +14884,10 @@ The first WCS axis is commonly used as right ascension
in images.
[Objects] Flux weighted declination 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, see
@option{--geox}.
The first WCS axis is commonly used as right ascension in images.
@@ -14857,6 +14896,11 @@ The first WCS axis is commonly used as right ascension
in images.
[Objects] Geometric center declination 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
magnitude}.
@@ -14968,7 +15012,15 @@ This is the square root of the mean variance under the
object, or the root mean
@item -a
@itemx --area
-The raw area (number of pixels) in any clump or object independent of what
pixel it lies over (if it is NaN/blank or unused for example).
+The raw area (number of pixels/voxels) in any clump or object independent of
what pixel it lies over (if it is NaN/blank or unused for example).
+
+@item --areaxy
+@cindex IFU: Integral Field Unit
+@cindex Integral Field Unit
+Similar to @option{--area}, when the clump or object is projected onto the
+first two dimensions. This is only available for 3-dimensional
+datasets. When working with Integral Field Unit (IFU) datasets, this
+projection onto the first two dimensions would be a narrow-band image.
@item --clumpsarea
[Objects] The total area of all the clumps in this object.
@@ -14981,6 +15033,12 @@ The area of all the pixels labeled with an object or
clump.
Note that unlike @option{--area}, pixel values are completely ignored in this
column.
For example, if a pixel value is blank, it won't be counted in
@option{--area}, but will be counted here.
+@item --geoareaxy
+Similar to @option{--geoarea}, when the clump or object is projected onto
+the first two dimensions. This is only available for 3-dimensional
+datasets. When working with Integral Field Unit (IFU) datasets, this
+projection onto the first two dimensions would be a narrow-band image.
+
@item -A
@itemx --semimajor
The pixel-value weighted root mean square (RMS) along the semi-major axis of
the profile (assuming it is an ellipse) in units of pixels.
@@ -15015,21 +15073,59 @@ The geometric (ignoring pixel values) angle of the
semi-major axis with the firs
@end table
+@cindex 3D data-cubes
+@cindex Cubes (3D data)
+@cindex IFU: Integral Field Unit
+@cindex Integral field unit (IFU)
+@cindex Spectrum (of astronomical source)
+Above, all of MakeCatalog's single-valued measurements were listed. As
+mentioned in the start of this section, MakeCatalog can also do
+multi-valued measurements per label. Currently the only such measurement is
+the creation of spectra from 3D data cubes as discussed below:
+
+@table @option
+@item --spectrum
+Generate a spectrum (measurement along the first two FITS dimensions) for
+each label when the input dataset is a 3D data cube. With this option, a
+seprate table/spectrum will be generated for every label. If the output is
+a FITS file, each label's spectrum will be written into an extension of
+that file with a standard name of @code{SPECTRUM_NN} (the label will be
+replaced with @code{NN}). If the output is a plain text file, each label's
+spectrum will be written into a separate file with the suffix
+@file{spec-NN.txt}. See @ref{MakeCatalog output} for more on specifying
+MakeCatalog's output file.
+
+The spectra will contain one row for every slice (third FITS dimension) of
+the cube. Since the physical nature of the third dimension is different,
+two types of spectra (along with their errors) are measured: 1) Sum of
+values in each slice that only have the requested label. 2) Sum of values
+on the 2D projection of the whole label (the area of this projection can be
+requested with the @option{--areaxy} column above).
+
+Labels can overlap when they are projected onto the first two FITS
+dimensions (the spatial domain). To help separate them, MakeCatalog does a
+third measurement on each slice: the area, sum of values and error of all
+pixels that belong to other labels but overlap with the 2D projection. This
+can be used to see how reliable the emission line measurement is (on the
+projected spectra) and also if multiple lines (labeled regions) belong to
+the same physical object.
+@end table
@node MakeCatalog output, , MakeCatalog measurements, Invoking astmkcatalog
@subsubsection MakeCatalog output
-
After it has completed all the requested measurements (see @ref{MakeCatalog
measurements}), MakeCatalog will store its measurements in table(s).
If an output filename is given (see @option{--output} in @ref{Input output
options}), the format of the table will be deduced from the name.
When it isn't given, the input name will be appended with a @file{_cat} suffix
(see @ref{Automatic output}) and its format will be determined from the
@option{--tableformat} option, which is also discussed in @ref{Input output
options}.
@option{--tableformat} is also necessary when the requested output name is a
FITS table (recall that FITS can accept ASCII and binary tables, see
@ref{Table}).
-By default only a single catalog/table will be created for ``objects'',
however, if @option{--clumpscat} is called, a secondary catalog/table will also
be created.
+By default (when @option{--spectrum} isn't called) only a single catalog/table
will be created for ``objects'', however, if @option{--clumpscat} is called, a
secondary catalog/table will also be created.
For more on ``objects'' and ``clumps'', see @ref{Segment}.
In short, if you only have one set of labeled images, you don't have to worry
about clumps (they are deactivated by default).
+When @option{--spectrum} is called, it is not mandatory to specify any
single-valued measurement columns. In this case, the output will only be the
spectra of each labeled region. See the description of @option{--spectrum} in
@ref{MakeCatalog measurements}.
+
The full list of MakeCatalog's output options are elaborated below.
@table @option
@@ -15128,6 +15224,10 @@ $ astmatch --aperture=1/3600,2/3600 in1.fits in2.txt
## and DEC_D columns of the second within a 0.5 arcseconds aperture.
$ astmatch --ccol1=RA,DEC --ccol2=RA_D,DEC_D --aperture=0.5/3600 \
in1.fits in2.fits
+
+## Match in 3D (RA, Dec and Wavelength).
+$ astmatch --ccol1=2,3,4 --ccol2=2,3,4 -a0.5/3600,0.5/3600,5e-10 \
+ in1.fits in2.txt
@end example
Match will find the rows that are nearest to each other in two catalogs (given
some coordinate columns).
@@ -15247,6 +15347,13 @@ Parameters of the aperture for matching.
The values given to this option can be fractions, for example when the
position columns are in units of degrees, @option{1/3600} can be used to ask
for one arcsecond.
The interpretation of the values depends on the requested dimensions
(determined from @option{--ccol1} and @code{--ccol2}) and how many values are
given to this option.
+When multiple objects are found within the aperture, the match is defined
+as the nearest one. In a multi-dimensional dataset, when the aperture is a
+general ellipse or ellipsoid (and not a circle or sphere), the distance is
+calculated in the elliptical space along the major axis. For the defintion
+of this distance, see @mymath{r_{el}} in @ref{Defining an ellipse and
+ellipsoid}.
+
@table @asis
@item 1D match
The aperture/interval can only take one value: half of the interval around
each point (maximum distance from each point).
@@ -15254,6 +15361,7 @@ The aperture/interval can only take one value: half of
the interval around each
@item 2D match
In a 2D match, the aperture can be a circle, an ellipse aligned in the axes or
an ellipse with a rotated major axis.
To simply the usage, you can determine the shape based on the number of free
parameters for each.
+
@table @asis
@item 1 number
For example @option{--aperture=2}.
@@ -15266,14 +15374,41 @@ The aperture will be an ellipse (if the two numbers
are different) with the resp
The numbers are in units of the first and second axis.
In the example above, the semi-axis value along the first axis will be 3 (in
units of the first coordinate) and along the second axis will be
@mymath{4\times10^{-10}} (in units of the second coordinate).
Such values can happen if you are comparing catalogs of a spectra for example.
-If more than one object exists in the aperture, the nearest will be found
along the major axis as described in @ref{Defining an ellipse}.
+If more than one object exists in the aperture, the nearest will be found
along the major axis as described in @ref{Defining an ellipse and ellipsoid}.
@item 3 numbers
For example @option{--aperture=2,0.6,30}.
The aperture will be an ellipse (if the second value is not 1).
The first number is the semi-major axis, the second is the axis ratio and the
third is the position angle (in degrees).
-If multiple matches are found within the ellipse, the distance (to find the
nearest) is calculated along the major axis in the elliptical space, see
@ref{Defining an ellipse}.
+If multiple matches are found within the ellipse, the distance (to find the
nearest) is calculated along the major axis in the elliptical space, see
@ref{Defining an ellipse and ellipsoid}.
@end table
+
+@item 3D match
+The aperture (matching volume) can be a sphere, an ellipsoid aligned on the
three axises or a genenral ellipsoid rotated in any direction.
+To simplifythe usage, the shape can be determined based on the number of
values given to this option.
+
+@table @asis
+@item 1 number
+For example @option{--aperture=3}.
+The matching volume will be a sphere of the given radius.
+The value is in the same units as the input coordinates.
+
+@item 3 numbers
+For example @option{--aperture=4,5,6e-10}.
+The aperture will be a general ellipsoid with the respective extent along each
dimension.
+The numbers must be in the same units as each axis.
+This is very similar to the two number case of 2D inputs.
+See there for more.
+
+@item 6 numbers
+For example @option{--aperture=4,0.5,0.6,10,20,30}.
+The numbers represent the full general ellipsoid definition (in any
orientation).
+For the definition of a general ellipsoid, see @ref{Defining an ellipse and
ellipsoid}.
+The first number is the semi-major axis.
+The second and third are the two axis ratios.
+The last three are the three Euler angles in units of degrees in the ZXZ order
as fully described in @ref{Defining an ellipse and ellipsoid}.
+@end table
+
@end table
@end table
@@ -15362,7 +15497,7 @@ In the subsections below, first a review of some very
basic information and conc
You can skip this subsection if you are already sufficiently familiar with
these concepts.
@menu
-* Defining an ellipse:: An ellipse on a pixelated grid.
+* Defining an ellipse and ellipsoid:: Definition of these important shapes.
* PSF:: Radial profiles for the PSF.
* Stars:: Making mock star profiles.
* Galaxies:: Radial profiles for galaxies.
@@ -15370,10 +15505,8 @@ You can skip this subsection if you are already
sufficiently familiar with these
* Oversampling:: Oversampling the model.
@end menu
-
-
-@node Defining an ellipse, PSF, Modeling basics, Modeling basics
-@subsubsection Defining an ellipse
+@node Defining an ellipse and ellipsoid, PSF, Modeling basics, Modeling basics
+@subsubsection Defining an ellipse and ellipsoid
@cindex Ellipse
@cindex Axis ratio
@@ -15399,6 +15532,33 @@ In that equation, the point is rotated, here the
coordinates are rotated and the
Hence, multiplying all elements of the ellipse definition with
@mymath{r_{el}^2} we get the elliptical distance at this point point located:
@mymath{r_{el}=\sqrt{i_r^2+(j_r/q)^2}}.
To place the radial profiles explained below over an ellipse,
@mymath{f(r_{el})} is calculated based on the functional radial profile desired.
+@cindex Ellipsoid
+@cindex Euler angles
+An ellipse in 3D, or an @url{https://en.wikipedia.org/wiki/Ellipsoid,
ellipsoid}, can be defined following similar principles as before.
+Labeling the major (largest) axis length as @mymath{a}, the second and third
(in a right-handed coordinate system) axis lengths can be labeled as @mymath{b}
and @mymath{c}.
+Hence we have two axis ratios: @mymath{q_1\equiv{b/a}} and
@mymath{q_2\equiv{c/a}}.
+The orientation of the ellipsoid can be defined from the orientation of its
major axis.
+There are many ways to define 3D orientation and order matters.
+So to be clear, here we use the ZXZ (or @mymath{Z_1X_2Z_3}) proper
@url{https://en.wikipedia.org/wiki/Euler_angles, Euler angles} to define the 3D
orientation.
+In short, when a point is rotated in this order, we first rotate it around the
Z axis (third axis) by @mymath{\alpha}, then about the (rotated) X axis by
@mymath{\beta} and finally about the (rotated) Z axis by @mymath{\gamma}.
+
+Following the discussion in @ref{Merging multiple warpings}, we can define the
full rotation with the following matrix multiplication.
+However, here we are rotating the coordinates, not the point.
+Therefore, both the rotation angles and rotation order are reversed.
+We are also not using homogeneous coordinates (see @ref{Warping basics}) since
we aren't concerned with translation in this context:
+
+@dispmath{\left[\matrix{i_r\cr j_r\cr k_r}\right] =
+ \left[\matrix{cos\gamma&sin\gamma&0\cr -sin\gamma&cos\gamma&0\cr
0&0&1}\right]
+ \left[\matrix{1&0&0\cr 0&cos\beta&sin\beta\cr
0&-sin\beta&cos\beta }\right]
+ \left[\matrix{cos\alpha&sin\alpha&0\cr -sin\alpha&cos\alpha&0\cr
0&0&1}\right]
+ \left[\matrix{i_c-i\cr j_c-j\cr k_c-k}\right] }
+
+@noindent
+Recall that an ellipsoid can be characterized with
+@mymath{(i_r/a)^2+(j_r/b)^2+(k_r/c)^2=1}, so similar to before
+(@mymath{r_{el}\equiv{a}}), we can find the ellipsoidal radius at pixel
+@mymath{(i,j,k)} as: @mymath{r_{el}=\sqrt{i_r^2+(j_r/q_1)^2+(k_r/q_2)^2}}.
+
@cindex Breadth first search
@cindex Inside-out construction
@cindex Making profiles pixel by pixel
@@ -15416,7 +15576,7 @@ Everything else after that (when the pixel index and
its radial profile have ent
-@node PSF, Stars, Defining an ellipse, Modeling basics
+@node PSF, Stars, Defining an ellipse and ellipsoid, Modeling basics
@subsubsection Point spread function
@cindex PSF
@@ -15568,7 +15728,7 @@ MacArthur et al.@footnote{MacArthur, L. A., S.
Courteau, and J. A. Holtzman (200
A pixel is the ultimate level of accuracy to gather data, we can't get any
more accurate in one image, this is known as sampling in signal processing.
However, the mathematical profiles which describe our models have infinite
accuracy.
Over a large fraction of the area of astrophysically interesting profiles (for
example galaxies or PSFs), the variation of the profile over the area of one
pixel is not too significant.
-In such cases, the elliptical radius (@mymath{r_{el}} of the center of the
pixel can be assigned as the final value of the pixel, see @ref{Defining an
ellipse}).
+In such cases, the elliptical radius (@mymath{r_{el}} of the center of the
pixel can be assigned as the final value of the pixel, see @ref{Defining an
ellipse and ellipsoid}).
@cindex Integration over pixel
@cindex Gradient over pixel area
@@ -15592,7 +15752,7 @@ Monte Carlo integration is first applied (which yields
@mymath{F_r}), then the c
If the fractional difference (@mymath{|F_r-F_c|/F_r}) is lower than a given
tolerance level (specified with @option{--tolerance}) MakeProfiles will stop
using Monte Carlo integration and only use the central pixel value.
@cindex Inside-out construction
-The ordering of the pixels in this inside-out construction is based on
@mymath{r=\sqrt{(i_c-i)^2+(j_c-j)^2}}, not @mymath{r_{el}}, see @ref{Defining
an ellipse}.
+The ordering of the pixels in this inside-out construction is based on
@mymath{r=\sqrt{(i_c-i)^2+(j_c-j)^2}}, not @mymath{r_{el}}, see @ref{Defining
an ellipse and ellipsoid}.
When the axis ratios are large (near one) this is fine.
But when they are small and the object is highly elliptical, it might seem
more reasonable to follow @mymath{r_{el}} not @mymath{r}.
The problem is that the gradient is stronger in pixels with smaller @mymath{r}
(and larger @mymath{r_{el}}) than those with smaller @mymath{r_{el}}.
@@ -15774,6 +15934,34 @@ The sections below discuss the options specific to
MakeProfiles based on context
Finally @ref{MakeProfiles output dataset} and @ref{MakeProfiles log file}
discuss the outputs of MakeProfiles and how you can configure them.
Besides these, MakeProfiles also supports all the common Gnuastro program
options that are discussed in @ref{Common options}, so please flip through them
is well for a more comfortable usage.
+When building 3D profiles, there are more degrees of freedom.
+Hence, more columns are necessary and all the values related to dimensions
(for example size of dataset in each dimension and the WCS properties) must
also have 3 values.
+To allow having an independent set of default values for creating 3D profiles,
MakeProfiles also installs a @file{astmkprof-3d.conf} configuration file (see
@ref{Configuration files}).
+You can use this for default 3D profile values.
+For example, if you installed Gnuastro with the prefix @file{/usr/local} (the
default location, see @ref{Installation directory}), you can benefit from this
configuration file by running MakeProfiles like the example below.
+As with all configuration files, if you want to customize a given option, call
it before the configuration file.
+
+@example
+$ astmkprof --config=/usr/local/etc/astmkprof-3d.conf catalog.txt
+@end example
+
+@cindex Shell alias
+@cindex Alias, shell
+@cindex Shell startup
+@cindex Startup, shell
+To further simplify the process, you can define a shell alias in any startup
file (for example @file{~/.bashrc}, see @ref{Installation directory}).
+Assuming that you installed Gnuastro in @file{/usr/local}, you can add this
line to the startup file (you may put it all in one line, it is broken into two
lines here for fitting within page limits).
+
+@example
+alias astmkprof-3d="astmkprof --config=/usr/local/etc/astmkprof-3d.conf"
+@end example
+
+@noindent
+Using this alias, you can call MakeProfiles with the name
@command{astmkprof-3d} (instead of @command{astmkprof}).
+It will automatically load the 3D specific configuration file first, and then
parse any other arguments, options or configuration files.
+You can change the default values in this 3D configuration file by calling
them on the command-line as you do with @command{astmkprof}@footnote{Recall
that for single-invocation options, the last command-line invocation takes
precedence over all previous invocations (including those in the 3D
configuration file).
+See the description of @option{--config} in @ref{Operating mode options}.}.
+
Please see @ref{Sufi simulates a detection} for a very complete tutorial
explaining how one could use MakeProfiles in conjunction with other Gnuastro's
programs to make a complete simulated image of a mock galaxy.
@menu
@@ -15850,7 +16038,7 @@ A fixed value will be used for all pixels less than or
equal to the truncation r
@item
Radial distance profile with `@code{distance}' or `@code{7}'.
-At the lowest level, each pixel only has an elliptical radial distance given
the profile's shape and orientation (see @ref{Defining an ellipse}).
+At the lowest level, each pixel only has an elliptical radial distance given
the profile's shape and orientation (see @ref{Defining an ellipse and
ellipsoid}).
When this profile is chosen, the pixel's elliptical radial distance from the
profile center is written as its value.
For this profile, the value in the magnitude column (@option{--mcol}) will be
ignored.
@@ -15867,9 +16055,30 @@ The S@'ersic index (@mymath{n}) or Moffat
@mymath{\beta}.
@item --pcol=STR/INT
The position angle (in degrees) of the profiles relative to the first FITS
axis (horizontal when viewed in SAO ds9).
+When building a 3D profile, this is the first Euler angle: first rotation of
the ellipsoid major axis from the first FITS axis (rotating about the third
axis).
+See @ref{Defining an ellipse and ellipsoid}.
+
+@item --p2col=STR/INT
+Second Euler angle (in degrees) when building a 3D ellipsoid.
+This is the second rotation of the ellipsoid major axis (following
@option{--pcol}) about the (rotated) X axis.
+See @ref{Defining an ellipse and ellipsoid}.
+This column is ignored when building a 2D profile.
+
+@item --p3col=STR/INT
+Third Euler angle (in degrees) when building a 3D ellipsoid.
+This is the third rotation of the ellipsoid major axis (following
@option{--pcol} and @option{--p2col}) about the (rotated) Z axis.
+See @ref{Defining an ellipse and ellipsoid}.
+This column is ignored when building a 2D profile.
@item --qcol=STR/INT
The axis ratio of the profiles (minor axis divided by the major axis in a 2D
ellipse).
+When building a 3D ellipse, this is the ratio of the major axis to the
semi-axis length of the second dimension (in a right-handed coordinate system).
+See @mymath{q1} in @ref{Defining an ellipse and ellipsoid}.
+
+@item --q2col=STR/INT
+The ratio of the ellipsoid major axis to the third semi-axis length (in a
right-handed coordinate system) of a 3D ellipsoid.
+See @mymath{q1} in @ref{Defining an ellipse and ellipsoid}.
+This column is ignored when building a 2D profile.
@item --mcol=STR/INT
The total pixelated magnitude of the profile within the truncation radius, see
@ref{Profile magnitude}.
@@ -16055,9 +16264,36 @@ Here are some examples:
A Moffat kernel with FWHM of 3 pixels, @mymath{\beta=2.8} which is truncated
at 5 times the FWHM.
@item --kernel=gaussian,2,3
-A Gaussian kernel with FWHM of 2 pixels and truncated at 3 times the FWHM.
+A circular Gaussian kernel with FWHM of 2 pixels and truncated at 3 times
+the FWHM.
@end table
+This option may also be used to create a 3D kernel.
+To do that, two small modifications are necessary: add a @code{-3d} (or
@code{-3D}) to the profile name (for example @code{moffat-3d}) and add a number
(axis-ratio along the third dimension) to the end of the parameters for all
profiles except @code{point}.
+The main reason behind providing an axis ratio in the third dimension is that
in 3D astronomical datasets, commonly the third dimension doesn't have the same
nature (units/sampling) as the first and second.
+
+For example in IFU datacubes, the first and second dimensions are
angularpositions (like RA and Dec) but the third is in units of Angstroms for
wavelength.
+Because of this different nature (which also affects theprocessing), it may be
necessary for the kernel to have a different extent in that direction.
+
+If the 3rd dimension axis ratio is equal to @mymath{1.0}, then the kernel will
be a spheroid.
+If its smaller than @mymath{1.0}, the kernel will be button-shaped: extended
less in the third dimension.
+However, when it islarger than @mymath{1.0}, the kernel will be bullet-shaped:
extended more in the third dimension.
+In the latter case, the radial parameter will correspond to the length along
the 3rd dimension.
+For example, let's have a look at the two examples above but in 3D:
+
+@table @option
+@item --kernel=moffat-3d,3,2.8,5,0.5
+An ellipsoid Moffat kernel with FWHM of 3 pixels, @mymath{\beta=2.8} which is
truncated at 5 times the FWHM.
+The ellipsoid is circular in the first two dimensions, but in the third
dimension its extent is half the first two.
+
+@item --kernel=gaussian-3d,2,3,1
+A spherical Gaussian kernel with FWHM of 2 pixels and truncated at 3 times
+the FWHM.
+@end table
+
+Ofcourse, if a specific kernel is needed that doesn't fit the constraints
imposed by this option, you can always use a catalog to define any arbitrary
kernel.
+Just call the @option{--individual} and @option{--nomerged} options to make
sure that it is built as a separate file (individually) and no ``merged'' image
of the input profiles is created.
+
@item -x INT,INT
@itemx --mergedsize=INT,INT
The number of pixels along each axis of the output, in FITS order.
@@ -22591,6 +22827,39 @@ in the center of the box, all the values in
@code{width} will be an odd
integer.
@end deftypefun
+@deftypefun void gal_box_bound_ellipsoid_extent (double @code{*semiaxes},
double @code{*euler_deg}, double @code{*extent})
+Return the maximum extent along each dimension of the given ellipsoid from
+its center. Therefore this is half the extent of the box in each
+dimension. The semi-axis lengths of the ellipsoid must be present in the 3
+element @code{semiaxis} array. The @code{euler_deg} array contains the
+three ellipsoid Euler angles in degrees. For a description of the Euler
+angles, see description of @code{gal_box_bound_ellipsoid} below. The extent
+in each dimension is in floating point format and stored in @code{extent}
+which must already be allocated before this function.
+@end deftypefun
+
+@deftypefun void gal_box_bound_ellipsoid (double @code{*semiaxes}, double
@code{*euler_deg}, long @code{*width})
+Any ellipsoid can be enclosed into a rectangular volume/box. The purpose of
+this function is to give the integer size/width of that box. The semi-axes
+lengths of the ellipse must be in the @code{semiaxes} array (with three
+elements). The major axis length must be the first element of
+@code{semiaxes}. The only other condition is that the next two semi-axes
+must both be smaller than the first. The orientation of the major axis is
+defined through three proper Euler angles (ZXZ order in degrees) that are
+given in the @code{euler_deg} array. The @code{width} array will contain
+the output size in long integer type (in FITS axis order). Since the
+ellipsoid center is assumed to be in the center of the box, all the values
+in @code{width} will be an odd integer.
+
+@cindex Euler angles
+The proper Euler angles can be defined in many ways (which axes to rotate
+about). For a full description of the Euler angles, please see
+@url{https://en.wikipedia.org/wiki/Euler_angles, Wikipedia}. Here we adopt
+the ZXZ (or @mymath{Z_1X_2Z_3}) proper Euler angles were the first rotation
+is done around the Z axis, the second one about the (rotated) X axis and
+the third about the (rotated) Z axis.
+@end deftypefun
+
@deftypefun void gal_box_border_from_center (double @code{center}, size_t
@code{ndim}, long @code{*width}, long @code{*fpixel}, long @code{*lpixel})
Given the center coordinates in @code{center} and the @code{width} (along
each dimension) of a box, return the coordinates of the first
@@ -22941,11 +23210,11 @@ points of one catalog lie within this aperture of a
point in the other, the
nearest is defined as the match. In a 2D situation (where the input lists
have two nodes), for the most generic case, it must have three elements:
the major axis length, axis ratio and position angle (see @ref{Defining an
-ellipse}). If @code{aperture[1]==1}, the aperture will be a circle of
-radius @code{aperture[0]} and the third value won't be used. When the
-aperture is an ellipse, distances between the points are also calculated in
-the respective elliptical distances (@mymath{r_{el}} in @ref{Defining an
-ellipse}).
+ellipse and ellipsoid}). If @code{aperture[1]==1}, the aperture will be a
+circle of radius @code{aperture[0]} and the third value won't be used. When
+the aperture is an ellipse, distances between the points are also
+calculated in the respective elliptical distances (@mymath{r_{el}} in
+@ref{Defining an ellipse and ellipsoid}).
To speed up the search, this function will sort the input coordinates by
their first column (first axis). If @emph{both} are already sorted by their
diff --git a/lib/binary.c b/lib/binary.c
index bfec6c7..a99461e 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -143,7 +143,7 @@ binary_erode_dilate_2d_8con(gal_data_t *input, unsigned
char dilate0_erode1)
/* Set the foreground and background values: */
if(dilate0_erode1==0) {f=1; b=0;}
- else {f=0; b=1;}
+ else {f=0; b=1;}
/* Check the 4 corners: */
if(byt[0]==b && (byt[1]==f
@@ -221,6 +221,44 @@ binary_erode_dilate_2d_8con(gal_data_t *input, unsigned
char dilate0_erode1)
+/* This is a general erosion and dilation function. It is less efficient
+ than the more specialized cases above. */
+static void
+binary_erode_dilate_general(gal_data_t *input, unsigned char dilate0_erode1,
+ int connectivity)
+{
+ uint8_t f, b, *pt, *fpt, *byt=input->array;
+ size_t i, *dinc=gal_dimension_increment(input->ndim, input->dsize);
+
+ /* Do a sanity check: */
+ if(dilate0_erode1!=1 && dilate0_erode1!=0)
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can fix "
+ "this problem. The value to dilate0_erode1 is %u while it should "
+ "be 0 or 1", __func__, PACKAGE_BUGREPORT, dilate0_erode1);
+
+ /* Set the foreground and background values: */
+ if(dilate0_erode1==0) {f=1; b=0;}
+ else {f=0; b=1;}
+
+ /* Go over the neighbors of each pixel. */
+ for(i=0;i<input->size;++i)
+ if(byt[i]==b)
+ GAL_DIMENSION_NEIGHBOR_OP(i, input->ndim, input->dsize, connectivity,
+ dinc,{
+ if(byt[i]!=GAL_BINARY_TMP_VALUE
+ && byt[nind]==f)
+ byt[i]=GAL_BINARY_TMP_VALUE;
+ });
+
+ /* Set all the changed pixels to the proper values: */
+ fpt=(pt=byt)+input->size;
+ do *pt = *pt==GAL_BINARY_TMP_VALUE ? f : *pt; while(++pt<fpt);
+}
+
+
+
+
+
/* Erode a binary dataset any number of times. If `inplace' is given a
value of `1', then do the erosion within the already allocated space,
otherwise, allocate a new array and save the result into that.
@@ -246,8 +284,8 @@ binary_erode_dilate(gal_data_t *input, size_t num, int
connectivity,
/* Set the dataset to work on. */
binary = ( (inplace && input->type==GAL_TYPE_UINT8)
- ? input :
- gal_data_copy_to_new_type(input, GAL_TYPE_UINT8) );
+ ? input
+ : gal_data_copy_to_new_type(input, GAL_TYPE_UINT8) );
/* Go over every element and do the erosion. */
switch(binary->ndim)
@@ -264,6 +302,11 @@ binary_erode_dilate(gal_data_t *input, size_t num, int
connectivity,
}
break;
+ case 3:
+ for(counter=0;counter<num;++counter)
+ binary_erode_dilate_general(binary, d0e1, connectivity);
+ break;
+
default:
error(EXIT_FAILURE, 0, "%s: currently doesn't work on %zu "
"dimensional datasets", __func__, binary->ndim);
diff --git a/lib/box.c b/lib/box.c
index be97e3c..1a48f71 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -39,10 +39,10 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/* Any ellipse can be enclosed into a rectangular box. The purpose of
- this function is to give the height and width of that box. The
- logic behind it is this: All the points on the circumference of an
- ellipse that is aligned on the x axis can be written as:
+/* Any ellipse can be enclosed into a rectangular box. The purpose of this
+ function is to give the height and width of that box. The logic behind
+ it is this: All the points on the circumference of an ellipse that is
+ aligned on the x axis can be written as:
(acos(t),bsin(t)) where 0<t<2\pi. (1)
@@ -100,6 +100,174 @@ gal_box_bound_ellipse(double a, double b, double
theta_deg, long *width)
+/* Find the bounding box of an ellipsoid. An ellipsoid is defined by its
+ three axises: the first (a) must be the major axis, the other two must
+ be smaller than `a' but no particular relation between them is
+ assumed. We will define the orientation of the ellipsoid from its major
+ axis and use the "Proper Euler angles" (ZXZ order) to define that
+ orientation.
+
+ The following derivation is taken from:
+
+ https://tavianator.com/exact-bounding-boxes-for-spheres-ellipsoids/
+
+ For the definition of the Euler angles (and the ZXZ order) see Wikipedia
+ in the link below. In short: first rotate around the Z axis, then the
+ (rotated) first axis, then the (rotated) Z axis.
+
+ https://en.wikipedia.org/wiki/Euler_angles
+
+ Defining the general point `p' as (the transpose of `p' with `p^T' and
+ its inverse with `p^-1'):
+
+ | x |
+ p = | y |
+ | z |
+ | 1 |
+
+ A sphere satisfies the following homogenous coordinate matrix and
+ general quadratic surface:
+
+ | 1 0 0 0 |
+ S = | 0 1 0 0 | --> p^T * S * p = 0
+ | 0 0 1 0 |
+ | 0 0 0 -1 |
+
+ The conversion from a sphere to an arbitrarily oriented ellipsoid will
+ be done by rotating the three semi-axes lengths and merging the three
+ vertical vectors into one matrix. The rotation can be written as the
+ following combined affine transformation (only in the three main axises
+ (since we aren't dealing with translation here). Here, we'll call
+ `sin(alpha)' as `s1', `cos(beta)' as `c2' and `sin(gamma)' as `s3'.
+
+ Rotate by Euler angles
+ ----------------------------
+ | c1 -s1 0 | | 1 0 0 | | c3 -s3 0 |
+ R = | s1 c1 0 | * | 0 c2 -s2 | * | s3 c3 0 |
+ | 0 0 1 | | 0 s2 c2 | | 0 0 1 |
+
+ Then `M' (rotation and scaling to obtain ellipsoid from sphere) will be:
+
+ | a | | 0 | | 0 | | A1 B1 C1 |
+ A = R * | 0 |, B = R * | b |, C = R * | 0 | --> M = | A2 B2 C2 |
+ | 0 | | 0 | | c | | A3 B3 C3 |
+
+ The final result is:
+
+ | a*c1*c3-a*s1*c2*s3 -b*c1*s3-b*s1*c2*c3 c*s1*s2 |
+ M = | a*s1*c3+a*c1*c2*s3 -b*s1*s3+b*c1*c2*c3 -c*c1*s2 |
+ | a*s2*s3 b*s2*c3 c*c2 |
+
+ The quadratic surface for this ellipsoid can now be written as:
+
+ (M^-1 * p)^T * S * (M^-1 * p) = 0 --> p^T * (M^-T * S * M^-1) * p = 0
+
+ Writing Q = M^-T * S * M^-1, we get: `p^T * Q * p = 0'. Now, we define a
+ plane with a horizontal vector `u = [a b c d ]', such that `u.p=0'. For
+ a point on the ellipsoid (at `p') we have: `u^T=p^T * Q'. This is
+ because: `u.p = u^T * p = p^T * Q * p = 0' (as we showed above).
+
+ Tangent planes will have the following useful property:
+
+ u^T * Q^-1 * u = p^T * Q * Q^-1 * Q * p = p^T * Q * p = 0
+
+ Now, the particular plane that is perpendicular to the X axis has the
+ general form: `u = [ 1 0 0 -x ]'. So, defining `R = Q^1', and using the
+ property above for tangential planes, we can find the X axis position.
+
+ However, getting to `R' from `M' as described above is not easy. So,
+ taking the following considerations into account, we can derive the
+ final values: [from that webpage] "Several details of the problem can
+ make computing the planes more efficient. The first is that `S' is
+ involutory, meaning `S^-1 = S'. This means that the product `M * S^-1'
+ can be computed implicitly: it is simply `M' with its last column
+ negated. The last column of `R = M * S^-1 * MT' is the same, because the
+ last column of `M^T' is `[ 0 0 0 1 ]'. In particular, `R[4,4]=-1'.
+
+ Not all values of RR are used; in fact, only values from the last column
+ and the diagonal appear in the formulae. We know the last column
+ implicitly, and the diagonal has the formula: "
+
+ R[i,i] = (sum_{j=1 to 3} M[i,j]^2 - M_[i,j]
+
+ So the bounding box lengths along each dimension are the
+ following. Recall that in homogenous coordinates, the last column is for
+ translation. So in the case of this function all the `M[i,4]' values are
+ zero.
+
+ x = M[1,4] \pm sqrt( M[1,1]^2 + M[1,2]^2 + M[1,3]^2 )
+ y = M[2,4] \pm sqrt( M[2,1]^2 + M[2,2]^2 + M[2,3]^2 )
+ z = M[3,4] \pm sqrt( M[3,1]^2 + M[3,2]^2 + M[3,3]^2 ) */
+void
+gal_box_bound_ellipsoid_extent(double *semiaxes, double *euler_deg,
+ double *extent)
+{
+ size_t i, j;
+ double a=semiaxes[0], b=semiaxes[1], c=semiaxes[2];
+ double c1=cos(euler_deg[0]*M_PI/180), s1=sin(euler_deg[0]*M_PI/180);
+ double c2=cos(euler_deg[1]*M_PI/180), s2=sin(euler_deg[1]*M_PI/180);
+ double c3=cos(euler_deg[2]*M_PI/180), s3=sin(euler_deg[2]*M_PI/180);
+ double R[9]={ a*c1*c3-a*s1*c2*s3, -1.0f*b*c1*s3-b*s1*c2*c3, c*s1*s2,
+ a*s1*c3+a*c1*c2*s3, -1.0f*b*s1*s3+b*c1*c2*c3, -1.0f*c*c1*s2,
+ a*s2*s3, b*s2*c3, c*c2 };
+
+ /* Sanity check. */
+ if(b>a || c>a)
+ error(EXIT_FAILURE, 0, "%s: the second and third semi-axes lengths "
+ "(%g, %g respectively) must both be smaller or equal to the "
+ "first (%g)", __func__, b, c, a);
+
+ /* Find the width along each dimension. */
+ for(i=0;i<3;++i)
+ {
+ extent[i]=0.0;
+ for(j=0;j<3;++j)
+ extent[i] += R[i*3+j] * R[i*3+j];
+ extent[i] = sqrt(extent[i]);
+ }
+
+ /* For a check:
+ printf("Axies: %g, %g, %g\n", a, b, c);
+ printf("Euler angles: %g, %g, %g\n", euler_deg[0], euler_deg[1],
+ euler_deg[2]);
+ printf("c1: %f, s1: %f\nc2: %f, s2: %f\nc3: %f, s3: %f\n",
+ c1, s1, c2, s2, c3, s3);
+ PRINT3BY3("R", R);
+ printf("extent: %ld, %ld, %ld\n", extent[0], extent[1], extent[2]);
+ exit(0);
+ */
+}
+
+
+
+
+
+/* Using `gal_box_bound_ellipsoid_extent', find the integer width of a box
+ that contains the ellipsoid. */
+#define PRINT3BY3(C, A){ \
+ printf("%s: | %-15g%-15g%-15g |\n" \
+ " | %-15g%-15g%-15g |\n" \
+ " | %-15g%-15g%-15g |\n\n", (C), (A)[0], (A)[1], (A)[2], \
+ (A)[3], (A)[4], (A)[5], (A)[6], (A)[7], (A)[8]); \
+ }
+void
+gal_box_bound_ellipsoid(double *semiaxes, double *euler_deg, long *width)
+{
+ size_t i;
+ double extent[3];
+
+ /* Find the extent of the ellipsoid in each axis. */
+ gal_box_bound_ellipsoid_extent(semiaxes, euler_deg, extent);
+
+ /* Find the width along each dimension. */
+ for(i=0;i<3;++i)
+ width[i] = 2 * ( (long)extent[i] ) + 1;
+}
+
+
+
+
+
/* We have the central pixel and box size of the crop box, find the
starting (first) and ending (last) pixels: */
void
diff --git a/lib/data.c b/lib/data.c
index d09c436..fd112cc 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -800,7 +800,6 @@ gal_data_copy_to_allocated(gal_data_t *in, gal_data_t *out)
gal_checkset_allocate_copy(in->unit, &out->unit);
gal_checkset_allocate_copy(in->comment, &out->comment);
-
/* Do the copying. */
switch(out->type)
{
diff --git a/lib/gnuastro-internal/options.h b/lib/gnuastro-internal/options.h
index 64c4785..dd5199c 100644
--- a/lib/gnuastro-internal/options.h
+++ b/lib/gnuastro-internal/options.h
@@ -55,6 +55,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
+
/* Standard groups of options. From the Argp manual (in GNU C Library): "In
a long help message, options are sorted alphabetically within each
group, and the groups presented in the order 0, 1, 2, ..., N, -M, ...,
diff --git a/lib/gnuastro/box.h b/lib/gnuastro/box.h
index c41a1a8..ee79a83 100644
--- a/lib/gnuastro/box.h
+++ b/lib/gnuastro/box.h
@@ -59,6 +59,13 @@ void
gal_box_bound_ellipse(double a, double b, double theta_deg, long *width);
void
+gal_box_bound_ellipsoid_extent(double *semiaxes, double *euler_deg,
+ double *extent);
+
+void
+gal_box_bound_ellipsoid(double *semiaxes, double *euler_deg, long *width);
+
+void
gal_box_border_from_center(double *center, size_t ndim, long *width,
long *fpixel, long *lpixel);
diff --git a/lib/label.c b/lib/label.c
index 791b8e2..895811d 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -557,8 +557,8 @@ label_clump_significance_sanity(gal_data_t *values,
gal_data_t *std,
gal_type_name(label->type, 1));
/* Dimentionality of the values dataset. */
- if( values->ndim>2 )
- error(EXIT_FAILURE, 0, "%s: currently only supports 1 or 2 "
+ if( values->ndim>3 )
+ error(EXIT_FAILURE, 0, "%s: currently only supports 1, 2 or 3 "
"dimensional datasets, but a %zu-dimensional dataset is given",
func, values->ndim);
@@ -660,7 +660,7 @@ label_clump_significance_raw(gal_data_t *values_d,
gal_data_t *std_d,
size_t ndim=values_d->ndim, *dsize=values_d->dsize;
double *row;
- size_t i, *a, *af, ii, coord[2];
+ size_t i, *a, *af, ii, coord[3];
size_t nngb=gal_dimension_num_neighbors(ndim);
int32_t nlab, *ngblabs, *label=label_d->array;
float *values=values_d->array, *std=std_d->array;
diff --git a/lib/match.c b/lib/match.c
index da59e41..a69fb86 100644
--- a/lib/match.c
+++ b/lib/match.c
@@ -166,10 +166,10 @@ match_coordinaes_sanity_check(gal_data_t *coord1,
gal_data_t *coord2,
gal_list_data_number(coord2));
- /* This function currently only works on 1 and 2 dimensional inputs. */
- if(ncoord1>2)
+ /* This function currently only works for less than 4 dimensions. */
+ if(ncoord1>3)
error(EXIT_FAILURE, 0, "%s: %zu dimension matching requested, this "
- "function currently only matches datasets with a maximum of 2 "
+ "function currently only matches datasets with a maximum of 3 "
"dimensions", __func__, ncoord1);
/* Check the column properties. */
@@ -180,11 +180,32 @@ match_coordinaes_sanity_check(gal_data_t *coord1,
gal_data_t *coord2,
if(aperture[0]<=0)
error(EXIT_FAILURE, 0, "%s: the first value in the aperture (%g) cannot "
"be zero or negative", __func__, aperture[0]);
- if(ncoord1==2)
- if(aperture[1]<=0 || aperture[1]>1)
- error(EXIT_FAILURE, 0, "%s: the second value in the aperture (%g) "
- "is the axis ratio, so it must be larger than zero and less "
- "than 1", __func__, aperture[1]);
+ switch(ncoord1)
+ {
+ case 1: /* We don't need any checks in a 1D match. */
+ break;
+
+ case 2:
+ if( (aperture[1]<=0 || aperture[1]>1))
+ error(EXIT_FAILURE, 0, "%s: the second value in the aperture (%g) "
+ "is the axis ratio, so it must be larger than zero and less "
+ "than 1", __func__, aperture[1]);
+ break;
+
+ case 3:
+ if(aperture[1]<=0 || aperture[1]>1 || aperture[2]<=0 || aperture[2]>1)
+ error(EXIT_FAILURE, 0, "%s: at least one of the second or third "
+ "values in the aperture (%g and %g respectively) is smaller "
+ "than zero or larger than one. In a 3D match, these are the "
+ "axis ratios, so they must be larger than zero and less than "
+ "1", __func__, aperture[1], aperture[2]);
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the issue. The value %zu not recognized for `ndim'", __func__,
+ PACKAGE_BUGREPORT, ncoord1);
+ }
}
@@ -297,6 +318,89 @@ match_coordinates_prepare(gal_data_t *coord1, gal_data_t
*coord2,
/********************************************************************/
/************* Coordinate matching *************/
/********************************************************************/
+/* Preparations for `match_coordinates_second_in_first'. */
+static void
+match_coordinates_sif_prepare(gal_data_t *A, gal_data_t *B,
+ double *aperture, size_t ndim, double **a,
+ double **b, double *dist, double *c,
+ double *s, int *iscircle)
+{
+ double semiaxes[3];
+
+ /* These two are common for all dimensions. */
+ a[0]=A->array;
+ b[0]=B->array;
+
+ /* See if the aperture is a circle or not. */
+ switch(ndim)
+ {
+ case 1:
+ *iscircle = 0; /* Irrelevant for 1D. */
+ dist[0]=aperture[0];
+ break;
+
+ case 2:
+ /* Set the main coordinate arrays. */
+ a[1]=A->next->array;
+ b[1]=B->next->array;
+
+ /* See if the aperture is circular. */
+ if( (*iscircle=(aperture[1]==1)?1:0)==0 )
+ {
+ /* Using the box that encloses the aperture, calculate the
+ distance along each axis. */
+ gal_box_bound_ellipse_extent(aperture[0], aperture[0]*aperture[1],
+ aperture[2], dist);
+
+ /* Calculate the sin and cos of the given ellipse if necessary
+ for ease of processing later. */
+ c[0] = cos( aperture[2] * M_PI/180.0 );
+ s[0] = sin( aperture[2] * M_PI/180.0 );
+ }
+ else
+ dist[0]=dist[1]=aperture[0];
+ break;
+
+ case 3:
+ /* Set the main coordinate arrays. */
+ a[1]=A->next->array;
+ b[1]=B->next->array;
+ a[2]=A->next->next->array;
+ b[2]=B->next->next->array;
+
+ if( (*iscircle=(aperture[1]==1 && aperture[2]==1)?1:0)==0 )
+ {
+ /* Using the box that encloses the aperture, calculate the
+ distance along each axis. */
+ semiaxes[0]=aperture[0];
+ semiaxes[1]=aperture[1]*aperture[0];
+ semiaxes[2]=aperture[2]*aperture[0];
+ gal_box_bound_ellipsoid_extent(semiaxes, &aperture[3], dist);
+
+ /* Calculate the sin and cos of the given ellipse if necessary
+ for ease of processing later. */
+ c[0] = cos( aperture[3] * M_PI/180.0 );
+ s[0] = sin( aperture[3] * M_PI/180.0 );
+ c[1] = cos( aperture[4] * M_PI/180.0 );
+ s[1] = sin( aperture[4] * M_PI/180.0 );
+ c[2] = cos( aperture[5] * M_PI/180.0 );
+ s[2] = sin( aperture[5] * M_PI/180.0 );
+ }
+ else
+ dist[0]=dist[1]=dist[2]=aperture[0];
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
+ "the problem. The value %zu is not recognized for ndim",
+ __func__, PACKAGE_BUGREPORT, ndim);
+ }
+}
+
+
+
+
+
static double
match_coordinates_elliptical_r_2d(double d1, double d2, double *ellipse,
double c, double s)
@@ -311,8 +415,29 @@ match_coordinates_elliptical_r_2d(double d1, double d2,
double *ellipse,
static double
+match_coordinates_elliptical_r_3d(double *delta, double *ellipsoid,
+ double *c, double *s)
+{
+ double Xr, Yr, Zr;
+ double c1=c[0], s1=s[0];
+ double c2=c[1], s2=s[1];
+ double c3=c[2], s3=s[2];
+ double q1=ellipsoid[1], q2=ellipsoid[2];
+ double x=delta[0], y=delta[1], z=delta[2];
+
+ Xr = x*( c3*c1 - s3*c2*s1 ) + y*( c3*s1 + s3*c2*c1) + z*( s3*s2 );
+ Yr = x*( -1*s3*c1 - c3*c2*s1 ) + y*(-1*s3*s1 + c3*c2*c1) + z*( c3*s2 );
+ Zr = x*( s1*s2 ) + y*(-1*s2*c1 ) + z*( c2 );
+ return sqrt( Xr*Xr + Yr*Yr/q1/q1 + Zr*Zr/q2/q2 );
+}
+
+
+
+
+
+static double
match_coordinates_distance(double *delta, int iscircle, size_t ndim,
- double *aperture, double c, double s)
+ double *aperture, double *c, double *s)
{
/* For more than one dimension, we'll need to calculate the distance from
the deltas (differences) in each dimension. */
@@ -325,7 +450,15 @@ match_coordinates_distance(double *delta, int iscircle,
size_t ndim,
return ( iscircle
? sqrt( delta[0]*delta[0] + delta[1]*delta[1] )
: match_coordinates_elliptical_r_2d(delta[0], delta[1],
- aperture, c, s) );
+ aperture, c[0], s[0]) );
+
+ case 3:
+ return ( iscircle
+ ? sqrt( delta[0]*delta[0]
+ + delta[1]*delta[1]
+ + delta[2]*delta[2] )
+ : match_coordinates_elliptical_r_3d(delta, aperture, c, s) );
+
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
"the problem. The value %zu is not recognized for ndim",
@@ -356,41 +489,17 @@ match_coordinates_second_in_first(gal_data_t *A,
gal_data_t *B,
redundant variables (those that equal a previous value) are only
defined to make it easy to read the code.*/
int iscircle=0;
- double r, c=NAN, s=NAN;
size_t i, ar=A->size, br=B->size;
size_t ai, bi, blow=0, prevblow=0;
size_t ndim=gal_list_data_number(A);
- double dist[2]={NAN,NAN}, delta[2]={NAN,NAN};
- double *a[2]={A->array, A->next ? A->next->array : NULL};
- double *b[2]={B->array, B->next ? B->next->array : NULL};
+ double r, c[3]={NAN, NAN, NAN}, s[3]={NAN, NAN, NAN};
+ double dist[3]={NAN, NAN, NAN}, delta[3]={NAN, NAN, NAN};
+ double *a[3]={NULL, NULL, NULL}, *b[3]={NULL, NULL, NULL};
- /* See if the aperture is a circle or not. */
- switch(ndim)
- {
- case 1: iscircle = 0; /* Irrelevant for 1D. */ break;
- case 2: iscircle = aperture[1]==1 ? 1 : 0; break;
- default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
- "the problem. The value %zu is not recognized for ndim",
- __func__, PACKAGE_BUGREPORT, ndim);
- }
-
- /* Preparations for the shape of the aperture. */
- if(ndim==1 || iscircle)
- dist[0]=dist[1]=aperture[0];
- else
- {
- /* Using the box that encloses the aperture, calculate the distance
- along each axis. */
- gal_box_bound_ellipse_extent(aperture[0], aperture[0]*aperture[1],
- aperture[2], dist);
-
- /* Calculate the sin and cos of the given ellipse if necessary for
- ease of processing later. */
- c = cos( aperture[2] * M_PI/180.0 );
- s = sin( aperture[2] * M_PI/180.0 );
- }
+ /* Necessary preperations. */
+ match_coordinates_sif_prepare(A, B, aperture, ndim, a, b, dist, c, s,
+ &iscircle);
/* For each row/record of catalog `a', make a list of the nearest records
in catalog b within the maximum distance. Note that both catalogs are
@@ -415,7 +524,7 @@ match_coordinates_second_in_first(gal_data_t *A, gal_data_t
*B,
a search here to change `blow' if necessary before doing further
searching.*/
for( blow=prevblow; blow<br && b[0][blow] < a[0][ai]-dist[0]; ++blow)
- {/* This can be blank, the `for' does all we need :-). */}
+ { /* This can be blank, the `for' does all we need :-). */ }
/* `blow' is now found for this `ai' and will be used unchanged to
@@ -440,7 +549,7 @@ match_coordinates_second_in_first(gal_data_t *A, gal_data_t
*B,
coordinates have EXACTLY the same floating point value as
each other to easily define an independent sorting in the
second axis. */
- if( ndim==1
+ if( ndim<2
|| ( b[1][bi] >= a[1][ai]-dist[1]
&& b[1][bi] <= a[1][ai]+dist[1] ) )
{
@@ -472,11 +581,16 @@ match_coordinates_second_in_first(gal_data_t *A,
gal_data_t *B,
The next two problems will be solved with the list after
parsing of the whole catalog is complete.*/
- for(i=0;i<ndim;++i) delta[i]=b[i][bi]-a[i][ai];
- r=match_coordinates_distance(delta, iscircle, ndim, aperture,
- c, s);
- if(r<aperture[0])
- match_coordinate_add_to_sfll(&bina[ai], bi, r);
+ if( ndim<3
+ || ( b[2][bi] >= a[2][ai]-dist[2]
+ && b[2][bi] <= a[2][ai]+dist[2] ) )
+ {
+ for(i=0;i<ndim;++i) delta[i]=b[i][bi]-a[i][ai];
+ r=match_coordinates_distance(delta, iscircle, ndim,
+ aperture, c, s);
+ if(r<aperture[0])
+ match_coordinate_add_to_sfll(&bina[ai], bi, r);
+ }
}
}
diff --git a/lib/options.c b/lib/options.c
index 4943d28..b71d8c3 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -921,7 +921,6 @@ gal_options_parse_sizes_reverse(struct argp_option *option,
char *arg,
"values to the `--%s' option must be positive",
arg, v[i], option->name);
-
if(ceil(v[i]) != v[i])
error_at_line(EXIT_FAILURE, 0, filename, lineno, "a given "
"value in `%s' (%g) is not an integer. The "
diff --git a/lib/wcs.c b/lib/wcs.c
index e1e8526..d0f9a77 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -990,8 +990,9 @@ gal_wcs_world_to_img(gal_data_t *coords, struct wcsprm
*wcs, int inplace)
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]);
+ 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]);
}
*/
@@ -1051,8 +1052,9 @@ gal_wcs_img_to_world(gal_data_t *coords, struct wcsprm
*wcs, int inplace)
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]);
+ 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]);
}
*/
diff --git a/tests/Makefile.am b/tests/Makefile.am
index c3ac5ec..f567a0f 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -96,8 +96,8 @@ if COND_COSMICCAL
cosmiccal/simpletest.sh: prepconf.sh.log
endif
if COND_CROP
- MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/imgcenter.sh \
- crop/imgcenternoblank.sh crop/section.sh crop/wcscenter.sh \
+ MAYBE_CROP_TESTS = crop/imgcat.sh crop/wcscat.sh crop/imgcenter.sh \
+ crop/imgcenternoblank.sh crop/section.sh crop/wcscenter.sh \
crop/imgpolygon.sh crop/imgoutpolygon.sh crop/wcspolygon.sh
crop/imgcat.sh: mkprof/mosaic1.sh.log
@@ -130,7 +130,7 @@ if COND_MATCH
match/merged-cols.sh: prepconf.sh.log
endif
if COND_MKCATALOG
- MAYBE_MKCATALOG_TESTS = mkcatalog/detections.sh \
+ MAYBE_MKCATALOG_TESTS = mkcatalog/detections.sh \
mkcatalog/objects-clumps.sh mkcatalog/aperturephot.sh
mkcatalog/objects-clumps.sh: segment/segment.sh.log
@@ -139,20 +139,24 @@ if COND_MKCATALOG
mkprof/clearcanvas.sh.log
endif
if COND_MKNOISE
- MAYBE_MKNOISE_TESTS = mknoise/addnoise.sh
+ MAYBE_MKNOISE_TESTS = mknoise/addnoise.sh mknoise/addnoise-3d.sh
mknoise/addnoise.sh: warp/warp_scale.sh.log
+ mknoise/addnoise-3d.sh: mkprof/3d-cat.sh.log
endif
if COND_MKPROF
- MAYBE_MKPROF_TESTS = mkprof/mosaic1.sh mkprof/mosaic2.sh \
- mkprof/mosaic3.sh mkprof/mosaic4.sh mkprof/radeccat.sh \
- mkprof/ellipticalmasks.sh mkprof/clearcanvas.sh
+ MAYBE_MKPROF_TESTS = mkprof/mosaic1.sh mkprof/mosaic2.sh \
+ mkprof/mosaic3.sh mkprof/mosaic4.sh mkprof/radeccat.sh \
+ mkprof/ellipticalmasks.sh mkprof/clearcanvas.sh mkprof/3d-cat.sh \
+ mkprof/3d-kernel.sh
+ mkprof/3d-cat.sh: prepconf.sh.log
mkprof/mosaic1.sh: prepconf.sh.log
mkprof/mosaic2.sh: prepconf.sh.log
mkprof/mosaic3.sh: prepconf.sh.log
mkprof/mosaic4.sh: prepconf.sh.log
mkprof/radeccat.sh: prepconf.sh.log
+ mkprof/3d-kernel.sh: prepconf.sh.log
mkprof/ellipticalmasks.sh: mknoise/addnoise.sh.log
mkprof/clearcanvas.sh: mknoise/addnoise.sh.log
endif
@@ -257,10 +261,10 @@ TESTS = prepconf.sh lib/multithread.sh $(MAYBE_CXX_TESTS)
\
# Files to distribute along with the tests.
EXTRA_DIST = $(TESTS) during-dev.sh buildprog/simpleio.c crop/cat.txt \
- match/positions-1.txt match/positions-2.txt mkprof/mkprofcat1.txt \
- mkprof/ellipticalmasks.txt mkprof/clearcanvas.txt mkprof/mkprofcat2.txt \
- mkprof/mkprofcat3.txt mkprof/mkprofcat4.txt mkprof/radeccat.txt \
- statistics/stdin-input.txt table/table.txt
+ mkprof/3d-cat.txt match/positions-1.txt match/positions-2.txt \
+ mkprof/mkprofcat1.txt mkprof/ellipticalmasks.txt mkprof/clearcanvas.txt \
+ mkprof/mkprofcat2.txt mkprof/mkprofcat3.txt mkprof/mkprofcat4.txt \
+ mkprof/radeccat.txt statistics/stdin-input.txt table/table.txt
diff --git a/tests/mkcatalog/simple-3d.sh b/tests/mkcatalog/simple-3d.sh
new file mode 100755
index 0000000..75e2986
--- /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 <address@hidden>
+# 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_detected_segmented.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
diff --git a/tests/mknoise/addnoise-3d.sh b/tests/mknoise/addnoise-3d.sh
new file mode 100755
index 0000000..ec010f4
--- /dev/null
+++ b/tests/mknoise/addnoise-3d.sh
@@ -0,0 +1,59 @@
+# Add noise to an input cube.
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <address@hidden>
+# 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).
+#
+# We will be adding noise to two images: the warped (smaller) and unwarped
+# (larger) mock images. The warped one will be used by programs that don't
+# care about the size of the image, but the larger one will be used by
+# those that do: for example SubtractSky and NoiseChisel will be better
+# tested on a larger image.
+prog=mknoise
+img=3d-cat.fits
+execname=../bin/$prog/ast$prog
+
+
+
+
+
+# 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 doesn't exist."; exit 77; fi
+if [ ! -f $img ]; then echo "$img does not exist."; exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+export GSL_RNG_SEED=1
+export GSL_RNG_TYPE=ranlxs2
+$execname --envseed $img
diff --git a/tests/mkprof/3d-cat.sh b/tests/mkprof/3d-cat.sh
new file mode 100755
index 0000000..459ce02
--- /dev/null
+++ b/tests/mkprof/3d-cat.sh
@@ -0,0 +1,48 @@
+# Create a 3D image with MakeProfiles
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <address@hidden>
+# 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=mkprof
+execname=../bin/$prog/ast$prog
+cat=$topsrc/tests/$prog/3d-cat.txt
+
+
+
+
+
+# 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).
+#
+# - Catalog doesn't exist (problem in tarball release).
+if [ ! -f $execname ]; then echo "$execname not created."; exit 77; fi
+if [ ! -f $cat ]; then echo "$cat does not exist."; exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+$execname $cat --config=.gnuastro/astmkprof-3d.conf --oversample=1
diff --git a/tests/mkprof/3d-cat.txt b/tests/mkprof/3d-cat.txt
new file mode 100644
index 0000000..ad54a01
--- /dev/null
+++ b/tests/mkprof/3d-cat.txt
@@ -0,0 +1,5 @@
+# Column 5: Profile [name, str8] Profile function.
+# x y z prof r n p1 p2 p3 q1 q2 m t
+1 10 5 1 sersic 20 2 45 0 0 0.5 0.25 -16 2
+2 50 50 160 sersic 10 4 20 30 40 1 0.5 -14 2
+3 80 50 80 sersic 30 1 80 50 20 0.4 0.8 -18 2
diff --git a/tests/mkprof/3d-kernel.sh b/tests/mkprof/3d-kernel.sh
new file mode 100755
index 0000000..d89f25e
--- /dev/null
+++ b/tests/mkprof/3d-kernel.sh
@@ -0,0 +1,47 @@
+# Create a 3D kernel with MakeProfiles.
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <address@hidden>
+# 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=mkprof
+execname=../bin/$prog/ast$prog
+
+
+
+
+
+# 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).
+#
+# - Catalog doesn't exist (problem in tarball release).
+if [ ! -f $execname ]; then echo "$execname not created."; exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+$execname $cat --kernel=gaussian-3d,2,3,0.5 --oversample=1 \
+ --output=3d-kernel.fits
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 63d062a: 3D features in Convolve, Match, MakeCatalog, MakeProfiles and library,
Mohammad Akhlaghi <=