[gnuastro-commits] master c70aa98: MakeCatalog: new column to measure the error in surface brightness

[Mohammad Akhlaghi]

branch: master
commit c70aa98219a7db563e03f11e0522074a31295dfe
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    MakeCatalog: new column to measure the error in surface brightness
    
    Until now there was no specific column to measure the error in surface
    brightness. Although it is very similar to the error in measuring the
    magnitude, it includes an area term that should also be taken into account.
    
    With this commit, two new columns are added to MakeCatalog:
    '--spatialresolution' and '--sberror'. The first is the error in measuring
    area (or the FWHM of the PSF) which should be specified by the user. The
    second is the measured error in measuring surface brightness that is a
    column in the output catalog.
    
    This task was suggested by Zahra Sharbaf and Juan Miro.
---
 NEWS                            |  5 ++++
 bin/mkcatalog/args.h            | 27 +++++++++++++++++++++
 bin/mkcatalog/astmkcatalog.conf |  3 +++
 bin/mkcatalog/columns.c         | 34 +++++++++++++++++++++++++-
 bin/mkcatalog/main.h            |  1 +
 bin/mkcatalog/ui.h              |  2 ++
 doc/announce-acknowledge.txt    |  1 +
 doc/gnuastro.texi               | 53 +++++++++++++++++++++++++++++++----------
 8 files changed, 113 insertions(+), 13 deletions(-)

diff --git a/NEWS b/NEWS
index 26403c4..7db74de 100644
--- a/NEWS
+++ b/NEWS
@@ -27,6 +27,11 @@ See the end of the file for license conditions.
      to the Hubble flow (which is the basis of the calculations here). The
      warning can be suppressed with the '--quiet' option.
 
+  MakeCatalog:
+   --areaerror: spatial resolution of image specified by user, used in
+     estimating the surface brightness error.
+   --sberror: error in measuring the surface brightness (mag/arcsec^2).
+
   Match:
    - When called with '--notmatched --outcols=AAA,BBB', Match will append
      non-matching rows of second table into first table's rows (for columns
diff --git a/bin/mkcatalog/args.h b/bin/mkcatalog/args.h
index a1c18c6..ef65086 100644
--- a/bin/mkcatalog/args.h
+++ b/bin/mkcatalog/args.h
@@ -427,6 +427,19 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET,
       gal_options_parse_csv_float64
     },
+    {
+      "spatialresolution",
+      UI_KEY_SPATIALRESOLUTION,
+      "FLT",
+      0,
+      "Spatial resolution (for surface brightness error).",
+      UI_GROUP_OTHERSETTINGS,
+      &p->spatialresolution,
+      GAL_TYPE_FLOAT32,
+      GAL_OPTIONS_RANGE_GE_0,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+    },
 
 
 
@@ -1519,6 +1532,20 @@ struct argp_option program_options[] =
       ui_column_codes_ll
     },
     {
+      "sberror",
+      UI_KEY_SBERROR,
+      0,
+      0,
+      "Surface brightness error (mag/arcsec^2).",
+      UI_GROUP_COLUMNS_MORPHOLOGY,
+      0,
+      GAL_TYPE_INVALID,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET,
+      ui_column_codes_ll
+    },
+    {
       "areaxy",
       UI_KEY_AREAXY,
       0,
diff --git a/bin/mkcatalog/astmkcatalog.conf b/bin/mkcatalog/astmkcatalog.conf
index 6815d2e..cbb2901 100644
--- a/bin/mkcatalog/astmkcatalog.conf
+++ b/bin/mkcatalog/astmkcatalog.conf
@@ -37,4 +37,7 @@
  upsigmaclip    3,0.2
  upnsigma           1
 
+# Settings for other columns:
+ spatialresolution  2
+
 # Catalog columns:
diff --git a/bin/mkcatalog/columns.c b/bin/mkcatalog/columns.c
index fb6ded7..b31853f 100644
--- a/bin/mkcatalog/columns.c
+++ b/bin/mkcatalog/columns.c
@@ -205,6 +205,7 @@ columns_wcs_preparation(struct mkcatalogparams *p)
             /* High-level. */
             case UI_KEY_RA:
             case UI_KEY_DEC:
+            case UI_KEY_SBERROR:
             case UI_KEY_HALFMAXSB:
             case UI_KEY_HALFSUMSB:
             case UI_KEY_AREAARCSEC2:
@@ -255,6 +256,7 @@ columns_wcs_preparation(struct mkcatalogparams *p)
         break;
 
       /* Calculate the pixel area if necessary. */
+      case UI_KEY_SBERROR:
       case UI_KEY_HALFMAXSB:
       case UI_KEY_HALFSUMSB:
       case UI_KEY_AREAARCSEC2:
@@ -485,6 +487,24 @@ columns_define_alloc(struct mkcatalogparams *p)
           oiflag[ OCOL_SUM ] = ciflag[ CCOL_SUM ] = 1;
           break;
 
+        case UI_KEY_SBERROR:
+          name           = "SB_ERROR";
+          unit           = "mag/arcsec^2";
+          ocomment       = "Error in measuring Surface brightness.";
+          ccomment       = ocomment;
+          otype          = GAL_TYPE_FLOAT32;
+          ctype          = GAL_TYPE_FLOAT32;
+          disp_fmt       = 0;
+          disp_width     = 6;
+          disp_precision = 0;
+          oiflag[ OCOL_NUM         ] = ciflag[ CCOL_NUM         ] = 1;
+          oiflag[ OCOL_SUM         ] = ciflag[ CCOL_SUM         ] = 1;
+          oiflag[ OCOL_SUM_VAR     ] = ciflag[ CCOL_SUM_VAR     ] = 1;
+          oiflag[ OCOL_SUM_VAR_NUM ] = ciflag[ CCOL_SUM_VAR_NUM ] = 1;
+                                       ciflag[ CCOL_RIV_SUM     ] = 1;
+                                       ciflag[ CCOL_RIV_SUM_VAR ] = 1;
+          break;
+
         case UI_KEY_AREAXY:
           name           = "AREAXY";
           unit           = "counter";
@@ -2262,7 +2282,11 @@ columns_xy_extrema(struct mkcatalog_passparams *pp, 
double *oi,
                                       ? columns_sn((P),(ROW),(O0C1))    \
                                       : NAN )                           \
                                     * log(10) ) )
