[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master fa3ef95: Pure rotation in Warp around pixel co

From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master fa3ef95: Pure rotation in Warp around pixel coordinate (0, 0)
Date: Sat, 1 Jul 2017 09:23:51 -0400 (EDT)

branch: master
commit fa3ef953118e6092ada99c045dc47db41193a8cc
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>

    Pure rotation in Warp around pixel coordinate (0,0)
    When doing rotation about pixel coordinate (0,0) by an angle theta, until
    now Warp would use the following wrong equation (which is correct for
    rotation of the coordinates, not rotation of the point):
      cos(theta)  sin(theta)
     -sin(theta)  cos(theta)
    When rotating the point (not the coordinates), the negative sign changes on
    the `sin' functions must change. However, this was not noticed! The reason
    was that in the standard aligned images, the first axis (RA) of the WCS
    points to the left, not the right. So the rotation would actually occur not
    around image coordinate (0,0), but image coordinate (X,0), where X is the
    length of the image along the first axis. This lead us to incorrectly
    associate rotation in the WCS axis directions (only when the RA and image X
    axis are in opposite directions).
    The problem became clear when trying rotation on an image where the RA and
    X axis were in the same direction. With this commit, the rotation matrix
    calculation is now corrected so like all warpings, the image coordinates
    are used, not WCS. The WCS information is just corrected.
    While working on this, a small issue in the `--align' option was noticed:
    when the first image axis and RA were in the same direction, `--align'
    wouldn't do any change! This has also been corrected with this commit.
    Also, as part of the fix to the `--align' option, it can now account for
    the possiblity that the pixels are not square (different pixel scales on
    different dimensions).
    This fixes bug #51353.
 bin/warp/main.h   |  2 ++
 bin/warp/ui.c     | 92 ++++++++++++++++++++++++++++---------------------------
 bin/warp/warp.c   | 69 +++++++++++++++++++++--------------------
 doc/gnuastro.texi | 76 +++++++++++++++++++++++++--------------------
 4 files changed, 127 insertions(+), 112 deletions(-)

diff --git a/bin/warp/main.h b/bin/warp/main.h
index 2d8c75e..6c544ea 100644
--- a/bin/warp/main.h
+++ b/bin/warp/main.h
@@ -55,6 +55,8 @@ struct warpparams
   gal_data_t      *matrix;  /* Warp/Transformation matrix.               */
   gal_data_t   *modularll;  /* List of modular warpings.                 */
   double         *inverse;  /* Inverse of the input matrix.              */
+  double     *inwcsmatrix;  /* Input WCS matrix.                         */
+  double      *pixelscale;  /* Pixel scale of input image.               */
   time_t          rawtime;  /* Starting time of the program.             */
   size_t       extinds[4];  /* Indexs of the minimum and maximum values. */
   size_t       ordinds[4];  /* Indexs of anticlockwise vertices.         */
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index bc36fdb..5000605 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -351,6 +351,11 @@ ui_check_options_and_arguments(struct warpparams *p)
       p->input->wcs=gal_wcs_read(p->inputname, p->cp.hdu, p->hstartwcs,
                                  p->hendwcs, &p->input->nwcs);
+      if(p->input->wcs)
+        {
+          p->pixelscale=gal_wcs_pixel_scale(p->input->wcs);
+          p->inwcsmatrix=gal_wcs_warp_matrix(p->input->wcs);
+        }
     error(EXIT_FAILURE, 0, "no input file is specified");
@@ -447,7 +452,7 @@ ui_matrix_prepare_raw(struct warpparams *p)
 static void
 ui_matrix_make_align(struct warpparams *p, double *tmatrix)
-  double A, *w, *ps, amatrix[4];
+  double A, *w, *P, x[4];
   /* Make sure the input image had a WCS structure. */
@@ -462,88 +467,83 @@ ui_matrix_make_align(struct warpparams *p, double 
           p->inputname, p->cp.hdu, p->input->wcs->naxis);
-  /* Find the pixel scale along the two dimensions. Note that we will be
-     using the scale along the image X axis for both values. */
-  w=gal_wcs_warp_matrix(p->input->wcs);
-  ps=gal_wcs_pixel_scale(p->input->wcs);
+  /* For help in reading, use aliases for the WCS matrix and pixel scale.*/
+  P=p->pixelscale;
+  w=p->inwcsmatrix;
   /* 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
-     value: `a' (which is the pixel scale):
+     to find as `X' and the final (aligned matrix) `P' (which is the pixel
+     scale):
-        x0  x1       w0  w1      -a  0
-        x2  x3   *   w2  w3   =   0  a
+        x0  x1       w0  w1      -P0   0
+        x2  x3   *   w2  w3   =   0    P1
      Let's open up the matrix multiplication, so we can find the `X'
      elements as function of the `W' elements and `a'.
-        x0*w0 + x1*w2 = -a                                         (1)
+        x0*w0 + x1*w2 = -P0                                        (1)
         x0*w1 + x1*w3 =  0                                         (2)
         x2*w0 + x3*w2 =  0                                         (3)
-        x2*w1 + x3*w3 =  a                                         (4)
+        x2*w1 + x3*w3 =  P1                                        (4)
      Let's bring the X with the smaller index in each equation to the left
-        x0 = (-w2/w0)*x1 - a/w0                                    (5)
+        x0 = (-w2/w0)*x1 - P0/w0                                   (5)
         x0 = (-w3/w1)*x1                                           (6)
         x2 = (-w2/w0)*x3                                           (7)
-        x2 = (-w3/w1)*x3 + a/w1                                    (8)
+        x2 = (-w3/w1)*x3 + P1/w1                                   (8)
     Using (5) and (6) we can find x0 and x1, by first eliminating x0:
-       (-w2/w0)*x1 - a/w0 = (-w3/w1)*x1 -> (w3/w1 - w2/w0) * x1 = a/w0
+       (-w2/w0)*x1 - P0/w0 = (-w3/w1)*x1  ->  (w3/w1 - w2/w0) * x1 = P0/w0
     For easy reading/writing, let's define: A = (w3/w1 - w2/w0)
-       --> x1 = a / w0 / A
+       --> x1 = P0 / w0 / A
        --> x0 = -1 * x1 * w3 / w1
     Similar to the above, we can find x2 and x3 from (7) and (8):
-       (-w2/w0)*x3 = (-w3/w1)*x3 + a/w1 -> (w3/w1 - w2/w0) * x3 = a/w1
+       (-w2/w0)*x3 = (-w3/w1)*x3 + P1/w1  ->  (w3/w1 - w2/w0) * x3 = P1/w1
-       --> x3 = a / w1 / A
+       --> x3 = P1 / w1 / A
        --> x2 = -1 * x3 * w2 / w0
-    Note that when the image is already aligned, a unity matrix should be
-    output.
-   */
+    Note that when the image is already aligned (off-diagonals are zero),
+    only the signs of the diagonal elements matter. */
   if( w[1]==0.0f && w[2]==0.0f )
-      amatrix[0]=1.0f;   amatrix[1]=0.0f;
-      amatrix[2]=0.0f;   amatrix[3]=1.0f;
+      x[0] = w[0]<0 ? 1.0f : -1.0f;  /* Has to be negative. */
+      x[1] = 0.0f;
+      x[2] = 0.0f;
+      x[3] = w[3]>0 ? 1.0f : 1.0f;   /* Has to be positive. */
       A = (w[3]/w[1]) - (w[2]/w[0]);
-      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];
+      x[1] = P[0] / w[0] / A;
+      x[3] = P[1] / w[1] / A;
+      x[0] = -1 * x[1] * w[3] / w[1];
+      x[2] = -1 * x[3] * w[2] / w[0];
   /* For a check:
-  printf("ps: %e\n", ps);
+  printf("ps: (%e, %e)\n", P[0], P[1]);
   printf("  %.8e    %.8e\n", w[0], w[1]);
   printf("  %.8e    %.8e\n", w[2], w[3]);
-  printf("  %.8e    %.8e\n", amatrix[0], amatrix[1]);
-  printf("  %.8e    %.8e\n", amatrix[2], amatrix[3]);
+  printf("  %.8e    %.8e\n", x[0], x[1]);
+  printf("  %.8e    %.8e\n", x[2], x[3]);
   /* 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);
+  tmatrix[0]=x[0];  tmatrix[1]=x[1]; tmatrix[2]=0.0f;
+  tmatrix[3]=x[2];  tmatrix[4]=x[3]; tmatrix[5]=0.0f;
+  tmatrix[6]=0.0f;  tmatrix[7]=0.0f; tmatrix[8]=1.0f;
@@ -584,7 +584,7 @@ ui_matrix_from_modular(struct warpparams *p)
   gal_data_t *pop;
   size_t dsize[]={3,3};
-  double s, c, v1, v2, *final, module[9];
+  double s, c, v1, v2, *final, module[9]={1,0,0,  0,1,0,  0,0,1};
   /* Reverse the list of modular warpings to be in the same order as the
      user specified.*/
