gnuastro-commits
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[gnuastro-commits] master ac1c619c 2/2: radial-profile: added WCS to pol


From: Mohammad Akhlaghi
Subject: [gnuastro-commits] master ac1c619c 2/2: radial-profile: added WCS to polar plot and bug fix in MakeProfile
Date: Fri, 31 May 2024 19:09:31 -0400 (EDT)

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

    radial-profile: added WCS to polar plot and bug fix in MakeProfile
    
    Until now the polar plot did not have any WCS to better guide the users of
    the output on what every pixel corresponds to. Also, while doing some
    tests, I noticed a bug in MakeProfiles (it does not write negative values
    when '--replace' is used).
    
    With this commit, the first draft of the polar plot code from the previous
    commit is edited and polished prior to merging with the 'master' branch and
    the bug in MakeProfiles has been fixed.
    
    This fixes bug #65822.
---
 NEWS                         |   2 +
 bin/mkprof/mkprof.c          |  56 +++++++----
 bin/mkprof/oneprofile.c      |   2 +-
 bin/script/radial-profile.sh | 217 +++++++++++++++++++++++++++----------------
 doc/gnuastro.texi            |  20 ++--
 5 files changed, 191 insertions(+), 106 deletions(-)

diff --git a/NEWS b/NEWS
index 2e3cc409..9ca75347 100644
--- a/NEWS
+++ b/NEWS
@@ -192,6 +192,8 @@ See the end of the file for license conditions.
   - bug #65743: Arithmetic's collapse-number mistakenly calculating
     collapse-median instead. Reported by Sepideh Eskandarlou.
 
+  - bug #65822: MakeProfiles does not write negative values with --replace.
+
 
 
 
diff --git a/bin/mkprof/mkprof.c b/bin/mkprof/mkprof.c
index 6dfcfd18..2000cb91 100644
--- a/bin/mkprof/mkprof.c
+++ b/bin/mkprof/mkprof.c
@@ -191,28 +191,28 @@ saveindividual(struct mkonthread *mkp)
                         0, NULL, 0);
   if( mkp->func==PROFILE_SERSIC || mkp->func==PROFILE_MOFFAT )
     gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PINDEX", 0,
-                          &p->r[id], 0, "Index (Sersic or Moffat) of profile"
-                          "in catalog", 0, NULL, 0);
+                          &p->r[id], 0, "Index (Sersic or Moffat) of "
+                          "profile in catalog", 0, NULL, 0);
   if(ndim==2)
     {
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA_DEG", 0,
                             &p->p1[id], 0, "Position angle of profile in "
                             "catalog", 0, "deg", 0);
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO", 0,
-                            &p->q1[id], 0, "Axis ratio of profile in catalog",
-                            0, NULL, 0);
+                            &p->q1[id], 0, "Axis ratio of profile in "
+                            "catalog", 0, NULL, 0);
     }
   else
     {
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA1_DEG", 0,
-                            &p->p1[id], 0, "First X-Z-X Euler angle in 3D", 0,
-                            "deg", 0);
+                            &p->p1[id], 0, "First X-Z-X Euler angle in "
+                            "3D", 0, "deg", 0);
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA2_DEG", 0,
-                            &p->p2[id], 0, "Second X-Z-X Euler angle in 3D", 0,
-                            "deg", 0);
+                            &p->p2[id], 0, "Second X-Z-X Euler angle in "
+                            "3D", 0, "deg", 0);
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "PA3_DEG", 0,
-                            &p->p3[id], 0, "Third X-Z-X Euler angle in 3D", 0,
-                            "deg", 0);
+                            &p->p3[id], 0, "Third X-Z-X Euler angle in "
+                            "3D", 0, "deg", 0);
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "AXISRATIO1", 0,
                             &p->q1[id], 0, "Axis ratio along second dim",
                             0, NULL, 0);
