[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 9738849 2/3: Wide kernel for NoiseChisel's qth
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 9738849 2/3: Wide kernel for NoiseChisel's qthresh estimation |
|
Date: |
Thu, 14 Sep 2017 19:27:23 -0400 (EDT) |
branch: master
commit 9738849c352334296b09cdfa6c4a9ef9215b794e
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Wide kernel for NoiseChisel's qthresh estimation
NoiseChisel uses the difference between the mode and median to identify if
a tile should be used for estimating the quantile thresholds. Until now,
NoiseChisel would convolve an image once and estimate the proper tiles for
quantile estimations on the convolved image and use the same convolved
image for quantile estimation. While a larger kernel does increase the
skewness (and thus difference between the mode and median), it disfigures
the shapes of the objects.
With this commit, there is a new `--widekernel' (and a corresponding
`--wkhdu' option to specify its HDU) option. When its given, the image will
be convolved with both the sharp 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.
Library functions that are changed with this commit:
- `gal_binary_fill_holes' now accepts a `maxsize' argument to only fill
holes less than a given size.
- `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.
---
bin/noisechisel/args.h | 43 ++++++++++++++++-
bin/noisechisel/astnoisechisel.conf | 44 +++++++++---------
bin/noisechisel/clumps.c | 14 ++++--
bin/noisechisel/detection.c | 93 ++++++++++++++++++++-----------------
bin/noisechisel/main.h | 9 +++-
bin/noisechisel/noisechisel.c | 27 ++++++++---
bin/noisechisel/sky.c | 6 ++-
bin/noisechisel/threshold.c | 93 ++++++++++++++++++++++++-------------
bin/noisechisel/ui.c | 68 +++++++++++++++++++--------
bin/noisechisel/ui.h | 5 +-
bin/statistics/sky.c | 16 +++----
bin/statistics/statistics.c | 2 +-
doc/gnuastro.texi | 79 +++++++++++++++++++++++++++----
lib/binary.c | 22 ++++++++-
lib/blank.c | 8 ++--
lib/gnuastro/binary.h | 2 +-
lib/gnuastro/tile.h | 6 +--
lib/statistics.c | 2 +-
lib/tile.c | 18 +++----
19 files changed, 388 insertions(+), 169 deletions(-)
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index dddb8a6..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,
@@ -378,6 +404,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "detgrowmaxholesize",
+ UI_KEY_DETGROWMAXHOLESIZE,
+ "INT",
+ 0,
+ "Max. area of holes after growth to fill.",
+ ARGS_GROUP_DETECTION,
+ &p->detgrowmaxholesize,
+ GAL_TYPE_SIZE_T,
+ GAL_OPTIONS_RANGE_GE_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"cleangrowndet",
UI_KEY_CLEANGROWNDET,
0,
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 9bbfb21..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,28 +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
- detgrowquant 0.51
+ 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 fc9fa75..df7b383 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -941,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 "
@@ -979,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 ];} );
@@ -1045,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);
@@ -1064,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). */
@@ -1295,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/detection.c b/bin/noisechisel/detection.c
index f1c3746..98bd855 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -243,7 +243,7 @@ detection_fill_holes_open(void *in_prm)
/* Fill the holes in this tile: holes with maximal connectivity means
that they are most strongly bounded. */
- gal_binary_fill_holes(copy, copy->ndim);
+ gal_binary_fill_holes(copy, copy->ndim, -1);
if(fho_prm->step==1)
{
detection_write_in_large(tile, copy);
@@ -746,7 +746,8 @@ detection_pseudo_real(struct noisechiselparams *p)
/* 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;
@@ -757,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");
@@ -769,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 )
{
@@ -810,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);
@@ -921,36 +912,35 @@ detection_remove_false_initial(struct noisechiselparams
*p,
static size_t
-detection_quantile_expand(struct noisechiselparams *p, gal_data_t **workbin_i)
+detection_quantile_expand(struct noisechiselparams *p, gal_data_t *workbin)
{
- 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;
+ 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;
- 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. */
+ 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);
+ 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 pixels that must be
- expanded and re-initialize the necessary pointers. */
+ /* 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;
- 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
@@ -975,12 +965,34 @@ detection_quantile_expand(struct noisechiselparams *p,
gal_data_t **workbin_i)
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);
+ workbin=gal_binary_dilate(workbin, 1, 1, 1);
+ gal_binary_fill_holes(workbin, 1, p->detgrowmaxholesize);
- /* Clean up and return. */
+ /* 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 gal_binary_connected_components(workbin, &p->olabel, 8);
+ return numexpanded;
}
@@ -1045,7 +1057,7 @@ detection(struct noisechiselparams *p)
final number of detections. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
if(p->detgrowquant!=1.0f)
- num_true_initial=detection_quantile_expand(p, &workbin);
+ num_true_initial=detection_quantile_expand(p, workbin);
if(!p->cp.quiet)
{
if(p->detgrowquant==1.0f)
@@ -1064,26 +1076,23 @@ detection(struct noisechiselparams *p)
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->cleangrowndet)
- 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 b60c2a0..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. */
@@ -73,6 +75,7 @@ struct noisechiselparams
uint8_t checkdetsn; /* Save pseudo-detection S/N values. */
float detquant; /* True detection quantile. */
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. */
@@ -99,8 +102,10 @@ 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. */
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/sky.c b/bin/noisechisel/sky.c
index 67fd346..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
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index cec4a72..d7da847 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -288,10 +288,12 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
(*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);
+ 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, filename, NULL, PROGRAM_NAME);
+ gal_tile_full_values_write(*third, tl, 1, filename, NULL,
+ PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
if(third) (*third)->name=NULL;
}
@@ -326,12 +328,12 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
(*first)->name="THRESH1_SMOOTH";
(*second)->name="THRESH2_SMOOTH";
if(third) (*third)->name="THRESH3_SMOOTH";
- gal_tile_full_values_write(*first, tl, filename, NULL,
+ 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, filename, NULL,
+ gal_tile_full_values_write(*third, tl, 1, filename, NULL,
PROGRAM_NAME);
(*first)->name = (*second)->name = NULL;
if(third) (*third)->name=NULL;
@@ -384,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;
@@ -409,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);
@@ -435,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);
@@ -449,10 +465,14 @@ qthresh_on_tile(void *in_prm)
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);
+ 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
{
@@ -460,8 +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);
- gal_blank_write(gal_data_ptr_increment(qprm->expand_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. */
@@ -498,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. */
@@ -536,21 +562,21 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
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->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;
if(qprm.expand_th)
{
qprm.expand_th->name="QTHRESH_EXPAND";
- gal_tile_full_values_write(qprm.expand_th, tl, p->qthreshname,
+ gal_tile_full_values_write(qprm.expand_th, tl, 1, p->qthreshname,
NULL, PROGRAM_NAME);
qprm.expand_th->name=NULL;
}
@@ -559,7 +585,8 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
/* Interpolate and smooth the derived values. */
threshold_interp_smooth(p, &qprm.erode_th, &qprm.noerode_th,
- &qprm.expand_th, p->qthreshname);
+ qprm.expand_th ? & qprm.expand_th : NULL,
+ p->qthreshname);
/* We now have a threshold for all tiles, apply it. */
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 3e72093..09f0d9d 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -245,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");
+ }
}
@@ -270,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");
- }
}
@@ -464,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);
}
@@ -572,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);
@@ -678,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);
}
}
@@ -760,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 10eaef4..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',
@@ -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,6 +73,7 @@ enum option_keys_enum
UI_KEY_OPENINGNGB,
UI_KEY_CHECKDETSKY,
UI_KEY_CHECKDETSN,
+ UI_KEY_DETGROWMAXHOLESIZE,
UI_KEY_CLEANGROWNDET,
UI_KEY_CHECKDETECTION,
UI_KEY_CHECKSKY,
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 b213e73..1e56e62 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -12188,6 +12188,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
@@ -12203,7 +12220,10 @@ 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.
+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
@option{--cleangrowndet}: A process to further clean/remove the possibility
@@ -12348,6 +12368,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
@@ -12669,9 +12706,26 @@ 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
+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
@@ -18125,7 +18179,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
@@ -20578,7 +20632,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
@@ -20599,6 +20653,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
@@ -20614,7 +20673,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})
@@ -20985,13 +21044,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})
diff --git a/lib/binary.c b/lib/binary.c
index f793057..4e5ba3e 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -641,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, int connectivity)
+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,7 +660,7 @@ gal_binary_fill_holes(gal_data_t *input, int connectivity)
/* Label the holes */
- gal_binary_connected_components(inv, &holelabs, connectivity);
+ numholes=gal_binary_connected_components(inv, &holelabs, connectivity);
/* Any pixel with a label larger than 1 is a hole in the input image and
@@ -668,6 +670,22 @@ gal_binary_fill_holes(gal_data_t *input, int connectivity)
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 7acc7aa..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, int connectivity);
+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);