[gnuastro-commits] master 0289aee: MakeProfiles: new ability to build any custom radial profile

[Mohammad Akhlaghi]

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

    MakeProfiles: new ability to build any custom radial profile
    
    Until now, MakeProfiles could only build profiles with a pre-defined
    parametric function. But in some scenarios users need to build custom
    profiles and using the output of the distance profile may not be enough.
    
    With this commit, MakeProfiles now has a new 'custom' profile that comes
    with the new '--customtable' and '--customtablehdu' options. For such
    profiles the a 3-column table will be necessary that has the interval's
    minimum radius, maximum radius and value to use.
---
 NEWS                      |   7 ++++
 bin/mkprof/args.h         |  29 ++++++++++++-
 bin/mkprof/astmkprof.conf |   1 +
 bin/mkprof/main.h         |   5 +++
 bin/mkprof/oneprofile.c   |   8 ++++
 bin/mkprof/profiles.c     |  38 ++++++++++++++++-
 bin/mkprof/profiles.h     |   3 ++
 bin/mkprof/ui.c           | 102 ++++++++++++++++++++++++++++++++++++++++++++++
 bin/mkprof/ui.h           |   2 +
 doc/gnuastro.texi         |  53 ++++++++++++++++++++++++
 10 files changed, 246 insertions(+), 2 deletions(-)

diff --git a/NEWS b/NEWS
index 7c72aaf..b4d14e9 100644
--- a/NEWS
+++ b/NEWS
@@ -83,6 +83,13 @@ See the end of the file for license conditions.
      (used to simulate Poisson noise) should be interpretted as brightness,
      not magnitude.
 
+  MakeProfiles:
+   - It is now possible to make any custom radial profile with the 'custom'
+     profile (with code '8'). A table should be given to the new
+     '--customtable' option which will define each radial interval and the
+     value to use for that radial interval. See the description of
+     '--customtable' for more.
+
   Table:
    - New '--noblank' option will remove all rows in output table that have
      atleast one blank value in the specified columns. For example if
diff --git a/bin/mkprof/args.h b/bin/mkprof/args.h
index 4654615..73f4140 100644
--- a/bin/mkprof/args.h
+++ b/bin/mkprof/args.h
@@ -83,6 +83,32 @@ struct argp_option program_options[] =
       GAL_OPTIONS_NOT_SET,
       ui_parse_kernel
     },
+    {
+      "customtable",
+      UI_KEY_CUSTOMTABLE,
+      "STR",
+      0,
+      "File for 'custom' profile.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->customname,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
+    {
+      "customtablehdu",
+      UI_KEY_CUSTOMTABLEHDU,
+      "INT/STR",
+      0,
+      "HDU of table given to '--customtable'.",
+      GAL_OPTIONS_GROUP_INPUT,
+      &p->customhdu,
+      GAL_TYPE_STRING,
+      GAL_OPTIONS_RANGE_ANY,
+      GAL_OPTIONS_NOT_MANDATORY,
+      GAL_OPTIONS_NOT_SET
+    },
 
 
 
@@ -364,7 +390,8 @@ struct argp_option program_options[] =
       "STR/INT",
       0,
       "sersic (1), moffat (2), gaussian (3), point (4), "
-      "flat (5), circumference (6), distance (7).",
+      "flat (5), circumference (6), distance (7),"
+      "radial-table (8)",
       UI_GROUP_CATALOG,
       &p->fcol,
       GAL_TYPE_STRING,
diff --git a/bin/mkprof/astmkprof.conf b/bin/mkprof/astmkprof.conf
index a707f82..6c330e5 100644
--- a/bin/mkprof/astmkprof.conf
+++ b/bin/mkprof/astmkprof.conf
@@ -21,6 +21,7 @@
 
 #input
  backhdu                   1
+ customtablehdu            1
 
 # Output:
  mergedsize        1000,1000
diff --git a/bin/mkprof/main.h b/bin/mkprof/main.h
index 87c97b2..12ccef5 100644
--- a/bin/mkprof/main.h
+++ b/bin/mkprof/main.h
@@ -66,6 +66,7 @@ enum profile_types
   PROFILE_FLAT,                 /* Flat profile.               */
   PROFILE_CIRCUMFERENCE,        /* Circumference profile.      */
   PROFILE_DISTANCE,             /* Elliptical radius of pixel. */
+  PROFILE_CUSTOM,          /* Radial prof. in file/table. */
 
   PROFILE_MAXIMUM_CODE,         /* Just for a sanity check.    */
 };
