[gnuastro-commits] master 5dbca4d: Warp checking for WCS when writing output

[Mohammad Akhlaghi]

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

    Warp checking for WCS when writing output
    
    Until now, when writing its output, warp was implicitly assuming that the
    dataset always has a WCS structure. Therefore when there was no WCS
    structure, Warp would given a segmentation fault. To fix this, the
    `correct_wcs_save_ouptut' function now checks if a WCS is present and only
    does WCS operations when it is.
    
    This bug was found thanks to a dataset provided by Bertrand Pain.
    
    This fixes bug #53580.
---
 NEWS                         |  1 +
 THANKS                       |  1 +
 bin/warp/ui.c                |  3 +-
 bin/warp/warp.c              | 69 +++++++++++++++++++++++---------------------
 doc/announce-acknowledge.txt |  1 +
 5 files changed, 41 insertions(+), 34 deletions(-)

diff --git a/NEWS b/NEWS
index e50e6e3..4dae408 100644
--- a/NEWS
+++ b/NEWS
@@ -142,6 +142,7 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
   bug #53312: Fits crash on keyword editing (except --delete).
   bug #53407: Instrumental noise in MakeNoise should be squared.
   bug #53424: Sigma-clipping seg-faults when there are no more elements.
+  bug #53580: Warp crash when no WCS present.
 
 
 
diff --git a/THANKS b/THANKS
index 219ff35..5103035 100644
--- a/THANKS
+++ b/THANKS
@@ -39,6 +39,7 @@ support in Gnuastro. The list is ordered alphabetically (by 
family name).
     Alan Lefor                           address@hidden
     Guillaume Mahler                     address@hidden
     Francesco Montanari                  address@hidden
+    Bertrand Pain                        address@hidden
     William Pence                        address@hidden
     Bob Proulx                           address@hidden
     Yahya Sefidbakht                     address@hidden
diff --git a/bin/warp/ui.c b/bin/warp/ui.c
index 085f828..e6d4f10 100644
--- a/bin/warp/ui.c
+++ b/bin/warp/ui.c
@@ -893,6 +893,7 @@ ui_preparations(struct warpparams *p)
 void
 ui_read_check_inputs_setup(int argc, char *argv[], struct warpparams *p)
 {
+  double *matrix;
   struct gal_options_common_params *cp=&p->cp;
 
 
@@ -940,7 +941,7 @@ ui_read_check_inputs_setup(int argc, char *argv[], struct 
warpparams *p)
   /* Everything is ready, notify the user of the program starting. */
   if(!p->cp.quiet)
     {
-      double *matrix=p->matrix->array;
+      matrix=p->matrix->array;
       printf(PROGRAM_NAME" started on %s", ctime(&p->rawtime));
       printf(" Using %zu CPU thread%s\n", p->cp.numthreads,
              p->cp.numthreads==1 ? "." : "s.");
diff --git a/bin/warp/warp.c b/bin/warp/warp.c
index 77ddd6c..7065ec6 100644
--- a/bin/warp/warp.c
+++ b/bin/warp/warp.c
@@ -414,7 +414,7 @@ correct_wcs_save_output(struct warpparams *p)
   double *m=p->matrix->array, diff;
   struct wcsprm *wcs=p->output->wcs;
   gal_fits_list_key_t *headers=NULL;
-  double *crpix=wcs->crpix, *w=p->inwcsmatrix;
+  double *crpix=wcs?wcs->crpix:NULL, *w=p->inwcsmatrix;
 
   /* `tinv' is the 2 by 2 inverse matrix. Recall that `p->inverse' is 3 by
      3 to account for homogeneous coordinates. */
@@ -422,28 +422,41 @@ correct_wcs_save_output(struct warpparams *p)
                   p->inverse[3]/p->inverse[8], p->inverse[4]/p->inverse[8]};
 
   /* Make the WCS corrections if necessary. */
-  if(p->keepwcs==0 && wcs)
+  if(wcs)
     {
-      /* 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;
+      if(p->keepwcs==0)
+        {
+          /* 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;
+        }
+
+      /* Due to floating point errors extremely small values of PC matrix
+         can 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( 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/p->pixelscale[0])<RELATIVEFLTERROR )
+        wcs->pc[3]=( (wcs->pc[3] < 0.0f ? -1.0f : 1.0f) * fabs(wcs->pc[0]) );
     }
 
   /* Add the appropriate headers: */
@@ -456,16 +469,6 @@ correct_wcs_save_output(struct warpparams *p)
                                 "Warp matrix element value", 0, NULL);
     }
 
-  /* Due to floating point errors extremely small values of PC matrix can
-     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( 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/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);
@@ -522,7 +525,7 @@ warp(struct warpparams *p)
   gal_threads_dist_in_threads(p->output->size, nt, &indexs, &thrdcols);
 
 
-  /* Start the convolution. */
+  /* Start the warp. */
   if(nt==1)
     {
       iwp[0].p=p;
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 28df3bd..c962921 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -5,6 +5,7 @@ Nima Dehdilani
 Antonio Diaz Diaz
 Lee Kelvin
 Guillaume Mahler
+Bertrand Pain
 Ole Streicher
 Michel Tallon
 Juan C. Tello