[gnuastro-commits] master fca5e80 117/125: First re-implementation of NoiseChisel's complete detection

[Mohammad Akhlaghi]

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

    First re-implementation of NoiseChisel's complete detection
    
    NoiseChisel now identifies false detections, removes them and dilates the
    final set of initial detections. For now, this whole process is done on one
    thread, but this may be changed in the next commits. So until now, there
    has been no need for large tile sizes.
    
    Several necessary library functions were also created, and slightly
    modified to complete this job, in particular these two library functions
    were added: `gal_tile_full_id_from_coord' (which will find the tile ID from
    a coordinate), and `gal_binary_fill_holes' (which identifies the holes in a
    binary image and removes them). The full two-layer tessellation structure
    now also saves the size of the first tile (currently for use in
    `gal_tile_full_id_from_coord').
---
 bin/noisechisel/args.h             |  18 +-
 bin/noisechisel/detection.c        | 506 ++++++++++++++++++++++++++++++++++---
 bin/noisechisel/detection.h        |   2 +-
 bin/noisechisel/main.h             |   9 +-
 bin/noisechisel/noisechisel.c      |   2 +-
 bin/noisechisel/threshold.c        |  96 ++++++-
 bin/noisechisel/threshold.h        |   4 +
 bin/noisechisel/ui.c               |  25 +-
 lib/binary.c                       | 121 ++++++++-
 lib/gnuastro-internal/commonopts.h |   2 +-
 lib/gnuastro/dimension.h           |   3 +
 lib/gnuastro/tile.h                |  15 +-
 lib/statistics.c                   |   5 +-
 lib/tile.c                         |  58 ++++-
 14 files changed, 803 insertions(+), 63 deletions(-)

diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 9d236b0..ee5f503 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -100,6 +100,22 @@ struct argp_option program_options[] =
     },
 
 
+    /* Tessellation.
+    {
+      "largetilesize",
+      ARGS_OPTION_KEY_LARGETILESIZE,
+      "INT[,INT]",
+      0,
+      "Regular tile size on dim.s (FITS order).",
+      GAL_OPTIONS_GROUP_TESSELLATION,
+      &cp->tl.tilesize,
+      GAL_TYPE_SIZE_T,
+      GAL_OPTIONS_RANGE_GT_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      gal_options_parse_sizes_reverse
+    },
+    */
 
 
     /* Output options. */
@@ -342,7 +358,7 @@ struct argp_option program_options[] =
       0,
       "Quantile in pseudo-det. to define true.",
       ARGS_GROUP_DETECTION,
-      &p->checkdetsn,
+      &p->detquant,
       GAL_TYPE_FLOAT32,
       GAL_OPTIONS_RANGE_GT_0_LT_1,
       GAL_OPTIONS_MANDATORY,
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 827d30e..d67fd8b 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -23,10 +23,14 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <config.h>
 
 #include <stdio.h>
+#include <errno.h>
+#include <error.h>
 #include <stdlib.h>
 
 #include <gnuastro/fits.h>
 #include <gnuastro/binary.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
 
 #include <gnuastro-internal/timing.h>
 
@@ -98,19 +102,13 @@ detection_initial(struct noisechiselparams *p)
               p->opening, p->openingngb==4 ? "4" : "8");
       gal_timing_report(&t1, msg, 2);
     }
-  if(p->detectionname)
-    {
-      p->binary->name="OPENED";
-      gal_fits_img_write(p->binary, p->detectionname, NULL, PROGRAM_STRING);
-      p->binary->name=NULL;
-    }
 
 
   /* Label the connected components. */
   p->numobjects=gal_binary_connected_components(p->binary, &p->olabel, 1);
   if(p->detectionname)
     {
-      p->olabel->name="LABELED";
+      p->olabel->name="OPENED-LABELED";
       gal_fits_img_write(p->olabel, p->detectionname, NULL, PROGRAM_STRING);
       p->olabel->name=NULL;
     }
@@ -145,14 +143,317 @@ detection_initial(struct noisechiselparams *p)
 
 
 /****************************************************************
- ************     Pseudo detection S/N threshold     ************
+ ************            Pseudo detections           ************
  ****************************************************************/
+/* Set all the pixels we don't need to Nan. */
+static void
+detection_pseudo_sky_or_det(struct noisechiselparams *p, uint8_t *w, int s0d1)
+{
+  uint32_t *l=p->olabel->array;
+  uint8_t *b=p->binary->array, *bf=b+p->binary->size;
+
+  if(s0d1)
+    /* Set all sky regions (label equal to zero) to zero. */
+    do *w++ = *l++ ? *b : 0; while(++b<bf);
+  else
+    /* Set all detected pixels (label larger than zero) to blank. */
+    do *w++ = *l++ ? GAL_BLANK_UINT8 : *b; while(++b<bf);
+}
+
+
+
+
+
 /* We have the thresholded image (with blank values for regions that should
    not be used). Find the pseudo-detections in those regions. */
-void
-detection_find_pseudos(struct noisechiselparams *p, gal_data_t workbin)
+static size_t
+detection_pseudo_find(struct noisechiselparams *p, gal_data_t *workbin,
+                      gal_data_t *worklab, int s0d1)
+{
+  uint8_t *b, *bf;
+
+
+  /* Set all the initial detected pixels to blank values. */
+  detection_pseudo_sky_or_det(p, workbin->array, s0d1);
+  if(p->detectionname)
+    {
+      workbin->name="DTHRESH-ON-SKY";
+      gal_fits_img_write(workbin, p->detectionname, NULL, PROGRAM_STRING);
+      workbin->name=NULL;
+    }
+
+
+  /* Fill the four-connected bounded holes. */
+  gal_binary_fill_holes(workbin);
+  if(p->detectionname)
+    {
+      workbin->name="HOLES-FILLED";
+      gal_fits_img_write(workbin, p->detectionname, NULL, PROGRAM_STRING);
+      workbin->name=NULL;
+    }
+
+
+  /* Open the image. */
+  gal_binary_open(workbin, 1, 4, 1);
+  if(p->detectionname)
+    {
+      workbin->name="OPENED";
+      gal_fits_img_write(workbin, p->detectionname, NULL, PROGRAM_STRING);
+      workbin->name=NULL;
+    }
+
+
+  /* Label all regions, but first, deal with the blank pixels in the
+     `workbin' dataset when working on the Sky. Recall that the blank
+     pixels are the other domain (detections if working on Sky and sky if
+     working on detections). On the Sky image, blank should be set to 1
+     (because we want the detected objects to have the same labels as the
+     pseudo-detections that cover them). This will allow us to later remove
+     these pseudo-detections. */
+  if(s0d1==0)
+    {
+      bf=(b=workbin->array)+workbin->size;
+      do if(*b==GAL_BLANK_UINT8) *b = !s0d1; while(++b<bf);
+    }
+  return gal_binary_connected_components(workbin, &worklab, 1);
+}
+
+
+
+
+
+#define PSN_EXTNAME "PSEUDOS-FOR-SN"
+static gal_data_t *
+detection_pseudo_sn(struct noisechiselparams *p, gal_data_t *worklab,
+                         size_t num, int s0d1)
+{
+  float *snarr;
+  uint8_t *flag;
+  size_t tablen=num+1;
+  gal_data_t *sn, *snind;
+  uint32_t *plabend, *indarr;
+  double ave, err, *xy, *brightness;
+  struct gal_linkedlist_stll *comments=NULL;
+  size_t ind, ndim=p->input->ndim, xyncols=1+ndim;
+  size_t i, *area, counter=0, *dsize=p->input->dsize;
+  size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+  float *img=p->input->array, *f=p->input->array, *ff=f+p->input->size;
+  uint32_t *plab = worklab->array, *dlab = s0d1 ? NULL : p->olabel->array;
+
+
+  /* Sanity check. */
+  if(p->input->type!=GAL_TYPE_FLOAT32)
+    error(EXIT_FAILURE, 0, "the input dataset to `detection_pseudo_sn' "
+          "must be float32 type, it is %s",
+          gal_type_to_string(p->input->type, 1));
+  if(!isnan(GAL_BLANK_FLOAT32))
+    error(EXIT_FAILURE, 0, "currently `detection_pseudo_sn' only "
+          "recognizes a NaN value for blank floating point data types, the "
+          "blank value is defined to be %f", GAL_BLANK_FLOAT32);
+  if(ndim!=2)
+    error(EXIT_FAILURE, 0, "currently `detection_pseudo_sn' only "
+          "works on 2D datasets, your input is %zu dimensions", ndim);
+
+
+  /* Allocate all the necessary arrays, note that since we want to put each
+     object's information into the same index, the number of allocated
+     spaces has to be `tablen=num+1'. */
+  area       = gal_data_calloc_array(GAL_TYPE_SIZE_T,  tablen          );
+  brightness = gal_data_calloc_array(GAL_TYPE_FLOAT64, tablen          );
+  xy         = gal_data_calloc_array(GAL_TYPE_FLOAT64, xyncols*tablen  );
+  flag       = s0d1==0 ? gal_data_calloc_array(GAL_TYPE_UINT8, tablen) : NULL;
+  sn         = gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &tablen, NULL, 1,
+                              p->cp.minmapsize, "SIGNAL-TO-NOISE", "ratio",
+                              NULL);
+  snind      = ( p->checkdetsn==0 ? NULL
+                 : gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &tablen, NULL, 1,
+                                  p->cp.minmapsize, "LABEL", "counter",
+                                  NULL) );
+
+
+  /* Go over all the pixels and get the necessary information. */
+  do
+    {
+      /* All this work os only necessary when we are actually on a
+         pseudo-detection label? */
+      if(*plab)
+        {
+          /* For Sky pseudo-detections we'll start to see if it has already
+             been determined that the object lies over a detected object or
+             not. If it does, then just ignore it. */
+          if(s0d1==0)
+            {
+              if( flag[*plab] ) { ++plab; ++dlab; continue; }
+              else if(*dlab)    /* We are on a detection. */
+                { flag[*plab]=1; area[*plab]=0; ++plab; ++dlab; continue; }
+            }
+
+          /* If we are on a blank pixel, ignore this pixel. */
+          if( isnan(*f) ) { ++plab; if(s0d1==0) ++dlab; continue; }
+
+          /* Save all the necessary values. */
+          ++area[*plab];
+          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;
+            }
+        }
+
+      /* Increment the other two labels. */
+      ++plab;
+      if(s0d1==0) ++dlab;
+    }
+  while(++f<ff);
+
+  /* A small sanity check.
+  {
+    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]);
+  }
+  */
+
+
+  /* If the user wants to see the steps (on the background), remove all the
+     pseudo-detections that will not be used in the final quantile
+     calcluation. */
+  if(p->detectionname)
+    {
+      plabend = (plab=worklab->array) + worklab->size;
+      do
+        if( *plab && ( area[*plab]<p->detsnminarea || brightness[*plab]<0) )
+          *plab=0;
+      while(++plab<plabend);
+      worklab->name=PSN_EXTNAME;
+      gal_fits_img_write(worklab, p->detectionname, NULL,
+                         PROGRAM_STRING);
+      worklab->name=NULL;
+    }
+
+
+  /* calculate the signal to noise for successful detections: */
+  snarr=sn->array;
+  if(snind) indarr=snind->array;
+  if(s0d1) { snarr[0]=NAN; if(snind) indarr[0]=GAL_BLANK_UINT32; }
+  for(i=1;i<num+1;++i)
+    {
+      ave=brightness[i]/area[i];
+      if( area[i]>p->detsnminarea && ave>0.0f && xy[i*xyncols]>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] );
+
+          /* Calculate the Sky and Sky standard deviation on this tile. */
+          ave -= ((float *)(p->sky->array))[
+                         gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
+          err  = ((float *)(p->std->array))[
+                         gal_tile_full_id_from_coord(&p->cp.tl, coord) ];
+
+          /* If the image was already sky subtracted, the second power of
+             the error needs to be doubled. */
+          err *= p->skysubtracted ? err : 2.0f*err;
+
+          /* Correct the index in the sn to store the Signal to noise
+             ratio. When we are dealing with the noise, we only want the
+             non-zero signal to noise values, so we will just use a
+             counter. But for initial detections, it is very important that
+             their Signal to noise ratio be placed in the same index as
+             their label. */
+          ind = s0d1 ? i : counter++;
+          if(snind) indarr[ind]=i;
+          snarr[ind] = ( sqrt( (float)(area[i])/p->cpscorr )
+                         * ave / sqrt(ave+err) );
+        }
+      else
+        /* In detection pseudo-detections, order matters, so we will set
+           all non-usable values to blank. */
+        if(s0d1)
+          {
+            snarr[i]=NAN;
+            if(snind) indarr[i]=GAL_BLANK_UINT32;;
+          }
+    }
+
+
+  /* If we are in Sky mode, the sizes have to be corrected */
+  if(s0d1==0)
+    {
+      sn->dsize[0]=sn->size=counter;
+      if(snind) snind->dsize[0]=snind->size=counter;
+    }
+
+
+  /* If the user wanted a list of S/N values for all pseudo-detections,
+     save it. */
+  if(snind)
+    {
+      /* Make the comments, then write the table. */
+      gal_linkedlist_add_to_stll(&comments, "See also: `"PSN_EXTNAME"' HDU "
+                                 "of output with `--checkdetection'", 1);
+      gal_linkedlist_add_to_stll(&comments, s0d1 ? "Pseudo-detection S/N "
+                                 "over initially detected regions."
+                                 : "Pseudo-detection S/N over initially "
+                                 "undetected regions.", 1);
+      threshold_write_sn_table(p, sn, snind, ( s0d1 ? p->detsn_d_name
+                                               : p->detsn_s_name ), comments);
+      gal_linkedlist_free_stll(comments, 1);
+    }
+
+
+  /* Clean up and return. */
+  free(xy);
+  free(area);
+  free(brightness);
+  if(flag) free(flag);
+  if(snind) gal_data_free(snind);
+  return sn;
+}
+
+
+
+
+
+static void
+detection_pseudo_remove_low_sn(struct noisechiselparams *p,
+                               gal_data_t *workbin, gal_data_t *worklab,
+                               gal_data_t *sn, float snthresh)
 {
+  size_t i;
+  float *snarr=sn->array;
+  uint8_t *b=workbin->array;
+  uint32_t *l=worklab->array, *lf=l+worklab->size;
+  uint8_t *keep=gal_data_calloc_array(GAL_TYPE_UINT8, sn->size);
+
+  /* Specify the new labels for those that must be kept/changed. Note that
+     when an object didn't have an S/N, its S/N was given a value of NaN
+     (which will fail on any condition), so it acts as if it had an S/N
+     lower than the required value. */
+  for(i=0;i<sn->size;++i)
+    if( snarr[i] > snthresh ) keep[i]=1;
+
+
+  /* Go over the pseudo-detection labels and only keep those that must be
+     kept (using the new labels). */
+  do *b++ = keep[ *l++ ] > 0; while(l<lf);
+
+
+  /* If the user wanted to see the steps. */
+  if(p->detectionname)
+    {
+      workbin->name="TRUE-PSEUDO-DETECTIONS";
+      gal_fits_img_write(workbin, p->detectionname, NULL,
+                         PROGRAM_STRING);
+      workbin->name=NULL;
+    }
 
+  /* Clean up: */
+  free(keep);
 }
 
 
