[gnuastro-commits] master 145070a2 21/23: color-faint-gray: new option --segment to define the color regions

[Mohammad Akhlaghi]

branch: master
commit 145070a202038bfddae13a143dc0cc81984334da
Author: Raul Infante-Sainz <infantesainz@gmail.com>
Commit: Mohammad Akhlaghi <mohammad@akhlaghi.org>

    color-faint-gray: new option --segment to define the color regions
    
    Until this commit, the way of separation the color, black, and gray regions
    were internally computed. These regions could be changed with two
    parameters, but it was always done by the script internally. However, in
    some situations it is useful to define these regions in a quantitative way.
    
    With this commit, a new option, '--segment' has been added. With this
    option it is possible to provide a segmentation image as an input that will
    be used for defining the different color regions.  Pixels equal to 2 in
    this image will be shown in color, pixels equal to 1 in pure black, and
    pixels equal to 0 will be shown in gray. Other minor modifications have
    been also done.
---
 bin/script/color-faint-gray.sh | 126 ++++++++++++++++++++++-------------------
 doc/gnuastro.texi              |  99 +++++++++++++++++++++++++++++++-
 2 files changed, 166 insertions(+), 59 deletions(-)

diff --git a/bin/script/color-faint-gray.sh b/bin/script/color-faint-gray.sh
index b2bc5ea1..124a3f22 100644
--- a/bin/script/color-faint-gray.sh
+++ b/bin/script/color-faint-gray.sh
@@ -44,6 +44,7 @@ export LANG=C
 
 # Default option values (can be changed with options on the command-line).
 hdu=""
+shdu=1
 globalhdu=""
 
 # Minimum, weights, and zeropoint values
@@ -51,13 +52,15 @@ weight=""
 minimum=""
 zeropoint=""
 
-# To control the asinh transformation
-qbright=1.0
-stretch=1.0
+# To control the asinh transformation, set both to 1 (scientific notation)
+qbright_default=$(astarithmetic 1.0 1.0 x --quiet)
+stretch_default=$(astarithmetic 1.0 1.0 x --quiet)
 
-# For color and gray background
+# For color, black, and gray regions
+segment=""
 grayval=""
 colorval=""
+coloronly=0
 graykernelfwhm=0
 colorkernelfwhm=0
 
@@ -69,7 +72,7 @@ colorkernelfwhm=0
 minvalrange=0.000
 maxvalrange=100.0
 
-# Transformation parameters to improve the contrast and bias
+# To enhance the image
 bias=0.0
 gamma=1.0
 contrast=1.0
@@ -78,9 +81,7 @@ quiet=""
 tmpdir=""
 keeptmp=0
 checkparams=0
-output="color-faint-gray.jpg"
-
-coloronly=0
+output="color-faint-gray.pdf"
 
 version=@VERSION@
 scriptname=@SCRIPT_NAME@
@@ -124,6 +125,7 @@ experienced Gnuastro users and developers. For more 
information, please run:
 $scriptname options:
  Input:
   -h, --hdu=STR           HDU/extension for the input channels.
+  -t, --shdu=STR          HDU/extension for the segment image.
   -g, --globalhdu=STR/INT Use this HDU for all inputs, ignore '--hdu'.
   -w, --weight=FLT        Relative weight for each input channel.
   -m, --minimum=FLT       Minimum value for each input channel.
@@ -139,7 +141,8 @@ $scriptname options:
   -G, --gamma             Gamma parameter (nonlinear, overrides bias/contrast).
 
  Color and gray parameters
-      --coloronly         No grayscale regions, background in black.
+      --coloronly         No grayscale regions, background in color (black).
+      --segment=STR       Segmentation image (color=2, black=1, gray=0).
       --grayval=FLT       Gray threshold (highest value to use grayscale).
       --colorval=FLT      Color threshold (lowest value to have color).
       --graykernelfwhm=FLT  Kernel FWHM for convolving the background image.
@@ -242,6 +245,9 @@ do
         -h|--hdu)            aux="$2";                                  
check_v "$1" "$aux"; hdu="$hdu $aux"; shift;shift;;
         -h=*|--hdu=*)        aux="${1#*=}";                             
check_v "$1" "$aux"; hdu="$hdu $aux"; shift;;
         -h*)                 aux="$(echo "$1"  | sed -e's/-h//')";      
check_v "$1" "$aux"; hdu="$hdu $aux"; shift;;
+        -t|--shdu)           shdu="$2";                                 
check_v "$1" "$shdu";  shift;shift;;
+        -t=*|--shdu=*)       shdu="${1#*=}";                            
check_v "$1" "$shdu";  shift;;
+        -t*)                 shdu=$(echo "$1" | sed -e's/-t//');        
check_v "$1" "$shdu";  shift;;
         -w|--weight)         aux="$2";                                  
