gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master 9285abe 057/113: Imported recent work from mas


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master 9285abe 057/113: Imported recent work from master with corrections
Date: Fri, 16 Apr 2021 10:33:46 -0400 (EDT)

branch: master
commit 9285abe85073f1d33d8396874b2abaf8c5bb3974
Merge: ec3c102 b4c2b5a
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    Imported recent work from master with corrections
    
    Prior to the last commit in master, there was no conflict. But with the
    last commit, we actually removed the extra layer of (zero-valued) pixels in
    the profiles created by MakeProfiles (through a correction in
    `gal_box_bound_ellipse'). This change also needed to be corrected in the
    `gal_box_bound_ellipsoid' function of this branch.
    
    So with this merge, we also made the respective correction here. Also,
    similar to the 2D kernel of the master branch, we also re-made the default
    3D kernel for NoiseChisel and Segment.
---
 NEWS                              |  10 +-
 bin/arithmetic/arithmetic.c       | 345 ++++++++++++++++++++++++--------------
 bin/arithmetic/arithmetic.h       |   3 +
 bin/mkprof/oneprofile.c           |  24 ++-
 doc/gnuastro.texi                 |  68 ++++++--
 lib/box.c                         |   6 +-
 lib/gnuastro-internal/kernel-2d.h |  50 +++---
 lib/gnuastro-internal/kernel-3d.h | 114 ++++++-------
 lib/interpolate.c                 |   4 +-
 9 files changed, 394 insertions(+), 230 deletions(-)

diff --git a/NEWS b/NEWS
index f0014f8..3ab0703 100644
--- a/NEWS
+++ b/NEWS
@@ -16,10 +16,13 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
       in charge of detection.
 
   Arithmetic:
-    - connected-components: label the connected elements of the input.
+    - erode: Erode the foreground of a binary dataset.
+    - dilate: Dilate the foreground of a binary dataset.
     - filter-sigclip-mean: sigma-clipped, mean filter operator.
     - filter-sigclip-median: sigma-clipped, median filter operator.
+    - connected-components: label the connected elements of the input.
     - invert: subtract the maximum of unsigned types (absorption to emission).
+    - interpolate-medianngb: Interpolate (only blanks) with nearest neighbors.
 
   ConvertType:
     - TIFF images can also be used as input.
@@ -105,6 +108,11 @@ GNU Astronomy Utilities NEWS                          -*- 
outline -*-
 
 ** Changed features
 
+  Arithmetic:
+    - filter-mean: a blank value in the input can be non-blank in the
+         output when non-blank elements present in filter.
+    - filter-median: similar to filter-mean.
+
   Fits:
     --history: can be called/written multiple times in one run.
     --comment: can be called/written multiple times in one run.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index baef282..ac1a183 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 #include <gnuastro/arithmetic.h>
+#include <gnuastro/interpolate.h>
 
 #include <gnuastro-internal/checkset.h>
 
@@ -179,121 +180,110 @@ arithmetic_filter(void *in_prm)
       /* For easy reading, put the index in `ind'. */
       ind=tprm->indexs[i];
 
-      /* If we are on a blank element, then just set the output to blank
-         also. */
-      if( afp->hasblank
-          && gal_blank_is(gal_data_ptr_increment(input->array, ind,
-                                                 input->type), input->type) )
-        gal_blank_write(gal_data_ptr_increment(afp->out->array, ind,
-                                               afp->out->type),
-                        afp->out->type);
-      else
-        {
-          /* Get the coordinate of the pixel. */
-          gal_dimension_index_to_coord(ind, ndim, dsize, coord);
+      /* Get the coordinate of the pixel. */
+      gal_dimension_index_to_coord(ind, ndim, dsize, coord);
 
-          /* See which dimensions need trimming. */
-          tile->size=1;
-          for(j=0;j<ndim;++j)
+      /* See which dimensions need trimming. */
+      tile->size=1;
+      for(j=0;j<ndim;++j)
+        {
+          /* Estimate the coordinate of the filter's starting point. Note
+             that we are dealing with size_t (unsigned int) type here, so
+             there are no negatives. A negative result will produce an
+             extremely large number, so instead of checking for negative,
+             we can just see if the result of a subtraction is less than
+             the width of the input. */
+          if( (coord[j] - hnfsize[j] > dsize[j])
+              || (coord[j] + hpfsize[j] >= dsize[j]) )
             {
-              /* Estimate the coordinate of the filter's starting
-                 point. Note that we are dealing with size_t (unsigned int)
-                 type here, so there are no negatives. A negative result
-                 will produce an extremely large number, so instead of
-                 checking for negative, we can just see if the result of a
-                 subtraction is less than the width of the input. */
-              if( (coord[j] - hnfsize[j] > dsize[j])
-                  || (coord[j] + hpfsize[j] >= dsize[j]) )
-                {
-                  start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
-                               ? 0 : coord[j] - hnfsize[j] );
-                  end[j]   = ( (coord[j] + hpfsize[j] >= dsize[j])
-                               ? dsize[j]
-                               : coord[j] + hpfsize[j] + 1);
-                  tsize[j] = end[j] - start[j];
-                }
-              else  /* NOT on the edge (given requested filter width). */
-                {
-                  tsize[j] = fsize[j];
-                  start[j] = coord[j] - hnfsize[j];
-                }
-              tile->size *= tsize[j];
+              start[j] = ( (coord[j] - hnfsize[j] > dsize[j])
+                           ? 0 : coord[j] - hnfsize[j] );
+              end[j]   = ( (coord[j] + hpfsize[j] >= dsize[j])
+                           ? dsize[j]
+                           : coord[j] + hpfsize[j] + 1);
+              tsize[j] = end[j] - start[j];
             }
-
-          /* For a test.
-             printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
-             printf("\tstart: %zu, %zu\n", start[1]+1, start[0]+1);
-             printf("\ttsize: %zu, %zu\n", tsize[1], tsize[0]);
-          */
-
-          /* Set the tile's starting pointer. */
-          index=gal_dimension_coord_to_index(ndim, dsize, start);
-          tile->array=gal_data_ptr_increment(input->array, index,
-                                             input->type);
-
-          /* Do the necessary calculation. */
-          switch(afp->operator)
+          else  /* NOT on the edge (given requested filter width). */
             {
-            case ARITHMETIC_OP_FILTER_MEDIAN:
-              result=gal_statistics_median(tile, 0);
-              break;
-
+              tsize[j] = fsize[j];
+              start[j] = coord[j] - hnfsize[j];
+            }
+          tile->size *= tsize[j];
+        }
 
-            case ARITHMETIC_OP_FILTER_MEAN:
-              result=gal_statistics_mean(tile);
-              break;
+      /* For a test.
+         printf("coord: %zu, %zu\n", coord[1]+1, coord[0]+1);
+         printf("\tstart: %zu, %zu\n", start[1]+1, start[0]+1);
+         printf("\ttsize: %zu, %zu\n", tsize[1], tsize[0]);
+      */
 
+      /* Set the tile's starting pointer. */
+      index=gal_dimension_coord_to_index(ndim, dsize, start);
+      tile->array=gal_data_ptr_increment(input->array, index,
+                                         input->type);
 
-            case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
-            case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
-              /* Find the sigma-clipped results. */
-              sigclip=gal_statistics_sigma_clip(tile, afp->sclip_multip,
-                                                afp->sclip_param, 0, 1);
+      /* Do the necessary calculation. */
+      switch(afp->operator)
+        {
+        case ARITHMETIC_OP_FILTER_MEDIAN:
+          result=gal_statistics_median(tile, 0);
+          break;
 
-              /* Set the required index. */
-              switch(afp->operator)
-                {
-                case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:   sind = 2; break;
-                case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN: sind = 1; break;
-                default:
-                  error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
-                        "%s to fix the problem. The `afp->operator' value "
-                        "%d is not recognized as sigma-clipped median or "
-                        "mean", __func__, PACKAGE_BUGREPORT, afp->operator);
-                }
 
-              /* Allocate the output and write the value into it. */
-              result=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
-                                    0, -1, NULL, NULL, NULL);
-              ((float *)(result->array))[0] =
-                ((float *)(sigclip->array))[sind];
+        case ARITHMETIC_OP_FILTER_MEAN:
+          result=gal_statistics_mean(tile);
+          break;
 
