[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 7e888ce 073/125: First implementation of MakeP
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 7e888ce 073/125: First implementation of MakeProfiles with gal_data_t |
|
Date: |
Sun, 23 Apr 2017 22:36:40 -0400 (EDT) |
branch: master
commit 7e888ce7e7a4c088e83165695957829a19c679a2
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
First implementation of MakeProfiles with gal_data_t
MakeProfiles now works with `gal_data_t', at least with the first set of
tests that were done with `make check'. `gal_data_t' can definitely be used
much more effectively than this implementation, but for now, this was what
the time constraints allowed. It will be improved in the future, and
possible bugs will be fixed.
Besides the major changes in MakeProfiles, the following side-note changes
were also made:
- The names for functions in different `.c' files of MakeProfiles have the
filename prefixed (like the libraries). This greatly helps in
readability.
- `gal_data_copy_to_new_type' now accounts for blank values in input
dataset and puts the output's blank value in those elements of the
output dataset. To do this, two macros were necessary, one for input
integer types and one for floats, because the `isnan' function couldn't
be used on integer types.
- The old `gal_options_initialize_numthreads' function was replaced with a
new `gal_threads_number()' function so library users who don't have
access to Gnulib can also easily/portably get the number of threads in
the host system.
- The new `gal_wcs_world_to_img' has replaced the old
`gal_wcs_radec_array_to_xy' function. It is suited for the new approach
of keeping columns as separate arrays, and also benefits from one call
to WCSLIB for all input coordinates, instead of the old function which
called WSCLIB for every input coordinate.
- The test catalogs for `make check' on MakeProfiles were updated to be
usable in the updated MakeProfiles.
---
bin/mkprof/main.c | 2 +-
bin/mkprof/main.h | 13 +-
bin/mkprof/mkprof.c | 173 +++++++--------
bin/mkprof/oneprofile.c | 93 ++++----
bin/mkprof/oneprofile.h | 8 +-
bin/mkprof/ui.c | 500 ++++++++++++++++++++++++++++++++++++++++----
bin/mkprof/ui.h | 2 +-
doc/gnuastro.texi | 135 +++++++-----
lib/data.c | 110 ++++++++--
lib/gnuastro/threads.h | 11 +-
lib/gnuastro/wcs.h | 4 +-
lib/options.c | 12 +-
lib/options.h | 3 -
lib/table.c | 4 +
lib/threads.c | 14 +-
lib/wcs.c | 74 +++++--
tests/mkprof/mkprofcat1.txt | 27 ++-
tests/mkprof/mkprofcat2.txt | 12 +-
tests/mkprof/mkprofcat3.txt | 12 +-
tests/mkprof/mkprofcat4.txt | 12 +-
tests/mkprof/radeccat.sh | 2 +-
tests/mkprof/radeccat.txt | 25 ++-
tmpfs-config-make | 4 +-
23 files changed, 915 insertions(+), 337 deletions(-)
diff --git a/bin/mkprof/main.c b/bin/mkprof/main.c
index a3a530a..58a2fd2 100644
--- a/bin/mkprof/main.c
+++ b/bin/mkprof/main.c
@@ -49,7 +49,7 @@ main (int argc, char *argv[])
mkprof(&p);
/* Free all non-freed allocations. */
- freeandreport(&p, &t1);
+ ui_free_report(&p, &t1);
/* Return successfully.*/
return EXIT_SUCCESS;
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index 2031a75..4114ff8 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -56,7 +56,7 @@ enum profile_types
PROFILE_FLAT, /* Flat profile. */
PROFILE_CIRCUMFERENCE, /* Circumference profile. */
- MAXIMUM_PROFILE_CODE, /* Just for a sanity check. */
+ PROFILE_MAXIMUM_CODE, /* Just for a sanity check. */
};
#define MINCIRCUMWIDTH 0.5f
@@ -137,16 +137,14 @@ struct mkprofparams
unsigned char mforflatpix; /* mcol is flat pixel value (f is 4 or 5). */
double crpix[2]; /* CRPIX FITS header keywords. */
double crval[2]; /* CRVAL FITS header keywords. */
- float resolution; /* PC1_1 and PC2_2 FITS header keywords. */
+ double resolution; /* For CDELTi FITS header keywords. */
/* Output */
+ gal_data_t *out; /* Output image. */
int type; /* User's desired output type. */
- char *basename; /* Merged image name with no directory. */
char *outdir; /* Output directory. */
- int anyblank; /* ==1: there are blanks in back. */
- int nwcs; /* Number of WCS. */
- struct wcsprm *wcs; /* WCSparam structure. */
+ char *basename; /* Merged image name with no directory. */
/* Processing parameters: */
@@ -162,9 +160,8 @@ struct mkprofparams
float *t; /* Truncation distance. */
gsl_rng *rng; /* Main instance of random number generator.*/
time_t rawtime; /* Starting time of the program. */
- gal_data_t *out; /* Output image. */
double *cat; /* Input catalog. */
- double *log; /* Log data to be printed. */
+ gal_data_t *log; /* Log data to be printed. */
struct builtqueue *bq; /* Top (last) elem of build queue. */
pthread_cond_t qready; /* bq is ready to be written. */
pthread_mutex_t qlock; /* Mutex lock to change builtq. */
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index edbe56e..3da39cb 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -213,17 +213,12 @@ saveindividual(struct mkonthread *mkp)
central pixel by one and completely ignore the fractional part.
*/
void *
-build(void *inparam)
+mkprof_build(void *inparam)
{
- printf("\n... build needs to be corrected ... \n");
- exit(0);
-#if 0
-
struct mkonthread *mkp=(struct mkonthread *)inparam;
struct mkprofparams *p=mkp->p;
- size_t i;
- double *cat;
+ size_t i, id;
int lockresult;
long lpixel_o[2];
pthread_mutex_t *qlock=&p->qlock;
@@ -241,26 +236,25 @@ build(void *inparam)
fbq->next=p->bq and then set p->bq to ibq.*/
builtqueue_addempty(&mkp->ibq);
ibq=mkp->ibq;
- ibq->id=mkp->indexs[i];
+ id=ibq->id=mkp->indexs[i];
if(fbq==NULL) fbq=ibq;
- cat=&p->cat[ibq->id*p->cs1];
/* Write the necessary parameters for this profile into mkp.*/
- setprofparams(mkp);
+ oneprof_set_prof_params(mkp);
/* Find the bounding box size (NOT oversampled). */
- if((int)cat[p->fcol]==PROFILE_POINT)
+ 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,
- cat[p->pcol]*DEGREESTORADIANS, mkp->width);
+ p->p[id]*DEGREESTORADIANS, mkp->width);
/* Get the overlapping pixels using the starting points (NOT
oversampled). */
- gal_box_border_from_center(cat[p->xcol], cat[p->ycol], mkp->width,
+ gal_box_border_from_center(p->x[id], p->y[id], mkp->width,
ibq->fpixel_i, ibq->lpixel_i);
mkp->fpixel_i[0]=ibq->fpixel_i[0];
mkp->fpixel_i[1]=ibq->fpixel_i[1];
@@ -283,7 +277,7 @@ build(void *inparam)
gsl_rng_set(mkp->rng, gal_timing_time_based_rng_seed());
/* Make the profile */
- makeoneprofile(mkp);
+ oneprofile_make(mkp);
if( p->individual || (ibq->ispsf && p->psfinimg==0))
{
saveindividual(mkp);
@@ -309,13 +303,12 @@ build(void *inparam)
p->bq=ibq;
/* If the list was empty when you locked the mutex, then
- either `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 `write`) or it just got the
+ 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. */
+ output). In either case, pthread_cond_signal will work. */
if(fbq->next==NULL)
pthread_cond_signal(qready);
pthread_mutex_unlock(qlock);
@@ -350,7 +343,6 @@ build(void *inparam)
pthread_barrier_wait(mkp->b);
return NULL;
-#endif
}
@@ -375,29 +367,38 @@ build(void *inparam)
/**************************************************************/
/************ The writer *************/
/**************************************************************/
-void
-writelog(struct mkprofparams *p)
+static void
+mkprof_write_log(struct mkprofparams *p)
{
- char comments[1000], gitdescribe[100], *gd=gal_git_describe();
+ char *comments, gitdescribe[100], *gd;
+
+ /* This function is only relevant if a log file was desired. */
+ if(!p->cp.log) return;
+ /* Get the Git description in the running folder. */
+ gd=gal_git_describe();
if(gd) sprintf(gitdescribe, " from %s,", gd);
else gitdescribe[0]='\0';
- sprintf(comments, "# Log file for "PROGRAM_STRING".\n"
- "# Created%s on %s"
- "# Column 0: Row number in catalog (starting from zero).\n"
- "# Column 1: Overlap magnitude with final image "
- "(zeropoint: %.3f).\n"
- "# Column 2: Number of Monte Carlo integration pixels.\n"
- "# Column 3: Fraction of brightness in Monte Carlo "
- "integrated pixels.\n"
- "# Column 4: An individual image was created.\n",
- ctime(&p->rawtime), gitdescribe, p->zeropoint);
-
- /*
- gal_txtarray_array_to_txt(p->log, p->cs0, LOGNUMCOLS, comments, int_cols,
- accu_cols, space, prec, 'f', LOGFILENAME);
- */
+ /* Make the comments. */
+ asprintf(&comments, "# %s\n# Created%s on %s# Zeropoint: %.4f",
+ PROGRAM_STRING, gitdescribe, ctime(&p->rawtime), p->zeropoint);
+
+ /* Write the log file to disk */
+ gal_table_write(p->log, comments, GAL_TABLE_FORMAT_TXT, LOGFILENAME,
+ p->cp.dontdelete);
+ free(comments);
+
+ /* In verbose mode, print the information. */
+ if(!p->cp.quiet)
+ {
+ asprintf(&comments, "%s created.", LOGFILENAME);
+ gal_timing_report(NULL, comments, 1);
+ free(comments);
+ }
+
+ /* Clean up. */
+ gal_data_free_ll(p->log);
}
@@ -405,34 +406,19 @@ writelog(struct mkprofparams *p)
void
-write(struct mkprofparams *p)
+mkprof_write(struct mkprofparams *p)
{
+ double sum;
char *jobname;
- double sum, *log;
struct timeval t1;
long os=p->oversample;
int replace=p->replace;
- gal_data_t *out, *towrite;
- size_t complete=0, num=p->num;
struct builtqueue *ibq=NULL, *tbq;
float *to, *from, *colend, *rowend;
+ size_t complete=0, num=p->num, clog;
+ gal_data_t *out=p->out, *log, *towrite;
size_t i, j, iw, jw, ii, jj, w=p->naxes[0], ow;
- /* Note that `naxes' is in the FITS standard, not C. */
- size_t dsize[2]={p->naxes[1], p->naxes[0]};
-
- /* The `out' data structure is the canvas on which all the profiles will
- be built. When there is no background image specified, this should be
- a cleared (all zeros) image. But the user might want to build the
- profiles on a canvas (of real data or other mock images). If so, the
- `p->out' was read-in/allocated in `ui.c'. */
- if(p->out)
- out=p->out;
- else
- out=p->out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT, 2, dsize, p->wcs, 1,
- p->cp.minmapsize, "Combined mock image", "Brightness",
- NULL);
-
/* Write each image into the output array. */
while(complete<p->num)
@@ -458,8 +444,9 @@ write(struct mkprofparams *p)
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 && p->nomerged==0)
+ 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);
@@ -484,7 +471,7 @@ write(struct mkprofparams *p)
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=out->array+i*w+j;
+ to = (float *)(out->array) + i*w+j;
from=ibq->img+ii*ow+jj;
rowend=to+iw*w;
do
@@ -512,39 +499,53 @@ write(struct mkprofparams *p)
}
/* Fill the log array. */
- log=&p->log[ibq->id*LOGNUMCOLS];
- log[0] = ibq->id;
- log[1] = sum>0.0f ? -2.5f*log10(sum)+p->zeropoint : NAN;
- log[2] = ibq->numaccu;
- log[3] = ibq->accufrac;
- log[4] = ibq->indivcreated;
+ if(p->cp.log)
+ {
+ clog=0;
+ for(log=p->log; log!=NULL; log=log->next)
+ switch(++clog)
+ {
+ case 5:
+ ((unsigned char *)(log->array))[ibq->id] = ibq->indivcreated;
+ break;
+ case 4:
+ ((float *)(log->array))[ibq->id] = ibq->accufrac;
+ break;
+ case 3:
+ ((unsigned long *)(log->array))[ibq->id]=ibq->numaccu;
+ break;
+ case 2:
+ ((float *)(log->array))[ibq->id] =
+ sum>0.0f ? -2.5f*log10(sum)+p->zeropoint : NAN;
+ break;
+ case 1:
+ ((unsigned long *)(log->array))[ibq->id]=ibq->id+1;
+ break;
+ }
+ }
/* Report if in verbose mode. */
++complete;
- if(!p->cp.quiet && p->nomerged==0)
+ if(!p->cp.quiet)
{
- errno=0;
- jobname=malloc(100*sizeof *jobname);
- if(jobname==NULL)
- error(EXIT_FAILURE, errno, "jobname in mkprof.c");
- sprintf(jobname, "row %zu complete, %zu left to go",
+ asprintf(&jobname, "row %zu complete, %zu left to go",
ibq->id, num-complete);
gal_timing_report(NULL, jobname, 2);
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 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);
tbq=ibq->next;
free(ibq);
ibq=tbq;
}
- /* Write the final array to the output FITS image. */
- if(p->nomerged==0)
+ /* Write the final array to the output FITS image if a merged image is to
+ be created. */
+ if(out->array)
{
/* Get the current time for verbose output. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
@@ -572,6 +573,11 @@ 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);
}
@@ -634,7 +640,7 @@ mkprof(struct mkprofparams *p)
mkp[0].onaxes=onaxes;
mkp[0].indexs=indexs;
mkp[0].rng=gsl_rng_clone(p->rng);
- build(&mkp[0]);
+ mkprof_build(&mkp[0]);
}
else
{
@@ -661,16 +667,17 @@ mkprof(struct mkprofparams *p)
mkp[i].onaxes=onaxes;
mkp[i].rng=gsl_rng_clone(p->rng);
mkp[i].indexs=&indexs[i*thrdcols];
- err=pthread_create(&t, &attr, build, &mkp[i]);
+ err=pthread_create(&t, &attr, mkprof_build, &mkp[i]);
if(err)
error(EXIT_FAILURE, 0, "can't create thread %zu", i);
}
}
/* Write the created arrays into the image. */
- write(p);
- if(p->cp.log)
- writelog(p);
+ mkprof_write(p);
+
+ /* Write the log file. */
+ mkprof_write_log(p);
/* If numthreads>1, then wait for all the jobs to finish and destroy
the attribute and barrier. */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 103bd87..33c3871 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -222,7 +222,7 @@ integ2d(struct mkonthread *mkp)
/************ Pixel by pixel building *************/
/********* Positions are in C not FITS *********/
/**************************************************************/
-void
+static void
makepixbypix(struct mkonthread *mkp)
{
float circ_r;
@@ -236,7 +236,7 @@ makepixbypix(struct mkonthread *mkp)
double tolerance=mkp->p->tolerance, pixfrac, junk;
double (*profile)(struct mkonthread *)=mkp->profile;
double xc=mkp->xc, yc=mkp->yc, os=mkp->p->oversample;
- double truncr=mkp->truncr, approx, hp=mkp->p->halfpixel;
+ double truncr=mkp->truncr, approx, hp=0.5f/mkp->p->oversample;
/* lQ: Largest. sQ: Smallest in queue */
struct gal_linkedlist_tosll *lQ=NULL, *sQ;
@@ -404,12 +404,9 @@ makepixbypix(struct mkonthread *mkp)
/************ Set profile parameters *************/
/**************************************************************/
int
-ispsf(double fcolvalue)
+oneprofile_ispsf(int fcode)
{
- int f=fcolvalue;
- if(f==PROFILE_MOFFAT || f==PROFILE_GAUSSIAN)
- return 1;
- else return 0;
+ return fcode==PROFILE_MOFFAT || fcode==PROFILE_GAUSSIAN;
}
@@ -418,28 +415,23 @@ ispsf(double fcolvalue)
/* About the shifts on the X column and y column:*/
void
-setprofparams(struct mkonthread *mkp)
+oneprof_set_prof_params(struct mkonthread *mkp)
{
- printf("\n ... setprofparams needs to be corrected ...\n");
- exit(0);
-
-#if 0
struct mkprofparams *p=mkp->p;
- double *cat, sigma;
+ double sigma;
int tp=p->tunitinp;
- size_t rcol=p->rcol, tcol=p->tcol;
+ size_t id=mkp->ibq->id;
/* Fill in the profile independant parameters. */
- cat=&p->cat[mkp->ibq->id*p->cs1];
- cat[p->xcol] += p->shift[0]/p->oversample;
- cat[p->ycol] += p->shift[1]/p->oversample;
- mkp->c = cos((90-cat[p->pcol])*DEGREESTORADIANS);
- mkp->s = sin((90-cat[p->pcol])*DEGREESTORADIANS);
- mkp->q = cat[p->qcol];
- mkp->brightness = pow( 10, (p->zeropoint - cat[p->mcol]) / 2.5f );
- mkp->ibq->ispsf = ispsf(cat[p->fcol]);
- mkp->func = mkp->ibq->func=cat[p->fcol];
+ p->x[id] += p->shift[0]/p->oversample; /* Shifts were multiplied by */
+ p->y[id] += p->shift[1]/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];
+ mkp->brightness = pow( 10, (p->zeropoint - p->m[id]) / 2.5f );
+ mkp->ibq->ispsf = oneprofile_ispsf(p->f[id]);
+ mkp->func = mkp->ibq->func = p->f[id];
/* Fill the profile dependent parameters. */
@@ -448,10 +440,10 @@ setprofparams(struct mkonthread *mkp)
case PROFILE_SERSIC:
mkp->correction = 1;
mkp->profile = &Sersic;
- mkp->sersic_re = cat[rcol];
- mkp->sersic_inv_n = 1.0f/cat[p->ncol];
- mkp->sersic_nb = -1.0f*sersic_b(cat[p->ncol]);
- mkp->truncr = tp ? cat[tcol] : cat[tcol]*cat[rcol];
+ 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->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
break;
@@ -459,15 +451,15 @@ setprofparams(struct mkonthread *mkp)
case PROFILE_MOFFAT:
mkp->correction = 1;
mkp->profile = &Moffat;
- mkp->moffat_nb = -1.0f*cat[p->ncol];
- mkp->moffat_alphasq = moffat_alpha(cat[rcol], cat[p->ncol]);
+ mkp->moffat_nb = -1.0f*p->n[id];
+ mkp->moffat_alphasq = moffat_alpha(p->r[id], p->n[id]);
mkp->moffat_alphasq *= mkp->moffat_alphasq;
- mkp->truncr = tp ? cat[tcol] : cat[tcol]*cat[rcol]/2;
+ mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id]/2;
if(p->psfinimg==0 && p->individual==0)
{
mkp->brightness = 1.0f; /* When the PSF is a separate image, */
- cat[p->xcol] = 0.0f; /* it should be centered and have a */
- cat[p->ycol] = 0.0f; /* total brightness of 1.0f. */
+ p->x[id] = 0.0f; /* it should be centered and have a */
+ p->y[id] = 0.0f; /* total brightness of 1.0f. */
}
break;
@@ -476,14 +468,14 @@ setprofparams(struct mkonthread *mkp)
case PROFILE_GAUSSIAN:
mkp->correction = 1;
mkp->profile = &Gaussian;
- sigma = cat[rcol]/2.35482f;
+ sigma = p->r[id]/2.35482f;
mkp->gaussian_c = -1.0f/(2.0f*sigma*sigma);
- mkp->truncr = tp ? cat[tcol] : cat[tcol]*cat[rcol]/2;
+ mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id]/2;
if(p->psfinimg==0 && p->individual==0)
{
mkp->brightness = 1.0f; /* Same as the explanations for */
- cat[p->xcol] = 0.0f; /* The Moffat profile. */
- cat[p->ycol] = 0.0f;
+ p->x[id] = 0.0f; /* The Moffat profile. */
+ p->y[id] = 0.0f;
}
break;
@@ -499,11 +491,11 @@ setprofparams(struct mkonthread *mkp)
case PROFILE_FLAT:
mkp->profile = &Flat;
- mkp->truncr = tp ? cat[tcol] : cat[tcol]*cat[rcol];
+ mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
if(p->mforflatpix)
{
mkp->correction = 0;
- mkp->fixedvalue = cat[p->mcol];
+ mkp->fixedvalue = p->m[id];
}
else
{
@@ -516,12 +508,12 @@ setprofparams(struct mkonthread *mkp)
case PROFILE_CIRCUMFERENCE:
mkp->profile = &Circumference;
- mkp->truncr = tp ? cat[tcol] : cat[tcol]*cat[rcol];
+ mkp->truncr = tp ? p->t[id] : p->t[id]*p->r[id];
mkp->intruncr = mkp->truncr - p->circumwidth;
if(p->mforflatpix)
{
mkp->correction = 0;
- mkp->fixedvalue = cat[p->mcol];
+ mkp->fixedvalue = p->m[id];
}
else
{
@@ -540,7 +532,6 @@ setprofparams(struct mkonthread *mkp)
"seen and reported prior to this step. Please contact us so "
"we can correct this");
}
-#endif
}
@@ -566,29 +557,24 @@ setprofparams(struct mkonthread *mkp)
/************ Outside functions *************/
/**************************************************************/
void
-makeoneprofile(struct mkonthread *mkp)
+oneprofile_make(struct mkonthread *mkp)
{
- printf("\n ... makeoneprofile needs to be corrected ...\n");
- exit(0);
-
-#if 0
struct mkprofparams *p=mkp->p;
float sum;
- size_t size;
long os=p->oversample;
double pixfrac, intpart;
- double *cat=&p->cat[ mkp->ibq->id*p->cs1 ];
+ size_t size, id=mkp->ibq->id;
- /* Find the profile center (see comments above). mkp->width
- is in the non-oversampled scale.*/
- pixfrac = modf(fabs(cat[p->xcol]), &intpart);
+ /* 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(cat[p->ycol]), &intpart);
+ 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;
@@ -607,7 +593,7 @@ makeoneprofile(struct mkonthread *mkp)
mkp->ibq->img=calloc(size, sizeof *mkp->ibq->img);
if(mkp->ibq->img==NULL)
error(EXIT_FAILURE, 0, "%zu bytes for object in row %zu of data in %s",
- size*sizeof *mkp->ibq->img, mkp->ibq->id, mkp->p->catname);
+ size*sizeof *mkp->ibq->img, mkp->ibq->id, p->catname);
/* Build the profile in the image. */
@@ -635,5 +621,4 @@ makeoneprofile(struct mkonthread *mkp)
gal_array_fmultip_const(mkp->ibq->img, size,
mkp->brightness/sum);
}
-#endif
}
diff --git a/bin/mkprof/oneprofile.h b/bin/mkprof/oneprofile.h
index f72d591..e04f85b 100644
--- a/bin/mkprof/oneprofile.h
+++ b/bin/mkprof/oneprofile.h
@@ -23,13 +23,15 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef ONEPROFILE_H
#define ONEPROFILE_H
+#include "mkprof.h"
+
int
-ispsf(double fcolvalue);
+oneprofile_ispsf(int fcolvalue);
void
-setprofparams(struct mkonthread *mkp);
+oneprof_set_prof_params(struct mkonthread *mkp);
void
-makeoneprofile(struct mkonthread *mkp);
+oneprofile_make(struct mkonthread *mkp);
#endif
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 8540475..2bcb4ca 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -28,6 +28,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <string.h>
+#include <gnuastro/wcs.h>
+#include <gnuastro/box.h>
#include <gnuastro/fits.h>
#include <gnuastro/table.h>
#include <gnuastro/linkedlist.h>
@@ -40,6 +42,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "main.h"
#include "ui.h"
+#include "oneprofile.h"
#include "authors-cite.h"
@@ -176,6 +179,9 @@ ui_initialize_options(struct mkprofparams *p,
cp->program_authors = PROGRAM_AUTHORS;
cp->coptions = gal_commonopts_options;
+ /* Default program parameters. */
+ p->type=GAL_DATA_TYPE_FLOAT;
+
/* Modify the common options for this program. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
@@ -202,14 +208,7 @@ ui_initialize_options(struct mkprofparams *p,
/* Read the number of threads available to the user, this should be done
before reading command-line and configuration file options, since they
can change it. */
- gal_options_initialize_numthreads(cp);
-
-
- /* Set the non-zero initial values, the structure was initialized to have
- a zero/NULL value for all elements. So, this is necessary only when
- `0' is meaningful in the context of the variable. */
- p->type=GAL_DATA_TYPE_INVALID;
- p->crpix[0]=p->crpix[1]=p->crval[0]=p->crval[1]=NAN;
+ cp->numthreads=gal_threads_number();
}
@@ -304,6 +303,10 @@ ui_read_check_only_options(struct mkprofparams *p)
p->racol?"racol":"deccol", p->racol?"Dec":"RA",
p->xcol?"deccol":"racol");
+ /* If a log file is to be created, make sure that one doesn't already
+ exist. */
+ if(p->cp.log)
+ gal_checkset_check_remove_file(LOGFILENAME, p->cp.dontdelete);
}
@@ -315,6 +318,9 @@ ui_read_check_only_options(struct mkprofparams *p)
static void
ui_check_options_and_arguments(struct mkprofparams *p)
{
+ int d0f1;
+ char *tmpname;
+
/* Make sure an input catalog is given, and if it is FITS, that the HDU
is also provided. */
if(p->catname)
@@ -329,6 +335,31 @@ ui_check_options_and_arguments(struct mkprofparams *p)
"you need to give a catalog/table containing the information of "
"the profiles");
}
+
+
+ /* If cp->output was not specified on the command line or in any of
+ the configuration files, then automatic output should be used, in
+ which case, cp->output should be the current directory. */
+ if(p->cp.output==NULL)
+ gal_checkset_allocate_copy("./", &p->cp.output);
+
+
+ /* Set the necessary output names. */
+ d0f1=gal_checkset_dir_0_file_1(p->cp.output, p->cp.dontdelete);
+ if(d0f1) /* --output is a file name. */
+ {
+ p->mergedimgname=p->cp.output;
+ p->outdir=gal_checkset_dir_part(p->mergedimgname);
+ }
+ else /* --output is a directory name. */
+ {
+ gal_checkset_allocate_copy(p->cp.output, &p->outdir);
+ gal_checkset_check_dir_write_add_slash(&p->outdir);
+ tmpname=gal_checkset_automatic_output(&p->cp, p->catname, ".fits");
+ p->mergedimgname=gal_checkset_malloc_cat(p->outdir, tmpname);
+ free(tmpname);
+ }
+ p->basename=gal_checkset_not_dir_part(p->mergedimgname);
}
@@ -363,16 +394,25 @@ ui_read_profile_function(struct mkprofparams *p, char
**strarr)
{
if( !strcmp("sersic", strarr[i]) )
p->f[i]=PROFILE_SERSIC;
+
else if ( !strcmp("moffat", strarr[i]) )
p->f[i]=PROFILE_MOFFAT;
+
else if ( !strcmp("gaussian", strarr[i]) )
p->f[i]=PROFILE_GAUSSIAN;
+
else if ( !strcmp("point", strarr[i]) )
p->f[i]=PROFILE_POINT;
+
else if ( !strcmp("flat", strarr[i]) )
p->f[i]=PROFILE_FLAT;
+
else if ( !strcmp("circumference", strarr[i]) )
p->f[i]=PROFILE_CIRCUMFERENCE;
+
+ else if ( !strcmp(GAL_DATA_BLANK_STRING, strarr[i]) )
+ error(EXIT_FAILURE, 0, "profile function column has blank values. "
+ "Input columns cannot contain blank values");
else
error(EXIT_FAILURE, 0, "`%s' not recognized as a profile function "
"name in row %zu", strarr[i], i);
@@ -386,7 +426,8 @@ ui_read_profile_function(struct mkprofparams *p, char
**strarr)
static void
ui_read_cols(struct mkprofparams *p)
{
- size_t counter=0;
+ char *colname;
+ size_t counter=0, i;
gal_data_t *cols, *tmp, *corrtype;
char *ax1col=p->racol?p->racol:p->xcol;
char *ax2col=p->deccol?p->deccol:p->ycol;
@@ -420,68 +461,88 @@ ui_read_cols(struct mkprofparams *p)
/* Pop out the top node. */
tmp=gal_data_pop_from_ll(&cols);
+
/* 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 9:
+ colname="first axis position";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_DOUBLE);
p->x=corrtype->array;
- corrtype->array=NULL;
break;
case 8:
+ colname="second axis position";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_DOUBLE);
p->y=corrtype->array;
- corrtype->array=NULL;
break;
case 7:
if(tmp->type==GAL_DATA_TYPE_STRING)
- ui_read_profile_function(p, tmp->array);
+ {
+ ui_read_profile_function(p, tmp->array);
+ 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_DATA_TYPE_INT);
p->f=corrtype->array;
- corrtype->array=NULL;
+
+ /* 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 "
+ "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 6:
+ colname="radius (`rcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->r=corrtype->array;
- corrtype->array=NULL;
break;
case 5:
+ colname="index (`ncol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->n=corrtype->array;
- corrtype->array=NULL;
break;
case 4:
+ colname="position angle (`pcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->p=corrtype->array;
- corrtype->array=NULL;
break;
case 3:
+ colname="axis ratio (`qcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->q=corrtype->array;
- corrtype->array=NULL;
break;
case 2:
+ colname="magnitude (`mcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->m=corrtype->array;
- corrtype->array=NULL;
break;
case 1:
+ colname="truncation (`tcol')";
corrtype=gal_data_copy_to_new_type_free(tmp, GAL_DATA_TYPE_FLOAT);
p->t=corrtype->array;
- corrtype->array=NULL;
break;
/* If the index isn't recognized, then it is larger, showing that
@@ -500,6 +561,24 @@ ui_read_cols(struct mkprofparams *p)
p->catname, p->catname);
}
+ /* 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( gal_data_has_blank(corrtype) )
+ 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;
+ free(corrtype);
+ }
}
}
@@ -508,17 +587,330 @@ ui_read_cols(struct mkprofparams *p)
static void
-ui_preparations(struct mkprofparams *p)
+ui_prepare_wcs(struct mkprofparams *p)
+{
+ int status;
+ struct wcsprm *wcs;
+ long os=p->oversample;
+
+ /* If the WCS structure is already set, then return. */
+ if(p->out->wcs) return;
+
+ /* Allocate the memory necessary for the wcsprm structure. */
+ errno=0;
+ wcs=p->out->wcs=malloc(sizeof *wcs);
+ if(wcs==NULL)
+ error(EXIT_FAILURE, errno, "%zu for wcs in preparewcs", sizeof *wcs);
+
+ /* Initialize the structure (allocate all its internal arrays). */
+ wcs->flag=-1;
+ if( (status=wcsini(1, 2, wcs)) )
+ error(EXIT_FAILURE, 0, "wcsinit error %d: %s",
+ status, wcs_errmsg[status]);
+
+ /* Correct the CRPIX values based on oversampling and shifting. */
+ p->crpix[0] = p->crpix[0]*os + p->shift[0] - os/2;
+ p->crpix[1] = p->crpix[1]*os + p->shift[1] - os/2;
+
+ /* Fill in all the important WCS structure parameters. */
+ wcs->equinox=2000.0f;
+ wcs->crpix[0]=p->crpix[0];
+ wcs->crpix[1]=p->crpix[1];
+ wcs->crval[0]=p->crval[0];
+ wcs->crval[1]=p->crval[1];
+ wcs->pc[0]=-1.0f;
+ wcs->pc[3]=1.0f;
+ wcs->pc[1]=wcs->pc[2]=0.0f;
+ wcs->cdelt[0]=wcs->cdelt[1]=p->resolution/3600;
+ strcpy(wcs->cunit[0], "deg");
+ strcpy(wcs->cunit[1], "deg");
+ strcpy(wcs->ctype[0], "RA---TAN");
+ strcpy(wcs->ctype[1], "DEC--TAN");
+
+ /* Set up the wcs structure with the constants defined above. */
+ status=wcsset(wcs);
+ if(status)
+ error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
+ wcs_errmsg[status]);
+}
+
+
+
+
+
+static void
+ui_prepare_canvas(struct mkprofparams *p)
+{
+ int status=0;
+ double truncr;
+ long width[2]={1,1};
+ size_t i, ndim, dsize[2];
+
+ /* If a background image is specified, then use that as the output
+ image to build the profiles over. */
+ if(p->backname)
+ {
+ /* Small sanity check. */
+ if(p->backhdu==NULL)
+ error(EXIT_FAILURE, 0, "no hdu specified for the background image "
+ "%s. Please run again `--backhdu' option", p->backname);
+
+ /* 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. */
+ if(p->nomerged)
+ p->out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT, 0, dsize, NULL, 1,
+ p->cp.minmapsize, NULL, NULL, NULL);
+ else
+ {
+ p->out=gal_fits_read_to_type(p->backname, p->backhdu,
+ GAL_DATA_TYPE_FLOAT, p->cp.minmapsize);
+ p->naxes[0]=p->out->dsize[1];
+ p->naxes[1]=p->out->dsize[0];
+ }
+
+ /* When a background image is specified, oversample must be 1 and
+ there is no shifts. */
+ p->oversample=1;
+ p->shift[0]=p->shift[1]=0;
+
+ /* Read the WCS structure of the background image. */
+ gal_fits_read_wcs(p->backname, p->backhdu, 0, 0, &p->out->nwcs,
+ &p->out->wcs);
+
+ }
+ else
+ {
+
+ /* If any of xshift or yshift is non-zero, the other should be too!
+ Note that conditional operators return 1 if true and 0 if false,
+ so if one is non-zero while the other is zero, then sum will be
+ 1. Otherwise the sum will either be 0 or 2.*/
+ switch ( (p->shift[0]!=0) + (p->shift[1]!=0) )
+ {
+ case 0:
+ /* `prepforconv' is only valid when xshift and yshift are both
+ zero. Also, a PSF profile should exist in the image. */
+ if(p->prepforconv)
+ {
+ /* Check if there is at least one Moffat or Gaussian profile. */
+ for(i=0;i<p->num;++i)
+ if( oneprofile_ispsf(p->f[i]) )
+ {
+ /* Calculate the size of the box holding the PSF. Note:
+
+ - For the Moffat and Gaussian profiles, the radius
+ columns is actually the FWHM which is actually the
+ diameter, not radius. So we have to divide it by
+ half.
+
+ - encloseellipse outputs the total width, we only want
+ half of it for the shift. */
+ 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);
+ p->shift[0] = (width[0]/2)*p->oversample;
+ p->shift[1] = (width[1]/2)*p->oversample;
+ }
+ }
+ break;
+
+ case 1:
+ error(EXIT_FAILURE, 0, "at least one of `--xshift` (`-X`) or "
+ "`--yshift` (`-Y`) are zero, they must either both be zero "
+ "or both given a positive value");
+ break;
+
+ case 2:
+ p->shift[0] *= p->oversample;
+ p->shift[1] *= p->oversample;
+ break;
+
+ default:
+ error(EXIT_FAILURE, 0, "a bug in ui_prepare_canvas! In checks "
+ "for shifts. Please contact us at %s so we can fix it",
+ PACKAGE_BUGREPORT);
+ }
+
+ /* Prepare the sizes of the final merged image (if it is to be
+ made). */
+ if(p->nomerged)
+ ndim=0;
+ else
+ {
+ /* Sanity check */
+ if(p->naxes[0]==0 || p->naxes[1]==0)
+ error(EXIT_FAILURE, 0, "No final merged image size is specified, "
+ "please use the `--naxis1', and `--naxis2' options to "
+ "specify the respective sizes");
+
+ /* Set the final merged image size. Note that even if we don't
+ want a merged image, we still need its WCS structure.*/
+ p->naxes[0] = (p->naxes[0] * p->oversample) + (2 * p->shift[0]);
+ p->naxes[1] = (p->naxes[1] * p->oversample) + (2 * p->shift[1]);
+ dsize[0] = p->naxes[1];
+ dsize[1] = p->naxes[0];
+ ndim = 2;
+ }
+
+ /* Make the output structure. */
+ p->out=gal_data_alloc(NULL, GAL_DATA_TYPE_FLOAT, ndim, 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);
+
+
+ /* 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)
+ {
+ if(p->out->name) free(p->out->name);
+ gal_checkset_allocate_copy("Mock profiles", &p->out->name);
+ if(p->out->unit==NULL)
+ gal_checkset_allocate_copy("Brightness", &p->out->unit);
+ }
+
+
+
+ /* When individual mode is requested, write the WCS structure to a header
+ 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)
+ {
+ status=wcshdo(WCSHDO_safe, p->out->wcs, &p->wcsnkeyrec, &p->wcsheader);
+ if(status)
+ error(EXIT_FAILURE, 0, "wcshdo error %d: %s", status,
+ wcs_errmsg[status]);
+ }
+}
+
+
+
+
+
+static void
+ui_finalize_coordinates(struct mkprofparams *p)
+{
+ size_t i;
+ double *x=NULL, *y=NULL;
+
+ /* 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
+ need to change them into actual image coordinates. */
+ if(p->racol)
+ {
+ /* 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);
+
+ /* If any conversions created a WCSLIB error, both the outputs will be
+ set to NaN. */
+ for(i=0;i<p->num;++i)
+ if( isnan(x[i]) )
+ error(EXIT_FAILURE, 0, "catalog row %zu: WCSLIB could not convert "
+ "(%f, %f) coordinates into image coordinates", i, p->x[i],
+ p->y[i]);
+
+ /* Free the RA and Dec arrays and put in the new image values. */
+ free(p->x);
+ free(p->y);
+ p->x=x;
+ p->y=y;
+ }
+
+
+ /* Correct the WCS scale. Note that when the WCS is read from a
+ 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. */
+ p->out->wcs->cdelt[0] /= p->oversample;
+ p->out->wcs->cdelt[1] /= p->oversample;
+
+
+
+ /* For a sanity check:
+ printf("\nui_finalize_coordinates sanity check:\n");
+ for(i=0;i<p->num;++i)
+ printf("%f, %f\n", p->x[i], p->y[i]);
+ */
+}
+
+
+
+
+/* Add all the columns of the log file. Just note that since this is a
+ linked list, we have to add them in the opposite order. */
+static void
+ui_make_log(struct mkprofparams *p)
{
+ char *comment;
+
+ /* Return if no long file is to be created. */
+ if(p->cp.log==0) return;
+
+ /* Individual created. */
+ gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_UCHAR, 1, &p->num, NULL,
+ 1, p->cp.minmapsize, "INDIV_CREATED", "bool",
+ "If an individual image was made (1) or not (0).");
+
+ /* Fraction of monte-carlo. */
+ gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_FLOAT, 1, &p->num, NULL,
+ 1, p->cp.minmapsize, "FRAC_MONTECARLO", "frac",
+ "Fraction of brightness in Monte-carlo integrated "
+ "pixels.");
+
+ /* Number of monte-carlo. */
+ gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_ULONG, 1, &p->num, NULL,
+ 1, p->cp.minmapsize, "NUM_MONTECARLO", "count",
+ "Number of Monte Carlo integrated pixels.");
+
+ /* Magnitude of profile overlap. */
+ gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_FLOAT, 1, &p->num, NULL,
+ 1, p->cp.minmapsize, "MAG_OVERLAP", "mag",
+ "Magnitude of profile's overlap with merged image.");
+
+ /* Row number in input catalog. */
+ if(gal_fits_name_is_fits(p->catname))
+ asprintf(&comment, "Row number of profile in %s (hdu: %s).",
+ p->catname, p->cp.hdu);
+ else
+ asprintf(&comment, "Row number of profile in %s.", p->catname);
+ gal_data_add_to_ll(&p->log, NULL, GAL_DATA_TYPE_ULONG, 1, &p->num, NULL,
+ 1, p->cp.minmapsize, "INPUT_ROW_NO", "count",
+ "Row number of profile in ");
+ free(comment);
+}
+
- /* Correct/set based on the given oversampling. */
- p->naxes[0] *= p->oversample;
- p->naxes[1] *= p->oversample;
- p->halfpixel = 0.5f/p->oversample;
+
+
+
+static void
+ui_preparations(struct mkprofparams *p)
+{
/* Read in all the columns. */
ui_read_cols(p);
+ /* Prepare the output canvas. */
+ ui_prepare_canvas(p);
+
+ /* Read the (possible) RA/Dec inputs into X and Y for the builder.*/
+ ui_finalize_coordinates(p);
+
+ /* Allocate the random number generator: */
+ gsl_rng_env_setup();
+ p->rng=gsl_rng_alloc(gsl_rng_default);
+
+ /* Make the log linked list. */
+ ui_make_log(p);
}
@@ -543,6 +935,52 @@ ui_preparations(struct mkprofparams *p)
/**************************************************************/
/************ Set the parameters *************/
/**************************************************************/
+static void
+ui_print_intro(struct mkprofparams *p)
+{
+ char *jobname;
+
+ if(p->cp.quiet) return;
+
+ printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
+
+ asprintf(&jobname, "%zu profile%sread from %s", p->num,
+ p->num>1?"s ":" ", p->catname);
+ gal_timing_report(NULL, jobname, 1);
+ free(jobname);
+
+ if(p->backname)
+ {
+ if(p->nomerged)
+ asprintf(&jobname, "WCS information read from %s", p->backname);
+ else
+ asprintf(&jobname, "%s is read and will be used as canvas",
+ p->backname);
+ gal_timing_report(NULL, jobname, 1);
+ free(jobname);
+ }
+
+ asprintf(&jobname, "Random number generator (RNG) type: %s",
+ gsl_rng_name(p->rng));
+ gal_timing_report(NULL, jobname, 1);
+ free(jobname);
+ if(p->envseed)
+ {
+ asprintf(&jobname, "RNG seed for all profiles: %zu",
+ gsl_rng_default_seed);
+ gal_timing_report(NULL, jobname, 1);
+ free(jobname);
+ }
+
+ asprintf(&jobname, "Using %zu threads.", p->cp.numthreads);
+ gal_timing_report(NULL, jobname, 1);
+ free(jobname);
+}
+
+
+
+
+
void
ui_read_check_inputs_setup(int argc, char *argv[], struct mkprofparams *p)
{
@@ -594,8 +1032,8 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
mkprofparams *p)
/* Read/allocate all the necessary starting arrays. */
ui_preparations(p);
- printf("\n... End of ui.c ...\n");
- exit(0);
+ /* Print introductory information. */
+ ui_print_intro(p);
}
@@ -621,16 +1059,13 @@ ui_read_check_inputs_setup(int argc, char *argv[],
struct mkprofparams *p)
/************ Free allocated, report *************/
/**************************************************************/
void
-freeandreport(struct mkprofparams *p, struct timeval *t1)
+ui_free_report(struct mkprofparams *p, struct timeval *t1)
{
- int status;
-
/* Free all the allocated arrays. */
free(p->cat);
free(p->cp.hdu);
free(p->outdir);
free(p->basename);
- if(p->individual==0) free(p->log);
/* p->cp.output might be equal to p->mergedimgname. In this case, if
we simply free them after each other, there will be a double free
@@ -648,11 +1083,6 @@ freeandreport(struct mkprofparams *p, struct timeval *t1)
if(p->individual)
free(p->wcsheader);
- /* Free the WCS structure. */
- if( (status=wcsvfree(&p->nwcs, &p->wcs)) )
- error(EXIT_FAILURE, 0, "wcsfree error %d: %s", status,
- wcs_errmsg[status]);
-
/* Free the random number generator: */
gsl_rng_free(p->rng);
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index fcabe30..c986541 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -27,6 +27,6 @@ void
ui_read_check_inputs_setup(int argc, char *argv[], struct mkprofparams *p);
void
-freeandreport(struct mkprofparams *p, struct timeval *t1);
+ui_free_report(struct mkprofparams *p, struct timeval *t1);
#endif
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6633eef..acaa76b 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5287,42 +5287,49 @@ each element interpretted as a @code{double} type (see
@ref{Data types}).
2 , 4.454 792 72.98348e7
@end example
-However, the example above has no other information about the columns. For
-example, Gnuastro's programs/libraries, you aren't limited to using the
-column's number. If the columns have names, units, or comments you can also
-select your columns based on searches/matches in these fields, for example
-see @ref{Table}. It is also bad for sending to a colleague, because they
-will find it hard to remember/use the columns properly. Also, in this
-manner, you can't guide the program reading the table on how to read the
-numbers. As an example, the first and third columns above can be read as
-integer types: the first column might be an ID and the third can be the
-number of pixels it occupies. So there is no need to read it as a
address@hidden type (which takes more memory, and is slower).
-
-In this bare-minimum example, you also can't use strings of characters, for
-example the names of filters, or some other identifier that includes
-non-numerical characters. In the absence of any information, only numbers
-can be read robustly. Assuming we read columns with non-numerical
+However, the example above has no other information about the columns (its
+raw data, with no meta-data). To use this table, you have to remember what
+each column was. Also, when you want to select columns, you have to count
+their position within the table. This can become frustrating and prone to
+bad errors (getting the columns wrong) especially as the number of columns
+increase. It is also bad for sending to a colleague, because they will find
+it hard to remember/use the columns properly.
+
+To solve these problems in Gnuastro's programs/libraries you aren't limited
+to using the column's number, see @ref{Selecting table columns}. If the
+columns have names, units, or comments you can also select your columns
+based on searches/matches in these fields, for example see @ref{Table}.
+Also, in this manner, you can't guide the program reading the table on how
+to read the numbers. As an example, the first and third columns above can
+be read as integer types: the first column might be an ID and the third can
+be the number of pixels an object occupies in an image. So there is no need
+to read these to columns as a @code{double} type (which takes more memory,
+and is slower).
+
+In the bare-minimum example above, you also can't use strings of
+characters, for example the names of filters, or some other identifier that
+includes non-numerical characters. In the absence of any information, only
+numbers can be read robustly. Assuming we read columns with non-numerical
characters as string, there would still be the problem that the strings
might contain space (or any delimiter) character for some rows. So, each
`word' in the string will be interpretted as a column and the program will
abort with an error that the rows don't have the same number of columns.
To correct for these limitations, Gnuastro defines the following convention
-for storing the table meta-data along with the plain text file to guide a
-human, or program, on how to read/interpret it. When the first non-white
-character in a line is @key{#}, or there are no non-white characters in it,
-then the line will be ignored (this is a pretty standard convention in many
+for storing the table meta-data along with the raw data in one plain text
+file. The format is primarily designed for ease of reading/writing by
+eye/fingers, but is also structured enough to be read by a program.
+
+When the first non-white character in a line is @key{#}, or there are no
+non-white characters in it, then the line will not be considered as a row
+of data in the table (this is a pretty standard convention in many
programs, and higher level languages). In the former case, the line is
interpretted as a @emph{comment}. If the comment line starts with
address@hidden
Column N:}', then it is assumed to contain information about column
@code{N} (a number, counting from 1). Comment lines that don't start with
this pattern are ignored and you can use them to include any further
-information you want to store with the table in the text file.
-
-The format is primarily defined for ease of reading/writing by eye/fingers,
-but is also structured enough to be read by a program. A column information
-comment is assumed to have the following format:
+information you want to store with the table in the text file. A column
+information comment is assumed to have the following format:
@example
# Column N: NAME [UNIT, TYPE, BLANK] COMMENT
@@ -5381,7 +5388,8 @@ with nothing but the column number is redundant):
# Column 3: mag_f160w [AB mag, f] Magnitude from the F160W filter
@end example
-The following type codes are recognized in Gnuastro's plain text tables:
+The data type of the column (see @ref{Data types}) should be specified with
+one of the following values:
@itemize
@item
@@ -5411,26 +5419,27 @@ The following type codes are recognized in Gnuastro's
plain text tables:
string (how many characters it has). The start of the string on each row is
the first non-delimiter character of the column that has the string
type. The next @code{N} characters will be interpretted as a string and all
-leading and trailing white space will be removed. If there is a new-line
-character (the string character is the last column, or we have reached the
-end of the line), then reading of the string will stop, even if the
address@hidden characters are not complete yet. See @file{tests/table/table.txt}
-for one example.
+leading and trailing white space will be removed.
-So the only time you have to pay attention to the positioning of the string
-column is not the last column. If the next column's characters, are closer
-than @code{N} characters to the start of the string column in that
-line/row, they will be considered part of the string column.
+If the next column's characters, are closer than @code{N} characters to the
+start of the string column in that line/row, they will be considered part
+of the string column. If there is a new-line character before the ending of
+the space given to the string column (in other words, the string column is
+the last column), then reading of the string will stop, even if the
address@hidden characters are not complete yet. See @file{tests/table/table.txt}
+for one example. Therefore, the only time you have to pay attention to the
+positioning and spaces given to the string column is when it is not the
+last column in the table.
The only limitation in this format is that trailing and leading white space
-characters will be removed from the columns that are read. So if trailing
-and leading white-spaces are critically important to your analysis, define
-your own starting and ending characters and remove them after the table has
-been read. For example in the sample table below, the two address@hidden|}'
-characters (which are arbitrary) will remain in the value of the second
-column and you can remove them manually later. If only one of the leading
-or trailing white spaces is important for your work, you can only use one
-of the address@hidden|}'s.
+characters will be removed from the columns that are read. In most cases,
+this is the desired behavior, but if trailing and leading white-spaces are
+critically important to your analysis, define your own starting and ending
+characters and remove them after the table has been read. For example in
+the sample table below, the two address@hidden|}' characters (which are
arbitrary)
+will remain in the value of the second column and you can remove them
+manually later. If only one of the leading or trailing white spaces is
+important for your work, you can only use one of the address@hidden|}'s.
@example
# Column 1: ID [label, uc]
@@ -13485,16 +13494,16 @@ Carlo integration). This is also, generally speaking,
what happens in
practice with the photons on the pixel. The number of random points
can be set with @option{--numrandom}.
-Unfortunately, repeating this Monte Carlo process would be extremely
-time and CPU consuming if it is to be applied to every pixel. In order
-to not loose too much accuracy, in MakeProfiles, the profile is built
-using both methods explained above. The building of the profile begins
-from its central pixel and continues outwards. Monte Carlo integration
-is first applied (which yields @mymath{F_r}), then the central pixel
-value (@mymath{F_c}) on the same pixel. If the fractional difference
-(@mymath{|F_r-F_c|/F_r}) is lower than a given tolerance level we will
-stop using Monte Carlo integration and only use the central pixel
-value.
+Unfortunately, repeating this Monte Carlo process would be extremely time
+and CPU consuming if it is to be applied to every pixel. In order to not
+loose too much accuracy, in MakeProfiles, the profile is built using both
+methods explained below. The building of the profile begins from its
+central pixel and continues (radially) outwards. Monte Carlo integration is
+first applied (which yields @mymath{F_r}), then the central pixel value
+(@mymath{F_c}) is calculated on the same pixel. If the fractional
+difference (@mymath{|F_r-F_c|/F_r}) is lower than a given tolerance level
+(specified with @option{--tolerance}) MakeProfiles will stop using Monte
+Carlo integration and only use the central pixel value.
@cindex Inside-out construction
The ordering of the pixels in this inside-out construction is based on
@@ -13755,6 +13764,10 @@ the command-line or in the configuration files. Note
that
@option{--oversample} will remain active even if a background image is
specified.
+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.
@@ -14056,11 +14069,17 @@ column number, where counting starts from zero.
@table @option
@item --fcol=INT/STR
-The functional form of the profile with one of the values below. The column
-can contain both string characters (for example address@hidden') and the
-numeric codes (address@hidden'). The string format can be easier when the
-catalog is to be read, or written by hand. For ways to benefit from string
-columns in a plain text table/catalog, see @ref{Gnuastro text table
+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 address@hidden') or string characters (for example
address@hidden'). 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
diff --git a/lib/data.c b/lib/data.c
index f853412..e7adcb0 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -561,7 +561,7 @@ gal_data_mmap(gal_data_t *data, int clear)
/* Write to the newly set file position so the space is allocated. */
- if( write(filedes, &uc, 1) == -1)
+ if( write(filedes, &uc, bsize) == -1)
error(EXIT_FAILURE, errno, "%s: unable to write one byte at the "
"%zu-th position", filename, bsize);
@@ -1849,11 +1849,64 @@ gal_data_flag_blank(gal_data_t *data)
************** Types and copying ***************
*************************************************************/
-/* gal_data_copy_to_new_type: Macro for all tyeps. */
-#define COPY_OTYPE_ITYPE_SET(otype, itype) { \
- itype *ia=in->array; \
- otype *oa=out->array, *of=oa+out->size; \
- do *oa=*ia++; while(++oa<of); \
+/* Copy to a new type for integers. */
+#define COPY_OTYPE_ITYPE_SET_INT(otype, itype) { \
+ itype *ia=in->array, iblank; \
+ otype *oa=out->array, *of=oa+out->size, oblank; \
+ \
+ /* Check if there are blank values in the input array and that */ \
+ /* the types of the two structures are different. */ \
+ if( in->type!=newtype && gal_data_has_blank(in) ) \
+ { \
+ /* Set the blank values */ \
+ gal_data_set_blank(&iblank, in->type); \
+ gal_data_set_blank(&oblank, newtype); \
+ \
+ /* Copy the input to the output. */ \
+ do { *oa = *ia==iblank ? oblank : *ia; ia++; } while(++oa<of); \
+ } \
+ \
+ /* There were no blank elements in the input, or the input and */ \
+ /* output have the same type. */ \
+ else \
+ do *oa=*ia++; while(++oa<of); \
+ }
+
+
+
+
+
+/* Copy to a new type for floating point values. */
+#define COPY_OTYPE_ITYPE_SET_FLT(otype, itype) { \
+ itype *ia=in->array, iblank; \
+ otype *oa=out->array, *of=oa+out->size, oblank; \
+ \
+ /* Check if there are blank values in the input array and that */ \
+ /* the types of the two structures are different. */ \
+ if( in->type!=newtype && gal_data_has_blank(in) ) \
+ { \
+ /* Set the blank values */ \
+ gal_data_set_blank(&iblank, in->type); \
+ gal_data_set_blank(&oblank, newtype); \
+ \
+ /* When the blank value isn't NaN, then we should use the */ \
+ /* equal operator to check for blank values. */ \
+ if( isnan(iblank) ) \
+ { \
+ do { *oa = isnan(*ia) ? oblank : *ia; ia++; } \
+ while(++oa<of); \
+ } \
+ else \
+ { \
+ do { *oa = *ia==iblank ? oblank : *ia; ia++; } \
+ while(++oa<of); \
+ } \
+ } \
+ \
+ /* There were no blank elements in the input, or the input and */ \
+ /* output have the same type. */ \
+ else \
+ do *oa=*ia++; while(++oa<of); \
}
@@ -1866,47 +1919,56 @@ gal_data_flag_blank(gal_data_t *data)
switch(in->type) \
{ \
case GAL_DATA_TYPE_UCHAR: \
- COPY_OTYPE_ITYPE_SET(otype, unsigned char); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, unsigned char); \
break; \
\
case GAL_DATA_TYPE_CHAR: \
- COPY_OTYPE_ITYPE_SET(otype, char); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, char); \
break; \
\
case GAL_DATA_TYPE_USHORT: \
- COPY_OTYPE_ITYPE_SET(otype, unsigned short); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, unsigned short); \
break; \
\
case GAL_DATA_TYPE_SHORT: \
- COPY_OTYPE_ITYPE_SET(otype, short); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, short); \
break; \
\
case GAL_DATA_TYPE_UINT: \
- COPY_OTYPE_ITYPE_SET(otype, unsigned int); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, unsigned int); \
break; \
\
case GAL_DATA_TYPE_INT: \
- COPY_OTYPE_ITYPE_SET(otype, int); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, int); \
break; \
\
case GAL_DATA_TYPE_ULONG: \
- COPY_OTYPE_ITYPE_SET(otype, unsigned long); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, unsigned long); \
break; \
\
case GAL_DATA_TYPE_LONG: \
- COPY_OTYPE_ITYPE_SET(otype, long); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, long); \
break; \
\
case GAL_DATA_TYPE_LONGLONG: \
- COPY_OTYPE_ITYPE_SET(otype, LONGLONG); \
+ COPY_OTYPE_ITYPE_SET_INT(otype, LONGLONG); \
break; \
\
case GAL_DATA_TYPE_FLOAT: \
- COPY_OTYPE_ITYPE_SET(otype, float); \
+ COPY_OTYPE_ITYPE_SET_FLT(otype, float); \
break; \
\
case GAL_DATA_TYPE_DOUBLE: \
- COPY_OTYPE_ITYPE_SET(otype, double); \
+ COPY_OTYPE_ITYPE_SET_FLT(otype, double); \
+ break; \
+ \
+ case GAL_DATA_TYPE_BIT: \
+ case GAL_DATA_TYPE_STRLL: \
+ case GAL_DATA_TYPE_COMPLEX: \
+ case GAL_DATA_TYPE_DCOMPLEX: \
+ error(EXIT_FAILURE, 0, "`gal_data_copy_to_new_type' currently " \
+ "doesn't support copying from %s type", \
+ gal_data_type_as_string(in->type, 1)); \
break; \
\
default: \
@@ -1981,6 +2043,20 @@ gal_data_copy_to_new_type(gal_data_t *in, int newtype)
COPY_OTYPE_SET(double);
break;
+ case GAL_DATA_TYPE_BIT:
+ case GAL_DATA_TYPE_STRLL:
+ case GAL_DATA_TYPE_COMPLEX:
+ case GAL_DATA_TYPE_DCOMPLEX:
+ error(EXIT_FAILURE, 0, "`gal_data_copy_to_new_type' currently doesn't "
+ "support copying to %s type",
+ gal_data_type_as_string(newtype, 1));
+ break;
+
+ case GAL_DATA_TYPE_STRING:
+ error(EXIT_FAILURE, 0, "`gal_data_copy_to_new_type' is only for "
+ "numeric data types, the desired new type is a string");
+ break;
+
default:
error(EXIT_FAILURE, 0, "type %d not recognized for "
"for newtype in gal_data_copy_to_new_type", newtype);
diff --git a/lib/gnuastro/threads.h b/lib/gnuastro/threads.h
index 48d3b1e..6c0b912 100644
--- a/lib/gnuastro/threads.h
+++ b/lib/gnuastro/threads.h
@@ -96,9 +96,14 @@ pthread_barrier_destroy(pthread_barrier_t *b);
-/*****************************************************************/
-/**************** gnuastro functions ******************/
-/*****************************************************************/
+
+
+/*******************************************************************/
+/************ Thread utilities **************/
+/*******************************************************************/
+size_t
+gal_threads_number();
+
void
gal_threads_dist_in_threads(size_t numactions, size_t numthreads,
size_t **outthrds, size_t *outthrdcols);
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index aa81d29..b2f4b11 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -53,8 +53,8 @@ gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy,
double *radec,
size_t number, size_t width);
void
-gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec, double *xy,
- size_t number, size_t width);
+gal_wcs_world_to_img(struct wcsprm *wcs, double *ra, double *dec,
+ double **x, double **y, size_t size);
double
gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
diff --git a/lib/options.c b/lib/options.c
index 954364c..4326c5a 100644
--- a/lib/options.c
+++ b/lib/options.c
@@ -34,7 +34,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/arithmetic.h>
#include <gnuastro/linkedlist.h>
-#include <nproc.h> /* from Gnulib, in Gnuastro's source */
+
#include <options.h>
#include <checkset.h>
@@ -49,16 +49,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/**********************************************************************/
/************ Option utilities ***************/
/**********************************************************************/
-void
-gal_options_initialize_numthreads(struct gal_options_common_params *cp)
-{
- cp->numthreads=num_processors(NPROC_CURRENT);
-}
-
-
-
-
-
int
gal_options_is_last(struct argp_option *option)
{
diff --git a/lib/options.h b/lib/options.h
index b8a6483..3deb04c 100644
--- a/lib/options.h
+++ b/lib/options.h
@@ -193,9 +193,6 @@ struct gal_options_common_params
/**********************************************************************/
/************ Option utilities ***************/
/**********************************************************************/
-void
-gal_options_initialize_numthreads(struct gal_options_common_params *cp);
-
int
gal_options_is_last(struct argp_option *option);
diff --git a/lib/table.c b/lib/table.c
index f012414..15fe7a8 100644
--- a/lib/table.c
+++ b/lib/table.c
@@ -808,6 +808,10 @@ gal_table_read(char *filename, char *hdu, struct
gal_linkedlist_stll *cols,
/************************************************************************/
/*************** Write a table ***************/
/************************************************************************/
+/* The input is a linked list of data structures and some comments. The
+ table will then be written into `filename' with a format that is
+ specified by `tableformat'. If it already exists, and `dontdelete' has a
+ value of 1, then it won't be deleted and an error will be printed. */
void
gal_table_write(gal_data_t *cols, char *comments, int tableformat,
char *filename, int dontdelete)
diff --git a/lib/threads.c b/lib/threads.c
index 96a55d0..099312e 100644
--- a/lib/threads.c
+++ b/lib/threads.c
@@ -30,7 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/threads.h>
-
+#include <nproc.h> /* from Gnulib, in Gnuastro's source */
@@ -166,8 +166,18 @@ pthread_barrier_destroy(pthread_barrier_t *b)
/*******************************************************************/
-/************ Distribute job indexs in threads **************/
+/************ Thread utilities **************/
/*******************************************************************/
+size_t
+gal_threads_number()
+{
+ return num_processors(NPROC_CURRENT);
+}
+
+
+
+
+
/* We have `numactions` jobs and we want their indexs to be divided
between `numthreads` CPU threads. This function will give each index to
a thread such that the maximum difference between the number of
diff --git a/lib/wcs.c b/lib/wcs.c
index 0c716df..e4cdd02 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -83,39 +83,67 @@ gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy,
double *radec,
-/* Similar to the gal_wcs_xyarray_to_radec, but the reverse: to convert an
- array of RA-Dec to X-Y. */
+/* Convert an array of world coordinates to image coordinates. Note that in
+ Gnuastro, each column is treated independently, so the inputs are
+ separate. If `*x==NULL', or `*y==NULL', then space will be allocated for
+ them, otherwise, it is assumed that space has already been
+ allocated. Note that they must be a 1 dimensional array (recall that in
+ Gnuastro columns are treated independently).*/
void
-gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec, double *xy,
- size_t number, size_t stride)
+gal_wcs_world_to_img(struct wcsprm *wcs, double *ra, double *dec,
+ double **x, double **y, size_t size)
{
size_t i;
- double imgcrd[2], phi, theta;
- int stat, status=0, ncoord=1, nelem=2;
-
- for(i=0;i<number;++i)
+ int status, *stat, ncoord=size, nelem=2;
+ double *phi, *theta, *world, *pixcrd, *imgcrd;
+
+ /* Allocate all the necessary arrays. */
+ stat=gal_data_calloc_array(GAL_DATA_TYPE_INT, size);
+ phi=gal_data_malloc_array(GAL_DATA_TYPE_DOUBLE, size);
+ theta=gal_data_malloc_array(GAL_DATA_TYPE_DOUBLE, size);
+ world=gal_data_malloc_array(GAL_DATA_TYPE_DOUBLE, 2*size);
+ imgcrd=gal_data_malloc_array(GAL_DATA_TYPE_DOUBLE, 2*size);
+ pixcrd=gal_data_malloc_array(GAL_DATA_TYPE_DOUBLE, 2*size);
+
+ /* Put in the values. */
+ for(i=0;i<size;++i) { world[i*2]=ra[i]; world[i*2+1]=dec[i]; }
+
+ /* Use WCSLIB's `wcss2p'. */
+ status=wcss2p(wcs, ncoord, nelem, world, phi, theta, imgcrd, pixcrd, stat);
+ if(status)
+ error(EXIT_FAILURE, 0, "wcss2p ERROR %d: %s", status,
+ wcs_errmsg[status]);
+
+ /* For a sanity check:
+ printf("\n\ngal_wcs_world_to_img sanity check:\n");
+ for(i=0;i<size;++i)
+ printf("world (%f, %f) --> pix (%f, %f), [stat: %d]\n",
+ world[i*2], world[i*2+1], pixcrd[i*2], pixcrd[i*2+1], stat[i]);
+ */
+
+ /* Allocate the output arrays if they were not already allocated. */
+ if(*x==NULL) *x=gal_data_calloc_array(GAL_DATA_TYPE_DOUBLE, size);
+ if(*y==NULL) *y=gal_data_calloc_array(GAL_DATA_TYPE_DOUBLE, size);
+
+ /* Put the values into the output arrays. */
+ for(i=0;i<size;++i)
{
- if(isnan(radec[i*stride]) || isnan(radec[i*stride+1]))
- radec[i*stride]=radec[i*stride+1]=NAN;
- else
- {
- status=wcss2p(wcs, ncoord, nelem, radec+i*stride, &phi, &theta,
- imgcrd, xy+i*stride, &stat);
- if(status)
- error(EXIT_FAILURE, 0, "wcss2p ERROR %d: %s", status,
- wcs_errmsg[status]);
-
- /* For a check:
- printf("(%f, %f) --> (%f, %f)\n", radec[i*stride], radec[i*stride+1],
- xy[i*stride], xy[i*stride+1]);
- */
- }
+ (*x)[i] = stat[i] ? NAN : pixcrd[i*2];
+ (*y)[i] = stat[i] ? NAN : pixcrd[i*2+1];
}
+
+ /* Clean up. */
+ free(phi);
+ free(stat);
+ free(theta);
+ free(world);
+ free(pixcrd);
}
+
/* The distance (along a great circle) on a sphere between two points
is calculated here. Since the pixel sides are usually very small,
we won't be using the direct formula:
diff --git a/tests/mkprof/mkprofcat1.txt b/tests/mkprof/mkprofcat1.txt
index c326154..92cdfa6 100644
--- a/tests/mkprof/mkprofcat1.txt
+++ b/tests/mkprof/mkprofcat1.txt
@@ -1,14 +1,13 @@
-# Random catalog generated by makerandomcat.py
-# Column 0: ID
-# Column 1: X
-# Column 2: Y
-# Column 3: Function (0: Sersic, 1: Moffat, 2: Gaussian, 3: Point).
-# Column 4: Effective radius
-# Column 5: Sersic index.
-# Column 6: Position angle
-# Column 7: Axis ratio
-# Column 8: Magnitude
-# Column 9: Truncation radius
-1 0.0000 0.0000 1 3.000 4.765 0.0000 1.000
0.000 5.000
-2 100.40 100.40 0 20.00 2.500 45.000 1.000
-18.0 5.000
-3 50.321 50.827 0 5.978 1.320 77.650 0.801
-15.0 5.000
+# Column 1: ID [count, uc] Object ID
+# Column 2: X [pixel, d] X axis position of profile center
+# Column 3: Y [pixel, d] Y axis position of profile center
+# Column 4: Function [name, str7] Profile's radial function
+# Column 5: Width [pixel, d] For Sersic: effective radius, for
Moffat, FWHM
+# Column 6: Sersic index [none, d] Sersic index, or Moffat beta
+# Column 7: Position angle [deg, d] Position angle of profile
+# Column 8: Axis ratio [frac, d] Axis ratio of profile
+# Column 9: Magnitude [ABmag, d] Magnitude of profile within
truncation radius
+# Column 10: Truncation radius [dist, d] Truncation radius to stop building
profile
+1 0.0000 0.0000 moffat 3.000 4.765 0.0000 1.000
0.000 5.000
+2 100.40 100.40 sersic 20.00 2.500 45.000 1.000
-18.0 5.000
+3 50.321 50.827 sersic 5.978 1.320 77.650 0.801
-15.0 5.000
diff --git a/tests/mkprof/mkprofcat2.txt b/tests/mkprof/mkprofcat2.txt
index d6d49db..aac3cff 100644
--- a/tests/mkprof/mkprofcat2.txt
+++ b/tests/mkprof/mkprofcat2.txt
@@ -1 +1,11 @@
-1 0.400 100.40 0 20.00 2.500 45.000 1.00
-18.0 5.000
+# Column 1: ID [count, uc] Object ID
+# Column 2: X [pixel, d] X axis position of profile center
+# Column 3: Y [pixel, d] Y axis position of profile center
+# Column 4: Function [code, i] Profile's radial function code
+# Column 5: Width [pixel, d] For Sersic: effective radius, for
Moffat, FWHM
+# Column 6: Sersic index [none, d] Sersic index, or Moffat beta
+# Column 7: Position angle [deg, d] Position angle of profile
+# Column 8: Axis ratio [frac, d] Axis ratio of profile
+# Column 9: Magnitude [ABmag, d] Magnitude of profile within
truncation radius
+# Column 10: Truncation radius [dist, d] Truncation radius to stop building
profile
+1 0.400 100.40 1 20.00 2.500 45.000 1.00
-18.0 5.000
diff --git a/tests/mkprof/mkprofcat3.txt b/tests/mkprof/mkprofcat3.txt
index e3fccdb..a6b05d6 100644
--- a/tests/mkprof/mkprofcat3.txt
+++ b/tests/mkprof/mkprofcat3.txt
@@ -1 +1,11 @@
-1 100.40 0.40 0 20.00 2.500 45.000 1.00
-18.0 5.000
+# Column 1: ID [count, uc] Object ID
+# Column 2: X [pixel, d] X axis position of profile center
+# Column 3: Y [pixel, d] Y axis position of profile center
+# Column 4: Function [code, i] Profile's radial function code
+# Column 5: Width [pixel, d] For Sersic: effective radius, for
Moffat, FWHM
+# Column 6: Sersic index [none, d] Sersic index, or Moffat beta
+# Column 7: Position angle [deg, d] Position angle of profile
+# Column 8: Axis ratio [frac, d] Axis ratio of profile
+# Column 9: Magnitude [ABmag, d] Magnitude of profile within
truncation radius
+# Column 10: Truncation radius [dist, d] Truncation radius to stop building
profile
+1 100.40 0.40 1 20.00 2.500 45.000 1.00
-18.0 5.000
diff --git a/tests/mkprof/mkprofcat4.txt b/tests/mkprof/mkprofcat4.txt
index 518f8bf..2f2b221 100644
--- a/tests/mkprof/mkprofcat4.txt
+++ b/tests/mkprof/mkprofcat4.txt
@@ -1 +1,11 @@
-1 0.40 0.40 0 20.00 2.500 45.000 1.00 -18.0
5.000
+# Column 1: ID [count, uc] Object ID
+# Column 2: X [pixel, d] X axis position of profile center
+# Column 3: Y [pixel, d] Y axis position of profile center
+# Column 4: Function [name, str7] Profile's radial function
+# Column 5: Width [pixel, d] For Sersic: effective radius, for
Moffat, FWHM
+# Column 6: Sersic index [none, d] Sersic index, or Moffat beta
+# Column 7: Position angle [deg, d] Position angle of profile
+# Column 8: Axis ratio [frac, d] Axis ratio of profile
+# Column 9: Magnitude [ABmag, d] Magnitude of profile within
truncation radius
+# Column 10: Truncation radius [dist, d] Truncation radius to stop building
profile
+1 0.40 0.40 sersic 20.00 2.500 45.000 1.00
-18.0 5.000
diff --git a/tests/mkprof/radeccat.sh b/tests/mkprof/radeccat.sh
index 023ee8e..a3a73aa 100755
--- a/tests/mkprof/radeccat.sh
+++ b/tests/mkprof/radeccat.sh
@@ -46,4 +46,4 @@ if [ ! -f $execname ]; then exit 77; fi
# Actual test script
# ==================
-$execname $cat --racol=1 --deccol=2 --naxis1=100 --naxis2=100
+$execname $cat --racol=RA --deccol=Dec --naxis1=100 --naxis2=100
diff --git a/tests/mkprof/radeccat.txt b/tests/mkprof/radeccat.txt
index 90505b8..7fc7b42 100644
--- a/tests/mkprof/radeccat.txt
+++ b/tests/mkprof/radeccat.txt
@@ -1,13 +1,12 @@
-# Random catalog generated by makerandomcat.py
-# Column 0: ID
-# Column 1: X
-# Column 2: Y
-# Column 3: Function (0: Sersic, 1: Moffat, 2: Gaussian, 3: Point).
-# Column 4: Effective radius
-# Column 5: Sersic index.
-# Column 6: Position angle
-# Column 7: Axis ratio
-# Column 8: Magnitude
-# Column 9: Truncation radius
-1 0.99987331 1.0001617 0 10.00 2.500 45.000 0.3
-6.0 2.000
-2 0.99929656 1.0006917 0 10.00 2.500 135.00 0.3
-6.0 2.000
+# Column 1: ID [count, uc] Object ID
+# Column 2: RA [deg, d] Right Ascension of profile center
+# Column 3: Dec [deg, d] Declination of profile center
+# Column 4: Function [name, str7] Profile's radial function
+# Column 5: Width [pixel, d] For Sersic: effective radius, for
Moffat, FWHM
+# Column 6: Sersic index [none, d] Sersic index, or Moffat beta
+# Column 7: Position angle [deg, d] Position angle of profile
+# Column 8: Axis ratio [frac, d] Axis ratio of profile
+# Column 9: Magnitude [ABmag, d] Magnitude of profile within
truncation radius
+# Column 10: Truncation radius [dist, d] Truncation radius to stop building
profile
+1 0.99987331 1.0001617 sersic 10.00 2.500 45.000 0.3
-6.0 2.000
+2 0.99929656 1.0006917 sersic 10.00 2.500 135.00 0.3
-6.0 2.000
diff --git a/tmpfs-config-make b/tmpfs-config-make
index e337626..a6a2f59 100755
--- a/tmpfs-config-make
+++ b/tmpfs-config-make
@@ -130,8 +130,8 @@ cd $build_dir
#
# ####################################
if [ ! -f Makefile ]; then
- $srcdir/configure --srcdir=$srcdir --disable-shared CFLAGS="-g -O0"\
- --enable-arithmetic --enable-table --enable-mkprof
+ $srcdir/configure --srcdir=$srcdir --disable-shared CFLAGS="-g -O0" \
+ --enable-arithmetic --enable-mkprof --enable-table
fi
- [gnuastro-commits] master 448d31e 092/125: Minor typo corrections in options of Header and Statistics, (continued)
- [gnuastro-commits] master 448d31e 092/125: Minor typo corrections in options of Header and Statistics, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c315491 075/125: Header using new option management system, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 581de34 084/125: Option strings directly converted to internal codes, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e2c7fff 054/125: Option values now usable in program, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 8f390d4 079/125: Convolve works with gal_data_t features, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 2e31ae4 103/125: Added -j8 to the suggestions after ./configure and make, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 18c9823 094/125: Image created on second HDU, common options to gnuastro.conf, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db9f983 065/125: Arithmetic uses new option management, no more mask images, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c1dba50 070/125: Column I/O options are now common options, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 5eb817a 052/125: Option values stored in argp_option, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 7e888ce 073/125: First implementation of MakeProfiles with gal_data_t,
Mohammad Akhlaghi <=
- [gnuastro-commits] master e8d9751 087/125: MakeNoise now uses the new gal_data_t infrastructure, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e978d44 068/125: Numeric option sanity checks, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e74eb25 093/125: Changed name of Header program to Fits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e353932 109/125: NaN magnitude will only blank the profile's area, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 96ab804 098/125: Working on tiles implemented in tile checking, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 31b3b76 099/125: Spatial convolution now using new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c0dd6db 102/125: The noedgecorrection option in spatial domain convolution, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dfdc34 069/125: Option values directly set in library, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 95ba655 100/125: Spatial convolution with new tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0dec9ac 110/125: Allow for floating point errors in Crop comparisons, Mohammad Akhlaghi, 2017/04/23