From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- pkg/proto/vol/src/i2sun/cnvimage.x | 142 +++++++ pkg/proto/vol/src/i2sun/i2sun.h | 46 +++ pkg/proto/vol/src/i2sun/mkpkg | 27 ++ pkg/proto/vol/src/i2sun/sigln.x | 783 +++++++++++++++++++++++++++++++++++++ pkg/proto/vol/src/i2sun/t_i2sun.x | 240 ++++++++++++ pkg/proto/vol/src/i2sun/trsetup.x | 32 ++ pkg/proto/vol/src/i2sun/trulut.x | 128 ++++++ pkg/proto/vol/src/i2sun/x_i2sun.x | 4 + 8 files changed, 1402 insertions(+) create mode 100644 pkg/proto/vol/src/i2sun/cnvimage.x create mode 100644 pkg/proto/vol/src/i2sun/i2sun.h create mode 100644 pkg/proto/vol/src/i2sun/mkpkg create mode 100644 pkg/proto/vol/src/i2sun/sigln.x create mode 100644 pkg/proto/vol/src/i2sun/t_i2sun.x create mode 100644 pkg/proto/vol/src/i2sun/trsetup.x create mode 100644 pkg/proto/vol/src/i2sun/trulut.x create mode 100644 pkg/proto/vol/src/i2sun/x_i2sun.x (limited to 'pkg/proto/vol/src/i2sun') diff --git a/pkg/proto/vol/src/i2sun/cnvimage.x b/pkg/proto/vol/src/i2sun/cnvimage.x new file mode 100644 index 00000000..59bd4655 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/cnvimage.x @@ -0,0 +1,142 @@ +include +include +include "i2sun.h" + + +# CNV_IMAGE -- Read each line of the input image, apply ztransform, convert +# to rasterfile form, and write to rasterfile. + +procedure cnv_image (im, slice, tr, uptr, rfd) +pointer im # input image +int slice # current slice if n>2 image +pointer tr # spatial & greyscale transform structure +pointer uptr # pointer to user-specified transfer lut +pointer rfd # output rasterfile + +real z1, z2, rz1, rz2 +int ztrans, row, xblk, yblk, ocols, olines +real px1, px2, py1, py2 # image coords in fractional image pixels +pointer sp, si, bufptr, out, rtemp, packed +short z1_s, z2_s, rz1_s, rz2_s +bool unitary_greyscale_transformation + +bool fp_equalr() +pointer siglns(), siglnr(), sigln_setup() +errchk siglns(), siglnr(), sigln_setup() + +begin + # Set up for scaled image input. + px1 = 1 + px2 = IM_LEN(im,COL) + py1 = 1 + py2 = IM_LEN(im,LINE) + ocols = TR_XE(tr) - TR_XS(tr) + 1 + olines = TR_YE(tr) - TR_YS(tr) + 1 + + # Round odd-numbered numbers of columns up due to restrictions on + # IRAF binary byte i/o (number of bytes of image data must match + # that specified in rasterfile header). + if (mod (ocols, 2) == 1) + ocols = ocols + 1 + + xblk = INDEFI + yblk = INDEFI + si = sigln_setup (im, px1,px2,ocols,xblk, py1,py2,olines,yblk, + TR_ORDER(tr)) + + # Type short pixels are treated as a special case to minimize vector + # operations for such images (which are common). If unity mapping is + # employed the data is simply copied, i.e., floor ceiling constraints + # are not applied. This is very fast and will produce a contoured + # image on the display which will be adequate for some applications. + + z1 = TR_Z1(tr) + z2 = TR_Z2(tr) + ztrans = TR_ZTRANS(tr) + rz1 = COLORSTART + rz2 = COLOREND + if (ztrans == Z_UNITARY) { + unitary_greyscale_transformation = true + } else if (ztrans == Z_LINEAR) { + unitary_greyscale_transformation = + ((fp_equalr(z1,rz1) && fp_equalr(z2,rz2)) || fp_equalr(z1,z2)) + } else + unitary_greyscale_transformation = false + + if (IM_PIXTYPE(im) == TY_SHORT && ztrans != Z_LOG) { + + call smark (sp) + call salloc (out, ocols, TY_SHORT) + call salloc (packed, ocols, TY_CHAR) + z1_s = z1; z2_s = z2 + + for (row=olines; row >= 1; row=row-1) { + bufptr = siglns (si, row, TR_SLICEAXIS(tr), slice) + + if (unitary_greyscale_transformation) { + call amovs (Mems[bufptr], Mems[out], ocols) + } else if (ztrans == Z_USER) { + rz1_s = U_Z1; rz2_s = U_Z2 + call amaps (Mems[bufptr], Mems[out], ocols, z1_s, z2_s, + rz1_s, rz2_s) + call aluts (Mems[out], Mems[out], ocols, Mems[uptr]) + } else { + rz1_s = rz1; rz2_s = rz2 + call amaps (Mems[bufptr], Mems[out], ocols, z1_s, z2_s, + rz1_s, rz2_s) + } + + # Pack to unsigned byte and write to file. + call achtsc (Mems[out], Memc[packed], ocols) + call chrpak (Memc[packed], 1, Memc[packed], 1, ocols) + call write (rfd, Memc[packed], ocols / SZB_CHAR) + } + call sfree (sp) + + } else if (ztrans == Z_USER) { + call smark (sp) + call salloc (rtemp, ocols, TY_REAL) + call salloc (out, ocols, TY_SHORT) + call salloc (packed, ocols, TY_CHAR) + + for (row=olines; row >= 1; row=row-1) { + bufptr = siglnr (si, row, TR_SLICEAXIS(tr), slice) + call amapr (Memr[bufptr], Memr[rtemp], ocols, z1, z2, + real(U_Z1), real(U_Z2)) + call achtrs (Memr[rtemp], Mems[out], ocols) + call aluts (Mems[out], Mems[out], ocols, Mems[uptr]) + call achtsc (Mems[out], Memc[packed], ocols) + call chrpak (Memc[packed], 1, Memc[packed], 1, ocols) + call write (rfd, Memc[packed], ocols / SZB_CHAR) + } + call sfree (sp) + + } else { + call smark (sp) + call salloc (rtemp, ocols, TY_REAL) + call salloc (packed, ocols, TY_CHAR) + + for (row=olines; row >= 1; row=row-1) { + bufptr = siglnr (si, row, TR_SLICEAXIS(tr), slice) + + if (unitary_greyscale_transformation) { + call amovr (Memr[bufptr], Memr[rtemp], ocols) + } else if (ztrans == Z_LOG) { + call amapr (Memr[bufptr], Memr[rtemp], ocols, + z1, z2, 1.0, 10.0 ** MAXLOG) + call alogr (Memr[rtemp], Memr[rtemp], ocols) + call amapr (Memr[rtemp], Memr[rtemp], ocols, + 1.0, real(MAXLOG), rz1, rz2) + } else + call amapr (Memr[bufptr], Memr[rtemp], ocols, z1, z2, + rz1, rz2) + call achtrc (Memr[rtemp], Memc[packed], ocols) + call chrpak (Memc[packed], 1, Memc[packed], 1, ocols) + call write (rfd, Memc[packed], ocols / SZB_CHAR) + } + call sfree (sp) + } + + call sigln_free (si) +end + diff --git a/pkg/proto/vol/src/i2sun/i2sun.h b/pkg/proto/vol/src/i2sun/i2sun.h new file mode 100644 index 00000000..73f2ea3f --- /dev/null +++ b/pkg/proto/vol/src/i2sun/i2sun.h @@ -0,0 +1,46 @@ +# I2SUNRAS.H -- Include file for IRAF to Sun rasterfile program i2sunras. + +define COL 1 +define LINE 2 +define BAND 3 +define Z_LINEAR 1 # linear ztransform +define Z_LOG 2 # log ztransform +define Z_UNITARY 3 # no ztransform +define Z_USER 4 # user-specified transform +define U_MAXPTS 4096 # max user-specified lut pairs (DISPLAY) +define U_Z1 0 # base user-specified transfer val +define U_Z2 4095 # upper user-specified transfer val +define MAXLOG 3 # if log, map to 1:10**MAXLOG b4 log10 +define DSP_MIN 0 # minimum display pixel value +define DSP_MAX 255 # maximum display pixel value +define RAS_HDR_INTS 8 # SunOS4.0 and earlier +define RMT_NONE 0 # SunOS4.0 and earlier +define RMT_EQUAL_RGB 1 # SunOS4.0 and earlier +define RMT_STANDARD 1 # SunOS4.0 and earlier +define RAS_MAGIC 1504078485 # SunOS4.0 and earlier +define NGREY 256 # SunOS4.0 and earlier, 8bit fb +define COLORSTART 1 # IMTOOL +define COLOREND 200 # IMTOOL +define COLORRANGE 200 # IMTOOL +define WHITE (NGREY-1) # IMTOOL +define BLACK 0 # IMTOOL +define NBITS_FB 8 +define wrapup_ 91 + +# Spatial and greyscale transformation structure. +define LEN_TR 20 +define TR_ZTRANS Memi[$1] # Greyscale transformation. +define TR_Z1 Memr[P2R($1+1)] # Minimum data z-value +define TR_Z2 Memr[P2R($1+2)] # Maximum data z-value +define TR_XSIZE Memi[$1+3] # Output rasterfile size in x +define TR_YSIZE Memi[$1+4] # Output rasterfile size in y +define TR_XMAG Memr[P2R($1+5)] # Magnification factor in x +define TR_YMAG Memr[P2R($1+6)] # Magnification factor in y +define TR_ORDER Memi[$1+7] # Interpolation order +define TR_XS Memi[$1+8] # Starting output x pixel index +define TR_XE Memi[$1+9] # Ending output x pixel index +define TR_YS Memi[$1+10] # Starting output y pixel index +define TR_YE Memi[$1+11] # Ending output y pixel index +define TR_SLICEAXIS Memi[$1+12] # Slice or frame axis when ndim>2 +define TR_SWAPBYTES Memb[$1+13] # Swap output bytes? +# # Reserved space diff --git a/pkg/proto/vol/src/i2sun/mkpkg b/pkg/proto/vol/src/i2sun/mkpkg new file mode 100644 index 00000000..b1a8c4f4 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/mkpkg @@ -0,0 +1,27 @@ +# Library for the I2SUN task. + +$checkout libpkg.a ../ +$update libpkg.a +$checkin libpkg.a ../ +$exit + +libpkg.a: + t_i2sun.x i2sun.h + trulut.x i2sun.h + trsetup.x i2sun.h + cnvimage.x i2sun.h + sigln.x + ; + +dbx: + $set XFLAGS = "-c -g -F -q" + $set LFLAGS = "-g -q" + $omake x_i2sun.x + $omake t_i2sun.x + $omake trulut.x + $omake trsetup.x + $omake cnvimage.x + $omake sigln.x + $link x_i2sun.o t_i2sun.o trulut.o trsetup.o cnvimage.o sigln.o \ + -o xx_i2sun.e + ; diff --git a/pkg/proto/vol/src/i2sun/sigln.x b/pkg/proto/vol/src/i2sun/sigln.x new file mode 100644 index 00000000..9d763b3f --- /dev/null +++ b/pkg/proto/vol/src/i2sun/sigln.x @@ -0,0 +1,783 @@ +include +include + +.help sigl2, sigl2_setup +.nf ___________________________________________________________________________ +SIGLN -- Get a line from a spatially scaled image of any dimensionality. +This procedure works like the regular IMIO get line procedure, but rescales +the input image in 1 or two axes upon input (for a resulting 2d output image). +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 = sigln_setup (im, x1,x2,nx,xblk, y1,y2,ny,yblk, order) + sigln_free (si) + ptr = sigln[sr] (si, linenumber) + +SIGLN_SETUP must be called to set up the transformations after mapping the +image and before performing any scaled i/o to the image. SIGLN_FREE must be +called when finished to return buffer space. +.endhelp ______________________________________________________________________ + +# Scaled image descriptor for 2-dim images + +define SI_LEN 16 +define SI_MAXDIM 2 # 2 dimensions of spatial scaling +define SI_NBUFS 3 # nbuffers used by SIGLN + +define SI_IM Memi[$1] # pointer to input image header +define SI_GRID Memi[$1+1+$2-1] # pointer to array of X coords +define SI_NPIX Memi[$1+3+$2-1] # number of X coords +define SI_BAVG Memi[$1+5+$2-1] # X block averaging factor +define SI_INTERP Memi[$1+7+$2-1] # interpolate X axis +define SI_BUF Memi[$1+9+$2-1] # line buffers +define SI_ORDER Memi[$1+12] # interpolator order, 0 or 1 +define SI_TYBUF Memi[$1+13] # buffer type +define SI_XOFF Memi[$1+14] # offset in input image to first X +define SI_INIT Memi[$1+15] # 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) + +# SIGLN_SETUP -- Set up the spatial transformation for SIGLN[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 sigln_setup (im, px1,px2,nx,xblk, py1,py2,ny,yblk, order) + +pointer im # the input image +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 + +begin + iferr (call calloc (si, SI_LEN, TY_STRUCT)) + call erract (EA_FATAL) + + SI_IM(si) = im + 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) && (blksize[i] != INDEFI)) { + 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) + } 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]) + 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]) + start = p1[1] - int (p1[i]) + 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) + } else { + do j = 0, npts[i]-1 + Memr[gp+j] = start + (j * tau[i]) + } + } + + return (si) +end + + +# SIGLN_FREE -- Free storage associated with an image opened for scaled +# input. This does not close and unmap the image. + +procedure sigln_free (si) + +pointer si +int i + +begin + # Free SIGLN 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 + + +# SIGLNS -- 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 siglns (si, lineno, slice_axis, slice) + +pointer si # pointer to SI descriptor +int lineno +int slice_axis # axis from which to slice section for ndim>2 images +int slice # current slice index + +pointer rawline, tempp, gp +int i, buf_y[2], new_y[2], tempi, curbuf, altbuf +int npix, nblks_y, ybavg, x1, x2 +real x, y, weight_1, weight_2 +pointer si_blkavgs() +errchk si_blkavgs + +begin + 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_blkavgs (SI_IM(si), x1, x2, int(y), + SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice)) + + # 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 + buf_y[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] == buf_y[i]) { + ; + } else if (new_y[i] == buf_y[altbuf]) { + SWAPP (SI_BUF(si,1), SI_BUF(si,2)) + SWAPI (buf_y[1], buf_y[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_blkavgs (SI_IM(si), x1, x2, new_y[i], + SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice) + + 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 { + call aluis (Mems[rawline], Mems[SI_BUF(si,i)], + Memr[SI_GRID(si,1)], npix) + } + + buf_y[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 buf_y[1] to buf_y[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 - buf_y[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 { + 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_BLKAVGS -- 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_blkavgs (im, x1, x2, y, xbavg, ybavg, slice_axis, slice) + +pointer im # input image +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 slice_axis # slice dimension if ndim>2 image +int slice # slice if ndim>2 image + +short temp_s +int nblks_x, nblks_y, ncols, nlines, xoff, i, j +int first_line, nlines_in_sum, npix, nfull_blks, count +long vs[IM_MAXDIM], ve[IM_MAXDIM] +real sum +pointer sp, a, b +pointer imgs2s(), imgs3s(), imggss() +errchk imgs2s, imgs3s, imggss + +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) + + if ((xbavg < 1) || (ybavg < 1)) + call error (1, "si_blkavg: illegal block size") + else if (x1 < 1 || x2 > ncols) + call error (2, "si_blkavg: column index out of bounds") + else if ((xbavg == 1) && (ybavg == 1)) + if (IM_NDIM(im) == 2) + return (imgs2s (im, xoff, xoff + npix - 1, y, y)) + else if (IM_NDIM(im) == 3) + return (imgs3s (im, xoff, xoff + npix - 1, y, y, slice, slice)) + else { + call amovkl (long(1), vs, IM_MAXDIM) + call amovkl (long(1), ve, IM_MAXDIM) + vs[1] = xoff + ve[1] = xoff + npix - 1 + vs[2] = y + ve[2] = y + vs[slice_axis] = slice + ve[slice_axis] = slice + return (imggss (im, vs, ve, 2)) + } + + nblks_x = (npix + xbavg-1) / xbavg + nblks_y = (nlines + ybavg-1) / ybavg + + if (y < 1 || y > nblks_y) + call error (2, "si_blkavg: block number out of range") + + call salloc (b, nblks_x, TY_SHORT) + + if (ybavg > 1) { + call aclrs (Mems[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. + if (IM_NDIM(im) == 2) + a = imgs2s (im, xoff, xoff + npix - 1, i, i) + else if (IM_NDIM(im) == 3) + a = imgs3s (im, xoff, xoff + npix - 1, i, i, slice, slice) + else { + call amovkl (long(1), vs, IM_MAXDIM) + call amovkl (long(1), ve, IM_MAXDIM) + vs[1] = xoff + ve[1] = xoff + npix - 1 + vs[2] = i + ve[2] = i + vs[slice_axis] = slice + ve[slice_axis] = slice + return (imggss (im, vs, ve, 2)) + } + + # Block average line in X. + if (xbavg > 1) { + # First block average only the full blocks. + nfull_blks = npix / xbavg + call abavs (Mems[a], Mems[a], nfull_blks, xbavg) + + # Now average the final partial block, if any. + if (nfull_blks < nblks_x) { + 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) { + call aadds (Mems[a], Mems[b], Mems[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) { + temp_s = nlines_in_sum + call adivks (Mems[b], temp_s, Mems[a], nblks_x) + } + + call sfree (sp) + return (a) +end + + +# SI_SAMPLES -- Resample a line via nearest neighbor, rather than linear +# interpolation (ALUI). The calling sequence is the same as for ALUIS. + +procedure si_samples (a, b, x, npix) + +short a[ARB], b[ARB] # input, output data arrays +real x[ARB] # sample grid +int npix, i + +begin + do i = 1, npix + b[i] = a[int(x[i])] +end + + +# SIGLNR -- 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 siglnr (si, lineno, slice_axis, slice) + +pointer si # pointer to SI descriptor +int lineno +int slice_axis # axis from which to slice section if ndim>2 +int slice # current slice index + +pointer rawline, tempp, gp +int i, buf_y[2], new_y[2], tempi, curbuf, altbuf +int npix, nblks_y, ybavg, x1, x2 +real x, y, weight_1, weight_2 +pointer si_blkavgr() +errchk si_blkavgr + +begin + 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_blkavgr (SI_IM(si), x1, x2, int(y), + SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice)) + + # 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 + buf_y[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] == buf_y[i]) { + ; + } else if (new_y[i] == buf_y[altbuf]) { + SWAPP (SI_BUF(si,1), SI_BUF(si,2)) + SWAPI (buf_y[1], buf_y[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_blkavgr (SI_IM(si), x1, x2, new_y[i], + SI_BAVG(si,1), SI_BAVG(si,2), slice_axis, slice) + + 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 { + call aluir (Memr[rawline], Memr[SI_BUF(si,i)], + Memr[SI_GRID(si,1)], npix) + } + + buf_y[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 buf_y[1] to buf_y[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 - buf_y[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 { + 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_BLKAVGR -- 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_blkavgr (im, x1, x2, y, xbavg, ybavg, slice_axis, slice) + +pointer im # input image +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 slice_axis # axis from which to slice section if ndim>2 +int slice # current slice + +int nblks_x, nblks_y, ncols, nlines, xoff, i, j +int first_line, nlines_in_sum, npix, nfull_blks, count +long vs[IM_MAXDIM], ve[IM_MAXDIM] +real sum +pointer sp, a, b +pointer imgs2r(), imgs3r(), imggsr() +errchk imgs2r, imgs3r, imggsr() + +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) + + if ((xbavg < 1) || (ybavg < 1)) + call error (1, "si_blkavg: illegal block size") + else if (x1 < 1 || x2 > ncols) + call error (2, "si_blkavg: column index out of bounds") + else if ((xbavg == 1) && (ybavg == 1)) + if (IM_NDIM(im) == 2) + return (imgs2r (im, xoff, xoff + npix - 1, y, y)) + else if (IM_NDIM(im) == 3) + return (imgs3r (im, xoff, xoff + npix - 1, y, y, slice, slice)) + else { + call amovkl (long(1), vs, IM_MAXDIM) + call amovkl (long(1), ve, IM_MAXDIM) + vs[1] = xoff + ve[1] = xoff + npix - 1 + vs[2] = y + ve[2] = y + vs[slice_axis] = slice + ve[slice_axis] = slice + return (imggsr (im, vs, ve, 2)) + } + + nblks_x = (npix + xbavg-1) / xbavg + nblks_y = (nlines + ybavg-1) / ybavg + + if (y < 1 || y > nblks_y) + call error (2, "si_blkavg: 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. + if (IM_NDIM(im) == 2) + a = imgs2r (im, xoff, xoff + npix - 1, i, i) + else if (IM_NDIM(im) == 3) + a = imgs3r (im, xoff, xoff + npix - 1, i, i, slice, slice) + else { + call amovkl (long(1), vs, IM_MAXDIM) + call amovkl (long(1), ve, IM_MAXDIM) + vs[1] = xoff + ve[1] = xoff + npix - 1 + vs[2] = i + ve[2] = i + vs[slice_axis] = slice + ve[slice_axis] = slice + return (imggsr (im, vs, ve, 2)) + } + + # Block average line in X. + if (xbavg > 1) { + # First block average only the full blocks. + nfull_blks = npix / xbavg + call abavr (Memr[a], Memr[a], nfull_blks, xbavg) + + # Now average the final partial block, if any. + if (nfull_blks < nblks_x) { + 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) { + 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) + call adivkr (Memr[b], real(nlines_in_sum), Memr[a], nblks_x) + + call sfree (sp) + return (a) +end + + +# SI_SAMPLER -- Resample a line via nearest neighbor, rather than linear +# interpolation (ALUI). The calling sequence is the same as for ALUIR. + +procedure si_sampler (a, b, x, npix) + +real a[ARB], b[ARB] # input, output data arrays +real x[ARB] # sample grid +int npix, i + +begin + do i = 1, npix + b[i] = a[int(x[i])] +end diff --git a/pkg/proto/vol/src/i2sun/t_i2sun.x b/pkg/proto/vol/src/i2sun/t_i2sun.x new file mode 100644 index 00000000..ab8119b7 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/t_i2sun.x @@ -0,0 +1,240 @@ +include +include +include +include "i2sun.h" + + +# I2SUN -- IRAF to Sun Rasterfile: convert either a list of IRAF images +# or all slices from a specified axis of a dimension>2 image into a series +# of Sun rasterfiles. This format-specific task is primarily used to make +# movies in the absence of a portable movie/filmloop utility, if a +# Sun-specific movie task is available. +# ** The format of the output Sun rasterfiles is hard-coded into this task, +# ** and thus could diverge from a future Sun format; we do not want to link +# ** with Sun libraries, as this task should be runnable on other machines. + +procedure t_i2sun + +pointer sp, tr, input, im, rfnames, clutfile, transform, cur_rf +pointer ulutfile, ulut, colormap, pk_colormap, lut +int list, lfd, rfd, nslices, stat, nimages +int rheader[RAS_HDR_INTS], ras_maptype, ras_maplength, frame, slice, i, j +short lut1, lut2 +bool use_clut, make_map + +pointer immap() +int open(), access(), clgeti(), imtopenp(), imtlen(), imtgetim(), read() +real clgetr() +bool streq(), clgetb() + +errchk open() + +begin + call smark (sp) + call salloc (tr, LEN_TR, TY_STRUCT) + call salloc (input, SZ_FNAME, TY_CHAR) + call salloc (rfnames, SZ_FNAME, TY_CHAR) + call salloc (clutfile, SZ_FNAME, TY_CHAR) + call salloc (cur_rf, SZ_FNAME, TY_CHAR) + call salloc (transform, SZ_LINE, TY_CHAR) + call salloc (pk_colormap, NGREY*3, TY_CHAR) + call salloc (colormap, NGREY*3, TY_SHORT) + lut = NULL + im = NULL + + # Input parameters. + list = imtopenp ("input") + call clgstr ("output", Memc[rfnames], SZ_FNAME) + call clgstr ("clutfile", Memc[clutfile], SZ_FNAME) + call clgstr ("ztrans", Memc[transform], SZ_LINE) + TR_Z1(tr) = clgetr ("z1") + TR_Z2(tr) = clgetr ("z2") + TR_ZTRANS(tr) = Z_LINEAR + if (streq (Memc[transform], "log")) + TR_ZTRANS(tr) = Z_LOG + else if (streq (Memc[transform], "none")) + TR_ZTRANS(tr) = Z_UNITARY + else if (streq (Memc[transform], "user")) { + + # Get user-specified transfer lookup table. + TR_ZTRANS(tr) = Z_USER + call salloc (ulutfile, SZ_FNAME, TY_CHAR) + call clgstr ("ulutfile", Memc[ulutfile], SZ_FNAME) + + # Borrowed from DISPLAY; mallocs storage for ulut: + call tr_ulut (Memc[ulutfile], TR_Z1(tr), TR_Z2(tr), ulut) + } + TR_XSIZE(tr) = clgeti ("xsize") + TR_YSIZE(tr) = clgeti ("ysize") + TR_ORDER(tr) = clgeti ("order") + TR_XMAG(tr) = clgetr ("xmag") + TR_YMAG(tr) = clgetr ("ymag") + + # Get input image axes to map to output frames. At present we + # can only traverse one slice axis. + TR_SLICEAXIS(tr) = clgeti ("sliceaxis") + + # Swap bytes in output rasterfile? (useful when I2SUN run on VAX etc.) + TR_SWAPBYTES(tr) = clgetb ("swap") + + # Check if there are no images. + nimages = imtlen (list) + if (nimages == 0) { + call eprintf (0, "No input images to convert") + goto wrapup_ + } + + # Open color lookup table file (an existing Sun rasterfile at present) + if (access (Memc[clutfile], READ_ONLY, BINARY_FILE) == YES) { + lfd = open (Memc[clutfile], READ_ONLY, BINARY_FILE) + use_clut = true + } else + use_clut = false + + # Read color lookup table. + make_map = false + if (use_clut) { + # Only the color table is used from the rasterfile; ignore all else. + stat = read (lfd, rheader, RAS_HDR_INTS * SZB_CHAR) + if (stat != RAS_HDR_INTS * SZB_CHAR) { + call eprintf ("Error reading header from file `%s'\n") + call pargstr (Memc[clutfile]) + goto wrapup_ + } + if (rheader[1] != RAS_MAGIC) { + call eprintf ("File `%s' not a valid Sun rasterfile\n") + call pargstr (Memc[clutfile]) + goto wrapup_ + } + ras_maptype = rheader[7] + ras_maplength = rheader[8] + if (ras_maptype != RMT_NONE && ras_maplength > 0) { + stat = read (lfd, Memc[colormap], ras_maplength / SZB_CHAR) + if (stat != ras_maplength / SZB_CHAR) { + call eprintf ("Error reading colormap from %s\n") + call pargstr (Memc[clutfile]) + goto wrapup_ + } + # Colormap was already packed on disk. + call achtsc (Mems[colormap], Memc[pk_colormap], ras_maplength) + } else { + make_map = true + call eprintf ("Invalid colormap in %s; using greyscale\n") + call pargstr (Memc[clutfile]) + } + } else + make_map = true + + if (make_map) { + # Construct a greyscale colormap of same range as IMTOOL. + ras_maptype = RMT_EQUAL_RGB + ras_maplength = NGREY * 3 + do i = 1, 3 { + Mems[colormap+(i-1)*NGREY] = WHITE + do j = COLORSTART+1, COLOREND + Mems[colormap+j-1+(i-1)*NGREY] = j * (WHITE+1) / + NGREY + Mems[colormap+COLOREND-1+1+(i-1)*NGREY] = WHITE + do j = COLOREND+2, NGREY + Mems[colormap+j-1+(i-1)*NGREY] = BLACK + } + call achtsc (Mems[colormap], Memc[pk_colormap], ras_maplength) + + # Pack to byte stream. + call chrpak (Memc[pk_colormap], 1, Memc[pk_colormap], 1, + ras_maplength) + } + + # For each IRAF image or band, construct and dispose of a rasterfile. + frame = 0 + while (imtgetim (list, Memc[input], SZ_FNAME) != EOF) { + im = immap (Memc[input], READ_ONLY, 0) + if (IM_NDIM(im) > 2 && TR_SLICEAXIS(tr) > IM_NDIM(im)) { + call eprintf ("Specified slice axis invalid for image %s\n") + call pargstr (Memc[input]) + goto wrapup_ + } + nslices = IM_LEN(im, TR_SLICEAXIS(tr)) + if (nslices < 1) + nslices = 1 + + # Set up spatial transformation (technically, could be different + # for each input image). + call tr_setup (im, tr) + + # We assume that if any n>2 images are present, the user wants + # all bands dumped out. + do slice = 1, nslices { + + # Construct next rasterfile name and open file; works in + # 'append' mode, next higher available frame number. + call sprintf (Memc[cur_rf], SZ_FNAME, Memc[rfnames]) + call pargi (frame) + while (access (Memc[cur_rf], READ_ONLY, BINARY_FILE) == YES) { + frame = frame + 1 + call sprintf (Memc[cur_rf], SZ_FNAME, Memc[rfnames]) + call pargi (frame) + } + iferr (rfd = open (Memc[cur_rf], NEW_FILE, BINARY_FILE)) { + call eprintf ("Cannot open output rasterfile `%s'\n") + call pargstr (Memc[cur_rf]) + goto wrapup_ + } + frame = frame + 1 + + # Write header to rasterfile: + rheader[1] = RAS_MAGIC + rheader[2] = TR_XE(tr) - TR_XS(tr) + 1 + rheader[3] = TR_YE(tr) - TR_YS(tr) + 1 + rheader[4] = NBITS_FB + rheader[5] = rheader[2] * rheader[3] + rheader[6] = RMT_STANDARD + rheader[7] = ras_maptype + rheader[8] = ras_maplength + if (TR_SWAPBYTES(tr)) + call bswap4 (rheader, 1, rheader, 1, RAS_HDR_INTS*4) + call write (rfd, rheader, RAS_HDR_INTS * SZB_CHAR) + + # Write colormap to rasterfile. + call write (rfd, Memc[pk_colormap], ras_maplength / SZB_CHAR) + + # Verify user-specified transfer function parameters. + if (TR_ZTRANS(tr) == Z_USER) { + call alims (Mems[ulut], U_MAXPTS, lut1, lut2) + if (lut2 < short(DSP_MIN) || lut1 > short(DSP_MAX)) { + call eprintf ("User specified greyscales <> range\n") + call eprintf ("ulut1=%D, dmin=%D; ulut2=%D, dmax=%D\n") + call pargi (lut1) + call pargi (DSP_MIN) + call pargi (lut2) + call pargi (DSP_MAX) + } + if (!IS_INDEF(TR_Z1(tr)) && !IS_INDEF(TR_Z2(tr)) && + TR_Z2(tr) < IM_MIN(im) || TR_Z1(tr) > IM_MAX(im)) { + call eprintf ("User specified intensities <> range\n") + call eprintf ("z1=%g, im_min=%g; z2=%g, im_max=%g\n") + call pargr (TR_Z1(tr)) + call pargr (IM_MIN(im)) + call pargr (TR_Z2(tr)) + call pargr (IM_MAX(im)) + call eprintf ("continuing anyway.\n") + } + } + + # Read image pixels and write to rasterfile. + call cnv_image (im, slice, tr, ulut, rfd) + + call close (rfd) + } + call imunmap (im) + } + +wrapup_ + if (im != NULL) + call imunmap(im) + call imtclose (list) + call close (rfd) + call sfree (sp) + if (ulut != NULL) + call mfree (ulut, TY_SHORT) +end diff --git a/pkg/proto/vol/src/i2sun/trsetup.x b/pkg/proto/vol/src/i2sun/trsetup.x new file mode 100644 index 00000000..1b14afb2 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/trsetup.x @@ -0,0 +1,32 @@ +include +include "i2sun.h" + + +# TR_SETUP -- Set up spatial transformation parameters. + +procedure tr_setup (im, tr) + +pointer im # An input image descriptor +pointer tr # Transformation structure + +int ncols, nlines + +begin + ncols = IM_LEN(im,COL) + nlines = IM_LEN(im,LINE) + + # Determine output raster dimensions. + TR_XS(tr) = 1 + TR_XE(tr) = ncols + if (!IS_INDEFI(TR_XSIZE(tr))) + TR_XE(tr) = max (1, TR_XSIZE(tr)) + else if (TR_XMAG(tr) != 1.0) + TR_XE(tr) = max (1, ncols * int(TR_XMAG(tr))) + + TR_YS(tr) = 1 + TR_YE(tr) = nlines + if (!IS_INDEFI(TR_YSIZE(tr))) + TR_YE(tr) = max (1, TR_YSIZE(tr)) + else if (TR_YMAG(tr) != 1.0) + TR_YE(tr) = max (1, nlines * int(TR_YMAG(tr))) +end diff --git a/pkg/proto/vol/src/i2sun/trulut.x b/pkg/proto/vol/src/i2sun/trulut.x new file mode 100644 index 00000000..4787b9b3 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/trulut.x @@ -0,0 +1,128 @@ +include +include +include "i2sun.h" + +# TR_ULUT -- Generates a look up table from data supplied by user. The +# data is read from a two column text file of intensity, greyscale values. +# The input data are sorted, then mapped to the x range [0-4095]. A +# piecewise linear look up table of 4096 values is then constructed from +# the (x,y) pairs given. A pointer to the look up table, as well as the z1 +# and z2 intensity endpoints, is returned. + +procedure tr_ulut (fname, z1, z2, lut) + +char fname[SZ_FNAME] # Name of file with intensity, greyscale values +real z1 # Intensity mapped to minimum gs value +real z2 # Intensity mapped to maximum gs value +pointer lut # Look up table - pointer is returned + +pointer sp, x, y +int nvalues, i, j, x1, x2, y1 +real delta_gs, delta_xv, slope +errchk ds_rlut, ds_sort, malloc + +begin + call smark (sp) + call salloc (x, U_MAXPTS, TY_REAL) + call salloc (y, U_MAXPTS, TY_REAL) + + # Read intensities and greyscales from the user's input file. The + # intensity range is then mapped into a standard range and the + # values sorted. + + call ds_rlut (fname, Memr[x], Memr[y], nvalues) + call alimr (Memr[x], nvalues, z1, z2) + call amapr (Memr[x], Memr[x], nvalues, z1, z2, real(U_Z1), real(U_Z2)) + call ds_sort (Memr[x], Memr[y], nvalues) + + # Fill lut in straight line segments - piecewise linear + call malloc (lut, U_MAXPTS, TY_SHORT) + do i = 1, nvalues-1 { + delta_gs = Memr[y+i] - Memr[y+i-1] + delta_xv = Memr[x+i] - Memr[x+i-1] + slope = delta_gs / delta_xv + x1 = int (Memr[x+i-1]) + x2 = int (Memr[x+i]) + y1 = int (Memr[y+i-1]) + do j = x1, x2 + Mems[lut+j] = y1 + slope * (j-x1) + } + Mems[lut+U_MAXPTS-1] = y1 + (slope * U_Z2) + + call sfree (sp) +end + + +# DS_RLUT -- Read text file of x, y, values. + +procedure ds_rlut (utab, x, y, nvalues) + +char utab[SZ_FNAME] # Name of list file +real x[U_MAXPTS] # Array of x values, filled on return +real y[U_MAXPTS] # Array of y values, filled on return +int nvalues # Number of values in x, y vectors - returned + +int n, fd +pointer sp, lbuf, ip +real xval, yval +int getline(), open() +errchk open, sscan, getline, salloc + +begin + call smark (sp) + call salloc (lbuf, SZ_LINE, TY_CHAR) + + iferr (fd = open (utab, READ_ONLY, TEXT_FILE)) + call error (1, "Error opening user lookup table") + + n = 0 + while (getline (fd, Memc[lbuf]) != EOF) { + # Skip comment lines and blank lines. + if (Memc[lbuf] == '#') + next + for (ip=lbuf; IS_WHITE(Memc[ip]); ip=ip+1) + ; + if (Memc[ip] == '\n' || Memc[ip] == EOS) + next + + # Decode the points to be plotted. + call sscan (Memc[ip]) + call gargr (xval) + call gargr (yval) + + n = n + 1 + if (n > U_MAXPTS) + call error (2, + "Intensity transformation table cannot exceed 4096 values") + + x[n] = xval + y[n] = yval + } + + nvalues = n + call close (fd) + call sfree (sp) +end + + +# DS_SORT -- Bubble sort of paired arrays. + +procedure ds_sort (xvals, yvals, nvals) + +real xvals[nvals] # Array of x values +real yvals[nvals] # Array of y values +int nvals # Number of values in each array + +int i, j +real temp +define swap {temp=$1;$1=$2;$2=temp} + +begin + for (i=nvals; i > 1; i=i-1) + for (j=1; j < i; j=j+1) + if (xvals[j] > xvals[j+1]) { + # Out of order; exchange y values + swap (xvals[j], xvals[j+1]) + swap (yvals[j], yvals[j+1]) + } +end diff --git a/pkg/proto/vol/src/i2sun/x_i2sun.x b/pkg/proto/vol/src/i2sun/x_i2sun.x new file mode 100644 index 00000000..20a36169 --- /dev/null +++ b/pkg/proto/vol/src/i2sun/x_i2sun.x @@ -0,0 +1,4 @@ +# X_I2SUN -- Task statement for I2SUN, used only for debugging (normally task +# resides in X_PVOL.E. + +task i2sun = t_i2sun -- cgit