@@ -114,6 +115,8 @@ struct mkprofparams
   char            *backname;  /* Name of background image file name.      */
   char             *catname;  /* Name of catalog of parameters.           */
   char             *backhdu;  /* HDU of background image.                 */
+  char          *customname;  /* Table to use for radial profile.         */
+  char           *customhdu;  /* HDU of table to use for radial profile.  */
   size_t             *dsize;  /* Size of the output image.                */
   uint8_t       clearcanvas;  /* Pixels in background image set to zero.  */
   gal_data_t        *kernel;  /* Parameters to define a kernel.           */
@@ -188,6 +191,8 @@ struct mkprofparams
   char           *wcsheader;  /* The WCS header information for main img. */
   int            wcsnkeyrec;  /* The number of keywords in the WCS header.*/
   char       *mergedimgname;  /* Name of merged image.                    */
+  gal_data_t        *custom;  /* Table containing custom values.          */
+  double   customregular[2];  /* Non-NaN if input table is regular.       */
   int                  nwcs;  /* for WCSLIB: no. coord. representations.  */
   struct wcsprm        *wcs;  /* WCS information for this dataset.        */
   size_t               ndim;  /* Number of dimensions (for 'nomerged').   */
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index b307744..a6fb713 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -694,6 +694,14 @@ oneprofile_set_prof_params(struct mkonthread *mkp)
 
 
 
+    case PROFILE_CUSTOM:
+      mkp->profile          = profiles_custom_table;
+      mkp->truncr           = tp ? p->t[id] : p->t[id]*p->r[id];
+      mkp->correction       = 0;
+      break;
+
+
+
     default:
       error(EXIT_FAILURE, 0, "%s: a bug! Please contact us so we can "
             "correct this problem. The profile code %u is not recognized.",
diff --git a/bin/mkprof/profiles.c b/bin/mkprof/profiles.c
index d2a00d9..b344181 100644
--- a/bin/mkprof/profiles.c
+++ b/bin/mkprof/profiles.c
@@ -41,7 +41,7 @@ along with Gnuastro. If not, see 
<http://www.gnu.org/licenses/>.
 /****************************************************************
  *****************         Profiles:         ********************
  ****************************************************************/
-/* The Gaussian function at a point. */
+/* The distance of this pixel. */
 double
 profiles_radial_distance(struct mkonthread *mkp)
 {
@@ -52,6 +52,42 @@ profiles_radial_distance(struct mkonthread *mkp)
 
 
 
+/* Read the values based on the distance from a table. */
+double
+profiles_custom_table(struct mkonthread *mkp)
+{
+  double out;
+  long i;  /* May become negative. */
+  double *reg=mkp->p->customregular;
+  double *min=mkp->p->custom->array;
+  double *max=mkp->p->custom->next->array;
+  double *value=mkp->p->custom->next->next->array;
+
+  /* If the table isn't regular ('reg[0]' isn't NaN), then we have to parse
+     over the whole table. However, if its regular, we can find the proper
+     value much more easily. */
+  if( isnan(reg[0]) )
+    {
+      out=0;
+      for(i=0;i<mkp->p->custom->size;++i)
+        if( mkp->r >= min[i] && mkp->r < max[i] )
+          { out=value[i]; break; }
+    }
+  else
+    {
+      i=(mkp->r - reg[0])/reg[1];
+      if(i<0 || i>mkp->p->custom->size) out=0;
+      else                                   out=value[i];
+    }
+
+  /* Return the output value. */
+  return isnan(out) ? 0 : out;
+}
+
+
+
+
+
 /* The integral of the Gaussian from -inf to +inf equals the square root of
    PI. So from zero to +inf it equals half of that.*/
 double
diff --git a/bin/mkprof/profiles.h b/bin/mkprof/profiles.h
index cdc6fa3..5f330eb 100644
--- a/bin/mkprof/profiles.h
+++ b/bin/mkprof/profiles.h
@@ -27,6 +27,9 @@ double
 profiles_radial_distance(struct mkonthread *mkp);
 
 double
+profiles_custom_table(struct mkonthread *mkp);
+
+double
 profiles_gaussian_total(double q);
 
 double
diff --git a/bin/mkprof/ui.c b/bin/mkprof/ui.c
index 823d289..5cf63f0 100644
--- a/bin/mkprof/ui.c
+++ b/bin/mkprof/ui.c
@@ -124,6 +124,9 @@ ui_profile_name_read(char *string, size_t row)
   else if ( !strcmp("distance", string) )
     return PROFILE_DISTANCE;
 
+  else if ( !strcmp("custom", string) )
+    return PROFILE_CUSTOM;
+
   else if ( !strcmp(GAL_BLANK_STRING, string) )
     error(EXIT_FAILURE, 0, "atleast one profile function is blank");
 
@@ -187,6 +190,9 @@ ui_initialize_options(struct mkprofparams *p,
   cp->numthreads         = gal_threads_number();
   cp->coptions           = gal_commonopts_options;
 
+  p->customregular[0]    = NAN;
+  p->customregular[1]    = NAN;
+
   /* Default program parameters. */
   p->zeropoint           = NAN;
   p->cp.type             = GAL_TYPE_FLOAT32;
@@ -1093,6 +1099,7 @@ ui_read_cols_3d(struct mkprofparams *p)
 static void
 ui_prepare_columns(struct mkprofparams *p)
 {
+  size_t i;
   double *karr;
   float r, n, t, q2;
 
@@ -1191,6 +1198,21 @@ ui_prepare_columns(struct mkprofparams *p)
                 __func__, PACKAGE_BUGREPORT, p->ndim);
         }
     }
+
+  /* If a custom profile is requested, make sure that a custom file is
+     given. */
+  for(i=0;i<p->num;++i)
+    if(p->f[i]==PROFILE_CUSTOM)
+      {
+        if(p->customname==NULL)
+          error(EXIT_FAILURE, 0, "at least one custom profile requested "
+                "(first occurance in row %zu), but no file/table was given "
+                "to the '--customtable' option. See the description of "
+                "this '--customtable' for more information on the "
+                "desired format", i+1);
+        break;
+      }
+
 }
 
 
