1 files changed, 160 insertions, 0 deletions
diff --git a/sys/imio/iki/fxf/fxfplread.x b/sys/imio/iki/fxf/fxfplread.x
new file mode 100644
index 00000000..4d4c3e83
--- /dev/null
+++ b/sys/imio/iki/fxf/fxfplread.x
@@ -0,0 +1,160 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+ 
+include <syserr.h>
+include <plset.h>
+include <imhdr.h>
+include <imio.h>
+include <mach.h>
+include "fxf.h"
+
+
+# FXF_PLREAD -- Read a PLIO mask stored in a FITS binary table extension
+# and load it into an image descriptor.
+#
+# There is a builtin assumption in this code (also in fxf_plwrite) that
+# masks will not be more than 3-dimensional.  This could be generalized
+# if necessary, but we have never seen a mask of dimensionality higher
+# than 3.  The dimensionality, size, and depth of the mask is preserved.
+
+procedure fxf_plread (im)
+
+pointer	im			#I image descriptor
+
+char	kwname[SZ_KEYWORD]
+pointer sp, fk, pl, lp, ip, ix
+long	data_offset, data_len, heap_offset, llen, loff
+int	naxes, axlen[IM_MAXDIM], depth, maxlen
+int	fd, i, j, nelem, nlines, v[PL_MAXDIM], maxoff, nbytes
+
+long	note()
+bool	streq()
+int	imgeti(), pl_create(), miireadi(), miireads()
+errchk	imgeti, pl_create, miireadi, miireads, seek, pl_update, syserrs
+
+begin
+	call smark (sp)
+
+	fk = IM_KDES(im)
+	fd = IM_HFD(im)
+
+	# The maximum encoded line list length is (normally) passed in via
+	# the binary table format keywords, and stored in FIT_PLMAXLEN.
+
+	maxlen = FIT_PLMAXLEN(fk)
+	if (maxlen <= 0)
+	    maxlen = DEF_PLMAXLEN
+
+	# Scratch buffer for encoded line lists.
+	call salloc (lp, maxlen, TY_SHORT)
+
+	# Get the dimensionality and size of the stored mask.
+	call amovki (1, axlen, IM_MAXDIM)
+	naxes = imgeti (im, "ZNAXIS")
+	call fxf_filter_keyw (im, "ZNAXIS")
+	do i = 1, naxes {
+	    call sprintf (kwname, LEN_CARD, "ZNAXIS%d")
+		call pargi(i)
+	    axlen[i] = imgeti (im, kwname)
+	    call fxf_filter_keyw (im, kwname)
+	    call sprintf (kwname, LEN_CARD, "ZTILE%d")
+		call pargi(i)
+	    call fxf_filter_keyw (im, kwname)
+	}
+
+	# Get the mask depth, passed as compression algorithm parameter
+	# number 1 for a PLIO-compressed image.
+
+	depth = DEF_PLDEPTH
+	ifnoerr (call imgstr (im, "ZNAME1", kwname, SZ_KEYWORD)) {
+	    if (streq (kwname, "depth"))
+		iferr (depth = imgeti (im, "ZVAL1"))
+		    depth = DEF_PLDEPTH
+	    call fxf_filter_keyw (im, "ZNAME1")
+	    call fxf_filter_keyw (im, "ZVAL1")
+	    call fxf_filter_keyw (im, "ZBITPIX")
+	    call fxf_filter_keyw (im, "ZIMAGE")
+	}
+
+	# Create an initially empty mask of the given size.
+	pl = pl_create (naxes, axlen, depth)
+
+	# Create a buffer for the line list index (maxdim 3 assumed).
+	nlines = axlen[3] * axlen[2]
+	call salloc (ix, nlines * 2, TY_INT)
+
+	# Compute the file offsets of the table data and heap areas.  The
+	# file position is assumed to be already positioned at the start
+	# of the data area of the file.
+
+	data_offset = note (fd)
+	data_len = FIT_LENAXIS(fk,3) * FIT_LENAXIS(fk,2) * FIT_LENAXIS(fk,1)
+	heap_offset = data_offset + data_len/SZB_CHAR
+
+	# Read the line list index from the input file.  The index contains
+	# one entry for every line in the (possibly multidimensional) image.
+	# Each entry consists of two integer values, the length of the
+	# stored line list, and the heap offset (in bytes) of the stored list.
+
+	nelem = miireadi (fd, Memi[ix], nlines * 2)
+	if (nelem != nlines * 2)
+	    call syserrs (SYS_FXFRMASK, IM_NAME(im))
+	
+	# Find out the maximum offset value to determine if they were
+	# written using the 2 byte units rather than the standard (byte unit)
+
+	maxoff = 0
+        ip = ix
+	do j = 1, axlen[3] {
+            do i = 1, axlen[2] {
+	       maxoff = max (maxoff, Memi[ip+1]+2*Memi[ip]) 
+   	       ip = ip + 2
+	    }
+	}
+
+	if (maxoff < (FIT_PCOUNT(fk) - maxoff/2)) {
+	    nbytes = 1
+        } else {	
+	    nbytes = 2
+	}
+
+	# Read the line list data and insert it into the PLIO mask.  
+	# pl_update will be called for each line of the mask even if multiple
+	# lines point to the same line list data, but pl_update will sort
+	# all this out and restore the multiple references as the mask is
+	# built.
+
+	ip = ix
+	v[1] = 1
+
+	do j = 1, axlen[3] {
+	    v[3] = j
+	    do i = 1, axlen[2] {
+		v[2] = i
+
+		llen = Memi[ip]
+
+		# This offset on the table data is in byte units, convert
+		# to short.
+
+		loff = Memi[ip+1] / nbytes
+
+		call seek (fd, heap_offset + loff)
+		nelem = miireads (fd, Mems[lp], llen)
+		if (nelem != llen)
+		    call syserrs (SYS_FXFRMASK, IM_NAME(im))
+
+		call pl_update (pl, v, Mems[lp])
+
+		ip = ip + 2
+	    }
+	}
+
+	# Set up IMIO descriptor.
+	call amovl (axlen, IM_LEN(im,1), IM_MAXDIM)
+	call amovl (axlen, IM_PHYSLEN(im,1), IM_MAXDIM)
+        IM_NDIM(im) = naxes
+	IM_PIXTYPE(im) = TY_INT
+	IM_PL(im) = pl
+
+	call sfree (sp)
+end