-              /* Clean up. */
-              gal_data_free(sigclip);
-              break;
 
+        case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:
+        case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN:
+          /* Find the sigma-clipped results. */
+          sigclip=gal_statistics_sigma_clip(tile, afp->sclip_multip,
+                                            afp->sclip_param, 0, 1);
 
+          /* Set the required index. */
+          switch(afp->operator)
+            {
+            case ARITHMETIC_OP_FILTER_SIGCLIP_MEAN:   sind = 2; break;
+            case ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN: sind = 1; break;
             default:
-              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
-                    "to fix the problem. `afp->operator' code %d is not "
-                    "recognized", PACKAGE_BUGREPORT, __func__,
-                    afp->operator);
+              error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at "
+                    "%s to fix the problem. The `afp->operator' value "
+                    "%d is not recognized as sigma-clipped median or "
+                    "mean", __func__, PACKAGE_BUGREPORT, afp->operator);
             }
 
-          /* Make sure the output array type and result's type are the
-             same. */
-          if(result->type!=afp->out->type)
-            result=gal_data_copy_to_new_type_free(result, afp->out->type);
+          /* Allocate the output and write the value into it. */
+          result=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL,
+                                0, -1, NULL, NULL, NULL);
+          ((float *)(result->array))[0] =
+            ((float *)(sigclip->array))[sind];
 
