[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9952509: Correction in alignment and getting p
From: |
Mohammad Akhlaghi |
Subject: |
[gnuastro-commits] master 9952509: Correction in alignment and getting pixel scale |
Date: |
Tue, 17 Jan 2017 22:50:36 +0000 (UTC) |
branch: master
commit 995250968ae307529b4a00c9111f4ceef8b8c977
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Correction in alignment and getting pixel scale
The conditionals to get the WCS rotation matrix in ImageWarp's
`makealignmatrix' function were missing an `else' before the `altlin |= 2'
check! Thus when both flags where on, both conditions were executed which
would lead the to input matrix begin reset to zero. Therefore no alignment
would be done.
The reading of WCSLIB's `wcsprm' structure into a proper matrix is an issue
that is needed in other contexts too. Hence, the new function
`gal_wcs_array_from_wcsprm' was defined in `lib/wcs.c' to be a central
place for all Gnuastro to get this matrix when necessary and avoid such
bugs in the future.
This issue was submitted by Lee Kelvin along with a FITS image as an
example. After this correction, the image would be rotated back to its
original aligned state. However, I noticed that the pixel scales had
changed! This would not happen for the FITS images that are produced by
running `make check'!
After going over each step, I noticed that the problem was in the
`gal_wcs_pixel_scale_deg' function, the pixel scale it returned was wrong
(even for the aligned image). Until now, the procedure to get the pixel
scale was very crude: we would use WCSLIB to get the coordinates of the
points (0,0), (0,1), and (1,0), then use `gal_wcs_angular_distance_deg' to
find the angular distance between the points. Since a very large number of
floating point calculations are necessary in the process, I guess it was
because of floating point errors that we got the wrong pixel scale for
Lee's image, but correct values for the others.
So, after reviewing my Linear Algebra and also some searches on the
internet, I found the robust solution to this problem: using `Singular
Value Decomposition'. Fortunately, the GNU Scientific Library had
implementations of the decomposition. Thus, as a side-effect to the bug
reported by Lee, `gal_wcs_pixel_scale_deg' now also returns the correct
pixel scale through a robust process.
---
THANKS | 1 +
bin/imgwarp/imgwarp.c | 9 ++--
bin/imgwarp/ui.c | 40 +++++-----------
lib/gnuastro/wcs.h | 6 ++-
lib/wcs.c | 125 ++++++++++++++++++++++++++++++++++++++++++-------
5 files changed, 130 insertions(+), 51 deletions(-)
diff --git a/THANKS b/THANKS
index 36d3ad5..036ae65 100644
--- a/THANKS
+++ b/THANKS
@@ -21,6 +21,7 @@ support in Gnuastro. The list is ordered alphabetically.
Antonio Diaz Diaz address@hidden
Takashi Ichikawa address@hidden
Brandon Invergo address@hidden
+ Lee Kelvin address@hidden
Mohammad-Reza Khellat address@hidden
Alan Lefor address@hidden
Francesco Montanari address@hidden
diff --git a/bin/imgwarp/imgwarp.c b/bin/imgwarp/imgwarp.c
index 1207f58..27c9533 100644
--- a/bin/imgwarp/imgwarp.c
+++ b/bin/imgwarp/imgwarp.c
@@ -430,8 +430,9 @@ correctwcssaveoutput(struct imgwarpparams *p)
{
size_t i;
void *array;
+ double *pixelscale;
+ double *m=p->matrix, diff;
char keyword[9*FLEN_KEYWORD];
- double *m=p->matrix, diff, dx, dy;
struct gal_fits_key_ll *headers=NULL;
double tpc[4], tcrpix[3], *crpix=p->wcs->crpix, *pc=p->wcs->pc;
double tinv[4]={p->inverse[0]/p->inverse[8], p->inverse[1]/p->inverse[8],
@@ -484,9 +485,9 @@ correctwcssaveoutput(struct imgwarpparams *p)
signs are usually different.*/
if( p->wcs->pc[1]<ABSOLUTEFLTERROR ) p->wcs->pc[1]=0.0f;
if( p->wcs->pc[2]<ABSOLUTEFLTERROR ) p->wcs->pc[2]=0.0f;
- gal_wcs_pixel_scale_deg(p->wcs, &dx, &dy);
+ pixelscale=gal_wcs_pixel_scale_deg(p->wcs);
diff=fabs(p->wcs->pc[0])-fabs(p->wcs->pc[3]);
- if( fabs(diff/dx)<RELATIVEFLTERROR )
+ if( fabs(diff/pixelscale[0])<RELATIVEFLTERROR )
p->wcs->pc[3] = ( (p->wcs->pc[3] < 0.0f ? -1.0f : 1.0f)
* fabs(p->wcs->pc[0]) );
@@ -495,6 +496,8 @@ correctwcssaveoutput(struct imgwarpparams *p)
p->onaxes[1], p->onaxes[0], p->numnul, p->wcs,
headers, SPACK_STRING);
+ /* Clean up. */
+ free(pixelscale);
if(array!=p->output)
free(array);
}
diff --git a/bin/imgwarp/ui.c b/bin/imgwarp/ui.c
index f8ba672..e3fbebf 100644
--- a/bin/imgwarp/ui.c
+++ b/bin/imgwarp/ui.c
@@ -541,9 +541,7 @@ read_matrix(struct imgwarpparams *p)
void
makealignmatrix(struct imgwarpparams *p, double *tmatrix)
{
- double A, dx, dy;
- double amatrix[4];
- double w[4]={0,0,0,0};
+ double A, *w, *ps, amatrix[4];
/* Check if there is only two WCS axises: */
if(p->wcs->naxis!=2)
@@ -551,30 +549,12 @@ makealignmatrix(struct imgwarpparams *p, double *tmatrix)
"axises. For the `--align' option to operate it must be 2",
p->up.inputname, p->cp.hdu, p->wcs->naxis);
- /* Depending on the type of data, make the input matrix. Note that
- wcs->altlin is actually bit flags, not integers, so we have to compare
- with powers of two. */
- if(p->wcs->altlin |= 1)
- {
- w[0]=p->wcs->cdelt[0]*p->wcs->pc[0];
- w[1]=p->wcs->cdelt[0]*p->wcs->pc[1];
- w[2]=p->wcs->cdelt[1]*p->wcs->pc[2];
- w[3]=p->wcs->cdelt[1]*p->wcs->pc[3];
- }
- if(p->wcs->altlin |= 2)
- {
- w[0]=p->wcs->cd[0];
- w[1]=p->wcs->cd[1];
- w[2]=p->wcs->cd[2];
- w[3]=p->wcs->cd[3];
- }
- else
- error(EXIT_FAILURE, 0, "currently the `--align' option only recognizes "
- "PCi_ja and CDi_ja keywords, not any others");
/* Find the pixel scale along the two dimensions. Note that we will be
using the scale along the image X axis for both values. */
- gal_wcs_pixel_scale_deg(p->wcs, &dx, &dy);
+ w=gal_wcs_array_from_wcsprm(p->wcs);
+ ps=gal_wcs_pixel_scale_deg(p->wcs);
+
/* Lets call the given WCS orientation `W', the rotation matrix we want
to find as `X' and the final (aligned matrix) to have just one useful
@@ -626,28 +606,32 @@ makealignmatrix(struct imgwarpparams *p, double *tmatrix)
else
{
A = (w[3]/w[1]) - (w[2]/w[0]);
- amatrix[1] = dx / w[0] / A;
- amatrix[3] = dx / w[1] / A;
+ amatrix[1] = ps[0] / w[0] / A;
+ amatrix[3] = ps[0] / w[1] / A;
amatrix[0] = -1 * amatrix[1] * w[3] / w[1];
amatrix[2] = -1 * amatrix[3] * w[2] / w[0];
}
/* For a check:
- printf("dx: %e\n", dx);
+ printf("ps: %e\n", ps);
printf("w:\n");
printf(" %.8e %.8e\n", w[0], w[1]);
printf(" %.8e %.8e\n", w[2], w[3]);
printf("x:\n");
printf(" %.8e %.8e\n", amatrix[0], amatrix[1]);
printf(" %.8e %.8e\n", amatrix[2], amatrix[3]);
+ exit(0);
*/
-
/* Put the matrix elements into the output array: */
tmatrix[0]=amatrix[0]; tmatrix[1]=amatrix[1]; tmatrix[2]=0.0f;
tmatrix[3]=amatrix[2]; tmatrix[4]=amatrix[3]; tmatrix[5]=0.0f;
tmatrix[6]=0.0f; tmatrix[7]=0.0f; tmatrix[8]=1.0f;
+
+ /* Clean up. */
+ free(w);
+ free(ps);
}
diff --git a/lib/gnuastro/wcs.h b/lib/gnuastro/wcs.h
index aa81d29..4adda93 100644
--- a/lib/gnuastro/wcs.h
+++ b/lib/gnuastro/wcs.h
@@ -47,6 +47,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
__BEGIN_C_DECLS /* From C++ preparations */
+double *
+gal_wcs_array_from_wcsprm(struct wcsprm *wcs);
void
gal_wcs_xy_array_to_radec(struct wcsprm *wcs, double *xy, double *radec,
@@ -59,8 +61,8 @@ gal_wcs_radec_array_to_xy(struct wcsprm *wcs, double *radec,
double *xy,
double
gal_wcs_angular_distance_deg(double r1, double d1, double r2, double d2);
-void
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy);
+double *
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs);
double
gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs);
diff --git a/lib/wcs.c b/lib/wcs.c
index 58967bd..2d053d4 100644
--- a/lib/wcs.c
+++ b/lib/wcs.c
@@ -31,11 +31,66 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <unistd.h>
#include <assert.h>
+#include <gsl/gsl_linalg.h>
+
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+/**************************************************************/
+/********** Utilities ************/
+/**************************************************************/
+double *
+gal_wcs_array_from_wcsprm(struct wcsprm *wcs)
+{
+ double *out;
+ size_t i, j, size=wcs->naxis*wcs->naxis;
+
+ /* Allocate the necessary array. */
+ errno=0;
+ out=malloc(size*sizeof *out);
+ if(out==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for `out' in "
+ "`gal_wcs_array_from_wcsprm'", size*sizeof *out);
+
+ /* Fill in the array. */
+ if(wcs->altlin |= 1) /* Has a PCi_j array. */
+ {
+ for(i=0;i<wcs->naxis;++i)
+ for(j=0;j<wcs->naxis;++j)
+ out[i*wcs->naxis+j] = wcs->cdelt[i] * wcs->pc[i*wcs->naxis+j];
+ }
+ else if(wcs->altlin |= 2) /* Has CDi_j array */
+ {
+ for(i=0;i<size;++i)
+ out[i]=wcs->cd[i];
+ }
+ else
+ error(EXIT_FAILURE, 0, "currently, `gal_wcs_pixel_scale_deg' only "
+ "recognizes PCi_ja and CDi_ja keywords");
+
+ /* Return the result */
+ return out;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
/**************************************************************/
@@ -148,18 +203,42 @@ gal_wcs_angular_distance_deg(double r1, double d1, double
r2, double d2)
/* Return the pixel scale of the image in both dimentions in degrees. */
-void
-gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx, double *dy)
+double *
+gal_wcs_pixel_scale_deg(struct wcsprm *wcs)
{
- double radec[6], xy[]={0,0,1,0,0,1};
-
- /* Get the RA and Dec of the bottom left, bottom right and top left
- sides of the first pixel in the image. */
- gal_wcs_xy_array_to_radec(wcs, xy, radec, 3, 2);
-
- /* Calculate the distances and convert back to degrees: */
- *dx = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[2], radec[3]);
- *dy = gal_wcs_angular_distance_deg(radec[0], radec[1], radec[4], radec[5]);
+ gsl_vector S;
+ gsl_matrix A, V;
+ size_t n=wcs->naxis;
+ double *a, *v, *pixscale;
+
+ /* Allocate space for the `v' and `pixscale' arrays. */
+ errno=0; v=malloc(n*n*sizeof *v);
+ if(v==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for `v' in "
+ "`gal_wcs_pixel_scale_deg'", n*n*sizeof *v);
+ errno=0; pixscale=malloc(n*sizeof *pixscale);
+ if(pixscale==NULL)
+ error(EXIT_FAILURE, errno, "%zu bytes for `pixscale' in "
+ "`gal_wcs_pixel_scale_deg'", n*sizeof *pixscale);
+
+ /* Write the full matrix into an array, irrespective of what type it was
+ stored in the wcsprm structure (`PCi_j' style or `CDi_j' style). */
+ a=gal_wcs_array_from_wcsprm(wcs);
+
+ /* Fill in the necessary GSL vector and Matrix structures. */
+ S.size=n; S.stride=1; S.data=pixscale;
+ V.size1=n; V.size2=n; V.tda=n; V.data=v;
+ A.size1=n; A.size2=n; A.tda=n; A.data=a;
+
+ /* Run GSL's Singular Value Decomposition, using one-sided Jacobi
+ orthogonalization which computes the singular (scale) values to a
+ higher relative accuracy.*/
+ gsl_linalg_SV_decomp_jacobi(&A, &V, &S);
+
+ /* Clean up and return. */
+ free(a);
+ free(v);
+ return pixscale;
}
@@ -174,11 +253,21 @@ gal_wcs_pixel_scale_deg(struct wcsprm *wcs, double *dx,
double *dy)
double
gal_wcs_pixel_area_arcsec2(struct wcsprm *wcs)
{
- double dx, dy;
-
- /* Get the pixel scales along each axis in degrees. */
- gal_wcs_pixel_scale_deg(wcs, &dx, &dy);
-
- /* Return the result */
- return dx * dy * 3600.0f * 3600.0f;
+ double out;
+ double *pixscale;
+
+ /* A small sanity check. Later, when higher dimensions are necessary, we
+ can find which ones correlate to RA and Dec and use them to find the
+ pixel area in arcsec^2. */
+ if(wcs->naxis!=2)
+ error(EXIT_FAILURE, 0, "`gal_wcs_pixel_area_arcsec2' can currently "
+ "calculate the area only when the image has 2 dimensions.");
+
+ /* Get the pixel scales along each axis in degrees, then multiply. */
+ pixscale=gal_wcs_pixel_scale_deg(wcs);
+
+ /* Clean up and return the result. */
+ out = pixscale[0] * pixscale[1] * 3600.0f * 3600.0f;
+ free(pixscale);
+ return out;
}
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 9952509: Correction in alignment and getting pixel scale,
Mohammad Akhlaghi <=