[gnuastro-commits] master 8348f1e 119/125: NoiseChisel pseudo-detections identified on threads

[Mohammad Akhlaghi]

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

    NoiseChisel pseudo-detections identified on threads
    
    Before implementing the new NoiseChisel (with `gal_data_t'), it would
    identify pseudo-detections (fill holes and open the result) on multiple
    threads and be much more efficient. But in the new implementation, it done
    everything on the whole image on one thread which was slow. So with this
    commit, it now applies the same strategy as before and will do the job on
    multiple threads. This was partly done in preparation for the segmentation,
    where a similar multi-threaded approach (when asking for a check image) is
    required.
    
    The tiles that the pseudo-detections are found on are larger tiles than the
    standard ones for sky or threshold definition. This was also necessary for
    segmentation (when trying to find the clump S/N value based on the
    undetected regions). In the process, a bug was found in the tessellation
    (which would appear when only one tile was necessary for a dimension), and
    it has been fixed.
---
 bin/noisechisel/Makefile.am                |   4 +-
 bin/noisechisel/args.h                     |  10 +-
 bin/noisechisel/astnoisechisel.conf        |  51 ++++----
 bin/noisechisel/detection.c                | 194 ++++++++++++++++++++++++++---
 bin/noisechisel/main.c                     |   2 +-
 bin/noisechisel/main.h                     |   6 +-
 bin/noisechisel/noisechisel.c              |  13 +-
 bin/noisechisel/{main.c => segmentation.c} |  29 +----
 bin/noisechisel/{main.c => segmentation.h} |  39 +-----
 bin/noisechisel/ui.c                       | 110 ++++++++++++----
 bin/noisechisel/ui.h                       |   5 +-
 lib/tile.c                                 |  57 +++++----
 12 files changed, 350 insertions(+), 170 deletions(-)

diff --git a/bin/noisechisel/Makefile.am b/bin/noisechisel/Makefile.am
index b5fe436..4918b8a 100644
--- a/bin/noisechisel/Makefile.am
+++ b/bin/noisechisel/Makefile.am
@@ -31,10 +31,10 @@ bin_PROGRAMS = astnoisechisel
 astnoisechisel_LDADD = -lgnuastro
 
 astnoisechisel_SOURCES = main.c ui.c detection.c noisechisel.c sky.c     \
-  threshold.c
+  segmentation.c threshold.c
 
 EXTRA_DIST = main.h authors-cite.h args.h ui.h detection.h noisechisel.h \
-  sky.h threshold.h
+  segmentation.h sky.h threshold.h
 
 
 
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 9c713ae..bc20412 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -100,22 +100,22 @@ struct argp_option program_options[] =
     },
 
 
-    /* Tessellation.
+    /* Tessellation. */
     {
       "largetilesize",
       ARGS_OPTION_KEY_LARGETILESIZE,
       "INT[,INT]",
       0,
-      "Regular tile size on dim.s (FITS order).",
+      "Sim. to --tilesize, but for larger tiles.",
       GAL_OPTIONS_GROUP_TESSELLATION,
-      &cp->tl.tilesize,
+      &p->ltl.tilesize,
       GAL_TYPE_SIZE_T,
       GAL_OPTIONS_RANGE_GT_0,
-      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_MANDATORY,
       GAL_OPTIONS_NOT_SET,
       gal_options_parse_sizes_reverse
     },
-    */
+
 
 
     /* Output options. */
diff --git a/bin/noisechisel/astnoisechisel.conf 
b/bin/noisechisel/astnoisechisel.conf
index 90d9556..573e3f5 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -18,33 +18,36 @@
 # without any warranty.
 
 # Input:
- khdu               1
- minbfrac         0.5
- minnumfalse      100
+ khdu                1
+ minbfrac          0.5
+ minnumfalse       100
+
+# Tessellation
+ largetilesize 200,200
 
 # Detection:
- mirrordist       1.5
- modmedqdiff     0.01
- qthresh          0.3
- smoothwidth        3
- erode              2
- erodengb           4
- noerodequant  0.9331
- opening            1
- openingngb         8
- sigmaclip      3,0.2
- dthresh          0.0
- detsnminarea      10
- detquant        0.95
- dilate             3
+ mirrordist        1.5
+ modmedqdiff      0.01
+ qthresh           0.3
+ smoothwidth         3
+ erode               2
+ erodengb            4
+ noerodequant   0.9331
+ opening             1
+ openingngb          8
+ sigmaclip       3,0.2
+ dthresh           0.0
+ detsnminarea       10
+ detquant         0.95
+ dilate              3
 
 # Segmentation
- segsnminarea      15
- keepmaxnearriver   0
- segquant        0.95
- gthresh          0.5
- minriverlength    15
- objbordersn        1
+ segsnminarea       15
+ keepmaxnearriver    0
+ segquant         0.95
+ gthresh           0.5
+ minriverlength     15
+ objbordersn         1
 
 # Operating mode
- continueaftercheck 0
\ No newline at end of file
+ continueaftercheck  0
\ No newline at end of file
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 7df73e5..a04f7c8 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -26,9 +26,11 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 #include <errno.h>
 #include <error.h>
 #include <stdlib.h>
+#include <string.h>
 
 #include <gnuastro/fits.h>
 #include <gnuastro/binary.h>
+#include <gnuastro/threads.h>
 #include <gnuastro/dimension.h>
 #include <gnuastro/statistics.h>
 
@@ -165,6 +167,105 @@ detection_pseudo_sky_or_det(struct noisechiselparams *p, 
uint8_t *w, int s0d1)
 
 
 
+/* Copy the space of this tile into the full/large array. */
+static void
+detection_write_in_large(gal_data_t *tile, gal_data_t *copy)
+{
+  uint8_t *c=copy->array;
+  GAL_TILE_PARSE_OPERATE({*i=*c++;}, tile, NULL, 0, 0);
+}
+
+
+
+
+
+/* Fill the holes and open on multiple threads to find the
+   pseudo-detections. Ideally both these should be done immediately after
+   each other on each large tile, but when the user wants to check the
+   steps, we need to break out of the threads at each step. */
+struct fho_params
+{
+  int                    step;
+  uint8_t          *copyspace;
+  gal_data_t         *workbin;
+  gal_data_t         *worklab;
+  struct noisechiselparams *p;
+};
+
+static void *
+detection_fill_holes_open(void *in_prm)
+{
+  struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+  struct fho_params *fho_prm=(struct fho_params *)(tprm->params);
+  struct noisechiselparams *p=fho_prm->p;
+
+  void *tarray;
+  gal_data_t *tile, *copy, *tblock;
+  size_t i, dsize[]={1,1,1,1,1,1,1,1,1,1}; /* For upto 10-Dimensions! */
+
+
+  /* A temporary data structure to wrap around the copy space. Note that
+     the initially allocated space for this tile is only 1 pixel! */
+  copy=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim, dsize,
+                      NULL, 0, -1, NULL, NULL, NULL);
+  free(copy->array);
+  copy->array=&fho_prm->copyspace[p->maxltcontig*tprm->id];
+
+
+  /* Go over all the tiles given to this thread. */
+  for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+    {
+      /* For easy reading. */
+      tile=&p->ltl.tiles[tprm->indexs[i]];
+
+      /* Change the tile pointers (temporarily). */
+      tarray=tile->array;
+      tblock=tile->block;
+      tile->array=gal_tile_block_relative_to_other(tile, fho_prm->workbin);
+      tile->block=fho_prm->workbin;
+
+      /* Copy the tile into the contiguous patch of memory to work on, but
+         first reset the size element so `gal_data_copy_to_allocated' knows
+         there is enough space. */
+      copy->size=p->maxltcontig;
+      gal_data_copy_to_allocated(tile, copy);
+
+      /* Fill the holes in this tile. */
+      gal_binary_fill_holes(copy);
+      if(fho_prm->step!=2)
+        {
+          detection_write_in_large(tile, copy);
+          tile->array=tarray;
+          tile->block=tblock;
+          continue;
+        }
+
+      /* Open all the regions. */
+      gal_binary_open(copy, 1, 4, 1);
+
+      /* Write the copied region back into the large input and AFTERWARDS,
+         correct the tile's pointers, the pointers must not be corrected
+         before writing the copy. */
+      detection_write_in_large(tile, copy);
+      tile->array=tarray;
+      tile->block=tblock;
+    }
+
+
+  /* Clean up. */
+  copy->array=NULL;
+  gal_data_free(copy);
+
+
+  /* Wait until all the threads finish and return. */
+  if(tprm->b) pthread_barrier_wait(tprm->b);
+  return NULL;
+}
+
+
+
+
+
 /* We have the thresholded image (with blank values for regions that should
    not be used). Find the pseudo-detections in those regions. */
 static size_t
