[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 6b59be2: Crop: added support for 3D crops usin
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 6b59be2: Crop: added support for 3D crops using --section and --center |
|
Date: |
Fri, 26 Jul 2019 06:45:24 -0400 (EDT) |
branch: master
commit 6b59be27d8d2bc0d388b4582cb4198fc159004e2
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Crop: added support for 3D crops using --section and --center
Until now, Crop only supported cropping of 2D images. With this commit it
also supports 3D crops, except for the `--polygon' mode (where it still
only works in 2D).
---
NEWS | 7 ++
bin/crop/astcrop.conf | 1 +
bin/crop/crop.c | 4 +-
bin/crop/main.h | 14 ++--
bin/crop/onecrop.c | 66 +++++++++++-----
bin/crop/onecrop.h | 32 ++++----
bin/crop/ui.c | 16 ++--
bin/crop/wcsmode.c | 207 ++++++++++++++++++++++++++++++++++++++++----------
8 files changed, 257 insertions(+), 90 deletions(-)
diff --git a/NEWS b/NEWS
index 022d54d..9756ca1 100644
--- a/NEWS
+++ b/NEWS
@@ -23,6 +23,13 @@ See the end of the file for license conditions.
dataset and returns a single-dimension output, containing only the
unique values in the dataset.
+ Crop:
+ - Can also crop 3D datasets (data cubes). A 3D crop has the same syntax
+ as the old 2D mode, only when the dataset is 3D, three coordinates
+ (values, ranges or catalog-columns) should be given to the relevant
+ option. Just note that `--polygon' crops are still not supported in
+ 3D.
+
CosmicCalculator:
--obsline: alternative way to set the used redshift. With this option
instead of explicitly giving the redshift, you can give a rest-frame
diff --git a/bin/crop/astcrop.conf b/bin/crop/astcrop.conf
index 8790a63..07f701d 100644
--- a/bin/crop/astcrop.conf
+++ b/bin/crop/astcrop.conf
@@ -31,6 +31,7 @@
# Crop by center (when a catalog is given)
coordcol 2
coordcol 3
+ coordcol 4
# Operating mode:
mode wcs
diff --git a/bin/crop/crop.c b/bin/crop/crop.c
index 64af887..3e19f71 100644
--- a/bin/crop/crop.c
+++ b/bin/crop/crop.c
@@ -420,8 +420,8 @@ crop(struct cropparams *p)
__func__, nt*sizeof *crp);
- /* Distribute the indexs into the threads (this is needed even if we
- only have one object where p->cs0 is not defined): */
+ /* Distribute the indexs into the threads (for clarity, this is needed
+ even if we only have one object). */
gal_threads_dist_in_threads(p->catname ? p->numout : 1, nt,
&indexs, &thrdcols);
diff --git a/bin/crop/main.h b/bin/crop/main.h
index 7b83426..528645e 100644
--- a/bin/crop/main.h
+++ b/bin/crop/main.h
@@ -40,7 +40,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Macros */
#define LOGFILENAME PROGRAM_EXEC".log"
#define FILENAME_BUFFER_IN_VERB 30
-#define MAXDIM 2
+#define MAXDIM 3
/* Modes to interpret coordinates. */
@@ -66,8 +66,8 @@ struct inputimgs
struct wcsprm *wcs; /* WCS structure of each input image. */
char *wcstxt; /* Text output of each WCS. */
int nwcskeys; /* Number of keywords in the header WCS. */
- double corners[8]; /* RA and Dec of this image corners (within). */
- double sized[2]; /* Width and height of image in degrees. */
+ double corners[24]; /* WCS of corners (24: for 3D, 8: for 2D). */
+ double sized[MAXDIM]; /* Width and height of image in degrees. */
double equatorcorr[2]; /* If image crosses the equator, see wcsmode.c.*/
};
@@ -93,7 +93,7 @@ struct cropparams
char *catname; /* Name of input catalog. */
char *cathdu; /* HDU of catalog if its a FITS file. */
char *namecol; /* Filename (without suffix) of crop col.*/
- gal_list_str_t *coordcol; /* Column in catalog with coordinates. */
+ gal_list_str_t *coordcol; /* Column in cat containing coordinates. */
char *section; /* Section string. */
char *polygon; /* Input string of polygon vertices. */
uint8_t outpolygon; /* ==1: Keep the inner polygon region. */
@@ -101,14 +101,14 @@ struct cropparams
/* Internal */
size_t numin; /* Number of input images. */
size_t numout; /* Number of output images. */
- double **centercoords; /* The center coordinates. */
+ double **centercoords; /* A 1D array for the center position. */
size_t checkcenter; /* width of a box to check for zeros */
char **name; /* filename of crop in row. */
double *wpolygon; /* Array of WCS polygon vertices. */
double *ipolygon; /* Array of image polygon vertices. */
size_t nvertices; /* Number of polygon vertices. */
- long iwidth[2]; /* Image mode width (in pixels). */
- double *pixscale; /* Resolution in each dimension. */
+ long iwidth[MAXDIM]; /* Image mode width (in pixels). */
+ double *pixscale; /* Raw resolution in each dimension. */
time_t rawtime; /* Starting time of the program. */
int outnameisfile; /* Output filename is a directory. */
int type; /* Type of output(s). */
diff --git a/bin/crop/onecrop.c b/bin/crop/onecrop.c
index bc5b30e..55c85cd 100644
--- a/bin/crop/onecrop.c
+++ b/bin/crop/onecrop.c
@@ -79,9 +79,10 @@ onecrop_parse_section(struct cropparams *p, size_t *dsize,
int add;
long read;
char *tailptr;
+ long naxes[MAXDIM];
char forl='f', *pt=p->section;
- long naxes[2]={dsize[1], dsize[0]};
- size_t i, dim=0, ndim=p->imgs->ndim;
+ size_t i, ndim=p->imgs->ndim, dim=0;
+
/* When the user asks for a section of the dataset, then the cropped
region is not defined by its center. So it makes no sense to later
@@ -89,6 +90,7 @@ onecrop_parse_section(struct cropparams *p, size_t *dsize,
zero. */
p->checkcenter=0;
+
/* Initialize the fpixel and lpixel arrays (note that `section' is only
defined in image mode, so there will only be one element in `imgs'. */
for(i=0;i<ndim;++i)
@@ -169,6 +171,7 @@ onecrop_parse_section(struct cropparams *p, size_t *dsize,
pt=tailptr;
}
+
/* Make sure the first pixel is located before/below the last pixel. */
for(i=0;i<ndim;++i)
if(fpixel[i]>lpixel[i])
@@ -298,7 +301,7 @@ onecrop_parse_polygon(struct cropparams *p)
-void
+static void
onecrop_ipolygon_fl(double *ipolygon, size_t nvertices, long *fpixel,
long *lpixel)
{
@@ -514,14 +517,14 @@ onecrop_flpixel(struct onecropparams *crp)
{
case IMGCROP_MODE_IMG:
- if(p->section) /* Defined by section. */
+ if(p->section) /* Defined by section. */
onecrop_parse_section(p, dsize, fpixel, lpixel);
- else if(p->polygon) /* Defined by polygon. */
+ else if(p->polygon) /* Defined by a polygon. */
{
if(p->outpolygon==0)
onecrop_ipolygon_fl(p->ipolygon, p->nvertices, fpixel, lpixel);
}
- else
+ else /* Defined by its center. */
{
for(i=0;i<ndim;++i) center[i] = p->centercoords[i][crp->out_ind];
gal_box_border_from_center(center, ndim, p->iwidth, fpixel, lpixel);
@@ -592,9 +595,9 @@ onecrop_make_array(struct onecropparams *crp, long
*fpixel_i,
long naxes[MAXDIM];
char *outname=crp->name;
int status=0, type=crp->p->type;
- size_t i, ndim=crp->p->imgs->ndim;
char **strarr, cpname[FLEN_KEYWORD];
gal_data_t *rkey=gal_data_array_calloc(1);
+ size_t i, ndim=crp->p->imgs->ndim, totsize;
char *cp, *cpf, blankrec[80], titlerec[80];
struct inputimgs *img=&crp->p->imgs[crp->in_ind];
@@ -671,7 +674,8 @@ onecrop_make_array(struct onecropparams *crp, long
*fpixel_i,
crp->p->blankptrwrite, "Pixels with no data.",
&status) )
gal_fits_io_error(status, "adding Blank");
- if(fits_write_null_img(ofp, 1, naxes[0]*naxes[1], &status))
+ totsize = naxes[0]*naxes[1] * (ndim==3?naxes[2]:1);
+ if(fits_write_null_img(ofp, 1, totsize, &status))
gal_fits_io_error(status, "writing null array");
@@ -738,7 +742,7 @@ onecrop(struct onecropparams *crp)
size_t i, j, cropsize=1, ndim=img->ndim;
char region[FLEN_VALUE], regionkey[FLEN_KEYWORD];
long fpixel_o[MAXDIM], lpixel_o[MAXDIM], inc[MAXDIM];
- long naxes[MAXDIM], fpixel_i[MAXDIM] , lpixel_i[MAXDIM];
+ long naxes[MAXDIM], fpixel_i[MAXDIM], lpixel_i[MAXDIM];
/* Fill the `naxes' and `inc' arrays. */
for(i=0;i<ndim;++i)
@@ -755,7 +759,6 @@ onecrop(struct onecropparams *crp)
memcpy(lpixel_i, crp->lpixel, ndim*sizeof *lpixel_i);
-
/* Find the overlap and apply it if there is any overlap. */
if( gal_box_overlap(naxes, fpixel_i, lpixel_i, fpixel_o, lpixel_o, ndim) )
{
@@ -767,7 +770,7 @@ onecrop(struct onecropparams *crp)
/* Allocate an array to keep the desired crop region, then read
- the desired pixels onto it. */
+ the desired pixels into it. */
status=0;
for(i=0;i<ndim;++i) cropsize *= ( lpixel_i[i] - fpixel_i[i] + 1 );
array=gal_pointer_allocate(p->type, cropsize, 0, __func__, "array");
@@ -790,6 +793,11 @@ onecrop(struct onecropparams *crp)
outside of it.*/
if(p->polygon)
{
+ /* A small sanity check until this part supports 3D. */
+ if(ndim!=2)
+ error(EXIT_FAILURE, 0, "%s: polygons only implemented in 2D",
+ __func__);
+
/* In WCS mode, crp->ipolygon was allocated and filled in
wcspolygonflpixel (wcsmode.c). */
if(p->mode==IMGCROP_MODE_IMG) crp->ipolygon=p->ipolygon;
@@ -868,15 +876,24 @@ onecrop_center_filled(struct onecropparams *crp)
fitsfile *ofp=crp->outfits;
int status=0, anynul=0, type;
long checkcenter=p->checkcenter;
- long naxes[2], fpixel[2], lpixel[2], inc[2]={1,1};
+ long naxes[3], fpixel[3], lpixel[3], inc[3]={1,1,1};
/* If checkcenter is zero, then don't check. */
if(checkcenter==0) return GAL_BLANK_UINT8;
/* Get the final size of the output image. */
gal_fits_img_info(ofp, &type, &ndim, &dsize, NULL, NULL);
- naxes[0]=dsize[1];
- naxes[1]=dsize[0];
+ if(ndim==2)
+ {
+ naxes[0]=dsize[1];
+ naxes[1]=dsize[0];
+ }
+ else
+ {
+ naxes[0]=dsize[2];
+ naxes[1]=dsize[1];
+ naxes[2]=dsize[0];
+ }
/* Get the size and range of the central region to check. The +1 is
because in FITS, counting begins from 1, not zero. It might happen
@@ -885,10 +902,23 @@ onecrop_center_filled(struct onecropparams *crp)
full image to check. */
size = ( (naxes[0]>checkcenter ? checkcenter : naxes[0])
* (naxes[1]>checkcenter ? checkcenter : naxes[1]) );
- fpixel[0] = naxes[0]>checkcenter ? (naxes[0]/2+1)-checkcenter/2 : 1;
- fpixel[1] = naxes[1]>checkcenter ? (naxes[1]/2+1)-checkcenter/2 : 1;
- lpixel[0] = naxes[0]>checkcenter ? (naxes[0]/2+1)+checkcenter/2 : naxes[0];
- lpixel[1] = naxes[1]>checkcenter ? (naxes[1]/2+1)+checkcenter/2 : naxes[1];
+ fpixel[0] = naxes[0]>checkcenter ? ((naxes[0]/2+1)-checkcenter/2) : 1;
+ fpixel[1] = naxes[1]>checkcenter ? ((naxes[1]/2+1)-checkcenter/2) : 1;
+ lpixel[0] = ( naxes[0]>checkcenter
+ ? ((naxes[0]/2+1)+checkcenter/2) : naxes[0] );
+ lpixel[1] = ( naxes[1]>checkcenter
+ ? ((naxes[1]/2+1)+checkcenter/2) : naxes[1] );
+
+
+ /* For the third dimension. */
+ if(ndim==3)
+ {
+ size *= (naxes[2]>checkcenter ? checkcenter : naxes[2]);
+ fpixel[2] = naxes[2]>checkcenter ? ((naxes[2]/2+1)-checkcenter/2) : 1;
+ lpixel[2] = ( naxes[2]>checkcenter
+ ? ((naxes[2]/2+1)+checkcenter/2) : naxes[2] );
+ }
+
/* For a check:
printf("naxes: %ld, %ld\nfpixel: (%ld, %ld)\nlpixel: (%ld, %ld)\n"
diff --git a/bin/crop/onecrop.h b/bin/crop/onecrop.h
index 7030c73..8c32101 100644
--- a/bin/crop/onecrop.h
+++ b/bin/crop/onecrop.h
@@ -35,28 +35,28 @@ struct onecropparams
struct cropparams *p;
/* About input image. */
- size_t in_ind; /* Index of this image in the input names. */
- fitsfile *infits; /* Pointer to the input FITS image. */
- long fpixel[2]; /* Position of first pixel in input image. */
- long lpixel[2]; /* Position of last pixel in input image. */
- double *ipolygon; /* Input image based polygon vertices. */
+ size_t in_ind; /* Index of this image in the input names. */
+ fitsfile *infits; /* Pointer to the input FITS image. */
+ long fpixel[MAXDIM]; /* Position of first pixel in input image. */
+ long lpixel[MAXDIM]; /* Position of last pixel in input image. */
+ double *ipolygon; /* Input image based polygon vertices. */
/* Output (cropped) image. */
- size_t out_ind; /* Index of this crop in the output list. */
- double world[2]; /* World coordinates of crop center. */
- double sized[2]; /* Width and height of image in degrees. */
- double corners[8]; /* RA and Dec of this crop's four sides. */
- double equatorcorr[2]; /* Crop crosses the equator, see wcsmode.c. */
- fitsfile *outfits; /* Pointer to the output FITS image. */
+ size_t out_ind; /* Index of this crop in the output list. */
+ double world[MAXDIM]; /* World coordinates of crop center. */
+ double sized[MAXDIM]; /* Width and height of image in degrees. */
+ double corners[24]; /* RA and Dec of this crop's corners. */
+ double equatorcorr[2]; /* Crop crosses the equator, see wcsmode.c. */
+ fitsfile *outfits; /* Pointer to the output FITS image. */
/* For log */
- char *name; /* Filename of crop. */
- size_t numimg; /* Number of images used to make this crop. */
- unsigned char centerfilled; /* ==1 if the center is filled. */
+ char *name; /* Filename of crop. */
+ size_t numimg; /* Number of images used to make this crop. */
+ unsigned char centerfilled; /* ==1 if the center is filled. */
/* Thread parameters. */
- size_t *indexs; /* Indexs to be used in this thread. */
- pthread_barrier_t *b; /* pthread barrier to keep threads waiting. */
+ size_t *indexs; /* Indexs to be used in this thread. */
+ pthread_barrier_t *b; /* pthread barrier to keep threads waiting. */
};
void
diff --git a/bin/crop/ui.c b/bin/crop/ui.c
index 6e8599c..939826b 100644
--- a/bin/crop/ui.c
+++ b/bin/crop/ui.c
@@ -66,7 +66,7 @@ const char *
argp_program_bug_address = PACKAGE_BUGREPORT;
static char
-args_doc[] = "[Crop-identifiers] ASTRdata ...";
+args_doc[] = "[Crop-Identifier] ASTRdata ...";
const char
doc[] = GAL_STRINGS_TOP_HELP_INFO PROGRAM_NAME" will create cutouts, "
@@ -359,9 +359,9 @@ ui_read_check_only_options(struct cropparams *p)
given, so we only need to check one of the two in each couple. */
if(p->coordcol==NULL)
error(EXIT_FAILURE, 0, "no crop center columns given to read from "
- "the input catalog (`%s'). Please use `--coordcol' two times "
- "to specify the column keeping the center position the "
- "respective dimension.\n\n"
+ "the input catalog (`%s'). Please use `--coordcol' several "
+ "times (depending on dimensionality) to specify the column "
+ "keeping the center position the respective dimension.\n\n"
"For more information on how to select columns in Gnuastro, "
"please run the following command:\n\n"
" $ info gnuastro \"Selecting table columns\"", p->catname);
@@ -797,7 +797,6 @@ ui_prepare_center(struct cropparams *p)
-
/* 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
@@ -938,6 +937,12 @@ ui_preparations(struct cropparams *p)
if(p->mode==IMGCROP_MODE_WCS) wcsmode_check_prepare(p, img);
}
+ /* Polygon cropping is currently only supported on 2D */
+ if(p->imgs->ndim!=2 && p->polygon)
+ error(EXIT_FAILURE, 0, "%s: polygon cropping is currently only "
+ "supported on 2D datasets (images), not %zuD datasets",
+ p->imgs->name, p->imgs->ndim);
+
/* Unify central crop methods into `p->centercoords'. */
if(p->catname || p->center)
@@ -1057,7 +1062,6 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct
cropparams *p)
gal_timing_report(NULL, msg, 1);
}
}
-
}
diff --git a/bin/crop/wcsmode.c b/bin/crop/wcsmode.c
index 8f6fe34..f57ef3d 100644
--- a/bin/crop/wcsmode.c
+++ b/bin/crop/wcsmode.c
@@ -56,9 +56,13 @@ wcsmode_check_prepare(struct cropparams *p, struct inputimgs
*img)
{
double *pixscale;
struct wcsprm *wcs=img->wcs;
- int i, status[4]={0,0,0,0}, ncorners=4;
size_t *dsize=img->dsize, ndim=img->ndim;
- double imgcrd[8], phi[4], theta[4], pixcrd[8];
+
+ /* For two dimensions, we have four corners (2 numbers for each) and for
+ three dimensions we have 8 corners (3 numbers for each), so we'll just
+ assume the largest. */
+ int i, status[8]={0,0,0,0,0,0,0,0}, ncorners=0;
+ double imgcrd[24], phi[8], theta[8], pixcrd[24];
/* Check if the image is aligned with the WCS coordinates. Note that
@@ -125,10 +129,30 @@ wcsmode_check_prepare(struct cropparams *p, struct
inputimgs *img)
/* Set the coordinates of the dataset's corners. Note that `dsize' is in
C order, while pixcrd is in FITS order.*/
- pixcrd[0] = 1; pixcrd[1] = 1;
- pixcrd[2] = dsize[1]; pixcrd[3] = 1;
- pixcrd[4] = 1; pixcrd[5] = dsize[0];
- pixcrd[6] = dsize[1]; pixcrd[7] = dsize[0];
+ switch(ndim)
+ {
+ case 2:
+ ncorners=4;
+ pixcrd[0] = 1; pixcrd[1] = 1;
+ pixcrd[2] = dsize[1]; pixcrd[3] = 1;
+ pixcrd[4] = 1; pixcrd[5] = dsize[0];
+ pixcrd[6] = dsize[1]; pixcrd[7] = dsize[0];
+ break;
+ case 3:
+ ncorners=8;
+ pixcrd[0] = 1; pixcrd[1] = 1; pixcrd[2] = 1;
+ pixcrd[3] = dsize[2]; pixcrd[4] = 1; pixcrd[5] = 1;
+ pixcrd[6] = 1; pixcrd[7] = dsize[1]; pixcrd[8] = 1;
+ pixcrd[9] = dsize[2]; pixcrd[10] = dsize[1]; pixcrd[11] = 1;
+ pixcrd[12] = 1; pixcrd[13] = 1; pixcrd[14] = dsize[0];
+ pixcrd[15] = dsize[2]; pixcrd[16] = 1; pixcrd[17] = dsize[0];
+ pixcrd[18] = 1; pixcrd[19] = dsize[1]; pixcrd[20] = dsize[0];
+ pixcrd[21] = dsize[2]; pixcrd[22] = dsize[1]; pixcrd[23] = dsize[0];
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: %zu dimensional datasets not supported",
+ __func__, ndim);
+ }
/* Get the coordinates of the corners of the dataset in WCS. */
@@ -146,9 +170,19 @@ wcsmode_check_prepare(struct cropparams *p, struct
inputimgs *img)
/* Fill in the size of the dataset in WCS from the first pixel in the
image. Note that `dsize' is in C axises, while the `pixscale',
`corners' and `sized' are in FITS axises. */
- img->sized[0] = ( img->dsize[1] * p->pixscale[0]
- / cos( img->corners[1] * M_PI / 180 ) );
- img->sized[1] = img->dsize[0] * p->pixscale[1];
+ if(ndim==2)
+ {
+ img->sized[0] = ( img->dsize[1] * p->pixscale[0]
+ / cos( img->corners[1] * M_PI / 180 ) );
+ img->sized[1] = img->dsize[0] * p->pixscale[1];
+ }
+ else /* 3D */
+ {
+ img->sized[0] = ( img->dsize[2] * p->pixscale[0] /* RA */
+ / cos( img->corners[1] * M_PI / 180 ) );
+ img->sized[1] = img->dsize[1] * p->pixscale[1]; /* Dec */
+ img->sized[2] = img->dsize[0] * p->pixscale[2]; /* 3D */
+ }
/* In case the image crosses the equator, we will calculate these values
@@ -170,14 +204,32 @@ wcsmode_check_prepare(struct cropparams *p, struct
inputimgs *img)
/* Just to check:
printf("\n\n%s:\n", img->name);
- printf("(%.10f, %.10f)\n"
- "(%.10f, %.10f)\n"
- "(%.10f, %.10f)\n"
- "(%.10f, %.10f)\n\n",
- img->corners[0], img->corners[1],
- img->corners[2], img->corners[3],
- img->corners[4], img->corners[5],
- img->corners[6], img->corners[7]);
+ if(ndim==2)
+ printf("(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n"
+ "(%.10f, %.10f)\n\n",
+ img->corners[0], img->corners[1],
+ img->corners[2], img->corners[3],
+ img->corners[4], img->corners[5],
+ img->corners[6], img->corners[7]);
+ else
+ printf("(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n"
+ "(%.10f, %.10f, %.10f)\n\n",
+ img->corners[0], img->corners[1], img->corners[2],
+ img->corners[3], img->corners[4], img->corners[5],
+ img->corners[6], img->corners[7], img->corners[8],
+ img->corners[9], img->corners[10], img->corners[11],
+ img->corners[12], img->corners[13], img->corners[14],
+ img->corners[15], img->corners[16], img->corners[17],
+ img->corners[18], img->corners[19], img->corners[20],
+ img->corners[21], img->corners[22], img->corners[23] );
exit(0);
*/
}
@@ -219,12 +271,18 @@ wcsmode_crop_corners(struct onecropparams *crp)
size_t i, ndim=p->imgs->ndim;
double minra=FLT_MAX, mindec=FLT_MAX;
double maxra=-FLT_MAX, maxdec=-FLT_MAX;
- double r, d, dr, h[MAXDIM], hr[MAXDIM];
+ double r, d, l, dr, h[MAXDIM], hr[MAXDIM];
size_t rmini=-1, rmaxi=-1, dmini=-1, dmaxi=-1;
/* Set the four corners of the WCS region. */
if(p->polygon)
{
+ /* A small sanity check. */
+ if(ndim!=2)
+ error(EXIT_FAILURE, 0, "%s: polygon crops are currently only "
+ "supported on 2D datasets, the input dataset is %zuD",
+ __func__, ndim);
+
/* Find their minimum and maximum values. */
for(i=0;i<p->nvertices;++i)
{
@@ -245,6 +303,7 @@ wcsmode_crop_corners(struct onecropparams *crp)
/* Set the RA and Dec to use as center. */
r=crp->world[0]=p->centercoords[0][crp->out_ind];
d=crp->world[1]=p->centercoords[1][crp->out_ind];
+ if(ndim==3) l=crp->world[2]=p->centercoords[2][crp->out_ind];
/* Calculate the declination in radians for easy readability. */
@@ -254,21 +313,56 @@ wcsmode_crop_corners(struct onecropparams *crp)
also calculate it in radians. */
hr[0] = ( h[0] = ((double *)(p->width->array))[0] / 2 ) * M_PI / 180;
hr[1] = ( h[1] = ((double *)(p->width->array))[1] / 2 ) * M_PI / 180;
+ if(ndim==3)
+ h[2] = ((double *)(p->width->array))[2] / 2;
/* Set the corners of this crop. */
- crp->corners[0] = r+h[0]/cos(dr-hr[1]);
- crp->corners[1] = d-h[1]; /* Bottom left. */
-
- crp->corners[2] = r-h[0]/cos(dr-hr[1]);
- crp->corners[3] = d-h[1]; /* Bottom Right. */
-
- crp->corners[4] = r+h[0]/cos(dr+hr[1]);
- crp->corners[5] = d+h[1]; /* Top Left. */
-
- crp->corners[6] = r-h[0]/cos(dr+hr[1]);
- crp->corners[7] = d+h[1]; /* Top Right. */
+ switch(ndim)
+ {
+ case 2:
+ crp->corners[0] = r+h[0]/cos(dr-hr[1]);
+ crp->corners[1] = d-h[1]; /* Bottom left. */
+
+ crp->corners[2] = r-h[0]/cos(dr-hr[1]);
+ crp->corners[3] = d-h[1]; /* Bottom Right. */
+
+ crp->corners[4] = r+h[0]/cos(dr+hr[1]);
+ crp->corners[5] = d+h[1]; /* Top Left. */
+
+ crp->corners[6] = r-h[0]/cos(dr+hr[1]);
+ crp->corners[7] = d+h[1]; /* Top Right. */
+ break;
+
+ case 3:
+ /* Note that the third dimension is assumed to be independent of
+ the first two. So the first two coordinates of its corners in
+ the front and back (on the two faces in the third dimension),
+ are equal. */
+ crp->corners[0] = crp->corners[12] = r+h[0]/cos(dr-hr[1]);
+ crp->corners[1] = crp->corners[13] = d-h[1];
+ crp->corners[2] = l-h[2]; /* Bottom left front. */
+
+ crp->corners[3] = crp->corners[15] = r-h[0]/cos(dr-hr[1]);
+ crp->corners[4] = crp->corners[16] = d-h[1];
+ crp->corners[5] = l-h[2]; /* Bottom right front.*/
+
+ crp->corners[6] = crp->corners[18] = r+h[0]/cos(dr+hr[1]);
+ crp->corners[7] = crp->corners[19] = d+h[1];
+ crp->corners[8] = l-h[2]; /* Top Left front. */
+
+ crp->corners[9] = crp->corners[21] = r-h[0]/cos(dr+hr[1]);
+ crp->corners[10] = crp->corners[22] = d+h[1];
+ crp->corners[11] = l-h[2]; /* Top right front. */
+
+ crp->corners[14] = l+h[2]; /* Bottom left back. */
+ crp->corners[17] = l+h[2]; /* Bottom right back. */
+ crp->corners[20] = l+h[2]; /* Top left back. */
+ crp->corners[23] = l+h[2]; /* Top right back. */
+ break;
+ }
}
+
/* Set the bottom width and height of the crop in degrees. Note that the
width changes as the height changes, so here we want the height and
the lowest declination. Note that in 2D on the bottom edge, corners[0]
@@ -282,10 +376,12 @@ wcsmode_crop_corners(struct onecropparams *crp)
rmini = ndim; /* First element in second corner. */
rmaxi = 0; /* First element. */
dmini = 1; /* Second element. */
- dmaxi = 5; /* Second element in third corner. */
+ dmaxi = ndim==2 ? 5 : 7; /* Second element in third corner. */
crp->sized[0]=( (crp->corners[rmaxi]-crp->corners[rmini])
/ cos(crp->corners[dmini]*M_PI/180) );
crp->sized[1]=crp->corners[dmaxi]-crp->corners[dmini];
+ if(ndim==3)
+ crp->sized[2] = crp->corners[14] - crp->corners[2];
/* In case the crop crosses the equator, then we need these two
@@ -300,16 +396,40 @@ wcsmode_crop_corners(struct onecropparams *crp)
crp->equatorcorr[1]=crp->sized[0]*cos(crp->corners[1]*M_PI/180);
}
+
/* Just to check:
- printf("\n\n%g, %g:\n", r, d);
- printf("\t(%.10f, %.10f)\n"
- "\t(%.10f, %.10f)\n"
- "\t(%.10f, %.10f)\n"
- "\t(%.10f, %.10f)\n\n",
- crp->corners[0], crp->corners[1],
- crp->corners[2], crp->corners[3],
- crp->corners[4], crp->corners[5],
- crp->corners[6], crp->corners[7]);
+ if(ndim==2)
+ {
+ printf("\n\n%g, %g:\n", r, d);
+ printf("\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n"
+ "\t(%.10f, %.10f)\n\n",
+ crp->corners[0], crp->corners[1],
+ crp->corners[2], crp->corners[3],
+ crp->corners[4], crp->corners[5],
+ crp->corners[6], crp->corners[7]);
+ }
+ else
+ {
+ printf("\n\n%g, %g, %g:\n", r, d, l);
+ printf("\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n"
+ "\t(%.10f, %.10f, %g)\n\n",
+ crp->corners[0], crp->corners[1], crp->corners[2],
+ crp->corners[3], crp->corners[4], crp->corners[5],
+ crp->corners[6], crp->corners[7], crp->corners[8],
+ crp->corners[9], crp->corners[10], crp->corners[11],
+ crp->corners[12], crp->corners[13], crp->corners[14],
+ crp->corners[15], crp->corners[16], crp->corners[17],
+ crp->corners[18], crp->corners[19], crp->corners[20],
+ crp->corners[21], crp->corners[22], crp->corners[23] );
+ }
exit(0);
*/
}
@@ -481,6 +601,11 @@ point_in_dataset(double *p, double *i, double *s, double
*c, size_t ndim)
{
double n;
+ /* If there is a third dimension, then first check that. Note that the
+ third dimension is assumed to be indendent of the first two. */
+ if(ndim==3 && ( p[2]<i[2] || p[2]>i[2]+s[2] ) )
+ return 0;
+
/* In the RA and Dec checks, first check the declination. If it is not in
range, you can safely return 0. */
if(p[1]>=i[1] && p[1]<=i[1]+s[1])
@@ -538,7 +663,7 @@ wcsmode_overlap(struct onecropparams *crp)
s=p->imgs[crp->in_ind].sized;
i=p->imgs[crp->in_ind].corners;
c=p->imgs[crp->in_ind].equatorcorr;
- fd=(d=crp->corners) + 8;
+ fd=(d=crp->corners) + (ndim==2 ? 8 : 24);
do
{
/* As long as one of the crop corners are in the image, we know there
@@ -554,7 +679,7 @@ wcsmode_overlap(struct onecropparams *crp)
s=crp->sized;
i=crp->corners;
c=crp->equatorcorr;
- fd=(d=p->imgs[crp->in_ind].corners) + 8;
+ fd=(d=p->imgs[crp->in_ind].corners) + (ndim==2 ? 8 : 24);
do
{
if( point_in_dataset(d, i, s, c, ndim) ) return 1;
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 6b59be2: Crop: added support for 3D crops using --section and --center,
Mohammad Akhlaghi <=