@@ -230,8 +230,8 @@ saveindividual(struct mkonthread *mkp)
                         (void *)(p->rng_name), 0,
                         "Name of random number generator", 0, NULL, 0);
   gal_fits_key_list_add(&keys, GAL_TYPE_ULONG, "RNGSEED", 0,
-                        &mkp->rng_seed, 0, "Seed of random number generator",
-                        0, NULL, 0);
+                        &mkp->rng_seed, 0, "Seed of random number "
+                        "generator", 0, NULL, 0);
   gal_fits_key_list_add(&keys, GAL_TYPE_SIZE_T, "NUMRANDOM", 0,
                         &p->numrandom, 0,
                         "Number of random points in central pixels", 0,
@@ -247,8 +247,8 @@ saveindividual(struct mkonthread *mkp)
                         &p->oversample, 0, "Oversampling factor", 0,
                         NULL, 0);
   gal_fits_key_list_add(&keys, GAL_TYPE_UINT8, "TUNITINP", 0,
-                        &p->tunitinp, 0, "Truncation is in units of pixels, "
-                        "not radius", 0, NULL, 0);
+                        &p->tunitinp, 0, "Truncation is in units of "
+                        "pixels, not radius", 0, NULL, 0);
   if( !isnan(p->zeropoint) )
       gal_fits_key_list_add(&keys, GAL_TYPE_FLOAT32, "ZEROPOINT", 0,
                             &p->zeropoint, 0, "Zeropoint magnitude", 0,
@@ -367,7 +367,6 @@ mkprof_build_single(struct mkonthread *mkp, long *fpixel_i, 
long *lpixel_i,
           if(dsize[i] != ibq->image->dsize[i]) needs_crop=1;
         }
 
-
       /* Define the individual overlap tile. */
       if(needs_crop)
         {
@@ -382,7 +381,6 @@ mkprof_build_single(struct mkonthread *mkp, long *fpixel_i, 
long *lpixel_i,
                                     NULL, 0, -1, 1, NULL, NULL, NULL);
       ibq->overlap_i->block=ibq->image;
 
-
       /* Define the merged overlap tile. */
       ind=gal_dimension_coord_to_index(ndim, p->out->dsize, start_mrg);
       ptr=gal_pointer_increment(p->out->array, ind, p->out->type);
@@ -618,10 +616,21 @@ mkprof_write(struct mkprofparams *p)
   double sum;
   char *jobname;
   struct timeval t1;
+  float *o, *of, oinit=NAN;
   gal_data_t *out=p->out, *log;
   struct builtqueue *ibq=NULL, *tbq;
   size_t complete=0, num=p->num, clog;
 
+  /* If we want to replace the pixel values, initialize the output array
+     with an impossible value (smallest possible 32-bit floating point
+     number). This is because we are always replacing the pixels with the
+     maximum to ensure reproducibility in multi-threaded operations. */
+  if(p->replace)
+    {
+      gal_type_min(GAL_TYPE_FLOAT32, &oinit);
+      of=(o=out->array)+out->size; do *o++=oinit; while(o<of);
+    }
+
   /* Write each image into the output array. */
   while(complete<p->num)
     {
@@ -649,7 +658,11 @@ mkprof_write(struct mkprofparams *p)
          array. */
       if(ibq->overlaps && out)
         GAL_TILE_PO_OISET(float,float,ibq->overlap_i,ibq->overlap_m,1,0, {
-            *o  = p->replace ? ( *i>*o ? *i : *o ) :  (*i + *o);
+            /* The '*i>*o' condition is there to have a reproducible output
+               in multi-threaded scenarios (where the order can differ). */
+            *o  = ( p->replace
+                    ? ( *o==oinit ? *i : ( *i>*o ? *i : *o ) )
+                    :  (*i + *o) );
             sum += *i;
           });
 
@@ -705,6 +718,15 @@ mkprof_write(struct mkprofparams *p)
     }
 
 
+  /* In case we had initialized the output pixels, set all the ones that
+     were not filled back to zero. */
+  if(p->replace)
+    {
+      of=(o=out->array)+out->size;
+      do *o = *o==oinit ? 0.0 : *o; while(++o<of);
+    }
+
+
   /* Write the final array to the output FITS image if a merged image is to
      be created. */
   if(out)
diff --git a/bin/mkprof/oneprofile.c b/bin/mkprof/oneprofile.c
index 47d97902..22bc9588 100644
--- a/bin/mkprof/oneprofile.c
+++ b/bin/mkprof/oneprofile.c
@@ -851,7 +851,7 @@ oneprofile_make(struct mkonthread *mkp)
     {
       /* First get the sum of all the pixels in the profile. */
       ff=(f=mkp->ibq->image->array) + mkp->ibq->image->size;
-      sum=0.0f; do sum+= isnan(*f) ? 0.0f : *f; while(++f<ff);
+      sum=0.0f; do sum += isnan(*f) ? 0.0f : *f; while(++f<ff);
 
       /* Correct the fraction of brightness that was calculated
          accurately (not using the pixel center). */
diff --git a/bin/script/radial-profile.sh b/bin/script/radial-profile.sh
index b3b433c4..6b6e641b 100755
--- a/bin/script/radial-profile.sh
+++ b/bin/script/radial-profile.sh
@@ -114,7 +114,7 @@ $scriptname options:
   -Q, --axis-ratio=FLT     Axis ratio for ellipse profiles (A/B).
   -a, --azimuth=FLT,FLT    Azimuthal angles range (from the major axis).
   -p, --position-angle=FLT Position angle for ellipse profiles.
-  -g, --polar              Polar grid (Radius vs Azimuth angle).
+  -g, --polar              Polar plot (azimuth angle vs. radius).
   -s, --sigmaclip=FLT,FLT  Sigma-clip multiple and tolerance.
   -z, --zeropoint=FLT      Zeropoint magnitude of input dataset.
   -Z, --zeroisnotblank     0.0 in float or double images are not blank.
@@ -124,7 +124,7 @@ $scriptname options:
   -t, --tmpdir             Directory to keep temporary files.
   -k, --keeptmp            Keep temporal/auxiliar files.
   -m, --measure=STR        Measurement operator (mean, sigclip-mean, etc.).
-  -P, --precision=INT      Number of digits after the decimal point for radius.
+  -P, --precision=INT      Number of digits after decimal point for radius.
   -v, --oversample=INT     Oversample for higher resolution radial profile.
   -u, --undersample=INT    Undersample for lower resolution radial profile.
   -i, --instd=FLT/STR      Sky standard deviation per pixel, as a single
@@ -436,20 +436,19 @@ fi
 # If no specific measurement has been requested, use the mean.
 if [ x"$measure" = x ]; then measure=mean; fi
 
-# For polar-plot azimuth option is mandatory. Because without azimuthal, we
-# can not create a polar-plot.
-if [ x"$azimuth" = x ] && [ x$polar != x0 ]; then
+# Without azimuthal angle, we can not create a polar-plot.
+if [ x"$azimuth" = x ] && [ x$polar = x1 ]; then
     azimuth="0,360"
 fi
 
 # If the user call the '--undersample' and '--oversample' option. On the
 # other hand, they use '--polar' option with these options
-if [ x$polar != x0 ] && [ x"$undersample" != x ] ; then
-    echo "If you want to use the --undersample, please remove the --polar 
option."; exit 1
+if [ x$polar = x1 ] && [ x"$undersample" != x ] ; then
+    echo "$scriptname: the polar plot does not yet support --undersample; 
please get in touch with us at 'bug-gnuastro@gnu.org' for its possible 
addition"; exit 1
 fi
 
-if [ x$polar != x0 ] && [ x"$oversample" != x ] ; then
-    echo "If you want to use the --oversample, please remove the --polar 
option."; exit 1
+if [ x$polar = x1 ] && [ x"$oversample" != x ] ; then
+    echo "$scriptname: the polar plot does not yet support --oversample; 
please get in touch with us at 'bug-gnuastro@gnu.org' for its possible 
addition"; exit 1
 fi
 
 
@@ -497,7 +496,7 @@ fi
 
 
 
-# Calculate the maximum radiusf
+# Calculate the maximum radius
 # ----------------------------
 #
 # If the user didn't set the '--rmax' parameter, then compute the maximum
@@ -697,7 +696,7 @@ echo "1 $xcenter $ycenter 7 $rmax 1 $positionangle 
$axisratio 1 1" \
                  --circumwidth=1 --replace --output=$radialaperturesholed
 astarithmetic $radialaperturesholed set-i \
               i 0 ne 1 fill-holes set-good \
-              i good i 1 + where $precfactor x uint32\
+              i good i 1 + where $precfactor x uint32 \
               $quiet --output $radialapertures
 rm $radialaperturesholed
 
@@ -773,10 +772,10 @@ if [ x"$azimuth" != x ]; then
                      | awk ' {if ($1<$2) print "and"; \
                               else       print "or"}')
     astarithmetic --output=$aperturesraw $quiet \
-                  $radialapertures  -h1 set-radial \
-                  $azimuthapertures -h1 set-azimuth \
-                  azimuth $azimuth_ini ge \
-                  azimuth $azimuth_fin le $condition \
+                  $radialapertures  -h1 uint32 set-radial \
+                  $azimuthapertures -h1 uint16 set-azimuth \
+                  azimuth $azimuth_ini uint16 ge \
+                  azimuth $azimuth_fin uint16 le $condition \
                   radial 1 eq or set-arc \
                   radial arc 1 uint8 ne 0 uint8 where
 else
@@ -800,8 +799,18 @@ else                           
finalzp="--zeropoint=$zeropoint";
 fi
 if [ x$polar = x0 ]; then polaropt=""
 else
+    # We are using 'raw_azimuth_fin' because it can be set to zero when the
+    # value was 360. Note also that we should convert to an integer at the
+    # start (before 'set-i'), otherwise, in floating point, the 'gt'
+    # condition is going to miss the outer edge.
+    condition=$(echo $azimuth_ini $azimuth_fin \
+                     | awk ' {if ($1<$2) print "or"; \
+                              else       print "and"}')
     intazimuthapertures=$tmpdir/azimuth-raw-int.fits
-    astarithmetic $azimuthapertures uint32 --output=$intazimuthapertures
+    astarithmetic $azimuthapertures uint16 set-i i \
+                  i $azimuth_ini     uint16 lt \
+                  i $raw_azimuth_fin uint16 gt $condition \
+                  0 where uint32 --output=$intazimuthapertures $quiet
     polaropt="--clumpsfile=$intazimuthapertures --clumpshdu=1 --clumpscat"
 fi
 
@@ -809,8 +818,6 @@ fi
 
 
 
-
-
 # Undersampling the aperture image
 # --------------------------------
 #
@@ -879,9 +886,9 @@ finalmeasure=$(echo "$measure" \
 #   the center, but after adding this, the '--sb-error' column was nicely
 #   consistent with the scatter in the mock profiles at central radii.
 cat=$tmpdir/catalog.fits
-astmkcatalog $apertures --hdu=1 --valuesfile=$values --valueshdu=1 --ids \
-             $finalinstd $finalmeasure $finalsigmaclip $finalzp $polaropt \
-            --spatialresolution=0.0 --output=$cat $quiet
+astmkcatalog $apertures --hdu=1 --valuesfile=$values --valueshdu=1 \
+             --ids $finalinstd $finalmeasure $finalsigmaclip $finalzp \
+             $polaropt --spatialresolution=0.0 --output=$cat $quiet
 
 
 
@@ -897,8 +904,9 @@ astmkcatalog $apertures --hdu=1 --valuesfile=$values 
--valueshdu=1 --ids \
 radraw=$tmpdir/radii.fits
 if [ x"$oversample" != x ]; then
 
-    asttable $cat -c'arith OBJ_ID float32 '$precfactor' / 1 - '$oversample' /' 
\
-             -o$radraw --colmetadata=1,RADIUS,pix,"Radial distance"
+    asttable $cat -o$radraw \
+             --colmetadata=1,RADIUS,pix,"Radial distance" \
+             -c'arith OBJ_ID float32 '$precfactor' / 1 - '$oversample' /'
 
 elif [ x"$undersample" != x ]; then
 
@@ -945,86 +953,131 @@ asttable $radraw --catcolumnfile=$cat $restcols 
--output=$outraw \
 # 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=RADIUS,0,$rmax --output=$output
+#
+# We are not using '$output' here, because the user may have asked to make
+# a polar plot and we may confront some bugs in the middle. It is better
+# that the final output only be built after everything is complete.
+outputraw=$tmpdir/output.fits
+asttable $outraw --range=RADIUS,0,$rmax --output=$outputraw
+
+
 
+
+
+# Polar plot
+# ----------
+#
 # If the polar option is called by the user then the second extension of
 # the catalog that is created based on the azimuth (as clump) and radial
 # (as objects) will be added to the final PSF.
 if [ x$polar != x0 ]; then
 
-    # Limit the radius and change the metadata.
+    # Set metadata on the polar-plot catalog for easy operations below (and
+    # for debugging). To include the maximum radius (99 in case of
+    # '--rmax=100'), because Table's '--range' option doesn't include the
+    # maximum number and we are dealing with integers, we need to add
+    # 'rmax' with one.
     polarcat=$tmpdir/polar-catalog.fits
-    asttable $cat -hCLUMPS --range=HOST_OBJ_ID,0,$rmax \
+    rmaxpone=$(astarithmetic $rmax uint32 1 + --quiet)
+    asttable $cat -hCLUMPS --range=HOST_OBJ_ID,0,$rmaxpone \
              --colmetadata=HOST_OBJ_ID,RADIUS,pix,"Radial distance" \
-             --colmetadata=ID_IN_HOST_OBJ,AZIMUTH,degree,"Azimuth angel" \
+             --colmetadata=ID_IN_HOST_OBJ,AZIMUTH,degree,"Azimuth angle" \
              --output=$polarcat
 
-    # If polar option call then check the output is FITS file or not. The
-    # output should be a FITS file because just in fits files we can write
-    # the outputs of the polar option in the second extenction.
-    if astfits $output &> /dev/null; then
-
-        # Find the maximum value of radius and azimuth angel.
-        maxrad=$(aststatistics $polarcat -h1 -cRADIUS --maximum)
-        maxazim=$(aststatistics $polarcat -h1 -cAZIMUTH --maximum)
-
-        # Check if we have zero value in the measurement column replace
-        # them with maximum value of measurment plus one.
-        numzero=$(asttable $polarcat --equal=3,0 | wc -l)
-        if ! [ $numzero = 0 ]; then
-            mv $polarcat $polarcat.fits
-            maxmesure=$(aststatistics $polarcat.fits -h1 -c3 --maximum -q)
-            placeholder=$(astarithmetic $maxmesure 1 + -q)
-            asttable $polarcat.fits -cRADIUS,AZIMUTH \
-                     -c'arith $3 set-m m m 0 eq '$placeholder' where' \
-                     --output=$polarcat
-        fi
-
-
-        # Add the polar-plot to the third extension of the output.
-        mkprof=$tmpdir/polar-plot-profile.fits
-        asttable $polarcat -h1 -cRADIUS,AZIMUTH -c3  \
-            | awk '{print(1, $2, $1, 4, 0, 0, 0, 0, $3, 0)}' \
-            | astmkprof --mforflatpix --mode=img --oversample=1 \
-                        --mergedsize=$maxazim,$maxrad --clearcanvas --replace 
--quiet \
-                        --output=$mkprof
-
-        # In polar-plot image replavce the zero pixels with NAN value and
-        # replace pixels that have maximum value plus one with zero value.
-        if  [ $numzero = 0 ]; then
-            mv $mkprof $mkprof.fits
-            astarithmetic $mkprof.fits set-p \
-                          p p 0 eq nan where \
-                          --output=$mkprof
-        else
-            mv $mkprof $mkprof.fits
-            astarithmetic $mkprof.fits                  set-p \
-                          p p 0            eq nan where set-z \
-                          z z $placeholder eq 0   where \
-                          --output=$mkprof
-        fi
-
-        # Add second extension to the final outpu that is the polar-plot
-        # table and in third extension is the polar-plot image.
-        astfits $polarcat --copy=1 --output=$output
-        astfits $mkprof --copy=1 --output=$output
+    # Find the maximum value of radius and azimuth angel.
+    maxrad=$(aststatistics  $polarcat -h1 -cRADIUS  --maximum)
+    maxazim=$(aststatistics $polarcat -h1 -cAZIMUTH --maximum)
+    minazim=$(aststatistics $polarcat -h1 -cAZIMUTH --minimum)
+
+    # In the end, many azimuth angles close to radius zero are not going to
+    # be filled at all! We need to set those to NaN within the radial
+    # profile script (because there was no data there and we want to
+    # preserve the zero-valued pixels in our radial profiles). So here
+    # (before calling MakeProfiles), we replace any possible zeros in the
+    # input table with a pre-defined number (the maximum value in the image
+    # plus one), and later, rever this specific number back to zero.
+    polarcatz=$tmpdir/polar-catalog-no-zeros.fits
+    numzero=$(asttable $polarcat --equal=3,0 | wc -l)
+    if [ $numzero = 0 ]; then
+        ln -sf $polarcat $polarcatz
+    else
+        maxmesure=$(aststatistics $polarcat -h1 -c3 --maximum -q)
+        placeholder=$(astarithmetic $maxmesure 1 + -q)
+        asttable $polarcat -cRADIUS,AZIMUTH --output=$polarcatz \
+                 -c'arith $3 set-m m m 0 eq '$placeholder' where'
+    fi
 
+    # Generate the raw polar plot. When the user sets a non-zero minimum
+    # azimuthal angle, we need the full range of azimuthal angles here, so
+    # the azimuthal width is effectively just the maximum azimuthal
+    # angle. If the minimum is non-zero, we will later crop the output to
+    # remove the all-zero columns of the output image.
+    ppraw=$tmpdir/polar-plot-raw.fits
+    asttable $polarcat -h1 -cAZIMUTH,RADIUS,3  \
+        | awk '{print(1, $1, $2, 4, 0, 0, 0, 0, $3, 0)}' \
+        | astmkprof --mforflatpix --mode=img --oversample=1 \
+                    --mergedsize=$maxazim,$maxrad --clearcanvas \
+                    --output=$ppraw $quiet \
+                    --crpix=1,1 \
+                    --pc=1,0,0,1 \
+                    --cdelt=1,1 \
+                    --cunit=deg,pix \
+                    --ctype=AZI,RAD \
+                    --crval=1,0
+
+    # If the minimum azimuthal angle is not zero, then we need to crop the
+    # zero-valued regions (as described above).
+    ppcrop=$tmpdir/polar-plot-crop.fits
+    if [ $minazim = 0 ]; then
+        ln -sf $ppraw $ppcrop
     else
+        astcrop $ppraw --mode=img --section=$minazim:$maxazim,: \
+                --output=$ppcrop $quiet
+    fi
 
-        # IF the output is not a FITS file save it in another files with
-        # name of the output that has polar in it name and save exatcly in
-        # the output directory.
-        suffix=$(echo $output | awk 'BEGIN{FS="."}{print $NF}')
-        polar_catalog=$(echo $output | sed 's|.'$suffix'|-polar.'$suffix'|')
-        asttable $polarcat --output=$polar_catalog
+    # Any zero-valued pixel in the output of MakeProfiles means that there
+    # was no input data there so we should set them to Nan! Recall that
+    # above, if there was any actual zeros, they were replaced and we set
+    # them back to zero after setting the NaNs.
+    ppout=$tmpdir/polar-plot.fits
+    if  [ $numzero = 0 ]; then
+        astarithmetic $ppcrop set-p p p 0 eq nan where \
+                      --output=$ppout $quiet
+    else
+        astarithmetic $ppcrop                       set-p \
+                      p p 0            eq nan where set-z \
+                      z z $placeholder eq 0   where \
+                      --output=$ppout $quiet
+    fi
 
+    # If the output is a FITS file, then add a new HDU. Otherwise (for
+    # example output is plain-text), we need to make a new file)
+    if astfits $outputraw &> /dev/null; then
+        exthdu=2
+        polarfile=$outputraw
+    else
+        exthdu=1
+        suffix=$(echo $outputraw | awk 'BEGIN{FS="."}{print $NF}')
+        polarfile=$(echo $outputraw \
+                        | sed 's|.'$suffix'|-polar.'$suffix'|')
     fi
+
+    # Copy the image into the desired output file and set its HDU name.
+    astfits $ppout --copy=1 --output=$polarfile
+    astfits $outputraw --hdu=$exthdu --update=EXTNAME,POLAR-PLOT
 fi
 
 
 
 
 
+# Write the final output.
+cp $outputraw $output
+
+
+
+
+
 
 # Remove temporary files
 # ----------------------
diff --git a/doc/gnuastro.texi b/doc/gnuastro.texi
index 3553b809..193dfe5f 100644
--- a/doc/gnuastro.texi
+++ b/doc/gnuastro.texi
@@ -33721,6 +33721,10 @@ $ astscript-radial-profile image.fits
 ## up to  a radial distance of 19 pixels, use the mean value.
 $ astscript-radial-profile image.fits --center=44,37 --rmax=19
 
+## Generate a 2D polar plot with the same properties as above:
+$ astscript-radial-profile image.fits --center=44,37 --rmax=19 \
+                           --polar
+
 ## Generate the radial profile centered at x=44 and y=37 (in pixels),
 ## up to a radial distance of 100 pixels, compute sigma clipped
 ## mean and standard deviation (sigclip-mean and sigclip-std) using
@@ -33943,12 +33947,16 @@ $ astscript-fits-view radial-tmp/values.fits \
 
 @item -g
 @itemx --polar
-To generate 2D polar-plot, as a table.
-Which is contained azimuthal angle, radius and mean vlues of each pixel.
-Using @option{--azimuth} option is mandatory with it.
-This option can not to be used with @option{--oversample} and 
@option{--undersample} options.
-If the output is not a FITS file, the polar output will be saved in another 
table with the same name of output with polar.
-Otherwise, if the output is a FITS, in second extension of the output of 
polar-plot is existed as a table and in third extension the image.
+Generate a 2D polar-plot image in the second HDU of the output (called 
@code{POLAR-PLOT}).
+A polar plot is a projection of the original pixel grid into polar coordinates 
(where the horizontal axis is the azimuthal angle and the vertical axis is the 
radius).
+
+By default it assumes the full azimuthal range (from 0 to 360 degrees); if a 
narrower azimuthal range is desired, use @option{--azimuth} (for example 
@option{--azimuth=30,50} to only generate the polar plot between 30 and 50 
degrees of azimuth).
+
+The output image contains WCS information to map the pixel coordinates into 
the polar coordinates.
+This is especially useful when the azimuthal range is not the full range: the 
first pixel in the horizonal axis is not 0 degrees.
+
+Currently, the polar plot cannot to be used with @option{--oversample} and 
@option{--undersample} options (please get in touch with us if you need it).
+Until it is implemented, you can use the @option{--scale} option of @ref{Warp} 
to do the oversampling of the input image yourself and generate the polar plot 
from that.
 
 
 @item -m STR



reply via email to

[Prev in Thread] Current Thread [Next in Thread]