+          /* Clean up. */
+          gal_data_free(sigclip);
+          break;
 
-          /* Copy the result into the output array. */
-          memcpy(gal_data_ptr_increment(afp->out->array, ind,
-                                        afp->out->type),
-                 result->array, gal_type_sizeof(afp->out->type));
 
-          /* Clean up for this pixel. */
-          gal_data_free(result);
+        default:
+          error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s "
+                "to fix the problem. `afp->operator' code %d is not "
+                "recognized", PACKAGE_BUGREPORT, __func__,
+                afp->operator);
         }
+
+      /* Make sure the output array type and result's type are the
+         same. */
+      if(result->type!=afp->out->type)
+        result=gal_data_copy_to_new_type_free(result, afp->out->type);
+
+
+      /* Copy the result into the output array. */
+      memcpy(gal_data_ptr_increment(afp->out->array, ind,
+                                    afp->out->type),
+             result->array, gal_type_sizeof(afp->out->type));
+
+      /* Clean up for this pixel. */
+      gal_data_free(result);
     }
 
 
@@ -513,54 +503,107 @@ wrapper_for_filter(struct arithmeticparams *p, char 
*token, int operator)
 /***************************************************************/
 /*************            Other functions          *************/
 /***************************************************************/
-static void
-arithmetic_connected_components(struct arithmeticparams *p, char *token)
+static int
+arithmetic_binary_conn_sanity_checks(gal_data_t *in, gal_data_t *conn,
+                                     char *operator)
 {
   int conn_int;
-  gal_data_t *out=NULL;
-
-  /* Pop the two necessary operands. */
-  gal_data_t *conn = operands_pop(p, token);
-  gal_data_t *in   = operands_pop(p, token);
 
   /* Do proper sanity checks on `conn'. */
   if(conn->size!=1)
-    error(EXIT_FAILURE, 0, "the first popped operand to "
-          "`connected-components' must be a single number. However, it has "
-          "%zu elements", conn->size);
+    error(EXIT_FAILURE, 0, "the first popped operand to `%s' must be a "
+          "single number. However, it has %zu elements", operator,
+          conn->size);
   if(conn->type==GAL_TYPE_FLOAT32 || conn->type==GAL_TYPE_FLOAT64)
-    error(EXIT_FAILURE, 0, "the first popped operand to "
-          "`connected-components' is the connectivity (a value between 1 "
-          "and the number of dimensions) therefore, it must NOT be a "
-          "floating point");
+    error(EXIT_FAILURE, 0, "the first popped operand to `%s' is the "
+          "connectivity (a value between 1 and the number of dimensions) "
+          "therefore, it must NOT be a floating point", operator);
 
   /* Convert the connectivity value to a 32-bit integer and read it in and
      make sure it is not larger than the number of dimensions. */
   conn=gal_data_copy_to_new_type_free(conn, GAL_TYPE_INT32);
   conn_int = *((int32_t *)(conn->array));
   if(conn_int>in->ndim)
-    error(EXIT_FAILURE, 0, "the first popped operand of "
-          "`connected-components' (%d) is larger than the number of "
-          "dimensions in the second-popped operand (%zu)", conn_int,
-          in->ndim);
+    error(EXIT_FAILURE, 0, "the first popped operand of `%s' (%d) is "
+          "larger than the number of dimensions in the second-popped "
+          "operand (%zu)", operator, conn_int, in->ndim);
 
   /* Make sure the array has an unsigned 8-bit type. */
   if(in->type!=GAL_TYPE_UINT8)
-    error(EXIT_FAILURE, 0, "the second popped operand of "
-          "`connected-components' doesn't have an 8-bit unsigned "
-          "integer type. It must be a binary dataset (only being equal "
-          "to zero is checked). You can use the `uint8' operator to "
-          "convert the type of this operand.");
+    error(EXIT_FAILURE, 0, "the second popped operand of `%s' doesn't "
+          "have an 8-bit unsigned integer type. It must be a binary "
+          "dataset (only being equal to zero is checked). You can use "
+          "the `uint8' operator for type conversion", operator);
+
+  /* Clean up and return the integer value of `conn'. */
+  gal_data_free(conn);
+  return conn_int;
+}
+
+
+
+
+
+static void
+arithmetic_erode_dilate(struct arithmeticparams *p, char *token, int op)
+{
+  int conn_int;
+
+  /* Pop the two necessary operands. */
+  gal_data_t *conn = operands_pop(p, token);
+  gal_data_t *in   = operands_pop(p, token);
+
+  /* Do the sanity checks and  */
+  switch(op)
+    {
+    case ARITHMETIC_OP_ERODE:
+      conn_int=arithmetic_binary_conn_sanity_checks(in, conn, "erode");
+      gal_binary_erode(in, 1, conn_int, 1);
+      break;
+
+    case ARITHMETIC_OP_DILATE:
+      conn_int=arithmetic_binary_conn_sanity_checks(in, conn, "dilate");
+      gal_binary_dilate(in, 1, conn_int, 1);
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix the "
+            "problem. The operator code %d not recognized", __func__,
+            PACKAGE_BUGREPORT, op);
+    }
+
+  /* Push the result onto the stack. */
+  operands_add(p, NULL, in);
+
+  /* Recall that`conn' was freed in the sanity check. */
+}
+
+
+
+
+
+static void
+arithmetic_connected_components(struct arithmeticparams *p, char *token)
+{
+  int conn_int;
+  gal_data_t *out=NULL;
+
+  /* Pop the two necessary operands. */
+  gal_data_t *conn = operands_pop(p, token);
+  gal_data_t *in   = operands_pop(p, token);
+
+  /* Basic sanity checks. */
+  conn_int=arithmetic_binary_conn_sanity_checks(in, conn,
+                                                "connected-components");
 
   /* Do the connected components labeling. */
-  gal_binary_connected_components(in, &out, 1);
+  gal_binary_connected_components(in, &out, conn_int);
 
   /* Push the result onto the stack. */
   operands_add(p, NULL, out);
 
-  /* Clean up. */
+  /* Clean up (`conn' was freed in the sanity check). */
   gal_data_free(in);
-  gal_data_free(conn);
 }
 
 