@@ -1637,6 +1659,83 @@ ui_make_log(struct mkprofparams *p)
 
 
 
+/* Read the input radial table. */
+static void
+ui_read_custom_table(struct mkprofparams *p)
+{
+  size_t i;
+  double diff;
+  int isregular;
+  gal_data_t *cols;
+  double *min, *max;
+
+  /* Read the input radial table. */
+  cols=gal_table_read(p->customname, p->customhdu,
+                      NULL, NULL, p->cp.searchin, p->cp.ignorecase,
+                      p->cp.minmapsize, p->cp.quietmmap, NULL);
+
+  /* Make sure the table only has three columns. */
+  if(gal_list_data_number(cols) != 3 )
+    error(EXIT_FAILURE, 0, "%s: has %zu columns, but it should only have "
+          "three columns. Column1: the radial interval's lower value. "
+          "Column 2: the radial interval's higher value. Column 3: the value "
+          "to use for pixels within that radius interval",
+          gal_fits_name_save_as_string(p->customname, p->customhdu),
+          gal_list_data_number(cols));
+
+  /* Make sure none of the three columns are string type. */
+  if( cols->type==GAL_TYPE_STRING
+      || cols->next->type==GAL_TYPE_STRING
+      || cols->next->next->type==GAL_TYPE_STRING )
+    error(EXIT_FAILURE, 0, "%s: the columns should only have numeric "
+          "data types", gal_fits_name_save_as_string(p->customname,
+                                                     p->customhdu));
+
+  /* Fill the final table as a double type. */
+  p->custom=gal_data_copy_to_new_type(cols, GAL_TYPE_FLOAT64);
+  p->custom->next=gal_data_copy_to_new_type(cols->next,
+                                                 GAL_TYPE_FLOAT64);
+  p->custom->next->next=gal_data_copy_to_new_type(cols->next->next,
+                                                       GAL_TYPE_FLOAT64);
+
+  /* Make sure the first column values are smaller than the second column's
+     values. */
+  min=p->custom->array;
+  max=p->custom->next->array;
+  for(i=0;i<p->custom->size;++i)
+    if(min[i]>=max[i])
+      error(EXIT_FAILURE, 0, "%s: the first column of row %zu (with value "
+            "%g) is larger or equal to the second column (with value %g). "
+            "However, the first column is the lower-limit of the radial "
+            "interval and the second column is the upper-limit. So the "
+            "first column must have a lower value",
+            gal_fits_name_save_as_string(p->customname,
+                                         p->customhdu), i+1,
+            min[i], max[i]);
+
+  /* Check if the input table is regular and sorted (which can greatly
+     speed up its usage). */
+  isregular=1;
+  diff=max[0]-min[0];
+  for(i=1;i<p->custom->size;++i)
+    if( min[i]<min[i-1]
+        || min[i] != max[i-1]
+        || max[i]-min[i] != diff )
+      isregular=0;
+  if(isregular)
+    {
+      p->customregular[0]=min[0];
+      p->customregular[1]=diff;
+    }
+
+  /* Clean up. */
+  gal_list_data_free(cols);
+}
+
+
+
+
+
 static void
 ui_read_ndim(struct mkprofparams *p)
 {
@@ -1721,6 +1820,9 @@ ui_preparations(struct mkprofparams *p)
      build the profiles). */
   ui_prepare_columns(p);
 
