[gnuastro-commits] master b4c2b5a: gal_box_bound_ellipse return width isn't extra

[Mohammad Akhlaghi]

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

    gal_box_bound_ellipse return width isn't extra
    
    Until now, the width returned by `gal_box_bound_ellipse' would be 2 pixels
    extra along each dimension, thus leaving some extra pixels on the edges
    when MakeProfiles used it. This was done to accommodate extreme cases,
    where the center of a profile was just on the edge of the pixel and the
    trunctation radius also had a fractional value.
    
    However, after some checks, we see that even in the most extreme cases,
    only a single (very distant) elliptical layer may be lost (not at all
    significant). But in most (normal cases), these extra layers will cause
    many extra pixels which can be a burden on MakeProfiles its self or
    higher-level programs (for example when building a kernel to convolve
    with).
    
    Also, while checking the `gal_box_bound_ellipse' functions, I noticed that
    in the documentation, we were mistakenly describing `a' and `b' as
    major/minor axis length. While in practice they are `semi-major' and
    `semi-minor'.
    
    Since there is no extra layer, we also don't need the call to Crop when
    making NoiseChisel's default internal kernel (in
    `lib/gnuastro-internal/kernel-2d.h'). A new kernel is now created and used
    (with a larger number of random points).
---
 bin/mkprof/oneprofile.c           | 23 +++++++++++++++++-
 doc/gnuastro.texi                 | 27 +++++++++++----------
 lib/box.c                         |  4 ++--
 lib/gnuastro-internal/kernel-2d.h | 50 +++++++++++++++++++--------------------
 4 files changed, 63 insertions(+), 41 deletions(-)

diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 865e63f..aab2524 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -172,7 +172,9 @@ oneprofile_r_circle(size_t index, struct mkonthread *mkp)
 float
 oneprofile_randompoints(struct mkonthread *mkp)
 {
+  double r_before=mkp->r;
   double range[2], sum=0.0f;
+  double coord_before[2]={mkp->coord[0], mkp->coord[1]};
   size_t i, j, numrandom=mkp->p->numrandom, ndim=mkp->p->ndim;
 
   /* Set the range in each dimension. */
@@ -188,6 +190,15 @@ oneprofile_randompoints(struct mkonthread *mkp)
       sum+=mkp->profile(mkp);
     }
 
+  /* Reset the original distance and coordinate of the pixel and return the
+     average random value. The resetting is mostly redundant (only useful
+     in checks), but since it has a very negligible cost (compared to the
+     random checks above) cost, its good to reset it to help in debugging
+     when necessary (avoid confusion when un-commenting the checks in
+     `oneprofile_pix_by_pix'). */
+  mkp->r=r_before;
+  mkp->coord[0]=coord_before[0];
+  mkp->coord[1]=coord_before[1];
   return sum/numrandom;
 }
 
@@ -376,7 +387,7 @@ oneprofile_pix_by_pix(struct mkonthread *mkp)
           oneprofile_r_el(mkp);
           if(mkp->r > truncr) continue;
 
-          /* Find the value for this pixel: */
+          /* Set the range for this pixel. */
           for(i=0;i<ndim;++i)
             {
               mkp->lower[i]  = mkp->coord[i] - hp;
@@ -389,6 +400,12 @@ oneprofile_pix_by_pix(struct mkonthread *mkp)
           if (fabs(array[p]-approx)/array[p] < tolerance)
             use_rand_points=0;
 
+          /* For a check:
+          printf("coord: %g, %g\n", mkp->coord[0], mkp->coord[1]);
+          printf("r_rand: %g (rand: %g, center: %g)\n\n", mkp->r, array[p],
+                 approx);
+          */
+
           /* Save the peak flux if this is the first pixel: */
           if(ispeak) { mkp->peakflux=array[p]; ispeak=0; }
 
@@ -440,6 +457,10 @@ oneprofile_pix_by_pix(struct mkonthread *mkp)
       /* Find the value for this pixel: */
       array[p]=profile(mkp);
 
+      /* For a check:
+      printf("r_center: %g\n", mkp->r);
+      */
+
       /* Save the peak flux if this is the first pixel: */
       if(ispeak) { mkp->peakflux=array[p]; ispeak=0; }
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 6b8759c..c850428 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -24460,22 +24460,23 @@ vertical).
 @deftypefun void gal_box_bound_ellipse_extent (double @code{a}, double 
@code{b}, double @code{theta_deg}, double @code{*extent})
 Return the maximum extent along each dimension of the given ellipse from
 the center of the ellipse. Therefore this is half the extent of the box in