@@ -160,6 +461,61 @@ detection_find_pseudos(struct noisechiselparams *p, 
gal_data_t workbin)
 
 
 
+/* Find and do the necessary work on pseudo-detections. */
+static gal_data_t *
+detection_pseudo_real(struct noisechiselparams *p)
+{
+  char *msg;
+  float snthresh;
+  size_t numpseudo;
+  struct timeval t1;
+  gal_data_t *sn, *quant, *workbin, *worklab;
+
+  /* Allocate the space for the working datasets. */
+  worklab=gal_data_copy(p->olabel);
+  workbin=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
+                         p->input->dsize, p->input->wcs, 0, p->cp.minmapsize,
+                         NULL, NULL, NULL);
+
+  /* Over the Sky: find the pseudo-detections and make the S/N table. */
+  if(!p->cp.quiet) gettimeofday(&t1, NULL);
+  numpseudo=detection_pseudo_find(p, workbin, worklab, 0);
+  sn=detection_pseudo_sn(p, worklab, numpseudo, 0);
+
+
+  /* Get the S/N quantile and report it if we are in non-quiet mode. */
+  quant=gal_statistics_quantile(sn, p->detquant, 1);
+  snthresh = *((float *)(quant->array));
+  if(!p->cp.quiet)
+    {
+      asprintf(&msg, "Pseudo-det S/N: %.2f (%.3f quant of %zu).",
+               snthresh, p->detquant, sn->size);
+      gal_timing_report(&t1, msg, 2);
+    }
+  gal_data_free(sn);
+  gal_data_free(quant);
+
+
+  /* Over the detections: find pseudo-detections and make S/N table. */
+  if(!p->cp.quiet) gettimeofday(&t1, NULL);
+  numpseudo=detection_pseudo_find(p, workbin, worklab, 1);
+  sn=detection_pseudo_sn(p, worklab, numpseudo, 1);
+
+
+  /* Remove the pseudo detections with a low S/N. */
+  detection_pseudo_remove_low_sn(p, workbin, worklab, sn, snthresh);
+
+
+  /* Clean up and return. */
+  gal_data_free(sn);
+  gal_data_free(worklab);
+  return workbin;
+}
+
+
+
+
+
 
 
 