@@ -622,9 +622,9 @@ ui_matrix_from_modular(struct warpparams *p)
         case UI_KEY_ROTATE:
           s = sin( v1 * M_PI / 180 );
           c = cos( v1 * M_PI / 180 );
-          module[0]=c;          module[1]=s;      module[2]=0.0f;
-          module[3]=-1.0f*s;    module[4]=c;      module[5]=0.0f;
-          module[6]=0.0f;       module[7]=0.0f;   module[8]=1.0f;
+          module[0]=c;          module[1]=-1.0f*s;  module[2]=0.0f;
+          module[3]=s;          module[4]=c;        module[5]=0.0f;
+          module[6]=0.0f;       module[7]=0.0f;     module[8]=1.0f;
         case UI_KEY_SCALE:
@@ -758,10 +758,10 @@ ui_matrix_finalize(struct warpparams *p)
     error(EXIT_FAILURE, 0, "the determinant of the given matrix "
           "is zero");
-  /* Check if the transformation is spatially invariant, in other words, if
-     it differs between differet regions of the output. If it doesn't we
-     can use this information for a more efficient processing. This is not
-     yet implemented. */
+  /* Note yet implemented: Check if the transformation is spatially
+     invariant, in other words, if it differs between differet regions of
+     the output. If it doesn't we can use this information for a more
+     efficient processing. */
    /* Make the inverse matrix: */
   inv=p->inverse=gal_data_malloc_array(GAL_TYPE_FLOAT64, 9, __func__,
@@ -993,6 +993,8 @@ ui_free_report(struct warpparams *p, struct timeval *t1)
+  if(p->pixelscale) free(p->pixelscale);
+  if(p->inwcsmatrix) free(p->inwcsmatrix);
   /* Report how long the operation took. */
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index 4584ec4..23fa4ee 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -111,8 +111,8 @@ along with Gnuastro. If not, see 
 /**************      Processing function      ******************/
-void *
-warponthread(void *inparam)
+static void *
+warp_onthread(void *inparam)
   struct iwpparams *iwp=(struct iwpparams*)inparam;
   struct warpparams *p=iwp->p;
@@ -291,8 +291,8 @@ warponthread(void *inparam)
    array to the input array. The order is fixed for all the pixels in
    the image altough the scale might change.
-warppreparations(struct warpparams *p)
+static void
+warp_preparations(struct warpparams *p)
   double is0=p->input->dsize[0], is1=p->input->dsize[1];
@@ -376,8 +376,6 @@ warppreparations(struct warpparams *p)
   p->opixarea=gal_polygon_area(forarea, 4);
   /* Find which index after transformation will have the minimum and
      maximum positions along the two axises. We can't use the starting
      loop because that is based on the input image which can be
@@ -405,38 +403,47 @@ warppreparations(struct warpparams *p)
 /* Correct the WCS coordinates (Multiply the 2x2 PC matrix of the WCS
-   structure by the INVERSE of the transform in 2x2 form that has been
-   converted from homogeneous coordinates). Then Multiply the crpix
-   array with the ACTUAL transformation matrix. */
+   structure by the INVERSE of the transform in 2x2). Then Multiply the
+   crpix array with the ACTUAL transformation matrix. */
 correct_wcs_save_output(struct warpparams *p)
   size_t i;
