[gnuastro-commits] master 5efb92b9: radial-profile: larger first azimuth angle will use outer part

[Mohammad Akhlaghi]

branch: master
commit 5efb92b9cc3d466d24bc05f6bca697448d4e5615
Author: Samane Raji <samaneraji@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    radial-profile: larger first azimuth angle will use outer part
    
    Until now, we were just able to generate an azimuthal profile with an
    initial angle smaller than the final angle. If we needed to have a
    symmetric wedge along the semi-major axis, in which the initial angle is
    bigger than the final angle for exmaple from 350 to 10, the radial profile
    script would crash.
    
    With this commit, this bug is fixed by adding a conditional statement in
    the call to 'astarithmetic'. In effect, this has also increased the
    functionality of the radial profile script: because the user can now ignore
    parts of the radial profile and use the reset of the azimuthal angles.
    
    In the process, the following issues have also been corrected:
    
      - When the azimuthal angles are used, the central pixel is effectively
        ignored because it is no longer a "hole" and gets connected to the
        regions outside the radial range. So the step to fill the central hole
        and increment the radial distances (that will be corrected in the final
        profile) is now done before the azimuthal correction.
    
     - The azimuthal profile (to select pixels in the desired azimuthal range)
       is now generated in float32 numerical type, not int16! This allows more
       accurate selection of the azimuthal range, especially in larger
       profiles.
    
    This fixes bug #62721
    
    The documentation and the parts relating to the central pixel have been
    added by Mohammad Akhlaghi.
---
 NEWS                         | 13 +++++++++
 bin/script/radial-profile.in | 67 +++++++++++++++++++++++++-------------------
 doc/gnuastro.texi            | 34 ++++++++++++++++++----
 3 files changed, 79 insertions(+), 35 deletions(-)

diff --git a/NEWS b/NEWS
index 381a4259..142a13bf 100644
--- a/NEWS
+++ b/NEWS
@@ -172,6 +172,15 @@ See the end of the file for license conditions.
      the old name was very confusing and would not directly convey this
      behavior.
 
+  Radial profile script:
+   --azimuth: if the first azimuthal angle is larger than the second, the
+     outer region between the two angles will be used. For example
+     '--azimuth=80,40' will use the full azimuthal range except for the
+     region with azimuthal angles between 40 and 80 degrees. This can be
+     useful to ignore parts of the profile. See the example in the
+     description of this option on how to visually check the pixels used
+     for the profile.
+
 ** Bugs fixed
   bug #62216: MakeProfiles crash when a 3D cube is requested and input
               catalog is from a pipe. Found by Irene Pintos Castro.
@@ -230,6 +239,10 @@ See the end of the file for license conditions.
               tables. Found by Hilderic Browne.
   bug #62720: Table not reading string columns that are shorter than the
               defined width in the metadata.
+  bug #62721: Radial profile script can't deal with angles around azimuth
+              zero (for example '--azimuth=355,5' to get only the azimuthal
+              range of 10 degrees around the major axis). Found and fixed
+              by Samane Raji.
   bug #62784: Conversion of sexagesimal RA or Dec ignores sign when first
               digit is zero (for example '-00d12m34'). Found by Manuel
               Sánchez-Benavente.
diff --git a/bin/script/radial-profile.in b/bin/script/radial-profile.in
index 0f783775..e98d56a5 100644
--- a/bin/script/radial-profile.in
+++ b/bin/script/radial-profile.in
@@ -568,12 +568,21 @@ fi
 # axisratio     -- axis ratio.
 # rmax          -- magnitude of the profile within the truncation radius 
(rmax).
 # 1             -- Truncation in radius unit.
+#
+# After making the apertures image, we will set the central pixel to have
+# the radius 1 (to avoid missing the central pixel in later phases.
 radialaperturesbase=radial-raw.fits
 radialapertures=$tmpdir/$radialaperturesbase
+radialaperturesholed=$tmpdir/radial-raw-with-hole.fits
 echo "1 $xcenter $ycenter 7 $rmax 1 $positionangle $axisratio 1 1" \
      | astmkprof --background=$values --backhdu=1 --mforflatpix \
                  --mode=img --clearcanvas --type=int16 $quiet \
-                 --circumwidth=1 --replace --output=$radialapertures
+                 --circumwidth=1 --replace --output=$radialaperturesholed
+astarithmetic $radialaperturesholed set-i \
+              i 0 ne 1 fill-holes set-good \
+              i good i 1 + where \
+              $quiet --output $radialapertures
+rm $radialaperturesholed
 
 
 
@@ -601,19 +610,34 @@ if [ x"$azimuth" != x ]; then
     azimuthapertures=$tmpdir/$azimuthaperturesbase
     echo "1 $xcenter $ycenter 9 $rmax 1 $positionangle $axisratio 1 1" \
          | astmkprof --background=$values --backhdu=1 --mforflatpix \
-                     --mode=img --clearcanvas --type=int16 $quiet \
+                     --mode=img --clearcanvas --type=float32 $quiet \
                      --circumwidth=1 --replace --output=$azimuthapertures
 
     # From the azimuthal aperture image, consider only that portion that
