[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 60c126f 120/125: New implementation of over-se
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 60c126f 120/125: New implementation of over-segmentation to find clumps |
|
Date: |
Sun, 23 Apr 2017 22:36:52 -0400 (EDT) |
branch: master
commit 60c126f9b79216e803b9fc2d98e02f4a2f679678
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
New implementation of over-segmentation to find clumps
The new implementation of over-segmentation, as part of NoiseChisel, has
been completed with this commit. Finding clumps is an important part of the
segmentation in NoiseChisel and arguably one of the most complicated!
Other changes that have been done along with this commit:
- The default value of `minbfrac' in NoiseChisel is now 0.8 as opposed to
0.5, this will allow much better estimation of noise properties (by
default). It may be slightly too high for a crowded field, but the
users can change it on the command-line (or in a configuration file)
for such datasets.
- Valgrind (with `--leak-check=full') was run and many cases of leaked
(non-freed) memory have been corrected.
- All the variables defined in `GAL_DIMENSION_NEIGHBOR_OP' now start with
a `gdn_' prefix so they don't get confused with those of the caller's
functions.
---
bin/noisechisel/Makefile.am | 8 +-
bin/noisechisel/args.h | 12 +-
bin/noisechisel/astnoisechisel.conf | 2 +-
bin/noisechisel/clumps.c | 825 +++++++++++++++++++++++++++
bin/noisechisel/{segmentation.c => clumps.h} | 16 +-
bin/noisechisel/detection.c | 20 +-
bin/noisechisel/main.h | 6 +-
bin/noisechisel/noisechisel.c | 7 +-
bin/noisechisel/segmentation.c | 49 ++
bin/noisechisel/sky.c | 19 +-
bin/noisechisel/threshold.c | 3 +-
bin/noisechisel/ui.c | 32 +-
bin/noisechisel/ui.h | 2 +-
lib/binary.c | 1 +
lib/data.c | 9 +-
lib/dimension.c | 3 +-
lib/gnuastro/dimension.h | 126 ++--
lib/interpolate.c | 14 +-
lib/statistics.c | 3 +-
lib/threads.c | 3 +-
lib/tile.c | 1 +
21 files changed, 1046 insertions(+), 115 deletions(-)
diff --git a/bin/noisechisel/Makefile.am b/bin/noisechisel/Makefile.am
index 4918b8a..64023c4 100644
--- a/bin/noisechisel/Makefile.am
+++ b/bin/noisechisel/Makefile.am
@@ -30,11 +30,11 @@ bin_PROGRAMS = astnoisechisel
astnoisechisel_LDADD = -lgnuastro
-astnoisechisel_SOURCES = main.c ui.c detection.c noisechisel.c sky.c \
- segmentation.c threshold.c
+astnoisechisel_SOURCES = main.c ui.c clumps.c detection.c noisechisel.c \
+ sky.c segmentation.c threshold.c
-EXTRA_DIST = main.h authors-cite.h args.h ui.h detection.h noisechisel.h \
- segmentation.h sky.h threshold.h
+EXTRA_DIST = main.h authors-cite.h args.h ui.h clumps.h detection.h \
+ noisechisel.h segmentation.h sky.h threshold.h
diff --git a/bin/noisechisel/args.h b/bin/noisechisel/args.h
index bc20412..b46cdf0 100644
--- a/bin/noisechisel/args.h
+++ b/bin/noisechisel/args.h
@@ -77,7 +77,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MINBFRAC,
"FLT",
0,
- "Min. fraction of undetected area in a tile.",
+ "Min. fraction of undetected area in tile.",
GAL_OPTIONS_GROUP_INPUT,
&p->minbfrac,
GAL_TYPE_FLOAT32,
@@ -160,7 +160,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_MIRRORDIST,
"FLT",
0,
- "Max. distance (error multip.) to find mode.",
+ "Max. dist. (error multip.) to find mode.",
ARGS_GROUP_DETECTION,
&p->mirrordist,
GAL_TYPE_FLOAT32,
@@ -290,7 +290,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_SIGMACLIP,
"FLT,FLT",
0,
- "Multiple of sigma and, tolerance or number.",
+ "Sigma multiple and, tolerance or number.",
ARGS_GROUP_DETECTION,
&p->sigmaclip,
GAL_TYPE_STRING,
@@ -430,7 +430,7 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_SEGQUANT,
"FLT",
0,
- "Quantile of clumps in sky to find true S/N.",
+ "S/N Quantile of true sky clumps.",
ARGS_GROUP_SEGMENTATION,
&p->segquant,
GAL_TYPE_FLOAT32,
@@ -495,10 +495,10 @@ struct argp_option program_options[] =
ARGS_OPTION_KEY_OBJBORDERSN,
"FLT",
0,
- "Min. S/N for grown clumps to be one object.",
+ "Min. S/N for grown clumps as one object.",
ARGS_GROUP_SEGMENTATION,
&p->objbordersn,
- GAL_TYPE_SIZE_T,
+ GAL_TYPE_FLOAT32,
GAL_OPTIONS_RANGE_GE_0,
GAL_OPTIONS_MANDATORY,
GAL_OPTIONS_NOT_SET
diff --git a/bin/noisechisel/astnoisechisel.conf
b/bin/noisechisel/astnoisechisel.conf
index 573e3f5..2d5cf8b 100644
--- a/bin/noisechisel/astnoisechisel.conf
+++ b/bin/noisechisel/astnoisechisel.conf
@@ -19,7 +19,7 @@
# Input:
khdu 1
- minbfrac 0.5
+ minbfrac 0.7
minnumfalse 100
# Tessellation
diff --git a/bin/noisechisel/clumps.c b/bin/noisechisel/clumps.c
new file mode 100644
index 0000000..a520ec4
--- /dev/null
+++ b/bin/noisechisel/clumps.c
@@ -0,0 +1,825 @@
+/*********************************************************************
+NoiseChisel - Detect and segment signal in a noisy dataset.
+NoiseChisel is part of GNU Astronomy Utilities (Gnuastro) package.
+
+Original author:
+ Mohammad Akhlaghi <address@hidden>
+Contributing author(s):
+Copyright (C) 2015, Free Software Foundation, Inc.
+
+Gnuastro is free software: you can redistribute it and/or modify it
+under the terms of the GNU General Public License as published by the
+Free Software Foundation, either version 3 of the License, or (at your
+option) any later version.
+
+Gnuastro is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
+**********************************************************************/
+#include <config.h>
+
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+#include <string.h>
+
+#include <gnuastro/fits.h>
+#include <gnuastro/qsort.h>
+#include <gnuastro/blank.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
+#include <gnuastro/linkedlist.h>
+
+#include <gnuastro-internal/checkset.h>
+
+#include "main.h"
+
+#include "ui.h"
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/***************** Basic structures and macros *****************/
+/**********************************************************************/
+/* Parameters for all threads. */
+struct clumps_params
+{
+ int sky0_det1; /* If working on the Sky or Detections. */
+ gal_data_t *sn; /* Array of clump S/N tables. */
+ gal_data_t *snind; /* Array of clump S/N index (for check). */
+ struct noisechiselparams *p; /* Pointer to main NoiseChisel parameters. */
+};
+
+/* Parameters for one thread. */
+struct clumps_thread_params
+{
+ size_t id; /* ID of this detection/tile over tile. */
+ size_t *topinds; /* Indexs of all local maxima. */
+ size_t numclumps; /* The number of clumps found in this run. */
+ gal_data_t *indexs; /* Array containing indexs of this object. */
+ gal_data_t *info; /* Information for all clumps. */
+ gal_data_t *sn; /* Signal-to-noise ratio for these clumps. */
+ gal_data_t *snind; /* Index of S/N for these clumps. */
+ struct clumps_params *clprm; /* Pointer to main structure. */
+};
+
+/* 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). */
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/****************************************************************
+ ***************** Over segmentation ********************
+ ****************************************************************/
+/* Over-segment the region specified by its indexs into peaks and their
+ respective regions (clumps). This is very similar to the immersion
+ method of Vincent & Soille(1991), but here, we will not separate the
+ image into layers, instead, we will work based on the ordered flux
+ values. If a certain pixel (at a certain level) has no neighbors, it is
+ a local maximum and will be assigned a new label. If it has a labeled
+ neighbor, it will take that label and if there is more than one
+ neighboring labeled region that pixel will be a `river` pixel. */
+void
+clumps_oversegment(struct clumps_thread_params *cltprm)
+{
+ struct noisechiselparams *p=cltprm->clprm->p;
+ size_t ndim=p->input->ndim;
+
+ float *arr=p->conv->array;
+ gal_data_t *indexs=cltprm->indexs;
+ 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;
+
+ /*********************************************
+ For checks and debugging:*
+ gal_data_t *crop;
+ size_t extcount=1;
+ uint32_t *cr, *crf;
+ size_t checkdsize[2]={10,10};
+ size_t checkstart[2]={50,145};
+ char *filename="clumpbuild.fits";
+ size_t checkstartind=gal_dimension_coord_to_index(2, dsize, checkstart);
+ gal_data_t *tile=gal_data_alloc(gal_data_ptr_increment(arr, checkstartind,
+ p->conv->type),
+ GAL_TYPE_INVALID, 2, checkdsize,
+ NULL, 0, 0, NULL, NULL, NULL);
+ tile->block=p->conv;
+ gal_checkset_check_remove_file(filename, 0, 0);
+ if(p->cp.numthreads!=1)
+ error(EXIT_FAILURE, 0, "in the debugging mode of `clumps_oversegment' "
+ "only one thread must be used");
+ crop=gal_data_copy(tile);
+ 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);
+ tile->array=gal_tile_block_relative_to_other(tile, p->clabel);
+ tile->block=p->clabel;
+ **********************************************/
+
+
+ /* Sort the given indexs based on their flux (`gal_qsort_index_arr' is
+ defined as static in `gnuastro/qsort.h') */
+ gal_qsort_index_arr=p->conv->array;
+ qsort(indexs->array, indexs->size, sizeof(size_t),
+ gal_qsort_index_float_decreasing);
+
+
+ /* Initialize the region we want to over-segment. */
+ af=(a=indexs->array)+indexs->size; do clabel[*a]=CLUMPS_INIT; while(++a<af);
+
+
+ /* Go over all the given indexs and pull out the clumps. */
+ af=(a=indexs->array)+indexs->size;
+ do
+ /* When regions of a constant flux or masked regions exist, some later
+ indexs (although they have same flux) will be filled before hand. If
+ they are done, there is no need to do them again. */
+ if(clabel[*a]==CLUMPS_INIT)
+ {
+ /* It might happen where one or multiple regions of the pixels
+ under study have the same flux. So two equal valued pixels of
+ two separate (but equal flux) regions will fall immediately
+ after each other in the sorted list of indexs and we have to
+ account for this.
+
+ Therefore, if we see that the next pixel in the index list has
+ the same flux as this one, it does not guarantee that it should
+ be given the same label. Similar to the breadth first search
+ algorithm for finding connected components, we will search all
+ the neighbours and the neighbours of those neighbours that have
+ the same flux of this pixel to see if they touch any label or
+ not and to finally give them all the same label. */
+ if( (a+1)<af && arr[*a]==arr[*(a+1)] )
+ {
+ /* Label of first neighbor found. */
+ n1=0;
+
+ /* A small sanity check. */
+ if(Q!=NULL || cleanup!=NULL)
+ error(EXIT_FAILURE, 0, "a bug! Please contact us at %s so we "
+ "can fix this problem. In `clumps_oversegment', `Q' and "
+ "`cleanup' should be NULL but while checking the equal"
+ "flux regions they aren't", PACKAGE_BUGREPORT);
+
+ /* Add this pixel to a queue. */
+ gal_linkedlist_add_to_sll(&Q, *a);
+ gal_linkedlist_add_to_sll(&cleanup, *a);
+ clabel[*a] = CLUMPS_TMPCHECK;
+
+ /* Find all the pixels that have the same flux and are
+ connected. */
+ while(Q!=NULL)
+ {
+ /* Pop an element from the queue. */
+ gal_linkedlist_pop_from_sll(&Q, &ind);
+
+ /* Look at the neighbors and see if we already have a
+ 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 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(&cleanup, nind);
+ gal_linkedlist_add_to_sll(&Q, 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)
+ {
+ if(n1==0) n1 = nlab;
+ else if(nlab!=n1) n1 = 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;
+ } );
+ }
+
+ /* Set the label that is to be given to this equal flux
+ region. If `n1' was set to any value, then that label should
+ be used for the whole region. Otherwise, this is a new
+ label, see the case for a non-flat region. */
+ if(n1) rlab = n1;
+ else
+ {
+ rlab = curlab++;
+ if( cltprm->topinds) /* This is a local maximum of */
+ cltprm->topinds[rlab]=*a; /* this region, save its index.*/
+ }
+
+ /* Give the same label to the whole connected equal flux
+ region, except those that might have been on the side of
+ the image and were a river pixel. */
+ while(cleanup!=NULL)
+ {
+ gal_linkedlist_pop_from_sll(&cleanup, &ind);
+ /* If it was on the sides of the image, it has been
+ changed to a river pixel. */
+ if( clabel[ ind ]==CLUMPS_TMPCHECK ) clabel[ ind ]=rlab;
+ }
+ }
+
+ /* The flux of this pixel is not the same as the next sorted
+ flux, so simply find the label for this object. */
+ else
+ {
+ /* `n1' is the label of the first labeled neighbor found, so
+ we'll initialize it to zero. */
+ 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
+ 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)
+ {
+ if(nlab==0) n1=CLUMPS_RIVER;
+ else
+ {
+ if(n1)
+ {
+ if(nlab!=n1) n1=CLUMPS_RIVER;
+ }
+ else n1=nlab;
+ }
+ }
+ });
+
+
+ /* 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
+ is either a river pixel (has more than one labeled neighbor
+ and has been set to `CLUMPS_RIVER' before) or all its
+ neighbors have the same label. In both such cases, rlab
+ should be set to n1.*/
+ if(n1) rlab = n1;
+ else
+ {
+ rlab = curlab++;
+ if( cltprm->topinds )
+ cltprm->topinds[ rlab ]=*a;
+ }
+
+ /* Put the found label in the pixel. */
+ clabel[ *a ] = rlab;
+ }
+
+
+ /*********************************************
+ For checks and debugging:
+ if( *a / dsize[1] >= checkstart[0]
+ && *a / dsize[1] < checkstart[0] + checkdsize[0]
+ && *a % dsize[1] >= checkstart[1]
+ && *a % dsize[1] < checkstart[1] + checkdsize[1] )
+ {
+ printf("%zu (%zu: %zu, %zu): %u\n", ++extcount, *a,
+ (*a%dsize[1])-checkstart[1], (*a/dsize[1])-checkstart[0],
+ clabel[*a]);
+ crop=gal_data_copy(tile);
+ crf=(cr=crop->array)+crop->size;
+ do if(*cr==CLUMPS_RIVER) *cr=0; while(++cr<crf);
+ gal_fits_img_write(crop, filename, NULL, PROGRAM_STRING);
+ gal_data_free(crop);
+ }
+ **********************************************/
+ }
+ while(++a<af);
+
+ /* Save the total number of clumps. */
+ cltprm->numclumps=curlab-1;
+
+
+ /* Set all the river pixels to zero, this is only necessary for the
+ clumps over detections, not the sky clumps (they will be converted
+ over all the image. */
+ if(cltprm->clprm->sky0_det1)
+ {
+ af=(a=indexs->array)+indexs->size;
+ do if( clabel[*a]==CLUMPS_RIVER ) clabel[*a]=0; while(++a<af);
+ }
+
+
+ /*********************************************
+ For checks and debugging:
+ tile->array=NULL;
+ gal_data_free(tile);
+ printf("Total number of clumps: %u\n", curlab-1);
+ **********************************************/
+
+ /* Clean up. */
+ free(dinc);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/***************** S/N threshold *****************/
+/**********************************************************************/
+/* In this function we want to find the general information for each clump
+ in an over-segmented labeled array. The signal in each clump is the
+ average signal inside it subtracted by the average signal in the river
+ pixels around it. So this function will go over all the pixels in the
+ object (already found in deblendclumps()) and add them appropriately.
+
+ The output is an array of size numclumps*INFO_NCOLS. as listed
+ below.*/
+enum infocols
+ {
+ INFO_X, /* Flux weighted X center col, 0 by C std. */
+ INFO_Y, /* Flux weighted Y center col. */
+ INFO_NFF, /* Number of non-negative pixels (for X,Y).*/
+ INFO_INFLUX, /* Tatal flux within clump. */
+ INFO_INAREA, /* Tatal area within clump. */
+ INFO_RIVFLUX, /* Tatal flux within rivers around clump. */
+ INFO_RIVAREA, /* Tatal area within rivers around clump. */
+ INFO_INSTD, /* Standard deviation at clump center. */
+
+ INFO_NCOLS, /* Total number of columns. */
+ };
+static void
+clumps_get_raw_info(struct clumps_thread_params *cltprm)
+{
+ struct noisechiselparams *p=cltprm->clprm->p;
+ size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+ size_t i, *a, *af, ii, coord[2];
+ double *row, *info=cltprm->info->array;
+ 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;
+
+ /* Allocate the array to keep the neighbor labels of river pixels. */
+ ngblabs=gal_data_malloc_array(GAL_TYPE_UINT32, nngb);
+
+ /* Go over all the pixels in this region. */
+ af=(a=cltprm->indexs->array)+cltprm->indexs->size;
+ do
+ if( !isnan(arr[ *a ]) )
+ {
+ /* This pixel belongs to a clump. */
+ if( clabel[ *a ] )
+ {
+ lab=clabel[*a];
+ ++info[ lab * INFO_NCOLS + INFO_INAREA ];
+ info[ lab * INFO_NCOLS + INFO_INFLUX ] += arr[*a];
+ if( arr[*a]>0.0f )
+ {
+ info[ lab * INFO_NCOLS + INFO_NFF ] += arr[*a];
+ info[ lab * INFO_NCOLS + INFO_X ] += arr[*a] * (*a/dsize[1]);
+ info[ lab * INFO_NCOLS + INFO_Y ] += arr[*a] * (*a%dsize[1]);
+ }
+ }
+
+ /* This pixel belongs to a river (has a value of zero and isn't
+ blank). */
+ else
+ {
+ /* We are on a river pixel. So the value of this pixel has to
+ be added to any of the clumps in touches. But since it might
+ touch a labeled region more than once, we use `ngblabs' to
+ keep track of which label we have already added its value
+ to. `ii` is the number of different labels this river pixel
+ has already been considered for. `ngblabs' will keep the list
+ labels. */
+ ii=0;
+ memset(ngblabs, 0, nngb*sizeof *ngblabs);
+
+ /* Look into the 8-connected neighbors (recall that a
+ connectivity of `ndim' means all pixels touching it (even on
+ one vertice). */
+ GAL_DIMENSION_NEIGHBOR_OP(*a, ndim, dsize, ndim, dinc, {
+ /* This neighbor's label. */
+ nlab=clabel[ nind ];
+
+ /* We only want those neighbors that are not rivers (>0) or
+ any of the flag values. */
+ if(nlab && nlab<CLUMPS_MAXLAB)
+ {
+ /* Go over all already checked labels and make sure
+ this clump hasn't already been considered. */
+ for(i=0;i<ii;++i) if(ngblabs[i]==nlab) break;
+
+ /* This neighbor clump hasn't been considered yet: */
+ if(i==ii)
+ {
+ ngblabs[ii++] = nlab;
+ ++info[ nlab * INFO_NCOLS + INFO_RIVAREA ];
+ info[ nlab * INFO_NCOLS + INFO_RIVFLUX ] += arr[*a];
+ }
+ }
+ } );
+ }
+ }
+ while(++a<af);
+
+
+ /* Do the final preparations. All the calculations are only necessary for
+ the clumps that satisfy the minimum area. So there is no need to waste
+ time on the smaller ones. */
+ for(lab=1; lab<=cltprm->numclumps; ++lab)
+ {
+ row = &info [ lab * INFO_NCOLS ];
+ if ( row[INFO_INAREA] > p->segsnminarea )
+ {
+ /* Especially over the undetected regions, it might happen that
+ none of the pixels were positive. In that case, set the total
+ area of the clump to zero so it is no longer considered.*/
+ if( row[INFO_NFF]==0.0f ) row[INFO_INAREA]=0;
+ else
+ {
+ coord[0]=GAL_DIMENSION_FLT_TO_INT(row[INFO_X]/row[INFO_NFF]);
+ coord[1]=GAL_DIMENSION_FLT_TO_INT(row[INFO_Y]/row[INFO_NFF]);
+ row[INFO_INSTD] = std[ gal_tile_full_id_from_coord(&p->cp.tl,
+ coord) ];
+ /* For a check
+ printf("---------\n");
+ printf("\t%f --> %zu\n", row[INFO_Y]/row[INFO_NFF], coord[1]);
+ printf("\t%f --> %zu\n", row[INFO_X]/row[INFO_NFF], coord[0]);
+ printf("%u: (%zu, %zu): %.3f\n", lab, coord[1]+1,
+ coord[0]+1, row[INFO_INSTD]);
+ */
+ }
+ }
+ }
+
+ /* Clean up. */
+ free(dinc);
+ free(ngblabs);
+}
+
+
+
+
+/* Make an S/N table for the clumps in a given region. */
+static void
+clumps_make_sn_table(struct clumps_thread_params *cltprm)
+{
+ struct noisechiselparams *p=cltprm->clprm->p;
+ size_t tablen=cltprm->numclumps+1;
+
+ float *snarr;
+ uint32_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};
+
+
+ /* Allocate the arrays to keep the final S/N table (and possibly S/N
+ index) for this object or tile. */
+ cltprm->sn = &cltprm->clprm->sn[ cltprm->id ];
+ cltprm->sn->ndim = 1; /* Depends on `cltprm->sn' */
+ cltprm->sn->type = GAL_TYPE_FLOAT32;
+ cltprm->sn->dsize = gal_data_malloc_array(GAL_TYPE_SIZE_T, 1);
+ cltprm->sn->array = gal_data_malloc_array(cltprm->sn->type, tablen);
+ cltprm->sn->size = cltprm->sn->dsize[0] = tablen; /* After dsize. */
+ if(p->checkclumpsn)
+ {
+ cltprm->snind = &cltprm->clprm->snind [ cltprm->id ];
+ cltprm->snind->ndim = 1; /* Depends on `cltprm->snind' */
+ cltprm->snind->type = GAL_TYPE_UINT32;
+ 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,
+ tablen);
+ }
+ else cltprm->snind=NULL;
+
+
+ /* Allocate the array to keep the raw information of each clump. Note the
+ `+1' in `infodsize', this is because the labels begin with 1 and we
+ want each label to have one row on the same label.*/
+ 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 detections */
+ snarr=cltprm->sn->array;
+ if(cltprm->snind) indarr=cltprm->snind->array;
+ for(i=1;i<tablen;++i)
+ {
+ /* For readability. */
+ row = &( ((double *)(cltprm->info->array))[ i * INFO_NCOLS ] );
+ Ni = row[ INFO_INAREA ];
+ I = row[ INFO_INFLUX ] / row[ INFO_INAREA ];
+ O = row[ INFO_RIVFLUX ] / row[ INFO_RIVAREA ];
+
+ /* If the inner flux is smaller than the outer flux (happens only in
+ noise cases) or the area is smaller than the minimum area to
+ calculate signal-to-noise, then set the S/N of this segment to
+ zero. */
+ if( Ni>p->segsnminarea && I>O ) /* This is O, not 0 (zero). */
+ {
+ /* Here we have done sky subtraction once. However, if the sky
+ was already subtracted (informed by the user), then the
+ varience should be multiplied by 2. */
+ var = ( (p->skysubtracted ? 2.0f : 1.0f)
+ * row[INFO_INSTD] * row[INFO_INSTD] );
+
+ /* Calculate the Signal to noise ratio, if we are on the
+ noise regions, we don't care about the IDs of the clumps
+ anymore, so store the Signal to noise ratios contiguously
+ (for easy sorting and etc). Note that counter will always
+ be smaller and equal to i. */
+ ind = sky0_det1 ? i : counter++;
+ if(cltprm->snind) indarr[ind]=i;
+ snarr[ind]=( sqrt(Ni/p->cpscorr)*(I-O)
+ / sqrt( (I>0?I:-1*I) + (O>0?O:-1*O) + var ) );
+ }
+ else
+ {
+ ind = sky0_det1 ? i : counter++;
+ if(cltprm->snind) indarr[ind]=i;
+ snarr[ind]=NAN;
+ }
+ }
+
+
+ /* If we are in Sky mode, the sizes have to be corrected */
+ if(sky0_det1==0)
+ {
+ cltprm->sn->dsize[0] = cltprm->sn->size = counter;
+ if(cltprm->snind) cltprm->snind->dsize[0] = cltprm->snind->size=counter;
+ }
+
+
+ /* Clean up. */
+ gal_data_free(cltprm->info);
+}
+
+
+
+
+
+static void *
+clumps_find_make_sn_table(void *in_prm)
+{
+ struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
+ struct clumps_params *clprm=(struct clumps_params *)(tprm->params);
+ struct noisechiselparams *p=clprm->p;
+ size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+ void *tarray;
+ double numdet;
+ 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 *scoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+ size_t *icoord=gal_data_malloc_array(GAL_TYPE_SIZE_T, ndim);
+
+
+ /* Initialize the parameters for this thread. */
+ cltprm.clprm = clprm;
+ cltprm.topinds = clprm->sky0_det1 ? NULL : NULL;
+
+
+ /* Go over all the tiles given to this thread. */
+ for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+ {
+ /* IDs. */
+ cltprm.id = tind = tprm->indexs[i];
+ tile = &p->ltl.tiles[tind];
+
+ /* Change the tile's pointers to the binary image (which has 1 for
+ detected pixels and 0 for un-detected regions). */
+ tarray=tile->array;
+ tblock=tile->block;
+ tile->array = gal_tile_block_relative_to_other(tile, p->binary);
+ tile->block = p->binary;
+
+
+ /* Find the number of detected pixels over this tile. Since this is
+ the binary image, this is just the sum of all the pixels. */
+ tmp=gal_statistics_sum(tile);
+ 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 )
+ {
+ /* Add the indexs of all undetected pixels in this tile into an
+ array. */
+ cltprm.indexs=gal_data_alloc(NULL, GAL_TYPE_SIZE_T, 1, &numsky,
+ NULL, 0, p->cp.minmapsize, NULL, NULL,
+ NULL);
+
+
+ /* Change the tile's block to the clump labels dataset (because
+ we'll need to set the labels of the rivers on the edge of the
+ tile here). */
+ tile->array = gal_tile_block_relative_to_other(tile, p->clabel);
+ tile->block = p->clabel;
+
+ /* We need to set all the pixels on the edge of the tile to
+ rivers and not include them in the list of indexs to set
+ clumps. To do that, we need this tile's starting
+ coordinates. */
+ gal_dimension_index_to_coord(gal_data_ptr_dist(p->clabel->array,
+ tile->array,
+ p->clabel->type),
+ ndim, dsize, scoord);
+
+ /* Add the index of every sky element to the array of indexs. */
+ c=0;
+ 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);
+ gal_dimension_index_to_coord(ind, ndim, dsize, icoord);
+
+ /* If the pixel is on the 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);
+ }, tile, NULL, 0, 1);
+
+ /* Correct the number of indexs. */
+ cltprm.indexs->size=cltprm.indexs->dsize[0]=c;
+
+ /* Generate the clumps over this region. */
+ clumps_oversegment(&cltprm);
+
+ /* Correct the river pixels */
+ GAL_TILE_PARSE_OPERATE({if(*i==CLUMPS_RIVER) *i=0;}, tile,
+ NULL, 0, 1);
+
+ /* Make the clump S/N table. */
+ clumps_make_sn_table(&cltprm);
+
+ /* Clean up. */
+ gal_data_free(cltprm.indexs);
+ }
+
+
+ /* Reset the tile's pointers back to what they were. */
+ tile->array=tarray;
+ tile->block=tblock;
+ }
+
+ /* Clean up. */
+ free(scoord);
+ free(icoord);
+
+ /* Wait for the all the threads to finish and return. */
+ if(tprm->b) pthread_barrier_wait(tprm->b);
+ return NULL;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/***************** High level functins *****************/
+/**********************************************************************/
+/* The job of this function is to find the best signal to noise value to
+ use as a threshold to detect real clumps.
+
+ 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
+ table into the sntablearr array (with the size of the number of
+ meshs). If no clumps could be found in a mesh, then
+ sntablearr[i]=NULL. Otherwise, it points to an array of the useful S/N
+ values in that clump. Note that we don't care about the order of S/N
+ values any more! There is also an accompanying array to keep the number
+ of elements in the final S/N array of each mesh: numclumpsarr.
+
+ Using these two arrays, after all the threads are finished, we can
+ concatenate all the S/N values into one array and send it to the main
+ findsnthresh function in thresh.c. */
+void
+clumps_on_undetected_sn(struct noisechiselparams *p)
+{
+ struct clumps_params clprm;
+
+ /* Initialize/allocate the clump parameters structure, */
+ clprm.p=p;
+ clprm.sky0_det1=0;
+ clprm.sn=gal_data_array_calloc(p->ltl.tottiles);
+ clprm.snind = ( p->checkclumpsn
+ ? gal_data_array_calloc(p->ltl.tottiles) : NULL );
+
+
+ /* Spin off the threads to start the work. */
+ gal_threads_spin_off(clumps_find_make_sn_table, &clprm, p->ltl.tottiles,
+ p->cp.numthreads);
+
+
+ /* If the user wanted to see the steps. */
+ if(p->segmentationname)
+ gal_fits_img_write(p->clabel, p->segmentationname, NULL, PROGRAM_STRING);
+
+ /* Clean up. */
+ gal_data_array_free(clprm.sn, p->ltl.tottiles, 1);
+ if(p->checkclumpsn) gal_data_array_free(clprm.snind, p->ltl.tottiles, 1);
+}
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/clumps.h
similarity index 86%
copy from bin/noisechisel/segmentation.c
copy to bin/noisechisel/clumps.h
index fe6e3f2..94092fd 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/clumps.h
@@ -20,18 +20,10 @@ General Public License for more details.
You should have received a copy of the GNU General Public License
along with Gnuastro. If not, see <http://www.gnu.org/licenses/>.
**********************************************************************/
-#include <config.h>
-
-#include <stdio.h>
-#include <errno.h>
-#include <error.h>
-#include <stdlib.h>
-
-#include "main.h"
-
+#ifndef CLUMPS_H
+#define CLUMPS_H
void
-segmentation(struct noisechiselparams *p)
-{
+clumps_on_undetected_sn(struct noisechiselparams *p);
-}
+#endif
diff --git a/bin/noisechisel/detection.c b/bin/noisechisel/detection.c
index a04f7c8..bcf1924 100644
--- a/bin/noisechisel/detection.c
+++ b/bin/noisechisel/detection.c
@@ -83,6 +83,7 @@ detection_initial(struct noisechiselparams *p)
asprintf(&msg, "Eroded %zu time%s (%zu-connectivity).", p->erode,
p->erode>1?"s":"", p->erodengb);
gal_timing_report(&t1, msg, 2);
+ free(msg);
}
if(p->detectionname)
{
@@ -104,6 +105,7 @@ detection_initial(struct noisechiselparams *p)
asprintf(&msg, "Opened (depth: %zu, %s connectivity).",
p->opening, p->openingngb==4 ? "4" : "8");
gal_timing_report(&t1, msg, 2);
+ free(msg);
}
@@ -386,7 +388,7 @@ detection_pseudo_sn(struct noisechiselparams *p, gal_data_t
*worklab,
uint8_t *flag;
size_t tablen=num+1;
gal_data_t *sn, *snind;
- uint32_t *plabend, *indarr;
+ uint32_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;
@@ -491,11 +493,11 @@ detection_pseudo_sn(struct noisechiselparams *p,
gal_data_t *worklab,
}
- /* calculate the signal to noise for successful detections: */
+ /* 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; }
- for(i=1;i<num+1;++i)
+ for(i=1;i<tablen;++i)
{
ave=brightness[i]/area[i];
if( area[i]>p->detsnminarea && ave>0.0f && xy[i*xyncols]>0.0f )
@@ -558,12 +560,18 @@ detection_pseudo_sn(struct noisechiselparams *p,
gal_data_t *worklab,
threshold_write_sn_table(p, sn, snind, ( s0d1 ? p->detsn_d_name
: p->detsn_s_name ), comments);
gal_linkedlist_free_stll(comments, 1);
+
+ /* Abort NoiseChisel if the user asked for it. */
+ if(s0d1 && !p->continueaftercheck)
+ ui_abort_after_check(p, p->detsn_s_name, p->detsn_d_name,
+ "pseudo-detection S/N values in a table");
}
/* Clean up and return. */
free(xy);
free(area);
+ free(coord);
free(brightness);
if(flag) free(flag);
if(snind) gal_data_free(snind);
@@ -632,6 +640,7 @@ detection_pseudo_real(struct noisechiselparams *p)
p->input->dsize, p->input->wcs, 0, p->cp.minmapsize,
NULL, NULL, NULL);
+
/* Over the Sky: find the pseudo-detections and make the S/N table. */
if(!p->cp.quiet) gettimeofday(&t1, NULL);
numpseudo=detection_pseudo_find(p, workbin, worklab, 0);
@@ -646,6 +655,7 @@ detection_pseudo_real(struct noisechiselparams *p)
asprintf(&msg, "Pseudo-det S/N: %.2f (%.3f quant of %zu).",
snthresh, p->detquant, sn->size);
gal_timing_report(&t1, msg, 2);
+ free(msg);
}
gal_data_free(sn);
gal_data_free(quant);
@@ -825,6 +835,7 @@ detection(struct noisechiselparams *p)
asprintf(&msg, "%zu detections after %zu dilation%s",
p->numobjects, p->dilate, p->dilate>1 ? "s." : ".");
gal_timing_report(&t0, msg, 1);
+ free(msg);
}
if(p->detectionname)
{
@@ -846,5 +857,6 @@ detection(struct noisechiselparams *p)
/* If the user wanted to check the threshold and hasn't called
`continueaftercheck', then stop NoiseChisel. */
if(p->detectionname && !p->continueaftercheck)
- ui_abort_after_check(p, p->detectionname, "showing all detection steps");
+ ui_abort_after_check(p, p->detectionname, NULL,
+ "showing all detection steps");
}
diff --git a/bin/noisechisel/main.h b/bin/noisechisel/main.h
index 4c9ea06..30d0cd3 100644
--- a/bin/noisechisel/main.h
+++ b/bin/noisechisel/main.h
@@ -92,13 +92,16 @@ struct noisechiselparams
char *detsn_d_name; /* Detection pseudo-detections S/N name. */
char *detectionname; /* Name of detection steps file. */
char *skyname; /* Name of Sky estimation steps file. */
- char *clumpsnname; /* Name of Clump S/N values file. */
+ char *clumpsn_s_name; /* Sky clump S/N name. */
+ char *clumpsn_d_name; /* Detection clumps S/N name. */
char *segmentationname; /* Name of segmentation steps file. */
+
gal_data_t *input; /* Input image. */
gal_data_t *kernel; /* Input image. */
gal_data_t *conv; /* Convolved input. */
gal_data_t *binary; /* For binary operations. */
gal_data_t *olabel; /* Labels of objects in the detection. */
+ gal_data_t *clabel; /* Labels of clumps in the detection. */
gal_data_t *sky; /* Mean of undetected pixels, per tile. */
gal_data_t *std; /* STD of undetected pixels, per tile. */
time_t rawtime; /* Starting time of the program. */
@@ -109,6 +112,7 @@ struct noisechiselparams
float cpscorr; /* Counts/second correction. */
size_t numobjects; /* Number of objects detected. */
+ size_t numclumps; /* Number of true clumps. */
size_t maxtcontig; /* Maximum contiguous space for a tile. */
size_t maxltcontig; /* Maximum contiguous space for a tile. */
size_t *maxtsize; /* Maximum size of a single small tile. */
diff --git a/bin/noisechisel/noisechisel.c b/bin/noisechisel/noisechisel.c
index be03d1d..85dd4d6 100644
--- a/bin/noisechisel/noisechisel.c
+++ b/bin/noisechisel/noisechisel.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "sky.h"
#include "detection.h"
#include "threshold.h"
+#include "segmentation.h"
@@ -85,8 +86,8 @@ noisechisel_find_sky_subtract(struct noisechiselparams *p)
/* Abort if the user only wanted to see until this point.*/
if(p->skyname && !p->continueaftercheck)
- ui_abort_after_check(p, p->skyname, "derivation of final Sky (and its "
- "STD) value");
+ ui_abort_after_check(p, p->skyname, NULL,
+ "derivation of final Sky (and its STD) value");
/* Subtract the Sky from the Input and Convolved (necessary for
segmentation) images. */
@@ -158,5 +159,5 @@ noisechisel(struct noisechiselparams *p)
noisechisel_convolve_correct_ch_edges(p);
/* Do the segmentation. */
-
+ segmentation(p);
}
diff --git a/bin/noisechisel/segmentation.c b/bin/noisechisel/segmentation.c
index fe6e3f2..3d83fa1 100644
--- a/bin/noisechisel/segmentation.c
+++ b/bin/noisechisel/segmentation.c
@@ -27,11 +27,60 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <error.h>
#include <stdlib.h>
+#include <gnuastro/fits.h>
+#include <gnuastro/blank.h>
+
#include "main.h"
+#include "ui.h"
+#include "clumps.h"
+#include "segmentation.h"
+
void
segmentation(struct noisechiselparams *p)
{
+ float *f;
+ uint32_t *l, *lf;
+
+
+ /* Start off the counter for the number of objects and clumps. The
+ value to these variables will be the label that is given to the
+ next clump or object found. Note that we stored a copy of the
+ initial number of objects in the numobjsinit variable above.*/
+ p->numclumps=1;
+ p->numobjects=1;
+
+
+ /* If a check segmentation image was requested, then put in the
+ inputs. */
+ if(p->segmentationname)
+ {
+ gal_fits_img_write(p->input, p->segmentationname, NULL, PROGRAM_STRING);
+ gal_fits_img_write(p->conv, p->segmentationname, NULL, PROGRAM_STRING);
+ gal_fits_img_write(p->olabel, p->segmentationname, NULL,
+ PROGRAM_STRING);
+ }
+
+
+ /* 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);
+
+
+ /* 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);
+
+
+ /* Find the clumps over the un-detected regions of the input. */
+ clumps_on_undetected_sn(p);
+
+ /* If the user wanted to check the segmentation and hasn't called
+ `continueaftercheck', then stop NoiseChisel. */
+ if(p->segmentationname && !p->continueaftercheck)
+ ui_abort_after_check(p, p->segmentationname, NULL,
+ "showing all segmentation steps");
}
diff --git a/bin/noisechisel/sky.c b/bin/noisechisel/sky.c
index 92dac5e..cfa009f 100644
--- a/bin/noisechisel/sky.c
+++ b/bin/noisechisel/sky.c
@@ -119,6 +119,9 @@ sky_mean_std_undetected(void *in_prm)
memcpy(gal_data_ptr_increment(p->std->array, tind, type),
gal_data_ptr_increment(meanstd->array, 1, type),
gal_type_sizeof(type));
+
+ /* Clean up. */
+ gal_data_free(meanstd);
}
else
{
@@ -175,19 +178,19 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
/* Get the basic information about the standard deviation
distribution. */
tmp=gal_statistics_median(p->std, 0);
- tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
memcpy(&p->medstd, tmp->array, sizeof p->medstd);
- free(tmp);
+ gal_data_free(tmp);
tmp=gal_statistics_minimum(p->std);
- tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
memcpy(&p->minstd, tmp->array, sizeof p->minstd);
- free(tmp);
+ gal_data_free(tmp);
tmp=gal_statistics_maximum(p->std);
- tmp=gal_data_copy_to_new_type(tmp, GAL_TYPE_FLOAT32);
+ tmp=gal_data_copy_to_new_type_free(tmp, GAL_TYPE_FLOAT32);
memcpy(&p->maxstd, tmp->array, sizeof p->maxstd);
- free(tmp);
+ gal_data_free(tmp);
/* In case the image is in electrons or counts per second, the standard
deviation of the noise will become smaller than unity, so we need to
@@ -201,8 +204,8 @@ sky_and_std(struct noisechiselparams *p, char *checkname)
/* If a check was requested, abort NoiseChisel. */
if(checkname && !p->continueaftercheck)
- ui_abort_after_check(p, checkname, "showing derivation of Sky value"
- "and its standard deviation, or STD");
+ ui_abort_after_check(p, checkname, NULL, "showing derivation of Sky "
+ "value and its standard deviation, or STD");
}
diff --git a/bin/noisechisel/threshold.c b/bin/noisechisel/threshold.c
index 9f9fc87..c510bc8 100644
--- a/bin/noisechisel/threshold.c
+++ b/bin/noisechisel/threshold.c
@@ -467,11 +467,12 @@ threshold_quantile_find_apply(struct noisechiselparams *p)
asprintf(&msg, "%.2f & %0.2f quantile thresholds applied.",
p->qthresh, p->noerodequant);
gal_timing_report(&t1, msg, 2);
+ free(msg);
}
/* If the user wanted to check the threshold and hasn't called
`continueaftercheck', then stop NoiseChisel. */
if(p->qthreshname && !p->continueaftercheck)
- ui_abort_after_check(p, p->qthreshname, "quantile threshold check");
+ ui_abort_after_check(p, p->qthreshname, NULL, "quantile threshold check");
}
diff --git a/bin/noisechisel/ui.c b/bin/noisechisel/ui.c
index 4927fdc..74ded96 100644
--- a/bin/noisechisel/ui.c
+++ b/bin/noisechisel/ui.c
@@ -360,13 +360,19 @@ ui_set_output_names(struct noisechiselparams *p)
/* Clump S/N values. */
if(p->checkclumpsn)
- p->clumpsnname=gal_checkset_automatic_output(&p->cp, basename,
- "_clumpsn.txt");
+ {
+ p->clumpsn_s_name=gal_checkset_automatic_output(&p->cp, basename,
+ ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
+ ? "_clumpsn_sky.txt" : "_clumpsn_sky.fits") );
+ p->clumpsn_d_name=gal_checkset_automatic_output(&p->cp, basename,
+ ( p->cp.tableformat==GAL_TABLE_FORMAT_TXT
+ ? "_clumpsn_det.txt" : "_clumpsn_det.fits") );
+ }
/* Segmentation steps. */
if(p->checksegmentation)
p->segmentationname=gal_checkset_automatic_output(&p->cp, basename,
- "_seg.txt");
+ "_seg.fits");
}
@@ -509,7 +515,7 @@ ui_prepare_tiles(struct noisechiselparams *p)
/* If `continueaftercheck' hasn't been called, abort NoiseChisel. */
if(!p->continueaftercheck)
- ui_abort_after_check(p, tl->tilecheckname,
+ ui_abort_after_check(p, tl->tilecheckname, NULL,
"showing all tiles over the image");
/* Free the name. */
@@ -667,21 +673,27 @@ ui_read_check_inputs_setup(int argc, char *argv[],
/**************************************************************/
void
ui_abort_after_check(struct noisechiselparams *p, char *filename,
- char *description)
+ char *file2name, char *description)
{
+ char *name;
+
+ if(file2name) asprintf(&name, "`%s' and `%s'", filename, file2name);
+ else asprintf(&name, "`%s'", filename);
+
/* Let the user know that NoiseChisel is aborting. */
fprintf(stderr,
"------------------------------------------------\n"
"%s aborted for a check\n"
"------------------------------------------------\n"
- "`%s' (%s) has been created.\n\n"
+ "%s (%s) has been created.\n\n"
"If you want %s to continue its processing AND save any "
"requested check(s), please run it again with "
"`--continueaftercheck'.\n"
"------------------------------------------------\n",
- PROGRAM_NAME, filename, description, PROGRAM_NAME);
+ PROGRAM_NAME, name, description, PROGRAM_NAME);
/* Clean up. */
+ free(name);
ui_free_report(p, NULL);
/* Abort. */
@@ -701,19 +713,23 @@ ui_free_report(struct noisechiselparams *p, struct
timeval *t1)
free(p->cp.output);
if(p->skyname) free(p->skyname);
if(p->detskyname) free(p->detskyname);
- if(p->clumpsnname) free(p->clumpsnname);
if(p->qthreshname) free(p->qthreshname);
if(p->detsn_s_name) free(p->detsn_s_name);
if(p->detsn_d_name) free(p->detsn_d_name);
if(p->detectionname) free(p->detectionname);
+ if(p->clumpsn_s_name) free(p->clumpsn_s_name);
+ if(p->clumpsn_d_name) free(p->clumpsn_d_name);
if(p->segmentationname) free(p->segmentationname);
/* Free the allocated datasets. */
+ gal_data_free(p->sky);
+ gal_data_free(p->std);
gal_data_free(p->conv);
gal_data_free(p->input);
gal_data_free(p->kernel);
gal_data_free(p->binary);
gal_data_free(p->olabel);
+ gal_data_free(p->clabel);
/* Clean up the tile structure. */
p->ltl.numchannels=NULL;
diff --git a/bin/noisechisel/ui.h b/bin/noisechisel/ui.h
index 785a1ea..513304a 100644
--- a/bin/noisechisel/ui.h
+++ b/bin/noisechisel/ui.h
@@ -87,7 +87,7 @@ ui_read_check_inputs_setup(int argc, char *argv[],
void
ui_abort_after_check(struct noisechiselparams *p, char *filename,
- char *description);
+ char *file2name, char *description);
void
ui_free_report(struct noisechiselparams *p, struct timeval *t1);
diff --git a/lib/binary.c b/lib/binary.c
index bb72446..a5bd201 100644
--- a/lib/binary.c
+++ b/lib/binary.c
@@ -519,6 +519,7 @@ binary_make_padded_inverse(gal_data_t *input, gal_data_t
**outtile)
/* Clean up and return. */
free(dsize);
+ free(startcoord);
return inv;
}
diff --git a/lib/data.c b/lib/data.c
index 59831bf..1c2bd0d 100644
--- a/lib/data.c
+++ b/lib/data.c
@@ -243,8 +243,8 @@ static void
gal_data_mmap(gal_data_t *data, int clear)
{
int filedes;
+ uint8_t uc=0;
char *filename;
- unsigned char uc=0;
size_t bsize=data->size*gal_type_sizeof(data->type);
/* Check if the .gnuastro folder exists, write the file there. If it
@@ -269,8 +269,11 @@ gal_data_mmap(gal_data_t *data, int clear)
"%zu bytes", filename, bsize);
- /* Write to the newly set file position so the space is allocated. */
- if( write(filedes, &uc, bsize) == -1)
+ /* Write to the newly set file position so the space is allocated. To do
+ this, we are simply writing `uc' (a byte with value 0) into the space
+ we identified by `lseek' (above). This will ensure that this space is
+ set a side for this array and prepare us to use `mmap'. */
+ if( write(filedes, &uc, 1) == -1)
error(EXIT_FAILURE, errno, "%s: unable to write one byte at the "
"%zu-th position", filename, bsize);
diff --git a/lib/dimension.c b/lib/dimension.c
index 5bc4987..03684ea 100644
--- a/lib/dimension.c
+++ b/lib/dimension.c
@@ -60,7 +60,8 @@ gal_dimension_increment(size_t ndim, size_t *dsize)
/* Along the fastest dimension, it is 1. */
out[ndim-1]=1;
- /* For the rest of the dimensions, it is */
+ /* For the rest of the dimensions, it is the multiple of the faster
+ dimension's length and the value for the previous dimension. */
if(ndim>1)
for(i=ndim-2;i>=0;--i)
out[i]=dsize[i+1]*out[i+1];
diff --git a/lib/gnuastro/dimension.h b/lib/gnuastro/dimension.h
index 9712bf2..61d6cd9 100644
--- a/lib/gnuastro/dimension.h
+++ b/lib/gnuastro/dimension.h
@@ -186,52 +186,52 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim);
*/
#define GAL_DIMENSION_NEIGHBOR_OP(index, ndim, dsize, connectivity, \
dinc, operation) { \
- uint32_t bitstr=0; \
- size_t nind, ind=index; \
- uint8_t D, *is_start, *is_end, *is_edge, one=1; \
+ uint32_t gdn_bitstr=0; \
+ size_t nind, gdn_ind=index; \
+ uint8_t gdn_D, *gdn_is_start, *gdn_is_end, *gdn_is_edge, gdn_one=1; \
\
/* Initialize the start/end. */ \
- is_start=(uint8_t *)(&bitstr); \
- is_end=(uint8_t *)(&bitstr)+1; \
+ gdn_is_start=(uint8_t *)(&gdn_bitstr); \
+ gdn_is_end=(uint8_t *)(&gdn_bitstr)+1; \
\
/* Start with the slowest dimension and see if it is on the edge */ \
/* or not, similar to `gal_dimension_index_to_coord'. In the */ \
/* process, also fill the `connectivity==1' neighbors. */ \
- for(D=0;D<ndim;++D) \
+ for(gdn_D=0;gdn_D<ndim;++gdn_D) \
{ \
/* If this dimension is only one element wide, no neighbors. */ \
- if( (dsize)[D] == 1 ) \
+ if( (dsize)[gdn_D] == 1 ) \
{ \
- *is_start |= 1<<D; \
- *is_end |= 1<<D; \
+ *gdn_is_start |= 1<<gdn_D; \
+ *gdn_is_end |= 1<<gdn_D; \
} \
else \
{ \
- if( ind / (dinc)[D] ) \
+ if( gdn_ind / (dinc)[gdn_D] ) \
{ \
/* We are at the end of this dimension. */ \
- if( ind / (dinc)[D] == (dsize)[D]-1 ) \
+ if( gdn_ind / (dinc)[gdn_D] == (dsize)[gdn_D]-1 ) \
{ \
- *is_end |= one<<D; \
- nind = (index) - (dinc)[D]; {operation;}; \
+ *gdn_is_end |= gdn_one<<gdn_D; \
+ nind = (index) - (dinc)[gdn_D]; {operation;}; \
} \
\
/* Middle of the dimension: both +1 and -1 possible. */ \
else \
{ \
- nind = (index) - (dinc)[D]; {operation;}; \
- nind = (index) + (dinc)[D]; {operation;}; \
+ nind = (index) - (dinc)[gdn_D]; {operation;}; \
+ nind = (index) + (dinc)[gdn_D]; {operation;}; \
} \
} \
else \
{ \
- *is_start |= one<<D; \
- nind = (index) + (dinc)[D]; {operation;}; \
+ *gdn_is_start |= gdn_one<<gdn_D; \
+ nind = (index) + (dinc)[gdn_D]; {operation;}; \
} \
} \
\
/* Change `ind' to the remainder of previous dimensions. */ \
- ind %= dinc[D]; \
+ gdn_ind %= dinc[gdn_D]; \
} \
\
/* We now know if the index is on the edge or not. During the */ \
@@ -240,19 +240,21 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim);
if(connectivity>1 && ndim>1) \
{ \
/* Finalize `is_edge' (bit value 1 for respective dim.). */ \
- is_edge=(uint8_t *)(&bitstr)+2; \
- *is_edge = *is_start | *is_end; \
+ gdn_is_edge=(uint8_t *)(&gdn_bitstr)+2; \
+ *gdn_is_edge = *gdn_is_start | *gdn_is_end; \
\
/* Shared between 2D and 3D datasets. */ \
- if(*is_edge) \
+ if(*gdn_is_edge) \
{ /* NOTE: these are bitwise operators, not conditionals. */ \
- if( !( *is_start & ( one | one<<1 ) ) ) \
+ if( !( *gdn_is_start & ( gdn_one | gdn_one<<1 ) ) ) \
{ nind=(index) - (dinc)[0] - (dinc)[1]; {operation;}; } \
- if( !( *is_start & one ) && !( *is_end & one<<1 ) ) \
+ if( !( *gdn_is_start & gdn_one ) \
+ && !( *gdn_is_end & gdn_one<<1 ) ) \
{ nind=(index) - (dinc)[0] + (dinc)[1]; {operation;}; } \
- if( !( *is_end & one ) && !( *is_start & one<<1 ) ) \
+ if( !( *gdn_is_end & gdn_one ) \
+ && !( *gdn_is_start & gdn_one<<1 ) ) \
{ nind=(index) + (dinc)[0] - (dinc)[1]; {operation;}; } \
- if( !( *is_end & ( one | one<<1 ) ) ) \
+ if( !( *gdn_is_end & ( gdn_one | gdn_one<<1 ) ) ) \
{ nind=(index) + (dinc)[0] + (dinc)[1]; {operation;}; } \
} \
else \
@@ -267,67 +269,79 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim);
if(ndim>2) \
{ \
/* Connectivity == 2. */ \
- if(*is_edge) \
- for(D=0;D<2;++D) \
+ if(*gdn_is_edge) \
+ for(gdn_D=0;gdn_D<2;++gdn_D) \
{ \
- if( !( *is_start & ( one<<D | one<<2 ) ) ) \
- { nind=(index) - (dinc)[D] - (dinc)[2]; {operation;}; } \
- if( !( *is_start & one<<D ) && !( *is_end & one<<2 ) ) \
- { nind=(index) - (dinc)[D] + (dinc)[2]; {operation;}; } \
- if( !( *is_end & one<<D ) && !( *is_start & one<<2 ) ) \
- { nind=(index) + (dinc)[D] - (dinc)[2]; {operation;}; } \
- if( !( *is_end & ( one<<D | one<<2 ) ) ) \
- { nind=(index) + (dinc)[D] + (dinc)[2]; {operation;}; } \
+ if( !( *gdn_is_start & ( gdn_one<<gdn_D | gdn_one<<2 ) ) ) \
+ { nind=(index) - (dinc)[gdn_D] - (dinc)[2]; \
+ {operation;}; } \
+ if( !( *gdn_is_start & gdn_one<<gdn_D ) \
+ && !( *gdn_is_end & gdn_one<<2 ) ) \
+ { nind=(index) - (dinc)[gdn_D] + (dinc)[2]; \
+ {operation;}; } \
+ if( !( *gdn_is_end & gdn_one<<gdn_D ) \
+ && !( *gdn_is_start & gdn_one<<2 ) ) \
+ { nind=(index) + (dinc)[gdn_D] - (dinc)[2]; \
+ {operation;}; } \
+ if( !( *gdn_is_end & ( gdn_one<<gdn_D | gdn_one<<2 ) ) ) \
+ { nind=(index) + (dinc)[gdn_D] + (dinc)[2]; \
+ {operation;}; } \
} \
else \
- for(D=0;D<2;++D) \
+ for(gdn_D=0;gdn_D<2;++gdn_D) \
{ \
- nind=(index) - (dinc)[D] - (dinc)[2]; {operation;}; \
- nind=(index) - (dinc)[D] + (dinc)[2]; {operation;}; \
- nind=(index) + (dinc)[D] - (dinc)[2]; {operation;}; \
- nind=(index) + (dinc)[D] + (dinc)[2]; {operation;}; \
+ nind=(index) - (dinc)[gdn_D] - (dinc)[2]; {operation;}; \
+ nind=(index) - (dinc)[gdn_D] + (dinc)[2]; {operation;}; \
+ nind=(index) + (dinc)[gdn_D] - (dinc)[2]; {operation;}; \
+ nind=(index) + (dinc)[gdn_D] + (dinc)[2]; {operation;}; \
} \
\
/* Connectivity == 3. */ \
if(connectivity>2) \
{ \
- if(*is_edge) \
+ if(*gdn_is_edge) \
{ \
- if( !*is_start ) \
+ if( !*gdn_is_start ) \
{ nind=(index) - (dinc)[0] - (dinc)[1] - (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_start & one) && !(*is_start & one<<1) \
- && !(*is_end & one<<2)) \
+ if( !(*gdn_is_start & gdn_one) \
+ && !(*gdn_is_start & gdn_one<<1) \
+ && !(*gdn_is_end & gdn_one<<2)) \
{ nind=(index) - (dinc)[0] - (dinc)[1] + (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_start & one) && !(*is_end & one<<1) \
- && !(*is_start & one<<2)) \
+ if( !(*gdn_is_start & gdn_one) \
+ && !(*gdn_is_end & gdn_one<<1) \
+ && !(*gdn_is_start & gdn_one<<2)) \
{ nind=(index) - (dinc)[0] + (dinc)[1] - (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_start & one) && !(*is_end & one<<1) \
- && !(*is_end & one<<2)) \
+ if( !(*gdn_is_start & gdn_one) \
+ && !(*gdn_is_end & gdn_one<<1) \
+ && !(*gdn_is_end & gdn_one<<2)) \
{ nind=(index) - (dinc)[0] + (dinc)[1] + (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_end & one) && !(*is_start & one<<1) \
- && !(*is_start & one<<2)) \
+ if( !(*gdn_is_end & gdn_one) \
+ && !(*gdn_is_start & gdn_one<<1) \
+ && !(*gdn_is_start & gdn_one<<2)) \
{ nind=(index) + (dinc)[0] - (dinc)[1] - (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_end & one) && !(*is_start & one<<1) \
- && !(*is_end & one<<2)) \
+ if( !(*gdn_is_end & gdn_one) \
+ && !(*gdn_is_start & gdn_one<<1) \
+ && !(*gdn_is_end & gdn_one<<2)) \
{ nind=(index) + (dinc)[0] - (dinc)[1] + (dinc)[2]; \
{operation;}; } \
\
- if( !(*is_end & one) && !(*is_end & one<<1) \
- && !(*is_start & one<<2)) \
+ if( !(*gdn_is_end & gdn_one) \
+ && !(*gdn_is_end & gdn_one<<1) \
+ && !(*gdn_is_start & gdn_one<<2)) \
{ nind=(index) + (dinc)[0] + (dinc)[1] - (dinc)[2]; \
{operation;}; } \
\
- if( !*is_end ) \
+ if( !*gdn_is_end ) \
{ nind=(index) + (dinc)[0] + (dinc)[1] + (dinc)[2]; \
{operation;}; } \
} \
@@ -356,7 +370,7 @@ gal_dimension_dist_manhattan(size_t *a, size_t *b, size_t
ndim);
\
/* For a check. */ \
/* printf("\nEdge bit flags: "); */ \
- /* gal_data_bit_print_stream(&bitstr, 3); printf("\n"); */ \
+ /* gal_data_bit_print_stream(&gdn_bitstr, 3); printf("\n"); */ \
}
diff --git a/lib/interpolate.c b/lib/interpolate.c
index eba3fbd..0b3f317 100644
--- a/lib/interpolate.c
+++ b/lib/interpolate.c
@@ -206,8 +206,13 @@ interpolate_close_neighbors_on_thread(void *in_prm)
tin=tin->next;
}
- /* If we have filled all the elements, break out. */
- if(++ngb_counter>=prm->numneighbors) break;
+ /* If we have filled all the elements clean up the linked
+ list and break out. */
+ if(++ngb_counter>=prm->numneighbors)
+ {
+ if(lQ) gal_linkedlist_tosll_free(lQ);
+ break;
+ }
}
/* Go over all the neighbors of this popped pixel and add them to
@@ -265,6 +270,7 @@ interpolate_close_neighbors_on_thread(void *in_prm)
gal_data_free_ll(nearest);
free(icoord);
free(ncoord);
+ free(dinc);
/* Wait for all the other threads to finish and return. */
@@ -386,8 +392,8 @@ gal_interpolate_close_neighbors(gal_data_t *input,
/* Clean up and return. */
- gal_data_free(prm.blanks);
free(prm.thread_flags);
- free(prm.ngb_vals);
+ gal_data_free(prm.blanks);
+ gal_linkedlist_free_vll(prm.ngb_vals, 1);
return prm.out;
}
diff --git a/lib/statistics.c b/lib/statistics.c
index 7476e9c..ca5a32e 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -892,7 +892,7 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
output values to NaN. */
if(oa[2]>GAL_STATISTICS_MODE_GOOD_SYM)
{
- b_val=gal_data_copy_to_new_type_free(mode, GAL_TYPE_FLOAT64);
+ b_val=gal_data_copy_to_new_type_free(b_val, GAL_TYPE_FLOAT64);
gal_data_copy_element_same_type(b_val, 0, &oa[3]);
}
else oa[0]=oa[1]=oa[2]=oa[3]=NAN;
@@ -905,6 +905,7 @@ gal_statistics_mode(gal_data_t *input, float mirrordist,
int inplace)
/* Clean up (if necessary), then return the output */
if(p.data!=input) gal_data_free(p.data);
+ gal_data_free(b_val);
gal_data_free(mode);
return out;
}
diff --git a/lib/threads.c b/lib/threads.c
index c1bfa46..faf78e5 100644
--- a/lib/threads.c
+++ b/lib/threads.c
@@ -281,10 +281,11 @@ gal_threads_attr_barrier_init(pthread_attr_t *attr,
pthread_barrier_t *b,
void *
run_on_thread(void *in_prm)
{
- size_t i;
struct gal_threads_params *tprm=(struct gal_threads_params *)in_prm;
struct my_params *prm=(struct my_params *)(tprm->params);
+ size_t i;
+
for(i=0; tprm->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
{
diff --git a/lib/tile.c b/lib/tile.c
index 1c2246a..3ad3c1e 100644
--- a/lib/tile.c
+++ b/lib/tile.c
@@ -825,6 +825,7 @@ gal_tile_full_two_layers(gal_data_t *input,
"`gal_tile_full_two_layers' don't match in dimension %zu, with "
"values of %zu and %zu respectively.", ndim-i,
tl->numchannels[i], junk[i]);
+ free(junk);
free(junk2);
/* Tile each channel. While tiling the first channel, we are also going
- [gnuastro-commits] master 1e961f5 067/125: Starting to use gal_data_t in MakeProfiles, (continued)
- [gnuastro-commits] master 1e961f5 067/125: Starting to use gal_data_t in MakeProfiles, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e58ed29 108/125: One value per tile in multi-channel tessellation, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0e344fe 086/125: Warp now uses gal_data_t and all its fruits, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 8348f1e 119/125: NoiseChisel pseudo-detections identified on threads, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master fca5e80 117/125: First re-implementation of NoiseChisel's complete detection, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master a1ff97b 111/125: Permutations and neighbor finding algorithms, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 0bc96c9 122/125: Clumps and their S/N tables over detections, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master 80cdd4f 118/125: Correction of convolution over channel borders, Mohammad Akhlaghi, 2017/04/23
- [gnuastro-commits] master e90b490 116/125: Internal headers moved to special directory, Mohammad Akhlaghi, 2017/04/23
- [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 <=
- [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, 2017/04/23
- [gnuastro-commits] master 3b26d5c 114/125: SubtractSky program removed, Statistics estimates Sky, Mohammad Akhlaghi, 2017/04/23