[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 73fdf0c 006/113: MakeProfiles builds 3D ellips
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 73fdf0c 006/113: MakeProfiles builds 3D ellipsoids |
|
Date: |
Fri, 16 Apr 2021 10:33:30 -0400 (EDT) |
branch: master
commit 73fdf0cdfd4b4db02d5a842f0744693e27c37929
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>
MakeProfiles builds 3D ellipsoids
3D ellipsoids can now be built with MakeProfiles with any shape and
orientation. This is a major commit unifying all the work that was
necessary to accomplish this feature. In the process, several library
fuctions were also corrected and several bugs were also found and
fixed. Bugs relating to the building of Gnuastro that were reported from
outside of this task were implemented in the master branch until now also,
so there is no need to merge master with this branch. All the major
changes/additions have been described in the NEWS file.
---
NEWS | 21 +
THANKS | 2 +-
bin/mkprof/Makefile.am | 2 +-
bin/mkprof/args.h | 151 +++++---
bin/mkprof/astmkprof-3d.conf | 57 +++
bin/mkprof/main.h | 25 +-
bin/mkprof/mkprof.c | 358 +++++++++--------
bin/mkprof/mkprof.h | 24 +-
bin/mkprof/oneprofile.c | 401 +++++++++++++------
bin/mkprof/oneprofile.h | 2 +-
bin/mkprof/profiles.c | 52 ++-
bin/mkprof/profiles.h | 23 +-
bin/mkprof/ui.c | 379 ++++++++++++++----
bin/mkprof/ui.h | 3 +
bin/noisechisel/threshold.c | 7 +-
bin/warp/ui.c | 2 +-
doc/gnuastro.texi | 900 +++++++++++++++++++++++++------------------
lib/box.c | 194 ++++++++--
lib/data.c | 2 +-
lib/gnuastro/box.h | 7 +-
lib/gnuastro/tile.h | 127 ++++--
lib/tile.c | 18 +-
lib/wcs.c | 77 +++-
tests/Makefile.am | 15 +-
tests/mkprof/3dcat.sh | 48 +++
tests/prepconf.sh | 7 +-
26 files changed, 1946 insertions(+), 958 deletions(-)
diff --git a/NEWS b/NEWS
index 926dc3b..dd65fa6 100644
--- a/NEWS
+++ b/NEWS
@@ -10,6 +10,14 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
Arithmetic: now has a new `--globalhdu' (`-g') option which can be used
once for all the input images.
+ MakeProfiles: can build 3D ellipsoids of any orientation and shape. The
+ following options are new to MakeProfiles for enabling 3D profiles:
+ `--p2col', `--p3col' (second two Euler angles), and `--q2col' (axis ratio
+ of third axis to major axis). To help in using the 3D features, a
+ `astmkprof-3d.conf' configuration file will also be installed with
+ Gnuastro. The `--config' option can be used to specify default values for
+ 3D profiles.
+
MakeProfiles: the new `--kernel' option can make a kernel image without
the need to define a catalog. With this option, a catalog (or
accompanying background image) must not be given.
@@ -18,6 +26,11 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
to specify the PC matrix, CUNIT and CTYPE world coordinate system
keywords of the output FITS file.
+ MakeProfiles: can build a new raidal profile called `distance' (code
+ 7). The value this profile gives to each pixel is equal to its distance
+ from the center. It can be used to build any other radial profile that is
+ not supported by MakeProfiles.
+
** Removed features
** Changed features
@@ -80,6 +93,14 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
NoiseChisel segfault when no usable region for sky clumps (bug #51372).
+ Pixel scale measurement when dimension scale isn't equal or doesn't
+ decrease (bug #51385).
+
+ Improper types for function code in MakeProfiles (bug #51467).
+
+ Crashes on 32-bit and big-endian systems (bug #51476).
+
+
diff --git a/THANKS b/THANKS
index bed7f05..144060a 100644
--- a/THANKS
+++ b/THANKS
@@ -26,8 +26,8 @@ support in Gnuastro. The list is ordered alphabetically.
Brandon Invergo brandon@gnu.org
Lee Kelvin l.s.kelvin@ljmu.ac.uk
Mohammad-Reza Khellat moha.khe@gmail.com
- Guillaume Mahler guillaume.mahler@univ-lyon1.fr
Alan Lefor alefor@astr.tohoku.ac.jp
+ Guillaume Mahler guillaume.mahler@univ-lyon1.fr
Francesco Montanari francesco.montanari@openmailbox.org
William Pence William.Pence@nasa.gov
Yahya Sefidbakht y.sefidbakht@gmail.com
diff --git a/bin/mkprof/Makefile.am b/bin/mkprof/Makefile.am
index 81453aa..6806258 100644
--- a/bin/mkprof/Makefile.am
+++ b/bin/mkprof/Makefile.am
@@ -39,4 +39,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 aaea679..0c7441b 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -178,6 +178,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "envseed",
+ UI_KEY_ENVSEED,
+ 0,
+ 0,
+ "Use GSL_RNG_SEED environment variable for seed.",
+ ARGS_GROUP_PROFILES,
+ &p->envseed,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"tolerance",
UI_KEY_TOLERANCE,
"FLT",
@@ -204,6 +217,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "mforflatpix",
+ UI_KEY_MFORFLATPIX,
+ 0,
+ 0,
+ "mcol is flat pixel value (when fcol is 5 or 6)",
+ ARGS_GROUP_CATALOG,
+ &p->mforflatpix,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"shift",
UI_KEY_SHIFT,
"INT[, ...]",
@@ -244,6 +270,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "magatpeak",
+ UI_KEY_MAGATPEAK,
+ 0,
+ 0,
+ "Magnitude is for peak pixel, not full profile.",
+ ARGS_GROUP_PROFILES,
+ &p->magatpeak,
+ GAL_OPTIONS_NO_ARG_TYPE,
+ GAL_OPTIONS_RANGE_0_OR_1,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"circumwidth",
UI_KEY_CIRCUMWIDTH,
"FLT",
@@ -269,32 +308,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_MANDATORY,
GAL_OPTIONS_NOT_SET
},
- {
- "magatpeak",
- UI_KEY_MAGATPEAK,
- 0,
- 0,
- "Magnitude is for peak pixel, not full profile.",
- ARGS_GROUP_PROFILES,
- &p->magatpeak,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "envseed",
- UI_KEY_ENVSEED,
- 0,
- 0,
- "Use GSL_RNG_SEED environment variable for seed.",
- ARGS_GROUP_PROFILES,
- &p->envseed,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
@@ -306,24 +319,11 @@ struct argp_option program_options[] =
ARGS_GROUP_CATALOG
},
{
- "ccol",
- UI_KEY_CCOL,
- "STR/INT",
- 0,
- "Coordinate columns (one call for each dimension).",
- ARGS_GROUP_CATALOG,
- &p->ccol,
- GAL_TYPE_STRLL,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
"mode",
UI_KEY_MODE,
"STR",
0,
- "Coordinate mode `img' or `wcs'.",
+ "Mode of `--ccol': `img' or `wcs'.",
GAL_OPTIONS_GROUP_INPUT,
&p->mode,
GAL_TYPE_STRING,
@@ -333,12 +333,25 @@ struct argp_option program_options[] =
ui_parse_coordinate_mode
},
{
+ "ccol",
+ UI_KEY_CCOL,
+ "STR/INT",
+ 0,
+ "Coordinate columns (one call for each dimension).",
+ ARGS_GROUP_CATALOG,
+ &p->ccol,
+ GAL_TYPE_STRLL,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"fcol",
UI_KEY_FCOL,
"STR/INT",
0,
"sersic (1), moffat (2), gaussian (3), point (4), "
- "flat (5), circumference (6).",
+ "flat (5), circumference (6), distance (7).",
ARGS_GROUP_CATALOG,
&p->fcol,
GAL_TYPE_STRING,
@@ -377,7 +390,7 @@ struct argp_option program_options[] =
UI_KEY_PCOL,
"STR/INT",
0,
- "Position angle.",
+ "Position angle (First X-Z-X Euler angle in 3D).",
ARGS_GROUP_CATALOG,
&p->pcol,
GAL_TYPE_STRING,
@@ -386,11 +399,37 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "p2col",
+ UI_KEY_P2COL,
+ "STR/INT",
+ 0,
+ "Second Euler angle (X-Z-X order).",
+ ARGS_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).",
+ ARGS_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 radius).",
ARGS_GROUP_CATALOG,
&p->qcol,
GAL_TYPE_STRING,
@@ -399,6 +438,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "q2col",
+ UI_KEY_Q2COL,
+ "STR/INT",
+ 0,
+ "Axis ratio (major/dim3 radius).",
+ ARGS_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",
@@ -424,19 +476,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
},
- {
- "mforflatpix",
- UI_KEY_MFORFLATPIX,
- 0,
- 0,
- "mcol is flat pixel value (when fcol is 5 or 6)",
- ARGS_GROUP_CATALOG,
- &p->mforflatpix,
- GAL_OPTIONS_NO_ARG_TYPE,
- GAL_OPTIONS_RANGE_0_OR_1,
- GAL_OPTIONS_NOT_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
diff --git a/bin/mkprof/astmkprof-3d.conf b/bin/mkprof/astmkprof-3d.conf
new file mode 100644
index 0000000..0bf19e5
--- /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:
+ naxis 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,1
+ 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 041ddeb..9d4e989 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -65,6 +65,7 @@ enum profile_types
PROFILE_POINT, /* Point profile. */
PROFILE_FLAT, /* Flat profile. */
PROFILE_CIRCUMFERENCE, /* Circumference profile. */
+ PROFILE_DISTANCE, /* Elliptical radius of pixel. */
PROFILE_MAXIMUM_CODE, /* Just for a sanity check. */
};
@@ -91,13 +92,10 @@ struct builtqueue
size_t id; /* ID of this object. */
int ispsf; /* This is a PSF profile. */
int overlaps; /* ==1: Overlaps with the image. */
- float *img; /* Array of this profile's image. */
- size_t imgwidth; /* Width of *img. */
- long fpixel_i[2]; /* First pixel in output image. */
- long lpixel_i[2]; /* Last pixel in output image. */
- long fpixel_o[2]; /* First pixel in this array. */
+ gal_data_t *image; /* Array of this profile's image. */
+ gal_data_t *overlap_i; /* Overlap tile over individual array. */
+ gal_data_t *overlap_m; /* Overlap tile over merged array. */
int func; /* Profile's radial function. */
-
int indivcreated; /* ==1: an individual file is created. */
size_t numaccu; /* Number of accurate pixels. */
double accufrac; /* Difference of accurate values. */
@@ -139,8 +137,11 @@ 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). */
@@ -162,11 +163,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 131519a..dbdf4ce 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -33,6 +33,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/git.h>
#include <gnuastro/fits.h>
#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
#include <gnuastro/statistics.h>
#include <gnuastro-internal/timing.h>
@@ -79,7 +80,7 @@ builtqueue_addempty(struct builtqueue **bq)
__func__, sizeof *tbq);
/* Initialize some of the values. */
- tbq->img=NULL;
+ tbq->image=NULL;
tbq->numaccu=0;
tbq->accufrac=0.0f;
tbq->indivcreated=0;
@@ -118,15 +119,11 @@ saveindividual(struct mkonthread *mkp)
struct mkprofparams *p=mkp->p;
double *crpix;
- gal_data_t *data;
long os=p->oversample;
size_t i, ndim=p->out->ndim;
struct builtqueue *ibq=mkp->ibq;
char *filename, *jobname, *outdir=p->outdir;
- /* Note that `width' is in FITS format, not C. */
- size_t dsize[2]={mkp->width[1], mkp->width[0]};
-
/* Write the name and remove a similarly named file when the `--kernel'
option wasn't called. If `--kernel' is called, then we should just use
@@ -140,25 +137,20 @@ saveindividual(struct mkonthread *mkp)
}
- /* Put the array into a data structure */
- data=gal_data_alloc(ibq->img, GAL_TYPE_FLOAT32, 2, dsize, NULL, 0,
- p->cp.minmapsize, "MockImage", "Brightness", NULL);
-
-
/* Write the array to file (a separately built PSF doesn't need WCS
coordinates). */
if(ibq->ispsf && p->psfinimg==0)
- gal_fits_img_write(data, filename, NULL, PROGRAM_STRING);
+ gal_fits_img_write(ibq->image, filename, NULL, PROGRAM_STRING);
else
{
/* Allocate space for the corrected crpix and fill it in. Both
`crpix' and `fpixel_i' are in FITS order. */
crpix=gal_data_malloc_array(GAL_TYPE_FLOAT64, ndim, __func__, "crpix");
for(i=0;i<ndim;++i)
- crpix[0] = ((double *)(p->crpix->array))[0] - os*(mkp->fpixel_i[0]-1);
+ crpix[i] = ((double *)(p->crpix->array))[i] - os*(mkp->fpixel_i[i]-1);
/* Write the image. */
- gal_fits_img_write_corr_wcs_str(data, filename, p->wcsheader,
+ gal_fits_img_write_corr_wcs_str(ibq->image, filename, p->wcsheader,
p->wcsnkeyrec, crpix, NULL,
PROGRAM_STRING);
}
@@ -199,6 +191,148 @@ saveindividual(struct mkonthread *mkp)
/**************************************************************/
/************ The builders *************/
/**************************************************************/
+/* High-level function to built a single profile and prepare it for the
+ next steps. */
+static void
+mkprof_build_single(struct mkonthread *mkp, long *fpixel_i, long *lpixel_i,
+ long *fpixel_o)
+{
+ struct mkprofparams *p = mkp->p;
+ struct builtqueue *ibq = mkp->ibq;
+
+ void *ptr;
+ int needs_crop=0;
+ size_t i, ind, fits_i, ndim=p->out->ndim;
+ 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). */
+ gsl_rng_memcpy(mkp->rng, p->rng);
+
+ /* Set the seed of the random number generator if the
+ environment is not to be used. */
+ if(mkp->p->envseed==0)
+ gsl_rng_set(mkp->rng, gal_timing_time_based_rng_seed());
+
+ /* Make the profile */
+ oneprofile_make(mkp);
+
+ /* Build an individual image if necessary. */
+ if( p->individual || (ibq->ispsf && p->psfinimg==0))
+ {
+ saveindividual(mkp);
+ if(ibq->ispsf && p->psfinimg==0)
+ ibq->overlaps=0;
+ }
+
+ /* If we want a merged image, then a tile needs to be defined over the
+ individual profile array and the output merged array to define the
+ overlapping region. */
+ if(p->out->array)
+ {
+ /* Note that `fpixel_o' and `lpixel_o' were in the un-oversampled
+ image, they are also in the FITS coordinates. */
+ for(i=0;i<ndim;++i)
+ {
+ /* Set the start and width of the overlap. */
+ fits_i = ndim-i-1;
+ start_indiv[i] = os * (fpixel_o[fits_i] - 1);
+ start_mrg[i] = os * (fpixel_i[fits_i] - 1);
+ dsize[i] = os * (lpixel_i[fits_i] - fpixel_i[fits_i] + 1);
+
+ /* Check if we need to crop the individual image or not. */
+ if(dsize[i] != ibq->image->dsize[i]) needs_crop=1;
+ }
+
+
+ /* Define the individual overlap tile. */
+ if(needs_crop)
+ {
+ /* If a crop is needed, set the starting pointer. */
+ ind=gal_dimension_coord_to_index(ndim, ibq->image->dsize,
+ start_indiv);
+ ptr=gal_data_ptr_increment(ibq->image->array, ind,
+ ibq->image->type);
+ }
+ else ptr=ibq->image->array;
+ ibq->overlap_i=gal_data_alloc(ptr, ibq->image->type, ndim, dsize, NULL,
+ 0, -1, NULL, NULL, NULL);
+ ibq->overlap_i->block=ibq->image;
+
+
+ /* Define the merged overlap tile. */
+ ind=gal_dimension_coord_to_index(ndim, p->out->dsize, start_mrg);
+ ptr=gal_data_ptr_increment(p->out->array, ind, p->out->type);
+ ibq->overlap_m=gal_data_alloc(ptr, p->out->type, ndim, dsize, NULL,
+ 0, -1, NULL, NULL, NULL);
+ ibq->overlap_m->block=p->out;
+ }
+}
+
+
+
+
+
+/* The profile has been built, now add it to the queue of profiles that
+ must be written into the final merged image. */
+static void
+mkprof_add_built_to_write_queue(struct mkonthread *mkp,
+ struct builtqueue *ibq,
+ struct builtqueue **fbq, size_t counter)
+{
+ struct mkprofparams *p = mkp->p;
+
+ int lockresult;
+ pthread_mutex_t *qlock=&p->qlock;
+ pthread_cond_t *qready=&p->qready;
+
+ /* Try locking the mutex so no thread can change the value of p->bq. If
+ you can lock it, then put the internal builtqueue on top of the
+ general builtqueue. If you can't, continue adding to the internal
+ builtqueue (make the next profiles) until you find a chance to lock
+ the mutex. */
+ lockresult=pthread_mutex_trylock(qlock);
+ if(lockresult==0) /* Mutex was successfully locked. */
+ {
+ /* Add this internal queue to system queue. */
+ (*fbq)->next=p->bq;
+ p->bq=ibq;
+
+ /* If the list was empty when you locked the mutex, then either
+ `mkprof_write` is waiting behind a condition variable for you to
+ fill it up or not (either it hasn't got to setting the condition
+ variable yet (this function locked the mutex before
+ `mkprof_write`) or it just got the list to be made and is busy
+ writing the arrays in the output). In either case,
+ pthread_cond_signal will work. */
+ if((*fbq)->next==NULL)
+ pthread_cond_signal(qready);
+ pthread_mutex_unlock(qlock);
+
+ /* Finally set both the internal queue and the first internal queue
+ element to NULL.*/
+ (*fbq)=NULL;
+ mkp->ibq=NULL;
+ }
+
+ /* The mutex couldn't be locked and there are no more objects for this
+ thread to build (giving a chance for this thread to add up its built
+ profiles). So we have to lock the mutex to pass on this built
+ structure to the builtqueue. */
+ else if (mkp->indexs[counter+1]==GAL_BLANK_SIZE_T)
+ {
+ pthread_mutex_lock(qlock);
+ (*fbq)->next=p->bq;
+ p->bq=ibq;
+ pthread_cond_signal(qready);
+ pthread_mutex_unlock(qlock);
+ }
+}
+
+
+
+
+
/* Build the profiles that are indexed in the indexs array of the
mkonthread structure that was assigned to it.
@@ -227,19 +361,17 @@ saveindividual(struct mkonthread *mkp)
fractional value >=0.5, we will just shift the integer part of the
central pixel by one and completely ignore the fractional part.
*/
-void *
+static void *
mkprof_build(void *inparam)
{
struct mkonthread *mkp=(struct mkonthread *)inparam;
struct mkprofparams *p=mkp->p;
- size_t i, id;
- int lockresult;
- double center[2];
- long lpixel_o[2];
- pthread_mutex_t *qlock=&p->qlock;
+ size_t i, id, ndim=p->out->ndim;
struct builtqueue *ibq, *fbq=NULL;
- pthread_cond_t *qready=&p->qready;
+ 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. */
for(i=0; mkp->indexs[i]!=GAL_BLANK_SIZE_T; ++i)
@@ -264,94 +396,44 @@ mkprof_build(void *inparam)
if( p->f[id] == PROFILE_POINT )
mkp->width[0]=mkp->width[1]=1;
else
- gal_box_ellipse_in_box(mkp->truncr, mkp->q*mkp->truncr,
- p->p[id]*DEGREESTORADIANS, mkp->width);
-
+ {
+ if(p->out->ndim==2)
+ gal_box_bound_ellipse(mkp->truncr, mkp->q[0]*mkp->truncr,
+ p->p1[id], mkp->width);
+ else
+ {
+ 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);
+ }
+ }
/* Get the overlapping pixels using the starting points (NOT
oversampled). */
center[0]=p->x[id];
center[1]=p->y[id];
- gal_box_border_from_center(center, p->out->ndim, mkp->width,
- ibq->fpixel_i, ibq->lpixel_i);
- mkp->fpixel_i[0]=ibq->fpixel_i[0];
- mkp->fpixel_i[1]=ibq->fpixel_i[1];
- ibq->overlaps = gal_box_overlap(mkp->onaxes, ibq->fpixel_i,
- ibq->lpixel_i, ibq->fpixel_o,
- lpixel_o, 2);
+ 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);
+ ibq->overlaps = gal_box_overlap(mkp->onaxes, fpixel_i, lpixel_i,
+ fpixel_o, lpixel_o, ndim);
- /* Build the profile if necessary, After this, the width is
- oversampled. */
+ /* Build the profile if necessary. */
if(ibq->overlaps || p->individual || (ibq->ispsf && p->psfinimg==0))
- {
- /* Put a copy of the main random number generator for this
- thread to use for this profile. */
- gsl_rng_memcpy(mkp->rng, p->rng);
-
- /* Set the seed of the random number generator if the
- environment is not to be used. */
- if(mkp->p->envseed==0)
- gsl_rng_set(mkp->rng, gal_timing_time_based_rng_seed());
-
- /* Make the profile */
- oneprofile_make(mkp);
- if( p->individual || (ibq->ispsf && p->psfinimg==0))
- {
- saveindividual(mkp);
- if(ibq->ispsf && p->psfinimg==0)
- ibq->overlaps=0;
- }
- }
+ mkprof_build_single(mkp, fpixel_i, lpixel_i, fpixel_o);
- /* Add ibq to bq if you can lock the mutex. */
- if(p->cp.numthreads>1)
- {
- /* Try locking the mutex so no thread can change the value
- of p->bq. If you can lock it, then put the internal
- builtqueue on top of the general builtqueue. If you
- can't, continue adding to the internal builtqueue (make
- the next profiles) until you find a chance to lock the
- mutex. */
- lockresult=pthread_mutex_trylock(qlock);
- if(lockresult==0) /* Mutex was successfully locked. */
- {
- /* Add this internal queue to system queue. */
- fbq->next=p->bq;
- p->bq=ibq;
-
- /* If the list was empty when you locked the mutex, then
- either `mkprof_write` is waiting behind a condition
- variable for you to fill it up or not (either it hasn't
- got to setting the condition variable yet (this function
- locked the mutex before `mkprof_write`) or it just got the
- list to be made and is busy writing the arrays in the
- output). In either case, pthread_cond_signal will work. */
- if(fbq->next==NULL)
- pthread_cond_signal(qready);
- pthread_mutex_unlock(qlock);
-
- /* Finally set both the internal queue and the first
- internal queue element to NULL.*/
- fbq=NULL;
- mkp->ibq=NULL;
- }
- /* The mutex couldn't be locked and there are no more
- objects for this thread to build (giving a chance for
- this thread to add up its built profiles). So we have to
- lock the mutex to pass on this built structure to the
- builtqueue. */
- else if (mkp->indexs[i+1]==GAL_BLANK_SIZE_T)
- {
- pthread_mutex_lock(qlock);
- fbq->next=p->bq;
- p->bq=ibq;
- pthread_cond_signal(qready);
- pthread_mutex_unlock(qlock);
- }
- }
+ /* Add this profile to the list of profiles that must be written onto
+ the final merged image with another thread. */
+ if(p->cp.numthreads>1)
+ mkprof_add_built_to_write_queue(mkp, ibq, &fbq, i);
}
/* Free the allocated space for this thread and wait until all other
@@ -387,22 +469,16 @@ mkprof_build(void *inparam)
/**************************************************************/
/************ The writer *************/
/**************************************************************/
-void
+static void
mkprof_write(struct mkprofparams *p)
{
double sum;
char *jobname;
struct timeval t1;
- long os=p->oversample;
- int replace=p->replace;
gal_data_t *out=p->out, *log;
- size_t i, j, iw, jw, ii, jj, ow;
- size_t w=p->dsize[ out->ndim - 1]; /* Number of rows. */
struct builtqueue *ibq=NULL, *tbq;
- float *to, *from, *colend, *rowend;
size_t complete=0, num=p->num, clog;
-
/* Write each image into the output array. */
while(complete<p->num)
{
@@ -423,66 +499,17 @@ mkprof_write(struct mkprofparams *p)
}
sum=0.0f;
- /* Write the array pointed to by ibq into the output image. Note
- that the FITS and C arrays have opposite axis orders and FITS
- counting starts from 1, not zero. Also fpixel is the first
- (inclusive) pixel and so is lpixel (it is inclusive). */
+
+ /* During the build process, we also defined the overlap tiles of
+ both the individual array and the final merged array, here we will
+ use those to put the required profile pixels into the final
+ array. */
if(ibq->overlaps && out->array)
- {
+ GAL_TILE_PO_OISET(float,float,ibq->overlap_i,ibq->overlap_m,1,0, {
+ *o = p->replace ? ( *i==0.0f ? *o : *i ) : (*i + *o);
+ sum += *i;
+ });
- /* Set the starting and ending points in the complete image. */
- i = os * (ibq->fpixel_i[1]-1);
- j = os * (ibq->fpixel_i[0]-1);
-
-
- /* Set the starting and ending points in the overlapping
- image. Note that oversampling has already been taken
- into account in ibq->imgwidth. */
- ow = ibq->imgwidth;
- ii = os * (ibq->fpixel_o[1]-1);
- jj = os * (ibq->fpixel_o[0]-1);
-
-
- /* Find the width of the overlapping region: */
- iw = os*(ibq->lpixel_i[1]-ibq->fpixel_i[1]+1);
- jw = os*(ibq->lpixel_i[0]-ibq->fpixel_i[0]+1);
-
-
- /* Write the overlap to the actual image. Instead of writing
- two for loops and summing all the row and column indexs
- for every pixel and each image, we use pointer arithmetic
- which is much more efficient. Just think of one pointer
- that is advancing over the final image (*to) and one that
- is advancing over the overlap image (*from). Since we
- know the images overlap, iw and jw are both smaller than
- the two image number of columns and number of rows, so
- w-jw and ow-jw will always be positive. */
- to = (float *)(out->array) + i*w+j;
- from=ibq->img+ii*ow+jj;
- rowend=to+iw*w;
- do
- {
- /* Go over all the pixels in this row and write this profile
- into the final output array. Just note that when
- replacing, we don't want to replace those pixels that have
- a zero value, since no profile existed there. */
- colend=to+jw;
- do
- {
- *to = ( replace
- ? ( *from==0.0f ? *to : *from )
- : *to + *from );
- sum += *from;
- ++from;
- }
- while(++to<colend);
-
- /* Go to the next row. */
- to += w-jw;
- from += ow-jw;
- }
- while(to<rowend);
- }
/* Fill the log array. */
if(p->cp.log)
@@ -510,6 +537,7 @@ mkprof_write(struct mkprofparams *p)
}
}
+
/* Report if in verbose mode. */
++complete;
if(!p->cp.quiet && p->num>1)
@@ -520,10 +548,13 @@ mkprof_write(struct mkprofparams *p)
free(jobname);
}
+
/* Free the array and the queue element and change it to the next one
and increment complete. Note that there is no problem to free a
NULL pointer (when the built array didn't overlap). */
- free(ibq->img);
+ gal_data_free(ibq->overlap_i);
+ gal_data_free(ibq->overlap_m);
+ gal_data_free(ibq->image);
tbq=ibq->next;
free(ibq);
ibq=tbq;
@@ -594,7 +625,6 @@ mkprof(struct mkprofparams *p)
size_t i, fi, *indexs, thrdcols;
size_t nb, ndim=p->out->ndim, nt=p->cp.numthreads;
-
/* Allocate the arrays to keep the thread and parameters for each
thread. Note that we only want nt-1 threads to do the
building. */
diff --git a/bin/mkprof/mkprof.h b/bin/mkprof/mkprof.h
index 8264317..03151b2 100644
--- a/bin/mkprof/mkprof.h
+++ b/bin/mkprof/mkprof.h
@@ -31,26 +31,22 @@ struct mkonthread
{
/* General parameters: */
double r; /* Elliptical radius at this point. */
- double x; /* x value of this point. */
- double xl; /* lower x boundary */
- double xh; /* higher x boundary. */
- double y; /* y value when integrating over x. */
- double yl; /* lower y boundary. */
- double yh; /* higher y boundary. */
- double c; /* Cosine of the position angle. */
- double s; /* Sine of the position angle. */
- double q; /* axis ratio of the position angle. */
- double xc; /* Center in C of created(oversampled). */
- double yc; /* Center in C of created(oversampled). */
+ 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. */
- int 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. */
/* Random number generator: */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 6996d3b..60698b8 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -51,30 +51,124 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/****************************************************************
- ************** Elliptical radius ******************
+ ************** Radial distance ******************
****************************************************************/
-/* Convert cartesian coordinates to the rotated elliptical radius. */
-void
-r_el(struct mkonthread *mkp)
+/* Set the center position of the profile in the oversampled image. Note
+ that `mkp->width' is in the non-oversampled scale. IMPORTANT: the
+ ordering is in FITS coordinate order. */
+static void
+oneprofile_center_oversampled(struct mkonthread *mkp)
{
- double c=mkp->c, s=mkp->s, q=mkp->q, x=mkp->x, y=mkp->y;
- mkp->r = sqrt( (x*c+y*s)*(x*c+y*s) + ((y*c-x*s)*(y*c-x*s)/q/q) );
+ struct mkprofparams *p=mkp->p;
+
+ double *dim;
+ long os=p->oversample;
+ double val, pixfrac, intpart;
+ size_t i, id=mkp->ibq->id, ndim=p->out->ndim;
+
+ for(i=0;i<ndim;++i)
+ {
+ 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) );
+ mkp->center[i] = round(val*100)/100;
+ }
}
-/* Calculate the cercular distance of a pixel to the profile center. */
-float
-r_circle(size_t p, struct mkonthread *mkp)
+static void
+oneprofile_set_coord(struct mkonthread *mkp, size_t index)
{
- double x, y;
+ size_t i, coord_c[3];
+ uint8_t os=mkp->p->oversample;
+ size_t ndim=mkp->ibq->image->ndim, *dsize=mkp->ibq->image->dsize;
+
+ /* Get the coordinates in C order. */
+ gal_dimension_index_to_coord(index, ndim, dsize, coord_c);
+
+ /* Convert these coordinates to one where the profile center is at the
+ center and the image is not over-sampled. Note that only `coord_c' is
+ in C order.*/
+ for(i=0;i<ndim;++i)
+ mkp->coord[i] = ( coord_c[ndim-i-1] - mkp->center[i] )/os;
+}
+
+
+
- x = p/mkp->width[0]; /* Note that width[0] is the First FITS */
- y = p%mkp->width[0]; /* axis, not first C axis. */
- return sqrt( (x-mkp->xc)*(x-mkp->xc) + (y-mkp->yc)*(y-mkp->yc) );
+/* Convert cartesian coordinates to the rotated elliptical radius. See the
+ "Defining an ellipse" section of the book for the full derivation. */
+static void
+oneprofile_r_el(struct mkonthread *mkp)
+{
+ /* ct: cos(theta) st: sin(theta)
+ cp: cos(phi) sp: sin(phi) */
+ 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->out->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->out->ndim'", __func__, PACKAGE_BUGREPORT,
+ mkp->p->out->ndim);
+ }
+}
+
+
+
+
+
+/* Calculate the circular/spherical distance of a pixel to the profile
+ center. This is just used to add pixels in the stack. Later, when the
+ pixels are popped from the stack, the elliptical radius will be used to
+ give them a value.*/
+static float
+oneprofile_r_circle(size_t index, struct mkonthread *mkp)
+{
+ size_t i, c[3];
+ double d, sum=0.0f;
+ size_t ndim=mkp->ibq->image->ndim, *dsize=mkp->ibq->image->dsize;
+
+ /* Convert the index into a coordinate. */
+ gal_dimension_index_to_coord(index, ndim, dsize, c);
+
+ /* Find the distance to the center along each dimension (in FITS
+ order). */
+ for(i=0;i<ndim;++i)
+ {
+ d = c[ndim-i-1] - mkp->center[i];
+ sum += d*d;
+ }
+
+ /* Return the distance. */
+ return sqrt(sum);
}
@@ -100,21 +194,21 @@ r_circle(size_t p, struct mkonthread *mkp)
****************************************************************/
/* Fill pixel with random values */
float
-randompoints(struct mkonthread *mkp)
+oneprofile_randompoints(struct mkonthread *mkp)
{
- double xrange, yrange, sum=0.0f;
- size_t i, numrandom=mkp->p->numrandom;
+ double range[3], sum=0.0f;
+ size_t i, j, numrandom=mkp->p->numrandom, ndim=mkp->p->out->ndim;
- /* Set the range of the x and y: */
- xrange=mkp->xh-mkp->xl;
- yrange=mkp->yh-mkp->yl;
+ /* Set the range in each dimension. */
+ for(i=0;i<ndim;++i)
+ range[i] = mkp->higher[i] - mkp->lower[i];
- /* Find the sum of the profile on the random positions */
+ /* Find the sum of the profile on the random positions. */
for(i=0;i<numrandom;++i)
{
- mkp->x = mkp->xl + gsl_rng_uniform(mkp->rng)*xrange;
- mkp->y = mkp->yl + gsl_rng_uniform(mkp->rng)*yrange;
- r_el(mkp);
+ for(j=0;j<ndim;++j)
+ mkp->coord[j] = mkp->lower[j] + gsl_rng_uniform(mkp->rng) * range[j];
+ oneprofile_r_el(mkp);
sum+=mkp->profile(mkp);
}
@@ -143,15 +237,16 @@ randompoints(struct mkonthread *mkp)
/****************************************************************
***************** 2D integration ********************
****************************************************************/
-/* This function is used in the integration of a profile. It
- assumes a fixed y and integrates over a range of x values. */
+/* This is an old implementation which we are not using now. But it is kept
+ here in case it can be useful */
+#if 0
double
twod_over_x(double x, void *params)
{
struct mkonthread *mkp=(struct mkonthread *) params;
mkp->x=x;
- r_el(mkp);
+ oneprofile_r_el(mkp);
return mkp->profile(mkp);
}
@@ -196,8 +291,7 @@ integ2d(struct mkonthread *mkp)
epsrel, &result, &abserr, &neval);
return result;
}
-
-
+#endif
@@ -220,49 +314,76 @@ integ2d(struct mkonthread *mkp)
/************ Pixel by pixel building *************/
/********* Positions are in C not FITS *********/
/**************************************************************/
+/* `oneprofile_center_oversampled' stored the center of the profile in
+ floating point coordinates. This function will convert that into a
+ pixel index. */
+static size_t
+oneprofile_center_pix_index(struct mkonthread *mkp)
+{
+ double pixfrac, intpart;
+ size_t *dsize=mkp->ibq->image->dsize;
+ size_t i, coord[3], ndim=mkp->ibq->image->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
+ `gal_dimension_coord_to_index'). */
+ for(i=0;i<ndim;++i)
+ {
+ pixfrac = modf(mkp->center[i], &intpart);
+ coord[ndim-i-1] = (long)(mkp->center[i]) + ( pixfrac<0.5f ? 0 : 1 );
+ }
+
+ /* Retun the pixel index of this coordinate. */
+ return gal_dimension_coord_to_index(ndim, dsize, coord);
+}
+
+
+
+
+
static void
-makepixbypix(struct mkonthread *mkp)
+oneprofile_pix_by_pix(struct mkonthread *mkp)
{
- size_t ndim=2, dsize[2]={mkp->width[1], mkp->width[0]};
+ struct builtqueue *ibq=mkp->ibq;
+ size_t ndim=ibq->image->ndim, *dsize=ibq->image->dsize;
uint8_t *byt;
gal_list_sizet_t *Q=NULL;
int use_rand_points=1, ispeak=1;
- struct builtqueue *ibq=mkp->ibq;
- float circ_r, *img=mkp->ibq->img;
- size_t *dinc=gal_dimension_increment(ndim, dsize);
- double tolerance=mkp->p->tolerance, pixfrac, junk;
+ double tolerance=mkp->p->tolerance;
+ float circ_r, *array=mkp->ibq->image->array;
double (*profile)(struct mkonthread *)=mkp->profile;
- double xc=mkp->xc, yc=mkp->yc, os=mkp->p->oversample;
- size_t p, x, y, is1=mkp->width[0], is0=mkp->width[1];
double truncr=mkp->truncr, approx, hp=0.5f/mkp->p->oversample;
+ size_t i, p, *dinc=gal_dimension_increment(ndim, dsize);
/* lQ: Largest. sQ: Smallest in queue */
gal_list_dosizet_t *lQ=NULL, *sQ;
/* Find the nearest pixel to the profile center and add it to the
queue. */
- pixfrac = modf(mkp->xc, &junk);
- x=(long)mkp->xc + ( pixfrac<0.5f ? 0 : 1 );
- pixfrac = modf(mkp->yc, &junk);
- y=(long)mkp->yc + ( pixfrac<0.5f ? 0 : 1 );
- p=x*mkp->width[0]+y;
+ p=oneprofile_center_pix_index(mkp);
- /* If this is a point source, just fill that one pixel and go. */
+ /* If this is a point source, just fill that one pixel and leave this
+ function. */
if(mkp->func==PROFILE_POINT)
- { img[p]=1; return; }
+ { array[p]=1; return; }
- /* Allocate the byt array to not repeat completed pixels. */
- byt = gal_data_calloc_array(GAL_TYPE_UINT8, is0*is1, __func__, "byt");
+ /* Allocate the `byt' array. It is used as a flag to make sure that we
+ don't re-calculate the profile value on a pixel more than once. */
+ byt = gal_data_calloc_array(GAL_TYPE_UINT8,
+ gal_dimension_total_size(ndim, dsize),
+ __func__, "byt");
/* Start the queue: */
byt[p]=1;
- gal_list_dosizet_add( &lQ, &sQ, p, r_circle(p, mkp) );
+ gal_list_dosizet_add( &lQ, &sQ, p, oneprofile_r_circle(p, mkp) );
/* If random points are necessary, then do it: */
- if(mkp->func==PROFILE_SERSIC || mkp->func==PROFILE_MOFFAT
- || mkp->func==PROFILE_GAUSSIAN)
+ switch(mkp->func)
{
+ case PROFILE_SERSIC:
+ case PROFILE_MOFFAT:
+ case PROFILE_GAUSSIAN:
while(sQ)
{
/* In case you want to see the status of the twosided ordered
@@ -270,39 +391,34 @@ makepixbypix(struct mkonthread *mkp)
line. Note that there will be a lot of lines printed! */
/*print_tossll(lQ, sQ);*/
- /* Pop the pixel from the queue and check if it is within the
- truncation radius. Note that `xc` and `p` both belong to the
- over sampled image. But all the profile parameters are in the
- non-oversampled image. So we divide the distance by os
- (p->oversample in double type) */
+ /* Pop a pixel from the queue, convert its index into coordinates
+ and use them to estimate the elliptical radius of the
+ pixel. If the pixel is outside the truncation radius, ignore
+ it. */
p=gal_list_dosizet_pop_smallest(&lQ, &sQ, &circ_r);
- mkp->x=(p/is1-xc)/os;
- mkp->y=(p%is1-yc)/os;
- r_el(mkp);
- if(mkp->r>truncr) continue;
+ oneprofile_set_coord(mkp, p);
+ oneprofile_r_el(mkp);
+ if(mkp->r > truncr) continue;
/* Find the value for this pixel: */
- mkp->xl=mkp->x-hp;
- mkp->xh=mkp->x+hp;
- mkp->yl=mkp->y-hp;
- mkp->yh=mkp->y+hp;
- /*
- printf("Center (%zu, %zu). r: %.4f. x: [%.4f--%.4f], "
- "y: [%.4f, %.4f]\n", p%is1+1, p/is1+1, mkp->r, mkp->xl,
- mkp->xh, mkp->yl, mkp->yh);
- */
+ for(i=0;i<ndim;++i)
+ {
+ mkp->lower[i] = mkp->coord[i] - hp;
+ mkp->higher[i] = mkp->coord[i] + hp;
+ }
+
/* Find the random points and profile center. */
- img[p]=randompoints(mkp);
+ array[p]=oneprofile_randompoints(mkp);
approx=profile(mkp);
- if (fabs(img[p]-approx)/img[p] < tolerance)
+ if (fabs(array[p]-approx)/array[p] < tolerance)
use_rand_points=0;
/* Save the peak flux if this is the first pixel: */
- if(ispeak) { mkp->peakflux=img[p]; ispeak=0; }
+ if(ispeak) { mkp->peakflux=array[p]; ispeak=0; }
/* For the log file: */
++ibq->numaccu;
- ibq->accufrac+=img[p];
+ ibq->accufrac+=array[p];
/* Go over the neighbors and add them to queue of elements to
check if they haven't been done already. */
@@ -311,7 +427,8 @@ makepixbypix(struct mkonthread *mkp)
if(byt[nind]==0)
{
byt[nind]=1;
- gal_list_dosizet_add( &lQ, &sQ, nind, r_circle(nind, mkp) );
+ gal_list_dosizet_add( &lQ, &sQ, nind,
+ oneprofile_r_circle(nind, mkp) );
}
} );
@@ -329,23 +446,26 @@ makepixbypix(struct mkonthread *mkp)
while(Q)
{
p=gal_list_sizet_pop(&Q);
- mkp->x=(p/is1-xc)/os;
- mkp->y=(p%is1-yc)/os;
- r_el(mkp);
+ oneprofile_set_coord(mkp, p);
+ oneprofile_r_el(mkp);
+
+ /* See if this pixel's radial distance is larger than the truncation
+ radius. If so, then don't add its neighbors to the queue and
+ continue to the next pixel in the queue. */
if(mkp->r>truncr)
{
/* For the circumference, if the profile is too elongated
and circumwidth is too small, then some parts of the
circumference will not be shown without this condition. */
- if(mkp->func==PROFILE_CIRCUMFERENCE) img[p]=profile(mkp);
+ if(mkp->func==PROFILE_CIRCUMFERENCE) array[p]=profile(mkp);
continue;
}
/* Find the value for this pixel: */
- img[p]=profile(mkp);
+ array[p]=profile(mkp);
/* Save the peak flux if this is the first pixel: */
- if(ispeak) { mkp->peakflux=img[p]; ispeak=0; }
+ if(ispeak) { mkp->peakflux=array[p]; ispeak=0; }
/* Go over the neighbours and add them to queue of elements
to check. */
@@ -386,7 +506,7 @@ makepixbypix(struct mkonthread *mkp)
/************ Set profile parameters *************/
/**************************************************************/
int
-oneprofile_ispsf(int fcode)
+oneprofile_ispsf(uint8_t fcode)
{
return fcode==PROFILE_MOFFAT || fcode==PROFILE_GAUSSIAN;
}
@@ -403,28 +523,59 @@ oneprof_set_prof_params(struct mkonthread *mkp)
double sigma;
int tp=p->tunitinp;
- size_t id=mkp->ibq->id;
-
- /* Fill in the profile independant parameters. */
- p->x[id] += p->shift[1]/p->oversample; /* Shifts were multiplied by */
- p->y[id] += p->shift[0]/p->oversample; /* `p->oversample' before. */
- mkp->c = cos( (90-p->p[id]) * DEGREESTORADIANS );
- mkp->s = sin( (90-p->p[id]) * DEGREESTORADIANS );
- mkp->q = p->q[id];
+ size_t id=mkp->ibq->id, ndim=p->out->ndim;
+
+ /* Fill the most basic dimension and profile agnostic parameters. */
mkp->brightness = pow( 10, (p->zeropoint - p->m[id]) / 2.5f );
mkp->ibq->ispsf = p->kernel ? 1 : oneprofile_ispsf(p->f[id]);
mkp->func = mkp->ibq->func = p->f[id];
- /* Fill the profile dependent parameters. */
+ /* 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 "
+ "`p->out->ndim'", __func__, PACKAGE_BUGREPORT, p->out->ndim);
+ }
+
+
+ /* Fill the profile-dependent parameters. */
switch (mkp->func)
{
case PROFILE_SERSIC:
mkp->correction = 1;
- mkp->profile = &Sersic;
+ mkp->profile = &profiles_sersic;
mkp->sersic_re = p->r[id];
mkp->sersic_inv_n = 1.0f/p->n[id];
- mkp->sersic_nb = -1.0f*sersic_b(p->n[id]);
+ mkp->sersic_nb = -1.0f*profiles_sersic_b(p->n[id]);
mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
break;
@@ -432,9 +583,9 @@ oneprof_set_prof_params(struct mkonthread *mkp)
case PROFILE_MOFFAT:
mkp->correction = 1;
- mkp->profile = &Moffat;
+ mkp->profile = &profiles_moffat;
mkp->moffat_nb = -1.0f*p->n[id];
- mkp->moffat_alphasq = moffat_alpha(p->r[id], p->n[id]);
+ mkp->moffat_alphasq = profiles_moffat_alpha(p->r[id], p->n[id]);
mkp->moffat_alphasq *= mkp->moffat_alphasq;
mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id]/2;
if(p->psfinimg==0 && p->individual==0)
@@ -442,6 +593,8 @@ oneprof_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;
@@ -449,7 +602,7 @@ oneprof_set_prof_params(struct mkonthread *mkp)
case PROFILE_GAUSSIAN:
mkp->correction = 1;
- mkp->profile = &Gaussian;
+ mkp->profile = &profiles_gaussian;
sigma = p->r[id]/2.35482f;
mkp->gaussian_c = -1.0f/(2.0f*sigma*sigma);
mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id]/2;
@@ -458,6 +611,8 @@ oneprof_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;
@@ -466,13 +621,13 @@ oneprof_set_prof_params(struct mkonthread *mkp)
case PROFILE_POINT:
mkp->correction = 1;
mkp->fixedvalue = 1.0f;
- mkp->profile = &Flat;
+ mkp->profile = &profiles_flat;
break;
case PROFILE_FLAT:
- mkp->profile = &Flat;
+ mkp->profile = &profiles_flat;
mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
if(p->mforflatpix)
{
@@ -489,7 +644,7 @@ oneprof_set_prof_params(struct mkonthread *mkp)
case PROFILE_CIRCUMFERENCE:
- mkp->profile = &Circumference;
+ mkp->profile = &profiles_circumference;
mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
mkp->intruncr = mkp->truncr - p->circumwidth;
if(p->mforflatpix)
@@ -508,9 +663,15 @@ oneprof_set_prof_params(struct mkonthread *mkp)
+ case PROFILE_DISTANCE:
+ mkp->profile = profiles_radial_distance;
+ mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
+ mkp->correction = 0;
+ break;
+
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us so we can correct "
- "this problem. The profile code %d is not recognized.", __func__,
+ "this problem. The profile code %u is not recognized.", __func__,
mkp->func);
}
}
@@ -542,44 +703,33 @@ oneprofile_make(struct mkonthread *mkp)
{
struct mkprofparams *p=mkp->p;
+ double sum;
float *f, *ff;
- long os=p->oversample;
- double sum, pixfrac, intpart;
- size_t size, id=mkp->ibq->id;
-
+ size_t i, dsize[3], ndim=p->out->ndim;
- /* Find the profile center (see comments above
- `mkprof_build'). mkp->width is in the non-oversampled scale.*/
- pixfrac = modf(fabs(p->x[id]), &intpart);
- mkp->yc = ( os * (mkp->width[0]/2 + pixfrac)
- + (pixfrac<0.50f ? os/2 : -1*os/2-1) );
- mkp->yc = round(mkp->yc*100)/100;
- pixfrac = modf(fabs(p->y[id]), &intpart);
- mkp->xc = ( os*(mkp->width[1]/2 + pixfrac)
- + (pixfrac<0.5f ? os/2 : -1*os/2-1) );
- mkp->xc = round(mkp->xc*100)/100;
+ /* Find the profile center in the over-sampled image in C
+ coordinates. IMPORTANT: width must not be oversampled.*/
+ oneprofile_center_oversampled(mkp);
- /* From this point on, the widths are the actual pixel
- widths (with oversampling). */
- mkp->width[0] *= os;
- mkp->width[1] *= os;
- mkp->ibq->imgwidth=mkp->width[0];
+ /* From this point on, the widths will be in the actual pixel widths
+ (with oversampling). */
+ for(i=0;i<ndim;++i)
+ {
+ mkp->width[i] *= p->oversample;
+ dsize[ndim-i-1] = mkp->width[i];
+ }
/* Allocate and clear the array for this one profile. */
- errno=0;
- size=mkp->width[0]*mkp->width[1];
- mkp->ibq->img=calloc(size, sizeof *mkp->ibq->img);
- if(mkp->ibq->img==NULL)
- error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for object in row %zu of "
- "data in %s", __func__,
- size*sizeof *mkp->ibq->img, mkp->ibq->id, p->catname);
+ mkp->ibq->image=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, ndim, dsize,
+ NULL, 1, p->cp.minmapsize, "MOCK",
+ "Brightness", NULL);
/* Build the profile in the image. */
- makepixbypix(mkp);
+ oneprofile_pix_by_pix(mkp);
/* Correct the sum of pixels in the profile so it has the fixed total
@@ -589,7 +739,8 @@ oneprofile_make(struct mkonthread *mkp)
if(mkp->correction)
{
/* First get the sum of all the pixels in the profile. */
- sum=0.0f; ff=(f=mkp->ibq->img)+size; do sum+=*f++; while(f<ff);
+ ff=(f=mkp->ibq->image->array) + mkp->ibq->image->size;
+ sum=0.0f; do sum+=*f++; while(f<ff);
/* Correct the fraction of brightness that was calculated
accurately (not using the pixel center). */
@@ -604,7 +755,7 @@ oneprofile_make(struct mkonthread *mkp)
value, then the whole profile's box is going to be NaN values,
which is inconvenient and with the simple check here we can avoid
it (only have the profile's pixels set to NaN. */
- ff=(f=mkp->ibq->img)+size;
+ ff = (f=mkp->ibq->image->array) + mkp->ibq->image->size;
if(isnan(mkp->brightness))
do *f = *f ? NAN : *f ; while(++f<ff);
else
diff --git a/bin/mkprof/oneprofile.h b/bin/mkprof/oneprofile.h
index e04f85b..08cc19c 100644
--- a/bin/mkprof/oneprofile.h
+++ b/bin/mkprof/oneprofile.h
@@ -26,7 +26,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "mkprof.h"
int
-oneprofile_ispsf(int fcolvalue);
+oneprofile_ispsf(uint8_t fcolvalue);
void
oneprof_set_prof_params(struct mkonthread *mkp);
diff --git a/bin/mkprof/profiles.c b/bin/mkprof/profiles.c
index 6ee0197..b869c1e 100644
--- a/bin/mkprof/profiles.c
+++ b/bin/mkprof/profiles.c
@@ -43,18 +43,19 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
****************************************************************/
/* The Gaussian function at a point. */
double
-Gaussian(struct mkonthread *mkp)
+profiles_radial_distance(struct mkonthread *mkp)
{
- return exp( mkp->gaussian_c * mkp->r * mkp->r );
+ return mkp->r;
}
+
/* The integral of the Gaussian from -inf to +inf equals the square
root of PI. So from zero to +inf it equals half of that.*/
double
-totgaussian(double q)
+profiles_gaussian_total(double q)
{
return q*sqrt(M_PI)/2;
}
@@ -63,6 +64,17 @@ totgaussian(double q)
+/* The Gaussian function at a point. */
+double
+profiles_gaussian(struct mkonthread *mkp)
+{
+ return exp( mkp->gaussian_c * mkp->r * mkp->r );
+}
+
+
+
+
+
/* This function will find the moffat function alpha value based on
the explantions here:
@@ -71,7 +83,7 @@ totgaussian(double q)
alpha=(FWHM/2)/(2^(1.0/beta)-1)^(0.5). Then the moffat
function at r is: (1.0 + (r/alpha)^2.0)^(-1.0*beta)*/
double
-moffat_alpha(double fwhm, double beta)
+profiles_moffat_alpha(double fwhm, double beta)
{
return (fwhm/2)/pow((pow(2, 1/beta)-1), 0.5f);
}
@@ -84,7 +96,7 @@ moffat_alpha(double fwhm, double beta)
from Pengetal 2010 (Galfit). In finding the profiles, I am assuming
\Sigma_0=1. So that is what I put here too.*/
double
-totmoffat(double alpha, double beta, double q)
+profiles_moffat_total(double alpha, double beta, double q)
{
return M_PI*alpha*alpha*q/(beta-1);
}
@@ -99,7 +111,7 @@ totmoffat(double alpha, double beta, double q)
This is done before hand to speed up the process. */
double
-Moffat(struct mkonthread *mkp)
+profiles_moffat(struct mkonthread *mkp)
{
return pow(1+mkp->r*mkp->r/mkp->moffat_alphasq, mkp->moffat_nb);
}
@@ -112,7 +124,7 @@ Moffat(struct mkonthread *mkp)
Courteau and Holtzman 2003:
http://adsabs.harvard.edu/abs/2003ApJ...582..689 */
double
-sersic_b(double n)
+profiles_sersic_b(double n)
{
if(n<=0.35f)
error(EXIT_FAILURE, 0, "the Sersic index cannot be smaller "
@@ -125,27 +137,27 @@ sersic_b(double n)
-/* Find the Sersic profile for a certain radius. rdre=r/re, inv_n=1/n,
- nb= -1*b. */
+/* Find the total brightness in a Sersic profile. From equation 4 in
+ Peng 2010. This assumes the surface brightness at the effective
+ radius is 1.*/
double
-Sersic(struct mkonthread *mkp)
+profiles_sersic_total(double n, double re, double b, double q)
{
- return exp( mkp->sersic_nb
- * ( pow(mkp->r/mkp->sersic_re, mkp->sersic_inv_n) -1 ) );
+ return (2*M_PI*re*re*exp(b)*n*pow(b, -2*n)*q*
+ gsl_sf_gamma(2*n));
}
-/* Find the total brightness in a Sersic profile. From equation 4 in
- Peng 2010. This assumes the surface brightness at the effective
- radius is 1.*/
+/* Find the Sersic profile for a certain radius. rdre=r/re, inv_n=1/n,
+ nb= -1*b. */
double
-totsersic(double n, double re, double b, double q)
+profiles_sersic(struct mkonthread *mkp)
{
- return (2*M_PI*re*re*exp(b)*n*pow(b, -2*n)*q*
- gsl_sf_gamma(2*n));
+ return exp( mkp->sersic_nb
+ * ( pow(mkp->r/mkp->sersic_re, mkp->sersic_inv_n) -1 ) );
}
@@ -154,7 +166,7 @@ totsersic(double n, double re, double b, double q)
/* Make a circumference (inner to the radius). */
double
-Circumference(struct mkonthread *mkp)
+profiles_circumference(struct mkonthread *mkp)
{
return mkp->r > mkp->intruncr ? mkp->fixedvalue : 0.0f;
}
@@ -165,7 +177,7 @@ Circumference(struct mkonthread *mkp)
/* Always returns a fixed value: */
double
-Flat(struct mkonthread *mkp)
+profiles_flat(struct mkonthread *mkp)
{
return mkp->fixedvalue;
}
diff --git a/bin/mkprof/profiles.h b/bin/mkprof/profiles.h
index 3f8cfbd..0efa236 100644
--- a/bin/mkprof/profiles.h
+++ b/bin/mkprof/profiles.h
@@ -24,33 +24,36 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#define PROFILES_H
double
-Gaussian(struct mkonthread *mkp);
+profiles_radial_distance(struct mkonthread *mkp);
double
-totgaussian(double q);
+profiles_gaussian_total(double q);
double
-Moffat(struct mkonthread *mkp);
+profiles_gaussian(struct mkonthread *mkp);
double
-moffat_alpha(double fwhm, double beta);
+profiles_moffat_alpha(double fwhm, double beta);
double
-totmoffat(double alpha, double beta, double q);
+profiles_moffat_total(double alpha, double beta, double q);
double
-Sersic(struct mkonthread *mkp);
+profiles_moffat(struct mkonthread *mkp);
double
-sersic_b(double n);
+profiles_sersic_b(double n);
double
-totsersic(double n, double re, double b, double q);
+profiles_sersic_total(double n, double re, double b, double q);
double
-Circumference(struct mkonthread *mkp);
+profiles_sersic(struct mkonthread *mkp);
double
-Flat(struct mkonthread *mkp);
+profiles_circumference(struct mkonthread *mkp);
+
+double
+profiles_flat(struct mkonthread *mkp);
#endif
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index ecd1cd4..2b9d2ec 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -109,7 +109,7 @@ enum program_args_groups
/**************************************************************/
/********* Initialize & Parse command-line **************/
/**************************************************************/
-static int
+static uint8_t
ui_profile_name_read(char *string, size_t row)
{
if( !strcmp("sersic", string) )
@@ -130,6 +130,9 @@ ui_profile_name_read(char *string, size_t row)
else if ( !strcmp("circum", string) )
return PROFILE_CIRCUMFERENCE;
+ else if ( !strcmp("distance", string) )
+ return PROFILE_DISTANCE;
+
else if ( !strcmp(GAL_BLANK_STRING, string) )
error(EXIT_FAILURE, 0, "atleast one profile function is blank");
@@ -453,14 +456,16 @@ ui_parse_coordinate_mode(struct argp_option *option, char
*arg,
/* We want to print the stored values. */
if(lineno==-1)
{
- gal_checkset_allocate_copy( *(int *)(option->value)==MKPROF_MODE_IMG
+ gal_checkset_allocate_copy( *(uint8_t *)(option->value)==MKPROF_MODE_IMG
? "img" : "wcs", &outstr );
return outstr;
}
else
{
- if (!strcmp(arg, "img")) *(int *)(option->value)=MKPROF_MODE_IMG;
- else if (!strcmp(arg, "wcs")) *(int *)(option->value)=MKPROF_MODE_WCS;
+ if(!strcmp(arg, "img"))
+ *(uint8_t *)(option->value)=MKPROF_MODE_IMG;
+ else if (!strcmp(arg, "wcs"))
+ *(uint8_t *)(option->value)=MKPROF_MODE_WCS;
else
error_at_line(EXIT_FAILURE, 0, filename, lineno, "`%s' (value to "
"`--mode') not recognized as a coordinate standard "
@@ -623,7 +628,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;
char *colname=NULL, **strarr;
@@ -631,16 +636,15 @@ ui_read_cols(struct mkprofparams *p)
gal_data_t *cols, *tmp, *corrtype=NULL;
size_t i, counter=0, coordcounter=0, ndim=p->out->ndim;
- /* Make sure the number of coordinate columns and number of dimensions in
- outputs are the same. There is no problem if it is more than
- `ndim'. In that case, the last values (possibly in configuration
- files) will be ignored.*/
- if( gal_list_str_number(p->ccol) < ndim )
- error(EXIT_FAILURE, 0, "%zu coordinate columns (calls to `--coordcol') "
- "given but output has %zu dimensions",
- gal_list_str_number(p->ccol), p->out->ndim);
+ /* The coordinate columns are a linked list of strings. */
+ ccol=p->ccol;
+ for(i=0;i<ndim;++i)
+ {
+ gal_list_str_add(&colstrs, ccol->v, 0);
+ ccol=ccol->next;
+ }
- /* Put the columns a specific order to read. */
+ /* Add the rest of the columns in a specific order. */
gal_list_str_add(&colstrs, p->fcol, 0);
gal_list_str_add(&colstrs, p->rcol, 0);
gal_list_str_add(&colstrs, p->ncol, 0);
@@ -649,15 +653,6 @@ ui_read_cols(struct mkprofparams *p)
gal_list_str_add(&colstrs, p->mcol, 0);
gal_list_str_add(&colstrs, p->tcol, 0);
- /* The coordinate columns might be different, So we'll add them to the
- end (top of the list). */
- ccol=p->ccol;
- for(i=0;i<ndim;++i)
- {
- gal_list_str_add(&colstrs, ccol->v, 0);
- ccol=ccol->next;
- }
-
/* Reverse the order to make the column orders correspond to how we added
them here and avoid possible bugs. */
gal_list_str_reverse(&colstrs);
@@ -669,8 +664,7 @@ ui_read_cols(struct mkprofparams *p)
/* Set the number of objects. */
p->num=cols->size;
- /* For a sanity check, make sure that the total number of columns read is
- the same as those that were wanted (it might be more). */
+ /* Put each column's data in the respective internal array. */
while(cols!=NULL)
{
/* Pop out the top column. */
@@ -686,9 +680,23 @@ ui_read_cols(struct mkprofparams *p)
switch(++counter)
{
case 1:
+ case 2:
+ ++coordcounter;
+ colname = ( counter==1
+ ? "first coordinate column (`--coordcol')"
+ : "second 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;
+ }
+ break;
+
+ case 3:
if(tmp->type==GAL_TYPE_STRING)
{
- p->f=gal_data_malloc_array(GAL_TYPE_INT32, p->num,
+ p->f=gal_data_malloc_array(GAL_TYPE_UINT8, p->num,
__func__, "p->f");
strarr=tmp->array;
for(i=0;i<p->num;++i)
@@ -700,14 +708,14 @@ ui_read_cols(struct mkprofparams *p)
{
/* Read the user's profile codes. */
colname="profile function code (`fcol')";
- corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_INT32);
+ 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: table row %zu, the function "
- "code is %d. It should be >%d and <%d. Please run "
+ 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 "
@@ -719,57 +727,254 @@ ui_read_cols(struct mkprofparams *p)
}
break;
- case 2:
+ case 4:
colname="radius (`rcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->r=corrtype->array;
+
+ /* Check if there is no zero-radius profile. */
+ for(i=0;i<p->num;++i)
+ if(p->r[i]<=0.0f)
+ error(EXIT_FAILURE, 0, "%s: row %zu, the radius value %g is "
+ "not acceptable. It has to be larger than 0", p->catname,
+ i+1, p->r[i]);
break;
- case 3:
+ case 5:
colname="index (`ncol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->n=corrtype->array;
break;
- case 4:
+ 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 5:
+ 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;
break;
- case 6:
+ case 8:
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 7:
+ case 9:
colname="truncation (`tcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->t=corrtype->array;
break;
- case 8:
- case 9:
+ /* 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);
+ }
+ }
+}
+
+
+
+
+
+/* Read the columns for a 3D profile. */
+static void
+ui_read_cols_3d(struct mkprofparams *p)
+{
+ int checkblank;
+ char *colname=NULL, **strarr;
+ gal_list_str_t *colstrs=NULL, *ccol;
+ gal_data_t *cols, *tmp, *corrtype=NULL;
+ size_t i, counter=0, coordcounter=0, ndim=p->out->ndim;
+
+ /* 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<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. */
+ cols=gal_table_read(p->catname, p->cp.hdu, colstrs, p->cp.searchin,
+ p->cp.ignorecase, p->cp.minmapsize);
+
+ /* 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;
+
+ /* Note that the input was a linked list, so the output order is the
+ inverse of the input order. For the position, we will store the
+ values into the `x' and `y' arrays even if they are RA/Dec. */
+ switch(++counter)
+ {
+ case 1:
+ case 2:
+ case 3:
++coordcounter;
- colname = ( counter==8
+ colname = ( counter==1
? "first coordinate column (`--coordcol')"
- : "second 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 8: p->x=corrtype->array; break;
- case 9: p->y=corrtype->array; break;
+ 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_data_malloc_array(GAL_TYPE_UINT8, p->num,
+ __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;
+ 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;
+ break;
+
+ case 11:
+ colname="axis ratio 2 (`q2col')";
+ corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
+ p->q2=corrtype->array;
+ 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;
+ break;
+
/* If the index isn't recognized, then it is larger, showing that
there was more than one match for the given criteria */
default:
@@ -804,13 +1009,15 @@ ui_read_cols(struct mkprofparams *p)
+
/* 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)
{
- float r, n, t;
double *karr;
+ float r, n, t;
+ size_t ndim=p->out->ndim;
/* If the kernel option was called, then we need to build a series of
single element columns to create an internal catalog. */
@@ -820,15 +1027,15 @@ ui_prepare_columns(struct mkprofparams *p)
p->num=1;
/* Allocate the necessary columns. */
- p->x=gal_data_malloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->x");
- p->y=gal_data_malloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->y");
- p->f=gal_data_malloc_array(GAL_TYPE_UINT8, 1, __func__, "p->f");
- p->r=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->r");
- p->n=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->n");
- p->p=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p");
- p->q=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->q");
- p->m=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->m");
- p->t=gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->t");
+ p->x = gal_data_malloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->x");
+ p->y = gal_data_malloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->y");
+ p->f = gal_data_malloc_array(GAL_TYPE_UINT8, 1, __func__, "p->f");
+ p->r = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->r");
+ p->n = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->n");
+ p->p1 = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p");
+ p->q1 = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->q");
+ p->m = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->m");
+ p->t = gal_data_malloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->t");
/* Set the values that need special consideration. */
if(p->kernel->size)
@@ -841,18 +1048,38 @@ ui_prepare_columns(struct mkprofparams *p)
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;
}
else
- ui_read_cols(p);
+ {
+ /* Make sure the number of coordinate columns and number of
+ dimensions in outputs are the same. There is no problem if it is
+ more than `ndim'. In that case, the last values (possibly in
+ configuration files) will be ignored.*/
+ if( gal_list_str_number(p->ccol) < ndim )
+ error(EXIT_FAILURE, 0, "%zu coordinate columns (calls to "
+ "`--coordcol') given but output has %zu dimensions",
+ gal_list_str_number(p->ccol), p->out->ndim);
+
+ /* Call the respective function. */
+ switch(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->out->ndim'",
+ __func__, PACKAGE_BUGREPORT, ndim);
+ }
+ }
}
@@ -985,7 +1212,6 @@ ui_prepare_wcs(struct mkprofparams *p)
}
for(i=0;i<ndim*ndim;++i) wcs->pc[i]=pc[i];
-
/* Set up the wcs structure with the constants defined above. */
status=wcsset(wcs);
if(status)
@@ -1001,10 +1227,10 @@ static void
ui_prepare_canvas(struct mkprofparams *p)
{
float *f, *ff;
- double truncr;
long width[2]={1,1};
int status=0, setshift=0;
size_t i, ndim, nshift=0, *dsize=NULL;
+ double truncr, semiaxes[3], euler_deg[3];
/* If a background image is specified, then use that as the output
image to build the profiles over. */
@@ -1035,10 +1261,6 @@ ui_prepare_canvas(struct mkprofparams *p)
p->cp.minmapsize);
ndim=p->out->ndim;
- /* Currently this only works on 2D images. */
- if(p->out->ndim!=2)
- error(EXIT_FAILURE, 0, "currently only works on 2D images");
-
/* If p->dsize was given as an option, free it. */
if( p->dsize ) free(p->dsize);
@@ -1067,10 +1289,6 @@ ui_prepare_canvas(struct mkprofparams *p)
{
/* Get the number of dimensions. */
ndim=0; for(i=0;p->dsize[i]!=GAL_BLANK_SIZE_T;++i) ++ndim;
- if(ndim!=2)
- error(EXIT_FAILURE, 0, "currently only 2D images are usable, you "
- "have provided %zu values for `--naxis'", ndim);
-
/* If any of the shift elements are zero, the others should be too!*/
if(p->shift && p->shift[0] && p->shift[1])
@@ -1109,8 +1327,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_ellipse_in_box(truncr, p->q[i]*truncr,
- p->p[i]*DEGREESTORADIANS, width);
+ if(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
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index e11cd3e..d75b774 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -65,7 +65,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/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 9c27c60..f58eedc 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -137,9 +137,10 @@ threshold_apply_on_thread(void *in_prm)
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can
"
- "address the problem. A value of %d had for `taprm->kind' "
- "is not valid", __func__, PACKAGE_BUGREPORT, taprm->kind);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we "
+ "can address the problem. A value of %d had for "
+ "`taprm->kind' is not valid", __func__, PACKAGE_BUGREPORT,
+ taprm->kind);
}
}
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 5000605..0e2df12 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -758,7 +758,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 487e911..ae66b75 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -478,7 +478,7 @@ MakeProfiles
Modeling basics
-* Defining an ellipse:: Defining an ellipse in 2D.
+* 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.
@@ -488,8 +488,9 @@ Modeling basics
Invoking MakeProfiles
* MakeProfiles catalog:: Required catalog properties.
-* MakeProfiles options:: Full list of MakeProfiles options.
-* MakeProfiles output:: The generated outputs.
+* MakeProfiles profile settings:: Configuration parameters for all profiles.
+* MakeProfiles output dataset:: The canvas/dataset to build profiles over.
+* MakeProfiles log file:: A description of the optional log file.
MakeNoise
@@ -1540,18 +1541,19 @@ suggestions during the GNU evaluation process. Bob
Proulx from Savannah,
has kindly supported Gnuastro's project webpage on Savannah and the
management of its version controlled source server there.
-@c To the developers: if you would like to thank someone, please change the
-@c `and' before the last name to a comma (','), then add " and Firstname
-@c Familyname'' after the last name.
-We would also like to gratefully thank Mohammad-reza Khellat, Alan Lefor,
-Yahya Sefidbakht, Francesco Montanari, Ole Streicher, Nicolas Bouché, Rosa
-Calvi, Lee Kelvin, Guillaume Mahler, William Pence, David Valls-Gabaud,
-Christopher Willmer, for their useful and constructive comments and
-suggestions. Finally we should thank all the (sometimes anonymous)
-developers in various online forums which patiently answered all our small
-(but important) technical questions. All work on Gnuastro has been
-voluntary, but we are most grateful to the following institutions (in
-chronological order) for hosting us in our research:
+@c To the developers: please keep this in the same order as the THANKS file
+@c (alphabetical, except for the names in the paragraph above).
+We would also like to gratefully thank (in alphabetical order) Roland
+Bacon, Nicolas Bouché, Adrian Bunk, Rosa Calvi, Antonio Diaz Diaz, Raúl
+Infante Sainz, Lee Kelvin, Mohammad-Reza Khellat, Alan Lefor, Guillaume
+Mahler, Francesco Montanari, William Pence, Yahya Sefidbakht, Ole
+Streicher, Ignacio Trujillo, David Valls-Gabaud and Christopher Willmer for
+their useful and constructive comments and suggestions. Finally we should
+thank all the (sometimes anonymous) people in various online forums which
+patiently answered all our small (but important) technical questions. All
+work on Gnuastro has been voluntary, but we are most grateful to the
+following institutions (in chronological order) for hosting us in our
+research:
@quotation
Ministry of education, culture, sports, science and technology (MEXT), Japan.@*
@@ -10651,12 +10653,12 @@ transpose of the one discussed here.}.
@cindex Multiplication, matrix
@cindex Non-commutative operations
@cindex Operations, non-commutative
-In @ref{Warping basics} we saw how one basic warping/transformation
-can be represented with a 3 by 3 matrix. To make more complex warpings
-these matrices have to be multiplied through matrix
-multiplication. However matrix multiplication is not commutative, so
-the order of the set of matrices you use for the multiplication is
-going to be very important.
+In @ref{Warping basics} we saw how a basic warp/transformation can be
+represented with a matrix. To make more complex warpings (for example to
+define a translattion, rotation and scale as one warp) the individual
+matrices have to be multiplied through matrix multiplication. However
+matrix multiplication is not commutative, so the order of the set of
+matrices you use for the multiplication is going to be very important.
The first warping should be placed as the left-most matrix. The second
warping to the right of that and so on. The second transformation is
@@ -14172,7 +14174,7 @@ given. You can skip this subsection if you are already
sufficiently
familiar with these concepts.
@menu
-* Defining an ellipse:: Defining an ellipse in 2D.
+* 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.
@@ -14180,8 +14182,8 @@ familiar with these concepts.
* 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
@@ -14204,34 +14206,70 @@ distance on the major axis to the center of an
ellipse which is located at
want to find @mymath{r_{el}} of a point located at @mymath{(i,j)} (in the
image coordinate system) from the center of the ellipse with axis ratio
@mymath{q} and position angle @mymath{\theta}. First the coordinate system
-is rotated@footnote{Do not confuse this with the rotation matrix of
-@ref{Warping basics}. In that equation, the point is rotated, here the
-coordinates are rotated and the point is fixed.} by @mymath{\theta} to get
-the new rotated coordinates of that point @mymath{(i_r,j_r)}:
+is rotated@footnote{Do not confuse the signs of @mymath{sin} with the
+rotation matrix defined in @ref{Warping basics}. In that equation, the
+point is rotated, here the coordinates are rotated and the point is fixed.}
+by @mymath{\theta} to get the new rotated coordinates of that point
+@mymath{(i_r,j_r)}:
-@dispmath{i_r(i,j)=+(i_c-i)\cos(\theta)+(j_c-j)\sin(\theta)}
-@dispmath{j_r(i,j)=-(i_c-i)\sin(\theta)+(j_c-j)\cos(\theta)}
+@dispmath{i_r(i,j)=+(i_c-i)\cos\theta+(j_c-j)\sin\theta}
+@dispmath{j_r(i,j)=-(i_c-i)\sin\theta+(j_c-j)\cos\theta}
@cindex Elliptical distance
@noindent Recall that an ellipse is defined by @mymath{(i_r/a)^2+(j_r/b)^2=1}
and that we defined @mymath{r_{el}\equiv{a}}. Hence, multiplying all
elements of the 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^2/q^2} }. To place the radial profiles
+@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
@cindex Pixel by pixel making of profiles
-MakeProfiles builds the profile starting from the nearest pixel in the
-image to the given profile center. The profile value is calculated for that
-central pixel using monte carlo integration, see @ref{Sampling from a
-function}. The next pixel is the next nearest neighbor to the central pixel
-as defined by @mymath{r_{el}}. This process goes on until the profile is
-fully built upto the trunctation radius. This is done fairly efficiently
-using a breadth first parsing
+MakeProfiles builds the profile starting from the nearest element (pixel in
+an image) in the dataset to the profile center. The profile value is
+calculated for that central pixel using monte carlo integration, see
+@ref{Sampling from a function}. The next pixel is the next nearest neighbor
+to the central pixel as defined by @mymath{r_{el}}. This process goes on
+until the profile is fully built upto the trunctation radius. This is done
+fairly efficiently using a breadth first parsing
strategy@footnote{@url{http://en.wikipedia.org/wiki/Breadth-first_search}}
which is implemented through an ordered linked list.
@@ -14246,7 +14284,7 @@ dimensions we are dealing with.
-@node PSF, Stars, Defining an ellipse, Modeling basics
+@node PSF, Stars, Defining an ellipse and ellipsoid, Modeling basics
@subsubsection Point Spread Function
@cindex PSF
@@ -14442,7 +14480,7 @@ 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}).
+pixel, see @ref{Defining an ellipse and ellipsoid}).
@cindex Integration over pixel
@cindex Gradient over pixel area
@@ -14479,15 +14517,15 @@ 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}. 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}}. In other words, the gradient is strongest
-along the minor axis. So if the next pixel is chosen based on
-@mymath{r_{el}}, the tolerance level will be reached sooner and lots
-of pixels with large fractional differences will be missed.
+@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}}. In other words, the gradient is strongest along the minor
+axis. So if the next pixel is chosen based on @mymath{r_{el}}, the
+tolerance level will be reached sooner and lots of pixels with large
+fractional differences will be missed.
Monte Carlo integration uses a random number of points. Thus,
every time you run it, by default, you will get a different
@@ -14723,20 +14761,21 @@ $ astmkprof --individual --oversample 3
--naxis=500,500 catalog.txt
@noindent
The parameters of the mock profiles can either be given through a catalog
-(which stores the parameters of many mock profiles), or the
-@option{--kernel} option (see @ref{MakeProfiles options}). The catalog can
-be in the FITS ASCII, FITS binary format, or plain text formats (see
-@ref{Tables}). The columns related to each parameter can be determined both
-by number, or by match/search criteria using the column names, units, or
-comments. with the options ending in @option{col}, see below.
+(which stores the parameters of many mock profiles, see @ref{MakeProfiles
+catalog}), or the @option{--kernel} option (see @ref{MakeProfiles output
+dataset}). The catalog can be in the FITS ASCII, FITS binary format, or
+plain text formats (see @ref{Tables}). The columns related to each
+parameter can be determined both by number, or by match/search criteria
+using the column names, units, or comments. with the options ending in
+@option{col}, see below.
Without any file given to the @option{--background} option, MakeProfiles
will make a zero-valued image and build the profiles on that (its size and
-main WCS parameters can also be defined through the options). Besides the
-main/merged image containing all the profiles in the catalog, it is also
-possible to build individual images for each profile (only enclosing one
-full profile to its truncation radius) with the @option{--individual}
-option.
+main WCS parameters can also be defined through the options described in
+@ref{MakeProfiles output dataset}). Besides the main/merged image
+containing all the profiles in the catalog, it is also possible to build
+individual images for each profile (only enclosing one full profile to its
+truncation radius) with the @option{--individual} option.
If an image is given to the @option{--background} option, the pixels of
that image are used as the background value for every pixel. The flux value
@@ -14746,41 +14785,57 @@ and WCS will be ignored if specified (for example
@option{--oversample},
@option{--naxis1}, @option{--naxis2} and @option{--prepforconv}) on the
command-line or in the configuration files.
+The sections below discuss the options specific to MakeProfiles based on
+context: the input catalog settings which can have many rows for different
+profiles are discussed in @ref{MakeProfiles catalog}, in @ref{MakeProfiles
+profile settings}, we discuss how you can set general profile settings
+(that are the same for all the profiles in the catalog). 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.
+
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
* MakeProfiles catalog:: Required catalog properties.
-* MakeProfiles options:: Full list of MakeProfiles options.
-* MakeProfiles output:: The generated outputs.
+* MakeProfiles profile settings:: Configuration parameters for all profiles.
+* MakeProfiles output dataset:: The canvas/dataset to build profiles over.
+* MakeProfiles log file:: A description of the optional log file.
@end menu
-@node MakeProfiles catalog, MakeProfiles options, Invoking astmkprof, Invoking
astmkprof
+@node MakeProfiles catalog, MakeProfiles profile settings, Invoking astmkprof,
Invoking astmkprof
@subsubsection MakeProfiles catalog
-The catalog can be in the FITS ASCII, FITS binary format, or plain text
-formats (see @ref{Tables}). Its columns can be ordered in any desired
-manner, you can specify which columns belong to which parameters using the
-set of options ending with @option{col} in @ref{MakeProfiles options}. For
-example through the @option{--xcol} and @option{--rcol} options, you can
-specify the column that contains the X axis position of the profile center
-and the radial parameter for the profile. See @ref{Selecting table columns}
-for a thorough discussion on how to identify a column in the input catalog.
-
-The value for the profile center in the catalog (in the @option{--xcol} and
-@option{--ycol} columns) can be a floating point number so the profile
-center can be on any sub-pixel position. Note that pixel positions in the
-FITS standard start from 1 and an integer is the pixel center. So a 2D
-image actually starts from the position (0.5, 0.5). When a
-@option{--background} image with WCS information is provided, you may also
-use RA and Dec to identify the center of each profile.
+The catalog containing information about each profile can be in the FITS
+ASCII, FITS binary, or plain text formats (see @ref{Tables}). Its columns
+can be ordered in any desired manner. You can specify which columns belong
+to which parameters using the set of options discussed below. For example
+through the @option{--rcol} and @option{--tcol} options, you can specify
+the column that contains the radial parameter for each profile and its
+truncation respectively. See @ref{Selecting table columns} for a thorough
+discussion on the values to these options.
+
+The value for the profile center in the catalog (the @option{--ccol}
+option) can be a floating point number so the profile center can be on any
+sub-pixel position. Note that pixel positions in the FITS standard start
+from 1 and an integer is the pixel center. So a 2D image actually starts
+from the position (0.5, 0.5), which is the bottom-left corner of the first
+pixel. When a @option{--background} image with WCS information is provided
+or you specify the WCS parameters with the respective options, you may also
+use RA and Dec to identify the center of each profile (see the
+@option{--mode} option below).
In MakeProfiles, profile centers do not have to be in (overlap with) the
final image. Even if only one pixel of the profile within the truncation
radius overlaps with the final image size, the profile is built and
included in the final image image. Profiles that are completely out of the
-image will not be created. You can use the output log file to see which
-profiles were within the image.
+image will not be created (unless you explicity ask for it with the
+@option{--individual} option). You can use the output log file (created
+with @option{--log} to see which profiles were within the image, see
+@ref{Common options}.
If PSF profiles (Moffat or Gaussian, see @ref{PSF}) are in the catalog and
the profiles are to be built in one image (when @option{--individual} is
@@ -14796,153 +14851,128 @@ unity, you have to set their magnitudes in the
catalog to the zero-point
magnitude and be sure that the central positions of the profiles don't have
any fractional part (the PSF center has to be in the center of the pixel).
-
-@node MakeProfiles options, MakeProfiles output, MakeProfiles catalog,
Invoking astmkprof
-@subsubsection MakeProfiles options
-The common options that are shared by Gnuastro programs, are fully
-explained in @ref{Common options} and are not repeated here. Since
-there are no image inputs, the@option{--hdu} option is ignored. The
-options can be classified into the following categories: Output,
-Profiles, Catalog and WCS. Below each one is reviewed.
-
-@noindent
-Output:
-
-@table @option
-
-@item -E STR/INT,FLT[,FLT,[...]]
-@itemx --kernel=STR/INT,FLT[,FLT,[...]]
-Only build one kernel profile with the parameters given as the values to
-this option. The different values must be separated by a comma
-(@key{,}). The first value identifies the radial function of the profile,
-either through a string or through a number (see description of
-@option{--fcol} below). Each radial profile needs a different total number
-of parameters: S@'ersic and Moffat functions need 3 parameters: radial,
-S@'ersic index or Moffat @mymath{\beta}, and truncation radius. The
-Gaussian function needs two parameters: radial and truncation radius. The
-point function doesn't need any parameters and flat and circumference
-profiles just need one parameter (truncation radius).
-
-The PSF or kernel is a unique (and highly constrained) type of profile: the
-sum of its pixels must be one, its center must be the center of the central
-pixel (in an image with an odd number of pixels on each side), and commonly
-it is circular, so its axis ratio and position angle are one and zero
-respectively. Kernels are commonly necessary for various data analysis and
-data manipulation steps (for example see @ref{Convolve}, and
-@ref{NoiseChisel}. Because of this it is inconvenient to define a catalog
-with one row and many zero valued columns (for all the non-necessary
-parameters). Hence, with this option, it is possible to create a kernel
-with MakeProfiles without the need to create a catalog. Here are some
-examples:
+The list of options directly related to the input catalog columns is shown
+below.
@table @option
-@item --kernel=moffat,3,2.8,5
-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.
-@end table
-
-@item -k STR
-@itemx --background=STR
-A background image FITS file to build the profiles on. The extension that
-contains the image should be specified with the @option{--backhdu} option,
-see below. When a background image is specified, it will be used to derive
-all the information about the output image. Hence, the following options
-will be ignored: @option{--naxis1}, @option{--naxis2}, @option{--crpix1},
-@option{--crpix2}, @option{--crval1}, @option{--crval2},
-@option{--resolution}, @option{--oversample}, and data type (see
-@option{--type} in @ref{Input output options}).
+@item --mode=STR
+Interpret the center position columns in image or WCS coordinates. This
+option thus accepts only two values: @option{img} and @option{wcs}. It is
+mandatory when a catalog is being used.
-The image will act like a canvas to build the profiles on: profile pixel
-values will be summed with the background image pixel values. With the
-@option{--replace} option you can disable this behavior and replace the
-profile pixels with the background pixels. If you want to use all the image
-information above, except for the pixel values (you want to have a blank
-canvas to build the profiles on, based on an input image), you can call
-@option{--clearcanvas}, to set all the input image's pixels to zero before
-starting to build the profiles over it (this is done in memory after
-reading the input, so nothing will happen to your input file).
+@item --ccol=STR/INT
+Center coordinate column for each dimension. This option must be called two
+times to define the center coordinates in an image. For example
+@option{--ccol=RA} and @option{--ccol=DEC} (along with @option{--mode=wcs})
+will inform MakeProfiles to look into the catalog columns named @option{RA}
+and @option{DEC} for the Right Ascension and Declination of the profile
+centers.
-@item -B STR/INT
-@itemx --backhdu=STR/INT
-The header data unit (HDU) of the file given to @option{--background}.
+@item --fcol=INT/STR
+The functional form of the profile with one of the values below depending
+on the desired profile. The column can contain either the numeric codes
+(for example `@code{1}') or string characters (for example
+`@code{sersic}'). The numeric codes are easier to use in scripts which
+generate catalogs with hundreds or thousands of profiles.
-@item -x INT,INT
-@itemx --naxis=INT,INT
-The number of pixels along each dimension axis of the output in FITS
-order. This is before over-sampling. For example if you call MakeProfiles
-with @option{--naxis=100,150 --oversample=5} (assuming no shift due for
-later convolution), then the final image size along the first axis will be
-500 by 750 pixels. Fractions are acceptable as values for each dimension,
-however, they must reduce to an integer, so @option{--naxis=150/3,300/3} is
-acceptable but @option{--naxis=150/4,300/4} is not.
+The string format can be easier when the catalog is to be written/checked
+by hand/eye before running MakeProfiles. It is much more readable and
+provides a level of documentation. All Gnuastro's recognized table formats
+(see @ref{Recognized table formats}) accept string type columns. To have
+string columns in a plain text table/catalog, see @ref{Gnuastro text table
+format}.
-When viewing a FITS image in DS9, the first FITS dimension is in the
-horizontal direction and the second is vertical. As an example, the image
-created with the example above will have 500 pixels horizontally and 750
-pixels vertically.
+@itemize
+@item
+S@'ersic profile with `@code{sersic}' or `@code{1}'.
+@item
+Moffat profile with `@code{moffat}' or `@code{2}'.
+@item
+Gaussian profile with `@code{gaussian}' or `@code{3}'.
+@item
+Point source with `@code{point}' or `@code{4}'.
+@item
+Flat profile with `@code{flat}' or `@code{5}'.
+@item
+Circumference profile with `@code{circum}' or `@code{6}'. A fixed
+value will be used for all pixels between the truncation radius
+(@mymath{r_t}) and @mymath{r_t-w} (@mymath{w} is the value to the
+@option{--circumwidth}).
+@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 orentiation (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. You can use this
+to define your own higher-level radial function. Note that for this
+profile, the value in the magnitude column (@option{--mcol}) will be
+ignored.
+@end itemize
-If a background image is specified, this option is ignored.
-@item -s INT
-@itemx --oversample=INT
-The scale to over-sample the profiles and final image. If not an odd
-number, will be added by one, see @ref{Oversampling}. Note that this
-@option{--oversample} will remain active even if an input image is
-specified. If your input catalog is based on the background image, be sure
-to set @option{--oversample=1}.
+@item --rcol=STR/INT
+The radius parameter of the profiles. Effective radius (@mymath{r_e}) if
+S@'ersic, FWHM if Moffat or Gaussian.
-@item --psfinimg
-Build the possibly existing PSF profiles (Moffat or Gaussian) in the
-catalog into the final image. By default they are built separately so
-you can convolve your images with them, thus their magnitude and
-positions are ignored. With this option, they will be built in the
-final image like every other galaxy profile. To have a final PSF in
-your image, make a point profile where you want the PSF and after
-convolution it will be the PSF.
+@item --ncol=STR/INT
+The S@'ersic index (@mymath{n}) or Moffat @mymath{\beta}.
-@item -i
-@itemx --individual
-@cindex Individual profiles
-@cindex Build individual profiles
-If this option is called, each profile is created in a separate FITS
-file within the same directory as the output and the row number of the
-profile (starting from zero) in the name. The file for each row's
-profile will be in the same directory as the final combined image of
-all the profiles and will have the final image's name as a suffix. So
-for example if the final combined image is named
-@file{./out/fromcatalog.fits}, then the first profile that will be
-created with this option will be named
-@file{./out/0_fromcatalog.fits}.
+@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.
-Since each image only has one full profile out to the truncation
-radius the profile is centered and so, only the sub-pixel position of
-the profile center is important for the outputs of this option. The
-output will have an odd number of pixels. If there is no oversampling,
-the central pixel will contain the profile center. If the value to
-@option{--oversample} is larger than unity, then the profile center is
-on any of the central @option{--oversample}'d pixels depending on the
-fractional value of the profile center.
+@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.
-If the fractional value is larger than half, it is on the bottom half
-of the central region. This is due to the FITS definition of a real
-number position: The center of a pixel has fractional value
-@mymath{0.00} so each pixel contains these fractions: .5 -- .75 -- .00
-(pixel center) -- .25 -- .5.
+@item --mcol=STR/INT
+The total pixelated magnitude of the profile within the truncation radius,
+see @ref{Profile magnitude}.
-@item -m
-@itemx --nomerged
-Don't make a merged image. By default after making the profiles, they
-are added to a final image with sides of @option{--naxis1} and
-@option{--naxis2} if they overlap with it.
+@item --tcol=STR/INT
+The truncation radius of this profile. By default it is in units of the
+radial parameter of the profile (the value in the @option{--rcol} of the
+catalog). If @option{--tunitinp} is given, this value is interpreted in
+units of pixels (prior to oversampling) irrespective of the profile.
@end table
-@noindent
-Profiles:
+@node MakeProfiles profile settings, MakeProfiles output dataset, MakeProfiles
catalog, Invoking astmkprof
+@subsubsection MakeProfiles profile settings
+
+The profile parameters that differ between each created profile are
+specified through the columns in the input catalog and described in
+@ref{MakeProfiles catalog}. Besides those there are general settings for
+some profiles that don't differ between one profile and another, they are a
+property of the general process. For example how many random points to use
+in the monte-carlo integration, this value is fixed for all the
+profiles. The options described in this section are for configuring such
+properties.
@table @option
@@ -14970,25 +15000,42 @@ default, the truncation column is considered to be in
units of the
radial parameters of the profile (@option{--rcol}). Read it as
`t-unit-in-p' for `truncation unit in pixels'.
-@item -X INT,INT
-@itemx --shift=INT,INT
-Shift all the profiles and enlarge the image along each dimension. To
-better understand this option, please see @mymath{n} in @ref{If convolving
-afterwards}. This is useful when you want to convolve the image
-afterwards. If you are using an external PSF, be sure to oversample it to
-the same scale used for creating the mock images. If a background image is
-specified, any possible value to this option is ignored.
+@item -f
+@itemx --mforflatpix
+When making fixed value profiles (flat and circumference, see
+`@option{--fcol}'), don't use the value in the column specified by
+`@option{--mcol}' as the magnitude. Instead use it as the exact value that
+all the pixels of these profiles should have. This option is irrelevant for
+other types of profiles. This option is very useful for creating masks, or
+labeled regions in an image. Any integer, or floating point value can used
+in this column with this option, including @code{NaN} (or `@code{nan}', or
+`@code{NAN}', case is irrelevant), and infinities (@code{inf}, @code{-inf},
+or @code{+inf}).
-@item -c
-@itemx --prepforconv
-Shift all the profiles and enlarge the image based on half the width
-of the first Moffat or Gaussian profile in the catalog, considering
-any possible oversampling see @ref{If convolving
-afterwards}. @option{--prepforconv} is only checked and possibly
-activated if @option{--xshift} and @option{--yshift} are both zero
-(after reading the command-line and configuration files). If a
-background image is specified, any possible value to this option is
-ignored.
+For example, with this option if you set the value in the magnitude column
+(@option{--mcol}) to @code{NaN}, you can create an elliptical or circular
+mask over an image (which can be given as the argument), see @ref{Blank
+pixels}. Another useful application of this option is to create labeled
+elliptical or circular apertures in an image. To do this, set the value in
+the magnitude column to the label you want for this profile. This labeled
+image can then be used in combination with NoiseChisel's output (see
+@ref{NoiseChisel output}) to do aperture photometry with MakeCatalog (see
+@ref{MakeCatalog}).
+
+Alternatively, if you want to mark regions of the image (for example with
+an elliptical circumference) and you don't want to use NaN values (as
+explained above) for some technical reason, you can get the minimum or
+maximum value in the image @footnote{The minimum will give a better result,
+because the maximum can be too high compared to most pixels in the image,
+making it harder to display.} using Arithmetic (see @ref{Arithmetic}), then
+use that value in the magnitude column along with this option for all the
+profiles.
+
+Please note that when using MakeProfiles on an already existing image, you
+have to set `@option{--oversample=1}'. Otherwise all the profiles will be
+scaled up based on the oversampling scale in your configuration files (see
+@ref{Configuration files}) unless you have accounted for oversampling in
+your catalog.
@item --magatpeak
The magnitude column in the catalog (see @ref{MakeProfiles catalog})
@@ -15019,6 +15066,36 @@ have to set @option{--oversample=1} otherwise after
resampling your profile
with Warp (see @ref{Warp}), the peak flux will be different.
@end cartouche
+@item -X INT,INT
+@itemx --shift=INT,INT
+Shift all the profiles and enlarge the image along each dimension. To
+better understand this option, please see @mymath{n} in @ref{If convolving
+afterwards}. This is useful when you want to convolve the image
+afterwards. If you are using an external PSF, be sure to oversample it to
+the same scale used for creating the mock images. If a background image is
+specified, any possible value to this option is ignored.
+
+@item -c
+@itemx --prepforconv
+Shift all the profiles and enlarge the image based on half the width
+of the first Moffat or Gaussian profile in the catalog, considering
+any possible oversampling see @ref{If convolving
+afterwards}. @option{--prepforconv} is only checked and possibly
+activated if @option{--xshift} and @option{--yshift} are both zero
+(after reading the command-line and configuration files). If a
+background image is specified, any possible value to this option is
+ignored.
+
+@item -z FLT
+@itemx --zeropoint=FLT
+The zero-point magnitude of the image.
+
+@item -w FLT
+@itemx --circumwidth=FLT
+The width of the circumference if the profile is to be an elliptical
+circumference or annulus. See the explanations for this type of profile in
+@option{--fcol}.
+
@item -R
@itemx --replace
Do not add the pixels of each profile over the background (possibly crowded
@@ -15043,6 +15120,44 @@ profile. However, when using flat profiles with the
`@option{--mforflatpix}' option, you should be careful not to give a
@code{0.0} value as the flat profile's pixel value.
+@end table
+
+@node MakeProfiles output dataset, MakeProfiles log file, MakeProfiles profile
settings, Invoking astmkprof
+@subsubsection MakeProfiles output dataset
+MakeProfiles takes an input catalog uses basic properties that are defined
+there to build a dataset, for example a 2D image containing the profiles in
+the catalog. In @ref{MakeProfiles catalog} and @ref{MakeProfiles profile
+settings}, the catalog and profile settings were discussed. The options of
+this section, allow you to configure the output dataset (or the canvas that
+will host the built profiles).
+
+@table @option
+
+@item -k STR
+@itemx --background=STR
+A background image FITS file to build the profiles on. The extension that
+contains the image should be specified with the @option{--backhdu} option,
+see below. When a background image is specified, it will be used to derive
+all the information about the output image. Hence, the following options
+will be ignored: @option{--naxis1}, @option{--naxis2}, @option{--crpix1},
+@option{--crpix2}, @option{--crval1}, @option{--crval2},
+@option{--resolution}, @option{--oversample}, and data type (see
+@option{--type} in @ref{Input output options}).
+
+The image will act like a canvas to build the profiles on: profile pixel
+values will be summed with the background image pixel values. With the
+@option{--replace} option you can disable this behavior and replace the
+profile pixels with the background pixels. If you want to use all the image
+information above, except for the pixel values (you want to have a blank
+canvas to build the profiles on, based on an input image), you can call
+@option{--clearcanvas}, to set all the input image's pixels to zero before
+starting to build the profiles over it (this is done in memory after
+reading the input, so nothing will happen to your input file).
+
+@item -B STR/INT
+@itemx --backhdu=STR/INT
+The header data unit (HDU) of the file given to @option{--background}.
+
@item -C
@itemx --clearcanvas
When an input image is specified (with the @option{--background} option,
@@ -15057,132 +15172,112 @@ catalog (see @ref{MakeCatalog}). In other cases,
you might have modeled the
objects in an image and want to create them on the same frame, but without
the original pixel values.
-@item -w FLT
-@itemx --circumwidth=FLT
-The width of the circumference if the profile is to be an elliptical
-circumference or annulus. See the explanations for this type of profile in
-@option{--fcol}.
-
-@item -z FLT
-@itemx --zeropoint=FLT
-The zero-point magnitude of the image.
-
-@end table
+@item -E STR/INT,FLT[,FLT,[...]]
+@itemx --kernel=STR/INT,FLT[,FLT,[...]]
+Only build one kernel profile with the parameters given as the values to
+this option. The different values must be separated by a comma
+(@key{,}). The first value identifies the radial function of the profile,
+either through a string or through a number (see description of
+@option{--fcol} below). Each radial profile needs a different total number
+of parameters: S@'ersic and Moffat functions need 3 parameters: radial,
+S@'ersic index or Moffat @mymath{\beta}, and truncation radius. The
+Gaussian function needs two parameters: radial and truncation radius. The
+point function doesn't need any parameters and flat and circumference
+profiles just need one parameter (truncation radius).
-@noindent
-Catalog: The value to all of these options is considered to be a
-column number, where counting starts from zero.
+The PSF or kernel is a unique (and highly constrained) type of profile: the
+sum of its pixels must be one, its center must be the center of the central
+pixel (in an image with an odd number of pixels on each side), and commonly
+it is circular, so its axis ratio and position angle are one and zero
+respectively. Kernels are commonly necessary for various data analysis and
+data manipulation steps (for example see @ref{Convolve}, and
+@ref{NoiseChisel}. Because of this it is inconvenient to define a catalog
+with one row and many zero valued columns (for all the non-necessary
+parameters). Hence, with this option, it is possible to create a kernel
+with MakeProfiles without the need to create a catalog. Here are some
+examples:
@table @option
+@item --kernel=moffat,3,2.8,5
+A Moffat kernel with FWHM of 3 pixels, @mymath{\beta=2.8} which is
+truncated at 5 times the FWHM.
-@item --mode=STR
-Interpret the center position columns in image or WCS coordinates. This
-option thus accepts only two values: @option{img} and @option{wcs}. It is
-mandatory when a catalog is being used.
-
-@item --ccol=STR/INT
-Center coordinate column for each dimension. This option must be called two
-times to define the center coordinates in an image. For example
-@option{--ccol=RA} and @option{--ccol=DEC} (along with @option{--mode=wcs})
-will inform MakeProfiles to look into the catalog columns named @option{RA}
-and @option{DEC} for the Right Ascension and Declination of the profile
-centers.
-
-@item --fcol=INT/STR
-The functional form of the profile with one of the values below depending
-on the desired profile. The column can contain either the numeric codes
-(for example `@code{1}') or string characters (for example
-`@code{sersic}'). The numeric codes are easier to use in scripts which
-generate catalogs with hundreds or thousands of profiles.
-
-The string format can be easier when the catalog is to be written/checked
-by hand/eye before running MakeProfiles. It is much more readable and
-provides a level of documentation. All Gnuastro's recognized table formats
-(see @ref{Recognized table formats}) accept string type columns. To have
-string columns in a plain text table/catalog, see @ref{Gnuastro text table
-format}.
-
-@itemize
-@item
-S@'ersic profile with `@code{sersic}' or `@code{1}'.
-@item
-Moffat profile with `@code{moffat}' or `@code{2}'.
-@item
-Gaussian profile with `@code{gaussian}' or `@code{3}'.
-@item
-Point source with `@code{point}' or `@code{4}'.
-@item
-Flat profile with `@code{flat}' or `@code{5}'.
-@item
-Circumference profile with `@code{circum}' or `@code{6}'. A fixed
-value will be used for all pixels between the truncation radius
-(@mymath{r_t}) and @mymath{r_t-w} (@mymath{w} is the value to the
-@option{--circumwidth}).
-@end itemize
-
-@item --rcol=STR/INT
-The radius parameter of the profiles. Effective radius (@mymath{r_e}) if
-S@'ersic, FWHM if Moffat or Gaussian.
+@item --kernel=gaussian,2,3
+A Gaussian kernel with FWHM of 2 pixels and truncated at 3 times the FWHM.
+@end table
-@item --ncol=STR/INT
-The S@'ersic index (@mymath{n}) or Moffat @mymath{\beta}.
+@item -x INT,INT
+@itemx --naxis=INT,INT
+The number of pixels along each dimension axis of the output in FITS
+order. This is before over-sampling. For example if you call MakeProfiles
+with @option{--naxis=100,150 --oversample=5} (assuming no shift due for
+later convolution), then the final image size along the first axis will be
+500 by 750 pixels. Fractions are acceptable as values for each dimension,
+however, they must reduce to an integer, so @option{--naxis=150/3,300/3} is
+acceptable but @option{--naxis=150/4,300/4} is not.
-@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 viewing a FITS image in DS9, the first FITS dimension is in the
+horizontal direction and the second is vertical. As an example, the image
+created with the example above will have 500 pixels horizontally and 750
+pixels vertically.
-@item --qcol=STR/INT
-The axis ratio of the profiles (minor axis divided by the major axis).
+If a background image is specified, this option is ignored.
-@item --mcol=STR/INT
-The total pixelated magnitude of the profile within the truncation radius,
-see @ref{Profile magnitude}.
+@item -s INT
+@itemx --oversample=INT
+The scale to over-sample the profiles and final image. If not an odd
+number, will be added by one, see @ref{Oversampling}. Note that this
+@option{--oversample} will remain active even if an input image is
+specified. If your input catalog is based on the background image, be sure
+to set @option{--oversample=1}.
-@item --tcol=STR/INT
-The truncation radius of this profile. By default it is in units of the
-radial parameter of the profile (the value in the @option{--rcol} of the
-catalog). If @option{--tunitinp} is given, this value is interpreted in
-units of pixels (prior to oversampling) irrespective of the profile.
+@item --psfinimg
+Build the possibly existing PSF profiles (Moffat or Gaussian) in the
+catalog into the final image. By default they are built separately so
+you can convolve your images with them, thus their magnitude and
+positions are ignored. With this option, they will be built in the
+final image like every other galaxy profile. To have a final PSF in
+your image, make a point profile where you want the PSF and after
+convolution it will be the PSF.
-@item -f
-@itemx --mforflatpix
-When making fixed value profiles (flat and circumference, see
-`@option{--fcol}'), don't use the value in the column specified by
-`@option{--mcol}' as the magnitude. Instead use it as the exact value that
-all the pixels of these profiles should have. This option is irrelevant for
-other types of profiles. This option is very useful for creating masks, or
-labeled regions in an image. Any integer, or floating point value can used
-in this column with this option, including @code{NaN} (or `@code{nan}', or
-`@code{NAN}', case is irrelevant), and infinities (@code{inf}, @code{-inf},
-or @code{+inf}).
+@item -i
+@itemx --individual
+@cindex Individual profiles
+@cindex Build individual profiles
+If this option is called, each profile is created in a separate FITS
+file within the same directory as the output and the row number of the
+profile (starting from zero) in the name. The file for each row's
+profile will be in the same directory as the final combined image of
+all the profiles and will have the final image's name as a suffix. So
+for example if the final combined image is named
+@file{./out/fromcatalog.fits}, then the first profile that will be
+created with this option will be named
+@file{./out/0_fromcatalog.fits}.
-For example, with this option if you set the value in the magnitude column
-(@option{--mcol}) to @code{NaN}, you can create an elliptical or circular
-mask over an image (which can be given as the argument), see @ref{Blank
-pixels}. Another useful application of this option is to create labeled
-elliptical or circular apertures in an image. To do this, set the value in
-the magnitude column to the label you want for this profile. This labeled
-image can then be used in combination with NoiseChisel's output (see
-@ref{NoiseChisel output}) to do aperture photometry with MakeCatalog (see
-@ref{MakeCatalog}).
+Since each image only has one full profile out to the truncation
+radius the profile is centered and so, only the sub-pixel position of
+the profile center is important for the outputs of this option. The
+output will have an odd number of pixels. If there is no oversampling,
+the central pixel will contain the profile center. If the value to
+@option{--oversample} is larger than unity, then the profile center is
+on any of the central @option{--oversample}'d pixels depending on the
+fractional value of the profile center.
-Alternatively, if you want to mark regions of the image (for example with
-an elliptical circumference) and you don't want to use NaN values (as
-explained above) for some technical reason, you can get the minimum or
-maximum value in the image @footnote{The minimum will give a better result,
-because the maximum can be too high compared to most pixels in the image,
-making it harder to display.} using Arithmetic (see @ref{Arithmetic}), then
-use that value in the magnitude column along with this option for all the
-profiles.
+If the fractional value is larger than half, it is on the bottom half
+of the central region. This is due to the FITS definition of a real
+number position: The center of a pixel has fractional value
+@mymath{0.00} so each pixel contains these fractions: .5 -- .75 -- .00
+(pixel center) -- .25 -- .5.
-Please note that when using MakeProfiles on an already existing image, you
-have to set `@option{--oversample=1}'. Otherwise all the profiles will be
-scaled up based on the oversampling scale in your configuration files (see
-@ref{Configuration files}) unless you have accounted for oversampling in
-your catalog.
+@item -m
+@itemx --nomerged
+Don't make a merged image. By default after making the profiles, they are
+added to a final image with side lengths specified by @option{--naxis}if
+they overlap with it.
@end table
+
@noindent
The options below can be used to define the world coordinate system (WCS)
properties of the MakeProfiles outputs. The option names are delibarately
@@ -15242,31 +15337,38 @@ later usage of them might cause trouble.
@end table
-@node MakeProfiles output, , MakeProfiles options, Invoking astmkprof
-@subsubsection MakeProfiles output
+@node MakeProfiles log file, , MakeProfiles output dataset, Invoking astmkprof
+@subsubsection MakeProfiles log file
-Besides the final merged image of all the profiles or individual
-profiles that can be built based on the input options, MakeProfiles
-will also create a log file in the current directory (where you run
-MockProfiles). The values for each column are explained in the first
-few commented (starting with @command{#} character). The log file
-includes the following information:
+Besides the final merged dataset of all the profiles, or the individual
+datasets (see @ref{MakeProfiles output dataset}), if the @option{--log}
+option is called MakeProfiles will also create a log file in the current
+directory (where you run MockProfiles). See @ref{Common options} for a full
+description of @option{--log} and other options that are shared between all
+Gnuastro programs. The values for each column are explained in the first
+few commented lines of the log file (starting with @command{#}
+character). Here is a more complete description.
@itemize
@item
-The total magnitude of the profile in the image. This will be
-different from your input magnitude if the profile was not completely
-in the image.
+An ID (row number of profile in input catalog).
+
+@item
+The total magnitude of the profile in the output dataset. When the profile
+does not completely overlap with the output dataset, this will be different
+from your input magnitude.
@item
The number of pixels (in the oversampled image) which used Monte Carlo
-integration and not the central pixel value.
+integration and not the central pixel value, see @ref{Sampling from a
+function}.
@item
The fraction of flux in the Monte Carlo integrated pixels.
@item
-If an individual image was created or not.
+If an individual image was created, this column will have a value of
+@code{1}, otherwise it will have a value of @code{0}.
@end itemize
@@ -18034,9 +18136,10 @@ the dataset, be sure to clean up all the allocated
spaces with
@end deftypefun
@deftypefun void gal_data_free_contents (gal_data_t @code{*data})
-Free all the non-@code{NULL} pointers in @code{gal_data_t}, for a complete
-description of the @code{gal_data_t} contents, see @ref{Generic data
-container}.
+Free all the non-@code{NULL} pointers in @code{gal_data_t}. If @code{data}
+is actually a tile (@code{data->block!=NULL}, see @ref{Tessellation
+library}), then @code{tile->array} is not freed. For a complete description
+of @code{gal_data_t} and its contents, see @ref{Generic data container}.
@end deftypefun
@deftypefun void gal_data_free (gal_data_t @code{*data})
@@ -20486,13 +20589,27 @@ operation will be done on @code{numthreads} threads.
@deffn {Function-like macro} GAL_TILE_PARSE_OPERATE (@code{IN}, @code{OTHER},
@code{PARSE_OTHER}, @code{CHECK_BLANK}, @code{OP})
-Parse over @code{IN} (which can be a tile or a fully allocated block of
-memory) and do the @code{OP} operation on it (can be any combination of C
-expressions). If @code{OTHER!=NULL}, this macro will allow access to its
-element(s) and it can optionally be parsed while parsing over
-@code{IN}. You can use this to parse the same region over two arrays. The
-input arguments are (the expected types of each argument are also written
-here):
+Parse @code{IN} (which can be a tile or a fully allocated block of memory)
+and do the @code{OP} operation on it. @code{OP} can be any combination of C
+expressions. If @code{OTHER!=NULL}, @code{OTHER} will be interpretted as a
+dataset and this macro will allow access to its element(s) and it can
+optionally be parsed while parsing over @code{IN}.
+
+If @code{OTHER} is a fully allocated block of memory (not a tile), then the
+same region that is covered by @code{IN} within its own block will be
+parsed (the same starting pixel with the same number of pixels in each
+dimension). Hence, in this case, the blocks of @code{OTHER} and @code{IN}
+must have the same size. When @code{OTHER} is a tile it must have the same
+size as @code{IN} and parsing will start from its starting element/pixel.
+Also, the respective allocated blocks of @code{OTHER} and @code{IN} (if
+different) may have different sizes. Using @code{OTHER} (along with
+@code{PARSE_OTHER}), this funciton-like macro will thus enable you to parse
+and define your own operation on two fixed size regions in one or two
+blocks of memory. In the latter case, they may have different numeric
+datatypes, see @ref{Numeric data types}).
+
+The input arguments to this macro are explained below, the expected type of
+each argument are also written following the argument name:
@table @code
@item IN (gal_data_t)
@@ -20502,8 +20619,12 @@ Input dataset, this can be a tile or an allocated
block of memory.
Dataset (@code{gal_data_t}) to parse along with @code{IN}. It can be
@code{NULL}. In that case, @code{o} (see description of @code{OP} below)
will be @code{NULL} and should not be used. If @code{PARSE_OTHER} is zero,
-the size of this dataset is irrelevant. Otherwise, it has to have the same
-size as the allocated block of @code{IN}.
+only its first element will be used and the size of this dataset is
+irrelevant.
+
+When @code{OTHER} is a block of memory, it has to have the same size as the
+allocated block of @code{IN}. When its a tile, it has to have the same size
+as @code{IN}.
@item PARSE_OTHER (int)
Parse the other dataset along with the input. When this is non-zero and
@@ -20555,14 +20676,14 @@ Blank value in the type of @code{INPUT}.
@end table
You can use a given tile (@code{tile} on a dataset that it was not
-initialized with (but has the same size, let's call it @code{new}) with the
+initialized with but has the same size, let's call it @code{new}) with the
following steps:
@example
void *tarray;
gal_data_t *tblock;
-/* `tile->block' must corrected AFTER `tile->array' */
+/* `tile->block' must be corrected AFTER `tile->array'. */
tarray = tile->array;
tblock = tile->block;
tile->array = gal_tile_block_relative_to_other(tile, new);
@@ -20602,15 +20723,15 @@ faketile->size = tile->size;
faketile->dsize = tile->dsize;
faketile->array = gal_tile_block_relative_to_other(tile, new);
-// Do your processing....
+/* Do your processing.... in a loop (over many tiles). */
GAL_TILE_PARSE_OPERATE(tile, faketile, 1, 1, @{
YOUR_PROCESSING_EXPRESSIONS;
@});
-// Clean up.
-bintile->array=NULL;
-bintile->dsize=NULL;
-gal_data_free(bintile);
+/* Clean up (outside the loop). */
+faketile->array=NULL;
+faketile->dsize=NULL;
+gal_data_free(faketile);
@end example
@end deffn
@@ -20858,13 +20979,38 @@ declared in @file{gnuastro/box.h}. All coordinates in
this header are in
the FITS format (first axis is the horizontal and the second axis is
vertical).
-@deftypefun void gal_box_ellipse_in_box (double @code{a}, double @code{b},
double @code{theta_rad}, long @code{*width})
+@deftypefun void gal_box_bound_ellipse (double @code{a}, double @code{b},
double @code{theta_deg}, long @code{*width})
Any ellipse can be enclosed into a rectangular box. The purpose of this
-function is to give the height and width of that box. @code{a} is the
-ellipse major axis, @code{b} is the minor axis, @code{theta_rad} is the
-position angle in radians. The @code{width} array will contain the output
-size in long integer type. @code{width[0]}, and @code{width[1]} are the
-number of pixels along the first and second FITS axis.
+function is to give the height and width of that box assuming the center of
+the ellipse is in the box center. @code{a} is the ellipse major axis,
+@code{b} is the minor axis, @code{theta_deg} is the position angle in
+degrees. The @code{width} array will contain the output size in long
+integer type. @code{width[0]}, and @code{width[1]} are the number of pixels
+along the first and second FITS axis. Since the ellipse center is assumed
+to be 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 (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 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})
diff --git a/lib/box.c b/lib/box.c
index ae5e708..8bfe0fb 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -1,5 +1,5 @@
/*********************************************************************
-overlap -- Find the overlap of a region and an image.
+Box -- Define bounding and overlapping boxes.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
@@ -39,10 +39,11 @@ 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 Find the width of
+ this box in each dimension. 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)
@@ -64,20 +65,20 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
and you will find the distance (about the center of the ellipse
that encloses the whole ellipse. */
void
-gal_box_ellipse_in_box(double a, double b, double theta_rad, long *width)
+gal_box_bound_ellipse(double a, double b, double theta_deg, long *width)
{
- double t_x, t_y, max_x, max_y;
- t_x=atan(b/a*tan(theta_rad));
- t_y=atan(-1*b/a/tan(theta_rad));
+ double t_r=theta_deg*M_PI/180;
+ double max_x, max_y, ct=cos(t_r), st=sin(t_r);
+ double t_x=atan(b/a*tan(t_r)), t_y=atan(-1.0f*b/a/tan(t_r));
- max_x=a*cos(t_x)*cos(theta_rad)+b*sin(t_x)*sin(theta_rad);
- max_y=-1*a*cos(t_y)*sin(theta_rad)+b*sin(t_y)*cos(theta_rad);
+ /* Calculate the maxima along each direction. */
+ max_x = a*cos(t_x)*ct + b*sin(t_x)*st;
+ max_y = -1*a*cos(t_y)*st + b*sin(t_y)*ct;
/* max_x and max_y are calculated from the center of the ellipse. We
want the final height and width of the box enclosing the
ellipse. So we have to multiply them by two, then take one from
- them (for the center).
- */
+ them (for the center). */
width[0]=2*( (size_t)fabs(max_x)+1 ) + 1;
width[1]=2*( (size_t)fabs(max_y)+1 ) + 1;
}
@@ -86,6 +87,151 @@ gal_box_ellipse_in_box(double a, double b, double
theta_rad, 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 ) */
+#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, j;
+ double sumsq, 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)
+ {
+ sumsq=0.0f; for(j=0;j<3;++j) sumsq += R[i*3+j] * R[i*3+j];
+ width[i] = 2 * ((long)sqrt(sumsq)+1) + 1;
+ }
+
+ /* 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("width: %ld, %ld, %ld\n", width[0], width[1], width[2]);
+ exit(0);
+ */
+}
+
+
+
+
+
/* We have the central pixel and box size of the crop box, find the
starting (first) and ending (last) pixels: */
void
@@ -129,19 +275,19 @@ gal_box_border_from_center(double *center, size_t ndim,
long *width,
-/* Problem to solve: We have set the first and last pixels in an input
- image (fpixel_i[2] and lpixel_i[2]). But those first and last
- pixels don't necessarily have to lie within the image, they can be
- outside of it or patially overlap with it. So the job of this
- function is to correct for such situations and find the starting
+/* Problem to solve: We have set the first and last pixels of a box in an
+ input image (fpixel_i[2] and lpixel_i[2]). But those first and last
+ pixels don't necessarily lie within the image's boundaries. They can be
+ outside of it or patially overlap with it (see examples below). The job
+ of this function is to corret for such situations and find the starting
and ending points of any overlap.
- It is assumed that your output (overlap) image's first pixel lies
- right ontop of the fpixel_i[0] in the input image. But since
- fpixel_i might be outside of the image, in this function we find
- the fpixel_o[2] and lpixel_o[2] in the overlap image coordinates
- that overlap with the input image. So the values of all four points
- might change after this function.
+ It is assumed that your output (overlap) image's first pixel lies right
+ ontop of the fpixel_i[0] in the input image. But since fpixel_i might be
+ outside of the image, in this function we find the fpixel_o[2] and
+ lpixel_o[2] in the overlap image coordinates that overlap with the input
+ image. So the values of all four points might change after this
+ function.
Before:
=======
diff --git a/lib/data.c b/lib/data.c
index 68178b2..ba7ff40 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -434,7 +434,7 @@ gal_data_free_contents(gal_data_t *data)
free(data->mmapname);
}
else
- if(data->array) free(data->array);
+ if(data->array && data->block==NULL) free(data->array);
data->array=NULL;
}
diff --git a/lib/gnuastro/box.h b/lib/gnuastro/box.h
index a1aa764..399d6a4 100644
--- a/lib/gnuastro/box.h
+++ b/lib/gnuastro/box.h
@@ -1,5 +1,5 @@
/*********************************************************************
-overlap -- Find the overlap of a region and an image.
+Box -- Define bounding and overlapping boxes.
This is part of GNU Astronomy Utilities (Gnuastro) package.
Original author:
@@ -52,7 +52,10 @@ __BEGIN_C_DECLS /* From C++ preparations */
*/
void
-gal_box_ellipse_in_box(double a, double b, double theta_rad, long *width);
+gal_box_bound_ellipse(double a, double b, double theta_deg, long *width);
+
+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,
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 0d17484..2809b23 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -71,6 +71,8 @@ gal_tile_series_from_minmax(gal_data_t *block, size_t
*minmax, size_t number);
+
+
/***********************************************************************/
/************** Allocated block **********************/
/***********************************************************************/
@@ -96,6 +98,8 @@ gal_tile_block_blank_flag(gal_data_t *tile_ll, size_t
numthreads);
+
+
/***********************************************************************/
/************** Tile full dataset ********************/
/***********************************************************************/
@@ -177,40 +181,90 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
`GAL_TILE_PARSE_OPERATE'), so some variables (basic definitions) that
are already defined in `GAL_TILE_PARSE_OPERATE' re-defined here. */
#define GAL_TILE_PO_OISET(IT, OT, IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \
+ IT *i=IN->array; \
+ gal_data_t *tpo_other=OTHER; /* `OTHER' may be NULL. */ \
+ gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \
+ \
+ size_t tpo_s_e_i_junk[2]={0,0}; \
+ IT b, *tpo_st=NULL, *tpo_f=i+IN->size; \
+ size_t tpo_i_increment=0, tpo_num_i_inc=1; \
+ size_t tpo_o_increment=0, tpo_num_o_inc=1; \
int tpo_parse_other=(OTHER && PARSE_OTHER); \
- size_t tpo_increment=0, tpo_num_increment=1; \
- gal_data_t *tpo_other_w=OTHER; /* `OTHER' may be NULL. */ \
gal_data_t *tpo_iblock = gal_tile_block(IN); \
- IT b, *tpo_st=NULL, *i=IN->array, *tpo_f=i+IN->size; \
- OT *tpo_ost=NULL, *o = tpo_other_w ? tpo_other_w->array : NULL; \
- gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \
+ OT *tpo_ost=NULL, *o = tpo_other ? tpo_other->array : NULL; \
int tpo_hasblank = CHECK_BLANK ? gal_blank_present(IN, 0) : 0; \
size_t tpo_s_e_i[2]={0,tpo_iblock->size-1}; /* -1: this is INCLUSIVE */ \
\
+ \
+ /* A small sanity check: if `OTHER' is given, and it is a block, */ \
+ /* then it must have the same size as `IN's block. On the other */ \
+ /* hand, when `OTHER' is a tile, its must have `IN's size. */ \
+ if( tpo_parse_other ) \
+ { \
+ if( OTHER==tpo_oblock ) /* `OTHER' is a block. */ \
+ { \
+ if(gal_data_dsize_is_different(tpo_iblock, tpo_oblock) ) \
+ { \
+ /* `error' function, is a GNU extension, see above. */ \
+ fprintf(stderr, "GAL_TILE_PO_OISET: when " \
+ "`PARSE_OTHER' is non-zero, the allocated " \
+ "block size of `IN' and `OTHER' must be " \
+ "equal, but they are not: %zu and %zu " \
+ "elements respectively\n", tpo_iblock->size, \
+ tpo_oblock->size); \
+ exit(EXIT_FAILURE); \
+ } \
+ } \
+ else \
+ if(gal_data_dsize_is_different(IN, OTHER) ) \
+ { \
+ /* The `error' function, is a GNU extension and this is */ \
+ /* a header, not a library which the user has to compile */ \
+ /* every time (on their own system). */ \
+ fprintf(stderr, "GAL_TILE_PO_OISET: when " \
+ "`PARSE_OTHER' is non-zero, the sizes of `IN' " \
+ "and `OTHER' must be equal (in all " \
+ "dimensions), but they are not: %zu and %zu " \
+ "elements respectively\n", IN->size, \
+ tpo_other->size); \
+ exit(EXIT_FAILURE); \
+ } \
+ } \
+ \
+ \
/* Write the blank value for the input type into `b'. */ \
gal_blank_write(&b, tpo_iblock->type); \
\
- /* If this is a tile, not a full block. */ \
+ \
+ /* If this is a tile, not a full block, then we need to set the */ \
+ /* starting pointers (`tpo_st' and `tpo_ost'). The latter needs */ \
+ /* special attention: if it is a block, then we will use the */ \
+ /* the same starting element as the input tile. If `OTHER' is a */ \
+ /* tile, then use its own starting position (recall that we have */ \
+ /* already made sure that `IN' and `OTHER' have the same size. */ \
if(IN!=tpo_iblock) \
{ \
tpo_st = gal_tile_start_end_ind_inclusive(IN, tpo_iblock, \
- tpo_s_e_i); \
+ tpo_s_e_i); \
if( tpo_parse_other ) \
- tpo_ost = ( (OT *)(tpo_oblock->array) \
- + ( tpo_st - (IT *)(tpo_iblock->array) ) ); \
+ tpo_ost = ( OTHER==tpo_oblock \
+ ? ( (OT *)(tpo_oblock->array) \
+ + ( tpo_st - (IT *)(tpo_iblock->array) ) ) \
+ : gal_tile_start_end_ind_inclusive(tpo_other, \
+ tpo_oblock, \
+ tpo_s_e_i_junk) ); \
} \
\
+ \
/* Go over contiguous patches of memory. */ \
- while( tpo_s_e_i[0] + tpo_increment <= tpo_s_e_i[1] ) \
+ while( tpo_s_e_i[0] + tpo_i_increment <= tpo_s_e_i[1] ) \
{ \
- \
- /* If we are on a tile, reset `i' and `f'. Also, if there */ \
- /* is more than one element in `OTHER', then set that. Note */ \
- /* that it is upto the caller to increment `o' in `OP'. */ \
+ /* If we are on a tile, reset `i' and `o'. */ \
if(IN!=tpo_iblock) \
{ \
- tpo_f = ( i = tpo_st + tpo_increment ) + IN->dsize[IN->ndim-1]; \
- if(tpo_parse_other) o = tpo_ost + tpo_increment; \
+ tpo_f = ( ( i = tpo_st + tpo_i_increment ) \
+ + IN->dsize[IN->ndim-1] ); \
+ if(tpo_parse_other) o = tpo_ost + tpo_o_increment; \
} \
\
/* Do the operation depending the nature of the blank value. */ \
@@ -232,17 +286,32 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
else \
do { {OP;} if(tpo_parse_other) ++o;} while(++i<tpo_f);\
\
- \
/* Set the incrementation. On a fully allocated iblock (when */ \
/* `IN==tpo_iblock'), we have already gone through the whole */ \
/* array, so we'll set the incrementation to the size of the */ \
- /* while block which will stop the `while' loop above. On a */ \
+ /* whole block. This will stop the `while' loop above. On a */ \
/* tile, we need to increment to the next contiguous patch */ \
/* of memory to continue parsing this tile. */ \
- tpo_increment += ( IN==tpo_iblock ? tpo_iblock->size \
- : gal_tile_block_increment(tpo_iblock, IN->dsize, \
- tpo_num_increment++, \
- NULL) ); \
+ tpo_i_increment += ( IN==tpo_iblock \
+ ? tpo_iblock->size \
+ : gal_tile_block_increment(tpo_iblock, \
+ IN->dsize, \
+ tpo_num_i_inc++, \
+ NULL) ); \
+ \
+ /* Similarly, increment the other array if necessary. Like */ \
+ /* the above, when `OTHER' is a full block, we'll just use */ \
+ /* the same increment as `IN'. Otherwise, when `OTHER' is a */ \
+ /* tile, calculate its increment based on its own block. */ \
+ if(tpo_parse_other) \
+ { \
+ if(OTHER==tpo_oblock) tpo_o_increment=tpo_i_increment; \
+ else \
+ tpo_o_increment += gal_tile_block_increment(tpo_oblock, \
+ tpo_other->dsize, \
+ tpo_num_o_inc++, \
+ NULL); \
+ } \
} \
\
/* This is done in case the caller doesn't need `o' to avoid */ \
@@ -302,23 +371,9 @@ gal_tile_full_free_contents(struct
gal_tile_two_layer_params *tl);
pixel/element of the input). See the documentation for more on this
macro and some examples. */
#define GAL_TILE_PARSE_OPERATE(IN, OTHER, PARSE_OTHER, CHECK_BLANK, OP) { \
- int tpo_parse_other=(OTHER && PARSE_OTHER); \
gal_data_t *tpo_iblock = gal_tile_block(IN); \
gal_data_t *tpo_oblock = OTHER ? gal_tile_block(OTHER) : NULL; \
\
- /* A small sanity check. */ \
- if( tpo_parse_other \
- && gal_data_dsize_is_different(tpo_iblock, tpo_oblock) ) \
- { \
- /* The `error' function, is a GNU extension. */ \
- fprintf(stderr, "GAL_TILE_PARSE_OPERATE: when `PARSE_OTHER' is "\
- "non-zero, the allocated block size of `IN' and " \
- "`OTHER' must be equal, but they are not: %zu and %zu " \
- "elements respectively)", tpo_iblock->size, \
- tpo_oblock->size); \
- exit(EXIT_FAILURE); \
- } \
- \
/* First set the OTHER type. */ \
if(OTHER) \
switch(tpo_oblock->type) \
diff --git a/lib/tile.c b/lib/tile.c
index 6ecbbad..5457375 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -348,9 +348,9 @@ size_t
gal_tile_block_increment(gal_data_t *block, size_t *tsize,
size_t num_increment, size_t *coord)
{
- size_t increment;
size_t n=block->ndim;
size_t *b=block->dsize, *t=tsize;
+ size_t increment=GAL_BLANK_SIZE_T;
if(n>3)
error(EXIT_FAILURE, 0, "%s: currently only implemented for at most 3 "
@@ -358,6 +358,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
switch(n)
{
+ /* A zero-dimensional dataset is not defined. */
case 0:
error(EXIT_FAILURE, 0, "%s: zero dimensional input is not acceptable",
__func__);
@@ -374,18 +375,17 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
if(coord) ++coord[0];
break;
- /* Higher dimensions. */
- default:
- if(num_increment % t[n-2])
+ /* 3D: The increment depends on which dimension we are reaching. */
+ case 3:
+ if(num_increment % t[1])
{
- increment=b[n-1];
- if(coord) ++coord[n-2];
+ increment = b[2];
+ if(coord) ++coord[1];
}
else
{
- increment=(b[n-2] * b[n-1]) - ( (t[n-2]-1) * b[n-1] );
- ++coord[n-3];
- if(coord) coord[n-2]=coord[n-1]=0;
+ increment=(b[1] * b[2]) - ( (t[1]-1) * b[2] );
+ if(coord) { ++coord[0]; coord[1]=coord[2]=0; }
}
break;
}
diff --git a/lib/wcs.c b/lib/wcs.c
index 14fdedc..1fe18da 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/tile.h>
#include <gnuastro/fits.h>
#include <gnuastro/dimension.h>
+#include <gnuastro/permutation.h>
@@ -431,44 +432,82 @@ gal_wcs_angular_distance_deg(double r1, double d1, double
r2, double d2)
-
+/* Return the pixel scale of the dataset in units of the WCS. */
/* Return the pixel scale of the dataset in units of the WCS. */
double *
gal_wcs_pixel_scale(struct wcsprm *wcs)
{
gsl_vector S;
+ double *a, *out;
+ int permute_set;
gsl_matrix A, V;
- size_t n=wcs->naxis;
- double *a, *v, *pixscale;
-
- /* Allocate space for the `v' and `pixscale' arrays. */
- errno=0; v=malloc(n*n*sizeof *v);
- if(v==NULL)
- error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `v'",
- __func__, n*n*sizeof *v);
- errno=0; pixscale=malloc(n*sizeof *pixscale);
- if(pixscale==NULL)
- error(EXIT_FAILURE, errno, "%s: allocating %zu bytes for `pixscale'",
- __func__, n*sizeof *pixscale);
-
- /* Write the full matrix into an array, irrespective of what type it was
- stored in the wcsprm structure (`PCi_j' style or `CDi_j' style). */
+ size_t i, j, n=wcs->naxis;
+ double *v=gal_data_malloc_array(GAL_TYPE_FLOAT64, n*n, __func__, "v");
+ size_t *permutation=gal_data_malloc_array(GAL_TYPE_SIZE_T, n, __func__,
+ "permutation");
+ gal_data_t *pixscale=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &n, NULL,
+ 0, -1, NULL, NULL, NULL);
+
+
+ /* Write the full WCS rotation matrix into an array, irrespective of what
+ style it was stored in the wcsprm structure (`PCi_j' style or `CDi_j'
+ style). */
a=gal_wcs_warp_matrix(wcs);
+
/* Fill in the necessary GSL vector and Matrix structures. */
- S.size=n; S.stride=1; S.data=pixscale;
+ S.size=n; S.stride=1; S.data=pixscale->array;
V.size1=n; V.size2=n; V.tda=n; V.data=v;
A.size1=n; A.size2=n; A.tda=n; A.data=a;
+
/* Run GSL's Singular Value Decomposition, using one-sided Jacobi
orthogonalization which computes the singular (scale) values to a
- higher relative accuracy.*/
+ higher relative accuracy. */
gsl_linalg_SV_decomp_jacobi(&A, &V, &S);
+
+ /* The raw pixel scale array produced from the singular value
+ decomposition above is ordered based on values, not the input. So when
+ the pixel scales in all the dimensions aren't the same (for example in
+ IFU datacubes), the order of the values in `pixelscale' will not
+ necessarily correspond to the input's dimensions.
+
+ To correct the order, we can use the `V' matrix to find the original
+ position of the pixel scale values and then use permutation to
+ re-order it correspondingly. This works when there is only one
+ non-zero element in each row of `V'. */
+ for(i=0;i<n;++i)
+ {
+ permute_set=0;
+ for(j=0;j<n;++j)
+ if(v[i*n+j])
+ {
+ /* Only works when each row only has one non-zero value. */
+ if(permute_set)
+ error(EXIT_FAILURE, 0, "%s: not able to find the proper "
+ "permutation for given rotation matrix", __func__);
+ else
+ {
+ permutation[i]=j;
+ permute_set=1;
+ }
+ }
+ }
+
+
+ /* Apply the permutation described above. */
+ gal_permutation_apply(pixscale, permutation);
+
+
/* Clean up and return. */
free(a);
free(v);
- return pixscale;
+ free(permutation);
+ out=pixscale->array;
+ pixscale->array=NULL;
+ gal_data_free(pixscale);
+ return out;
}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 4050a2b..f0ee158 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -127,10 +127,11 @@ if COND_MKNOISE
mknoise/addnoise.sh: warp/warp_scale.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/3dcat.sh
+ mkprof/3dcat.sh: prepconf.sh.log
mkprof/mosaic1.sh: prepconf.sh.log
mkprof/mosaic2.sh: prepconf.sh.log
mkprof/mosaic3.sh: prepconf.sh.log
@@ -224,10 +225,10 @@ TESTS = prepconf.sh lib/versioncpp.sh lib/multithread.sh
\
# Files to distribute along with the tests.
-EXTRA_DIST = $(TESTS) during-dev.sh buildprog/simpleio.c crop/cat.txt \
- mkprof/mkprofcat1.txt mkprof/ellipticalmasks.txt mkprof/clearcanvas.txt \
- mkprof/mkprofcat2.txt mkprof/mkprofcat3.txt mkprof/mkprofcat4.txt \
- mkprof/radeccat.txt table/table.txt
+EXTRA_DIST = $(TESTS) during-dev.sh buildprog/simpleio.c crop/cat.txt \
+ mkprof/3dcat.txt mkprof/mkprofcat1.txt mkprof/ellipticalmasks.txt \
+ mkprof/clearcanvas.txt mkprof/mkprofcat2.txt mkprof/mkprofcat3.txt \
+ mkprof/mkprofcat4.txt mkprof/radeccat.txt table/table.txt
diff --git a/tests/mkprof/3dcat.sh b/tests/mkprof/3dcat.sh
new file mode 100755
index 0000000..4d39abd
--- /dev/null
+++ b/tests/mkprof/3dcat.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 <akhlaghi@gnu.org>
+# Contributing author(s):
+#
+# Copying and distribution of this file, with or without modification,
+# are permitted in any medium without royalty provided the copyright
+# notice and this notice are preserved. This file is offered as-is,
+# without any warranty.
+
+
+# Preliminaries
+# =============
+#
+# Set the variables (The executable is in the build tree). Do the
+# basic checks to see if the executable is made or if the defaults
+# file exists (basicchecks.sh is in the source tree).
+prog=mkprof
+execname=../bin/$prog/ast$prog
+cat=$topsrc/tests/$prog/3dcat.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/prepconf.sh b/tests/prepconf.sh
index 0c9bafb..664cc0f 100755
--- a/tests/prepconf.sh
+++ b/tests/prepconf.sh
@@ -68,10 +68,13 @@ rm addedoptions.txt
# ---------------------------------
#
# Each utility's configuration file is copied in the `tests' directory for
-# easy readability.
+# easy readability. We are just using a wildcard to copy any configuration
+# file in each program's source. This is because some programs come with
+# some supplementary configuration files to help the users (like
+# MakeProfiles to operate in 3D).
for prog in arithmetic convertt convolve cosmiccal crop fits \
mkcatalog mknoise mkprof noisechisel statistics \
table warp
do
- cp $topsrc/bin/$prog/ast$prog.conf .gnuastro/ast$prog.conf
+ cp $topsrc/bin/$prog/ast*.conf .gnuastro/
done
- [gnuastro-commits] master e47f8db 024/113: Merged recent work in master, (continued)
- [gnuastro-commits] master e47f8db 024/113: Merged recent work in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8fa5ff1 026/113: Minor edit in book (part added in last commit), Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 40f0a56 013/113: Minor corrections to MakeProfiles continued, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 8372486 020/113: NoiseChisel's detection complete in 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6cc3d25 027/113: No 3D projections in function to inspect NoiseChisel outputs, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master e990860 029/113: Merged recent updates in master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6bdc5d6 030/113: Oversegmentation connectivity one less for 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master be3bfd8 033/113: NoiseChisel configuration file in 3D updated, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 3b6c15d 036/113: Merged with recent work in master, no conflicts, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master fb3660f 037/113: MakeCatalog works in 3D, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 73fdf0c 006/113: MakeProfiles builds 3D ellipsoids,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 03d73fe 011/113: Some minor corrections in code and comments, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 87ab805 014/113: Identifiers for integer constants in BZERO comparisons, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 1483201 009/113: Noised cube created in make check, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master bdeaba9 012/113: Corrected name of 3D catalog for tarball inclusion, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 6b53c5e 015/113: Bug fixes in master: commits 1ff1c25 and 540af65, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master adeb1c6 016/113: NoiseChisel's convolution step in 3D completed, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master ace4160 017/113: New --mcolisbrightness option for MakeProfiles, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master f06acc5 018/113: Merged recent work from master, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 61aec7b 019/113: Convolve name and dimensionality checks, Mohammad Akhlaghi, 2021/04/16
- [gnuastro-commits] master 4627ddf 021/113: NoiseChisel segmentation working in 3D, test added, Mohammad Akhlaghi, 2021/04/16