diff options
Diffstat (limited to 'pkg/proto/ringavg.cl')
| -rw-r--r-- | pkg/proto/ringavg.cl | 172 |
1 files changed, 172 insertions, 0 deletions
diff --git a/pkg/proto/ringavg.cl b/pkg/proto/ringavg.cl new file mode 100644 index 00000000..a7504642 --- /dev/null +++ b/pkg/proto/ringavg.cl @@ -0,0 +1,172 @@ +# RINGAVG (Nov02) proto RINGAVG (Nov02) +# +# +# NAME +# ringavg -- compute pixel averages in concentric rings about given +# center +# +# +# USAGE +# ringavg image xc yc +# +# +# PARAMETERS +# +# image +# Image to be used. +# +# xc, yc +# Pixel coordinate for center of rings. +# +# r1 = 0, r2 = 10, dr = 1 +# Rings to be measured. R1 is the inner radius of the first ring, +# R2 is the outer radius of the last bin, and DR is the widths of +# the rings. The values are in units of pixels. +# +# labels = yes +# Print column labels for the output? +# +# vebar = no +# If VEBAR is yes then the standard deviation and standard error +# will be printed as negative values for use with GRAPH. +# +# +# DESCRIPTION +# Pixels are binned into a series of concentric rings centered on a +# given position in the input image. The rings are defined by an +# inner radius for the first ring, an outer radius for the last ring, +# and the width of the rings. The statistics of the pixel values in +# each ring are then computed and list to the standard output. The +# output lines consist of the inner and outer ring radii, the number +# of pixels, the average value, the standard deviation of the value +# (corrected for population size), and the standard error. The +# parameter LABEL selects whether to include column labels. +# +# If the ring average are to be plotted with the task GRAPH using the +# option to plot error bars based on the standard deviation or +# standard error then the VEBAR parameter may be set to write the +# values as negative values are required by that task. +# +# This task is a script and so users may copy it and modify it as +# desired. Because it is a script it will be very slow if r2 becomes +# large. +# +# +# EXAMPLES +# 1. Compute the ring averages with labels and output to the terminal. +# +# cl> ringavg pwimage 17 17 +# # R min R max Npix Average Std Dev Std Err +# 0.00 1.00 5 7.336 9.16 4.096 +# 1.00 2.00 8 0.2416 0.2219 0.07844 +# 2.00 3.00 16 0.3994 0.5327 0.1332 +# 3.00 4.00 20 0.06211 0.05491 0.01228 +# 4.00 5.00 32 0.0987 0.08469 0.01497 +# 5.00 6.00 32 0.06983 0.06125 0.01083 +# 6.00 7.00 36 0.0641 0.0839 0.01398 +# 7.00 8.00 48 0.06731 0.05373 0.007755 +# 8.00 9.00 56 0.06146 0.07601 0.01016 +# 9.00 10.00 64 0.05626 0.05846 0.007308 +# +# 2. Plot the ring averages with standard errors used for error bars. +# +# cl> ringavg pwimage 17 17 label- vebar+ | fields STDIN 2,4,6 | +# >>> graph point+ marker=vebar +# +# 3. Plot ring averages for galaxy in dev$pix. +# +# cl> ringavg dev$pix 256 256 r2=100 dr=5 label- | fields STDIN 2,4 | +# >>> graph logy+ +# +# +# +# SEE ALSO +# pradprof, psfmeasure, radprof +# +# +# To install: +# +# Copy to your home or other personal directory. Enter the command +# "task ringavg = home$ringavg.cl" interactively, in login.cl or in +# your loginuser.cl. Substitute the host or logical path for home$ +# if the script is placed in a different directory. + + +procedure ringavg (image, xc, yc) + +file image {prompt="Input image"} +real xc {prompt="X center"} +real yc {prompt="Y center"} + +real r1 = 0 {prompt="Inner radius of first bin"} +real r2 = 10 {prompt="Outer radius of last bin"} +real dr = 1 {prompt="Radial bin width"} + +bool labels = yes {prompt="Print column labels?"} +bool vebars = no {prompt="Format for error bars in GRAPH?"} + +struct *fd + +begin + file temp + real n, r, val, ra, rb, ravg, rstddev, rmean + + # Extract the pixel values sorted by radius. + temp = mktemp ("temp") + pradprof (image, xc, yc, radius=r2, center=no, list=yes) | + sort ("STDIN", column=1, ignore_white+, numeric+, reverse-, > temp) + + if (label) + printf ("# %6s %8s %8s %10s %10s %10s\n", "R min", "R max", "Npix", + "Average", "Std Dev", "Std Err") + + # Read through the file. Skip the first two comment lines. + fd = temp + i = fscan (fd) + i = fscan (fd) + n = 0 + rb = -1 + while (fscan (fd, r, val) != EOF) { + if (r < r1) + next + if (r > r2) + break + if (r > rb) { + if (n > 0) { + ravg = ravg / n + rstddev = sqrt (rstddev / n - ravg ** 2) + if (vebar) + rstddev = -rstddev + if (n > 1) + rstddev = rstddev * sqrt (n / (n - 1.)) + rmean = rstddev / sqrt (n) + printf ("%8.2f %8.2f %8d %10.4g %10.4g %10.4g\n", + ra, rb, n, ravg, rstddev, rmean) + } + ravg = 0. + rstddev = 0. + n = 0 + ra = int (r / dr) * dr + rb = ra + dr + } + + ravg = ravg + val + rstddev = rstddev + val * val + n = n + 1 + } + + if (n > 0) { + ravg = ravg / n + rstddev = sqrt (rstddev / n - ravg ** 2) + if (vebar) + rstddev = -rstddev + if (n > 1) + rstddev = rstddev * sqrt (n / (n - 1.)) + rmean = rstddev / sqrt (n) + printf ("%8.2f %8.2f %8d %10.4g %10.4g %10.4g\n", + ra, rb, n, ravg, rstddev, rmean) + } + + fd = "" + delete (temp, verify-) +end |