Initial commit

1 files changed, 305 insertions, 0 deletions

diff --git a/pkg/plot/t_hafton.x b/pkg/plot/t_hafton.x
new file mode 100644
index 00000000..07e70393
--- /dev/null
+++ b/pkg/plot/t_hafton.x
@@ -0,0 +1,305 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<error.h>
+include	<gset.h>
+include	<config.h>
+include	<mach.h>
+include	<imhdr.h>
+include	<xwhen.h>
+include	<fset.h>
+
+define	DUMMY	6
+define	SAMPLE_SIZE	1000
+define	LEN_STDLINE	40
+
+# HAFTON -- Draw a half tone plot of an image section.  This is an
+# interface to the NCAR HAFTON routine.  
+
+procedure t_hafton()
+
+int	sign
+bool	sub, pre
+pointer	im, subras, gp
+int	tcojmp[LEN_JUMPBUF]
+char	imsect[SZ_FNAME], mapping_function[SZ_FNAME]
+char	device[SZ_FNAME], title[SZ_LINE], system_id[SZ_LINE]
+int	ncols, nlines, epa, status, wkid, mode, old_onint
+int	nlevels, nprm, nopt, xres, yres, nfunction, nx, ny
+real	z1, z2, wx1, wx2, wy1, wy2, contrast
+real	xs, xe, ys, ye, vx1, vx2, vy1, vy2
+
+real	clgetr()
+extern	hf_tco_onint()
+int	clgeti(),  strncmp()
+pointer	gopen(), plt_getdata(), immap()
+bool	clgetb(), fp_equalr(), streq()
+common	/tcocom/ tcojmp
+
+begin
+	# Get image section string and output device.
+	call clgstr ("image", imsect, SZ_FNAME)
+	call clgstr ("device", device, SZ_FNAME)
+
+	# Map image.
+	im = immap (imsect, READ_ONLY, 0)
+
+	z1 = clgetr ("z1")
+	z2 = clgetr ("z2")
+
+	# Gaurantee that image min/max is up to date
+	if (IM_LIMTIME(im) < IM_MTIME(im))
+	    call hf_minmax (im, IM_MIN(im), IM_MAX(im))
+
+	if (fp_equalr (z1, z2)) {
+	    z1 = IM_MIN(im)
+	    z2 = IM_MAX(im)
+	}
+
+	# User can specify the type of mapping function used, and whether
+	# the contrast is negative or positive.  
+
+	nlevels = clgeti ("nlevels")
+	contrast = clgetr ("contrast")
+	call clgstr ("mapping_function", mapping_function, SZ_FNAME)
+
+	# Assign integer code to specified mapping function
+	if (strncmp (mapping_function, "linear", 2) == 0)
+	    nfunction = 1
+	else if (strncmp (mapping_function, "exponential", 1) == 0)
+	    nfunction = 2
+	else if (strncmp (mapping_function, "logarithmic", 2) == 0)
+	    nfunction = 3
+	else if (strncmp (mapping_function, "sinusoidal",  1) == 0)
+	    nfunction = 4
+	else if (strncmp (mapping_function, "arcsine",     1) == 0)
+	    nfunction = 5
+	else if (strncmp (mapping_function, "crtpict",     1) == 0)
+	    nfunction = 6
+	else
+	    call error (0, "Hafton: unknown mapping function")
+
+	sign = 1.0
+	if (contrast < 0.0)
+	    sign = -1.0
+	nopt = sign * nfunction
+
+	mode = NEW_FILE
+	if (clgetb ("append"))
+	    mode = APPEND
+
+	# Read in subraster.  Image resolution can be decreased by
+	# subsampling or block averaging.
+
+	xres = clgeti ("xres")
+	yres = clgeti ("yres")
+	sub = clgetb ("subsample")
+	pre = clgetb ("preserve")
+
+	# Retrieve values from image header that will be needed.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	if (streq (title, "imtitle")) {
+	    call strcpy (imsect, title, SZ_LINE)
+	    call strcat (": ", title, SZ_LINE)
+	    call strcat (IM_TITLE(im), title, SZ_LINE)
+	}
+
+	xs = 1.0
+	xe = real (ncols)
+	ys = 1.0
+	ye = real (nlines)
+
+	# Get data with proper resolution.  Procedure plt_getdata returns
+	# a pointer to the data matrix to be contoured.  The resolution
+	# is decreased by the specified mathod in this procedure.  The
+	# dimensions of the data array are also returned.  The image
+	# header pointer can be unmapped after plt_getdata is called.
+
+	nx = 0
+	ny = 0
+	subras = plt_getdata (im, sub, pre, xres, yres, nx, ny)
+
+	if (nfunction == 6) {
+	    # User wants crtpict automatic algorithm - linear mapping
+	    # between calculated z1, z2 using possible non-integer contrast.
+	    # Get z1, z2 as if positive contrast.  Set nopt later to negative
+	    # if necessary.
+
+	    call zscale (im, z1, z2, abs(contrast), SAMPLE_SIZE, LEN_STDLINE)
+	}
+
+	call eprintf ("Intensities from z1=%.2f to z2=%.2f mapped with a")
+    	    call pargr (z1)
+    	    call pargr (z2)
+
+        switch (nfunction) {
+        case (1):
+	    call eprintf (" linear function\n")
+        case (2):
+	    call eprintf ("n exponential function\n")
+        case (3):
+	    call eprintf (" logarithmic function\n")
+        case (4):
+	    call eprintf (" sinusodial function\n")
+        case (5):
+	    call eprintf ("n arcsine function\n")
+        case (6):
+	    call eprintf (" CRTPICT function\n")
+	    if (nopt > 0) {
+		# Positive contrast.  Set nopt to positive linear mapping.
+	        nopt = 1
+	    } else  {
+		# Negative contrast. Set nopt to negative linear mapping.
+		nopt = -1
+	    }
+        }
+
+	vx1 = clgetr ("vx1")
+	vx2 = clgetr ("vx2")
+	vy1 = clgetr ("vy1")
+	vy2 = clgetr ("vy2")
+
+	# Open device and make contour plot.
+	call gopks (STDERR)
+	wkid = 1
+	gp = gopen (device, mode, STDGRAPH)
+	call gopwk (wkid, DUMMY, gp)
+	call gacwk (wkid)
+
+	call pl_map_viewport (gp,
+	    ncols, nlines, vx1, vx2, vy1, vy2, clgetb ("fill"), true)
+	nprm = -1
+
+	# Install interrupt exception handler.
+	call zlocpr (hf_tco_onint, epa)
+	call xwhen (X_INT, epa, old_onint)
+
+	# Make the hafton plot.  If an interrupt occurs ZSVJMP is reentered
+	# with an error status.
+
+	call zsvjmp (tcojmp, status)
+	if (status == OK) {
+	    call hafton (Memr[subras], nx, nx, ny, z1, z2,
+	        nlevels, nopt, nprm, 0, 0.)
+	} else {
+	    call gcancel (gp)
+	    call fseti (STDOUT, F_CANCEL, OK)
+	}
+
+	# Should a fancy (crtpict like) perimeter be drawn around the plot?
+	if (clgetb ("perimeter")) {
+	    call gswind (gp, xs, xe, ys, ye)
+	    call draw_perimeter (gp)
+	} else
+	    call perim (1, ncols - 1, nlines - 1, 1)
+
+	# Now find window and output text string title.  The window is
+	# set to the full image coordinates for labelling.
+
+	call ggview (gp, wx1, wx2, wy1, wy2)
+	call gseti (gp, G_WCS, 0)
+	call gtext (gp, (wx1 + wx2) / 2.0, wy2 + .03, title, "h=c;v=b;f=b;s=.7")
+
+	# Add system id banner to plot.
+	call gseti (gp, G_CLIP, NO)
+	call sysid (system_id, SZ_LINE)
+	call gtext (gp, (wx1+wx2)/2.0, wy1-0.07, system_id, "h=c;v=b;s=.5")
+
+	call gdawk (wkid)
+	call gclwk (wkid)
+	call gclks ()
+	call imunmap (im)
+
+	# Free space used for scaled input routines.
+	call mfree (subras, TY_REAL)
+end
+
+
+# HF_TCO_ONINT -- Interrupt handler for the task hafton.  Branches back to 
+# ZSVJMP in the main routine to permit shutdown without an error message.
+
+procedure hf_tco_onint (vex, next_handler)
+
+int	vex		# virtual exception
+int	next_handler	# not used
+
+int	tcojmp[LEN_JUMPBUF]
+common	/tcocom/ tcojmp
+
+begin
+	call xer_reset()
+	call zdojmp (tcojmp, vex)
+end
+
+
+# HF_MINMAX -- Compute the minimum and maximum pixel values of an image.
+# Works for images of any dimensionality, size, or datatype, although
+# the min and max values can currently only be stored in the image header
+# as real values.
+
+procedure hf_minmax (im, min_value, max_value)
+
+pointer	im				# image descriptor
+real	min_value			# minimum pixel value in image (out)
+real	max_value			# maximum pixel value in image (out)
+
+pointer	buf
+bool	first_line
+long	v[IM_MAXDIM]
+short	minval_s, maxval_s
+long	minval_l, maxval_l
+real	minval_r, maxval_r
+int	imgnls(), imgnll(), imgnlr()
+errchk	amovkl, imgnls, imgnll, imgnlr, alims, aliml, alimr
+
+begin
+	call amovkl (long(1), v, IM_MAXDIM)		# start vector
+	first_line = true
+	min_value = INDEF
+	max_value = INDEF
+
+	switch (IM_PIXTYPE(im)) {
+	case TY_SHORT:
+	    while (imgnls (im, buf, v) != EOF) {
+		call alims (Mems[buf], IM_LEN(im,1), minval_s, maxval_s)
+		if (first_line) {
+		    min_value = minval_s
+		    max_value = maxval_s
+		    first_line = false
+		} else {
+		    if (minval_s < min_value)
+			min_value = minval_s
+		    if (maxval_s > max_value)
+			max_value = maxval_s
+		}
+	    }
+	case TY_USHORT, TY_INT, TY_LONG:
+	    while (imgnll (im, buf, v) != EOF) {
+		call aliml (Meml[buf], IM_LEN(im,1), minval_l, maxval_l)
+		if (first_line) {
+		    min_value = minval_l
+		    max_value = maxval_l
+		    first_line = false
+		} else {
+		    if (minval_l < min_value)
+			min_value = minval_l
+		    if (maxval_l > max_value)
+			max_value = maxval_l
+		}
+	    }
+	default:
+	    while (imgnlr (im, buf, v) != EOF) {
+		call alimr (Memr[buf], IM_LEN(im,1), minval_r, maxval_r)
+		if (first_line) {
+		    min_value = minval_r
+		    max_value = maxval_r
+		    first_line = false
+		} else {
+		    if (minval_r < min_value)
+			min_value = minval_r
+		    if (maxval_r > max_value)
+			max_value = maxval_r
+		}
+	    }
+	}
+end