@@ -172,6 +273,8 @@ detection_pseudo_find(struct noisechiselparams *p, 
gal_data_t *workbin,
                       gal_data_t *worklab, int s0d1)
 {
   uint8_t *b, *bf;
+  gal_data_t *bin;
+  struct fho_params fho_prm={0, NULL, workbin, worklab, p};
 
 
   /* Set all the initial detected pixels to blank values. */
@@ -184,33 +287,84 @@ detection_pseudo_find(struct noisechiselparams *p, 
gal_data_t *workbin,
     }
 
 
-  /* 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;
-    }
+  /* Allocate the space necessary to work on each tile (to avoid having to
+     allocate it it separately for each tile and within each
+     thread. `maxltcontig' is the maximum contiguous patch of memory needed
+     to store all tiles. Finally, since we are working on a `uint8_t' type,
+     the size of each element is only 1 byte. */
+  fho_prm.copyspace=gal_data_malloc_array(GAL_TYPE_UINT8,
+                                          p->cp.numthreads*p->maxltcontig);
 
 
-  /* Open the image. */
-  gal_binary_open(workbin, 1, 4, 1);
-  if(p->detectionname)
+  /* Fill the holes and open on each large tile. When no check image is
+     requested, the two steps can be done independently on each tile, but
+     when a check image is requested, we need to break out of the thread
+     spinning function to save the full image then continue it. */
+  if( p->detectionname )
     {
-      workbin->name="OPENED";
-      gal_fits_img_write(workbin, p->detectionname, NULL, PROGRAM_STRING);
-      workbin->name=NULL;
+      /* Necessary initializations. */
+      bin=gal_data_copy(workbin); /*  - Temporary array for demonstration.  */
+      fho_prm.workbin=bin;        /*  - To pass onto the thread.            */
+      fho_prm.step=1;             /*  - So we can break out of the threads. */
+
+      /* Do each step. */
+      while(fho_prm.step<3)
+        {
+          /* Put a copy of `workbin' into `bin' for every step (only
+             necessary for the second step and after). For the first time
+             it was already copied.*/
+          if(fho_prm.step>1)
+            memcpy(bin->array, workbin->array, workbin->size);
+
+          /* Do the respective step. */
+          gal_threads_spin_off(detection_fill_holes_open, &fho_prm,
+                               p->ltl.tottiles, p->cp.numthreads);
+
+          /* Set the extension name based on the step. */
+          switch(fho_prm.step)
+            {
+            case 1:
+              bin->name="HOLES-FILLED";
+              break;
+            case 2:
+              bin->name="OPENED";
+              break;
+            default:
+              error(EXIT_FAILURE, 0, "a bug! the value %d is not recognized "
+                    "in `detection_pseudo_find'. Please contact us at %s so "
+                    "we can address the issue", fho_prm.step,
+                    PACKAGE_BUGREPORT);
+            }
+
+          /* Write the temporary array into the check image. */
+          gal_fits_img_write(bin, p->detectionname, NULL, PROGRAM_STRING);
+
+          /* Increment the step counter. */
+          ++fho_prm.step;
+        }
+
+      /* Clean up: the array in `bin' should just be replaced with that in
+         `workbin' because it is used in later steps. */
+      free(workbin->array);
+      workbin->array=bin->array;
+      bin->name=bin->array=NULL;
+      gal_data_free(bin);
     }
+  else
+    gal_threads_spin_off(detection_fill_holes_open, &fho_prm,
+                         p->ltl.tottiles, p->cp.numthreads);
+
+
+  /* Clean up. */
+  free(fho_prm.copyspace);
 
 
   /* 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. */
+     `workbin' dataset when working on the Sky. Recall that in this case,
+     the blank pixels are the 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;
diff --git a/bin/noisechisel/main.c b/bin/noisechisel/main.c
index 53565a4..dcdb43e 100644
--- a/bin/noisechisel/main.c
+++ b/bin/noisechisel/main.c
@@ -38,7 +38,7 @@ int
 main (int argc, char *argv[])
 {
   struct timeval t1;
-  struct noisechiselparams p={{{0},0},0};
+  struct noisechiselparams p={{{0},0},{0},0};
 
   /* Set they starting time. */
   time(&p.rawtime);
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 36920e1..4c9ea06 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -44,7 +44,7 @@ struct noisechiselparams
 {
   /* From command-line */
   struct gal_options_common_params cp; /* Common parameters.              */
-  /*struct gal_tile_two_layer_params ltl;*/ /* Large tessellation.        */
+  struct gal_tile_two_layer_params ltl;/* Large tessellation.             */
   char             *inputname;  /* Input filename.                        */
   char            *kernelname;  /* Input kernel filename.                 */
   char                  *khdu;  /* Kernel HDU.                            */
@@ -110,7 +110,9 @@ struct noisechiselparams
 
   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.         */
+  size_t          maxltcontig;  /* Maximum contiguous space for a tile.   */
+  size_t            *maxtsize;  /* Maximum size of a single small tile.   */
+  size_t           *maxltsize;  /* Maximum size of a single large tile.   */
 };
 
 #endif
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 1fb90d1..be03d1d 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -105,15 +105,12 @@ noisechisel_convolve_correct_ch_edges(struct 
noisechiselparams *p)
 {
   struct gal_tile_two_layer_params *tl=&p->cp.tl;
 
-  gal_checkset_check_remove_file("conv-correct.fits", 0, 0);
-  gal_fits_img_write(p->conv, "conv-correct.fits", NULL, PROGRAM_STRING);
-
+  /* Correct the convolved image if necessary. */
   if( tl->totchannels>1 && tl->workoverch==0 )
     gal_convolve_spatial_correct_ch_edge(tl->tiles, p->kernel,
                                          p->cp.numthreads, 1, p->conv);
 
-  gal_fits_img_write(p->conv, "conv-correct.fits", NULL, PROGRAM_STRING);
-
+  /* Inform the user. */
   if(!p->cp.quiet)
     gal_timing_report(NULL, "Corrected convolution of touching channel "
                       "edges", 1);
@@ -153,9 +150,13 @@ noisechisel(struct noisechiselparams *p)
   /* Remove false detections. */
   detection(p);
 
-  /* Find the Sky value and subtract it. */
+  /* Find the Sky value and subtract it from the input and convolved
+     images. */
   noisechisel_find_sky_subtract(p);
 
   /* Correct the convolved image channel edges if necessary. */
   noisechisel_convolve_correct_ch_edges(p);
+
+  /* Do the segmentation. */
+
 }
diff --git a/bin/noisechisel/main.c b/bin/noisechisel/segmentation.c
similarity index 67%
copy from bin/noisechisel/main.c
copy to bin/noisechisel/segmentation.c
index 53565a4..fe6e3f2 100644
--- a/bin/noisechisel/main.c
+++ b/bin/noisechisel/segmentation.c
@@ -23,36 +23,15 @@ 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-internal/timing.h>
-
 #include "main.h"
 
-#include "ui.h"
-#include "noisechisel.h"
-
 
-/* Main function */
-int
-main (int argc, char *argv[])
+void
+segmentation(struct noisechiselparams *p)
 {
-  struct timeval t1;
-  struct noisechiselparams p={{{0},0},0};
-
-  /* Set they starting time. */
-  time(&p.rawtime);
-  gettimeofday(&t1, NULL);
-
-  /* Read the input parameters. */
-  ui_read_check_inputs_setup(argc, argv, &p);
-
-  /* Run MakeProfiles */
-  noisechisel(&p);
-
-  /* Free all non-freed allocations. */
-  ui_free_report(&p, &t1);
 
-  /* Return successfully.*/
-  return EXIT_SUCCESS;
 }
diff --git a/bin/noisechisel/main.c b/bin/noisechisel/segmentation.h
similarity index 62%
copy from bin/noisechisel/main.c
copy to bin/noisechisel/segmentation.h
index 53565a4..a9f4c38 100644
--- a/bin/noisechisel/main.c
+++ b/bin/noisechisel/segmentation.h
@@ -20,39 +20,10 @@ General Public License for more details.
 You should have received a copy of the GNU General Public License
 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
 **********************************************************************/
-#include <config.h>
+#ifndef SEGMENTATION_H
+#define SEGMENTATION_H
 
-#include <stdio.h>
-#include <stdlib.h>
+void
+segmentation(struct noisechiselparams *p);
 
-#include <gnuastro-internal/timing.h>
-
-#include "main.h"
-
-#include "ui.h"
-#include "noisechisel.h"
-
-
-/* Main function */
-int
-main (int argc, char *argv[])
-{
-  struct timeval t1;
-  struct noisechiselparams p={{{0},0},0};
-
-  /* Set they starting time. */
-  time(&p.rawtime);
-  gettimeofday(&t1, NULL);
-
-  /* Read the input parameters. */
-  ui_read_check_inputs_setup(argc, argv, &p);
-
-  /* Run MakeProfiles */
-  noisechisel(&p);
-
-  /* Free all non-freed allocations. */
-  ui_free_report(&p, &t1);
-
-  /* Return successfully.*/
-  return EXIT_SUCCESS;
-}
+#endif
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index b8a03ec..4927fdc 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -145,7 +145,7 @@ ui_initialize_options(struct noisechiselparams *p,
 
 
 /* Parse a single option: */
-error_t
+static error_t
 parse_opt(int key, char *arg, struct argp_state *state)
 {
   struct noisechiselparams *p = state->input;
@@ -323,6 +323,11 @@ ui_set_output_names(struct noisechiselparams *p)
     p->cp.output=gal_checkset_automatic_output(&p->cp, p->inputname,
                                                "_labeled.fits");
 
+  /* Tile check. */
+  if(p->cp.tl.checktiles)
+    p->cp.tl.tilecheckname=gal_checkset_automatic_output(&p->cp, basename,
+                                                         "_tiles.fits");
+
   /* Quantile threshold. */
   if(p->checkqthresh)
     p->qthreshname=gal_checkset_automatic_output(&p->cp, basename,
@@ -449,51 +454,97 @@ ui_prepare_kernel(struct noisechiselparams *p)
 
 
 
-void
-ui_preparations(struct noisechiselparams *p)
+static void
+ui_prepare_tiles(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;
+  struct gal_tile_two_layer_params *tl=&p->cp.tl, *ltl=&p->ltl;
 
-  /* Prepare the names of the outputs. */
-  ui_set_output_names(p);
-
-  /* Read the input as a single precision floating point dataset. */
-  p->input = gal_fits_img_read_to_type(p->inputname, p->cp.hdu,
-                                       GAL_TYPE_FLOAT32,
-                                       p->cp.minmapsize);
-  gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs,
-               &p->input->wcs);
-  if(p->input->name==NULL)
-    gal_checkset_allocate_copy("INPUT", &p->input->name);
 
-  /* Read in the kernel for convolution. */
-  ui_prepare_kernel(p);
-
-  /* Check the tile parameters and make the tile structure. We will also
-     need the dimensions of the tile with the maximum required memory.  */
+  /* Check the tile parameters for the small tile sizes and make the tile
+     structure. We will also need the dimensions of the tile with the
+     maximum required memory. */
   p->maxtsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->input->ndim);
   gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, tl);
   gal_tile_full_two_layers(p->input, tl);
   gal_tile_full_permutation(tl);
   for(check=tl->tiles; check!=NULL; check=check->next)
-    if( check->size > p->maxtcontig )/* p->maxtcontig was initialized to 0. */
+    if( check->size > p->maxtcontig )/* p->maxtcontig initialized to 0. */
       {
         p->maxtcontig=check->size;
         memcpy(p->maxtsize, check->dsize, tl->ndim*sizeof *p->maxtsize);
       }
 
+
+  /* Make the large tessellation, except for the size, the rest of the
+     parameters are the same as the small tile sizes. */
+  ltl->numchannels    = tl->numchannels;
+  ltl->remainderfrac  = tl->remainderfrac;
+  ltl->workoverch     = tl->workoverch;
+  ltl->checktiles     = tl->checktiles;
+  ltl->oneelempertile = tl->oneelempertile;
+  p->maxltsize=gal_data_malloc_array(GAL_TYPE_SIZE_T, p->input->ndim);
+  gal_tile_full_sanity_check(p->inputname, p->cp.hdu, p->input, ltl);
+  gal_tile_full_two_layers(p->input, ltl);
+  gal_tile_full_permutation(ltl);
+  for(check=ltl->tiles; check!=NULL; check=check->next)
+    if( check->size > p->maxltcontig )/* p->maxltcontig initialized to 0. */
+      {
+        p->maxltcontig=check->size;
+        memcpy(p->maxltsize, check->dsize, ltl->ndim*sizeof *p->maxltsize);
+      }
+
   /* Make the tile check image if requested. */
   if(tl->checktiles)
     {
-      tl->tilecheckname=gal_checkset_automatic_output(&p->cp, p->inputname,
-                                                      "_tiled.fits");
+      /* Large tiles. */
+      check=gal_tile_block_check_tiles(ltl->tiles);
+      gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
+      gal_data_free(check);
+
+      /* Small tiles. */
       check=gal_tile_block_check_tiles(tl->tiles);
       gal_fits_img_write(check, tl->tilecheckname, NULL, PROGRAM_NAME);
       gal_data_free(check);
+
+      /* If `continueaftercheck' hasn't been called, abort NoiseChisel. */
+      if(!p->continueaftercheck)
+        ui_abort_after_check(p, tl->tilecheckname,
+                             "showing all tiles over the image");
+
+      /* Free the name. */
       free(tl->tilecheckname);
     }