@@ -600,6 +643,49 @@ arithmetic_invert(struct arithmeticparams *p, char *token)
 
 
 
+static void
+arithmetic_interpolate(struct arithmeticparams *p, char *token)
+{
+  int num_int;
+  gal_data_t *interpolated;
+
+  /* First pop the number of nearby neighbors.*/
+  gal_data_t *num = operands_pop(p, token);
+
+  /* Then pop the actual dataset to interpolate. */
+  gal_data_t *in = operands_pop(p, token);
+
+  /* Do proper sanity checks on `num'. */
+  if(num->size!=1)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`interpolate-medianngb' must be a single number. However, "
+          "it has %zu elements", num->size);
+  if(num->type==GAL_TYPE_FLOAT32 || num->type==GAL_TYPE_FLOAT64)
+    error(EXIT_FAILURE, 0, "the first popped operand to "
+          "`interpolate-medianngb' is the number of nearby neighbors (a "
+          "counter, an integer). It must NOT be a floating point.\n\n"
+          "If its already an integer, but in a floating point container, "
+          "you can use the `int32' operator to convert it to a 32-bit "
+          "integer for example");
+
+  /* Convert the given number to a 32-bit integer and read it in. */
+  num=gal_data_copy_to_new_type_free(num, GAL_TYPE_INT32);
+  num_int = *((int32_t *)(num->array));
+
+  /* Call the interpolation function. */
+  interpolated=gal_interpolate_close_neighbors(in, NULL, num_int,
+                                               p->cp.numthreads, 1, 0);
+
+  /* Clean up and push the interpolated array onto the stack. */
+  gal_data_free(in);
+  gal_data_free(num);
+  operands_add(p, NULL, interpolated);
+}
+
+
+
+
+
 
 
 
