[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master f5a1bfc 1/2: MakeProfiles also accepts WCS pos
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master f5a1bfc 1/2: MakeProfiles also accepts WCS positions |
|
Date: |
Wed, 10 Aug 2016 21:18:04 +0000 (UTC) |
branch: master
commit f5a1bfc3f665e92916f85fe0dfca098429444592
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
MakeProfiles also accepts WCS positions
Until now the profile centers in MakeProfiles could only be specified by
image coordinates (X and Y). However, for most practical purposes, the
World Coordinate System (WCS, using right ascension and declination) is
needed. So MakeProfiles was modified to also be able to read WCS positions
in the input catalog when the user asks so. To do this, two options were
added for the user to specify the RA and Dec columns. When these columns
are given, the X and Y coordinates are no longer used.
This process mainly involved the starting of MakeProfiles (mostly in
`ui.c'). The construction of the WCS structure (when a background image is
not provided) was moved to `ui.c' (from `mkprof.c'). Having the WCS
structure, the given RA and Decs are converted to X and Y, replaced with
the RA and Dec and the `xcol' and `ycol' values set to the old RA and Dec
columns. Reading the WCS structure in `ui.c' also simplified the output
writing parts of `mkprof.c'.
While doing this I noticed that if the given oversampling factor is not
odd, then it is incremented to become odd! This is very bad behavior (I
honestly don't know what I was thinking when I put this line!). If the user
is expecting to oversample by 2, we can't secretly oversample by 3 for
example, we will just greatly confuse the user! Instead, now an error
message is printed which explains that oversampling must be an odd number.
While doing the work, in some places, the indentation were also corrected
(mainly following the refactoring of function names in commit c63305c).
A simple test was also written to build profiles based on their input WCS
coordinates.
This finishes task #13566.
---
doc/gnuastro.texi | 20 +++-
lib/fits.c | 9 +-
lib/wcs.c | 27 +++---
src/mkcatalog/columns.c | 3 +-
src/mkprof/Makefile.am | 8 +-
src/mkprof/args.h | 28 +++++-
src/mkprof/main.h | 4 +
src/mkprof/mkprof.c | 116 ++++++-----------------
src/mkprof/ui.c | 222 ++++++++++++++++++++++++++++++++++++++++++---
tests/Makefile.am | 27 +++---
tests/mkprof/radeccat.sh | 49 ++++++++++
tests/mkprof/radeccat.txt | 13 +++
12 files changed, 385 insertions(+), 141 deletions(-)
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index ba9727d..88b7e67 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12652,11 +12652,27 @@ used for making an elliptical annulus (with a width
of 1 or 2 pixels).
@item --xcol
(@option{=INT}) The center of the profiles along the first FITS axis
-(horizontal when viewed in SAO ds9).
+(horizontal when viewed in SAO ds9). See the explanations for
address@hidden for precedence when both image and WCS coordinate columns
+are given.
@item --ycol
(@option{=INT}) The center of the profiles along the second FITS axis
-(vertical when viewed in SAO ds9).
+(vertical when viewed in SAO ds9). Similar to @option{--xcol}.
+
address@hidden --racol
+(@option{=INT}) The profile center's right ascension. Along with
address@hidden, these WCS coordinate columns are not mandatory. If they
+are not given, the @option{--xcol} and @option{--ycol} options will be used
+to specify the profile's central position and vice-versa. However, if image
+coordinate columns (@option{--xcol} and @option{--ycol}) and WCS coordinate
+columns (@option{--racol} and @option{--deccol}) are given, the WCS
+coordinate columns take precedence and image coordinate columns will be
+ignored.
+
address@hidden --deccol
+(@option{=INT}) The profile center's declination. Similar to
address@hidden
@item --rcol
(@option{=INT}) The radius parameter of the profiles. Effective radius
diff --git a/lib/fits.c b/lib/fits.c
index 5030001..afe8ffd 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -1417,11 +1417,14 @@ gal_fits_array_to_file(char *filename, char *hdu, int
bitpix,
+/* `atof' stands for "array to file". This is essentially the same as
+ `gal_fits_array_to_file' except that the WCS structure's CRPIX values
+ have changed. */
void
gal_fits_atof_correct_wcs(char *filename, char *hdu, int bitpix,
- void *array, size_t s0, size_t s1,
- char *wcsheader, int wcsnkeyrec,
- double *crpix, char *spack_string)
+ void *array, size_t s0, size_t s1,
+ char *wcsheader, int wcsnkeyrec,
+ double *crpix, char *spack_string)
{
fitsfile *fptr;
int status=0, datatype;
diff --git a/lib/wcs.c b/lib/wcs.c
index bb2fc4e..fd20e3b 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -41,17 +41,15 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/**************************************************************/
/********** XY to RADEC ************/
/**************************************************************/
-/* Convert the X and y coordinates in a larger array to ra and dec
- coordinates in that same array. xy points to the first element in
- the X column and radec points to the first element in the RA
- column. The columns for Y and Dec have to be immediately after X
- and RA.
-
- It appears that WCSLIB can only deal with static allocation. At
- least in its tests it only uses static allocation. I tried dynamic
- allocation, but it didn't work. So I can't use the vector
- functionalities of WCSLIB and have to translate each point
- separately.
+/* Use the X and Y columns in a larger array to fill the RA and Dec columns
+ in that same array. `xy' points to the first element in the X column and
+ `radec' points to the first element in the RA column. The columns for Y
+ and Dec have to be immediately after X and RA.
+
+ It appears that WCSLIB can only deal with static allocation. At least in
+ its tests it only uses static allocation. I tried dynamic allocation,
+ but it didn't work. So I can't use the vector functionalities of WCSLIB
+ and have to translate each point separately.
*/
void
gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy, double *radec,
@@ -85,6 +83,8 @@ 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. */
void
gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec, double *xy,
size_t number, size_t width)
@@ -104,9 +104,10 @@ gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double
*radec, double *xy,
if(status)
error(EXIT_FAILURE, 0, "wcss2p ERROR %d: %s", status,
wcs_errmsg[status]);
+
/* For a check:
- printf("(%f, %f) --> (%f, %f)\n", xy[i*width], xy[i*width+1],
- radec[i*width], radec[i*width+1]);
+ printf("(%f, %f) --> (%f, %f)\n", radec[i*width], radec[i*width+1],
+ xy[i*width], xy[i*width+1]);
*/
}
}
diff --git a/src/mkcatalog/columns.c b/src/mkcatalog/columns.c
index d3cdd99..5ef2532 100644
--- a/src/mkcatalog/columns.c
+++ b/src/mkcatalog/columns.c
@@ -541,8 +541,7 @@ preparewcs(struct mkcatalogparams *p, size_t col)
the first row is not used by any object or colump (since
their indexes begin from 1).*/
gal_wcs_xy_array_to_radec(p->wcs, p->info+p->icols+xc,
- p->info+p->icols+rc,
- p->num, p->icols);
+ p->info+p->icols+rc, p->num, p->icols);
/* Set the flag of the converted columns to 1.0f, so the
diff --git a/src/mkprof/Makefile.am b/src/mkprof/Makefile.am
index 7ba9815..a40d592 100644
--- a/src/mkprof/Makefile.am
+++ b/src/mkprof/Makefile.am
@@ -28,10 +28,10 @@ bin_PROGRAMS = astmkprof
astmkprof_SOURCES = main.c main.h ui.c ui.h mkprof.c mkprof.h args.h \
oneprofile.c oneprofile.h profiles.c profiles.h cite.h
-astmkprof_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la \
--lgalconfigfiles -lgalfits -lgaltxtarray -lgalcheckset -lgallinkedlist \
--lgaltiming -lgalthreads -lgalbox -lgalstatistics -lgalarraymanip \
--lgalqsort
+astmkprof_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la \
+-lgalconfigfiles -lgalfits -lgalwcs -lgaltxtarray -lgalcheckset \
+-lgallinkedlist -lgaltiming -lgalthreads -lgalbox -lgalstatistics \
+-lgalarraymanip -lgalqsort
diff --git a/src/mkprof/args.h b/src/mkprof/args.h
index e316892..05f18dc 100644
--- a/src/mkprof/args.h
+++ b/src/mkprof/args.h
@@ -63,7 +63,7 @@ const char doc[] =
a d f g j k l u v
B C E F G H I J L O Q T U W Z
- Maximum integer used so far: 515.
+ Maximum integer used so far: 517
*/
static struct argp_option options[] =
{
@@ -279,6 +279,22 @@ static struct argp_option options[] =
4
},
{
+ "racol",
+ 516,
+ "INT",
+ 0,
+ "Center right ascension.",
+ 4
+ },
+ {
+ "deccol",
+ 517,
+ "INT",
+ 0,
+ "Center declination.",
+ 4
+ },
+ {
"fcol",
502,
"INT",
@@ -532,6 +548,16 @@ parse_opt(int key, char *arg, struct argp_state *state)
NULL, 0);
p->up.ycolset=1;
break;
+ case 516:
+ gal_checkset_sizet_el_zero(arg, &p->racol, "racol", ' ', p->cp.spack,
+ NULL, 0);
+ p->up.racolset=1;
+ break;
+ case 517:
+ gal_checkset_sizet_el_zero(arg, &p->deccol, "deccol", ' ', p->cp.spack,
+ NULL, 0);
+ p->up.deccolset=1;
+ break;
case 502:
gal_checkset_sizet_el_zero(arg, &p->fcol, "fcol", ' ', p->cp.spack,
NULL, 0);
diff --git a/src/mkprof/main.h b/src/mkprof/main.h
index ac08fcb..da853d5 100644
--- a/src/mkprof/main.h
+++ b/src/mkprof/main.h
@@ -113,6 +113,8 @@ struct uiparams
int fcolset;
int xcolset;
int ycolset;
+ int racolset;
+ int deccolset;
int rcolset;
int ncolset;
int pcolset;
@@ -158,6 +160,8 @@ struct mkprofparams
size_t fcol; /* Column specifying profile function. */
size_t xcol; /* X column of profile center. */
size_t ycol; /* Y column of profile center. */
+ size_t racol; /* RA column of profile center. */
+ size_t deccol; /* Dec column of profile center. */
size_t rcol; /* Effective radius of profile. */
size_t ncol; /* Sersic index column of profile. */
size_t pcol; /* Position angle column of profile. */
diff --git a/src/mkprof/mkprof.c b/src/mkprof/mkprof.c
index b50bed6..5ad6fe6 100644
--- a/src/mkprof/mkprof.c
+++ b/src/mkprof/mkprof.c
@@ -59,73 +59,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
-/**************************************************************/
-/************ Prepare WCS parameters *************/
-/**************************************************************/
-void
-preparewcs(struct mkprofparams *p)
-{
- int status;
- struct wcsprm wcs;
- long os=p->oversample;
-
- /* Initialize the structure (allocate all the 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. */
- 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 input array values. */
- 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*p->resolution/3600/p->oversample;
- wcs.pc[3]=p->resolution/3600/p->oversample;
- wcs.pc[1]=wcs.pc[2]=0.0f;
- wcs.cdelt[0]=wcs.cdelt[1]=1.0f;
- 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: */
- if( (status=wcsset(&wcs)) )
- error(EXIT_FAILURE, 0, "wcsset error %d: %s", status,
- wcs_errmsg[status]);
-
- /* Write the WCS structure to a header string. */
- if( (status=wcshdo(WCSHDO_safe, &wcs, &p->wcsnkeyrec, &p->wcsheader)) )
- error(EXIT_FAILURE, 0, "wcshdo error %d: %s", status,
- wcs_errmsg[status]);
-
- /* Free the allocated spaces by wcsini/wcsset: */
- if( (status=wcsfree(&wcs)) )
- error(EXIT_FAILURE, 0, "wcsfree error %d: %s", status,
- wcs_errmsg[status]);
-}
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -203,7 +136,7 @@ saveindividual(struct mkonthread *mkp)
crpix[0] = p->crpix[0] - os*(mkp->fpixel_i[0]-1);
crpix[1] = p->crpix[1] - os*(mkp->fpixel_i[1]-1);
- /* Write the name and save the FITS image: */
+ /* Write the name and remove a similarly named file. */
sprintf(outname, "%s%lu_%s", outdir, ibq->id, p->basename);
gal_checkset_check_remove_file(outname, p->cp.dontdelete);
@@ -614,27 +547,41 @@ write(struct mkprofparams *p)
ibq=tbq;
}
- /* Write the final array to the final FITS image. */
+ /* Write the final array to the output FITS image. */
if(p->nomerged==0)
{
+ /* Get the current time for verbose output. */
if(verb) gettimeofday(&t1, NULL);
+
+ /* If the background image had a special file type, MakeProfile's
+ output should also have the same type. So here we check the type
+ and if its not float, we change to that type. */
if(p->up.backname)
{
- if(bitpix==FLOAT_IMG) array=out;
- else gal_fits_change_type(p->out, FLOAT_IMG,
- p->naxes[1]*p->naxes[0],
- p->anyblank, &array, bitpix);
- gal_fits_array_to_file(p->mergedimgname, "MockImg on back",
- bitpix, array, p->naxes[1],
- p->naxes[0], p->anyblank, p->wcs,
- NULL, SPACK_STRING);
- if(bitpix!=FLOAT_IMG) free(array);
+ if(bitpix==FLOAT_IMG)
+ array=out;
+ else
+ gal_fits_change_type(out, FLOAT_IMG, p->naxes[1]*p->naxes[0],
+ p->anyblank, &array, bitpix);
}
else
- gal_fits_atof_correct_wcs(p->mergedimgname, "MockImg", FLOAT_IMG,
- out, p->naxes[1], p->naxes[0],
- p->wcsheader, p->wcsnkeyrec,
- NULL, SPACK_STRING);
+ {
+ array=out;
+ bitpix=FLOAT_IMG;
+ }
+
+ /* Write the array inside the output FITS file. */
+ gal_fits_array_to_file(p->mergedimgname, "MockImg on back",
+ bitpix, array, p->naxes[1],
+ p->naxes[0], p->anyblank, p->wcs,
+ NULL, SPACK_STRING);
+
+ /* Free `array' when it was not a simple assignment but an
+ allocation,*/
+ if(p->up.backname && bitpix!=FLOAT_IMG)
+ free(array);
+
+ /* In verbose mode, print the information. */
if(verb)
{
errno=0;
@@ -684,10 +631,6 @@ mkprof(struct mkprofparams *p)
size_t nt=p->cp.numthreads, nb;
long onaxes[2], os=p->oversample;
- /* Get the WCS header strings ready to put into the FITS
- image(s). */
- if(p->up.backname==NULL)
- preparewcs(p);
/* Allocate the arrays to keep the thread and parameters for each
thread. Note that we only want nt-1 threads to do the
@@ -766,5 +709,4 @@ mkprof(struct mkprofparams *p)
/* Free the allocated spaces. */
free(mkp);
free(indexs);
- free(p->wcsheader);
}
diff --git a/src/mkprof/ui.c b/src/mkprof/ui.c
index 08ad2bf..7efa550 100644
--- a/src/mkprof/ui.c
+++ b/src/mkprof/ui.c
@@ -29,9 +29,10 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdio.h>
#include <stdlib.h>
-#include <nproc.h> /* From Gnulib. */
+#include <nproc.h> /* From Gnulib. */
#include <gnuastro/box.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
#include <gnuastro/timing.h> /* Includes time.h and sys/time.h */
#include <gnuastro/checkset.h>
@@ -211,6 +212,20 @@ readconfig(char *filename, struct mkprofparams *p)
filename, lineno);
up->ycolset=1;
}
+ else if(strcmp(name, "racol")==0)
+ {
+ if(up->racolset) continue;
+ gal_checkset_sizet_el_zero(value, &p->racol, name, key, SPACK,
+ filename, lineno);
+ up->racolset=1;
+ }
+ else if(strcmp(name, "deccol")==0)
+ {
+ if(up->deccolset) continue;
+ gal_checkset_sizet_el_zero(value, &p->deccol, name, key, SPACK,
+ filename, lineno);
+ up->deccolset=1;
+ }
else if(strcmp(name, "fcol")==0)
{
if(up->fcolset) continue;
@@ -360,6 +375,10 @@ printvalues(FILE *fp, struct mkprofparams *p)
fprintf(fp, CONF_SHOWFMT"%lu\n", "xcol", p->xcol);
if(up->ycolset)
fprintf(fp, CONF_SHOWFMT"%lu\n", "ycol", p->ycol);
+ if(up->racolset)
+ fprintf(fp, CONF_SHOWFMT"%lu\n", "racol", p->racol);
+ if(up->deccolset)
+ fprintf(fp, CONF_SHOWFMT"%lu\n", "deccol", p->deccol);
if(up->fcolset)
fprintf(fp, CONF_SHOWFMT"%lu\n", "fcol", p->fcol);
if(up->rcolset)
@@ -414,10 +433,6 @@ checkifset(struct mkprofparams *p)
GAL_CONFIGFILES_REPORT_NOTSET("tolerance");
if(up->zeropointset==0)
GAL_CONFIGFILES_REPORT_NOTSET("zeropoint");
- if(up->xcolset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("xcol");
- if(up->ycolset==0)
- GAL_CONFIGFILES_REPORT_NOTSET("ycol");
if(up->fcolset==0)
GAL_CONFIGFILES_REPORT_NOTSET("fcol");
if(up->rcolset==0)
@@ -448,6 +463,27 @@ checkifset(struct mkprofparams *p)
GAL_CONFIGFILES_REPORT_NOTSET("crval2");
if(up->resolutionset==0)
GAL_CONFIGFILES_REPORT_NOTSET("resolution");
+
+
+ /* The X and Y columns are only needed when the RA and Dec columns have
+ not been given. */
+ if(up->racolset==0 && up->deccolset==0)
+ {
+ if(up->xcolset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("xcol");
+ if(up->ycolset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("ycol");
+ }
+ /* The `if' statment above made sure that at least one of the RA and Dec
+ columns have been specified. So make sure that both are specified. */
+ else if(up->racolset!=up->deccolset)
+ {
+ if(up->racolset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("racol");
+ if(up->deccolset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("deccol");
+ }
+
GAL_CONFIGFILES_END_OF_NOTSET_REPORT;
}
@@ -478,18 +514,23 @@ sanitycheck(struct mkprofparams *p)
{
int d0f1;
long width[2]={1,1};
- char *tmpname=NULL;
+ struct uiparams *up=&p->up;
double truncr, *cat=p->cat, *row;
size_t i, j, columns[9], cs1=p->cs1;
+ char *tmpname=NULL, *xcolstr, *ycolstr;
/* Make sure the input file exists. */
gal_checkset_check_file(p->up.catname);
- /* Check if over-sampling is an odd number, then convert the
- oversampling rate into the double type. */
- if(p->oversample%2==0) ++p->oversample;
+ /* Check if over-sampling is an odd number, then set/modify the
+ respective values.*/
+ if(p->oversample%2==0)
+ error(EXIT_FAILURE, 0, "the value to the `--oversample' (`-s') option "
+ "must be an odd number. Please run the following command for a "
+ "complete explanation:\n\n info gnuastro \"Oversampling\"\n\nOr "
+ "See the \"Oversampling\" section of the Gnuastro book.");
p->halfpixel = 0.5f/p->oversample;
p->naxes[0] *= p->oversample;
p->naxes[1] *= p->oversample;
@@ -503,6 +544,29 @@ sanitycheck(struct mkprofparams *p)
"one value. So these two options cannot be called together");
+
+ /* When the RA and Dec columns have been given use them for the profile
+ positions instead of the X and Y columns. In the next step we are
+ going to convert the RAs and Decs to Xs and Ys and until then, we are
+ just dealing with the columns, not the actual values, so it is safe
+ (and greatfly simplifies the sanity checks below) to set xcol to ra
+ column and ycol to deccol. Also use this check to set the string that
+ should be printed if the column is not within the catalog's number of
+ columsn*/
+ if(up->racolset)
+ {
+ p->xcol=p->racol;
+ p->ycol=p->deccol;
+ xcolstr="racol";
+ ycolstr="deccol";
+ }
+ else
+ {
+ xcolstr="xcol";
+ ycolstr="ycol";
+ }
+
+
/* If the column numbers are not equal. */
columns[0]=p->xcol; columns[1]=p->ycol; columns[2]=p->fcol;
columns[3]=p->rcol; columns[4]=p->ncol; columns[5]=p->pcol;
@@ -517,8 +581,8 @@ sanitycheck(struct mkprofparams *p)
/* If all the columns are within the catalog: */
- GAL_CHECKSET_CHECK_COL_IN_CAT(p->xcol, "xcol");
- GAL_CHECKSET_CHECK_COL_IN_CAT(p->ycol, "ycol");
+ GAL_CHECKSET_CHECK_COL_IN_CAT(p->xcol, xcolstr);
+ GAL_CHECKSET_CHECK_COL_IN_CAT(p->ycol, ycolstr);
GAL_CHECKSET_CHECK_COL_IN_CAT(p->fcol, "fcol");
GAL_CHECKSET_CHECK_COL_IN_CAT(p->rcol, "rcol");
GAL_CHECKSET_CHECK_COL_IN_CAT(p->ncol, "ncol");
@@ -644,10 +708,71 @@ sanitycheck(struct mkprofparams *p)
/*************** Preparations *******************/
/**************************************************************/
void
+preparewcs(struct mkprofparams *p)
+{
+ int status;
+ struct wcsprm *wcs;
+ long os=p->oversample;
+
+ /* Allocate the memory necessary for the wcsprm structure. */
+ errno=0;
+ wcs=p->wcs=malloc(sizeof *wcs);
+ if(wcs==NULL)
+ error(EXIT_FAILURE, errno, "%lu for wcs in preparewcs", sizeof *wcs);
+
+ /* Initialize the structure (allocate all the 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*p->resolution/3600/p->oversample;
+ wcs->pc[3]=p->resolution/3600/p->oversample;
+ wcs->pc[1]=wcs->pc[2]=0.0f;
+ wcs->cdelt[0]=wcs->cdelt[1]=1.0f;
+ 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]);
+
+ /* 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. */
+ if(p->individual)
+ {
+ status=wcshdo(WCSHDO_safe, wcs, &p->wcsnkeyrec, &p->wcsheader);
+ if(status)
+ error(EXIT_FAILURE, 0, "wcshdo error %d: %s", status,
+ wcs_errmsg[status]);
+ }
+}
+
+
+
+
+
+void
preparearrays(struct mkprofparams *p)
{
- void *array;
- size_t naxes[2];
+ void *array=NULL;
+ size_t i, naxes[2];
+ double *wcstoimg=NULL;
/* Allocate space for the log file: */
p->log=malloc(p->cs0*LOGNUMCOLS*sizeof *p->log);
@@ -693,6 +818,63 @@ preparearrays(struct mkprofparams *p)
}
+ /* If there is no WCS structure (either no background image given or the
+ background image didn't have any WCS structure), then prepare the WCS
+ based on the given options. */
+ if(p->wcs==NULL)
+ preparewcs(p);
+
+
+ /* Convert the RA and Dec to X and Y. We will make a temporary array with
+ four rows and the same number of columns, then use that array and the
+ WCS structure to get X and Y values. Those X and Y values will then be
+ replaced with the RAs and Decs of the original catalog. Recall that
+ this is only done in the RAM, not the actual input file, so there is
+ no problem. One of the reasons we are doing this is that WCSLIB needs
+ the X and Ys and the RA and Decs to be in touching pieces of memory
+ and we can't guarantee that (the user might have these columns in
+ any order). */
+ if(p->up.racolset)
+ {
+ /* Allocate the space for the temporary array. */
+ errno=0;
+ wcstoimg=malloc(4*p->cs0*sizeof *wcstoimg);
+ if(wcstoimg==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for wcstoimg",
+ 4*p->cs0*sizeof *wcstoimg);
+
+
+ /* Fill in the first two columns with the RA and Dec. */
+ for(i=0;i<p->cs0;++i)
+ {
+ wcstoimg[ i*4 ] = p->cat[ i*p->cs1+p->racol ];
+ wcstoimg[ i*4+1 ] = p->cat[ i*p->cs1+p->deccol ];
+ wcstoimg[ i*4+2 ] = wcstoimg[ i*4+3 ] = NAN;
+ }
+
+
+ /* Convert the X and Y to RA and Dec. */
+ gal_wcs_radec_array_to_xy(p->wcs, wcstoimg, wcstoimg+2, p->cs0, 4);
+
+
+ /* Write the produced X and Y into the input catalog, note that in
+ `sanitycheck', xcol became identical to racol and ycol to
+ deccol. Note that oversampling has been applied to the WCS
+ structure. However, when X and Y are given, oversampling is not
+ applied at this point, so we have to correct for the WCS's
+ oversampling.*/
+ for(i=0;i<p->cs0;++i)
+ {
+ p->cat[ i*p->cs1+p->xcol ] = wcstoimg[ i*4+2 ] / p->oversample;
+ p->cat[ i*p->cs1+p->ycol ] = wcstoimg[ i*4+3 ] / p->oversample;
+ }
+
+
+ /* Free the temporary array space. */
+ free(wcstoimg);
+ }
+
+
/* If the constant is to be NaN, then set it: */
if(p->setconsttonan) p->constant=CONSTFORNAN;
@@ -738,6 +920,7 @@ setparams(int argc, char *argv[], struct mkprofparams *p)
cp->removedirinfo = 1;
p->constant = 1;
+ p->wcs = NULL;
/* Read the arguments. */
errno=0;
@@ -828,6 +1011,8 @@ setparams(int argc, char *argv[], struct mkprofparams *p)
void
freeandreport(struct mkprofparams *p, struct timeval *t1)
{
+ int status;
+
/* Free all the allocated arrays. */
free(p->cat);
free(p->cp.hdu);
@@ -846,8 +1031,15 @@ freeandreport(struct mkprofparams *p, struct timeval *t1)
free(p->cp.output);
free(p->mergedimgname);
}
- if(p->up.backname)
- wcsvfree(&p->nwcs, &p->wcs);
+
+ /* Free the WCS headers string that was defined for individual mode. */
+ 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/tests/Makefile.am b/tests/Makefile.am
index e29a4b0..8720303 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -30,20 +30,18 @@ export topsrc=$(top_srcdir); \
export haslibjpeg=$(MAYBE_HASLIBJPEG); \
export hasghostscript=$(MAYBE_HASGHOSTSCRIPT);
-TESTS = prepconf.sh mkprof/mosaic1.sh mkprof/mosaic2.sh
\
-mkprof/mosaic3.sh mkprof/mosaic4.sh imgcrop/imgcat.sh \
-imgcrop/wcscat.sh imgcrop/xcyc.sh imgcrop/xcycnoblank.sh \
-imgcrop/section.sh imgcrop/radec.sh imgcrop/imgpolygon.sh \
-imgcrop/imgoutpolygon.sh imgcrop/wcspolygon.sh convertt/fitstotxt.sh \
-convertt/fitstojpeg.sh convertt/blankch.sh convertt/jpegtotxt.sh \
-convertt/fitstojpegcmyk.sh convertt/jpegtofits.sh \
-convertt/fitstopdf.sh convolve/spatial.sh convolve/frequency.sh
\
-imgwarp/imgwarp_scale.sh imgwarp/homographic.sh mknoise/addnoise.sh \
-mkprof/ellipticalmasks.sh header/write.sh header/print.sh \
-header/update.sh header/delete.sh imgstat/basicstats.sh
\
-subtractsky/subtractsky.sh noisechisel/noisechisel.sh \
-mkcatalog/simple.sh arithmetic/snimage.sh arithmetic/onlynumbers.sh \
-cosmiccal/simpletest.sh
+TESTS = prepconf.sh mkprof/mosaic1.sh mkprof/mosaic2.sh mkprof/mosaic3.sh \
+mkprof/mosaic4.sh mkprof/radeccat.sh imgcrop/imgcat.sh imgcrop/wcscat.sh \
+imgcrop/xcyc.sh imgcrop/xcycnoblank.sh imgcrop/section.sh imgcrop/radec.sh \
+imgcrop/imgpolygon.sh imgcrop/imgoutpolygon.sh imgcrop/wcspolygon.sh \
+convertt/fitstotxt.sh convertt/fitstojpeg.sh convertt/blankch.sh \
+convertt/jpegtotxt.sh convertt/fitstojpegcmyk.sh convertt/jpegtofits.sh \
+convertt/fitstopdf.sh convolve/spatial.sh convolve/frequency.sh
\
+imgwarp/imgwarp_scale.sh imgwarp/homographic.sh mknoise/addnoise.sh \
+mkprof/ellipticalmasks.sh header/write.sh header/print.sh header/update.sh \
+header/delete.sh imgstat/basicstats.sh subtractsky/subtractsky.sh \
+noisechisel/noisechisel.sh mkcatalog/simple.sh arithmetic/snimage.sh \
+arithmetic/onlynumbers.sh cosmiccal/simpletest.sh
EXTRA_DIST = $(TESTS) mkprof/mkprofcat1.txt \
mkprof/mkprofcat1_mask.txt mkprof/mkprofcat2.txt \
@@ -77,6 +75,7 @@ mkprof/mosaic1.sh: prepconf.sh.log
mkprof/mosaic2.sh: prepconf.sh.log
mkprof/mosaic3.sh: prepconf.sh.log
mkprof/mosaic4.sh: prepconf.sh.log
+mkprof/radeccat.sh: prepconf.sh.log
imgcrop/imgcat.sh: mkprof/mosaic1.sh.log
imgcrop/wcscat.sh: mkprof/mosaic1.sh.log mkprof/mosaic2.sh.log \
mkprof/mosaic3.sh.log mkprof/mosaic4.sh.log
diff --git a/tests/mkprof/radeccat.sh b/tests/mkprof/radeccat.sh
new file mode 100755
index 0000000..16c0686
--- /dev/null
+++ b/tests/mkprof/radeccat.sh
@@ -0,0 +1,49 @@
+# Create a mock image from radeccat.txt
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+# Mohammad Akhlaghi <address@hidden>
+# Contributing author(s):
+#
+# Copying and distribution of this file, with or without modification,
+# are permitted in any medium without royalty provided the copyright
+# notice and this notice are preserved. This file is offered as-is,
+# without any warranty.
+
+
+
+
+
+# Preliminaries
+# =============
+#
+# Set the variabels (The executable is in the build tree). Do the
+# basic checks to see if the executable is made or if the defaults
+# file exists (basicchecks.sh is in the source tree).
+prog=mkprof
+execname=../src/$prog/ast$prog
+cat=$topsrc/tests/$prog/radeccat.txt
+
+
+
+
+
+# Skip?
+# =====
+#
+# If the dependencies of the test don't exist, then skip it. There are two
+# types of dependencies:
+#
+# - The executable was not made (for example due to a configure option),
+#
+if [ ! -f $execname ]; then exit 77; fi
+
+
+
+
+
+# Actual test script
+# ==================
+$execname $cat --racol=1 --deccol=2 --naxis1=100 --naxis2=100
diff --git a/tests/mkprof/radeccat.txt b/tests/mkprof/radeccat.txt
new file mode 100644
index 0000000..90505b8
--- /dev/null
+++ b/tests/mkprof/radeccat.txt
@@ -0,0 +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.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