-
+#define SB_ERROR(P,ROW,O0C1) ( MAG_ERROR(P,ROW,O0C1)                    \
+                               + ( 2.5/log(10)*p->spatialresolution     \
+                                   / ( ROW[   O0C1?CCOL_NUM:OCOL_NUM ]  \
+                                       ? ROW[ O0C1?CCOL_NUM:OCOL_NUM ]  \
+                                       : NAN ) ) )
 
 
 
@@ -2340,6 +2364,10 @@ columns_fill(struct mkcatalog_passparams *pp)
           ((float *)colarr)[oind] = MKC_SB(oi[OCOL_SUM], oi[OCOL_NUM]);
           break;
 
+        case UI_KEY_SBERROR:
+          ((float *)colarr)[oind] = SB_ERROR(p, oi, 0);
+          break;
+
         case UI_KEY_AREAXY:
           ((int32_t *)colarr)[oind] = oi[OCOL_NUMXY];
           break;
@@ -2749,6 +2777,10 @@ columns_fill(struct mkcatalog_passparams *pp)
             ((float *)colarr)[cind]=MKC_SB(ci[CCOL_SUM], ci[CCOL_NUM]);
             break;
 
+          case UI_KEY_SBERROR:
+            ((float *)colarr)[cind]=SB_ERROR(p, ci, 1);
+            break;
+
           case UI_KEY_AREAXY:
             ((int32_t *)colarr)[cind]=ci[CCOL_NUMXY];
             break;
diff --git a/bin/mkcatalog/main.h b/bin/mkcatalog/main.h
index aa2dd34..954cec2 100644
--- a/bin/mkcatalog/main.h
+++ b/bin/mkcatalog/main.h
@@ -244,6 +244,7 @@ struct mkcatalogparams
   int32_t       checkuplim[2];  /* Object & clump ID to check dist.     */
 
   gal_data_t         *fracmax;  /* Fractions to use in --fracsumarea.   */
+  float     spatialresolution;  /* Error in area (used in SB error).    */
 
   /* Internal. */
   char           *relabclumps;  /* Name of new file for clump labels.   */
diff --git a/bin/mkcatalog/ui.h b/bin/mkcatalog/ui.h
index 3d390da..eca4d3e 100644
--- a/bin/mkcatalog/ui.h
+++ b/bin/mkcatalog/ui.h
@@ -105,11 +105,13 @@ enum option_keys_enum
   UI_KEY_CHECKUPLIM,
   UI_KEY_NOCLUMPSORT,
   UI_KEY_FRACMAX,
+  UI_KEY_SPATIALRESOLUTION,
 
   UI_KEY_OBJID,                         /* Catalog columns. */
   UI_KEY_IDINHOSTOBJ,
   UI_KEY_AREAARCSEC2,
   UI_KEY_SURFACEBRIGHTNESS,
