[gnuastro-commits] master 8372486 020/113: NoiseChisel's detection complete in 3D

[Mohammad Akhlaghi]

branch: master
commit 83724869b01742b41f40dd71009e87dc726ec438
Author: Mohammad Akhlaghi <akhlaghi@gnu.org>
Commit: Mohammad Akhlaghi <akhlaghi@gnu.org>

    NoiseChisel's detection complete in 3D
    
    NoiseChisel's full detection process is now complete for 3D datacubes. For
    the process the following changes have been made:
    
     - A new NoiseChisel configuration specifically for 3D datacubes.
    
     - Some corrections in the extension names of check images.
    
     - Documentation was updated.
    
     - A general (static/internal) function `binary_erode_dilate_general' was
       defined in the `lib/binary.c' to do erosion and dilation on any
       dimensionality (based on `GAL_DIMENSION_NEIGHBOR_OP').
    
     - Several checks/updates were made in the user interface steps (checks of
       inputs).
---
 bin/noisechisel/astnoisechisel-3d.conf |  4 ++
 bin/noisechisel/detection.c            | 65 +++++++++++++++-------
 bin/noisechisel/noisechisel.c          | 16 +++---
 bin/noisechisel/sky.c                  |  6 ++-
 bin/noisechisel/ui.c                   | 89 ++++++++++++++++++++++++++----
 doc/gnuastro.texi                      | 98 ++++++++++++++++++----------------
 lib/binary.c                           | 49 +++++++++++++++--
 7 files changed, 240 insertions(+), 87 deletions(-)

diff --git a/bin/noisechisel/astnoisechisel-3d.conf 
b/bin/noisechisel/astnoisechisel-3d.conf
index f7cb3a2..58026c6 100644
--- a/bin/noisechisel/astnoisechisel-3d.conf
+++ b/bin/noisechisel/astnoisechisel-3d.conf
@@ -21,3 +21,7 @@
  numchannels    1,1,1
  tilesize       20,20,20
  largetilesize  50,50,50
+
+# Detection
+ erodengb       6
+ openingngb     18
\ No newline at end of file
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index a8df736..ff25c1e 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -54,6 +54,7 @@ detection_initial(struct noisechiselparams *p)
 {
   char *msg;
   uint8_t *b, *bf;
+  int connectivity=0;
   struct timeval t0, t1;
 
 
@@ -75,9 +76,25 @@ detection_initial(struct noisechiselparams *p)
     }
 
 
+  /* Set the connectivity of the erosion. */
+  switch(p->binary->ndim)
+    {
+    case 2:
+      connectivity = p->erodengb==4 ? 1 : 2;
+      break;
+    case 3:
+      connectivity = p->erodengb==6 ? 1 : (p->erodengb==18 ? 2 : 3);
+      break;
+    default:
+      error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
+            "the problem. %zu not a valid value to p->binary->ndim",
+            __func__, PACKAGE_BUGREPORT, p->binary->ndim);
+    }
+
+
   /* Erode the image. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
-  gal_binary_erode(p->binary, p->erode, p->erodengb==4 ? 1 : 2, 1);
+  gal_binary_erode(p->binary, p->erode, connectivity, 1);
   if(!p->cp.quiet)
     {
       asprintf(&msg, "Eroded %zu time%s (%zu-connectivity).", p->erode,
@@ -114,7 +131,7 @@ detection_initial(struct noisechiselparams *p)
   p->numinitialdets=gal_binary_connected_components(p->binary, &p->olabel, 1);
   if(p->detectionname)
     {
-      p->olabel->name="OPENED_AND_LABELED";
+      p->olabel->name="OPENED-AND-LABELED";
       gal_fits_img_write(p->olabel, p->detectionname, NULL, PROGRAM_STRING);
       p->olabel->name=NULL;
     }
@@ -439,9 +456,9 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
   size_t tablen=num+1;
   gal_data_t *sn, *snind;
   int32_t *plabend, *indarr=NULL;
-  double ave, err, *xy, *brightness;
-  size_t ind, ndim=p->input->ndim, xyncols=1+ndim;
-  size_t i, *area, counter=0, *dsize=p->input->dsize;
+  double ave, err, *pos, *brightness;
+  size_t ind, ndim=p->input->ndim, pcols=1+ndim;
+  size_t i, j, *area, counter=0, *dsize=p->input->dsize;
   float *img=p->input->array, *f=p->input->array, *ff=f+p->input->size;
   int32_t *plab = worklab->array, *dlab = s0d1D2 ? NULL : p->olabel->array;
   size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim, __func__,
@@ -456,9 +473,9 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
     error(EXIT_FAILURE, 0, "%s: only a NaN value is recognized for blank "
           "floating point data types, the blank value is defined to be %f",
           __func__, GAL_BLANK_FLOAT32);
-  if(ndim!=2)
-    error(EXIT_FAILURE, 0, "%s: only 2D datasets are acceptable, your input "
-          "is %zu dimensions", __func__, ndim);
+  if(ndim!=2 && ndim!=3)
+    error(EXIT_FAILURE, 0, "%s: only 2D images or 3D datacubes are "
+          "acceptable, but input has %zu dimensions", __func__, ndim);
 
 
   /* Allocate all the necessary arrays, note that since we want to put each
@@ -468,8 +485,8 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
                                      "area");
   brightness = gal_data_calloc_array(GAL_TYPE_FLOAT64, tablen, __func__,
                                      "brightness");
-  xy         = gal_data_calloc_array(GAL_TYPE_FLOAT64, xyncols*tablen,
-                                     __func__, "xy");
+  pos        = gal_data_calloc_array(GAL_TYPE_FLOAT64, pcols*tablen,
+                                     __func__, "pos");
   flag       = ( s0d1D2==0
                  ? gal_data_calloc_array(GAL_TYPE_UINT8, tablen, __func__,
                                          "flag")
@@ -508,9 +525,19 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
           brightness[*plab] += *f;
           if( *f > 0.0f )  /* For calculatiing the approximate center, */
             {              /* necessary for calculating Sky and STD.   */
-              xy[*plab*xyncols  ] += *f;
-              xy[*plab*xyncols+1] += (double)((f-img)/dsize[1]) * *f;
-              xy[*plab*xyncols+2] += (double)((f-img)%dsize[1]) * *f;
+              pos[*plab*pcols  ] += *f;
+              switch(ndim)
+                {
+                case 2:
+                  pos[*plab*pcols+1] += ( (f-img)/dsize[1] ) * *f;
+                  pos[*plab*pcols+2] += ( (f-img)%dsize[1] ) * *f;
+                  break;
+                case 3:
+                  pos[*plab*pcols+1] += ( (f-img) % dsize[1]*dsize[2] ) * *f;
+                  pos[*plab*pcols+2] += ( (f-img)/dsize[2] ) * *f;
+                  pos[*plab*pcols+3] += ( (f-img)%dsize[2] ) * *f;
+                  break;
+                }
             }
         }
 
@@ -526,8 +553,8 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
     size_t i;
     for(i=1;i<num+1;++i)
       printf("%zu (%u): %-5zu %-13.3f %-13.3f %-13.3f %-13.3f\n", i, flag[i],
-             area[i], brightness[i], xy[i*xyncols], xy[i*xyncols+1],
-             xy[i*xyncols+2]);
+             area[i], brightness[i], pos[i*pcols], pos[i*pcols+1],
+             pos[i*pcols+2]);
   }
   */
 
@@ -559,11 +586,11 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
   for(i=1;i<tablen;++i)
     {
       ave=brightness[i]/area[i];
-      if( area[i]>p->detsnminarea && ave>0.0f && xy[i*xyncols]>0.0f )
+      if( area[i]>p->detsnminarea && ave>0.0f && pos[i*pcols]>0.0f )
         {
           /* Get the flux weighted center coordinates. */
-          coord[0]=GAL_DIMENSION_FLT_TO_INT( xy[i*xyncols+1]/xy[i*xyncols] );
-          coord[1]=GAL_DIMENSION_FLT_TO_INT( xy[i*xyncols+2]/xy[i*xyncols] );
+          for(j=0;j<ndim;++j)
+            coord[j]=GAL_DIMENSION_FLT_TO_INT(pos[i*pcols+j+1]/pos[i*pcols]);
 
           /* Calculate the Sky and Sky standard deviation on this tile. */
           ave -= ((float *)(p->sky->array))[
@@ -611,7 +638,7 @@ detection_sn(struct noisechiselparams *p, gal_data_t 
*worklab, size_t num,
 
 
   /* Clean up and return. */
-  free(xy);
+  free(pos);
   free(area);
   free(coord);
   free(brightness);
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index b6ba506..71b808c 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -62,7 +62,7 @@ noisechisel_convolve(struct noisechiselparams *p)
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
   p->conv=gal_convolve_spatial(tl->tiles, p->kernel, p->cp.numthreads,
                                1, tl->workoverch);
-  gal_checkset_allocate_copy("INPUT_CONVOLVED", &p->conv->name);
+  gal_checkset_allocate_copy("INPUT-CONVOLVED", &p->conv->name);
 
 
   if(!p->cp.quiet) gal_timing_report(&t1, "Convolved with kernel.", 1);
@@ -71,12 +71,6 @@ noisechisel_convolve(struct noisechiselparams *p)
       gal_fits_img_write(p->input, p->detectionname, NULL, PROGRAM_STRING);
       gal_fits_img_write(p->conv, p->detectionname, NULL, PROGRAM_STRING);
     }
-
-  if(p->conv->ndim==3)
-    {
-      printf("\n... end of %s ...\n", __func__);
-      exit(0);
-    }
 }
 
 
@@ -295,6 +289,14 @@ noisechisel(struct noisechiselparams *p)
      images. */
   noisechisel_find_sky_subtract(p);
 
+  /***********************************************/
+  if(p->conv->ndim==3)
+    {
+      printf("\n... end of %s ...\n", __func__);
+      exit(0);
+    }
+  /***********************************************/
+
   /* If the user only wanted detection, ignore the segmentation steps. */
   if(p->onlydetection==0)
     {
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 7f369e4..efb5961 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -155,7 +155,11 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
   /* When the check image has the same resolution as the input, write the
      binary array as a reference to help in the comparison. */
   if(checkname && !tl->oneelempertile)
-    gal_fits_img_write(p->binary, checkname, NULL, PROGRAM_STRING);
+    {
+      p->binary->name="DETECTED";
+      gal_fits_img_write(p->binary, checkname, NULL, PROGRAM_STRING);
+      p->binary->name=NULL;
+    }
 
 
   /* Allocate space for the mean and standard deviation. */
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 7c53eee..d35b3ff 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -221,13 +221,17 @@ parse_opt(int key, char *arg, struct argp_state *state)
 static void
 ui_read_check_only_options(struct noisechiselparams *p)
 {
-  /* Make sure the connectivity is defined. */
-  if(p->erodengb!=4 && p->erodengb!=8)
-    error(EXIT_FAILURE, 0, "%zu not acceptable for `--erodengb'. It must "
-          "be 4 or 8 (specifying the type of connectivity)", p->erodengb);
-  if(p->openingngb!=4 && p->openingngb!=8)
-    error(EXIT_FAILURE, 0, "%zu not acceptable for `--openingngb'. It must "
-          "be 4 or 8 (specifying the type of connectivity)", p->openingngb);
+  /* A general check on the neighbor connectivity values. */
+  if(p->erodengb!=4 && p->erodengb!=8 && p->erodengb!=6 && p->erodengb!=18
+     && p->erodengb!=26)
+    error(EXIT_FAILURE, 0, "%zu is not an acceptable value for "
+          "`--erodengb'. Acceptable values are 4 or 8 (for 2D inputs) and "
+          "6, 18, or 26 (for 3D inputs)", p->erodengb);
+  if(p->openingngb!=4 && p->openingngb!=8 && p->openingngb!=6
+     && p->openingngb!=18 && p->openingngb!=26)
+    error(EXIT_FAILURE, 0, "%zu is not an acceptable value for "
+          "`--openingngb'. Acceptable values are 4 or 8 (for 2D inputs) and "
+          "6, 18, or 26 (for 3D inputs)", p->openingngb);
 
   /* Make sure that the no-erode-quantile is not smaller or equal to
      qthresh. */
