[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master f50cc13: Query: new --overlapwith option to au
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master f50cc13: Query: new --overlapwith option to automatically determine region |
|
Date: |
Sat, 19 Sep 2020 23:39:21 -0400 (EDT) |
branch: master
commit f50cc1366072ae43bae53d55f469e9a617b98066
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Query: new --overlapwith option to automatically determine region
Until now, if you wanted to query a certain database based on an image, you
would have to run 'astfits' with '--skycoverage', then pass the values. But
this is very inconvenient (especially when you are on the command-line and
want a fast result).
With this commit, the main workhorse of the '--skycoverage' option of
'astfits' has been moved inside the WCS library (and is used by
'astfits'). But because its in the library, it can now be automatically
used by 'astquery' too with the new '--overlapwith' option.
---
NEWS | 1 +
bin/fits/fits.c | 148 ++++++++------------------------------------------
bin/query/args.h | 13 +++++
bin/query/main.h | 1 +
bin/query/query.c | 47 +++++++++-------
bin/query/ui.c | 15 ++++--
bin/query/ui.h | 3 +-
doc/gnuastro.texi | 16 ++++++
lib/fits.c | 1 -
lib/gnuastro/wcs.h | 5 +-
lib/wcs.c | 154 +++++++++++++++++++++++++++++++++++++++++++++++++++++
11 files changed, 252 insertions(+), 152 deletions(-)
diff --git a/NEWS b/NEWS
index 4449490..e33e2bf 100644
--- a/NEWS
+++ b/NEWS
@@ -38,6 +38,7 @@ See the end of the file for license conditions.
within the image.
Library:
+ - gal_wcs_coverage: Return the sky coverage of given image HDU.
- gal_wcs_dimension_name: return the name of the requested WCS dim.
** Removed features
diff --git a/bin/fits/fits.c b/bin/fits/fits.c
index ebd000e..90a9348 100644
--- a/bin/fits/fits.c
+++ b/bin/fits/fits.c
@@ -357,131 +357,24 @@ fits_pixelscale(struct fitsparams *p)
static void
fits_skycoverage(struct fitsparams *p)
{
- fitsfile *fptr;
+ int nwcs=0;
+ size_t i, ndim;
struct wcsprm *wcs;
- int nwcs=0, type, status=0;
- char *name=NULL, *unit=NULL;
- gal_data_t *tmp, *coords=NULL;
- double *x, *y, *z, min[3], max[3];
- size_t i, ndim, numrows, *dsize=NULL;
-
- /* Read the desired WCS. */
- wcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
-
- /* If a WCS doesn't exist, let the user know and return. */
- if(wcs==NULL)
- error(EXIT_FAILURE, 0, "%s (hdu %s): no WCS could be read by WCSLIB, "
- "hence the sky coverage cannot be determined", p->filename,
- p->cp.hdu);
-
- /* Make sure the input HDU is an image. */
- if( gal_fits_hdu_format(p->filename, p->cp.hdu) != IMAGE_HDU )
- error(EXIT_FAILURE, 0, "%s (hdu %s): is not an image HDU, the "
- "'--skycoverage' option only applies to image extensions",
- p->filename, p->cp.hdu);
-
- /* Get the array information of the image. */
- fptr=gal_fits_hdu_open(p->filename, p->cp.hdu, READONLY);
- gal_fits_img_info(fptr, &type, &ndim, &dsize, &name, &unit);
- fits_close_file(fptr, &status);
+ double *center, *width, *min, *max;
- /* Abort if we have more than 3 dimensions. */
- if(ndim==1 || ndim>3)
- error(EXIT_FAILURE, 0, "%s (hdu: %s): has %zu dimensions. Currently "
- "'--skycoverage' only supports 2 or 3 dimensional datasets",
- p->filename, p->cp.hdu, ndim);
-
- /* Now that we have the number of dimensions in the image, allocate the
- space needed for the coordinates. */
- numrows = (ndim==2) ? 5 : 9;
- for(i=0;i<ndim;++i)
- {
- tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numrows, NULL, 0,
- p->cp.minmapsize, p->cp.quietmmap, NULL,
- NULL, NULL);
- tmp->next=coords;
- coords=tmp;
- }
-
- /* Fill in the coordinate arrays, Note that 'dsize' is ordered in C
- dimensions, for the WCS conversion, we need to have the dimensions
- ordered in FITS/Fortran order. */
- switch(ndim)
- {
- case 2:
- x=coords->array; y=coords->next->array;
- x[0] = 1; y[0] = 1;
- x[1] = dsize[1]; y[1] = 1;
- x[2] = 1; y[2] = dsize[0];
- x[3] = dsize[1]; y[3] = dsize[0];
- x[4] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
- y[4] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
- break;
- case 3:
- x=coords->array; y=coords->next->array; z=coords->next->next->array;
- x[0] = 1; y[0] = 1; z[0]=1;
- x[1] = dsize[2]; y[1] = 1; z[1]=1;
- x[2] = 1; y[2] = dsize[1]; z[2]=1;
- x[3] = dsize[2]; y[3] = dsize[1]; z[3]=1;
- x[4] = 1; y[4] = 1; z[4]=dsize[0];
- x[5] = dsize[2]; y[5] = 1; z[5]=dsize[0];
- x[6] = 1; y[6] = dsize[1]; z[6]=dsize[0];
- x[7] = dsize[2]; y[7] = dsize[1]; z[7]=dsize[0];
- x[8] = dsize[2]/2 + (dsize[2]%2 ? 1 : 0.5f);
- y[8] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
- z[8] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
- break;
- default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
- "fix the problem. 'ndim' of %zu is not recognized.",
- __func__, PACKAGE_BUGREPORT, ndim);
- }
-
- /* For a check:
- printf("IMAGE COORDINATES:\n");
- for(i=0;i<numrows;++i)
- if(ndim==2)
- printf("%-15g%-15g\n", x[i], y[i]);
- else
- printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
- */
-
- /* Convert to the world coordinate system. */
- gal_wcs_img_to_world(coords, wcs, 1);
-
- /* For a check:
- printf("\nWORLD COORDINATES:\n");
- for(i=0;i<numrows;++i)
- if(ndim==2)
- printf("%-15g%-15g\n", x[i], y[i]);
- else
- printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
- */
-
- /* Get the minimum and maximum values in each dimension. */
- tmp=gal_statistics_minimum(coords);
- min[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- tmp=gal_statistics_maximum(coords);
- max[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- tmp=gal_statistics_minimum(coords->next);
- min[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- tmp=gal_statistics_maximum(coords->next);
- max[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- if(ndim>2)
- {
- tmp=gal_statistics_minimum(coords->next->next);
- min[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- tmp=gal_statistics_maximum(coords->next->next);
- max[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
- }
+ /* Find the coverage. */
+ if( gal_wcs_coverage(p->filename, p->cp.hdu, &ndim,
+ ¢er, &width, &min, &max)==0 )
+ error(EXIT_FAILURE, 0, "%s (hdu %s): is not usable for finding "
+ "sky coverage (either doesn't have a WCS, or isn't an image "
+ "or cube HDU with 2 or 3 dimensions", p->filename, p->cp.hdu);
/* Inform the user. */
if(p->cp.quiet)
{
/* First print the center and full-width. */
- for(tmp=coords; tmp!=NULL; tmp=tmp->next)
- printf("%-15.10g ", ((double *)(tmp->array))[ndim==2 ? 4 : 8]);
- for(i=0;i<ndim;++i) printf("%-15.10g ", max[i]-min[i]);
+ for(i=0;i<ndim;++i) printf("%-15.10g ", center[i]);
+ for(i=0;i<ndim;++i) printf("%-15.10g ", width[i]);
printf("\n");
/* Then print the range in coordinates. */
@@ -495,13 +388,14 @@ fits_skycoverage(struct fitsparams *p)
switch(ndim)
{
case 2:
- printf(" Center: %-15.10g%-15.10g\n", x[4], y[4]);
- printf(" Width: %-15.10g%-15.10g\n", max[0]-min[0], max[1]-min[1]);
+ printf(" Center: %-15.10g%-15.10g\n", center[0], center[1]);
+ printf(" Width: %-15.10g%-15.10g\n", width[0], width[1]);
break;
case 3:
- printf(" Center: %-15.10g%-15.10g%-15.10g\n", x[8], y[8], z[8]);
- printf(" width: %-15.10g%-15.10g%-15.10g\n", max[0]-min[0],
- max[1]-min[1], max[2]-min[2]);
+ printf(" Center: %-15.10g%-15.10g%-15.10g\n", center[0],
+ center[1], center[2]);
+ printf(" width: %-15.10g%-15.10g%-15.10g\n", width[0],
+ width[1], width[2]);
break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
@@ -510,15 +404,19 @@ fits_skycoverage(struct fitsparams *p)
}
/* For the range type of coverage. */
+ wcs=gal_wcs_read(p->filename, p->cp.hdu, 0, 0, &nwcs);
printf("\nSky coverage by range along dimensions:\n");
for(i=0;i<ndim;++i)
printf(" %-8s %-15.10g%-15.10g\n", gal_wcs_dimension_name(wcs, i),
min[i], max[i]);
+ wcsfree(wcs);
}
/* Clean up. */
- free(dsize);
- wcsfree(wcs);
+ free(min);
+ free(max);
+ free(width);
+ free(center);
}
diff --git a/bin/query/args.h b/bin/query/args.h
index 4832af9..b18a02e 100644
--- a/bin/query/args.h
+++ b/bin/query/args.h
@@ -84,6 +84,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET,
},
{
+ "overlapwith",
+ UI_KEY_OVERLAPWITH,
+ "STR",
+ 0,
+ "Set query region to overlap with this image.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->overlapwith,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ },
+ {
"center",
UI_KEY_CENTER,
"FLT[,...]",
diff --git a/bin/query/main.h b/bin/query/main.h
index f763300..f716ea3 100644
--- a/bin/query/main.h
+++ b/bin/query/main.h
@@ -46,6 +46,7 @@ struct queryparams
struct gal_options_common_params cp; /* Common parameters. */
int database; /* ID of database to use. */
char *datasetstr; /* ID of dataset in database to use. */
+ char *overlapwith; /* Image to use instead of center. */
gal_data_t *center; /* Center position of query. */
gal_data_t *radius; /* Radius around center. */
gal_data_t *range; /* Range of magnitudes to query. */
diff --git a/bin/query/query.c b/bin/query/query.c
index 2a86bce..94a488e 100644
--- a/bin/query/query.c
+++ b/bin/query/query.c
@@ -28,6 +28,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <string.h>
+#include <gnuastro/wcs.h>
#include <gnuastro/pointer.h>
#include <gnuastro-internal/checkset.h>
@@ -82,15 +83,16 @@ void
query_gaia_sanitycheck(struct queryparams *p)
{
/* Make sure that atleast one type of constraint is specified. */
- if(p->center==NULL && p->query==NULL)
- error(EXIT_FAILURE, 0, "no '--center' or '--query' specified. At least "
- "one of these options are necessary in the Gaia dataset");
+ if(p->query==NULL && p->center==NULL && p->overlapwith==NULL)
+ error(EXIT_FAILURE, 0, "no '--query', '--center' or '--overlapwith' "
+ "specified. At least one of these options are necessary in the "
+ "Gaia dataset");
/* If '--center' is given, '--radius' is also necessary. */
- if(p->center)
+ if(p->center || p->overlapwith)
{
/* Make sure the radius is given, and that it isn't zero. */
- if( p->radius==NULL && p->width==NULL)
+ if(p->overlapwith==NULL && p->radius==NULL && p->width==NULL)
error(EXIT_FAILURE, 0, "the '--radius' ('-r') or '--width' ('-w') "
"options are necessary with the '--center' ('-C') option");
@@ -99,7 +101,6 @@ query_gaia_sanitycheck(struct queryparams *p)
error(EXIT_FAILURE, 0, "the '--dataset' ('-s') option is necessary "
"with the '--center' ('-C') option");
-
/* Use simpler names for the commonly used datasets. */
if( !strcmp(p->datasetstr, "dr2") )
{
@@ -131,10 +132,12 @@ query_gaia_sanitycheck(struct queryparams *p)
void
query_gaia(struct queryparams *p)
{
+ size_t ndim;
gal_data_t *tmp;
- double *center, *darray;
+ double width2, *center, *darray;
char *tmpstr, *regionstr, *rangestr=NULL;
char *command, *columns, allcols[]="*", *querystr;
+ double *ocenter=NULL, *owidth=NULL, *omin=NULL, *omax=NULL;
/* Make sure everything is fine. */
query_gaia_sanitycheck(p);
@@ -145,16 +148,18 @@ query_gaia(struct queryparams *p)
querystr=p->query;
else
{
- /* For easy reading. */
- center=p->center->array;
-
/* If certain columns have been requested use them, otherwise
- download all existing columns.
-
- columns="source_id,ra,dec,phot_g_mean_mag";
- */
+ download all existing columns.*/
columns = p->columns ? query_strlist_to_str(p->columns) : allcols;
+ /* If the user wanted an overlap with an image, then calculate it. */
+ if(p->overlapwith)
+ gal_wcs_coverage(p->overlapwith, p->cp.hdu, &ndim, &ocenter,
+ &owidth, &omin, &omax);
+
+ /* For easy reading. */
+ center = p->overlapwith ? ocenter : p->center->array;
+
/* Write the region. */
if(p->radius)
{
@@ -164,12 +169,13 @@ query_gaia(struct queryparams *p)
error(EXIT_FAILURE, 0, "%s: asprintf allocation ('regionstr')",
__func__);
}
- else if(p->width)
+ else if(p->width || p->overlapwith)
{
- darray=p->width->array;
+ darray = p->overlapwith ? owidth : p->width->array;
+ width2 = ( (p->overlapwith || p->width->size==2)
+ ? darray[1] : darray[0] );
if( asprintf( ®ionstr, "BOX('ICRS', %.8f, %.8f, %.8f, %.8f)",
- center[0], center[1], darray[0],
- p->width->size==1 ? darray[0] : darray[1] )<0 )
+ center[0], center[1], darray[0], width2 )<0 )
error(EXIT_FAILURE, 0, "%s: asprintf allocation ('regionstr')",
__func__);
}
@@ -195,11 +201,14 @@ query_gaia(struct queryparams *p)
"WHERE 1=CONTAINS( POINT('ICRS', ra, dec), %s ) %s",
columns, p->datasetstr, regionstr,
rangestr ? rangestr : "")<0 )
- error(EXIT_FAILURE, 0, "%s: asprintf allocation ('querystr')",
__func__);
+ error(EXIT_FAILURE, 0, "%s: asprintf allocation ('querystr')",
+ __func__);
/* Clean up. */
free(regionstr);
if(columns!=allcols) free(columns);
+ if(p->overlapwith)
+ {free(ocenter); free(owidth); free(omin); free(omax);}
}
diff --git a/bin/query/ui.c b/bin/query/ui.c
index 27ff1ea..13962b8 100644
--- a/bin/query/ui.c
+++ b/bin/query/ui.c
@@ -116,7 +116,6 @@ ui_initialize_options(struct queryparams *p,
/* Select individually. */
switch(cp->coptions[i].key)
{
- case GAL_OPTIONS_KEY_HDU:
case GAL_OPTIONS_KEY_LOG:
case GAL_OPTIONS_KEY_TYPE:
case GAL_OPTIONS_KEY_SEARCHIN:
@@ -269,9 +268,15 @@ ui_read_check_only_options(struct queryparams *p)
"command) for the current databases");
/* Make sure that '--query' and '--center' are not called together. */
- if(p->center && p->query)
- error(EXIT_FAILURE, 0, "the '--center' and '--query' options cannot be "
- "called together (they are parallel ways to define a query)");
+ if(p->query && (p->center || p->overlapwith) )
+ error(EXIT_FAILURE, 0, "the '--query' option cannot be called together "
+ "together with '--center' or '--overlapwith'");
+
+ /* Overlapwith cannot be called with the manual query. */
+ if( p->overlapwith && (p->center || p->width || p->radius) )
+ error(EXIT_FAILURE, 0, "the '--overlapwith' option cannot be called "
+ "with the manual region specifiers ('--center', '--width' or "
+ "'--radius')");
/* The radius and width cannot be called together. */
if(p->radius && p->width)
@@ -332,7 +337,7 @@ ui_read_check_only_options(struct queryparams *p)
"FITS file outputs are supported", p->cp.output);
}
else
- gal_checkset_automatic_output(&p->cp, "./query.fits", ".fits");
+ p->cp.output=gal_checkset_automatic_output(&p->cp, "./query.fits",
".fits");
}
diff --git a/bin/query/ui.h b/bin/query/ui.h
index 145bfb7..599f952 100644
--- a/bin/query/ui.h
+++ b/bin/query/ui.h
@@ -42,7 +42,7 @@ enum program_args_groups
/* Available letters for short options:
- a b e f i j k m n p t u v x y z
+ a b e f i j k m n p t u x y z
A B E G H J L R W X Y
*/
enum option_keys_enum
@@ -52,6 +52,7 @@ enum option_keys_enum
UI_KEY_QUERY = 'Q',
UI_KEY_DATASET = 's',
UI_KEY_CENTER = 'C',
+ UI_KEY_OVERLAPWITH = 'v',
UI_KEY_RADIUS = 'r',
UI_KEY_RANGE = 'g',
UI_KEY_COLUMN = 'c',
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 7e071ac..959c7a3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9994,6 +9994,10 @@ $ astquery --database=gaia --dataset=dr2
--output=my-gaia.fits \
--center=113.8729761,31.9027152 --radius=0.1 \
--column=source_id,ra,dec,phot_g_mean_mag
+## Find the ID, RA and Dec of all Gaia sources within an image.
+$ astquery --database=gaia --dataset=dr2 --overlapwith=image.fits
+ -csource_id,ra,dec
+
## Use a custom query to extract entries in the Gaia DR2 catalog.
## The 'XXXX YYYY' can be a query of any size on the command-line.
$ astquery --database=gaia --query="XXXX YYYY" --output=my-gaia.fits
@@ -10059,6 +10063,13 @@ Below is a list of the simplified names for the
databases that have them.
@end itemize
@end table
+@item -v STR
+@itemx --overlapwith=STR
+File name containing image (in the HDU given by @option{--hdu}) to use for
identifying the region to query in the give database and dataset.
+Based on the image's WCS, the area it covers is estimated and values to the
@option{--center}, @option{--width} will be calculated internally.
+Hence this option cannot be used with @code{--center}, @code{--width} or
@code{--radius}.
+Also, since it internally generates the query, it can't be used with
@code{--query}.
+
@item -C FLT,FLT
@itemx --center=FLT,FLT
The center to use for the automatically generated query (not compatible with
@option{--query}).
@@ -23217,6 +23228,11 @@ structure is not two dimensional and the units
(@code{CUNIT} keywords) are
not @code{deg} (for degrees), then this function will return a NaN.
@end deftypefun
+@deftypefun int gal_wcs_coverage (char @code{*filename}, char @code{*hdu},
size_t @code{*ondim}, double @code{**ocenter}, double @code{**owidth}, double
@code{**omin}, double @code{**omax})
+Find the sky coverage of the image HDU (@code{hdu}) within @file{filename}.
+The the number of dimensions is written into @code{ndim}, and space for the
various output arrays is internally allocated and filled with the respective
values.
+@end deftypefun
+
@deftypefun {gal_data_t *} gal_wcs_world_to_img (gal_data_t @code{*coords},
struct wcsprm @code{*wcs}, int @code{inplace})
Convert the linked list of world coordinates in @code{coords} to a linked
list of image coordinates given the input WCS structure. @code{coords} must
diff --git a/lib/fits.c b/lib/fits.c
index ef4cb76..36eaa28 100644
--- a/lib/fits.c
+++ b/lib/fits.c
@@ -821,7 +821,6 @@ gal_fits_hdu_open_format(char *filename, char *hdu, int
img0_tab1)
-
/**************************************************************/
/********** Header keywords ************/
/**************************************************************/
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index b4c0aac..7457411 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -154,7 +154,10 @@ gal_wcs_pixel_scale(struct wcsprm *wcs);
double
gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
-
+int
+gal_wcs_coverage(char *filename, char *hdu, size_t *ndim,
+ double **center, double **width, double **min,
+ double **max);
diff --git a/lib/wcs.c b/lib/wcs.c
index 492eb64..ef4e79d 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -39,6 +39,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/pointer.h>
#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
#include <gnuastro/permutation.h>
#include <gnuastro-internal/checkset.h>
@@ -1431,6 +1432,159 @@ gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
+int
+gal_wcs_coverage(char *filename, char *hdu, size_t *ondim,
+ double **ocenter, double **owidth, double **omin,
+ double **omax)
+{
+ fitsfile *fptr;
+ struct wcsprm *wcs;
+ int nwcs=0, type, status=0;
+ char *name=NULL, *unit=NULL;
+ gal_data_t *tmp, *coords=NULL;
+ size_t i, ndim, *dsize=NULL, numrows;
+ double *x, *y, *z, *min, *max, *center, *width;
+
+ /* Read the desired WCS. */
+ wcs=gal_wcs_read(filename, hdu, 0, 0, &nwcs);
+
+ /* If a WCS doesn't exist, return NULL. */
+ if(wcs==NULL) return 0;
+
+ /* Make sure the input HDU is an image. */
+ if( gal_fits_hdu_format(filename, hdu) != IMAGE_HDU )
+ error(EXIT_FAILURE, 0, "%s (hdu %s): is not an image HDU, the "
+ "'--skycoverage' option only applies to image extensions",
+ filename, hdu);
+
+ /* Get the array information of the image. */
+ fptr=gal_fits_hdu_open(filename, hdu, READONLY);
+ gal_fits_img_info(fptr, &type, ondim, &dsize, &name, &unit);
+ fits_close_file(fptr, &status);
+ ndim=*ondim;
+
+ /* Abort if we have more than 3 dimensions. */
+ if(ndim==1 || ndim>3) return 0;
+
+ /* Allocate the output datasets. */
+ center=*ocenter=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
+ "ocenter");
+ width=*owidth=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
+ "owidth");
+ min=*omin=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
+ "omin");
+ max=*omax=gal_pointer_allocate(GAL_TYPE_FLOAT64, ndim, 0, __func__,
+ "omax");
+
+ /* Now that we have the number of dimensions in the image, allocate the
+ space needed for the coordinates. */
+ numrows = (ndim==2) ? 5 : 9;
+ for(i=0;i<ndim;++i)
+ {
+ tmp=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &numrows, NULL, 0,
+ -1, 1, NULL, NULL, NULL);
+ tmp->next=coords;
+ coords=tmp;
+ }
+
+ /* Fill in the coordinate arrays, Note that 'dsize' is ordered in C
+ dimensions, for the WCS conversion, we need to have the dimensions
+ ordered in FITS/Fortran order. */
+ switch(ndim)
+ {
+ case 2:
+ x=coords->array; y=coords->next->array;
+ x[0] = 1; y[0] = 1;
+ x[1] = dsize[1]; y[1] = 1;
+ x[2] = 1; y[2] = dsize[0];
+ x[3] = dsize[1]; y[3] = dsize[0];
+ x[4] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
+ y[4] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
+ break;
+ case 3:
+ x=coords->array; y=coords->next->array; z=coords->next->next->array;
+ x[0] = 1; y[0] = 1; z[0]=1;
+ x[1] = dsize[2]; y[1] = 1; z[1]=1;
+ x[2] = 1; y[2] = dsize[1]; z[2]=1;
+ x[3] = dsize[2]; y[3] = dsize[1]; z[3]=1;
+ x[4] = 1; y[4] = 1; z[4]=dsize[0];
+ x[5] = dsize[2]; y[5] = 1; z[5]=dsize[0];
+ x[6] = 1; y[6] = dsize[1]; z[6]=dsize[0];
+ x[7] = dsize[2]; y[7] = dsize[1]; z[7]=dsize[0];
+ x[8] = dsize[2]/2 + (dsize[2]%2 ? 1 : 0.5f);
+ y[8] = dsize[1]/2 + (dsize[1]%2 ? 1 : 0.5f);
+ z[8] = dsize[0]/2 + (dsize[0]%2 ? 1 : 0.5f);
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to "
+ "fix the problem. 'ndim' of %zu is not recognized.",
+ __func__, PACKAGE_BUGREPORT, ndim);
+ }
+
+ /* For a check:
+ printf("IMAGE COORDINATES:\n");
+ for(i=0;i<numrows;++i)
+ if(ndim==2)
+ printf("%-15g%-15g\n", x[i], y[i]);
+ else
+ printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
+ */
+
+ /* Convert to the world coordinate system. */
+ gal_wcs_img_to_world(coords, wcs, 1);
+
+ /* For a check:
+ printf("\nWORLD COORDINATES:\n");
+ for(i=0;i<numrows;++i)
+ if(ndim==2)
+ printf("%-15g%-15g\n", x[i], y[i]);
+ else
+ printf("%-15g%-15g%-15g\n", x[i], y[i], z[i]);
+ */
+
+ /* Get the minimum and maximum values in each dimension. */
+ tmp=gal_statistics_minimum(coords);
+ min[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ tmp=gal_statistics_maximum(coords);
+ max[0] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ tmp=gal_statistics_minimum(coords->next);
+ min[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ tmp=gal_statistics_maximum(coords->next);
+ max[1] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ if(ndim>2)
+ {
+ tmp=gal_statistics_minimum(coords->next->next);
+ min[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ tmp=gal_statistics_maximum(coords->next->next);
+ max[2] = ((double *)(tmp->array))[0]; gal_data_free(tmp);
+ }
+
+ /* Write the center and width. */
+ switch(ndim)
+ {
+ case 2:
+ center[0]=x[4]; center[1]=y[4];
+ width[0]=max[0]-min[0]; width[1]=max[1]-min[1];
+ break;
+ case 3:
+ center[0]=x[8]; center[1]=y[8]; center[2]=z[8];
+ width[0]=max[0]-min[0]; width[1]=max[1]-min[1]; width[2]=max[2]-min[2];
+ break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to solve the "
+ "problem. The value %zu is not a recognized dimension", __func__,
+ PACKAGE_BUGREPORT, ndim);
+ }
+
+ /* Clean up and return success. */
+ free(dsize);
+ wcsfree(wcs);
+ return 1;
+}
+
+
+
+
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master f50cc13: Query: new --overlapwith option to automatically determine region,
Mohammad Akhlaghi <=