-    # are in between the two specified angles. Set the rest of the pixels
-    # to zero values (on the radial apertures image)
-    astarithmetic $radialapertures -h1 set-radial \
+    # are in between (or outside) the two specified angles. Set the rest of
+    # the pixels to zero values (on the radial apertures image).
+    #
+    # There are two ways that the angle range can be defined:
+    #
+    #   - The first is smaller than the second (for example
+    #     '--azimuth=10,20'). In this case, we will use an 'and' between
+    #     the regions of each angle so the user gets the region between the
+    #     two angles.
+    #
+    #   - The first is smaller than the second (for example
+    #     '--azimuth=355,5') In this case, we assume that the user wants an
+    #     azimuthal range outside the two angles. In the example above, its
+    #     the 10 degrees around the azimuthal angle 0. For this case, we
+    #     should use 'or' between the two regions.
+    condition=$(echo $azimuth_ini $azimuth_fin \
+                     | awk ' {if ($1<$2) print "and"; \
+                              else       print "or"}')
+    astarithmetic --output=$aperturesraw $quiet \
+                  $radialapertures  -h1 set-radial \
                   $azimuthapertures -h1 set-azimuth \
-                  azimuth $azimuth_ini gt \
-                  azimuth $azimuth_fin lt \
-                  and int16 set-arc \
-                  radial arc 1 ne 0 where \
-                  $quiet --output $aperturesraw
+                  azimuth $azimuth_ini ge \
+                  azimuth $azimuth_fin le $condition set-arc \
+                  radial arc 1 uint8 ne 0 uint8 where
 else
     cd $tmpdir; ln -fs $radialaperturesbase $aperturesrawbase; cd $curdir
 fi
@@ -629,34 +653,19 @@ fi
 # number of pixels on each aperture (at each radius) will be larger. Most
 # of times this option is good to average over a larger number of pixels
 # and increase the signal-to-noise ratio of the measurement.
-aperturesrebinnedbase=apertures-rebin.fits
-aperturesrebinned=$tmpdir/$aperturesrebinnedbase
+aperturesbase=apertures.fits
+apertures=$tmpdir/$aperturesbase
 if [ x"$undersample" != x ]; then
     astarithmetic $aperturesraw $undersample / \
-                  $quiet --output $aperturesrebinned
+                  $quiet --output $apertures
 else
-    cd $tmpdir; ln -fs $aperturesrawbase $aperturesrebinnedbase; cd $curdir
+    cd $tmpdir; ln -fs $aperturesrawbase $aperturesbase; cd $curdir
 fi
 
 
 
 
 
-# Fill the central pixel(s)
-# -------------------------
-#
-# The central pixel(s) have a distance of 0! So we need to add a single
-# value to all the profile pixels (but keep the outer parts at 0).
-apertures=$tmpdir/apertures.fits
-astarithmetic $aperturesrebinned set-i \
-              i 0 ne 1 fill-holes set-good \
-              i good i 1 + where \
-              $quiet --output $apertures
-
-
-
-
-
 # Extract each measurement column(s)
 # ----------------------------------
 #
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index c40a5741..5776c6b9 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -25498,17 +25498,39 @@ The axis ratio of the apertures (minor axis divided 
by the major axis in a 2D el
 By default (when this option isn't given), the radial profile will be circular 
(axis ratio of 1).
 This parameter is used as the option @option{--qcol} in the generation of the 
apertures with @command{astmkprof}.
 
-@item -a FLT[,FLT[,...]]
-@itemx --azimuth=FLT[,FLT[,...]]
-Azimuthal angles range (in degrees, from 0 to 360) from the major axis.
-This parameter is used as the option @option{--qcol=azimuth} in the generation 
of the apertures with @command{astmkprof}.
-For example, if you want to obtain the radial profile only in between 55 and 
78 degrees (taking as the reference the major axis), use 
@command{--azimuth=55,78}.
-
 @item -p FLT
 @itemx --positionangle=FLT
 The position angle (in degrees) of the profiles relative to the first FITS 
axis (horizontal when viewed in SAO DS9).
 By default, it is @option{--positionangle=0}, which means that the semi-major 
axis of the profiles will be parallel to the first FITS axis.
 
+@item -a FLT,FLT
+@itemx --azimuth=FLT,FLT
+@cindex Wedge (radial profile)
+@cindex Azimuthal range (radial profile)
+Limit the profile to the given azimuthal angle range (two numbers given to 
this option, in degrees, from 0 to 360) from the major axis (defined by 
@option{--positionangle}).
+The radial profile will therefore be created on a wedge-like shape, not the 
full circle/ellipse.
+
+If the first angle is @emph{smaller} than the second (for example 
@option{--azimuth=10,80}), the region between, or @emph{inside}, the two angles 
will be used.
+Otherwise (for example @option{--azimuth=80,10}), the region @emph{outside} 
the two angles will be used.
+The latter case can be useful when you want to ignore part of the 2D shape 
(for example due to a bright star that can be contaminating it).
+
+You can visually see the shape of the region used by running this script with 
@option{--keeptmp} and viewing the @file{values.fits} and @file{apertures.fits} 
files of the temporary directory with a FITS image viewer like @ref{SAO DS9}.
+You can use @ref{Viewing FITS file contents with DS9 or TOPCAT} to open them 
together in one instance of DS9, with both frames matched and locked (for easy 
comparison in case you want to zoom-in or out).
+For example see the commands below (based on your target object, just change 
the image name, center, position angle and etc):
+
+@example
+## Generate the radial profile
+$ astscript-radial-profile image.fits --center=1.234,6.789 \
+            --mode=wcs --rmax=50 --positionangle=20 \
+            --axisratio=0.8 --azimuth=95,150 --keeptmp \
+            --tmpdir=radial-tmp
+
+## Visually check the values and apertures used.
+$ astscript-fits-view radial-tmp/values.fits \
+                      radial-tmp/apertures.fits
+@end example
+
+
 @item -m STR
 @itemx --measure=STR
 The operator for measuring the values over each radial distance.