-  double *m=p->matrix->array, diff;
+  double tcrpix[3];
   char keyword[9*FLEN_KEYWORD];
+  double *m=p->matrix->array, diff;
   struct wcsprm *wcs=p->output->wcs;
   gal_fits_list_key_t *headers=NULL;
-  double tpc[4], tcrpix[3], *pixelscale;
-  double *crpix=wcs->crpix, *pc=wcs->pc;
+  double *crpix=wcs->crpix, *w=p->inwcsmatrix;
+  /* `tinv' is the 2 by 2 inverse matrix. Recall that `p->inverse' is 3 by
+     3 to account for homogeneous coordinates. */
   double tinv[4]={p->inverse[0]/p->inverse[8], p->inverse[1]/p->inverse[8],
                   p->inverse[3]/p->inverse[8], p->inverse[4]/p->inverse[8]};
+  /* Make the WCS corrections if necessary. */
   if(p->keepwcs==0 && wcs)
-      /* Correct the PC matrix: */
-      tpc[0]=pc[0]*tinv[0]+pc[1]*tinv[2];
-      tpc[1]=pc[0]*tinv[1]+pc[1]*tinv[3];
-      tpc[2]=pc[2]*tinv[0]+pc[3]*tinv[2];
-      tpc[3]=pc[2]*tinv[1]+pc[3]*tinv[3];
-      pc[0]=tpc[0]; pc[1]=tpc[1]; pc[2]=tpc[2]; pc[3]=tpc[3];
+      /* Correct the input WCS matrix. Since we are re-writing the PC
+         matrix from the full rotation matrix (including pixel scale),
+         we'll also have to set the CDELT fields to 1. Just to be sure that
+         the PC matrix is used in the end by WCSLIB, we'll also set altlin
+         to 1.*/
+      wcs->altlin=1;
+      wcs->cdelt[0] = wcs->cdelt[1] = 1.0f;
+      wcs->pc[0] = w[0]*tinv[0] + w[1]*tinv[2];
+      wcs->pc[1] = w[0]*tinv[1] + w[1]*tinv[3];
+      wcs->pc[2] = w[2]*tinv[0] + w[3]*tinv[2];
+      wcs->pc[3] = w[2]*tinv[1] + w[3]*tinv[3];
       /* Correct the CRPIX point. The +1 in the end of the last two
          lines is because FITS counts from 1. */