@@ -778,10 +864,16 @@ reversepolish(struct arithmeticparams *p)
             { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEAN;   nop=0;  }
           else if (!strcmp(token->v, "filter-sigclip-median"))
             { op=ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN; nop=0;  }
+          else if (!strcmp(token->v, "erode"))
+            { op=ARITHMETIC_OP_ERODE;                 nop=0;  }
+          else if (!strcmp(token->v, "dilate"))
+            { op=ARITHMETIC_OP_DILATE;                nop=0;  }
           else if (!strcmp(token->v, "connected-components"))
             { op=ARITHMETIC_OP_CONNECTED_COMPONENTS;  nop=0;  }
           else if (!strcmp(token->v, "invert"))
             { op=ARITHMETIC_OP_INVERT;                nop=0;  }
+          else if (!strcmp(token->v, "interpolate-medianngb"))
+            { op=ARITHMETIC_OP_INTERPOLATE_MEDIANNGB; nop=0;  }
 
 
           /* Finished checks with known operators */
@@ -856,6 +948,11 @@ reversepolish(struct arithmeticparams *p)
                   wrapper_for_filter(p, token->v, op);
                   break;
 
+                case ARITHMETIC_OP_ERODE:
+                case ARITHMETIC_OP_DILATE:
+                  arithmetic_erode_dilate(p, token->v, op);
+                  break;
+
                 case ARITHMETIC_OP_CONNECTED_COMPONENTS:
                   arithmetic_connected_components(p, token->v);
                   break;
@@ -864,6 +961,10 @@ reversepolish(struct arithmeticparams *p)
                   arithmetic_invert(p, token->v);
                   break;
 
+                case ARITHMETIC_OP_INTERPOLATE_MEDIANNGB:
+                  arithmetic_interpolate(p, token->v);
+                  break;
+
                 default:
                   error(EXIT_FAILURE, 0, "%s: a bug! please contact us at "
                         "%s to fix the problem. The code %d is not "
diff --git a/bin/arithmetic/arithmetic.h b/bin/arithmetic/arithmetic.h
index 59b2170..e5eca69 100644
--- a/bin/arithmetic/arithmetic.h
+++ b/bin/arithmetic/arithmetic.h
@@ -35,8 +35,11 @@ enum arithmetic_prog_operators
   ARITHMETIC_OP_FILTER_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEAN,
   ARITHMETIC_OP_FILTER_SIGCLIP_MEDIAN,
+  ARITHMETIC_OP_ERODE,
+  ARITHMETIC_OP_DILATE,
   ARITHMETIC_OP_CONNECTED_COMPONENTS,
   ARITHMETIC_OP_INVERT,
+  ARITHMETIC_OP_INTERPOLATE_MEDIANNGB,
 };
 
 
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 3b1fc28..79b3d29 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -195,8 +195,10 @@ oneprofile_r_circle(size_t index, struct mkonthread *mkp)
 float
 oneprofile_randompoints(struct mkonthread *mkp)
 {
+  double r_before=mkp->r;
   double range[3], sum=0.0f;
   size_t i, j, numrandom=mkp->p->numrandom, ndim=mkp->p->ndim;
+  double coord_before[3]={mkp->coord[0], mkp->coord[1], mkp->coord[2]};
 
   /* Set the range in each dimension. */
   for(i=0;i<ndim;++i)
@@ -211,6 +213,16 @@ 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];
+  mkp->coord[2]=coord_before[2];
   return sum/numrandom;
 }
 
@@ -399,7 +411,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;
@@ -412,6 +424,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; }
 
@@ -463,6 +481,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 70c318d..cbca5a0 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -10116,6 +10116,13 @@ The final effect of mean filtering is to smooth the 
input image, it is
 essentially a convolution with a kernel that has identical values for all
 its pixels (is flat), see @ref{Convolution process}.
 