@@ -599,12 +603,12 @@ ui_prepare_tiles(struct noisechiselparams *p)
 
 
 
+/* Read the input image and check the dimensions and neighbors. */
 static void
-ui_preparations(struct noisechiselparams *p)
+ui_preparations_read_input(struct noisechiselparams *p)
 {
-  /* Prepare the names of the outputs. */
-  ui_set_output_names(p);
-
+  size_t value;
+  char *option_name=NULL, *good_values=NULL;
 
   /* Read the input as a single precision floating point dataset. */
   p->input = gal_fits_img_read_to_type(p->inputname, p->cp.hdu,
@@ -615,6 +619,69 @@ ui_preparations(struct noisechiselparams *p)
     gal_checkset_allocate_copy("INPUT", &p->input->name);
 
 
+  /* Check dimensionality and neighbors. */
+  switch(p->input->ndim)
+    {
+    case 2:
+      if(p->erodengb!=4 && p->erodengb!=8)
+        {
+          value=p->erodengb;
+          good_values="4 or 8";
+          option_name="erodengb";
+        }
+      else if(p->openingngb!=4 && p->openingngb!=8)
+        {
+          value=p->openingngb;
+          good_values="4 or 8";
+          option_name="openingngb";
+        }
+      break;
+
+    case 3:
+      if(p->erodengb!=6 && p->erodengb!=18 && p->erodengb!=26)
+        {
+          value=p->erodengb;
+          good_values="6, 18, or 26";
+          option_name="erodengb";
+        }
+      else if(p->openingngb!=6 && p->openingngb!=18 && p->openingngb!=26)
+        {
+          value=p->openingngb;
+          good_values="6, 18, or 26";
+          option_name="openingngb";
+        }
+      break;
+
+    default:
+      error(EXIT_FAILURE, 0, "%s (hdu %s) is a %zu-dimensional dataset "
+            "which is currently not supported", p->inputname, p->cp.hdu,
+            p->input->ndim);
+    }
+
+
+  /* Abort with an error message if necessary. */
+  if(option_name)
+    error(EXIT_FAILURE, 0, "%zu not acceptable for `--%s' for the "
+          "input %zu-dimensional dataset in `%s' (hdu %s). It must be %s "
+          "(specifying the type of connectivity)", value, option_name,
+          p->input->ndim, p->inputname, p->cp.hdu, good_values);
+}
+
+
+
+
+
+static void
+ui_preparations(struct noisechiselparams *p)
+{
+  /* Prepare the names of the outputs. */
+  ui_set_output_names(p);
+
+
+  /* Read the input dataset and check the dimensions. */
+  ui_preparations_read_input(p);
+
+
   /* Check for blank values to help later processing. AFTERWARDS, set the
      USE_ZERO flag, so the zero-bit (if the input doesn't have any blank
      value) will be meaningful. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index bcd7dcf..52bbb43 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -5624,14 +5624,14 @@ also contain non-numerical columns.
 @node Gnuastro text table format, Selecting table columns, Recognized table 
formats, Tables
 @subsection Gnuastro text table format
 
-Plain text files are most generic, portable, and easiest way to (manually)
-create, (visually) inspect, or (manually) edit a table. In this format, the
-ending of a row is defined by the new-line character (a line on a text
-editor). So when you view it on a text editor, every row will occupy one
-line. The delimiters (or characters separating the columns) are white space
-characters (space, horizontal tab, vertical tab) and a comma (@key{,}). The
-only further requirement is that all rows/lines must have the same number
-of columns.
+Plain text files are the most generic, portable, and easiest way to
+(manually) create, (visually) inspect, or (manually) edit a table. In this
+format, the ending of a row is defined by the new-line character (a line on
+a text editor). So when you view it on a text editor, every row will occupy
+one line. The delimiters (or characters separating the columns) are white
+space characters (space, horizontal tab, vertical tab) and a comma
+(@key{,}). The only further requirement is that all rows/lines must have
+the same number of columns.
 
 The columns don't have to be exactly under each other and the rows can be
 arbitrarily long with different lengths. For example the following contents
@@ -5644,13 +5644,14 @@ types}).
 2 , 4.454        792     72.98348e7
 @end example
 
-However, the example above has no other information about the columns (its
-raw data, with no meta-data). To use this table, you have to remember what
-each column was. Also, when you want to select columns, you have to count
-their position within the table. This can become frustrating and prone to
-bad errors (getting the columns wrong) especially as the number of columns
-increase. It is also bad for sending to a colleague, because they will find
-it hard to remember/use the columns properly.
+However, the example above has no other information about the columns (it
+is just raw data, with no meta-data). To use this table, you have to
+remember what the numbers in each column represent. Also, when you want to
+select columns, you have to count their position within the table. This can
+become frustrating and prone to bad errors (getting the columns wrong)
+especially as the number of columns increase. It is also bad for sending to
+a colleague, because they will find it hard to remember/use the columns
+properly.
 
 To solve these problems in Gnuastro's programs/libraries you aren't limited
 to using the column's number, see @ref{Selecting table columns}. If the
@@ -5720,24 +5721,24 @@ space characters will be stripped from all of the 
elements. For example in
 this line:
 
 @example
-# Column 5:  column name   [km/s,    f,-99] Redshift as speed
+# Column 5:  column name   [km/s,    f32,-99] Redshift as speed
 @end example
 
-The @code{NAME} field will be `@code{column name}', or @code{TYPE} will be
-`@code{f}'. Note how all the white space characters before and after
-strings are not used, but those in the middle remained. Also, white space
-characters aren't mandatory, so in the example above @code{BLANK} will be
-`@code{-99}'.
+The @code{NAME} field will be `@code{column name}' and the @code{TYPE}
+field will be `@code{f32}'. Note how all the white space characters before
+and after strings are not used, but those in the middle remained. Also,
+white space characters aren't mandatory. Hence, in the example above, the
+@code{BLANK} field will be given the value of `@code{-99}'.
 
 Except for the column number (@code{N}), the rest of the fields are
-optional and the column information comments don't have to be in order. In
-other words, the information for column @mymath{N+m} (@mymath{m>0}) can be
-given before column @mymath{N}. Also, you don't have to specify information
-for all columns. Those columns that don't have this information will be
-interpreted with the default settings (like the case above: values are
-double precision floating point, and the column has no name, unit, or
-comment). So these lines are all acceptable for any table (the first one,
-with nothing but the column number is redundant):
+optional. Alao, the column information comments don't have to be in
+order. In other words, the information for column @mymath{N+m}
+(@mymath{m>0}) can be given in a line before column @mymath{N}. Also, you
+don't have to specify information for all columns. Those columns that don't
+have this information will be interpreted with the default settings (like
+the case above: values are double precision floating point, and the column
+has no name, unit, or comment). So these lines are all acceptable for any
+table (the first one, with nothing but the column number is redundant):
 
 @example
 # Column 5:
@@ -12535,26 +12536,31 @@ output will have the same pixel size as the input, 
but with the
 The number of erosions to apply to the binary thresholded image. Erosion is
 simply the process of flipping (from 1 to 0) any of the foreground pixels
 that neighbor a background pixel. In a 2D image, there are two kinds of
-neighbors, 4-connected and 8-connected neighbors. You can specify which
-type of neighbors should be used for erosion with the @option{--erodengb}
-option, see below.
-
-Erosion has the effect of shrinking the foreground pixels. To put it
-another way, it expands the holes. This is a founding principle in
-NoiseChisel: it exploits the fact that with very low thresholds, the
-holes in the very low surface brightness regions of an image will be
-smaller than regions that have no signal. Therefore by expanding those
-holes, we are able to separate the regions harboring signal.
+neighbors, 4-connected and 8-connected neighbors. In a 3D dataset, there
+are three: 6-connected, 18-connected, and 26-connected. You can specify
+which class of neighbors should be used for erosion with the
+@option{--erodengb} option, see below.
+
+Erosion has the effect of shrinking the foreground (above threshold)
+pixels. To put it another way, it expands the holes. This is a founding
+principle in NoiseChisel: it exploits the fact that with very low
+thresholds, the holes in the very low surface brightness regions of an
+image will be smaller than regions that have no signal. Therefore by
+expanding those holes, we are able to separate the regions harboring
+signal.
 
 @item --erodengb=INT
 The type of neighborhood (structuring element) used in erosion, see
-@option{--erode} for an explanation on erosion. Only two integer values are
-acceptable: 4 or 8. In 4-connectivity, the neighbors of a pixel are defined
-as the four pixels on the top, bottom, right and left of a pixel that share
-an edge with it. The 8-connected neighbors on the other hand include the
+@option{--erode} for an explanation on erosion. If the input is a 2D image,
+only two integer values are acceptable: 4 or 8. For a 3D input datacube,
+the acceptable values are: 6, 18 and 26.
+
+In 2D 4-connectivity, the neighbors of a pixel are defined as the four
+pixels on the top, bottom, right and left of a pixel that share an edge
+with it. The 8-connected neighbors on the other hand include the
 4-connected neighbors along with the other 4 pixels that share a corner
 with this pixel. See Figure 6 (a) and (b) in Akhlaghi and Ichikawa (2015)
-for a demonstration.
+for a demonstration. A similar argument applies to 3D datacubes.
 
 @item --noerodequant
 Pure erosion is going to carve off sharp and small objects completely out
@@ -12588,8 +12594,8 @@ separate foreground `islands' (detections) thereby 
completely isolating
 them. Once opening is complete, we have @emph{initial} detections.
 
 @item --openingngb=INT
-The structuring element used for opening, see @option{--erodengb} for more
-information about a structuring element.
+The class of neighbors used for opening, see @option{--erodengb} for more
+information on acceptable values.
 
 @item -s FLT,FLT
 @itemx --sigmaclip=FLT,FLT
diff --git a/lib/binary.c b/lib/binary.c
index a1adfaf..8d7fe21 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -142,7 +142,7 @@ binary_erode_dilate_2d_8con(gal_data_t *input, unsigned 
char dilate0_erode1)
 
   /* Set the foreground and background values: */
   if(dilate0_erode1==0) {f=1; b=0;}
-  else         {f=0; b=1;}
+  else                  {f=0; b=1;}
 
   /* Check the 4 corners: */
   if(byt[0]==b && (byt[1]==f
@@ -220,6 +220,44 @@ binary_erode_dilate_2d_8con(gal_data_t *input, unsigned 
char dilate0_erode1)
 
 
 
+/* This is a general erosion and dilation function. It is less efficient
+   than the more specialized cases above. */
+static void
+binary_erode_dilate_general(gal_data_t *input, unsigned char dilate0_erode1,
+                            int connectivity)
+{
+  uint8_t f, b, *pt, *fpt, *byt=input->array;
+  size_t i, *dinc=gal_dimension_increment(input->ndim, input->dsize);
+
+  /* Do a sanity check: */
+  if(dilate0_erode1!=1 && dilate0_erode1!=0)
+    error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can fix "
+          "this problem. The value to dilate0_erode1 is %u while it should "
+          "be 0 or 1", __func__, PACKAGE_BUGREPORT, dilate0_erode1);
+
+  /* Set the foreground and background values: */
+  if(dilate0_erode1==0) {f=1; b=0;}
+  else                  {f=0; b=1;}
+
+  /* Go over the neighbors of each pixel. */
+  for(i=0;i<input->size;++i)
+    if(byt[i]==b)
+      GAL_DIMENSION_NEIGHBOR_OP(i, input->ndim, input->dsize, connectivity,
+                                dinc,{
+                                  if(byt[i]!=GAL_BINARY_TMP_VALUE
+                                     && byt[nind]==f)
+                                    byt[i]=GAL_BINARY_TMP_VALUE;
+                                });
+
+  /* Set all the changed pixels to the proper values: */
+  fpt=(pt=byt)+input->size;
+  do *pt = *pt==GAL_BINARY_TMP_VALUE ? f : *pt; while(++pt<fpt);
+}
+
+
+
+
+
 /* Erode a binary dataset any number of times. If `inplace' is given a
    value of `1', then do the erosion within the already allocated space,
    otherwise, allocate a new array and save the result into that.
@@ -245,8 +283,8 @@ binary_erode_dilate(gal_data_t *input, size_t num, int 
connectivity,
 
   /* Set the dataset to work on. */
   binary = ( (inplace && input->type==GAL_TYPE_UINT8)
-             ? input :
-             gal_data_copy_to_new_type(input, GAL_TYPE_UINT8) );
+             ? input
+             : gal_data_copy_to_new_type(input, GAL_TYPE_UINT8) );
 
   /* Go over every element and do the erosion. */
   switch(binary->ndim)
@@ -263,6 +301,11 @@ binary_erode_dilate(gal_data_t *input, size_t num, int 
connectivity,
           }
       break;
 
+    case 3:
+      for(counter=0;counter<num;++counter)
+        binary_erode_dilate_general(binary, d0e1, connectivity);
+      break;
+
     default:
       error(EXIT_FAILURE, 0, "%s: currently doesn't work on %zu "
             "dimensional datasets", __func__, binary->ndim);