[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 42f0405 2/3: Position in upperlimit distributi
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 42f0405 2/3: Position in upperlimit distribution (sigma and quantile) |
|
Date: |
Tue, 6 Mar 2018 09:20:21 -0500 (EST) |
branch: master
commit 42f040597d9748c083b0fde3d6fa336029ed470a
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Position in upperlimit distribution (sigma and quantile)
Three new options are now added to MakeCatalog: two will allow identifying
the position of an object or clump within the upper-limit distribution, and
another can be used to identify the upper limit distribution irrespective
of the `--upnsigma' value.
Also, while working on the quantile position, I noticed that the fact that
`gal_statistics_quantile_function' only returned a NaN when the value was
out of range was not too meaningful (is it smaller than the minimum or
larger than the maximum for example). So now, it returns `inf' or `-inf' to
specify in which direction the value is out of range.
---
NEWS | 11 +++++
bin/mkcatalog/args.h | 42 +++++++++++++++++
bin/mkcatalog/columns.c | 114 +++++++++++++++++++++++++++++++++++---------
bin/mkcatalog/main.h | 8 +++-
bin/mkcatalog/mkcatalog.c | 2 +-
bin/mkcatalog/mkcatalog.h | 4 +-
bin/mkcatalog/ui.h | 3 ++
bin/mkcatalog/upperlimit.c | 115 +++++++++++++++++++++++++++++++++++++--------
doc/gnuastro.texi | 33 +++++++++++--
lib/statistics.c | 88 +++++++++++++++++++++++++++-------
10 files changed, 350 insertions(+), 70 deletions(-)
diff --git a/NEWS b/NEWS
index 24a1aaa..0fd4fd6 100644
--- a/NEWS
+++ b/NEWS
@@ -8,6 +8,13 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
MakeCatalog: The new `--mean' and `--median' options will calculate the
mean and median pixel value within an object or clump respectively.
+ MakeCatalog: The new `--upperlimitsigma' and `--upperlimitquantile'
+ columns will report the postion of the object's brightness with respect
+ to the distribution of randomly measured values. The former as a multiple
+ of sigma and the latter as a quantile. Also, `--upperlimitonesigma' will
+ return the 1-sigma value of the randomly placed upper limit magnitudes
+ (irrespective of the value given to `--upnsigma').
+
NoiseChisel: a value of `none' to the `--kernel' option will disable
convolution.
@@ -19,6 +26,10 @@ GNU Astronomy Utilities NEWS -*-
outline -*-
** Changed features
+ Libraries: `gal_statistics_quantile_function' returns `inf' or `-inf' if
+ the given value is smaller than the minimum or larger than the maximum of
+ the input dataset's range.
+
** Bug fixes
Many unused result warnings for asprintf in some compilers (bug #52979).
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index c3a4a04..3a547fe 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -801,6 +801,48 @@ struct argp_option program_options[] =
ui_column_codes_ll
},
{
+ "upperlimitonesigma",
+ UI_KEY_UPPERLIMITONESIGMA,
+ 0,
+ 0,
+ "Upper-limit one sigma value.",
+ UI_GROUP_COLUMNS_BRIGHTNESS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "upperlimitsigma",
+ UI_KEY_UPPERLIMITSIGMA,
+ 0,
+ 0,
+ "Place in upperlimit distribution (sigma multiple).",
+ UI_GROUP_COLUMNS_BRIGHTNESS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
+ "upperlimitquantile",
+ UI_KEY_UPPERLIMITQUANTILE,
+ 0,
+ 0,
+ "Quantile in random distribution (max 1).",
+ UI_GROUP_COLUMNS_BRIGHTNESS,
+ 0,
+ GAL_TYPE_INVALID,
+ GAL_OPTIONS_RANGE_ANY,
+ GAL_OPTIONS_NOT_MANDATORY,
+ GAL_OPTIONS_NOT_SET,
+ ui_column_codes_ll
+ },
+ {
"riverave",
UI_KEY_RIVERAVE,
0,
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index bb539a7..aa354e5 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -275,9 +275,9 @@ columns_define_alloc(struct mkcatalogparams *p)
smaller domain of raw measurements. So to avoid having to calculate
something multiple times, each parameter will flag the intermediate
parameters it requires in these arrays. */
- oiflag = p->oiflag = gal_data_malloc_array(GAL_TYPE_UINT8, OCOL_NUMCOLS,
+ oiflag = p->oiflag = gal_data_calloc_array(GAL_TYPE_UINT8, OCOL_NUMCOLS,
__func__, "oiflag");
- ciflag = p->ciflag = gal_data_malloc_array(GAL_TYPE_UINT8, CCOL_NUMCOLS,
+ ciflag = p->ciflag = gal_data_calloc_array(GAL_TYPE_UINT8, CCOL_NUMCOLS,
__func__, "ciflag");
/* Allocate the columns. */
@@ -770,6 +770,45 @@ columns_define_alloc(struct mkcatalogparams *p)
p->hasmag = 1; /* Doesn't need per-pixel calculations. */
break;
+ case UI_KEY_UPPERLIMITONESIGMA:
+ name = "UPPERLIMIT_ONE_SIGMA";
+ unit = p->input->unit;
+ ocomment = "One sigma value of all random measurements.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 8;
+ disp_precision = 3;
+ p->upperlimit = 1;
+ break;
+
+ case UI_KEY_UPPERLIMITSIGMA:
+ name = "UPPERLIMIT_SIGMA";
+ unit = "frac";
+ ocomment = "Place in upperlimit distribution (sigma
multiple).";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 8;
+ disp_precision = 3;
+ p->upperlimit = 1;
+ break;
+
+ case UI_KEY_UPPERLIMITQUANTILE:
+ name = "UPPERLIMIT_QUANTILE";
+ unit = "quantile";
+ ocomment = "Quantile of brightness in random distribution.";
+ ccomment = ocomment;
+ otype = GAL_TYPE_FLOAT32;
+ ctype = GAL_TYPE_FLOAT32;
+ disp_fmt = GAL_TABLE_DISPLAY_FMT_GENERAL;
+ disp_width = 8;
+ disp_precision = 3;
+ p->upperlimit = 1;
+ break;
+
case UI_KEY_RIVERAVE:
name = "RIVER_AVE";
unit = p->input->unit ? p->input->unit : "pixelunit";
@@ -1227,6 +1266,27 @@ columns_second_order(struct mkcatalog_passparams *pp,
double *row,
+/* The clump brightness is needed in several places, so we've defined this
+ function to have an easier code. */
+static double
+columns_clump_brightness(double *ci)
+{
+ double tmp;
+ /* Calculate the river flux over the clump area. But only when rivers are
+ present. When grown clumps are requested, the clumps can fully cover a
+ detection (that has one or no clumps). */
+ tmp = ( ci[ CCOL_RIV_NUM ]>0.0f
+ ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
+ : 0 );
+
+ /* Subtract it from the clump's brightness. */
+ return ci[ CCOL_NUM ]>0.0f ? (ci[ CCOL_SUM ] - tmp) : NAN;
+}
+
+
+
+
+
/* The magnitude error is directly derivable from the S/N:
To derive the error in measuring the magnitude from the S/N, let's take
@@ -1425,6 +1485,20 @@ columns_fill(struct mkcatalog_passparams *pp)
((float *)colarr)[oind] = MKC_MAG(oi[ OCOL_UPPERLIMIT_B ]);
break;
+ case UI_KEY_UPPERLIMITONESIGMA:
+ ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_S ];
+ break;
+
+ case UI_KEY_UPPERLIMITSIGMA:
+ ((float *)colarr)[oind] = ( ( oi[ OCOL_NUM ]>0.0f
+ ? oi[ OCOL_SUM ] : NAN )
+ / oi[ OCOL_UPPERLIMIT_S ] );
+ break;
+
+ case UI_KEY_UPPERLIMITQUANTILE:
+ ((float *)colarr)[oind] = oi[ OCOL_UPPERLIMIT_Q ];
+ break;
+
case UI_KEY_SN:
((float *)colarr)[oind] = columns_sn(p, oi, 0);
break;
@@ -1545,17 +1619,7 @@ columns_fill(struct mkcatalog_passparams *pp)
break;
case UI_KEY_BRIGHTNESS:
- /* Calculate the river flux over the clump area. But only when
- rivers are present. When grown clumps are requested, the
- clumps can fully cover a detection (that has one or no
- clumps). */
- tmp = ( ci[ CCOL_RIV_NUM ]>0.0f
- ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]*ci[ CCOL_NUM ]
- : 0 );
-
- /* Subtract it from the clump's brightness. */
- ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
- ? (ci[ CCOL_SUM ] - tmp) : NAN );
+ ((float *)colarr)[cind] = columns_clump_brightness(ci);
break;
case UI_KEY_NORIVERBRIGHTNESS:
@@ -1564,16 +1628,8 @@ columns_fill(struct mkcatalog_passparams *pp)
break;
case UI_KEY_MEAN:
- /* Similar to brightness. */
- tmp = ( ci[ CCOL_RIV_NUM ]>0.0f
- ? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM ]
- : 0 );
-
- /* Subtract it from the clump's mean. */
- ((float *)colarr)[cind] = ( ci[ CCOL_NUM ]>0.0f
- ? (ci[CCOL_SUM]/ci[CCOL_NUM] - tmp)
- : NAN );
-
+ ((float *)colarr)[cind] = ( columns_clump_brightness(ci)
+ /ci[CCOL_NUM] );
break;
case UI_KEY_MEDIAN:
@@ -1600,6 +1656,18 @@ columns_fill(struct mkcatalog_passparams *pp)
((float *)colarr)[cind] = MKC_MAG(ci[ CCOL_UPPERLIMIT_B ]);
break;
+ case UI_KEY_UPPERLIMITONESIGMA:
+ ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_S ];
+ break;
+
+ case UI_KEY_UPPERLIMITSIGMA:
+ ((float *)colarr)[cind] = ( columns_clump_brightness(ci)
+ / ci[ CCOL_UPPERLIMIT_S ] );
+
+ case UI_KEY_UPPERLIMITQUANTILE:
+ ((float *)colarr)[cind] = ci[ CCOL_UPPERLIMIT_Q ];
+ break;
+
case UI_KEY_RIVERAVE:
((float *)colarr)[cind] = ( ci[ CCOL_RIV_NUM]
? ci[ CCOL_RIV_SUM ]/ci[ CCOL_RIV_NUM]
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 0792db4..5fb5ed0 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -86,7 +86,9 @@ enum objectcols
OCOL_GXX, /* Second order geometric variable: X*X. */
OCOL_GYY, /* Second order geometric variable: Y*Y. */
OCOL_GXY, /* Second order geometric variable: X*Y. */
- OCOL_UPPERLIMIT_B, /* Upper limit magnitude. */
+ OCOL_UPPERLIMIT_B, /* Upper limit brightness. */
+ OCOL_UPPERLIMIT_S, /* Upper limit one-sigma value. */
+ OCOL_UPPERLIMIT_Q, /* Quantile of object in random distribution.*/
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. */
@@ -122,7 +124,9 @@ enum clumpcols
CCOL_GXX, /* Second order geometric moment. */
CCOL_GYY, /* Second order geometric moment. */
CCOL_GXY, /* Second order geometric moment. */
- CCOL_UPPERLIMIT_B, /* Upper limit magnitude. */
+ CCOL_UPPERLIMIT_B, /* Upper limit brightness. */
+ CCOL_UPPERLIMIT_S, /* Upper limit one-sigma value. */
+ CCOL_UPPERLIMIT_Q, /* Quantile of object in random distribution.*/
CCOL_NUMCOLS, /* SHOULD BE LAST: total number of columns. */
};
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 517e8ce..29e0ae7 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -612,7 +612,7 @@ mkcatalog_single_object(void *in_prm)
mkcatalog_first_pass(&pp);
/* Currently the second pass is only necessary when there is a clumps
- image or the median is requested. */
+ image. */
if(p->clumps)
{
/* Allocate space for the properties of each clump. */
diff --git a/bin/mkcatalog/mkcatalog.h b/bin/mkcatalog/mkcatalog.h
index 959b309..4e93c35 100644
--- a/bin/mkcatalog/mkcatalog.h
+++ b/bin/mkcatalog/mkcatalog.h
@@ -33,11 +33,11 @@ struct mkcatalog_passparams
gal_data_t *tile; /* The tile to pass-over. */
float *st_i; /* Starting pointer for input image. */
int32_t *st_o; /* Starting pointer for objects image. */
- int32_t *st_c; /* Starting pointer for objects image. */
+ int32_t *st_c; /* Starting pointer for clumps image. */
float *st_sky; /* Starting pointer for input image. */
float *st_std; /* Starting pointer for input image. */
size_t start_end_inc[2]; /* Starting and ending indexs. */
- size_t *shift; /* Shift coordinates for coordinates. */
+ size_t *shift; /* Shift coordinates. */
gsl_rng *rng; /* Random number generator. */
size_t clumpstartindex; /* Clump starting row in final catalog. */
gal_data_t *up_vals; /* Container for upper-limit values. */
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 4586527..c161970 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -119,6 +119,9 @@ enum option_keys_enum
UI_KEY_MEDIAN,
UI_KEY_CLUMPSMAGNITUDE,
UI_KEY_UPPERLIMIT,
+ UI_KEY_UPPERLIMITONESIGMA,
+ UI_KEY_UPPERLIMITSIGMA,
+ UI_KEY_UPPERLIMITQUANTILE,
UI_KEY_RIVERAVE,
UI_KEY_RIVERNUM,
UI_KEY_SKY,
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
index c767888..645f695 100644
--- a/bin/mkcatalog/upperlimit.c
+++ b/bin/mkcatalog/upperlimit.c
@@ -34,6 +34,8 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/statistics.h>
#include "main.h"
+
+#include "ui.h"
#include "mkcatalog.h"
@@ -263,17 +265,103 @@ upperlimit_random_position(struct mkcatalog_passparams
*pp, gal_data_t *tile,
-static double
+/* Given the distribution of values, do the upper-limit calculations. */
+static void
+upperlimit_measure(struct mkcatalog_passparams *pp, int32_t clumplab,
+ int do_measurement)
+{
+ float *scarr;
+ gal_data_t *column;
+ size_t init_size, col, one=1;
+ struct mkcatalogparams *p=pp->p;
+ gal_data_t *sigclip=NULL, *sum, *qfunc;
+ double *o = ( clumplab
+ ? &pp->ci[ (clumplab-1) * CCOL_NUMCOLS ]
+ : pp->oi );
+
+ /* If the random distribution exsits, then fill it in. */
+ if(do_measurement)
+ {
+ /* These columns are for both objects and clumps, so if they are
+ requested in objects, they will also be written for clumps here
+ (the order is irrelevant here). */
+ for(column=p->objectcols; column!=NULL; column=column->next)
+ {
+ switch(column->status)
+ {
+ /* Columns that depend on the sigma of the distribution. */
+ case UI_KEY_UPPERLIMIT:
+ case UI_KEY_UPPERLIMITMAG:
+ case UI_KEY_UPPERLIMITONESIGMA:
+
+ /* We only need to do this once. */
+ if(sigclip==NULL)
+ {
+ /* Calculate the sigma-clipped standard deviation. Since
+ it is done in place, the size will change, so we'll
+ keep the size here and put it back after we are
+ done. */
+ init_size=pp->up_vals->size;
+ sigclip=gal_statistics_sigma_clip(pp->up_vals,
+ p->upsigmaclip[0],
+ p->upsigmaclip[1], 1, 1);
+ pp->up_vals->size=pp->up_vals->dsize[0]=init_size;
+ scarr=sigclip->array;
+
+ /* Write the raw sigma. */
+ col = clumplab ? CCOL_UPPERLIMIT_S : OCOL_UPPERLIMIT_S;
+ o[col] = scarr[3];
+
+ /* Write the multiple of `upnsigma'. */
+ col = clumplab ? CCOL_UPPERLIMIT_B : OCOL_UPPERLIMIT_B;
+ o[col] = scarr[3] * p->upnsigma;
+
+ /* Clean up. */
+ gal_data_free(sigclip);
+ }
+ break;
+
+ /* Quantile column. */
+ case UI_KEY_UPPERLIMITQUANTILE:
+
+ /* Similar to the case for sigma-clipping, we'll need to keep
+ the size here also. */
+ init_size=pp->up_vals->size;
+ sum=gal_data_alloc(NULL, GAL_TYPE_FLOAT32, 1, &one, NULL, 0,
+ -1, NULL, NULL, NULL);
+ ((float *)(sum->array))[0]=o[clumplab ? CCOL_SUM : OCOL_SUM];
+ qfunc=gal_statistics_quantile_function(pp->up_vals, sum, 1);
+
+ /* Fill in the column. */
+ col = clumplab ? CCOL_UPPERLIMIT_Q : OCOL_UPPERLIMIT_Q;
+ pp->up_vals->size=pp->up_vals->dsize[0]=init_size;
+ o[col] = ((double *)(qfunc->array))[0];
+ break;
+ }
+ }
+ }
+ else
+ {
+ o[ clumplab ? CCOL_UPPERLIMIT_S : OCOL_UPPERLIMIT_S ] = NAN;
+ o[ clumplab ? CCOL_UPPERLIMIT_B : OCOL_UPPERLIMIT_B ] = NAN;
+ o[ clumplab ? CCOL_UPPERLIMIT_Q : OCOL_UPPERLIMIT_Q ] = NAN;
+ }
+}
+
+
+
+
+
+static void
upperlimit_one_tile(struct mkcatalog_passparams *pp, gal_data_t *tile,
unsigned long seed, int32_t clumplab)
{
struct mkcatalogparams *p=pp->p;
size_t ndim=p->input->ndim, *dsize=p->input->dsize;
+ double sum;
void *tarray;
- double sum, out;
int continueparse;
- gal_data_t *sigclip;
uint8_t *M=NULL, *st_m=NULL;
float *uparr=pp->up_vals->array;
float *I, *II, *SK, *st_i, *st_sky;
@@ -375,19 +463,12 @@ upperlimit_one_tile(struct mkcatalog_passparams *pp,
gal_data_t *tile,
++tcounter;
}
- /* Calculate the standard deviation of this distribution. */
- if(counter==p->upnum)
- {
- sigclip=gal_statistics_sigma_clip(pp->up_vals, p->upsigmaclip[0],
- p->upsigmaclip[1], 1, 1);
- out = ((float *)(sigclip->array))[3] * p->upnsigma;
- }
- else out=NAN;
+ /* Do the measurement on the random distribution. */
+ upperlimit_measure(pp, clumplab, counter==p->upnum);
/* Reset the tile's array pointer, clean up and return. */
tile->array=tarray;
free(rcoord);
- return out;
}
@@ -416,16 +497,14 @@ void
upperlimit_calculate(struct mkcatalog_passparams *pp)
{
size_t i;
- double *ci;
unsigned long seed;
gal_data_t *clumptiles;
struct mkcatalogparams *p=pp->p;
/* First find the upper limit magnitude for this object. */
- pp->oi[OCOL_UPPERLIMIT_B] = upperlimit_one_tile(pp, pp->tile,
- p->seed+pp->object, 0);
+ upperlimit_one_tile(pp, pp->tile, p->seed+pp->object, 0);
- /* If a clumps image is present (a clump catalog is requested( and this
+ /* If a clumps image is present (a clump catalog is requested) and this
object has clumps, then find the upper limit magnitude for the clumps
within this object. */
if(p->clumps && pp->clumpsinobj)
@@ -440,10 +519,8 @@ upperlimit_calculate(struct mkcatalog_passparams *pp)
IDs. */
for(i=0;i<pp->clumpsinobj;++i)
{
- ci=&pp->ci[ i * CCOL_NUMCOLS ];
seed = p->seed + p->numobjects + p->numclumps * pp->object + i;
- ci[CCOL_UPPERLIMIT_B] = upperlimit_one_tile(pp, &clumptiles[i],
- seed, i+1);
+ upperlimit_one_tile(pp, &clumptiles[i], seed, i+1);
}
/* Clean up the clump tiles. */
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 137b2ba..c7c15fd 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -15669,13 +15669,34 @@ magnitude settings} for a complete explanation. This
is very important for
the fainter and smaller objects in the image where the measured magnitudes
are not reliable.
-
@item --upperlimitmag
The upper limit magnitude for this object or clump. See @ref{Quantifying
measurement limits} and @ref{Upper-limit magnitude settings} for a complete
explanation. This is very important for the fainter and smaller objects in
the image where the measured magnitudes are not reliable.
address@hidden --upperlimitonesigma
+The @mymath{1\sigma} upper limit value (in units of the input image) for
+this object or clump. See @ref{Quantifying measurement limits} and
address@hidden magnitude settings} for a complete explanation. When
address@hidden, this column's values will be the same as
address@hidden
+
address@hidden --upperlimitsigma
+The position of the total brightness measured within the distribution of
+randomly placed upperlimit measurements in units of the distribution's
address@hidden or standard deviation. See @ref{Quantifying measurement
+limits} and @ref{Upper-limit magnitude settings} for a complete
+explanation.
+
address@hidden --upperlimitquantile
+The position of the total brightness measured within the distribution of
+randomly placed upperlimit measurements as a quantile (value between 0 or
+1). See @ref{Quantifying measurement limits} and @ref{Upper-limit magnitude
+settings} for a complete explanation. If the object is brighter than the
+brightest randomly placed profile, a value of @code{inf} is returned. If it
+is less than the minimum, a value of @code{-inf} is reported.
+
@item --riverave
[Clumps] The average brightness of the river pixels around this
clump. River pixels were defined in Akhlaghi and Ichikawa 2015. In
@@ -23551,16 +23572,18 @@ datatype of the output is the same as @code{input}.
See
Return the index of the quantile function (inverse quantile) of
@code{input} at @code{value}. In other words, this function will return the
index of the nearest element (of a sorted and non-blank) @code{input} to
address@hidden See @code{gal_statistics_median} for a description of
address@hidden
address@hidden If the value is outside the range of the input, then this
+function will return @code{GAL_BLANK_SIZE_T}.
@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_quantile_function (gal_data_t
@code{*input}, gal_data_t @code{*value}, int @code{inplace})
Return a single-element (@code{double} or @code{float64}) dataset
containing the quantile function of the non-blank values in @code{input} at
@code{value}. In other words, this function will return the quantile of
address@hidden in @code{input}. See @code{gal_statistics_median} for a
-description of @code{inplace}.
address@hidden in @code{input}. If the value is smaller than the input's
+smallest element, the returned value will be zero. If the value is larger
+than the input's largest element, then the returned value is 1. See
address@hidden for a description of @code{inplace}.
@end deftypefun
@deftypefun {gal_data_t *} gal_statistics_mode (gal_data_t @code{*input},
float @code{mirrordist}, int @code{inplace})
diff --git a/lib/statistics.c b/lib/statistics.c
index ce71667..6c19f36 100644
--- a/lib/statistics.c
+++ b/lib/statistics.c
@@ -353,29 +353,44 @@ gal_statistics_quantile(gal_data_t *input, double
quantile, int inplace)
/* Return the index of the (first) point in the sorted dataset that has the
closest value to `value' (which has to be the same type as the `input'
dataset). */
-#define STATS_QFUNC(IT) { \
+#define STATS_QFUNC_IND(IT) { \
IT *r, *a=nbs->array, *af=a+nbs->size, v=*((IT *)(value->array)); \
\
- /* For a reference. Since we are comparing with the previous. */ \
+ /* For a reference. Since we are comparing with the previous */ \
/* element, we need to start with the second element.*/ \
r=a++; \
\
/* Increasing array: */ \
if( *a < *(a+1) ) \
- do if(*a>v) { if( v - *(a-1) < *a - v ) --a; break; } while(++a<af); \
+ { \
+ if( v>=*r ) \
+ { \
+ do if(*a>v) { if( v - *(a-1) < *a - v ) --a; break; } \
+ while(++a<af); \
+ parsed=1; \
+ } \
+ } \
\
/* Decreasing array. */ \
else \
- do if(*a<v) { if( *(a-1) - v < v - *a ) --a; break; } while(++a<af); \
+ { \
+ if(v<=*r) \
+ { \
+ do if(*a<v) { if( *(a-1) - v < v - *a ) --a; break; } \
+ while(++a<af); \
+ parsed=1; \
+ } \
+ } \
\
- /* Set the difference. */ \
- if(a<af) index=a-r; \
+ /* Set the difference if the value is actually in the range. */ \
+ if(parsed && a<af) index = a-r; \
}
size_t
gal_statistics_quantile_function_index(gal_data_t *input, gal_data_t *value,
int inplace)
{
- size_t index=-1;
+ int parsed=0;
+ size_t index=GAL_BLANK_SIZE_T;
gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
/* A sanity check. */
@@ -386,16 +401,16 @@ gal_statistics_quantile_function_index(gal_data_t *input,
gal_data_t *value,
/* Find the result: */
switch(nbs->type)
{
- case GAL_TYPE_UINT8: STATS_QFUNC( uint8_t ); break;
- case GAL_TYPE_INT8: STATS_QFUNC( int8_t ); break;
- case GAL_TYPE_UINT16: STATS_QFUNC( uint16_t ); break;
- case GAL_TYPE_INT16: STATS_QFUNC( int16_t ); break;
- case GAL_TYPE_UINT32: STATS_QFUNC( uint32_t ); break;
- case GAL_TYPE_INT32: STATS_QFUNC( int32_t ); break;
- case GAL_TYPE_UINT64: STATS_QFUNC( uint64_t ); break;
- case GAL_TYPE_INT64: STATS_QFUNC( int64_t ); break;
- case GAL_TYPE_FLOAT32: STATS_QFUNC( float ); break;
- case GAL_TYPE_FLOAT64: STATS_QFUNC( double ); break;
+ case GAL_TYPE_UINT8: STATS_QFUNC_IND( uint8_t ); break;
+ case GAL_TYPE_INT8: STATS_QFUNC_IND( int8_t ); break;
+ case GAL_TYPE_UINT16: STATS_QFUNC_IND( uint16_t ); break;
+ case GAL_TYPE_INT16: STATS_QFUNC_IND( int16_t ); break;
+ case GAL_TYPE_UINT32: STATS_QFUNC_IND( uint32_t ); break;
+ case GAL_TYPE_INT32: STATS_QFUNC_IND( int32_t ); break;
+ case GAL_TYPE_UINT64: STATS_QFUNC_IND( uint64_t ); break;
+ case GAL_TYPE_INT64: STATS_QFUNC_IND( int64_t ); break;
+ case GAL_TYPE_FLOAT32: STATS_QFUNC_IND( float ); break;
+ case GAL_TYPE_FLOAT64: STATS_QFUNC_IND( double ); break;
default:
error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
__func__, nbs->type);
@@ -411,12 +426,24 @@ gal_statistics_quantile_function_index(gal_data_t *input,
gal_data_t *value,
/* Return the quantile function of the given value as float64. */
+#define STATS_QFUNC(IT) { \
+ IT *a=nbs->array, v=*((IT *)(value->array)); \
+ \
+ /* Increasing array: */ \
+ if( *a < *(a+1) ) \
+ d[0] = v<*a ? -INFINITY : INFINITY; \
+ \
+ /* Decreasing array. */ \
+ else \
+ d[0] = v>*a ? INFINITY : -INFINITY; \
+ }
gal_data_t *
gal_statistics_quantile_function(gal_data_t *input, gal_data_t *value,
int inplace)
{
double *d;
size_t dsize=1;
+ gal_data_t *nbs=gal_statistics_no_blank_sorted(input, inplace);
gal_data_t *out=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, &dsize,
NULL, 1, -1, NULL, NULL, NULL);
size_t ind=gal_statistics_quantile_function_index(input, value, inplace);
@@ -424,7 +451,32 @@ gal_statistics_quantile_function(gal_data_t *input,
gal_data_t *value,
/* Note that counting of the index starts from 0, so for the quantile we
should divided by (size - 1). */
d=out->array;
- d[0] = ( ind==-1 ? NAN : ((double)ind) / ((double)(input->size - 1)) );
+ if(ind==GAL_BLANK_SIZE_T)
+ {
+ /* See if the value is larger or smaller than the input's minimum or
+ maximum. */
+ switch(nbs->type)
+ {
+ case GAL_TYPE_UINT8: STATS_QFUNC( uint8_t ); break;
+ case GAL_TYPE_INT8: STATS_QFUNC( int8_t ); break;
+ case GAL_TYPE_UINT16: STATS_QFUNC( uint16_t ); break;
+ case GAL_TYPE_INT16: STATS_QFUNC( int16_t ); break;
+ case GAL_TYPE_UINT32: STATS_QFUNC( uint32_t ); break;
+ case GAL_TYPE_INT32: STATS_QFUNC( int32_t ); break;
+ case GAL_TYPE_UINT64: STATS_QFUNC( uint64_t ); break;
+ case GAL_TYPE_INT64: STATS_QFUNC( int64_t ); break;
+ case GAL_TYPE_FLOAT32: STATS_QFUNC( float ); break;
+ case GAL_TYPE_FLOAT64: STATS_QFUNC( double ); break;
+ default:
+ error(EXIT_FAILURE, 0, "%s: type code %d not recognized",
+ __func__, nbs->type);
+ }
+ }
+ else
+ d[0] = (double)ind / ((double)(nbs->size - 1));
+
+ /* Clean up and return. */
+ if(nbs!=input) gal_data_free(nbs);
return out;
}