[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[gnuastro-commits] master ea015c3b 2/2: Table: new sorted-to-interval op
|
From: |
Mohammad Akhlaghi |
|
Subject: |
[gnuastro-commits] master ea015c3b 2/2: Table: new sorted-to-interval operator in column arithmetic |
|
Date: |
Sat, 10 Dec 2022 19:49:08 -0500 (EST) |
branch: master
commit ea015c3bae5720ba3912c4c41b2f07b002f2172d
Author: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>
Table: new sorted-to-interval operator in column arithmetic
Until now, the radial profile script had the '--radinterval' option to
generate an input to the '--customtable' option of MakeProfile. It had
multiple problems:
1. When the output was generated in this format, the raw profile was no
longer created.
2. It was hard to do some operations on the radial profile (like merging
multiple profiles) and generate an input to the '--customtable' option
of MakeProfiles.
With this commit, that option has been removed from the radial profile
script and instead, a new column arithmetic operator has been added that is
called 'sorted-to-interval'. This allows users to easily generate an input
to the '--customtable' option of MakeProfiles from any radial profile
easily. This fixes both the problems above and is much more modular and
generic.
---
NEWS | 17 +++---
bin/script/radial-profile.in | 36 +----------
bin/table/arithmetic.c | 77 ++++++++++++++++++++++++
bin/table/arithmetic.h | 1 +
doc/gnuastro.texi | 140 +++++++++++++++++++++++++------------------
5 files changed, 173 insertions(+), 98 deletions(-)
diff --git a/NEWS b/NEWS
index a43f3990..12b26728 100644
--- a/NEWS
+++ b/NEWS
@@ -77,19 +77,22 @@ See the end of the file for license conditions.
Statistics:
--outliernumngb: see description of same option in NoiseChisel.
+ Table:
+ - New column arithmetic operator:
+ - sorted-to-interval: return two columns from a single (sorted) input,
+ containing the minimum and maximum values of an interval. The
+ intervals are defined such that the intervals of every row fully
+ cover the full range of values in the given column. This can be
+ useful in scenarios where you need to build the inputs of the
+ '--customtable' feature in MakeProfiles (to build a 2D image a
+ custom profile).
+
astscript-radial-profile:
--precision: sample the radial profile at precisions less than one
pixel. This is useful when you need to sample the profile within the
central few pixels more accurately than integer pixels. See the
documentation of this option in the book for a complete explanation
with examples.
- --radinterval: output will contain two radius columns (showing the
- radial interval of each measurement). Without this option, the output
- only contains a single radius column (showing the central radius of
- the measurement in each row). This greatly simplifies scenarios where
- you want to build a 2D image of the profile, see the description of
- this option in the book for a more complete description and an example
- usage.
Library:
- GAL_ARITHMETIC_OP_SWAP: swap the top two operands.
diff --git a/bin/script/radial-profile.in b/bin/script/radial-profile.in
index d3e0ce8b..2a12bf64 100644
--- a/bin/script/radial-profile.in
+++ b/bin/script/radial-profile.in
@@ -56,7 +56,6 @@ zeropoint=""
sigmaclip=""
measuretmp=""
oversample=""
-radinterval=0
undersample=""
positionangle=0
zeroisnotblank=""
@@ -129,8 +128,6 @@ $scriptname options:
-i, --instd=FLT/STR Sky standard deviation per pixel, as a single
number or as the filename.
-d, --stdhdu=STR HDU/extension of the sky standard deviation image.
- -I, --radinterval Format output for the '--customtable' feature
- of MakeProfile.
Operating mode:
-?, --help Print this help list.
@@ -275,8 +272,6 @@ do
-d|--stdhdu) stdhdu="$2"; check_v "$1"
"$stdhdu"; shift;shift;;
-d=*|--stdhdu=*) stdhdu="${1#*=}"; check_v "$1"
"$stdhdu"; shift;;
-d*) stdhdu=$(echo "$1" | sed -e's/-d//'); check_v "$1"
"$stdhdu"; shift;;
- -F|--radinterval) radinterval=1; shift;;
- -F*|--radinterval=*) on_off_option_error --radinterval -I;;
# Non-operating options.
-q|--quiet) quiet="--quiet"; shift;;
@@ -819,33 +814,8 @@ outraw=$tmpdir/out-raw.fits
restcols=$(astfits $cat -h1 \
| awk '/^TFIELDS/{for(i=2;i<=$3;++i) \
printf "--catcolumns=%d ", i}')
-if [ x"$radinterval" = x0 ]; then
- checkcol=RADIUS
- asttable $radraw --catcolumnfile=$cat $restcols --output=$outraw \
- --colmetadata=1,RADIUS,pix,"Radial distance"
-else
- # Since the larger interval value will be outside of the '--rmax' by
- # default, we'll use the smaller interval value for the final check in
- # the next step.
- checkcol=RADIUS_MIN
-
- # Build the radial intervals from the radii. Note that order important
- # in the AWK commands below.
- radintervalraw=$tmpdir/radii-interval.fits
- asttable $radraw -c1 \
- | awk 'NR==2 {print 0, ($1-prev)/2} \
- NR>2 {print prev-(prev-pprev)/2, prev+($1-prev)/2} \
- {pprev=prev; prev=$1} \
- END{print prev-(prev-pprev)/2, prev+(prev-pprev)/2}' \
- | asttable -c'arith $1 float32','arith $2 float32' \
- --colmetadata=1,RADIUS_MIN,pix,"Minimum radius of this
interval" \
- --colmetadata=2,RADIUS_MAX,pix,"Maximum radius of this
interval" \
- -o$radintervalraw
-
- # Add the rest of the catalog columns to the interval.
- asttable $radintervalraw --catcolumnfile=$cat $restcols --output=$outraw
- rm $radintervalraw
-fi
+asttable $radraw --catcolumnfile=$cat $restcols --output=$outraw \
+ --colmetadata=1,RADIUS,pix,"Radial distance"
@@ -859,7 +829,7 @@ fi
# MakeCatalog). Therefore, once the final radii are set (accounting for
# over/under sampling), we should make sure that the final output doesn't
# contain radii larger than what the user asked for.
-asttable $outraw --range=$checkcol,0,$rmax --output=$output
+asttable $outraw --range=RADIUS,0,$rmax --output=$output
diff --git a/bin/table/arithmetic.c b/bin/table/arithmetic.c
index afe60ce0..ac9f0636 100644
--- a/bin/table/arithmetic.c
+++ b/bin/table/arithmetic.c
@@ -31,6 +31,7 @@ along with Gnuastro. If not, see
<http://www.gnu.org/licenses/>.
#include <gnuastro/wcs.h>
#include <gnuastro/type.h>
#include <gnuastro/pointer.h>
+#include <gnuastro/statistics.h>
#include <gnuastro-internal/checkset.h>
#include <gnuastro-internal/arithmetic-set.h>
@@ -139,6 +140,7 @@ arithmetic_operator_name(int operator)
case ARITHMETIC_TABLE_OP_DISTANCEFLAT: out="distance-flat"; break;
case ARITHMETIC_TABLE_OP_DATETOMILLISEC: out="date-to-millisec"; break;
case ARITHMETIC_TABLE_OP_DISTANCEONSPHERE: out="distance-on-sphere";
break;
+ case ARITHMETIC_TABLE_OP_SORTEDTOINTERVAL: out="sorted-to-interval";
break;
default:
error(EXIT_FAILURE, 0, "%s: a bug! Please contact us at %s to fix "
"the problem. %d is not a recognized operator code", __func__,
@@ -203,6 +205,8 @@ arithmetic_set_operator(struct tableparams *p, char *string,
{ op=ARITHMETIC_TABLE_OP_DISTANCEFLAT; *num_operands=0; }
else if( !strcmp(string, "distance-on-sphere"))
{ op=ARITHMETIC_TABLE_OP_DISTANCEONSPHERE; *num_operands=0; }
+ else if( !strcmp(string, "sorted-to-interval"))
+ { op=ARITHMETIC_TABLE_OP_SORTEDTOINTERVAL; *num_operands=0; }
else
{ op=GAL_ARITHMETIC_OP_INVALID; *num_operands=GAL_BLANK_INT; }
}
@@ -834,6 +838,75 @@ arithmetic_datetosec(struct tableparams *p, gal_data_t
**stack,
+/* Convert the ISO date format to seconds since Unix time. */
+static void
+arithmetic_sortedtointerval(struct tableparams *p, gal_data_t **stack,
+ int operator)
+{
+ int f32out;
+ size_t i, size;
+ double *in, *min, *max;
+ gal_data_t *min_d, *max_d;
+
+ /* Input dataset. */
+ gal_data_t *in_d=arithmetic_stack_pop(stack, operator, NULL);
+
+ /* Make sure the input is sorted and is numeric */
+ if(in_d->type==GAL_TYPE_STRING)
+ error(EXIT_FAILURE, 0, "the operand given to 'sortedtointerval' "
+ "should have a numeric data type");
+ if(gal_statistics_is_sorted(in_d, 1)==0)
+ error(EXIT_FAILURE, 0, "the operand given to 'sortedtointerval' "
+ "should be sorted");
+
+ /* Convert the input to double precision for internal processing. In the
+ end we'll convert the outputs to float32 unless the input was
+ originally float64. */
+ f32out = in_d->type!=GAL_TYPE_FLOAT64;
+ in_d=gal_data_copy_to_new_type_free(in_d, GAL_TYPE_FLOAT64);
+
+ /* Allocate the output columns (minimum and maximum of each interval). */
+ min_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, in_d->dsize, NULL, 0,
+ in_d->minmapsize, in_d->quietmmap,
+ NULL, NULL, NULL);
+ max_d=gal_data_alloc(NULL, GAL_TYPE_FLOAT64, 1, in_d->dsize, NULL, 0,
+ in_d->minmapsize, in_d->quietmmap,
+ NULL, NULL, NULL);
+
+ /* Fill the minimum and maximum intervals. Just note that the minimum
+ value of the first interval and maximum of the last element need
+ special attention (we'll extrapolate from the after/before elements to
+ preserve constant spacing in case that is the case). */
+ in=in_d->array;
+ size=in_d->size;
+ min=min_d->array;
+ max=max_d->array;
+ min[0] = in[0] - (in[1] - in[0] )/2;
+ max[0] = in[0] + (in[1] - in[0] )/2;
+ min[size-1] = in[size-1] - (in[size-1] - in[size-2])/2;
+ max[size-1] = in[size-1] + (in[size-1] - in[size-2])/2;
+ for(i=1;i<size-1;++i)
+ {
+ min[i] = max[i-1];
+ max[i] = in[i] + (in[i+1]-in[i])/2;
+ }
+
+ /* Convert the output to float32, unless the input was float64. */
+ if(f32out)
+ {
+ min_d=gal_data_copy_to_new_type_free(min_d, GAL_TYPE_FLOAT32);
+ max_d=gal_data_copy_to_new_type_free(max_d, GAL_TYPE_FLOAT32);
+ }
+
+ /* Clean up and put the resulting calculation back on the stack. */
+ max_d->next=min_d;
+ gal_data_free(in_d);
+ gal_list_data_add(stack, max_d);
+}
+
+
+
+
@@ -941,6 +1014,10 @@ arithmetic_operator_run(struct tableparams *p,
arithmetic_distance(p, stack, token->operator);
break;
+ case ARITHMETIC_TABLE_OP_SORTEDTOINTERVAL:
+ arithmetic_sortedtointerval(p, stack, token->operator);
+ break;
+
case ARITHMETIC_TABLE_OP_SET:
gal_arithmetic_set_name(setprm, token->name_def);
break;
diff --git a/bin/table/arithmetic.h b/bin/table/arithmetic.h
index 091241b5..47a121b2 100644
--- a/bin/table/arithmetic.h
+++ b/bin/table/arithmetic.h
@@ -42,6 +42,7 @@ enum arithmetic_operators
ARITHMETIC_TABLE_OP_DISTANCEFLAT,
ARITHMETIC_TABLE_OP_DATETOMILLISEC,
ARITHMETIC_TABLE_OP_DISTANCEONSPHERE,
+ ARITHMETIC_TABLE_OP_SORTEDTOINTERVAL,
};
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 0fdbc12b..7f89b1fb 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -14303,6 +14303,39 @@ The returned value is not a floating point type
because for large numbers, float
Other than the units of the output, this operator behaves similarly to
@code{date-to-sec}.
See the description of that operator for an example.
+
+@item sorted-to-interval
+Given a single column (which must be already sorted and have a numeric data
type), return two columns: the first returned column is the minimum and the
second returned column is the maximum value of the interval of each row row.
+The maximum of each row is equal to the minimum of the previous row; thus
creating a contiguous interval coverage of the input column's range in all rows.
+
+The minimum value of the first row and maximum of the last row will be
smaller/larger than the respective row of the input (based on the distance to
the next/previous element).
+This is done to ensure that if your input has a fixed interval length between
all elements, the first and last intervals also have that fixed length.
+
+For example, with the command below, we'll use this operator on a hypothetical
radial profile.
+Note how the intervals are contiguous even though the radius values are not
equally distant (if the row with a radius of 2.5 didn't exist, the intervals
would all be the same length).
+For another example of the usage of this operator, see the example in the
description of @option{--customtable} in @ref{MakeProfiles profile settings}.
+
+@example
+$ cat radial-profile.txt
+# Column 1: RADIUS [pix,f32,] Distance to center in pixels.
+# Column 2: MEAN [ADU,f32,] Mean value at that radius.
+0 100
+1 40
+2 30
+2.5 25
+3 20
+
+$ asttable radial-profile.txt --txtf32f=fixed --txtf32p=3 \
+ -c'arith RADIUS sorted-to-interval',MEAN
+-0.500 0.500 100.000
+0.500 1.500 40.000
+1.500 2.250 30.000
+2.250 2.750 25.000
+2.750 3.250 20.000
+@end example
+
+Such intervals can be useful in scenarios like generating the input to
@option{--customtable} in MakeProfiles (see @ref{MakeProfiles profile
settings}) from a radial profile (see @ref{Generate radial profile}).
+
@end table
@node Operation precedence in Table, Invoking asttable, Column arithmetic,
Table
@@ -26578,35 +26611,59 @@ The interval's maximum radius.
The value to be used for pixels within the given interval (including
NaN/blank).
@end table
-If you are using Gnuastro's radial profile script to build the profile, you
can use its @option{--radinterval} option to make the properly formatted input
for a custom table in MakeProfiles.
-For more, see the description of @option{--radinterval}, and its example usage
in @ref{Invoking astscript-radial-profile}.
-
-For example, let's assume you have the radial profile below in a file called
@file{radial.txt}, and that the first column is the @emph{larger} interval
radius (in units of pixels) 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
-In this case, you can construct the table to give to @option{--customtable}
with the command below (using Gnuastro's @ref{Column arithmetic}).
-
-@example
-asttable radial.txt -c'arith $1 1 -' -c1,2 -ocustom.fits
-@end example
-
-@noindent
-In case the intervals are different from 1 (for example, 0.5), change the
@code{$1 1 -} to @code{$1 0.5 -}.
-On a side-note, Gnuastro has features to extract the radial profile of an
object from the image, see @ref{Generate radial profile}.
+Gnuastro's column arithmetic in the Table program has the
@code{sorted-to-interval} operator that will generate the first two columns
from a single column (your radial profile).
+See the description of that operator in @ref{Column arithmetic} and the
example below.
By default, once a 2D image is constructed for the radial profile, it will be
scaled such that its total magnitude corresponds to the value in the magnitude
column (@option{--mcol}) of the main input catalog.
If you want to disable the scaling and use the raw values in your custom
profile (in other words: you want to ignore the magnitude column) you need to
call @option{--mcolnocustprof} (see above).
+In the example below, we'll start with a certain radial profile, and use this
option to build its 2D representation in an image (recall that you can build
radial profiles with @ref{Generate radial profile}).
+But first, we will need to use the @code{sorted-to-interval} to build the
necessary input format (see @ref{Column arithmetic}).
+
+@example
+$ cat radial.txt
+# Column 1: RADIUS [pix ,f32,] Radial distance
+# Column 2: MEAN [input-units,f32,] Mean of values.
+0.0 1.00000
+1.0 0.50184
+1.4 0.37121
+2.0 0.26414
+2.2 0.23427
+2.8 0.17868
+3.0 0.16627
+3.1 0.15567
+3.6 0.13132
+4.0 0.11404
+
+## Convert the radius in each row to an interval
+$ asttable radial.txt --output=interval.fits \
+ -c'arith RADIUS sorted-to-interval',MEAN
+
+## Inspect the table containing intervals
+$ asttable interval.fits -ffixed
+-0.500000 0.500000 1.000000
+0.500000 1.200000 0.501840
+1.200000 1.700000 0.371210
+1.700000 2.100000 0.264140
+2.100000 2.500000 0.234270
+2.500000 2.900000 0.178680
+2.900000 3.050000 0.166270
+3.050000 3.350000 0.155670
+3.350000 3.800000 0.131320
+3.800000 4.200000 0.114040
+
+## Build the 2D image of the profile from the interval.
+$ echo "1 7 7 8 10 2.5 0 1 1 2" \
+ | astmkprof --mergedsize=13,13 --oversample=1 \
+ --customtable=interval.fits \
+ --output=image.fits
+
+## View the created FITS image.
+$ astscript-fits-view image.fits --ds9scale=minmax
+@end example
+
+Recall that if you want your image pixels to have the same values as the
@code{MEAN} column in your profile, you should run MakeProfiles with
@option{--mcolnocustprof}.
+
@item --customtablehdu INT/STR
The HDU/extension in the FITS file given to @option{--customtable}.
@@ -28514,39 +28571,6 @@ You can disable the deletion of the temporary
directory with the @option{--keept
Do Not delete the temporary directory (see description of @option{--tmpdir}
above).
This option is useful for debugging.
For example, to check that the profiles generated for obtaining the radial
profile have the desired center, shape and orientation.
-
-@item -I
-@itemx --radinterval
-Output will contain the radial interval of each measurement (in two radius
columns: minimum and maximum radii of each measurement/row).
-Without this option, only the central radius value is used (so the default
output only has a single radius column).
-
-With this option, the output can be directly fed to the @option{--customtable}
option of MakeProfiles (see @ref{MakeProfiles profile settings}).
-In this format, the radial interval of each row is needed, not just the
central value.
-Therefore when this option is called, the radial profile script output will
have two radius columns showing the minimum (inclusive) and maximum (exclusive)
radii that the measurement was done on.
-
-In the example below, we are first measuring the radial profile of a
hypothetical object (at RA,Dec of 1.234,6.789), while activating the
@option{--radinterval} feature.
-The resulting radial profile is then given to MakeProfiles.
-Note that this is a minimal call to MakeProfiles just to show this feature,
you can heavily customize MakeProfiles (for example build the profile over a
background image, or have many rows with many different types of other profiles
and etc.).
-For more on MakeProfiles see @ref{MakeProfiles}.
-
-@example
-$ astscript-radial-profile image.fits --center=1.234,6.789 \
- --mode=wcs --rmax=10 --precision=1 --radinterval \
- --measure=mean --output=rad.fits
-
-$ echo "1 50 50 8 20 0 0 1 1 1" \
- | astmkprof --oversample=1 --mergedsize=99,99 \
- --customtable=rad.fits --mcolnocustprof \
- --output=prof2d.fits
-
-$ astscript-fits-view prof2d.fits
-@end example
-
-@cartouche
-@noindent
-@strong{Use a single measurement to feed to MakeProfiles:} The table that is
given to the @option{--customtable} option must have three columns.
-Therefore if you plan to feed the output to MakeProfiles, only use a single
measurement for the radial profile script.
-@end cartouche
@end table