-each dimension. @code{a} is the ellipse major axis, @code{b} is the minor
-axis, @code{theta_deg} is the position angle in degrees. The extent in each
-dimension is in floating point format and stored in @code{extent} which
-must already be allocated before this function.
+each dimension. @code{a} is the ellipse semi-major axis, @code{b} is the
+semi-minor axis, @code{theta_deg} is the position angle in degrees. The
+extent in each dimension is in floating point format and stored in
address@hidden which must already be allocated before this function.
 @end deftypefun
 
 @deftypefun void gal_box_bound_ellipse (double @code{a}, double @code{b}, 
double @code{theta_deg}, long @code{*width})
-Any ellipse can be enclosed into a rectangular box. The purpose of this
-function is to give the height and width of that box assuming the center of
-the ellipse is in the box center. @code{a} is the ellipse major axis,
address@hidden is the minor axis, @code{theta_deg} is the position angle in
-degrees. The @code{width} array will contain the output size in long
-integer type. @code{width[0]}, and @code{width[1]} are the number of pixels
-along the first and second FITS axis. Since the ellipse center is assumed
-to be in the center of the box, all the values in @code{width} will be an
-odd integer.
+Any ellipse can be enclosed into a rectangular box. This function will
+write the height and width of that box where @code{width} points to. It
+assumes the center of the ellipse is located within the central pixel of
+the box. @code{a} is the ellipse semi-major axis length, @code{b} is the
+semi-minor axis, @code{theta_deg} is the position angle in degrees. The
address@hidden array will contain the output size in long integer
+type. @code{width[0]}, and @code{width[1]} are the number of pixels along
+the first and second FITS axis. Since the ellipse center is assumed to be
+in the center of the box, all the values in @code{width} will be an odd
+integer.
 @end deftypefun
 
 @deftypefun void gal_box_border_from_center (double @code{center}, size_t 
@code{ndim}, long @code{*width}, long @code{*fpixel}, long @code{*lpixel})
diff --git a/lib/box.c b/lib/box.c
index 9d80d37..55f77c6 100644
--- a/lib/box.c
+++ b/lib/box.c
@@ -92,8 +92,8 @@ gal_box_bound_ellipse(double a, double b, double theta_deg, 
long *width)
      want the final height and width of the box enclosing the
      ellipse. So we have to multiply them by two, then take one from
      them (for the center). */
-  width[0] = 2 * ( (long)extent[0] + 1 ) + 1;
-  width[1] = 2 * ( (long)extent[1] + 1 ) + 1;
+  width[0] = 2 * ( (long)extent[0] ) + 1;
+  width[1] = 2 * ( (long)extent[1] ) + 1;
 }
 
 
diff --git a/lib/gnuastro-internal/kernel-2d.h 
b/lib/gnuastro-internal/kernel-2d.h
index ccc2d49..3462467 100644
--- a/lib/gnuastro-internal/kernel-2d.h
+++ b/lib/gnuastro-internal/kernel-2d.h
@@ -34,22 +34,21 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    Make the kernel
    ---------------
 
-   We'll first make the kernel using MakeProfiles, then crop the outer
-   layers that are zero. Put these lines into a script and run it.
+   We'll first make the kernel using MakeProfiles with the following
+   commands. IMPORTANT NOTE: because the kernel is so sharp, random
+   sampling is going to be used for all the pixels (the center won't be
+   used). So it is important to have a large number of random points to
+   make the very slight differences between symmetric parts of the profile
+   even less significant.
 
-     set -o errexit           # Stop if a program returns false.
-     export GSL_RNG_TYPE=ranlxs2
      export GSL_RNG_SEED=1
-     astmkprof --kernel=gaussian,2,5 --oversample=1 --envseed -ok2.fits
-     astcrop k2.fits --section=2:*-1,2:*-1 --zeroisnotblank    \
-             --mode=img --output=fwhm2.fits
-     rm k2.fits
+     export GSL_RNG_TYPE=ranlxs2
+     astmkprof --kernel=gaussian,2,5 --oversample=1 --envseed 
--numrandom=100000
 
    Convert it to C code
    --------------------
 
-   We'll use this tiny C program to convert the FITS file into the two
-   important C variables:
+   Put the following C program into a file called `kernel.c'.
 
      #include <stdio.h>
      #include <stdlib.h>
@@ -60,7 +59,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
      {
        size_t i;
        float *arr;
-       gal_data_t *img=gal_fits_img_read_to_type("fwhm2.fits", "1",
+       gal_data_t *img=gal_fits_img_read_to_type("kernel.fits", "1",
                                                  GAL_TYPE_FLOAT32, -1);
 
        arr=img->array;
@@ -97,33 +96,34 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
    -----------------
 
    We can now compile and run that C program and put the outputs in
-   `ready.c'. Once its created, copy the contents of `ready.c' after these
-   comments.
+   `kernel.c'. Once its created, copy the contents of `kernel-2d.h' after
+   these comments.
 