-      tcrpix[0]=m[0]*crpix[0]+m[1]*crpix[1]+m[2];
-      tcrpix[1]=m[3]*crpix[0]+m[4]*crpix[1]+m[5];
-      tcrpix[2]=m[6]*crpix[0]+m[7]*crpix[1]+m[8];
-      crpix[0]=tcrpix[0]/tcrpix[2]-p->outfpixval[0]+1;
-      crpix[1]=tcrpix[1]/tcrpix[2]-p->outfpixval[1]+1;
+      tcrpix[0] = m[0]*crpix[0]+m[1]*crpix[1]+m[2];
+      tcrpix[1] = m[3]*crpix[0]+m[4]*crpix[1]+m[5];
+      tcrpix[2] = m[6]*crpix[0]+m[7]*crpix[1]+m[8];
+      crpix[0] = tcrpix[0]/tcrpix[2] - p->outfpixval[0] + 1;
+      crpix[1] = tcrpix[1]/tcrpix[2] - p->outfpixval[1] + 1;
   /* Add the appropriate headers: */
@@ -453,20 +460,16 @@ correct_wcs_save_output(struct warpparams *p)
      be set to zero and extremely small differences between PC1_1 and PC2_2
      can be ignored. The reason for all the `fabs' functions is because the
      signs are usually different.*/
-  if( wcs->pc[1]<ABSOLUTEFLTERROR ) wcs->pc[1]=0.0f;
-  if( wcs->pc[2]<ABSOLUTEFLTERROR ) wcs->pc[2]=0.0f;
-  pixelscale=gal_wcs_pixel_scale(wcs);
+  if( fabs(wcs->pc[1])<ABSOLUTEFLTERROR ) wcs->pc[1]=0.0f;
+  if( fabs(wcs->pc[2])<ABSOLUTEFLTERROR ) wcs->pc[2]=0.0f;
-  if( fabs(diff/pixelscale[0])<RELATIVEFLTERROR )
+  if( fabs(diff/p->pixelscale[0])<RELATIVEFLTERROR )
     wcs->pc[3] =  ( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f) * fabs(wcs->pc[0]) );
   /* Save the output into the proper type and write it. */
     p->output=gal_data_copy_to_new_type_free(p->output, p->cp.type);
   gal_fits_img_write(p->output, p->cp.output, headers, PROGRAM_STRING);
-  /* Clean up. */
-  free(pixelscale);
@@ -512,7 +515,7 @@ warp(struct warpparams *p)
   /* Prepare the output array and all the necessary things: */
-  warppreparations(p);
+  warp_preparations(p);
   /* Distribute the output pixels into the threads: */
@@ -524,7 +527,7 @@ warp(struct warpparams *p)
-      warponthread(&iwp[0]);
+      warp_onthread(&iwp[0]);
@@ -543,7 +546,7 @@ warp(struct warpparams *p)
-            err=pthread_create(&t, &attr, warponthread, &iwp[i]);
+            err=pthread_create(&t, &attr, warp_onthread, &iwp[i]);
               error(EXIT_FAILURE, 0, "%s: can't create thread %zu",
                     __func__, i);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index b1ca151..fdd0cb5 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14188,57 +14188,65 @@ familiar with these concepts.
 @subsubsection Defining an ellipse
 @cindex Ellipse
address@hidden Axis ratio
address@hidden Position angle
 The PSF, see @ref{PSF}, and galaxy radial profiles are generally defined on
 an ellipse so in this section first defining an ellipse on a pixelated 2D
 surface is discussed. Labeling the major axis of an ellipse @mymath{a}, and
-its minor axis with @mymath{b}, the axis ratio is defined as:
+its minor axis with @mymath{b}, the @emph{axis ratio} is defined as:
 @mymath{q\equiv b/a}. The major axis of an ellipse can be aligned in any
 direction, therefore the angle of the major axis with respect to the
-horizontal axis of the image is defined to be the position angle of the
-ellipse and in this book, we show it with @mymath{\theta}.
+horizontal axis of the image is defined to be the @emph{position angle} of
+the ellipse and in this book, we show it with @mymath{\theta}.
 @cindex Radial profile on ellipse
-Our aim is to put a radial profile of any functional form
address@hidden(r)} over an ellipse. Let's define the radial distance
address@hidden as the distance on the major axis to the center of the
-ellipse which is located at @mymath{x_c} and @mymath{y_c}. We want to
-find the elliptical distance of a point located at @mymath{(i,j)}, in
-the image coordinate system, from the center of the ellipse. First the
-coordinate system is rotated by @mymath{\theta} to get the new rotated
-coordinates of that point @mymath{(i_r,j_r)}:
+Our aim is to put a radial profile of any functional form @mymath{f(r)}
+over an ellipse. Hence we need to associate a radius/distance to every
+point in space. Let's define the radial distance @mymath{r_{el}} as the
+distance on the major axis to the center of an ellipse which is located at
address@hidden and @mymath{i_c} (in other words @mymath{r_{el}\equiv{a}}). We
+want to find @mymath{r_{el}} of a point located at @mymath{(i,j)} (in the
+image coordinate system) from the center of the ellipse with axis ratio
address@hidden and position angle @mymath{\theta}. First the coordinate system
+is address@hidden not confuse this with the rotation matrix of
address@hidden basics}. In that equation, the point is rotated, here the
+coordinates are rotated and the point is fixed.} by @mymath{\theta} to get
+the new rotated coordinates of that point @mymath{(i_r,j_r)}:
 @cindex Elliptical distance
address@hidden The elliptical distance of a point located at @mymath{(i,j)}
-can now be defined as: @mymath{r_{el}^2=\sqrt{i_r^2+j_r^2/q^2} }. To
-place the radial profiles explained below over an ellipse,
address@hidden(r_{el}(i,j))} is calculated based on the functional radial
-profile desired.
address@hidden Recall that an ellipse is defined by 
+and that we defined @mymath{r_{el}\equiv{a}}. Hence, multiplying all
+elements of the the ellipse definition with @mymath{r_{el}^2} we get the
+elliptical distance at this point point located:
address@hidden/q^2} }. To place the radial profiles
+explained below over an ellipse, @mymath{f(r_{el})} is calculated based on
+the functional radial profile desired.
 @cindex Breadth first search
 @cindex Inside-out construction
 @cindex Making profiles pixel by pixel
 @cindex Pixel by pixel making of profiles
-The way MakeProfiles builds the profile is that the nearest pixel in the
-image to the given profile center is found and the profile value is
-calculated for it, see @ref{Sampling from a function}. The next pixel
-which the profile value is calculated on is the next nearest neighbor
-of the initial pixel to the profile center (as defined by
address@hidden). This is done fairly efficiently using a breadth
-first parsing
+MakeProfiles builds the profile starting from the nearest pixel in the
+image to the given profile center. The profile value is calculated for that
+central pixel using monte carlo integration, see @ref{Sampling from a
+function}. The next pixel is the next nearest neighbor to the central pixel
+as defined by @mymath{r_{el}}. This process goes on until the profile is
+fully built upto the trunctation radius. This is done fairly efficiently
+using a breadth first parsing
 which is implemented through an ordered linked list.
-Using this approach, we only go over one layer of pixels on the
-circumference of the profile to build the profile. Not one more extra
-pixel has to be checked. Another consequence of this strategy is that
-extending MakeProfiles to three dimensions becomes very simple: only
-the neighbors of each pixel have to be changed. Everything else after
-that (when the pixel index and its radial profile have entered the
-linked list) is the same, no matter the number of dimensions we are
-dealing with.
+Using this approach, we build the profile by expanding the
+circumference. Not one more extra pixel has to be checked (the calculation
+of @mymath{r_{el}} from above is not cheap in CPU terms). Another
+consequence of this strategy is that extending MakeProfiles to three
+dimensions becomes very simple: only the neighbors of each pixel have to be
+changed. Everything else after that (when the pixel index and its radial
+profile have entered the linked list) is the same, no matter the number of
+dimensions we are dealing with.

reply via email to

[Prev in Thread] Current Thread [Next in Thread]