[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 5c1860f 1/3: MakeCatalog can calculate the med
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 5c1860f 1/3: MakeCatalog can calculate the median value |
|
Date: |
Fri, 23 Feb 2018 10:00:24 -0500 (EST) |
branch: master
commit 5c1860f969e3988962c55a41e907847163e26743
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
MakeCatalog can calculate the median value
For calculating the surface brightness limit of extended and diffuse
objects, the average (brightness/area) can be too strongly affected by
outliers. So MakeCatalogs will now also calculate the median of pixel
values.
---
NEWS | 3 ++
bin/mkcatalog/args.h | 16 +++++-
bin/mkcatalog/columns.c | 29 +++++++++++
bin/mkcatalog/main.h | 6 ++-
bin/mkcatalog/mkcatalog.c | 125 +++++++++++++++++++++++++++++++++++++++++++++-
bin/mkcatalog/ui.h | 1 +
doc/gnuastro.texi | 37 +++++++-------
7 files changed, 196 insertions(+), 21 deletions(-)
diff --git a/NEWS b/NEWS
index 19c3ca3..b0ae1d5 100644
--- a/NEWS
+++ b/NEWS
@@ -9,6 +9,9 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
cumulative frequency plots to be set outside the minimum or maximum
values of the dataset.
+ MakeCatalog: The new `--median' option will calculate the median pixel
+ value within an object or clump.
+
** Removed features
** Changed features
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 53b67ae..7bd1cb1 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -703,6 +703,20 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "median",
+ UI_KEY_MEDIAN,
+ 0,
+ 0,
+ "Median of values in object/clump.",
+ UI_GROUP_COLUMNS_BRIGHTNESS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"magnitude",
UI_KEY_MAGNITUDE,
0,
@@ -846,7 +860,7 @@ struct argp_option program_options[] =
- /* Morphology/shapre related columns. */
+ /* Morphology/shape related columns. */
{
0, 0, 0, 0,
"Morphology/shape related columns",
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index fcfa30c..3cdb618 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -660,6 +660,24 @@ columns_define_alloc(struct mkcatalogparams *p)
ciflag[ CCOL_SUM ] = 1;
break;
+ case UI_KEY_MEDIAN:
+ name = "MEDIAN";
+ unit = p->input->unit ? p->input->unit : "pixelunit";
+ ocomment = "Median of sky subtracted values.";
+ ccomment = "Median of pixels subtracted by rivers.";
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 10;
+ disp_precision = 4;
+ oiflag[ OCOL_MEDIAN ] = 1;
+ oiflag[ OCOL_NUMALL ] = 1;
+ ciflag[ CCOL_MEDIAN ] = 1;
+ ciflag[ CCOL_NUMALL ] = 1;
+ ciflag[ CCOL_RIV_NUM ] = 1;
+ ciflag[ CCOL_RIV_SUM ] = 1;
+ break;
+
case UI_KEY_MAGNITUDE:
name = "MAGNITUDE";
unit = "log";
@@ -1357,6 +1375,12 @@ columns_fill(struct mkcatalog_passparams *pp)
: NAN );
break;
+ case UI_KEY_MEDIAN:
+ ((float *)colarr)[oind] = ( oi[ OCOL_NUM ]>0.0f
+ ? oi[ OCOL_MEDIAN ]
+ : NAN );
+ break;
+
case UI_KEY_MAGNITUDE:
((float *)colarr)[oind] = MKC_MAG(oi[ OCOL_SUM ]);
break;
@@ -1515,6 +1539,11 @@ columns_fill(struct mkcatalog_passparams *pp)
? ci[ CCOL_SUM ] : NAN );
break;
+ case UI_KEY_MEDIAN:
+ ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
+ ? ci[ CCOL_MEDIAN ] : NAN );
+ break;
+
case UI_KEY_MAGNITUDE: /* Similar: brightness for clumps */
tmp = ( ci[ CCOL_RIV_NUM ]
? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 477cc7f..0792db4 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -69,6 +69,7 @@ enum objectcols
OCOL_NUMALL, /* Number of all pixels with this label. */
OCOL_NUM, /* Number of pixels above threshold. */
OCOL_SUM, /* Sum of (value-sky) in object. */
+ OCOL_MEDIAN, /* Median of value in object. */
OCOL_VX, /* Sum of (value-sky) * x. */
OCOL_VY, /* Sum of (value-sky) * y. */
OCOL_SX, /* Shift along X axis. */
@@ -88,7 +89,7 @@ enum objectcols
OCOL_UPPERLIMIT_B, /* Upper limit magnitude. */
OCOL_C_NUMALL, /* Value independent no. of pixels in clumps.*/
OCOL_C_NUM, /* Area of clumps in this object. */
- OCOL_C_SUM, /* Brightness in object clumps. */
+ OCOL_C_SUM, /* Brightness in object clumps. */
OCOL_C_VX, /* Sum of (value-sky)*x on clumps. */
OCOL_C_VY, /* Sum of (value-sky)*y on obj. clumps. */
OCOL_C_GX, /* Geometric center of clumps in object X. */
@@ -109,6 +110,7 @@ enum clumpcols
CCOL_VYY, /* Sum of flux*y*y of this clump. */
CCOL_VXY, /* Sum of flux*x*y of this clump. */
CCOL_SUM, /* River subtracted brightness. */
+ CCOL_MEDIAN, /* Median of values in clump. */
CCOL_RIV_SUM, /* Sum of rivers around clump. */
CCOL_RIV_NUM, /* Num river pixels around this clump. */
CCOL_SUMSKY, /* Sum of sky value on this object. */
@@ -156,7 +158,7 @@ struct mkcatalogparams
char *upmaskfile; /* Name of upper limit mask file. */
char *upmaskhdu; /* HDU of upper limit mask file. */
size_t upnum; /* Number of upper-limit random samples.*/
- size_t *uprange; /* Range of random positions about target. */
+ size_t *uprange; /* Range of random positions about target.*/
uint8_t envseed; /* Use the environment for random seed. */
double upsigmaclip[2]; /* Sigma clip to measure upper limit. */
float upnsigma; /* Multiple of sigma to define up-lim. */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 08cdded..3f9ea99 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -37,6 +37,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/fits.h>
#include <gnuastro/threads.h>
#include <gnuastro/dimension.h>
+#include <gnuastro/statistics.h>
#include <gnuastro-internal/timing.h>
@@ -405,6 +406,124 @@ mkcatalog_second_pass(struct mkcatalog_passparams *pp)
+static void
+mkcatalog_median_pass(struct mkcatalog_passparams *pp)
+{
+ struct mkcatalogparams *p=pp->p;
+ size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+
+ double *ci;
+ gal_data_t *median;
+ int32_t *O, *C=NULL;
+ gal_data_t **clumpsmed=NULL;
+ float ss, *I, *II, *SK, *ST;
+ size_t i, increment=0, num_increment=1;
+ size_t counter=0, *ccounter=NULL, tsize=pp->oi[OCOL_NUMALL];
+ gal_data_t *objmed=gal_data_alloc(NULL, p->input->type, 1, &tsize, NULL, 0,
+ p->cp.minmapsize, NULL, NULL, NULL);
+
+ /* A small sanity check. */
+ if(p->input->type!=GAL_TYPE_FLOAT32)
+ error(EXIT_FAILURE, 0, "%s: the input has to be float32 type", __func__);
+
+
+ /* Allocate space for the clump medians. */
+ if(p->clumps)
+ {
+ errno=0;
+ clumpsmed=malloc(pp->clumpsinobj * sizeof *clumpsmed);
+ if(clumpsmed==NULL)
+ error(EXIT_FAILURE, errno, "%s: couldn't allocate `clumpsmed' for "
+ "%zu clumps", __func__, pp->clumpsinobj);
+
+
+ /* Allocate the array necessary to keep the values of each clump. */
+ ccounter=gal_data_calloc_array(GAL_TYPE_SIZE_T, pp->clumpsinobj,
+ __func__, "ccounter");
+ for(i=0;i<pp->clumpsinobj;++i)
+ {
+ tsize=pp->ci[ i * CCOL_NUMCOLS + CCOL_NUMALL ];
+ clumpsmed[i]=gal_data_alloc(NULL, p->input->type, 1, &tsize, NULL,
+ 0, p->cp.minmapsize, NULL, NULL, NULL);
+ }
+ }
+
+
+ /* Parse each contiguous patch of memory covered by this object. */
+ while( pp->start_end_inc[0] + increment <= pp->start_end_inc[1] )
+ {
+ /* Set the contiguous range to parse, we will check the count
+ over the `I' pointer and just increment the rest. */
+ O = pp->st_o + increment;
+ SK = pp->st_sky + increment;
+ ST = pp->st_std + increment;
+ if(p->clumps) C = pp->st_c + increment;
+ II = ( I = pp->st_i + increment ) + pp->tile->dsize[ndim-1];
+
+ /* Parse the next contiguous region of this tile. */
+ do
+ {
+ /* If this pixel belongs to the requested object, is a clumps and
+ isn't NAN, then do the processing. `hasblank' is constant, so
+ when the input doesn't have any blank values, the `isnan' will
+ never be checked. */
+ if( *O==pp->object && !isnan(*I) )
+ {
+ /* Copy the value for the whole object. */
+ ss = *I - *SK;
+ memcpy( gal_data_ptr_increment(objmed->array, counter++,
+ p->input->type),
+ &ss, gal_type_sizeof(p->input->type) );
+
+ /* We are also on a clump. */
+ if(p->clumps && *C>0)
+ memcpy( gal_data_ptr_increment(clumpsmed[*C-1]->array,
+ ccounter[*C-1]++,
+ p->input->type),
+ &ss, gal_type_sizeof(p->input->type) );
+ }
+
+ /* Increment the other pointers. */
+ ++O; ++SK; ++ST; if(p->clumps) ++C;
+ }
+ while(++I<II);
+
+ /* Increment to the next contiguous region of this tile. */
+ increment += ( gal_tile_block_increment(p->input, dsize,
+ num_increment++, NULL) );
+ }
+
+
+ /* Calculate the final medians for objects. */
+ median=gal_data_copy_to_new_type_free(gal_statistics_median(objmed, 1),
+ GAL_TYPE_FLOAT64);
+ pp->oi[OCOL_MEDIAN]=*((double *)(median->array));
+ gal_data_free(objmed);
+ gal_data_free(median);
+
+
+ /* Calculate the median for clumps. */
+ if(p->clumps)
+ {
+ for(i=0;i<pp->clumpsinobj;++i)
+ {
+ ci=&pp->ci[ i * CCOL_NUMCOLS ];
+ median=gal_statistics_median(clumpsmed[i], 1);
+ median=gal_data_copy_to_new_type_free(median, GAL_TYPE_FLOAT64);
+ ci[ CCOL_MEDIAN ] = ( *((double *)(median->array))
+ - (ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]) );
+ gal_data_free(clumpsmed[i]);
+ gal_data_free(median);
+ }
+ free(clumpsmed);
+ free(ccounter);
+ }
+}
+
+
+
+
+
@@ -492,7 +611,7 @@ mkcatalog_single_object(void *in_prm)
mkcatalog_first_pass(&pp);
/* Currently the second pass is only necessary when there is a clumps
- image. */
+ image or the median is requested. */
if(p->clumps)
{
/* Allocate space for the properties of each clump. */
@@ -509,6 +628,10 @@ mkcatalog_single_object(void *in_prm)
mkcatalog_second_pass(&pp);
}
+ /* If the median is requested, another pass is necessary. */
+ if( p->oiflag[ OCOL_MEDIAN ] )
+ mkcatalog_median_pass(&pp);
+
/* Calculate the upper limit magnitude (if necessary). */
if(p->upperlimit) upperlimit_calculate(&pp);
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index e5f6d15..27da86c 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -115,6 +115,7 @@ enum option_keys_enum
UI_KEY_CLUMPSGEOW2,
UI_KEY_CLUMPSBRIGHTNESS,
UI_KEY_NORIVERBRIGHTNESS,
+ UI_KEY_MEDIAN,
UI_KEY_CLUMPSMAGNITUDE,
UI_KEY_UPPERLIMIT,
UI_KEY_RIVERAVE,
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a4b380d..58d61f3 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15046,7 +15046,7 @@ variable to the @code{mkcatalogparams} structure.
@item ui.h
The @code{option_keys_enum} associates a unique value for each option to
-MakeProfiles. The options that have a short option version, the single
+MakeCatalog. The options that have a short option version, the single
character short comment is used for the value. Those that don't have a
short option version, get a large integer automatically. You should add a
variable here to identify your desired column.
@@ -15066,15 +15066,14 @@ sanity checks and preparations for it here.
Otherwise, you can ignore this
file.
@item columns.c
-This file will contain the main definition and high-level calculation of
-your new column through the @code{columns_define_alloc} and
address@hidden functions. In the first, you specify the basic
-information about the column: its name, units, comments, type (see
address@hidden data types}) and how it should be printed if the output is a
-text file. You should also specify the raw/internal columns that are
-necessary for this column here as the many existing examples show. Through
-the types for objects and rows, you can specify if this column is only for
-clumps, objects or both.
+This file contains the main definition and high-level calculation of your
+new column through the @code{columns_define_alloc} and @code{columns_fill}
+functions. In the first, you specify the basic information about the
+column: its name, units, comments, type (see @ref{Numeric data types}) and
+how it should be printed if the output is a text file. You should also
+specify the raw/internal columns that are necessary for this column here as
+the many existing examples show. Through the types for objects and rows,
+you can specify if this column is only for clumps, objects or both.
The second main function (@code{columns_fill}) writes the final value into
the appropriate column for each object and clump. As you can see in the
@@ -15595,6 +15594,11 @@ usable pixels (blank or below the threshold) are
present over the clump or
object, the stored value will be NaN, because zero (note that zero is
meaningful).
address@hidden --median
+The median sky subtracted value of pixels within the object or clump. For
+clumps, the average river flux is subtracted from the sky subtracted
+median.
+
@item --noriverbrightness
[Clumps] The Sky (not river) subtracted clump brightness. By definition,
for the clumps, the average brightness of the rivers surrounding it are
@@ -22707,15 +22711,14 @@ use to access an element of the @code{OTHER} dataset
(if
(similar to @code{itype} and @code{INPUT} discussed above). If
@code{PARSE_OTHER} is non-zero, then @code{o} will also be incremented to
the same index element but in the other array. You can use these along with
-any other variable you define before this macro to process the input and
-store it in the other.
+any other variable you define before this macro to process the input and/or
+the other.
All variables within this function-like macro begin with @code{tpo_} except
-for the three variables listed above. So as long as you don't start the
-names of your variables with this prefix everything will be fine. Note that
address@hidden (and possibly @code{o}) will be incremented once by this
-function-like macro, so don't increment them within @code{OP}. As a
-summary, the three variables you have access to within @code{OP} are:
+for the three variables listed below. Therefore, as long as you don't start
+the names of your variables with this prefix everything will be fine. Note
+that @code{i} (and possibly @code{o}) will be incremented once by this
+function-like macro, so don't increment them within @code{OP}.
@table @code
@item i