aboutsummaryrefslogtreecommitdiff
path: root/pkg/proto/t_bscale.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/proto/t_bscale.x')
-rw-r--r--pkg/proto/t_bscale.x581
1 files changed, 581 insertions, 0 deletions
diff --git a/pkg/proto/t_bscale.x b/pkg/proto/t_bscale.x
new file mode 100644
index 00000000..c3133454
--- /dev/null
+++ b/pkg/proto/t_bscale.x
@@ -0,0 +1,581 @@
+include <imhdr.h>
+include <error.h>
+include <ctype.h>
+include <mach.h>
+
+define OPTIONS "|mean|median|mode|"
+define MEAN 1 # mean of image
+define MEDIAN 2 # median of image
+define MODE 3 # mode of image
+
+define BINWIDTH 0.1 # bin width (in sigmas) for mode
+define BINSEP 0.01 # bin separation (in sigmas) for mode
+
+# T_BSCALE -- Linearly transform the intensity scales of a list of images
+# using the following expression.
+#
+# out = (in - bzero) / bscale
+
+procedure t_bscale ()
+
+pointer inlist # list of input images
+pointer outlist # list of output images
+real bzero # zero point
+real bscale # scale factor
+real lower # lower limit for mean, median, or mode computation
+real upper # upper limit for mean, median, or mode computation
+pointer section # image section for statistics
+int step # default sample step
+int verbose # verbose mode
+
+double temp
+int i, bz, bs
+real mean, median, mode, sigma, tlower, tupper
+pointer sp, str, image1, image2, imtemp, inim, outim
+
+bool clgetb()
+int btoi(), imtopenp(), strdic(), gctod(), clgeti(), imtgetim(), imtlen()
+pointer immap()
+real clgetr()
+
+begin
+ # Allocate working space.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (image1, SZ_FNAME, TY_CHAR)
+ call salloc (image2, SZ_FNAME, TY_CHAR)
+ call salloc (imtemp, SZ_FNAME, TY_CHAR)
+ call salloc (section, SZ_FNAME, TY_CHAR)
+
+ # Open the input and output image lists.
+ inlist = imtopenp ("input")
+ outlist = imtopenp ("output")
+ if (imtlen(inlist) != imtlen(outlist)) {
+ call sfree (sp)
+ call imtclose (inlist)
+ call imtclose (outlist)
+ call error (0, "Length of input and output lists not equal.")
+ }
+
+ # Get the zero point algorithm.
+ call clgstr ("bzero", Memc[str], SZ_LINE)
+ bz = strdic (Memc[str], Memc[str], SZ_LINE, OPTIONS)
+ if (bz == 0) {
+ i = 1
+ if (gctod (Memc[str], i, temp) == 0)
+ call error (0, "Invalid `bzero' parameter")
+ bzero = temp
+ }
+
+ # Get the scale algorithm.
+ call clgstr ("bscale", Memc[str], SZ_LINE)
+ bs = strdic (Memc[str], Memc[str], SZ_LINE, OPTIONS)
+ if (bs == 0) {
+ i = 1
+ if (gctod (Memc[str], i, temp) == 0)
+ call error (0, "Invalid `bscale' parameter")
+ bscale = temp
+ }
+
+ # Get the section to be used for statistics computation.
+ call clgstr ("section", Memc[section], SZ_FNAME)
+ step = max (1, clgeti ("step"))
+
+ # Get the upper and lower good data limits.
+ lower = clgetr ("lower")
+ if (IS_INDEFR(lower))
+ tlower = -MAX_REAL
+ else
+ tlower = lower
+ upper = clgetr ("upper")
+ if (IS_INDEFR(upper))
+ tupper = MAX_REAL
+ else
+ tupper = upper
+
+ verbose = btoi (clgetb ("verbose"))
+
+ # Loop over the input and output image lists.
+ while ((imtgetim (inlist, Memc[image1], SZ_FNAME) != EOF) &&
+ (imtgetim (outlist, Memc[image2], SZ_FNAME) != EOF)) {
+
+ # Open the input and output images.
+ call xt_mkimtemp (Memc[image1], Memc[image2], Memc[imtemp],
+ SZ_FNAME)
+ iferr (inim = immap (Memc[image1], READ_ONLY, 0)) {
+ call erract (EA_WARN)
+ next
+ }
+ iferr (outim = immap (Memc[image2], NEW_COPY, inim)) {
+ call imunmap (inim)
+ call erract (EA_WARN)
+ next
+ }
+
+ # Compute the required statistics.
+ if ((bz != 0) || (bs != 0))
+ call bs_imstats (inim, Memc[section], step, BINWIDTH, BINSEP,
+ mean, median, mode, sigma, tupper, tlower)
+ else {
+ mean = INDEF
+ median = INDEF
+ mode = INDEF
+ }
+
+ switch (bz) {
+ case MODE:
+ bzero = mode
+ case MEAN:
+ bzero = mean
+ case MEDIAN:
+ bzero = median
+ }
+
+ switch (bs) {
+ case MODE:
+ bscale = mode
+ case MEAN:
+ bscale = mean
+ case MEDIAN:
+ bscale = median
+ }
+
+ # Log the output.
+ if (verbose == YES) {
+ call bs_log (Memc[image1], Memc[imtemp], mean, median, mode,
+ bzero, bscale, upper, lower)
+ call flush (STDOUT)
+ }
+
+ # Scale the image.
+ call bs_scale (inim, outim, bzero, bscale)
+
+ call imunmap (inim)
+ call imunmap (outim)
+ call xt_delimtemp (Memc[image2], Memc[imtemp])
+ }
+
+ call imtclose (inlist)
+ call imtclose (outlist)
+ call sfree (sp)
+end
+
+
+# BSCALE -- Scale the image brightness.
+
+procedure bs_scale (inim, outim, bzero, bscale)
+
+pointer inim # pointer to the input image
+pointer outim # pointer to the output image
+real bzero # zero point
+real bscale # scale
+
+int nc
+long v1[IM_MAXDIM], v2[IM_MAXDIM]
+real bz, bs
+
+pointer in, out
+int imgnlr(), impnlr(), imgnlx(), impnlx(), imgnld(), impnld()
+
+begin
+ call amovkl (long(1), v1, IM_MAXDIM)
+ call amovkl (long(1), v2, IM_MAXDIM)
+
+ bz = -bzero
+ bs = 1. / bscale
+ nc = IM_LEN(inim,1)
+
+ switch (IM_PIXTYPE(inim)) {
+ case TY_DOUBLE :
+ while ((imgnld (inim, in, v1) != EOF) && (impnld (outim,
+ out, v2) != EOF))
+ call altad (Memd[in], Memd[out], nc, double(bz), double(bs))
+
+ case TY_COMPLEX:
+ while ((imgnlx (inim, in, v1) != EOF) && (impnlx (outim,
+ out, v2) != EOF))
+ call altax (Memx[in], Memx[out], nc, bz, bs)
+
+ default:
+ while ((imgnlr (inim, in, v1) != EOF) && (impnlr (outim,
+ out, v2) != EOF))
+ call altar (Memr[in], Memr[out], nc, bz, bs)
+ }
+end
+
+
+# BS_LOG -- Log the scaling operation.
+
+procedure bs_log (image1, image2, mean, median, mode, bzero, bscale, upper,
+ lower)
+
+char image1[ARB] # input image name
+char image2[ARB] # output image name
+real mean # input image mean
+real median # input image median
+real mode # input image mode
+real bzero, bscale # the computed scale values
+real upper # upper limit for mean
+real lower # lower limit for mean
+
+begin
+ call printf ("%s -> %s using bzero: %g and bscale: %g\n")
+ call pargstr (image1)
+ call pargstr (image2)
+ call pargr (bzero)
+ call pargr (bscale)
+
+ if (! IS_INDEF(mean)) {
+ call printf (" mean: %g median: %g mode: %g ")
+ call pargr (mean)
+ call pargr (median)
+ call pargr (mode)
+ call printf (" upper: %g lower: %g\n")
+ call pargr (upper)
+ call pargr (lower)
+ }
+end
+
+
+# BS_IMSTATS -- Compute the image statistics within a section of an image.
+# This routine parses the image section and samples the image. The actual
+# statistics are evaluated by BS_STATS.
+
+procedure bs_imstats (im, section, step, binwidth, binsep, mean, median, mode,
+ sigma, upper, lower)
+
+pointer im # input image
+char section[ARB] # sample section
+int step # default sample step
+real binwidth # bin width
+real binsep # separation between bins
+real mean # mean
+real median # median
+real mode # mode
+real sigma # sigma
+real upper # upper limit for statistics
+real lower # lower limit for statistics
+
+int i, n, nx, ndim
+pointer sp, x1, x2, xs, v, v1, dv, data, ptr1, ptr2
+int imgnlr()
+
+begin
+ call smark (sp)
+ call salloc (x1, IM_MAXDIM, TY_INT)
+ call salloc (x2, IM_MAXDIM, TY_INT)
+ call salloc (xs, IM_MAXDIM, TY_INT)
+ call salloc (v, IM_MAXDIM, TY_LONG)
+ call salloc (v1, IM_MAXDIM, TY_LONG)
+ call salloc (dv, IM_MAXDIM, TY_LONG)
+
+ # Initialize the section.
+ ndim = IM_NDIM(im)
+ do i = 1, ndim {
+ Memi[x1+i-1] = 1
+ Memi[x2+i-1] = IM_LEN(im,i)
+ Memi[xs+i-1] = 0
+ }
+
+ # Parse the sample section.
+ call bs_section (section, Memi[x1], Memi[x2], Memi[xs], ndim)
+
+ # Check the step sizes.
+ do i = 1, ndim {
+ if (Memi[xs+i-1] == 0)
+ Memi[xs+i-1] = step
+ }
+
+ # Define the region of the image to be extracted.
+ n = 1
+ do i = 1, ndim {
+ nx = (Memi[x2+i-1] - Memi[x1+i-1]) / Memi[xs+i-1] + 1
+ Meml[v+i-1] = Memi[x1+i-1]
+ if (nx == 1)
+ Meml[dv+i-1] = 1
+ else
+ Meml[dv+i-1] = (Memi[x2+i-1] - Memi[x1+i-1]) / (nx - 1)
+ n = n * nx
+ }
+
+ # Accumulate the pixel values within the section.
+ call salloc (data, n, TY_REAL)
+ Meml[v] = 1
+ ptr1 = data
+ call amovl (Meml[v], Meml[v1], IM_MAXDIM)
+ while (imgnlr (im, ptr2, Meml[v1]) != EOF) {
+
+ ptr2 = ptr2 + Memi[x1] - 1
+ do i = Memi[x1], Memi[x2], Meml[dv] {
+ Memr[ptr1] = Memr[ptr2]
+ ptr1 = ptr1 + 1
+ ptr2 = ptr2 + Meml[dv]
+ }
+
+ for (i=2; i<=ndim; i=i+1) {
+ Meml[v+i-1] = Meml[v+i-1] + Meml[dv+i-1]
+ if (Meml[v+i-1] <= Memi[x2+i-1])
+ break
+ Meml[v+i-1] = Memi[x1+i-1]
+ }
+
+ if (i > ndim)
+ break
+
+ call amovl (Meml[v], Meml[v1], IM_MAXDIM)
+ }
+
+ # Compute the statistics.
+ call bs_stats (Memr[data], n, binwidth, binsep, mean, median, mode,
+ sigma, upper, lower)
+
+ call sfree (sp)
+end
+
+
+# BS_STATS -- Compute the vector statistics.
+#
+# 1. Sort the data
+# 2. Exclude the extreme points
+# 3. The median is at the midpoint of the sorted data
+# 4. Compute the mean
+# 5. Compute the sigmas about the mean
+# 6. Scale the bin width and separations by the sigma
+# 7. Find the mode over all the bins (which may overlap)
+
+procedure bs_stats (data, npts, binwidth, binsep, mean, median, mode, sigma,
+ upper, lower)
+
+real data[npts] # sata array which will be sorted.
+int npts # number of data points
+real binwidth # bin width
+real binsep # separation between bins
+real mean # mean
+real median # median
+real mode # mode
+real sigma # sigma
+real upper # upper limit for mean
+real lower # lower limit for mean
+
+int x1, x2, x3, n, nmax
+real width, sep, y1, y2
+int bs_awvgr()
+
+begin
+ # Initialize.
+ if (npts <= 0) {
+ mean = INDEFR
+ median = INDEFR
+ mode = INDEFR
+ sigma = INDEFR
+ return
+ }
+
+ # Sort the data.
+ call asrtr (data, data, npts)
+
+ # Find the array indices for the lower and upper data bounds.
+ x1 = 1
+ while (data[x1] < lower)
+ x1 = x1 + 1
+ x3 = npts
+ while (data[x3] > upper)
+ x3 = x3 - 1
+
+ # Assign number of elements within the bounds.
+ n = x3 - x1 + 1
+
+ # Compute the median.
+ median = data[x1 + n/2 - 1]
+ mode = median
+
+ # Compute the mean and sigma.
+ if (bs_awvgr (data[x1], n, mean, sigma, 0.0, 0.0) <= 0)
+ return
+
+ # Check for no dispersion in the data.
+ if (sigma <= 0.0)
+ return
+
+ width = binwidth * sigma
+ sep = binsep * sigma
+
+ # Compute the mode.
+ nmax = 0
+ x2 = x1
+ for (y1 = data[x1]; x2 < x3; y1 = y1 + sep) {
+ for (; data[x1] < y1; x1 = x1 + 1)
+ ;
+ y2 = y1 + width
+ for (; (x2 < x3) && (data[x2] < y2); x2 = x2 + 1)
+ ;
+ if (x2 - x1 > nmax) {
+ nmax = x2 - x1
+ mode = data[(x2+x1)/2]
+ }
+ }
+end
+
+
+# BS_AWVGR -- Compute the mean and standard deviation (sigma) of a sample.
+# Pixels whose value lies outside the specified lower and upper limits are
+# not used. If the upper and lower limits have the same value (e.g., zero),
+# no limit checking is performed. The number of pixels in the sample is
+# returned as the function value.
+
+int procedure bs_awvgr (a, npix, mean, sigma, lcut, hcut)
+
+real a[ARB] # input array of data
+int npix # the number of data points
+real mean # the computed mean
+real sigma # the computed standard deviation
+real lcut, hcut # lower and upper cutoff for statistics calculation
+
+int i, ngpix
+real value, sum, sumsq, temp
+
+begin
+ sum = 0.0
+ sumsq = 0.0
+ ngpix = 0
+
+ # Accumulate sum, sum of squares. The test to disable limit checking
+ # requires numerical equality of two floating point numbers; this
+ # should be ok since they are used as flags not as numbers (they are
+ # not used in computations).
+
+ if (hcut == lcut) {
+ do i = 1, npix {
+ value = a[i]
+ sum = sum + value
+ sumsq = sumsq + value ** 2
+ }
+ ngpix = npix
+
+ } else {
+ do i = 1, npix {
+ value = a[i]
+ if (value >= lcut && value <= hcut) {
+ ngpix = ngpix + 1
+ sum = sum + value
+ sumsq = sumsq + value ** 2
+ }
+ }
+ }
+
+ switch (ngpix) { # compute mean and sigma
+ case 0:
+ mean = INDEFR
+ sigma = INDEFR
+ case 1:
+ mean = sum
+ sigma = INDEFR
+ default:
+ mean = sum / ngpix
+ temp = (sumsq - mean * sum) / (ngpix - 1)
+ if (temp <= 0.0)
+ sigma = 0.0
+ else
+ sigma = sqrt (temp)
+ }
+
+ return (ngpix)
+end
+
+
+# BS_SECTION -- Parse an image section into its elements.
+#
+# 1. The default values must be set by the caller.
+# 2. A null image section is OK.
+# 3. The first nonwhitespace character must be '['.
+# 4. The last interpreted character must be ']'.
+#
+# This procedure should be replaced with an IMIO procedure at some point.
+
+procedure bs_section (section, x1, x2, xs, ndim)
+
+char section[ARB] # Image section
+int x1[ndim] # Starting pixel
+int x2[ndim] # Ending pixel
+int xs[ndim] # Step size in pixels
+int ndim # Number of dimensions
+
+int i, ip, a, b, c, temp, ctoi()
+define error_ 99
+
+begin
+ # Decode the section string.
+ ip = 1
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == '[')
+ ip = ip + 1
+ else if (section[ip] == EOS)
+ return
+ else
+ goto error_
+
+ do i = 1, ndim {
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == ']')
+ break
+
+ # Default values
+ a = x1[i]
+ b = x2[i]
+ c = xs[i]
+
+ # Get a:b:c. Allow notation such as "-*:c"
+ # (or even "-:c") where the step is obviously negative.
+
+ if (ctoi (section, ip, temp) > 0) { # a
+ a = temp
+ if (section[ip] == ':') {
+ ip = ip + 1
+ if (ctoi (section, ip, b) == 0) # a:b
+ goto error_
+ } else
+ b = a
+ c = 1
+ } else if (section[ip] == '-') { # -*
+ temp = a
+ a = b
+ b = temp
+ c = 1
+ ip = ip + 1
+ if (section[ip] == '*')
+ ip = ip + 1
+ } else if (section[ip] == '*') { # *
+ c = 1
+ ip = ip + 1
+ }
+
+ if (section[ip] == ':') { # :step
+ ip = ip + 1
+ if (ctoi (section, ip, c) == 0)
+ goto error_
+ else if (c == 0)
+ goto error_
+ }
+
+ if ((a > b) && (c > 0))
+ c = -c
+
+ x1[i] = a
+ x2[i] = b
+ xs[i] = c
+
+ while (IS_WHITE(section[ip]))
+ ip = ip + 1
+ if (section[ip] == ',')
+ ip = ip + 1
+ }
+
+ if (section[ip] != ']')
+ goto error_
+
+ return
+error_
+ call error (0, "Error in image section specification")
+end
generated by cgit (git 2.34.1) at 2025-09-20 05:12:47 -0400