[gnuastro-commits] master 2136295: Minimum sharp PSF spectrum is now an option not macro

[Mohammad Akhlaghi]

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

    Minimum sharp PSF spectrum is now an option not macro
    
    The minimum sharp PSF's spectrum value to use in Convolve's deconvolution
    was a macro until now, making it very hard to play with. So now Convolve
    accepts this value as an option. For the job, a new
    `gal_checkset_double_l_0_s_1' was defined in the internal Gnuastro library
    to directly read a double type value and check if it is between 0 and 1.
---
 bin/convolve/args.h           |   33 +++++++++++++++++++++++++--------
 bin/convolve/astconvolve.conf |    1 +
 bin/convolve/convolve.c       |   13 +++++++------
 bin/convolve/main.h           |    3 ++-
 bin/convolve/ui.c             |   21 ++++++++++++++++-----
 doc/gnuastro.texi             |   39 ++++++++++++++++++++++-----------------
 lib/checkset.c                |   26 +++++++++++++++++++++++++-
 lib/checkset.h                |    4 ++++
 8 files changed, 102 insertions(+), 38 deletions(-)

diff --git a/bin/convolve/args.h b/bin/convolve/args.h
index 6e8f335..53317b5 100644
--- a/bin/convolve/args.h
+++ b/bin/convolve/args.h
@@ -55,7 +55,7 @@ const char doc[] =
 
 /* Free letters for options:
 
-   c d e g i j l n r t u w x y z
+   d e g i j l n r t u w x y z
    A B C E F G H I J O Q R T W X Y Z
 
    Free numbers: >=504
@@ -115,6 +115,15 @@ static struct argp_option options[] =
       "Do not normalize the kernel image.",
       1
     },
+    {
+      "minsharpspec",
+      'c',
+      "FLT",
+      0,
+      "Deconvolution: min spectrum of sharp img.",
+      1
+    },
+
 
 
 
@@ -273,7 +282,8 @@ parse_opt(int key, char *arg, struct argp_state *state)
 
     /* Inputs: */
     case 'M':
-      gal_checkset_allocate_copy_set(arg, &p->up.maskname, &p->up.masknameset);
+      gal_checkset_allocate_copy_set(arg, &p->up.maskname,
+                                     &p->up.masknameset);
       break;
     case 'H':
       gal_checkset_allocate_copy_set(arg, &p->up.mhdu, &p->up.mhduset);
@@ -291,6 +301,11 @@ parse_opt(int key, char *arg, struct argp_state *state)
     case 501:
       p->kernelnorm=0;
       break;
+    case 'c':
+      gal_checkset_double_l_0_s_1(arg, &p->minsharpspec, "minsharpspec",
+                                 key, SPACK, NULL, 0);
+      p->up.minsharpspecset=1;
+      break;
 
 
     /* Output: */
@@ -298,21 +313,23 @@ parse_opt(int key, char *arg, struct argp_state *state)
 
    /* Mesh grid: */
     case 's':
-      gal_checkset_sizet_l_zero(arg, &p->mp.meshsize, "meshsize", key, SPACK,
-                                NULL, 0);
+      gal_checkset_sizet_l_zero(arg, &p->mp.meshsize, "meshsize", key,
+                                SPACK, NULL, 0);
       p->up.meshsizeset=1;
       break;
     case 'a':
-      gal_checkset_sizet_l_zero(arg, &p->mp.nch1, "nch1", key, SPACK, NULL, 0);
+      gal_checkset_sizet_l_zero(arg, &p->mp.nch1, "nch1", key, SPACK,
+                                NULL, 0);
       p->up.nch1set=1;
       break;
     case 'b':
-      gal_checkset_sizet_l_zero(arg, &p->mp.nch2, "nch2", key, SPACK, NULL, 0);
+      gal_checkset_sizet_l_zero(arg, &p->mp.nch2, "nch2", key, SPACK,
+                                NULL, 0);
       p->up.nch2set=1;
       break;
     case 'L':
-      gal_checkset_float_l_0_s_1(arg, &p->mp.lastmeshfrac, "lastmeshfrac", key,
-                                 SPACK, NULL, 0);
+      gal_checkset_float_l_0_s_1(arg, &p->mp.lastmeshfrac, "lastmeshfrac",
+                                 key, SPACK, NULL, 0);
       p->up.lastmeshfracset=1;
       break;
     case 503:
diff --git a/bin/convolve/astconvolve.conf b/bin/convolve/astconvolve.conf
index 79169c9..fb06348 100644
--- a/bin/convolve/astconvolve.conf
+++ b/bin/convolve/astconvolve.conf
@@ -21,6 +21,7 @@
  spatial            0
  frequency          1
  makekernel         0
+ minsharpspec       0.005
 
 # Input:
  hdu                0
diff --git a/bin/convolve/convolve.c b/bin/convolve/convolve.c
index e372239..45284ec 100644
--- a/bin/convolve/convolve.c
+++ b/bin/convolve/convolve.c
@@ -118,9 +118,8 @@ complexarraymultiply(double *a, double *b, size_t size)
 
 
 
-/* Divide the elements of the first array by the elements of the
-   second array and put the result into the elements of the first
-   array.
+/* Divide the elements of the first array by the elements of the second
+   array and put the result into the elements of the first array.
 
    (a+ib)/(c+id)=[(a+ib)*(c-id)]/[(c+id)*(c-id)]
                 =(ac-iad+ibc+bd)/(c^2+d^2)
@@ -130,14 +129,14 @@ complexarraymultiply(double *a, double *b, size_t size)
    on the loop.
  */
 void
-complexarraydivide(double *a, double *b, size_t size)
+complexarraydivide(double *a, double *b, size_t size, double minsharpspec)
 {
   double r, *af;
 
   af=a+2*size;
   do
     {
-      if (sqrt(*b**b + *(b+1)**(b+1))>MINGOODDIVSPEC)
+      if (sqrt(*b**b + *(b+1)**(b+1))>minsharpspec)
         {
           r      = ( ( (*a * *b) + (*(a+1) * *(b+1)) )
                      / ( *b * *b + *(b+1) * *(b+1) ) );
@@ -145,6 +144,8 @@ complexarraydivide(double *a, double *b, size_t size)
                      / ( *b * *b + *(b+1) * *(b+1) ) );
           *a=r;
 
+          /* Just as a sanity check (the result should never be larger than
+             one. */
           if(sqrt(*a**a + *(a+1)**(a+1))>1.00001f)
             *a=*(a+1)=0.0f;
         }
@@ -679,7 +680,7 @@ frequencyconvolve(struct convolveparams *p)
   if(verb) gettimeofday(&t1, NULL);
   if(p->makekernel)
     {
-      complexarraydivide(p->pimg, p->pker, p->ps0*p->ps1);
+      complexarraydivide(p->pimg, p->pker, p->ps0*p->ps1, p->minsharpspec);
       if(verb)
         gal_timing_report(&t1, "Divided in the frequency domain.", 1);
     }
diff --git a/bin/convolve/main.h b/bin/convolve/main.h
index 7cb25d0..d5158e5 100644
--- a/bin/convolve/main.h
+++ b/bin/convolve/main.h
@@ -39,7 +39,6 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 
 
 
-#define MINGOODDIVSPEC       0.005
 #define CONVFLOATINGPOINTERR 1e-10
 #define COMPLEXTOREALSPEC    1  /* Spectrum of complex number.  */
 #define COMPLEXTOREALPHASE   2  /* Phase of the complex number. */
@@ -63,6 +62,7 @@ struct uiparams
   int            mhduset;
   int      kernelnameset;
   int            khduset;
+  int    minsharpspecset;
   int        meshsizeset;
   int            nch1set;
   int            nch2set;
@@ -90,6 +90,7 @@ struct convolveparams
   size_t             is1;   /* Input image size along C's second axis.  */
   size_t             ks0;   /* Kernel size along C's first axis.        */
   size_t             ks1;   /* Kernel size along C's second axis.       */
+  double    minsharpspec;   /* Deconvolution: min spectrum of sharp img.*/
   int         kernelflip;   /* ==1: Flip the kernel.                    */
   int         kernelnorm;   /* ==1: Normalize the kernel.               */
   int               nwcs;   /* Number of WCS headers.                   */
diff --git a/bin/convolve/ui.c b/bin/convolve/ui.c
index 6b57113..0d65275 100644
--- a/bin/convolve/ui.c
+++ b/bin/convolve/ui.c
@@ -110,7 +110,13 @@ readconfig(char *filename, struct convolveparams *p)
                                        &up->kernelnameset);
       else if (strcmp(name, "khdu")==0)
         gal_checkset_allocate_copy_set(value, &up->khdu, &up->khduset);
-
+      else if(strcmp(name, "minsharpspec")==0)
+        {
+          if(up->minsharpspecset) continue;
+          gal_checkset_double_l_0_s_1(value, &p->minsharpspec, "minsharpspec",
+                                     key, SPACK, filename, lineno);
+          up->minsharpspecset=1;
+        }
 
 
       /* Outputs: */
@@ -243,6 +249,8 @@ printvalues(FILE *fp, struct convolveparams *p)
     GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("kernel", up->kernelname);
   if(up->khduset)
     GAL_CHECKSET_PRINT_STRING_MAYBE_WITH_SPACE("khdu", up->khdu);
+  if(up->minsharpspecset)
+    fprintf(fp, CONF_SHOWFMT"%f\n", "minsharpspec", p->minsharpspec);
 
 
 
@@ -298,6 +306,8 @@ checkifset(struct convolveparams *p)
     GAL_CONFIGFILES_REPORT_NOTSET("kernel");
   if(up->khduset==0)
     GAL_CONFIGFILES_REPORT_NOTSET("khdu");
+  if(p->makekernel && up->minsharpspecset==0)
+    GAL_CONFIGFILES_REPORT_NOTSET("minsharpspec");
 
 
   /* Mesh grid: */
@@ -447,14 +457,15 @@ preparearrays(struct convolveparams *p)
             "pixels' section of the Gnuastro manual for more information.");
 
 
-  /* Read the file specified by --kernel. If makekernel is specified,
-     then this is actually the low resolution image. */
+  /* Read the file specified by --kernel. If makekernel is specified, then
+     this is actually the sharper image and the input image (given as an
+     argument) is the blurry image. */
   if(p->makekernel)
     {
       /* Read in the kernel array: */
       gal_fits_file_to_float(up->kernelname, NULL, up->khdu, NULL,
-                                  &p->kernel, &bitpix, &anyblank,
-                                  &p->ks0, &p->ks1);
+                             &p->kernel, &bitpix, &anyblank,
+                             &p->ks0, &p->ks1);
       if(p->ks0!=p->is0 || p->ks1!=p->is1)
         error(EXIT_FAILURE, 0, "with the `--makekernel' (`-m') option, "
               "the input image and the image specified with the kernel "
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 28b6d2b..be99ad3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -8659,7 +8659,7 @@ One line examples:
 
 @example
 $ astconvolve --kernel=psf.fits mockimg.fits