+  /* Read the radial table. */
+  if(p->customname) ui_read_custom_table(p);
+
   /* If the kernel option was given, some parameters need to be
      over-written: */
   if(p->kernel)
diff --git a/bin/mkprof/ui.h b/bin/mkprof/ui.h
index fb873fc..35283dd 100644
--- a/bin/mkprof/ui.h
+++ b/bin/mkprof/ui.h
@@ -93,6 +93,8 @@ enum option_keys_enum
   UI_KEY_PC,
   UI_KEY_CUNIT,
   UI_KEY_CTYPE,
+  UI_KEY_CUSTOMTABLE,
+  UI_KEY_CUSTOMTABLEHDU,
 };
 
 
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index a8715d9..f7c674e 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -17607,6 +17607,10 @@ For this profile, the value in the magnitude column 
(@option{--mcol}) will be ig
 
 You can use this for checks or as a first approximation to define your own 
higher-level radial function.
 In the latter case, just note that the central values are going to be 
incorrect (see @ref{Sampling from a function}).
+
+@item
+Custom profile with `@code{custom}' or `@code{8}'.
+The values to use for each radial interval should be in the table given to 
@option{--customtable}. For more, see @ref{MakeProfiles profile settings}.
 @end itemize
 
 @item --rcol=STR/INT
@@ -17734,6 +17738,55 @@ Particularly in sharper profiles, the flux in the peak 
pixel is strongly decreas
 Also note that in such cases, besides de-convolution, you will have to set 
@option{--oversample=1} otherwise after resampling your profile with Warp (see 
@ref{Warp}), the peak flux will be different.
 @end cartouche
 
+@item --customtable STR
+The filename of the table to use in the custom profiles (see description of 
@option{--fcol} in @ref{MakeProfiles catalog}.
+This can be a plain-text table, or FITS table, see @ref{Tables}, if its a FITS 
table, you can use @option{--customtablehdu} to specify which HDU should be 
used (described below).
+
+A custom profile can have any value you want for a given radial profile.
+Each interval is defined by its minimum (inclusive) and maximum (exclusive) 
radius, when a pixel falls within this radius the value specified for that 
interval will be used.
+If a pixel is not in the given intervals, a value of 0 will be used for it.
+
+The table should have 3 columns as shown below.
+If the intervals are contiguous (the maximum value of the previous interval is 
equal to the minimum value of an interval) and the intervals all have the same 
size (difference between minimum and maximum values) the creation of these 
profiles will be fast.
+However, if the intervals are not sorted and contiguous, Makeprofiles will 
parse the intervals from the top of the table and use the first interval that 
contains the pixel center.
+
+@table @asis
+@item Column 1:
+The interval's minimum radius.
+@item Column 2:
+The interval's maximum radius.
+@item Column 3:
+The value to be used for pixels within the given interval.
+@end table
+
+For example let's assume you have the radial profile below in a file called 
@file{radial.txt}.
+The first column is the larger interval radius and the second column is the 
value in that interval:
+
+@example
+1    100
+2    90
+3    50
+4    10
+5    2
+6    0.1
+7    0.05
+@end example
+
+@noindent
+You can construct the table to give to @option{--customtable} with either of 
the commands below: the first one with Gnuastro's @ref{Column arithmetic} which 
can also work on FITS tables, and the second one with an AWK command that only 
works on plain-text tables..
+
+@example
+asttable radial.fits -c'arith $1 1 -' -c1,2 -ocustom.fits
+awk '@{print $1-1, $1, $2@}' radial.txt > custom.txt
+@end example
+
+@noindent
+In case the intervals are different from 1 (for example 0.5), change them 
respectively: for Gnuastro's table change @code{$1 1 -} to @code{$1 0.5 -} and 
for AWK change  @code{$1-1} to @code{$1-0.5}.
+
+
+@item --customtablehdu INT/STR
+The HDU/extension in the FITS file given to @option{--customtable}.
+
 @item -X INT,INT
 @itemx --shift=INT,INT
 Shift all the profiles and enlarge the image along each dimension.