[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 42e1dda 2/2: gal_label_clump_significance can
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 42e1dda 2/2: gal_label_clump_significance can work on 1D or 2D |
|
Date: |
Thu, 17 May 2018 20:11:20 -0400 (EDT) |
branch: master
commit 42e1ddaa5545b888dbdeba4eb36f529e961106e8
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
gal_label_clump_significance can work on 1D or 2D
Until now, `gal_label_clump_significance' was only defined for 2D datasets,
but since moving to the library, it may be used on 1D datasets too. So it
was modified to work on 1D datasets is well.
---
lib/label.c | 396 ++++++++++++++++++++++++++++++------------------------------
1 file changed, 197 insertions(+), 199 deletions(-)
diff --git a/lib/label.c b/lib/label.c
index af76ca0..0d0528d 100644
--- a/lib/label.c
+++ b/lib/label.c
@@ -507,133 +507,6 @@ gal_label_watershed(gal_data_t *values, gal_data_t
*indexs,
-/* Grow the given labels without creating new ones. */
-void
-gal_label_grow_indexs(gal_data_t *labels, gal_data_t *indexs, int withrivers,
- int connectivity)
-{
- int searchngb;
- size_t *iarray=indexs->array;
- int32_t n1, nlab, *olabel=labels->array;
- size_t *s, *sf, thisround, ninds=indexs->size;
- size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
-
- /* Some basic sanity checks: */
- label_check_type(indexs, GAL_TYPE_SIZE_T, "indexs", __func__);
- label_check_type(labels, GAL_TYPE_INT32, "labels", __func__);
- if(indexs->ndim!=1)
- error(EXIT_FAILURE, 0, "%s: `indexs' has to be a 1D array, but it is "
- "%zuD", __func__, indexs->ndim);
-
- /* The basic idea is this: after growing, not all the blank pixels are
- necessarily filled, for example the pixels might belong to two regions
- above the growth threshold. So the pixels in between them (which are
- below the threshold will not ever be able to get a label). Therefore,
- the safest way we can terminate the loop of growing the objects is to
- stop it when the number of pixels left to fill in this round
- (thisround) equals the number of blanks.
-
- To start the loop, we set `thisround' to one more than the number of
- indexed pixels. Note that it will be corrected immediately after the
- loop has started, it is just important to pass the `while'. */
- thisround=ninds+1;
- while( thisround > ninds )
- {
- /* `thisround' will keep the number of pixels to be inspected in this
- round. `ninds' will count the number of pixels left without a
- label by the end of this round. Since `ninds' comes from the
- previous loop (or outside, for the first round) it has to be saved
- in `thisround' to begin counting a fresh. */
- thisround=ninds;
- ninds=0;
-
- /* Go over all the available indexs. NOTE: while the `indexs->array'
- pointer remains unchanged, `indexs->size' can/will change (get
- smaller) in every loop. */
- sf = (s=indexs->array) + indexs->size;
- do
- {
- /* We'll begin by assuming the nearest neighbor of this pixel
- has no label (has a value of 0). */
- n1=0;
-
- /* Check the neighbors of this pixel. Note that since this
- macro has multiple loops within it, we can't use
- break. We'll use the `searchngb' variable instead. */
- searchngb=1;
- GAL_DIMENSION_NEIGHBOR_OP(*s, labels->ndim, labels->dsize,
- connectivity, dinc,
- {
- if(searchngb)
- {
- /* For easy reading. */
- nlab = olabel[nind];
-
- /* This neighbor's label is meaningful. */
- if(nlab>0) /* This is a real label. */
- {
- if(n1) /* A prev. ngb label has been found. */
- {
- if( n1 != nlab ) /* Different label from */
- { /* prevously found ngb for this pixel. */
- n1=GAL_LABEL_RIVER;
- searchngb=0;
- }
- }
- else
- { /* This is the first labeld neighbor found. */
- n1=nlab;
-
- /* If we want to completely fill in the region
- (`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) searchngb=0;
- }
- }
- }
- } );
-
- /* The loop over neighbors (above) finishes with three
- possibilities:
-
- n1==0 --> No labeled neighbor was found.
- n1==GAL_LABEL_RIVER --> Connecting two labeled regions.
- n1>0 --> Only has one neighbouring label.
-
- 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)
- {
- /* Set the label. */
- olabel[*s]=n1;
-
- /* If this pixel is a river (can only happen when
- `withrivers' is zero), keep it in the loop, because we
- want the `indexs' dataset to contain all non-positive
- (non-labeled) pixels, including rivers. */
- if(n1==GAL_LABEL_RIVER)
- iarray[ ninds++ ] = *s;
- }
- else
- iarray[ ninds++ ] = *s;
-
- /* Correct the size of the `indexs' dataset. */
- indexs->size = indexs->dsize[0] = ninds;
- }
- while(++s<sf);
- }
-
- /* Clean up. */
- free(dinc);
-}
-
-
-
-
-
@@ -650,7 +523,7 @@ gal_label_grow_indexs(gal_data_t *labels, gal_data_t
*indexs, int withrivers,
/**********************************************************************/
-/************* Clump peak intensity *************/
+/************* Clump significance *************/
/**********************************************************************/
static int
label_clump_significance_sanity(gal_data_t *values, gal_data_t *std,
@@ -719,7 +592,7 @@ label_clump_significance_sanity(gal_data_t *values,
gal_data_t *std,
"significance values must have NULL pointers for its `array' "
"and `dsize' pointers (they will be allocated here)", func);
- /* See if the clumps are to be build starting from local maxima or local
+ /* See if the clumps are to be built starting from local maxima or local
minima. */
af=(a=indexs->array)+indexs->size;
do
@@ -730,11 +603,14 @@ label_clump_significance_sanity(gal_data_t *values,
gal_data_t *std,
first=f[*a];
else
{
- /* Note that the elements may have equal values, so for
- `second', we want the first non-blank AND different
- value. */
- if( isnan(second) && f[*a]!=first )
- second=f[*a];
+ if( isnan(second) )
+ {
+ /* Note that the elements may have equal values, so for
+ `second', we want the first non-blank AND different
+ value. */
+ if( f[*a]!=first )
+ second=f[*a];
+ }
else
break;
}
@@ -766,7 +642,7 @@ enum infocols
INFO_X, /* Flux weighted X center col, 0 by C std. */
INFO_Y, /* Flux weighted Y center col. */
INFO_SFF, /* Sum of non-negative pixels (for X,Y). */
- INFO_INSTD, /* Standard deviation at clump center. */
+ INFO_STD, /* Standard deviation. */
INFO_INAREA, /* Tatal area within clump. */
INFO_RIVAREA, /* Tatal area within rivers around clump. */
INFO_PEAK_RIVER, /* Peak (min or max) river value around a clump. */
@@ -778,17 +654,16 @@ static void
label_clump_significance_raw(gal_data_t *values_d, gal_data_t *std_d,
gal_data_t *label_d, gal_data_t *indexs,
struct gal_tile_two_layer_params *tl,
- double *info, size_t numclumps, size_t minarea,
- int variance)
+ double *info)
{
size_t ndim=values_d->ndim, *dsize=values_d->dsize;
double *row;
size_t i, *a, *af, ii, coord[2];
size_t nngb=gal_dimension_num_neighbors(ndim);
+ int32_t nlab, *ngblabs, *label=label_d->array;
float *values=values_d->array, *std=std_d->array;
size_t *dinc=gal_dimension_increment(ndim, dsize);
- int32_t lab, nlab, *ngblabs, *label=label_d->array;
/* Allocate the array to keep the neighbor labels of river pixels. */
ngblabs=gal_pointer_allocate(GAL_TYPE_INT32, nngb, 0, __func__, "ngblabs");
@@ -808,9 +683,11 @@ label_clump_significance_raw(gal_data_t *values_d,
gal_data_t *std_d,
++row[ INFO_INAREA ];
if( values[*a]>0.0f )
{
- row[ INFO_SFF ] += values[*a];
- row[ INFO_X ] += ( values[*a] * (*a/dsize[1]) );
- row[ INFO_Y ] += ( values[*a] * (*a%dsize[1]) );
+ gal_dimension_index_to_coord(*a, ndim, dsize, coord);
+ row[ INFO_SFF ] += values[*a];
+ row[ INFO_X ] += values[*a] * coord[0];
+ if(ndim>1)
+ row[ INFO_Y ] += values[*a] * coord[1];
}
/* In the loop `INFO_INAREA' is just the pixel counter of this
@@ -858,7 +735,21 @@ label_clump_significance_raw(gal_data_t *values_d,
gal_data_t *std_d,
++row[INFO_RIVAREA];
if( row[INFO_RIVAREA]==1.0f )
- row[INFO_PEAK_RIVER] = values[*a];
+ {
+ /* Get the maximum river value. */
+ row[INFO_PEAK_RIVER] = values[*a];
+
+ /* Get the standard deviation. */
+ if(std_d->size==1 || std_d->size==values_d->size)
+ row[INFO_STD]=std_d->size==1?std[0]:std[*a];
+ else
+ {
+ gal_dimension_index_to_coord(*a, ndim, dsize,
+ coord);
+ row[INFO_STD]=
+ std[gal_tile_full_id_from_coord(tl,coord)];
+ }
+ }
}
}
} );
@@ -866,51 +757,6 @@ label_clump_significance_raw(gal_data_t *values_d,
gal_data_t *std_d,
}
while(++a<af);
- /* Based on the position of each clump, find a representative standard
- deviation. */
- for(lab=1; lab<=numclumps; ++lab)
- {
- /* To help in reading. */
- row = &info [ lab * INFO_NCOLS ];
-
- /* The calculations are only necessary for the clumps that satisfy
- the minimum area. There is no need to waste time on the smaller
- ones. */
- if ( row[INFO_INAREA] > minarea )
- {
- /* It might happen that none of the pixels were positive
- (especially over the undetected regions). In that case, set
- the total area of the clump to zero so it is no longer
- considered.*/
- if( row[INFO_SFF]==0.0f )
- row[INFO_INAREA]=0;
- else
- {
- /* Find the coordinates of the clump's weighted center. */
- coord[0]=GAL_DIMENSION_FLT_TO_INT(row[INFO_X]/row[INFO_SFF]);
- coord[1]=GAL_DIMENSION_FLT_TO_INT(row[INFO_Y]/row[INFO_SFF]);
-
- /* Find the corresponding standard deviation. */
- row[INFO_INSTD]=( std_d->size>1
- ? ( std_d->size==values_d->size
- ? std[gal_dimension_coord_to_index(ndim,
- dsize, coord)]
- : std[gal_tile_full_id_from_coord(tl,
- coord)] )
- : std[0] );
- if(variance) row[INFO_INSTD] = sqrt(row[INFO_INSTD]);
-
- /* For a check
- printf("---------\n");
- printf("\t%f --> %zu\n", row[INFO_Y]/row[INFO_SFF], coord[1]);
- printf("\t%f --> %zu\n", row[INFO_X]/row[INFO_SFF], coord[0]);
- printf("%u: (%zu, %zu): %.3f\n", lab, coord[1]+1,
- coord[0]+1, row[INFO_INSTD]);
- */
- }
- }
- }
-
/* Clean up. */
free(dinc);
free(ngblabs);
@@ -931,9 +777,9 @@ gal_label_clump_significance(gal_data_t *values, gal_data_t
*std,
double *info;
int max1_min0;
float *sigarr;
+ double C, R, S, *row;
int32_t *indarr=NULL;
size_t i, ind, counter=0;
- double C, R, S, Ni, *row;
size_t tablen=numclumps+1;
/* If there were no initial clumps, then ignore this function. */
@@ -968,8 +814,7 @@ gal_label_clump_significance(gal_data_t *values, gal_data_t
*std,
/* First, get the raw information necessary for making the S/N table. */
- label_clump_significance_raw(values, std, label, indexs, tl, info,
- numclumps, minarea, variance);
+ label_clump_significance_raw(values, std, label, indexs, tl, info);
/* Calculate the signal to noise ratio for successful clumps */
@@ -982,19 +827,27 @@ gal_label_clump_significance(gal_data_t *values,
gal_data_t *std,
/* If we have a sufficient area and any rivers were actually found
for this clump, then do the measurement. */
- Ni=row[ INFO_INAREA ];
- if( Ni>minarea && row[ INFO_RIVAREA ])
+ if( row[ INFO_INAREA ]>minarea && row[ INFO_RIVAREA ])
{
- /* Calculate the significance ratio, if `keepsmall' is not
+ /* Set the index to write the values. If `keepsmall' is not
called, 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. */
- S = row[ INFO_INSTD ];
- R = row[ INFO_PEAK_RIVER ];
- C = row[ INFO_PEAK_CENTER ];
+ store the signal-to-noise ratios contiguously. Note that
+ counter will always be smaller and equal to i. */
ind = keepsmall ? i : counter++;
+ /* For easy reading. */
+ R = row[ INFO_PEAK_RIVER ];
+ C = row[ INFO_PEAK_CENTER ];
+ S = variance ? sqrt(row[ INFO_STD ]) : row[ INFO_STD ];
+
+ /* NEGATIVE VALUES: Rivers are also defined on the edges of the
+ image and on pixels touching blank pixels. In a strong
+ gradient, such sitations can cause the river to be
+ larger/smaller than the minimum/maximum within the clump. This
+ can only happen in very strong gradients, so for now, I think
+ it is safe to ignore that clump (its negative value will
+ automatically discard it). Later, if we find a problem with
+ this, we'll have to figure out a better solution. */
if(sigind) indarr[ind]=i;
sigarr[ind] = ( max1_min0 ? (C-R) : (R-C) ) / S;
}
@@ -1023,3 +876,148 @@ gal_label_clump_significance(gal_data_t *values,
gal_data_t *std,
/* Clean up. */
free(info);
}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/**********************************************************************/
+/************* Growing labels *************/
+/**********************************************************************/
+/* Grow the given labels without creating new ones. */
+void
+gal_label_grow_indexs(gal_data_t *labels, gal_data_t *indexs, int withrivers,
+ int connectivity)
+{
+ int searchngb;
+ size_t *iarray=indexs->array;
+ int32_t n1, nlab, *olabel=labels->array;
+ size_t *s, *sf, thisround, ninds=indexs->size;
+ size_t *dinc=gal_dimension_increment(labels->ndim, labels->dsize);
+
+ /* Some basic sanity checks: */
+ label_check_type(indexs, GAL_TYPE_SIZE_T, "indexs", __func__);
+ label_check_type(labels, GAL_TYPE_INT32, "labels", __func__);
+ if(indexs->ndim!=1)
+ error(EXIT_FAILURE, 0, "%s: `indexs' has to be a 1D array, but it is "
+ "%zuD", __func__, indexs->ndim);
+
+ /* The basic idea is this: after growing, not all the blank pixels are
+ necessarily filled, for example the pixels might belong to two regions
+ above the growth threshold. So the pixels in between them (which are
+ below the threshold will not ever be able to get a label). Therefore,
+ the safest way we can terminate the loop of growing the objects is to
+ stop it when the number of pixels left to fill in this round
+ (thisround) equals the number of blanks.
+
+ To start the loop, we set `thisround' to one more than the number of
+ indexed pixels. Note that it will be corrected immediately after the
+ loop has started, it is just important to pass the `while'. */
+ thisround=ninds+1;
+ while( thisround > ninds )
+ {
+ /* `thisround' will keep the number of pixels to be inspected in this
+ round. `ninds' will count the number of pixels left without a
+ label by the end of this round. Since `ninds' comes from the
+ previous loop (or outside, for the first round) it has to be saved
+ in `thisround' to begin counting a fresh. */
+ thisround=ninds;
+ ninds=0;
+
+ /* Go over all the available indexs. NOTE: while the `indexs->array'
+ pointer remains unchanged, `indexs->size' can/will change (get
+ smaller) in every loop. */
+ sf = (s=indexs->array) + indexs->size;
+ do
+ {
+ /* We'll begin by assuming the nearest neighbor of this pixel
+ has no label (has a value of 0). */
+ n1=0;
+
+ /* Check the neighbors of this pixel. Note that since this
+ macro has multiple loops within it, we can't use
+ break. We'll use the `searchngb' variable instead. */
+ searchngb=1;
+ GAL_DIMENSION_NEIGHBOR_OP(*s, labels->ndim, labels->dsize,
+ connectivity, dinc,
+ {
+ if(searchngb)
+ {
+ /* For easy reading. */
+ nlab = olabel[nind];
+
+ /* This neighbor's label is meaningful. */
+ if(nlab>0) /* This is a real label. */
+ {
+ if(n1) /* A prev. ngb label has been found. */
+ {
+ if( n1 != nlab ) /* Different label from */
+ { /* prevously found ngb for this pixel. */
+ n1=GAL_LABEL_RIVER;
+ searchngb=0;
+ }
+ }
+ else
+ { /* This is the first labeld neighbor found. */
+ n1=nlab;
+
+ /* If we want to completely fill in the region
+ (`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) searchngb=0;
+ }
+ }
+ }
+ } );
+
+ /* The loop over neighbors (above) finishes with three
+ possibilities:
+
+ n1==0 --> No labeled neighbor was found.
+ n1==GAL_LABEL_RIVER --> Connecting two labeled regions.
+ n1>0 --> Only has one neighbouring label.
+
+ 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)
+ {
+ /* Set the label. */
+ olabel[*s]=n1;
+
+ /* If this pixel is a river (can only happen when
+ `withrivers' is zero), keep it in the loop, because we
+ want the `indexs' dataset to contain all non-positive
+ (non-labeled) pixels, including rivers. */
+ if(n1==GAL_LABEL_RIVER)
+ iarray[ ninds++ ] = *s;
+ }
+ else
+ iarray[ ninds++ ] = *s;
+
+ /* Correct the size of the `indexs' dataset. */
+ indexs->size = indexs->dsize[0] = ninds;
+ }
+ while(++s<sf);
+ }
+
+ /* Clean up. */
+ free(dinc);
+}