[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 76eae52 1/3: Dilation of true detections repla
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 76eae52 1/3: Dilation of true detections replaced with growth |
|
Date: |
Thu, 14 Sep 2017 19:27:22 -0400 (EDT) |
branch: master
commit 76eae5284ecc89715dd475e8a6ad6f537655fb38
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Dilation of true detections replaced with growth
Until now, dilation was used to grow the true detections (because we
initially carved layers of pixels off). But dilation is a blind process, it
doesn't follow the signal, it is just based on the shape of the initial
detected region. As a result to go down to low surface brightness, the
objects would get boxy (if 8-connectivity was used) or diamond shaped (if
4-connectivity was used).
With this commit, a process like the growth of clumps (to define objects)
is used. It is indeed slower than dilation, but it is able to identify the
shape of very low surface brightness signal accurately. See the manual for
a more detailed discussion.
In the libraries, it is now possible to define the connectivity of the
holes to fill in the `gal_binary_fill_holes' function.
---
bin/noisechisel/args.h | 35 +++------
bin/noisechisel/astnoisechisel.conf | 3 +-
bin/noisechisel/clumps.c | 35 ++++++---
bin/noisechisel/clumps.h | 2 +-
bin/noisechisel/detection.c | 138 +++++++++++++++++++++++++++---------
bin/noisechisel/main.h | 6 +-
bin/noisechisel/segmentation.c | 5 +-
bin/noisechisel/sky.c | 2 +-
bin/noisechisel/threshold.c | 86 +++++++++++++++++++---
bin/noisechisel/threshold.h | 3 +-
bin/noisechisel/ui.c | 6 +-
bin/noisechisel/ui.h | 5 +-
doc/gnuastro.texi | 92 ++++++++++++++----------
lib/binary.c | 13 ++--
lib/gnuastro/binary.h | 2 +-
tests/noisechisel/noisechisel.sh | 4 +-
16 files changed, 294 insertions(+), 143 deletions(-)
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index 772b09a..dddb8a6 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -365,39 +365,26 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "dilate",
- UI_KEY_DILATE,
- "INT",
- 0,
- "Number of times to dilate true detections.",
- ARGS_GROUP_DETECTION,
- &p->dilate,
- GAL_TYPE_SIZE_T,
- GAL_OPTIONS_RANGE_ANY,
- GAL_OPTIONS_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
- "dilatengb",
- UI_KEY_DILATENGB,
- "INT",
+ "detgrowquant",
+ UI_KEY_DETGROWQUANT,
+ "FLT",
0,
- "4 or 8 connectivity in final dilation.",
+ "Minimum quant. to expand true detections.",
ARGS_GROUP_DETECTION,
- &p->dilatengb,
- 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
},
{
- "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..9bbfb21 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -39,8 +39,7 @@
dthresh 0.0
detsnminarea 10
detquant 0.95
- dilate 3
- dilatengb 8
+ detgrowquant 0.51
# Segmentation
segsnminarea 15
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index e2a5d9c..fc9fa75 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)
{
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 839e8ad..f1c3746 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);
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,7 +744,7 @@ 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)
{
@@ -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);
@@ -887,13 +888,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 +909,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 +920,73 @@ detection_remove_false_initial(struct noisechiselparams *p,
+static size_t
+detection_quantile_expand(struct noisechiselparams *p, gal_data_t **workbin_i)
+{
+ gal_data_t *workbin=*workbin_i; /* workbin comes from and is later used */
+ /* outside this function. */
+ int32_t *o, *of;
+ size_t *d, counter=0;
+ gal_data_t *exp_thresh_full, *diffuseindexs;
+ uint8_t *b=workbin->array, *bf=b+workbin->size;
+ float *e_th, *arr=p->conv ? p->conv->array : p->input->array;
+
+ /* Expand the threshold values (from one value for each tile) to the
+ whole dataset. */
+ exp_thresh_full=gal_tile_block_write_const_value(p->expand_thresh,
+ p->cp.tl.tiles, 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 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;
+ d=diffuseindexs->array;
+ e_th=exp_thresh_full->array;
+ of=(o=p->olabel->array)+p->olabel->size;
+ arr=p->conv ? p->conv->array : p->input->array;
+ 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, 0);
+ gal_binary_fill_holes(workbin, 1);
+
+ /* Clean up and return. */
+ gal_data_free(exp_thresh_full);
+ return gal_binary_connected_components(workbin, &p->olabel, 8);
+}
+
+
+
+
+
/* The initial detection has been done, now we want to remove false
detections. */
void
@@ -963,6 +1030,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)
{
@@ -973,31 +1041,33 @@ detection(struct noisechiselparams *p)
}
- /* If the user asked for dilation, then apply it. */
- 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 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->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(&t0, 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)
+ if(p->cleangrowndet)
p->numdetections = detection_final_remove_small_sn(p, num_true_initial);
else
+
{
p->numdetections = num_true_initial;
if(p->detectionname)
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index bdd7396..b60c2a0 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -72,9 +72,8 @@ 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. */
+ 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. */
@@ -105,6 +104,7 @@ struct noisechiselparams
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/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..67fd346 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -202,7 +202,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..cec4a72 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,43 @@ 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";
+ if(third) (*third)->name="THRESH3_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)
+ gal_tile_full_values_write(*third, tl, filename, NULL, PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
+ if(third) (*third)->name=NULL;
}
/* Smooth the threshold if requested. */
@@ -292,16 +311,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";
+ if(third) (*third)->name="THRESH3_SMOOTH";
gal_tile_full_values_write(*first, tl, filename, NULL,
PROGRAM_NAME);
gal_tile_full_values_write(*second, tl, filename, NULL,
PROGRAM_NAME);
+ if(third)
+ gal_tile_full_values_write(*third, tl, filename, NULL,
+ PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
+ if(third) (*third)->name=NULL;
}
}
}
@@ -332,6 +365,7 @@ struct qthreshparams
{
gal_data_t *erode_th;
gal_data_t *noerode_th;
+ gal_data_t *expand_th;
void *usage;
struct noisechiselparams *p;
};
@@ -408,11 +442,17 @@ 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. */
+ 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 +460,8 @@ qthresh_on_tile(void *in_prm)
tind, type), type);
gal_blank_write(gal_data_ptr_increment(qprm->noerode_th->array,
tind, type), type);
+ gal_blank_write(gal_data_ptr_increment(qprm->expand_th->array,
+ tind, type), type);
}
/* Clean up and fix the tile's pointers. */
@@ -467,6 +509,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,8 +531,12 @@ 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;
+ qprm.expand_th->flag |= GAL_DATA_FLAG_BLANK_CH;
if(p->qthreshname)
{
qprm.erode_th->name="QTHRESH_ERODE";
@@ -493,13 +545,21 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
PROGRAM_NAME);
gal_tile_full_values_write(qprm.noerode_th, tl, 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, 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,
- p->qthreshname);
+ &qprm.expand_th, p->qthreshname);
/* We now have a threshold for all tiles, apply it. */
@@ -512,6 +572,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..3e72093 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. */
@@ -360,7 +356,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. */
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 416cd76..10eaef4 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -48,7 +48,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',
@@ -71,8 +71,7 @@ enum option_keys_enum
UI_KEY_OPENINGNGB,
UI_KEY_CHECKDETSKY,
UI_KEY_CHECKDETSN,
- UI_KEY_DILATENGB,
- UI_KEY_CLEANDILATED,
+ UI_KEY_CLEANGROWNDET,
UI_KEY_CHECKDETECTION,
UI_KEY_CHECKSKY,
UI_KEY_CHECKCLUMPSN,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 66e703c..b213e73 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12195,29 +12195,20 @@ 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.
@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
@@ -12660,18 +12651,46 @@ 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.
+
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
@@ -21730,10 +21749,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..f793057 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,7 +641,7 @@ 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)
{
uint8_t *in;
gal_data_t *inv, *tile, *holelabs=NULL;
@@ -658,8 +657,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 */
+ gal_binary_connected_components(inv, &holelabs, connectivity);
/* Any pixel with a label larger than 1 is a hole in the input image and
diff --git a/lib/gnuastro/binary.h b/lib/gnuastro/binary.h
index cb19650..7acc7aa 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);
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