-     $ astbuildprog -q prepare.c > ready.c
+     $ astbuildprog -q kernel.c > kernel-2d.h
  */
 
 size_t kernel_2d_dsize[2]={11, 11};
-float kernel_2d[121]={0, 0, 0, 0, 0, 2.570758e-08, 0, 0, 0, 0, 0,
+float kernel_2d[121]={0, 0, 0, 0, 0, 2.599797e-08, 0, 0, 0, 0, 0,
+
+0, 0, 3.008479e-08, 6.938075e-07, 4.493532e-06, 8.276223e-06, 4.515019e-06, 
6.947793e-07, 3.04628e-08, 0, 0,
 
-0, 0, 2.981546e-08, 7.249833e-07, 4.468747e-06, 8.409227e-06, 4.554846e-06, 
7.034199e-07, 3.002102e-08, 0, 0,
+0, 3.009687e-08, 2.556034e-06, 5.936867e-05, 0.0003808578, 0.0007126221, 
0.0003827095, 5.902729e-05, 2.553342e-06, 2.978137e-08, 0,
 
-0, 3.054498e-08, 2.614154e-06, 5.891601e-05, 0.0003810036, 0.000708165, 
0.0003842406, 5.963722e-05, 2.618934e-06, 2.990584e-08, 0,
+0, 7.021852e-07, 5.912285e-05, 0.00137637, 0.008863639, 0.01648383, 
0.008855942, 0.001365171, 5.925718e-05, 7.021184e-07, 0,
 
-0, 7.199899e-07, 5.801019e-05, 0.001365485, 0.009023659, 0.01638159, 
0.008892864, 0.001345278, 5.920425e-05, 6.984741e-07, 0,
+0, 4.490787e-06, 0.0003826718, 0.008857355, 0.05742518, 0.1062628, 0.05727194, 
0.008880079, 0.0003826067, 4.478989e-06, 0,
 
-0, 4.584869e-06, 0.0003830431, 0.008917402, 0.05743425, 0.1061489, 0.05746412, 
0.008902563, 0.0003849257, 4.448404e-06, 0,
+2.595735e-08, 8.31301e-06, 0.0007113572, 0.01640853, 0.1061298, 0.1971036, 
0.1062611, 0.01647962, 0.000708363, 8.379878e-06, 2.593496e-08,
 
-2.572769e-08, 8.414041e-06, 0.0007008284, 0.0164456, 0.1055995, 0.19753, 
0.1061855, 0.01653461, 0.0007141303, 8.41643e-06, 2.550312e-08,
+0, 4.516684e-06, 0.0003846966, 0.008860709, 0.05739478, 0.1062216, 0.05725683, 
0.00881713, 0.000383981, 4.473017e-06, 0,
 
-0, 4.582525e-06, 0.0003775396, 0.00898499, 0.05741534, 0.1062144, 0.05700329, 
0.008838926, 0.0003822096, 4.543726e-06, 0,
+0, 6.950547e-07, 5.920586e-05, 0.00137483, 0.00887785, 0.0164709, 0.008855232, 
0.001372743, 5.939038e-05, 7.016624e-07, 0,
 
-0, 6.883925e-07, 6.09077e-05, 0.001339333, 0.008817007, 0.01636454, 
0.008995386, 0.001407854, 6.004799e-05, 7.203602e-07, 0,
+0, 3.006322e-08, 2.587011e-06, 5.92911e-05, 0.0003843824, 0.0007118155, 
0.000386519, 5.974654e-05, 2.585581e-06, 3.048036e-08, 0,
 
-0, 3.095966e-08, 2.575403e-06, 5.89859e-05, 0.0003804447, 0.0007091904, 
0.0003810006, 5.903253e-05, 2.575202e-06, 2.934356e-08, 0,
+0, 0, 3.041056e-08, 7.05225e-07, 4.497418e-06, 8.388542e-06, 4.478833e-06, 
7.018358e-07, 2.995504e-08, 0, 0,
 
-0, 0, 3.040937e-08, 7.018197e-07, 4.543086e-06, 8.296753e-06, 4.434901e-06, 
6.659026e-07, 3.066215e-08, 0, 0,
+0, 0, 0, 0, 0, 2.567377e-08, 0, 0, 0, 0, 0};
 
-0, 0, 0, 0, 0, 2.603901e-08, 0, 0, 0, 0, 0};
 
 #endif