@@ -177,19 +533,69 @@ detection_find_pseudos(struct noisechiselparams *p, 
gal_data_t workbin)
 /****************************************************************
  ************        Removing false detections       ************
  ****************************************************************/
-/*  */
-static void
-detection_keep_sky_or_det(struct noisechiselparams *p, uint8_t *w, int s0d1)
+static size_t
+detection_remove_false_initial(struct noisechiselparams *p,
+                               gal_data_t *workbin)
 {
-  uint32_t *l=p->olabel->array;
-  uint8_t *b=p->binary->array, *bf=b+p->binary->size;
-
-  if(s0d1)
-    /* Set all sky regions (label equal to zero) to blank. */
-    do *w++ = *l++ ? *b : GAL_BLANK_UINT8; while(++b<bf);
+  size_t i;
+  uint8_t *b=workbin->array;
+  uint32_t *l=p->olabel->array, *lf=l+p->olabel->size, curlab=1;
+  uint32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32, p->numobjects+1);
+
+  /* Find the new labels for all the existing labels. Recall that
+     `newlabels' was initialized to zero, so any label that is not given a
+     new label here will be automatically removed. After the first pixel of
+     a label overlaps with dbyt[i], we don't need to check the rest of that
+     object's pixels. At this point, tokeep is only binary: 0 or 1.
+
+     Note that the zeroth element of tokeep can also be non zero, this is
+     because the holes of the labeled regions might be filled during
+     filling the holes, but have not been filled in the original labeled
+     array. They are not important so you can just ignore them. */
+  do
+    {
+      if(*l)
+        {
+          newlabels[ *l ] =
+            newlabels[ *l ]     /* Have we already checked this label?    */
+            ? 1                 /* Yes we have. Just set it to 1.         */
+            : *b;               /* No we haven't, check pseudo-detection. */
+        }
+      ++b;
+    }
+  while(++l<lf);
+  newlabels[0]=0;
+
+
+  /* Now that we know which labels to keep, set the new labels for the
+     detections that must be kept. */
+  for(i=0;i<p->numobjects;++i) if(newlabels[i]) newlabels[i] = curlab++;
+
+
+  /* Replace the byt and olab values with their proper values. If the
+     user doesn't want to dilate, then change the labels in `lab'
+     too. Otherwise, you don't have to worry about the label
+     array. After dilation a new labeling will be done and the whole
+     labeled array will be re-filled.*/
+  b=workbin->array; l=p->olabel->array;
+  if(p->dilate)
+    do
+      {
+        if(*l!=GAL_BLANK_UINT32)
+          *b = newlabels[ *l ] > 0;
+        ++b;
+      }
+    while(++l<lf);
   else