-$ astconvolve --kernel=sharperimage.fits --makekernel\
+$ astconvolve --kernel=sharperimage.fits --makekernel=10 \
               blurryimage.fits
 @end example
 
@@ -8764,21 +8764,21 @@ of the image), we also remove any pixel with a value 
less than
 
 @item -m
 @itemx --makekernel
-(@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. The two images have to be the same size. Some notes to take
-into account for a good result:
-
-Noise has large frequencies which can make the result less reliable
-for the higher frequencies of the kernel. So all the frequencies which
-have a spectrum smaller than 0.005 in the frequency domain 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. 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.
+(@option{=INT}) If this option is called, Convolve will do de-convolution
+(see @ref{Convolution theorem}). The image specified by the
address@hidden option is assumed to be the sharper (less blurry) image
+and the input image 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:
 @itemize
 
 @item
@@ -8798,8 +8798,13 @@ pixels is one) and then take their average to decrease 
this effect.
 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
+
address@hidden -c
address@hidden --minsharpspec
+(@option{=FLT}) The minimum frequency spectrum (or coefficient, or pixel
+value in the frequency domain image) to use in deconvolution, see the
+explanations under the @option{--makekernel} option for more information.
 @end table
 
 
diff --git a/lib/checkset.c b/lib/checkset.c
index cacacee..7141f45 100644
--- a/lib/checkset.c
+++ b/lib/checkset.c
@@ -402,7 +402,7 @@ gal_checkset_any_float(char *optarg, float *var, char *lo, 
char so,
 
 void
 gal_checkset_double_l_0(char *optarg, double *var, char *lo, char so,
-                        char* spack, char *filename, size_t lineno)
+                        char *spack, char *filename, size_t lineno)
 {
   double tmp;
   char *tailptr;
@@ -425,6 +425,30 @@ gal_checkset_double_l_0(char *optarg, double *var, char 
*lo, char so,
 
 
 void
+gal_checkset_double_l_0_s_1(char *optarg, double *var, char *lo, char so,
+                            char *spack, char *filename, size_t lineno)
+{
+  double tmp;
+  char *tailptr;
+  *var=tmp=strtod(optarg, &tailptr);
+
+  CHECKFULLNUMBER;
+  if(tmp>1.0f || tmp<0)
+    {
+      if(filename)
+        error_at_line(EXIT_FAILURE, 0, filename, lineno,
+                      FIXEDFORFILE" "NOTEMSG_SMALLERONE, lo, optarg);
+      else
+        error(EXIT_FAILURE, 0, FIXEDFOROPTION" "NOTEMSG_SMALLERONE,
+              lo, so, optarg);
+    }
+}
+
+
+
+
+
+void
 gal_checkset_double_el_0(char *optarg, double *var, char *lo, char so,
                          char* spack, char *filename, size_t lineno)
 {
diff --git a/lib/checkset.h b/lib/checkset.h
index 7954a7e..16471ec 100644
--- a/lib/checkset.h
+++ b/lib/checkset.h
@@ -163,6 +163,10 @@ gal_checkset_double_l_0(char *optarg, double *var, char 
*lo, char so,
                         char *spack, char *filename, size_t lineno);
 
 void
+gal_checkset_double_l_0_s_1(char *optarg, double *var, char *lo, char so,
+                            char *spack, char *filename, size_t lineno);
+
+void
 gal_checkset_double_el_0(char *optarg, double *var, char *lo, char so,
                          char* spack, char *filename, size_t lineno);