+Note that blank pixels will also be affected by this operator: if there are
+any non-blank elements in the box surrounding a blank pixel, in the
+flitered image, it will have the mean of the non-blank elements, therefore
+it won't be blank any more. If blank elements are important for your
+analysis, you can use the @code{isblank} with the @code{where} operator to
+set them back to blank after filtering.
+
 @item filter-median
 Apply @url{https://en.wikipedia.org/wiki/Median_filter, median filtering}
 on the input dataset. This is very similar to @command{filter-mean}, except
@@ -10160,6 +10167,35 @@ dataset. This operator and its necessary operands are 
almost identical to
 median value (which is less affected by outliers than the mean) is added
 back to the stack.
 
+@item interpolate-medianngb
+Interpolate all the blank elements of the second popped operand with the
+median of its nearest non-blank neighbors. The number of the nearest
+non-blank neighbors used to calculate the median is given by the first
+popped operand. Note that the distance of the nearest non-blank neighbors
+is irrelevant in this interpolation.
+
+@item erode
+@cindex Erosion
+Erode the foreground pixels (with value @code{1}) of the input dataset
+(second popped operand). The first popped operand is the connectivity (see
+description in @command{connected-components}). Erosion is simply a
+flipping of all foreground pixels (to background; with value @code{0}) that
+are ``touching'' background pixels. ``Touching'' is defined by the
+connectivity. In effect, this carves off the outer borders of the
+foreground, making them thinner. This operator assumes a binary dataset
+(all pixels are @code{0} and @code{1}).
+
+@item dilate
+@cindex Dilation
+Dilate the foreground pixels (with value @code{1}) of the input dataset
+(second popped operand). The first popped operand is the connectivity (see
+description in @command{connected-components}). Erosion is simply a
+flipping of all background pixels (with value @code{0}) to foreground that
+are ``touching'' foreground pixels. ``Touching'' is defined by the
+connectivity. In effect, this expands the outer borders of the
+foreground. This operator assumes a binary dataset (all pixels are @code{0}
+and @code{1}).
+
 @item connected-components
 @cindex Connected components
 Find the connected components in the input dataset (second popped
@@ -24717,22 +24753,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
+@code{extent} 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 integer height and width of that box assuming the
-center of the ellipse is in the box center. @code{a} is the ellipse major
-axis, @code{b} 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
+@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.
 @end deftypefun
 
 @deftypefun void gal_box_bound_ellipsoid_extent (double @code{*semiaxes}, 
double @code{*euler_deg}, double @code{*extent})
@@ -25679,7 +25716,7 @@ nature. The produced values are also always within the 
range of the known
 values and strong outliers do not get created. We will hopefully implement
 other methods too (wrappers around the GNU Scientific Library's very
 complete set of functions), but currently the developers are too busy. So
-if you do have the chance to help your contribution would be very welcome
+if you do have the chance to help, your contribution would be very welcome
 and appreciated.
 
 @deftypefun {gal_data_t *} gal_interpolate_close_neighbors (gal_data_t 
@code{*input}, struct gal_tile_two_layer_params @code{*tl}, size_t 
@code{numneighbors}, size_t @code{numthreads}, int @code{onlyblank}, int 
@code{aslinkedlist})
@@ -26042,7 +26079,8 @@ main(void)
   /* Read the image into memory as a float32 data type. We are using
    * `-1' for `minmapsize' to ensure that the image is read into
    * memory. */
-  p.image=gal_fits_img_read_to_type(filename, hdu, GAL_TYPE_FLOAT32, -1);
+  p.image=gal_fits_img_read_to_type(filename, hdu, GAL_TYPE_FLOAT32,
+                                    -1);
 
 
   /* Print some basic information before the actual contents: */
diff --git a/lib/box.c b/lib/box.c
index d42c37c..407700b 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;
 }
 
 
@@ -261,7 +261,7 @@ gal_box_bound_ellipsoid(double *semiaxes, double 
*euler_deg, long *width)
 
   /* Find the width along each dimension. */
   for(i=0;i<3;++i)
-    width[i] = 2 * ((long)extent[i] + 1) + 1;
+    width[i] = 2 * ( (long)extent[i] ) + 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
diff --git a/lib/gnuastro-internal/kernel-3d.h 
b/lib/gnuastro-internal/kernel-3d.h
index 30ca97e..5b10cf8 100644
--- a/lib/gnuastro-internal/kernel-3d.h
+++ b/lib/gnuastro-internal/kernel-3d.h
@@ -23,17 +23,32 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #ifndef __GAL_KERNEL3D_H__
 #define __GAL_KERNEL3D_H__
 
-/* How to build this kernel.
+/* How to produce this kernel
+   ==========================
 
-   Run MakeProfiles with the following commands:
+   Below, the steps necessary to easily create the C contents of this file
+   are described. The first step can be modified to change the default
+   kernel properties and put the new contents into this file.
 
-       $ export GSL_RNG_SEED=1
-       $ export GSL_RNG_TYPE=ranlxs2
-       $ astmkprof --kernel=gaussian-3d,1.5,5,0.5 --oversample=1 --envseed
+   Make the kernel
+   ---------------
 
-   The resulting fits file is converted to plain text with the this C
-   program to generate the `kernel-3d.h' header file which is then put
-   under these comments.
+   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.
+
+     export GSL_RNG_SEED=1
+     export GSL_RNG_TYPE=ranlxs2
+     astmkprof --kernel=gaussian-3d,1.5,5,0.5 --oversample=1 --envseed \
+               --numrandom=100000
+
+   Convert it to C code
+   --------------------
+
+   Put the following C program into a file called `kernel.c'.
 
        #include <stdio.h>
        #include <stdlib.h>
@@ -77,62 +92,39 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
          return EXIT_SUCCESS;
        }
 
-   Assuming this C program is in a file named `kernel.c', it can be
-   compiled, run and saved into `kernel-3d.h' with the command below:
+   Run the C program
+   -----------------
 
-       $ astbuildprog -q kernel.c > kernel-3d.h
+   We can now compile and run that C program and put the outputs in
+   `kernel.c'. Once its created, copy the contents of `kernel-3d.h' after
+   these comments.
 
+     $ astbuildprog -q kernel.c > kernel-3d.h
  */
 
-size_t kernel_3d_dsize[3]={5, 9, 9};
-float kernel_3d[405]={0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 5.293481e-07, 1.391654e-06, 5.283476e-07, 0, 0, 0,
-0, 0, 5.023404e-06, 0.0001093707, 0.0003028791, 0.0001103754, 5.206822e-06, 0, 
0,
-0, 5.27038e-07, 0.0001140605, 0.002517939, 0.006957374, 0.002468025, 
0.0001131454, 5.31773e-07, 0,
-0, 1.481738e-06, 0.0002994444, 0.006708808, 0.01894199, 0.006776441, 
0.0003116253, 1.493134e-06, 0,
-0, 5.293459e-07, 0.0001081432, 0.002473192, 0.006704316, 0.002578867, 
0.0001105525, 4.866729e-07, 0,
-0, 0, 4.985141e-06, 0.000110672, 0.0002885369, 0.0001184523, 4.688276e-06, 0, 
0,
-0, 0, 0, 5.610833e-07, 1.336052e-06, 5.216995e-07, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 3.462786e-07, 7.853481e-06, 2.153742e-05, 7.564613e-06, 3.611804e-07, 0, 
0,
-0, 3.52808e-07, 7.161698e-05, 0.001676089, 0.004618002, 0.001656355, 
7.414426e-05, 3.43867e-07, 0,
-0, 7.990315e-06, 0.001704389, 0.03857445, 0.1031222, 0.03779064, 0.001703231, 
7.975499e-06, 0,
-0, 2.231775e-05, 0.00473238, 0.1043504, 0.2858114, 0.10241, 0.004618072, 
2.193002e-05, 0,
-0, 7.621142e-06, 0.001647367, 0.03755038, 0.1036048, 0.03796006, 0.001672143, 
7.758124e-06, 0,
-0, 3.68421e-07, 7.659229e-05, 0.001695716, 0.004599371, 0.001670764, 
7.450273e-05, 3.481919e-07, 0,
-0, 0, 3.521386e-07, 7.924394e-06, 2.207155e-05, 7.889852e-06, 3.53245e-07, 0, 
0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 5.182329e-07, 1.397928e-06, 5.17901e-07, 0, 0, 0,
-0, 0, 4.812169e-06, 0.0001135803, 0.0003043477, 0.000109481, 4.974718e-06, 0, 
0,
-0, 5.377616e-07, 0.0001073098, 0.002424047, 0.006817353, 0.002492242, 
0.0001147361, 5.33439e-07, 0,
-0, 1.44089e-06, 0.0002933802, 0.006909565, 0.01873765, 0.006651767, 
0.0003042257, 1.366647e-06, 0,
-0, 5.166781e-07, 0.0001119049, 0.002454791, 0.006940499, 0.002578837, 
0.0001061706, 5.329928e-07, 0,
-0, 0, 4.656712e-06, 0.0001124511, 0.0003010639, 0.0001131462, 4.786908e-06, 0, 
0,
-0, 0, 0, 5.303103e-07, 1.363214e-06, 5.122773e-07, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0,
-0, 0, 0, 0, 0, 0, 0, 0, 0};
+size_t kernel_3d_dsize[3]={3, 7, 7};
+float kernel_3d[147]={0, 0, 5.219212e-07, 1.43663e-06, 5.288961e-07, 0, 0,
+0, 4.890969e-06, 0.0001122713, 0.0003051736, 0.0001101779, 4.908836e-06, 0,
+5.281005e-07, 0.0001115758, 0.002482525, 0.006792462, 0.002476586, 
0.0001099998, 5.217535e-07,
+1.451268e-06, 0.0003033727, 0.006850741, 0.01871265, 0.006743792, 
0.0003061892, 1.43298e-06,
+5.209756e-07, 0.000110545, 0.002490561, 0.006798376, 0.002503618, 
0.0001109581, 5.295083e-07,
+0, 4.884327e-06, 0.0001122016, 0.0003074409, 0.0001103678, 4.984748e-06, 0,
+0, 0, 5.224761e-07, 1.425723e-06, 5.33545e-07, 0, 0,
+
+0, 3.538044e-07, 7.891201e-06, 2.166429e-05, 7.979216e-06, 3.570717e-07, 0,
+3.536664e-07, 7.511148e-05, 0.001693119, 0.004623666, 0.001691994, 
7.563465e-05, 3.541075e-07,
+7.921932e-06, 0.001687588, 0.03788469, 0.1038283, 0.03776859, 0.001679332, 
7.979274e-06,
+2.200398e-05, 0.004617298, 0.1038528, 0.2849551, 0.1039089, 0.004635326, 
2.165271e-05,
+7.955934e-06, 0.001697289, 0.0379024, 0.1037254, 0.03768558, 0.001666544, 
7.878737e-06,
+3.523428e-07, 7.45342e-05, 0.001689582, 0.00461954, 0.001684659, 7.568914e-05, 
3.499421e-07,
+0, 3.56209e-07, 7.915472e-06, 2.17837e-05, 7.854093e-06, 3.555503e-07, 0,
+
+0, 0, 5.263017e-07, 1.431613e-06, 5.12058e-07, 0, 0,
+0, 4.907304e-06, 0.0001125078, 0.00030509, 0.000111018, 5.084412e-06, 0,
+5.256628e-07, 0.0001103357, 0.002493239, 0.006835639, 0.002478384, 
0.0001131122, 5.260213e-07,
+1.442301e-06, 0.0003005426, 0.006833205, 0.01862562, 0.006816466, 0.000302246, 
1.436046e-06,
+5.298979e-07, 0.0001107397, 0.002476231, 0.006837885, 0.00252707, 
0.0001124105, 5.305631e-07,
+0, 5.012419e-06, 0.0001110889, 0.0003017522, 0.0001124517, 4.790862e-06, 0,
+0, 0, 5.104474e-07, 1.45196e-06, 5.191671e-07, 0, 0};
 
 #endif
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 1276aaa..8d66deb 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -333,8 +333,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
   int permute=(tl && tl->totchannels>1 && tl->workoverch);
 
 
-  /* If there are no blank values in the array we should only fill blank
-     values, then simply copy the input and abort. */
+  /* If there are no blank values in the array, AND we should only fill
+     blank values, then simply copy the input and abort. */
   if( (input->flag | GAL_DATA_FLAG_BLANK_CH)     /* Zero bit is meaningful.*/
       && !(input->flag | GAL_DATA_FLAG_HASBLANK) /* There are no blanks.   */
       && onlyblank )                             /* Only interpolate blank.*/



reply via email to

[Prev in Thread] Current Thread [Next in Thread]