-    /* Set all detected pixels (label larger than zero) to blank. */
-    do *w++ = *l++ ? GAL_BLANK_UINT8 : *b; while(++b<bf);
+    do
+      if(*l!=GAL_BLANK_UINT32)
+        *l = newlabels[ *l ];
+    while(++l<lf);
+
+
+  /* Clean up and return. */
+  free(newlabels);
+  return curlab-1;
 }
 
 
@@ -199,11 +605,20 @@ detection_keep_sky_or_det(struct noisechiselparams *p, 
uint8_t *w, int s0d1)
 /* The initial detection has been done, now we want to remove false
    detections. */
 void
-detection_remove_false(struct noisechiselparams *p)
+detection(struct noisechiselparams *p)
 {
   char *msg;
-  struct timeval t1;
   gal_data_t *workbin;
+  struct timeval t0, t1;
+  size_t num_true_initial;
+
+  /* Report for the user. */
+  if(!p->cp.quiet)
+    {
+      gal_timing_report(NULL, "Starting to find/remove false detections.", 1);
+      gettimeofday(&t0, NULL);
+    }
+
 
   /* Find the Sky and its Standard Deviation from the initial detectios. */
   if(!p->cp.quiet) gettimeofday(&t1, NULL);
@@ -224,24 +639,47 @@ detection_remove_false(struct noisechiselparams *p)
     }
 
 
-  /* Allocate the space for the pseudo-detection threshold. */
-  workbin=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
-                         p->input->dsize, p->input->wcs, 0, p->cp.minmapsize,
-                         NULL, NULL, NULL);
+  /* Find the real pseudo-detections. */
+  workbin=detection_pseudo_real(p);
 
 
-  /* Set all the initial detected pixels to blank values. */
-  detection_keep_sky_or_det(p, workbin->array, 0);
+  /* Only keep the initial detections that overlap with the real
+     pseudo-detections. */
+  num_true_initial=detection_remove_false_initial(p, workbin);
+  if(!p->cp.quiet)
+    {
+      asprintf(&msg, "%zu false initial detections removed.",
+               p->numobjects - num_true_initial);
+      gal_timing_report(&t1, msg, 2);
+      free(msg);
+    }
+
+
+
+  /* If the user asked for dilation, then apply it. */
+  if(p->dilate)
+    {
+      gal_binary_dilate(workbin, p->dilate, 8, 1);
+      p->numobjects = gal_binary_connected_components(workbin, &p->olabel, 8);
+    }
+  else p->numobjects=num_true_initial;
+  if(!p->cp.quiet)
+    {
+      asprintf(&msg, "%zu detections after %zu dilation%s",
+              p->numobjects, p->dilate, p->dilate>1 ? "s." : ".");
+      gal_timing_report(&t0, msg, 1);
+    }
   if(p->detectionname)
     {
-      workbin->name="DTHRESH-ON-SKY";
-      gal_fits_img_write(workbin, p->detectionname, NULL, PROGRAM_STRING);
-      workbin->name=NULL;
+      p->olabel->name="TRUE-INITIAL-DETECTIONS";
+      gal_fits_img_write(p->olabel, p->detectionname, NULL,
+                         PROGRAM_STRING);
+      p->olabel->name=NULL;
     }
 
 
-  /* Clean up. */
-  gal_data_free(workbin);
+
+  /*  */
 
 
   /* If the user wanted to check the threshold and hasn't called
diff --git a/bin/noisechisel/detection.h b/bin/noisechisel/detection.h
index b785f9f..ed2a100 100644
--- a/bin/noisechisel/detection.h
+++ b/bin/noisechisel/detection.h
@@ -27,6 +27,6 @@ void
 detection_initial(struct noisechiselparams *p);
 
 void
-detection_remove_false(struct noisechiselparams *p);
+detection(struct noisechiselparams *p);
 
 #endif
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 5bf608c..36920e1 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -44,6 +44,7 @@ struct noisechiselparams
 {
   /* From command-line */
   struct gal_options_common_params cp; /* Common parameters.              */
+  /*struct gal_tile_two_layer_params ltl;*/ /* Large tessellation.        */
   char             *inputname;  /* Input filename.                        */
   char            *kernelname;  /* Input kernel filename.                 */
   char                  *khdu;  /* Kernel HDU.                            */
@@ -87,7 +88,8 @@ struct noisechiselparams
   /* Internal. */
   char           *qthreshname;  /* Name of Quantile threshold check image.*/
   char            *detskyname;  /* Name of Initial det sky check image.   */
-  char             *detsnname;  /* Name of pseudo-detections S/N values.  */
+  char          *detsn_s_name;  /* Sky pseudo-detections S/N name.        */
+  char          *detsn_d_name;  /* Detection pseudo-detections S/N name.  */
   char         *detectionname;  /* Name of detection steps file.          */
   char               *skyname;  /* Name of Sky estimation steps file.     */
   char           *clumpsnname;  /* Name of Clump S/N values file.         */
