aboutsummaryrefslogtreecommitdiff
path: root/pkg/images/tv/display/sigm2.x
diff options
context:
space:
mode:
Diffstat (limited to 'pkg/images/tv/display/sigm2.x')
-rw-r--r--pkg/images/tv/display/sigm2.x1110
1 files changed, 1110 insertions, 0 deletions
diff --git a/pkg/images/tv/display/sigm2.x b/pkg/images/tv/display/sigm2.x
new file mode 100644
index 00000000..41a3b5da
--- /dev/null
+++ b/pkg/images/tv/display/sigm2.x
@@ -0,0 +1,1110 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <imhdr.h>
+include <error.h>
+
+.help sigm2, sigm2_setup
+.nf ___________________________________________________________________________
+SIGM2 -- Get a line from a spatially scaled 2-dimensional image. This procedure
+works like the regular IMIO get line procedure, but rescales the input
+2-dimensional image in either or both axes upon input. If the magnification
+ratio required is greater than 0 and less than 2 then linear interpolation is
+used to resample the image. If the magnification ratio is greater than or
+equal to 2 then the image is block averaged by the smallest factor which
+reduces the magnification to the range 0-2 and then interpolated back up to
+the desired size. In some cases this will smooth the data slightly, but the
+operation is efficient and avoids aliasing effects.
+
+ si = sigm2_setup (im,pm, x1,x2,nx,xblk, y1,y2,ny,yblk, order)
+ sigm2_free (si)
+ ptr = sigm2[sr] (si, linenumber)
+
+SIGM2_SETUP must be called to set up the transformations after mapping the
+image and before performing any scaled i/o to the image. SIGM2_FREE must be
+called when finished to return buffer space.
+
+The SIGM routines are like SIGL routines except for the addition of
+interpolation over bad pixels and order=-1 takes the maximum rather
+than the average when doing block averaging or interpolation.
+.endhelp ______________________________________________________________________
+
+# Scaled image descriptor for 2-dim images
+
+define SI_LEN 19
+define SI_MAXDIM 2 # images of 2 dimensions supported
+define SI_NBUFS 3 # nbuffers used by SIGL2
+
+define SI_IM Memi[$1] # pointer to input image header
+define SI_FP Memi[$1+1] # pointer to fixpix structure
+define SI_GRID Memi[$1+2+$2-1] # pointer to array of X coords
+define SI_NPIX Memi[$1+4+$2-1] # number of X coords
+define SI_BAVG Memi[$1+6+$2-1] # X block averaging factor
+define SI_INTERP Memi[$1+8+$2-1] # interpolate X axis
+define SI_BUF Memi[$1+10+$2-1]# line buffers
+define SI_BUFY Memi[$1+13+$2-1]# Y values of buffers
+define SI_ORDER Memi[$1+15] # interpolator order
+define SI_TYBUF Memi[$1+16] # buffer type
+define SI_XOFF Memi[$1+17] # offset in input image to first X
+define SI_INIT Memi[$1+18] # YES until first i/o is done
+
+define OUTBUF SI_BUF($1,3)
+
+define SI_TOL (1E-5) # close to a pixel
+define INTVAL (abs ($1 - nint($1)) < SI_TOL)
+define SWAPI {tempi=$2;$2=$1;$1=tempi}
+define SWAPP {tempp=$2;$2=$1;$1=tempp}
+define NOTSET (-9999)
+
+# SIGM2_SETUP -- Set up the spatial transformation for SIGL2[SR]. Compute
+# the block averaging factors (1 if no block averaging is required) and
+# the sampling grid points, i.e., pixel coordinates of the output pixels in
+# the input image.
+
+pointer procedure sigm2_setup (im, pm, px1,px2,nx,xblk, py1,py2,ny,yblk, order)
+
+pointer im # the input image
+pointer pm # pixel mask
+real px1, px2 # range in X to be sampled on an even grid
+int nx # number of output pixels in X
+int xblk # blocking factor in x
+real py1, py2 # range in Y to be sampled on an even grid
+int ny # number of output pixels in Y
+int yblk # blocking factor in y
+int order # interpolator order (0=replicate, 1=linear)
+
+int npix, noldpix, nbavpix, i, j
+int npts[SI_MAXDIM] # number of output points for axis
+int blksize[SI_MAXDIM] # block averaging factor (npix per block)
+real tau[SI_MAXDIM] # tau = p(i+1) - p(i) in fractional pixels
+real p1[SI_MAXDIM] # starting pixel coords in each axis
+real p2[SI_MAXDIM] # ending pixel coords in each axis
+real scalar, start
+pointer si, gp, xt_fpinit()
+
+begin
+ iferr (call calloc (si, SI_LEN, TY_STRUCT))
+ call erract (EA_FATAL)
+
+ SI_IM(si) = im
+ SI_FP(si) = xt_fpinit (pm, 1, INDEFI)
+ SI_NPIX(si,1) = nx
+ SI_NPIX(si,2) = ny
+ SI_ORDER(si) = order
+ SI_INIT(si) = YES
+
+ p1[1] = px1 # X = index 1
+ p2[1] = px2
+ npts[1] = nx
+ blksize[1] = xblk
+
+ p1[2] = py1 # Y = index 2
+ p2[2] = py2
+ npts[2] = ny
+ blksize[2] = yblk
+
+ # Compute block averaging factors if not defined.
+ # If there is only one pixel then the block average is the average
+ # between the first and last point.
+
+ do i = 1, SI_MAXDIM {
+ if ((blksize[i] >= 1) && !IS_INDEFI (blksize[i])) {
+ if (npts[i] == 1)
+ tau[i] = 0.
+ else
+ tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+ } else {
+ if (npts[i] == 1) {
+ tau[i] = 0.
+ blksize[i] = int (p2[i] - p1[i] + 1 + SI_TOL)
+ } else {
+ tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+ if (tau[i] >= 2.0) {
+
+ # If nx or ny is not an integral multiple of the block
+ # averaging factor, noldpix is the next larger number
+ # which is an integral multiple. When the image is
+ # block averaged pixels will be replicated as necessary
+ # to fill the last block out to this size.
+
+ blksize[i] = int (tau[i] + SI_TOL)
+ npix = p2[i] - p1[i] + 1
+ noldpix = (npix+blksize[i]-1) / blksize[i] * blksize[i]
+ nbavpix = noldpix / blksize[i]
+ scalar = real (nbavpix - 1) / real (noldpix - 1)
+ p1[i] = (p1[i] - 1.0) * scalar + 1.0
+ p2[i] = (p2[i] - 1.0) * scalar + 1.0
+ tau[i] = (p2[i] - p1[i]) / (npts[i] - 1)
+ } else
+ blksize[i] = 1
+ }
+ }
+ }
+
+ SI_BAVG(si,1) = blksize[1]
+ SI_BAVG(si,2) = blksize[2]
+
+# if (IS_INDEFI (xblk))
+# xblk = blksize[1]
+# if (IS_INDEFI (yblk))
+# yblk = blksize[2]
+
+ # Allocate and initialize the grid arrays, specifying the X and Y
+ # coordinates of each pixel in the output image, in units of pixels
+ # in the input (possibly block averaged) image.
+
+ do i = 1, SI_MAXDIM {
+ # The X coordinate is special. We do not want to read entire
+ # input image lines if only a range of input X values are needed.
+ # Since the X grid vector passed to ALUI (the interpolator) must
+ # contain explicit offsets into the vector being interpolated,
+ # we must generate interpolator grid points starting near 1.0.
+ # The X origin, used to read the block averaged input line, is
+ # given by XOFF.
+
+ if (i == 1) {
+ SI_XOFF(si) = int (p1[i] + SI_TOL)
+ start = p1[1] - int (p1[i] + SI_TOL) + 1.0
+ } else
+ start = p1[i]
+
+ # Do the axes need to be interpolated?
+ if (INTVAL(start) && INTVAL(tau[i]))
+ SI_INTERP(si,i) = NO
+ else
+ SI_INTERP(si,i) = YES
+
+ # Allocate grid buffer and set the grid points.
+ iferr (call malloc (gp, npts[i], TY_REAL))
+ call erract (EA_FATAL)
+ SI_GRID(si,i) = gp
+ if (SI_ORDER(si) <= 0) {
+ do j = 0, npts[i]-1
+ Memr[gp+j] = int (start + (j * tau[i]) + 0.5 + SI_TOL)
+ } else {
+ do j = 0, npts[i]-1
+ Memr[gp+j] = start + (j * tau[i])
+ }
+ }
+
+ return (si)
+end
+
+
+# SIGM2_FREE -- Free storage associated with an image opened for scaled
+# input. This does not close and unmap the image.
+
+procedure sigm2_free (si)
+
+pointer si
+int i
+
+begin
+ # Free fixpix structure.
+ call xt_fpfree (SI_FP(si))
+
+ # Free SIGM2 buffers.
+ do i = 1, SI_NBUFS
+ if (SI_BUF(si,i) != NULL)
+ call mfree (SI_BUF(si,i), SI_TYBUF(si))
+
+ # Free GRID buffers.
+ do i = 1, SI_MAXDIM
+ if (SI_GRID(si,i) != NULL)
+ call mfree (SI_GRID(si,i), TY_REAL)
+
+ call mfree (si, TY_STRUCT)
+end
+
+
+# SIGM2S -- Get a line of type short from a scaled image. Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure sigm2s (si, lineno)
+
+pointer si # pointer to SI descriptor
+int lineno
+
+pointer rawline, tempp, gp
+int i, new_y[2], tempi, curbuf, altbuf
+int nraw, npix, nblks_y, ybavg, x1, x2
+real x, y, weight_1, weight_2
+pointer si_blmavgs()
+errchk si_blmavgs
+
+begin
+ nraw = IM_LEN(SI_IM(si),1)
+ npix = SI_NPIX(si,1)
+
+ # Determine the range of X (in pixels on the block averaged input image)
+ # required for the interpolator.
+
+ gp = SI_GRID(si,1)
+ x1 = SI_XOFF(si)
+ x = Memr[gp+npix-1]
+ x2 = x1 + int(x)
+ if (INTVAL(x))
+ x2 = x2 - 1
+ x2 = max (x1 + 1, x2)
+
+ gp = SI_GRID(si,2)
+ y = Memr[gp+lineno-1]
+
+ # The following is an optimization provided for the case when it is
+ # not necessary to interpolate in either X or Y. Block averaging is
+ # permitted.
+
+ if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+ return (si_blmavgs (SI_IM(si), SI_FP(si), x1, x2, int(y),
+ SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si)))
+
+ # If we are interpolating in Y two buffers are required, one for each
+ # of the two input image lines required to interpolate in Y. The lines
+ # stored in these buffers are interpolated in X to the output grid but
+ # not in Y. Both buffers are not required if we are not interpolating
+ # in Y, but we use them anyhow to simplify the code.
+
+ if (SI_INIT(si) == YES) {
+ do i = 1, 2 {
+ if (SI_BUF(si,i) != NULL)
+ call mfree (SI_BUF(si,i), SI_TYBUF(si))
+ call malloc (SI_BUF(si,i), npix, TY_SHORT)
+ SI_TYBUF(si) = TY_SHORT
+ SI_BUFY(si,i) = NOTSET
+ }
+ if (OUTBUF(si) != NULL)
+ call mfree (OUTBUF(si), SI_TYBUF(si))
+ call malloc (OUTBUF(si), npix, TY_SHORT)
+ SI_INIT(si) = NO
+ }
+
+ # If the Y value of the new line is not in range of the contents of the
+ # current line buffers, refill one or both buffers. To refill we must
+ # read a (possibly block averaged) input line and interpolate it onto
+ # the X grid. The X and Y values herein are in the coordinate system
+ # of the (possibly block averaged) input image.
+
+ new_y[1] = int(y)
+ new_y[2] = int(y) + 1
+
+ # Get the pair of lines whose integral Y values form an interval
+ # containing the fractional Y value of the output line. Sometimes the
+ # desired line will happen to be in the other buffer already, in which
+ # case we just have to swap buffers. Often the new line will be the
+ # current line, in which case nothing is done. This latter case occurs
+ # frequently when the magnification ratio is large.
+
+ curbuf = 1
+ altbuf = 2
+
+ do i = 1, 2 {
+ if (new_y[i] == SI_BUFY(si,i)) {
+ ;
+ } else if (new_y[i] == SI_BUFY(si,altbuf)) {
+ SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+ SWAPI (SI_BUFY(si,1), SI_BUFY(si,2))
+
+ } else {
+ # Get line and interpolate onto output grid. If interpolation
+ # is not required merely copy data out. This code is set up
+ # to always use two buffers; in effect, there is one buffer of
+ # look ahead, even when Y[i] is integral. This means that we
+ # will go out of bounds by one line at the top of the image.
+ # This is handled by copying the last line.
+
+ ybavg = SI_BAVG(si,2)
+ nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+ if (new_y[i] <= nblks_y)
+ rawline = si_blmavgs (SI_IM(si), SI_FP(si), x1, x2,
+ new_y[i], SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si))
+
+ if (SI_INTERP(si,1) == NO) {
+ call amovs (Mems[rawline], Mems[SI_BUF(si,i)], npix)
+ } else if (SI_ORDER(si) == 0) {
+ call si_samples (Mems[rawline], Mems[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ } else if (SI_ORDER(si) == -1) {
+ call si_maxs (Mems[rawline], nraw,
+ Memr[SI_GRID(si,1)], Mems[SI_BUF(si,i)], npix)
+ } else {
+ call aluis (Mems[rawline], Mems[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ }
+
+ SI_BUFY(si,i) = new_y[i]
+ }
+
+ SWAPI (altbuf, curbuf)
+ }
+
+ # We now have two line buffers straddling the output Y value,
+ # interpolated to the X grid of the output line. To complete the
+ # bilinear interpolation operation we take a weighted sum of the two
+ # lines. If the range from SI_BUFY(si,1) to SI_BUFY(si,2) is repeatedly
+ # interpolated in Y no additional i/o occurs and the linear
+ # interpolation operation (ALUI) does not have to be repeated (only the
+ # weighted sum is required). If the distance of Y from one of the
+ # buffers is zero then we do not even have to take a weighted sum.
+ # This is not unusual because we may be called with a magnification
+ # of 1.0 in Y.
+
+ weight_1 = 1.0 - (y - SI_BUFY(si,1))
+ weight_2 = 1.0 - weight_1
+
+ if (weight_1 < SI_TOL)
+ return (SI_BUF(si,2))
+ else if (weight_2 < SI_TOL || SI_ORDER(si) == 0)
+ return (SI_BUF(si,1))
+ else if (SI_ORDER(si) == -1) {
+ call amaxs (Mems[SI_BUF(si,1)], Mems[SI_BUF(si,2)],
+ Mems[OUTBUF(si)], npix)
+ return (OUTBUF(si))
+ } else {
+ call awsus (Mems[SI_BUF(si,1)], Mems[SI_BUF(si,2)],
+ Mems[OUTBUF(si)], npix, weight_1, weight_2)
+ return (OUTBUF(si))
+ }
+end
+
+
+# SI_BLMAVGS -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on. If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blmavgs (im, fp, x1, x2, y, xbavg, ybavg, order)
+
+pointer im # input image
+pointer fp # fixpix structure
+int x1, x2 # range of x blocks to be read
+int y # y block to be read
+int xbavg, ybavg # X and Y block averaging factors
+int order # averaging option
+
+real sum
+short blkmax
+pointer sp, a, b
+int nblks_x, nblks_y, ncols, nlines, xoff, blk1, blk2, i, j, k
+int first_line, nlines_in_sum, npix, nfull_blks, count
+pointer xt_fps()
+errchk xt_fps
+
+begin
+ call smark (sp)
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ xoff = (x1 - 1) * xbavg + 1
+ npix = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1) - xoff + 1
+
+ if ((xbavg < 1) || (ybavg < 1))
+ call error (1, "si_blmavg: illegal block size")
+ else if (x1 < 1 || x2 > ncols)
+ call error (2, "si_blmavg: column index out of bounds")
+ else if ((xbavg == 1) && (ybavg == 1))
+ return (xt_fps (fp, im, y, NULL) + xoff - 1)
+
+ nblks_x = (npix + xbavg-1) / xbavg
+ nblks_y = (nlines + ybavg-1) / ybavg
+
+ if (y < 1 || y > nblks_y)
+ call error (2, "si_blmavg: block number out of range")
+
+ if (ybavg > 1) {
+ call salloc (b, nblks_x, TY_LONG)
+ call aclrl (Meml[b], nblks_x)
+ nlines_in_sum = 0
+ }
+
+ # Read and accumulate all input lines in the block.
+ first_line = (y - 1) * ybavg + 1
+
+ do i = first_line, min (nlines, first_line + ybavg - 1) {
+ # Get line from input image.
+ a = xt_fps (fp, im, i, NULL) + xoff - 1
+
+ # Block average line in X.
+ if (xbavg > 1) {
+ # First block average only the full blocks.
+ nfull_blks = npix / xbavg
+ if (order == -1) {
+ blk1 = a
+ do j = 1, nfull_blks {
+ blk2 = blk1 + xbavg
+ blkmax = Mems[blk1]
+ do k = blk1+1, blk2-1
+ blkmax = max (blkmax, Mems[k])
+ Mems[a+j-1] = blkmax
+ blk1 = blk2
+ }
+ } else
+ call abavs (Mems[a], Mems[a], nfull_blks, xbavg)
+
+ # Now average the final partial block, if any.
+ if (nfull_blks < nblks_x) {
+ if (order == -1) {
+ blkmax = Mems[blk1]
+ do k = blk1+1, a+npix-1
+ blkmax = max (blkmax, Mems[k])
+ Mems[a+j-1] = blkmax
+ } else {
+ sum = 0.0
+ count = 0
+ do j = nfull_blks * xbavg + 1, npix {
+ sum = sum + Mems[a+j-1]
+ count = count + 1
+ }
+ Mems[a+nblks_x-1] = sum / count
+ }
+ }
+ }
+
+ # Add line into block sum. Keep track of number of lines in sum
+ # so that we can compute block average later.
+
+ if (ybavg > 1) {
+ if (order == -1) {
+ do j = 0, nblks_x-1
+ Meml[b+j] = max (Meml[b+j], long (Mems[a+j]))
+ } else {
+ do j = 0, nblks_x-1
+ Meml[b+j] = Meml[b+j] + Mems[a+j]
+ nlines_in_sum = nlines_in_sum + 1
+ }
+ }
+ }
+
+ # Compute the block average in Y from the sum of all lines block
+ # averaged in X. Overwrite buffer A, the buffer returned by IMIO.
+ # This is kosher because the block averaged line is never longer
+ # than an input line.
+
+ if (ybavg > 1) {
+ if (order == -1) {
+ do i = 0, nblks_x-1
+ Mems[a+i] = Meml[b+i]
+ } else {
+ do i = 0, nblks_x-1
+ Mems[a+i] = Meml[b+i] / real(nlines_in_sum)
+ }
+ }
+
+ call sfree (sp)
+ return (a)
+end
+
+
+# SI_MAXS -- Resample a line via maximum value.
+
+procedure si_maxs (a, na, x, b, nb)
+
+short a[na] # input array
+int na # input size
+real x[nb] # sample grid
+short b[nb] # output arrays
+int nb # output size
+
+int i
+
+begin
+ do i = 1, nb
+ b[i] = max (a[int(x[i])], a[min(na,int(x[i]+1))])
+end
+
+
+# SIGM2I -- Get a line of type short from a scaled image. Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure sigm2i (si, lineno)
+
+pointer si # pointer to SI descriptor
+int lineno
+
+pointer rawline, tempp, gp
+int i, new_y[2], tempi, curbuf, altbuf
+int nraw, npix, nblks_y, ybavg, x1, x2
+real x, y, weight_1, weight_2
+pointer si_blmavgi()
+errchk si_blmavgi
+
+begin
+ nraw = IM_LEN(SI_IM(si),1)
+ npix = SI_NPIX(si,1)
+
+ # Determine the range of X (in pixels on the block averaged input image)
+ # required for the interpolator.
+
+ gp = SI_GRID(si,1)
+ x1 = SI_XOFF(si)
+ x = Memr[gp+npix-1]
+ x2 = x1 + int(x)
+ if (INTVAL(x))
+ x2 = x2 - 1
+ x2 = max (x1 + 1, x2)
+
+ gp = SI_GRID(si,2)
+ y = Memr[gp+lineno-1]
+
+ # The following is an optimization provided for the case when it is
+ # not necessary to interpolate in either X or Y. Block averaging is
+ # permitted.
+
+ if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+ return (si_blmavgi (SI_IM(si), SI_FP(si), x1, x2, int(y),
+ SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si)))
+
+ # If we are interpolating in Y two buffers are required, one for each
+ # of the two input image lines required to interpolate in Y. The lines
+ # stored in these buffers are interpolated in X to the output grid but
+ # not in Y. Both buffers are not required if we are not interpolating
+ # in Y, but we use them anyhow to simplify the code.
+
+ if (SI_INIT(si) == YES) {
+ do i = 1, 2 {
+ if (SI_BUF(si,i) != NULL)
+ call mfree (SI_BUF(si,i), SI_TYBUF(si))
+ call malloc (SI_BUF(si,i), npix, TY_INT)
+ SI_TYBUF(si) = TY_INT
+ SI_BUFY(si,i) = NOTSET
+ }
+ if (OUTBUF(si) != NULL)
+ call mfree (OUTBUF(si), SI_TYBUF(si))
+ call malloc (OUTBUF(si), npix, TY_INT)
+ SI_INIT(si) = NO
+ }
+
+ # If the Y value of the new line is not in range of the contents of the
+ # current line buffers, refill one or both buffers. To refill we must
+ # read a (possibly block averaged) input line and interpolate it onto
+ # the X grid. The X and Y values herein are in the coordinate system
+ # of the (possibly block averaged) input image.
+
+ new_y[1] = int(y)
+ new_y[2] = int(y) + 1
+
+ # Get the pair of lines whose integral Y values form an interval
+ # containing the fractional Y value of the output line. Sometimes the
+ # desired line will happen to be in the other buffer already, in which
+ # case we just have to swap buffers. Often the new line will be the
+ # current line, in which case nothing is done. This latter case occurs
+ # frequently when the magnification ratio is large.
+
+ curbuf = 1
+ altbuf = 2
+
+ do i = 1, 2 {
+ if (new_y[i] == SI_BUFY(si,i)) {
+ ;
+ } else if (new_y[i] == SI_BUFY(si,altbuf)) {
+ SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+ SWAPI (SI_BUFY(si,1), SI_BUFY(si,2))
+
+ } else {
+ # Get line and interpolate onto output grid. If interpolation
+ # is not required merely copy data out. This code is set up
+ # to always use two buffers; in effect, there is one buffer of
+ # look ahead, even when Y[i] is integral. This means that we
+ # will go out of bounds by one line at the top of the image.
+ # This is handled by copying the last line.
+
+ ybavg = SI_BAVG(si,2)
+ nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+ if (new_y[i] <= nblks_y)
+ rawline = si_blmavgi (SI_IM(si), SI_FP(si), x1, x2,
+ new_y[i], SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si))
+
+ if (SI_INTERP(si,1) == NO) {
+ call amovi (Memi[rawline], Memi[SI_BUF(si,i)], npix)
+ } else if (SI_ORDER(si) == 0) {
+ call si_samplei (Memi[rawline], Memi[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ } else if (SI_ORDER(si) == -1) {
+ call si_maxi (Memi[rawline], nraw,
+ Memr[SI_GRID(si,1)], Memi[SI_BUF(si,i)], npix)
+ } else {
+ call aluii (Memi[rawline], Memi[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ }
+
+ SI_BUFY(si,i) = new_y[i]
+ }
+
+ SWAPI (altbuf, curbuf)
+ }
+
+ # We now have two line buffers straddling the output Y value,
+ # interpolated to the X grid of the output line. To complete the
+ # bilinear interpolation operation we take a weighted sum of the two
+ # lines. If the range from SI_BUFY(si,1) to SI_BUFY(si,2) is repeatedly
+ # interpolated in Y no additional i/o occurs and the linear
+ # interpolation operation (ALUI) does not have to be repeated (only the
+ # weighted sum is required). If the distance of Y from one of the
+ # buffers is zero then we do not even have to take a weighted sum.
+ # This is not unusual because we may be called with a magnification
+ # of 1.0 in Y.
+
+ weight_1 = 1.0 - (y - SI_BUFY(si,1))
+ weight_2 = 1.0 - weight_1
+
+ if (weight_1 < SI_TOL)
+ return (SI_BUF(si,2))
+ else if (weight_2 < SI_TOL || SI_ORDER(si) == 0)
+ return (SI_BUF(si,1))
+ else if (SI_ORDER(si) == -1) {
+ call amaxi (Memi[SI_BUF(si,1)], Memi[SI_BUF(si,2)],
+ Memi[OUTBUF(si)], npix)
+ return (OUTBUF(si))
+ } else {
+ call awsui (Memi[SI_BUF(si,1)], Memi[SI_BUF(si,2)],
+ Memi[OUTBUF(si)], npix, weight_1, weight_2)
+ return (OUTBUF(si))
+ }
+end
+
+
+# SI_BLMAVGI -- Get a line from a block averaged image of type integer.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on. If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blmavgi (im, fp, x1, x2, y, xbavg, ybavg, order)
+
+pointer im # input image
+pointer fp # fixpix structure
+int x1, x2 # range of x blocks to be read
+int y # y block to be read
+int xbavg, ybavg # X and Y block averaging factors
+int order # averaging option
+
+real sum
+int blkmax
+pointer sp, a, b
+int nblks_x, nblks_y, ncols, nlines, xoff, blk1, blk2, i, j, k
+int first_line, nlines_in_sum, npix, nfull_blks, count
+pointer xt_fpi()
+errchk xt_fpi
+
+begin
+ call smark (sp)
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ xoff = (x1 - 1) * xbavg + 1
+ npix = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1) - xoff + 1
+
+ if ((xbavg < 1) || (ybavg < 1))
+ call error (1, "si_blmavg: illegal block size")
+ else if (x1 < 1 || x2 > ncols)
+ call error (2, "si_blmavg: column index out of bounds")
+ else if ((xbavg == 1) && (ybavg == 1))
+ return (xt_fpi (fp, im, y, NULL) + xoff - 1)
+
+ nblks_x = (npix + xbavg-1) / xbavg
+ nblks_y = (nlines + ybavg-1) / ybavg
+
+ if (y < 1 || y > nblks_y)
+ call error (2, "si_blmavg: block number out of range")
+
+ if (ybavg > 1) {
+ call salloc (b, nblks_x, TY_LONG)
+ call aclrl (Meml[b], nblks_x)
+ nlines_in_sum = 0
+ }
+
+ # Read and accumulate all input lines in the block.
+ first_line = (y - 1) * ybavg + 1
+
+ do i = first_line, min (nlines, first_line + ybavg - 1) {
+ # Get line from input image.
+ a = xt_fpi (fp, im, i, NULL) + xoff - 1
+
+ # Block average line in X.
+ if (xbavg > 1) {
+ # First block average only the full blocks.
+ nfull_blks = npix / xbavg
+ if (order == -1) {
+ blk1 = a
+ do j = 1, nfull_blks {
+ blk2 = blk1 + xbavg
+ blkmax = Memi[blk1]
+ do k = blk1+1, blk2-1
+ blkmax = max (blkmax, Memi[k])
+ Memi[a+j-1] = blkmax
+ blk1 = blk2
+ }
+ } else
+ call abavi (Memi[a], Memi[a], nfull_blks, xbavg)
+
+ # Now average the final partial block, if any.
+ if (nfull_blks < nblks_x) {
+ if (order == -1) {
+ blkmax = Memi[blk1]
+ do k = blk1+1, a+npix-1
+ blkmax = max (blkmax, Memi[k])
+ Memi[a+j-1] = blkmax
+ } else {
+ sum = 0.0
+ count = 0
+ do j = nfull_blks * xbavg + 1, npix {
+ sum = sum + Memi[a+j-1]
+ count = count + 1
+ }
+ Memi[a+nblks_x-1] = sum / count
+ }
+ }
+ }
+
+ # Add line into block sum. Keep track of number of lines in sum
+ # so that we can compute block average later.
+
+ if (ybavg > 1) {
+ if (order == -1) {
+ do j = 0, nblks_x-1
+ Meml[b+j] = max (Meml[b+j], long (Memi[a+j]))
+ } else {
+ do j = 0, nblks_x-1
+ Meml[b+j] = Meml[b+j] + Memi[a+j]
+ nlines_in_sum = nlines_in_sum + 1
+ }
+ }
+ }
+
+ # Compute the block average in Y from the sum of all lines block
+ # averaged in X. Overwrite buffer A, the buffer returned by IMIO.
+ # This is kosher because the block averaged line is never longer
+ # than an input line.
+
+ if (ybavg > 1) {
+ if (order == -1) {
+ do i = 0, nblks_x-1
+ Memi[a+i] = Meml[b+i]
+ } else {
+ do i = 0, nblks_x-1
+ Memi[a+i] = Meml[b+i] / real(nlines_in_sum)
+ }
+ }
+
+ call sfree (sp)
+ return (a)
+end
+
+
+# SI_MAXI -- Resample a line via maximum value.
+
+procedure si_maxi (a, na, x, b, nb)
+
+int a[na] # input array
+int na # input size
+real x[nb] # sample grid
+int b[nb] # output arrays
+int nb # output size
+
+int i
+
+begin
+ do i = 1, nb
+ b[i] = max (a[int(x[i])], a[min(na,int(x[i]+1))])
+end
+
+
+# SIGM2R -- Get a line of type real from a scaled image. Block averaging is
+# done by a subprocedure; this procedure gets a line from a possibly block
+# averaged image and if necessary interpolates it to the grid points of the
+# output line.
+
+pointer procedure sigm2r (si, lineno)
+
+pointer si # pointer to SI descriptor
+int lineno
+
+pointer rawline, tempp, gp
+int i, new_y[2], tempi, curbuf, altbuf
+int nraw, npix, nblks_y, ybavg, x1, x2
+real x, y, weight_1, weight_2
+pointer si_blmavgr()
+errchk si_blmavgr
+
+begin
+ nraw = IM_LEN(SI_IM(si))
+ npix = SI_NPIX(si,1)
+
+ # Deterine the range of X (in pixels on the block averaged input image)
+ # required for the interpolator.
+
+ gp = SI_GRID(si,1)
+ x1 = SI_XOFF(si)
+ x = Memr[gp+npix-1]
+ x2 = x1 + int(x)
+ if (INTVAL(x))
+ x2 = x2 - 1
+ x2 = max (x1 + 1, x2)
+
+ gp = SI_GRID(si,2)
+ y = Memr[gp+lineno-1]
+
+ # The following is an optimization provided for the case when it is
+ # not necessary to interpolate in either X or Y. Block averaging is
+ # permitted.
+
+ if (SI_INTERP(si,1) == NO && SI_INTERP(si,2) == NO)
+ return (si_blmavgr (SI_IM(si), SI_FP(si), x1, x2, int(y),
+ SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si)))
+
+ # If we are interpolating in Y two buffers are required, one for each
+ # of the two input image lines required to interpolate in Y. The lines
+ # stored in these buffers are interpolated in X to the output grid but
+ # not in Y. Both buffers are not required if we are not interpolating
+ # in Y, but we use them anyhow to simplify the code.
+
+ if (SI_INIT(si) == YES) {
+ do i = 1, 2 {
+ if (SI_BUF(si,i) != NULL)
+ call mfree (SI_BUF(si,i), SI_TYBUF(si))
+ call malloc (SI_BUF(si,i), npix, TY_REAL)
+ SI_TYBUF(si) = TY_REAL
+ SI_BUFY(si,i) = NOTSET
+ }
+ if (OUTBUF(si) != NULL)
+ call mfree (OUTBUF(si), SI_TYBUF(si))
+ call malloc (OUTBUF(si), npix, TY_REAL)
+ SI_INIT(si) = NO
+ }
+
+ # If the Y value of the new line is not in range of the contents of the
+ # current line buffers, refill one or both buffers. To refill we must
+ # read a (possibly block averaged) input line and interpolate it onto
+ # the X grid. The X and Y values herein are in the coordinate system
+ # of the (possibly block averaged) input image.
+
+ new_y[1] = int(y)
+ new_y[2] = int(y) + 1
+
+ # Get the pair of lines whose integral Y values form an interval
+ # containing the fractional Y value of the output line. Sometimes the
+ # desired line will happen to be in the other buffer already, in which
+ # case we just have to swap buffers. Often the new line will be the
+ # current line, in which case nothing is done. This latter case occurs
+ # frequently when the magnification ratio is large.
+
+ curbuf = 1
+ altbuf = 2
+
+ do i = 1, 2 {
+ if (new_y[i] == SI_BUFY(si,i)) {
+ ;
+ } else if (new_y[i] == SI_BUFY(si,altbuf)) {
+ SWAPP (SI_BUF(si,1), SI_BUF(si,2))
+ SWAPI (SI_BUFY(si,1), SI_BUFY(si,2))
+
+ } else {
+ # Get line and interpolate onto output grid. If interpolation
+ # is not required merely copy data out. This code is set up
+ # to always use two buffers; in effect, there is one buffer of
+ # look ahead, even when Y[i] is integral. This means that we
+ # will go out of bounds by one line at the top of the image.
+ # This is handled by copying the last line.
+
+ ybavg = SI_BAVG(si,2)
+ nblks_y = (IM_LEN (SI_IM(si), 2) + ybavg-1) / ybavg
+ if (new_y[i] <= nblks_y)
+ rawline = si_blmavgr (SI_IM(si), SI_FP(si), x1, x2,
+ new_y[i], SI_BAVG(si,1), SI_BAVG(si,2), SI_ORDER(si))
+
+ if (SI_INTERP(si,1) == NO) {
+ call amovr (Memr[rawline], Memr[SI_BUF(si,i)], npix)
+ } else if (SI_ORDER(si) == 0) {
+ call si_sampler (Memr[rawline], Memr[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ } else if (SI_ORDER(si) == -1) {
+ call si_maxr (Memr[rawline], nraw,
+ Memr[SI_GRID(si,1)], Memr[SI_BUF(si,i)], npix)
+ } else {
+ call aluir (Memr[rawline], Memr[SI_BUF(si,i)],
+ Memr[SI_GRID(si,1)], npix)
+ }
+
+ SI_BUFY(si,i) = new_y[i]
+ }
+
+ SWAPI (altbuf, curbuf)
+ }
+
+ # We now have two line buffers straddling the output Y value,
+ # interpolated to the X grid of the output line. To complete the
+ # bilinear interpolation operation we take a weighted sum of the two
+ # lines. If the range from SI_BUFY(si,1) to SI_BUFY(si,2) is repeatedly
+ # interpolated in Y no additional i/o occurs and the linear
+ # interpolation operation (ALUI) does not have to be repeated (only the
+ # weighted sum is required). If the distance of Y from one of the
+ # buffers is zero then we do not even have to take a weighted sum.
+ # This is not unusual because we may be called with a magnification
+ # of 1.0 in Y.
+
+ weight_1 = 1.0 - (y - SI_BUFY(si,1))
+ weight_2 = 1.0 - weight_1
+
+ if (weight_1 < SI_TOL)
+ return (SI_BUF(si,2))
+ else if (weight_2 < SI_TOL || SI_ORDER(si) == 0)
+ return (SI_BUF(si,1))
+ else if (SI_ORDER(si) == -1) {
+ call amaxr (Memr[SI_BUF(si,1)], Memr[SI_BUF(si,2)],
+ Memr[OUTBUF(si)], npix)
+ return (OUTBUF(si))
+ } else {
+ call awsur (Memr[SI_BUF(si,1)], Memr[SI_BUF(si,2)],
+ Memr[OUTBUF(si)], npix, weight_1, weight_2)
+ return (OUTBUF(si))
+ }
+end
+
+
+# SI_BLMAVGR -- Get a line from a block averaged image of type short.
+# For example, block averaging by a factor of 2 means that pixels 1 and 2
+# are averaged to produce the first output pixel, 3 and 4 are averaged to
+# produce the second output pixel, and so on. If the length of an axis
+# is not an integral multiple of the block size then the last pixel in the
+# last block will be replicated to fill out the block; the average is still
+# defined even if a block is not full.
+
+pointer procedure si_blmavgr (im, fp, x1, x2, y, xbavg, ybavg, order)
+
+pointer im # input image
+pointer fp # fixpix structure
+int x1, x2 # range of x blocks to be read
+int y # y block to be read
+int xbavg, ybavg # X and Y block averaging factors
+int order # averaging option
+
+int nblks_x, nblks_y, ncols, nlines, xoff, blk1, blk2, i, j, k
+int first_line, nlines_in_sum, npix, nfull_blks, count
+real sum, blkmax
+pointer sp, a, b
+pointer xt_fpr()
+errchk xt_fpr
+
+begin
+ call smark (sp)
+
+ ncols = IM_LEN(im,1)
+ nlines = IM_LEN(im,2)
+ xoff = (x1 - 1) * xbavg + 1
+ npix = min (ncols, xoff + (x2 - x1 + 1) * xbavg - 1) - xoff + 1
+
+ if ((xbavg < 1) || (ybavg < 1))
+ call error (1, "si_blmavg: illegal block size")
+ else if (x1 < 1 || x2 > ncols)
+ call error (2, "si_blmavg: column index out of bounds")
+ else if ((xbavg == 1) && (ybavg == 1))
+ return (xt_fpr (fp, im, y, NULL) + xoff - 1)
+
+ nblks_x = (npix + xbavg-1) / xbavg
+ nblks_y = (nlines + ybavg-1) / ybavg
+
+ if (y < 1 || y > nblks_y)
+ call error (2, "si_blmavg: block number out of range")
+
+ call salloc (b, nblks_x, TY_REAL)
+
+ if (ybavg > 1) {
+ call aclrr (Memr[b], nblks_x)
+ nlines_in_sum = 0
+ }
+
+ # Read and accumulate all input lines in the block.
+ first_line = (y - 1) * ybavg + 1
+
+ do i = first_line, min (nlines, first_line + ybavg - 1) {
+ # Get line from input image.
+ a = xt_fpr (fp, im, i, NULL) + xoff - 1
+
+ # Block average line in X.
+ if (xbavg > 1) {
+ # First block average only the full blocks.
+ nfull_blks = npix / xbavg
+ if (order == -1) {
+ blk1 = a
+ do j = 1, nfull_blks {
+ blk2 = blk1 + xbavg
+ blkmax = Memr[blk1]
+ do k = blk1+1, blk2-1
+ blkmax = max (blkmax, Memr[k])
+ Memr[a+j-1] = blkmax
+ blk1 = blk2
+ }
+ } else
+ call abavr (Memr[a], Memr[a], nfull_blks, xbavg)
+
+ # Now average the final partial block, if any.
+ if (nfull_blks < nblks_x) {
+ if (order == -1) {
+ blkmax = Memr[blk1]
+ do k = blk1+1, a+npix-1
+ blkmax = max (blkmax, Memr[k])
+ Memr[a+j-1] = blkmax
+ } else {
+ sum = 0.0
+ count = 0
+ do j = nfull_blks * xbavg + 1, npix {
+ sum = sum + Memr[a+j-1]
+ count = count + 1
+ }
+ Memr[a+nblks_x-1] = sum / count
+ }
+ }
+ }
+
+ # Add line into block sum. Keep track of number of lines in sum
+ # so that we can compute block average later.
+ if (ybavg > 1) {
+ if (order == -1)
+ call amaxr (Memr[a], Memr[b], Memr[b], nblks_x)
+ else {
+ call aaddr (Memr[a], Memr[b], Memr[b], nblks_x)
+ nlines_in_sum = nlines_in_sum + 1
+ }
+ }
+ }
+
+ # Compute the block average in Y from the sum of all lines block
+ # averaged in X. Overwrite buffer A, the buffer returned by IMIO.
+ # This is kosher because the block averaged line is never longer
+ # than an input line.
+
+ if (ybavg > 1) {
+ if (order == -1)
+ call amovr (Memr[b], Memr[a], nblks_x)
+ else
+ call adivkr (Memr[b], real(nlines_in_sum), Memr[a], nblks_x)
+ }
+
+ call sfree (sp)
+ return (a)
+end
+
+
+# SI_MAXR -- Resample a line via maximum value.
+
+procedure si_maxr (a, na, x, b, nb)
+
+real a[na] # input array
+int na # input size
+real x[nb] # sample grid
+real b[nb] # output arrays
+int nb # output size
+
+int i
+
+begin
+ do i = 1, nb
+ b[i] = max (a[int(x[i])], a[min(na,int(x[i]+1))])
+end
generated by cgit (git 2.34.1) at 2025-09-20 11:21:24 -0400