check_v "$1" "$aux"; weight="$weight $aux"; shift;shift;;
         -w=*|--weight=*)     aux="${1#*=}";                             
check_v "$1" "$aux"; weight="$weight $aux"; shift;;
         -w*)                 aux="$(echo "$1"  | sed -e's/-w//')";      
check_v "$1" "$aux"; weight="$weight $aux"; shift;;
@@ -270,6 +276,8 @@ do
         -b*)                 bias=$(echo "$1"  | sed -e's/-b//');       
check_v "$1" "$bias";  shift;;
 
         --coloronly)        coloronly=1; shift;;
+        --segment)          segment="$2";                              check_v 
"$1" "$segment";  shift;shift;;
+        --segment=*)        segment="${1#*=}";                         check_v 
"$1" "$segment";  shift;;
         --grayval)          grayval="$2";                              check_v 
"$1" "$grayval";  shift;shift;;
         --grayval=*)        grayval="${1#*=}";                         check_v 
"$1" "$grayval";  shift;;
         --colorval)         colorval="$2";                             check_v 
"$1" "$colorval";  shift;shift;;
@@ -431,6 +439,19 @@ else
 fi
 
 
+# Bright. If the user provides --qbright, make sure it is not equal to zero
+# (with 8 decimals), that would crash the asinh transformation.
+if [ x$qbright != x ]; then
+    qbright_check=$(echo "$qbright" | awk 'BEGIN{FS=","} {printf "%.8f", $1}')
+    if [ x$qbright_check = x0.00000000 ]; then
+        echo "$scriptname: --qbright value ($qbright) cannot be zero (8 
decimals)."
+        exit 1
+    fi
+else
+    qbright=$qbright_default
+fi
+
+
 # Stretch. If the user provides --stretch, make sure it is not equal to
 # zero (with 8 decimals), that would crash the asinh transformation.
 if [ x$stretch != x ]; then
@@ -439,19 +460,11 @@ if [ x$stretch != x ]; then
         echo "$scriptname: --stretch value ($stretch) cannot be zero (8 
decimals)."
         exit 1
     fi
+else
+    stretch=$stretch_default
 fi
 
 
-# Bright. If the user provides --qbright, make sure it is not equal to zero
-# (with 8 decimals), that would crash the asinh transformation.
-if [ x$qbright != x ]; then
-    qbright_check=$(echo "$qbright" | awk 'BEGIN{FS=","} {printf "%.8f", $1}')
-    if [ x$qbright_check = x0.00000000 ]; then
-        echo "$scriptname: --qbright value ($qbright) cannot be zero (8 
decimals)."
-        exit 1
-    fi
-fi
-
 
 
 
@@ -773,9 +786,6 @@ else
     # channels.
 
 
-
-
-
     # Background image
     # ----------------
     #
@@ -832,42 +842,43 @@ else
     # asinh-transformed image. If the user does not provide a value then use
     # ghe computed one (guessed). If the user provide a value, then use it
     # directly. Note that the guessed value is computed in any case.
-    colorval_guessed=$(aststatistics $I_COLORGRAY_threshold --median --quiet)
-    grayval_guessed=$colorval_guessed
+    colorval_estimated=$(aststatistics $I_COLORGRAY_threshold --median --quiet)
     if [ x$colorval = x"" ]; then
-      colorval=$colorval_guessed
+      colorval=$colorval_estimated
     fi
 
+    grayval_estimated=$(aststatistics $I_COLORGRAY_threshold --median --quiet)
     if [ x$grayval = x"" ]; then
-      grayval=$colorval
+      grayval=$grayval_estimated
     fi
 
 
 
 
 
-    # Mask images
-    # -----------
+    # Segmentation image
+    # ------------------
     #
