[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 15c524d 124/125: New implementation of NoiseCh
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete |
|
Date: |
Sun, 23 Apr 2017 22:36:52 -0400 (EDT) |
branch: master
commit 15c524db1964853cea6de3fd9902effa7422d557
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
New implementation of NoiseChisel complete
The new implementation of NoiseChisel is now complete (with the
documentation and also several improvements compared to the last
commit). The list of changes are:
- With the `--onlydetection' option, no segmentation will be done and
there will also be no clumps extension. This is also a new name to the
old `--onlydetect' option.
- As a result of the change above, in MakeCatalog, the default `skyhdu'
and `stdhdu' have been changed to the extension names instead of
numbers.
- The old `--minbfrac' option in NoiseChisel is now called `--minskyfrac'
to be more clear, `b' might be confused with `blank'.
- A major change compared to the previous commit is that the labeled
images (for clumps and objects) now have the type `int32_t' (before they
were `uint32_t'). This change was made because it is much more easier to
use negative values for internal values and also for the diffuse region
(for the clumps image). Also, the signed 32-bit integer is a standard
recognized CFITSIO type. All the implementation of segmentation were
also corrected.
- A new bit-wise `flag' element was added to `gal_data_t' along with the
necessary bit-wise macros for having blank values and if the dataset is
sorted. Although currently on the blank values flag are
implemented. `gal_blank_present' will now use these flags (if they are
present) and this can greatly help in not re-checking a dataset. This
benefited many library functions that NoiseChisel currently uses now and
will benefit many more in the future.
- NoiseChisel's output functions have been written.
- `gal_type_bit_string' has replaced `gal_data_bit_print_stream', since we
usually want the bit-stream when thinking in a specific type.
- The Arithmetic program tests that depended on NoiseChisel have been
corrected and now pass.
- `tests/during-dev.sh' now also accepts the number of threads to build
Gnuastro with.
- `tmpfs-config-make' now also builds with no debugging and shared
libraries by default. With the next commit (when work on MakeCatalog
starts), all the separate `--enable-XXX' for programs will also be
removed it.
---
bin/mkcatalog/astmkcatalog.conf | 4 +-
bin/noisechisel/args.h | 40 +-
bin/noisechisel/astnoisechisel.conf | 2 +-
bin/noisechisel/clumps.c | 386 ++++++++++++-------
bin/noisechisel/clumps.h | 7 +-
bin/noisechisel/detection.c | 83 +++--
bin/noisechisel/main.h | 6 +-
bin/noisechisel/noisechisel.c | 157 +++++++-
bin/noisechisel/segmentation.c | 93 ++---
bin/noisechisel/sky.c | 2 +-
bin/noisechisel/threshold.c | 115 ++++--
bin/noisechisel/ui.c | 35 +-
bin/noisechisel/ui.h | 4 +-
bin/statistics/ui.c | 2 +
doc/gnuastro.texi | 720 ++++++++++++++++++++----------------
lib/binary.c | 28 +-
lib/blank.c | 184 +++++----
lib/convolve.c | 12 +-
lib/data.c | 80 ++--
lib/gnuastro/data.h | 40 +-
lib/gnuastro/tile.h | 5 +
lib/gnuastro/type.h | 3 +-
lib/interpolate.c | 63 +++-
lib/statistics.c | 34 +-
lib/tile.c | 61 ++-
lib/type.c | 34 +-
tests/arithmetic/snimage.sh | 2 +-
tests/arithmetic/where.sh | 2 +-
tests/during-dev.sh | 10 +-
tests/noisechisel/noisechisel.sh | 2 +-
tmpfs-config-make | 2 +-
31 files changed, 1442 insertions(+), 776 deletions(-)
diff --git a/bin/mkcatalog/astmkcatalog.conf b/bin/mkcatalog/astmkcatalog.conf
index 89e5e96..cdd17c8 100644
--- a/bin/mkcatalog/astmkcatalog.conf
+++ b/bin/mkcatalog/astmkcatalog.conf
@@ -20,8 +20,8 @@
# Input:
objhdu 2
clumphdu 3
- skyhdu 4
- stdhdu 5
+ skyhdu SKY
+ stdhdu SKY_STD
zeropoint 0.0
skysubtracted 0
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index b46cdf0..b265b2c 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -73,13 +73,13 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "minbfrac",
- ARGS_OPTION_KEY_MINBFRAC,
+ "minskyfrac",
+ ARGS_OPTION_KEY_MINSKYFRAC,
"FLT",
0,
"Min. fraction of undetected area in tile.",
GAL_OPTIONS_GROUP_INPUT,
- &p->minbfrac,
+ &p->minskyfrac,
GAL_TYPE_FLOAT32,
GAL_OPTIONS_RANGE_GE_0_LE_1,
GAL_OPTIONS_MANDATORY,
@@ -120,13 +120,13 @@ struct argp_option program_options[] =
/* Output options. */
{
- "onlydetect",
- ARGS_OPTION_KEY_ONLYDETECT,
+ "onlydetection",
+ ARGS_OPTION_KEY_ONLYDETECTION,
0,
0,
- "Don't do any segmentation.",
+ "Stop at the end of detection.",
GAL_OPTIONS_GROUP_OUTPUT,
- &p->onlydetect,
+ &p->onlydetection,
GAL_OPTIONS_NO_ARG_TYPE,
GAL_OPTIONS_RANGE_0_OR_1,
GAL_OPTIONS_NOT_MANDATORY,
@@ -426,19 +426,6 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
- "segquant",
- ARGS_OPTION_KEY_SEGQUANT,
- "FLT",
- 0,
- "S/N Quantile of true sky clumps.",
- ARGS_GROUP_SEGMENTATION,
- &p->segquant,
- GAL_TYPE_FLOAT32,
- GAL_OPTIONS_RANGE_GT_0,
- GAL_OPTIONS_MANDATORY,
- GAL_OPTIONS_NOT_SET
- },
- {
"checkclumpsn",
ARGS_OPTION_KEY_CHECKCLUMPSN,
0,
@@ -452,6 +439,19 @@ struct argp_option program_options[] =
GAL_OPTIONS_NOT_SET
},
{
+ "segquant",
+ ARGS_OPTION_KEY_SEGQUANT,
+ "FLT",
+ 0,
+ "S/N Quantile of true sky clumps.",
+ ARGS_GROUP_SEGMENTATION,
+ &p->segquant,
+ GAL_TYPE_FLOAT32,
+ GAL_OPTIONS_RANGE_GT_0,
+ GAL_OPTIONS_MANDATORY,
+ GAL_OPTIONS_NOT_SET
+ },
+ {
"keepmaxnearriver",
ARGS_OPTION_KEY_KEEPMAXNEARRIVER,
0,
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 2d5cf8b..23a78e9 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -19,7 +19,7 @@
# Input:
khdu 1
- minbfrac 0.7
+ minskyfrac 0.7
minnumfalse 100
# Tessellation
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
index ad20722..2a5541b 100644
--- a/bin/noisechisel/clumps.c
+++ b/bin/noisechisel/clumps.c
@@ -76,13 +76,13 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
size_t *a, *af, ind, *dsize=p->input->dsize;
struct gal_linkedlist_sll *Q=NULL, *cleanup=NULL;
size_t *dinc=gal_dimension_increment(ndim, dsize);
- uint32_t n1, nlab, rlab, curlab=1, *clabel=p->clabel->array;
+ int32_t n1, nlab, rlab, curlab=1, *clabel=p->clabel->array;
/*********************************************
For checks and debugging:*
gal_data_t *crop;
size_t extcount=1;
- uint32_t *cr, *crf;
+ int32_t *cr, *crf;
size_t checkdsize[2]={10,10};
size_t checkstart[2]={50,145};
char *filename="clumpbuild.fits";
@@ -100,9 +100,9 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
gal_fits_img_write(crop, filename, NULL, PROGRAM_STRING);
gal_data_free(crop);
printf("blank: %u\nriver: %u\ntmpcheck: %u\ninit: %u\nmaxlab: %u\n",
- (uint32_t)GAL_BLANK_UINT32, (uint32_t)CLUMPS_RIVER,
- (uint32_t)CLUMPS_TMPCHECK, (uint32_t)CLUMPS_INIT,
- (uint32_t)CLUMPS_MAXLAB);
+ (int32_t)GAL_BLANK_INT32, (int32_t)CLUMPS_RIVER,
+ (int32_t)CLUMPS_TMPCHECK, (int32_t)CLUMPS_INIT,
+ (int32_t)CLUMPS_MAXLAB);
tile->array=gal_tile_block_relative_to_other(tile, p->clabel);
tile->block=p->clabel;
**********************************************/
@@ -168,44 +168,64 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
label. */
GAL_DIMENSION_NEIGHBOR_OP(ind, ndim, dsize, ndim, dinc,
{
- /* For easy reading. */
- nlab=clabel[ nind ];
-
- /* This neighbor is within the search region. */
- if(nlab)
+ /* If it is already decided to be a river, then stop
+ looking at the neighbors. */
+ if(n1!=CLUMPS_RIVER)
{
- /* If this neighbour has not been labeled yet
- and has an equal flux, add it to the queue
- to expand the studied region.*/
- if( nlab==CLUMPS_INIT && arr[nind]==arr[*a] )
- {
- clabel[nind]=CLUMPS_TMPCHECK;
- gal_linkedlist_add_to_sll(&Q, nind);
- gal_linkedlist_add_to_sll(&cleanup, nind);
- }
+ /* For easy reading. */
+ nlab=clabel[ nind ];
- /* If this neighbour has a positive nlab, it
- belongs to another object, so if `n1' has not
- been set for the whole region, put `nlab' into
- `n1'. If `n1' has been set and is different
- from `nlab' then this whole equal flux region
- should be a wide river because it is connecting
- two connected regions.*/
- else if(nlab<CLUMPS_MAXLAB)
+ /* This neighbor's label isn't zero. */
+ if(nlab)
{
- if(n1==0) n1 = nlab;
- else if(nlab!=n1) n1 = CLUMPS_RIVER;
+ /* If this neighbor has not been labeled yet
+ and has an equal flux, add it to the queue
+ to expand the studied region.*/
+ if( nlab==CLUMPS_INIT && arr[nind]==arr[*a] )
+ {
+ clabel[nind]=CLUMPS_TMPCHECK;
+ gal_linkedlist_add_to_sll(&Q, nind);
+ gal_linkedlist_add_to_sll(&cleanup, nind);
+ }
+ else
+ n1=( nlab>0
+
+ /* If this neighbor has a positive
+ nlab, it belongs to another object,
+ so if `n1' has not been set for the
+ whole region (n1==0), put `nlab'
+ into `n1'. If `n1' has been set and
+ is different from `nlab' then this
+ whole equal flux region should be a
+ wide river because it is connecting
+ two connected regions.*/
+ ? ( n1
+ ? (n1==nlab ? n1 : CLUMPS_RIVER)
+ : nlab )
+
+ /* If the data has blank pixels (recall
+ that blank in int32 is negative),
+ see if the neighbor is blank and if
+ so, set the label to a river. Since
+ the flag checking can be done
+ outside this loop, for datasets with
+ no blank element this last step will
+ be completley ignored. */
+ : ( ( (p->input->flag
+ & GAL_DATA_FLAG_HASBLANK)
+ && nlab==GAL_BLANK_INT32 )
+ ? CLUMPS_RIVER : n1 ) );
}
- }
- /* If this neigbour has a label of zero, then we are
- on the edge of the indexed region (the neighbor is
- not in the initial list of pixels to segment). When
- over-segmenting the noise and the detections,
- `clabel' is zero for the parts of the image that we
- are not interested in here. */
- else
- clabel[*a]=CLUMPS_RIVER;
+ /* If this neigbour has a label of zero, then we
+ are on the edge of the indexed region (the
+ neighbor is not in the initial list of pixels
+ to segment). When over-segmenting the noise and
+ the detections, `clabel' is zero for the parts
+ of the image that we are not interested in
+ here. */
+ else clabel[*a]=CLUMPS_RIVER;
+ }
} );
}
@@ -242,30 +262,57 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
n1=0;
/* Go over all the fully connected neighbors of this pixel and
- see if all the neighbors (connectivity equals the number of
- dimensions) that have a non-macro value (less than
+ see if all the neighbors (with maximum connectivity: the
+ number of dimensions) that have a non-macro value (less than
CLUMPS_MAXLAB) belong to one label or not. If the pixel is
neighboured by more than one label, set it as a river
pixel. Also if it is touching a zero valued pixel (which
does not belong to this object), set it as a river pixel.*/
GAL_DIMENSION_NEIGHBOR_OP(*a, ndim, dsize, ndim, dinc,
{
- nlab=clabel[ nind ];
- if(nlab<CLUMPS_MAXLAB)
+ /* When `n1' has already been set as a river, there is no
+ point in looking at the other neighbors. */
+ if(n1!=CLUMPS_RIVER)
{
- if(nlab==0) n1=CLUMPS_RIVER;
- else
- {
- if(n1)
- {
- if(nlab!=n1) n1=CLUMPS_RIVER;
- }
- else n1=nlab;
- }
+ /* For easy reading. */
+ nlab=clabel[ nind ];
+
+ /* If this neighbor is on a non-processing label, then
+ set the first neighbor accordingly. Note that we
+ also want the zero valued neighbors (detections if
+ working on sky, and sky if working on detection):
+ we want rivers between the two domains. */
+ n1 = ( nlab
+
+ /* nlab is non-zero. */
+ ? ( nlab>0
+
+ /* Neighbor has a meaningful label, so
+ check with any previously found labeled
+ neighbors. */
+ ? ( n1
+ ? ( nlab==n1 ? n1 : CLUMPS_RIVER )
+ : nlab )
+
+ /* If the dataset has blank values and this
+ neighbor is blank, then the pixel should
+ be a river. Note that the blank checking
+ can be optimized out, so if the input
+ doesn't have blank values,
+ `nlab==GAL_BLANK_INT32' will never be
+ checked. */
+ : ( (p->input->flag & GAL_DATA_FLAG_HASBLANK)
+ && nlab==GAL_BLANK_INT32
+ ? CLUMPS_RIVER : n1 ) )
+
+ /* `nlab==0' (the neighbor lies in the other
+ domain (sky or detections). To avoid the
+ different domains touching, this pixel
+ should be a river. */
+ : CLUMPS_RIVER );
}
});
-
/* Either assign a new label to this pixel, or give it the one
of its neighbors. If n1 equals zero, then this is a new
peak, and a new label should be created. But if n1!=0, it
@@ -285,7 +332,6 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
clabel[ *a ] = rlab;
}
-
/*********************************************
For checks and debugging:
if( *a / dsize[1] >= checkstart[0]
@@ -311,13 +357,13 @@ clumps_oversegment(struct clumps_thread_params *cltprm)
/* Set all the river pixels to zero. Note that this is only necessary for
the detected clumps. When finding clumps over the Sky, we will be
- going over the full tile and removing rivers. This is because, we set
- the borders of the tile to a river value and didn't include them in
- the list of indexs. */
+ going over the full tile and removing rivers after this function. This
+ is because, we set the borders of the tile to a river value and didn't
+ include them in the list of indexs. */
if(cltprm->clprm->sky0_det1)
{
af=(a=indexs->array)+indexs->size;
- do if( clabel[*a]==CLUMPS_RIVER ) clabel[*a]=0; while(++a<af);
+ do if( clabel[*a]==CLUMPS_RIVER ) clabel[*a]=CLUMPS_INIT; while(++a<af);
}
/*********************************************
@@ -368,7 +414,7 @@ clumps_grow_prepare_initial(struct clumps_thread_params
*cltprm)
size_t ndiffuse=0, coord[2], *dindexs;
double wcoord[2]={0.0f,0.0f}, brightness=0.0f;
float glimit, *imgss=input->array, *std=p->std->array;
- uint32_t *olabel=p->olabel->array, *clabel=p->clabel->array;
+ int32_t *olabel=p->olabel->array, *clabel=p->clabel->array;
/* Find the flux weighted center (meaningful only for positive valued
@@ -423,20 +469,17 @@ clumps_grow_prepare_initial(struct clumps_thread_params
*cltprm)
sf=(s=indexs->array)+indexs->size;
do
{
- if( clabel[*s] ) olabel[*s] = clabel[*s];
- else
- {
- olabel[*s] = CLUMPS_INIT;
- if( imgss[*s]>glimit ) dindexs[ ndiffuse++ ] = *s;
- }
+ olabel[*s] = clabel[*s];
+ if( clabel[*s]==CLUMPS_INIT )
+ if( imgss[*s]>glimit ) dindexs[ ndiffuse++ ] = *s;
}
while(++s<sf);
/* When a user just wants clumps and doesn't want multiple clumps to be
considered part of one object, they are going to set an impossibly
- high `--gthresh' so no growth actually happens. In that case, we don't
- need the allocated array of indexs. So just free it. */
+ high `--gthresh', so no growth actually happens. In that case, we
+ don't need the allocated array of indexs. So just free it. */
if(ndiffuse==0)
{
free(cltprm->diffuseindexs->array);
@@ -467,14 +510,14 @@ clumps_grow_prepare_final(struct clumps_thread_params
*cltprm)
{
size_t ndiffuse=0;
size_t *dindexs=cltprm->diffuseindexs->array;
- uint32_t *olabel=cltprm->clprm->p->olabel->array;
+ int32_t *olabel=cltprm->clprm->p->olabel->array;
size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
/* Recall that we initially allocated `diffuseindexs' to have the same
size as the indexs. So there is no problem if there are more pixels in
this final round compared to the initial round. */
do
- if( olabel[*s] > CLUMPS_MAXLAB )
+ if( olabel[*s] < 0 )
dindexs[ ndiffuse++ ] = *s;
while(++s<sf);
@@ -505,9 +548,9 @@ clumps_grow(struct clumps_thread_params *cltprm, int
withrivers)
gal_data_t *diffuseindexs=cltprm->diffuseindexs;
size_t ndim=p->input->ndim, *dsize=p->input->dsize;
- int stopngbsearch;
+ int searchngb;
size_t *diarray=cltprm->diffuseindexs->array;
- uint32_t n1, nlab, *olabel=p->olabel->array;
+ int32_t n1, nlab, *olabel=p->olabel->array;
size_t *dinc=gal_dimension_increment(ndim, dsize);
size_t *s, *sf, thisround, ndiffuse=diffuseindexs->size;
@@ -544,20 +587,23 @@ clumps_grow(struct clumps_thread_params *cltprm, int
withrivers)
/* Check the very closest neighbors of each pixel (4-connectivity
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. */
- stopngbsearch=0;
+ searchngb=1;
GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, 1, dinc,
{
- if(!stopngbsearch)
+ if(searchngb)
{
+ /* For easy reading. */
nlab = olabel[nind];
- if(nlab && nlab<CLUMPS_MAXLAB) /* This is a real label. */
+
+ /* This neighbor's label is meaningful. */
+ if(nlab>0) /* This is a real label. */
{
if(n1) /* A prev. neighboring label has been found. */
{
if( n1 != nlab ) /* Different label from */
{ /* prevously found neighbor for this pixel. */
n1=CLUMPS_RIVER;
- stopngbsearch=1;
+ searchngb=0;
}
}
else
@@ -568,7 +614,7 @@ clumps_grow(struct clumps_thread_params *cltprm, int
withrivers)
(`withrivers==0'), then there is no point in
looking in other neighbors, the first neighbor
we find is the one we'll use. */
- if(!withrivers) stopngbsearch=1;
+ if(!withrivers) searchngb=0;
}
}
}
@@ -579,20 +625,17 @@ clumps_grow(struct clumps_thread_params *cltprm, int
withrivers)
n1==0 --> No labeled neighbor was found.
n1==CLUMPS_RIVER --> Connecting two labeled regions.
- n1>0 && n1<CLUMPS_MAXLAB --> Only has one neighbouring label.
+ n1>0 --> Only has one neighbouring label.
- The first one means that no neighbors was found and this pixel
- should be kept for the next loop (we'll be growing the objects
- pixel-layer by pixel-layer). In the other two cases, we just
- need to write in the value of `n1'. */
+ The first one means that no neighbors were found and this
+ pixel should be kept for the next loop (we'll be growing the
+ objects pixel-layer by pixel-layer). In the other two cases,
+ we just need to write in the value of `n1'. */
if(n1)
{
olabel[*s]=n1;
- if(n1>CLUMPS_MAXLAB) /* We want to keep river */
- {
- if(withrivers==0) printf("%zu\n", *s);
- diarray[ ndiffuse++ ] = *s; /* pixels the diffuse list.*/
- }
+ if(withrivers && n1==CLUMPS_RIVER) /* To keep rivers in */
+ diarray[ ndiffuse++ ] = *s; /* the diffuse list. */
}
else
diarray[ ndiffuse++ ] = *s;
@@ -659,10 +702,10 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
size_t nngb=gal_dimension_num_neighbors(ndim);
float *arr=p->input->array, *std=p->std->array;
size_t *dinc=gal_dimension_increment(ndim, dsize);
- uint32_t lab, nlab, *ngblabs, *clabel=p->clabel->array;
+ int32_t lab, nlab, *ngblabs, *clabel=p->clabel->array;
/* Allocate the array to keep the neighbor labels of river pixels. */
- ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb);
+ ngblabs=gal_data_malloc_array(GAL_TYPE_INT32, nngb);
/* Go over all the pixels in this region. */
af=(a=cltprm->indexs->array)+cltprm->indexs->size;
@@ -670,7 +713,7 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
if( !isnan(arr[ *a ]) )
{
/* This pixel belongs to a clump. */
- if( clabel[ *a ] )
+ if( clabel[ *a ]>0 )
{
lab=clabel[*a];
++info[ lab * INFO_NCOLS + INFO_INAREA ];
@@ -706,7 +749,7 @@ clumps_get_raw_info(struct clumps_thread_params *cltprm)
/* We only want those neighbors that are not rivers (>0) or
any of the flag values. */
- if(nlab && nlab<CLUMPS_MAXLAB)
+ if(nlab>0)
{
/* Go over all already checked labels and make sure
this clump hasn't already been considered. */
@@ -771,7 +814,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
size_t tablen=cltprm->numinitclumps+1;
float *snarr;
- uint32_t *indarr=NULL;
+ int32_t *indarr=NULL;
double I, O, Ni, var, *row;
int sky0_det1=cltprm->clprm->sky0_det1;
size_t i, ind, counter=0, infodsize[2]={tablen, INFO_NCOLS};
@@ -789,7 +832,7 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
{
cltprm->snind = &cltprm->clprm->snind [ cltprm->id ];
cltprm->snind->ndim = 1; /* Depends on `cltprm->snind' */
- cltprm->snind->type = GAL_TYPE_UINT32;
+ cltprm->snind->type = GAL_TYPE_INT32;
cltprm->snind->dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, 1);
cltprm->snind->size = cltprm->snind->dsize[0]=tablen;/* After dsize */
cltprm->snind->array = gal_data_malloc_array(cltprm->snind->type,
@@ -804,9 +847,11 @@ clumps_make_sn_table(struct clumps_thread_params *cltprm)
cltprm->info=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 2, infodsize,
NULL, 1, p->cp.minmapsize, NULL, NULL, NULL);
+
/* First get the raw information necessary for making the S/N table. */
clumps_get_raw_info(cltprm);
+
/* Calculate the signal to noise for successful clumps */
snarr=cltprm->sn->array;
if(cltprm->snind) indarr=cltprm->snind->array;
@@ -875,9 +920,9 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
gal_data_t *tile)
{
gal_data_t *newinds;
+ int32_t *ninds, curlab, *l, *lf;
size_t len=cltprm->numinitclumps+1;
struct noisechiselparams *p=cltprm->clprm->p;
- uint32_t *ninds, curlab, *l=cltprm->snind->array, *lf=l+cltprm->snind->size;
/* A small sanity check. */
if(gal_tile_block(tile)!=p->clabel)
@@ -889,7 +934,7 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
/* Allocate a dataset with the new indexs, note that it will need to have
one element for each initial label (the excluded clumps need to be set
to zero). So we also need to clear the allocated space. */
- newinds=gal_data_alloc(NULL, p->clabel->type, 1, &len, NULL, 1,
+ newinds=gal_data_alloc(NULL, p->clabel->type, 1, &len, NULL, 0,
p->cp.minmapsize, NULL, NULL, NULL);
@@ -904,16 +949,23 @@ clumps_correct_sky_labels_for_check(struct
clumps_thread_params *cltprm,
if(p->cp.numthreads>1) pthread_mutex_unlock(&cltprm->clprm->labmutex);
+ /* Initialize the newinds array to CLUMPS_INIT (which be used as a new
+ label for all the clumps that must be removed. */
+ lf = (l=newinds->array) + newinds->size; do *l++=CLUMPS_INIT; while(l<lf);
+
+
/* The new indexs array has been initialized to zero. So we just need to
go over the labels in `cltprm->sninds' and give them a value of
`curlab++'. */
ninds=newinds->array;
+ lf = (l=cltprm->snind->array) + cltprm->snind->size;
do { ninds[*l]=curlab++; *l=ninds[*l]; } while(++l<lf);
- /* Go over this tile and correct the values. */
- GAL_TILE_PARSE_OPERATE({*i = ninds[*(uint32_t *)i];}, tile, NULL, 0, 1);
+ /* Go over this tile and correct the values. */
+ GAL_TILE_PARSE_OPERATE({if(*i>0) *i=ninds[ *(int32_t *)i ];},
+ tile, NULL, 0, 1);
/* Clean up. */
gal_data_free(newinds);
@@ -936,7 +988,7 @@ clumps_find_make_sn_table(void *in_prm)
gal_data_t *tile, *tblock, *tmp;
uint8_t *binary=p->binary->array;
struct clumps_thread_params cltprm;
- size_t i, c, ind, tind, numsky, *indarr;
+ size_t i, c, ind, tind, num, numsky, *indarr;
size_t *scoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
size_t *icoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
@@ -946,7 +998,7 @@ clumps_find_make_sn_table(void *in_prm)
cltprm.topinds = NULL;
- /* Go over all the tiles given to this thread. */
+ /* Go over all the tiles/detections given to this thread. */
for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
{
/* IDs. */
@@ -960,6 +1012,15 @@ clumps_find_make_sn_table(void *in_prm)
tile->array = gal_tile_block_relative_to_other(tile, p->binary);
tile->block = p->binary;
+ /* Get the number of usable elements in this tile (note that tiles
+ can have blank pixels), so we can't simply use `tile->size'. */
+ if(tile->flag & GAL_DATA_FLAG_HASBLANK)
+ {
+ tmp=gal_statistics_number(tile);
+ num=*((size_t *)(tmp->array));
+ gal_data_free(tmp);
+ }
+ else num=tile->size;
/* Find the number of detected pixels over this tile. Since this is
the binary image, this is just the sum of all the pixels. */
@@ -967,11 +1028,12 @@ clumps_find_make_sn_table(void *in_prm)
numdet=*((double *)(tmp->array));
gal_data_free(tmp);
-
/* See if this tile should be used or not (has enough undetected
- pixels). */
- numsky=tile->size-numdet;
- if( (float)numsky/(float)tile->size > p->minbfrac )
+ pixels). Note that it might happen that some tiles are fully
+ blank. In such cases, it is important to first check the number of
+ detected pixels. */
+ numsky=num-numdet;
+ if( num && (float)numsky/(float)num > p->minskyfrac )
{
/* Add the indexs of all undetected pixels in this tile into an
array. */
@@ -999,23 +1061,44 @@ clumps_find_make_sn_table(void *in_prm)
indarr=cltprm.indexs->array;
GAL_TILE_PARSE_OPERATE({
/* This pixel's index over all the image. */
- ind = (uint32_t *)i - (uint32_t *)(p->clabel->array);
+ ind = (int32_t *)i - (int32_t *)(p->clabel->array);
gal_dimension_index_to_coord(ind, ndim, dsize, icoord);
- /* If the pixel is on the edge, set it as river and
+ /* If the pixel is on the tile edge, set it as river and
don't include it in the indexs. */
if( icoord[0]==scoord[0]
|| icoord[0]==scoord[0]+tile->dsize[0]-1
|| icoord[1]==scoord[1]
|| icoord[1]==scoord[1]+tile->dsize[1]-1 )
- *(uint32_t *)i=CLUMPS_RIVER;
-
- /* This pixel is not on the edge, check if it had a value
- of `0' in the binary image (is not detected), and if
- so, then add it to the list of indexs. */
- else if(binary[ind]==0)
- indarr[c++]=gal_data_ptr_dist(tile->block->array, i,
- p->clabel->type);
+ *(int32_t *)i=CLUMPS_RIVER;
+
+ /* This pixel is not on the edge, check if it had a value of
+ `0' in the binary image (is not detected) then add it to
+ the list of indexs (note that the binary image also
+ contains the blank pixels, so only sky regions have a
+ value of 0 in the binary image). */
+ else if( binary[ind]==0 )
+ {
+ /*
+ if(c!=cltprm.indexs->size)
+ {
+ if(cltprm.id==282) *i+=2;
+ */
+ indarr[c++]=gal_data_ptr_dist(p->clabel->array, i,
+ p->clabel->type);
+ /*
+ }
+ else
+ if(cltprm.id==282)
+ {
+ int32_t *clabel=p->clabel->array;
+ size_t kjd=gal_data_ptr_dist(p->clabel->array, i,
+ p->clabel->type);
+ printf("%zu, %zu: %u\n", kjd%dsize[1]+1,
+ kjd/dsize[1]+1, clabel[kjd]);
+ }
+ */
+ }
}, tile, NULL, 0, 1);
/* Correct the number of indexs. */
@@ -1024,9 +1107,10 @@ clumps_find_make_sn_table(void *in_prm)
/* Generate the clumps over this region. */
clumps_oversegment(&cltprm);
- /* Set all river pixels to zero. */
- GAL_TILE_PARSE_OPERATE({if(*i==CLUMPS_RIVER) *i=0;}, tile,
- NULL, 0, 1);
+ /* Set all river pixels to CLUMPS_INIT (to be distinguishable
+ from the detected regions). */
+ GAL_TILE_PARSE_OPERATE({if(*i==CLUMPS_RIVER) *i=CLUMPS_INIT;},
+ tile, NULL, 0, 1);
/* For a check, the step variable will be set. */
if(clprm->step==1)
@@ -1063,8 +1147,8 @@ clumps_find_make_sn_table(void *in_prm)
-/* The job of this function is to find the best signal to noise value to
- use as a threshold to detect real clumps.
+/* Find the true clump signal to noise value from the clumps in the sky
+ region.
Each thread will find the useful signal to noise values for the tiles
that have been assigned to it. It will then store the pointer to the S/N
@@ -1082,13 +1166,11 @@ void
clumps_true_find_sn_thresh(struct noisechiselparams *p)
{
char *msg;
- uint8_t *b;
- uint32_t *l, *lf;
struct timeval t1;
size_t i, j, c, numsn=0;
struct clumps_params clprm;
struct gal_linkedlist_stll *comments=NULL;
- gal_data_t *sn, *tmp, *snind, *quant, *claborig;
+ gal_data_t *sn, *snind, *quant, *claborig;
/* Get starting time for later reporting if necessary. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
@@ -1158,13 +1240,8 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
default values are hard to view, so we'll make a copy of the
demo, set all Sky regions to blank and all clump macro values
to zero. */
- tmp=gal_data_copy(p->clabel);
- b=p->binary->array; lf=(l=tmp->array)+tmp->size;
- do *l = ( *b++
- ? GAL_BLANK_UINT32
- : ( *l > CLUMPS_MAXLAB ? 0 : *l ) ); while(++l<lf);
- gal_fits_img_write(tmp, p->segmentationname, NULL, PROGRAM_STRING);
- gal_data_free(tmp);
+ gal_fits_img_write(p->clabel, p->segmentationname, NULL,
+ PROGRAM_STRING);
/* Increment the step counter. */
++clprm.step;
@@ -1181,7 +1258,6 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
p->cp.numthreads);
}
-
/* Destroy the mutex if it was initialized. */
if( p->cp.numthreads>1 && (p->checksegmentation || p->checkclumpsn) )
pthread_mutex_destroy(&clprm.labmutex);
@@ -1194,7 +1270,7 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
/* Allocate the space to keep all the S/N values. */
- sn=gal_data_alloc(NULL, clprm.sn->type, 1, &numsn, NULL, 0,
+ sn=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &numsn, NULL, 0,
p->cp.minmapsize, "CLUMP_S/N", "ratio",
"Signal-to-noise ratio");
snind = ( p->checkclumpsn
@@ -1212,12 +1288,17 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
{
((float *)(sn->array))[c] = ((float *)(clprm.sn[i].array))[j];
if(snind)
- ((uint32_t *)(snind->array))[c] =
- ((uint32_t *)(clprm.snind[i].array))[j];
+ ((int32_t *)(snind->array))[c] =
+ ((int32_t *)(clprm.snind[i].array))[j];
++c;
}
+ /* The S/N array of sky clumps is desiged to have no blank values, so set
+ the flags accordingly to avoid a redundant blank search. */
+ sn->flag = GAL_DATA_FLAG_BLANK_CH | GAL_DATA_FLAG_HASBLANK;
+
+
/* If the user wanted to see the S/N table, then save it. */
if(p->checkclumpsn)
{
@@ -1237,7 +1318,7 @@ clumps_true_find_sn_thresh(struct noisechiselparams *p)
/* Find the desired quantile from the full S/N distribution. */
- quant=gal_statistics_quantile(sn, p->segquant, 1);
+ quant = gal_statistics_quantile(sn, p->segquant, 1);
p->clumpsnthresh = *((float *)(quant->array));
if(!p->cp.quiet)
{
@@ -1283,18 +1364,32 @@ gal_data_t *
clumps_det_label_indexs(struct noisechiselparams *p)
{
size_t i, *areas;
- uint32_t *a, *l, *lf;
+ int32_t *a, *l, *lf;
gal_data_t *labindexs=gal_data_array_calloc(p->numdetections+1);
/* Find the area in each detected objects (to see how much space we need
to allocate). */
areas=gal_data_calloc_array(GAL_TYPE_SIZE_T, p->numdetections+1);
- lf=(l=p->olabel->array)+p->olabel->size; do ++areas[*l++]; while(l<lf);
+ if(p->input->flag & GAL_DATA_FLAG_HASBLANK)
+ {
+ lf=(l=p->olabel->array)+p->olabel->size; /* Blank pixels have a */
+ do if(*l>0) ++areas[*l]; while(++l<lf); /* negative value in int32. */
+ }
+ else
+ {
+ lf=(l=p->olabel->array)+p->olabel->size; do ++areas[*l]; while(++l<lf);
+ areas[0]=0;
+ }
+
+ /* For a check.
+ for(i=0;i<p->numdetections+1;++i)
+ printf("detection %zu: %zu\n", i, areas[i]);
+ exit(0);
+ */
/* Allocate/Initialize the dataset containing the indexs of each
object. We don't want the labels of the non-detected regions
(areas[0]). So we'll set that to zero.*/
- areas[0]=0;
for(i=1;i<p->numdetections+1;++i)
gal_data_initialize(&labindexs[i], NULL, GAL_TYPE_SIZE_T, 1,
&areas[i], NULL, 0, p->cp.minmapsize, NULL, NULL,
@@ -1305,7 +1400,7 @@ clumps_det_label_indexs(struct noisechiselparams *p)
memset(areas, 0, (p->numdetections+1)*sizeof *areas);
lf=(a=l=p->olabel->array)+p->olabel->size;
do
- if(*l) /* We don't want the undetected regions (*l==0) */
+ if(*l>0) /* No undetected regions (*l==0), or blank (<0) */
((size_t *)(labindexs[*l].array))[ areas[*l]++ ] = l-a;
while(++l<lf);
@@ -1328,10 +1423,15 @@ clumps_det_keep_true_relabel(struct
clumps_thread_params *cltprm)
int istouching;
size_t *s, *sf;
float *sn=cltprm->sn->array;
- uint32_t *clabel=p->clabel->array;
+ int32_t *l, *lf, curlab=1, *clabel=p->clabel->array;
size_t i, *dinc=gal_dimension_increment(ndim, dsize);
- uint32_t curlab=1, *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32,
- cltprm->numinitclumps+1);
+ int32_t *newlabs=gal_data_malloc_array(GAL_TYPE_INT32,
+ cltprm->numinitclumps+1);
+
+ /* Initialize the new labels with CLUMPS_INIT (so the diffuse area can be
+ distinguished from the clumps). */
+ lf=(l=newlabs)+cltprm->numinitclumps+1; do *l++=CLUMPS_INIT; while(l<lf);
+
/* Set the new labels. Here we will also be removing clumps with a peak
that touches a river pixel. */
@@ -1357,12 +1457,18 @@ clumps_det_keep_true_relabel(struct
clumps_thread_params *cltprm)
}
}
+
/* Save the total number of true clumps in this detection. */
cltprm->numtrueclumps=curlab-1;
- /* Correct the clump labels. */
+
+ /* Correct the clump labels. Note that the non-clumpy regions over the
+ detections (rivers) have already been initialized to CLUMPS_INIT
+ (which is negative). So we'll just need to correct the ones with a
+ value larger than 0. */
sf=(s=cltprm->indexs->array)+cltprm->indexs->size;
- do if(clabel[*s]) clabel[*s]=newlabs[ clabel[*s] ]; while(++s<sf);
+ do if(clabel[*s]>0) clabel[*s] = newlabs[ clabel[*s] ]; while(++s<sf);
+
/* Clean up. */
free(dinc);
diff --git a/bin/noisechisel/clumps.h b/bin/noisechisel/clumps.h
index 7e17c90..85e5d57 100644
--- a/bin/noisechisel/clumps.h
+++ b/bin/noisechisel/clumps.h
@@ -25,10 +25,9 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/* Constants for the clump over-segmentation. */
-#define CLUMPS_RIVER UINT32_MAX-2
-#define CLUMPS_TMPCHECK UINT32_MAX-3
-#define CLUMPS_INIT UINT32_MAX-4
-#define CLUMPS_MAXLAB UINT32_MAX-5 /* Largest clump label (unsigned). */
+#define CLUMPS_INIT -1
+#define CLUMPS_RIVER -2
+#define CLUMPS_TMPCHECK -3
/* Parameters for all threads. */
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index 1aebc48..a61646a 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -92,9 +92,9 @@ detection_initial(struct noisechiselparams *p)
p->binary->name=NULL;
}
-
/* Correct the no-erode values. */
- bf=(b=p->binary->array)+p->binary->size; do *b = *b>0; while(++b<bf);
+ bf=(b=p->binary->array)+p->binary->size;
+ do *b = *b==THRESHOLD_NO_ERODE_VALUE ? 1 : *b; while(++b<bf);
/* Do the opening. */
@@ -108,12 +108,11 @@ detection_initial(struct noisechiselparams *p)
free(msg);
}
-
/* Label the connected components. */
p->numinitialdets=gal_binary_connected_components(p->binary, &p->olabel, 1);
if(p->detectionname)
{
- p->olabel->name="OPENED-LABELED";
+ p->olabel->name="OPENED_AND_LABELED";
gal_fits_img_write(p->olabel, p->detectionname, NULL, PROGRAM_STRING);
p->olabel->name=NULL;
}
@@ -154,15 +153,22 @@ detection_initial(struct noisechiselparams *p)
static void
detection_pseudo_sky_or_det(struct noisechiselparams *p, uint8_t *w, int s0d1)
{
- uint32_t *l=p->olabel->array;
+ int32_t *l=p->olabel->array;
uint8_t *b=p->binary->array, *bf=b+p->binary->size;
if(s0d1)
- /* Set all sky regions (label equal to zero) to zero. */
+ /* Set all sky regions (label equal to zero) to zero, since a blank
+ pixel is also non-zero, we don't need to check for blanks in this
+ case. */
do *w++ = *l++ ? *b : 0; while(++b<bf);
else
- /* Set all detected pixels (label larger than zero) to blank. */
- do *w++ = *l++ ? GAL_BLANK_UINT8 : *b; while(++b<bf);
+ /* Set all detected pixels to 1. */
+ do
+ {
+ *w++ = *l ? *l==GAL_BLANK_INT32 ? GAL_BLANK_UINT8 : 1 : *b;
+ ++l;
+ }
+ while(++b<bf);
}
@@ -274,7 +280,6 @@ static size_t
detection_pseudo_find(struct noisechiselparams *p, gal_data_t *workbin,
gal_data_t *worklab, int s0d1)
{
- uint8_t *b, *bf;
gal_data_t *bin;
struct fho_params fho_prm={0, NULL, workbin, worklab, p};
@@ -356,7 +361,6 @@ detection_pseudo_find(struct noisechiselparams *p,
gal_data_t *workbin,
gal_threads_spin_off(detection_fill_holes_open, &fho_prm,
p->ltl.tottiles, p->cp.numthreads);
-
/* Clean up. */
free(fho_prm.copyspace);
@@ -366,12 +370,13 @@ detection_pseudo_find(struct noisechiselparams *p,
gal_data_t *workbin,
the blank pixels are the detections. On the Sky image, blank should be
set to 1 (because we want the detected objects to have the same labels
as the pseudo-detections that cover them). This will allow us to later
- remove these pseudo-detections. */
+ remove these pseudo-detections.
if(s0d1==0)
{
bf=(b=workbin->array)+workbin->size;
do if(*b==GAL_BLANK_UINT8) *b = !s0d1; while(++b<bf);
}
+ */
return gal_binary_connected_components(workbin, &worklab, 1);
}
@@ -379,23 +384,23 @@ detection_pseudo_find(struct noisechiselparams *p,
gal_data_t *workbin,
-#define PSN_EXTNAME "PSEUDOS-FOR-SN"
+
static gal_data_t *
detection_pseudo_sn(struct noisechiselparams *p, gal_data_t *worklab,
- size_t num, int s0d1)
+ size_t num, int s0d1)
{
float *snarr;
uint8_t *flag;
size_t tablen=num+1;
gal_data_t *sn, *snind;
- uint32_t *plabend, *indarr=NULL;
+ int32_t *plabend, *indarr=NULL;
double ave, err, *xy, *brightness;
struct gal_linkedlist_stll *comments=NULL;
size_t ind, ndim=p->input->ndim, xyncols=1+ndim;
size_t i, *area, counter=0, *dsize=p->input->dsize;
size_t *coord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
float *img=p->input->array, *f=p->input->array, *ff=f+p->input->size;
- uint32_t *plab = worklab->array, *dlab = s0d1 ? NULL : p->olabel->array;
+ int32_t *plab = worklab->array, *dlab = s0d1 ? NULL : p->olabel->array;
/* Sanity check. */
@@ -423,17 +428,17 @@ detection_pseudo_sn(struct noisechiselparams *p,
gal_data_t *worklab,
p->cp.minmapsize, "SIGNAL-TO-NOISE", "ratio",
NULL);
snind = ( p->checkdetsn==0 ? NULL
- : gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &tablen, NULL, 1,
+ : gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &tablen, NULL, 1,
p->cp.minmapsize, "LABEL", "counter",
NULL) );
-
/* Go over all the pixels and get the necessary information. */
do
{
- /* All this work os only necessary when we are actually on a
- pseudo-detection label? */
- if(*plab)
+ /* All this work is only necessary when we are actually on a
+ pseudo-detection label: it is non-zero and not blank. */
+ if(*plab && ( (p->input->flag | GAL_DATA_FLAG_HASBLANK)
+ && *plab!=GAL_BLANK_INT32 ) )
{
/* For Sky pseudo-detections we'll start to see if it has already
been determined that the object lies over a detected object or
@@ -465,6 +470,7 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t
*worklab,
}
while(++f<ff);
+
/* A small sanity check.
{
size_t i;
@@ -483,10 +489,11 @@ detection_pseudo_sn(struct noisechiselparams *p,
gal_data_t *worklab,
{
plabend = (plab=worklab->array) + worklab->size;
do
- if( *plab && ( area[*plab]<p->detsnminarea || brightness[*plab]<0) )
+ if( *plab!=GAL_BLANK_INT32
+ && ( area[*plab]<p->detsnminarea || brightness[*plab]<0) )
*plab=0;
while(++plab<plabend);
- worklab->name=PSN_EXTNAME;
+ worklab->name="PSEUDOS-FOR-SN";
gal_fits_img_write(worklab, p->detectionname, NULL,
PROGRAM_STRING);
worklab->name=NULL;
@@ -496,7 +503,7 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t
*worklab,
/* Calculate the signal to noise for successful detections: */
snarr=sn->array;
if(snind) indarr=snind->array;
- if(s0d1) { snarr[0]=NAN; if(snind) indarr[0]=GAL_BLANK_UINT32; }
+ if(s0d1) { snarr[0]=NAN; if(snind) indarr[0]=GAL_BLANK_INT32; }
for(i=1;i<tablen;++i)
{
ave=brightness[i]/area[i];
@@ -533,7 +540,7 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t
*worklab,
if(s0d1)
{
snarr[i]=NAN;
- if(snind) indarr[i]=GAL_BLANK_UINT32;;
+ if(snind) indarr[i]=GAL_BLANK_INT32;;
}
}
@@ -551,8 +558,8 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t
*worklab,
if(snind)
{
/* Make the comments, then write the table. */
- gal_linkedlist_add_to_stll(&comments, "See also: `"PSN_EXTNAME"' HDU "
- "of output with `--checkdetection'", 1);
+ gal_linkedlist_add_to_stll(&comments, "See also: `PSEUDOS-FOR-SN' "
+ "HDU of output with `--checkdetection'", 1);
gal_linkedlist_add_to_stll(&comments, s0d1 ? "Pseudo-detection S/N "
"over initially detected regions."
: "Pseudo-detection S/N over initially "
@@ -590,7 +597,7 @@ detection_pseudo_remove_low_sn(struct noisechiselparams *p,
size_t i;
float *snarr=sn->array;
uint8_t *b=workbin->array;
- uint32_t *l=worklab->array, *lf=l+worklab->size;
+ int32_t *l=worklab->array, *lf=l+worklab->size;
uint8_t *keep=gal_data_calloc_array(GAL_TYPE_UINT8, sn->size);
/* Specify the new labels for those that must be kept/changed. Note that
@@ -602,8 +609,13 @@ detection_pseudo_remove_low_sn(struct noisechiselparams *p,
/* Go over the pseudo-detection labels and only keep those that must be
- kept (using the new labels). */
- do *b++ = keep[ *l++ ] > 0; while(l<lf);
+ kept (using the new labels) in the binary array.. */
+ if( p->input->flag & GAL_DATA_FLAG_HASBLANK )
+ do
+ *b++ = *l == GAL_BLANK_INT32 ? GAL_BLANK_UINT8 : keep[ *l ] > 0;
+ while(++l<lf);
+ else
+ do *b++ = keep[ *l ] > 0; while(++l<lf);
/* If the user wanted to see the steps. */
@@ -638,6 +650,7 @@ detection_pseudo_real(struct noisechiselparams *p)
workbin=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
p->input->dsize, p->input->wcs, 0, p->cp.minmapsize,
NULL, NULL, NULL);
+ workbin->flag=p->input->flag;
/* Over the Sky: find the pseudo-detections and make the S/N table. */
@@ -651,7 +664,7 @@ detection_pseudo_real(struct noisechiselparams *p)
p->detsnthresh = *((float *)(quant->array));
if(!p->cp.quiet)
{
- asprintf(&msg, "Pseudo-det S/N: %.2f (%.3f quant of %zu).",
+ asprintf(&msg, "Pseudo-det S/N: %.2f (%.2f quant of %zu).",
p->detsnthresh, p->detquant, sn->size);
gal_timing_report(&t1, msg, 2);
free(msg);
@@ -703,8 +716,8 @@ detection_remove_false_initial(struct noisechiselparams *p,
{
size_t i;
uint8_t *b=workbin->array;
- uint32_t *l=p->olabel->array, *lf=l+p->olabel->size, curlab=1;
- uint32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32,
+ int32_t *l=p->olabel->array, *lf=l+p->olabel->size, curlab=1;
+ int32_t *newlabels=gal_data_calloc_array(GAL_TYPE_UINT32,
p->numinitialdets+1);
/* Find the new labels for all the existing labels. Recall that
@@ -719,7 +732,7 @@ detection_remove_false_initial(struct noisechiselparams *p,
array. They are not important so you can just ignore them. */
do
{
- if(*l)
+ if( *l && *l!=GAL_BLANK_INT32 )
{
newlabels[ *l ] =
newlabels[ *l ] /* Have we already checked this label? */
@@ -746,7 +759,7 @@ detection_remove_false_initial(struct noisechiselparams *p,
if(p->dilate)
do
{
- if(*l!=GAL_BLANK_UINT32)
+ if(*l!=GAL_BLANK_INT32)
*b = newlabels[ *l ] > 0;
++b;
}
@@ -754,7 +767,7 @@ detection_remove_false_initial(struct noisechiselparams *p,
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_UINT32) /* Sky and its STD. */
+ if(*l!=GAL_BLANK_INT32) /* Sky and its STD. */
*b = ( *l = newlabels[ *l ] ) > 0;
++b;
}
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 367048f..b22aa60 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -49,10 +49,10 @@ struct noisechiselparams
char *kernelname; /* Input kernel filename. */
char *khdu; /* Kernel HDU. */
uint8_t skysubtracted; /* Input has been Sky subtracted before. */
- float minbfrac; /* Undetected area min. frac. in tile. */
+ float minskyfrac; /* Undetected area min. frac. in tile. */
size_t minnumfalse; /* Min No. of det/seg for true quantile. */
- uint8_t onlydetect; /* Do not do any segmentation. */
+ uint8_t onlydetection; /* Do not do any segmentation. */
uint8_t grownclumps; /* Save grown clumps instead of original. */
uint8_t continueaftercheck; /* Don't abort after the check steps. */
@@ -77,8 +77,8 @@ struct noisechiselparams
uint8_t checksky; /* Check the Sky value estimation. */
size_t segsnminarea; /* Minimum area for segmentation. */
- float segquant; /* Quantile of clumps in sky for true S/N.*/
uint8_t checkclumpsn; /* Save the clump S/N values to a file. */
+ float segquant; /* Quantile of clumps in sky for true S/N.*/
uint8_t keepmaxnearriver; /* Keep clumps with a peak near a river. */
float gthresh; /* Multiple of STD to stop growing clumps.*/
size_t minriverlength; /* Min, len of good grown clump rivers. */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index 5520a3c..b9cd4d0 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -26,6 +26,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <stdlib.h>
#include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
#include <gnuastro/convolve.h>
#include <gnuastro-internal/timing.h>
@@ -61,6 +62,8 @@ noisechisel_convolve(struct noisechiselparams *p)
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("INPUT_CONVOLVED", &p->conv->name);
+
if(!p->cp.quiet) gal_timing_report(&t1, "Convolved with kernel.", 1);
if(p->detectionname)
@@ -140,6 +143,145 @@ noisechisel_convolve_correct_ch_edges(struct
noisechiselparams *p)
/***********************************************************************/
+/************* Output ***************/
+/***********************************************************************/
+/* The input image has been sky subtracted for further processing. So we'll
+ need to copy the input image directly into the output. */
+static void
+noisechisel_output_copy_input(struct noisechiselparams *p)
+{
+ int status=0;
+ fitsfile *in, *out;
+ char card[FLEN_CARD];
+
+
+ /* Create/open the output file. */
+ out=gal_fits_open_to_write(p->cp.output);
+
+
+ /* Open the input FITS file in the proper extension. */
+ in=gal_fits_hdu_open(p->inputname, p->cp.hdu, READWRITE);
+
+
+ /* Copy the input HDU into the output. */
+ if( fits_copy_hdu(in, out, 0, &status) )
+ gal_fits_io_error(status, "copying input hdu into first output hdu");
+
+
+ /* If an extension name exists in the input HDU, then don't touch it. If
+ the input doesn't have any, then make an EXTNAME keyword for it. Note
+ that `fits_read_card' will return a non-zero if it doesn't find the
+ keyword. */
+ if( fits_read_card(out, "EXTNAME", card, &status) )
+ {
+ status=0;
+ fits_write_key(out, TSTRING, "EXTNAME", "INPUT", "", &status);
+ }
+
+
+ /* Close the two files. */
+ fits_close_file(in, &status);
+ fits_close_file(out, &status);
+}
+
+
+
+
+
+/* Write the output file. */
+static void
+noisechisel_output(struct noisechiselparams *p)
+{
+ struct gal_fits_key_ll *keys=NULL;
+
+ /* Copy the input image into the first extension. */
+ noisechisel_output_copy_input(p);
+
+
+ /* Write the object labels and useful information into it's header. */
+ if(p->onlydetection)
+ {
+ p->olabel->name = "OBJECTS";
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_SIZE_T, "NDETS", 0,
+ &p->numobjects, 0, "Total number of detections",
+ 0, "counter");
+ }
+ else
+ {
+ p->olabel->name = "OBJECTS";
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_STRING, "WCLUMPS", 0, "yes", 0,
+ "Generate catalog with clumps?", 0, "bool");
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_SIZE_T, "NOBJS", 0,
+ &p->numobjects, 0, "Total number of objects",
+ 0, "counter");
+ }
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_FLOAT32, "DETSN", 0, &p->detsnthresh,
+ 0, "Minimum S/N of true pseudo-detections", 0,
+ "ratio");
+ gal_fits_img_write(p->olabel, p->cp.output, keys, PROGRAM_STRING);
+ p->olabel->name=NULL;
+ keys=NULL;
+
+
+ /* Write the clumps labels and useful information into it's header. Note
+ that to make the clumps image more easily viewable, we will set all
+ sky pixels to blank. Only clump pixels that have an overlapping object
+ pixel will be use anyway, so the sky pixels are irrelevant. */
+ if(p->onlydetection==0)
+ {
+ p->clabel->name="CLUMPS";
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_SIZE_T, "NCLUMPS", 0,
+ &p->numclumps, 0, "Total number of clumps", 0,
+ "counter");
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_FLOAT32, "CLUMPSN", 0,
+ &p->clumpsnthresh, 0,
+ "Minimum S/N of true clumps", 0, "ratio");
+ gal_fits_img_write(p->clabel, p->cp.output, keys, PROGRAM_STRING);
+ p->clabel->name=NULL;
+ keys=NULL;
+ }
+
+
+ /* 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, PROGRAM_STRING);
+ p->sky->name=NULL;
+
+
+ /* Write the Sky standard deviation into the output. */
+ p->std->name="SKY_STD";
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_FLOAT32, "MAXSTD", 0,
+ &p->maxstd, 0, "Maximum raw tile standard deviation",
+ 0, p->input->unit);
+ gal_fits_key_add_to_ll(&keys, GAL_TYPE_FLOAT32, "MINSTD", 0,
+ &p->minstd, 0, "Minimum raw tile standard deviation",
+ 0, p->input->unit);
+ gal_fits_key_add_to_ll(&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, PROGRAM_STRING);
+ p->std->name=NULL;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/***********************************************************************/
/************* High level function ***************/
/***********************************************************************/
void
@@ -158,9 +300,16 @@ noisechisel(struct noisechiselparams *p)
images. */
noisechisel_find_sky_subtract(p);
- /* Correct the convolved image channel edges if necessary. */
- noisechisel_convolve_correct_ch_edges(p);
+ /* If the user only wanted detection, ignore the segmentation steps. */
+ if(p->onlydetection==0)
+ {
+ /* Correct the convolved image channel edges if necessary. */
+ noisechisel_convolve_correct_ch_edges(p);
+
+ /* Do the segmentation. */
+ segmentation(p);
+ }
- /* Do the segmentation. */
- segmentation(p);
+ /* Write the output. */
+ noisechisel_output(p);
}
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index 6856742..5d4968b 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -56,7 +56,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
static void
segmentation_relab_noseg(struct clumps_thread_params *cltprm)
{
- uint32_t *olabel=cltprm->clprm->p->olabel->array;
+ int32_t *olabel=cltprm->clprm->p->olabel->array;
size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
do olabel[ *s ] = 1; while(++s<sf);
}
@@ -92,11 +92,11 @@ segmentation_relab_to_objects(struct clumps_thread_params
*cltprm)
float *imgss=p->input->array;
uint8_t *adjacency=adjacency_d->array;
size_t nngb=gal_dimension_num_neighbors(ndim);
- uint32_t *clumptoobj, *olabel=p->olabel->array;
+ int32_t *clumptoobj, *olabel=p->olabel->array;
size_t *dinc=gal_dimension_increment(ndim, dsize);
size_t *s, *sf, i, j, ii, rpnum, *nums=nums_d->array;
double ave, rpsum, c=sqrt(1/p->cpscorr), *sums=sums_d->array;
- uint32_t *ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb);
+ int32_t *ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb);
double err=cltprm->std*cltprm->std*(p->skysubtracted?1.0f:2.0f);
@@ -113,7 +113,7 @@ segmentation_relab_to_objects(struct clumps_thread_params
*cltprm)
if( olabel[ *s ]==CLUMPS_RIVER )
{
/* Initialize the values. */
- ii=0;
+ i=ii=0;
rpnum=1; /* River-pixel number of points used. */
rpsum=imgss[*s]; /* River-pixel sum of values used. */
memset(ngblabs, 0, nngb*sizeof *ngblabs);
@@ -121,7 +121,7 @@ segmentation_relab_to_objects(struct clumps_thread_params
*cltprm)
/* Check all the fully-connected neighbors of this pixel and see if
it touches a label or not */
GAL_DIMENSION_NEIGHBOR_OP(*s, ndim, dsize, ndim, dinc, {
- if( olabel[nind] && olabel[nind] < CLUMPS_MAXLAB )
+ if( olabel[nind] > 0 )
{
/* Add this neighbor's value and increment the number. */
if( !isnan(imgss[nind]) ) { ++rpnum; rpsum+=imgss[nind]; }
@@ -189,23 +189,26 @@ segmentation_relab_to_objects(struct clumps_thread_params
*cltprm)
/* For a check:
- printf("=====================\n");
- printf("%zu:\n--------\n", cltprm->id);
- for(i=1;i<amwidth;++i)
+ if(cltprm->id==XXX)
{
- printf(" %zu...\n", i);
- for(j=1;j<amwidth;++j)
+ printf("=====================\n");
+ printf("%zu:\n--------\n", cltprm->id);
+ for(i=1;i<amwidth;++i)
{
- ii=i*amwidth+j;
- if(nums[ii])
+ printf(" %zu...\n", i);
+ for(j=1;j<amwidth;++j)
{
- ave=sums[ii]/nums[ii];
- printf(" ...%zu: N:%-4zu S:%-10.2f S/N: %-10.2f "
- "--> %u\n", j, nums[ii], sums[ii], c*ave/sqrt(ave+err),
- adjacency[ii]);
+ ii=i*amwidth+j;
+ if(nums[ii])
+ {
+ ave=sums[ii]/nums[ii];
+ printf(" ...%zu: N:%-4zu S:%-10.2f S/N: %-10.2f "
+ "--> %u\n", j, nums[ii], sums[ii],
+ c*ave/sqrt(ave+err), adjacency[ii]);
+ }
}
+ printf("\n");
}
- printf("\n");
}
*/
@@ -215,11 +218,21 @@ segmentation_relab_to_objects(struct clumps_thread_params
*cltprm)
&cltprm->numobjects);
clumptoobj = cltprm->clumptoobj->array;
+ /* For a check
+ if(cltprm->id==XXXX)
+ {
+ printf("NUMTRUECLUMPS: %zu\n----------\n", cltprm->numtrueclumps);
+ for(i=0;i<cltprm->numtrueclumps+1;++i)
+ printf("\t%zu --> %d\n", i, clumptoobj[i]);
+ printf("=== numobjects: %zu====\n", cltprm->numobjects);
+ exit(0);
+ }
+ */
/* Correct all the labels. */
sf=(s=cltprm->indexs->array)+cltprm->indexs->size;
do
- if( olabel[*s]<CLUMPS_MAXLAB )
+ if( olabel[*s] > 0 )
olabel[*s] = clumptoobj[ olabel[*s] ];
while(++s<sf);
@@ -243,18 +256,18 @@ segmentation_relab_clumps_in_objects(struct
clumps_thread_params *cltprm)
{
size_t numobjects=cltprm->numobjects, numtrueclumps=cltprm->numtrueclumps;
- uint32_t *clumptoobj=cltprm->clumptoobj->array;
- uint32_t *clabel=cltprm->clprm->p->clabel->array;
+ int32_t *clumptoobj=cltprm->clumptoobj->array;
+ int32_t *clabel=cltprm->clprm->p->clabel->array;
size_t i, *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
size_t *nclumpsinobj=gal_data_calloc_array(GAL_TYPE_SIZE_T, numobjects+1);
- uint32_t *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32, numtrueclumps+1);
+ int32_t *newlabs=gal_data_calloc_array(GAL_TYPE_UINT32, numtrueclumps+1);
/* Fill both arrays. */
for(i=1;i<numtrueclumps+1;++i)
newlabs[i] = ++nclumpsinobj[ clumptoobj[i] ];
/* Reset the clump labels over the detection region. */
- do if(clabel[*s]) clabel[*s] = newlabs[ clabel[*s] ]; while(++s<sf);
+ do if(clabel[*s]>0) clabel[*s] = newlabs[ clabel[*s] ]; while(++s<sf);
/* Clean up. */
free(newlabs);
@@ -276,7 +289,7 @@ static void
segmentation_relab_overall(struct clumps_thread_params *cltprm)
{
struct clumps_params *clprm=cltprm->clprm;
- uint32_t startinglab, *olabel=clprm->p->olabel->array;
+ int32_t startinglab, *olabel=clprm->p->olabel->array;
size_t *s=cltprm->indexs->array, *sf=s+cltprm->indexs->size;
/* Lock the mutex if we are working on more than one thread. NOTE: it is
@@ -332,7 +345,7 @@ segmentation_on_threads(void *in_prm)
size_t i, *s, *sf;
gal_data_t *topinds;
struct clumps_thread_params cltprm;
- uint32_t *clabel=p->clabel->array, *olabel=p->olabel->array;
+ int32_t *clabel=p->clabel->array, *olabel=p->olabel->array;
/* Initialize the general parameters for this thread. */
cltprm.clprm = clprm;
@@ -433,7 +446,7 @@ segmentation_on_threads(void *in_prm)
{
sf=(s=cltprm.indexs->array)+cltprm.indexs->size;
do
- if(olabel[*s]<CLUMPS_MAXLAB) clabel[*s]=olabel[*s];
+ if(olabel[*s]>0) clabel[*s]=olabel[*s];
while(++s<sf);
}
@@ -474,7 +487,7 @@ segmentation_on_threads(void *in_prm)
if(cltprm.numobjects>1)
segmentation_relab_clumps_in_objects(&cltprm);
gal_data_free(cltprm.clumptoobj);
- if(clprm->step==6) continue;
+ if(clprm->step==6) {continue;}
}
/* Convert the object labels to their final value */
@@ -497,7 +510,7 @@ segmentation_save_sn_table(struct clumps_params *clprm)
{
char *msg;
float *sarr;
- uint32_t *oiarr, *cioarr;
+ int32_t *oiarr, *cioarr;
size_t i, j, c=0, totclumps=0;
gal_data_t *sn, *objind, *clumpinobj;
struct noisechiselparams *p=clprm->p;
@@ -514,10 +527,10 @@ segmentation_save_sn_table(struct clumps_params *clprm)
sn=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &totclumps, NULL, 0,
p->cp.minmapsize, "CLUMP_S/N", "ratio",
"Signal-to-noise ratio.");
- objind=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &totclumps, NULL, 0,
+ objind=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &totclumps, NULL, 0,
p->cp.minmapsize, "HOST_DET_ID", "counter",
"ID of detection hosting this clump.");
- clumpinobj=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &totclumps, NULL, 0,
+ clumpinobj=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &totclumps, NULL, 0,
p->cp.minmapsize, "CLUMP_ID_IN_OBJ", "counter",
"ID of clump in host detection.");
@@ -576,10 +589,8 @@ static void
segmentation_detections(struct noisechiselparams *p)
{
char *msg;
- uint8_t *b;
- uint32_t *l, *lf;
struct clumps_params clprm;
- gal_data_t *tmp, *demo, *labindexs, *claborig;
+ gal_data_t *labindexs, *claborig, *demo=NULL;
/* Get the indexs of all the pixels in each label. */
@@ -726,13 +737,7 @@ segmentation_detections(struct noisechiselparams *p)
default values are hard to view, so we'll make a copy of the
demo, set all Sky regions to blank and all clump macro values
to zero. */
- tmp=gal_data_copy(demo);
- b=p->binary->array; lf=(l=tmp->array)+tmp->size;
- do *l = ( *b++
- ? ( *l > CLUMPS_MAXLAB ? 0 : *l )
- : GAL_BLANK_UINT32 ); while(++l<lf);
- gal_fits_img_write(tmp, p->segmentationname, NULL, PROGRAM_STRING);
- gal_data_free(tmp);
+ gal_fits_img_write(demo, p->segmentationname, NULL, PROGRAM_STRING);
/* If the user wanted to check the clump S/N values, then break
out of the loop, we don't need the rest of the process any
@@ -796,7 +801,7 @@ segmentation(struct noisechiselparams *p)
{
float *f;
char *msg;
- uint32_t *l, *lf;
+ int32_t *l, *lf;
struct timeval t1;
/* To keep the user up to date. */
@@ -807,7 +812,6 @@ segmentation(struct noisechiselparams *p)
1);
}
-
/* If a check segmentation image was requested, then put in the
inputs. */
if(p->segmentationname)
@@ -820,16 +824,16 @@ segmentation(struct noisechiselparams *p)
p->olabel->name=NULL;
}
-
/* Allocate the clump labels image. */
p->clabel=gal_data_alloc(NULL, p->olabel->type, p->olabel->ndim,
p->olabel->dsize, p->olabel->wcs, 1,
p->cp.minmapsize, NULL, NULL, NULL);
+ p->clabel->flag=p->input->flag;
/* Set any possibly existing NaN values to blank. */
f=p->input->array; lf=(l=p->clabel->array)+p->clabel->size;
- do if(isnan(*f++)) *l=GAL_BLANK_UINT32; while(++l<lf);
+ do if(isnan(*f++)) *l=GAL_BLANK_INT32; while(++l<lf);
/* Find the clump S/N threshold. */
@@ -838,12 +842,13 @@ segmentation(struct noisechiselparams *p)
/* Reset the clabel array to find true clumps in objects. */
f=p->input->array; lf=(l=p->clabel->array)+p->clabel->size;
- do *l = isnan(*f++) ? GAL_BLANK_UINT32 : 0; while(++l<lf);
+ do *l = isnan(*f++) ? GAL_BLANK_INT32 : 0; while(++l<lf);
/* Find true clumps over the detected regions. */
segmentation_detections(p);
+
/* Report the results and timing to the user. */
if(!p->cp.quiet)
{
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index cfa009f..cf97484 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -93,7 +93,7 @@ sky_mean_std_undetected(void *in_prm)
/* Only continue, if the fraction of Sky values are less than the
requested fraction. */
- if( (float)(numsky)/(float)(tile->size) > p->minbfrac)
+ if( (float)(numsky)/(float)(tile->size) > p->minskyfrac)
{
/* Calculate the mean and STD over this tile. */
s=s2=0.0f;
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index c510bc8..6d38004 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -53,26 +53,39 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/**********************************************************************/
/*************** Apply a given threshold. *****************/
/**********************************************************************/
-void
-threshold_apply(struct noisechiselparams *p, float *value1, float *value2,
- int type)
+struct threshold_apply_p
{
- size_t tid;
+ float *value1;
+ float *value2;
+ int kind;
+ struct noisechiselparams *p;
+};
+
+
+
+
+/* Apply the threshold on the tiles given to this thread. */
+static void *
+threshold_apply_on_thread(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct threshold_apply_p *taprm=(struct threshold_apply_p *)(tprm->params);
+ struct noisechiselparams *p=taprm->p;
+
+ size_t i, tid;
void *tarray=NULL;
gal_data_t *tile, *tblock=NULL;
+ float *value1=taprm->value1, *value2=taprm->value2;
- /* Clear the binary array (this is mainly because the input may contain
- blank values and we won't be doing the thresholding no those
- pixels. */
- memset(p->binary->array, 0, p->binary->size);
-
- /* Go over all the tiles. */
- for(tid=0; tid<p->cp.tl.tottiles; ++tid)
+ /* Go over all the tiles assigned to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
{
/* For easy reading. */
+ tid=tprm->indexs[i];
tile=&p->cp.tl.tiles[tid];
- switch(type)
+ /* Based on the kind of threshold. */
+ switch(taprm->kind)
{
/* This is a quantile threshold. */
@@ -86,12 +99,23 @@ threshold_apply(struct noisechiselparams *p, float *value1,
float *value2,
tile->block=p->conv;
}
- /* Apply the threshold. */
+ /* Apply the threshold: When the `>' comparison fails, it can be
+ either because the pixel was actually smaller than the
+ threshold, or that it was a NaN value. In the first case,
+ return 0, in the second, return a blank value.
+
+ We already know if a tile contains a blank value (which is a
+ constant over the whole loop). So before checking if the value
+ is blank, see if the tile actually has a blank value. This
+ will help in efficiency, because the compiler can move this
+ check out of the loop and only check for NaN values when we
+ know the tile has blank pixels. */
GAL_TILE_PARSE_OPERATE({
*o = ( *i > value1[tid]
? ( *i > value2[tid] ? THRESHOLD_NO_ERODE_VALUE : 1 )
- : 0 );
- }, tile, p->binary, 1, 1);
+ : ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
+ ? GAL_BLANK_UINT8 : 0 ) );
+ }, tile, p->binary, 1, 0);
/* Revert the tile's pointers back to what they were. */
if(p->conv) { tile->array=tarray; tile->block=tblock; }
@@ -101,19 +125,14 @@ threshold_apply(struct noisechiselparams *p, float
*value1, float *value2,
/* This is a Sky and Sky STD threshold. */
case THRESHOLD_SKY_STD:
- /* The threshold is always low. So for the majority of non-NaN
- pixels in the image, the condition above will be true. If we
- come over a NaN pixel, then by definition of NaN, all
- conditionals will fail.
-
- If an image doesn't have any NaN pixels, only the pixels below
- the threshold have to be checked for a NaN which are by
- definition a very small fraction of the total pixels. And if
- there are NaN pixels in the image. */
+ /* See the explanation above the same step in the quantile
+ threshold for an explanation. */
GAL_TILE_PARSE_OPERATE({
*o = ( ( *i - value1[tid] > p->dthresh * value2[tid] )
- ? 1 : *i==*i ? 0 : GAL_BLANK_UINT8 );
- }, tile, p->binary, 1, 1);
+ ? 1
+ : ( (tile->flag & GAL_DATA_FLAG_HASBLANK) && !(*i==*i)
+ ? GAL_BLANK_UINT8 : 0 ) );
+ }, tile, p->binary, 1, 0);
break;
@@ -121,15 +140,32 @@ threshold_apply(struct noisechiselparams *p, float
*value1, float *value2,
error(EXIT_FAILURE, 0, "a bug! please contact us at %s so we can "
"address the problem. For some reason a value of %d had "
"been given to `type' in `threshold_apply'",
- PACKAGE_BUGREPORT, type);
+ PACKAGE_BUGREPORT, taprm->kind);
}
}
+
+ /* Wait until all the other threads finish. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
}
+/* Apply a given threshold threshold on the tiles. */
+void
+threshold_apply(struct noisechiselparams *p, float *value1,
+ float *value2, int kind)
+{
+ struct threshold_apply_p taprm={value1, value2, kind, p};
+ gal_threads_spin_off(threshold_apply_on_thread, &taprm, p->cp.tl.tottiles,
+ p->cp.numthreads);
+}
+
+
+
+
@@ -235,11 +271,13 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
(*first)->next=(*second)->next=NULL;
if(filename)
{
+ (*first)->name="THRESH1_INTERP";
+ (*second)->name="THRESH2_INTERP";
gal_tile_full_values_write(*first, tl, filename, PROGRAM_STRING);
gal_tile_full_values_write(*second, tl, filename, PROGRAM_STRING);
+ (*first)->name = (*second)->name = NULL;
}
-
/* Smooth the threshold if requested. */
if(p->smoothwidth>1)
{
@@ -258,8 +296,11 @@ threshold_interp_smooth(struct noisechiselparams *p,
gal_data_t **first,
/* Add them to the check image. */
if(filename)
{
+ (*first)->name="THRESH1_SMOOTH";
+ (*second)->name="THRESH2_SMOOTH";
gal_tile_full_values_write(*first, tl, filename, PROGRAM_STRING);
gal_tile_full_values_write(*second, tl, filename, PROGRAM_STRING);
+ (*first)->name = (*second)->name = NULL;
}
}
}
@@ -421,27 +462,37 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
/* Allocate space for the quantile threshold values. */
qprm.erode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
- "QTHRESH-ERODE", p->input->unit, NULL);
+ NULL, p->input->unit, NULL);
qprm.noerode_th=gal_data_alloc(NULL, p->input->type, p->input->ndim,
tl->numtiles, NULL, 0, cp->minmapsize,
- "QTHRESH-NOERODE", p->input->unit, NULL);
+ NULL, p->input->unit, NULL);
/* Allocate temporary space for processing in each tile. */
qprm.usage=gal_data_malloc_array(p->input->type,
cp->numthreads * p->maxtcontig);
- /* Find the threshold on each tile, then clean up the temporary space. */
+
+ /* Find the threshold on each tile, free the temporary processing space
+ and set the blank flag. */
qprm.p=p;
gal_threads_spin_off(qthresh_on_tile, &qprm, tl->tottiles, cp->numthreads);
+ free(qprm.usage);
+ if( gal_blank_present(qprm.erode_th) )
+ {
+ qprm.erode_th->flag |= GAL_DATA_FLAG_HASBLANK;
+ qprm.noerode_th->flag |= GAL_DATA_FLAG_HASBLANK;
+ }
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,
PROGRAM_STRING);
gal_tile_full_values_write(qprm.noerode_th, tl, p->qthreshname,
PROGRAM_STRING);
+ qprm.erode_th->name = qprm.noerode_th->name = NULL;
}
- free(qprm.usage);
/* Interpolate and smooth the derived values. */
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 25f739c..e8d2a8e 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -30,6 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
#include <gnuastro/threads.h>
#include <gnuastro/dimension.h>
@@ -470,7 +471,7 @@ ui_prepare_kernel(struct noisechiselparams *p)
static void
ui_prepare_tiles(struct noisechiselparams *p)
{
- gal_data_t *check;
+ gal_data_t *check, *tile;
struct gal_tile_two_layer_params *tl=&p->cp.tl, *ltl=&p->ltl;
@@ -507,6 +508,21 @@ ui_prepare_tiles(struct noisechiselparams *p)
memcpy(p->maxltsize, check->dsize, ltl->ndim*sizeof *p->maxltsize);
}
+
+ /* If the input has blank elements, then go over all the tiles and set
+ their blank flag to 1 if they have blank values within
+ them. Afterwards, set the use-zero flag for the tiles in any case.*/
+ if( p->input->flag & GAL_DATA_FLAG_HASBLANK )
+ {
+ gal_tile_full_blank_flag(tl->tiles, p->cp.numthreads);
+ gal_tile_full_blank_flag(ltl->tiles, p->cp.numthreads);
+ }
+ for(tile=tl->tiles;tile!=NULL;tile=tile->next)
+ tile->flag |= GAL_DATA_FLAG_BLANK_CH;
+ for(tile=ltl->tiles;tile!=NULL;tile=tile->next)
+ tile->flag |= GAL_DATA_FLAG_BLANK_CH;
+
+
/* Make the tile check image if requested. */
if(tl->checktiles)
{
@@ -551,6 +567,20 @@ ui_preparations(struct noisechiselparams *p)
gal_checkset_allocate_copy("INPUT", &p->input->name);
+ /* NoiseChisel currently only works on 2D datasets (images). */
+ if(p->input->ndim!=2)
+ error(EXIT_FAILURE, 0, "%s (hdu: %s) has %zu dimensions but NoiseChisel "
+ "can only operate on 2D datasets (images)", p->inputname, p->cp.hdu,
+ 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. */
+ if( gal_blank_present(p->input) ) p->input->flag |= GAL_DATA_FLAG_HASBLANK;
+ p->input->flag |= GAL_DATA_FLAG_BLANK_CH;
+
+
/* Read in the kernel for convolution. */
ui_prepare_kernel(p);
@@ -563,9 +593,10 @@ ui_preparations(struct noisechiselparams *p)
p->binary=gal_data_alloc(NULL, GAL_TYPE_UINT8, p->input->ndim,
p->input->dsize, p->input->wcs, 0,
p->cp.minmapsize, NULL, "binary", NULL);
- p->olabel=gal_data_alloc(NULL, GAL_TYPE_UINT32, p->input->ndim,
+ p->olabel=gal_data_alloc(NULL, GAL_TYPE_INT32, p->input->ndim,
p->input->dsize, p->input->wcs, 0,
p->cp.minmapsize, NULL, "labels", NULL);
+ p->binary->flag = p->olabel->flag = p->input->flag;
}
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 513304a..ca28276 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -38,7 +38,7 @@ enum option_keys_enum
ARGS_OPTION_KEY_LARGETILESIZE = 'L',
ARGS_OPTION_KEY_KERNEL = 'k',
ARGS_OPTION_KEY_SKYSUBTRACTED = 'E',
- ARGS_OPTION_KEY_MINBFRAC = 'B',
+ ARGS_OPTION_KEY_MINSKYFRAC = 'B',
ARGS_OPTION_KEY_MIRRORDIST = 'r',
ARGS_OPTION_KEY_MODMEDQDIFF = 'Q',
ARGS_OPTION_KEY_QTHRESH = 't',
@@ -62,7 +62,7 @@ enum option_keys_enum
automatically). */
ARGS_OPTION_KEY_KHDU = 1000,
ARGS_OPTION_KEY_MINNUMFALSE,
- ARGS_OPTION_KEY_ONLYDETECT,
+ ARGS_OPTION_KEY_ONLYDETECTION,
ARGS_OPTION_KEY_GROWNCLUMPS,
ARGS_OPTION_KEY_SMOOTHWIDTH,
ARGS_OPTION_KEY_CHECKQTHRESH,
diff --git a/bin/statistics/ui.c b/bin/statistics/ui.c
index ca8a6f0..46d6607 100644
--- a/bin/statistics/ui.c
+++ b/bin/statistics/ui.c
@@ -841,6 +841,8 @@ ui_preparations(struct statisticsparams *p)
{
/* Only keep the elements we want. */
gal_blank_remove(p->input);
+ p->input->flag &= ~GAL_DATA_FLAG_HASBLANK ;
+ p->input->flag |= GAL_DATA_FLAG_BLANK_CH ;
/* Make sure there is data remaining: */
if(p->input->size==0)
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a9fa45d..6cb04e2 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -441,8 +441,10 @@ NoiseChisel
Invoking NoiseChisel
-* NoiseChisel options:: The options specific to NoiseChisel.
-* NoiseChisel output:: NoiseChisel's output. Why no catalog?
+* General NoiseChisel options:: General NoiseChisel preprocessing.
+* Detection options:: Configure detection in NoiseChisel.
+* Segmentation options:: Configure segmentation in NoiseChisel.
+* NoiseChisel output:: NoiseChisel's output format.
MakeCatalog
@@ -11908,46 +11910,73 @@ one for the Sky standard deviation).
@node NoiseChisel, MakeCatalog, Statistics, Data analysis
@section NoiseChisel
-Once raw data have gone through the initial reduction process (through the
-programs in @ref{Data manipulation}). We are ready to derive scientific
-results out of them. Unfortunately in most cases, the scientifically
-interesting targets are deeply drowned in a sea of noise. NoiseChisel is
-Gnuastro's tool to detect signal in noise. In fact, NoiseChisel was the
-motivation behind creating Gnuastro and has a journal article devoted to
-its techniques: @url{http://arxiv.org/abs/1505.01664, arXiv:1505.01664},
-published in 2015 by the Astrophysical Journal Supplement Series
-220.
-
-Following the explanations for the options in @ref{NoiseChisel options}
-should provide a relatively good idea of how NoiseChisel works, but it is
-recommended to see the paper since it does a very thorough job at
-explaining the concepts and methods of NoiseChisel with abundant
-demonstrations for each step. However, the paper cannot undergo any further
-updates, but NoiseChisel will evolve: better algorithms or steps will be
-found, thus options will be added or removed. So this section is the
-definitive guide. Please follow the @file{NEWS} file with each release for
-such updates.
-
address@hidden Labeling
@cindex Detection
-Detection is the process of separating signal from noise. In other
-words, after detection is complete, one set of data elements (pixels
-in an image) will be distinguished as signal and another set of the
-data elements will be noise. Segmentation is the process of labeling
-the detected pixels into possibly multiple components (objects). For
-example when two galaxies lie sufficiently close to each other to be
-detected as one object.
-
address@hidden Noise-based detection
-NoiseChisel was the first software to make use of a noise-based concept to
-detection and segmentation. In this method, instead of emphasizing on the
-signal and trying to guess the properties of the to-be-detected targets
-prior to detection (for example assuming that it is an ellipse), the
-emphasis is put on the noise in the image and it imposes statistically
-negligible requirements on the signal. The name of NoiseChisel is derived
-from the first thing it does after thresholding the image: to erode it. In
-mathematical morphology, erosion on pixels can be pictured as carving off
-boundary pixels. So what NoiseChisel does is similar to what a wood chisel
-or stone chisel do. It is just not a hardware but software.
address@hidden Segmentation
+Once instrumental signatures are removed from the raw data in the initial
+reduction process (see @ref{Data manipulation}). We are ready to derive
+scientific results out of them. But we can't do anything special with a raw
+dataset, for example an image is just an array of values. Every pixel just
+has one value and its position within the image. Therefore, the first step
+of your high-level analysis will be to classify/label the dataset
+elements/pixels into two classes: signal and noise. This process is
+formally known as @emph{detection}. Afterwards, you want to separate the
+detections into multiple components (for example when two detected regions
+aren't touching, they should be treated independently as two distant
+galaxies for example). This higher level classification of the detections
+is known as @emph{segmentation}. NoiseChisel is Gnuastro's program for
+detection and segmentation.
+
+NoiseChisel works based on a new noise-based approach to signal detection
+and was introduced to the astronomical community in
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa
+[2015]}. NoiseChisel's primary output is an array (image) with the same
+size as the input but containing labels: those pixels with a label of 0 are
+noise/sky while those pixels with labels larger than 0 are detections
+(separate segments will be given positive integers, starting from 1). For
+more on NoiseChisel's particular output format and its benefits (especially
+in conjunction with @ref{MakeCatalog}), please see
address@hidden://arxiv.org/abs/1611.06387, Akhlaghi [2016]}.
+
+Data is inherently mixed with noise: only mock/simulated datasets are free
+of noise. So this process of separating signal from noise is not
+trivial. In particular, most scientifically interesting astronomical
+targets are faint, can have a large variety of morphologies along with a
+large distribution in brightness and size which are all drowned in a ocean
+of noise. So detection is a uniquely vital aspect of any scientific work
+and even more so for astronomical research. This is such a fundamental step
+that designing of NoiseChisel was the primary motivation behind creating
+Gnuastro: the first generation of Gnuastro's programs were all first part
+of what later became NoiseChisel, afterwards they spinned-off into separate
+programs.
+
+Before using NoiseChisel it is strongly recommended to read
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]} to
+gain a good understanding of what it does and how each parameter influences
+the output. Thanks to that paper, there is no more need to continue this
+introduction any further and we can just dive into the details of running
+NoiseChisel in the following sections. However, the paper cannot undergo
+any further updates, but NoiseChisel will evolve: better algorithms or
+steps will be found, thus options will be added or removed. So this section
+is the definitive guide to the options. Please follow the
address@hidden@footnote{The @file{NEWS} file is in the released Gnuastro
+tarball (see @ref{Release tarball}). You can also see it online at
address@hidden://git.savannah.gnu.org/cgit/gnuastro.git/plain/NEWS}.} file with
+each release to see how they have changed since the publication of that
+paper.
+
address@hidden Erosion
+The name of NoiseChisel is derived from the first thing it does after
+thresholding the dataset: to erode it. In mathematical morphology, erosion
+on pixels can be pictured as carving off boundary pixels. Hence, what
+NoiseChisel does is similar to what a wood chisel or stone chisel do. It is
+just not a hardware, but a software. Infact looking at it as a chisel and
+your dataset as a solid cube of rock will greatly help in best using it:
+with NoiseChisel you literally carve the galaxies/stars/comets out of the
+noise. Try running it with the @option{--checkdetection} option to see each
+step of the carving process on your input dataset. You can then change a
+specific option to carve out your signal out of the noise more
+successfully.
@menu
* Invoking astnoisechisel:: Options and arguments for NoiseChisel.
@@ -11956,7 +11985,9 @@ or stone chisel do. It is just not a hardware but
software.
@node Invoking astnoisechisel, , NoiseChisel, NoiseChisel
@subsection Invoking NoiseChisel
-NoiseChisel will detect signal in noise. The executable name is
+NoiseChisel will detect and segment signal in noise producing a
+multi-extension labeled image, ready for input into @ref{MakeCatalog} to
+generate a catalog or other processing. The executable name is
@file{astnoisechisel} with the following general template
@example
@@ -11968,142 +11999,208 @@ One line examples:
@example
$ astnoisechisel input.fits
-$ astnoisechisel --nch1=4 --lmesh=256 input.fits
+$ astnoisechisel --numchannels=4,1 --tilesize=100,100 input.fits
@end example
@cindex Gaussian
@noindent
-If NoiseChisel is to do processing, an input image should be provided with
-the recognized extensions as input data, see @ref{Arguments}. The options
-that are shared by all Gnuastro programs can be seen in @ref{Common
-options}. In order to ignore some pixels during the analysis, you can use
-Gnuastro's @ref{Arithmetic} program to set them as blank (see @ref{Blank
-pixels}). Any blank pixel is ignored by all the steps of NoiseChisel.
-
-A convolution kernel can also be optionally given. If a value (file
-name) is given to @option{--kernel} on the command-line or in a
-configuration file (see @ref{Configuration files}), then that file
-will be used to convolve the image prior to thresholding. Otherwise a
-default kernel will be used. The default kernel is a 2D Gaussian with
-a FWHM of 2 pixels truncated at 5 times the FWHM. See Section 3.1.1 of
-Akhlaghi and Ichikawa (2015) to learn why this particular kernel was
-chosen as default. See @ref{Convolution kernel} for kernel related
-options.
-
-NoiseChisel uses a mesh grid to tile the image. This enables it to deal
-with possible gradients, see @ref{Tessellation}. The mesh grid options are
-common to all the programs using it, and are listed in @ref{Processing
-options}. In particular, NoiseChisel uses two mesh grids: a large and a
-small one. The sizes of the meshs in each grid can be specified with the
-following two options (the @option{--meshsize} option is not recognized by
-NoiseChisel). The rest of the options explained in @ref{Processing options}
-apply to both grids.
-
address@hidden @option
-
address@hidden -s
address@hidden --smeshsize
-Similar to @option{--meshsize} in @ref{Processing options}, but for the
-smaller mesh grid, which is most dependent on the gradients in the
-image.
-
address@hidden -l
address@hidden --lmeshsize
-Similar to @option{--meshsize} in @ref{Processing options}, but for the
-larger mesh grid, used for detection and segmentation Signal to noise
-ratio analysis.
-
address@hidden table
+If NoiseChisel is to do processing (for example you don't want to get help,
+or see the values to each input parameter), an input image should be
+provided with the recognized extensions (see @ref{Arguments}). NoiseChisel
+shares a large set of common operations with other Gnuastro programs,
+mainly regarding input/output, general processing steps, and general
+operating modes. To help in a unified experience between all of Gnuastro's
+programs, these operations have the same command-line options, see
address@hidden options} for a full list. Since the common options are
+thoroughly discussed there, they are no longer reviewed here. You can see
+all the options with a short description on the command-line with the
address@hidden option, see @ref{Getting help}.
+
+NoiseChisel's input image may contain blank elements (see @ref{Blank
+pixels}). Blank elements will be ignored in all steps of NoiseChisel. Hence
+if your dataset has bad pixels which should be masked with a mask image,
+please use Gnuastro's @ref{Arithmetic} program (in particular its
address@hidden operator) to convert those pixels to blank pixels before
+running NoiseChisel. Gnuastro's Arithmetic program has bitwise operators
+helping you select specific kinds of bad-pixels when necessary.
+
+A convolution kernel can also be optionally given. If a value (file name)
+is given to @option{--kernel} on the command-line or in a configuration
+file (see @ref{Configuration files}), then that file will be used to
+convolve the image prior to thresholding. Otherwise a default kernel will
+be used. The default kernel is a 2D Gaussian with a FWHM of 2 pixels
+truncated at 5 times the FWHM. This choice of the default kernel is
+discussed in Section 3.1.1 of @url{https://arxiv.org/abs/1505.01664,
+Akhlaghi and Ichikawa [2015]}. See @ref{Convolution kernel} for kernel
+related options.
+
+NoiseChisel defines two tessellations over the input (see
address@hidden). This enables it to deal with possible gradients in the
+input dataset and also significantly improve speed by processing each tile
+on different threads. The tessellation related options are discussed in
address@hidden options}. In particular, NoiseChisel uses two tessellations
+(with everything between them identical except the tile sizes): a
+fine-grained one with smaller tiles (mainly used in detection) and a more
+larger tiled one which is used for multi-threaded processing. The common
+Tessellation options discribed in @ref{Processing options} define all
+parameters of both tessellations, only the large tile size for the latter
+tessellation is set through the @option{--largetilesize} option. To inspect
+the tessellations on your input dataset, run NoiseChisel with
address@hidden
address@hidden
@noindent
-Please run NoiseChisel with the @option{--help} option to list all the
-recognized options with a short explanation, irrespective of which
-part of the Gnuastro book they are fully explained in, see
address@hidden help}.
address@hidden TIP:} Frequently use the options starting with
address@hidden Depending on what you want to detect in the data, you can
+often play with the paramters/options for a better result than the default
+parameters. You can start with @option{--checkdetection} and
address@hidden for the main steps. For their full list please
+run:
address@hidden
+$ astnoisechisel --help | grep check
address@hidden example
address@hidden cartouche
+In the sections below, NoiseChisel's options are classified into three
+general classes to help in easy navigation. @ref{General NoiseChisel
+options} mainly discusses the options relating to input and those that are
+shared in both detection and segmentation. Options to configure the
+detection are described in @ref{Detection options} and @ref{Segmentation
+options} we discuss how you can fine-tune the segmentation of the
+detections. Finally in @ref{NoiseChisel output} the format of NoiseChisel's
+output is discussed. The order of options here follow the same logical
+order that the respective action takes place within NoiseChisel (note that
+the output of @option{--help} is sorted alphabetically).
@menu
-* NoiseChisel options:: The options specific to NoiseChisel.
-* NoiseChisel output:: NoiseChisel's output. Why no catalog?
+* General NoiseChisel options:: General NoiseChisel preprocessing.
+* Detection options:: Configure detection in NoiseChisel.
+* Segmentation options:: Configure segmentation in NoiseChisel.
+* NoiseChisel output:: NoiseChisel's output format.
@end menu
address@hidden NoiseChisel options, NoiseChisel output, Invoking
astnoisechisel, Invoking astnoisechisel
address@hidden NoiseChisel options
-The options particular to NoiseChisel are listed below. They are
-classified by context and also sorted in the same order that the
-operations are done on the image. See Akhlaghi and Ichikawa (2015) for
-a very complete, detailed and illustrated explanation of each
-step. Reading through the option explanations should be enough to
-obtain a general idea of how NoiseChisel works. Before the procedures
-explained by these options begin, the image is convolved with a
-kernel. The first group of options listed below are those that apply
-to both the detection and the segmentation processes.
address@hidden General NoiseChisel options, Detection options, Invoking
astnoisechisel, Invoking astnoisechisel
address@hidden General NoiseChisel options
+
+The options discussed in this section are mainly regarding the input(s),
+output, and some general processing options that are shared between both
+detection and segmentation. Recall that you can always see the full list of
+Gnuastro's options with the @option{--help} option.
+
@table @option
address@hidden -k STR
address@hidden --kernel=STR
+File name of kernel to smooth the image before applying the threshold, see
address@hidden kernel}. The first step of NoiseChisel is to
+convolve/smooth the image and use the convolved image in multiple steps
+during the processing. It will be used to define (and later apply) the
+quantile threshold (see @option{--qthresh}). The convolved image is also
+used to define the clumps (see Section 3.2.1 and Figure 8 of
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa [2015]}).
+
+The @option{--kernel} option is not mandatory. If no kernel is provided, a
+2D Gaussian profile with a FWHM of 2 pixels truncated at 5 times the FWHM
+is used. This choice of the default kernel is discussed in Section 3.1.1 of
+Akhlaghi and Ichikawa [2015].
+
address@hidden --khdu=STR
+HDU containing the kernel in the file given to the @option{--kernel}
+option.
+
@item -E
@itemx --skysubtracted
-If this option is called, it is assumed that the image has already
-been sky subtracted once. Knowing if the sky has already been
-subtracted once or not is very important in estimating the Signal to
-noise ratio of the detections and clumps. In short an extra
address@hidden must be added in the error (noise or
-denominator in the Signal to noise ratio) for every flux value that is
-present in the calculation. This can be interpreted as the error in
-measuring that sky value when it was subtracted by any other
-program. See Section 3.3 in Akhlaghi and Ichikawa (2015) for a
-complete explanation.
+If this option is called, it is assumed that the image has already been sky
+subtracted once. Knowing if the sky has already been subtracted once or not
+is very important in estimating the Signal to noise ratio of the detections
+and clumps. In short an extra @mymath{\sigma_{sky}^2} must be added in the
+error (noise or denominator in the Signal to noise ratio) for every flux
+value that is present in the calculation. This can be interpreted as the
+error in measuring that sky value when it was subtracted by any other
+program. See Section 3.3 in @url{https://arxiv.org/abs/1505.01664, Akhlaghi
+and Ichikawa [2015]}) for a complete explanation.
@item -B FLT
address@hidden --minbfrac=FLT
-Minimum fraction (value between 0 and 1) of blank (undetected) area in a
-mesh for it to be considered in measuring the following properties.
address@hidden --minskyfrac=FLT
+Minimum fraction (value between 0 and 1) of sky (undetected) areas in a
+tile for it to be considered in measuring the following detection and
+segmentation properties.
@itemize
@item
-Measuring the Signal to noise ratio of false detections during the
-false detection removal.
+Measuring the Signal to noise ratio of false detections during the false
+detection removal on small tiles.
@item
-Measuring the sky value (average of undetected pixels). Both before
-the removal of false detections and after it.
+Measuring the sky value (average of undetected pixels) on small tiles. Both
+before the removal of false detections and after it.
@item
-Clump Signal to noise ratio in the blank regions.
+Clump Signal to noise ratio in the sky regions of large files.
@end itemize
-Because of the PSF, astronomical objects, other than cosmic rays,
-never have a clear cutoff and commonly sink into the noise very
-slowly. Even below the very low thresholds used by NoiseChisel. So
-when a large fraction of the area of one mesh is covered by
-detections, it is very probable that their faint wings are present in
-the undetected regions. Therefore, to get an accurate measurement of
-the above parameters over the full mesh grid, meshs that harbor too
-many detected regions should be excluded.
+Because of the PSF and their intrinsic amorphous properties, astronomical
+objects (except cosmic rays) never have a clear cutoff and commonly sink
+into the noise very slowly. Even below the very low thresholds used by
+NoiseChisel. So when a large fraction of the area of one mesh is covered by
+detections, it is very plausible that their faint wings are present in the
+undetected regions (hence causing a bias in any measurement). To get an
+accurate measurement of the above parameters over the tessellation, tiles
+that harbor too many detected regions should be excluded. The used tiles
+are visible in the respective @option{--check} option of the given step.
@item --minnumfalse=INT
-The minimum number of `psudo-detections' (in identifying false detections)
-or clumps (in identifying false clumps) in each large mesh grid. If their
-number is less than this value, this mesh will be left blank and filled
-during mesh interpolation.
-
-The Signal to noise ratio of false detections and clumps in each mesh
-is found using the quantile of the Signal to noise ratio distribution
-of the psudo-detections and clumps over the undetected pixels in each
-mesh. If the number of Signal to noise ratio measurements in each mesh
-is not enough, the quantile will not be accurate. For example if you
-set @option{--detquant=0.99} (or the top 1 percent), then it is best
-to have at least 100 Signal to noise ratio measurements.
-
address@hidden table
-
+The minimum number of `pseudo-detections' (to identify false initial
+detections) or clumps (to identifying false clumps) found over the
+un-detected regions to identify a Signal-to-Noise ratio threshold.
+
+The Signal to noise ratio (S/N) of false pseudo-detections and clumps in
+each tile is found using the quantile of the S/N distribution of the
+psudo-detections and clumps over the undetected pixels in each mesh. If the
+number of S/N measurements is not large enough, the quantile will not be
+accurate (can have large scatter). For example if you set
address@hidden (or the top 1 percent), then it is best to have at
+least 100 S/N measurements.
+
address@hidden -L INT[,INT]
address@hidden --largetilesize=INT[,INT]
+The size of each tile for the tessellation with the larger tile
+sizes. Except for the tile size, all the other parameters for this
+tessellation are taken from the common options described in @ref{Processing
+options}. The format is identical to that of the @option{--tilesize} option
+that is discussed in that section.
+
address@hidden --onlydetection
+If this option is called, no segmentation will be done and the output will
+only have four extensions (no clumps extension, see @ref{NoiseChisel
+output}). The second extension of the output is not going to be objects but
+raw detections (a large region will be given one label): labeling is only
+done based on connectivity. The last two extensions of the output will be
+the Sky and its Standard deviation.
+
+This option can result in faster processing when only the noise properties
+of the image are desired for a catalog using another image's labels for
+example. A common case is when you want to measure colors or SEDs in
+several images. Let's say you have images in two colors: A and B. For
+simplicity also assume that they are exactly on the same position in the
+sky with the same pixel scale.
+
+You choose to set A as a reference, so you run the NoiseChisel fully on
+A. Then you run NoiseChisel on B with @option{--onlydetection} since you
+only need the noise properties of B (for the signal to noise column in its
+catalog). You can then run MakeCatalog on A normally, see
address@hidden To run MakeCatalog on B, you simply set the object and
+clump labels images to those that NoiseChisel produced for A, see
address@hidden astmkcatalog}.
address@hidden
-Detection is the process of separating the pixels in the image into
-two groups: 1) Signal and 2) Noise. Through the parameters below, you
-can customize the detection process in NoiseChisel.
address@hidden @option
address@hidden --grownclumps
+In the output (see @ref{NoiseChisel output}) store the grown clumps (or
+full detected region if only one clump was present in that detection). By
+default the original clumps are stored as the third extension of the
+output, if this option is called, it is replaced with the grown clump
+labels.
@item --continueaftercheck
Continue NoiseChisel after any of the options starting with
@@ -12111,12 +12208,26 @@ Continue NoiseChisel after any of the options
starting with
are many checks, allowing to inspect the status of the processing. The
results of each step affect the next steps of processing, so, when you are
want to check the status of the processing at one step, the time spent to
-complete NoiseChisel is just wasted time. As a result, by default, when you
-use any of the NoiseChisel options that start with @option{--check},
-NoiseChisel will abort once all the desired extensions (steps of the
-processign) to the file have been written into it. With this option, you
-can disable this behavior and ask NoiseChisel to continue with the rest of
-the processing.
+complete NoiseChisel is just wasted/distracting time.
+
+To encourage easier experimentation with the option values, when you use
+any of the NoiseChisel options that start with @option{--check},
+NoiseChisel will abort once all the desired check file(s) is (are)
+completed. If you call the @option{--continueaftercheck} option, you can
+disable this behavior and ask NoiseChisel to continue with the rest of the
+processing after completeting the check file(s).
+
address@hidden table
+
address@hidden Detection options, Segmentation options, General NoiseChisel
options, Invoking astnoisechisel
address@hidden Detection options
+
+Detection is the process of separating the pixels in the image into two
+groups: 1) Signal and 2) Noise. Through the parameters below, you can
+customize the detection process in NoiseChisel. Recall that you can always
+see the full list of Gnuastro's options with the @option{--help} option.
+
address@hidden @option
@item -r FLT
@itemx --mirrordist=FLT
@@ -12134,32 +12245,39 @@ tiles that satisfy this mode and median difference.
@item -t FLT
@itemx --qthresh=FLT
The quantile threshold to apply to the convolved image. The detection
-process begins with applying a quantile threshold to each of the small mesh
-grid elements, see @ref{Tessellation}. The quantile is only calculated for
-those meshs that don't have any significant signal within them, see
address@hidden signal in a tile}.
+process begins with applying a quantile threshold to each of the tiles in
+the small tessellation. The quantile is only calculated for tiles that
+don't have any significant signal within them, see @ref{Quantifying signal
+in a tile}. Interpolation is then used to give a value to the un-successful
+tiles and it is finally smoothed.
@cindex Quantile
@cindex Binary image
@cindex Foreground pixels
@cindex Background pixels
-The quantile value is a floating point value between 0 and 1. Assume
-that we have sorted the @mymath{N} data elements of a distribution
-(the pixels in each mesh on the convolved image). The quantile
-(@mymath{q}) of this distribution is the value of the element with an
-index of (the nearest integer to) @mymath{q\times{N}} in the sorted
-data set. After thresholding is complete, we will have a binary (two
-valued) image. The pixels above the threshold are known as foreground
-pixels (have a value of 1) while those which lie below the threshold
-are known as background (have a value of 0).
+The quantile value is a floating point value between 0 and 1. Assume that
+we have sorted the @mymath{N} data elements of a distribution (the pixels
+in each mesh on the convolved image). The quantile (@mymath{q}) of this
+distribution is the value of the element with an index of (the nearest
+integer to) @mymath{q\times{N}} in the sorted data set. After thresholding
+is complete, we will have a binary (two valued) image. The pixels above the
+threshold are known as foreground pixels (have a value of 1) while those
+which lie below the threshold are known as background (have a value of 0).
+
address@hidden --smoothwidth=INT
+Width of flat kernel used to smooth the interpolated quantile thresholds,
+see @option{--qthresh} for more.
@cindex NaN
@item --checkqthresh
Check the quantile threshold values on the mesh grid. A file suffixed with
address@hidden will be created, see @ref{Tessellation}. With this
-option, NoiseChisel will abort as soon as quantile estimation has been
-completed, allowing you to inspect the status. This behavior can be
-disabled with @option{--continueaftercheck}.
address@hidden will be created showing each step. With this option,
+NoiseChisel will abort as soon as quantile estimation has been completed,
+allowing you to inspect the steps leading to the final quantile threshold,
+this can be disabled with @option{--continueaftercheck}. By default the
+output will have the same pixel size as the input, but with the
address@hidden option, only one pixel will be used for each tile
+(see @ref{Processing options}).
@item -e INT
@itemx --erode=INT
@@ -12237,22 +12355,25 @@ must be an integer.
@item --checkdetsky
Check the initial approximation of the sky value and its standard deviation
-in a FITS file ending with @file{_detsky.fits}. See @ref{Tessellation} for
-more information. If @option{--oneelempertile} is not called (see
address@hidden options}), then the first extension will be the the binary
-image with initial detections labeled one and background labeled zero. The
-mesh values will be in the subsequent extensions.
+in a FITS file ending with @file{_detsky.fits}. With this option,
+NoiseChisel will abort as soon as the sky value used for defining
+pseudo-detections is complete. This allows you to inspect the steps leading
+to the final quantile threshold, this behavior can be disabled with
address@hidden By default the output will have the same
+pixel size as the input, but with the @option{--oneelempertile} option,
+only one pixel will be used for each tile (see @ref{Processing options}).
@item -R FLT
@itemx --dthresh=FLT
The detection threshold: a multiple of the initial sky standard deviation
-added with the initial sky approximation. This flux threshold is applied to
-the initially undetected regions on the unconvolved image. The background
-pixels that are completely engulfed in a 4-connected foreground region are
-converted to background (holes are filled) and one opening (depth of 1) is
-applied over both the initially detected and undetected regions. The Signal
-to noise ratio of the resulting `psudo-detections' are used to identify
-true vs. false detections. See Section 3.1.5 and Figure 7 in Akhlaghi and
+added with the initial sky approximation (which you can inspect with
address@hidden). This flux threshold is applied to the initially
+undetected regions on the unconvolved image. The background pixels that are
+completely engulfed in a 4-connected foreground region are converted to
+background (holes are filled) and one opening (depth of 1) is applied over
+both the initially detected and undetected regions. The Signal to noise
+ratio of the resulting `psudo-detections' are used to identify true
+vs. false detections. See Section 3.1.5 and Figure 7 in Akhlaghi and
Ichikawa (2015) for a very complete explanation.
@item -i INT
@@ -12266,21 +12387,20 @@ psudo-detection that is smaller than this area. Use
@option{--detsnhistnbins} to check if this value is reasonable or not.
@item --checkdetsn
-Save the S/N values of the pseudo-detections into a file ending with
address@hidden You can use these to inspect the values/distribution,
-you can use the histogram creating feature of @ref{Statistics} to make a
-histogram of the distribution (ready for plotting in a text file, or a
-crude ASCII-art demonstration on the command-line).
-
-You can use this for an empirical way to estimate the best
address@hidden value for your dataset. A good minimum area results
-in the histogram having a sharp drop towards the higher S/Ns. In other
-words, when there is a prominent peak in the histogram and the last few
-bins have less than 10 (an arbitrary number, meaning very few!)
-pseudo-detections. When the minimum area is too large or too small, the
-slope of the histogram in the higher S/N bins will not be too sharp with
-the very high S/N bins having a large number of objects and the peak being
-less significant.
+Save the S/N values of the pseudo-detections into two files ending with
address@hidden and @file{_detsn_det.XXX}. The @file{.XXX} is
+determined from the @option{--tableformat} option (see @ref{Input output
+options}, for example @file{.txt} or @file{.fits}). You can use these to
+inspect the S/N values and their distribution (in combination with the
address@hidden option to see where the pseudo-detections are).
+You can use Gnuastro's @ref{Statistics} to make a histogram of the
+distribution (ready for plotting in a text file, or a crude ASCII-art
+demonstration on the command-line).
+
+With this option, NoiseChisel will abort as soon as the two tables are
+created. This allows you to inspect the steps leading to the final quantile
+threshold, this behavior can be disabled with
address@hidden
@item -c FLT
@itemx --detquant=FLT
@@ -12300,114 +12420,93 @@ never have a clear cutoff, so all the 8-pixels
connected to the border
pixels of a detection might harbor data.
@item --checkdetection
-Every step of the detection process will be added as an extension to a
-file with the suffix @file{_det.fits}. Going through each would just
-be a repeat of the explanations above and also of those in Akhlaghi
-and Ichikawa (2015). The extension label should be sufficient to
-recognize which step you are observing. Viewing all the steps can be
-the best guide in choosing the best set of parameters. Note that
-calling this function will significantly slow NoiseChisel. Normally
-the detection steps are done in parallel, but to show you each step
-individually, the parallel processing has to be halted and restarted
-multiple times. Below are some notes that might be useful in
-interpreting certain steps, beyond the paper.
-
address@hidden
address@hidden
-Going from the first ``Labeled'' extension (for the background
-pseudo-detections, which shows the labeled pseudo-detections) to the
-next extension (``For S/N''), with a relatively low @option{--dthresh}
-value, you will notice that a lot of the large area pseudo-detections
-are removed and not used in calculating the S/N threshold. The main
-reason for this is that they overlap with possible detections. You can
-check by going to the next extension and seeing how there are
-detections there. The filled holes have been covering initial
-detections. If you have many such cases, it might be a good sign that
-your @option{--dthresh} value is too low.
address@hidden itemize
+Every step of the detection process will be added as an extension to a file
+with the suffix @file{_det.fits}. Going through each would just be a repeat
+of the explanations above and also of those in Akhlaghi and Ichikawa
+(2015). The extension label should be sufficient to recognize which step
+you are observing. Viewing all the steps can be the best guide in choosing
+the best set of parameters. With this option, NoiseChisel will abort as
+soon as a snapshot of all the detection process is saved. This behavior can
+be disabled with @option{--continueaftercheck}.
@item --checksky
-Check the final sky and its standard deviation values on the mesh
-grid. Similar to @option{--checkdetectionsky}.
-
address@hidden --checkmaskdet
-Mask (set to NaN) the undetected pixels in one extension and the
-detected pixels in the next extension of a file ending with
address@hidden
+Check the derivation of the final sky and its standard deviation values on
+the mesh grid. With this option, NoiseChisel will abort as soon as the sky
+value is estimated over the image (on each tile). This behavior can be
+disabled with @option{--continueaftercheck}. By default the output will
+have the same pixel size as the input, but with the
address@hidden option, only one pixel will be used for each tile
+(see @ref{Processing options}).
@end table
address@hidden
-Segmentation is the process of possibly breaking up a detection into
-multiple objects and clumps. In deep surveys segmentation becomes
-particularly important since galaxies might fall along the same line
-of sight or they might be merging. It is thus very important to be
-able separate the pixels within a detection if it is necessary. After
-segmentation, such a detected region will get different labels.
-
-In NoiseChisel, segmentation is done by first finding the `true'
-clumps over a detection and then expanding those clumps to a certain
-flux limit. True clumps are found in a process very similar to the
-true detections explained above, see Akhlaghi and Ichikawa (2015) for
-more information. If the connections between the grown clumps are
-weaker than a given threshold, the grown clumps are considered to be
-separate objects. Otherwise, the are considered to be part of the same
-object.
address@hidden @option
address@hidden Segmentation options, NoiseChisel output, Detection options,
Invoking astnoisechisel
address@hidden Segmentation options
address@hidden --onlydetect
-If this option is called, no segmentation will be done. The object
-labels extension in the output will simply be the detection (connected
-components) labels and the clumps image will be blank (see
address@hidden output}).
-
-This option can result in faster processing when only the noise
-properties of the image are desired for a catalog using another
-image's labels. A common case is when you want to measure colors or
-SEDs in several images. Let's say you have images in two colors: A and
-B. For simplicity also assume that they are exactly on the same
-position in the sky with the same pixel scale.
-
-You choose to set A as a reference, so you run the full NoiseChisel on
-A. Then you run NoiseChisel on B with this option since you only need
-the noise properties of B (for the signal to noise column in its
-catalog). You can then run MakeCatalog on A normally, see
address@hidden To run MakeCatalog on B, you simply set the object
-and clump labels images to those that NoiseChisel produced for A, see
address@hidden astmkcatalog}.
+Segmentation is the process of (possibly) breaking up a detection into
+multiple segments (technically called @emph{objects} and @emph{clumps} in
+NoiseChisel). In deep surveys segmentation becomes particularly important
+because we will be detecting more diffuse flux so galaxy images are going
+to overlap more. It is thus very important to be able separate the pixels
+within a detection.
+
+In NoiseChisel, segmentation is done by first finding the `true' clumps
+over a detection and then expanding those clumps to a certain flux
+limit. True clumps are found in a process very similar to the true
+detections explained in @ref{Detection options}, see
address@hidden://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 2015} for more
+information. If the connections between the grown clumps are weaker than a
+given threshold, the grown clumps are considered to be separate objects.
+
address@hidden @option
@item -m INT
@itemx --segsnminarea=INT
The minimum area which a clump in the undetected regions should have in
order to be considered in the clump Signal to noise ratio measurement. If
this size is set to a small value, the Signal to noise ratio of false
-clumps will not be accurately found. Note that this has to be larger than
-the value to @option{--detsnminarea}. Because the clumps are found on the
-convolved (smoothed) image while the psudo-detections are found on the
-input image. Use @option{--clumpsnhistnbins} to check if this value is
-reasonable or not.
-
address@hidden --clumpsnhistnbins=INT
-Similar to @option{--detsnhistnbins} (see there for more explanations), but
-for the distribution of the S/N of clumps over the undetected regions. The
-relevant parameter to check here is be @option{--segsnminarea}.
+clumps will not be accurately found. It is recommended that this value be
+larger than the value to @option{--detsnminarea}. Because the clumps are
+found on the convolved (smoothed) image while the psudo-detections are
+found on the input image. You can use @option{--checkclumpsn} and
address@hidden to see if your chosen value is reasonable or
+not.
+
address@hidden --checkclumpsn
+Save the S/N values of the clumps into two files ending with
address@hidden and @file{_clumpsn_det.XXX}. The @file{.XXX} is
+determined from the @option{--tableformat} option (see @ref{Input output
+options}, for example @file{.txt} or @file{.fits}). You can use these to
+inspect the S/N values and their distribution (in combination with the
address@hidden option to see where the clumps are). You can
+use Gnuastro's @ref{Statistics} to make a histogram of the distribution
+(ready for plotting in a text file, or a crude ASCII-art demonstration on
+the command-line).
+
+With this option, NoiseChisel will abort as soon as the two tables are
+created. This allows you to inspect the steps leading to the final S/N
+quantile threshold, this behavior can be disabled with
address@hidden
@item -g FLT
@itemx --segquant=FLT
The quantile of the noise clump Signal to noise ratio distribution. This
-value is used to identify true clumps over the detected regions.
+value is used to identify true clumps over the detected regions. You can
+get the full distribution of clumps S/Ns over the undetected areas with the
address@hidden option and see them with
address@hidden
@item -v
@itemx --keepmaxnearriver
-Keep a clump whose maximum flux is 8-connected to a river pixel. By
-default such clumps over detections are considered to be noise and are
-removed irrespective of their brightness (see @ref{Flux Brightness and
-magnitude}). Over large profiles, that sink into the noise very
-slowly, noise can cause part of the profile (which was flat without
-noise) to become a very large and with a very high Signal to noise
-ratio. In such cases, the pixel with the maximum flux in the clump
-will be immediately touching a river pixel.
+Keep a clump whose maximum flux is 8-connected to a river pixel. By default
+such clumps over detections are considered to be noise and are removed
+irrespective of their brightness (see @ref{Flux Brightness and
+magnitude}). Over large profiles, that sink into the noise very slowly,
+noise can cause part of the profile (which was flat without noise) to
+become a very large and with a very high Signal to noise ratio. In such
+cases, the pixel with the maximum flux in the clump will be immediately
+touching a river pixel.
@item -G FLT
@itemx --gthresh=FLT
@@ -12416,13 +12515,6 @@ stop growing true clumps. Once true clumps are found,
they are set as the
basis to segment the detected region. They are grown until the threshold
specified by this option.
address@hidden --grownclumps
-In the output (see @ref{NoiseChisel output}) store the grown clumps
-(or full detected region if only one clump was present in that
-detection). By default the original clumps are stored as the third
-extension of the output, if this option is called, it is replaced with
-the grown clump labels.
-
@item -y INT
@itemx --minriverlength=INT
The minimum length of a river between two grown clumps for it to be
@@ -12458,18 +12550,22 @@ identified from the noise characteristics of the
image, so you don't
have to give any hand input Signal to noise ratio.
@item --checksegmentation
-A file with the suffix @file{_seg.fits} will be created. This file
-keeps all the relevant steps in finding true clumps and segmenting the
-detections in various extensions. Having read the paper or the steps
-above, the extension name should be enough to understand which step
-each extension is showing. Examining this file can be an excellent guide
-in choosing the best set of parameters. Note that calling this
-function will significantly slow NoiseChisel.
+A file with the suffix @file{_seg.fits} will be created. This file keeps
+all the relevant steps in finding true clumps and segmenting the detections
+into multiple objects in various extensions. Having read the paper or the
+steps above. Examining this file can be an excellent guide in choosing the
+best set of parameters. Note that calling this function will significantly
+slow NoiseChisel. In verbose mode (without the @option{--quiet} option, see
address@hidden mode options}) the important steps (along with their
+extension names) will also be reported.
+
+With this option, NoiseChisel will abort as soon as the two tables are
+created. This behavior can be disabled with @option{--continueaftercheck}.
@end table
address@hidden NoiseChisel output, , NoiseChisel options, Invoking
astnoisechisel
address@hidden NoiseChisel output, , Segmentation options, Invoking
astnoisechisel
@subsubsection NoiseChisel output
The default name and directory of the outputs are explained in
@@ -12497,18 +12593,22 @@ very easy.
@end itemize
@item
-The object labels. Each pixel in the input image is given a label in
-this extension, the labels start from one. The total number of labels
-is stored as the value to the @code{NOBJS} keyword in the header of
-this extension. It is also printed in verbose mode.
+The object/detection labels. Each pixel in the input image is given a label
+in this extension, the labels start from one. If the
address@hidden option is given, each large connected part of the
+image has one label. Without that option, this extension is going to show
+the labels of the objects that are found after segmentation. The total
+number of labels is stored as the value to the @code{NOBJS}/@code{NDETS}
+keyword in the header of this extension. This number is also printed in
+verbose mode.
@item
-The clump labels. All the pixels in the input image that belong to a
-true clump are given a positive label in this extension. The detected
-regions that were not a clump are given a negative value to clearly
-identify the noise from the detections. The total number of clumps in
-this image is stored in the @code{NCLUMPS} keyword of this extension
-and printed in verbose output.
+The clump labels when @option{--onlydetection} is not called. All the
+pixels in the input image that belong to a true clump are given a positive
+label in this extension. The detected regions that were not a clump are
+given a negative value to clearly identify the sky noise from the diffuse
+detections. The total number of clumps in this image is stored in the
address@hidden keyword of this extension and printed in verbose output.
If the @option{--grownclumps} option is called, or a value of @code{1}
is given to it in any of the configuration files, then instead of the
diff --git a/lib/binary.c b/lib/binary.c
index f283450..da6ce4c 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -345,15 +345,15 @@ gal_binary_open(gal_data_t *input, size_t num, int
connectivity,
through the breadth first search algorithm. `binary' has to have an
`uint8' datatype and only zero and non-zero values in it will be
distinguished. The output dataset (which will contain a label on each
- pixel) maybe already allocated (with type `uint32'). If `*out!=NULL',
- the labels will be reset to zero before the start and the labels will be
+ pixel) maybe already allocated (with type `int32'). If `*out!=NULL', the
+ labels will be reset to zero before the start and the labels will be
written into it. If `*out==NULL', the necessary dataset will be
allocated here and put into it. */
size_t
gal_binary_connected_components(gal_data_t *binary, gal_data_t **out,
int connectivity)
{
- uint32_t *l;
+ int32_t *l;
uint8_t *b, *bf;
gal_data_t *lab;
size_t p, i, curlab=1;
@@ -380,10 +380,10 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
error(EXIT_FAILURE, 0, "the `binary' and `out' datasets must have "
"the same size in `gal_binary_connected_components'");
- /* Make sure it has a `uint32' type. */
- if( lab->type!=GAL_TYPE_UINT32 )
+ /* Make sure it has a `int32' type. */
+ if( lab->type!=GAL_TYPE_INT32 )
error(EXIT_FAILURE, 0, "the `out' dataset in "
- "`gal_binary_connected_components' must have `uint32' type"
+ "`gal_binary_connected_components' must have `int32' type"
"but the array you have given is `%s' type",
gal_type_to_string(lab->type, 1));
@@ -391,7 +391,7 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
memset(lab->array, 0, lab->size * gal_type_sizeof(lab->type));
}
else
- lab=*out=gal_data_alloc(NULL, GAL_TYPE_UINT32, binary->ndim,
+ lab=*out=gal_data_alloc(NULL, GAL_TYPE_INT32, binary->ndim,
binary->dsize, binary->wcs, 1,
binary->minmapsize, NULL, "labels", NULL);
@@ -402,7 +402,7 @@ gal_binary_connected_components(gal_data_t *binary,
gal_data_t **out,
l=lab->array;
bf=(b=binary->array)+binary->size;
if( gal_blank_present(binary) )
- do *l++ = *b==GAL_BLANK_UINT8 ? GAL_BLANK_UINT32 : 0; while(++b<bf);
+ do *l++ = *b==GAL_BLANK_UINT8 ? GAL_BLANK_INT32 : 0; while(++b<bf);
/* Go over all the pixels and do a breadth-first: any pixel that is not
@@ -479,11 +479,11 @@ gal_data_t *
gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
size_t *numconnected)
{
- uint32_t *newlabs;
gal_data_t *newlabs_d;
+ int32_t *newlabs, curlab=1;
+ uint8_t *adj=adjacency->array;
struct gal_linkedlist_sll *Q=NULL;
size_t i, j, p, num=adjacency->dsize[0];
- uint8_t *adj=adjacency->array, curlab=1;
/* Some small sanity checks. */
if(adjacency->type != GAL_TYPE_UINT8)
@@ -503,7 +503,7 @@ gal_binary_connected_adjacency_matrix(gal_data_t *adjacency,
/* Allocate (and clear) the output datastructure. */
- newlabs_d=gal_data_alloc(NULL, GAL_TYPE_UINT32, 1, &num, NULL, 1,
+ newlabs_d=gal_data_alloc(NULL, GAL_TYPE_INT32, 1, &num, NULL, 1,
adjacency->minmapsize, NULL, NULL, NULL);
newlabs=newlabs_d->array;
@@ -622,6 +622,10 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t
**outtile)
tile->block=inv;
+ /* Put the input's flags into the inverted array and the tile. */
+ inv->flag = tile->flag = input->flag;
+
+
/* Fill the central regions. */
in=input->array;
GAL_TILE_PARSE_OPERATE({*i = *in==GAL_BLANK_UINT8 ? *in : !*in; ++in;},
@@ -689,7 +693,7 @@ 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'. */
GAL_TILE_PARSE_OPERATE({
- *in = *i>1 && *i!=GAL_BLANK_UINT32 ? 1 : *in;
+ *in = *i>1 && *i!=GAL_BLANK_INT32 ? 1 : *in;
++in;
}, tile, NULL, 0, 0);
diff --git a/lib/blank.c b/lib/blank.c
index ee6351a..05889c6 100644
--- a/lib/blank.c
+++ b/lib/blank.c
@@ -111,7 +111,23 @@ gal_blank_alloc_write(uint8_t type)
-/* Return 1 if the dataset has a blank value and zero if it doesn't. */
+/* Return 1 if the dataset has a blank value and zero if it doesn't. Before
+ checking the dataset, this function will look at its flags. If the
+ `GAL_DATA_FLAG_HASBLANK' or `GAL_DATA_FLAG_DONT_CHECK_ZERO' bits of
+ `input->flag' are set to 1, this function will not do any check and will
+ just use the information in the flags.
+
+ If you want to re-check a dataset which has non-zero flags, then
+ explicitly set the appropriate flag to zero before calling this
+ function. When there are no other flags, you can just set `input->flags'
+ to zero, otherwise you can use this expression:
+
+ input->flags &= ~ (GAL_DATA_FLAG_HASBLANK | GAL_DATA_FLAG_USE_ZERO);
+
+ This function has no side-effects on the dataset: it will not toggle the
+ flags, to avoid repeating parsing of the full dataset multiple times
+ (when it occurs), please toggle the flags your self after the first
+ check. */
#define HAS_BLANK(IT) { \
IT b, *a=input->array, *af=a+input->size, *start; \
gal_blank_write(&b, block->type); \
@@ -123,10 +139,15 @@ gal_blank_alloc_write(uint8_t type)
/* Go over all the elements. */ \
while( start_end_inc[0] + increment <= start_end_inc[1] ) \
{ \
+ /* Necessary when we are on a tile. */ \
if(input!=block) \
af = ( a = start + increment ) + input->dsize[input->ndim-1]; \
+ \
+ /* Check for blank values. */ \
if(b==b) do if(*a==b) return 1; while(++a<af); \
else do if(*a!=*a) return 1; while(++a<af); \
+ \
+ /* Necessary when we are on a tile. */ \
if(input!=block) \
increment += gal_tile_block_increment(block, input->dsize, \
num_increment++, NULL); \
@@ -145,6 +166,12 @@ gal_blank_present(gal_data_t *input)
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 ) )
+ return input->flag & GAL_DATA_FLAG_HASBLANK;
+
/* Go over the pixels and check: */
switch(block->type)
{
@@ -201,60 +228,68 @@ gal_blank_present(gal_data_t *input)
type that has a value of 1 for data that are blank and 0 for those that
aren't. */
#define FLAG_BLANK(IT) { \
- IT b, *a=data->array; \
- gal_blank_write(&b, data->type); \
+ IT b, *a=input->array; \
+ gal_blank_write(&b, input->type); \
if(b==b) /* Blank value can be checked with the equal comparison */ \
do { *o = *a==b; ++a; } while(++o<of); \
else /* Blank value will fail with the equal comparison */ \
do { *o = *a!=*a; ++a; } while(++o<of); \
}
gal_data_t *
-gal_blank_flag(gal_data_t *data)
+gal_blank_flag(gal_data_t *input)
{
uint8_t *o, *of;
gal_data_t *out;
- char **str=data->array, **strf=str+data->size;
+ char **str=input->array, **strf=str+input->size;
- /* Allocate the output array. */
- out=gal_data_alloc(NULL, GAL_TYPE_UINT8, data->ndim, data->dsize,
- data->wcs, 0, data->minmapsize, data->name, data->unit,
- data->comment);
-
- /* Set the pointers for easy looping. */
- of=(o=out->array)+data->size;
-
- /* Go over the pixels and set the output values. */
- switch(data->type)
+ if( gal_blank_present(input) )
{
- /* Numeric types */
- case GAL_TYPE_UINT8: FLAG_BLANK( uint8_t ); break;
- case GAL_TYPE_INT8: FLAG_BLANK( int8_t ); break;
- case GAL_TYPE_UINT16: FLAG_BLANK( uint16_t ); break;
- case GAL_TYPE_INT16: FLAG_BLANK( int16_t ); break;
- case GAL_TYPE_UINT32: FLAG_BLANK( uint32_t ); break;
- case GAL_TYPE_INT32: FLAG_BLANK( int32_t ); break;
- case GAL_TYPE_UINT64: FLAG_BLANK( uint64_t ); break;
- case GAL_TYPE_INT64: FLAG_BLANK( int64_t ); break;
- case GAL_TYPE_FLOAT32: FLAG_BLANK( float ); break;
- case GAL_TYPE_FLOAT64: FLAG_BLANK( double ); break;
-
- /* String. */
- case GAL_TYPE_STRING:
- do *o++ = !strcmp(*str,GAL_BLANK_STRING); while(++str<strf);
- break;
-
- /* Currently unsupported types. */
- case GAL_TYPE_BIT:
- case GAL_TYPE_COMPLEX32:
- case GAL_TYPE_COMPLEX64:
- error(EXIT_FAILURE, 0, "%s type not yet supported in `gal_blank_flag'",
- gal_type_to_string(data->type, 1));
-
- /* Bad input. */
- default:
- error(EXIT_FAILURE, 0, "type value (%d) not recognized "
- "in `gal_blank_flag'", data->type);
+ /* Allocate a non-cleared output array, we are going to parse the
+ input and fill in each element. */
+ out=gal_data_alloc(NULL, GAL_TYPE_UINT8, input->ndim, input->dsize,
+ input->wcs, 0, input->minmapsize, NULL, "bool",
+ NULL);
+
+ /* Set the pointers for easy looping. */
+ of=(o=out->array)+input->size;
+
+ /* Go over the pixels and set the output values. */
+ switch(input->type)
+ {
+ /* Numeric types */
+ case GAL_TYPE_UINT8: FLAG_BLANK( uint8_t ); break;
+ case GAL_TYPE_INT8: FLAG_BLANK( int8_t ); break;
+ case GAL_TYPE_UINT16: FLAG_BLANK( uint16_t ); break;
+ case GAL_TYPE_INT16: FLAG_BLANK( int16_t ); break;
+ case GAL_TYPE_UINT32: FLAG_BLANK( uint32_t ); break;
+ case GAL_TYPE_INT32: FLAG_BLANK( int32_t ); break;
+ case GAL_TYPE_UINT64: FLAG_BLANK( uint64_t ); break;
+ case GAL_TYPE_INT64: FLAG_BLANK( int64_t ); break;
+ case GAL_TYPE_FLOAT32: FLAG_BLANK( float ); break;
+ case GAL_TYPE_FLOAT64: FLAG_BLANK( double ); break;
+
+ /* String. */
+ case GAL_TYPE_STRING:
+ do *o++ = !strcmp(*str,GAL_BLANK_STRING); while(++str<strf);
+ break;
+
+ /* Currently unsupported types. */
+ case GAL_TYPE_BIT:
+ case GAL_TYPE_COMPLEX32:
+ case GAL_TYPE_COMPLEX64:
+ error(EXIT_FAILURE, 0, "%s type not yet supported in "
+ "`gal_blank_flag'", gal_type_to_string(input->type, 1));
+
+ /* Bad input. */
+ default:
+ error(EXIT_FAILURE, 0, "type value (%d) not recognized "
+ "in `gal_blank_flag'", input->type);
+ }
}
+ else
+ /* Allocate a CLEAR data structure (all zeros). */
+ out=gal_data_alloc(NULL, GAL_TYPE_UINT8, input->ndim, input->dsize,
+ input->wcs, 1, input->minmapsize, NULL, "bool", NULL);
/* Return */
return out;
@@ -269,43 +304,50 @@ gal_blank_flag(gal_data_t *data)
the input array, all it does is to shift the blank eleemnts to the end
and adjust the size elements of the `gal_data_t'. */
#define BLANK_REMOVE(IT) { \
- IT b, *a=data->array, *af=a+data->size, *o=data->array; \
- if( gal_blank_present(data) ) \
- { \
- gal_blank_write(&b, data->type); \
- if(b==b) /* Blank value can be be checked with equal. */ \
- do if(*a!=b) { ++num; *o++=*a; } while(++a<af); \
- else /* Blank value will fail on equal comparison. */ \
- do if(*a==*a) { ++num; *o++=*a; } while(++a<af); \
- } \
- else num=data->size; \
+ IT b, *a=input->array, *af=a+input->size, *o=input->array; \
+ gal_blank_write(&b, input->type); \
+ if(b==b) /* Blank value can be be checked with equal. */ \
+ do if(*a!=b) { ++num; *o++=*a; } while(++a<af); \
+ else /* Blank value will fail on equal comparison. */ \
+ do if(*a==*a) { ++num; *o++=*a; } while(++a<af); \
}
void
-gal_blank_remove(gal_data_t *data)
+gal_blank_remove(gal_data_t *input)
{
size_t num=0;
- /* Shift all non-blank elements to the start of the array. */
- switch(data->type)
+ /* This function currently assumes a contiguous patch of memory. */
+ if(input->block && input->ndim!=1 )
+ error(EXIT_FAILURE, 0, "`gal_blank_remove' doesn't currently work on "
+ "tiles in datasets with more dimensions than 1, your input has "
+ "%zu dimensions", input->ndim);
+
+ /* If the dataset doesn't have blank values, then just get the size. */
+ if( gal_blank_present(input) )
{
- case GAL_TYPE_UINT8: BLANK_REMOVE( uint8_t ); break;
- case GAL_TYPE_INT8: BLANK_REMOVE( int8_t ); break;
- case GAL_TYPE_UINT16: BLANK_REMOVE( uint16_t ); break;
- case GAL_TYPE_INT16: BLANK_REMOVE( int16_t ); break;
- case GAL_TYPE_UINT32: BLANK_REMOVE( uint32_t ); break;
- case GAL_TYPE_INT32: BLANK_REMOVE( int32_t ); break;
- case GAL_TYPE_UINT64: BLANK_REMOVE( uint64_t ); break;
- case GAL_TYPE_INT64: BLANK_REMOVE( int64_t ); break;
- case GAL_TYPE_FLOAT32: BLANK_REMOVE( float ); break;
- case GAL_TYPE_FLOAT64: BLANK_REMOVE( double ); break;
- default:
- error(EXIT_FAILURE, 0, "type code %d not recognized in "
- "`gal_blank_remove'", data->type);
+ /* Shift all non-blank elements to the start of the array. */
+ switch(input->type)
+ {
+ case GAL_TYPE_UINT8: BLANK_REMOVE( uint8_t ); break;
+ case GAL_TYPE_INT8: BLANK_REMOVE( int8_t ); break;
+ case GAL_TYPE_UINT16: BLANK_REMOVE( uint16_t ); break;
+ case GAL_TYPE_INT16: BLANK_REMOVE( int16_t ); break;
+ case GAL_TYPE_UINT32: BLANK_REMOVE( uint32_t ); break;
+ case GAL_TYPE_INT32: BLANK_REMOVE( int32_t ); break;
+ case GAL_TYPE_UINT64: BLANK_REMOVE( uint64_t ); break;
+ case GAL_TYPE_INT64: BLANK_REMOVE( int64_t ); break;
+ case GAL_TYPE_FLOAT32: BLANK_REMOVE( float ); break;
+ case GAL_TYPE_FLOAT64: BLANK_REMOVE( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "type code %d not recognized in "
+ "`gal_blank_remove'", input->type);
+ }
}
+ else num=input->size;
/* Adjust the size elements of the dataset. */
- data->ndim=1;
- data->dsize[0]=data->size=num;
+ input->ndim=1;
+ input->dsize[0]=input->size=num;
}
diff --git a/lib/convolve.c b/lib/convolve.c
index 556eefe..699932f 100644
--- a/lib/convolve.c
+++ b/lib/convolve.c
@@ -516,7 +516,6 @@ gal_convolve_spatial_general(gal_data_t *tiles, gal_data_t
*kernel,
size_t numthreads, int edgecorrection,
int convoverch, gal_data_t *tocorrect)
{
- char *name;
struct spatial_params params;
gal_data_t *out, *block=gal_tile_block(tiles);
@@ -535,12 +534,15 @@ gal_convolve_spatial_general(gal_data_t *tiles,
gal_data_t *kernel,
if(tocorrect) out=tocorrect;
else
{
- name = ( block->name
- ? gal_checkset_malloc_cat(block->name, "_CONVOLVED") : NULL );
+ /* Allocate the space for the convolved image. */
out=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, block->ndim, block->dsize,
- block->wcs, 0, block->minmapsize, name,
+ block->wcs, 0, block->minmapsize, NULL,
block->unit, NULL);
- if(name) free(name);
+
+ /* Spatial convolution won't change the blank bit-flag, so use the
+ block structure's blank bit flag. */
+ out->flag = ( block->flag
+ | ( GAL_DATA_FLAG_BLANK_CH | GAL_DATA_FLAG_HASBLANK ) );
}
diff --git a/lib/data.c b/lib/data.c
index f389042..83aed5b 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -64,27 +64,6 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
/*********************************************************************/
/************* Size and allocation *******************/
/*********************************************************************/
-/* Print the contents of any place in memory that `in' points to for `size'
- bytes. This solution was inspired from:
-
-
http://stackoverflow.com/questions/111928/is-there-a-printf-converter-to-print-in-binary-format
*/
-void
-gal_data_bit_print_stream(void *in, size_t size)
-{
- size_t i;
- char *byte=in;
- for(i=0;i<size;++i)
- printf("%c%c%c%c%c%c%c%c ",
- (byte[i] & 0x80 ? '1' : '0'), (byte[i] & 0x40 ? '1' : '0'),
- (byte[i] & 0x20 ? '1' : '0'), (byte[i] & 0x10 ? '1' : '0'),
- (byte[i] & 0x08 ? '1' : '0'), (byte[i] & 0x04 ? '1' : '0'),
- (byte[i] & 0x02 ? '1' : '0'), (byte[i] & 0x01 ? '1' : '0') );
-}
-
-
-
-
-
int
gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second)
{
@@ -325,13 +304,14 @@ gal_data_initialize(gal_data_t *data, void *array,
uint8_t type,
/* Do the simple copying cases. For the display elements, set them all to
impossible (negative) values so if not explicitly set by later steps,
the default values are used if/when printing.*/
- data->status=0;
- data->next=NULL;
- data->ndim=ndim;
- data->type=type;
- data->block=NULL;
- data->mmapname=NULL;
- data->minmapsize=minmapsize;
+ data->flag = 0;
+ data->status = 0;
+ data->next = NULL;
+ data->ndim = ndim;
+ data->type = type;
+ data->block = NULL;
+ data->mmapname = NULL;
+ data->minmapsize = minmapsize;
gal_checkset_allocate_copy(name, &data->name);
gal_checkset_allocate_copy(unit, &data->unit);
gal_checkset_allocate_copy(comment, &data->comment);
@@ -1134,22 +1114,10 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
{
gal_data_t *out;
- /* Allocate space for the output type */
+ /* Allocate the output datastructure. */
out=gal_data_alloc(NULL, newtype, in->ndim, in->dsize, in->wcs,
0, in->minmapsize, in->name, in->unit, in->comment);
- /* Set the rest of the values. */
- out->next=in->next;
- out->status=in->status;
- out->disp_width=in->disp_width;
- out->disp_precision=in->disp_precision;
-
- /* For debugging.
- printf("in: %d (%s)\nout: %d (%s)\n\n", in->type,
- gal_data_type_as_string(in->type, 1), out->type,
- gal_data_type_as_string(out->type, 1));
- */
-
/* Fill in the output array: */
gal_data_copy_to_new_type_to_allocated(in, out, newtype);
@@ -1161,11 +1129,21 @@ gal_data_copy_to_new_type(gal_data_t *in, uint8_t
newtype)
-/* Copy an array into the already allocated space of `out'. Note this
- function won't re-allocate the space if the required size is larger. It
- will abort with an error. If the output's size parameters are larger
- than the input, then this function will update them to be the same as
- the input. */
+/* Copy a given dataset (`in') into an already allocated dataset `out' (the
+ actual dataset and its `array' element). The meta-data of `in' will be
+ fully copied into `out' also. `out->size' will be used to find the
+ available space in the allocated space.
+
+ When `in->size != out->size' this function will behave as follows:
+
+ `out->size < in->size': it won't re-allocate the necessary space, it
+ will abort with an error, so please check before calling this
+ function.
+
+ `out->size > in->size': it will update `out->size' and `out->dsize'
+ to be the same as the input. So if you want to re-use a
+ pre-allocated space with varying input sizes, be sure to reset
+ `out->size' before every call to this function. */
void
gal_data_copy_to_new_type_to_allocated(gal_data_t *in, gal_data_t *out,
uint8_t newtype)
@@ -1185,6 +1163,12 @@ gal_data_copy_to_new_type_to_allocated(gal_data_t *in,
gal_data_t *out,
"dimensions "
"are %zu and %zu respectively", out->ndim, in->ndim);
+ /* Write the constant values. */
+ out->flag=in->flag;
+ out->next=in->next;
+ out->status=in->status;
+ out->disp_width=in->disp_width;
+ out->disp_precision=in->disp_precision;
/* Do the copying. */
switch(out->type)
@@ -1215,7 +1199,6 @@ gal_data_copy_to_new_type_to_allocated(gal_data_t *in,
gal_data_t *out,
"for `out->type' in gal_data_copy_to_new_type", out->type);
}
-
/* Correct the sizes of the output to be the same as the input. If it is
equal, there is no problem, if not, the size information will be
changed, so if you want to use this allocated space again, be sure to
@@ -1228,7 +1211,8 @@ gal_data_copy_to_new_type_to_allocated(gal_data_t *in,
gal_data_t *out,
-/* Copy the input data structure into a new type */
+/* Copy the input data structure into a new type and free the allocated
+ space. */
gal_data_t *
gal_data_copy_to_new_type_free(gal_data_t *in, uint8_t type)
{
diff --git a/lib/gnuastro/data.h b/lib/gnuastro/data.h
index 741f4ef..1383efe 100644
--- a/lib/gnuastro/data.h
+++ b/lib/gnuastro/data.h
@@ -55,6 +55,40 @@ __BEGIN_C_DECLS /* From C++ preparations */
+/* Flag values for the dataset. Note that these are bit-values, so to be
+ more clear, we'll use hexadecimal notation: `0x1' (=1), `0x2' (=2),
+ `0x4' (=4), `0x8' (=8), `0x10' (=16), `0x20' (=32) and so on. */
+
+/* Number of bytes in the unsigned integer hosting the bit-flags (`flag'
+ element) of `gal_data_t'. */
+#define GAL_DATA_FLAG_SIZE 1
+
+/* Bit 0: The has-blank flag has been checked, so a flag value of 0 for the
+ blank flag is strustable. This can be very useful to avoid
+ repetative checks when the necessary value of the bit is 0. */
+#define GAL_DATA_FLAG_BLANK_CH 0x1
+
+/* Bit 1: Dataset contains blank values. */
+#define GAL_DATA_FLAG_HASBLANK 0x2
+
+/* Bit 2: Sorted flags have been checked, see GAL_DATA_FLAG_BLANK_CH. */
+#define GAL_DATA_FLAG_SORT_CH 0x4
+
+/* Bit 3: Dataset is sorted and increasing. */
+#define GAL_DATA_FLAG_SORTED_I 0x8
+
+/* Bit 4: Dataset is sorted and decreasing. */
+#define GAL_DATA_FLAG_SORTED_D 0x10
+
+/* Maximum internal flag value, using bitwise shift operators, this can
+ be used to define internal flags for libraries/programs that depend
+ on Gnuastro without causing any conflict with the internal flags or
+ having to check the values manually on every release. */
+#define GAL_DATA_FLAG_MAXFLAG GAL_DATA_FLAG_SORTED_D
+
+
+
+
/* Main data structure.
mmap (keep data outside of RAM)
@@ -164,7 +198,8 @@ typedef struct gal_data_t
struct wcsprm *wcs; /* WCS information for this dataset. */
/* Content descriptions. */
- int status; /* Any context-specific status value. */
+ uint8_t flag; /* Flags: currently 8-bits are enough. */
+ int status; /* Context specific value for the dataset. */
char *name; /* e.g., EXTNAME, or column, or keyword. */
char *unit; /* Units of the data. */
char *comment; /* A more detailed description of the data. */
@@ -186,9 +221,6 @@ typedef struct gal_data_t
/*********************************************************************/
/************* Size and allocation *******************/
/*********************************************************************/
-void
-gal_data_bit_print_stream(void *in, size_t size);
-
int
gal_data_dsize_is_different(gal_data_t *first, gal_data_t *second);
diff --git a/lib/gnuastro/tile.h b/lib/gnuastro/tile.h
index f198e97..1856c4d 100644
--- a/lib/gnuastro/tile.h
+++ b/lib/gnuastro/tile.h
@@ -85,6 +85,8 @@ gal_tile_block_relative_to_other(gal_data_t *tile, gal_data_t
*other);
+
+
/***********************************************************************/
/************** Tile full dataset ********************/
/***********************************************************************/
@@ -148,6 +150,9 @@ gal_tile_full_id_from_coord(struct
gal_tile_two_layer_params *tl,
size_t *coord);
void
+gal_tile_full_blank_flag(gal_data_t *tile_ll, size_t numthreads);
+
+void
gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl);
diff --git a/lib/gnuastro/type.h b/lib/gnuastro/type.h
index f3b6fa2..4454e77 100644
--- a/lib/gnuastro/type.h
+++ b/lib/gnuastro/type.h
@@ -126,7 +126,8 @@ gal_type_is_linked_list(uint8_t type);
int
gal_type_out(int first_type, int second_type);
-
+char *
+gal_type_bit_string(void *in, size_t size);
__END_C_DECLS /* From C++ preparations */
diff --git a/lib/interpolate.c b/lib/interpolate.c
index 0b3f317..d543a7f 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -282,6 +282,40 @@ interpolate_close_neighbors_on_thread(void *in_prm)
+/* When no interpolation is needed, then we can just copy the input into
+ the output. */
+static gal_data_t *
+interpolate_copy_input(gal_data_t *input, int aslinkedlist)
+{
+ gal_data_t *tin, *tout;
+
+ /* Make a copy of the first input. */
+ tout=gal_data_copy(input);
+ tout->next=NULL;
+
+ /* If we have a linked list, copy each element. */
+ if(aslinkedlist)
+ {
+ /* Note that we have already copied the first input. */
+ for(tin=input->next; tin!=NULL; tin=tin->next)
+ {
+ /* Copy this dataset (will also copy flags). */
+ tout->next=gal_data_copy(tin);
+ tout=tout->next;
+ }
+
+ /* Output is the reverse of the input, so reverse it. */
+ gal_data_reverse_ll(&tout);
+ }
+
+ /* Return the copied list. */
+ return tout;
+}
+
+
+
+
+
/* Interpolate blank values in an array. If the `tl!=NULL', then it is
assumed that the tile values correspond to given tessellation. Such that
`input[i]' corresponds to `tiles[i]' in the tessellation. */
@@ -291,13 +325,20 @@ gal_interpolate_close_neighbors(gal_data_t *input,
size_t numneighbors, size_t numthreads,
int onlyblank, int aslinkedlist)
{
- char *name;
gal_data_t *tin, *tout;
struct interpolate_params prm;
size_t ngbvnum=numthreads*numneighbors;
int permute=(tl && tl->totchannels>1 && tl->workoverch);
+ /* If there are no blank values in the array we should only fill blank
+ values, then simply copy the input and abort. */
+ if( (input->flag | GAL_DATA_FLAG_BLANK_CH) /* Zero bit is meaningful.*/
+ && !(input->flag | GAL_DATA_FLAG_HASBLANK) /* There are no blanks. */
+ && onlyblank ) /* Only interpolate blank.*/
+ return interpolate_copy_input(input, aslinkedlist);
+
+
/* Initialize the constant parameters. */
prm.tl = tl;
prm.ngb_vals = NULL;
@@ -330,11 +371,9 @@ gal_interpolate_close_neighbors(gal_data_t *input,
/* Allocate space for the (first) output. */
- name = input->name ? gal_checkset_malloc_cat("INTRP_", input->name) : NULL;
prm.out=gal_data_alloc(NULL, input->type, input->ndim, input->dsize,
- input->wcs, 0, input->minmapsize, name, input->unit,
- NULL);
- free(name);
+ input->wcs, 0, input->minmapsize, NULL,
+ input->unit, NULL);
gal_linkedlist_add_to_vll(&prm.ngb_vals,
gal_data_malloc_array(input->type, ngbvnum));
@@ -354,12 +393,9 @@ gal_interpolate_close_neighbors(gal_data_t *input,
"`gal_interpolate_close_neighbors'");
/* Allocate the output array for this node. */
- name = ( tin->name
- ? gal_checkset_malloc_cat("INTRP_", tin->name) : NULL );
gal_data_add_to_ll(&prm.out, NULL, tin->type, tin->ndim, tin->dsize,
- tin->wcs, 0, tin->minmapsize, name, tin->unit,
+ tin->wcs, 0, tin->minmapsize, NULL, tin->unit,
NULL);
- if(name) free(name);
/* Allocate the space for the neighbor values of this input. */
gal_linkedlist_add_to_vll(&prm.ngb_vals,
@@ -391,6 +427,15 @@ gal_interpolate_close_neighbors(gal_data_t *input,
}
+ /* The interpolated array doesn't have blank values. So set the blank
+ flag to 0 and set the use-zero to 1. */
+ for(tout=prm.out; tout!=NULL; tout=tout->next)
+ {
+ tout->flag |= GAL_DATA_FLAG_BLANK_CH;
+ tout->flag &= ~GAL_DATA_FLAG_HASBLANK;
+ }
+
+
/* Clean up and return. */
free(prm.thread_flags);
gal_data_free(prm.blanks);
diff --git a/lib/statistics.c b/lib/statistics.c
index ca5a32e..eae030b 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -705,7 +705,7 @@ mode_golden_section(struct statistics_mode_params *p)
/* Once the mode is found, we need to do a quality control. This quality
- control is the measure of its symmetricity. Lets assume the mode index
+ control is the measure of its symmetricity. Let's assume the mode index
is at `m', since an index is just a count, from the Poisson
distribution, the error in `m' is sqrt(m).
@@ -719,7 +719,7 @@ mode_golden_section(struct statistics_mode_params *p)
a completly symmetric mode, this should be 1. Note that the search for
`b` only goes to the 95% of the distribution. */
#define MODE_SYM(IT) { \
- IT *a=p->data->array, af, bf, mf, fi; \
+ IT *a=p->data->array, af=0, bf=0, mf=0, fi; \
\
/* Set the values at the mirror and at `a' (see above). */ \
mf=a[m]; \
@@ -758,9 +758,10 @@ mode_golden_section(struct statistics_mode_params *p)
if(bi==0) bi=topi; \
\
bf = *(IT *)b_val = a[bi]; \
- /* printf("%f, %f, %f\n", af, mf, bf); */ \
+ /*printf("%zu: %f,%f,%f\n", m, (double)af, (double)mf, (double)bf);*/ \
\
- return (bf-mf)/(mf-af); \
+ /* For a bad result, return 0 (which will not output any mode). */ \
+ return bf==af ? 0 : (bf-mf)/(mf-af); \
}
static double
mode_symmetricity(struct statistics_mode_params *p, size_t m, void *b_val)
@@ -847,6 +848,12 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
p.data=gal_statistics_no_blank_sorted(input, inplace);
+ /* It can happen that the whole array is blank. In such cases,
+ `p.data->size==0', so set all output elements to NaN and return. */
+ oa=out->array;
+ if(p.data->size==0) { oa[0]=oa[1]=oa[2]=oa[3]=NAN; return out; }
+
+
/* Basic constants. */
p.tolerance = 0.01;
p.mirrordist = mirrordist;
@@ -876,7 +883,6 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
/* Do the golden-section search iteration, read the mode value from the
input array and save it as the first element of the output dataset's
array, then free the `mode' structure. */
- oa=out->array;
modeindex = mode_golden_section(&p);
gal_data_copy_element_same_type(p.data, modeindex, mode->array);
mode=gal_data_copy_to_new_type_free(mode, GAL_TYPE_FLOAT64);
@@ -1195,16 +1201,14 @@ gal_statistics_no_blank_sorted(gal_data_t *input, int
inplace)
`noblank'.*/
if( gal_blank_present(contig) )
{
- if(inplace) /* If we can free the input, then */
- { /* just remove the blank elements */
- gal_blank_remove(contig); /* in place. */
- noblank=contig;
- }
- else
- {
- noblank=gal_data_copy(contig); /* We aren't allowed to touch */
- gal_blank_remove(noblank); /* the input, so make a copy. */
- }
+ /* See if we should allocate a new dataset to remove blanks or if we
+ can use the actual contiguous patch of memory. */
+ noblank = inplace ? contig : gal_data_copy(contig);
+ gal_blank_remove(noblank);
+
+ /* If we are working in place, then remove the blank flag. */
+ if(inplace)
+ noblank->flag &= ~GAL_DATA_FLAG_HASBLANK;
}
else noblank=contig;
diff --git a/lib/tile.c b/lib/tile.c
index 3ad3c1e..2e2c406 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -358,8 +358,22 @@ gal_tile_block_write_const_value(gal_data_t *tilevalues,
gal_data_t *tilesll,
0, block->minmapsize, tilevalues->name,
tilevalues->unit, tilevalues->comment);
- /* If requested, initialize `tofill'. */
+ /* 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);
+ else
+ {
+ /* Copy the flags. */
+ tofill->flag=tilevalues->flag;
+
+ /* If we have more than one dimension, then remove the possibly
+ sorted flags. */
+ if(block->ndim>1)
+ {
+ tofill->flag &= ~GAL_DATA_FLAG_SORTED_I;
+ tofill->flag &= ~GAL_DATA_FLAG_SORTED_D;
+ }
+ }
/* Go over the tiles and write the values in. Note Recall that `tofill'
has the same type as `tilevalues'. So we are using memcopy. */
@@ -440,6 +454,8 @@ gal_tile_block_relative_to_other(gal_data_t *tile,
gal_data_t *other)
+
+
/***********************************************************************/
/************** Tile full dataset ********************/
/***********************************************************************/
@@ -657,7 +673,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
tiles[i].array=gal_data_ptr_increment(block->array, tind, block->type);
/* Set the sizes of the tile. */
- tiles[i].size=1;
+ tiles[i].size=1; /* Just an initializer, will be changed. */
tiles[i].ndim=input->ndim;
tiles[i].minmapsize=input->minmapsize;
tiles[i].dsize=gal_data_malloc_array(GAL_TYPE_SIZE_T,input->ndim);
@@ -683,6 +699,7 @@ gal_tile_full(gal_data_t *input, size_t *regular,
/* Set the block structure for this tile to the `input', and set the
next pointer as the next tile. Note that only when we are dealing
with the last tile should the `next' pointer be set to NULL.*/
+ tiles[i].flag = 0;
tiles[i].block = input;
tiles[i].next = i==numtiles-1 ? NULL : &tiles[i+1];
@@ -1118,6 +1135,46 @@ gal_tile_full_id_from_coord(struct
gal_tile_two_layer_params *tl,
+/* To use within `gal_tile_full_blank_bit'. */
+static void *
+tile_full_blank_flag(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ gal_data_t *tile_ll=(gal_data_t *)(tprm->params);
+
+ size_t i;
+ gal_data_t *tile;
+
+ /* Check all the tiles given to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+ {
+ tile=&tile_ll[ tprm->indexs[i] ];
+ if( gal_blank_present(tile) )
+ tile->flag |= GAL_DATA_FLAG_HASBLANK;
+ }
+
+ /* Wait for all the other threads to finish. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+/* Update the blank flag on the tiles within the list of input tiles. */
+void
+gal_tile_full_blank_flag(gal_data_t *tile_ll, size_t numthreads)
+{
+ /* Go over all the tiles and update their blank flag. */
+ gal_threads_spin_off(tile_full_blank_flag, tile_ll,
+ gal_data_num_in_ll(tile_ll), numthreads);
+}
+
+
+
+
+
/* Clean up the allocated spaces in the parameters. */
void
gal_tile_full_free_contents(struct gal_tile_two_layer_params *tl)
diff --git a/lib/type.c b/lib/type.c
index 11c39a3..3e4eaf7 100644
--- a/lib/type.c
+++ b/lib/type.c
@@ -30,7 +30,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <string.h>
#include <gnuastro/type.h>
-
+#include <gnuastro/data.h>
@@ -295,3 +295,35 @@ gal_type_out(int first_type, int second_type)
{
return first_type > second_type ? first_type : second_type;
}
+
+
+
+
+
+/* Write the bit (0 or 1) contents of `in' into a string ready for
+ printing. `size' is used to determine the number of bytes to print. The
+ output string will be dynamically allocated within this function. This
+ can be useful for easy checking of bit flag values, for example in an
+ expression like below:
+
+ printf("flag: %s\n", gal_type_bit_string(&flag, sizeof flag) ); */
+char *
+gal_type_bit_string(void *in, size_t size)
+{
+ size_t i;
+ char *byte=in;
+ char *str=gal_data_malloc_array(GAL_TYPE_UINT8, 8*size+1);
+
+ /* Print the bits into the allocated string. This was inspired from
+
+
http://stackoverflow.com/questions/111928/is-there-a-printf-converter-to-print-in-binary-format
*/
+ for(i=0;i<size;++i)
+ sprintf(str+i*8, "%c%c%c%c%c%c%c%c ",
+ (byte[i] & 0x80 ? '1' : '0'), (byte[i] & 0x40 ? '1' : '0'),
+ (byte[i] & 0x20 ? '1' : '0'), (byte[i] & 0x10 ? '1' : '0'),
+ (byte[i] & 0x08 ? '1' : '0'), (byte[i] & 0x04 ? '1' : '0'),
+ (byte[i] & 0x02 ? '1' : '0'), (byte[i] & 0x01 ? '1' : '0') );
+
+ /* Return the allocated and filled string. */
+ return str;
+}
diff --git a/tests/arithmetic/snimage.sh b/tests/arithmetic/snimage.sh
index 361465c..309d254 100755
--- a/tests/arithmetic/snimage.sh
+++ b/tests/arithmetic/snimage.sh
@@ -48,5 +48,5 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img $img - $img / --hdu=0 --hdu=3 --hdu=4 \
+$execname $img $img - $img / --hdu=1 --hdu=4 --hdu=5 \
--output=snimage.fits
diff --git a/tests/arithmetic/where.sh b/tests/arithmetic/where.sh
index 3b82f8f..97a4178 100755
--- a/tests/arithmetic/where.sh
+++ b/tests/arithmetic/where.sh
@@ -48,4 +48,4 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img $img 0 eq nan where -h0 -h1 --output=where.fits
+$execname $img $img 0 eq nan where -h1 -h2 --output=where.fits
diff --git a/tests/during-dev.sh b/tests/during-dev.sh
index eaa2399..4c2a798 100755
--- a/tests/during-dev.sh
+++ b/tests/during-dev.sh
@@ -68,6 +68,7 @@
# symbolic link in `tmpfs-config-make', it is also assumed in the version
# controlled version of this script. Note, if your directory names have
# space characters in them, quote the full value
+numthreads=8
builddir=build
outdir=
@@ -91,9 +92,10 @@ options=
# First, make sure the variables are set. Note that arguments and options
# are not absolutly vital! so they are not checked here. The utility will
# warn and halt if it needs them.
-if [ x"$outdir" = x ]; then echo "outdir is not set."; exit 1; fi
-if [ x$utilname = x ]; then echo "utilname is not set."; exit 1; fi
-if [ x"$builddir" = x ]; then echo "builddir is not set."; exit 1; fi
+if [ x"$outdir" = x ]; then echo "outdir is not set."; exit 1; fi
+if [ x$utilname = x ]; then echo "utilname is not set."; exit 1; fi
+if [ x"$builddir" = x ]; then echo "builddir is not set."; exit 1; fi
+if [ x"$numthreads" = x ]; then echo "numthreads is not set."; exit 1; fi
# If builddir is relative, then append the current directory to make it
@@ -123,7 +125,7 @@ if [ -f "$utility" ]; then rm "$utility"; fi
# Before actually running put a copy of the configuration file in the
# output directory and also add the onlydirconf option so user or system
# wide configuration files don't interfere.
-if make -C "$builddir"; then
+if make -j$numthreads -C "$builddir"; then
# Change to the output directory.
cd "$outdir"
diff --git a/tests/noisechisel/noisechisel.sh b/tests/noisechisel/noisechisel.sh
index a1fb2cc..05d0b05 100755
--- a/tests/noisechisel/noisechisel.sh
+++ b/tests/noisechisel/noisechisel.sh
@@ -48,4 +48,4 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img --checkdetection --checksegmentation
+$execname $img --checkdetection --checksegmentation --continueaftercheck
diff --git a/tmpfs-config-make b/tmpfs-config-make
index 0eacca9..cec5f2a 100755
--- a/tmpfs-config-make
+++ b/tmpfs-config-make
@@ -135,7 +135,7 @@ if [ ! -f Makefile ]; then
--enable-cosmiccal --enable-crop --enable-fits \
--enable-mknoise --enable-mkprof --enable-noisechisel \
--enable-statistics --enable-table --enable-warp \
- CFLAGS="-g -O0" --disable-shared
+ #CFLAGS="-g -O0" --disable-shared
fi
- [gnuastro-commits] master 29d7e46 091/125: Statistics now successfully uses gal_data_t, (continued)
- [gnuastro-commits] master 29d7e46 091/125: Statistics now successfully uses gal_data_t, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master c250a6d 123/125: NoiseChisel's clumps grown and objects defined, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master ffaa059 101/125: Tessellation options now common to all programs, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master dec2b22 078/125: ImageCrop works with gal_data_t features, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master f8b1820 061/125: Other configuration files read immediately when seen, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 46b4651 085/125: Image prefix removed from some program names, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master df716bc 097/125: Tesselation with the new data structure: work started, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master db39cee 088/125: Statistics returns basic results, no array lib, new blank lib, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 60cd76e 090/125: Statistics mostly complete, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 15c524d 124/125: New implementation of NoiseChisel complete,
Mohammad Akhlaghi <=
- [gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 99ac7dc 083/125: Gnuastro now uses fixed width (from stdint.h) types, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master eaf19a3 115/125: NoiseChisel's initial detection step implemented, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 3339b92 125/125: New implementation of MakeCatalog now complete, Mohammad Akhlaghi, 2017/04/23