+}
+
+
+
+
+
+static void
+ui_preparations(struct noisechiselparams *p)
+{
+  /* Prepare the names of the outputs. */
+  ui_set_output_names(p);
+
+
+  /* Read the input as a single precision floating point dataset. */
+  p->input = gal_fits_img_read_to_type(p->inputname, p->cp.hdu,
+                                       GAL_TYPE_FLOAT32,
+                                       p->cp.minmapsize);
+  gal_wcs_read(p->inputname, p->cp.hdu, 0, 0, &p->input->nwcs,
+               &p->input->wcs);
+  if(p->input->name==NULL)
+    gal_checkset_allocate_copy("INPUT", &p->input->name);
+
+
+  /* Read in the kernel for convolution. */
+  ui_prepare_kernel(p);
+
+
+  /* Prepare the tessellation. */
+  ui_prepare_tiles(p);
+
 
   /* Allocate space for the over-all necessary arrays. */
   p->binary=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
@@ -523,11 +574,11 @@ ui_preparations(struct noisechiselparams *p)
 
 
 /**************************************************************/
-/************         Set the parameters          *************/
+/************     High level reading function     *************/
 /**************************************************************/
-
 void
-ui_read_check_inputs_setup(int argc, char *argv[], struct noisechiselparams *p)
+ui_read_check_inputs_setup(int argc, char *argv[],
+                           struct noisechiselparams *p)
 {
   struct gal_options_common_params *cp=&p->cp;
 
@@ -664,6 +715,11 @@ ui_free_report(struct noisechiselparams *p, struct timeval 
*t1)
   gal_data_free(p->binary);
   gal_data_free(p->olabel);
 
+  /* Clean up the tile structure. */
+  p->ltl.numchannels=NULL;
+  gal_tile_full_free_contents(&p->ltl);
+  gal_tile_full_free_contents(&p->cp.tl);
+
   /* Print the final message. */
   if(!p->cp.quiet && t1)
     gal_timing_report(t1, PROGRAM_NAME" finished in: ", 0);
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 6eea6c4..785a1ea 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -30,11 +30,12 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /* Available letters for short options:
 
    a b f j l n u w x z
-   A H J L W X Y
+   A H J W X Y
 */
 enum option_keys_enum
 {
   /* With short-option version. */
+  ARGS_OPTION_KEY_LARGETILESIZE      = 'L',
   ARGS_OPTION_KEY_KERNEL             = 'k',
   ARGS_OPTION_KEY_SKYSUBTRACTED      = 'E',
   ARGS_OPTION_KEY_MINBFRAC           = 'B',
@@ -59,7 +60,7 @@ enum option_keys_enum
 
   /* Only with long version (start with a value 1000, the rest will be set
      automatically). */
-  ARGS_OPTION_KEY_KHDU              = 1000,
+  ARGS_OPTION_KEY_KHDU               = 1000,
   ARGS_OPTION_KEY_MINNUMFALSE,
   ARGS_OPTION_KEY_ONLYDETECT,
   ARGS_OPTION_KEY_GROWNCLUMPS,
diff --git a/lib/tile.c b/lib/tile.c
index e72dcf5..1c2246a 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -487,30 +487,44 @@ gal_tile_full_regular_first(gal_data_t *parent, size_t 
*regular,
   /* For each dimension, set the size of the first tile. */
   for(i=0;i<parent->ndim;++i)
     {
-      /* Calculate the remainder in this dimension. */
-      remainder=dsize[i] % regular[i];
-
-      /* Depending on the remainder, set the first tile size and number. */
-      if(remainder)
+      /* It might happen that the tile size is bigger than the parent size
+         in a dimension, in that case the analysis in the comments above
+         are useless and only one tile should cover this dimension with the
+         size of the parent. */
+      if( regular[i] >= dsize[i] )
         {
-          if( remainder > remainderfrac * regular[i] )
+          tsize[i]=1;
+          first[i]=last[i]=dsize[i];
+        }
+      else
+        {
+          /* Calculate the remainder in this dimension. */
+          remainder=dsize[i] % regular[i];
+
+          /* Depending on the remainder, set the first tile size and
+             number. */
+          if(remainder)
             {
-              first[i]  = ( remainder + regular[i] )/2;
-              tsize[i] = dsize[i]/regular[i] + 1 ;
-              last[i]   = dsize[i] - ( first[i] + regular[i]*(tsize[i]-2) );
+              if( remainder > remainderfrac * regular[i] )
+                {
+                  first[i]  = ( remainder + regular[i] )/2;
+                  tsize[i]  = dsize[i]/regular[i] + 1 ;
+                  last[i]   = ( dsize[i]
+                                - ( first[i] + regular[i]*(tsize[i]-2) ) );
+                }
+              else
+                {
+                  first[i]  = remainder + regular[i];
+                  tsize[i]  = dsize[i]/regular[i];
+                  last[i]   = regular[i];
+                }
             }
           else
             {
-              first[i]  = remainder + regular[i];
-              tsize[i]  = dsize[i]/regular[i];
-              last[i]   = regular[i];
+              first[i]  = last[i] = regular[i];
+              tsize[i] = dsize[i]/regular[i];
             }
         }
-      else
-        {
-          first[i]  = last[i] = regular[i];
-          tsize[i] = dsize[i]/regular[i];
-        }
     }
 }
 
@@ -592,7 +606,6 @@ gal_tile_full(gal_data_t *input, size_t *regular,
                               first, last, tsize);
   numtiles=gal_dimension_total_size(input->ndim, tsize);
 
-
   /* Allocate the necessary space for all the tiles (if necessary). */
   if(*out)        tiles = *out;
   else     *out = tiles = gal_data_array_calloc(numtiles*multiple);
@@ -671,7 +684,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
          next pointer as the next tile. Note that only when we are dealing
          with the last tile should the `next' pointer be set to NULL.*/
       tiles[i].block = input;
-      tiles[i].next = i==numtiles-1 ? NULL : &tiles[i+1];
+      tiles[i].next  = i==numtiles-1 ? NULL : &tiles[i+1];
 
       /* For a check:
       printf("%zu:\n\tStart index: %zu\n\tsize: %zu x %zu\n", i, tind,
@@ -687,9 +700,9 @@ gal_tile_full(gal_data_t *input, size_t *regular,
   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'.           */
-}
+  tsize[input->ndim]=-1; /* `tsize' had ndim+1 values, we will mark the  */
+  return tsize;          /* extra space with the largest possible value: */
+}                        /* -1, see `gal_tile_full_sanity_check'.        */