[gnuastro-commits] master 09a9212b: Convolve: kernel is read before its sanity checks in PSF-matching

[Mohammad Akhlaghi]

branch: master
commit 09a9212b60783936ff814b601b557f6e1faab50c
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Convolve: kernel is read before its sanity checks in PSF-matching
    
    Until now, when requested to do PSF-matching, Convolve would crash with a
    segmentation fault. The issue can be reproduced with the commands below.
    
        astmkprof --kernel=moffat,5,7,5 --oversample=1 --output=wide.fits
        astmkprof --kernel=gaussian,3,5 --oversample=1 --output=sharp.fits
        astconvolve wide.fits --kernel=sharp.fits --makekernel=10 \
                    --output=matched-psf.fits
    
    This happened because the check for the dimensionality of the kernel was
    done before actually reading it!
    
    With this commit, the issue has been addressed by bringing up the step to
    read the kernel before the part where we check its dimensionality. Also the
    documentation of the '--makekernel' option was slightly edited to be more
    clear. To help catch PSF-matching crashes, a PSF-matchin test has been
    added in the 'make check' phase.
    
    This bug was found with the help of Raul Infante-Sainz.
    
    This fixes bug #62250.
---
 NEWS                        |  6 ++--
 bin/convolve/ui.c           |  6 ++--
 doc/gnuastro.texi           | 22 ++++++++----
 tests/Makefile.am           |  4 ++-
 tests/convolve/psf-match.sh | 88 +++++++++++++++++++++++++++++++++++++++++++++
 tests/mkprof/mkprofcat1.txt |  3 +-
 tests/mkprof/mosaic1.sh     |  3 +-
 7 files changed, 117 insertions(+), 15 deletions(-)

diff --git a/NEWS b/NEWS
index 2cc27cb4..13846c8d 100644
--- a/NEWS
+++ b/NEWS
@@ -14,8 +14,10 @@ See the end of the file for license conditions.
 ** Bugs fixed
   bug #62216: MakeProfiles crash when a 3D cube is requested and input
               catalog is from a pipe. Found by Irene Pintos Castro.
