[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master da0249e 3/3: Merged recent additions/changes i
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master da0249e 3/3: Merged recent additions/changes in NoiseChisel |
|
Date: |
Thu, 14 Sep 2017 19:27:23 -0400 (EDT) |
branch: master
commit da0249eaff04bf43fc7ce75f89d10fb36a22dea4
Merge: 36f9ad2 9738849
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Merged recent additions/changes in NoiseChisel
Several new features were added and changed in NoiseChisel on a separate
branch (for testing). They appear to give reasonable results, so the
changes are now merged with master.
---
NEWS | 22 ++++
bin/noisechisel/args.h | 62 ++++++++----
bin/noisechisel/astnoisechisel.conf | 45 +++++----
bin/noisechisel/clumps.c | 49 ++++++---
bin/noisechisel/clumps.h | 2 +-
bin/noisechisel/detection.c | 194 +++++++++++++++++++++++++-----------
bin/noisechisel/main.h | 15 ++-
bin/noisechisel/noisechisel.c | 27 +++--
bin/noisechisel/segmentation.c | 5 +-
bin/noisechisel/sky.c | 8 +-
bin/noisechisel/threshold.c | 155 ++++++++++++++++++++++------
bin/noisechisel/threshold.h | 3 +-
bin/noisechisel/ui.c | 74 +++++++++-----
bin/noisechisel/ui.h | 10 +-
bin/statistics/sky.c | 16 +--
bin/statistics/statistics.c | 2 +-
doc/gnuastro.texi | 163 ++++++++++++++++++++++--------
lib/binary.c | 31 ++++--
lib/blank.c | 8 +-
lib/gnuastro/binary.h | 2 +-
lib/gnuastro/tile.h | 6 +-
lib/statistics.c | 2 +-
lib/tile.c | 18 ++--
tests/noisechisel/noisechisel.sh | 4 +-
24 files changed, 657 insertions(+), 266 deletions(-)
diff --git a/NEWS b/NEWS
index e3b992d..9f9922f 100644
--- a/NEWS
+++ b/NEWS
@@ -5,6 +5,20 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** New features
+ NoiseChisel: with the new `--widekernel' option it is now possible to use
+ a wider kernel to identify which tiles contain signal. The rest of the
+ steps (identifying the quantile threshold on the selected tiles and etc)
+ are done on the dataset convolved with `--kernel' as they were
+ before. Since it is time consuming, this is an optional feature.
+
+ NoiseChisel: the finally selected true detections are now grown based on
+ signal contiguity, not by blind dilation. The growth process is the same
+ as the growing of clumps to define objects. Only for true detections, the
+ growth occurs in the noise. You can configure this growth with the
+ `--detgrowquant' and `--detgrowmaxholesize'. With this new feature it is
+ now possible to detect signal out to much lower surface brightness limits
+ and the detections don't look boxy any more.
+
`gal_fits_key_img_blank': returns the value that must be used in the
BLANK keyword for the given type as defined by the FITS standard.
@@ -12,6 +26,14 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** Changed features
+ `gal_binary_fill_holes' now accepts a `connectivity' and `maxsize'
+ argument to specify the connectivity of the holes and the maximum size of
+ acceptable holes to fill.
+
+ `gal_tile_block_write_const_value' and `gal_tile_full_values_write' now
+ accept a new `withblank' option to set all pixels that are blank in the
+ tile's block to be blank in the check image.
+
** Bug fixes
ConvertType crash when changing values (bug #52010).
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 772b09a..80de2ce 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -38,7 +38,7 @@ struct argp_option program_options[] =
UI_KEY_KERNEL,
"STR",
0,
- "Filename of Kernel to convolve with input",
+ "Filename of kernel to convolve with input",
GAL_OPTIONS_GROUP_INPUT,
&p->kernelname,
GAL_TYPE_STRING,
@@ -51,7 +51,7 @@ struct argp_option program_options[] =
UI_KEY_KHDU,
"STR",
0,
- "HDU containing Kernel image.",
+ "HDU containing kernel image.",
GAL_OPTIONS_GROUP_INPUT,
&p->khdu,
GAL_TYPE_STRING,
@@ -60,6 +60,32 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "widekernel",
+ UI_KEY_WIDEKERNEL,
+ "STR",
+ 0,
+ "Filename of wider kernel for better qthresh",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->widekernelname,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
+ "wkhdu",
+ UI_KEY_WKHDU,
+ "STR",
+ 0,
+ "HDU containing wide kernel image.",
+ GAL_OPTIONS_GROUP_INPUT,
+ &p->wkhdu,
+ GAL_TYPE_STRING,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"skysubtracted",
UI_KEY_SKYSUBTRACTED,
0,
@@ -365,39 +391,39 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "dilate",
- UI_KEY_DILATE,
- "INT",
+ "detgrowquant",
+ UI_KEY_DETGROWQUANT,
+ "FLT",
0,
- "Number of times to dilate true detections.",
+ "Minimum quant. to expand true detections.",
ARGS_GROUP_DETECTION,
- &p->dilate,
- GAL_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_ANY,
+ &p->detgrowquant,
+ GAL_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GE_0_LE_1,
GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
},
{
- "dilatengb",
- UI_KEY_DILATENGB,
+ "detgrowmaxholesize",
+ UI_KEY_DETGROWMAXHOLESIZE,
"INT",
0,
- "4 or 8 connectivity in final dilation.",
+ "Max. area of holes after growth to fill.",
ARGS_GROUP_DETECTION,
- &p->dilatengb,
+ &p->detgrowmaxholesize,
GAL_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_RANGE_GE_0,
GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
},
{
- "cleandilated",
- UI_KEY_CLEANDILATED,
+ "cleangrowndet",
+ UI_KEY_CLEANGROWNDET,
0,
0,
- "Remove small S/N dilated objects.",
+ "Remove small S/N grown detections.",
ARGS_GROUP_DETECTION,
- &p->cleandilated,
+ &p->cleangrowndet,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_ANY,
GAL_OPTIONS_NOT_MANDATORY,
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 6217e43..a037355 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -19,6 +19,7 @@
# Input:
khdu 1
+ wkhdu 1
minskyfrac 0.7
minnumfalse 100
@@ -26,29 +27,29 @@
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
- dilatengb 8
+ 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
+ detgrowquant 0.65
+ detgrowmaxholesize 100
# 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
+ continueaftercheck 0
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index e2a5d9c..df7b383 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -532,19 +532,34 @@ clumps_grow_prepare_final(struct clumps_thread_params
*cltprm)
This function is going to be used before identifying objects and also
after it (to completely fill in the diffuse area). The distinguishing
point between these two steps is the presence of rivers, so you can use
- the `withrivers' argument. */
+ the `withrivers' argument.
+
+
+ Input:
+
+ labels: The labels array that must be operated on. The pixels that
+ must be "grown" must have the value `CLUMPS_INIT' (negative).
+
+ diffuseindexs: The indexs of the pixels that must be grown.
+
+ withrivers: as described abvoe.
+*/
void
-clumps_grow(struct clumps_thread_params *cltprm, int withrivers)
+clumps_grow(gal_data_t *labels, gal_data_t *diffuseindexs, int withrivers)
{
- struct noisechiselparams *p=cltprm->clprm->p;
- gal_data_t *diffuseindexs=cltprm->diffuseindexs;
- size_t ndim=p->input->ndim, *dsize=p->input->dsize;
-
int searchngb;
- size_t *diarray=cltprm->diffuseindexs->array;
- int32_t n1, nlab, *olabel=p->olabel->array;
- size_t *dinc=gal_dimension_increment(ndim, dsize);
+ size_t *diarray=diffuseindexs->array;
+ int32_t n1, nlab, *olabel=labels->array;
size_t *s, *sf, thisround, ndiffuse=diffuseindexs->size;
+ size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
+
+ /* A small sanity check: */
+ if(labels->type!=GAL_TYPE_INT32)
+ error(EXIT_FAILURE, 0, "%s: `labels' has to have type of int32_t",
+ __func__);
+ if(diffuseindexs->type!=GAL_TYPE_SIZE_T)
+ error(EXIT_FAILURE, 0, "%s: `diffuseindexs' has to have type of size_t",
+ __func__);
/* The basic idea is this: after growing, not all the blank pixels are
necessarily filled, for example the pixels might belong to two regions
@@ -580,7 +595,7 @@ clumps_grow(struct clumps_thread_params *cltprm, int
withrivers)
in a 2D image). Note that since this macro has multiple loops
within it, we can't use break. We'll use a variable instead. */
searchngb=1;
- GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, 1, dinc,
+ GAL_DIMENSION_NEIGHBOR_OP(*s, labels->ndim, labels->dsize, 1, dinc,
{
if(searchngb)
{
@@ -926,6 +941,11 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
size_t len=cltprm->numinitclumps+1;
struct noisechiselparams *p=cltprm->clprm->p;
+ /* If there are no clumps in this tile, then this function can be
+ ignored. */
+ if(cltprm->snind->size==0) return;
+
+
/* A small sanity check. */
if(gal_tile_block(tile)!=p->clabel)
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to address "
@@ -964,7 +984,6 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
do { ninds[*l]=curlab++; *l=ninds[*l]; } while(++l<lf);
-
/* Go over this tile and correct the values. */
GAL_TILE_PARSE_OPERATE( tile, NULL, 0, 1,
{if(*i>0) *i=ninds[ *(int32_t *)i ];} );
@@ -1030,7 +1049,11 @@ clumps_find_make_sn_table(void *in_prm)
/* Find the number of detected pixels over this tile. Since this is
- the binary image, this is just the sum of all the pixels. */
+ the binary image, this is just the sum of all the pixels.
+
+ Note that `numdet' can be `nan' when the whole tile is blank and
+ so there was no values to sum. Recall that in summing, when there
+ is not input, the output is `nan'. */
tmp=gal_statistics_sum(tile);
numdet=*((double *)(tmp->array));
gal_data_free(tmp);
@@ -1049,6 +1072,7 @@ clumps_find_make_sn_table(void *in_prm)
NULL, 0, p->cp.minmapsize, NULL, NULL,
NULL);
+
/* Change the tile's block to the clump labels dataset (because
we'll need to set the labels of the rivers on the edge of the
tile here). */
@@ -1280,7 +1304,6 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
p->ltl.tottiles, p->cp.numthreads);
}
-
/* Destroy the mutex if it was initialized. */
if( p->cp.numthreads>1 && (p->checksegmentation || p->checkclumpsn) )
pthread_mutex_destroy(&clprm.labmutex);
diff --git a/bin/noisechisel/clumps.h b/bin/noisechisel/clumps.h
index 85e5d57..1557955 100644
--- a/bin/noisechisel/clumps.h
+++ b/bin/noisechisel/clumps.h
@@ -79,7 +79,7 @@ void
clumps_grow_prepare_final(struct clumps_thread_params *cltprm);
void
-clumps_grow(struct clumps_thread_params *cltprm, int withrivers);
+clumps_grow(gal_data_t *labels, gal_data_t *diffuseindexs, int withrivers);
void
clumps_true_find_sn_thresh(struct noisechiselparams *p);
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index cffe814..2c5b88b 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -40,6 +40,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "ui.h"
#include "sky.h"
+#include "clumps.h"
#include "threshold.h"
@@ -240,8 +241,9 @@ detection_fill_holes_open(void *in_prm)
copy->size=p->maxltcontig;
gal_data_copy_to_allocated(tile, copy);
- /* Fill the holes in this tile. */
- gal_binary_fill_holes(copy);
+ /* Fill the holes in this tile: holes with maximal connectivity means
+ that they are most strongly bounded. */
+ gal_binary_fill_holes(copy, copy->ndim, -1);
if(fho_prm->step==1)
{
detection_write_in_large(tile, copy);
@@ -404,7 +406,7 @@ detection_sn_write_to_file(struct noisechiselparams *p,
gal_data_t *sn,
/* Description comment. */
str = ( s0d1D2
? ( s0d1D2==2
- ? "S/N of dilated detections."
+ ? "S/N of grown detections."
: "Pseudo-detection S/N over initial detections." )
: "Pseudo-detection S/N over initial undetections.");
gal_list_str_add(&comments, str, 1);
@@ -421,8 +423,8 @@ detection_sn_write_to_file(struct noisechiselparams *p,
gal_data_t *sn,
/* Abort NoiseChisel if the user asked for it. */
if(s0d1D2==2 && !p->continueaftercheck)
ui_abort_after_check(p, p->detsn_s_name, p->detsn_d_name,
- "pseudo-detection and dilated S/N values in "
- "a table");
+ "pseudo-detection and grown/final detection S/N "
+ "values in a table");
}
@@ -742,9 +744,10 @@ detection_pseudo_real(struct noisechiselparams *p)
-/* This is for the final detections (dilated) detections. */
+/* This is for the final detections (grown) detections. */
static size_t
-detection_final_remove_small_sn(struct noisechiselparams *p, size_t num)
+detection_final_remove_small_sn(struct noisechiselparams *p,
+ gal_data_t *workbin, size_t num)
{
size_t i;
int8_t *b;
@@ -755,7 +758,6 @@ detection_final_remove_small_sn(struct noisechiselparams
*p, size_t num)
int32_t *newlabs=gal_data_calloc_array(GAL_TYPE_INT32, num+1, __func__,
"newlabs");
-
/* Get the Signal to noise ratio of all detections. */
sn=detection_sn(p, p->olabel, num, 2, "DILATED");
@@ -767,7 +769,7 @@ detection_final_remove_small_sn(struct noisechiselparams
*p, size_t num)
/* Go over the labeled image and correct the labels. */
- b=p->binary->array;
+ b=workbin->array;
lf=(l=p->olabel->array)+p->olabel->size;
if( p->input->flag & GAL_DATA_FLAG_HASBLANK )
{
@@ -799,8 +801,7 @@ detection_final_remove_small_sn(struct noisechiselparams
*p, size_t num)
/* Make the comments, then write the table. */
gal_list_str_add(&comments, "See also: `DILATED' "
"HDU of output with `--checkdetection'.", 1);
- gal_list_str_add(&comments, "S/N of finally dilated "
- "detections.", 1);
+ gal_list_str_add(&comments, "S/N of finally grown detections.", 1);
threshold_write_sn_table(p, sn, snind, p->detsn_D_name, comments);
@@ -809,15 +810,6 @@ detection_final_remove_small_sn(struct noisechiselparams
*p, size_t num)
}
- /* For a check image. */
- if(p->detectionname)
- {
- p->olabel->name="DETECTION-FINAL";
- gal_fits_img_write(p->olabel, p->detectionname, NULL,
- PROGRAM_NAME);
- p->olabel->name=NULL;
- }
-
/* Clean up and return. */
free(newlabs);
gal_data_free(sn);
@@ -887,13 +879,20 @@ detection_remove_false_initial(struct noisechiselparams
*p,
for(i=0;i<p->numinitialdets;++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.*/
+ /* Replace the byt and olab values with their proper values. If the user
+ doesn't want to grow, 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)
+ if(p->detgrowquant==1.0f) /* We need the binary array even when */
+ do /* there is no growth: the binary */
+ { /* array is used for estimating the */
+ if(*l!=GAL_BLANK_INT32) /* Sky and its STD. */
+ *b = ( *l = newlabels[ *l ] ) > 0;
+ ++b;
+ }
+ while(++l<lf);
+ else
do
{
if(*l!=GAL_BLANK_INT32)
@@ -901,14 +900,6 @@ detection_remove_false_initial(struct noisechiselparams *p,
++b;
}
while(++l<lf);
- else /* We need the binary array even when */
- do /* there is no dilation: the binary */
- { /* array is used for estimating the */
- if(*l!=GAL_BLANK_INT32) /* Sky and its STD. */
- *b = ( *l = newlabels[ *l ] ) > 0;
- ++b;
- }
- while(++l<lf);
/* Clean up and return. */
@@ -920,6 +911,94 @@ detection_remove_false_initial(struct noisechiselparams *p,
+static size_t
+detection_quantile_expand(struct noisechiselparams *p, gal_data_t *workbin)
+{
+ int32_t *o, *of;
+ size_t *d, counter=0, numexpanded;
+ float *i, *e_th, *arr=p->conv->array;
+ gal_data_t *exp_thresh_full, *diffuseindexs;
+ uint8_t *b=workbin->array, *bf=b+workbin->size;
+
+ /* Expand the threshold values (from one value for each tile) to the
+ whole dataset. Since we know the tiles cover the whole image, we don't
+ neeed to initialize or check for blanks.*/
+ exp_thresh_full=gal_tile_block_write_const_value(p->expand_thresh,
+ p->cp.tl.tiles, 0, 0);
+
+ /* Count the pixels that must be expanded. */
+ e_th=exp_thresh_full->array;
+ do { if(*b++==0 && *arr>*e_th) ++counter; ++arr; ++e_th; } while(b<bf);
+
+ /* Allocate the space necessary to keep all the index of all the pixels
+ that must be expanded and re-initialize the necessary pointers. */
+ diffuseindexs=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &counter, NULL, 0,
+ p->cp.minmapsize, NULL, NULL, NULL);
+
+ /* Fill in the diffuse indexs and initialize the objects dataset. */
+ b=workbin->array;
+ arr=p->conv->array;
+ d=diffuseindexs->array;
+ e_th=exp_thresh_full->array;
+ of=(o=p->olabel->array)+p->olabel->size;
+ do
+ {
+ /* If the binary value is 1, then we want an initial label of 1 (the
+ object is already detected). If it isn't, then we only want it if
+ it is above the threshold. */
+ *o = *b==1 ? 1 : ( *arr>*e_th ? CLUMPS_INIT : 0);
+ if(*b==0 && *arr>*e_th)
+ *d++ = o - (int32_t *)(p->olabel->array);
+
+ /* Increment the points and go onto the next pixel. */
+ ++b;
+ ++arr;
+ ++e_th;
+ }
+ while(++o<of);
+
+ /* Expand the detections. */
+ clumps_grow(p->olabel, diffuseindexs, 0);
+
+ /* Only keep the 1 valued pixels in the binary array and fill its
+ holes. */
+ o=p->olabel->array;
+ bf=(b=workbin->array)+workbin->size;
+ do *b = (*o++ == 1); while(++b<bf);
+ workbin=gal_binary_dilate(workbin, 1, 1, 1);
+ gal_binary_fill_holes(workbin, 1, p->detgrowmaxholesize);
+
+ /* Get the labeled image. */
+ numexpanded=gal_binary_connected_components(workbin, &p->olabel, 8);
+
+ /* Set all the input's blank pixels to blank in the labeled and binary
+ arrays. */
+ if( gal_blank_present(p->input, 1) )
+ {
+ b=workbin->array;
+ i=p->input->array;
+ of=(o=p->olabel->array)+p->olabel->size;
+ do
+ {
+ if(isnan(*i++))
+ {
+ *o=GAL_BLANK_INT32;
+ *b=GAL_BLANK_UINT8;
+ }
+ ++b;
+ }
+ while(++o<of);
+ }
+
+ /* Clean up and return */
+ gal_data_free(exp_thresh_full);
+ return numexpanded;
+}
+
+
+
+
+
/* The initial detection has been done, now we want to remove false
detections. */
void
@@ -964,6 +1043,7 @@ detection(struct noisechiselparams *p)
/* Only keep the initial detections that overlap with the real
pseudo-detections. */
+ if(!p->cp.quiet) gettimeofday(&t1, NULL);
num_true_initial=detection_remove_false_initial(p, workbin);
if(!p->cp.quiet)
{
@@ -974,48 +1054,46 @@ detection(struct noisechiselparams *p)
}
- /* If the user asked for dilation, then apply it. */
+ /* If the user asked for dilation/expansion, then apply it and report the
+ final number of detections. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
- if(p->dilate)
- {
- gal_binary_dilate(workbin, p->dilate, p->dilatengb==4 ? 1 : 2 , 1);
- num_true_initial = gal_binary_connected_components(workbin, &p->olabel,
- 8);
- }
+ if(p->detgrowquant!=1.0f)
+ num_true_initial=detection_quantile_expand(p, workbin);
if(!p->cp.quiet)
{
- asprintf(&msg, "%zu detections after %zu dilation%s",
- num_true_initial, p->dilate, p->dilate>1 ? "s." : ".");
+ if(p->detgrowquant==1.0f)
+ asprintf(&msg, "%zu detections with no growth.", num_true_initial);
+ else
+ asprintf(&msg, "%zu detections after growth to %.3f quantile.",
+ num_true_initial, p->detgrowquant);
gal_timing_report(&t1, msg, 2);
free(msg);
}
- /* When the final (dilated or over-all object) detection's S/N is less
- than the pseudo-detection's S/N limit, the object is false. For a real
+ /* When the final (grown or over-all object) detection's S/N is less than
+ the pseudo-detection's S/N limit, the object is false. For a real
detection, the actual object S/N should be higher than any of its
pseudo-detection because it has a much larger area (and possibly more
flux under it). So when the final S/N is smaller than the minimum
acceptable S/N threshold, we have a false pseudo-detection. */
- if(p->cleandilated)
- p->numdetections = detection_final_remove_small_sn(p, num_true_initial);
- else
- {
- p->numdetections = num_true_initial;
- if(p->detectionname)
- {
- p->olabel->name="DETECTION-FINAL";
- gal_fits_img_write(p->olabel, p->detectionname, NULL,
- PROGRAM_NAME);
- p->olabel->name=NULL;
- }
- }
+ p->numdetections = ( p->cleangrowndet
+ ? detection_final_remove_small_sn(p, workbin,
+ num_true_initial)
+ : num_true_initial );
if(!p->cp.quiet)
{
asprintf(&msg, "%zu final true detections.", p->numdetections);
gal_timing_report(&t0, msg, 1);
free(msg);
}
+ if(p->detectionname)
+ {
+ p->olabel->name="DETECTION-FINAL";
+ gal_fits_img_write(p->olabel, p->detectionname, NULL,
+ PROGRAM_NAME);
+ p->olabel->name=NULL;
+ }
/* p->binary was used to keep the initial pseudo-detection threshold. But
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index bdd7396..5d08e59 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -47,7 +47,9 @@ struct noisechiselparams
struct gal_tile_two_layer_params ltl;/* Large tessellation. */
char *inputname; /* Input filename. */
char *kernelname; /* Input kernel filename. */
+ char *widekernelname; /* Name of wider kernel to be used. */
char *khdu; /* Kernel HDU. */
+ char *wkhdu; /* Wide kernel HDU. */
uint8_t skysubtracted; /* Input has been Sky subtracted before. */
float minskyfrac; /* Undetected area min. frac. in tile. */
size_t minnumfalse; /* Min No. of det/seg for true quantile. */
@@ -72,9 +74,9 @@ struct noisechiselparams
size_t detsnminarea; /* Minimum pseudo-detection area for S/N. */
uint8_t checkdetsn; /* Save pseudo-detection S/N values. */
float detquant; /* True detection quantile. */
- size_t dilate; /* Number of times to dilate true dets. */
- size_t dilatengb; /* Connectivity for final dilation. */
- uint8_t cleandilated; /* Remove dilated objects with small S/N. */
+ float detgrowquant; /* Quantile to grow true detections. */
+ size_t detgrowmaxholesize; /* Max. size of holes to fill in growth. */
+ uint8_t cleangrowndet; /* Remove grown objects with small S/N. */
uint8_t checkdetection; /* Save all detection steps to a file. */
uint8_t checksky; /* Check the Sky value estimation. */
@@ -100,11 +102,14 @@ struct noisechiselparams
char *segmentationname; /* Name of segmentation steps file. */
gal_data_t *input; /* Input image. */
- gal_data_t *kernel; /* Input image. */
- gal_data_t *conv; /* Convolved input. */
+ gal_data_t *kernel; /* Sharper kernel. */
+ gal_data_t *widekernel; /* Wider kernel. */
+ gal_data_t *conv; /* Convolved wth sharper kernel. */
+ gal_data_t *wconv; /* Convolved with wider kernel. */
gal_data_t *binary; /* For binary operations. */
gal_data_t *olabel; /* Labels of objects in the detection. */
gal_data_t *clabel; /* Labels of clumps in the detection. */
+ gal_data_t *expand_thresh; /* Quantile threshold to expand per tile. */
gal_data_t *sky; /* Mean of undetected pixels, per tile. */
gal_data_t *std; /* STD of undetected pixels, per tile. */
size_t maxtcontig; /* Maximum contiguous space for a tile. */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index eb9c222..3417f9d 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -58,19 +58,34 @@ noisechisel_convolve(struct noisechiselparams *p)
struct timeval t1;
struct gal_tile_two_layer_params *tl=&p->cp.tl;
- /* Do the convolution. */
+ /* Convovle with sharper kernel. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
p->conv=gal_convolve_spatial(tl->tiles, p->kernel, p->cp.numthreads,
1, tl->workoverch);
- gal_checkset_allocate_copy("CONVOLVED", &p->conv->name);
+ gal_checkset_allocate_copy(p->widekernel?"CONVOLVED-SHARPER":"CONVOLVED",
+ &p->conv->name);
-
- if(!p->cp.quiet) gal_timing_report(&t1, "Convolved with kernel.", 1);
+ /* Report and write check images if necessary. */
+ if(!p->cp.quiet)
+ gal_timing_report(&t1, "Convolved with sharper kernel.", 1);
if(p->detectionname)
{
gal_fits_img_write(p->input, p->detectionname, NULL, PROGRAM_NAME);
gal_fits_img_write(p->conv, p->detectionname, NULL, PROGRAM_NAME);
}
+
+ /* Convolve with wider kernel (if requested). */
+ if(p->widekernel)
+ {
+ if(!p->cp.quiet) gettimeofday(&t1, NULL);
+ p->wconv=gal_convolve_spatial(tl->tiles, p->widekernel,
+ p->cp.numthreads, 1, tl->workoverch);
+ gal_checkset_allocate_copy("CONVOLVED-WIDER", &p->wconv->name);
+
+ /* Report the status: */
+ if(!p->cp.quiet)
+ gal_timing_report(&t1, "Convolved with wider kernel.", 1);
+ }
}
@@ -232,7 +247,7 @@ noisechisel_output(struct noisechiselparams *p)
/* Write the Sky image into the output */
if(p->sky->name) free(p->sky->name);
p->sky->name="SKY";
- gal_tile_full_values_write(p->sky, &p->cp.tl, p->cp.output,
+ gal_tile_full_values_write(p->sky, &p->cp.tl, 1, p->cp.output,
NULL, PROGRAM_NAME);
p->sky->name=NULL;
@@ -248,7 +263,7 @@ noisechisel_output(struct noisechiselparams *p)
gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "MEDSTD", 0, &p->medstd, 0,
"Median raw tile standard deviation", 0,
p->input->unit);
- gal_tile_full_values_write(p->std, &p->cp.tl, p->cp.output, keys,
+ gal_tile_full_values_write(p->std, &p->cp.tl, 1, p->cp.output, keys,
PROGRAM_NAME);
p->std->name=NULL;
}
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index eb3e706..0b6bd65 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -469,7 +469,8 @@ segmentation_on_threads(void *in_prm)
{
/* Grow the true clumps over the detection. */
clumps_grow_prepare_initial(&cltprm);
- if(cltprm.diffuseindexs->size) clumps_grow(&cltprm, 1);
+ if(cltprm.diffuseindexs->size)
+ clumps_grow(p->olabel, cltprm.diffuseindexs, 1);
if(clprm->step==3)
{ gal_data_free(cltprm.diffuseindexs); continue; }
@@ -508,7 +509,7 @@ segmentation_on_threads(void *in_prm)
clumps_grow_prepare_final(&cltprm);
/* Cover the whole area. */
- clumps_grow(&cltprm, 0);
+ clumps_grow(p->olabel, cltprm.diffuseindexs, 0);
}
gal_data_free(cltprm.diffuseindexs);
if(clprm->step==5) { gal_data_free(cltprm.clumptoobj); continue; }
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 5cc9340..23ac8da 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -174,8 +174,10 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
cp->numthreads);
if(checkname)
{
- gal_tile_full_values_write(p->sky, tl, checkname, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(p->std, tl, checkname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(p->sky, tl, 1, checkname, NULL,
+ PROGRAM_NAME);
+ gal_tile_full_values_write(p->std, tl, 1, checkname, NULL,
+ PROGRAM_NAME);
}
/* Get the basic information about the standard deviation
@@ -202,7 +204,7 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
p->cpscorr = p->minstd>1 ? 1.0f : p->minstd;
/* Interpolate and smooth the derived values. */
- threshold_interp_smooth(p, &p->sky, &p->std, checkname);
+ threshold_interp_smooth(p, &p->sky, &p->std, NULL, checkname);
/* If a check was requested, abort NoiseChisel. */
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 1604942..d7da847 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -137,9 +137,10 @@ threshold_apply_on_thread(void *in_prm)
default:
- error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we can
"
- "address the problem. A value of %d had for `taprm->kind' "
- "is not valid", __func__, PACKAGE_BUGREPORT, taprm->kind);
+ error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s so we "
+ "can address the problem. A value of %d had for "
+ "`taprm->kind' is not valid", __func__, PACKAGE_BUGREPORT,
+ taprm->kind);
}
}
@@ -247,7 +248,8 @@ threshold_write_sn_table(struct noisechiselparams *p,
gal_data_t *insn,
/* Interpolate and smooth the values for each tile over the whole image. */
void
threshold_interp_smooth(struct noisechiselparams *p, gal_data_t **first,
- gal_data_t **second, char *filename)
+ gal_data_t **second, gal_data_t **third,
+ char *filename)
{
gal_data_t *tmp;
struct gal_options_common_params *cp=&p->cp;
@@ -255,26 +257,45 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
/* A small sanity check. */
if( (*first)->next )
- error(EXIT_FAILURE, 0, "%s: the `first' argument to must not have any "
- "`next' pointer.", __func__);
+ error(EXIT_FAILURE, 0, "%s: `first' must not have any `next' pointer.",
+ __func__);
+ if( (*second)->next )
+ error(EXIT_FAILURE, 0, "%s: `second' must not have any `next' pointer.",
+ __func__);
+ if( third && (*third)->next )
+ error(EXIT_FAILURE, 0, "%s: `third' must not have any `next' pointer.",
+ __func__);
/* Do the interpolation of both arrays. */
(*first)->next = *second;
+ if(third) (*second)->next = *third;
tmp=gal_interpolate_close_neighbors(*first, tl, cp->interpnumngb,
cp->numthreads, cp->interponlyblank, 1);
gal_data_free(*first);
gal_data_free(*second);
+ if(third) gal_data_free(*third);
*first=tmp;
- *second=tmp->next;
+ *second=(*first)->next;
+ if(third)
+ {
+ *third=(*second)->next;
+ (*third)->next=NULL;
+ }
(*first)->next=(*second)->next=NULL;
if(filename)
{
(*first)->name="THRESH1_INTERP";
(*second)->name="THRESH2_INTERP";
- gal_tile_full_values_write(*first, tl, filename, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(*second, tl, filename, NULL, PROGRAM_NAME);
+ if(third) (*third)->name="THRESH3_INTERP";
+ gal_tile_full_values_write(*first, tl, 1, filename, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(*second, tl, 1, filename, NULL,
+ PROGRAM_NAME);
+ if(third)
+ gal_tile_full_values_write(*third, tl, 1, filename, NULL,
+ PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
+ if(third) (*third)->name=NULL;
}
/* Smooth the threshold if requested. */
@@ -292,16 +313,30 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
gal_data_free(*second);
*second=tmp;
+ /* Smooth the third */
+ if(third)
+ {
+ tmp=gal_tile_full_values_smooth(*third, tl, p->smoothwidth,
+ p->cp.numthreads);
+ gal_data_free(*third);
+ *third=tmp;
+ }
+
/* Add them to the check image. */
if(filename)
{
(*first)->name="THRESH1_SMOOTH";
(*second)->name="THRESH2_SMOOTH";
- gal_tile_full_values_write(*first, tl, filename, NULL,
+ if(third) (*third)->name="THRESH3_SMOOTH";
+ gal_tile_full_values_write(*first, tl, 1, filename, NULL,
PROGRAM_NAME);
- gal_tile_full_values_write(*second, tl, filename, NULL,
+ gal_tile_full_values_write(*second, tl, 1, filename, NULL,
PROGRAM_NAME);
+ if(third)
+ gal_tile_full_values_write(*third, tl, 1, filename, NULL,
+ PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
+ if(third) (*third)->name=NULL;
}
}
}
@@ -332,6 +367,7 @@ struct qthreshparams
{
gal_data_t *erode_th;
gal_data_t *noerode_th;
+ gal_data_t *expand_th;
void *usage;
struct noisechiselparams *p;
};
@@ -350,6 +386,7 @@ qthresh_on_tile(void *in_prm)
double *darr;
void *tarray=NULL;
int type=qprm->erode_th->type;
+ gal_data_t *modeconv = p->wconv ? p->wconv : p->conv;
gal_data_t *tile, *mode, *qvalue, *usage, *tblock=NULL;
size_t i, tind, twidth=gal_type_sizeof(type), ndim=p->input->ndim;
@@ -375,23 +412,18 @@ qthresh_on_tile(void *in_prm)
tile = &p->cp.tl.tiles[tind];
- /* If we have a convolved image, temporarily change the tile's
- pointers so we can do the work on the convolved image, then copy
- the desired contents into the already allocated `usage' array. */
- if(p->conv)
- {
- tarray=tile->array; tblock=tile->block;
- tile->array=gal_tile_block_relative_to_other(tile, p->conv);
- tile->block=p->conv;
- }
+ /* Temporarily change the tile's pointers so we can do the work on
+ the convolved image, then copy the desired contents into the
+ already allocated `usage' array. */
+ tarray=tile->array; tblock=tile->block;
+ tile->array=gal_tile_block_relative_to_other(tile, modeconv);
+ tile->block=modeconv;
gal_data_copy_to_allocated(tile, usage);
- if(p->conv) { tile->array=tarray; tile->block=tblock; }
+ tile->array=tarray; tile->block=tblock;
- /* Find the mode on this dataset, note that we have set the `inplace'
- flag to `1'. So as a byproduct of finding the mode, `usage' is not
- going to have any blank elements and will be sorted (thus ready to
- be used by the quantile functions). */
+ /* Find the mode on this tile, note that we have set the `inplace'
+ flag to `1' to avoid extra allocation. */
mode=gal_statistics_mode(usage, p->mirrordist, 1);
@@ -401,6 +433,24 @@ qthresh_on_tile(void *in_prm)
darr=mode->array;
if( fabs(darr[1]-0.5f) < p->modmedqdiff )
{
+ /* The mode was found on the wider convolved image, but the
+ qthresh values have to be found on the sharper convolved
+ images. This is because the distribution becomes more skewed
+ with a wider kernel, helping us find tiles with no data more
+ easily. But for the quantile threshold, we want to use the
+ sharper convolved image to loose less of the spatial
+ information. */
+ if(modeconv!=p->conv)
+ {
+ tarray=tile->array; tblock=tile->block;
+ tile->array=gal_tile_block_relative_to_other(tile, p->conv);
+ tile->block=p->conv;
+ usage->ndim=ndim; /* Since usage was modified in */
+ usage->size=p->maxtcontig; /* place, it needs to be */
+ gal_data_copy_to_allocated(tile, usage); /* re-initialized. */
+ tile->array=tarray; tile->block=tblock;
+ }
+
/* Get the erosion quantile for this tile and save it. Note that
the type of `qvalue' is the same as the input dataset. */
qvalue=gal_statistics_quantile(usage, p->qthresh, 1);
@@ -408,11 +458,21 @@ qthresh_on_tile(void *in_prm)
qvalue->array, twidth);
gal_data_free(qvalue);
- /* Do the same for the no-erode quantile. */
+ /* Same for the no-erode quantile. */
qvalue=gal_statistics_quantile(usage, p->noerodequant, 1);
memcpy(gal_data_ptr_increment(qprm->noerode_th->array, tind, type),
qvalue->array, twidth);
gal_data_free(qvalue);
+
+ /* Same for the expansion quantile. */
+ if(p->detgrowquant!=1.0f)
+ {
+ qvalue=gal_statistics_quantile(usage, p->detgrowquant, 1);
+ memcpy(gal_data_ptr_increment(qprm->expand_th->array, tind,
+ type),
+ qvalue->array, twidth);
+ gal_data_free(qvalue);
+ }
}
else
{
@@ -420,6 +480,9 @@ qthresh_on_tile(void *in_prm)
tind, type), type);
gal_blank_write(gal_data_ptr_increment(qprm->noerode_th->array,
tind, type), type);
+ if(p->detgrowquant!=1.0f)
+ gal_blank_write(gal_data_ptr_increment(qprm->expand_th->array,
+ tind, type), type);
}
/* Clean up and fix the tile's pointers. */
@@ -456,8 +519,13 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
as the input, making it hard to inspect visually. So we'll only put
the full input when `oneelempertile' isn't requested. */
if(p->qthreshname && !tl->oneelempertile)
- gal_fits_img_write(p->conv ? p->conv : p->input, p->qthreshname, NULL,
- PROGRAM_NAME);
+ {
+ gal_fits_img_write(p->conv ? p->conv : p->input, p->qthreshname, NULL,
+ PROGRAM_NAME);
+ if(p->wconv)
+ gal_fits_img_write(p->wconv ? p->wconv : p->input, p->qthreshname,
+ NULL, PROGRAM_NAME);
+ }
/* Allocate space for the quantile threshold values. */
@@ -467,6 +535,12 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
qprm.noerode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
NULL, p->input->unit, NULL);
+ if(p->detgrowquant!=1.0f)
+ qprm.expand_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
+ tl->numtiles, NULL, 0, cp->minmapsize,
+ NULL, p->input->unit, NULL);
+ else
+ qprm.expand_th=NULL;
/* Allocate temporary space for processing in each tile. */
@@ -483,22 +557,35 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
gal_threads_spin_off(qthresh_on_tile, &qprm, tl->tottiles, cp->numthreads);
free(qprm.usage);
if( gal_blank_present(qprm.erode_th, 1) )
- qprm.noerode_th->flag |= GAL_DATA_FLAG_HASBLANK;
+ {
+ qprm.noerode_th->flag |= GAL_DATA_FLAG_HASBLANK;
+ if(qprm.expand_th) qprm.expand_th->flag |= GAL_DATA_FLAG_HASBLANK;
+ }
qprm.noerode_th->flag |= GAL_DATA_FLAG_BLANK_CH;
+ if(p->detgrowquant!=1.0f) qprm.expand_th->flag |= GAL_DATA_FLAG_BLANK_CH;
if(p->qthreshname)
{
qprm.erode_th->name="QTHRESH_ERODE";
qprm.noerode_th->name="QTHRESH_NOERODE";
- gal_tile_full_values_write(qprm.erode_th, tl, p->qthreshname, NULL,
+ gal_tile_full_values_write(qprm.erode_th, tl, 1, p->qthreshname, NULL,
PROGRAM_NAME);
- gal_tile_full_values_write(qprm.noerode_th, tl, p->qthreshname, NULL,
+ gal_tile_full_values_write(qprm.noerode_th, tl, 1, p->qthreshname, NULL,
PROGRAM_NAME);
- qprm.erode_th->name = qprm.noerode_th->name = NULL;
+ qprm.erode_th->name=qprm.noerode_th->name=NULL;
+
+ if(qprm.expand_th)
+ {
+ qprm.expand_th->name="QTHRESH_EXPAND";
+ gal_tile_full_values_write(qprm.expand_th, tl, 1, p->qthreshname,
+ NULL, PROGRAM_NAME);
+ qprm.expand_th->name=NULL;
+ }
}
/* Interpolate and smooth the derived values. */
threshold_interp_smooth(p, &qprm.erode_th, &qprm.noerode_th,
+ qprm.expand_th ? & qprm.expand_th : NULL,
p->qthreshname);
@@ -512,6 +599,10 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
gal_fits_img_write(p->binary, p->qthreshname, NULL, PROGRAM_NAME);
+ /* Set the expansion quantile if necessary. */
+ p->expand_thresh = qprm.expand_th ? qprm.expand_th : NULL;
+
+
/* Clean up and report duration if necessary. */
gal_data_free(qprm.erode_th);
gal_data_free(qprm.noerode_th);
diff --git a/bin/noisechisel/threshold.h b/bin/noisechisel/threshold.h
index 75f88c0..595fa11 100644
--- a/bin/noisechisel/threshold.h
+++ b/bin/noisechisel/threshold.h
@@ -44,7 +44,8 @@ threshold_write_sn_table(struct noisechiselparams *p,
gal_data_t *sntable,
void
threshold_interp_smooth(struct noisechiselparams *p, gal_data_t **first,
- gal_data_t **second, char *filename);
+ gal_data_t **second, gal_data_t **third,
+ char *filename);
void
threshold_quantile_find_apply(struct noisechiselparams *p);
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 960eaac..09f0d9d 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -119,7 +119,6 @@ ui_initialize_options(struct noisechiselparams *p,
cp->numthreads = gal_threads_number();
cp->coptions = gal_commonopts_options;
-
/* Modify common options. */
for(i=0; !gal_options_is_last(&cp->coptions[i]); ++i)
{
@@ -228,9 +227,6 @@ ui_read_check_only_options(struct noisechiselparams *p)
if(p->openingngb!=4 && p->openingngb!=8)
error(EXIT_FAILURE, 0, "%zu not acceptable for `--openingngb'. It must "
"be 4 or 8 (specifying the type of connectivity)", p->openingngb);
- if(p->dilatengb!=4 && p->dilatengb!=8)
- error(EXIT_FAILURE, 0, "%zu not acceptable for `--dilatengb'. It must "
- "be 4 or 8 (specifying the type of connectivity)", p->dilatengb);
/* Make sure that the no-erode-quantile is not smaller or equal to
qthresh. */
@@ -249,6 +245,36 @@ ui_read_check_only_options(struct noisechiselparams *p)
"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\"");
+
+ /* Kernel checks. */
+ if(p->kernelname)
+ {
+ /* Check if it exists. */
+ gal_checkset_check_file(p->kernelname);
+
+ /* If its FITS, see if a HDU has been provided. */
+ if( gal_fits_name_is_fits(p->kernelname) && p->khdu==NULL )
+ error(EXIT_FAILURE, 0, "no HDU specified for kernel. When the "
+ "kernel is a FITS file, a HDU must also be specified. You "
+ "can use the `--khdu' option and give it the HDU number "
+ "(starting from zero), extension name, or anything "
+ "acceptable by CFITSIO");
+ }
+
+ /* Wide kernel checks. */
+ if(p->widekernelname)
+ {
+ /* Check if it exists. */
+ gal_checkset_check_file(p->widekernelname);
+
+ /* If its FITS, see if a HDU has been provided. */
+ if( gal_fits_name_is_fits(p->widekernelname) && p->wkhdu==NULL )
+ error(EXIT_FAILURE, 0, "no HDU specified for wide kernel. When the "
+ "wide kernel is a FITS file, a HDU must also be specified. You "
+ "can use the `--khdu' option and give it the HDU number "
+ "(starting from zero), extension name, or anything "
+ "acceptable by CFITSIO");
+ }
}
@@ -274,21 +300,6 @@ ui_check_options_and_arguments(struct noisechiselparams *p)
}
else
error(EXIT_FAILURE, 0, "no input file is specified");
-
- /* Basic Kernel file checks. */
- if(p->kernelname)
- {
- /* Check if it exists. */
- gal_checkset_check_file(p->kernelname);
-
- /* If its FITS, see if a HDU has been provided. */
- if( gal_fits_name_is_fits(p->kernelname) && p->khdu==NULL )
- error(EXIT_FAILURE, 0, "no HDU specified for kernel. When the "
- "kernel is a FITS file, a HDU must also be specified, you "
- "can use the `--khdu' option and give it the HDU number "
- "(starting from zero), extension name, or anything "
- "acceptable by CFITSIO");
- }
}
@@ -360,7 +371,7 @@ ui_set_output_names(struct noisechiselparams *p)
? "_detsn_det.txt" : "_detsn_det.fits") );
p->detsn_D_name=gal_checkset_automatic_output(&p->cp, basename,
( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
- ? "_detsn_dilated.txt" : "_detsn_dilated.fits") );
+ ? "_detsn_grown.txt" : "_detsn_grown.fits") );
}
/* Detection steps. */
@@ -468,6 +479,13 @@ ui_prepare_kernel(struct noisechiselparams *p)
ff=(f=kernel_2d)+gal_dimension_total_size(2, p->kernel->dsize);
do *k++=*f; while(++f<ff);
}
+
+
+ /* If a wide kernel is given, then read it into memory. Otherwise, just
+ ignore it. */
+ if(p->widekernelname)
+ p->widekernel=gal_fits_img_read_kernel(p->widekernelname, p->wkhdu,
+ p->cp.minmapsize);
}
@@ -576,9 +594,7 @@ ui_preparations(struct noisechiselparams *p)
p->input->ndim);
- /* Check for blank values to help later processing. AFTERWARDS, set the
- USE_ZERO flag, so the zero-bit (if the input doesn't have any blank
- value) will be meaningful. */
+ /* Check for blank values to help later processing. */
gal_blank_present(p->input, 1);
@@ -682,9 +698,15 @@ ui_read_check_inputs_setup(int argc, char *argv[],
p->cp.numthreads==1 ? "." : "s.");
printf(" - Input: %s (hdu: %s)\n", p->inputname, p->cp.hdu);
if(p->kernelname)
- printf(" - Kernel: %s (hdu: %s)\n", p->kernelname, p->khdu);
+ printf(" - %s: %s (hdu: %s)\n",
+ p->widekernelname ? "Sharp Kernel" : "Kernel",
+ p->kernelname, p->khdu);
else
- printf(" - Kernel: FWHM=2 pixel Gaussian.\n");
+ printf(" - %s: FWHM=2 pixel Gaussian.\n",
+ p->widekernelname ? "Sharp Kernel" : "Kernel");
+ if(p->widekernelname)
+ printf(" - Wide Kernel: %s (hdu: %s)\n", p->widekernelname,
+ p->wkhdu);
}
}
@@ -764,11 +786,13 @@ ui_free_report(struct noisechiselparams *p, struct
timeval *t1)
gal_data_free(p->sky);
gal_data_free(p->std);
gal_data_free(p->conv);
+ gal_data_free(p->wconv);
gal_data_free(p->input);
gal_data_free(p->kernel);
gal_data_free(p->binary);
gal_data_free(p->olabel);
gal_data_free(p->clabel);
+ gal_data_free(p->widekernel);
/* Clean up the tile structure. */
p->ltl.numchannels=NULL;
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 416cd76..bfe3b22 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -29,7 +29,7 @@ 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 b f j l n u x z
A H J W X Y
*/
enum option_keys_enum
@@ -37,6 +37,7 @@ enum option_keys_enum
/* With short-option version. */
UI_KEY_LARGETILESIZE = 'L',
UI_KEY_KERNEL = 'k',
+ UI_KEY_WIDEKERNEL = 'w',
UI_KEY_SKYSUBTRACTED = 'E',
UI_KEY_MINSKYFRAC = 'B',
UI_KEY_MIRRORDIST = 'r',
@@ -48,7 +49,7 @@ enum option_keys_enum
UI_KEY_DTHRESH = 'R',
UI_KEY_DETSNMINAREA = 'i',
UI_KEY_DETQUANT = 'c',
- UI_KEY_DILATE = 'd',
+ UI_KEY_DETGROWQUANT = 'd',
UI_KEY_SEGSNMINAREA = 'm',
UI_KEY_SEGQUANT = 'g',
UI_KEY_KEEPMAXNEARRIVER = 'v',
@@ -61,6 +62,7 @@ enum option_keys_enum
/* Only with long version (start with a value 1000, the rest will be set
automatically). */
UI_KEY_KHDU = 1000,
+ UI_KEY_WKHDU,
UI_KEY_MINNUMFALSE,
UI_KEY_ONLYDETECTION,
UI_KEY_GROWNCLUMPS,
@@ -71,8 +73,8 @@ enum option_keys_enum
UI_KEY_OPENINGNGB,
UI_KEY_CHECKDETSKY,
UI_KEY_CHECKDETSN,
- UI_KEY_DILATENGB,
- UI_KEY_CLEANDILATED,
+ UI_KEY_DETGROWMAXHOLESIZE,
+ UI_KEY_CLEANGROWNDET,
UI_KEY_CHECKDETECTION,
UI_KEY_CHECKSKY,
UI_KEY_CHECKCLUMPSN,
diff --git a/bin/statistics/sky.c b/bin/statistics/sky.c
index f68678c..098338f 100644
--- a/bin/statistics/sky.c
+++ b/bin/statistics/sky.c
@@ -192,9 +192,9 @@ sky(struct statisticsparams *p)
}
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->sky_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->std_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
}
@@ -213,9 +213,9 @@ sky(struct statisticsparams *p)
gal_timing_report(&t1, "All blank tiles filled (interplated).", 1);
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->sky_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->std_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
}
@@ -237,9 +237,9 @@ sky(struct statisticsparams *p)
1);
if(p->checksky)
{
- gal_tile_full_values_write(p->sky_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->sky_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, p->checkskyname, NULL,
+ gal_tile_full_values_write(p->std_t, tl, 1, p->checkskyname, NULL,
PROGRAM_NAME);
}
}
@@ -250,8 +250,8 @@ sky(struct statisticsparams *p)
( p->cp.output
? p->cp.output
: p->inputname ), "_sky.fits");
- gal_tile_full_values_write(p->sky_t, tl, outname, NULL, PROGRAM_NAME);
- gal_tile_full_values_write(p->std_t, tl, outname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(p->sky_t, tl, 1, outname, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(p->std_t, tl, 1, outname, NULL, PROGRAM_NAME);
if(!cp->quiet)
printf(" - Written to `%s'.\n", outname);
diff --git a/bin/statistics/statistics.c b/bin/statistics/statistics.c
index 7b27985..e4a2e2b 100644
--- a/bin/statistics/statistics.c
+++ b/bin/statistics/statistics.c
@@ -254,7 +254,7 @@ statistics_interpolate_and_write(struct statisticsparams *p,
}
/* Write the values. */
- gal_tile_full_values_write(values, &cp->tl, output, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(values, &cp->tl, 1, output, NULL, PROGRAM_NAME);
}
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index cdde8ae..76b94dd 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12191,6 +12191,23 @@ major changes since that paper was published.
@itemize
@item
address@hidden: NoiseChisel uses the difference between the mode and
+median to identify if a tile should be used for estimating the quantile
+thresholds (see @ref{Quantifying signal in a tile}). Until now, NoiseChisel
+would convolve an image once and estimate the proper tiles for quantile
+estimations on the convolved image. The same convolved image would later be
+used for quantile estimation. A larger kernel does increase the skewness
+(and thus difference between the mode and median), however, it disfigures
+the shapes/morphology of the objects.
+
+This new @option{--widekernel} option (and a corresponding @option{--wkhdu}
+option to specify its HDU) option are added to solve such cases. When its
+given, the input will be convolved with both the sharp (given through the
address@hidden option) and wide kernels. The mode and median are
+calculated on the dataset that is convolved with the wider kernel, then the
+quantiles are estimated on the image convolved with the sharper kernel.
+
address@hidden
@option{--noerodequant}: to specify a quantile threshold where erosion
will not apply. This is useful to detect sharper point-like sources that
will be missed due to too much erosion. To completely ignore this features
@@ -12198,29 +12215,23 @@ give this option a value of 1 (only the largest
valued pixel in the input
will not be eroded).
@item
address@hidden: to specify the final dilation connectivity. Prior to
-Gnuastro 0.4, it was assumed to be 8 (maximum connectivity in a 2D image
-where each pixel has 8 neighbors).
address@hidden: is used to grow the final true detections until a
+given quantile in the same way that clumps are grown during segmentation
+(compare columns 2 and 3 in Figure 10 of the paper). It replaces the old
address@hidden option in the paper and older versions of
+Gnuastro. Dilation is a blind growth method which causes objects to be boxy
+or diamond shaped when too many layers are added. However, with the growth
+method that is defined now, we can follow the signal into the noise with
+any shape. The appropriate quantile depends on your dataset's correlated
+noise properties and how cleanly it was Sky subtracted. The new
address@hidden can also be used to specify the maximum hole
+size to fill as part of this growth, see the description in @ref{Detection
+options} for more details.
@item
address@hidden: After dilation, if the signal-to-noise ratio of a
-detection is less than the derived pseudo-detection S/N limit, that
-detection will be discarded. In an ideal/clean noise, a true detection's
-S/N should be larger than its constituent pseudo-detections because its
-area is larger and it also covers more signal. However, on a false
-detections (especially at lower @option{--detquant} values), the increase
-in size can cause a decrease in S/N below that threshold.
-
-This will improve purity and not change completeness (a true detection will
-not be discarded). Because a true detection has flux in its vicinity and
-dilation will catch more of that flux and increase the S/N. So on a true
-detection, the final S/N cannot be less than pseudo-detections.
-
-However, in many real images bad processing creates artifacts that cannot
-be accurately removed by the Sky subtraction. In such cases, this option
-will decrease the completeness (will artificially discard true
-detections). So this feature is not default and should to be explicitly
-called when you know the noise is clean.
address@hidden: A process to further clean/remove the possibility
+of false detections, see the descriptions under this option in
address@hidden options}.
@end itemize
@@ -12360,6 +12371,23 @@ Akhlaghi and Ichikawa [2015].
HDU containing the kernel in the file given to the @option{--kernel}
option.
address@hidden -w STR
address@hidden --widekernel=STR
+File name of a wider kernel to use in estimating the difference of the mode
+and median in a tile (this difference is used to identify the significance
+of signal in that tile, see @ref{Quantifying signal in a tile}). This
+convolved image is only used for this purpose. Once the mode is found to be
+sufficiently close to the median, the input convolved with the sharper
+kernel (@option{--kernel}) is used to identify the quantile threshold (see
address@hidden).
+
+Since this convolution will significantly slow down the processing, it is
+optional. If no image is given, the mode and median will be found on the
+image that is convolved with the dataset in @option{--kernel}.
+
address@hidden --wkhdu=STR
+HDU containing the kernel file given to the @option{--widekernel} option.
+
@item -E
@itemx --skysubtracted
If this option is called, it is assumed that the image has already been sky
@@ -12663,18 +12691,63 @@ that this is only calculated for the large mesh grids
that satisfy the
minimum fraction of undetected pixels (value of @option{--minbfrac}) and
minimum number of psudo-detections (value of @option{--minnumfalse}).
address@hidden -d INT
address@hidden --dilate=INT
-Number of times to dilate the final true detections. See the explanations
-in @option{--opening} for more information on dilation. The structuring
-element for this final dilation is fixed to an 8-connected
-neighborhood. This is because astronomical objects, except cosmic rays,
-never have a clear cutoff, so all the 8-pixels connected to the border
-pixels of a detection might harbor data.
address@hidden -d FLT
address@hidden --detgrowquant=FLT
+Quantile limit to ``grow'' the final detections. As discussed in the
+previous options, after applying the initial quantile threshold, layers of
+pixels are carved off the objects to identify true signal. With this step
+you can return those low surface brightness layers that were carved off
+back to the detections. To disable growth, set the value of this option to
address@hidden
+
+The process is as follows: after the true detections are found, all the
+non-detected pixels above this quantile will be put in a list and used to
+``grow'' the true detections (seeds of the growth). Like all quantile
+thresholds, this threshold is defined and applied to the convolved
+dataset. Afterwards, the dataset is dilated once (with minimum
+connectivity) to connect very thin regions on the boundary: imagine
+building a dam at the point rivers spill into an open sea/ocean. Finally,
+all holes are filled. In the geography metaphor, holes can be seen as the
+closed (by the dams) rivers and lakes, so this process is like turning the
+water in all such rivers and lakes into soil. See
address@hidden for configuring the hole filling.
+
address@hidden --detgrowmaxholesize=INT
+The maximum hole size to fill during the final expansion of the true
+detections as described in @option{--detgrowquant}. This is necessary when
+the input contains many smaller objects and can be used to avoid marking
+blank sky regions as detections. For example multiple galaxies can be
+positioned such that they surround an empty region of sky. If all the holes
+are filled, the Sky region in between them will be taken as a detection
+which is not desired. To avoid such cases, the integer given to this option
+must be smaller than the hole between the objects.
+
+On the other hand, if you have a very large (and extended) galaxy, the
+diffuse wings of the galaxy may create very large holes. In such cases, a
+larger value to this option will cause the whole region to be detected as
+part of the large galaxy and thus detect it to extremely low surface
+brightness limits.
+
address@hidden --cleangrowndet
+After dilation, if the signal-to-noise ratio of a detection is less than
+the derived pseudo-detection S/N limit, that detection will be
+discarded. In an ideal/clean noise, a true detection's S/N should be larger
+than its constituent pseudo-detections because its area is larger and it
+also covers more signal. However, on a false detections (especially at
+lower @option{--detquant} values), the increase in size can cause a
+decrease in S/N below that threshold.
+
+This will improve purity and not change completeness (a true detection will
+not be discarded). Because a true detection has flux in its vicinity and
+dilation will catch more of that flux and increase the S/N. So on a true
+detection, the final S/N cannot be less than pseudo-detections.
+
+However, in many real images bad processing creates artifacts that cannot
+be accurately removed by the Sky subtraction. In such cases, this option
+will decrease the completeness (will artificially discard true
+detections). So this feature is not default and should to be explicitly
+called when you know the noise is clean.
address@hidden --dilatengb=INT
-The connectivity used for the final dilation, see @option{--erodengb} for
-more information about connectivity or a structuring element.
@item --checkdetection
Every step of the detection process will be added as an extension to a file
@@ -18109,7 +18182,7 @@ to manually copy and paste the function name or worry
about it changing
later (@code{__func__} was standardized in C99).
@end deftypefun
address@hidden {void *} gal_data_calloc_array (uint8_t @code{type}, size_t
@code{size)}, const char @code{*funcname}, const char @code{*varname})
address@hidden {void *} gal_data_calloc_array (uint8_t @code{type}, size_t
@code{size}, const char @code{*funcname}, const char @code{*varname})
Similar to @code{gal_data_malloc_array}, but the space is cleared (set to
0) after allocation.
@end deftypefun
@@ -20576,7 +20649,7 @@ function, it should be @code{1}. In subsequent calls
(while you are parsing
a tile), it should be increased by one.
@end deftypefun
address@hidden {gal_data_t *} gal_tile_block_write_const_value (gal_data_t
@code{*tilevalues}, gal_data_t @code{*tilesll}, int @code{initialize})
address@hidden {gal_data_t *} gal_tile_block_write_const_value (gal_data_t
@code{*tilevalues}, gal_data_t @code{*tilesll}, int @code{withblank}, int
@code{initialize})
Write a constant value for each tile over the area it covers in an
allocated dataset that is the size of @code{tile}'s allocated block of
memory (found through @code{gal_tile_block} described above). The arguments
@@ -20597,6 +20670,11 @@ doesn't care, it will parse the @code{next} elements
to go to the next
tile. This function will not pop-from or free the @code{tilesll}, it will
only parse it from start to end.
address@hidden withblank
+If the block containing the tiles has blank elements, those blank elements
+will be blank in the output of this function also, hence the array will be
+initialzied with blank values when this option is called (see below).
+
@item initialize
Initialize the allocated space with blank values before writing in the
constant values. This can be useful when the tiles don't cover the full
@@ -20612,7 +20690,7 @@ values. The output dataset will have a type of
@code{GAL_TYPE_INT32} (see
This function can be used when you want to check the coverage of each tile
over the allocated block of memory. It is just a wrapper over the
address@hidden
address@hidden (with @code{withblank} set to zero).
@end deftypefun
@deftypefun {void *} gal_tile_block_relative_to_other (gal_data_t
@code{*tile}, gal_data_t @code{*other})
@@ -20983,13 +21061,17 @@ inverse: IN_ALL[ perm[i] ] = IN_MEMORY[ i
]
@end example
@end deftypefun
address@hidden void gal_tile_full_values_write (gal_data_t @code{*tilevalues},
struct gal_tile_two_layer_params @code{*tl}, char @code{*filename},
gal_fits_list_key_t @code{*keys}, char @code{*program_string})
address@hidden void gal_tile_full_values_write (gal_data_t @code{*tilevalues},
struct gal_tile_two_layer_params @code{*tl}, int @code{withblank}, char
@code{*filename}, gal_fits_list_key_t @code{*keys}, char @code{*program_string})
Write one value for each tile into a file. It is important to note that the
values in @code{tilevalues} must be ordered in the same manner as the
tiles, so @code{tilevalues->array[i]} is the value that should be given to
@code{tl->tiles[i]}. The @code{tl->permutation} array must have been
initialized before calling this function with
@code{gal_tile_full_permutation}.
+
+If @code{withblank} is non-zero, then block structure of the tiles will be
+checked and all blank pixels in the block will be blank in the final output
+file also.
@end deftypefun
@deftypefun {gal_data_t *} gal_tile_full_values_smooth (gal_data_t
@code{*tilevalues}, struct gal_tile_two_layer_params @code{*tl}, size_t
@code{width}, size_t @code{numthreads})
@@ -21747,10 +21829,11 @@ Although the adjacency matrix as used here is
symmetric, currently this
function assumes that it is filled on both sides of the diagonal.
@end deftypefun
address@hidden void gal_binary_fill_holes (gal_data_t @code{*input})
-Fill all the holes that are bounded within a 4-connected region of the
-binary @code{input} dataset. This function currently only works on a 2D
-dataset.
address@hidden void gal_binary_fill_holes (gal_data_t @code{*input}, int
@code{connectivity})
+Fill all the holes (0 valued pixels surrounded by 1 valued pixels) within a
+region of the binary @code{input} dataset. The connectivity of the holes
+can be set with @code{connectivity}. This function currently only works on
+a 2D dataset.
@end deftypefun
@node Convolution functions, Interpolation, Binary datasets, Gnuastro library
diff --git a/lib/binary.c b/lib/binary.c
index a1adfaf..4e5ba3e 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -368,8 +368,8 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
/* Make sure the given dataset has the same size as the input. */
if( gal_data_dsize_is_different(binary, lab) )
- error(EXIT_FAILURE, 0, "%s: the `binary' and `out' datasets must have "
- "the same size", __func__);
+ error(EXIT_FAILURE, 0, "%s: the `binary' and `out' datasets must "
+ "have the same size", __func__);
/* Make sure it has a `int32' type. */
if( lab->type!=GAL_TYPE_INT32 )
@@ -618,8 +618,7 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t
**outtile)
-/* Fill all the holes in an input unsigned char array that are bounded
- within a 4-connected region.
+/* Fill all the holes in an input unsigned char array.
The basic method is this:
@@ -642,9 +641,11 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t
**outtile)
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)
+gal_binary_fill_holes(gal_data_t *input, int connectivity, size_t maxsize)
{
uint8_t *in;
+ uint32_t *i, *fi;
+ size_t numholes, *sizes;
gal_data_t *inv, *tile, *holelabs=NULL;
/* A small sanity check. */
@@ -658,8 +659,8 @@ gal_binary_fill_holes(gal_data_t *input)
inv=binary_make_padded_inverse(input, &tile);
- /* Label the 8-connected (connectivity==2) holes. */
- gal_binary_connected_components(inv, &holelabs, 2);
+ /* Label the holes */
+ numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
/* Any pixel with a label larger than 1 is a hole in the input image and
@@ -669,6 +670,22 @@ gal_binary_fill_holes(gal_data_t *input)
tile->array=gal_tile_block_relative_to_other(tile, holelabs);
tile->block=holelabs; /* has to be after correcting `tile->array'. */
+ /* If the user wants to only fill holes to a certain size, then remove
+ those with a larger size. */
+ if(maxsize<-1)
+ {
+ /* Allocate space to keep the size of each hole: */
+ sizes=gal_data_calloc_array(GAL_TYPE_SIZE_T, numholes+1, __func__,
+ "sizes");
+ fi=(i=holelabs->array)+holelabs->size; do ++sizes[*i]; while(++i<fi);
+
+ /* Set those labels with a larger size to 1 (treat it as
+ background). */
+ fi=(i=holelabs->array)+holelabs->size;
+ do
+ if(*i!=GAL_BLANK_INT32) *i = sizes[*i]>maxsize ? 1 : *i;
+ while(++i<fi);
+ }
/* The type of the tile is already known (it is `int32_t') and we have no
output, so we'll just put `int' as a place-holder. In this way we can
diff --git a/lib/blank.c b/lib/blank.c
index 44dee32..b3550a9 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -168,10 +168,10 @@ gal_blank_present(gal_data_t *input, int updateflag)
blank is present. */
if(input->size==0) return 0;
- /* From the user's flags, it might not be necessary to go through the
- dataset any more. */
- if( ( input->flag & GAL_DATA_FLAG_HASBLANK )
- || ( input->flag & GAL_DATA_FLAG_BLANK_CH ) )
+ /* From the user's flags, you can tell if the dataset has already been
+ checked for blank values or not. If it has, then just return the
+ checked result. */
+ if( input->flag & GAL_DATA_FLAG_BLANK_CH )
return input->flag & GAL_DATA_FLAG_HASBLANK;
/* Go over the pixels and check: */
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index cb19650..1f25915 100644
--- a/lib/gnuastro/binary.h
+++ b/lib/gnuastro/binary.h
@@ -98,7 +98,7 @@ gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
/***************** Fill holes ********************/
/*********************************************************************/
void
-gal_binary_fill_holes(gal_data_t *input);
+gal_binary_fill_holes(gal_data_t *input, int connectivity, size_t maxsize);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index e33f1d8..02ae328 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -85,7 +85,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
gal_data_t *
gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll,
- int initialize);
+ int withblank, int initialize);
gal_data_t *
gal_tile_block_check_tiles(gal_data_t *tiles);
@@ -151,8 +151,8 @@ gal_tile_full_permutation(struct gal_tile_two_layer_params
*tl);
void
gal_tile_full_values_write(gal_data_t *tilevalues,
struct gal_tile_two_layer_params *tl,
- char *filename, gal_fits_list_key_t *keys,
- char *program_string);
+ int withblank, char *filename,
+ gal_fits_list_key_t *keys, char *program_string);
gal_data_t *
gal_tile_full_values_smooth(gal_data_t *tilevalues,
diff --git a/lib/statistics.c b/lib/statistics.c
index 7e77964..6d5e013 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -58,8 +58,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
gal_data_t *
gal_statistics_number(gal_data_t *input)
{
- uint64_t counter=0;
size_t dsize=1;
+ uint64_t counter=0;
gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_UINT64, 1, &dsize,
NULL, 1, -1, NULL, NULL, NULL);
diff --git a/lib/tile.c b/lib/tile.c
index cbff41e..05c8eac 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -423,7 +423,7 @@ gal_tile_block_increment(gal_data_t *block, size_t *tsize,
don't cover the full allocated block. */
gal_data_t *
gal_tile_block_write_const_value(gal_data_t *tilevalues, gal_data_t *tilesll,
- int initialize)
+ int withblank, int initialize)
{
void *in;
int type=tilevalues->type;
@@ -443,7 +443,7 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
/* If requested, initialize `tofill', otherwise it is assumed that the
full area of the output is covered by the tiles. */
- if(initialize) gal_blank_initialize(tofill);
+ if(withblank || initialize) gal_blank_initialize(tofill);
else
{
/* Copy the flags. */
@@ -465,9 +465,9 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
{
/* Set the pointer to use as input. The `if(o)' statement is set
because GCC 7.1.1 complained about the possiblity of the first
- argument of `memcpy' being NULL. */
+ argument of `memcpy' being NULL. Recall that `o' is a pointer. */
in=gal_data_ptr_increment(tilevalues->array, tile_ind++, type);;
- GAL_TILE_PARSE_OPERATE( tile, tofill, 1, 0, {
+ GAL_TILE_PARSE_OPERATE( tile, tofill, 1, withblank, {
if(o) memcpy(o, in, gal_type_sizeof(type));
} );
}
@@ -497,7 +497,7 @@ gal_tile_block_check_tiles(gal_data_t *tilesll)
arr=ids->array; for(i=0;i<dsize;++i) arr[i]=i;
/* Make the output. */
- out=gal_tile_block_write_const_value(ids, tilesll, 1);
+ out=gal_tile_block_write_const_value(ids, tilesll, 0, 1);
/* Clean up and return. */
gal_data_free(ids);
@@ -1088,8 +1088,8 @@ gal_tile_full_permutation(struct
gal_tile_two_layer_params *tl)
void
gal_tile_full_values_write(gal_data_t *tilevalues,
struct gal_tile_two_layer_params *tl,
- char *filename, gal_fits_list_key_t *keys,
- char *program_string)
+ int withblank, char *filename,
+ gal_fits_list_key_t *keys, char *program_string)
{
gal_data_t *disp;
@@ -1113,8 +1113,8 @@ gal_tile_full_values_write(gal_data_t *tilevalues,
else disp = tilevalues;
}
else
- disp=gal_tile_block_write_const_value(tilevalues, tl->tiles, 0);
-
+ disp=gal_tile_block_write_const_value(tilevalues, tl->tiles,
+ withblank, 0);
/* Write the array as a file and then clean up (if necessary). */
gal_fits_img_write(disp, filename, keys, program_string);
diff --git a/tests/noisechisel/noisechisel.sh b/tests/noisechisel/noisechisel.sh
index 14578c4..c60aebf 100755
--- a/tests/noisechisel/noisechisel.sh
+++ b/tests/noisechisel/noisechisel.sh
@@ -49,5 +49,5 @@ if [ ! -f $img ]; then echo "$img does not exist.";
exit 77; fi
# Actual test script
# ==================
-$execname $img --cleandilated --checkdetection --checksegmentation \
- --continueaftercheck
+$execname $img --cleangrowndet --checkdetection --checksegmentation \
+ --continueaftercheck --tilesize=100,100