[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 8161798: MakeCatalog: positions of min and max
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 8161798: MakeCatalog: positions of min and max value account for many pixels |
|
Date: |
Mon, 26 Oct 2020 15:27:18 -0400 (EDT) |
branch: master
commit 816179883ef6d726b88cd3cc0177b46614c1133e
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
MakeCatalog: positions of min and max value account for many pixels
Until now, in the recently added columns like '--maxvx' and '--maxvy' we
were mistakely reporting the postiion of the first pixel that has the
maximum value. But in scenarios where the centers of saturated stars may be
saturated (and they are replaced by the maximum value), this would cause
significantly different positions (when there is a saturation spike in
particular).
With this commit, the mean position of all the maximum pixels is returned,
not just the first one. This fixes the problem above.
Also, I added to new columns: '--areaminv' and '--areamaxv' to report the
number of pixels that have the minimum/maximum value. This can be useful in
identifying saturated stars in the scenarios above.
---
NEWS | 2 +
bin/arithmetic/arithmetic.c | 11 ++-
bin/mkcatalog/args.h | 28 ++++++++
bin/mkcatalog/columns.c | 139 +++++++++++++++++++++++++++++--------
bin/mkcatalog/main.h | 4 ++
bin/mkcatalog/parse.c | 164 +++++++++++++++++++++++++-------------------
bin/mkcatalog/ui.c | 4 ++
bin/mkcatalog/ui.h | 2 +
doc/gnuastro.texi | 6 ++
9 files changed, 261 insertions(+), 99 deletions(-)
diff --git a/NEWS b/NEWS
index a29ad64..60a7367 100644
--- a/NEWS
+++ b/NEWS
@@ -70,6 +70,8 @@ See the end of the file for license conditions.
--halfradiusobs: radius derived from '--halfsumarea', underestimates r_e.
--fracradiusobs1: similar to --halfradiusobs, but for custom fraction 1.
--fracradiusobs2: similar to --halfradiusobs, but for custom fraction 2.
+ --areaminv: the number of pixels that are equal to the minimum value.
+ --areamaxv: the number of pixels that are equal to the maximum value.
- New columns to return the position of pixel with minimum or maximum
value: '--minvx', '--maxvx', '--minvy', '--maxvy', '--minvz',
'--maxvz'.
diff --git a/bin/arithmetic/arithmetic.c b/bin/arithmetic/arithmetic.c
index 09209e7..02b698e 100644
--- a/bin/arithmetic/arithmetic.c
+++ b/bin/arithmetic/arithmetic.c
@@ -703,14 +703,21 @@ arithmetic_invert(struct arithmeticparams *p, char *token)
#define INTERPOLATE_REGION(TYPE,OP,FUNC) { \
- TYPE mm, *a=in->array, *m=minmax->array; \
+ TYPE mm, b, *a=in->array, *m=minmax->array; \
FUNC(in->type, &mm); \
+ gal_blank_write(&b, in->type); \
for(i=0;i<minmax->size;++i) m[i]=mm; \
for(i=0;i<in->size;++i) \
{ \
if( l[i]>0 ) \
GAL_DIMENSION_NEIGHBOR_OP(i, in->ndim, in->dsize, in->ndim, \
- dinc, { if( a[nind] OP m[l[i]] ) m[l[i]]=a[nind]; }); \
+ dinc, \
+ { \
+ if(in->type==GAL_TYPE_FLOAT32 || in->type==GAL_TYPE_FLOAT64) \
+ { if( a[nind] OP m[l[i]] ) m[l[i]]=a[nind]; } \
+ else \
+ { if( a[nind]!=b && a[nind] OP m[l[i]] ) m[l[i]]=a[nind]; } \
+ }); \
} \
for(i=0;i<in->size;++i) if( l[i]>0 ) { a[i]=m[l[i]]; } \
}
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 03554fb..3a2bb8b 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -1477,6 +1477,34 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "areaminv",
+ UI_KEY_MINVNUM,
+ 0,
+ 0,
+ "Number of pixels with minimum value.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "areamaxv",
+ UI_KEY_MAXVNUM,
+ 0,
+ 0,
+ "Number of pixels with maximum value.",
+ UI_GROUP_COLUMNS_POSITION_PIXEL,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"surfacebrightness",
UI_KEY_SURFACEBRIGHTNESS,
0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index c31bf7c..0442fe0 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -726,12 +726,13 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Minimum value's X pixel position.";
ccomment = ocomment;
- otype = GAL_TYPE_UINT32;
- ctype = GAL_TYPE_UINT32;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MINVX ] = ciflag[ CCOL_MINVX ] = 1;
+ oiflag[ OCOL_MINVX ] = ciflag[ CCOL_MINVX ] = 1;
+ oiflag[ OCOL_MINVNUM ] = ciflag[ CCOL_MINVNUM ] = 1;
break;
case UI_KEY_MAXVX:
@@ -739,12 +740,13 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Maximum value's X pixel position.";
ccomment = ocomment;
- otype = GAL_TYPE_UINT32;
- ctype = GAL_TYPE_UINT32;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MAXVX ] = ciflag[ CCOL_MAXVX ] = 1;
+ oiflag[ OCOL_MAXVX ] = ciflag[ CCOL_MAXVX ] = 1;
+ oiflag[ OCOL_MAXVNUM ] = ciflag[ CCOL_MAXVNUM ] = 1;
break;
case UI_KEY_MINVY:
@@ -752,12 +754,13 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Minimum value's Y pixel position.";
ccomment = ocomment;
- otype = GAL_TYPE_UINT32;
- ctype = GAL_TYPE_UINT32;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MINVY ] = ciflag[ CCOL_MINVY ] = 1;
+ oiflag[ OCOL_MINVY ] = ciflag[ CCOL_MINVY ] = 1;
+ oiflag[ OCOL_MINVNUM ] = ciflag[ CCOL_MINVNUM ] = 1;
break;
case UI_KEY_MAXVY:
@@ -765,12 +768,13 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Maximum value's Y pixel position.";
ccomment = ocomment;
- otype = GAL_TYPE_UINT32;
- ctype = GAL_TYPE_UINT32;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MAXVY ] = ciflag[ CCOL_MAXVY ] = 1;
+ oiflag[ OCOL_MAXVY ] = ciflag[ CCOL_MAXVY ] = 1;
+ oiflag[ OCOL_MAXVNUM ] = ciflag[ CCOL_MAXVNUM ] = 1;
break;
case UI_KEY_MINVZ:
@@ -778,12 +782,13 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Minimum value's Z pixel position.";
ccomment = ocomment;
- otype = GAL_TYPE_UINT32;
- ctype = GAL_TYPE_UINT32;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MINVZ ] = ciflag[ CCOL_MINVZ ] = 1;
+ oiflag[ OCOL_MINVZ ] = ciflag[ CCOL_MINVZ ] = 1;
+ oiflag[ OCOL_MINVNUM ] = ciflag[ CCOL_MINVNUM ] = 1;
break;
case UI_KEY_MAXVZ:
@@ -791,12 +796,39 @@ columns_define_alloc(struct mkcatalogparams *p)
unit = "pixel";
ocomment = "Maximum value's Z pixel position.";
ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = 0;
+ disp_width = 10;
+ disp_precision = 0;
+ oiflag[ OCOL_MAXVZ ] = ciflag[ CCOL_MAXVZ ] = 1;
+ oiflag[ OCOL_MAXVNUM ] = ciflag[ CCOL_MAXVNUM ] = 1;
+ break;
+
+ case UI_KEY_MINVNUM:
+ name = "MIN_V_NUM";
+ unit = "counter";
+ ocomment = "Number of pixels with the minimum value.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_UINT32;
+ ctype = GAL_TYPE_UINT32;
+ disp_fmt = 0;
+ disp_width = 10;
+ disp_precision = 0;
+ oiflag[ OCOL_MINVNUM ] = ciflag[ CCOL_MINVNUM ] = 1;
+ break;
+
+ case UI_KEY_MAXVNUM:
+ name = "MAX_V_NUM";
+ unit = "counter";
+ ocomment = "Number of pixels with the maximum value..";
+ ccomment = ocomment;
otype = GAL_TYPE_UINT32;
ctype = GAL_TYPE_UINT32;
disp_fmt = 0;
disp_width = 10;
disp_precision = 0;
- oiflag[ OCOL_MAXVZ ] = ciflag[ CCOL_MAXVZ ] = 1;
+ oiflag[ OCOL_MAXVNUM ] = ciflag[ CCOL_MAXVNUM ] = 1;
break;
case UI_KEY_MINX:
@@ -2325,12 +2357,37 @@ columns_fill(struct mkcatalog_passparams *pp)
((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_C_GZ],
oi[OCOL_C_NUMALL] );
- case UI_KEY_MINVX: ((uint32_t *)colarr)[oind] = oi[OCOL_MINVX]; break;
- case UI_KEY_MAXVX: ((uint32_t *)colarr)[oind] = oi[OCOL_MAXVX]; break;
- case UI_KEY_MINVY: ((uint32_t *)colarr)[oind] = oi[OCOL_MINVY]; break;
- case UI_KEY_MAXVY: ((uint32_t *)colarr)[oind] = oi[OCOL_MAXVY]; break;
- case UI_KEY_MINVZ: ((uint32_t *)colarr)[oind] = oi[OCOL_MINVZ]; break;
- case UI_KEY_MAXVZ: ((uint32_t *)colarr)[oind] = oi[OCOL_MAXVZ]; break;
+ case UI_KEY_MINVX:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MINVX],
oi[OCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVX:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MAXVX],
oi[OCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVY:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MINVY],
oi[OCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVY:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MAXVY],
oi[OCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVZ:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MINVZ],
oi[OCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVZ:
+ ((float *)colarr)[oind] = MKC_RATIO( oi[OCOL_MAXVZ],
oi[OCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVNUM:
+ ((uint32_t *)colarr)[oind] = oi[OCOL_MINVNUM];
+ break;
+
+ case UI_KEY_MAXVNUM:
+ ((uint32_t *)colarr)[oind] = oi[OCOL_MAXVNUM];
+ break;
case UI_KEY_MINX:
case UI_KEY_MAXX:
@@ -2670,18 +2727,44 @@ columns_fill(struct mkcatalog_passparams *pp)
((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_GZ],
ci[CCOL_NUMALL] );
+ case UI_KEY_MINVX:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MINVX],
ci[CCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVX:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MAXVX],
ci[CCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVY:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MINVY],
ci[CCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVY:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MAXVY],
ci[CCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVZ:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MINVZ],
ci[CCOL_MINVNUM] );
+ break;
+
+ case UI_KEY_MAXVZ:
+ ((float *)colarr)[cind] = MKC_RATIO( ci[CCOL_MAXVZ],
ci[CCOL_MAXVNUM] );
+ break;
+
+ case UI_KEY_MINVNUM:
+ ((uint32_t *)colarr)[cind] = ci[CCOL_MINVNUM];
+ break;
+
+ case UI_KEY_MAXVNUM:
+ ((uint32_t *)colarr)[cind] = ci[CCOL_MAXVNUM];
+ break;
+
case UI_KEY_MINX: ((uint32_t *)colarr)[cind] = ci[CCOL_MINX];
break;
case UI_KEY_MAXX: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXX];
break;
case UI_KEY_MINY: ((uint32_t *)colarr)[cind] = ci[CCOL_MINY];
break;
case UI_KEY_MAXY: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXY];
break;
case UI_KEY_MINZ: ((uint32_t *)colarr)[cind] = ci[CCOL_MINZ];
break;
case UI_KEY_MAXZ: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXZ];
break;
- case UI_KEY_MINVX: ((uint32_t *)colarr)[cind] = ci[CCOL_MINVX];
break;
- case UI_KEY_MAXVX: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXVX];
break;
- case UI_KEY_MINVY: ((uint32_t *)colarr)[cind] = ci[CCOL_MINVY];
break;
- case UI_KEY_MAXVY: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXVY];
break;
- case UI_KEY_MINVZ: ((uint32_t *)colarr)[cind] = ci[CCOL_MINVZ];
break;
- case UI_KEY_MAXVZ: ((uint32_t *)colarr)[cind] = ci[CCOL_MAXVZ];
break;
case UI_KEY_W1:
case UI_KEY_W2:
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index e05ffc4..d4c9547 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -94,6 +94,8 @@ enum objectcols
OCOL_MAXVY, /* Y of maximum pixel in values file. */
OCOL_MINVZ, /* Z of minimum pixel in values file. */
OCOL_MAXVZ, /* Z of maximum pixel in values file. */
+ OCOL_MINVNUM, /* Number of pixels with minimum value. */
+ OCOL_MAXVNUM, /* Number of pixels with maximum value. */
OCOL_SUMSKY, /* Sum of sky value on this object. */
OCOL_NUMSKY, /* Number of sky value on this object. */
OCOL_SUMVAR, /* Sum of sky variance value on this object. */
@@ -159,6 +161,8 @@ enum clumpcols
CCOL_MAXVY, /* Y of maximum pixel in values array. */
CCOL_MINVZ, /* Z of minimum pixel in values array. */
CCOL_MAXVZ, /* Z of maximum pixel in values array. */
+ CCOL_MINVNUM, /* Number of pixels with minimum value. */
+ CCOL_MAXVNUM, /* Number of pixels with maximum value. */
CCOL_SUMSKY, /* Sum of sky value on this clump. */
CCOL_NUMSKY, /* Number of sky value on this clump. */
CCOL_SUMVAR, /* Sum of sky variance value on this clump. */
diff --git a/bin/mkcatalog/parse.c b/bin/mkcatalog/parse.c
index ab7744a..903147d 100644
--- a/bin/mkcatalog/parse.c
+++ b/bin/mkcatalog/parse.c
@@ -440,16 +440,13 @@ parse_objects(struct mkcatalog_passparams *pp)
double *oi=pp->oi;
gal_data_t *xybin=NULL;
- size_t maxima_c[3]={0,0,0};
size_t *tsize=pp->tile->dsize;
uint8_t *u, *uf, goodvalue, *xybinarr=NULL;
+ double minima_v=FLT_MAX, maxima_v=-FLT_MAX;
size_t d, pind=0, increment=0, num_increment=1;
- double minima_v[3]={ FLT_MAX, FLT_MAX, FLT_MAX};
- double maxima_v[3]={-FLT_MAX, -FLT_MAX, -FLT_MAX};
int32_t *O, *OO, *C=NULL, *objarr=p->objects->array;
float var, sval, varval, skyval, *V=NULL, *SK=NULL, *ST=NULL;
float *std=p->std?p->std->array:NULL, *sky=p->sky?p->sky->array:NULL;
- size_t minima_c[3]={GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T, GAL_BLANK_SIZE_T};
/* If tile processing isn't necessary, set 'tid' to a blank value. */
size_t tid = ( ( (p->sky && p->sky->size>1 && pp->st_sky == NULL )
@@ -480,6 +477,8 @@ parse_objects(struct mkcatalog_passparams *pp)
|| oif[ OCOL_MAXVY ]
|| oif[ OCOL_MINVZ ]
|| oif[ OCOL_MAXVZ ]
+ || oif[ OCOL_MINVNUM ]
+ || oif[ OCOL_MAXVNUM ]
|| sc
/* When the sky and its STD are tiles, we'll also need
the coordinate to find which tile a pixel belongs
@@ -589,19 +588,49 @@ parse_objects(struct mkcatalog_passparams *pp)
if(oif[ OCOL_C_SUM ]) oi[ OCOL_C_SUM ] += *V;
}
- /* Get the extrema of the values. */
- if( oif[ OCOL_MINVX ] && *V<minima_v[0] )
- { minima_v[0] = *V; minima_c[0] = c[ ndim-1 ]; }
- if( oif[ OCOL_MAXVX ] && *V>maxima_v[0] )
- { maxima_v[0] = *V; maxima_c[0] = c[ ndim-1 ]; }
- if( oif[ OCOL_MINVY ] && *V<minima_v[1] )
- { minima_v[1] = *V; minima_c[1] = c[ ndim-2 ]; }
- if( oif[ OCOL_MAXVY ] && *V>maxima_v[1] )
- { maxima_v[1] = *V; maxima_c[1] = c[ ndim-2 ]; }
- if( oif[ OCOL_MINVZ ] && *V<minima_v[2] )
- { minima_v[2] = *V; minima_c[2] = c[ ndim-3 ]; }
- if( oif[ OCOL_MAXVZ ] && *V>maxima_v[2] )
- { maxima_v[2] = *V; maxima_c[2] = c[ ndim-3 ]; }
+ /* Get the extrema of the values. Note that if the minima
+ or maxima value's coordinates are requested in any
+ dimension, then 'OCOL_MINVNUM' or 'OCOL_MAXVNUM' will
+ be activated). */
+ if( oif[ OCOL_MINVNUM ] && *V<=minima_v )
+ {
+ /* If the value is smaller than the smallest found so
+ far, reset the counter to one, and reset the sum
+ of positions this one's position. */
+ if( *V<minima_v )
+ {
+ minima_v = *V;
+ oi[ OCOL_MINVNUM ]=1;
+ if(oif[OCOL_MINVX]) oi[ OCOL_MINVX ] = c[ ndim-1 ]+1;
+ if(oif[OCOL_MINVY]) oi[ OCOL_MINVY ] = c[ ndim-2 ]+1;
+ if(oif[OCOL_MINVZ]) oi[ OCOL_MINVZ ] = c[ ndim-3 ]+1;
+ }
+ else
+ {
+ oi[ OCOL_MINVNUM ]++;
+ if(oif[OCOL_MINVX]) oi[ OCOL_MINVX ] += c[ ndim-1
]+1;
+ if(oif[OCOL_MINVY]) oi[ OCOL_MINVY ] += c[ ndim-2
]+1;
+ if(oif[OCOL_MINVZ]) oi[ OCOL_MINVZ ] += c[ ndim-3
]+1;
+ }
+ }
+ if( oif[ OCOL_MAXVNUM ] && *V>=maxima_v )
+ {
+ if( *V>maxima_v )
+ {
+ maxima_v = *V;
+ oi[ OCOL_MAXVNUM ]=1;
+ if(oif[OCOL_MAXVX]) oi[ OCOL_MAXVX ] = c[ ndim-1 ]+1;
+ if(oif[OCOL_MAXVY]) oi[ OCOL_MAXVY ] = c[ ndim-2 ]+1;
+ if(oif[OCOL_MAXVZ]) oi[ OCOL_MAXVZ ] = c[ ndim-3 ]+1;
+ }
+ else
+ {
+ oi[ OCOL_MAXVNUM ]++;
+ if(oif[OCOL_MAXVX]) oi[ OCOL_MAXVX ] += c[ ndim-1
]+1;
+ if(oif[OCOL_MAXVY]) oi[ OCOL_MAXVY ] += c[ ndim-2
]+1;
+ if(oif[OCOL_MAXVZ]) oi[ OCOL_MAXVZ ] += c[ ndim-3
]+1;
+ }
+ }
/* For flux weighted centers, we can only use positive
values, so do those measurements here. */
@@ -696,14 +725,6 @@ parse_objects(struct mkcatalog_passparams *pp)
pind=0;
}
- /* Write the extrema. */
- if( oif[ OCOL_MINVX ] ) oi[ OCOL_MINVX ] = minima_c[0] + 1;
- if( oif[ OCOL_MAXVX ] ) oi[ OCOL_MAXVX ] = maxima_c[0] + 1;
- if( oif[ OCOL_MINVY ] ) oi[ OCOL_MINVY ] = minima_c[1] + 1;
- if( oif[ OCOL_MAXVY ] ) oi[ OCOL_MAXVY ] = maxima_c[1] + 1;
- if( oif[ OCOL_MINVZ ] ) oi[ OCOL_MINVZ ] = minima_c[2] + 1;
- if( oif[ OCOL_MAXVZ ] ) oi[ OCOL_MAXVZ ] = maxima_c[2] + 1;
-
/* Write the projected area columns. */
if(xybin)
{
@@ -798,7 +819,6 @@ parse_clumps(struct mkcatalog_passparams *pp)
gal_data_t *xybin=NULL;
int32_t *O, *OO, *C=NULL, nlab;
size_t cind, *tsize=pp->tile->dsize;
- size_t *minima_c=NULL, *maxima_c=NULL;
double *minima_v=NULL, *maxima_v=NULL;
uint8_t *u, *uf, goodvalue, *cif=p->ciflag;
size_t nngb=gal_dimension_num_neighbors(ndim);
@@ -837,6 +857,8 @@ parse_clumps(struct mkcatalog_passparams *pp)
|| cif[ CCOL_MAXVY ]
|| cif[ CCOL_MINVZ ]
|| cif[ CCOL_MAXVZ ]
+ || cif[ CCOL_MINVNUM ]
+ || cif[ CCOL_MAXVNUM ]
|| sc
|| tid==GAL_BLANK_SIZE_T )
? gal_pointer_allocate(GAL_TYPE_SIZE_T, ndim, 0, __func__,
@@ -865,16 +887,12 @@ parse_clumps(struct mkcatalog_passparams *pp)
}
/* For the extrema columns. */
- if( cif[ CCOL_MINVX ] || cif[ CCOL_MINVX ] || cif[ CCOL_MINVZ ] )
- {
- minima_c=parse_init_extrema(cif, GAL_TYPE_SIZE_T, ndim*pp->clumpsinobj,
0);
- minima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, ndim*pp->clumpsinobj,
0);
- }
- if( cif[ CCOL_MAXVX ] || cif[ CCOL_MAXVY ] || cif[ CCOL_MAXVZ ] )
- {
- maxima_c=parse_init_extrema(cif, GAL_TYPE_SIZE_T, ndim*pp->clumpsinobj,
1);
- maxima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, ndim*pp->clumpsinobj,
1);
- }
+ if( cif[ CCOL_MINVNUM ] || cif[ CCOL_MINVX ]
+ || cif[ CCOL_MINVX ] || cif[ CCOL_MINVZ ] )
+ minima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, pp->clumpsinobj, 0);
+ if( cif[ CCOL_MAXVNUM ] || cif[ CCOL_MAXVX ]
+ || cif[ CCOL_MAXVY ] || cif[ CCOL_MAXVZ ] )
+ maxima_v=parse_init_extrema(cif, GAL_TYPE_FLOAT64, pp->clumpsinobj, 1);
/* Parse each contiguous patch of memory covered by this object. */
while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
@@ -967,25 +985,43 @@ parse_clumps(struct mkcatalog_passparams *pp)
if(cif[ CCOL_NUMXY ])
((uint8_t *)(xybin[cind].array))[ pind ] = 2;
- /* Extrema columns. */
- if( cif[ CCOL_MINVX ] && *V<minima_v[ cind*ndim+0 ] )
- { minima_v[ cind*ndim+0 ] = *V;
- minima_c[ cind*ndim+0 ] = c[ ndim-1 ]; }
- if( cif[ CCOL_MAXVX ] && *V>maxima_v[ cind*ndim+0 ] )
- { maxima_v[ cind*ndim+0 ] = *V;
- maxima_c[ cind*ndim+0 ] = c[ ndim-1 ]; }
- if( cif[ CCOL_MINVY ] && *V<minima_v[ cind*ndim+1 ] )
- { minima_v[ cind*ndim+1 ] = *V;
- minima_c[ cind*ndim+1 ] = c[ ndim-2 ]; }
- if( cif[ CCOL_MAXVY ] && *V>maxima_v[ cind*ndim+1 ] )
- { maxima_v[ cind*ndim+1 ] = *V;
- maxima_c[ cind*ndim+1 ] = c[ ndim-2 ]; }
- if( cif[ CCOL_MINVZ ] && *V<minima_v[ cind*ndim+2 ] )
- { minima_v[ cind*ndim+2 ] = *V;
- minima_c[ cind*ndim+2 ] = c[ ndim-3 ]; }
- if( cif[ CCOL_MAXVZ ] && *V>maxima_v[ cind*ndim+2 ] )
- { maxima_v[ cind*ndim+2 ] = *V;
- maxima_c[ cind*ndim+2 ] = c[ ndim-3 ]; }
+ /* Minimum/maximum pixel positions. */
+ if( cif[ CCOL_MINVNUM ] && *V<=minima_v[cind] )
+ {
+ if( *V<minima_v[cind] )
+ {
+ minima_v[cind] = *V;
+ ci[ CCOL_MINVNUM ]=1;
+ if(cif[CCOL_MINVX]) ci[ CCOL_MINVX ] = c[ ndim-1
]+1;
+ if(cif[CCOL_MINVY]) ci[ CCOL_MINVY ] = c[ ndim-2
]+1;
+ if(cif[CCOL_MINVZ]) ci[ CCOL_MINVZ ] = c[ ndim-3
]+1;
+ }
+ else
+ {
+ ci[ CCOL_MINVNUM ]++;
+ if(cif[CCOL_MINVX]) ci[ CCOL_MINVX ] += c[
ndim-1 ]+1;
+ if(cif[CCOL_MINVY]) ci[ CCOL_MINVY ] += c[
ndim-2 ]+1;
+ if(cif[CCOL_MINVZ]) ci[ CCOL_MINVZ ] += c[
ndim-3 ]+1;
+ }
+ }
+ if( cif[ CCOL_MAXVNUM ] && *V>=maxima_v[cind] )
+ {
+ if( *V>maxima_v[cind] )
+ {
+ maxima_v[cind] = *V;
+ ci[ CCOL_MAXVNUM ]=1;
+ if(cif[CCOL_MAXVX]) ci[ CCOL_MAXVX ] = c[ ndim-1
]+1;
+ if(cif[CCOL_MAXVY]) ci[ CCOL_MAXVY ] = c[ ndim-2
]+1;
+ if(cif[CCOL_MAXVZ]) ci[ CCOL_MAXVZ ] = c[ ndim-3
]+1;
+ }
+ else
+ {
+ ci[ CCOL_MAXVNUM ]++;
+ if(cif[CCOL_MAXVX]) ci[ CCOL_MAXVX ] += c[
ndim-1 ]+1;
+ if(cif[CCOL_MAXVY]) ci[ CCOL_MAXVY ] += c[
ndim-2 ]+1;
+ if(cif[CCOL_MAXVZ]) ci[ CCOL_MAXVZ ] += c[
ndim-3 ]+1;
+ }
+ }
/* Columns that need positive values. */
if( *V > 0.0f )
@@ -1155,27 +1191,17 @@ parse_clumps(struct mkcatalog_passparams *pp)
gal_fits_img_write(&xybin[i], "xybin.fits", NULL, NULL);
}
-
- /* Write the position of the maximum values. */
- if( cif[ CCOL_MINVX ] ) ci[ CCOL_MINVX ] = minima_c[ i*ndim + 0 ] + 1;
- if( cif[ CCOL_MAXVX ] ) ci[ CCOL_MAXVX ] = maxima_c[ i*ndim + 0 ] + 1;
- if( cif[ CCOL_MINVY ] ) ci[ CCOL_MINVY ] = minima_c[ i*ndim + 1 ] + 1;
- if( cif[ CCOL_MAXVY ] ) ci[ CCOL_MAXVY ] = maxima_c[ i*ndim + 1 ] + 1;
- if( cif[ CCOL_MINVZ ] ) ci[ CCOL_MINVZ ] = minima_c[ i*ndim + 2 ] + 1;
- if( cif[ CCOL_MAXVZ ] ) ci[ CCOL_MAXVZ ] = maxima_c[ i*ndim + 2 ] + 1;
}
/* Clean up. */
- if(c) free(c);
- if(sc) free(sc);
- if(dinc) free(dinc);
+ if(c) free(c);
+ if(sc) free(sc);
+ if(dinc) free(dinc);
if(ngblabs) free(ngblabs);
if(minima_v) free(minima_v);
- if(minima_c) free(minima_c);
if(maxima_v) free(maxima_v);
- if(maxima_c) free(maxima_c);
- if(xybin) gal_data_array_free(xybin, pp->clumpsinobj, 1);
+ if(xybin) gal_data_array_free(xybin, pp->clumpsinobj, 1);
}
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index e144028..4dbe04a 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -984,6 +984,8 @@ ui_necessary_inputs(struct mkcatalogparams *p, int *values,
int *sky,
case OCOL_MAXVY: *values = 1; break;
case OCOL_MINVZ: *values = 1; break;
case OCOL_MAXVZ: *values = 1; break;
+ case OCOL_MINVNUM: *values = 1; break;
+ case OCOL_MAXVNUM: *values = 1; break;
case OCOL_SUMSKY: *sky = 1; break;
case OCOL_SUMVAR: *std = 1; break;
case OCOL_SUMWHT: *values = 1; break;
@@ -1054,6 +1056,8 @@ ui_necessary_inputs(struct mkcatalogparams *p, int
*values, int *sky,
case CCOL_MAXVX: *values = 1; break;
case CCOL_MAXVY: *values = 1; break;
case CCOL_MAXVZ: *values = 1; break;
+ case CCOL_MINVNUM: *values = 1; break;
+ case CCOL_MAXVNUM: *values = 1; break;
case CCOL_SUMSKY: *sky = 1; break;
case CCOL_SUMVAR: *std = 1; break;
case CCOL_SUMWHT: *values = 1; break;
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 9db5c3d..18bd57c 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -124,6 +124,8 @@ enum option_keys_enum
UI_KEY_MAXVY,
UI_KEY_MINVZ,
UI_KEY_MAXVZ,
+ UI_KEY_MINVNUM,
+ UI_KEY_MAXVNUM,
UI_KEY_CLUMPSX,
UI_KEY_CLUMPSY,
UI_KEY_CLUMPSZ,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 5d21df9..45ff007 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -16426,6 +16426,12 @@ The used (non-blank in values image) area of the
labeled region in units of arcs
This column is just the value of the @option{--area} column, multiplied by the
area of each pixel in the input image (in units of arcsec^2).
Similar to the @option{--ra} or @option{--dec} columns, for this option to
work, the objects extension used has to have a WCS structure.
+@item --areaminv
+The number of pixels that are equal to the minimum value of the labeled region
(clump or object).
+
+@item --areamaxv
+The number of pixels that are equal to the maximum value of the labeled region
(clump or object).
+
@item --surfacebrightness
The surface brightness (in units of mag/arcsec@mymath{^2}) of the labeled
region (objects or clumps).
For more on the definition of the surface brightness, see @ref{Brightness flux
magnitude}.
| [Prev in Thread] |
Current Thread |
[Next in Thread] |
- [gnuastro-commits] master 8161798: MakeCatalog: positions of min and max value account for many pixels,
Mohammad Akhlaghi <=