-  bug #62220: astscript-psf-stamp crashes when coordinate not in image,
-              found and fixed by Sepideh Eskandarlou.
+  bug #62220: astscript-psf-stamp crashes when coordinate not in image.
+              Found and fixed by Sepideh Eskandarlou.
+  bug #62250: Convolve crashes when doing PSF matching. Found with help of
+              Raul Infante-Sainz.
 
 
 
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index f26fe45a..f0eba4fa 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -561,15 +561,15 @@ ui_preparations(struct convolveparams *p)
      argument) is the blurry image. */
   if(p->makekernel)
     {
+      /* Read the kernel. */
+      ui_read_kernel(p);
+
       /* Currently this is not implemented in 1D. */
       if(p->kernel->ndim==1)
         error(EXIT_FAILURE, 0, "'--makekernel' is currently not available "
               "on 1D datasets");
       else
         {
-          /* Read in the kernel array. */
-          ui_read_kernel(p);
-
           /* Make sure the size of the kernel is the same as the input */
           if( p->input->dsize[0]!=p->kernel->dsize[0]
               || p->input->dsize[1]!=p->kernel->dsize[1] )
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a07947e6..eabf5574 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -16907,8 +16907,9 @@ $ astconvolve --kernel=psf.fits mockimg.fits
 ## Convolve in the spatial domain:
 $ astconvolve observedimg.fits --kernel=psf.fits --domain=spatial
 
-## Find the kernel to match sharper and blurry PSF images:
-$ astconvolve --kernel=sharperimage.fits --makekernel=10           \
+## Find the kernel to match sharper and blurry PSF images (they both
+## have to have the same pixel size).
+$ astconvolve --kernel=sharperimage.fits --makekernel=10 \
               blurryimage.fits
 
 ## Convolve a Spectrum (column 14 in the FITS table below) with a
@@ -16998,21 +16999,29 @@ For a full discussion, please see @ref{Edges in the 
spatial domain}.
 
 @item -m INT
 @itemx --makekernel=INT
-(@option{=INT}) If this option is called, Convolve will do de-convolution (see 
@ref{Convolution theorem}).
-The image specified by the @option{--kernel} option is assumed to be the 
sharper (less blurry) image and the input image is assumed to be the more 
blurry image.
+If this option is called, Convolve will do PSF-matching: the output will be 
the kernel that you should convolve with the sharper image to obtain the blurry 
one (see @ref{Convolution theorem}).
+The two images must have the same size (number of pixels).
+This option is not yet supported in 1-dimensional datasets.
+In effect, it is only necessary to give the two PSFs of your two datasets, 
find the matching kernel based on that, then apply that kernel to the 
higher-resolution (sharper image).
+
+The image given to the @option{--kernel} option is assumed to be the sharper 
(less blurry) image and the input image (with no option) is assumed to be the 
more blurry image.
 The value given to this option will be used as the maximum radius of the 
kernel.
 Any pixel in the final kernel that is larger than this distance from the 
center will be set to zero.
-The two images must have the same size.
 
 Noise has large frequencies which can make the result less reliable for the 
higher frequencies of the final result.
 So all the frequencies which have a spectrum smaller than the value given to 
the @option{minsharpspec} option in the sharper input image are set to zero and 
not divided.
 This will cause the wings of the final kernel to be flatter than they would 
ideally be which will make the convolved image result unreliable if it is too 
high.
-Some notes to take into account for a good result:
+
+Some notes to on how to prepare your two input PSFs.
+Note that these (and several other issues that relate to an accurate 
estimation of the PSF) are practically described in the following tutorial: 
@ref{Building the extended PSF}.
 
 @itemize
 @item
 Choose a bright (unsaturated) star and use a region box (with Crop for 
example, see @ref{Crop}) that is sufficiently above the noise.
 
+@item
+Mask all background sources that may be nearby (you can use Segment's clumps, 
see @ref{Segment}).
+
 @item
 Use Warp (see @ref{Warp}) to warp the pixel grid so the star's center is 
exactly on the center of the central pixel in the cropped image.
 This will certainly slightly degrade the result, however, it is necessary.
@@ -17022,7 +17031,6 @@ If there are multiple good stars, you can shift all of 
them, then normalize them
 The shifting might move the center of the star by one pixel in any direction, 
so crop the central pixel of the warped image to have a clean image for the 
de-convolution.
 @end itemize
 
-Note that this feature is not yet supported in 1-dimensional datasets.
 
 @item -c
 @itemx --minsharpspec
diff --git a/tests/Makefile.am b/tests/Makefile.am
index b4de5606..0c7be0fe 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -86,10 +86,11 @@ if COND_CONVERTT
 endif
 if COND_CONVOLVE
   MAYBE_CONVOLVE_TESTS = convolve/spatial.sh convolve/frequency.sh \
-                         convolve/spectrum-1d.sh
+                         convolve/psf-match.sh convolve/spectrum-1d.sh
 
   convolve/spectrum-1d.sh: prepconf.sh.log
   convolve/spatial.sh: mkprof/mosaic1.sh.log
+  convolve/psf-match.sh: mkprof/mosaic1.sh.log
   convolve/frequency.sh: mkprof/mosaic1.sh.log
 endif
 if COND_COSMICCAL
@@ -227,6 +228,7 @@ script/sort-by-night.sh: mkcatalog/aperturephot.sh.log
 # Environment variables for the test scripts.
 AM_TESTS_ENVIRONMENT = \
 export LANG=C; \
+export AWK=$(AWK); \
 export mkdir_p="$(MKDIR_P)"; \
 export progbdir=programs-built; \
 export topsrc=$(abs_top_srcdir); \
diff --git a/tests/convolve/psf-match.sh b/tests/convolve/psf-match.sh
new file mode 100755
index 00000000..d3cf3948
--- /dev/null
+++ b/tests/convolve/psf-match.sh
@@ -0,0 +1,88 @@
+# Find the matching kernel for two PSFs
+#
+# See the Tests subsection of the manual for a complete explanation
+# (in the Installing gnuastro section).
+#
+# Original author:
+#     Mohammad Akhlaghi <mohammad@akhlaghi.org>
+# Contributing author(s):
+# Copyright (C) 2022 Free Software Foundation, Inc.
+#
+# Copying and distribution of this file, with or without modification,
+# are permitted in any medium without royalty provided the copyright
+# notice and this notice are preserved.  This file is offered as-is,
+# without any warranty.
+
+
+
+
+
+# Preliminaries
+# =============
+#
+# Set the variables (The executable is in the build tree). Do the
+# basic checks to see if the executable is made or if the defaults
+# file exists (basicchecks.sh is in the source tree).
+prog=convolve
+psfwide=psf.fits
+psfsharp=psf-sharp.fits
+fitsprog=$progbdir/astfits
+cropprog=$progbdir/astcrop
+execname=$progbdir/ast$prog
+arithprog=$progbdir/astarithmetic
+
+
+
+
+
+# Skip?
+# =====
+#
+# If the dependencies of the test don't exist, then skip it. There are two
+# types of dependencies:
+#
+#   - The executable was not made (for example due to a configure option),
+#
+#   - The input data was not made (for example the test that created the
+#     data file failed).
+if [ ! -f $execname  ]; then echo "$execname not created."; exit 77; fi
+if [ ! -f $psfwide   ]; then echo "$psfwide does not exist."; exit 77; fi
+if [ ! -f $psfsharp  ]; then echo "$psfsharp does not exist."; exit 77; fi
+if [ ! -f $fitsprog  ]; then echo "$fitsprog does not exist."; exit 77; fi
+if [ ! -f $cropprog  ]; then echo "$cropprog does not exist."; exit 77; fi
+if [ ! -f $arithprog ]; then echo "$arithprog does not exist."; exit 77; fi
+
+
+
+
+
+# Prepare Sharper PSF
+# -------------------
+#
+# It is necessary that the two PSFs be the same width. So before finding
+# the matching PSF, it is necessary to make the sharper (smaller) PSF have
+# the same pixel size as the wider one.
+psfsharpwidth=psf-sharp-width-set.fits
+psfsharpnoblank=psf-sharp-no-blank.fits
+wwidth=$($fitsprog $psfwide --keyvalue=naxis1 --quiet)
+center=$($fitsprog $psfsharp --keyvalue=naxis1 --quiet \
+             | $AWK '{print int($1/2)+1}')
+$cropprog $psfsharp --center=$center,$center --width=$wwidth \
+          --mode=img --output=$psfsharpwidth
+$arithprog $psfsharpwidth set-i i i isblank 0 where \
+           --output=$psfsharpnoblank
+
+
+
+
+
+# Actual test script
+# ==================
+#
+# 'check_with_program' can be something like Valgrind or an empty
+# string. Such programs will execute the command if present and help in
+# debugging when the developer doesn't have access to the user's system.
+$check_with_program $execname $psfwide --kernel=$psfsharpnoblank \
+                              --makekernel=10   \
+                              --output=psf-match.fits \
+    && rm $psfsharpwidth $psfsharpnoblank
diff --git a/tests/mkprof/mkprofcat1.txt b/tests/mkprof/mkprofcat1.txt
index b94a0728..3f463ac1 100644
--- a/tests/mkprof/mkprofcat1.txt
+++ b/tests/mkprof/mkprofcat1.txt
@@ -1,7 +1,7 @@
 # Column 1:  ID                [count, u8]   Object ID
 # Column 2:  X                 [pixel, f64]  X axis position of profile center
 # Column 3:  Y                 [pixel, f64]  Y axis position of profile center
-# Column 4:  Function          [name, str7]  Profile's radial function
+# Column 4:  Function          [name, str9]  Profile's radial function
 # Column 5:  Width             [pixel, f64]  For Sersic: effective radius, for 
Moffat, FWHM
 # Column 6:  Sersic index      [none, f64]   Sersic index, or Moffat beta
 # Column 7:  Position angle    [deg, f64]    Position angle of profile
@@ -16,5 +16,6 @@
 # this notice are preserved.  This file is offered as-is, without any
 # warranty.
 1     0.0000     0.0000     moffat     3.000     4.765      0.0000     1.000   
   0.000     5.000
+1     0.0000     0.0000     gaussian   2.000     0.000      0.0000     1.000   
   0.000     5.000
 2     100.40     100.40     sersic     20.00     2.500      45.000     1.000   
   -18.0     5.000
 3     50.321     50.827     sersic     5.978     1.320      77.650     0.801   
   -15.0     5.000
diff --git a/tests/mkprof/mosaic1.sh b/tests/mkprof/mosaic1.sh
index 2269b3fd..c028ba5b 100755
--- a/tests/mkprof/mosaic1.sh
+++ b/tests/mkprof/mosaic1.sh
@@ -57,4 +57,5 @@ if [ ! -f $cat      ]; then echo "$cat does not exist.";   
exit 77; fi
 # debugging when the developer doesn't have access to the user's system.
 $check_with_program $execname $cat --mergedsize=100,100 \
                     --output=mkprofcat1.fits \
-    && mv 0_mkprofcat1.fits psf.fits
+    && mv 0_mkprofcat1.fits psf.fits \
+    && mv 1_mkprofcat1.fits psf-sharp.fits