@@ -101,6 +103,11 @@ struct noisechiselparams
   gal_data_t             *std;  /* STD of undetected pixels, per tile.    */
   time_t              rawtime;  /* Starting time of the program.          */
 
+  float                medstd;  /* Median STD before interpolation.       */
+  float                minstd;  /* Minimum STD before interpolation.      */
+  float                maxstd;  /* Maximum STD before interpolation.      */
+  float               cpscorr;  /* Counts/second correction.              */
+
   size_t           numobjects;  /* Number of objects detected.            */
   size_t           maxtcontig;  /* Maximum contiguous space for a tile.   */
   size_t            *maxtsize;  /* Maximum size of a single tile.         */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 39ec6ec..2244460 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -69,5 +69,5 @@ noisechisel(struct noisechiselparams *p)
   detection_initial(p);
 
   /* Remove false detections */
-  detection_remove_false(p);
+  detection(p);
 }
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 8b97bfa..e47a794 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -146,6 +146,67 @@ threshold_apply(struct noisechiselparams *p, float 
*value1, float *value2,
 
 
 /**********************************************************************/
+/***************            Write S/N values          *****************/
+/**********************************************************************/
+void
+threshold_write_sn_table(struct noisechiselparams *p, gal_data_t *insn,
+                         gal_data_t *inind, char *filename,
+                         struct gal_linkedlist_stll *comments)
+{
+  gal_data_t *sn, *ind, *cols;
+
+  /* Remove all blank elements. The index and sn values must have the same
+     set of blank elements, but checking on the integer array is faster. */
+  if( gal_blank_present(inind) )
+    {
+      ind=gal_data_copy(inind);
+      sn=gal_data_copy(insn);
+      gal_blank_remove(ind);
+      gal_blank_remove(sn);
+    }
+  else
+    {
+      sn  = insn;
+      ind = inind;
+    }
+
+  /* Set the columns. */
+  cols       = ind;
+  cols->next = sn;
+
+
+  /* Prepare the comments. */
+  gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
+
+
+  /* write the table. */
+  gal_table_write(cols, comments, p->cp.tableformat, filename, 1);
+
+
+  /* Clean up (if necessary). */
+  if(sn!=insn) gal_data_free(sn);
+  if(ind==inind) ind->next=NULL; else gal_data_free(ind);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
 /***************     Interpolation and smoothing      *****************/
 /**********************************************************************/
 /* Interpolate and smooth the values for each tile over the whole image. */
@@ -477,7 +538,7 @@ threshold_mean_std_undetected(void *in_prm)
       bintile->size=tile->size;
       bintile->dsize=tile->dsize;
       bintile->array=gal_tile_block_relative_to_other(tile, p->binary);
-      GAL_TILE_PARSE_OPERATE({if(!*i) numsky++;}, bintile, NULL, 0, 0);
+      GAL_TILE_PARSE_OPERATE({if(!*o) numsky++;}, tile, bintile, 1, 1);
 
       /* Only continue, if the fraction of Sky values are less than the
          requested fraction. */
@@ -485,8 +546,14 @@ threshold_mean_std_undetected(void *in_prm)
         {
           /* Calculate the mean and STD over this tile. */
           s=s2=0.0f;
-          GAL_TILE_PARSE_OPERATE({ if(!*o) {s += *i; s2 += *i * *i;} },
-                                 tile, bintile, 1, 1);
+          GAL_TILE_PARSE_OPERATE(
+                                 {
+                                   if(!*o)
+                                     {
+                                       s  += *i;
+                                       s2 += *i * *i;
+                                     }
+                                 }, tile, bintile, 1, 1);
           darr[0]=s/numsky;
           darr[1]=sqrt( (s2-s*s/numsky)/numsky );
 
@@ -527,6 +594,7 @@ threshold_mean_std_undetected(void *in_prm)
 void
 threshold_sky_and_std(struct noisechiselparams *p)
 {
+  gal_data_t *tmp;
   struct gal_options_common_params *cp=&p->cp;
   struct gal_tile_two_layer_params *tl=&cp->tl;
 
@@ -553,6 +621,28 @@ threshold_sky_and_std(struct noisechiselparams *p)
       gal_tile_full_values_write(p->std, tl, p->detskyname, PROGRAM_STRING);
     }
 
+  /* Get the basic information about the standard deviation
+     distribution. */
+  tmp=gal_statistics_median(p->std, 0);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->medstd, tmp->array, sizeof p->medstd);
+  free(tmp);
+
+  tmp=gal_statistics_minimum(p->std);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->minstd, tmp->array, sizeof p->minstd);
+  free(tmp);
+
+  tmp=gal_statistics_maximum(p->std);
+  tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+  memcpy(&p->maxstd, tmp->array, sizeof p->maxstd);
+  free(tmp);
+
+  /* In case the image is in electrons or counts per second, the standard
+     deviation of the noise will become smaller than unity, so we need to
+     correct it in the S/N calculation. So, we'll calculate the correction
+     factor here. */
+  p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
 
   /* Interpolate and smooth the derived values. */
   threshold_interp_smooth(p, &p->sky, &p->std, p->detskyname);
diff --git a/bin/noisechisel/threshold.h b/bin/noisechisel/threshold.h
index 9c9a1fb..159699f 100644
--- a/bin/noisechisel/threshold.h
+++ b/bin/noisechisel/threshold.h
@@ -36,6 +36,10 @@ enum threshold_type
 void
 threshold_quantile_find_apply(struct noisechiselparams *p);
 
+void
+threshold_write_sn_table(struct noisechiselparams *p, gal_data_t *sntable,
+                         gal_data_t *snind, char *filename,
+                         struct gal_linkedlist_stll *comments);
 
 void
 threshold_sky_and_std(struct noisechiselparams *p);
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 7d8f413..b8a03ec 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -127,7 +127,6 @@ ui_initialize_options(struct noisechiselparams *p,
         {
         case GAL_OPTIONS_KEY_SEARCHIN:
         case GAL_OPTIONS_KEY_IGNORECASE:
-        case GAL_OPTIONS_KEY_TABLEFORMAT:
           cp->coptions[i].flags=OPTION_HIDDEN;
           break;
 
@@ -229,6 +228,16 @@ ui_read_check_only_options(struct noisechiselparams *p)
           "must be larger than the base quantile threshold (`--qthresh', "
           "or `-t'). You have provided %.4f and %.4f for the former and "
           "latter, respectively", p->noerodequant, p->qthresh);
+
+  /* For the options that make tables, the table formation option is
+     mandatory. */
+  if( (p->checkdetsn || p->checkclumpsn) && p->cp.tableformat==0 )
+    error(EXIT_FAILURE, 0, "`--tableformat' is necessary with the "
+          "`--checkdetsn' and `--checkclumpsn' options.\n"
+          "Please see description for `--tableformat' after running the "
+          "following command for more information (use `SPACE' to go down "
+          "the page and `q' to return to the command-line):\n\n"
+          "    $ info gnuastro \"Input Output options\"");
 }
 
 
@@ -326,8 +335,14 @@ ui_set_output_names(struct noisechiselparams *p)
 
   /* Pseudo-detection S/N values. */
   if(p->checkdetsn)
-    p->detsnname=gal_checkset_automatic_output(&p->cp, basename,
-                                               "_detsn.txt");
+    {
+      p->detsn_s_name=gal_checkset_automatic_output(&p->cp, basename,
+                 ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
+                   ? "_detsn_sky.txt" : "_detsn_sky.fits") );
+      p->detsn_d_name=gal_checkset_automatic_output(&p->cp, basename,
+                 ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
+                   ? "_detsn_det.txt" : "_detsn_det.fits") );
+    }
 
   /* Detection steps. */
   if(p->checkdetection)
