[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master d4eedae: Radial distance profile added to Make
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master d4eedae: Radial distance profile added to MakeProfiles |
|
Date: |
Fri, 21 Jul 2017 18:34:23 -0400 (EDT) |
branch: master
commit d4eedae062a012a509bb83b6eebb7316064fcf9b
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Radial distance profile added to MakeProfiles
A new radial distance profile is added to MakeProfiles. To make it more
clear, MakeProfile's source was also cleaned to be easier to read and
debug.
---
NEWS | 4 +
bin/mkprof/args.h | 108 ++++++-------
bin/mkprof/main.h | 17 +-
bin/mkprof/mkprof.c | 406 ++++++++++++++++++++++++++----------------------
bin/mkprof/mkprof.h | 18 +--
bin/mkprof/oneprofile.c | 341 +++++++++++++++++++++++++---------------
bin/mkprof/profiles.c | 56 ++++---
bin/mkprof/profiles.h | 23 +--
bin/mkprof/ui.c | 255 +++++++++++++++++-------------
doc/gnuastro.texi | 12 ++
10 files changed, 721 insertions(+), 519 deletions(-)
diff --git a/NEWS b/NEWS
index dac6478..43e68bf 100644
--- a/NEWS
+++ b/NEWS
@@ -18,6 +18,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
to specify the PC matrix, CUNIT and CTYPE world coordinate system
keywords of the output FITS file.
+ MakeProfiles: the new `distance' profile will save the radial distance of
+ each pixel. This may be used to define your own profiles that are not
+ currently supported in MakeProfiles.
+
** Removed features
** Changed features
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index aaea679..7f3a824 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -165,6 +165,20 @@ struct argp_option program_options[] =
ARGS_GROUP_PROFILES
},
{
+ "mode",
+ UI_KEY_MODE,
+ "STR",
+ 0,
+ "Mode of `--ccol': `img' or `wcs'.",
+ ARGS_GROUP_PROFILES,
+ &p->mode,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_parse_coordinate_mode
+ },
+ {
"numrandom",
UI_KEY_NUMRANDOM,
"INT",
@@ -178,6 +192,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 +231,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_PROFILES,
+ &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 +284,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 +322,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
- },
@@ -319,26 +346,12 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "mode",
- UI_KEY_MODE,
- "STR",
- 0,
- "Coordinate mode `img' or `wcs'.",
- GAL_OPTIONS_GROUP_INPUT,
- &p->mode,
- GAL_TYPE_STRING,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_MANDATORY,
- GAL_OPTIONS_NOT_SET,
- ui_parse_coordinate_mode
- },
- {
"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,
@@ -424,19 +437,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/main.h b/bin/mkprof/main.h
index 041ddeb..5f58583 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. */
@@ -165,7 +163,7 @@ struct mkprofparams
uint8_t *f; /* Profile function code. */
float *r; /* Radius of profile. */
float *n; /* Index of profile. */
- float *p; /* Position angle of profile */
+ float *p; /* Position angle of profile. */
float *q; /* Axis ratio of profile. */
float *m; /* Magnitude of profile. */
float *t; /* Truncation distance. */
@@ -180,6 +178,11 @@ struct mkprofparams
char *wcsheader; /* The WCS header information for main img. */
int wcsnkeyrec; /* The number of keywords in the WCS header.*/
char *mergedimgname; /* Name of merged image. */
+ int nwcs; /* for WCSLIB: no. coord. representations. */
+ struct wcsprm *wcs; /* WCS information for this dataset. */
+ size_t ndim; /* Number of dimensions (for `nomerged'). */
+ /* We can't put it in `out' because it is */
+ /* meaning ful there. */
};
#endif
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 5a46b58..2821bd6 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>
@@ -78,11 +79,17 @@ builtqueue_addempty(struct builtqueue **bq)
error(EXIT_FAILURE, 0, "%s: allocating %zu bytes for `tbq'",
__func__, sizeof *tbq);
- /* Initialize some of the values. */
- tbq->img=NULL;
- tbq->numaccu=0;
- tbq->accufrac=0.0f;
- tbq->indivcreated=0;
+ /* Initialize the values (same order as in structure definition). */
+ tbq->id = GAL_BLANK_SIZE_T;
+ tbq->ispsf = 0;
+ tbq->overlaps = 0;
+ tbq->image = NULL;
+ tbq->overlap_i = NULL;
+ tbq->overlap_m = NULL;
+ tbq->func = PROFILE_MAXIMUM_CODE;
+ tbq->indivcreated = 0;
+ tbq->numaccu = 0;
+ tbq->accufrac = 0.0f;
/* Set its next element to the input bq and re-set the input bq. */
tbq->next=*bq;
@@ -118,15 +125,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;
+ size_t i, ndim=p->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 +143,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 +197,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->ndim;
+ size_t start_indiv[2], start_mrg[2], dsize[2], 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)
+ {
+ /* 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 +367,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->ndim;
struct builtqueue *ibq, *fbq=NULL;
- pthread_cond_t *qready=&p->qready;
+ long fpixel_i[2], lpixel_i[2], fpixel_o[2], lpixel_o[2];
+
/* Make each profile that was specified for this thread. */
for(i=0; mkp->indexs[i]!=GAL_BLANK_SIZE_T; ++i)
@@ -264,94 +402,33 @@ mkprof_build(void *inparam)
if( p->f[id] == PROFILE_POINT )
mkp->width[0]=mkp->width[1]=1;
else
- gal_box_bound_ellipse(mkp->truncr, mkp->q*mkp->truncr, p->p[id],
- mkp->width);
-
+ gal_box_bound_ellipse(mkp->truncr, mkp->q[0]*mkp->truncr,
+ p->p[id], 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);
-
-
- /* Build the profile if necessary, After this, the width is
- oversampled. */
- if(ibq->overlaps || p->individual || (ibq->ispsf && p->psfinimg==0))
+ if(p->out)
{
- /* 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;
- }
+ center[0]=p->x[id];
+ center[1]=p->y[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);
}
- /* 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);
- }
- }
+ /* Build the profile if necessary. */
+ if(ibq->overlaps || p->individual || (ibq->ispsf && p->psfinimg==0))
+ mkprof_build_single(mkp, fpixel_i, lpixel_i, fpixel_o);
+
+
+ /* 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 +464,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 +494,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). */
- if(ibq->overlaps && out->array)
- {
- /* 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);
- }
+ /* 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)
+ 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;
+ });
+
/* Fill the log array. */
if(p->cp.log)
@@ -510,6 +532,7 @@ mkprof_write(struct mkprofparams *p)
}
}
+
/* Report if in verbose mode. */
++complete;
if(!p->cp.quiet && p->num>1)
@@ -520,25 +543,33 @@ 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;
}
+
/* Write the final array to the output FITS image if a merged image is to
be created. */
- if(out->array)
+ if(out)
{
/* Get the current time for verbose output. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
- /* Write the final image into a FITS file with the requested type. */
+ /* Write the final image into a FITS file with the requested
+ type. Until now, we were using `p->wcs' for the WCS, but from now
+ on, will put it in `out' to also free it while freeing `out'. */
+ out->wcs=p->wcs;
gal_fits_img_write_to_type(out, p->mergedimgname, NULL,
PROGRAM_STRING, p->cp.type);
+ p->wcs=NULL;
/* Clean up */
gal_data_free(out);
@@ -551,11 +582,6 @@ mkprof_write(struct mkprofparams *p)
free(jobname);
}
}
-
- /* Even with no merged image, there might still be pointers in `out' that
- need to be freed. */
- else
- gal_data_free(p->out);
}
@@ -590,10 +616,9 @@ mkprof(struct mkprofparams *p)
pthread_barrier_t b;
struct mkonthread *mkp;
gal_list_str_t *comments=NULL;
- long *onaxes, os=p->oversample;
size_t i, fi, *indexs, thrdcols;
- size_t nb, ndim=p->out->ndim, nt=p->cp.numthreads;
-
+ long *onaxes=NULL, os=p->oversample;
+ size_t nb, ndim=p->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
@@ -611,13 +636,18 @@ mkprof(struct mkprofparams *p)
gal_threads_dist_in_threads(p->num, nt, &indexs, &thrdcols);
- /* onaxes are sides of the image without over-sampling in FITS order. */
- onaxes=gal_data_malloc_array(GAL_TYPE_LONG, ndim, __func__, "onaxes");
- for(fi=0; fi < ndim; ++fi)
+ /* `onaxes' are size of the merged output image without over-sampling or
+ shifting in FITS order. When no output merged image is needed, we can
+ ignore it. */
+ if(p->out)
{
- i=ndim-fi-1;
- onaxes[fi] = ( ( p->dsize[i] - 2 * p->shift[i] ) / os
- + 2 * p->shift[i]/os );
+ onaxes=gal_data_malloc_array(GAL_TYPE_LONG, ndim, __func__, "onaxes");
+ for(fi=0; fi < ndim; ++fi)
+ {
+ i=ndim-fi-1;
+ onaxes[fi] = ( ( p->dsize[i] - 2 * p->shift[i] ) / os
+ + 2 * p->shift[i]/os );
+ }
}
@@ -694,5 +724,5 @@ mkprof(struct mkprofparams *p)
/* Clean up. */
free(mkp);
free(indexs);
- free(onaxes);
+ if(onaxes) free(onaxes);
}
diff --git a/bin/mkprof/mkprof.h b/bin/mkprof/mkprof.h
index b786d78..1f2d758 100644
--- a/bin/mkprof/mkprof.h
+++ b/bin/mkprof/mkprof.h
@@ -31,17 +31,13 @@ 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[2]; /* Pixel coordinate. */
+ double lower[2]; /* Coordinates of lower pixel position. */
+ double higher[2]; /* Coordinates of higher pixel position. */
+ double c[1]; /* Cosine of position angle(s). */
+ double s[1]; /* Sine of position angle(s). */
+ double q[1]; /* Axis ratio(s). */
+ double center[2]; /* Center (in FITS) in oversampled image.*/
double (*profile)(struct mkonthread *); /* Function to use. */
double truncr; /* Truncation radius in pixels. */
double intruncr; /* Inner truncation radius in pixels. */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index e36385a..d756e42 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -51,30 +51,100 @@ 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;
+ size_t i, id=mkp->ibq->id;
+ double val, pixfrac, intpart;
+
+ for(i=0;i<p->ndim;++i)
+ {
+ dim = i==0 ? p->x : p->y;
+
+ 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)
+{
+ size_t i, coord_c[2];
+ 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;
+}
+
+
+
+
+
+/* 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)
{
- double x, y;
+ double Xr, Yr; /* Rotated x and y. */
+ double q=mkp->q[0];
+ double c=mkp->c[0], s=mkp->s[0];
+ double x=mkp->coord[0], y=mkp->coord[1];
+
+ Xr = x * ( c ) + y * ( s );
+ Yr = x * ( -1.0f*s ) + y * ( c );
+ mkp->r = sqrt( Xr*Xr + Yr*Yr/q/q );
+}
+
+
+
- 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) );
+/* 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[2];
+ 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 +170,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[2], sum=0.0f;
+ size_t i, j, numrandom=mkp->p->numrandom, ndim=mkp->p->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 +213,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 +267,7 @@ integ2d(struct mkonthread *mkp)
epsrel, &result, &abserr, &neval);
return result;
}
-
-
+#endif
@@ -220,49 +290,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[2], ndim=mkp->p->ndim;
+
+ /* Find the coordinates of the center point. Note `mkp->center' is in
+ FITS coordinates, while coord must be in C coordinates (to be used in
+ `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 +367,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 +403,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 +422,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. */
@@ -395,7 +491,7 @@ oneprofile_ispsf(uint8_t fcode)
-/* About the shifts on the X column and y column:*/
+/* Prepare all the parameters for any type of profile. */
void
oneprof_set_prof_params(struct mkonthread *mkp)
{
@@ -405,26 +501,31 @@ oneprof_set_prof_params(struct mkonthread *mkp)
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];
+ /* 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. */
+ /* Shifts were already multiplied with oversample. Just note that
+ p->x and p->y are in the FITS ordering, while p->shift is in C
+ ordering. */
+ mkp->q[0] = p->q[id];
+ p->x[id] += p->shift[1]/p->oversample;
+ p->y[id] += p->shift[0]/p->oversample;
+ mkp->c[0] = cos( p->p[id] * DEGREESTORADIANS );
+ mkp->s[0] = sin( p->p[id] * DEGREESTORADIANS );
+
+
+ /* Fill 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 +533,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)
@@ -449,7 +550,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;
@@ -466,13 +567,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 +590,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,6 +609,12 @@ 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 %u is not recognized.", __func__,
@@ -542,44 +649,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->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;
+ /* Find the profile center in the over-sampled image in C
+ coordinates. IMPORTANT: width must not be oversampled.*/
+ oneprofile_center_oversampled(mkp);
- 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;
-
- /* 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 +685,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 +701,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/profiles.c b/bin/mkprof/profiles.c
index 6ee0197..39e38f4 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.*/
+
+/* 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 5e2714d..41fa790 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -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");
@@ -161,6 +164,7 @@ ui_profile_name_write(int profile_code)
case PROFILE_POINT: return "point";
case PROFILE_FLAT: return "flat";
case PROFILE_CIRCUMFERENCE: return "circum";
+ case PROFILE_DISTANCE: return "distance";
default:
error(EXIT_FAILURE, 0, "%s: %d not recognized as a profile code",
__func__, profile_code);
@@ -286,7 +290,7 @@ ui_parse_kernel(struct argp_option *option, char *arg,
long profcode;
double *darray;
gal_data_t *kernel;
- size_t i, nc, numneeded=0;
+ size_t i, nc, need=0;
char *c, *profile, *tailptr;
char *str, sstr[GAL_OPTIONS_STATIC_MEM_FOR_VALUES];
@@ -299,7 +303,8 @@ ui_parse_kernel(struct argp_option *option, char *arg,
/* First write the profile function code into the output string. */
nc=0;
- nc += sprintf(sstr+nc, "%s,", ui_profile_name_write(kernel->status));
+ profile=ui_profile_name_write(kernel->status);
+ nc += sprintf(sstr+nc, "%s,", profile);
/* Write the values into a string. */
for(i=0;i<kernel->size;++i)
@@ -326,13 +331,25 @@ ui_parse_kernel(struct argp_option *option, char *arg,
the rest.*/
c=arg;while(*c!='\0' && *c!=',') ++c;
profile=arg;
- arg = (*c=='\0') ? NULL : c+1; /* `point' doesn't need any numbers. */
- *c='\0';
+ arg = (*c=='\0') ? NULL : c+1; /* the `point' profile doesn't need */
+ *c='\0'; /* any numbers. */
+
/* Read the parameters. */
kernel=gal_options_parse_list_of_numbers(arg, filename, lineno);
*(gal_data_t **)(option->value) = kernel;
+
+ /* All parameters must be positive. */
+ darray=kernel->array;
+ for(i=0;i<kernel->size;++i)
+ if(darray[i]<=0)
+ error(EXIT_FAILURE, 0, "value number %zu (%g) in the given list "
+ "of kernel parameters (`%s') is not acceptable. All "
+ "parameters to the `--kernel' option must be non-zero and "
+ "positive", i+1, darray[i], arg);
+
+
/* Write the profile type code into `kernel->status'. If it starts
with a digit, then the user might have given the code of the
profile directly. In that case, parse the number. Otherwise,
@@ -353,15 +370,17 @@ ui_parse_kernel(struct argp_option *option, char *arg,
else
kernel->status=ui_profile_name_read(profile, 0);
+
/* Make sure the number of parameters conforms with the profile. */
switch(kernel->status)
{
- case PROFILE_SERSIC: numneeded=3; break;
- case PROFILE_MOFFAT: numneeded=3; break;
- case PROFILE_GAUSSIAN: numneeded=2; break;
- case PROFILE_POINT: numneeded=0; break;
- case PROFILE_FLAT: numneeded=1; break;
- case PROFILE_CIRCUMFERENCE: numneeded=1; break;
+ case PROFILE_SERSIC: need = 3; break;
+ case PROFILE_MOFFAT: need = 3; break;
+ case PROFILE_GAUSSIAN: need = 2; break;
+ case PROFILE_POINT: need = 0; break;
+ case PROFILE_FLAT: need = 1; break;
+ case PROFILE_CIRCUMFERENCE: need = 1; break;
+ case PROFILE_DISTANCE: need = 1; break;
default:
error_at_line(EXIT_FAILURE, 0, filename, lineno, "%s: a bug! "
"Please contact us at %s to correct the issue. "
@@ -369,13 +388,16 @@ ui_parse_kernel(struct argp_option *option, char *arg,
PACKAGE_BUGREPORT, kernel->status);
}
+
/* Make sure the number of parameters given are the same number that
are needed. */
- if( kernel->size != numneeded )
- error_at_line(EXIT_FAILURE, 0, filename, lineno, "as a kernel, a "
- "`%s' profile needs %zu parameters, but %zu is given",
- ui_profile_name_write(kernel->status), numneeded,
- kernel->size);
+ if( kernel->size != need )
+ error_at_line(EXIT_FAILURE, 0, filename, lineno, "as a kernel, "
+ "a `%s' profile needs %zu parameters, but %zu "
+ "parameter%s given to `--kernel'",
+ ui_profile_name_write(kernel->status), need,
+ kernel->size, kernel->size>1?"s are":" is");
+
/* Our job is done, return NULL. */
return NULL;
@@ -459,8 +481,10 @@ ui_parse_coordinate_mode(struct argp_option *option, char
*arg,
}
else
{
- if (!strcmp(arg, "img")) *(uint8_t
*)(option->value)=MKPROF_MODE_IMG;
- else if (!strcmp(arg, "wcs")) *(uint8_t
*)(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 "
@@ -626,23 +650,14 @@ static void
ui_read_cols(struct mkprofparams *p)
{
int checkblank;
+ size_t i, counter=0;
char *colname=NULL, **strarr;
gal_list_str_t *colstrs=NULL, *ccol;
gal_data_t *cols, *tmp, *corrtype=NULL;
- size_t i, counter=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)
+ for(i=0; i<p->ndim; ++i)
{
gal_list_str_add(&colstrs, ccol->v, 0);
ccol=ccol->next;
@@ -668,8 +683,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. */
@@ -679,7 +693,7 @@ ui_read_cols(struct mkprofparams *p)
turned off for some columns. */
checkblank=1;
- /* Put each column's data in the respective internal pointer. */
+ /* See which column we are currently reading. */
switch(++counter)
{
case 1:
@@ -695,6 +709,7 @@ ui_read_cols(struct mkprofparams *p)
}
break;
+
case 3:
if(tmp->type==GAL_TYPE_STRING)
{
@@ -716,7 +731,7 @@ ui_read_cols(struct mkprofparams *p)
/* 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 "
+ 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 "
@@ -729,6 +744,7 @@ ui_read_cols(struct mkprofparams *p)
}
break;
+
case 4:
colname="radius (`rcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
@@ -742,18 +758,21 @@ ui_read_cols(struct mkprofparams *p)
i+1, p->r[i]);
break;
+
case 5:
colname="index (`ncol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->n=corrtype->array;
break;
+
case 6:
colname="position angle (`pcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
p->p=corrtype->array;
break;
+
case 7:
colname="axis ratio (`qcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
@@ -765,9 +784,9 @@ ui_read_cols(struct mkprofparams *p)
error(EXIT_FAILURE, 0, "%s: row %zu, the axis ratio value %g "
"is not acceptable. It has to be >0 and <=1", p->catname,
i+1, p->q[i]);
-
break;
+
case 8:
colname="magnitude (`mcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
@@ -775,6 +794,7 @@ ui_read_cols(struct mkprofparams *p)
checkblank=0; /* Magnitude can be NaN: to mask regions. */
break;
+
case 9:
colname="truncation (`tcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
@@ -786,9 +806,9 @@ ui_read_cols(struct mkprofparams *p)
error(EXIT_FAILURE, 0, "%s: row %zu, the truncation radius "
"value %g is not acceptable. It has to be larger than 0",
p->catname, i+1, p->t[i]);
-
break;
+
/* If the index isn't recognized, then it is larger, showing that
there was more than one match for the given criteria */
default:
@@ -828,8 +848,8 @@ ui_read_cols(struct mkprofparams *p)
static void
ui_prepare_columns(struct mkprofparams *p)
{
- float r, n, t;
double *karr;
+ float r, n, t;
/* If the kernel option was called, then we need to build a series of
single element columns to create an internal catalog. */
@@ -839,20 +859,21 @@ 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");
-
- /* Set the values that need special consideration. */
+ p->x = gal_data_calloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->x");
+ p->y = gal_data_calloc_array(GAL_TYPE_FLOAT64, 1, __func__, "p->y");
+ p->f = gal_data_calloc_array(GAL_TYPE_UINT8, 1, __func__, "p->f");
+ p->r = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->r");
+ p->n = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->n");
+ p->p = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->p");
+ p->q = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->q");
+ p->m = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->m");
+ p->t = gal_data_calloc_array(GAL_TYPE_FLOAT32, 1, __func__, "p->t");
+
+ /* For profiles that need a different number of input values. Note
+ that when a profile doesn't need a value, it will be ignored. */
+ karr=p->kernel->array;
if(p->kernel->size)
{
- karr=p->kernel->array;
r = karr[0];
n = p->kernel->size==2 ? 0.0f : karr[1];
t = p->kernel->size==1 ? 1.0f : karr[ p->kernel->size - 1 ];
@@ -871,7 +892,19 @@ ui_prepare_columns(struct mkprofparams *p)
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) < p->ndim )
+ error(EXIT_FAILURE, 0, "%zu coordinate columns (calls to "
+ "`--coordcol') given but output has %zu dimensions",
+ gal_list_str_number(p->ccol), p->ndim);
+
+ /* Call the column-reading function. */
+ ui_read_cols(p);
+ }
}
@@ -885,7 +918,7 @@ ui_prepare_columns(struct mkprofparams *p)
static int
ui_wcs_sanity_check(struct mkprofparams *p)
{
- size_t ndim=p->out->ndim;
+ size_t ndim=p->ndim;
if(p->crpix)
{
@@ -953,14 +986,10 @@ ui_prepare_wcs(struct mkprofparams *p)
int status;
struct wcsprm *wcs;
char **cunit, **ctype;
- size_t i, ndim=p->out->ndim;
+ size_t i, ndim=p->ndim;
double *crpix, *crval, *cdelt, *pc;
- /* If the WCS structure is already set, then return. */
- if(p->out->wcs) return;
-
-
/* Check and initialize the WCS information. If any of the necessary WCS
parameters are missing, then don't build any WCS. */
if( ui_wcs_sanity_check(p) ) return;
@@ -974,7 +1003,7 @@ ui_prepare_wcs(struct mkprofparams *p)
/* Allocate the memory necessary for the wcsprm structure. */
errno=0;
- wcs=p->out->wcs=malloc(sizeof *wcs);
+ wcs=p->wcs=malloc(sizeof *wcs);
if(wcs==NULL)
error(EXIT_FAILURE, errno, "%zu for wcs in preparewcs", sizeof *wcs);
@@ -1004,7 +1033,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)
@@ -1021,9 +1049,10 @@ ui_prepare_canvas(struct mkprofparams *p)
{
float *f, *ff;
double truncr;
+ gal_data_t *keysll;
long width[2]={1,1};
int status=0, setshift=0;
- size_t i, ndim, nshift=0, *dsize=NULL;
+ size_t i, nshift=0, *dsize=NULL, ndim_counter;
/* If a background image is specified, then use that as the output
image to build the profiles over. */
@@ -1041,55 +1070,68 @@ ui_prepare_canvas(struct mkprofparams *p)
/* Read in the background image and its coordinates, note that when
no merged image is desired, we just need the WCS information of
- the background image. So `ndim==0' and what `dsize' points to is
- irrelevant. */
+ the background image and the number of its dimensions. So
+ `ndim==0' and what `dsize' points to is irrelevant. */
if(p->nomerged)
- p->out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 0, dsize, NULL,
- 1, p->cp.minmapsize, NULL, NULL, NULL);
+ {
+ /* Get the number of the background image's dimensions. */
+ keysll=gal_data_array_calloc(1);
+ keysll->name="NAXIS"; keysll->type=GAL_TYPE_SIZE_T;
+ gal_fits_key_read(p->backname, p->backhdu, keysll, 0, 0);
+ p->ndim = *(size_t *)(keysll->array);
+ keysll->name=NULL;
+ gal_data_array_free(keysll, 1, 1);
+ }
else
{
/* Read the image. */
p->out=gal_fits_img_read_to_type(p->backname, p->backhdu,
GAL_TYPE_FLOAT32,
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");
+ p->ndim=p->out->ndim;
/* If p->dsize was given as an option, free it. */
if( p->dsize ) free(p->dsize);
/* Write the size of the background image into `dsize'. */
- p->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->out->ndim,
- __func__, "p->dsize");
- for(i=0;i<p->out->ndim;++i) p->dsize[i] = p->out->dsize[i];
+ p->dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
+ "p->dsize");
+ for(i=0;i<p->ndim;++i) p->dsize[i] = p->out->dsize[i];
/* Set all pixels to zero if the user wanted a clear canvas. */
if(p->clearcanvas)
{ff=(f=p->out->array)+p->out->size; do *f++=0.0f; while(f<ff);}
}
+
+ /* Currently, things are only implemented for 2D. */
+ if(p->ndim!=2)
+ error(EXIT_FAILURE, 0, "%s (hdu %s) has %zu dimensions. Currently "
+ "only a 2 dimensional background image is acceptable",
+ p->backname, p->backhdu, p->ndim);
+
+
/* When a background image is specified, oversample must be 1 and
there is no shifts. */
p->oversample=1;
if(p->shift) free(p->shift);
- p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->out->ndim,
- __func__, "p->shift (1)");
+ p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
+ "p->shift (1)");
/* Read the WCS structure of the background image. */
- p->out->wcs=gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->out->nwcs);
-
+ p->wcs=gal_wcs_read(p->backname, p->backhdu, 0, 0, &p->nwcs);
}
else
{
/* 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);
+ ndim_counter=0;
+ for(i=0;p->dsize[i]!=GAL_BLANK_SIZE_T;++i) ++ndim_counter;
+ p->ndim=ndim_counter;
+ /* Currently, things are only implemented for 2D. */
+ if(p->ndim!=2)
+ error(EXIT_FAILURE, 0, "%zu numbers given to `--naxis', only 2 "
+ "values may be given", p->ndim);
/* If any of the shift elements are zero, the others should be too!*/
if(p->shift && p->shift[0] && p->shift[1])
@@ -1102,10 +1144,10 @@ ui_prepare_canvas(struct mkprofparams *p)
}
/* Make sure it has the same number of elements as naxis. */
- if(ndim!=nshift)
+ if(p->ndim!=nshift)
error(EXIT_FAILURE, 0, "%zu and %zu elements given to `--ndim' "
"and `--shift' respectively. These two numbers must be the "
- "same", ndim, nshift);
+ "same", p->ndim, nshift);
}
else
{
@@ -1137,7 +1179,7 @@ ui_prepare_canvas(struct mkprofparams *p)
shifts (from zero). So, we'll just free it and reset
it. */
if(p->shift) free(p->shift);
- p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->out->ndim,
+ p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim,
__func__, "p->shift (2)");
if(setshift)
{
@@ -1149,43 +1191,40 @@ ui_prepare_canvas(struct mkprofparams *p)
/* If shift has not been set until now, set it. */
if(p->shift==NULL)
- p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
+ p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim, __func__,
"p->shift (3)");
/* Prepare the sizes of the final merged image (if it is to be
made). Note that even if we don't want a merged image, we still
need its WCS structure. */
- if(p->nomerged)
- ndim=0;
- else
+ if(p->nomerged==0)
{
- ndim=0;
+ ndim_counter=0;
for(i=0;p->dsize[i]!=GAL_BLANK_SIZE_T;++i)
{
/* Count the number of dimensions. */
- ++ndim;
+ ++ndim_counter;
/* Correct dsize. */
p->dsize[i] = (p->dsize[i]*p->oversample) + (2*p->shift[i]);
}
dsize = p->dsize;
- }
- /* Make the output structure. */
- p->out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, ndim, dsize, NULL, 1,
- p->cp.minmapsize, NULL, NULL, NULL);
+ /* Make the output structure. */
+ p->out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, ndim_counter, dsize,
+ NULL, 1, p->cp.minmapsize, NULL, NULL, NULL);
+ }
}
/* Make the WCS structure of the output data structure if it has not
been set yet. */
- ui_prepare_wcs(p);
+ if(p->wcs==NULL)
+ ui_prepare_wcs(p);
- /* Set the name, comments and units of the output structure/file. Note
- that when no merged image is to be created, the array in `p->out' will
- be NULL.*/
- if(p->out->array)
+ /* Set the name, comments and units of the final merged output. */
+ if(p->out)
{
if(p->out->name) free(p->out->name);
gal_checkset_allocate_copy("Mock profiles", &p->out->name);
@@ -1199,9 +1238,9 @@ ui_prepare_canvas(struct mkprofparams *p)
string to speed up the process: if we don't do it here, this process
will be necessary on every individual profile's output. So it is much
more efficient done once here. */
- if(p->individual && p->out->wcs)
+ if(p->individual && p->wcs)
{
- status=wcshdo(WCSHDO_safe, p->out->wcs, &p->wcsnkeyrec, &p->wcsheader);
+ status=wcshdo(WCSHDO_safe, p->wcs, &p->wcsnkeyrec, &p->wcsheader);
if(status)
error(EXIT_FAILURE, 0, "wcshdo error %d: %s", status,
wcs_errmsg[status]);
@@ -1215,10 +1254,10 @@ ui_prepare_canvas(struct mkprofparams *p)
static void
ui_finalize_coordinates(struct mkprofparams *p)
{
+ size_t i, ndim=p->ndim;
double *x=NULL, *y=NULL;
uint8_t os=p->oversample;
- size_t i, ndim=p->out->ndim;
- double *cdelt=p->out->wcs->cdelt, *crpix=p->out->wcs->crpix;
+ double *cdelt=p->wcs->cdelt, *crpix=p->wcs->crpix;
/* When the user specified RA and Dec columns, the respective values
where stored in the `p->x' and `p->y' arrays. So before proceeding, we
@@ -1228,7 +1267,7 @@ ui_finalize_coordinates(struct mkprofparams *p)
/* Note that we read the RA and Dec columns into the `p->x' and `p->y'
arrays temporarily before. Here, we will convert them, free the old
ones and replace them with the proper X and Y values. */
- gal_wcs_world_to_img(p->out->wcs, p->x, p->y, &x, &y, p->num);
+ gal_wcs_world_to_img(p->wcs, p->x, p->y, &x, &y, p->num);
/* If any conversions created a WCSLIB error, both the outputs will be
set to NaN. */
@@ -1249,7 +1288,7 @@ ui_finalize_coordinates(struct mkprofparams *p)
background image, oversample is set to 1. This is done here because
the conversion of WCS to pixel coordinates needs to be done with the
non-over-sampled image.*/
- for(i=0;i<p->out->ndim;++i)
+ for(i=0;i<p->ndim;++i)
{
/* Oversampling has already been applied in `p->shift'. Also note
that shift is in the C dimension ordring, while crpix is in FITS
@@ -1322,19 +1361,25 @@ ui_preparations(struct mkprofparams *p)
over-written: */
if(p->kernel)
{
+ /* Set the necessary constants. */
+ p->ndim=2;
p->nomerged=1;
p->psfinimg=0;
p->individual=1;
- }
- /* Prepare the output canvas. */
- ui_prepare_canvas(p);
+ /* Set the shift array. */
+ p->shift=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->ndim,
+ __func__, "p->shift");
+ }
+ else
+ ui_prepare_canvas(p);
/* Read in all the columns. */
ui_prepare_columns(p);
/* Read the (possible) RA/Dec inputs into X and Y for the builder.*/
- ui_finalize_coordinates(p);
+ if(p->wcs)
+ ui_finalize_coordinates(p);
/* Allocate the random number generator: */
gsl_rng_env_setup();
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 62ede17..cbf3e87 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14860,6 +14860,18 @@ Circumference profile with address@hidden' or
address@hidden'. 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}).
address@hidden
+Radial distance profile with address@hidden' or address@hidden'. At the lowest
+level, each pixel only has an elliptical radial distance given the
+profile's shape and orentiation (see @ref{Defining an ellipse}). When this
+profile is chosen, the pixel's elliptical radial distance from the profile
+center is written as its value. For this profile, the value in the
+magnitude column (@option{--mcol}) will be ignored.
+
+You can use this for checks or as a first approximation to define your own
+higher-level radial function. In the latter case, just note that the
+central values are going to be incorrect (see @ref{Sampling from a
+function}).
@end itemize
@item --rcol=STR/INT
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master d4eedae: Radial distance profile added to MakeProfiles,
Mohammad Akhlaghi <=