+  UI_KEY_SBERROR,
   UI_KEY_AREAXY,
   UI_KEY_CLUMPSAREA,
   UI_KEY_WEIGHTAREA,
diff --git a/doc/announce-acknowledge.txt b/doc/announce-acknowledge.txt
index 89a4c5d..0a6e709 100644
--- a/doc/announce-acknowledge.txt
+++ b/doc/announce-acknowledge.txt
@@ -2,6 +2,7 @@ Alphabetically ordered list to acknowledge in the next release.
 
 Mark Calabretta
 Leslie Hunt
+Juan Miro
 Juan Molina Tobar
 Zahra Sharbaf
 Ignacio Trujillo
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 157c491..af397db 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -409,7 +409,7 @@ ConvertType
 
 Annotations for figure in paper
 
-* Full script of annotations on figure::
+* Full script of annotations on figure::  All the steps in one script
 
 Table
 
@@ -540,6 +540,7 @@ MakeCatalog
 Quantifying measurement limits
 
 * Magnitude measurement error of each detection::  Derivation of mag error 
equation
+* Surface brightness error of each detection::  Error in measuring the Surface 
brightness.
 * Completeness limit of each detection::  Possibility of detecting similar 
objects?
 * Upper limit magnitude of each detection::  How reliable is your magnitude?
 * Surface brightness limit of image::  How deep is your data?
@@ -10375,7 +10376,7 @@ PGFPlots is also highly customizable, you can make a 
lot of changes and customiz
 Both 