-    # Different mask are computed here: color, black, and gray images. They
-    # correspond to the regions that are represented in those colormaps.
-    # At the end, for debugging and finding the best paramters pourpouses,
-    # a total mask is computed.
-    GRAY_MASK="$tmpdir/GRAY_mask.fits"
-    COLOR_MASK="$tmpdir/COLOR_mask.fits"
-    BLACK_MASK="$tmpdir/BLACK_mask.fits"
-    TOTAL_MASK="$tmpdir/TOTAL_mask-2color-1black-0gray.fits"
-    astarithmetic $I_COLORGRAY_threshold set-i \
-                  i $colorval gt --output $COLOR_MASK
-    astarithmetic $I_COLORGRAY_threshold set-i \
-                  i $colorval lt i $grayval gt and \
-                  --output $BLACK_MASK
-    astarithmetic $I_COLORGRAY_threshold set-i \
-                  i $colorval lt --output $GRAY_MASK
-    astarithmetic $COLOR_MASK -h1 f32 2 x set-c \
-                  $BLACK_MASK -h1 f32 1 x set-b \
-                  $GRAY_MASK  -h1 f32 0 x set-g \
-                  c b g 3 sum uint8 --output $TOTAL_MASK
+    # The TOTAL_MASK consists in a segmentation image whose pixel values
+    # correspond to the three regions. Here it is defined as follow:
+    #   pixels=2 will be shown in color
+    #   pixels=1 will be shown in pure black
+    #   pixels=0 will be shown in gray
+    # This image is computed by default from the colorval and grayval
+    # parameters. Alternatively, it can be provided by the user.
+    if [ x$segment = x ]; then
+        TOTAL_MASK="$tmpdir/TOTAL_mask-2color-1black-0gray.fits"
+        astarithmetic $I_COLORGRAY_threshold                     set-i \
+                      i $colorval gt                   2 uint8 x set-c \
+                      i $colorval lt i $grayval gt and 1 uint8 x set-b \
+                      i $colorval lt                   0 uint8 x set-g \
+                      c b g 3 sum uint8 --output $TOTAL_MASK
+        shdu=1
+   else
+       grayval="$segment"
+       colorval="$segment"
+       TOTAL_MASK=$segment
+   fi
 
 
 
@@ -878,7 +889,8 @@ else
     #
     # We use the specified gray background image after some modifications to
     # obtain the GRAY part. For futher highlight the interesting diffuse flux,
-    # we need to invert it the brigher pixels are darker.
+    # we need to invert it the brigher pixels are darker. Here, the color
+    # region is masked (pixels=2 in TOTAL_MASK)
     #
     # Here is the logic of the steps we take:
     #   1. We call the convolved image 'values'.
@@ -896,8 +908,8 @@ else
     #   function is specified. E.g., log, sqrt, asinh, etc.)
     grayscale=""
     I_GRAY_colormasked="$tmpdir/GRAY_colormasked.fits"
-    astarithmetic $I_BACK_convolved -h1 set-values \
-                  $COLOR_MASK       -h1 set-mask \
+    astarithmetic $I_BACK_convolved -h1      set-values \
+                  $TOTAL_MASK -h$shdu 2 uint8 eq set-mask \
                   values mask nan where $grayscale set-masked \
                   masked minvalue set-oldmin \
                   masked maxvalue set-oldmax \
@@ -913,14 +925,14 @@ else
 
 
 
-    # Set the black region to zero pixel values
-    # -----------------------------------------
+    # Set the pure black region to zero pixel values
+    # ----------------------------------------------
     #
-    # Put the black pixels equal to zero. By doing this, those pixels will
-    # be set to pure black color in the final image.
+    # Put black pixels (pixels=1 in TOTAL_MASK) equal to zero. By doing
+    # this, those pixels will be set to pure black color.
     I_GRAY_colormasked_zeroblack="$tmpdir/GRAY_colormasked_zeroblack.fits"
     astarithmetic $I_GRAY_colormasked -h1 set-i \
-                  $BLACK_MASK         -h1 set-b \
+                  $TOTAL_MASK -h$shdu 1 uint8 eq -h1 set-b \
                   i b 0.0 where float32 \
                   --output=$I_GRAY_colormasked_zeroblack
 
@@ -1012,8 +1024,8 @@ TIPS:
       '$I_COLORGRAY_threshold'
 
 PARAMETERS:
-  Estimated:  --colorval=$colorval_guessed --grayval=$grayval_guessed
-  Used     :  --colorval=$colorval --grayval=$grayval --qbright=$qbright 
--stretch=$stretch
+  Default:  --qbright=$qbright_default --stretch=$stretch_default 
--colorval=$colorval_estimated --grayval=$grayval_estimated
+  Used   :  --qbright=$qbright --stretch=$stretch --colorval=$colorval 
--grayval=$grayval
 
 Output written to '$output'.
 EOF
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index fa641b5f..d41be600 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -9165,11 +9165,106 @@ $ astscript-fits-view 
aux/TOTAL-mask-2color-1black-0gray.fits
 In this segmentation image, pixels equal to 2 will be shown in color, pixels 
equal to 1 will be shown as pure black, and pixels equal to zero are shown in 
gray.
 By default, the script sets the same value for both thresholds.
 That means that there is not pure black pixels unless they are equal to zero 
in the three channels.
-
 By adjusting these two parameters, you can obtain an optimal result to show 
the bright and faint parts of your data within one printable image.
 The values used here are somewhat extreme to illustrate the logic of the 
procedure, but we encourage you to experiment with values close to the 
estimated by default in order to have a smooth transition between the three 
regions (color, black, and gray).
 The script can provide additional information about the pixel value 