@@ -438,6 +453,7 @@ void
 ui_preparations(struct noisechiselparams *p)
 {
   gal_data_t *check;
+  /*struct gal_tile_two_layer_params *ltl=&p->ltl;*/
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
   /* Prepare the names of the outputs. */
@@ -633,10 +649,11 @@ ui_free_report(struct noisechiselparams *p, struct 
timeval *t1)
   free(p->maxtsize);
   free(p->cp.output);
   if(p->skyname)          free(p->skyname);
-  if(p->detsnname)        free(p->detsnname);
   if(p->detskyname)       free(p->detskyname);
   if(p->clumpsnname)      free(p->clumpsnname);
   if(p->qthreshname)      free(p->qthreshname);
+  if(p->detsn_s_name)     free(p->detsn_s_name);
+  if(p->detsn_d_name)     free(p->detsn_d_name);
   if(p->detectionname)    free(p->detectionname);
   if(p->segmentationname) free(p->segmentationname);
 
diff --git a/lib/binary.c b/lib/binary.c
index faac02f..bb72446 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -28,7 +28,9 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <stdlib.h>
 #include <string.h>
 
+#include <gnuastro/fits.h>
 #include <gnuastro/tile.h>
+#include <gnuastro/blank.h>
 #include <gnuastro/binary.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/linkedlist.h>
@@ -400,7 +402,7 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
   b=binary->array;
   for(i=0;i<binary->size;++i)
     /* Check if this pixel is already labeled. */
-    if( b[i] && !l[i] )
+    if( b[i] && l[i]==0 )
       {
         /* This is the first pixel of this connected region that we have
            got to. */
@@ -421,7 +423,7 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
             GAL_DIMENSION_NEIGHBOR_OP(p, binary->ndim, binary->dsize,
                                       connectivity, dinc,
               {
-                if( b[ nind ] && !l[ nind ] )
+                if( b[ nind ] && l[ nind ]==0 )
                   {
                     l[ nind ] = curlab;
                     gal_linkedlist_add_to_sll(&Q, nind);
@@ -462,12 +464,127 @@ gal_binary_connected_components(gal_data_t *binary, 
gal_data_t **out,
 /*********************************************************************/
 /*****************            Fill holes          ********************/
 /*********************************************************************/
+/* Make the array that is the inverse of the input byt of fill
+   holes. The inverse array will also be 4 pixels larger in both
+   dimensions. This is because we might also want to fill those holes
+   that are touching the side of the image. One pixel for a pixel that
+   is one pixel away from the image border. Another pixel for those
+   objects that are touching the image border. */
+static gal_data_t *
+binary_make_padded_inverse(gal_data_t *input, gal_data_t **outtile)
+{
+  uint8_t *in;
+  size_t i, startind;
+  gal_data_t *inv, *tile;
+  size_t *dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim);
+  size_t *startcoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, input->ndim);
+
+
+  /* Set the size of the padded inverse image and the coordinates of the
+     start. We want the inverse to be padded on the edges of each dimension
+     by 2 pixels, so each dimension should be padded by 4 pixels. */
+  for(i=0;i<input->ndim;++i)
+    {
+      startcoord[i]=2;
+      dsize[i]=input->dsize[i]+4;
+    }
+
+
+  /* Allocate the inverse dataset and initialize it to 1 (mainly for the
+     edges, the inner region will be set afterwards).
+
+     PADDING MUST BE INITIALIZED WITH 1: This is done so the connected body
+     of 1 valued pixels (after inversion) gets a label of 1 after labeling
+     the connected components and any hole, will get a value larger than
+     1. */
+  inv=gal_data_alloc(NULL, GAL_TYPE_UINT8, input->ndim, dsize, NULL, 0,
+                     input->minmapsize, "INVERSE", "binary", NULL);
+  memset(inv->array, 1, inv->size);
+
+
+  /* Define a tile to fill the central regions of the inverse. */
+  startind=gal_dimension_coord_to_index(input->ndim, inv->dsize, startcoord);
+  tile=gal_data_alloc(gal_data_ptr_increment(inv->array, startind, inv->type),
+                      inv->type, input->ndim, input->dsize, NULL, 0, 0, NULL,
+                      NULL, NULL);
+  *outtile=tile;
+  tile->block=inv;
+
+
+  /* Fill the central regions. */
+  in=input->array;
+  GAL_TILE_PARSE_OPERATE({*i = *in==GAL_BLANK_UINT8 ? *in : !*in; ++in;},
+                         tile, NULL, 0, 0);
+
+
+  /* Clean up and return. */
+  free(dsize);
+  return inv;
+}
+
+
+
+
+
+/* Fill all the holes in an input unsigned char array that are bounded
+   within a 4-connected region.
+
+   The basic method is this:
+
+   1. An inverse image is created:
+
+        * For every pixel in the input that is 1, the inverse is 0.
+
+        * The inverse image has two extra pixels on each edge to
+          ensure that all the inv[i]==1 pixels around the image are
+          touching each other and a diagonal object passing through
+          the image does not cause the inv[i]==1 pixels on the edges
+          of the image to get a different label.
+
+   2. The 8 connected regions in this inverse image are found.
+
+   3. Since we had a 2 pixel padding on the edges of the image, we
+      know for sure that all labeled regions with a label of 1 are
+      actually connected `holes' in the input image.
+
+      Any pixel with a label larger than 1, is therefore a bounded
+      hole that is not 8-connected to the rest of the holes.  */
 void
 gal_binary_fill_holes(gal_data_t *input)
 {
+  uint8_t *in;
+  gal_data_t *inv, *tile, *holelabs=NULL;
+
   /* A small sanity check. */
   if( input->type != GAL_TYPE_UINT8 )
     error(EXIT_FAILURE, 0, "input to `gal_binary_fill_holes' must have "
           "`uint8' type, but its input dataset has `%s' type",
           gal_type_to_string(input->type, 1));
+
+
+  /* Make the inverse image. */
+  inv=binary_make_padded_inverse(input, &tile);
+
+
+  /* Label the 8-connected (connectivity==2) holes. */
+  gal_binary_connected_components(inv, &holelabs, 2);
+
+
+  /* Any pixel with a label larger than 1 is a hole in the input image and
+     we should invert the respective pixel. To do it, we'll use the tile
+     that was defined before, just change its block and array.*/
+  in=input->array;
+  tile->array=gal_tile_block_relative_to_other(tile, holelabs);
+  tile->block=holelabs; /* has to be after correcting `tile->array'. */
+  GAL_TILE_PARSE_OPERATE({
+      *in = *i>1 && *i!=GAL_BLANK_UINT32 ? 1 : *in;
+      ++in;
+    }, tile, NULL, 0, 0);
+
+
+  /* Clean up and return. */
+  tile->array=NULL;
+  gal_data_free(inv);
+  gal_data_free(tile);
+  gal_data_free(holelabs);
 }
diff --git a/lib/gnuastro-internal/commonopts.h 
b/lib/gnuastro-internal/commonopts.h
index e4ac0ad..78ad0d1 100644
--- a/lib/gnuastro-internal/commonopts.h
+++ b/lib/gnuastro-internal/commonopts.h
@@ -315,7 +315,7 @@ struct argp_option gal_commonopts_options[] =
       GAL_OPTIONS_KEY_MINMAPSIZE,
       "INT",
       0,
-      "Minimum no. bytes to map arrays to hdd/ssd.",
+      "Minimum bytes in array to not use ram RAM.",
       GAL_OPTIONS_GROUP_OPERATING_MODE,
       &cp->minmapsize,
       GAL_TYPE_SIZE_T,
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 6bd4b8c..9712bf2 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -65,6 +65,9 @@ gal_dimension_num_neighbors(size_t ndim);
 /************************************************************************/
 /********************          Coordinates         **********************/
 /************************************************************************/
+#define GAL_DIMENSION_FLT_TO_INT(FLT) (FLT)-(int)(FLT) > 0.5f    \
+  ? (int)(FLT)+1 : (int)(FLT);
+
 void
 gal_dimension_add_coords(size_t *c1, size_t *c2, size_t *out, size_t ndim);
 
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index 9b5c829..f198e97 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -109,15 +109,18 @@ struct gal_tile_two_layer_params
   size_t         *numtilesinch; /* Tile no. in each dim. on one channel.  */
   char          *tilecheckname; /* Name of file to check tiles.           */
   size_t          *permutation; /* Tile pos. in memory --> pos. overall.  */
+  size_t           *firsttsize; /* See `gal_tile_full_regular_first'.     */
 
   /* Actual tile and channel data structures. */
   gal_data_t            *tiles; /* Tiles array (also linked with `next'). */
   gal_data_t         *channels; /* Channels array (linked with `next').   */
 };
 
+
 size_t *
 gal_tile_full(gal_data_t *input, size_t *regular,
-              float remainderfrac, gal_data_t **out, size_t multiple);
+              float remainderfrac, gal_data_t **out, size_t multiple,
+              size_t **firsttsize);
 
 void
 gal_tile_full_sanity_check(char *filename, char *hdu, gal_data_t *input,
@@ -140,6 +143,10 @@ gal_tile_full_values_smooth(gal_data_t *tilevalues,
                             struct gal_tile_two_layer_params *tl,
                             size_t width, size_t numthreads);
 
+size_t
+gal_tile_full_id_from_coord(struct gal_tile_two_layer_params *tl,
+                            size_t *coord);
+
 void
 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
 
@@ -191,10 +198,10 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl);
         /* itself or not. */                                            \
         if(hasblank)                                                    \
           {                                                             \
-            if(b==b) do if(*i!=b)  {OP; if(parse_out) ++o;} while(++i<f); \
-            else     do if(*i==*i) {OP; if(parse_out) ++o;} while(++i<f); \
+            if(b==b) do{if(*i!=b) {OP;}if(parse_out) ++o;} while(++i<f); \
+            else     do{if(*i==*i){OP;}if(parse_out) ++o;} while(++i<f); \
           }                                                             \
-        else         do            {OP; if(parse_out) ++o;} while(++i<f); \
+        else         do          {{OP;}if(parse_out) ++o;} while(++i<f); \
                                                                         \
                                                                         \
         /* Set the incrementation. On a fully allocated iblock (when */ \
diff --git a/lib/statistics.c b/lib/statistics.c
index 33f6996..7476e9c 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -263,9 +263,10 @@ gal_data_t *
 gal_statistics_median(gal_data_t *input, int inplace)
 {
   size_t dsize=1;
-  gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
+  gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);;
   gal_data_t *out=gal_data_alloc(NULL, nbs->type, 1, &dsize, NULL, 1, -1,
                                  NULL, NULL, NULL);
+
   /* Write the median. */
   statistics_median_in_sorted_no_blank(nbs, out->array);
 
@@ -1201,7 +1202,7 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int 
inplace)
       else
         {
           noblank=gal_data_copy(contig);   /* We aren't allowed to touch */
-          gal_blank_remove(contig);        /* the input, so make a copy. */
+          gal_blank_remove(noblank);       /* the input, so make a copy. */
         }
     }
   else noblank=contig;