TiKZ@footnote{@url{http://mirrors.ctan.org/graphics/pgf/base/doc/pgfmanual.pdf}}
 and 
PGFPLots@footnote{@url{http://mirrors.ctan.org/graphics/pgf/contrib/pgfplots/doc/pgfplots.pdf}}
 have wonderful manuals, so have a look trough them.
 
 @menu
-* Full script of annotations on figure:: All the steps in one script
+* Full script of annotations on figure::  All the steps in one script
 @end menu
 
 @node Full script of annotations on figure,  , Annotations for figure in 
paper, Annotations for figure in paper
@@ -17942,13 +17943,14 @@ Therefore the measurements discussed here are 
commonly used in units of magnitud
 
 @menu
 * Magnitude measurement error of each detection::  Derivation of mag error 
equation
+* Surface brightness error of each detection::  Error in measuring the Surface 
brightness.
 * Completeness limit of each detection::  Possibility of detecting similar 
objects?
 * Upper limit magnitude of each detection::  How reliable is your magnitude?
 * Surface brightness limit of image::  How deep is your data?
 * Upper limit magnitude of image::  How deep is your data for certain 
footprint?
 @end menu
 
-@node Magnitude measurement error of each detection, Completeness limit of 
each detection, Quantifying measurement limits, Quantifying measurement limits
+@node Magnitude measurement error of each detection, Surface brightness error 
of each detection, Quantifying measurement limits, Quantifying measurement 
limits
 @subsubsection Magnitude measurement error of each detection
 The raw error in measuring the magnitude is only meaningful when the object's 
magnitude is brighter than the upper-limit magnitude (see below).
 As discussed in @ref{Brightness flux magnitude}, the magnitude (@mymath{M}) of 
an object with brightness @mymath{B} and zero point magnitude @mymath{z} can be 
written as:
@@ -17973,7 +17975,22 @@ But, @mymath{\Delta{B}/B} is just the inverse of the 
Signal-to-noise ratio (@mym
 MakeCatalog uses this relation to estimate the magnitude errors.
 The signal-to-noise ratio is calculated in different ways for clumps and 
objects (see @url{https://arxiv.org/abs/1505.01664, Akhlaghi and Ichikawa 
[2015]}), but this single equation can be used to estimate the measured 
magnitude error afterwards for any type of target.
 
-@node Completeness limit of each detection, Upper limit magnitude of each 
detection, Magnitude measurement error of each detection, Quantifying 
measurement limits
+@node Surface brightness error of each detection, Completeness limit of each 
detection, Magnitude measurement error of each detection, Quantifying 
measurement limits
+@subsubsection Surface brightness error of each detection
+
+@cindex Surface brightness error
+@cindex Error in surface brightness
+We can derive the error in measuring the surface brightness based on the 
surface brightness (SB) equation of @ref{Brightness flux magnitude} and the 
generic magnitude error (@mymath{\Delta{M}}) of @ref{Magnitude measurement 
error of each detection}.
+Let's set @mymath{A} to represent the area and @mymath{\Delta{A}} to represent 
the error in measuring the area.
+For more on @mymath{\Delta{A}}, see the description of 
@option{--spatialresolution} in @ref{MakeCatalog inputs and basic settings}.
+
+@dispmath{\Delta{(SB)} = \Delta{M} + \left|{-2.5\over 
ln(10)}\right|\times{\Delta{A}\over{A}}}
+
+In the surface brightness equation mentioned above, @mymath{A} is in units of 
arcsecond squared and the conversion between arcseconds to pixels is a 
multiplication factor.
+Therefore as long as @mymath{A} and @mymath{\Delta{A}} have the same units, it 
doesn't matter if they are in arcseconds or pixels.
+Since the measure of spatial resolution (or area error) is the FWHM of the PSF 
which is usually defined in terms of pixels, its more intuitive to use pixles 
for @mymath{A} and @mymath{\Delta{A}}.
+
+@node Completeness limit of each detection, Upper limit magnitude of each 
detection, Surface brightness error of each detection, Quantifying measurement 
limits
 @subsubsection Completeness limit of each detection
 @cindex Completeness
 As the surface brightness of the objects decreases, the ability to detect them 
will also decrease.
@@ -18507,9 +18524,18 @@ For the maximum value, see the description of 
@option{--maximum} column below.
 The value(s) of this option must be larger than 0 and smaller than 1 (they are 
a fraction).
 When only @option{--fracmaxarea1} or @option{--fracmaxsum1} is requested, one 
value must be given to this option, but if @option{--fracmaxarea2} or 
@option{--fracmaxsum2} are also requested, two values must be given to this 
option.
 The values can be written as simple floating point numbers, or as fractions, 
for example @code{0.25,0.75} and @code{0.25,3/4} are the same.
+
+@item --spatialresolution=FLT
+The error in measuring spatial properties (for example the area) in units of 
pixels.
+You can think of this as the FWHM of the dataset's PSF and is used in 
measurements like the error in surface brightness (@option{--sberror}, see 
@ref{MakeCatalog measurements}).
+Ideally, images are taken in the optimal Nyquist sampling @ref{Sampling 
theorem}, so the default value for this option is 2.
+But in practice real images my be over-sampled (usually ground-based images, 
where you will need to increase the default value) or undersampled (some 
space-based images, where you will need to decrease the default value).
 @end table
 
 
+
+
+
 @node Upper-limit settings, MakeCatalog measurements, MakeCatalog inputs and 
basic settings, Invoking astmkcatalog
 @subsubsection Upper-limit settings
 
@@ -18954,13 +18980,17 @@ The number of pixels that are equal to the maximum 
value of the labeled region (
 The surface brightness (in units of mag/arcsec@mymath{^2}) of the labeled 
region (objects or clumps).
 For more on the definition of the surface brightness, see @ref{Brightness flux 
magnitude}.
 
+@item --sberror
+Error in measuring the surface brightness (the @option{--surfacebrightness} 
column).
+This column will use the value given to @option{--spatialresolution} in the 
processing (in pixels).
+For more on @option{--spatialresolution}, see see @ref{MakeCatalog inputs and 
basic settings} and for the equation used to derive the surface brightness 
error, see @ref{Surface brightness error of each detection}.
+
 @item --areaxy
 @cindex IFU: Integral Field Unit
 @cindex Integral Field Unit
-Similar to @option{--area}, when the clump or object is projected onto the
-first two dimensions. This is only available for 3-dimensional
-datasets. When working with Integral Field Unit (IFU) datasets, this
-projection onto the first two dimensions would be a narrow-band image.
+Similar to @option{--area}, when the clump or object is projected onto the 
first two dimensions.
+This is only available for 3-dimensional datasets.
+When working with Integral Field Unit (IFU) datasets, this projection onto the 
first two dimensions would be a narrow-band image.
 
 @item --fwhm
 @cindex FWHM
@@ -19065,10 +19095,9 @@ Note that unlike @option{--area}, pixel values are 
completely ignored in this co
 For example, if a pixel value is blank, it won't be counted in 
@option{--area}, but will be counted here.
 
 @item --geoareaxy
-Similar to @option{--geoarea}, when the clump or object is projected onto
-the first two dimensions. This is only available for 3-dimensional
-datasets. When working with Integral Field Unit (IFU) datasets, this
-projection onto the first two dimensions would be a narrow-band image.
+Similar to @option{--geoarea}, when the clump or object is projected onto the 
first two dimensions.
+This is only available for 3-dimensional datasets.
+When working with Integral Field Unit (IFU) datasets, this projection onto the 
first two dimensions would be a narrow-band image.
 
 @item -A
 @itemx --semimajor