"Fossies" - the Fresh Open Source Software Archive

Member "gnuastro-0.9/bin/noisechisel/threshold.c" (6 Mar 2019, 24994 Bytes) of package /linux/privat/gnuastro-0.9.tar.lz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "threshold.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 0.8_vs_0.9.

    1 /*********************************************************************
    2 NoiseChisel - Detect signal in a noisy dataset.
    3 NoiseChisel is part of GNU Astronomy Utilities (Gnuastro) package.
    4 
    5 Original author:
    6      Mohammad Akhlaghi <mohammad@akhlaghi.org>
    7 Contributing author(s):
    8 Copyright (C) 2015-2019, Free Software Foundation, Inc.
    9 
   10 Gnuastro is free software: you can redistribute it and/or modify it
   11 under the terms of the GNU General Public License as published by the
   12 Free Software Foundation, either version 3 of the License, or (at your
   13 option) any later version.
   14 
   15 Gnuastro is distributed in the hope that it will be useful, but
   16 WITHOUT ANY WARRANTY; without even the implied warranty of
   17 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   18 General Public License for more details.
   19 
   20 You should have received a copy of the GNU General Public License
   21 along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
   22 **********************************************************************/
   23 #include <config.h>
   24 
   25 #include <stdio.h>
   26 #include <errno.h>
   27 #include <error.h>
   28 #include <string.h>
   29 #include <stdlib.h>
   30 
   31 #include <gnuastro/fits.h>
   32 #include <gnuastro/blank.h>
   33 #include <gnuastro/threads.h>
   34 #include <gnuastro/pointer.h>
   35 #include <gnuastro/statistics.h>
   36 #include <gnuastro/interpolate.h>
   37 
   38 #include <gnuastro-internal/timing.h>
   39 #include <gnuastro-internal/checkset.h>
   40 #include <gnuastro-internal/tile-internal.h>
   41 
   42 #include "main.h"
   43 
   44 #include "ui.h"
   45 #include "threshold.h"
   46 
   47 
   48 
   49 
   50 
   51 
   52 
   53 
   54 
   55 
   56 /**********************************************************************/
   57 /***************        Apply a given threshold.      *****************/
   58 /**********************************************************************/
   59 struct threshold_apply_p
   60 {
   61   float               *value1;
   62   float               *value2;
   63   int                    kind;
   64   struct noisechiselparams *p;
   65 };
   66 
   67 
   68 
   69 
   70 /* Apply the threshold on the tiles given to this thread. */
   71 static void *
   72 threshold_apply_on_thread(void *in_prm)
   73 {
   74   struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
   75   struct threshold_apply_p *taprm=(struct threshold_apply_p *)(tprm->params);
   76   struct noisechiselparams *p=taprm->p;
   77 
   78   size_t i, tid;
   79   void *tarray=NULL;
   80   gal_data_t *tile, *tblock=NULL;
   81   float *value1=taprm->value1, *value2=taprm->value2;
   82 
   83   /* Go over all the tiles assigned to this thread. */
   84   for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
   85     {
   86       /* For easy reading. */
   87       tid=tprm->indexs[i];
   88       tile=&p->cp.tl.tiles[tid];
   89 
   90       /* Based on the kind of threshold. */
   91       switch(taprm->kind)
   92         {
   93 
   94         /* This is a quantile threshold. */
   95         case THRESHOLD_QUANTILES:
   96           /* Correct the tile's pointers to apply the threshold on the
   97              convolved image. */
   98           if(p->conv)
   99             {
  100               tarray=tile->array; tblock=tile->block;
  101               tile->array=gal_tile_block_relative_to_other(tile, p->conv);
  102               tile->block=p->conv;
  103             }
  104 
  105           /* Apply the threshold: When the `>' comparison fails, it can be
  106              either because the pixel was actually smaller than the
  107              threshold, or that it was a NaN value. In the first case,
  108              return 0, in the second, return a blank value.
  109 
  110              We already know if a tile contains a blank value (which is a
  111              constant over the whole loop). So before checking if the value
  112              is blank, see if the tile actually has a blank value. This
  113              will help in efficiency, because the compiler can move this
  114              check out of the loop and only check for NaN values when we
  115              know the tile has blank pixels. */
  116           GAL_TILE_PO_OISET(float, uint8_t, tile, p->binary, 1, 0, {
  117               *o = ( *i > value1[tid]
  118                      ? ( *i > value2[tid] ? THRESHOLD_NO_ERODE_VALUE : 1 )
  119                      : ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
  120                          ? GAL_BLANK_UINT8 : 0 ) );
  121             });
  122 
  123           /* Revert the tile's pointers back to what they were. */
  124           if(p->conv) { tile->array=tarray; tile->block=tblock; }
  125           break;
  126 
  127 
  128         /* This is a Sky and Sky STD threshold. */
  129         case THRESHOLD_SKY_STD:
  130 
  131           /* See the explanation above the same step in the quantile
  132              threshold for an explanation. */
  133           GAL_TILE_PO_OISET(float, uint8_t, tile, p->binary, 1, 0, {
  134               *o = ( ( *i - value1[tid] > p->dthresh * value2[tid] )
  135                      ? 1
  136                      : ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
  137                          ? GAL_BLANK_UINT8 : 0 ) );
  138             });
  139           break;
  140 
  141 
  142         default:
  143           error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we "
  144                 "can address the problem. A value of %d had for "
  145                 "`taprm->kind' is not valid", __func__, PACKAGE_BUGREPORT,
  146                 taprm->kind);
  147         }
  148     }
  149 
  150   /* Wait until all the other threads finish. */
  151   if(tprm->b) pthread_barrier_wait(tprm->b);
  152   return NULL;
  153 }
  154 
  155 
  156 
  157 
  158 
  159 /* Apply a given threshold threshold on the tiles. */
  160 void
  161 threshold_apply(struct noisechiselparams *p, float *value1,
  162                 float *value2, int kind)
  163 {
  164   struct threshold_apply_p taprm={value1, value2, kind, p};
  165   gal_threads_spin_off(threshold_apply_on_thread, &taprm, p->cp.tl.tottiles,
  166                        p->cp.numthreads);
  167 }
  168 
  169 
  170 
  171 
  172 
  173 
  174 
  175 
  176 
  177 
  178 
  179 
  180 
  181 
  182 
  183 
  184 
  185 
  186 
  187 /**********************************************************************/
  188 /***************            Write S/N values          *****************/
  189 /**********************************************************************/
  190 void
  191 threshold_write_sn_table(struct noisechiselparams *p, gal_data_t *insn,
  192                          gal_data_t *inind, char *filename,
  193                          gal_list_str_t *comments, char *extname)
  194 {
  195   gal_data_t *sn, *ind, *cols;
  196 
  197   /* Remove all blank elements. The index and sn values must have the same
  198      set of blank elements, but checking on the integer array is faster. */
  199   if( gal_blank_present(inind, 1) )
  200     {
  201       /* Remove blank elements. */
  202       ind=gal_data_copy(inind);
  203       sn=gal_data_copy(insn);
  204       gal_blank_remove(ind);
  205       gal_blank_remove(sn);
  206 
  207       /* A small sanity check. */
  208       if(ind->size==0 || sn->size==0)
  209         error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
  210               "the problem. For some reason, all the elements in `ind' or "
  211               "`sn' are blank", __func__, PACKAGE_BUGREPORT);
  212     }
  213   else
  214     {
  215       sn  = insn;
  216       ind = inind;
  217     }
  218 
  219   /* Set the columns. */
  220   cols       = ind;
  221   cols->next = sn;
  222 
  223 
  224   /* Prepare the comments. */
  225   gal_table_comments_add_intro(&comments, PROGRAM_STRING, &p->rawtime);
  226 
  227 
  228   /* Write the table. Note that we'll set the `dontdelete' argument to 0
  229      because when the output is a FITS table, we want all the tables in one
  230      FITS file. We have already deleted any existing file with the same
  231      name in `ui_set_output_names'.*/
  232   gal_table_write(cols, comments, p->cp.tableformat, filename, extname, 0);
  233 
  234 
  235   /* Clean up (if necessary). */
  236   if(sn!=insn) gal_data_free(sn);
  237   if(ind==inind) ind->next=NULL; else gal_data_free(ind);
  238 }
  239 
  240 
  241 
  242 
  243 
  244 
  245 
  246 
  247 
  248 
  249 
  250 
  251 
  252 
  253 
  254 
  255 
  256 
  257 
  258 /**********************************************************************/
  259 /***************     Interpolation and smoothing      *****************/
  260 /**********************************************************************/
  261 /* Interpolate and smooth the values for each tile over the whole image. */
  262 void
  263 threshold_interp_smooth(struct noisechiselparams *p, gal_data_t **first,
  264                         gal_data_t **second, gal_data_t **third,
  265                         char *filename)
  266 {
  267   gal_data_t *tmp;
  268   struct gal_options_common_params *cp=&p->cp;
  269   struct gal_tile_two_layer_params *tl=&cp->tl;
  270 
  271   /* A small sanity check. */
  272   if( (*first)->next )
  273     error(EXIT_FAILURE, 0, "%s: `first' must not have any `next' pointer.",
  274           __func__);
  275   if( (*second)->next )
  276     error(EXIT_FAILURE, 0, "%s: `second' must not have any `next' pointer.",
  277           __func__);
  278   if( third && (*third)->next )
  279     error(EXIT_FAILURE, 0, "%s: `third' must not have any `next' pointer.",
  280           __func__);
  281 
  282   /* Do the interpolation of both arrays. */
  283   (*first)->next = *second;
  284   if(third) (*second)->next = *third;
  285   tmp=gal_interpolate_close_neighbors(*first, tl, cp->interpmetric,
  286                                       cp->interpnumngb, cp->numthreads,
  287                                       cp->interponlyblank, 1);
  288   gal_data_free(*first);
  289   gal_data_free(*second);
  290   if(third) gal_data_free(*third);
  291   *first=tmp;
  292   *second=(*first)->next;
  293   if(third)
  294     {
  295       *third=(*second)->next;
  296       (*third)->next=NULL;
  297     }
  298   (*first)->next=(*second)->next=NULL;
  299   if(filename)
  300     {
  301       (*first)->name="THRESH1_INTERP";
  302       (*second)->name="THRESH2_INTERP";
  303       if(third) (*third)->name="THRESH3_INTERP";
  304       gal_tile_full_values_write(*first, tl, !p->ignoreblankintiles,
  305                                  filename, NULL, PROGRAM_NAME);
  306       gal_tile_full_values_write(*second, tl, !p->ignoreblankintiles,
  307                                  filename, NULL, PROGRAM_NAME);
  308       if(third)
  309         gal_tile_full_values_write(*third, tl, !p->ignoreblankintiles,
  310                                    filename, NULL, PROGRAM_NAME);
  311       (*first)->name = (*second)->name = NULL;
  312       if(third) (*third)->name=NULL;
  313     }
  314 
  315   /* Smooth the threshold if requested. */
  316   if(p->smoothwidth>1)
  317     {
  318       /* Smooth the first. */
  319       tmp=gal_tile_full_values_smooth(*first, tl, p->smoothwidth,
  320                                       p->cp.numthreads);
  321       gal_data_free(*first);
  322       *first=tmp;
  323 
  324       /* Smooth the second */
  325       tmp=gal_tile_full_values_smooth(*second, tl, p->smoothwidth,
  326                                       p->cp.numthreads);
  327       gal_data_free(*second);
  328       *second=tmp;
  329 
  330       /* Smooth the third */
  331       if(third)
  332         {
  333           tmp=gal_tile_full_values_smooth(*third, tl, p->smoothwidth,
  334                                           p->cp.numthreads);
  335           gal_data_free(*third);
  336           *third=tmp;
  337         }
  338 
  339       /* Add them to the check image. */
  340       if(filename)
  341         {
  342           (*first)->name="THRESH1_SMOOTH";
  343           (*second)->name="THRESH2_SMOOTH";
  344           if(third) (*third)->name="THRESH3_SMOOTH";
  345           gal_tile_full_values_write(*first, tl, !p->ignoreblankintiles,
  346                                      filename, NULL, PROGRAM_NAME);
  347           gal_tile_full_values_write(*second, tl, !p->ignoreblankintiles,
  348                                      filename, NULL, PROGRAM_NAME);
  349           if(third)
  350             gal_tile_full_values_write(*third, tl, !p->ignoreblankintiles,
  351                                        filename, NULL, PROGRAM_NAME);
  352           (*first)->name = (*second)->name = NULL;
  353           if(third) (*third)->name=NULL;
  354         }
  355     }
  356 }
  357 
  358 
  359 
  360 
  361 
  362 
  363 
  364 
  365 
  366 
  367 
  368 
  369 
  370 
  371 
  372 
  373 
  374 
  375 
  376 
  377 /****************************************************************
  378  ************           Quantile threshold           ************
  379  ****************************************************************/
  380 struct qthreshparams
  381 {
  382   gal_data_t        *erode_th;
  383   gal_data_t      *noerode_th;
  384   gal_data_t       *expand_th;
  385   void                 *usage;
  386   struct noisechiselparams *p;
  387 };
  388 
  389 
  390 
  391 
  392 
  393 static void *
  394 qthresh_on_tile(void *in_prm)
  395 {
  396   struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
  397   struct qthreshparams *qprm=(struct qthreshparams *)tprm->params;
  398   struct noisechiselparams *p=qprm->p;
  399 
  400   void *tarray=NULL;
  401   int type=qprm->erode_th->type;
  402   gal_data_t *meanconv = p->wconv ? p->wconv : p->conv;
  403   size_t i, tind, twidth=gal_type_sizeof(type), ndim=p->input->ndim;
  404   gal_data_t *tile, *mean, *num, *meanquant, *qvalue, *usage, *tblock=NULL;
  405 
  406   /* Put the temporary usage space for this thread into a data set for easy
  407      processing. */
  408   usage=gal_data_alloc(gal_pointer_increment(qprm->usage,
  409                                              tprm->id*p->maxtcontig, type),
  410                        type, ndim, p->maxtsize, NULL, 0, p->cp.minmapsize,
  411                        NULL, NULL, NULL);
  412 
  413   /* Go over all the tiles given to this thread. */
  414   for(i=0; tprm->indexs[i] != GAL_BLANK_SIZE_T; ++i)
  415     {
  416       /* Re-initialize the usage array's space (will be changed in
  417          `gal_data_copy_to_allocated' for each tile). */
  418       usage->ndim=ndim;
  419       usage->size=p->maxtcontig;
  420       memcpy(usage->dsize, p->maxtsize, ndim*sizeof *p->maxtsize);
  421 
  422 
  423       /* For easy reading. */
  424       tind = tprm->indexs[i];
  425       tile = &p->cp.tl.tiles[tind];
  426 
  427 
  428       /* Temporarily change the tile's pointers so we can do the work on
  429          the convolved image, then copy the desired contents into the
  430          already allocated `usage' array. */
  431       tarray=tile->array; tblock=tile->block;
  432       tile->array=gal_tile_block_relative_to_other(tile, meanconv);
  433       tile->block=meanconv;
  434       gal_data_copy_to_allocated(tile, usage);
  435       tile->array=tarray;
  436       tile->block=tblock;
  437 
  438 
  439       /* Find the mean's quantile on this tile, note that we have already
  440          copied the tile's dataset to a newly allocated place. So we have
  441          set the `inplace' flag to `1' to avoid extra allocation. */
  442       mean=gal_statistics_mean(usage);
  443       num=gal_statistics_number(usage);
  444       mean=gal_data_copy_to_new_type_free(mean, usage->type);
  445       meanquant = ( *(size_t *)(num->array)
  446                     ? gal_statistics_quantile_function(usage, mean, 1)
  447                     : NULL );
  448 
  449       /* Only continue if the mean's quantile is close enough to the
  450          median.  */
  451       if( meanquant
  452           && fabs( *(double *)(meanquant->array)-0.5f) < p->meanmedqdiff )
  453         {
  454           /* The mean was found on the wider convolved image, but the
  455              qthresh values have to be found on the sharper convolved
  456              images. This is because the distribution becomes more skewed
  457              with a wider kernel, helping us find tiles with no data more
  458              easily. But for the quantile threshold, we want to use the
  459              sharper convolved image to loose less of the spatial
  460              information. */
  461           if(meanconv!=p->conv)
  462             {
  463               tarray=tile->array; tblock=tile->block;
  464               tile->array=gal_tile_block_relative_to_other(tile, p->conv);
  465               tile->block=p->conv;
  466               usage->ndim=ndim;             /* Since usage was modified in */
  467               usage->size=p->maxtcontig;    /* place, it needs to be       */
  468               gal_data_copy_to_allocated(tile, usage);  /* re-initialized. */
  469               tile->array=tarray; tile->block=tblock;
  470             }
  471 
  472           /* Get the erosion quantile for this tile and save it. Note that
  473              the type of `qvalue' is the same as the input dataset. */
  474           qvalue=gal_statistics_quantile(usage, p->qthresh, 1);
  475           memcpy(gal_pointer_increment(qprm->erode_th->array, tind, type),
  476                  qvalue->array, twidth);
  477           gal_data_free(qvalue);
  478 
  479           /* Same for the no-erode quantile. */
  480           qvalue=gal_statistics_quantile(usage, p->noerodequant, 1);
  481           memcpy(gal_pointer_increment(qprm->noerode_th->array, tind, type),
  482                  qvalue->array, twidth);
  483           gal_data_free(qvalue);
  484 
  485           /* Same for the expansion quantile. */
  486           if(qprm->expand_th)
  487             {
  488               qvalue=gal_statistics_quantile(usage, p->detgrowquant, 1);
  489               memcpy(gal_pointer_increment(qprm->expand_th->array, tind,
  490                                             type),
  491                      qvalue->array, twidth);
  492               gal_data_free(qvalue);
  493             }
  494         }
  495       else
  496         {
  497           gal_blank_write(gal_pointer_increment(qprm->erode_th->array,
  498                                                  tind, type), type);
  499           gal_blank_write(gal_pointer_increment(qprm->noerode_th->array,
  500                                                  tind, type), type);
  501           if(qprm->expand_th)
  502             gal_blank_write(gal_pointer_increment(qprm->expand_th->array,
  503                                                    tind, type), type);
  504         }
  505 
  506       /* Clean up and fix the tile's pointers. */
  507       gal_data_free(num);
  508       gal_data_free(mean);
  509       gal_data_free(meanquant);
  510     }
  511 
  512   /* Clean up and wait for the other threads to finish, then return. */
  513   usage->array=NULL;  /* Not allocated here. */
  514   gal_data_free(usage);
  515   if(tprm->b) pthread_barrier_wait(tprm->b);
  516   return NULL;
  517 }
  518 
  519 
  520 
  521 
  522 
  523 static void
  524 threshold_good_error(size_t number, int before0_after1, size_t interpnumngb)
  525 {
  526   before0_after1=1;
  527 
  528   /* Set the differing strings. */
  529   char *in1 = ( before0_after1
  530                 ? "after removing outliers"
  531                 : "for defining a quantile threshold" );
  532   char *in2 = ( before0_after1
  533                 ? ""
  534                 : "NOTE that this is happening *BEFORE* outlier rejection "
  535                   "(where the number may decrease even further).");
  536   char *in3 = ( before0_after1
  537                 ? "\n"
  538                   "  - (slightly) Increase `--outliersclip' to reject less "
  539                   "as outliers.\n"
  540                   "  - (slightly) Increase `--outliersigma' to reject less "
  541                   "as outliers.\n"
  542                 : "\n");
  543 
  544   /* Print the error message and abort. */
  545   error(EXIT_FAILURE, 0, "%zu tiles usable %s!\n\n"
  546 
  547         "This is smaller than the requested minimum value of %zu (value to "
  548         "the `--interpnumngb' option). %s\n\n"
  549 
  550         "There are several ways to address the problem. The best and most "
  551         "highly recommended is to use a larger input if possible (when the "
  552         "input is a crop from a larger dataset). If this is not the case, "
  553         "or it doesn't solve the problem, you need to loosen the "
  554         "parameters mentioned below in the respective order (and therefore "
  555         "cause scatter/inaccuracy in the final result). Hence its best to "
  556         "not loosen them too much (recall that you can see all the option "
  557         "values to Gnuastro's programs by appending `-P' to the end of your "
  558         "command).\n"
  559         "  - (slightly) Decrease `--tilesize' so your tile-grid has more "
  560         "tiles.\n"
  561         "  - (slightly) Increase `--meanmedqdiff' to accept more tiles.%s"
  562         "  - (slightly) Decrease `--interpnumngb' to be less than %zu.\n\n"
  563 
  564         "---- Tip ----\n"
  565         "Append your command with `--checkqthresh' to see the "
  566         "successful tiles in relation with this dataset's contents "
  567         "before this crash. A visual inspection will greatly help in "
  568         "finding the cause/solution for this particular dataset (note "
  569         "that the output of `--checkqthresh' is a multi-extension FITS "
  570         "file).\n\n"
  571         "To better understand this important step, please run the "
  572         "following command (press `SPACE'/arrow-keys to navigate and "
  573         "`Q' to return back to the command-line):\n\n"
  574         "    $ info gnuastro \"Quantifying signal in a tile\"\n", number,
  575         in1, interpnumngb, in2, in3, number);
  576 }
  577 
  578 
  579 
  580 
  581 
  582 void
  583 threshold_quantile_find_apply(struct noisechiselparams *p)
  584 {
  585   char *msg;
  586   gal_data_t *num;
  587   struct timeval t1;
  588   size_t nval, nblank;
  589   struct qthreshparams qprm;
  590   struct gal_options_common_params *cp=&p->cp;
  591   struct gal_tile_two_layer_params *tl=&cp->tl;
  592 
  593 
  594   /* Get the starting time if necessary. */
  595   if(!p->cp.quiet) gettimeofday(&t1, NULL);
  596 
  597 
  598   /* Add image to check image if requested. If the user has asked for
  599      `oneelempertile', then the size of values is not going to be the same
  600      as the input, making it hard to inspect visually. So we'll only put
  601      the full input when `oneelempertile' isn't requested. */
  602   if(p->qthreshname && !tl->oneelempertile)
  603     {
  604       gal_fits_img_write(p->conv ? p->conv : p->input, p->qthreshname, NULL,
  605                          PROGRAM_NAME);
  606       if(p->wconv)
  607         gal_fits_img_write(p->wconv ? p->wconv : p->input, p->qthreshname,
  608                            NULL, PROGRAM_NAME);
  609     }
  610 
  611 
  612   /* Allocate space for the quantile threshold values. */
  613   qprm.erode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
  614                                tl->numtiles, NULL, 0, cp->minmapsize,
  615                                NULL, p->input->unit, NULL);
  616   qprm.noerode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
  617                                  tl->numtiles, NULL, 0, cp->minmapsize,
  618                                  NULL, p->input->unit, NULL);
  619   qprm.expand_th = ( p->detgrowquant!=1.0f
  620                      ? gal_data_alloc(NULL, p->input->type, p->input->ndim,
  621                                       tl->numtiles, NULL, 0, cp->minmapsize,
  622                                       NULL, p->input->unit, NULL)
  623                      : NULL );
  624 
  625 
  626   /* Allocate temporary space for processing in each tile. */
  627   qprm.usage=gal_pointer_allocate(p->input->type,
  628                                   cp->numthreads * p->maxtcontig, 0,
  629                                   __func__, "qprm.usage");
  630 
  631 
  632   /* Find the threshold on each tile, free the temporary processing space
  633      and set the blank flag on both. Since they have the same blank
  634      elements, it is only necessary to check one (with the `updateflag'
  635      value set to 1), then update the next. */
  636   qprm.p=p;
  637   gal_threads_spin_off(qthresh_on_tile, &qprm, tl->tottiles, cp->numthreads);
  638   free(qprm.usage);
  639   if( gal_blank_present(qprm.erode_th, 1) )
  640     {
  641       qprm.noerode_th->flag |= GAL_DATA_FLAG_HASBLANK;
  642       if(qprm.expand_th) qprm.expand_th->flag  |= GAL_DATA_FLAG_HASBLANK;
  643     }
  644   qprm.noerode_th->flag |= GAL_DATA_FLAG_BLANK_CH;
  645   if(qprm.expand_th) qprm.expand_th->flag  |= GAL_DATA_FLAG_BLANK_CH;
  646   if(p->qthreshname)
  647     {
  648       qprm.erode_th->name="QTHRESH_ERODE";
  649       qprm.noerode_th->name="QTHRESH_NOERODE";
  650       gal_tile_full_values_write(qprm.erode_th, tl,
  651                                  !p->ignoreblankintiles,
  652                                  p->qthreshname, NULL, PROGRAM_NAME);
  653       gal_tile_full_values_write(qprm.noerode_th, tl,
  654                                  !p->ignoreblankintiles,
  655                                  p->qthreshname, NULL, PROGRAM_NAME);
  656       qprm.erode_th->name=qprm.noerode_th->name=NULL;
  657 
  658       if(qprm.expand_th)
  659         {
  660           qprm.expand_th->name="QTHRESH_EXPAND";
  661           gal_tile_full_values_write(qprm.expand_th, tl,
  662                                      !p->ignoreblankintiles,
  663                                      p->qthreshname, NULL, PROGRAM_NAME);
  664           qprm.expand_th->name=NULL;
  665         }
  666     }
  667 
  668   /* A small sanity check. */
  669   nblank=gal_blank_number(qprm.erode_th, 1);
  670   if( nblank > qprm.erode_th->size-cp->interpnumngb )
  671     threshold_good_error(qprm.erode_th->size-nblank, 0, cp->interpnumngb);
  672 
  673   /* Remove outliers if requested. */
  674   if(p->outliersigma!=0.0)
  675     gal_tileinternal_no_outlier(qprm.erode_th, qprm.noerode_th,
  676                                 qprm.expand_th, &p->cp.tl, p->outliersclip,
  677                                 p->outliersigma, p->qthreshname);
  678 
  679 
  680   /* Check if the number of acceptable tiles is more than the minimum
  681      interpolated number. Since this is a common problem for users, it is
  682      much more useful to do the check here rather than printing multiple
  683      errors in parallel. */
  684   num=gal_statistics_number(qprm.erode_th);
  685   nval=((size_t *)(num->array))[0];
  686   if( nval < cp->interpnumngb )
  687     threshold_good_error(nval, 1, cp->interpnumngb);
  688 
  689 
  690   /* Interpolate and smooth the derived values. */
  691   threshold_interp_smooth(p, &qprm.erode_th, &qprm.noerode_th,
  692                           qprm.expand_th ? &qprm.expand_th : NULL,
  693                           p->qthreshname);
  694 
  695 
  696   /* We now have a threshold for all tiles, apply it. */
  697   threshold_apply(p, qprm.erode_th->array, qprm.noerode_th->array,
  698                   THRESHOLD_QUANTILES);
  699 
  700 
  701   /* Write the binary image if check is requested. */
  702   if(p->qthreshname && !tl->oneelempertile)
  703     {
  704       p->binary->name="QTHRESH-APPLIED";
  705       gal_fits_img_write(p->binary, p->qthreshname, NULL, PROGRAM_NAME);
  706       p->binary->name=NULL;
  707     }
  708 
  709 
  710   /* Set the expansion quantile if necessary. */
  711   p->expand_thresh = qprm.expand_th ? qprm.expand_th : NULL;
  712 
  713 
  714   /* Clean up and report duration if necessary. */
  715   gal_data_free(qprm.erode_th);
  716   gal_data_free(qprm.noerode_th);
  717   if(!p->cp.quiet)
  718     {
  719       if( asprintf(&msg, "%.2f & %0.2f quantile thresholds applied.",
  720                    p->qthresh, p->noerodequant)<0 )
  721         error(EXIT_FAILURE, 0, "%s: asprintf allocation", __func__);
  722       gal_timing_report(&t1, msg, 2);
  723       free(msg);
  724     }
  725 
  726 
  727   /* If the user wanted to check the threshold and hasn't called
  728      `continueaftercheck', then stop NoiseChisel. */
  729   if(p->qthreshname && !p->continueaftercheck)
  730     ui_abort_after_check(p, p->qthreshname, NULL, "quantile threshold check");
  731 }