diff --git a/lib/tile.c b/lib/tile.c
index 4fe885d..e72dcf5 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -479,7 +479,7 @@ gal_tile_block_relative_to_other(gal_data_t *tile, 
gal_data_t *other)
    input's length. */
 static void
 gal_tile_full_regular_first(gal_data_t *parent, size_t *regular,
-                            float significance, size_t *first, size_t *last,
+                            float remainderfrac, size_t *first, size_t *last,
                             size_t *tsize)
 {
   size_t i, remainder, *dsize=parent->dsize;;
@@ -493,7 +493,7 @@ gal_tile_full_regular_first(gal_data_t *parent, size_t 
*regular,
       /* Depending on the remainder, set the first tile size and number. */
       if(remainder)
         {
-          if( remainder > significance * regular[i] )
+          if( remainder > remainderfrac * regular[i] )
             {
               first[i]  = ( remainder + regular[i] )/2;
               tsize[i] = dsize[i]/regular[i] + 1 ;
@@ -574,7 +574,8 @@ gal_tile_full_regular_first(gal_data_t *parent, size_t 
*regular,
 */
 size_t *
 gal_tile_full(gal_data_t *input, size_t *regular,
-              float remainderfrac, gal_data_t **out, size_t multiple)
+              float remainderfrac, gal_data_t **out, size_t multiple,
+              size_t **firsttsize)
 {
   size_t i, d, tind, numtiles, *start=NULL;
   gal_data_t *tiles, *block=gal_tile_block(input);
@@ -588,7 +589,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
   /* Set the first tile size and total number of tiles along each
      dimension, then allocate the array of tiles. */
   gal_tile_full_regular_first(input, regular, remainderfrac,
-                             first, last, tsize);
+                              first, last, tsize);
   numtiles=gal_dimension_total_size(input->ndim, tsize);
 
 
@@ -682,9 +683,9 @@ gal_tile_full(gal_data_t *input, size_t *regular,
 
   /* Clean up and return. */
   free(last);
-  free(first);
   free(coord);
   free(tcoord);
+  *firsttsize=first;
   if(start) free(start);
   tsize[input->ndim]=-1; /* So it can be used for another tessellation, */
   return tsize;          /* see `gal_tile_full_sanity_check'.           */
@@ -796,14 +797,14 @@ gal_tile_full_two_layers(gal_data_t *input,
                          struct gal_tile_two_layer_params *tl)
 {
   gal_data_t *t;
-  size_t i, *junk, ndim=tl->ndim=input->ndim;
+  size_t i, *junk, *junk2, ndim=tl->ndim=input->ndim;
 
   /* Initialize.  */
   tl->channels=tl->tiles=NULL;
 
   /* Initialize necessary values and do the channels tessellation. */
   junk = gal_tile_full(input, tl->channelsize, tl->remainderfrac,
-                       &tl->channels, 1);
+                       &tl->channels, 1, &junk2);
   tl->totchannels = gal_dimension_total_size(ndim, tl->numchannels);
   for(i=0;i<ndim;++i)
     if(junk[i]!=tl->numchannels[i])
@@ -811,13 +812,14 @@ gal_tile_full_two_layers(gal_data_t *input,
             "`gal_tile_full_two_layers' don't match in dimension %zu, with "
             "values of %zu and %zu respectively.", ndim-i,
             tl->numchannels[i], junk[i]);
+  free(junk2);
 
   /* Tile each channel. While tiling the first channel, we are also going
      to allocate the space for the other channels. Then pass those pointers
      when we want to fill in each tile of the other channels. */
   tl->numtilesinch = gal_tile_full(tl->channels, tl->tilesize,
                                    tl->remainderfrac, &tl->tiles,
-                                   tl->totchannels);
+                                   tl->totchannels, &tl->firsttsize);
   tl->tottilesinch = gal_dimension_total_size(ndim, tl->numtilesinch);
   for(i=1; i<tl->totchannels; ++i)
     {
@@ -830,8 +832,9 @@ gal_tile_full_two_layers(gal_data_t *input,
       /* Fill in the information for all the tiles in this channel. Note
          that we already have the returned value, so it isn't important.*/
       junk=gal_tile_full(&tl->channels[i], tl->tilesize, tl->remainderfrac,
-                         &t, 1);
+                         &t, 1, &junk2);
       free(junk);
+      free(junk2);
     }
 
   /* Multiply the number of tiles along each dimension OF ONE CHANNEL by
@@ -1065,6 +1068,42 @@ gal_tile_full_values_smooth(gal_data_t *tilevalues,
 
 
 
+size_t
+gal_tile_full_id_from_coord(struct gal_tile_two_layer_params *tl,
+                            size_t *coord)
+{
+  /* This function only works for 10 dimensions. */
+  size_t i, tr, chid, tile[10];
+
+
+  /* Host channel's ID. */
+  for(i=0;i<tl->ndim;++i)
+    tile[i] = tl->totchannels == 1 ? 0 : coord[i] / tl->channelsize[i];
+  chid=gal_dimension_coord_to_index(tl->ndim, tl->numchannels, tile);
+
+
+  /* Find the tile within the channel. */
+  for(i=0;i<tl->ndim;++i)
+    {
+      tr=coord[i] % tl->channelsize[i];
+      if( tl->firsttsize[i] != tl->tilesize[i] )
+        tile[i] = ( tr <= tl->firsttsize[i]
+                    ? 0
+                    : 1 + (tr - tl->firsttsize[i]) / tl->tilesize[i] );
+      else
+        tile[i] = tr / tl->tilesize[i];
+    }
+
+
+  /* Return the tile ID. */
+  return ( chid * tl->tottilesinch
+           + gal_dimension_coord_to_index(tl->ndim, tl->numtilesinch, tile) );
+}
+
+
+
+
+
 /* Clean up the allocated spaces in the parameters. */
 void
 gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl)
@@ -1077,6 +1116,7 @@ gal_tile_full_free_contents(struct 
gal_tile_two_layer_params *tl)
   if(tl->numtilesinch)  free(tl->numtilesinch);
   if(tl->tilecheckname) free(tl->tilecheckname);
   if(tl->permutation)   free(tl->permutation);
+  if(tl->firsttsize)    free(tl->firsttsize);
 
   /* Free the arrays of `gal_data_c' for each tile and channel. */
   if(tl->tiles)    gal_data_array_free(tl->tiles,    tl->tottiles,    0);