[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master 3600a5e 3/6: Upper limit magnitude calculation
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master 3600a5e 3/6: Upper limit magnitude calculation in MakeCatalog |
|
Date: |
Sun, 2 Oct 2016 22:23:45 +0000 (UTC) |
branch: master
commit 3600a5eb4258d253fb227ca75c128b4d4dc4c143
Author: Mohammad Akhlaghi <address@hidden>
Commit: Mohammad Akhlaghi <address@hidden>
Upper limit magnitude calculation in MakeCatalog
The upper limit magnitude is a very important factor when doing photometry
on the fainter objects in an image (as explained more fully in the
book). With this commit, MakeCatalog can now produce the upper limit
magnitude for each object.
In MakeCatalog, the `upperlimit.c' and `upperlimit.h' source files were
added that are fully devoted to measuing the upper limit magnitude for each
object. The external function in `upperlimit.c' was written to not depend
on the MakeCatalog parameters structure. In this way, later it can easily
be moved with the libraries if necessary. Seven new parameters were added
to MakeCatalog for the job: `envseed', `upmask', `upmaskhdu', `upnum',
`upsclipmultip', `upsclipaccu', and `upnsigma'. In the `--help' output,
these are all classified in their own category to be more clear. Since it
involves random numbers, all the input parameters are also writen to the
comments of the output catalog for reproducibility.
In the book, the old "Depth and limiting magnitude" section was renamed to
"Quantifying data limits" and the explanations were updated to be more
consistent.
This finishes task #14167.
---
bin/mkcatalog/Makefile.am | 2 +-
bin/mkcatalog/args.h | 203 +++++++++++---
bin/mkcatalog/astmkcatalog.conf | 6 +
bin/mkcatalog/columns.c | 37 +++
bin/mkcatalog/columns.h | 3 +
bin/mkcatalog/main.h | 24 ++
bin/mkcatalog/mkcatalog.c | 47 ++++
bin/mkcatalog/ui.c | 120 +++++++-
bin/mkcatalog/upperlimit.c | 571 +++++++++++++++++++++++++++++++++++++++
bin/mkcatalog/upperlimit.h | 33 +++
doc/gnuastro.texi | 458 +++++++++++++++++++------------
tests/mkcatalog/simple.sh | 2 +-
12 files changed, 1289 insertions(+), 217 deletions(-)
diff --git a/bin/mkcatalog/Makefile.am b/bin/mkcatalog/Makefile.am
index 161a944..fcc7bcd 100644
--- a/bin/mkcatalog/Makefile.am
+++ b/bin/mkcatalog/Makefile.am
@@ -29,7 +29,7 @@ AM_CPPFLAGS = -I\$(top_srcdir)/bootstrapped/lib
bin_PROGRAMS = astmkcatalog
astmkcatalog_SOURCES = main.c main.h cite.h ui.c ui.h args.h \
-mkcatalog.c mkcatalog.h columns.c columns.h
+mkcatalog.c mkcatalog.h columns.c columns.h upperlimit.c upperlimit.h
astmkcatalog_LDADD = $(top_builddir)/bootstrapped/lib/libgnu.la \
-lgnuastro
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index 5778b85..0eb25d1 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -74,7 +74,7 @@ const char doc[] =
f g k l u v w
F G J L Q R U W X Y Z
- Number keys used: <=533
+ Number keys used: <= 541
Options with keys (second structure element) larger than 500 do not
have a short version.
@@ -193,6 +193,7 @@ static struct argp_option options[] =
+
{
0, 0, 0, 0,
"Output:",
@@ -251,16 +252,82 @@ static struct argp_option options[] =
{
0, 0, 0, 0,
- "Catalog columns:",
+ "Upper limit magnitude:",
+ 3
+ },
+ {
+ "upmask",
+ 535,
+ "STR",
+ 0,
+ "Mask image file name only for upper limit.",
+ 3
+ },
+ {
+ "upmaskhdu",
+ 536,
+ "STR",
+ 0,
+ "Upper limit mask image header name or number.",
+ 3
+ },
+ {
+ "upnum",
+ 537,
+ "INT",
+ 0,
+ "Number of randomly positioned samples.",
+ 3
+ },
+ {
+ "envseed",
+ 540,
+ 0,
+ 0,
+ "Random number generator info from environment.",
+ 3
+ },
+ {
+ "upsclipmultip",
+ 538,
+ "FLT",
+ 0,
+ "Sigma clipping multiple of std.",
+ 3
+ },
+ {
+ "upsclipaccu",
+ 539,
+ "FLT",
+ 0,
+ "Acceptable relative diff to stop clipping.",
3
},
{
+ "upnsigma",
+ 541,
+ "FLT",
+ 0,
+ "Multiple of final sigma to use.",
+ 3
+ },
+
+
+
+
+
+ {
+ 0, 0, 0, 0,
+ "Catalog columns:",
+ 4
+ },
+ {
"id",
'i',
0,
0,
"Overall ID of this object or clump.",
- 3
+ 4
},
{
"hostobjid",
@@ -268,7 +335,7 @@ static struct argp_option options[] =
0,
0,
"ID of object hosting this clump.",
- 3
+ 4
},
{
"idinhostobj",
@@ -276,7 +343,7 @@ static struct argp_option options[] =
0,
0,
"ID of clump in host object.",
- 3
+ 4
},
{
"numclumps",
@@ -284,7 +351,7 @@ static struct argp_option options[] =
0,
0,
"Number of clumps in this object.",
- 3
+ 4
},
{
"area",
@@ -292,7 +359,7 @@ static struct argp_option options[] =
0,
0,
"Number of pixels.",
- 3
+ 4
},
{
"clumpsarea",
@@ -300,7 +367,7 @@ static struct argp_option options[] =
0,
0,
"Area of clumps in an object.",
- 3
+ 4
},
{
"x",
@@ -308,7 +375,7 @@ static struct argp_option options[] =
0,
0,
"All obj. flux weighted center (first FITS axis).",
- 3
+ 4
},
{
"y",
@@ -316,7 +383,7 @@ static struct argp_option options[] =
0,
0,
"All obj. flux weighted center (second FITS axis).",
- 3
+ 4
},
{
"geox",
@@ -324,7 +391,7 @@ static struct argp_option options[] =
0,
0,
"All obj. geometric center (first FITS axis).",
- 3
+ 4
},
{
"geoy",
@@ -332,7 +399,7 @@ static struct argp_option options[] =
0,
0,
"All obj. geometric center (second FITS axis).",
- 3
+ 4
},
{
"clumpsx",
@@ -340,7 +407,7 @@ static struct argp_option options[] =
0,
0,
"Clumps flux weighted center (first FITS axis).",
- 3
+ 4
},
{
"clumpsy",
@@ -348,7 +415,7 @@ static struct argp_option options[] =
0,
0,
"Clumps flux weighted center (second FITS axis).",
- 3
+ 4
},
{
"clumpsgeox",
@@ -356,7 +423,7 @@ static struct argp_option options[] =
0,
0,
"Clumps geometric center (first FITS axis).",
- 3
+ 4
},
{
"clumpsgeoy",
@@ -364,7 +431,7 @@ static struct argp_option options[] =
0,
0,
"Clumps geometric center (second FITS axis).",
- 3
+ 4
},
{
"ra",
@@ -372,7 +439,7 @@ static struct argp_option options[] =
0,
0,
"All object flux weighted center right ascension.",
- 3
+ 4
},
{
"dec",
@@ -380,7 +447,7 @@ static struct argp_option options[] =
0,
0,
"All object flux weighted center declination.",
- 3
+ 4
},
{
"geora",
@@ -388,7 +455,7 @@ static struct argp_option options[] =
0,
0,
"All object geometric center right ascension.",
- 3
+ 4
},
{
"geodec",
@@ -396,7 +463,7 @@ static struct argp_option options[] =
0,
0,
"All object geometric center declination.",
- 3
+ 4
},
{
"clumpsra",
@@ -404,7 +471,7 @@ static struct argp_option options[] =
0,
0,
"Clumps flux weighted center right ascension.",
- 3
+ 4
},
{
"clumpsdec",
@@ -412,7 +479,7 @@ static struct argp_option options[] =
0,
0,
"Clumps flux weighted center declination.",
- 3
+ 4
},
{
"clumpsgeora",
@@ -420,7 +487,7 @@ static struct argp_option options[] =
0,
0,
"Clumps geometric center right ascension.",
- 3
+ 4
},
{
"clumpsgeodec",
@@ -428,7 +495,7 @@ static struct argp_option options[] =
0,
0,
"Clumps geometric center declination.",
- 3
+ 4
},
{
"brightness",
@@ -436,7 +503,7 @@ static struct argp_option options[] =
0,
0,
"Brightness (sum of pixel values).",
- 3
+ 4
},
{
"clumpbrightness",
@@ -444,7 +511,7 @@ static struct argp_option options[] =
0,
0,
"Brightness in clumps of an object.",
- 3
+ 4
},
{
"noriverbrightness",
@@ -452,7 +519,7 @@ static struct argp_option options[] =
0,
0,
"Sky (not river) subtracted clump brightness.",
- 3
+ 4
},
{
"magnitude",
@@ -460,7 +527,7 @@ static struct argp_option options[] =
0,
0,
"Total magnitude.",
- 3
+ 4
},
{
"magnitudeerr",
@@ -468,7 +535,7 @@ static struct argp_option options[] =
0,
0,
"Total magnitude error.",
- 3
+ 4
},
{
"clumpsmagnitude",
@@ -476,7 +543,15 @@ static struct argp_option options[] =
0,
0,
"Total magnitude of clumps in this object.",
- 3
+ 4
+ },
+ {
+ "upperlimitmag",
+ 534,
+ 0,
+ 0,
+ "Upper limit magnitude for each clump or object.",
+ 4
},
{
"riverave",
@@ -484,7 +559,7 @@ static struct argp_option options[] =
0,
0,
"Average river value surrounding this clump.",
- 3
+ 4
},
{
"rivernum",
@@ -492,7 +567,7 @@ static struct argp_option options[] =
0,
0,
"Number of river pixels surrounding this clump.",
- 3
+ 4
},
{
"sn",
@@ -500,7 +575,7 @@ static struct argp_option options[] =
0,
0,
"Signal to noise ratio column.",
- 3
+ 4
},
{
"sky",
@@ -508,7 +583,7 @@ static struct argp_option options[] =
0,
0,
"Sky value.",
- 3
+ 4
},
{
"std",
@@ -516,7 +591,7 @@ static struct argp_option options[] =
0,
0,
"Sky standard deviation.",
- 3
+ 4
},
{
"semimajor",
@@ -524,7 +599,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Semi-major axis.",
- 3
+ 4
},
{
"semiminor",
@@ -532,7 +607,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Semi-minor axis.",
- 3
+ 4
},
{
"positionangle",
@@ -540,7 +615,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Position angle.",
- 3
+ 4
},
{
"geosemimajor",
@@ -548,7 +623,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Semi-major axis.",
- 3
+ 4
},
{
"geosemiminor",
@@ -556,7 +631,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Semi-minor axis.",
- 3
+ 4
},
{
"geopositionangle",
@@ -564,7 +639,7 @@ static struct argp_option options[] =
0,
0,
"Flux weighted Position angle.",
- 3
+ 4
},
@@ -613,7 +688,8 @@ parse_opt(int key, char *arg, struct argp_state *state)
/* Input: */
case 'M':
- gal_checkset_allocate_copy_set(arg, &p->up.maskname, &p->up.masknameset);
+ gal_checkset_allocate_copy_set(arg, &p->up.maskname,
+ &p->up.masknameset);
break;
case 'H':
gal_checkset_allocate_copy_set(arg, &p->up.mhdu, &p->up.mhduset);
@@ -630,7 +706,8 @@ parse_opt(int key, char *arg, struct argp_state *state)
&p->up.clumplabsnameset);
break;
case 502:
- gal_checkset_allocate_copy_set(arg, &p->up.clumphdu, &p->up.clumphduset);
+ gal_checkset_allocate_copy_set(arg, &p->up.clumphdu,
+ &p->up.clumphduset);
break;
case 's':
gal_checkset_allocate_copy_set(arg, &p->up.skyname, &p->up.skynameset);
@@ -693,6 +770,42 @@ parse_opt(int key, char *arg, struct argp_state *state)
break;
+ /* Upper limit magnitude */
+ case 535:
+ gal_checkset_allocate_copy_set(arg, &p->up.upmaskname,
+ &p->up.upmasknameset);
+ break;
+ case 536:
+ gal_checkset_allocate_copy_set(arg, &p->up.upmaskhdu,
+ &p->up.upmaskhduset);
+ break;
+ case 537:
+ gal_checkset_sizet_l_zero(arg, &p->upnum, "upnum", key, SPACK,
+ NULL, 0);
+ p->up.upnumset=1;
+ break;
+ case 540:
+ p->envseed=1;
+ p->up.envseedset=1;
+ break;
+ case 538:
+ gal_checkset_float_l_0(arg, &p->upsclipmultip, "upsclipmultip",
+ key, SPACK, NULL, 0);
+ p->up.upsclipmultipset=1;
+ break;
+ case 539:
+ gal_checkset_float_l_0_s_1(arg, &p->upsclipaccu, "upsclipaccu",
+ key, SPACK, NULL, 0);
+ p->up.upsclipaccuset=1;
+ break;
+ case 541:
+ gal_checkset_float_l_0(arg, &p->upnsigma, "upnsigma", key,
+ SPACK, NULL, 0);
+ p->up.upnsigmaset=1;
+ break;
+
+
+
/* Catalog columns: */
case 'i':
gal_linkedlist_add_to_sll(&p->allcolsll, CATID);
@@ -806,6 +919,10 @@ parse_opt(int key, char *arg, struct argp_state *state)
gal_linkedlist_add_to_sll(&p->allcolsll, CATCLUMPSMAGNITUDE);
p->up.clumpsmagnitudeset=1;
break;
+ case 534:
+ gal_linkedlist_add_to_sll(&p->allcolsll, CATUPPERLIMITMAG);
+ p->up.upperlimitmagset=1;
+ break;
case 514:
gal_linkedlist_add_to_sll(&p->allcolsll, CATRIVERAVE);
p->up.riveraveset=1;
diff --git a/bin/mkcatalog/astmkcatalog.conf b/bin/mkcatalog/astmkcatalog.conf
index 061bdad..c552bcc 100644
--- a/bin/mkcatalog/astmkcatalog.conf
+++ b/bin/mkcatalog/astmkcatalog.conf
@@ -34,6 +34,12 @@
floatprecision 3
accuprecision 8
+# Upper limit magnitude:
+ upnum 200
+ upsclipmultip 3
+ upsclipaccu 0.2
+ upnsigma 1
+
# Catalog columns:
# sn 1
# std 1
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index ed41bc5..0febae3 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -43,6 +43,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include "columns.h"
#include "mkcatalog.h"
+#include "upperlimit.h"
@@ -889,6 +890,42 @@ brightnessmag(struct mkcatalogparams *p, size_t col, char
*target,
void
+upperlimitcol(struct mkcatalogparams *p)
+{
+ size_t i;
+ float *std;
+ double *ptr;
+
+ /* For the comments: */
+ p->unitp = CATUNITMAG;
+ sprintf(p->description, "%lu: Upper limit magnitude for this %s.",
+ p->curcol, p->name);
+
+
+ /* Correct the raw values (divide them by area) if not already
+ done. */
+ std=upperlimit(p->img, p->sky, p->objects, p->upmask, p->s0, p->s1,
+ p->upnum, p->cp.numthreads, p->envseed, p->upsclipmultip,
+ p->upsclipaccu);
+
+
+ /* Write the standard deviations values in the final catalog as
+ magnitudes. */
+ for(i=0;i<p->num;++i)
+ {
+ /* `ptr' is defined for a short/readable line. */
+ ptr = &p->cat[i * p->numcols + p->curcol ];
+ *ptr = -2.5f * log10( p->upnsigma*std[i+1] ) + p->zeropoint;
+ }
+
+ free(std);
+}
+
+
+
+
+
+void
skystd(struct mkcatalogparams *p, size_t col)
{
size_t i;
diff --git a/bin/mkcatalog/columns.h b/bin/mkcatalog/columns.h
index 91cce4b..2f4a5d3 100644
--- a/bin/mkcatalog/columns.h
+++ b/bin/mkcatalog/columns.h
@@ -73,6 +73,9 @@ brightnessmag(struct mkcatalogparams *p, size_t col, char
*target,
char *scale);
void
+upperlimitcol(struct mkcatalogparams *p);
+
+void
skystd(struct mkcatalogparams *p, size_t col);
void
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index 483e207..5ce8f30 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -23,6 +23,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#ifndef MAIN_H
#define MAIN_H
+#include <gsl/gsl_rng.h>
#include <gnuastro/fits.h>
#include <gnuastro/commonparams.h>
@@ -86,6 +87,7 @@ enum objectcols
OFlxWhtDec, /* Dec of (OFlxWhtX, OFlxWhtY). */
OAREAC, /* Area of clumps in this object. */
OBrightnessC, /* Brightness in object clumps. */
+ OUpperBright, /* Upper limit brightness for this object. */
OFlxWhtCX, /* Sum of (flux-sky)*x on object clumps. */
OFlxWhtCY, /* Sum of (flux-sky)*y on obj. clumps. */
OFlxWhtCXX, /* Sum of (flux-sky)*x*x on object clumps. */
@@ -134,6 +136,7 @@ enum clumpcols
CFlxWhtDec, /* Dec of (CFlxWhtX, CFlxWhtY). */
CBrightness, /* River subtracted brightness. */
CNoRiverBrightness, /* Sky (not river) subtracted brightness. */
+ CUpperBright, /* Upper limit brightness for this clump. */
CRivAve, /* Average value in rivers around clump. */
CRivArea, /* Num river pixels around this clump. */
CSKY, /* Sum of sky value on this object. */
@@ -192,6 +195,7 @@ enum outcols
CATMAGNITUDE,
CATMAGNITUDEERR,
CATCLUMPSMAGNITUDE,
+ CATUPPERLIMITMAG,
CATRIVERAVE,
CATRIVERNUM,
CATSN,
@@ -222,6 +226,8 @@ struct uiparams
char *skyhdu; /* Sky HDU name. */
char *stdname; /* Sky STD value image file name. */
char *stdhdu; /* Sky STD HDU name. */
+ char *upmaskname; /* Name of mask for upper magnitude. */
+ char *upmaskhdu; /* HDU Name of mask for upper mag. */
int masknameset;
int mhduset;
@@ -233,6 +239,8 @@ struct uiparams
int skyhduset;
int stdnameset;
int stdhduset;
+ int envseedset;
+
int zeropointset;
int skysubtractedset;
int thresholdset;
@@ -244,6 +252,13 @@ struct uiparams
int accuprecisionset;
int nsigmagset;
+ int upmasknameset;
+ int upmaskhduset;
+ int upnumset;
+ int upsclipmultipset;
+ int upsclipaccuset;
+ int upnsigmaset;
+
int idset;
int hostobjidset;
int idinhostobjset;
@@ -272,6 +287,7 @@ struct uiparams
int magnitudeset;
int magnitudeerrset;
int clumpsmagnitudeset;
+ int upperlimitmagset;
int riveraveset;
int rivernumset;
int snset;
@@ -309,6 +325,7 @@ struct mkcatalogparams
int skysubtracted; /* Input is already sky subtracted. */
double nsigmag; /* Multiple of Sky STD to report mag. */
double threshold; /* Only pixels larger than this *STD. */
+ int envseed; /* Use GSL_RNG_SEED for random seed. */
/* Output: */
char *ocatname; /* File name of object catalog. */
@@ -319,6 +336,13 @@ struct mkcatalogparams
int floatprecision; /* Floating point precision. */
int accuprecision; /* Accurate floating point precision. */
+ /* Upper limit: */
+ long *upmask; /* Mask array, only for upper liimit. */
+ size_t upnum; /* Number of random samples. */
+ float upsclipmultip; /* Sigma clip multiple for std calc. */
+ float upsclipaccu; /* Sigma clip accuracy for std calc. */
+ float upnsigma; /* Multiple of sigma for upper-limit. */
+
/* Operating mode: */
/* Internal: */
diff --git a/bin/mkcatalog/mkcatalog.c b/bin/mkcatalog/mkcatalog.c
index 67a419b..88b63eb 100644
--- a/bin/mkcatalog/mkcatalog.c
+++ b/bin/mkcatalog/mkcatalog.c
@@ -476,6 +476,14 @@ makeoutput(struct mkcatalogparams *p)
sprintf(p->line, "# Sky STD %s (hdu: %s)\n", p->up.stdname,
p->up.stdhdu);
strcat(comment, p->line);
+ if(p->obj0clump1==0 && p->upmask)
+ {
+ sprintf(p->line, "# Upper limit magnitude mask %s (hdu: %s)\n",
+ p->up.upmaskname, p->up.upmaskhdu);
+ strcat(comment, p->line);
+ }
+
+
/* If a magnitude is also desired, print the zero point
@@ -523,6 +531,41 @@ makeoutput(struct mkcatalogparams *p)
strcat(comment, p->line);
}
+ /* If an upper limit was output then report its parameters: */
+ if(p->obj0clump1==0 && p->up.upperlimitmagset)
+ {
+ if(p->upmask)
+ {
+ sprintf(p->line, "# Upper limit magnitude mask %s (hdu: %s)\n",
+ p->up.upmaskname, p->up.upmaskhdu);
+ strcat(comment, p->line);
+ }
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%lu\n", "Number of upper "
+ "limit magnitude samples", p->upnum);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%lu\n", "Number of threads "
+ "used for upper limit magnitude",
+ p->cp.numthreads);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%s\n", "Random number "
+ "generator type for upper limit magnitude",
+ gsl_rng_default->name);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%lu\n", "Random number "
+ "generator seed for upper limit magnitude",
+ gsl_rng_default_seed);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%.3f\n", "STD multiple for "
+ "upper limit magnitude sigma-clip", p->upsclipmultip);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%.3f\n", "STD accuracy "
+ "to stop upper limit magnitude sigma-clip", p->upsclipaccu);
+ strcat(comment, p->line);
+ sprintf(p->line, "# "CATDESCRIPTLENGTH"%.3f\n", "Multiple of "
+ "sigma for final upper limit magnitude", p->upnsigma);
+ strcat(comment, p->line);
+ }
+
/* Prepare for printing the columns: */
strcat(comment, "#\n# Columns:\n# --------\n");
@@ -660,6 +703,10 @@ makeoutput(struct mkcatalogparams *p)
brightnessmag(p, OBrightnessC, MKCATCINO, MKCATMAG);
break;
+ case CATUPPERLIMITMAG:
+ upperlimitcol(p);
+ break;
+
case CATRIVERAVE:
brightnessmag(p, CRivAve, MKRIVERSSUR, MKCATBRIGHT);
break;
diff --git a/bin/mkcatalog/ui.c b/bin/mkcatalog/ui.c
index 881561e..2cd9e1b 100644
--- a/bin/mkcatalog/ui.c
+++ b/bin/mkcatalog/ui.c
@@ -157,6 +157,7 @@ readconfig(char *filename, struct mkcatalogparams *p)
}
+
/* Outputs */
else if(strcmp(name, "output")==0)
gal_checkset_allocate_copy_set(value, &cp->output,
@@ -205,6 +206,51 @@ readconfig(char *filename, struct mkcatalogparams *p)
}
+ /* Upper limit magnitude */
+ else if (strcmp(name, "upmaskname")==0)
+ gal_checkset_allocate_copy_set(value, &up->upmaskname,
+ &up->upmasknameset);
+ else if (strcmp(name, "upmaskhdu")==0)
+ gal_checkset_allocate_copy_set(value, &up->upmaskhdu,
+ &up->upmaskhduset);
+ else if(strcmp(name, "upnum")==0)
+ {
+ if(up->upnumset) continue;
+ gal_checkset_sizet_l_zero(value, &p->upnum, name,
+ key, SPACK, filename, lineno);
+ up->upnumset=1;
+ }
+ else if(strcmp(name, "envseed")==0)
+ {
+ if(up->envseedset) continue;
+ gal_checkset_int_zero_or_one(value, &p->envseed, name, key,
+ SPACK, filename, lineno);
+ up->envseedset=1;
+ }
+ else if(strcmp(name, "upsclipmultip")==0)
+ {
+ if(up->upsclipmultipset) continue;
+ gal_checkset_float_l_0(value, &p->upsclipmultip, name,
+ key, SPACK, filename, lineno);
+ up->upsclipmultipset=1;
+ }
+ else if(strcmp(name, "upsclipaccu")==0)
+ {
+ if(up->upsclipaccuset) continue;
+ gal_checkset_float_l_0_s_1(value, &p->upsclipaccu, name,
+ key, SPACK, filename, lineno);
+ up->upsclipaccuset=1;
+ }
+ else if(strcmp(name, "upnsigma")==0)
+ {
+ if(up->upnsigmaset) continue;
+ gal_checkset_float_l_0(value, &p->upnsigma, name,
+ key, SPACK, filename, lineno);
+ up->upnsigmaset=1;
+ }
+
+
+
/* Catalog columns */
else if(strcmp(name, "id")==0)
{
@@ -458,6 +504,15 @@ readconfig(char *filename, struct mkcatalogparams *p)
gal_linkedlist_add_to_sll(&p->allcolsll, CATCLUMPSMAGNITUDE);
up->clumpsmagnitudeset=1;
}
+ else if(strcmp(name, "upperlimitmag")==0)
+ {
+ if(up->upperlimitmagset) continue;
+ gal_checkset_int_zero_or_one(value, &yes, name, key, SPACK,
+ filename, lineno);
+ if(!yes) continue;
+ gal_linkedlist_add_to_sll(&p->allcolsll, CATUPPERLIMITMAG);
+ up->upperlimitmagset=1;
+ }
else if(strcmp(name, "riverave")==0)
{
if(up->riveraveset) continue;
@@ -640,6 +695,24 @@ printvalues(FILE *fp, struct mkcatalogparams *p)
if(up->accuprecisionset)
fprintf(fp, CONF_SHOWFMT"%d\n", "accuprecision", p->accuprecision);
+ /* Upper limit magnitude */
+ fprintf(fp, "\n# Upper limit magnitude:\n");
+ if(up->upmasknameset)
+ fprintf(fp, CONF_SHOWFMT"%s\n", "upmaskname", up->upmaskname);
+ if(up->upmaskhduset)
+ fprintf(fp, CONF_SHOWFMT"%s\n", "upmaskhdu", up->upmaskhdu);
+ if(up->upnumset)
+ fprintf(fp, CONF_SHOWFMT"%lu\n", "upnum", p->upnum);
+ if(up->envseedset)
+ fprintf(fp, CONF_SHOWFMT"%d\n", "envseed", p->envseed);
+ if(up->upsclipmultipset)
+ fprintf(fp, CONF_SHOWFMT"%.3f\n", "upsclipmultip", p->upsclipmultip);
+ if(up->upsclipaccuset)
+ fprintf(fp, CONF_SHOWFMT"%.3f\n", "upsclipaccu", p->upsclipaccu);
+ if(up->upnsigmaset)
+ fprintf(fp, CONF_SHOWFMT"%.3f\n", "upnsigma", p->upnsigma);
+
+
/* Catalog columns, since order is important. Notice that they have
to be printed in opposite order (because of the way they are read
through a simple linked list). */
@@ -719,6 +792,9 @@ printvalues(FILE *fp, struct mkcatalogparams *p)
case CATCLUMPSBRIGHTNESS:
fprintf(fp, CONF_SHOWFMT"%d\n", "clumpsbrightness", 1);
break;
+ case CATUPPERLIMITMAG:
+ fprintf(fp, CONF_SHOWFMT"%d\n", "upperlimitmag", 1);
+ break;
case CATNORIVERBRIGHTNESS:
fprintf(fp, CONF_SHOWFMT"%d\n", "noriverbrightness", 1);
break;
@@ -821,6 +897,19 @@ checkifset(struct mkcatalogparams *p)
if(up->accuprecisionset==0)
GAL_CONFIGFILES_REPORT_NOTSET("accuprecision");
+ /* Upper limit magnitude. */
+ if(p->up.upperlimitmagset)
+ {
+ if(p->up.upnumset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("upnum");
+ if(p->up.upsclipmultipset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("upsclipmultip");
+ if(p->up.upsclipaccuset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("upsclipaccu");
+ if(p->up.upnsigmaset==0)
+ GAL_CONFIGFILES_REPORT_NOTSET("upnsigma");
+ }
+
GAL_CONFIGFILES_END_OF_NOTSET_REPORT;
}
@@ -1199,6 +1288,9 @@ preparearrays(struct mkcatalogparams *p)
case CATCLUMPSBRIGHTNESS:
p->objcols[p->objncols++] = p->allcols[i];
break;
+ case CATUPPERLIMITMAG:
+ p->objcols[p->objncols++] = p->allcols[i];
+ break;
case CATNORIVERBRIGHTNESS:
p->clumpcols[p->clumpncols++] = p->allcols[i];
break;
@@ -1285,8 +1377,11 @@ preparearrays(struct mkcatalogparams *p)
checksetlong(p, p->up.objlabsname, p->up.objhdu, &p->objects);
if(p->up.clumplabsname)
checksetlong(p, p->up.clumplabsname, p->up.clumphdu, &p->clumps);
- else
- p->clumps=NULL;
+ else p->clumps=NULL;
+ if(p->up.upmasknameset)
+ checksetlong(p, p->up.upmaskname, p->up.upmaskhdu, &p->upmask);
+ else p->upmask=NULL;
+
/* Read the necessary keywords. */
readkeywords(p);
@@ -1336,6 +1431,14 @@ preparearrays(struct mkcatalogparams *p)
for(i=1;i<=p->numclumps;++i)
p->cinfo[i*CCOLUMNS+CPOSSHIFTX]=p->cinfo[i*CCOLUMNS+CPOSSHIFTY]=NAN;
+ /* Allocate the random number generator: */
+ if(p->up.upperlimitmagset)
+ {
+ gsl_rng_env_setup();
+ if(p->envseed==0)
+ gsl_rng_default_seed=gal_timing_time_based_rng_seed();
+ }
+
/* Clean up: */
gal_linkedlist_free_sll(p->allcolsll);
}
@@ -1380,7 +1483,6 @@ setparams(int argc, char *argv[], struct mkcatalogparams
*p)
if(argp_parse(&thisargp, argc, argv, 0, 0, p))
error(EXIT_FAILURE, errno, "parsing arguments");
-
/* Add the user default values and save them if asked. */
GAL_CONFIGFILES_CHECK_SET_CONFIG;
@@ -1394,7 +1496,6 @@ setparams(int argc, char *argv[], struct mkcatalogparams
*p)
if(p->up.inputname)
sanitycheck(p);
-
/* Make the array of input images. */
preparearrays(p);
@@ -1416,6 +1517,16 @@ setparams(int argc, char *argv[], struct mkcatalogparams
*p)
p->up.clumphdu);
printf(" - Sky %s (hdu: %s)\n", p->up.skyname, p->up.skyhdu);
printf(" - Sky STD %s (hdu: %s)\n", p->up.stdname, p->up.stdhdu);
+ if(p->up.upperlimitmagset)
+ {
+ if(p->up.upmasknameset)
+ printf(" - Upper limit mask: %s (hdu: %s)\n", p->up.inputname,
+ p->cp.hdu);
+ printf(" - RNG type for upper-limit magnitude: %s\n",
+ gsl_rng_default->name);
+ printf(" - RNG seed for upper-limit magnitude: %ld\n",
+ gsl_rng_default_seed);
+ }
}
}
@@ -1462,6 +1573,7 @@ freeandreport(struct mkcatalogparams *p, struct timeval
*t1)
free(p->up.skyhdu);
free(p->up.stdhdu);
free(p->up.clumphdu);
+ if(p->upmask) free(p->upmask);
if(p->up.mhduset) free(p->up.mhdu);
if(p->wcs) wcsvfree(&p->nwcs, &p->wcs);
if(p->up.skynameset) free(p->up.skyname);
diff --git a/bin/mkcatalog/upperlimit.c b/bin/mkcatalog/upperlimit.c
new file mode 100644
index 0000000..886f577
--- /dev/null
+++ b/bin/mkcatalog/upperlimit.c
@@ -0,0 +1,571 @@
+/*********************************************************************
+MakeCatalog - Make a catalog from an input and labeled image.
+MakeCatalog 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 <math.h>
+#include <stdio.h>
+#include <errno.h>
+#include <error.h>
+#include <stdlib.h>
+#include <sys/time.h>
+
+#include <gsl/gsl_rng.h>
+#include <gnuastro/fits.h>
+#include <gnuastro/timing.h>
+#include <gnuastro/threads.h>
+#include <gnuastro/txtarray.h>
+#include <gnuastro/statistics.h>
+#include <gsl/gsl_statistics_double.h>
+
+#include "upperlimit.h"
+
+
+
+
+/*********************************************************************/
+/************ Structures and macros ****************/
+/*********************************************************************/
+/* These structures and macros are defined in the C file since they are
+ only necessary in the processing here, they do not need to be given to
+ the user through the header, it will just complicate their
+ name-space. */
+struct upperlimitparams
+{
+ size_t s0; /* Number of rows in image. */
+ size_t s1; /* Number of columns in image. */
+ int envseed; /* ==1: evironment for rng setup. */
+ long minlab; /* Minimum label in image. */
+ long numlabs; /* Number of labels in image. */
+ float *img; /* The input image. */
+ size_t *box; /* The box parameters of each label. */
+ unsigned char **pix; /* The pixels of each label. */
+ float *std; /* Standard deviation of the sums. */
+ size_t upnum; /* Number of random samples. */
+ size_t numthreads; /* Number of threads to use. */
+ float sclipmultip; /* Multiple of STD for sigma clip. */
+ float sclipaccu; /* Accuracy to stop sigma clipping. */
+};
+
+
+
+
+
+/* Information for each thread. */
+struct tupperlimitparams
+{
+ size_t id; /* Id of this thread. */
+ size_t *indexs; /* Indexes for this thread. */
+ pthread_barrier_t *b; /* Barrier for all threads. */
+ struct upperlimitparams *p; /* Pointer to main program struct. */
+};
+
+
+
+
+
+/* Columns of labeled region information. Note that the width columns
+ (WIDCOL) are actually the same as the max columns (MAXCOL). The
+ width values are written over the maximum columns because we don't
+ need the maximum values any more. */
+#define XMINCOL 0
+#define YMINCOL 1
+#define XMAXCOL 2
+#define YMAXCOL 3
+#define XWIDCOL XMAXCOL
+#define YWIDCOL YMAXCOL
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/************ Preparations ****************/
+/*********************************************************************/
+static void
+fillseginfo(struct upperlimitparams *p, long *seg)
+{
+ float *f, *ff;
+ long *l, *ll, maxlab;
+ unsigned char *u, *uu;
+ size_t xw, yw, xmin, ymin;
+ size_t i, j, asize, label, size=p->s0*p->s1;
+
+
+ /* Find the minimum and maximum labels: */
+ maxlab = INT32_MIN;
+ p->minlab = INT32_MAX;
+ ll=(l=seg)+size;
+ do
+ /* We don't want to look at pixels with value zero or blank pixels, so
+ ignore them. */
+ if( *l!=0 && *l!=GAL_FITS_LONG_BLANK )
+ {
+ if(*l>maxlab) maxlab=*l;
+ if(*l<p->minlab) p->minlab=*l;
+ }
+ while(++l<ll);
+
+
+ /* Do a small sanity check: */
+ if(maxlab<0 || p->minlab<0)
+ error(EXIT_FAILURE, 0, "the labeled image must not contain negative "
+ "pixels");
+
+
+ /* Allocate the space to keep the object information. This array will
+ keep the minimum and maximum coordinate for each labeled region. Note
+ that max is in the image, and we want to use the labels as indexs (row
+ numbers), so we must allocate max+1 rows.
+
+ Note: we are not assuming that the first index is 1. Because many
+ cases can arise where the lowest index can be a large number (for
+ example when the user can be working on a sub-set of objects selected
+ from a large catalog), so we don't want to allocate a huge amount of
+ memory and waste thread resources on them. We have stored the minimum
+ label and will use that in reporting the final output to make a full
+ array. */
+ p->numlabs = maxlab - p->minlab + 1;
+ errno=0;
+ p->box=malloc(p->numlabs * 4 * sizeof *p->box);
+ if(p->box==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for p->box",
+ p->numlabs * 4 * sizeof *p->box);
+
+ /* Initialize the information array*/
+ for(i=0;i<p->numlabs;++i)
+ {
+ p->box[ i*4 + XMINCOL ] = p->s0;
+ p->box[ i*4 + YMINCOL ] = p->s1;
+ p->box[ i*4 + XMAXCOL ] = p->box[ i*4 + YMAXCOL ] = 0;
+ }
+
+ /* Fill it in with the minimum and maximum positions. */
+ for(i=0;i<p->s0;++i)
+ for(j=0;j<p->s1;++j)
+ if(seg[i*p->s1+j]>0)
+ {
+ label = seg[i*p->s1+j] - p->minlab;
+ if(i<p->box[ label*4 + XMINCOL ]) p->box[ label*4 + XMINCOL ] = i;
+ if(j<p->box[ label*4 + YMINCOL ]) p->box[ label*4 + YMINCOL ] = j;
+ if(i>p->box[ label*4 + XMAXCOL ]) p->box[ label*4 + XMAXCOL ] = i;
+ if(j>p->box[ label*4 + YMAXCOL ]) p->box[ label*4 + YMAXCOL ] = j;
+ }
+
+ /* For a check. Note that in DS9 the coordinates are inverted and
+ start from 1, not zero.
+ {
+ i=0;
+ printf("%lu: (%lu, %lu) --> (%lu, %lu)\n", i+p->minlab,
+ p->box[i*4+YMINCOL]+1, p->box[i*4+XMINCOL]+1,
+ p->box[i*4+YMAXCOL]+1, p->box[i*4+XMAXCOL]+1);
+ }
+ */
+
+
+ /* Allocate the pointers to the object's pixels and allocate space
+ for the area of each detection. In the meantime, Change the xmax
+ and ymax to width in X and width in Y*/
+ errno=0;
+ p->pix=malloc(p->numlabs * sizeof *p->pix);
+ if(p->pix==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for p->pix",
+ p->numlabs * sizeof *p->pix);
+ for(i=0; i<p->numlabs; ++i)
+ {
+ /* Replace the max values with the width of each box. */
+ p->box[i*4+XWIDCOL] = p->box[i*4+XMAXCOL] - p->box[i*4+XMINCOL] + 1;
+ p->box[i*4+YWIDCOL] = p->box[i*4+YMAXCOL] - p->box[i*4+YMINCOL] + 1;
+
+ /* Allocate space for the pixels of this label. */
+ errno=0;
+ asize = p->box[i*4+XWIDCOL] * p->box[i*4+YWIDCOL] * sizeof **p->pix;
+ p->pix[i] = malloc(asize);
+ if(p->pix[i]==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for p->pix[i] (i=%lu)",
+ asize, i);
+
+ /* Set the array values. Any pixel that is equal to the object's
+ label over the box is given a value of 1.0f, other pixels are
+ given a value of 0.0f. We set the starting pointers for the
+ larger and smaller arrays in this row, then go along the row
+ until it finshes. `j' will take us to the next row after this
+ row finishes. */
+ xmin = p->box[ i*4+XMINCOL ];
+ ymin = p->box[ i*4+YMINCOL ];
+ xw = p->box[ i*4+XWIDCOL ];
+ yw = p->box[ i*4+YWIDCOL ];
+ for(j=0; j<xw; ++j)
+ {
+ uu = ( u = p->pix[i] + yw*j ) + yw;
+ l = seg + (xmin+j)*p->s1 + ymin;
+ do *u = (*l++ == i+p->minlab); while(++u<uu);
+ }
+ }
+
+
+ /* For a check:
+ i=0;
+ gal_fits_array_to_file("tmpf.fits", "check", BYTE_IMG, p->pix[i],
+ p->box[i*4+XWIDCOL], p->box[i*4+YWIDCOL], 0,
+ NULL, NULL, "myprog");
+ exit(0);
+ */
+
+ /* Allocate an array for the output values, note that this is for the
+ full list of IDs, not just those in the image. */
+ p->std = malloc( (maxlab+1) * sizeof *p->std );
+ if(p->std==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for p->std (upperlimit.c)",
+ (maxlab+1) * sizeof *p->std);
+ ff=(f=p->std)+(maxlab+1); do *f++=NAN; while(f<ff);
+}
+
+
+
+
+
+/* Subtract the Sky and set all the segmentation and mask pixels to
+ NaN. Note that the mask and sky images are optional. If the caller
+ didn't want them, they have to be set to NULL.
+
+ We could also be doing this during the actual flux calculation, but
+ since many iterations of many objects will be done, it is more efficient
+ to do all the preparations (setting blank pixels and subtracting the
+ sky. */
+void
+prepareimg(struct upperlimitparams *p, float *img, float *sky, long *seg,
+ long *mask)
+{
+ float *f, *ff;
+ size_t size=p->s0*p->s1;
+
+ /* Allocate the space for the image. */
+ errno=0;
+ p->img=malloc(size * sizeof *p->img);
+ if(p->img==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for p->img (upperlimit.c)",
+ size * sizeof *p->img);
+
+ /* Set the values. */
+ ff=(f=p->img)+size;
+ do
+ {
+ /* Note that the Sky and mask are optional (can be NULL), so we are
+ checking their existance (not NULL) first. */
+ if(mask)
+ {
+ *f = *img++ - (( *seg || *mask ) ? NAN : ( sky ? *sky : 0 ));
+ ++mask;
+ }
+ else
+ *f = *img++ - (*seg ? NAN : ( sky ? *sky : 0 ));
+
+ /* Increment the pointers (when they exist) to the next pixel. */
+ ++seg;
+ if(sky) ++sky;
+ }
+ while(++f<ff);
+
+ /* For a check:
+ gal_fits_array_to_file("in.fits", "check", FLOAT_IMG, p->img,
+ p->s0, p->s1, 0, NULL, NULL, "myprog");
+ exit(0);
+ */
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/************ Actual operation on threads ****************/
+/*********************************************************************/
+static void *
+upperlimit_on_thread(void *inparam)
+{
+ /* The first thing to do is to say what the input pointer actually is. */
+ struct tupperlimitparams *tp = (struct tupperlimitparams *)inparam;
+ struct upperlimitparams *p = tp->p;
+
+ /* Now you can go onto do defining the function like any other
+ function: first you define the variables and so on... */
+ double *sum;
+ gsl_rng *rng;
+ unsigned char *u, *uu;
+ float *l, ave, med, *fsum;
+ size_t i, j, c, xw, yw, lab, xmin, ymin;
+
+ /* Allcate space to keep the total sums */
+ errno=0;
+ sum=malloc(p->upnum*sizeof *sum);
+ if(sum==NULL)
+ error(EXIT_FAILURE, errno, "%lu bytes for sum (upperlimit.c)",
+ p->upnum*sizeof *sum);
+
+ /* Allocate the random number generator for this thread. Note that when
+ envseed is non-zero, then we need to use the given environment
+ value. However, we cannot use the same value for all threads because
+ all the positions of all threads will be the same, so we add the
+ thread-id (which is also a known paramter) to the seed. A given number
+ of segments will be distributed between a given number of threads in a
+ reproducible manner, so with this seed value for this thread, we will
+ get reproducible results.*/
+ rng=gsl_rng_alloc(gsl_rng_default);
+ gsl_rng_set(rng, gsl_rng_default_seed+tp->id);
+
+
+ /* Go over the jobs indexed for this thread: */
+ for(i=0; tp->indexs[i] != GAL_THREADS_NON_THRD_INDEX; ++i)
+ {
+ /* Setup this labels parameters */
+ c=0;
+ lab=tp->indexs[i];
+ xw = p->box[ lab*4 + XWIDCOL ];
+ yw = p->box[ lab*4 + YWIDCOL ];
+
+ /* The starting point should be placed between 0 and Image_width
+ - box_width. */
+ while(c<p->upnum)
+ {
+ /* Get a randomly placed start for this label */
+ sum[c] = 0.0f;
+ xmin = gsl_rng_uniform_int(rng, p->s0-xw);
+ ymin = gsl_rng_uniform_int(rng, p->s1-yw);
+
+ /* Go over the pixels and get the sum of the pixel values,
+ this is very similar to the loop in fillseginfo. */
+ for(j=0; j<xw; ++j)
+ {
+ uu = ( u = p->pix[lab] + yw*j ) + yw;
+ l = p->img + (xmin+j)*p->s1 + ymin;
+
+ /* We can't simply multiply, since NaN values might be
+ involved over regions that we don't need.*/
+ do {sum[c] += *u==1.0f ? *l : 0.0f; ++l;} while(++u<uu);
+ }
+
+ /* Only use this sum if it is not a NaN. */
+ /*printf("%lu: %f\n", c, sum[c]);*/
+ if(!isnan(sum[c])) ++c;
+ }
+
+ /* Get the standard deviation by sigma-clipping */
+ gal_fits_change_type (sum, DOUBLE_IMG, c, 0, (void **)&fsum, FLOAT_IMG);
+ gal_statistics_sigma_clip_converge(fsum, 0, c, p->sclipmultip,
+ p->sclipaccu, &ave, &med,
+ &p->std[p->minlab+lab], 0);
+
+ /*
+ printf("%f\n", p->std[ lab*p->numoutcols + p->nimg ]);
+ exit(0);
+
+ p->std[ lab*p->numoutcols + p->nimg ] = gsl_stats_sd(sum, 1, c);
+ */
+ }
+
+ /* Free the random number generator space. */
+ free(sum);
+ gsl_rng_free(rng);
+
+ /* Wait until all other threads finish. When there was only one thread,
+ we explicitly set the pointer to the barrier structure to NULL, so
+ only wait when a barrier is actually defined.*/
+ if(tp->b)
+ pthread_barrier_wait(tp->b);
+
+ /* Return the NULL pointer. */
+ return NULL;
+}
+
+
+
+
+
+/* This is the thread spinner function. For each image we will divide
+ all the objects between the threads and the */
+void
+upperlimit_manager(struct upperlimitparams *p)
+{
+ int err;
+ pthread_t t; /* All thread ids saved in this, not used. */
+ size_t numbarriers;
+ pthread_attr_t attr;
+ pthread_barrier_t b;
+ size_t i, *indexs, thrdcols;
+ struct tupperlimitparams *tp;
+
+ /* Allocate the array of `param' structures. Note that in most cases, the
+ number of threads will not be a constant like this simple case, it
+ will be a variable passed to the thread-spinner. So we are using
+ dynamic allocation for more general use as a tutorial. */
+ errno=0;
+ tp = malloc(p->numthreads*sizeof *tp);
+ if( tp==NULL )
+ error(EXIT_FAILURE, errno, "%lu bytes for tp (upperlimit.c)",
+ p->numthreads*sizeof *tp);
+
+ /* Distribute the actions into the threads: */
+ gal_threads_dist_in_threads(p->numlabs, p->numthreads, &indexs, &thrdcols);
+
+ /* Do the job: when only one thread is necessary, there is no need to
+ spin off one thread, just call the function directly (spinning off
+ threads is expensive). This is for the generic thread spinner
+ function, not this simple function where `p->numthreads' is a
+ constant. */
+ if(p->numthreads==1)
+ {
+ tp[0].p=p;
+ tp[0].id=0;
+ tp[0].b=NULL;
+ tp[0].indexs=indexs;
+ upperlimit_on_thread(&tp[0]);
+ }
+ else
+ {
+ /* Initialize the attributes. Note that this running thread
+ (that spinns off the nt threads) is also a thread, so the
+ number the barriers should be one more than the number of
+ threads spinned off. */
+ numbarriers = ( (p->numlabs<p->numthreads ? p->numlabs : p->numthreads)
+ + 1 );
+ gal_threads_attr_barrier_init(&attr, &b, numbarriers);
+
+ /* Spin off the threads: */
+ for(i=0;i<p->numthreads;++i)
+ if(indexs[i*thrdcols]!=GAL_THREADS_NON_THRD_INDEX)
+ {
+ tp[i].p=p;
+ tp[i].id=i;
+ tp[i].b=&b;
+ tp[i].indexs=&indexs[i*thrdcols];
+ err=pthread_create(&t, &attr, upperlimit_on_thread, &tp[i]);
+ if(err)
+ {
+ fprintf(stderr, "can't create thread %lu", i);
+ exit(EXIT_FAILURE);
+ }
+ }
+
+ /* Wait for all threads to finish and free the spaces. */
+ pthread_barrier_wait(&b);
+ pthread_attr_destroy(&attr);
+ pthread_barrier_destroy(&b);
+ }
+
+ /* Clean up. */
+ free(tp);
+ free(indexs);
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+/*********************************************************************/
+/************ External function ****************/
+/*********************************************************************/
+float *
+upperlimit(float *img, float *sky, long *seg, long *mask, size_t s0,
+ size_t s1, size_t upnum, size_t numthreads, int envseed,
+ float sclipmultip, float sclipaccu)
+{
+ size_t i;
+ struct upperlimitparams p;
+
+ /* Put the input parameters in the structure. */
+ p.s0=s0;
+ p.s1=s1;
+ p.upnum=upnum;
+ p.envseed=envseed;
+ p.sclipaccu=sclipaccu;
+ p.numthreads=numthreads;
+ p.sclipmultip=sclipmultip;
+
+ /* Fill the information for each segment. */
+ fillseginfo(&p, seg);
+
+ /* Prepare the image. */
+ prepareimg(&p, img, sky, seg, mask);
+
+ /* Find the upper limit for all the objects on a thread. */
+ upperlimit_manager(&p);
+
+ /* For a check:
+ for(i=1;i<p.minlab+p.numlabs;++i)
+ printf("i=%lu: %-10.5f\n", i, p.std[i]);
+ exit(0);
+ */
+
+ /* Clean up and return. */
+ for(i=0; i<p.numlabs; ++i)
+ free(p.pix[i]);
+ free(p.img);
+ free(p.pix);
+ free(p.box);
+ return p.std;
+}
diff --git a/bin/mkcatalog/upperlimit.h b/bin/mkcatalog/upperlimit.h
new file mode 100644
index 0000000..e1dff06
--- /dev/null
+++ b/bin/mkcatalog/upperlimit.h
@@ -0,0 +1,33 @@
+/*********************************************************************
+MakeCatalog - Make a catalog from an input and labeled image.
+MakeCatalog 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/>.
+**********************************************************************/
+#ifndef __MKCATALOG_UPPERLIMIT_H__
+#define __MKCATALOG_UPPERLIMIT_H__
+
+
+/* Functions: */
+float *
+upperlimit(float *img, float *sky, long *seg, long *mask, size_t s0,
+ size_t s1, size_t upnum, size_t numthreads, int envseed,
+ float sclipmultip, float sclipaccu);
+
+#endif
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index f6eb667..60743f7 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -439,7 +439,7 @@ Invoking NoiseChisel
MakeCatalog
-* Depth and limiting magnitude:: For comparing different catalogs.
+* Quantifying data limits:: For comparing different catalogs.
* Measuring elliptical parameters:: Estimating elliptical parameters.
* Invoking astmkcatalog:: Options and arguments to MakeCatalog.
* Adding new columns to MakeCatalog:: How to add new columns.
@@ -7024,7 +7024,7 @@ the maximum of the respective pixels in all operands in
the stack.
@item average
Similar to @command{min}, but the pixels of the output will contain
-the average of the respective pixels in all operands in the stack
+the average of the respective pixels in all operands in the stack.
@item median
Similar to @command{min}, but the pixels of the output will contain
@@ -7736,24 +7736,22 @@ the following sum of periodic functions:
f(l)=\displaystyle\sum_{n=-\infty}^{\infty}c_ne^{i{2{\pi}n\over L}l} }
@noindent
-Note that the different frequencies (@mymath{2{\pi}n/L}, in units of
-cycles per meters for example) are not arbitrary. They are all integer
-multiples of the fundamental frequency of
address@hidden/L}. Recall that @mymath{L} was the length of the
-signal we want to model. Therefore, we see that the smallest possible
-frequency (or the frequency resolution) in the end, depends on the
-length we observed the signal or @mymath{L}. In the case of each
-dimension on an image, this is the size of the image in the respective
-dimension. The frequencies have been defined in this ``harmonic''
-fashion to insure that the final sum is periodic outside of the
address@hidden, l_0+L]} interval too. At this point, you might be
-thinking that the sky is not periodic with the same period as my
-camera's view angle. You are absolutely right! The important thing is
-that since your camera's observed region is the only region we are
-``observing'' and will be using, the rest of the sky is irrelevant; so
-we can safely assume the is periodic outside of it. However, this
-working assumption will haunt us later in @ref{Edges in the frequency
-domain}.
+Note that the different frequencies (@mymath{2{\pi}n/L}, in units of cycles
+per meters for example) are not arbitrary. They are all integer multiples
+of the fundamental frequency of @mymath{\omega_0=2\pi/L}. Recall that
address@hidden was the length of the signal we want to model. Therefore, we see
+that the smallest possible frequency (or the frequency resolution) in the
+end, depends on the length we observed the signal or @mymath{L}. In the
+case of each dimension on an image, this is the size of the image in the
+respective dimension. The frequencies have been defined in this
+``harmonic'' fashion to insure that the final sum is periodic outside of
+the @mymath{[l_0, l_0+L]} interval too. At this point, you might be
+thinking that the sky is not periodic with the same period as my camera's
+view angle. You are absolutely right! The important thing is that since
+your camera's observed region is the only region we are ``observing'' and
+will be using, the rest of the sky is irrelevant; so we can safely assume
+the sky is periodic outside of it. However, this working assumption will
+haunt us later in @ref{Edges in the frequency domain}.
The frequencies are thus determined by definition. So all we need to
do is to find the coefficients (@mymath{c_n}), or magnitudes, or radii
@@ -11281,7 +11279,7 @@ transformations on the labeled image.
@menu
-* Depth and limiting magnitude:: For comparing different catalogs.
+* Quantifying data limits:: For comparing different catalogs.
* Measuring elliptical parameters:: Estimating elliptical parameters.
* Invoking astmkcatalog:: Options and arguments to MakeCatalog.
* Adding new columns to MakeCatalog:: How to add new columns.
@@ -11290,132 +11288,161 @@ transformations on the labeled image.
address@hidden Depth and limiting magnitude, Measuring elliptical parameters,
MakeCatalog, MakeCatalog
address@hidden Depth and limiting magnitude
address@hidden Quantifying data limits, Measuring elliptical parameters,
MakeCatalog, MakeCatalog
address@hidden Quantifying data limits
@cindex Depth
@cindex Clump magnitude limit
@cindex Object magnitude limit
@cindex Limit, object/clump magnitude
@cindex Magnitude, object/clump detection limit
-Different observations have different noise properties and different
-detection methods (or one method with a different set of parameters)
-will have different abilities to detect certain objects in an
-image. Therefore it is very important that there be a scale on which
-we can compare different observations (images) and detection methods
-to objectively quantify the noise and our ability to detect signal in
-it.
-
-Due to the presence of correlated noise in processed images (which are
-used for scientific deductions), we cannot simply deduce the limiting
-signal properties from those of the noise. Hence a different measure
-is needed for each. To quantify the level of noise, we define
address@hidden and to quantify the ability to reliably detect/study
-objects with that methodology we define the magnitude limit. In
-astronomy, it is common to use the magnitude (a unit-less scale) and
-not physical units, see @ref{Flux Brightness and magnitude}. Therefore
-both these measures will be in units of magnitudes, but since
-magnitudes don't have units, we are just showing them like units as a
-place holder for clarity.
+Different datasets (images in the case of MakeCatalog) have different noise
+properties and different detection methods (or one method with a different
+set of parameters) will have different abilities to detect or measure
+certain kinds of objects in an image. Therefore it is very important to
+quantify our ability to detect and measure signal in noise. In this section
+we discuss these limits that are very important in any analysis.
+
+In astronomy, it is common to use the magnitude (a unit-less scale) and
+physical units, see @ref{Flux Brightness and magnitude}. Therefore all the
+measured discussed here are commonly defined in units of magnitudes.
@table @asis
@item Depth
-[magnitude] As we make more observations on one region of the sky and
-add the images over each other, we are able to decrease the standard
-deviation of the noise in each address@hidden is true for any
-noisy data, not just astronomical images}. Qualitatively, this
-decrease manifests its self by making fainter (per pixel) parts of the
-objects in the image more visible. Quantitatively, it increases the
-Signal to noise ratio, since the signal is fixed but the noise is
-decreased. It is very important to have in mind that noise is defined
-per pixel (or independent data measurement unit), not per object.
-
-You can think of the noise as muddy water that is completely covering
-a flat address@hidden ground is the sky value in this analogy,
-see @ref{Sky value}. Note that this analogy only holds for a flat sky
-value across the surface of the image or ground.} with some
address@hidden hills are astronomical objects in this
-analogy. Brighter parts of the object are higher from the ground.} in
-it. Let's assume that in your first observation the muddy water has
-just been stirred and you can't see anything through it. As you wait
-and make more observations, the mud settles down and the @emph{depth}
-of the transparent water increases as you wait. The summits of hills
-begin to appear. As the depth of clear water increases, the parts of
-the hills with lower heights can be seen more clearly.
+As we make more observations on one region of the sky and add the images
+over each other, we are able to decrease the standard deviation of the
+noise in each address@hidden is true for any noisy data, not just
+astronomical images}. Qualitatively, this decrease manifests its self by
+making fainter (per pixel) parts of the objects in the image more
+visible. Quantitatively, it increases the Signal to noise ratio, since the
+signal is fixed but the noise is decreased. It is very important to have in
+mind that here, noise is defined per pixel (or in the units of our data
+measurement), not per object.
+
+You can think of the noise as muddy water that is completely covering a
+flat address@hidden ground is the sky value in this analogy, see
address@hidden value}. Note that this analogy only holds for a flat sky value
+across the surface of the image or ground.} with some regions higher than
+the address@hidden peaks are (parts of) astronomical objects in this
+analogy.} in it. In this analogy, height is brightness. Let's assume that
+in your first observation the muddy water has just been stirred and you
+can't see anything through it. As you wait and make more observations, the
+mud settles down and the @emph{depth} of the transparent water increases as
+you wait. The summits of hills begin to appear. As the depth of clear water
+increases, the parts of the hills with lower heights (less surface
+brightness) can be seen more clearly.
@cindex Depth
The outputs of NoiseChisel include the Sky standard deviation
-(@mymath{\sigma}) on every group of pixels (a mesh) that were
-calculated from the undetected pixels in that mesh, see @ref{Tiling an
-image} and @ref{NoiseChisel output}. Let's take @mymath{\sigma_m} as
-the median @mymath{\sigma} over the successful meshes in the image
-(prior to interpolation or smoothing, see @ref{Grid interpolation and
-smoothing}). Note that even though on different instruments, pixels
-have different physical sizes (for example in @mymath{\mu}m),
-nevertheless, a pixel is the unit of data collection. Therefore, as
-far as noise is considered, the physical or projected size of the
-pixels is irrelevant. We thus define the @emph{depth} of each data set
-as the magnitude of @mymath{\sigma_m}.
-
-As an example, the XDF survey covers part of the sky that the Hubble
-space telescope has observed the most (for 85 orbits) and is
-consequently very small (@mymath{\sim4} address@hidden). On the
-other hand, the CANDELS survey, is one of the widest multi-color
-surveys covering several fields (about 720 address@hidden) but its
-deepest fields have only 9 orbits observation. The depth of the XDF
-and CANDELS-deep surveys in the near infrared WFC3/F160W filter are
-respectively 34.40 and 32.45 magnitudes. In a single orbit image, this
-same field has a depth of 31.32. Note that a larger magnitude
-corresponds to less brightness, see @ref{Flux Brightness and
-magnitude}.
-
address@hidden Magnitude limit
address@hidden False detections
address@hidden Detections false
-[magnitude] Because of noise and methodology, no detection algorithm
-can be perfect and this parameter specifies the limits of the combined
-data and methodology used for detection. Assuming a fixed shape for
-all the targets, for any algorithm and its accompanying set of input
-parameters, there is always a certain magnitude limit below which the
-detections (pseudo-detections or clumps in NoiseChisel, see
address@hidden) are no longer usable/reliable.
+(@mymath{\sigma}) on every group of pixels (a mesh) that were calculated
+from the undetected pixels in that mesh, see @ref{Tiling an image} and
address@hidden output}. Let's take @mymath{\sigma_m} as the median
address@hidden over the successful meshes in the image (prior to
+interpolation or smoothing, see @ref{Grid interpolation and smoothing}).
+
+On different instruments pixels have different physical sizes (for example
+in micro-meters, or spatial angle over the sky), nevertheless, a pixel is
+our unit of data collection. In other words, while quantifying the noise,
+the physical or projected size of the pixels is irrelevant. We thus define
+the @emph{depth} of each dataset (or image) as the magnitude of
address@hidden
+
address@hidden XDF
address@hidden CANDELS
+As an example, the XDF survey covers part of the sky that the Hubble space
+telescope has observed the most (for 85 orbits) and is consequently very
+small (@mymath{\sim4} address@hidden). On the other hand, the CANDELS
+survey, is one of the widest multi-color surveys covering several fields
+(about 720 address@hidden) but its deepest fields have only 9 orbits
+observation. The depth of the XDF and CANDELS-deep surveys in the near
+infrared WFC3/F160W filter are respectively 34.40 and 32.45 magnitudes. In
+a single orbit image, this same field has a depth of 31.32. Recall that a
+larger magnitude corresponds to less brightness, see @ref{Flux Brightness
+and magnitude}.
+
address@hidden Upper limit magnitude
+Due to the noisy nature of data, it is possible to get arbitrarily low
+values for a faint object's brightness (or arbitrarily high
+magnitudes). Given the scatter caused by the noise, such small values are
+meaningless: another similar depth observation will give a radically
+different value. This problem is most common when you use one image/filter
+to generate target labels (which specify which pixels belong to which
+object, see @ref{NoiseChisel output} and @ref{MakeCatalog}) and another
+image/filter to generate a catalog for measuring colors.
+
+The object might not be visible in the filter used for the latter image, or
+the image @emph{depth} (see above) might be much shallower. So you will get
+unreasonbly faint magnitudes. For example when the depth of the image is 32
+magnitudes, a measurement that gives a magnitude of 36 for a
address@hidden pixel object is clearly unreliable. In another similar
+depth image, we might measure a magnitude of 30 for it, and yet another
+might give 33. Furthermore, due to the noise scatter so close to the depth
+of the dataset, the total brightness might actually get measured as a
+negative value, so no magnitude can be defined (recall that a magnitude is
+a base-10 logarithm).
+
address@hidden Upper limit magnitude
address@hidden Magnitude, upper limit
+Using such unreliable measurements will directly affect our analysis, so we
+must not use them. However, all is not lost! Given our limited depth, there
+is one thing we can deduce about the object's magnitude: we can say that if
+something actually exists here (possibly burried deep under the noise), it
+must have a magnitude that is fainter than an @emph{upper limit
+magnitude}. To find this upper limit magnitude, we place the object's
+footprint (segmentation map) over random parts of the image where there are
+no detections, so we only have pure (possibly correlated) noise and
+undetected objects. Doing this a large number of times will give us a
+distribution of brightness values. The standard deviation (@mymath{\sigma})
+of that distribution can be used to quantify the upper limit magnitude.
@cindex Correlated noise
-While adding more data sets does have the advantage of decreasing the
-standard deviation of the noise, it also produces correlated
-noise. Correlated noise is produced because the raw data sets are
-warped (rotated, shifted or resampled in general, see
address@hidden) before they are added with each other. This
-correlated noise manifests as a `smoothing' or `blurring' over the
-image. Therefore pixels in added images are no longer separate or
-independently measured data elements, they are `correlated' and this
-produces a hurdle in our ability to detect objects in them.
-
address@hidden Number count
-To find the limiting magnitude, you have to use the output of
-MakeCatalog and plot the number of objects as a function of magnitude
-with your favorite plotting tool, this is called a ``number count''
-plot. It is simply a histogram of the catalog in each magnitude
-bin. This histogram can be used in many ways to specify a magnitude
-limit, for example see Akhlaghi et al. (2015, in preparation) for one
-method of using multiple depth images in order to find this limit.
+Traditionally, faint/small object photometry was done using fixed circular
+apertures (for example with a diameter of @mymath{N} arcseconds). In this
+way, the upper limit was like the depth discussed above: one value for the
+whole image. But with the much more advanced hardware and software of
+today, we can make customized segmentation maps for each object. The number
+of pixels (are of the object) used directly affects the final distribution
+and thus magnitude. Also the image correlated noise might actually create
+certain patters, so the shape of the object can also affect the result. So
+in MakeCatalog, the upper limit magnitude is found for each object in the
+image separately. Not one value for the whole image.
+
address@hidden Completeness limit
address@hidden Completeness
+As the surface brightness of the objects decreases, the ability to detect
+them will also decrease. An important statistic is thus the fraction of
+objects of similar morphology and brightness that will be identified with
+our detection algorithm/parameters in the given image. This fraction is
+known as completeness. For brighter objects, completeness is 1: all bright
+objects that might exist over the imagde will be detected. However, as we
+go to lower surface brightness objects, we fail to detect some and
+gradually we are not able to detect anything any more. For a given profile,
+the magnitude where the completeness drops below a certain level usually
+above @mymath{90\%} is known as the completeness limit.
+
address@hidden Purity
address@hidden False detections
address@hidden Detections false
+Another important parameter in measuring completeness is purity: the
+fraction of true detections to all true detections. In effect purity is the
+measure of contamination by false detections: the higher the purity, the
+lower the contamination. Completeness and purity are anti-correlated: if we
+can allow a large number of false detections (that we might be able to
+remove by other means), we can significantly increase theh completeness
+limit.
+
+One traditional way to measure the completeness and purity of a given
+sample is by embedding mock profiles in regions of the image with no
+detection. However in such a study we must be really careful to choose
+model profiles as similar to the target of interest as possible.
@end table
address@hidden Zero-point magnitude
address@hidden
-For any data-set and detection algorithm (with a specific set of
-parameters), the depth and limiting magnitudes can differ. The first
-is reported in the comments section of the catalog plain text
-file. Note that the accuracy in estimating the zero-point magnitude is
-a very important factor in an accurate comparison between magnitudes
-measured for different images, on possibly different instruments and
-cameras, see @ref{Flux Brightness and magnitude}.
address@hidden Measuring elliptical parameters, Invoking astmkcatalog, Depth
and limiting magnitude, MakeCatalog
address@hidden Measuring elliptical parameters, Invoking astmkcatalog,
Quantifying data limits, MakeCatalog
@subsection Measuring elliptical parameters
One of the most important class of parameters that are important for
@@ -11807,7 +11834,92 @@ in more accurate floating point display.
@end table
address@hidden
+The upper limit magnitude was discussed in @ref{Quantifying data
+limits}. Unlike other measured values/columns in MakeCatalog, the upper
+limit magnitude needs several defined parameters which are discussed
+here. All the upper limit magnitude specific options start with @option{up}
+for upper limit, except for @option{--envseed} that is also present in
+other programs and is general for any job requiring random number
+generation.
+
address@hidden Reproducibility
+One very important consideration in Gnuastro is reproducibility. Therefore,
+the values to all of these parameters along with others (like the random
+number generator type and seed) are also reported in the comments of the
+final catalog when the upper limit magnitude column is desired. One of
+those parameters is the number of threads used in this run of MakeCatalog.
+This might seem irrelevant.
+
+The job of sampling a large number of positions over the image is an
+``embarrasingly parallel'' kind of problem which can greatly benefit from
+threads. For maximum efficiency, the objects in the catalog are divided
+between the threads and the random number generator is used independently
+on each thread. So even when the random number generator type and seed are
+identical, if the number of threads differ, you will get different
+results. Ofcourse, the difference is due to scatter and is statistically
+negligible, but it is not exactly reproducible (which is important in some
+contexts).
+
address@hidden @option
+
address@hidden --upmask
+(@option{=STR}) File name of mask image to use for upper-limit
+calculation. In some cases (especially when doing matched photometry), the
+object labels specified in the main input and mask image might not be
+adequate. In other words they do not necessarily have to cover @emph{all}
+detected objects: the user might have selected only a few of the objects in
+their labeled image. All the non-zero pixels of the image specified by this
+option (in the @option{--upmaskhdu} extension) will be set as blank for the
+upper limit magnitude calculation.
+
+For example, when you are using labels from another image, you can give
+NoiseChisel's objects image output for this image as the value to this
+option. In this way, you can be sure that regions with data do not harm
+your distribution. See @ref{Quantifying data limits} for more on the upper
+limit magnitude.
+
address@hidden --upmaskhdu
+(@option{=STR}) The extension in the file specified by @option{--upmask}.
+
address@hidden --upnum
+(@option{=INT}) The number of random samples to take for all the objects. A
+larger value to this option will give a more accurate result
+(assymptotically), but it will also slow down the process. When a randomly
+positioned sample overlaps with a detected/masked pixel it is not counted
+and another random position is found until the object completely lies over
+an undetected region. So you can be sure that for each object, this many
+samples over undetected objects are made. See the upper limit magnitude
+discussion in @ref{Quantifying data limits} for more.
+
address@hidden --envseed
+Read the random number generator type and seed value from the environment
+(see @ref{Generating random numbers}). Random numbers are used in
+calculating the random positions of different samples of each object.
+
address@hidden --upsclipmultip
+(@option{=FLT}) The multiple of @mymath{\sigma} to use in
address@hidden the final distribution (see @ref{Sigma
+clipping}). The images might have some artifacts or some regions with
+signal in the image might not have been removed correctly. If a random
+sampling falls over such regions, they can significantly bias the final
+standard deviation. To avoid the effect of such outliers, after the
+distribution is found, it is @mymath{\sigma}-clipped (by convergence).
+
address@hidden --upsclipaccu
+(@option{=FLT}) Relative difference between two successive @mymath{\sigma}
+measurements to halt the @mymath{\sigma}-clipping process by convergence,
+see the explanations for @option{--upsclipmultip} for more.
+
address@hidden --upnsigma
+(@option{=FLT}) The multiple of the standard deviation (or @mymath{\sigma})
+used to measure the magnitude. Note that this is the final @mymath{\sigma}
+produced after the @mymath{\sigma}-clipping.
+
address@hidden table
+
+
+
+
The final group of options particular to MakeCatalog are those that
specify which columns should be displayed in the output catalogs. For
each column there is an option, if it has been called on the command
@@ -12004,6 +12116,14 @@ defined for work on adding these sources of error too.
[Objects] The magnitude of all clumps in this object, see
@option{--clumpbrightness}.
address@hidden --upperlimitmag
+[Objects] The upper limit magnitude for this object. The object's footprint
+is used over randomly positioned parts of the image that do not cover a
+detected object. The standard deviation of the final distribution is then
+the upper limit magnitude, see @ref{Quantifying data limits} for a complete
+explanation. This is very important for the fainter objects in the image
+where the measured magnitudes are not reliable.
+
@item --riverave
[Clumps] The average brightness of the river pixels around this
clump. River pixels were defined in Akhlaghi and Ichikawa 2015. In
@@ -12094,26 +12214,26 @@ will need to also modify the respective
@code{firstpass} and
@code{secondpass} functions and define new information table columns
in @file{main.h}.
-For all the steps below, it is easiest to just copy and paste an
-existing option and change the variables in each step. In all these
-different places, the options are sorted in the same order (same order
-as @ref{Invoking astmkcatalog}). This allows a particular option to be
-easily found in all steps. Therefore in adding your new option, be
-sure to keep it in the same relative place in the list (it doesn't
-necessarily have to be in the end), and near conceptually similar
-options.
+For all the steps below, it is easiest to just copy and paste an existing
+option and change the variables in each step. In all these different
+places, the options are sorted in the same order (same order as
address@hidden astmkcatalog}). This allows a particular column/option to be
+easily found in all steps. Therefore in adding your new option, be sure to
+keep it in the same relative place in the list in all the separate places
+(it doesn't necessarily have to be in the end), and near conceptually
+similar options.
@table @file
@item main.h
The @code{objectcols} and @code{clumpcols} enumerated variables
-(@code{enum}) define the labels for the separate information table
-columns (the raw pixel calculations). If your new calculation requires
-new raw columns, the add them. The @code{outcols} enumerated variables
-define the recognized output columns. You will need to add a new name
-for your desired column here, then add a row for your new parameter in
-the @code{uiparams} structure, the name should be the same as the one
-you used in @code{outcols} (in small caps).
+(@code{enum}) define the labels for the separate information table columns
+(the raw pixel calculations). If your new calculation requires new raw
+columns, the add new macros for them. The @code{outcols} enumerated
+variables define the recognized output columns. You will need to add a new
+name for your desired column here, then add a row for your new parameter in
+the @code{uiparams} structure, the name should be the same as the one you
+used in @code{outcols} (in small caps).
@cindex GNU C library
@item args.h
@@ -12123,26 +12243,28 @@ that the option name be the same as the name you put
in @code{outcols}
(in @file{main.h}), like the other options that are already
defined.
-Define a new option in @code{argp_option} (rows between the @address@hidden
-and @address@hidden). The first is the name of the option, the second is its
-short option name (see @ref{Options}), it is best to use a number (so
-no short option is defined for the user) which in Gnuastro is defined
-to be a number larger than 500. All the used numbers are shown in the
-comment immediately above @code{argp_option} structure, increase that
-number by one and use that number for your new option. The next two
-values should be @code{0} and the last should be a short description
-of your option.
+Define a new option in @code{argp_option} (rows between the @address@hidden and
address@hidden@}}). The first is the name of the option, the second is its short
+option name (see @ref{Options}). The next two values should be @code{0} and
+the last should be a short description of your option/column. For the
+second element is best to use a number (so no short option is defined for
+the user) which in Gnuastro is defined to be a number larger than 500. The
+range of used numbers (and short option names) are shown in the comment
+immediately above @code{argp_option} structure, increase that number by one
+and use that number for your new option.
-Go down to the @code{parse_opt} function and add a @code{case} for
-your specific option (using the same number you specified above), then
-update the contents.
+Go down to the @code{parse_opt} function and add a @code{case} for your
+specific option (using the same number you specified above), then update
+the contents, based on the new macros you defined @code{main.h}.
@item ui.c
-In the @code{readconfig} function, add an @code{else if} for your
-option (easiest to just copy and paste, then update). In the
address@hidden function's @code{switch} structure, add a
address@hidden for your option. Finally, in the @code{preparearrays}
-function, add a @code{case} for your new output column.
+In the @code{readconfig} function, add an @code{else if} for your option
+(easiest to just copy and paste, then update). In the @code{printvalues}
+function's @code{switch} structure, add a @code{case} for your
+option. Finally, in the @code{preparearrays} function, add a @code{case}
+for your new output column. Note that if the column is to be used for both
+objects and clumps, you need to set both @code{p->objcols} and
address@hidden>clumpcols}.
@item mkcatalog.c
If your procedure needs new raw calculations, add a new line to
@@ -12366,14 +12488,14 @@ namely the pixels, also produces a very small
`spread' in the image.
@cindex Convolution
@cindex Image blurring
@cindex PSF image size
-Convolution is the mathematical process by which we can apply a
-`spread' to an image, or in other words blur the image, see
address@hidden process}. The Brightness of an object should remain
-unchanged after convolution, see @ref{Flux Brightness and
-magnitude}. Therefore, it is important that the sum of all the pixels
-of the PSF be unity. The image also has to have an odd number of
-pixels on its sides so one pixel can be defined as the center. In
-MakeProfiles, the PSF can be set by the two methods explained below.
+Convolution is the mathematical process by which we can apply a `spread' to
+an image, or in other words blur the image, see @ref{Convolution
+process}. The Brightness of an object should remain unchanged after
+convolution, see @ref{Flux Brightness and magnitude}. Therefore, it is
+important that the sum of all the pixels of the PSF be unity. The PSF image
+also has to have an odd number of pixels on its sides so one pixel can be
+defined as the center. In MakeProfiles, the PSF can be set by the two
+methods explained below.
@table @asis
@@ -15095,7 +15217,7 @@ defines two major formats: images and tables. An image
has a single type
for all its data elements (a pixel in a 2D image). A table has a single
type for every one of its columns. In both cases, the FITS header stores
this information. All the FITS data types correspond to C types (for
-example @code{short}, @code{int}, @code{float}, or @code{long}.
+example @code{short}, @code{int}, @code{float}, or @code{long}).
@cindex CFITSIO
@cindex @code{BITPIX}
diff --git a/tests/mkcatalog/simple.sh b/tests/mkcatalog/simple.sh
index ca3cb6a..0ff05bd 100755
--- a/tests/mkcatalog/simple.sh
+++ b/tests/mkcatalog/simple.sh
@@ -48,4 +48,4 @@ if [ ! -f $execname ] || [ ! -f $img ]; then exit 77; fi
# Actual test script
# ==================
-$execname $img --sn --magnitude --dec --ra --y --x
+$execname $img --sn --upperlimitmag --magnitude --dec --ra --y --x
- [gnuastro-commits] master updated (958018d -> c074b09), Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master 8644305 5/6: Correctly report conf files checked, Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master c074b09 6/6: Duplicate filename in MakeCatalog's output removed, Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master 008df03 2/6: Minor random small corrections, Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master b9c3d30 1/6: envseed as a conf. option in MakeNoise and MakeProfiles, Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master 708090b 4/6: Non-installed headers moved to lib/, Mohammad Akhlaghi, 2016/10/02
- [gnuastro-commits] master 3600a5e 3/6: Upper limit magnitude calculation in MakeCatalog,
Mohammad Akhlaghi <=