diff options
Diffstat (limited to 'pkg/images/tv/display')
69 files changed, 7645 insertions, 0 deletions
diff --git a/pkg/images/tv/display/README b/pkg/images/tv/display/README new file mode 100644 index 00000000..f31a6aa4 --- /dev/null +++ b/pkg/images/tv/display/README @@ -0,0 +1,15 @@ +DISPLAY -- Prototype routines for loading and controlling the image display. +The lower level code is device dependent. + + display loads the display + dcontrol adjusts the display (frame select, window, etc.) + +The basic strategy is that the image display device is interfaced to IRAF +file i/o as a binary file. IMIO is then used to access the image or graphics +planes of the device as a disk resident imagefile would be referenced. +Each image plane of each image device is a separate "imagefile", and has a +distinct image header file in the directory "dev$". + +This package uses the ZFIOGD (binary graphics device) device driver, the +source for which is in host$gdev. It is this driver which implements physical +i/o to the device (actually, to the host system device driver for the device). diff --git a/pkg/images/tv/display/ace.h b/pkg/images/tv/display/ace.h new file mode 100755 index 00000000..4c4f40bf --- /dev/null +++ b/pkg/images/tv/display/ace.h @@ -0,0 +1,38 @@ +define NUMSTART 11 # First object number + +# Mask Flags. +define MASK_NUM 000777777B # Mask number +define MASK_GRW 001000000B # Grow pixel +define MASK_SPLIT 002000000B # Split flag +define MASK_BNDRY 004000000B # Boundary flag +define MASK_BP 010000000B # Bad pixel +define MASK_BPFLAG 020000000B # Bad pixel flag +define MASK_DARK 040000000B # Dark flag + +define MSETFLAG ori($1,$2) +define MUNSETFLAG andi($1,noti($2)) + +define MNUM (andi($1,MASK_NUM)) +define MNOTGRW (andi($1,MASK_GRW)==0) +define MGRW (andi($1,MASK_GRW)!=0) +define MNOTBP (andi($1,MASK_BP)==0) +define MBP (andi($1,MASK_BP)!=0) +define MNOTBPFLAG (andi($1,MASK_BPFLAG)==0) +define MBPFLAG (andi($1,MASK_BPFLAG)!=0) +define MNOTBNDRY (andi($1,MASK_BNDRY)==0) +define MBNDRY (andi($1,MASK_BNDRY)!=0) +define MNOTSPLIT (andi($1,MASK_SPLIT)==0) +define MSPLIT (andi($1,MASK_SPLIT)!=0) +define MNOTDARK (andi($1,MASK_DARK)==0) +define MDARK (andi($1,MASK_DARK)!=0) + +# Output object masks types. +define OM_TYPES "|boolean|numbers|colors|all|\ + |bboolean|bnumbers|bcolors|" +define OM_BOOL 1 # Boolean (0=sky, 1=object+bad+grow) +define OM_ONUM 2 # Object number only +define OM_COLORS 3 # Bad=1, Objects=2-9 +define OM_ALL 4 # All values +define OM_BBOOL 6 # Boolean (0=sky, 1=object+bad+grow) +define OM_BONUM 7 # Object number only +define OM_BCOLORS 8 # Bad=1, Objects=2-9 diff --git a/pkg/images/tv/display/display.h b/pkg/images/tv/display/display.h new file mode 100644 index 00000000..fa89a479 --- /dev/null +++ b/pkg/images/tv/display/display.h @@ -0,0 +1,42 @@ +# Display modes: + +define RGB 1 # True color mode +define FRAME 2 # Single frame mode + +# Color selections: + +define BLUE 1B # BLUE Select +define GREEN 2B # GREEN Select +define RED 4B # RED Select +define MONO 7B # RED + GREEN + BLUE + +# Size limiting parameters. + +define MAXCHAN 2 +define SAMPLE_SIZE 600 + +# If a logarithmic greyscale transformation is desired, the input range Z1:Z2 +# will be mapped into the range 1.0 to 10.0 ** MAXLOG before taking the log +# to the base 10. + +define MAXLOG 3 + +# The following parameter is used to compare display pixel coordinates for +# equality. It determines the maximum permissible magnification. The machine +# epsilon is not used because the computations are nontrivial and accumulation +# of error is a problem. + +define DS_TOL (1E-4) + +# These parameters are needed for user defined transfer functions. + +define U_MAXPTS 4096 +define U_Z1 0 +define U_Z2 4095 + +# BPDISPLAY options: + +define BPDISPLAY "|none|overlay|interpolate|" +define BPDNONE 1 # Ignore bad pixel mask +define BPDOVRLY 2 # Overlay bad pixels +define BPDINTERP 3 # Interpolate bad pixels diff --git a/pkg/images/tv/display/dsmap.x b/pkg/images/tv/display/dsmap.x new file mode 100644 index 00000000..4a5f7e9c --- /dev/null +++ b/pkg/images/tv/display/dsmap.x @@ -0,0 +1,33 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <imset.h> +include <fset.h> + +# DSMAP -- Map the display, i.e., open the display device as an imagefile. + +pointer procedure dsmap (frame, mode, color, chan) + +int frame +int mode +int color +int chan[ARB] + +pointer ds +char device[SZ_FNAME] + +int imstati(), fstati(), envgets(), imdopen() +extern imdopen() +pointer imdmap() +errchk imdmap + +begin + if (envgets ("stdimage", device, SZ_FNAME) == 0) + call error (1, "variable `stdimage' not defined in environment") + + ds = imdmap (device, mode, imdopen) + chan[1] = fstati (imstati (ds, IM_PIXFD), F_CHANNEL) + chan[2] = color + + return (ds) +end diff --git a/pkg/images/tv/display/dspmmap.x b/pkg/images/tv/display/dspmmap.x new file mode 100644 index 00000000..e20689f1 --- /dev/null +++ b/pkg/images/tv/display/dspmmap.x @@ -0,0 +1,20 @@ +# DS_PMMAP -- Open a pixel mask READ_ONLY. + +pointer procedure ds_pmmap (pmname, refim) + +char pmname[ARB] #I Pixel mask name +pointer refim #I Reference image pointer + +pointer sp, mname +pointer im, yt_mappm() +errchk yt_mappm + +begin + call smark (sp) + call salloc (mname, SZ_FNAME, TY_CHAR) + + im = yt_mappm (pmname, refim, "pmmatch", Memc[mname], SZ_FNAME) + + call sfree (sp) + return (im) +end diff --git a/pkg/images/tv/display/dsulut.x b/pkg/images/tv/display/dsulut.x new file mode 100644 index 00000000..2069bd68 --- /dev/null +++ b/pkg/images/tv/display/dsulut.x @@ -0,0 +1,141 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <error.h> +include <ctype.h> +include "display.h" + +# DS_ULUTALLOC -- 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. + +pointer procedure ds_ulutalloc (fname, z1, z2) + +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, sp, x, y +int nvalues, i, j, x1, x2, y1 +real delta_gs, delta_xv, slope +errchk ds_ulutread, ds_ulutsort, 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_ulutread (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_ulutsort (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) + return (lut) +end + + +# DS_ULUTFREE -- Free the lookup table allocated by DS_ULUT. + +procedure ds_ulutfree (lut) + +pointer lut + +begin + call mfree (lut, TY_SHORT) +end + + +# DS_ULUTREAD -- Read text file of x, y, values. + +procedure ds_ulutread (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_ULUTSORT -- Bubble sort of paired arrays. + +procedure ds_ulutsort (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/images/tv/display/findz.x b/pkg/images/tv/display/findz.x new file mode 100644 index 00000000..e1f0f73e --- /dev/null +++ b/pkg/images/tv/display/findz.x @@ -0,0 +1,62 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include "iis.h" + +# FINDZ -- Estimate the range of greylevels Z1 to Z2 containing a specified +# fraction of the greylevels in the image. The technique is to sample the +# image at some interval, computing the values of the greylevels a fixed +# distance either side of the median. Since it is not necessary to compute +# the full histogram we do not need to know the image zmin, zmax in advance. +# Works for images of any dimensionality, size, or datatype. + +procedure findz (im, z1, z2, zfrac, maxcols, nsample_lines) + +pointer im +real z1, z2, zfrac +int maxcols, nsample_lines + +real rmin, rmax +real frac +int imin, imax, ncols, nlines +int i, n, step, sample_size, imlines + +pointer sp, buf +pointer imgl2r() +include "iis.com" + +begin + call smark (sp) + call salloc (buf, ncols, TY_REAL) + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + + # Try to include a constant number of pixels in the sample + # regardless of the image size. The entire image is used if we + # have a small image, and at least sample_lines lines are read + # if we have a large image. + + sample_size = iis_ydim * nsample_lines + imlines = min(nlines, max(nsample_lines, sample_size / ncols)) + step = nlines / (imlines + 1) + + frac = (1.0 - zfrac) / 2. + imin = frac * (ncols - 1) + imax = (1.0 - frac) * (ncols - 1) + rmin = 0.0 + rmax = 0.0 + n = 0 + + do i = 1 + step, nlines, max (1, step) { + call asrtr (Memr[imgl2r (im, i)], Memr[buf], ncols) + rmin = rmin + Memr[buf + imin] + rmax = rmax + Memr[buf + imax] + n = n + 1 + } + + z1 = rmin / n + z2 = rmax / n + + call sfree (sp) +end diff --git a/pkg/images/tv/display/gwindow.h b/pkg/images/tv/display/gwindow.h new file mode 100644 index 00000000..ae91e2ea --- /dev/null +++ b/pkg/images/tv/display/gwindow.h @@ -0,0 +1,49 @@ +# Window descriptor structure. + +define LEN_WDES (210+(W_MAXWC+1)*LEN_WC) +define LEN_WC 10 # 4=[XbXeYbYe]+2=tr_type[xy] +define W_MAXWC 5 # max world coord systems +define W_SZSTRING 99 # size of strings +define W_SZIMSECT W_SZSTRING # image section string + +define W_DEVICE Memi[$1] +define W_FRAME Memi[$1+1] # device frame number +define W_XRES Memi[$1+2] # device resolution, x +define W_YRES Memi[$1+3] # device resolution, y +define W_BPDISP Memi[$1+4] # bad pixel display option +define W_BPCOLORS Memi[$1+5] # overlay colors +define W_OCOLORS Memi[$1+6] # badpixel colors +define W_IMSECT Memc[P2C($1+10)] # image section +define W_OVRLY Memc[P2C($1+60)] # overlay mask +define W_BPM Memc[P2C($1+110)] # bad pixel mask +define W_ZPM Memc[P2C($1+160)] # Z scaling pixel mask +define W_WC ($1+$2*LEN_WC+210) # ptr to coord descriptor + +# Fields of the WC coordinate descriptor, a substructure of the window +# descriptor. "W_XB(W_WC(w,0))" is the XB field of wc 0 of window W. + +define W_XS Memr[P2R($1)] # starting X value +define W_XE Memr[P2R($1+1)] # ending X value +define W_XT Memi[$1+2] # X transformation type +define W_YS Memr[P2R($1+3)] # starting Y value +define W_YE Memr[P2R($1+4)] # ending Y value +define W_YT Memi[$1+5] # Y transformation type +define W_ZS Memr[P2R($1+6)] # starting Z value (greyscale) +define W_ZE Memr[P2R($1+7)] # ending Z value +define W_ZT Memi[$1+8] # Z transformation type +define W_UPTR Memi[$1+9] # LUT when ZT=USER + +# WC types. + +define W_NWIN 0 # Display window in NDC coordinates +define W_DWIN 1 # Display window in image pixel coordinates +define W_WWIN 2 # Display window in image world coordinates +define W_IPIX 3 # Image pixel coordinates (in pixels) +define W_DPIX 4 # Display pixel coordinates (in pixels) + +# Types of coordinate and greyscale transformations. + +define W_UNITARY 0 # values map without change +define W_LINEAR 1 # linear mapping +define W_LOG 2 # logarithmic mapping +define W_USER 3 # user specifies transformation diff --git a/pkg/images/tv/display/iis.com b/pkg/images/tv/display/iis.com new file mode 100644 index 00000000..8b367132 --- /dev/null +++ b/pkg/images/tv/display/iis.com @@ -0,0 +1,25 @@ +# Common for IIS display + +int iischan # the device channel used by FIO +int iisnopen # number of times the display has been opened +int iisframe # frame number at iisopn time (kludge). +int iis_xdim, iis_ydim # frame size, pixels +int iis_config # frame size configuration +int iis_server # device is actually a display server +bool packit # byte pack data for i/o +bool swap_bytes # byte swap the IIS header +short hdr[LEN_IISHDR] # header + +int iis_version # WCS version +int iis_valid # valid mapping info flag +char iis_region[SZ_FNAME] # region name +real iis_sx, iis_sy # source raster offset +int iis_snx, iis_sny # source raster size +int iis_dx, iis_dy # dest raster offset +int iis_dnx, iis_dny # dest raster size +char iis_objref[SZ_FNAME] # object reference + +common /iiscom/ iischan, iisnopen, iisframe, iis_xdim, iis_ydim, iis_config, + iis_server, packit, swap_bytes, hdr, iis_version, iis_valid, + iis_region, iis_sx, iis_sy, iis_snx, iis_sny, + iis_dx, iis_dy, iis_dnx, iis_dny, iis_objref diff --git a/pkg/images/tv/display/iis.h b/pkg/images/tv/display/iis.h new file mode 100644 index 00000000..bdd4f33a --- /dev/null +++ b/pkg/images/tv/display/iis.h @@ -0,0 +1,121 @@ +# This file contains the hardware definitions for the iis model 70/f +# at Kitt Peak. + +# Define header +define LEN_IISHDR 8 # Length of IIS header + +define XFERID $1[1] # transfer id +define THINGCT $1[2] # thing count +define SUBUNIT $1[3] # subuint select +define CHECKSUM $1[4] # check sum +define XREG $1[5] # x register +define YREG $1[6] # y register +define ZREG $1[7] # z register +define TREG $1[8] # t register + + +# Transfer ID definitions +define IREAD 100000B +define IWRITE 0B +define PACKED 40000B +define SAMPLE 40000B +define BYPASSIFM 20000B +define BYTE 10000B +define ADDWRITE 4000B +define ACCUM 2000B +define BLOCKXFER 1000B +define VRETRACE 400B +define MUX32 200B +define IMT800 100B # [IMTOOL SPECIAL] + +# Subunits +define REFRESH 1 +define LUT 2 +define OFM 3 +define IFM 4 +define FEEDBACK 5 +define SCROLL 6 +define VIDEOM 7 +define SUMPROC 8 +define GRAPHICS 9 +define CURSOR 10 +define ALU 11 +define ZOOM 12 +define IMCURSOR 20B +define WCS 21B + +# Command definitions +define COMMAND 100000B +define ADVXONTC 100000B # Advance x on thing count +define ADVXONYOV 40000B # Advance x on y overflow +define ADVYONXOV 100000B # Advance y on x overflow +define ADVYONTC 40000B # Advance y on thing count +define ERASE 100000B # Erase + +# 4 - Button Trackball +define PUSH 40000B +define BUTTONA 400B +define BUTTONB 1000B +define BUTTONC 2000B +define BUTTOND 4000B + +# Display channels +define CHAN1 1B +define CHAN2 2B +define CHAN3 4B +define CHAN4 10B +define CHAN5 20B +define CHAN6 40B +define CHAN7 100B +define CHAN8 200B +define CHAN9 400B +define CHAN10 1000B +define CHAN11 2000B +define CHAN12 4000B +define CHAN13 10000B +define CHAN14 20000B +define CHAN15 40000B +define CHAN16 100000B +define GRCHAN 100000B + +define LEN_IISFRAMES 16 +define IISFRAMES CHAN1, CHAN2, CHAN3, CHAN4, CHAN5, CHAN6, CHAN7, CHAN8, CHAN9, CHAN10, CHAN11, CHAN12, CHAN13, CHAN14, CHAN15, CHAN16 + +# Colors + +define BLUE 1B +define GREEN 2B +define RED 4B +define MONO 7B + +# Bit plane selections +define BITPL0 1B +define BITPL1 2B +define BITPL2 4B +define BITPL3 10B +define BITPL4 20B +define BITPL5 40B +define BITPL6 100B +define BITPL7 200B +define ALLBITPL 377B + +# IIS Sizes +define IIS_XDIM 512 +define IIS_YDIM 512 +define MCXSCALE 64 # metacode x scale +define MCYSCALE 64 # metacode y scale +define SZB_IISHDR 16 # size of IIS header in bytes +define SZB_IMCURVAL 160 # size of imcursor value buffer, bytes +define LEN_ZOOM 3 # zoom parameters +define LEN_CURSOR 3 # cursor parameters +define LEN_SPLIT 12 # split screen +define LEN_LUT 256 # look up table +define LEN_OFM 1024 # output function look up table +define SZ_OLD_WCSTEXT 320 # old max WCS text chars +define SZ_WCSTEXT 1024 # max WCS text chars + +# IIS Status Words +define IIS_FILSIZE (IIS_XDIM * IIS_YDIM * SZB_CHAR) +define IIS_BLKSIZE 1024 +define IIS_OPTBUFSIZE (IIS_XDIM * SZB_CHAR) +define IIS_MAXBUFSIZE 32768 diff --git a/pkg/images/tv/display/iisblk.x b/pkg/images/tv/display/iisblk.x new file mode 100644 index 00000000..1ff81d49 --- /dev/null +++ b/pkg/images/tv/display/iisblk.x @@ -0,0 +1,40 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISBLK -- Blink IIS display frames at millisecond time resolution. + +procedure iisblk (chan1, chan2, chan3, chan4, nframes, rate) + +int chan1[ARB] +int chan2[ARB] +int chan3[ARB] +int chan4[ARB] +int nframes +real rate + +int msec, status, xcur, ycur +int and() + +begin + status = 0 + msec = int (rate * 1000.) + + while (and (status, PUSH) == 0) { + call zwmsec (msec) + call iisrgb (chan1, chan1, chan1) + call zwmsec (msec) + call iisrgb (chan2, chan2, chan2) + if (nframes >= 3) { + call zwmsec (msec) + call iisrgb (chan3, chan3, chan3) + } + if (nframes == 4) { + call zwmsec (msec) + call iisrgb (chan4, chan4, chan4) + } + call iisrcr (status, xcur, ycur) + } +end diff --git a/pkg/images/tv/display/iiscls.x b/pkg/images/tv/display/iiscls.x new file mode 100644 index 00000000..71da6c35 --- /dev/null +++ b/pkg/images/tv/display/iiscls.x @@ -0,0 +1,24 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <knet.h> +include "zdisplay.h" +include "iis.h" + +# IISCLS -- Close IIS display. + +procedure iiscls (chan, status) + +int chan[ARB] +int status +include "iis.com" + +begin + if (iisnopen == 1) { + call zclsgd (iischan, status) + iisnopen = 0 + } else if (iisnopen > 1) { + iisnopen = iisnopen - 1 + } else + iisnopen = 0 +end diff --git a/pkg/images/tv/display/iisers.x b/pkg/images/tv/display/iisers.x new file mode 100644 index 00000000..de276a99 --- /dev/null +++ b/pkg/images/tv/display/iisers.x @@ -0,0 +1,28 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISERS -- Erase IIS frame. + +procedure iisers (chan) + +int chan[ARB] +short erase + +int status, tid +int iisflu(), andi() +include "iis.com" + +begin + call achtiu (andi (ERASE, 0177777B), erase, 1) + + # IMTOOL special - IIS frame bufrer configuration code. + tid = IWRITE+BYPASSIFM+BLOCKXFER + tid = tid + max (0, iis_config - 1) + + call iishdr (tid, 1, FEEDBACK, ADVXONTC, ADVYONXOV, iisflu(chan), + ALLBITPL) + call iisio (erase, SZB_CHAR, status) +end diff --git a/pkg/images/tv/display/iisflu.x b/pkg/images/tv/display/iisflu.x new file mode 100644 index 00000000..3fee9d63 --- /dev/null +++ b/pkg/images/tv/display/iisflu.x @@ -0,0 +1,24 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISFLU -- IIS frame look up table. + +int procedure iisflu (chan) + +int chan[ARB] +int frame +int iisframe[LEN_IISFRAMES] +data iisframe/IISFRAMES/ + +begin + frame = chan[1] - IIS_CHAN * DEVCODE + if (frame < 1) + return (iisframe[1]) + else if (frame > LEN_IISFRAMES) + return (GRCHAN) + else + return (iisframe[frame]) +end diff --git a/pkg/images/tv/display/iisgop.x b/pkg/images/tv/display/iisgop.x new file mode 100644 index 00000000..c33f21d2 --- /dev/null +++ b/pkg/images/tv/display/iisgop.x @@ -0,0 +1,14 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "iis.h" + +# IISGOP -- Open IIS graphics display. + +procedure iisgop (frame, mode, chan) + +int frame, mode, chan[ARB] + +begin + call iisopn (frame + LEN_IISFRAMES, mode, chan) +end diff --git a/pkg/images/tv/display/iishdr.x b/pkg/images/tv/display/iishdr.x new file mode 100644 index 00000000..38ea733d --- /dev/null +++ b/pkg/images/tv/display/iishdr.x @@ -0,0 +1,30 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" +# IISHDR -- Form IIS header. + +procedure iishdr (id, count, subunit, x, y, z, t) + +int id, count, subunit, x, y, z, t +int i, sum +include "iis.com" + +begin + call achtiu (id, XFERID(hdr), 1) + call achtiu (count, THINGCT(hdr), 1) + call achtiu (subunit, SUBUNIT(hdr), 1) + call achtiu (x, XREG(hdr), 1) + call achtiu (y, YREG(hdr), 1) + call achtiu (z, ZREG(hdr), 1) + call achtiu (t, TREG(hdr), 1) + CHECKSUM(hdr) = 1 + + if (THINGCT(hdr) > 0) + THINGCT(hdr) = -THINGCT(hdr) + sum = 0 + for (i = 1; i <= LEN_IISHDR; i = i + 1) + sum = sum + hdr[i] + call achtiu (-sum, CHECKSUM(hdr), 1) +end diff --git a/pkg/images/tv/display/iisio.x b/pkg/images/tv/display/iisio.x new file mode 100644 index 00000000..ad3902ed --- /dev/null +++ b/pkg/images/tv/display/iisio.x @@ -0,0 +1,43 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <knet.h> +include "zdisplay.h" +include "iis.h" + +# IISIO -- Synchronous i/o to the IIS. + +procedure iisio (buf, nbytes, status) + +short buf[ARB] +int nbytes +int status + +int xferid +int and() +include "iis.com" + +begin + call iiswt (iischan, status) + xferid = XFERID(hdr) + + if (swap_bytes) + call bswap2 (hdr, 1, hdr, 1, SZB_IISHDR) + call zawrgd (iischan, hdr, SZB_IISHDR, 0) + call iiswt (iischan, status) + + if (and (xferid, IREAD) != 0) { + call zardgd (iischan, buf, nbytes, 0) + call iiswt (iischan, status) + if (swap_bytes && and(xferid,PACKED) == 0) + call bswap2 (buf, 1, buf, 1, nbytes) + } else { + if (swap_bytes && and(xferid,PACKED) == 0) + call bswap2 (buf, 1, buf, 1, nbytes) + call zawrgd (iischan, buf, nbytes, 0) + call iiswt (iischan, status) + } + + if (status <= 0) + status = EOF +end diff --git a/pkg/images/tv/display/iismtc.x b/pkg/images/tv/display/iismtc.x new file mode 100644 index 00000000..2d6eb2cf --- /dev/null +++ b/pkg/images/tv/display/iismtc.x @@ -0,0 +1,21 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISMTC -- Match channel lut to frame2. + +procedure iismtc (chan1, chan2) + +int chan1[ARB], chan2[ARB] +short lut[LEN_LUT] + +int iisflu() + +begin + if (iisflu (chan2) == GRCHAN) + return + call iisrlt (chan1, lut) + call iiswlt (chan2, lut) +end diff --git a/pkg/images/tv/display/iisofm.x b/pkg/images/tv/display/iisofm.x new file mode 100644 index 00000000..24259fd3 --- /dev/null +++ b/pkg/images/tv/display/iisofm.x @@ -0,0 +1,183 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <math.h> +include "zdisplay.h" +include "iis.h" + +# These procedures have been modified to limit the maximum output level. + +define NIN 256 # Number of input levels +define NOUT 1024 # Number of output levels + +# IISOFM -- Output color mapping. + +procedure iisofm (map) + +char map[ARB] # type of mapping + +int i +short lutr[LEN_OFM] +short lutg[LEN_OFM] +short lutb[LEN_OFM] + +begin + if (map[1] == 'm') { # MONO + do i = 1, LEN_OFM + lutr[i] = min ((i - 1) * NOUT / NIN, NOUT) + call iiswom (MONO, lutr) + return + } + + call aclrs (lutr, LEN_OFM) + call aclrs (lutg, LEN_OFM) + call aclrs (lutb, LEN_OFM) + + if (map[1] == 'l') { # LINEAR + call iislps (lutb, lutg, lutr) + + } else if (map[1] == '8') { # 8COLOR + do i = 33, 64 { + lutb[i] = NOUT - 1 + lutr[i] = NOUT - 1 + } + do i = 65, 96 + lutb[i] = NOUT - 1 + do i = 97, 128 { + lutb[i] = NOUT - 1 + lutg[i] = NOUT - 1 + } + do i = 129, 160 + lutg[i] = NOUT - 1 + do i = 161, 192 { + lutg[i] = NOUT - 1 + lutr[i] = NOUT - 1 + } + do i = 193, 224 + lutr[i] = NOUT - 1 + do i = 225, 256 { + lutr[i] = NOUT - 1 + lutg[i] = NOUT - 1 + lutb[i] = NOUT - 1 + } + do i = 257, LEN_OFM { + lutr[i] = NOUT - 1 + lutg[i] = NOUT - 1 + lutb[i] = NOUT - 1 + } + + } else if (map[1] == 'r') { # RANDOM + do i = 2, LEN_OFM, 8 { + lutr[i] = NOUT - 1 + lutb[i] = NOUT - 1 + } + do i = 3, LEN_OFM, 8 + lutb[i] = NOUT - 1 + do i = 4, LEN_OFM, 8 { + lutb[i] = NOUT - 1 + lutg[i] = NOUT - 1 + } + do i = 5, LEN_OFM, 8 + lutg[i] = NOUT - 1 + do i = 6, LEN_OFM, 8 { + lutg[i] = NOUT - 1 + lutr[i] = NOUT - 1 + } + do i = 7, LEN_OFM, 8 + lutr[i] = NOUT - 1 + do i = 8, LEN_OFM, 8 { + lutr[i] = NOUT - 1 + lutg[i] = NOUT - 1 + lutb[i] = NOUT - 1 + } + } + + call iiswom (RED, lutr) + call iiswom (GREEN, lutg) + call iiswom (BLUE, lutb) +end + + +# IISWOM -- Write output color look up table. + +procedure iiswom (color, lut) + +int color +short lut[ARB] +int status + +begin + call iishdr (IWRITE+VRETRACE, LEN_OFM, OFM, ADVXONTC, ADVYONXOV, + color, 0) + call iisio (lut, LEN_OFM * SZB_CHAR, status) +end + + +# IISROM -- Read color look up table. + +procedure iisrom (color, lut) + +int color +short lut[ARB] +int status + +begin + call iishdr (IREAD+VRETRACE, LEN_OFM, LUT, ADVXONTC, ADVYONXOV, + color, 0) + call iisio (lut, LEN_OFM * SZB_CHAR, status) +end + + +# Linear Pseudocolor Modelling code. + +define BCEN 64 +define GCEN 128 +define RCEN 196 + +# IISLPS -- Load the RGB luts for linear pseudocolor. + +procedure iislps (lutb, lutg, lutr) + +short lutb[ARB] # blue lut +short lutg[ARB] # green lut +short lutr[ARB] # red lut + +begin + # Set the mappings for the primary color bands. + call iislps_curve (lutb, NIN, BCEN, NOUT - 1, NIN/2) + call iislps_curve (lutg, NIN, GCEN, NOUT - 1, NIN/2) + call iislps_curve (lutr, NIN, RCEN, NOUT - 1, NIN/2) + + # Add one half band of white color at the right. + call iislps_curve (lutb, NIN, NIN, NOUT - 1, NIN/2) + call iislps_curve (lutg, NIN, NIN, NOUT - 1, NIN/2) + call iislps_curve (lutr, NIN, NIN, NOUT - 1, NIN/2) +end + + +# IISLPS_CURVE -- Compute the lookup table for a single color. + +procedure iislps_curve (y, npts, xc, height, width) + +short y[npts] # output curve +int npts # number of points +int xc # x center +int height, width + +int i +real dx, dy, hw + +begin + hw = width / 2.0 + dy = height / hw * 2.0 + + do i = 1, npts { + dx = abs (i - xc) + if (dx > hw) + ; + else if (dx > hw / 2.0) + y[i] = max (int(y[i]), min (height, int((hw - dx) * dy))) + else + y[i] = height + } +end diff --git a/pkg/images/tv/display/iisopn.x b/pkg/images/tv/display/iisopn.x new file mode 100644 index 00000000..a310e168 --- /dev/null +++ b/pkg/images/tv/display/iisopn.x @@ -0,0 +1,76 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <knet.h> +include "zdisplay.h" +include "iis.h" + +# ---------------------------------------------------------------------- +# MODIFIED VERSION OF IISOPN.X FOR IMTOOL -- DO NOT DELETE. +# Referenced by the Sun/IRAF special file list: see hlib$mkpkg.sf. +# ---------------------------------------------------------------------- + +# IISOPN -- Open IIS display. + +procedure iisopn (devinfo, mode, chan) + +char devinfo[ARB] # device info for zopen (packed) +int mode # access mode +int chan[ARB] # receives IIS descriptor + +int delim +char resource[SZ_FNAME] +char node[SZ_FNAME] +bool first_time +data first_time /true/ +int ki_gnode(), strncmp() +include "iis.com" +include "imd.com" +define quit_ 91 + +begin + if (first_time) { + iisnopen = 0 + iis_version = 0 + first_time = false + } + + # We permit multiple opens but only open the physical device once. + if (iisnopen == 0) { + call zopngd (devinfo, mode, iischan) + + # Initialize imd_gcur. + call strcpy (devinfo, imd_devinfo, SZ_LINE) + imd_mode = mode + imd_magic = -1 + } + + if (iischan != ERR) { + iisnopen = iisnopen + 1 + chan[1] = FRTOCHAN(iisframe) + + # The following code is DEVICE DEPENDENT (horrible kludge, but + # it simplifies things and this is throw away code). + + # Byte pack i/o if the device is on a remote node since the i/o + # bandwidth is the limiting factor; do not bytepack if on local + # node since cpu time is the limiting factor. + + call strupk (devinfo, resource, SZ_FNAME) + packit = (ki_gnode (resource, node, delim) != 0) + if (!packit) + packit = (strncmp (resource[delim+1], "imt", 3) == 0) + + # Enable byte swapping if the device is byte swapped but the + # local host is not (assumes that if there is an IIS it is on + # a byte swapped VAX - this should be done in graphcap instead). + + swap_bytes = (strncmp (resource[delim+1], "iis", 3) == 0 && + BYTE_SWAP2 == NO) + + # Initialize zoom. + call iiszm (1, 0, 0) + + } else + chan[1] = ERR +end diff --git a/pkg/images/tv/display/iispio.x b/pkg/images/tv/display/iispio.x new file mode 100644 index 00000000..81e2512d --- /dev/null +++ b/pkg/images/tv/display/iispio.x @@ -0,0 +1,97 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <knet.h> +include "zdisplay.h" +include "iis.h" + +# IISPIO -- Asynchronous pixel i/o to the IIS. + +procedure iispio (buf, nx, ny) + +short buf[nx,ny] # Cell array +int nx, ny # length, number of image lines + +pointer iobuf +bool first_time +int xferid, status, nbytes, szline, i +int and() +include "iis.com" +data first_time /true/ + +begin + if (first_time) { + if (packit) + i = IIS_MAXBUFSIZE + else + i = IIS_MAXBUFSIZE * (SZ_SHORT * SZB_CHAR) + call malloc (iobuf, i, TY_SHORT) + first_time = false + } + + # Wait for the last i/o transfer. + call iiswt (iischan, status) + if (status == ERR) + return + + # Disable interrupts while transmitting to or receiving data from + # the display, to avoid loss of synch on the datastream and resulting + # loss of communications with the device. + + call intr_disable() + xferid = XFERID(hdr) + + # Transmit the packet header. + if (swap_bytes) + call bswap2 (hdr, 1, hdr, 1, SZB_IISHDR) + call zawrgd (iischan, hdr, SZB_IISHDR, 0) + call iiswt (iischan, status) + if (status == ERR) { + call intr_enable() + return + } + + # Read or write the data block. + nbytes = ny * iis_xdim + szline = iis_xdim + + if (packit) + szline = szline / (SZ_SHORT * SZB_CHAR) + else + nbytes = nbytes * (SZ_SHORT * SZB_CHAR) + + # Transmit the data byte-packed to increase the i/o bandwith + # when using network i/o. + + if (and (xferid, IREAD) != 0) { + # Read from the IIS. + + call zardgd (iischan, Mems[iobuf], nbytes, 0) + call iiswt (iischan, status) + + # Unpack and line flip the packed data. + if (packit) { + do i = 0, ny-1 + call achtbs (Mems[iobuf+i*szline], buf[1,ny-i], iis_xdim) + } else { + do i = 0, ny-1 + call amovs (Mems[iobuf+i*szline], buf[1,ny-i], szline) + } + + } else { + # Write to the IIS. + + # Bytepack the image lines, doing a line flip in the process. + if (packit) { + do i = 0, ny-1 + call achtsb (buf[1,ny-i], Mems[iobuf+i*szline], iis_xdim) + } else { + do i = 0, ny-1 + call amovs (buf[1,ny-i], Mems[iobuf+i*szline], szline) + } + + call zawrgd (iischan, Mems[iobuf], nbytes, 0) + } + + call intr_enable() +end diff --git a/pkg/images/tv/display/iisrcr.x b/pkg/images/tv/display/iisrcr.x new file mode 100644 index 00000000..53119d06 --- /dev/null +++ b/pkg/images/tv/display/iisrcr.x @@ -0,0 +1,32 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +define DELAY 30 # milliseconds between cursor reads + + +# IISRCR -- Read cursor from display. Note that the position is 1 indexed. + +procedure iisrcr (status, xcur, ycur) + +int status, xcur, ycur +short cursor[LEN_CURSOR] +include "iis.com" + +begin + call iishdr(IREAD+VRETRACE, LEN_CURSOR, COMMAND+CURSOR, ADVXONTC, 0,0,0) + + call zwmsec (DELAY) + + call iisio (cursor, LEN_CURSOR * SZB_CHAR, status) + if (status <= 0) { + status = EOF + return + } + + status = cursor[1] + xcur = MCXSCALE * mod (cursor[2] + 31, iis_xdim) + ycur = MCYSCALE * mod (cursor[3] + 31, iis_ydim) +end diff --git a/pkg/images/tv/display/iisrd.x b/pkg/images/tv/display/iisrd.x new file mode 100644 index 00000000..3421a71f --- /dev/null +++ b/pkg/images/tv/display/iisrd.x @@ -0,0 +1,42 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISRD -- Read data from IIS. + +procedure iisrd (chan, buf, nbytes, offset) + +int chan[ARB] +short buf[ARB] +int nbytes +long offset + +long off1, off2 +int nchars, thing_count, tid, y1, y2, x +int or(), iisflu() +include "iis.com" + +begin + # Convert to chars and clip at the top of the display. + off1 = (offset - 1) / SZB_CHAR + 1 + off2 = min (iis_xdim * iis_ydim, (offset + nbytes - 1) / SZB_CHAR) + 1 + nchars = off2 - off1 + + x = 0 + y1 = (off1-1 ) / iis_xdim + y2 = (off2-1 - iis_xdim) / iis_xdim + y2 = max (y1, y2) + + if (packit) + tid = IREAD+PACKED + else + tid = IREAD + thing_count = nchars + + call iishdr (tid, thing_count, REFRESH, or(x,ADVXONTC), + or(iis_ydim-y2-1, ADVYONXOV), iisflu(chan), ALLBITPL) + + call iispio (buf, iis_xdim, y2 - y1 + 1) +end diff --git a/pkg/images/tv/display/iisrgb.x b/pkg/images/tv/display/iisrgb.x new file mode 100644 index 00000000..9dcc38cd --- /dev/null +++ b/pkg/images/tv/display/iisrgb.x @@ -0,0 +1,32 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISRGB -- Enable RGB display. + +procedure iisrgb (red_chan, green_chan, blue_chan) + +int red_chan[ARB], green_chan[ARB], blue_chan[ARB] + +int i, frm, status +short split[LEN_SPLIT] +int iisflu() + +begin + frm = iisflu (blue_chan) + do i = 1, 4 + split[i] = frm + + frm = iisflu (green_chan) + do i = 5, 8 + split[i] = frm + + frm = iisflu (red_chan) + do i = 9, 12 + split[i] = frm + + call iishdr (IWRITE+VRETRACE, LEN_SPLIT, COMMAND+LUT, ADVXONTC, 0, 0, 0) + call iisio (split, LEN_SPLIT * SZB_CHAR, status) +end diff --git a/pkg/images/tv/display/iissfr.x b/pkg/images/tv/display/iissfr.x new file mode 100644 index 00000000..f6e92013 --- /dev/null +++ b/pkg/images/tv/display/iissfr.x @@ -0,0 +1,15 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "iis.h" + +# IIS_SETFRAME -- Set the frame number for IISOPN. This is a kludge to pass +# this number to IISOPN via the iis common. + +procedure iis_setframe (frame) + +int frame +include "iis.com" + +begin + iisframe = frame +end diff --git a/pkg/images/tv/display/iisstt.x b/pkg/images/tv/display/iisstt.x new file mode 100644 index 00000000..86474d25 --- /dev/null +++ b/pkg/images/tv/display/iisstt.x @@ -0,0 +1,29 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <fio.h> +include "zdisplay.h" +include "iis.h" + +# IISSTT -- IIS status. +# [OBSOLETE - NO LONGER USED (see zsttim)] + +procedure iisstt (chan, what, lvalue) + +int chan[ARB], what +long lvalue + +begin + switch (what) { + case FSTT_FILSIZE: + lvalue = IIS_FILSIZE + case FSTT_BLKSIZE: + lvalue = IIS_BLKSIZE + case FSTT_OPTBUFSIZE: + lvalue = IIS_OPTBUFSIZE + case FSTT_MAXBUFSIZE: + lvalue = IIS_MAXBUFSIZE + default: + lvalue = ERR + } +end diff --git a/pkg/images/tv/display/iiswcr.x b/pkg/images/tv/display/iiswcr.x new file mode 100644 index 00000000..3970f230 --- /dev/null +++ b/pkg/images/tv/display/iiswcr.x @@ -0,0 +1,20 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISWCR -- Write cursor to display. Note that the position is 1 indexed. + +procedure iiswcr (status, xcur, ycur) + +int status, xcur, ycur +short cursor[LEN_CURSOR] +include "iis.com" + +begin + call iishdr (IWRITE+VRETRACE, 2, COMMAND+CURSOR, 1+ADVXONTC, 0,0,0) + cursor[2] = mod (xcur / MCXSCALE - 32, iis_xdim) + cursor[3] = mod (ycur / MCYSCALE - 32, iis_ydim) + call iisio (cursor[2], 2 * SZB_CHAR, status) +end diff --git a/pkg/images/tv/display/iiswnd.x b/pkg/images/tv/display/iiswnd.x new file mode 100644 index 00000000..e906cc1f --- /dev/null +++ b/pkg/images/tv/display/iiswnd.x @@ -0,0 +1,117 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISWND -- Window IIS display frame with the trackball. + +procedure iiswnd3 (chan1, chan2, chan3) + +int chan1[ARB], chan2[ARB], chan3[ARB] + +int i, j +real x, y +short lut[LEN_LUT] +int status, xcur, ycur, lutval +int iisflu(), and() + +begin + if (iisflu(chan1) == GRCHAN) + return + call iisrlt (chan1, lut) + + # Starting point at lut[2] because lut[1] is background + for (i=3; (i < 257) && (lut[i] == lut[2]); i=i+1) + ; + i = i - 1 + + for (j=255; (j > i) && (lut[j] == lut[256]); j=j-1) + ; + j = j + 1 + + if ((i == j) || (lut[i] == lut[j])) { + xcur = 256 + ycur = 384 + } else { + y = real (lut[j] - lut[i]) / (j - i) + xcur = 2 * (i - 1) - (2 * lut[i] - 256) / y + 1 + if (y > 1) + y = 2 - (1 / y) + if (y < -1) + y = -2 - (1 / y) + ycur = 128 * y + 256.5 + } + + xcur = xcur * MCXSCALE + ycur = ycur * MCYSCALE + call iiswcr (status, xcur, ycur) + status = 0 + + while (and (status, PUSH) == 0) { + call iisrcr (status, xcur, ycur) + if (status == EOF) + break + + xcur = xcur / MCXSCALE + ycur = ycur / MCYSCALE + x = xcur / 2 + y = (ycur - 255.5) / 128. + + if (y > 1) + y = 1. / (2 - y) + if (y < - 1) + y = -1. / (2 + y) + do i = 1, 256 { + lutval = y * (i - 1 - x) + 127.5 + lut[i] = max (0, min (255, lutval)) + } + + lut[1] = 0 # Make background black + if ((chan1[1] == chan2[1]) && (chan1[1] == chan3[1])) + call iiswlt (chan1, lut) + else { + call iiswlt (chan1, lut) + call iiswlt (chan2, lut) + call iiswlt (chan3, lut) + } + } +end + + +# IISWLT -- Write monochrome look up table. + +procedure iiswlt (chan, lut) + +int chan[ARB] +short lut[ARB] + +int status +int iisflu() + +begin + if (iisflu (chan) == GRCHAN) + return + call iishdr (IWRITE+VRETRACE, LEN_LUT, LUT, ADVXONTC, 0, chan[2], + iisflu (chan)) + call iisio (lut, LEN_LUT * SZB_CHAR, status) +end + + +# IISRLT -- Read monochrome look up table. + +procedure iisrlt (chan, lut) + +int chan[ARB] +short lut[ARB] + +int status +int iisflu() + +begin + if (iisflu (chan) == GRCHAN) + return + call iishdr (IREAD+VRETRACE, LEN_LUT, LUT, ADVXONTC, 0, 0, + iisflu (chan)) + call iisio (lut, LEN_LUT * SZB_CHAR, status) +end diff --git a/pkg/images/tv/display/iiswr.x b/pkg/images/tv/display/iiswr.x new file mode 100644 index 00000000..68a1a583 --- /dev/null +++ b/pkg/images/tv/display/iiswr.x @@ -0,0 +1,48 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISWR -- Write pixel data to IIS. Writes are limited to entire display lines. +# The data is line-flipped, causing the first line to be displayed at the bottom +# of the screen. + +procedure iiswr (chan, buf, nbytes, offset) + +int chan[ARB] # io channel +short buf[ARB] # pixels +int nbytes # length of pixel array in bytes +long offset # pixel offset in image display + +long off1, off2 +int nchars, thing_count, tid, y1, y2, x +int or(), iisflu() +include "iis.com" + +begin + # Convert to chars and clip at the top of the display. + off1 = (offset - 1) / SZB_CHAR + 1 + off2 = min (iis_xdim * iis_ydim, (offset + nbytes - 1) / SZB_CHAR) + 1 + nchars = off2 - off1 + + x = 0 + y1 = (off1-1 ) / iis_xdim + y2 = (off2-1 - iis_xdim) / iis_xdim + y2 = max (y1, y2) + +#call eprintf ("iiswr: %d bytes at %d, x=%d, y=[%d:%d]\n") +#call pargi(nbytes); call pargi(offset) +#call pargi(x); call pargi(y1); call pargi(y2) + + if (packit) + tid = IWRITE+BYPASSIFM+BLOCKXFER+BYTE+PACKED + else + tid = IWRITE+BYPASSIFM + thing_count = nchars + + call iishdr (tid, thing_count, REFRESH, or(x,ADVXONTC), + or(iis_ydim-y2-1, ADVYONXOV), iisflu(chan), ALLBITPL) + + call iispio (buf, iis_xdim, y2 - y1 + 1) +end diff --git a/pkg/images/tv/display/iiswt.x b/pkg/images/tv/display/iiswt.x new file mode 100644 index 00000000..ae18ebff --- /dev/null +++ b/pkg/images/tv/display/iiswt.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <knet.h> +include "zdisplay.h" +include "iis.h" + +# IISWT -- Wait for IIS display. + +procedure iiswt (chan, nbytes) + +int chan[ARB], nbytes +include "iis.com" + +begin + call zawtgd (iischan, nbytes) + if (packit) + nbytes = nbytes * (SZ_SHORT * SZB_CHAR) +end diff --git a/pkg/images/tv/display/iiszm.x b/pkg/images/tv/display/iiszm.x new file mode 100644 index 00000000..d207f47a --- /dev/null +++ b/pkg/images/tv/display/iiszm.x @@ -0,0 +1,38 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include "zdisplay.h" +include "iis.h" + +# IISZM -- Zoom IIS window. + +procedure iiszm (zfactor, x, y) + +int zfactor, x, y +short zoom[LEN_ZOOM] +int status + +begin + call iishdr (IWRITE+VRETRACE, LEN_ZOOM, ZOOM, ADVXONTC, 0, 0, 0) + zoom[1] = zfactor - 1 + zoom[2] = x / MCXSCALE + zoom[3] = y / MCYSCALE + call iisio (zoom, LEN_ZOOM * SZB_CHAR, status) +end + + +# IISRM -- Roam IIS display. + +procedure iisrm (zfactor) + +int zfactor +int status, xcur, ycur +int and() + +begin + status = 0 + while (status != EOF && and (status, PUSH) == 0) { + call iisrcr (status, xcur, ycur) + call iiszm (zfactor, xcur, ycur) + } +end diff --git a/pkg/images/tv/display/imd.com b/pkg/images/tv/display/imd.com new file mode 100644 index 00000000..9738e89b --- /dev/null +++ b/pkg/images/tv/display/imd.com @@ -0,0 +1,7 @@ +# IMD.COM -- Common for the IMD routines. + +int imd_magic # set to -1 when initialized +int imd_mode # display access mode +char imd_devinfo[SZ_LINE] # device information for zopngd + +common /imdcom/ imd_magic, imd_mode, imd_devinfo diff --git a/pkg/images/tv/display/imdgcur.x b/pkg/images/tv/display/imdgcur.x new file mode 100644 index 00000000..0f8cf658 --- /dev/null +++ b/pkg/images/tv/display/imdgcur.x @@ -0,0 +1,37 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <knet.h> +include "iis.h" + +# IMD_GCUR -- This is functionally equivalent to CLGCUR and should be used in +# place of the latter routine in programs which directly map the display. +# Its function is to close off the display at a low level in order to free +# the display device for access by the CL process for the cursor read. + +int procedure imd_gcur (param, wx, wy, wcs, key, strval, maxch) + +char param[ARB] # parameter to be read [not used] +real wx, wy # cursor coordinates +int wcs # wcs to which coordinates belong +int key # keystroke value of cursor event +char strval[ARB] # string value, if any +int maxch + +int status +bool devopen +int clgcur() +include "iis.com" +include "imd.com" + +begin + devopen = (iisnopen > 0) + if (imd_magic == -1 && devopen) + call zclsgd (iischan, status) + + status = clgcur (param, wx, wy, wcs, key, strval, maxch) + + if (imd_magic == -1 && devopen) + call zopngd (imd_devinfo, imd_mode, iischan) + + return (status) +end diff --git a/pkg/images/tv/display/imdgetwcs.x b/pkg/images/tv/display/imdgetwcs.x new file mode 100644 index 00000000..57f432bc --- /dev/null +++ b/pkg/images/tv/display/imdgetwcs.x @@ -0,0 +1,188 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <ctype.h> +include "zdisplay.h" +include "iis.h" + +# IMD_GETWCS -- Get the saved WCS for the given frame of the given display +# device. (No great attempt at generality here). +# [INTERNAL ROUTINE - RESTRICTED USE]. +# +# Example: +# +# dev$pix - m51 B 600s +# 1. 0. 0. -1. 1. 512. 36. 320.0713 1 +# +# The file format is the image title, followed by a line specifying the +# coordinate transformation matrix (6 numbers: a b c d tx ty) and the +# greyscale transformation (z1 z2 zt). +# +# The procedure returns OK if the WCS for the frame is sucessfully accessed, +# or ERR if the WCS cannot be read. In the latter case the output WCS will +# be the default unitary WCS. + +int procedure imd_getwcs (frame, server, image, sz_image, title, sz_title, + a, b, c, d, tx, ty) + +int frame #I frame (wcs) number of current device +int server #I device is a display server +char image[ARB] #O image name +int sz_image #I max image name length +char title[ARB] #O image title string +int sz_title #I max image title length +real a, d #O x, y scale factors +real b, c #O cross terms (rotations) +real tx, ty #O x, y offsets + +char ch +int fd, chan, status, wcs_status, zt +real z1, z2 +pointer sp, dir, device, fname, wcstext +int envfind(), strncmp(), open(), fscan(), nscan(), stropen(), iisflu() + +include "iis.com" + +begin + call smark (sp) + call salloc (dir, SZ_PATHNAME, TY_CHAR) + call salloc (fname, SZ_PATHNAME, TY_CHAR) + call salloc (device, SZ_FNAME, TY_CHAR) + call salloc (wcstext, SZ_WCSTEXT, TY_CHAR) + + wcs_status = OK + + # Retrieve the WCS text and open a file descriptor on it. + + if (server == YES) { + # Retrieve the WCS information from a display server. + chan = iisflu(FRTOCHAN(frame)) + + # Cannot use iisio here as the data is byte packed and cannot be + # swapped (while the header still has to be swapped). + + if (iis_version > 0) { + iis_valid = NO + call iishdr (IREAD+PACKED, SZ_WCSTEXT, WCS, 1, 0, chan, 0) + call iisio (Memc[wcstext], SZ_WCSTEXT, status) + if (status > 0) + call strupk (Memc[wcstext], Memc[wcstext], SZ_WCSTEXT) + + iferr (fd = stropen (Memc[wcstext], SZ_WCSTEXT, READ_ONLY)) + fd = NULL + + } else { + call iishdr (IREAD+PACKED, SZ_OLD_WCSTEXT, WCS, 0, 0, chan, 0) + call iisio (Memc[wcstext], SZ_OLD_WCSTEXT, status) + if (status > 0) + call strupk (Memc[wcstext], Memc[wcstext], SZ_OLD_WCSTEXT) + + iferr (fd = stropen (Memc[wcstext], SZ_OLD_WCSTEXT, READ_ONLY)) + fd = NULL + } + + } else { + # Construct the WCS filename, "dir$device_frame.wcs". (Copied from + # the make-WCS code in t_display.x). + + if (envfind ("wcsdir", Memc[dir], SZ_PATHNAME) <= 0) + if (envfind ("WCSDIR", Memc[dir], SZ_PATHNAME) <= 0) + if (envfind ("uparm", Memc[dir], SZ_PATHNAME) <= 0) + call strcpy ("tmp$", Memc[dir], SZ_PATHNAME) + + if (envfind ("stdimage", Memc[device], SZ_FNAME) <= 0) + call strcpy ("display", Memc[device], SZ_FNAME) + + # Get the WCS file filename. + call sprintf (Memc[fname], SZ_PATHNAME, "%s%s_%d.wcs") + call pargstr (Memc[dir]) + if (strncmp (Memc[device], "imt", 3) == 0) + call pargstr ("imtool") + else + call pargstr (Memc[device]) + call pargi (frame) + + if (sz_image > 0) + image[1] = EOS + if (sz_title > 0) + title[1] = EOS + + # Get the saved WCS. + iferr (fd = open (Memc[fname], READ_ONLY, TEXT_FILE)) + fd = NULL + } + + # Decode the WCS from the WCS text. + if (fd != NULL) { + image[1] = EOS + title[1] = EOS + + if (fscan (fd) != EOF) { + # Decode "image - title". + if (sz_image > 0) + call gargwrd (image, sz_image) + if (sz_title > 0) { + call gargwrd (title, sz_title) + repeat { + call gargc (ch) + } until (!IS_WHITE(ch)) + title[1] = ch + call gargstr (title[2], sz_title - 1) + } + + # Decode the WCS information. + if (fscan (fd) != EOF) { + call gargr (a) + call gargr (b) + call gargr (c) + call gargr (d) + call gargr (tx) + call gargr (ty) + call gargr (z1) + call gargr (z2) + call gargi (zt) + if (nscan() == 9) + wcs_status = OK + + if (iis_version > 0) { + if (fscan (fd) != EOF) { + call gargstr (iis_region, SZ_FNAME) + call gargr (iis_sx) + call gargr (iis_sy) + call gargi (iis_snx) + call gargi (iis_sny) + call gargi (iis_dx) + call gargi (iis_dy) + call gargi (iis_dnx) + call gargi (iis_dny) + } + if (nscan() == 9) { + if (fscan (fd) != EOF) + call gargstr (iis_objref, SZ_FNAME) + if (nscan() == 1) + iis_valid = YES + } else + iis_valid = NO + } else { + if (nscan() != 9) { + # Set up the unitary transformation if we + # cannot retrieve the real one. + a = 1.0 + b = 0.0 + c = 0.0 + d = 1.0 + tx = 1.0 + ty = 1.0 + wcs_status = ERR + } + } + } + } + } + + + if (fd != NULL) + call close (fd) + call sfree (sp) + + return (wcs_status) +end diff --git a/pkg/images/tv/display/imdmapfr.x b/pkg/images/tv/display/imdmapfr.x new file mode 100644 index 00000000..745febe2 --- /dev/null +++ b/pkg/images/tv/display/imdmapfr.x @@ -0,0 +1,108 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imset.h> +include <imhdr.h> +include <mach.h> +include <fset.h> +include "display.h" +include "iis.h" + +# IMD_MAPFRAME -- Open the given frame of the stdimage display device on an +# IMIO image descriptor. + +pointer procedure imd_mapframe (frame, mode, select_frame) + +int frame #I frame to be opened [1:N] +int mode #I access mode +int select_frame #I make frame the display frame + +pointer ds +int chan[MAXCHAN] +char device[SZ_FNAME] + +pointer imdmap() +extern imdopen() +int imstati(), fstati(), envgets() +errchk imdmap, imseti +include "iis.com" + +begin + if (envgets ("stdimage", device, SZ_FNAME) == 0) + call error (1, "variable `stdimage' not defined in environment") + + # Pass frame number into IIS code. + call iis_setframe (frame) + + # Map the frame onto an image descriptor. + ds = imdmap (device, mode, imdopen) + # call imseti (ds, IM_CLOSEFD, YES) + chan[1] = fstati (imstati (ds, IM_PIXFD), F_CHANNEL) + chan[2] = MONO + + # Pick up the frame size. + iis_xdim = IM_LEN(ds,1) + iis_ydim = IM_LEN(ds,2) + iis_config = IM_LEN(ds,3) + + # Optimize for sequential i/o. + call imseti (ds, IM_ADVICE, SEQUENTIAL) + + # Display frame being loaded? + if (select_frame == YES) + call zfrmim (chan) + + return (ds) +end + +# IMD_MAPFRAME1 -- Open the given frame of the stdimage display device on an +# IMIO image descriptor. +# This differs from imd_mapframe only in the addition of the erase option. + +pointer procedure imd_mapframe1 (frame, mode, select_frame, erase) + +int frame #I frame to be opened [1:N] +int mode #I access mode +int select_frame #I make frame the display frame +int erase #I erase frame + +pointer ds +int chan[MAXCHAN] +char device[SZ_FNAME] + +pointer imdmap() +extern imdopen() +int imstati(), fstati(), envgets() +errchk imdmap, imseti +include "iis.com" + +begin + if (envgets ("stdimage", device, SZ_FNAME) == 0) + call error (1, "variable `stdimage' not defined in environment") + + # Pass frame number into IIS code. + call iis_setframe (frame) + + # Map the frame onto an image descriptor. + ds = imdmap (device, mode, imdopen) + # call imseti (ds, IM_CLOSEFD, YES) + chan[1] = fstati (imstati (ds, IM_PIXFD), F_CHANNEL) + chan[2] = MONO + + # Pick up the frame size. + iis_xdim = IM_LEN(ds,1) + iis_ydim = IM_LEN(ds,2) + iis_config = IM_LEN(ds,3) + + # Optimize for sequential i/o. + call imseti (ds, IM_ADVICE, SEQUENTIAL) + + # Display frame being loaded? + if (select_frame == YES) + call zfrmim (chan) + + # Erase frame being loaded? + if (erase == YES) + call zersim (chan) + + return (ds) +end diff --git a/pkg/images/tv/display/imdmapping.x b/pkg/images/tv/display/imdmapping.x new file mode 100644 index 00000000..049bef1b --- /dev/null +++ b/pkg/images/tv/display/imdmapping.x @@ -0,0 +1,194 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <ctype.h> +include "iis.h" +include "zdisplay.h" + +.help imd_setmapping, imd_getmapping, imd_query_map +.nf ____________________________________________________________________________ + + Interface routines for setting and getting display server mappings. + + imd_setmapping (region, sx,sy,snx,sny, dx,dy,dnx,dny, objref) + status = imd_getmapping (region, sx,sy,snx,sny, dx,dy,dnx,dny, objref) + status = imd_query_map (wcs, region, sx,sy,snx,sny, dx,dy,dnx,dny, objref) + +The imd_setmapping() procedure should be called prior to an imd_putwcs() +if the mapping information is to be sent with the next WCS write. The +imd_getmapping() function returns a non-zero status if the last WCS query +returned valid mapping information during the read. Both routines depend +upon a previous call to imd_wcsver() (imdmapping.x) to initialize the common +to query the server for this new capability. The imd_query_map() function +returns a non-zero status if a valid mapping is available for the given WCS +number (e.g. the wcs number returned by a cursor read can be entered and +information such as the image name can be returned for the associated mapping). + +.endhelp _______________________________________________________________________ + + +# IMD_SETMAPPING -- Set the mapping information to be sent with the next +# SETWCS command. + +procedure imd_setmapping (reg, sx, sy, snx, sny, dx, dy, dnx, dny, objref) + +char reg[SZ_FNAME] #i region name +real sx, sy #i source raster +int snx, sny +int dx, dy #i destination raster +int dnx, dny +char objref[SZ_FNAME] #i object reference + +bool streq() + +include "iis.com" + +begin + call strcpy (reg, iis_region, SZ_FNAME) + iis_sx = sx + iis_sy = sy + iis_snx = snx + iis_sny = sny + iis_dx = dx + iis_dy = dy + iis_dnx = dnx + iis_dny = dny + + if (streq (objref, "dev$pix")) + call fpathname ("dev$pix.imh", iis_objref, SZ_FNAME) + else + call strcpy (objref, iis_objref, SZ_FNAME) + + iis_valid = YES +end + + +# IMD_GETMAPPING -- Get the mapping information returned with the last +# GETWCS command. + +int procedure imd_getmapping (reg, sx, sy, snx, sny, dx, dy, dnx, dny, objref) + +char reg[SZ_FNAME] #o region name +real sx, sy #o source raster +int snx, sny +int dx, dy #o destination raster +int dnx, dny +char objref[SZ_FNAME] #o object reference + +include "iis.com" + +begin + if (iis_valid == YES) { + call strcpy (iis_region, reg, SZ_FNAME) + sx = iis_sx + sy = iis_sy + snx = iis_snx + sny = iis_sny + dx = iis_dx + dy = iis_dy + dnx = iis_dnx + dny = iis_dny + call strcpy (iis_objref, objref, SZ_FNAME) + } + return (iis_valid) +end + + +# IMD_QUERY_MAP -- Return the mapping information in the server for the +# specified WCS number. + +int procedure imd_query_map (wcs, reg, sx,sy,snx,sny, dx,dy,dnx,dny, objref) + +int wcs #i WCS number of request +char reg[SZ_FNAME] #o region name +real sx, sy #o source raster +int snx, sny +int dx, dy #o destination raster +int dnx, dny +char objref[SZ_FNAME] #o object reference + +pointer sp, wcstext, ip, ds +int fd, frame, chan, status, wcs_status, nl + +int fscan(), stropen(), iisflu() +pointer imd_mapframe1() + +include "iis.com" +define done_ 91 + +begin + call smark (sp) + call salloc (wcstext, SZ_WCSTEXT, TY_CHAR) + call aclrc (Memc[wcstext], SZ_WCSTEXT) + + wcs_status = ERR + iis_valid = NO + frame = wcs / 100 + ds = NULL + + if (iis_version > 0) { + + # If the channel isn't currently open, map the frame temporarily + # so we get a valid read. + if (iisnopen == 0) + ds = imd_mapframe1 (frame, READ_ONLY, NO, NO) + + # Retrieve the WCS information from a display server. + chan = iisflu(FRTOCHAN(frame)) + + # Query the server using the X register to indicate this is + # a "new form" of the WCS query, and pass the requested WCS in + # the T register (which is normally zero). + + call iishdr (IREAD+PACKED, SZ_WCSTEXT, WCS, 1, 0, chan, wcs) + call iisio (Memc[wcstext], SZ_WCSTEXT, status) + if (status > 0) + call strupk (Memc[wcstext], Memc[wcstext], SZ_WCSTEXT) + else + goto done_ + + + # Skip the wcs part of the string, we only want the mapping. + nl = 0 + for (ip=wcstext ; Memc[ip] != NULL; ip=ip+1) { + if (Memc[ip] == '\n') + nl = nl + 1 + if (nl == 2) + break + } + ip = ip + 1 + + # Open the string for reading. + iferr (fd = stropen (Memc[ip], SZ_WCSTEXT, READ_ONLY)) + fd = NULL + + # Decode the Mapping from the WCS text. + if (fd != NULL) { + if (fscan (fd) != EOF) { + call gargwrd (reg, SZ_FNAME) + call gargr (sx) + call gargr (sy) + call gargi (snx) + call gargi (sny) + call gargi (dx) + call gargi (dy) + call gargi (dnx) + call gargi (dny) + + if (fscan (fd) != EOF) { + call gargstr (objref, SZ_FNAME) + wcs_status = OK + iis_valid = YES + } + } + } + + # Close any temporary connection to the server. + if (ds != NULL) + call imunmap (ds) + } + +done_ if (fd != NULL) + call close (fd) + call sfree (sp) + return (wcs_status) +end diff --git a/pkg/images/tv/display/imdopen.x b/pkg/images/tv/display/imdopen.x new file mode 100644 index 00000000..85950270 --- /dev/null +++ b/pkg/images/tv/display/imdopen.x @@ -0,0 +1,16 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <knet.h> + +# IMDOPEN -- Open the image display device as a binary file. + +int procedure imdopen (fname, access_mode) + +char fname[ARB] +int access_mode, fopnbf() +extern zopnim(), zclsim(), zardim(), zawrim(), zawtim(), zsttim() + +begin + return (fopnbf (fname, access_mode, + zopnim, zardim, zawrim, zawtim, zsttim, zclsim)) +end diff --git a/pkg/images/tv/display/imdputwcs.x b/pkg/images/tv/display/imdputwcs.x new file mode 100644 index 00000000..a7b55c8c --- /dev/null +++ b/pkg/images/tv/display/imdputwcs.x @@ -0,0 +1,139 @@ +include <imhdr.h> +include <error.h> +include <imset.h> +include <fset.h> +include "display.h" +include "iis.h" + + +# IMD_PUTWCS -- Write WCS. + +procedure imd_putwcs (ds, frame, str1, str2, a, b, c, d, tx, ty, z1, z2, ztr) +pointer ds #I IMIO descriptor for image display. +int frame #I Frame number for which WCS is to be set. +char str1[ARB] #I First title string (image name). +char str2[ARB] #I Second title string (image title). +real a, d #I x, y scale factors. +real b, c #I cross terms (rotations). +real tx, ty #I x, y offsets. +real z1, z2 #I min and maximum grey scale values. +int ztr #I greyscale transformation code. + +pointer sp, old_wcs, mapping, wcstext, dir, fname, ftemp, device +int wcsfile, server, chan[MAXCHAN] +int fstati(), imstati(), envfind(), open(), strncmp() + +include "iis.com" + +begin + call smark (sp) + call salloc (old_wcs, SZ_WCSTEXT, TY_CHAR) + call salloc (mapping, SZ_WCSTEXT, TY_CHAR) + call salloc (wcstext, SZ_WCSTEXT, TY_CHAR) + + # Format the WCS text. + call sprintf (Memc[old_wcs], SZ_WCSTEXT, + "%s - %s\n%g %g %g %g %g %g %g %g %d\n") + call pargstr (str1) + call pargstr (str2) + call pargr (a) + call pargr (b) + call pargr (c) + call pargr (d) + call pargr (tx) + call pargr (ty) + call pargr (z1) + call pargr (z2) + call pargi (ztr) + + # Add the mapping information if it's valid and we have a capable + # server. + if (iis_version > 0 && iis_valid == YES) { + call sprintf (Memc[mapping], SZ_WCSTEXT, + "%s %g %g %d %d %d %d %d %d\n%s\n") + call pargstr (iis_region) + call pargr (iis_sx) + call pargr (iis_sy) + call pargi (iis_snx) + call pargi (iis_sny) + call pargi (iis_dx) + call pargi (iis_dy) + call pargi (iis_dnx) + call pargi (iis_dny) + call pargstr (iis_objref) + + call sprintf (Memc[wcstext], SZ_WCSTEXT, "%s%s") + call pargstr (Memc[old_wcs]) + call pargstr (Memc[mapping]) + } else + call strcpy (Memc[old_wcs], Memc[wcstext], SZ_OLD_WCSTEXT) + + + # If we are writing to a display server (device has the logical + # cursor capability), output the WCS text via the datastream, + # else use a text file. The datastream set-WCS is also used to + # pass the frame buffer configuration to server devices. + + server = IM_LEN (ds, 4) + + if (server == YES) { + chan[1] = fstati (imstati (ds, IM_PIXFD), F_CHANNEL) + chan[2] = MONO + call imd_setwcs (chan, Memc[wcstext]) + + # Invalidate the mapping once it's been sent. + iis_valid = NO + + } else { + # Construct the WCS filename, "dir$device_frame.wcs". + call salloc (dir, SZ_PATHNAME, TY_CHAR) + call salloc (fname, SZ_PATHNAME, TY_CHAR) + call salloc (ftemp, SZ_PATHNAME, TY_CHAR) + call salloc (device, SZ_FNAME, TY_CHAR) + + if (envfind ("wcsdir", Memc[dir], SZ_PATHNAME) <= 0) + if (envfind ("WCSDIR", Memc[dir], SZ_PATHNAME) <= 0) + if (envfind ("uparm", Memc[dir], SZ_PATHNAME) <= 0) + call strcpy ("tmp$", Memc[dir], SZ_PATHNAME) + + if (envfind ("stdimage", Memc[device], SZ_FNAME) <= 0) + call strcpy ("display", Memc[device], SZ_FNAME) + + # Get a temporary file in the WCS directory. + call sprintf (Memc[ftemp], SZ_PATHNAME, "%swcs") + call pargstr (Memc[dir]) + call mktemp (Memc[ftemp], Memc[ftemp], SZ_PATHNAME) + + # Make the final WCS file filename. + call sprintf (Memc[fname], SZ_PATHNAME, "%s%s_%d.wcs") + call pargstr (Memc[dir]) + if (strncmp (Memc[device], "imt", 3) == 0) + call pargstr ("imtool") + else + call pargstr (Memc[device]) + call pargi (frame) + + # Update the WCS file. + iferr (wcsfile = open (Memc[ftemp], TEMP_FILE, TEXT_FILE)) + call erract (EA_WARN) + else { + # Now delete the old file, if any, and write the new one. + # To avoid process race conditions, create the new file as an + # atomic operation, first writing a new file and then renaming + # it to create the WCS file. + + iferr (call delete (Memc[fname])) + ; + + # Output the file version. + call putline (wcsfile, Memc[wcstext]) + call close (wcsfile) + + # Install the new file. + iferr (call rename (Memc[ftemp], Memc[fname])) + call erract (EA_WARN) + } + } + + call sfree (sp) +end diff --git a/pkg/images/tv/display/imdrcur.x b/pkg/images/tv/display/imdrcur.x new file mode 100644 index 00000000..34148b5b --- /dev/null +++ b/pkg/images/tv/display/imdrcur.x @@ -0,0 +1,117 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <ctype.h> + +# IMDRCUR -- Read the logical image cursor of the named image display device. +# opened with IMDOPEN). This is a high level cursor read, returning image +# pixel coordinates and relying upon the display server to use the keyboard or +# mouse to terminate the cursor read. Nonblocking reads and frame buffer +# coordinates are available as options. The user is expected to select the +# frame for which coordintes are to be returned; the frame number is returned +# in the encoded WCS. The cursor key is returned as the function value. + +int procedure imdrcur (device, x, y, wcs, key, strval, maxch, in_wcs, pause) + +char device[ARB] #I image display device +real x, y #O cursor coords given WCS +int wcs #O WCS of coordinates (frame*100+in_wcs) +int key #O keystroke which triggered cursor read +char strval[maxch] #O optional string value +int maxch #I max chars out +int in_wcs #I desired wcs: 0=frame, 1=image +int pause #I blocking cursor read? (YES|NO) + +char ch +int fd, op +pointer sp, curval, devname, tty, dd, ip + +bool streq() +pointer ttygdes() +int imdopen(), ttygets(), envgets(), nscan(), stg_getline() + +string eof "EOF\n" +string stdimage "stdimage" +errchk ttygdes, imdopen, imdrcuro + +begin + call smark (sp) + call salloc (devname, SZ_FNAME, TY_CHAR) + call salloc (curval, SZ_LINE, TY_CHAR) + call salloc (dd, SZ_LINE, TY_CHAR) + + # Get the logical device name. + if (streq (device, stdimage)) { + if (envgets (stdimage, Memc[devname], SZ_FNAME) <= 0) + call strcpy (device, Memc[devname], SZ_FNAME) + } else + call strcpy (device, Memc[devname], SZ_FNAME) + + # Get the DD kernel driver string for the device. + tty = ttygdes (Memc[devname]) + if (ttygets (tty, "DD", Memc[dd], SZ_LINE) <= 0) + call strcpy (Memc[devname], Memc[dd], SZ_FNAME) + + # Open the device and read the logical image cursor. + fd = imdopen (Memc[dd], READ_WRITE) + call imdrcuro (tty, Memc[curval], SZ_LINE, in_wcs, pause) + + # Decode the formatted cursor value string. + if (streq (Memc[curval], eof)) { + key = EOF + } else { + call sscan (Memc[curval]) + call gargr (x) + call gargr (y) + call gargi (wcs) + call gargc (ch) + call gargstr (Memc[curval], SZ_LINE) + + key = ch + if (nscan() < 4) + key = ERR + + ip = curval + if (nscan() < 5) + Memc[curval] = EOS + else { + while (IS_WHITE(Memc[ip]) || Memc[ip] == '\n') + ip = ip + 1 + } + } + + # In this implementation, string input for colon commands is via the + # terminal to avoid the complexities of character i/o to the display. + # Note that the lower level code can return the string value if it + # chooses to (must be a nonnull string). + + strval[1] = EOS + if (key == ':') { + # String value not already set by imdrcuro? + if (Memc[ip] == EOS) { + call stg_putline (STDOUT, ":") + if (stg_getline (STDIN, Memc[curval]) == EOF) + Memc[curval] = EOS + else + for (ip=curval; IS_WHITE (Memc[ip]); ip=ip+1) + ; + } + + # Copy to the output string argument. + op = 1 + while (Memc[ip] != '\n' && Memc[ip] != EOS) { + strval[op] = Memc[ip] + op = min (op + 1, maxch) + ip = ip + 1 + } + strval[op] = EOS + } + + # Map ctrl/d and ctrl/z onto EOF. + if (key == '\004' || key == '\032') + key = EOF + + call close (fd) + call ttycdes (tty) + + return (key) +end diff --git a/pkg/images/tv/display/imdrcuro.x b/pkg/images/tv/display/imdrcuro.x new file mode 100644 index 00000000..2296fd03 --- /dev/null +++ b/pkg/images/tv/display/imdrcuro.x @@ -0,0 +1,206 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <chars.h> +include <ctype.h> +include "zdisplay.h" +include "iis.h" + +define NEXT_FRAME '\006' +define PREV_FRAME '\022' +define TOGGLE_MARK '\015' + +# IMDRCURO -- Read the logical image cursor from an already opened image +# display device (opened with IMDOPEN). This is a high level cursor read, +# returning image pixel coordinates and relying upon the display server to use +# the keyboard or mouse to terminate the cursor read. Nonblocking reads and +# frame buffer coordinates are available as options. The cursor value is +# returned as an ascii string encoded as follows: +# +# wx wy wcs key [strval] +# +# where WX,WY are the cursor coordinates in the coordinate system defined by +# WCS (= framenumber*100 + wcs, wcs=0 for frame buffer coordinates, wcs=1 for +# image pixel coordinates, the default), KEY is the keystroke used to terminate +# the cursor read, and STRVAL is the string value of the cursor, if key=':' +# (a colon command). Nonprintable keys are returned as octal escapes. + +procedure imdrcuro (tty, outstr, maxch, wcs, pause) + +pointer tty #I graphcap descriptor for device +char outstr[maxch] #O formatted output cursor value +int maxch #I max chars out +int wcs #I desired wcs: 0=framecoords, 1=imagecoords +int pause #I blocking cursor read? (YES|NO) + +short cursor[3] +char key, str[1] +short split[LEN_SPLIT] +pointer sp, strval, imcurval +real a, b, c, d, tx, ty, wx, wy +int status, frame, tid, z, n, keystat, sx, sy, ip, chan, i + +bool mark_cursor +data mark_cursor /false/ + +bool ttygetb() +int rdukey(), ttygeti(), cctoc(), iisflu(), imd_getwcs() +define again_ 91 +include "iis.com" + +begin + call smark (sp) + call salloc (strval, SZ_LINE, TY_CHAR) + call salloc (imcurval, SZB_IMCURVAL, TY_CHAR) + + if (ttygetb (tty, "LC")) { + # Logical image cursor read; the display server supports the + # logical image cursor read as an atomic operation, via the + # logical subunit IMCURSOR (an IRAF special extension to the + # regular IIS datastream protocol). + + if (pause == NO) + tid = IREAD + SAMPLE + else + tid = IREAD + + call iishdr (tid, SZB_IMCURVAL, COMMAND+IMCURSOR, 0,0, wcs, 0) + + call iisio (Memc[imcurval], SZB_IMCURVAL, status) + if (status <= 0) + call strcpy ("EOF\n", outstr, maxch) + else + call strupk (Memc[imcurval], outstr, maxch) + + } else { + # IIS compatible cursor read. Implement the logical cursor read + # using only the primitive IIS cursor functions and the terminal + # driver, accessing the WCS file directly to get the coordinate + # transformation from IIS device coords to image pixel coords. + + # Pick up the frame size and configuration number. + iis_xdim = ttygeti (tty, "xr") + iis_ydim = ttygeti (tty, "yr") + iis_config = ttygeti (tty, "cn") +again_ + if (pause == YES) { + # Enable cursor blink to indicate cursor read in progress. + call iishdr (IWRITE+VRETRACE,1,COMMAND+CURSOR, ADVXONTC, 0,0,0) + cursor[1] = 57B + call iisio (cursor, SZ_SHORT * SZB_CHAR, status) + + # Wait for the user to type a key on the keyboard. The value + # is returned as a newline delimited string. + + keystat = rdukey (Memc[strval], SZ_LINE) + + } else { + Memc[strval] = '\n' + Memc[strval+1] = EOS + keystat = 1 + } + + # Sample the cursor position. + call iisrcr (status, sx, sy) + sx = sx / MCXSCALE + sy = sy / MCYSCALE + + # Determine which frame was being displayed. + call iishdr (IREAD, LEN_SPLIT, COMMAND+LUT, ADVXONTC, 0,0,0) + call iisio (split, LEN_SPLIT * SZB_CHAR, status) + + z = split[1] + if (z == 0) + z = 1 + for (n=1; and(z,1) == 0; z = z / 2) + n = n + 1 + frame = max(1, min(4, n)) + chan = FRTOCHAN(frame) + + if (pause == YES) { + # Turn off cursor blink. + call iishdr (IWRITE+VRETRACE,1,COMMAND+CURSOR, ADVXONTC, 0,0,0) + cursor[1] = 47B + call iisio (cursor, SZ_SHORT * SZB_CHAR, status) + } + + # Decode the trigger keystroke. + ip = 1 + if (cctoc (Memc[strval], ip, key) <= 0) + key = 0 + + # Check for the builtin pseudo "cursor mode" commands. + switch (key) { + case NEXT_FRAME: + # Display the next frame in sequence. + frame = frame + 1 + if (frame > 4) + frame = 1 + chan = IIS_CHAN * DEVCODE + frame + call iisrgb (chan, chan, chan) + goto again_ + case PREV_FRAME: + # Display the previous frame. + frame = frame - 1 + if (frame <= 0) + frame = 1 + chan = IIS_CHAN * DEVCODE + frame + call iisrgb (chan, chan, chan) + goto again_ + case TOGGLE_MARK: + # Toggle the mark cursor enable. + mark_cursor = !mark_cursor + goto again_ + } + + # Mark the cursor position by editing the frame buffer. + if (mark_cursor && keystat > 1 && key != '\004' && key != '\032') { + do i = 1, 3 + cursor[i] = 1 + call achtsb (cursor, cursor, 3) + + call iishdr (IWRITE+BYPASSIFM+PACKED+VRETRACE, 3, REFRESH, + or(sx-1,ADVXONTC), or(sy-1,ADVYONXOV), + iisflu(chan), ALLBITPL) + call iisio (cursor, 3, status) + + call iishdr (IWRITE+BYPASSIFM+PACKED+VRETRACE, 3, REFRESH, + or(sx-1,ADVXONTC), or(sy,ADVYONXOV), + iisflu(chan), ALLBITPL) + call iisio (cursor, 3, status) + + call iishdr (IWRITE+BYPASSIFM+PACKED+VRETRACE, 3, REFRESH, + or(sx-1,ADVXONTC), or(sy+1,ADVYONXOV), + iisflu(chan), ALLBITPL) + call iisio (cursor, 3, status) + } + + # Perform the transformation to image pixel coordinates. + if (wcs != 0) { + if (imd_getwcs (frame,NO, str,0,str,0, a,b,c,d,tx,ty) == ERR) { + call eprintf ("Warning: cannot retrieve WCS for frame %d\n") + call pargi (frame) + } + if (abs(a) > .001) + wx = sx * a + tx + if (abs(d) > .001) + wy = sy * d + ty + } else { + wx = sx + wy = sy + } + + # Format the output cursor value string. + if (keystat == EOF) + call strcpy ("EOF\n", outstr, maxch) + else { + call sprintf (outstr, maxch, "%.6g %.6g %d %s") + call pargr (wx) + call pargr (wy) + call pargi (frame * 100 + wcs) + call pargstr (Memc[strval]) + } + } + + call sfree (sp) +end diff --git a/pkg/images/tv/display/imdsetwcs.x b/pkg/images/tv/display/imdsetwcs.x new file mode 100644 index 00000000..98e8afdc --- /dev/null +++ b/pkg/images/tv/display/imdsetwcs.x @@ -0,0 +1,32 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <knet.h> +include <mach.h> +include "iis.h" + +# IMD_SETWCS -- Pass the WCS for the indicated reference frame to a display +# server. The frame buffer configuration is also passed. + +procedure imd_setwcs (chan, wcstext) + +int chan #I display channel code (frame) +char wcstext[ARB] #I wcs text + +pointer sp, pkwcs +int status, count +int strlen(), iisflu() +include "iis.com" + +begin + count = strlen (wcstext) + 1 + + call smark (sp) + call salloc (pkwcs, count, TY_CHAR) + call strpak (wcstext, Memc[pkwcs], count) + + call iishdr (IWRITE+PACKED, count, WCS, iis_version, 0, iisflu(chan), + max(0,iis_config-1)) + call iisio (Memc[pkwcs], count, status) + + call sfree (sp) +end diff --git a/pkg/images/tv/display/imdwcs.x b/pkg/images/tv/display/imdwcs.x new file mode 100644 index 00000000..66d6b4b5 --- /dev/null +++ b/pkg/images/tv/display/imdwcs.x @@ -0,0 +1,118 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +.help imdwcs +.nf ------------------------------------------------------------------------- +IMDWCS -- Simple interim WCS package for the display interface. This is a +restricted use interface which will be obsoleted by a future interface. + + iw = iw_open (ds, frame, imname, sz_imname, status) + iw_fb2im (iw, fb_x,fb_y, im_x,im_y) + iw_im2fb (iw, im_x,im_y, fb_x,fb_y) + iw_close (iw) + + +This facility uses the WCSDIR file mechanism to retrieve the WCS information +for a display frame. The display name is given by the current value of the +'stdimage' environment variable. Although the WCSDIR info supports a full +2D rotation matrix we recognize only scale and shift terms here. + +NOTE -- The frame buffer coordinates used here are defined in the coordinate +system of the DISPLAY program, IMD_MAPFRAME, etc., i.e., the origin is at the +lower left corner of the frame, and the system is one-indexed. The WCS file, +on the other hand, stores device frame buffer coordinates, which are zero +indexed with the origin at the upper left. +.endhelp -------------------------------------------------------------------- + +define LEN_IWDES 6 + +define IW_A Memr[P2R($1)] # x scale +define IW_B Memr[P2R($1+1)] # cross term (not used) +define IW_C Memr[P2R($1+2)] # cross term (not used) +define IW_D Memr[P2R($1+3)] # y scale +define IW_TX Memr[P2R($1+4)] # x shift +define IW_TY Memr[P2R($1+5)] # y shift + + +# IW_OPEN -- Retrieve the WCS information for the given frame of the stdimage +# display device. If the WCS for the frame cannot be accessed for any reason +# a unitary transformation is returned and wcs_status is set to ERR. Note that +# this is not a hard error, i.e., a valid descriptor is still returned. + +pointer procedure iw_open (ds, frame, imname, sz_imname, wcs_status) + +pointer ds #I display image descriptor +int frame #I frame number for which WCS is desired +char imname[ARB] #O receives name of image loaded into frame (if any) +int sz_imname #I max chars out to imname[]. +int wcs_status #O ERR if WCS cannot be accessed, OK otherwise + +pointer iw +int server +char junk[1] +int imd_getwcs() +errchk calloc + +begin + call calloc (iw, LEN_IWDES, TY_STRUCT) + + # Get the WCS. + server = IM_LEN(ds,4) + wcs_status = imd_getwcs (frame, server, imname, sz_imname, junk,0, + IW_A(iw), IW_B(iw), IW_C(iw), IW_D(iw), IW_TX(iw), IW_TY(iw)) + + # Avoid divide by zero problems if invalid WCS. + if (abs(IW_A(iw)) < .0001 || abs(IW_D(iw)) < .0001) { + + IW_A(iw) = 1.0; IW_D(iw) = 1.0 + IW_TX(iw) = 0.0; IW_TY(iw) = 0.0 + wcs_status = ERR + + } else { + # Convert hardware FB to display FB coordinates. + IW_TY(iw) = IW_TY(iw) + (IW_D(iw) * (IM_LEN(ds,2)-1)) + IW_D(iw) = -IW_D(iw) + } + + return (iw) +end + + +# IW_FB2IM -- Convert frame buffer coordinates to image pixel coordinates. + +procedure iw_fb2im (iw, fb_x,fb_y, im_x,im_y) + +pointer iw #I imd wcs descriptor +real fb_x,fb_y #I frame buffer X,Y coordinates +real im_x,im_y #O image pixel X,Y coordinates + +begin + im_x = (fb_x - 1) * IW_A(iw) + IW_TX(iw) + im_y = (fb_y - 1) * IW_D(iw) + IW_TY(iw) +end + + +# IW_IM2FB -- Convert image pixel coordinates to frame buffer coordinates. + +procedure iw_im2fb (iw, im_x,im_y, fb_x,fb_y) + +pointer iw #I imd wcs descriptor +real im_x,im_y #I image pixel X,Y coordinates +real fb_x,fb_y #O frame buffer X,Y coordinates + +begin + fb_x = (im_x - IW_TX(iw)) / IW_A(iw) + 1 + fb_y = (im_y - IW_TY(iw)) / IW_D(iw) + 1 +end + + +# IW_CLOSE -- Close the IW descriptor. + +procedure iw_close (iw) + +pointer iw #I imd wcs descriptor + +begin + call mfree (iw, TY_STRUCT) +end diff --git a/pkg/images/tv/display/imdwcsver.x b/pkg/images/tv/display/imdwcsver.x new file mode 100644 index 00000000..f8fd9a08 --- /dev/null +++ b/pkg/images/tv/display/imdwcsver.x @@ -0,0 +1,65 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "iis.h" +include "zdisplay.h" + +# IMD_WCSVER -- Query the server for the WCS version supported. A zero +# will be returned for the "old" wcs format used, otherwise the server +# will return a version identifier. + +int procedure imd_wcsver () + +pointer ds +int chan, status, frame, ip +char wcstext[SZ_OLD_WCSTEXT] + +int strncmp(), ctoi(), iisflu() +pointer imd_mapframe1() +bool envgetb() + +include "iis.com" + +begin + iis_valid = NO # initialize + + # Check the environment for a flag to disable the new WCS info. + if (envgetb ("disable_wcs_maps")) { + iis_version = 0 + return (iis_version) + } + + # Open a temporary connection to the server if needed. + ds = NULL + if (iisnopen == 0) + ds = imd_mapframe1 (1, READ_ONLY, NO, NO) + + # Send a WCS query with the X and Y register set. This tells a + # knowledgeable server to reply with a WCS version string, + # otherwise it is a no-op and we get the normal WCS response + # indicating the old format. + + frame = 1 + chan = iisflu (FRTOCHAN(frame)) + call aclrc (wcstext, SZ_OLD_WCSTEXT) + call iishdr (IREAD+PACKED, SZ_OLD_WCSTEXT, WCS, 1, 1, chan, 0) + call iisio (wcstext, SZ_OLD_WCSTEXT, status) + if (status > 0) + call strupk (wcstext, wcstext, SZ_OLD_WCSTEXT) + else { + iis_version = 0 + call imunmap (ds) + return (iis_version) + } + + # Decode the version from the WCS text. + if (strncmp (wcstext, "version=", 8) == 0) { + ip = 9 + status = ctoi (wcstext, ip, iis_version) + } else + iis_version = 0 + + + if (ds != NULL) + call imunmap (ds) + return (iis_version) +end diff --git a/pkg/images/tv/display/maskcolor.x b/pkg/images/tv/display/maskcolor.x new file mode 100644 index 00000000..aa78d77b --- /dev/null +++ b/pkg/images/tv/display/maskcolor.x @@ -0,0 +1,478 @@ +include <ctotok.h> +include <evvexpr.h> +include "ace.h" + +define COLORS "|black|white|red|green|blue|yellow|cyan|magenta|transparent|" +define DEFCOLOR 203 + + +# MASKCOLOR_MAP -- Create the mask colormap object. + +pointer procedure maskcolor_map (colorstring) + +char colorstring #I Color specification string +pointer colors #O Mask colormap object + +int i, j, ip, ncolors, token, lasttoken, maskval1, maskval2, color, offset +int strdic(), ctoi(), nowhite() +pointer sp, str, op + +int coltrans[9] +data coltrans/202,203,204,205,206,207,208,209,-1/ + +define err_ 10 + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + + # If the colorstring is an expression just save the string + # and set the number of colors to 0. + i = nowhite (colorstring, Memc[str], SZ_LINE) + if (Memc[str] == '(') { + call malloc (colors, SZ_LINE, TY_INT) + call malloc (op, LEN_OPERAND, TY_STRUCT) + Memi[colors] = 0 + Memi[colors+1] = op + call strcpy (colorstring, Memc[P2C(colors+2)], SZ_LINE) + O_TYPE(op) = TY_INT + O_VALP(op) = NULL + O_FLAGS(op) = O_FREEOP + # Check expression here. + return (colors) + } + + # Allocate memory for the colormap object. + call malloc (colors, 4*10, TY_INT) + + # Initialize + ncolors = 1 + maskval1 = INDEFI + maskval2 = INDEFI + color = DEFCOLOR + offset = NO + + Memi[colors] = ncolors + Memi[colors+2] = color + Memi[colors+3] = offset + + # Parse the color specification. + token = 0 + call sscan (colorstring) + repeat { + lasttoken = token + call gargtok (token, Memc[str], SZ_LINE) + switch (token) { + case TOK_IDENTIFIER: + call strlwr (Memc[str]) + i = strdic (Memc[str], Memc[str], SZ_LINE, COLORS) + if (i == 0) + goto err_ + color = coltrans[i] + case TOK_NUMBER: + if (lasttoken == TOK_NUMBER) { + if (Memc[str] != '-') + goto err_ + ip = 2 + if (ctoi (Memc[str], ip, maskval2) == 0) + goto err_ + } else { + if (Memc[str] == '+') { + offset = YES + ip = 2 + } else if (Memc[str] == '-') { + offset = YES + ip = 1 + } else + ip = 1 + if (ctoi (Memc[str], ip, color) == 0) + goto err_ + if (lasttoken != TOK_OPERATOR) + maskval2 = color + } + case TOK_OPERATOR: + if (Memc[str] != '=' || lasttoken != TOK_NUMBER) + goto err_ + maskval1 = min (color, maskval2) + maskval2 = max (color, maskval2) + + if (Memc[str+1] == '+') { + call gargtok (token, Memc[str+2], SZ_LINE) + offset = YES + ip = 3 + if (ctoi (Memc[str], ip, color) == 0) + goto err_ + } else if (Memc[str+1] == '-') { + call gargtok (token, Memc[str+2], SZ_LINE) + offset = YES + ip = 2 + if (ctoi (Memc[str], ip, color) == 0) + goto err_ + } + case TOK_PUNCTUATION, TOK_EOS: + if (Memc[str] != ',' && Memc[str] != EOS) + goto err_ + if (!IS_INDEFI(maskval1)) { + do i = 2, ncolors { + j = 4 * i - 4 + if (Memi[colors+j] == maskval1 && + Memi[colors+j+1] == maskval2) + break + } + if (i > ncolors) { + if (mod (ncolors, 10) == 0) + call realloc (colors, 4*(ncolors+10), TY_INT) + ncolors = ncolors + 1 + } + j = 4 * i - 4 + Memi[colors+j] = maskval1 + Memi[colors+j+1] = maskval2 + Memi[colors+j+2] = color + Memi[colors+j+3] = offset + } else { + Memi[colors+2] = color + Memi[colors+3] = offset + } + if (token == TOK_EOS) + break + maskval1 = INDEFI + maskval2 = INDEFI + offset = NO + default: + goto err_ + } + } + + Memi[colors] = ncolors + call sfree (sp) + return (colors) + +err_ + call mfree (colors, TY_INT) + call sfree (sp) + call error (1, "Error in color specifications") +end + + +# MASKCOLOR_FREE -- Free the mask color object. + +procedure maskcolor_free (colors) + +pointer colors #I Mask colormap object + +begin + if (Memi[colors] == 0) + call evvfree (Memi[colors+1]) + call mfree (colors, TY_INT) +end + + +# MASKCOLOR -- Return a color for a mask value. + +int procedure maskcolor (colors, maskval) + +pointer colors #I Mask colormap object +int maskval #I Mask value +int color #O Color value + +int i, j, offset + +begin + # If there is no color array return the mask value. + if (Memi[colors] == 0) + return (maskval) + + color = Memi[colors+2] + offset = Memi[colors+3] + do i = 2, Memi[colors] { + j = 4 * i - 4 + if (maskval >= Memi[colors+j] && maskval <= Memi[colors+j+1]) { + color = Memi[colors+j+2] + offset = Memi[colors+j+3] + break + } + } + + if (offset == YES) + color = maskval + color + return (color) +end + + +procedure maskexprn (colors, maskvals, nmaskvals) + +pointer colors #I Mask colormap object +pointer maskvals #O Pointer to mask values (TY_INT) +int nmaskvals #I Number of mask values + +int i +pointer op, o, evvexpr() +errchk evvexpr + +int locpr +extern maskoperand, maskfunc + +begin + if (Memi[colors] > 0) + return + + op = Memi[colors+1] + O_LEN(op) = nmaskvals + O_VALP(op) = maskvals + + o = evvexpr (Memc[P2C(colors+2)], locpr(maskoperand), op, + locpr(maskfunc), NULL, O_FREEOP) + + #call amovi (Memi[O_VALP(o)], Memi[maskvals], nmaskvals) + switch (O_TYPE(o)) { + case TY_SHORT: + do i = 0, O_LEN(o) { + if (Memi[maskvals+i] > 0) + Memi[maskvals+i] = max (0, Mems[O_VALP(o)+i]) + } + case TY_BOOL, TY_INT: + do i = 0, O_LEN(o) { + if (Memi[maskvals+i] > 0) + Memi[maskvals+i] = max (0, Memi[O_VALP(o)+i]) + } + case TY_REAL: + do i = 0, O_LEN(o) { + if (Memi[maskvals+i] > 0) + Memi[maskvals+i] = max (0, nint(Memr[O_VALP(o)+i])) + } + case TY_DOUBLE: + do i = 0, O_LEN(o) { + if (Memi[maskvals+i] > 0) + Memi[maskvals+i] = max (0, nint(Memd[O_VALP(o)+i])) + } + } + + call evvfree (o) +end + + +# MASKOPERAND -- Handle mask expression operands. + +procedure maskoperand (op, operand, o) + +pointer op #I Input operand pointer +char operand[ARB] #I Operand name +pointer o #O Operand object + +char str[10] +int i, coltrans[9], strdic() +data coltrans/202,203,204,205,206,207,208,209,-1/ + +begin + if (operand[1] == '$') { + call xvv_initop (o, O_LEN(op), O_TYPE(op)) + call amovi (Memi[O_VALP(op)], Memi[O_VALP(o)], O_LEN(op)) + return + } + + call strcpy (operand, str, 11) + call strlwr (str) + i = strdic (str, str, 11, COLORS) + if (i > 0) { + call xvv_initop (o, 0, TY_INT) + O_VALI(o) = coltrans[i] + return + } + + call xvv_error1 ("Unknown mask operand %s", operand) +end + + +define KEYWORDS "|acenum|colors|" + +define F_ACENUM 1 # acenum (maskcodes,[flags]) +define F_COLORS 2 # colors (maskcodes,[col1,col2,col3]) + +# MASKFUNC -- Special processing functions. + +procedure maskfunc (data, func, args, nargs, val) + +pointer data #I client data +char func[ARB] #I function to be called +pointer args[ARB] #I pointer to arglist descriptor +int nargs #I number of arguments +pointer val #O output operand (function value) + +char str[12] +int i, j, c1, c2, c3 +int iresult, optype, oplen, opcode, v_nargs +double dresult + +bool strne() +int strdic(), btoi(), andi() +errchk malloc + +begin + # Lookup the function name in the dictionary. An exact match is + # required (strdic permits abbreviations). Abort if the function + # is not known. + + opcode = strdic (func, str, 12, KEYWORDS) + if (strne (func, str)) + call xvv_error1 ("unknown function `%s' called", func) + + # Verify correct number of arguments. + switch (opcode) { + case F_ACENUM, F_COLORS: + v_nargs = -1 + default: + v_nargs = 1 + } + + if (v_nargs > 0 && nargs != v_nargs) + call xvv_error2 ("function `%s' requires %d arguments", + func, v_nargs) + else if (v_nargs < 0 && nargs < abs(v_nargs)) + call xvv_error2 ("function `%s' requires at least %d arguments", + func, abs(v_nargs)) + + # Group some common operations. + switch (opcode) { + case F_ACENUM: + # Check types of arguments. + if (O_TYPE(args[1]) != TY_INT) + call xvv_error1 ("error in argument types for function `%s'", + func) + if (nargs > 1) { + if (O_TYPE(args[2]) != TY_CHAR) + call xvv_error1 ( + "error in argument types for function `%s'", func) + } + optype = TY_INT + oplen = O_LEN(args[1]) + if (oplen > 0) + call malloc (iresult, oplen, TY_INT) + case F_COLORS: + # Check types of arguments. + do i = 1, nargs { + if (O_TYPE(args[i]) != TY_INT) + call xvv_error1 ("function `%s' requires integer arguments", + func) + } + optype = TY_INT + oplen = O_LEN(args[1]) + if (oplen > 0) + call malloc (iresult, oplen, TY_INT) + } + + # Evaluate the function. + switch (opcode) { + case F_ACENUM: + if (nargs == 1) + call strcpy ("BDEG", str, 12) + else + call strcpy (O_VALC(args[2]), str, 12) + call strupr (str) + c1 = 0; c2 = 0 + for (i=1; str[i]!=EOS; i=i+1) { + switch (str[i]) { + case 'B': + c1 = c1 + MASK_BP + case 'D': + c2 = c2 + MASK_GRW + MASK_SPLIT + case 'E': + c1 = c1 + MASK_BNDRY + case 'F': + c1 = c1 + MASK_BPFLAG + case 'G': + c1 = c1 + MASK_GRW + case 'S': + c1 = c1 + MASK_SPLIT + } + } + + if (oplen == 0) { + i = O_VALI(args[1]) + if (i > 10) { + if (andi(i,c1)!=0 && andi(i,c2)==0) + i = MNUM(i) + else + i = -1 + } else + i = 0 + iresult = i + } else { + do j = 0, oplen-1 { + i = Memi[O_VALP(args[1])+j] + if (i > 10) { + if (andi(i,c1)!=0) + i = MNUM(i) + else if (c2 != 0 && i <= MASK_NUM) + i = MNUM(i) + else + i = -1 + } else + i = 0 + Memi[iresult+j] = i + } + } + case F_COLORS: + c1 = 0; c2 = 204; c3 = 217 + if (nargs > 1) + c1 = O_VALI(args[2]) + if (nargs > 2) { + c2 = O_VALI(args[3]) + c3 = c2 + } + if (nargs > 3) + c3 = O_VALI(args[4]) + if (c3 < c2) { + i = c2; c2 = c3; c3 = i + } + c3 = c3 - c2 + 1 + + optype = TY_INT + oplen = O_LEN(args[1]) + if (oplen == 0) { + i = O_VALI(args[1]) + if (i == 0) + i = c1 + else if (i > 0) + i = c2 + mod (i-1, c3) + iresult = i + } else { + do j = 0, oplen-1 { + i = Memi[O_VALP(args[1])+j] + if (i == 0) + i = c1 + else if (i > 0) + i = c2 + mod (i-1, c3) + Memi[iresult+j] = i + } + } + } + + # Write the result to the output operand. Bool results are stored in + # iresult as an integer value, string results are stored in iresult as + # a pointer to the output string, and integer and real/double results + # are stored in iresult and dresult without any tricks. + + call xvv_initop (val, oplen, optype) + if (oplen == 0) { + switch (optype) { + case TY_BOOL: + O_VALI(val) = btoi (iresult != 0) + case TY_CHAR: + O_VALP(val) = iresult + case TY_INT: + O_VALI(val) = iresult + case TY_REAL: + O_VALR(val) = dresult + case TY_DOUBLE: + O_VALD(val) = dresult + } + } else { + O_VALP(val) = iresult + O_FLAGS(val) = O_FREEVAL + } + + # Free any storage used by the argument list operands. + do i = 1, nargs + call xvv_freeop (args[i]) + +end diff --git a/pkg/images/tv/display/maxmin.x b/pkg/images/tv/display/maxmin.x new file mode 100644 index 00000000..30f281f7 --- /dev/null +++ b/pkg/images/tv/display/maxmin.x @@ -0,0 +1,54 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <imhdr.h> +include "iis.h" + +# MAXMIN -- Get the minimum and maximum pixel values of an image. If valid +# header values are available they are used, otherwise the image is sampled +# on an even grid and the min and max values of this sample are returned. + +procedure maxmin (im, zmin, zmax, nsample_lines) + +pointer im +real zmin, zmax # min and max intensity values +int nsample_lines # amount of image to sample + +int step, ncols, nlines, sample_size, imlines, i +real minval, maxval +pointer imgl2r() +include "iis.com" + +begin + # Only calculate minimum, maximum pixel values if the current + # values are unknown, or if the image was modified since the + # old values were computed. + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + + if (IM_LIMTIME(im) >= IM_MTIME(im)) { + # Use min and max values in image header if they are up to date. + zmin = IM_MIN(im) + zmax = IM_MAX(im) + + } else { + zmin = MAX_REAL + zmax = -MAX_REAL + + # Try to include a constant number of pixels in the sample + # regardless of the image size. The entire image is used if we + # have a small image, and at least sample_lines lines are read + # if we have a large image. + + sample_size = iis_xdim * nsample_lines + imlines = min(nlines, max(nsample_lines, sample_size / ncols)) + step = nlines / (imlines + 1) + + do i = 1 + step, nlines, max (1, step) { + call alimr (Memr[imgl2r(im,i)], ncols, minval, maxval) + zmin = min (zmin, minval) + zmax = max (zmax, maxval) + } + } +end diff --git a/pkg/images/tv/display/mkpkg b/pkg/images/tv/display/mkpkg new file mode 100644 index 00000000..4d6d8885 --- /dev/null +++ b/pkg/images/tv/display/mkpkg @@ -0,0 +1,79 @@ +# Make the DISPLAY libraries. + +$checkout libds.a lib$ +$update libds.a +$checkin libds.a lib$ +$exit + +zzdebug: +zzdebug.e: + $omake zzdebug.x <imhdr.h> + $link zzdebug.o -lds -lstg -o zzdebug.e + ; + +libds.a: + dsmap.x <fset.h> <imset.h> <mach.h> + dspmmap.x <ctype.h> <error.h> <imhdr.h> <imset.h> <mach.h> \ + <pmset.h> + dsulut.x <ctype.h> display.h <error.h> + findz.x iis.com iis.h <imhdr.h> + iisblk.x iis.h <mach.h> zdisplay.h + iiscls.x iis.com iis.h <knet.h> <mach.h> zdisplay.h + iisers.x iis.com iis.h <mach.h> zdisplay.h + iisflu.x iis.h <mach.h> zdisplay.h + iisgop.x iis.h <mach.h> + iishdr.x iis.com iis.h <mach.h> zdisplay.h + iisio.x iis.com iis.h <knet.h> <mach.h> zdisplay.h + iismtc.x iis.h <mach.h> zdisplay.h + iisofm.x iis.h <mach.h> <math.h> zdisplay.h + iisopn.x iis.com iis.h imd.com <knet.h> <mach.h> zdisplay.h + iispio.x iis.com iis.h <knet.h> <mach.h> zdisplay.h + iisrcr.x iis.com iis.h <mach.h> zdisplay.h + iisrd.x iis.com iis.h <mach.h> zdisplay.h + iisrgb.x iis.h <mach.h> zdisplay.h + iissfr.x iis.com iis.h + iisstt.x <fio.h> iis.h <mach.h> zdisplay.h + iiswcr.x iis.com iis.h <mach.h> zdisplay.h + iiswnd.x iis.h <mach.h> zdisplay.h + iiswr.x iis.com iis.h <mach.h> zdisplay.h + iiswt.x iis.com iis.h <knet.h> <mach.h> zdisplay.h + iiszm.x iis.h <mach.h> zdisplay.h + imdgcur.x iis.com iis.h imd.com <knet.h> + imdgetwcs.x <ctype.h> iis.com iis.h zdisplay.h + imdmapfr.x display.h <fset.h> iis.com iis.h <imhdr.h> <imset.h> \ + <mach.h> + imdmapping.x <ctype.h> iis.com iis.h zdisplay.h + imdopen.x <knet.h> + imdputwcs.x display.h <error.h> <fset.h> iis.com iis.h <imhdr.h> \ + <imset.h> + imdrcuro.x <chars.h> <ctype.h> iis.com iis.h <mach.h> zdisplay.h + imdrcur.x <ctype.h> + imdsetwcs.x iis.com iis.h <knet.h> <mach.h> + imdwcsver.x iis.com iis.h zdisplay.h + imdwcs.x <imhdr.h> + maskcolor.x ace.h <ctotok.h> <evvexpr.h> + maxmin.x iis.com iis.h <imhdr.h> <mach.h> + sigl2.x <error.h> <imhdr.h> + sigm2.x <error.h> <imhdr.h> + t_dcontrol.x display.h <fset.h> iis.com iis.h zdisplay.h + t_display.x display.h <error.h> gwindow.h iis.h \ + <imhdr.h> <imset.h> <mach.h> <pmset.h> + zardim.x zdisplay.h + zawrim.x zdisplay.h + zawtim.x zdisplay.h + zblkim.x zdisplay.h + zclrim.x zdisplay.h + zclsim.x zdisplay.h + zersim.x zdisplay.h + zfrmim.x zdisplay.h + zmapim.x zdisplay.h + zmtcim.x zdisplay.h + zopnim.x zdisplay.h + zrcrim.x zdisplay.h + zrgbim.x zdisplay.h + zrmim.x zdisplay.h + zscale.x <ctype.h> <imhdr.h> <imio.h> <imset.h> <pmset.h> + zsttim.x <fio.h> iis.com iis.h <knet.h> + zwndim.x zdisplay.h + zzdebug.x <imhdr.h> + ; diff --git a/pkg/images/tv/display/sigl2.x b/pkg/images/tv/display/sigl2.x new file mode 100644 index 00000000..cbc465ec --- /dev/null +++ b/pkg/images/tv/display/sigl2.x @@ -0,0 +1,976 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> +include <error.h> + +.help sigl2, sigl2_setup +.nf ___________________________________________________________________________ +SIGL2 -- 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 = sigl2_setup (im, x1,x2,nx,xblk, y1,y2,ny,yblk, order) + sigl2_free (si) + ptr = sigl2[sr] (si, linenumber) + +SIGL2_SETUP must be called to set up the transformations after mapping the +image and before performing any scaled i/o to the image. SIGL2_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 # 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_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) + +# SIGL2_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 sigl2_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) && !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) + } 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 + + +# SIGL2_FREE -- Free storage associated with an image opened for scaled +# input. This does not close and unmap the image. + +procedure sigl2_free (si) + +pointer si +int i + +begin + # Free SIGL2 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 + + +# SIGL2S -- 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 sigl2s (si, lineno) + +pointer si # pointer to SI descriptor +int lineno + +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))) + + # 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)) + + 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) + +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 + +real sum +pointer sp, a, b +int nblks_x, nblks_y, ncols, nlines, xoff, i, j +int first_line, nlines_in_sum, npix, nfull_blks, count +pointer imgs2s() +errchk imgs2s + +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)) + return (imgs2s (im, xoff, xoff + npix - 1, y, y)) + + 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") + + 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 = imgs2s (im, xoff, xoff + npix - 1, i, i) + + # 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) { + 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) { + do i = 0, nblks_x-1 + Mems[a+i] = Meml[b+i] / real(nlines_in_sum) + } + + 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 + + +# SIGL2I -- Get a line of type int 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 sigl2i (si, lineno) + +pointer si # pointer to SI descriptor +int lineno + +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_blkavgi() +errchk si_blkavgi + +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_blkavgi (SI_IM(si), x1, x2, int(y), + SI_BAVG(si,1), SI_BAVG(si,2))) + + # 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 + buf_y[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] == 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_blkavgi (SI_IM(si), x1, x2, new_y[i], + SI_BAVG(si,1), SI_BAVG(si,2)) + + 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 { + call aluii (Memi[rawline], Memi[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 awsui (Memi[SI_BUF(si,1)], Memi[SI_BUF(si,2)], + Memi[OUTBUF(si)], npix, weight_1, weight_2) + return (OUTBUF(si)) + } +end + + +# SI_BLKAVGI -- 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_blkavgi (im, x1, x2, y, xbavg, ybavg) + +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 + +real sum +pointer sp, a, b +int nblks_x, nblks_y, ncols, nlines, xoff, i, j +int first_line, nlines_in_sum, npix, nfull_blks, count +pointer imgs2i() +errchk imgs2i + +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)) + return (imgs2i (im, xoff, xoff + npix - 1, y, y)) + + 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") + + 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 = imgs2i (im, xoff, xoff + npix - 1, i, i) + + # Block average line in X. + if (xbavg > 1) { + # First block average only the full blocks. + nfull_blks = npix / xbavg + call abavi (Memi[a], Memi[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 + 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) { + 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) { + do i = 0, nblks_x-1 + Memi[a+i] = Meml[b+i] / real(nlines_in_sum) + } + + call sfree (sp) + return (a) +end + + +# SI_SAMPLEI -- Resample a line via nearest neighbor, rather than linear +# interpolation (ALUI). The calling sequence is the same as for ALUII. + +procedure si_samplei (a, b, x, npix) + +int 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 + + +# SIGL2R -- 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 sigl2r (si, lineno) + +pointer si # pointer to SI descriptor +int lineno + +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))) + + # 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)) + + 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 real. +# 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) + +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 nblks_x, nblks_y, ncols, nlines, xoff, i, j +int first_line, nlines_in_sum, npix, nfull_blks, count +real sum +pointer sp, a, b +pointer imgs2r() +errchk imgs2r + +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)) + return (imgs2r (im, xoff, xoff + npix - 1, y, y)) + + 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. + a = imgs2r (im, xoff, xoff + npix - 1, i, i) + + # 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/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 diff --git a/pkg/images/tv/display/t_dcontrol.x b/pkg/images/tv/display/t_dcontrol.x new file mode 100644 index 00000000..8b68a66b --- /dev/null +++ b/pkg/images/tv/display/t_dcontrol.x @@ -0,0 +1,193 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <syserr.h> +include <fset.h> +include "display.h" +include "zdisplay.h" +include "iis.h" + +# DCONTROL -- Control functions for the image display device. This has been +# cleaned up to eliminate unecessary operations and make it more efficient, +# but is only a throwaway program which breaks a few rules. This file contains +# some explicitly IIS dependent code. + +procedure t_dcontrol() + +real rate +int zoom, type, status +pointer sp, device, devinfo, tty +bool erase, window, rgb_window, blink, match, roam +int red_frame, green_frame, blue_frame, prim_frame, alt_frame, nframes +int red_chan[2], green_chan[2], blue_chan[2], prim_chan[2], alt_chan[2] +char type_string[SZ_FNAME], map_string[SZ_FNAME] +int chan[2], alt1[2], alt2[2] alt3[2] alt4[2] + +real clgetr() +pointer ttygdes() +bool clgetb(), streq(), ttygetb() +int clgeti(), clscan(), nscan(), envgets(), ttygets(), ttygeti(), btoi() +string stdimage "stdimage" +include "iis.com" +define err_ 91 + +begin + call smark (sp) + call salloc (device, SZ_FNAME, TY_CHAR) + call salloc (devinfo, SZ_LINE, TY_CHAR) + + # Get display parameters. + + call clgstr ("type", type_string, SZ_FNAME) + call clgstr ("map", map_string, SZ_FNAME) + + red_frame = clgeti ("red_frame") + green_frame = clgeti ("green_frame") + blue_frame = clgeti ("blue_frame") + prim_frame = clgeti ("frame") + alt_frame = clgeti ("alternate") + + zoom = clgeti ("zoom") + rate = clgetr ("rate") + erase = clgetb ("erase") + window = clgetb ("window") + rgb_window = clgetb ("rgb_window") + blink = clgetb ("blink") + match = clgetb ("match") + roam = clgetb ("roam") + + # Remember current frame. + call clputi ("frame", prim_frame) + call iis_setframe (prim_frame) + + # Get device information. + call clgstr ("device", Memc[device], SZ_FNAME) + if (streq (device, stdimage)) { + if (envgets (stdimage, Memc[device], SZ_FNAME) <= 0) + call syserrs (SYS_ENVNF, stdimage) + } + tty = ttygdes (Memc[device]) + if (ttygets (tty, "DD", Memc[devinfo], SZ_LINE) <= 0) + call error (1, "no `DD' entry in graphcap entry for device") + + # Pick up the frame size and configuration number. + iis_xdim = ttygeti (tty, "xr") + iis_ydim = ttygeti (tty, "yr") + iis_config = ttygeti (tty, "cn") + iis_server = btoi (ttygetb (tty, "LC")) + + # Verify operation is legal on device. + if (iis_server == YES) { + if (!streq (type_string, "frame")) + goto err_ + if (!streq (map_string, "mono")) + goto err_ + if (erase) + ; + if (roam) + goto err_ + if (window) + goto err_ + if (rgb_window) + goto err_ + if (blink) + goto err_ + if (match) { +err_ call eprintf ("operation not supported for display device %s\n") + call pargstr (Memc[device]) + call ttycdes (tty) + call sfree (sp) + return + } + } + + # Access display. + call strpak (Memc[devinfo], Memc[devinfo], SZ_LINE) + call iisopn (Memc[devinfo], READ_WRITE, chan) + if (chan[1] == ERR) + call error (2, "cannot open display") + + call fseti (STDOUT, F_FLUSHNL, YES) + + red_chan[1] = FRTOCHAN(red_frame) + green_chan[1] = FRTOCHAN(green_frame) + blue_chan[1] = FRTOCHAN(blue_frame) + prim_chan[1] = FRTOCHAN(prim_frame) + alt_chan[1] = FRTOCHAN(alt_frame) + + red_chan[2] = MONO + green_chan[2] = MONO + blue_chan[2] = MONO + prim_chan[2] = MONO + alt_chan[2] = MONO + + # Execute the selected control functions. + if (streq (type_string, "rgb")) { + type = RGB + call zrgbim (red_chan, green_chan, blue_chan) + } else if (streq (type_string, "frame")) { + type = FRAME + call zfrmim (prim_chan) + } else + call error (3, "unknown display type") + + # Set display mapping. + call zmapim (prim_chan, map_string) + + if (erase) { + switch (type) { + case RGB: + call zersim (red_chan) + call zersim (green_chan) + call zersim (blue_chan) + case FRAME: + call zersim (prim_chan) + } + + } else { + if (roam) { + call printf ("Roam display and exit by pushing any button\n") + call zrmim (prim_chan, zoom) + } + + if (window) { + call printf ("Window display and exit by pushing any button\n") + call zwndim (prim_chan) + } + + if (rgb_window) { + call printf ("Window display and exit by pushing any button\n") + call zwndim3 (red_chan, green_chan, blue_chan) + } + + if (match) + call zmtcim (alt_chan, prim_chan) + + if (blink) { + if (clscan ("alternate") != EOF) { + call gargi (alt1[1]) + call gargi (alt2[1]) + call gargi (alt3[1]) + call gargi (alt4[1]) + nframes = nscan() + + alt1[1] = FRTOCHAN(alt1[1]) + alt2[1] = FRTOCHAN(alt2[1]) + alt3[1] = FRTOCHAN(alt3[1]) + alt4[1] = FRTOCHAN(alt4[1]) + + alt1[2] = MONO + alt2[2] = MONO + alt3[2] = MONO + alt4[2] = MONO + + call printf ("Exit by pushing any button\n") + call zblkim (alt1, alt2, alt3, alt4, nframes, rate) + } + } + } + + # Close display. + call zclsim (chan[1], status) + call ttycdes (tty) + call sfree (sp) +end diff --git a/pkg/images/tv/display/t_display.x b/pkg/images/tv/display/t_display.x new file mode 100644 index 00000000..f4156f39 --- /dev/null +++ b/pkg/images/tv/display/t_display.x @@ -0,0 +1,885 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <mach.h> +include <imset.h> +include <imhdr.h> +include <error.h> +include <pmset.h> +include "display.h" +include "gwindow.h" +include "iis.h" + +# DISPLAY - Display an image. The specified image section is mapped into +# the specified section of an image display frame. The mapping involves +# a linear transformation in X and Y and a linear or logarithmic transformation +# in Z (greyscale). Images of all pixel datatypes are supported, and there +# no upper limit on the size of an image. The display device is interfaced +# to FIO as a file and is accessed herein via IMIO as just another imagefile. +# The physical characteristics of the display (i.e., X, Y, and Z resolution) +# are taken from the image header. The display frame buffer is the pixel +# storage "file". + +procedure t_display() + +char image[SZ_FNAME] # Image to display +int frame # Display frame +int erase # Erase frame? + +int i +pointer sp, wdes, im, ds + +bool clgetb() +int clgeti(), btoi(), imd_wcsver(), imtlen(), imtgetim() +pointer immap(), imd_mapframe1(), imtopenp() +errchk immap, imd_mapframe1 +errchk ds_getparams, ds_setwcs, ds_load_display, ds_erase_border + +begin + call smark (sp) + call salloc (wdes, LEN_WDES, TY_STRUCT) + call aclri (Memi[wdes], LEN_WDES) + + # Open input imagefile. + im = imtopenp ("image") + if (imtlen (im) != 1) + call error (1, "Only one image may be displayed") + i = imtgetim (im, image, SZ_FNAME) + call imtclose (im) + #call clgstr ("image", image, SZ_FNAME) + im = immap (image, READ_ONLY, 0) + if (IM_NDIM(im) <= 0) + call error (1, "image has no pixels") + + # Query server to get the WCS version, this also tells us whether + # we can use the all 16 supported frames. + if (imd_wcsver() == 0) + call clputi ("display.frame.p_max", 4) + else + call clputi ("display.frame.p_max", 16) + + + # Open display device as an image. + frame = clgeti ("frame") + W_FRAME(wdes) = frame + + erase = btoi (clgetb ("erase")) + if (erase == YES) + ds = imd_mapframe1 (frame, WRITE_ONLY, + btoi (clgetb ("select_frame")), erase) + else + ds = imd_mapframe1 (frame, READ_WRITE, + btoi (clgetb ("select_frame")), erase) + + # Get display parameters and set up transformation. + call ds_getparams (im, ds, wdes) + + # Compute and output the screen to image pixel WCS. + call ds_setwcs (im, ds, wdes, image, frame) + + # Display the image and zero the border if necessary. + call ds_load_display (im, ds, wdes) + if (!clgetb ("erase") && clgetb ("border_erase")) + call ds_erase_border (im, ds, wdes) + + # Free storage. + call maskcolor_free (W_OCOLORS(wdes)) + call maskcolor_free (W_BPCOLORS(wdes)) + do i = 0, W_MAXWC + if (W_UPTR(W_WC(wdes,i)) != NULL) + call ds_ulutfree (W_UPTR(W_WC(wdes,i))) + call imunmap (ds) + call imunmap (im) + + call sfree (sp) +end + + +# DS_GETPARAMS -- Get the parameters controlling how the image is mapped +# into the display frame. Set up the transformations and save in the graphics +# descriptor file. If "repeat" mode is enabled, read the graphics descriptor +# file and reuse the transformations therein. + +procedure ds_getparams (im, ds, wdes) + +pointer im, ds, wdes #I Image, display, and graphics descriptors + +bool fill, zscale_flag, zrange_flag, zmap_flag +real xcenter, ycenter, xsize, ysize +real xmag, ymag, xscale, yscale, pxsize, pysize +real z1, z2, contrast +int nsample, ncols, nlines +pointer wnwin, wdwin, wwwin, wipix, wdpix, zpm, bpm +pointer sp, str, ztrans, lutfile + +int clgeti(), clgwrd(), nowhite() +real clgetr() +pointer maskcolor_map(), ds_pmmap(), zsc_pmsection() +pointer ds_ulutalloc() +bool streq(), clgetb() +errchk maskcolor_map, ds_pmmap, zsc_pmsection, mzscale + +begin + call smark (sp) + call salloc (str, SZ_LINE, TY_CHAR) + call salloc (ztrans, SZ_FNAME, TY_CHAR) + + # Get overlay mask and colors. + call clgstr ("overlay", W_OVRLY(wdes), W_SZSTRING) + call clgstr ("ocolors", Memc[str], SZ_LINE) + W_OCOLORS(wdes) = maskcolor_map (Memc[str]) + + # Get bad pixel mask. + call clgstr ("bpmask", W_BPM(wdes), W_SZSTRING) + W_BPDISP(wdes) = clgwrd ("bpdisplay", Memc[str], SZ_LINE, BPDISPLAY) + call clgstr ("bpcolors", Memc[str], SZ_LINE) + W_BPCOLORS(wdes) = maskcolor_map (Memc[str]) + + # Determine the display window into which the image is to be mapped + # in normalized device coordinates. + + xcenter = max(0.0, min(1.0, clgetr ("xcenter"))) + ycenter = max(0.0, min(1.0, clgetr ("ycenter"))) + xsize = max(0.0, min(1.0, clgetr ("xsize"))) + ysize = max(0.0, min(1.0, clgetr ("ysize"))) + + # Set up a new graphics descriptor structure defining the coordinate + # transformation used to map the image into the display frame. + + wnwin = W_WC(wdes,W_NWIN) + wdwin = W_WC(wdes,W_DWIN) + wwwin = W_WC(wdes,W_WWIN) + wipix = W_WC(wdes,W_IPIX) + wdpix = W_WC(wdes,W_DPIX) + + # Determine X and Y scaling ratios required to map the image into the + # normalized display window. If spatial scaling is not desired filling + # must be disabled and XMAG and YMAG must be set to 1.0 in the + # parameter file. Fill mode will always produce an aspect ratio of 1; + # if nonequal scaling is required then the magnification ratios must + # be set explicitly by the user. + + fill = clgetb ("fill") + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + if (fill) { + # Compute scale in units of window coords per data pixel required + # to scale image to fit window. + + xmag = (IM_LEN(ds,1) * xsize) / ncols + ymag = (IM_LEN(ds,2) * ysize) / nlines + + if (xmag > ymag) + xmag = ymag + else + ymag = xmag + + } else { + # Compute scale required to provide image magnification ratios + # specified by the user. Magnification is specified in units of + # display pixels, i.e, a magnification ratio of 1.0 means that + # image pixels will map to display pixels without scaling. + + xmag = clgetr ("xmag") + ymag = clgetr ("ymag") + } + + xscale = 1.0 / (IM_LEN(ds,1) / xmag) + yscale = 1.0 / (IM_LEN(ds,2) / ymag) + + # Set device window limits in normalized device coordinates. + # World coord system 0 is used for the device window. + + W_XS(wnwin) = xcenter - xsize / 2.0 + W_XE(wnwin) = xcenter + xsize / 2.0 + W_YS(wnwin) = ycenter - ysize / 2.0 + W_YE(wnwin) = ycenter + ysize / 2.0 + + # Set pixel coordinates of window. + # If the image is too large to fit in the window given the scaling + # factors XSCALE and YSCALE, the following will set starting and ending + # pixel coordinates in the interior of the image. If the image is too + # small to fill the window then the pixel coords will reference beyond + # the bounds of the image. Note that the 0.5 is because NDC has + # the screen corner at 0 while screen pixels have the corner at 0.5. + + pxsize = xsize / xscale + pysize = ysize / yscale + + W_XS(wdwin) = (ncols / 2.0) - (pxsize / 2.0) + 0.5 + W_XE(wdwin) = W_XS(wdwin) + pxsize + W_YS(wdwin) = (nlines / 2.0) - (pysize / 2.0) + 0.5 + W_YE(wdwin) = W_YS(wdwin) + pysize + + # Compute X and Y magnification ratios required to map image into + # the device window in device pixel units. + + xmag = (W_XE(wnwin)-W_XS(wnwin))*IM_LEN(ds,1)/(W_XE(wdwin)-W_XS(wdwin)) + ymag = (W_YE(wnwin)-W_YS(wnwin))*IM_LEN(ds,2)/(W_YE(wdwin)-W_YS(wdwin)) + + # Compute the coordinates of the image section to be displayed. + # Round down if upper pixel is exactly at one-half. + + W_XS(wipix) = max (1, nint(W_XS(wdwin))) + W_XE(wipix) = min (ncols, nint(W_XE(wdwin)-1.01)) + W_YS(wipix) = max (1, nint(W_YS(wdwin))) + W_YE(wipix) = min (nlines, nint(W_YE(wdwin)-1.01)) + + # Now compute the image and display pixels to be used. + # The image may be truncated to fit in the display window. + # These are integer coordinates at the pixel centers. + + pxsize = W_XE(wipix) - W_XS(wipix) + 1 + pysize = W_YE(wipix) - W_YS(wipix) + 1 + xcenter = (W_XE(wnwin) + W_XS(wnwin)) / 2.0 * IM_LEN(ds,1) + 0.5 + ycenter = (W_YE(wnwin) + W_YS(wnwin)) / 2.0 * IM_LEN(ds,2) + 0.5 + + #W_XS(wdpix) = max (1, nint (xcenter - (pxsize/2.0*xmag) + 0.5)) + W_XS(wdpix) = max (1, int (xcenter - (pxsize/2.0*xmag) + 0.5)) + W_XE(wdpix) = min (IM_LEN(ds,1), nint (W_XS(wdpix)+pxsize*xmag - 1.01)) + #W_YS(wdpix) = max (1, nint (ycenter - (pysize/2.0*ymag) + 0.5)) + W_YS(wdpix) = max (1, int (ycenter - (pysize/2.0*ymag) + 0.5)) + W_YE(wdpix) = min (IM_LEN(ds,2), nint (W_YS(wdpix)+pysize*ymag - 1.01)) + + # Now adjust the display window to be consistent with the image and + # display pixels to be used. + + W_XS(wdwin) = W_XS(wnwin) * IM_LEN(ds,1) + 0.5 + W_XE(wdwin) = W_XE(wnwin) * IM_LEN(ds,1) + 0.5 + W_YS(wdwin) = W_YS(wnwin) * IM_LEN(ds,2) + 0.5 + W_YE(wdwin) = W_YE(wnwin) * IM_LEN(ds,2) + 0.5 + W_XS(wdwin) = (W_XS(wipix)-0.5) + (W_XS(wdwin)-(W_XS(wdpix)-0.5))/xmag + W_XE(wdwin) = (W_XS(wipix)-0.5) + (W_XE(wdwin)-(W_XS(wdpix)-0.5))/xmag + W_YS(wdwin) = (W_YS(wipix)-0.5) + (W_YS(wdwin)-(W_YS(wdpix)-0.5))/ymag + W_YE(wdwin) = (W_YS(wipix)-0.5) + (W_YE(wdwin)-(W_YS(wdpix)-0.5))/ymag + + # Order of interpolator used for spatial transformation. + W_XT(wdwin) = max(0, min(1, clgeti ("order"))) + W_YT(wdwin) = W_XT(wdwin) + + # Determine the greyscale transformation. + call clgstr ("ztrans", Memc[ztrans], SZ_FNAME) + if (streq (Memc[ztrans], "log")) + W_ZT(wdwin) = W_LOG + else if (streq (Memc[ztrans], "linear")) + W_ZT(wdwin) = W_LINEAR + else if (streq (Memc[ztrans], "none")) + W_ZT(wdwin) = W_UNITARY + else if (streq (Memc[ztrans], "user")) { + W_ZT(wdwin) = W_USER + call salloc (lutfile, SZ_FNAME, TY_CHAR) + call clgstr ("lutfile", Memc[lutfile], SZ_FNAME) + W_UPTR(wdwin) = ds_ulutalloc (Memc[lutfile], z1, z2) + } else { + call eprintf ("Bad greylevel transformation '%s'\n") + call pargstr (Memc[ztrans]) + W_ZT(wdwin) = W_LINEAR + } + + # The zscale, and zrange parameters determine the algorithms for + # determining Z1 and Z2, the range of input z values to be mapped + # into the fixed range of display greylevels. If sampling and no + # sample mask is given then create one as a subsampled image section. + # If greyscale mapping is disabled the zscale and zrange options are + # disabled. Greyscale mapping can also be disabled by turning off + # zscale and zrange and setting Z1 and Z2 to the device greyscale min + # and max values, producing a unitary transformation. + + if (W_ZT(wdwin) == W_UNITARY || W_ZT(wdwin) == W_USER) { + zscale_flag = false + zrange_flag = false + zmap_flag = false + } else { + zmap_flag = true + zscale_flag = clgetb ("zscale") + if (!zscale_flag) + zrange_flag = clgetb ("zrange") + } + + if (zscale_flag || (zrange_flag && IM_LIMTIME(im) < IM_MTIME(im))) { + call clgstr ("zmask", W_ZPM(wdes), W_SZSTRING) + nsample = max (100, clgeti ("nsample")) + if (nowhite (W_ZPM(wdes), W_ZPM(wdes), W_SZSTRING) > 0) { + if (W_ZPM(wdes) == '[') + zpm = zsc_pmsection (W_ZPM(wdes), im) + else + zpm = ds_pmmap (W_ZPM(wdes), im) + } else + zpm = NULL + iferr (bpm = ds_pmmap (W_BPM(wdes), im)) { + call erract (EA_WARN) + bpm = NULL + } + } + + if (zscale_flag) { + # Autoscaling is desired. Compute Z1 and Z2 which straddle the + # median computed by sampling a portion of the image. + + contrast = clgetr ("contrast") + call mzscale (im, zpm, bpm, contrast, nsample, z1, z2) + if (zpm != NULL) + call imunmap (zpm) + if (bpm != NULL) + call imunmap (bpm) + + } else if (zrange_flag) { + # Use the limits in the header if current otherwise get the + # minimum and maximum of the sample mask. + if (IM_LIMTIME(im) >= IM_MTIME(im)) { + z1 = IM_MIN(im) + z2 = IM_MAX(im) + } else { + call mzscale (im, zpm, bpm, 0., nsample, z1, z2) + if (zpm != NULL) + call imunmap (zpm) + if (bpm != NULL) + call imunmap (bpm) + } + + } else if (zmap_flag) { + z1 = clgetr ("z1") + z2 = clgetr ("z2") + } else { + z1 = IM_MIN(ds) + z2 = IM_MAX(ds) + } + + W_ZS(wdwin) = z1 + W_ZE(wdwin) = z2 + + call printf ("z1=%g z2=%g\n") + call pargr (z1) + call pargr (z2) + call flush (STDOUT) + + # The user world coordinate system should be set from the CTRAN + # structure in the image header, but for now we just make it equal + # to the pixel coordinate system. + + call amovi (Memi[wdwin], Memi[wwwin], LEN_WC) + W_UPTR(wwwin) = NULL # should not copy pointers!! + call sfree (sp) +end + + +# DS_SETWCS -- Compute the rotation matrix needed to convert screen coordinates +# (zero indexed, y-flipped) to image pixel coordinates, allowing both for the +# transformation from screen space to the image section being displayed, and +# from the image section to the physical input image. +# +# NOTE -- This code assumes that the display device is zero-indexed and +# y-flipped; this is usually the case, but should be parameterized in the +# graphcap. This code also assumes that the full device screen is being used, +# and that we are not assigning multiple WCS to different regions of the screen. + +procedure ds_setwcs (im, ds, wdes, image, frame) + +pointer im, ds, wdes # image, display, and coordinate descriptors +char image[SZ_FNAME] # image section name +int frame # frame + +real a, b, c, d, tx, ty +int ip, i, j, axis[2] +real sx, sy +int dx, dy, snx, sny, dnx, dny +pointer sp, imname, title, wnwin, wdwin +pointer src, dest, region, objref +long lv[IM_MAXDIM], pv1[IM_MAXDIM], pv2[IM_MAXDIM] + +bool streq() + +begin + call smark (sp) + call salloc (imname, SZ_FNAME, TY_CHAR) + call salloc (title, SZ_LINE, TY_CHAR) + call salloc (region, SZ_FNAME, TY_CHAR) + call salloc (objref, SZ_FNAME, TY_CHAR) + + # Compute the rotation matrix needed to transform screen pixel coords + # to image section coords. + + wnwin = W_WC(wdes,W_NWIN) + wdwin = W_WC(wdes,W_DWIN) + + # X transformation. + a = (W_XE(wdwin)-W_XS(wdwin))/((W_XE(wnwin)-W_XS(wnwin))*IM_LEN(ds,1)) + c = 0.0 # not rotated, cross term is zero + tx = W_XS(wdwin) - a * (W_XS(wnwin) * IM_LEN(ds,1)) + + # Y transformation. + b = 0.0 # not rotated, cross term is zero + d = (W_YE(wdwin)-W_YS(wdwin))/((W_YE(wnwin)-W_YS(wnwin))*IM_LEN(ds,2)) + ty = W_YS(wdwin) - d * (W_YS(wnwin) * IM_LEN(ds,2)) + + # Now allow for the Y-flip (origin at upper left in display window). + d = -d + ty = W_YE(wdwin) - d * ((1.0 - W_YE(wnwin)) * IM_LEN(ds,2)) + + # Now translate the screen corner to the center of the screen pixel. + tx = tx + 0.5 * a + ty = ty + 0.5 * d + + # Determine the logical to physical mapping by evaluating two points. + # and determining the axis reduction if any. pv1 will be the + # offset and pv2-pv1 will be the scale. + + call aclrl (pv1, IM_MAXDIM) + call aclrl (lv, IM_MAXDIM) + call imaplv (im, lv, pv1, 2) + call amovkl (long(1), lv, IM_MAXDIM) + call aclrl (pv2, IM_MAXDIM) + call imaplv (im, lv, pv2, 2) + + i = 1 + axis[1] = 1; axis[2] = 2 + do j = 1, IM_MAXDIM + if (pv1[j] != pv2[j]) { + axis[i] = j + i = i + 1 + } + + pv2[axis[1]] = (pv2[axis[1]] - pv1[axis[1]]) + pv2[axis[2]] = (pv2[axis[2]] - pv1[axis[2]]) + + # These imply a new rotation matrix which we won't bother to work out + # separately here. Multiply the two rotation matrices and add the + # translation vectors to get the overall transformation from screen + # coordinates to image coordinates. + a = a * pv2[axis[1]] + d = d * pv2[axis[2]] + tx = tx * pv2[axis[1]] + pv1[axis[1]] + ty = ty * pv2[axis[2]] + pv1[axis[2]] + + # Get the image name (minus image section) and + # title string (minus any newline. + call ds_gimage (im, image, Memc[imname], SZ_FNAME) + call strcpy (IM_TITLE(im), Memc[title], SZ_LINE) + for (ip=title; Memc[ip] != '\n' && Memc[ip] != EOS; ip=ip+1) + ; + Memc[ip] = EOS + + + # Define the mapping from the image pixels to frame buffer pixels. + src = W_WC(wdes,W_IPIX) + sx = W_XS(src) + sy = W_YS(src) + snx = (W_XE(src) - W_XS(src) + 1) + sny = (W_YE(src) - W_YS(src) + 1) + + dest = W_WC(wdes,W_DPIX) + dx = W_XS(dest) + dy = W_YS(dest) + dnx = (W_XE(dest) - W_XS(dest) + 1) + dny = (W_YE(dest) - W_YS(dest) + 1) + + # For a single image display the 'region' is fixed. The object ref + # is the fully defined image node!prefix path, including any sections. + # We need a special kludge to keep backward compatability with the + # use of "dev$pix" as the standard test image name. + call strcpy ("image", Memc[region], SZ_FNAME) + if (streq (image, "dev$pix")) + call fpathname ("dev$pix.imh", Memc[objref], SZ_PATHNAME) + else + call fpathname (image, Memc[objref], SZ_PATHNAME) + + # Add the mapping info to be written with the WCS. + call imd_setmapping (Memc[region], sx, sy, snx, sny, + dx, dy, dnx, dny, Memc[objref]) + + # Write the WCS. + call imd_putwcs (ds, frame, Memc[imname], Memc[title], + a, b, c, d, tx, ty, W_ZS(wdwin), W_ZE(wdwin), W_ZT(wdwin)) + + call sfree (sp) +end + + +# DS_GIMAGE -- Convert input image section name to a 2D physical image section. + +procedure ds_gimage (im, input, output, maxchar) + +pointer im #I IMIO pointer +char input[ARB] #I Input image name +char output[maxchar] #O Output image name +int maxchar #I Maximum characters in output name. + +int i, fd +pointer sp, section, lv, pv1, pv2 + +int stropen(), strlen() +bool streq() + +begin + call smark (sp) + call salloc (section, SZ_FNAME, TY_CHAR) + call salloc (lv, IM_MAXDIM, TY_LONG) + call salloc (pv1, IM_MAXDIM, TY_LONG) + call salloc (pv2, IM_MAXDIM, TY_LONG) + + # Get endpoint coordinates in original image. + call amovkl (long(1), Meml[lv], IM_MAXDIM) + call aclrl (Meml[pv1], IM_MAXDIM) + call imaplv (im, Meml[lv], Meml[pv1], 2) + call amovl (IM_LEN(im,1), Meml[lv], IM_NDIM(im)) + call aclrl (Meml[pv2], IM_MAXDIM) + call imaplv (im, Meml[lv], Meml[pv2], 2) + + # Set image section. + fd = stropen (Memc[section], SZ_FNAME, NEW_FILE) + call fprintf (fd, "[") + do i = 1, IM_MAXDIM { + if (Meml[pv1+i-1] != Meml[pv2+i-1]) + call fprintf (fd, "*") + else if (Meml[pv1+i-1] != 0) { + call fprintf (fd, "%d") + call pargi (Meml[pv1+i-1]) + } else + break + call fprintf (fd, ",") + } + call close (fd) + i = strlen (Memc[section]) + Memc[section+i-1] = ']' + + if (streq ("[*,*]", Memc[section])) + Memc[section] = EOS + + # Strip existing image section and add new section. +# call imgimage (input, output, maxchar) +# call strcat (Memc[section], output, maxchar) + + if (Memc[section] == EOS) + call imgimage (input, output, maxchar) + else + call strcpy (input, output, maxchar) + + call sfree (sp) +end + + +# DS_LOAD_DISPLAY -- Map an image into the display window. In general this +# involves independent linear transformations in the X, Y, and Z (greyscale) +# dimensions. If a spatial dimension is larger than the display window then +# the image is block averaged. If a spatial dimension or a block averaged +# dimension is smaller than the display window then linear interpolation is +# used to expand the image. Both the input image and the output device appear +# to us as images, accessed via IMIO. All spatial scaling is +# handled by the "scaled input" package, i.e., SIGM2[SR]. Our task is to +# get lines from the scaled input image, transform the greyscale if necessary, +# and write the lines to the output device. + +procedure ds_load_display (im, ds, wdes) + +pointer im # input image +pointer ds # output image +pointer wdes # graphics window descriptor + +real z1, z2, dz1, dz2, px1, px2, py1, py2 +int i, order, zt, wx1, wx2, wy1, wy2, wy, nx, ny, xblk, yblk, color +pointer wdwin, wipix, wdpix, ovrly, bpm, pm, uptr +pointer in, out, si, si_ovrly, si_bpovrly, ocolors, bpcolors, rtemp +bool unitary_greyscale_transformation +short lut1, lut2, dz1_s, dz2_s, z1_s, z2_s + +bool fp_equalr() +int imstati(), maskcolor() +pointer ds_pmmap(), imps2s(), imps2r() +pointer sigm2s(), sigm2i(), sigm2r(), sigm2_setup() +errchk ds_pmmap, imps2s, imps2r, sigm2s, sigm2i, sigm2r, sigm2_setup +errchk maskexprn + +begin + wdwin = W_WC(wdes,W_DWIN) + wipix = W_WC(wdes,W_IPIX) + wdpix = W_WC(wdes,W_DPIX) + + # Set image and display pixels. + px1 = nint (W_XS(wipix)) + px2 = nint (W_XE(wipix)) + py1 = nint (W_YS(wipix)) + py2 = nint (W_YE(wipix)) + wx1 = nint (W_XS(wdpix)) + wx2 = nint (W_XE(wdpix)) + wy1 = nint (W_YS(wdpix)) + wy2 = nint (W_YE(wdpix)) + + z1 = W_ZS(wdwin) + z2 = W_ZE(wdwin) + zt = W_ZT(wdwin) + uptr = W_UPTR(wdwin) + order = max (W_XT(wdwin), W_YT(wdwin)) + + # Setup scaled input and masks. + si = NULL + si_ovrly = NULL + si_bpovrly = NULL + nx = wx2 - wx1 + 1 + ny = wy2 - wy1 + 1 + xblk = INDEFI + yblk = INDEFI + + ocolors = W_OCOLORS(wdes) + iferr (ovrly = ds_pmmap (W_OVRLY(wdes), im)) { + call erract (EA_WARN) + ovrly = NULL + } + if (ovrly != NULL) { + xblk = INDEFI + yblk = INDEFI + si_ovrly = sigm2_setup (ovrly, NULL, px1,px2,nx,xblk, + py1,py2,ny,yblk, -1) + } + + bpcolors = W_BPCOLORS(wdes) + switch (W_BPDISP(wdes)) { + case BPDNONE: + si = sigm2_setup (im, NULL, px1,px2,nx,xblk, py1,py2,ny,yblk, order) + case BPDOVRLY: + si = sigm2_setup (im, NULL, px1,px2,nx,xblk, py1,py2,ny,yblk, order) + iferr (bpm = ds_pmmap (W_BPM(wdes), im)) + bpm = NULL + if (bpm != NULL) + si_bpovrly = sigm2_setup (bpm, NULL, px1,px2,nx,xblk, + py1,py2,ny,yblk, -1) + case BPDINTERP: + iferr (bpm = ds_pmmap (W_BPM(wdes), im)) + bpm = NULL + if (bpm != NULL) + pm = imstati (bpm, IM_PMDES) + else + pm = NULL + si = sigm2_setup (im, pm, px1,px2,nx,xblk, py1,py2,ny,yblk, order) + } + + # The device IM_MIN and IM_MAX parameters define the acceptable range + # of greyscale values for the output device (e.g., 0-255 for most 8-bit + # display devices). Values Z1 and Z2 are mapped linearly or + # logarithmically into IM_MIN and IM_MAX. + + dz1 = IM_MIN(ds) + dz2 = IM_MAX(ds) + if (fp_equalr (z1, z2)) { + z1 = z1 - 1 + z2 = z2 + 1 + } + + # If the user specifies the transfer function, verify that the + # intensity and greyscale are in range. + + if (zt == W_USER) { + call alims (Mems[uptr], U_MAXPTS, lut1, lut2) + dz1_s = short (dz1) + dz2_s = short (dz2) + if (lut2 < dz1_s || lut1 > dz2_s) + call eprintf ("User specified greyscales out of range\n") + if (z2 < IM_MIN(im) || z1 > IM_MAX(im)) + call eprintf ("User specified intensities out of range\n") + } + + # Type short pixels are treated as a special case to minimize vector + # operations for such images (which are common). If the image pixels + # are either short or real then only the ALTR (greyscale transformation) + # vector operation is required. The ALTR operator linearly maps + # greylevels in the range Z1:Z2 to DZ1:DZ2, and does a floor ceiling + # of DZ1:DZ2 on all pixels outside the range. 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. + + if (zt == W_UNITARY) { + unitary_greyscale_transformation = true + } else if (zt == W_LINEAR) { + unitary_greyscale_transformation = + (fp_equalr(z1,dz1) && fp_equalr(z2,dz2)) + } else + unitary_greyscale_transformation = false + + if (IM_PIXTYPE(im) == TY_SHORT && zt != W_LOG) { + z1_s = z1; z2_s = z2 + if (z1_s == z2_s) { + z1_s = z1_s - 1 + z2_s = z2_s + 1 + } + + for (wy=wy1; wy <= wy2; wy=wy+1) { + in = sigm2s (si, wy - wy1 + 1) + out = imps2s (ds, wx1, wx2, wy, wy) + + if (unitary_greyscale_transformation) { + call amovs (Mems[in], Mems[out], nx) + } else if (zt == W_USER) { + dz1_s = U_Z1; dz2_s = U_Z2 + call amaps (Mems[in],Mems[out],nx, z1_s,z2_s, dz1_s,dz2_s) + call aluts (Mems[out], Mems[out], nx, Mems[uptr]) + } else { + dz1_s = dz1; dz2_s = dz2 + call amaps (Mems[in],Mems[out],nx, z1_s,z2_s, dz1_s,dz2_s) + } + + if (si_ovrly != NULL) { + in = sigm2i (si_ovrly, wy - wy1 + 1) + call maskexprn (ocolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (ocolors, Memi[in+i]) + if (color >= 0) + Mems[out+i] = color + } + } + } + if (si_bpovrly != NULL) { + in = sigm2i (si_bpovrly, wy - wy1 + 1) + call maskexprn (bpcolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (bpcolors, Memi[in+i]) + if (color >= 0) + Mems[out+i] = color + } + } + } + } + + } else if (zt == W_USER) { + call salloc (rtemp, nx, TY_REAL) + + for (wy=wy1; wy <= wy2; wy=wy+1) { + in = sigm2r (si, wy - wy1 + 1) + out = imps2s (ds, wx1, wx2, wy, wy) + + call amapr (Memr[in], Memr[rtemp], nx, z1, z2, + real(U_Z1), real(U_Z2)) + call achtrs (Memr[rtemp], Mems[out], nx) + call aluts (Mems[out], Mems[out], nx, Mems[uptr]) + + if (si_ovrly != NULL) { + in = sigm2i (si_ovrly, wy - wy1 + 1) + call maskexprn (ocolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (ocolors, Memi[in+i]) + if (color >= 0) + Mems[out+i] = color + } + } + } + if (si_bpovrly != NULL) { + in = sigm2i (si_bpovrly, wy - wy1 + 1) + call maskexprn (bpcolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (bpcolors, Memi[in+i]) + if (color >= 0) + Mems[out+i] = color + } + } + } + } + + } else { + for (wy=wy1; wy <= wy2; wy=wy+1) { + in = sigm2r (si, wy - wy1 + 1) + out = imps2r (ds, wx1, wx2, wy, wy) + + if (unitary_greyscale_transformation) { + call amovr (Memr[in], Memr[out], nx) + } else if (zt == W_LOG) { + call amapr (Memr[in], Memr[out], nx, + z1, z2, 1.0, 10.0 ** MAXLOG) + do i = 0, nx-1 + Memr[out+i] = log10 (Memr[out+i]) + call amapr (Memr[out], Memr[out], nx, + 0.0, real(MAXLOG), dz1, dz2) + } else + call amapr (Memr[in], Memr[out], nx, z1, z2, dz1, dz2) + + if (si_ovrly != NULL) { + in = sigm2i (si_ovrly, wy - wy1 + 1) + call maskexprn (ocolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (ocolors, Memi[in+i]) + if (color >= 0) + Memr[out+i] = color + } + } + } + if (si_bpovrly != NULL) { + in = sigm2i (si_bpovrly, wy - wy1 + 1) + call maskexprn (bpcolors, in, nx) + do i = 0, nx-1 { + if (Memi[in+i] != 0) { + color = maskcolor (bpcolors, Memi[in+i]) + if (color >= 0) + Memr[out+i] = color + } + } + } + } + } + + call sigm2_free (si) + if (si_ovrly != NULL) + call sigm2_free (si_ovrly) + if (si_bpovrly != NULL) + call sigm2_free (si_bpovrly) + if (ovrly != NULL) + call imunmap (ovrly) + if (bpm != NULL) + call imunmap (bpm) +end + + +# DS_ERASE_BORDER -- Zero the border of the window if the frame has not been +# erased, and if the displayed section does not occupy the full window. +# It would be more efficient to do this while writing the greyscale data to +# the output image, but that would complicate the display procedures and frames +# are commonly erased before displaying an image. + +procedure ds_erase_border (im, ds, wdes) + +pointer im # input image +pointer ds # output image (display) +pointer wdes # window descriptor + +int wx1,wx2,wy1,wy2 # section of display window filled by image data +int dx1,dx2,dy1,dy2 # coords of full display window in device pixels +int i, nx +pointer wdwin, wdpix +pointer imps2s() +errchk imps2s + +begin + wdwin = W_WC(wdes,W_DWIN) + wdpix = W_WC(wdes,W_DPIX) + + # Set display pixels and display window pixels. + wx1 = nint (W_XS(wdpix)) + wx2 = nint (W_XE(wdpix)) + wy1 = nint (W_YS(wdpix)) + wy2 = nint (W_YE(wdpix)) + dx1 = max (1, nint (W_XS(wdwin))) + dx2 = min (IM_LEN(ds,1), nint (W_XE(wdwin) - 0.01)) + dy1 = max (1, nint (W_YS(wdwin))) + dy2 = min (IM_LEN(ds,2), nint (W_YE(wdwin) - 0.01)) + nx = dx2 - dx1 + 1 + + # Erase lower margin. + for (i=dy1; i < wy1; i=i+1) + call aclrs (Mems[imps2s (ds, dx1, dx2, i, i)], nx) + + # Erase left and right margins. By doing the right margin of a line + # immediately after the left margin we have a high liklihood that the + # display line will still be in the FIO buffer. + + for (i=wy1; i <= wy2; i=i+1) { + if (dx1 < wx1) + call aclrs (Mems[imps2s (ds, dx1, wx1-1, i, i)], wx1 - dx1) + if (wx2 < dx2) + call aclrs (Mems[imps2s (ds, wx2+1, dx2, i, i)], dx2 - wx2) + } + + # Erase upper margin. + for (i=wy2+1; i <= dy2; i=i+1) + call aclrs (Mems[imps2s (ds, dx1, dx2, i, i)], nx) +end diff --git a/pkg/images/tv/display/zardim.x b/pkg/images/tv/display/zardim.x new file mode 100644 index 00000000..e09c4b10 --- /dev/null +++ b/pkg/images/tv/display/zardim.x @@ -0,0 +1,21 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZARDIM -- Read data from a binary file display device. + +procedure zardim (chan, buf, nbytes, offset) + +int chan[ARB] +short buf[ARB] +int nbytes +long offset +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisrd (chan, buf, nbytes, offset) + } +end diff --git a/pkg/images/tv/display/zawrim.x b/pkg/images/tv/display/zawrim.x new file mode 100644 index 00000000..a7219b07 --- /dev/null +++ b/pkg/images/tv/display/zawrim.x @@ -0,0 +1,21 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZAWRIM -- Write data to a binary file display device. + +procedure zawrim (chan, buf, nbytes, offset) + +int chan[ARB] +short buf[ARB] +int nbytes +long offset +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iiswr (chan, buf, nbytes, offset) + } +end diff --git a/pkg/images/tv/display/zawtim.x b/pkg/images/tv/display/zawtim.x new file mode 100644 index 00000000..13756adc --- /dev/null +++ b/pkg/images/tv/display/zawtim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZAWTIM -- Wait for an image display frame which is addressable as +# a binary file. + +procedure zawtim (chan, nbytes) + +int chan[ARB], nbytes +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iiswt (chan, nbytes) + } +end diff --git a/pkg/images/tv/display/zblkim.x b/pkg/images/tv/display/zblkim.x new file mode 100644 index 00000000..55041809 --- /dev/null +++ b/pkg/images/tv/display/zblkim.x @@ -0,0 +1,23 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZBLKIM -- Blink binary file display device (millisecond time resolution). + +procedure zblkim (chan1, chan2, chan3, chan4, nframes, rate) + +int chan1[ARB] +int chan2[ARB] +int chan3[ARB] +int chan4[ARB] +int nframes +real rate +int device + +begin + device = chan1[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisblk (chan1, chan2, chan3, chan4, nframes, rate) + } +end diff --git a/pkg/images/tv/display/zclrim.x b/pkg/images/tv/display/zclrim.x new file mode 100644 index 00000000..268123cc --- /dev/null +++ b/pkg/images/tv/display/zclrim.x @@ -0,0 +1,18 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZCLRIM -- Color window binary file display device. + +procedure zclrim (chan) + +int chan[ARB] +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisclr (chan) + } +end diff --git a/pkg/images/tv/display/zclsim.x b/pkg/images/tv/display/zclsim.x new file mode 100644 index 00000000..8f3f34b0 --- /dev/null +++ b/pkg/images/tv/display/zclsim.x @@ -0,0 +1,22 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZCLSIM -- Close an image display frame which is addressable as +# a binary file. + +procedure zclsim (chan, status) + +int chan[ARB] +int status +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iiscls (chan, status) + default: + status = ERR + } +end diff --git a/pkg/images/tv/display/zdisplay.h b/pkg/images/tv/display/zdisplay.h new file mode 100644 index 00000000..b55b94dc --- /dev/null +++ b/pkg/images/tv/display/zdisplay.h @@ -0,0 +1,6 @@ +# Display devices defined by OS + +define IIS "/dev/iis" # IIS display device +define IIS_CHAN 1 # Device channel identifier +define DEVCODE 100 # Channel = DEVCODE * DEVCHAN +define FRTOCHAN (IIS_CHAN*DEVCODE+($1)) diff --git a/pkg/images/tv/display/zersim.x b/pkg/images/tv/display/zersim.x new file mode 100644 index 00000000..c1b280e4 --- /dev/null +++ b/pkg/images/tv/display/zersim.x @@ -0,0 +1,18 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZERSIM -- Erase binary file display device. + +procedure zersim (chan) + +int chan[ARB] +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisers (chan) + } +end diff --git a/pkg/images/tv/display/zfrmim.x b/pkg/images/tv/display/zfrmim.x new file mode 100644 index 00000000..de2bfee2 --- /dev/null +++ b/pkg/images/tv/display/zfrmim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZFRMIM -- Set FRAME display. + +procedure zfrmim (chan) + +int chan[ARB] + +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisrgb (chan, chan, chan) + } +end diff --git a/pkg/images/tv/display/zmapim.x b/pkg/images/tv/display/zmapim.x new file mode 100644 index 00000000..5c3e663a --- /dev/null +++ b/pkg/images/tv/display/zmapim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZMAPIM -- Set display map. + +procedure zmapim (chan, maptype) + +int chan[ARB] +char maptype[ARB] +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisofm (maptype) + } +end diff --git a/pkg/images/tv/display/zmtcim.x b/pkg/images/tv/display/zmtcim.x new file mode 100644 index 00000000..11dddb65 --- /dev/null +++ b/pkg/images/tv/display/zmtcim.x @@ -0,0 +1,18 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZMTCIM -- Match lut to frame. + +procedure zmtcim (chan1, chan2) + +int chan1[ARB], chan2[ARB] +int device + +begin + device = chan1[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iismtc (chan1, chan2) + } +end diff --git a/pkg/images/tv/display/zopnim.x b/pkg/images/tv/display/zopnim.x new file mode 100644 index 00000000..ddd18d3a --- /dev/null +++ b/pkg/images/tv/display/zopnim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZOPNIM -- Open an image display frame which is addressable as +# a binary file. + +procedure zopnim (devinfo, mode, chan) + +char devinfo[ARB] # packed devinfo string +int mode # access mode +int chan + +int iischan[2] # Kludge + +begin + call iisopn (devinfo, mode, iischan) + chan = iischan[1] +end diff --git a/pkg/images/tv/display/zrcrim.x b/pkg/images/tv/display/zrcrim.x new file mode 100644 index 00000000..3f4f939b --- /dev/null +++ b/pkg/images/tv/display/zrcrim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZRCRIM -- Read Cursor from binary file display device. + +procedure zrcrim (chan, xcur, ycur) + +int chan[ARB] +int status, xcur, ycur +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisrcr (status, xcur, ycur) + } +end diff --git a/pkg/images/tv/display/zrgbim.x b/pkg/images/tv/display/zrgbim.x new file mode 100644 index 00000000..04c0e147 --- /dev/null +++ b/pkg/images/tv/display/zrgbim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZRGBIM -- Set RGB display. + +procedure zrgbim (red_chan, green_chan, blue_chan) + +int red_chan[ARB], green_chan[ARB], blue_chan[ARB] + +int device + +begin + device = red_chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisrgb (red_chan, green_chan, blue_chan) + } +end diff --git a/pkg/images/tv/display/zrmim.x b/pkg/images/tv/display/zrmim.x new file mode 100644 index 00000000..f26ee6ef --- /dev/null +++ b/pkg/images/tv/display/zrmim.x @@ -0,0 +1,19 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZRMIM -- Zoom and roam display. + +procedure zrmim (chan, zfactor) + +int chan[ARB] +int zfactor +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iisrm (zfactor) + } +end diff --git a/pkg/images/tv/display/zscale.x b/pkg/images/tv/display/zscale.x new file mode 100644 index 00000000..abbf2ecb --- /dev/null +++ b/pkg/images/tv/display/zscale.x @@ -0,0 +1,623 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <ctype.h> +include <imhdr.h> +include <imset.h> +include <pmset.h> +include <imio.h> + +# User callable routines. +# ZSCALE -- Sample an image and compute greyscale limits. +# MZSCALE -- Sample an image with pixel masks and compute greyscale limits. +# ZSC_PMSECTION -- Create a pixel mask from an image section. +# ZSC_ZLIMITS -- Compute Z transform limits from a sample of pixels. + + +# ZSCALE -- Sample an image and compute greyscale limits. +# A sample mask is created based on the input parameters and then +# MZSCALE is called. + +procedure zscale (im, z1, z2, contrast, optimal_sample_size, len_stdline) + +pointer im # image to be sampled +real z1, z2 # output min and max greyscale values +real contrast # adj. to slope of transfer function +int optimal_sample_size # desired number of pixels in sample +int len_stdline # optimal number of pixels per line + +int nc, nl +pointer sp, section, zpm, zsc_pmsection() +errchk zsc_pmsection, mzscale + +begin + call smark (sp) + call salloc (section, SZ_FNAME, TY_CHAR) + + # Make the sample image section. + switch (IM_NDIM(im)) { + case 1: + call sprintf (Memc[section], SZ_FNAME, "[*]") + default: + nc = max (1, min (IM_LEN(im,1), len_stdline)) + nl = max (1, min (IM_LEN(im,2), optimal_sample_size / nc)) + call sprintf (Memc[section], SZ_FNAME, "[*:%d,*:%d]") + call pargi (IM_LEN(im,1) / nc) + call pargi (IM_LEN(im,2) / nl) + } + + # Make a mask and compute the greyscale limits. + zpm = zsc_pmsection (Memc[section], im) + call mzscale (im, zpm, NULL, contrast, optimal_sample_size, z1, z2) + call imunmap (zpm) + call sfree (sp) +end + + +# MZSCALE -- Sample an image with pixel masks and compute greyscale limits. +# The image is sampled through a pixel mask. If no pixel mask is given +# a uniform sample mask is generated. If a bad pixel mask is given +# bad pixels in the sample are eliminated. Once the sample is obtained +# the greyscale limits are obtained using the ZSC_ZLIMITS algorithm. + +procedure mzscale (im, zpm, bpm, contrast, maxpix, z1, z2) + +pointer im #I image to be sampled +pointer zpm #I pixel mask for sampling +pointer bpm #I bad pixel mask +real contrast #I contrast parameter +int maxpix #I maximum number of pixels in sample +real z1, z2 #O output min and max greyscale values + +int i, ndim, nc, nl, npix, nbp, imstati() +pointer sp, section, v, sample, zmask, bp, zim, pmz, pmb, buf +pointer zsc_pmsection(), imgnlr() +bool pm_linenotempty() +errchk zsc_pmsection, zsc_zlimits + +begin + call smark (sp) + call salloc (section, SZ_FNAME, TY_CHAR) + call salloc (v, IM_MAXDIM, TY_LONG) + call salloc (sample, maxpix, TY_REAL) + zmask = NULL + bp = NULL + + ndim = min (2, IM_NDIM(im)) + nc = IM_LEN(im,1) + nl = IM_LEN(im,2) + + # Generate a uniform sample mask if none is given. + if (zpm == NULL) { + switch (IM_NDIM(im)) { + case 1: + call sprintf (Memc[section], SZ_FNAME, "[*]") + default: + i = max (1., sqrt ((nc-1)*(nl-1) / real (maxpix))) + call sprintf (Memc[section], SZ_FNAME, "[*:%d,*:%d]") + call pargi (i) + call pargi (i) + } + zim = zsc_pmsection (Memc[section], im) + pmz = imstati (zim, IM_PMDES) + } else + pmz = imstati (zpm, IM_PMDES) + + # Set bad pixel mask. + if (bpm != NULL) + pmb = imstati (bpm, IM_PMDES) + else + pmb = NULL + + # Get the sample up to maxpix pixels. + npix = 0 + nbp = 0 + call amovkl (long(1), Memi[v], IM_MAXDIM) + repeat { + if (pm_linenotempty (pmz, Meml[v])) { + if (zmask == NULL) + call salloc (zmask, nc, TY_INT) + call pmglpi (pmz, Meml[v], Memi[zmask], 0, nc, 0) + if (pmb != NULL) { + if (pm_linenotempty (pmb, Meml[v])) { + if (bp == NULL) + call salloc (bp, nc, TY_INT) + call pmglpi (pmb, Meml[v], Memi[bp], 0, nc, 0) + nbp = nc + } else + nbp = 0 + + } + if (imgnlr (im, buf, Meml[v]) == EOF) + break + do i = 0, nc-1 { + if (Memi[zmask+i] == 0) + next + if (nbp > 0) + if (Memi[bp+i] != 0) + next + Memr[sample+npix] = Memr[buf+i] + npix = npix + 1 + if (npix == maxpix) + break + } + if (npix == maxpix) + break + } else { + do i = 2, ndim { + Meml[v+i-1] = Meml[v+i-1] + 1 + if (Meml[v+i-1] <= IM_LEN(im,i)) + break + else if (i < ndim) + Meml[v+i-1] = 1 + } + } + } until (Meml[v+ndim-1] > IM_LEN(im,ndim)) + + if (zpm == NULL) + call imunmap (zim) + + # Compute greyscale limits. + call zsc_zlimits (Memr[sample], npix, contrast, z1, z2) + + call sfree (sp) +end + + +# ZSC_PMSECTION -- Create a pixel mask from an image section. +# This only applies the mask to the first plane of the image. + +pointer procedure zsc_pmsection (section, refim) + +char section[ARB] #I Image section +pointer refim #I Reference image pointer + +int i, j, ip, ndim, temp, a[2], b[2], c[2], rop, ctoi() +pointer pm, im, mw, dummy, pm_newmask(), im_pmmapo(), imgl1i(), mw_openim() +define error_ 99 + +begin + # Decode the section string. + call amovki (1, a, 2) + call amovki (1, b, 2) + call amovki (1, c, 2) + ndim = min (2, IM_NDIM(refim)) + do i = 1, ndim + b[i] = IM_LEN(refim,i) + + ip = 1 + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (section[ip] == '[') { + ip = ip + 1 + + do i = 1, ndim { + while (IS_WHITE(section[ip])) + ip = ip + 1 + + # 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[i] = temp + if (section[ip] == ':') { + ip = ip + 1 + if (ctoi (section, ip, b[i]) == 0) # a:b + goto error_ + } else + b[i] = a[i] + } else if (section[ip] == '-') { # -* + temp = a[i] + a[i] = b[i] + b[i] = temp + ip = ip + 1 + if (section[ip] == '*') + ip = ip + 1 + } else if (section[ip] == '*') # * + ip = ip + 1 + if (section[ip] == ':') { # ..:step + ip = ip + 1 + if (ctoi (section, ip, c[i]) == 0) + goto error_ + else if (c[i] == 0) + goto error_ + } + if (a[i] > b[i] && c[i] > 0) + c[i] = -c[i] + + while (IS_WHITE(section[ip])) + ip = ip + 1 + if (i < ndim) { + if (section[ip] != ',') + goto error_ + } else { + if (section[ip] != ']') + goto error_ + } + ip = ip + 1 + } + } + + # In this case make the values be increasing only. + do i = 1, ndim + if (c[i] < 0) { + temp = a[i] + a[i] = b[i] + b[i] = temp + c[i] = -c[i] + } + + # Make the mask. + pm = pm_newmask (refim, 16) + + rop = PIX_SET+PIX_VALUE(1) + if (c[1] == 1 && c[2] == 1) + call pm_box (pm, a[1], a[2], b[1], b[2], rop) + + else if (c[1] == 1) + for (i=a[2]; i<=b[2]; i=i+c[2]) + call pm_box (pm, a[1], i, b[1], i, rop) + + else + for (i=a[2]; i<=b[2]; i=i+c[2]) + for (j=a[1]; j<=b[1]; j=j+c[1]) + call pm_point (pm, j, i, rop) + + i = IM_NPHYSDIM(refim) + IM_NPHYSDIM(refim) = ndim + im = im_pmmapo (pm, refim) + IM_NPHYSDIM(refim) = i + dummy = imgl1i (im) # Force I/O to set header + ifnoerr (mw = mw_openim (refim)) { # Set WCS + call mw_saveim (mw, im) + call mw_close (mw) + } + + return (im) + +error_ + call error (1, "Error in image section specification") +end + + +.help zsc_zlimits +.nf ___________________________________________________________________________ +ZSC_ZLIMITS -- Compute limits for a linear transform that best samples the +the histogram about the median value. This is often called to compute +greyscale limits from a sample of pixel values. + +If the number of pixels is too small an error condition is returned. If +the contrast parameter value is zero the limits of the sample are +returned. Otherwise the sample is sorted and the median is found from the +central value(s). A straight line is fitted to the sorted sample with +interative rejection. If more than half the pixels are rejected the full +range is returned. The contrast parameter is used to adjust the transfer +slope about the median. The final limits are the extension of the fitted +line to the first and last array index. +.endhelp ______________________________________________________________________ + +define MIN_NPIXELS 5 # smallest permissible sample +define MAX_REJECT 0.5 # max frac. of pixels to be rejected +define GOOD_PIXEL 0 # use pixel in fit +define BAD_PIXEL 1 # ignore pixel in all computations +define REJECT_PIXEL 2 # reject pixel after a bit +define KREJ 2.5 # k-sigma pixel rejection factor +define MAX_ITERATIONS 5 # maximum number of fitline iterations + + +# ZSC_ZLIMITS -- Compute Z transform limits from a sample of pixels. + +procedure zsc_zlimits (sample, npix, contrast, z1, z2) + +real sample[ARB] #I Sample of pixel values (possibly resorted) +int npix #I Number of pixels +real contrast #I Contrast algorithm parameter +real z1, z2 #O Z transform limits + +int center_pixel, minpix, ngoodpix, ngrow, zsc_fit_line() +real zmin, zmax, median +real zstart, zslope + +begin + # Check for a sufficient sample. + if (npix < MIN_NPIXELS) + call error (1, "Insufficient sample pixels found") + + # If contrast is zero return the range. + if (contrast == 0.) { + call alimr (sample, npix, z1, z2) + return + } + + # Sort the sample, compute the range, and median pixel values. + # The median value is the average of the two central values if there + # are an even number of pixels in the sample. + + call asrtr (sample, sample, npix) + zmin = sample[1] + zmax = sample[npix] + + center_pixel = (npix + 1) / 2 + if (mod (npix, 2) == 1) + median = sample[center_pixel] + else + median = (sample[center_pixel] + sample[center_pixel+1]) / 2 + + # Fit a line to the sorted sample vector. If more than half of the + # pixels in the sample are rejected give up and return the full range. + # If the user-supplied contrast factor is not 1.0 adjust the scale + # accordingly and compute Z1 and Z2, the y intercepts at indices 1 and + # npix. + + minpix = max (MIN_NPIXELS, int (npix * MAX_REJECT)) + ngrow = max (1, nint (npix * .01)) + ngoodpix = zsc_fit_line (sample, npix, zstart, zslope, + KREJ, ngrow, MAX_ITERATIONS) + + if (ngoodpix < minpix) { + z1 = zmin + z2 = zmax + } else { + if (contrast > 0) + zslope = zslope / contrast + z1 = max (zmin, median - (center_pixel - 1) * zslope) + z2 = min (zmax, median + (npix - center_pixel) * zslope) + } +end + + +# ZSC_FIT_LINE -- Fit a straight line to a data array of type real. This is +# an iterative fitting algorithm, wherein points further than ksigma from the +# current fit are excluded from the next fit. Convergence occurs when the +# next iteration does not decrease the number of pixels in the fit, or when +# there are no pixels left. The number of pixels left after pixel rejection +# is returned as the function value. + +int procedure zsc_fit_line (data, npix, zstart, zslope, krej, ngrow, maxiter) + +real data[npix] # data to be fitted +int npix # number of pixels before rejection +real zstart # Z-value of pixel data[1] (output) +real zslope # dz/pixel (output) +real krej # k-sigma pixel rejection factor +int ngrow # number of pixels of growing +int maxiter # max iterations + +int i, ngoodpix, last_ngoodpix, minpix, niter +real xscale, z0, dz, x, z, mean, sigma, threshold +double sumxsqr, sumxz, sumz, sumx, rowrat +pointer sp, flat, badpix, normx +int zsc_reject_pixels(), zsc_compute_sigma() + +begin + call smark (sp) + + if (npix <= 0) + return (0) + else if (npix == 1) { + zstart = data[1] + zslope = 0.0 + return (1) + } else + xscale = 2.0 / (npix - 1) + + # Allocate a buffer for data minus fitted curve, another for the + # normalized X values, and another to flag rejected pixels. + + call salloc (flat, npix, TY_REAL) + call salloc (normx, npix, TY_REAL) + call salloc (badpix, npix, TY_SHORT) + call aclrs (Mems[badpix], npix) + + # Compute normalized X vector. The data X values [1:npix] are + # normalized to the range [-1:1]. This diagonalizes the lsq matrix + # and reduces its condition number. + + do i = 0, npix - 1 + Memr[normx+i] = i * xscale - 1.0 + + # Fit a line with no pixel rejection. Accumulate the elements of the + # matrix and data vector. The matrix M is diagonal with + # M[1,1] = sum x**2 and M[2,2] = ngoodpix. The data vector is + # DV[1] = sum (data[i] * x[i]) and DV[2] = sum (data[i]). + + sumxsqr = 0 + sumxz = 0 + sumx = 0 + sumz = 0 + + do i = 1, npix { + x = Memr[normx+i-1] + z = data[i] + sumxsqr = sumxsqr + (x ** 2) + sumxz = sumxz + z * x + sumz = sumz + z + } + + # Solve for the coefficients of the fitted line. + z0 = sumz / npix + dz = sumxz / sumxsqr + + # Iterate, fitting a new line in each iteration. Compute the flattened + # data vector and the sigma of the flat vector. Compute the lower and + # upper k-sigma pixel rejection thresholds. Run down the flat array + # and detect pixels to be rejected from the fit. Reject pixels from + # the fit by subtracting their contributions from the matrix sums and + # marking the pixel as rejected. + + ngoodpix = npix + minpix = max (MIN_NPIXELS, int (npix * MAX_REJECT)) + + for (niter=1; niter <= maxiter; niter=niter+1) { + last_ngoodpix = ngoodpix + + # Subtract the fitted line from the data array. + call zsc_flatten_data (data, Memr[flat], Memr[normx], npix, z0, dz) + + # Compute the k-sigma rejection threshold. In principle this + # could be more efficiently computed using the matrix sums + # accumulated when the line was fitted, but there are problems with + # numerical stability with that approach. + + ngoodpix = zsc_compute_sigma (Memr[flat], Mems[badpix], npix, + mean, sigma) + threshold = sigma * krej + + # Detect and reject pixels further than ksigma from the fitted + # line. + ngoodpix = zsc_reject_pixels (data, Memr[flat], Memr[normx], + Mems[badpix], npix, sumxsqr, sumxz, sumx, sumz, threshold, + ngrow) + + # Solve for the coefficients of the fitted line. Note that after + # pixel rejection the sum of the X values need no longer be zero. + + if (ngoodpix > 0) { + rowrat = sumx / sumxsqr + z0 = (sumz - rowrat * sumxz) / (ngoodpix - rowrat * sumx) + dz = (sumxz - z0 * sumx) / sumxsqr + } + + if (ngoodpix >= last_ngoodpix || ngoodpix < minpix) + break + } + + # Transform the line coefficients back to the X range [1:npix]. + zstart = z0 - dz + zslope = dz * xscale + + call sfree (sp) + return (ngoodpix) +end + + +# ZSC_FLATTEN_DATA -- Compute and subtract the fitted line from the data array, +# returned the flattened data in FLAT. + +procedure zsc_flatten_data (data, flat, x, npix, z0, dz) + +real data[npix] # raw data array +real flat[npix] # flattened data (output) +real x[npix] # x value of each pixel +int npix # number of pixels +real z0, dz # z-intercept, dz/dx of fitted line +int i + +begin + do i = 1, npix + flat[i] = data[i] - (x[i] * dz + z0) +end + + +# ZSC_COMPUTE_SIGMA -- Compute the root mean square deviation from the +# mean of a flattened array. Ignore rejected pixels. + +int procedure zsc_compute_sigma (a, badpix, npix, mean, sigma) + +real a[npix] # flattened data array +short badpix[npix] # bad pixel flags (!= 0 if bad pixel) +int npix +real mean, sigma # (output) + +real pixval +int i, ngoodpix +double sum, sumsq, temp + +begin + sum = 0 + sumsq = 0 + ngoodpix = 0 + + # Accumulate sum and sum of squares. + do i = 1, npix + if (badpix[i] == GOOD_PIXEL) { + pixval = a[i] + ngoodpix = ngoodpix + 1 + sum = sum + pixval + sumsq = sumsq + pixval ** 2 + } + + # Compute mean and sigma. + switch (ngoodpix) { + case 0: + mean = INDEF + sigma = INDEF + case 1: + mean = sum + sigma = INDEF + default: + mean = sum / ngoodpix + temp = sumsq / (ngoodpix - 1) - sum**2 / (ngoodpix * (ngoodpix - 1)) + if (temp < 0) # possible with roundoff error + sigma = 0.0 + else + sigma = sqrt (temp) + } + + return (ngoodpix) +end + + +# ZSC_REJECT_PIXELS -- Detect and reject pixels more than "threshold" greyscale +# units from the fitted line. The residuals about the fitted line are given +# by the "flat" array, while the raw data is in "data". Each time a pixel +# is rejected subtract its contributions from the matrix sums and flag the +# pixel as rejected. When a pixel is rejected reject its neighbors out to +# a specified radius as well. This speeds up convergence considerably and +# produces a more stringent rejection criteria which takes advantage of the +# fact that bad pixels tend to be clumped. The number of pixels left in the +# fit is returned as the function value. + +int procedure zsc_reject_pixels (data, flat, normx, badpix, npix, + sumxsqr, sumxz, sumx, sumz, threshold, ngrow) + +real data[npix] # raw data array +real flat[npix] # flattened data array +real normx[npix] # normalized x values of pixels +short badpix[npix] # bad pixel flags (!= 0 if bad pixel) +int npix +double sumxsqr,sumxz,sumx,sumz # matrix sums +real threshold # threshold for pixel rejection +int ngrow # number of pixels of growing + +int ngoodpix, i, j +real residual, lcut, hcut +double x, z + +begin + ngoodpix = npix + lcut = -threshold + hcut = threshold + + do i = 1, npix + if (badpix[i] == BAD_PIXEL) + ngoodpix = ngoodpix - 1 + else { + residual = flat[i] + if (residual < lcut || residual > hcut) { + # Reject the pixel and its neighbors out to the growing + # radius. We must be careful how we do this to avoid + # directional effects. Do not turn off thresholding on + # pixels in the forward direction; mark them for rejection + # but do not reject until they have been thresholded. + # If this is not done growing will not be symmetric. + + do j = max(1,i-ngrow), min(npix,i+ngrow) { + if (badpix[j] != BAD_PIXEL) { + if (j <= i) { + x = normx[j] + z = data[j] + sumxsqr = sumxsqr - (x ** 2) + sumxz = sumxz - z * x + sumx = sumx - x + sumz = sumz - z + badpix[j] = BAD_PIXEL + ngoodpix = ngoodpix - 1 + } else + badpix[j] = REJECT_PIXEL + } + } + } + } + + return (ngoodpix) +end diff --git a/pkg/images/tv/display/zsttim.x b/pkg/images/tv/display/zsttim.x new file mode 100644 index 00000000..dc6c91f6 --- /dev/null +++ b/pkg/images/tv/display/zsttim.x @@ -0,0 +1,26 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <knet.h> +include <fio.h> +include "iis.h" + +# ZSTTIM -- Return status on binary file display device. + +procedure zsttim (chan, what, lvalue) + +int chan[ARB], what +long lvalue + +include "iis.com" + +begin + call zsttgd (iischan, what, lvalue) + + if (what == FSTT_MAXBUFSIZE) { + # Return the maximum transfer size in bytes. + if (lvalue == 0) + lvalue = FSTT_MAXBUFSIZE + if (!packit) + lvalue = min (IIS_MAXBUFSIZE, lvalue) * 2 + } +end diff --git a/pkg/images/tv/display/zwndim.x b/pkg/images/tv/display/zwndim.x new file mode 100644 index 00000000..d27027cf --- /dev/null +++ b/pkg/images/tv/display/zwndim.x @@ -0,0 +1,31 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include "zdisplay.h" + +# ZWNDIM -- Window binary file display device. + +procedure zwndim (chan) + +int chan[ARB] +int device + +begin + device = chan[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iiswnd3 (chan, chan, chan) + } +end + +procedure zwndim3 (chan1, chan2, chan3) + +int chan1[ARB], chan2[ARB], chan3[ARB] +int device + +begin + device = chan1[1] / DEVCODE + switch (device) { + case IIS_CHAN: + call iiswnd3 (chan1, chan2, chan3) + } +end diff --git a/pkg/images/tv/display/zzdebug.x b/pkg/images/tv/display/zzdebug.x new file mode 100644 index 00000000..eb642d42 --- /dev/null +++ b/pkg/images/tv/display/zzdebug.x @@ -0,0 +1,165 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <imhdr.h> + +task mktest = t_mktest, + sigl2 = t_sigl2, + wrimage = t_wrimage, + zscale = t_zscale, + rcur = t_rcur + +define TWOPI 6.23 + + +# MKTEST -- Make a test image containing a circularly symetric sinusoid. + +procedure t_mktest() + +char imname[SZ_FNAME] +int nx, ny +int i, j +real period, xcen, ycen, radius +pointer im, line + +int clgeti() +real clgetr() +pointer immap(), impl2r() + +begin + call clgstr ("imname", imname, SZ_FNAME) + im = immap (imname, NEW_IMAGE, 0) + + nx = clgeti ("nx") + ny = clgeti ("ny") + period = clgetr ("period") + + IM_LEN(im,1) = nx + IM_LEN(im,2) = ny + + xcen = (nx + 1) / 2.0 + ycen = (ny + 1) / 2.0 + + do j = 1, ny { + line = impl2r (im, j) + do i = 1, nx { + radius = sqrt ((i - xcen) ** 2 + (j - ycen) ** 2) + Memr[line+i-1] = sin ((radius / period) * TWOPI) * 255.0 + } + } + + call imunmap (im) +end + + +# READ -- Benchmark scaled input procedure. + +procedure t_sigl2 () + +char imname[SZ_FNAME] +pointer im, si, buf +int i, nx, ny, xblk, yblk +pointer sigl2_setup(), sigl2s(), immap() + +begin + call clgstr ("imname", imname, SZ_FNAME) + im = immap (imname, READ_ONLY, 0) + + nx = IM_LEN(im,1) + ny = IM_LEN(im,2) + + xblk = INDEFI + yblk = INDEFI + si = sigl2_setup (im, 1.0,real(nx),nx,xblk, 1.0,real(ny),ny,yblk,0) + + do i = 1, ny + buf = sigl2s (si, i) + + call sigl2_free (si) + call imunmap (im) +end + + +# WRIMAGE -- Benchmark image output as used in the display program. + +procedure t_wrimage () + +char imname[SZ_FNAME] +int i, ncols, nlines +pointer im, buf +int clgeti() +pointer immap(), imps2s() + +begin + call clgstr ("imname", imname, SZ_FNAME) + im = immap (imname, NEW_IMAGE, 0) + + ncols = clgeti ("ncols") + nlines = clgeti ("nlines") + + IM_LEN(im,1) = ncols + IM_LEN(im,2) = nlines + IM_PIXTYPE(im) = TY_SHORT + + do i = 1, nlines + buf = imps2s (im, 1, ncols, i, i) + + call imunmap (im) +end + + +# ZSCALE -- Test the zscale procedure, used to determine the smallest range of +# greyscale values which preserves the most information in an image. + +procedure t_zscale() + +char imname[SZ_FNAME] +int sample_size, len_stdline +real z1, z2, contrast +int clgeti() +real clgetr() +pointer im, immap() + +begin + call clgstr ("imname", imname, SZ_FNAME) + im = immap (imname, READ_ONLY, 0) + + sample_size = clgeti ("npix") + len_stdline = clgeti ("stdline") + contrast = clgetr ("contrast") + + call zscale (im, z1, z2, contrast, sample_size, len_stdline) + call printf ("z1=%g, z2=%g\n") + call pargr (z1) + call pargr (z2) +end + + +# RCUR -- Try reading the image cursor. + +procedure t_rcur() + +real x, y +int wcs, key +int wci, pause +char device[SZ_FNAME] +char strval[SZ_LINE] + +bool clgetb() +int btoi(), clgeti(), imdrcur() + +begin + call clgstr ("device", device, SZ_FNAME) + wci = clgeti ("wcs") + pause = btoi (clgetb ("pause")) + + while (imdrcur (device, x,y,wcs,key,strval,SZ_LINE, wci,pause) != EOF) { + call printf ("%8.2f %8.2f %d %o %s\n") + call pargr (x) + call pargr (y) + call pargi (wcs) + call pargi (key) + call pargstr (strval) + if (key == 'q') + break + } +end |