distributions used to estimate the parameters by using the 
@option{--checkparams} option.
-We will use @option{--colorval=5} and @option{--grayval=1} in what follows.
+
+A very interesting characteristic is the possibility of providing the 
segmentation image as an input.
+This feature expand the possibilities of generating color images in a more 
quantitative way.
+For example, let's use this for showing the different surface brightness 
intervals in the color image.
+First, let's generate a surface brightness image from the R channel.
+That is, the value of each pixel will be in the units of surface brightness 
(mag/arcsec@mymath{^2}).
+To to that, we need obtain the pixel area in arcsec and use the zeropoint 
value of the image.
+Then, the @option{counts-to-sb} operator of @command{astarithmetic} is used.
+
+@example
+$ Rzpoint=23.43
+$ Rchannel=aligned/i-jplus.fits
+$ pixarea=$(astfits $Rchannel --pixelareaarcsec2 --quiet)
+
+# Compute the SB image
+$ astarithmetic $Rchannel $Rzpoint $pixarea counts-to-sb \
+                --output=R-sb.fits
+
+# Have a look at the image
+$ astscript-fits-view R-sb.fits
+@end example
+
+Remember that because @file{R-sb.fits} is a surface brightness image, lower 
pixel values mean brighter regions, and higher pixel values mean fainter 
regions.
+Note also that there are some NaN pixels, they correspond to negative values 
on the original image, we will take care of them later.
+
+Now let's define the segmentation image (@file{R-sb-seg.fits}) step-by-step 
with the following criteria:
+
+- Bright parts, SB < 21, in color: @file{R-sb-seg.fits}@mymath{=2}
+- Faint parts, 21 < SB < 25, in black: @file{R-sb-seg.fits}@mymath{=1}
+- Background parts, 25 < SB, in gray: @file{R-sb-seg.fits}@mymath{=0}
+
+@example
+# SB thresholds (low and high)
+$ sbh=21
+$ sbl=25
+
+# Bright parts in color
+$ astarithmetic R-sb.fits set-i \
+                i $sbh lt 2 uint8 x set-c \
+                c --output=color.fits
+
+# Faint parts in black
+$ astarithmetic R-sb.fits set-i \
+                i $sbh gt \
+                i $sbl lt \
+                and 1 uint8 x set-b \
+                b --output=black.fits
+
+# Background in gray
+$ astarithmetic R-sb.fits set-i \
+                i $sbl gt 0 uint8 x set-g \
+                g --output=gray.fits
+
+# Check the images
+$ astscript-fits-view color.fits black.fits gray.fits
+@end example
+
+Have a look at the three images and note check that they define the three 
regions according to the surface brightness levels from the @file{R-sb.fits} 
image.
+This is great!
+Now we can sum these images to obtain the segmentation image and use it as an 
input so the color image will consider it to generate the different colormap 
parts.
+But before this, we need to account for the NaN pixels.
+Essentially, we need to consider those NaN pixels to be shown in gray because 
they are noise or background (negative values on the original image that become 
NaN in the SB image).
+Let's summarize the above commands, the sum of the images, and the NaN 
conversion, in a single command.
+
+@example
+$ astarithmetic sb.fits set-i \
+                i isblank set-n \
+                i $sbh lt 2 uint8 x set-c \
+                i $sbl gt 0 uint8 x set-g \
+                i $sbh gt i $sbl lt and 1 uint8 x set-b \
+                c b g 3 sum uint8 set-seg \
+                seg n 0 where set-segnzero \
+                segnzero --output=R-sb-seg.fits
+
+# Check the image
+$ astscript-fits-view R-sb-seg.fits
+@end example
+
+Let's say in words what this mean.
+Pixels on this image equal to 2 will be shown in color, pixels equal to 1 will 
be shown in pure black, pixels equal to 0 will be shown in gray.
+Now, let's use this segmentation image with the @option{--segment} option for 
obtaining the color image.
+
+@example
+$ astscript-color-faint-gray $params --output=m51-sb.pdf \
+                             --qbright=0.01 --stretch=100 \
+                             --segment=R-sb-seg.fits
+@end example
+
+Open @file{R-sb.pdf} and have a look.
+Do you see how the different regions (SB intervals) have been coloured 
differently?
+They come from the SB levels we defined, and because of that, this is actually 
a quantitative image that tells you brightness information and not only a nice 
color image!
+The separation between the color and black regions is not smooth because we 
wanted to show these diferent regions, but playing and modifying the parameters 
it is possible to obtain better results.
+This is really interesting because now it is possible to use color images for 
detecting low surface brightness features at the same time they provide 
quantitative measurements.
+Of course, here we have defined this segmentation image just using two SB 
values, but it is possible to define any other segmentation image that you may 
need for your particular science.
+
+In what follow we will not use the segmentation image.
+To make things simple, let's use @option{--colorval=5} and 
@option{--grayval=1} on the internally computed image for the thresholds.
 
 @example
 params="$R $G $B \