[gnuastro-commits] master 17d6f65 004/113: Merged recent corrections in master

[Mohammad Akhlaghi]

branch: master
commit 17d6f651288bf90aebc3834daa69a12df4080c1d
Merge: a11e71c fa3ef95
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    Merged recent corrections in master
---
 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->cp.minmapsize);
       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);
+        }
     }
   else
     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. */
   if(p->input->wcs==NULL)
@@ -462,88 +467,83 @@ ui_matrix_make_align(struct warpparams *p, double 
*tmatrix)
           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
      side:
 
-        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. */
     }
   else
     {
       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("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]);
+  printf("  %.8e    %.8e\n", x[0], x[1]);
+  printf("  %.8e    %.8e\n", x[2], x[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);
+  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;
           break;
 
         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)
   free(p->cp.output);
   gal_data_free(p->input);
   gal_data_free(p->matrix);
+  if(p->pixelscale) free(p->pixelscale);
+  if(p->inwcsmatrix) free(p->inwcsmatrix);
 
   /* Report how long the operation took. */
   if(!p->cp.quiet)
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 
<http://www.gnu.org/licenses/>.
 /***************************************************************/
 /**************      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.
 */
-void
-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. */
 void
 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;
   diff=fabs(wcs->pc[0])-fabs(wcs->pc[3]);
-  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. */
   if(p->cp.type!=p->output->type)
     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)
     {
       iwp[0].p=p;
       iwp[0].indexs=indexs;
-      warponthread(&iwp[0]);
+      warp_onthread(&iwp[0]);
     }
   else
     {
@@ -543,7 +546,7 @@ warp(struct warpparams *p)
             iwp[i].p=p;
             iwp[i].b=&b;
             iwp[i].indexs=&indexs[i*thrdcols];
-            err=pthread_create(&t, &attr, warponthread, &iwp[i]);
+            err=pthread_create(&t, &attr, warp_onthread, &iwp[i]);
             if(err)
               error(EXIT_FAILURE, 0, "%s: can't create thread %zu",
                     __func__, i);
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c36dac0..faad1b4 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14186,57 +14186,65 @@ familiar with these concepts.
 @subsubsection Defining an ellipse
 
 @cindex Ellipse
+@cindex Axis ratio
+@cindex 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
-@mymath{f(r)} over an ellipse. Let's define the radial distance
-@mymath{r_{el}} 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)}:
-
-@dispmath{i_r(i,j)=(i_c-i)\cos(\theta)+(j_c-j)\sin(\theta)}
-@dispmath{j_r(i,j)=(j_c-j)\cos(\theta)-(i_c-i)\sin(\theta)}
+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
+@mymath{i_c} 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
+@mymath{q} and position angle @mymath{\theta}. First the coordinate system
+is rotated@footnote{Do not confuse this with the rotation matrix of
+@ref{Warping 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)}:
+
+@dispmath{i_r(i,j)=+(i_c-i)\cos(\theta)+(j_c-j)\sin(\theta)}
+@dispmath{j_r(i,j)=-(i_c-i)\sin(\theta)+(j_c-j)\cos(\theta)}
 
 @cindex Elliptical distance
-@noindent 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,
-@mymath{f(r_{el}(i,j))} is calculated based on the functional radial
-profile desired.
+@noindent Recall that an ellipse is defined by @mymath{(i_r/a)^2+(j_r/b)^2=1}
+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:
+@mymath{r_{el}=\sqrt{i_r^2+j_r^2/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
-@mymath{r_{el}}). 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
 strategy@footnote{@url{http://en.wikipedia.org/wiki/Breadth-first_search}}
 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.