Initial commit

1 files changed, 447 insertions, 0 deletions

diff --git a/noao/imred/quadred/src/quad/quadsections.x b/noao/imred/quadred/src/quad/quadsections.x
new file mode 100644
index 00000000..1d21d94c
--- /dev/null
+++ b/noao/imred/quadred/src/quad/quadsections.x
@@ -0,0 +1,447 @@
+include <imhdr.h>
+include "quadgeom.h"
+
+# QUADSECTIONS.X -- File header comment.
+
+define	OPTION_DICT	"|biassec|datasec|trimsec|reflect|duplicate|"
+
+define	OPT_BIASSEC	1
+define	OPT_DATASEC	2
+define	OPT_TRIMSEC	3
+define	OPT_REFLECT	4
+define	OPT_DUPLICATE	5
+
+define	DEFAULT_MACRO	"$I$S\\n"
+
+
+# QUADSECTIONS -- Quadsections task.
+
+procedure t_quadsections ()
+
+pointer	inlist			#TI List of input image name.
+char	instrument[SZ_FNAME]	#TI Instrument translation file
+char	section[SZ_LINE]	#TI Section of CCD required
+int	option			#TI Type of section required
+char	format[SZ_LINE]		#TI strmac macro string for building output
+int	xtrim1			#TI X pixels to trim at start of data
+int	xtrim2			#TI X pixels to trim at end   of data
+int	ytrim1			#TI Y pixels to trim at start of data
+int	ytrim2			#TI Y pixels to trim at end   of data
+int	xskip1			#TI X pixels to skip at start of overscan
+int	xskip2			#TI X pixels to skip at end   of overscan
+ 
+char	buffer[SZ_LINE], input[SZ_LINE]
+
+int	clgwrd(), imtopenp(), imtgetim()
+
+begin
+	# Open input image
+	inlist = imtopenp ("images")
+
+	# Open instrument file
+	call clgstr    ("instrument",  instrument,  SZ_FNAME)
+	call hdmopen   (instrument)
+
+	# Get option.
+	option = clgwrd ("window", buffer, SZ_LINE, OPTION_DICT)
+
+	# Get default section
+	call clgstr ("section", section, SZ_LINE)
+
+	# Get user defined trim and overscan margins
+	#xtrim1 = clgeti ("xtrim1")
+	#xtrim2 = clgeti ("xtrim2")
+	#ytrim1 = clgeti ("ytrim1")
+	#ytrim2 = clgeti ("ytrim2")
+	#xskip1 = clgeti ("xskip1")
+	#xskip2 = clgeti ("xskip2")
+	xtrim1 = INDEFI
+	xtrim2 = INDEFI
+	ytrim1 = INDEFI
+	ytrim2 = INDEFI
+	xskip1 = INDEFI
+	xskip2 = INDEFI
+
+	# Get format string and convert to a strmac macro string.
+	call clgstr ("template", format, SZ_LINE)
+	if (format[1] == EOS)
+	    call sprintf (format, SZ_LINE, "%s")
+		call pargstr (DEFAULT_MACRO)
+	call qsmkmacro (format, format, SZ_LINE)
+
+	while (imtgetim (inlist, input, SZ_LINE) != EOF)  {
+	    call quadsections (input, section, xtrim1, xtrim2, ytrim1, ytrim2,
+	    xskip1, xskip2, option, format)
+	}
+
+	# Tidy up
+	call hdmclose ()
+	call imtclose (inlist)
+end
+
+
+# QUADSECTIONS -- ???
+
+procedure quadsections (input, section, xtrim1, xtrim2, ytrim1, ytrim2,
+	xskip1, xskip2, option, format)
+
+char	input[SZ_FNAME]		#I Input image name.
+char	section[SZ_LINE]	#I Default section specification
+int	xtrim1			#I X pixels to trim at start of data
+int	xtrim2			#I X pixels to trim at end   of data
+int	ytrim1			#I Y pixels to trim at start of data
+int	ytrim2			#I Y pixels to trim at end   of data
+int	xskip1			#I X pixels to skip at start of overscan
+int	xskip2			#I X pixels to skip at end   of overscan
+int	option			#I Type of section required
+char	format[SZ_LINE]		#I strmac macro string for building output
+
+char	image[SZ_LINE], argstr[SZ_LINE], buffer[SZ_LINE]
+char	insection[SZ_LINE], outsection[SZ_LINE]
+int	amp, arg[9], len
+int	i, x1,  x2,  y1,  y2
+
+pointer	in, qg
+bool	reduced
+
+pointer	immap()
+int	hdmaccf(), gstrcpy(), strlen(), strmac()
+bool	quadsect()
+
+begin
+	# Parse input into an image name and an image section.
+	call imgimage   (input, image,     SZ_LINE)
+	call imgsection (input, insection, SZ_LINE)
+	# if no section was supplied in the image name use the default section.
+	if (insection[1] == EOS) {
+	    call strcpy (section, insection, SZ_LINE)
+	}
+
+	# Open input image
+	in = immap  (image, READ_ONLY, 0)
+
+	# Determine if image has been trimmed or not
+	reduced = (hdmaccf (in, "trim") == YES)
+
+	if (reduced) { 
+	    # OPT_BIASSEC does not make sense for reduced images.
+	    if (option == OPT_BIASSEC)
+		return
+	    # Trimsec and datasec are identical for reduced images
+	    if (option == OPT_TRIMSEC)
+		option = OPT_DATASEC
+	}
+
+	# Set-up quadgeom structure
+	call quadalloc (qg)
+        if (hdmaccf (in, "HDR_REV") == NO) {
+	    if (reduced) {
+		call quadgeomred (in, qg)
+	    } else {
+		call quadgeom (in, qg, "", "")
+	    }
+        } else {
+            call qghdr2 (in, qg)
+        }
+
+#	call quaddump (qg)
+
+	# Adjust quadgeom structure for user trim and overscan margins
+	if (! reduced) {
+	    call qguser (qg, xtrim1, xtrim2, ytrim1, ytrim2, xskip1, xskip2)
+	}
+#	call quaddump (qg)
+
+
+	# Store image name as first argument in macro argument string "argstr"
+	arg[1] = 1
+	arg[2] = 1 + arg[1] + gstrcpy (image, argstr, SZ_LINE)
+
+	# Blank output string
+	buffer[1] = EOS
+
+	# Determine the intersection of the specified section with the portion
+	# of the image read through each readout.
+	do amp = 1, QG_NAMPS (qg) {
+
+	    # skip any phantoms in raw images
+	    if (QG_PHANTOM (qg, amp) == NO) {
+
+		if (quadsect (qg, insection, option, amp, x1, x2, y1, y2)) {
+		    # Build resulting section string ...
+		    call sprintf (outsection, SZ_LINE, "[%d:%d,%d:%d]")
+			call pargi (x1)
+			call pargi (x2)
+			call pargi (y1)
+			call pargi (y2)
+
+		    # ... and save it as second argument
+		    arg[3] = 1 + arg[2] + gstrcpy (outsection, argstr[arg[2]],
+		    SZ_LINE-arg[2]+1)
+
+		    # Save Ampid as third argument
+		    call strcpy (Memc[QG_AMPID(qg, amp)], argstr[arg[3]],
+		    SZ_LINE-arg[3]+1)
+
+		    # Process macro string
+		    i = strmac (format, argstr, buffer, SZ_LINE)
+		    call printf (buffer)
+		}
+	    }
+	}
+
+	# Output <lf> if format does not explicitly include one.
+	len = strlen (buffer)
+	if ((len > 2) && !(buffer[len-1]=='\\' && buffer[len]=='n')) {
+	    call printf ("\n")
+	}
+
+	call flush (STDOUT)
+
+	# Tidy up
+	call quadfree (qg)
+	call imunmap (in)
+end
+
+
+# QSMKMACRO -- Perform the following substitutions on the given macro string
+#
+# $I --> $1
+# $S --> $2
+# $A --> $3
+# $? --> $?
+
+procedure qsmkmacro (instr, outstr, maxchars)
+
+char	instr[ARB]			#I Input macro string.
+char	outstr[maxchars]		#O Output macro string.
+int	maxchars			#I Maximum length of outstr
+
+char	ch
+int	ip, op
+
+begin
+
+	op = 1
+	for (ip=1; instr[ip] != EOS; ip=ip+1) {
+	    ch = instr[ip]
+	    outstr[op] = ch
+	    op = op + 1
+	    if (op > maxchars)
+		call error (0, "qsmkmacro: Output buffer overflow")
+	    
+	    if (ch == '$') {
+		ip = ip + 1
+		ch = instr[ip]
+
+		if (ch == 'I') {
+		    outstr (op) = '1'
+		    op = op + 1
+		} else if (ch == 'S') {
+		    outstr (op) = '2'
+		    op = op + 1
+		} else if (ch == 'A') {
+		    outstr (op) = '3'
+		    op = op + 1
+		} else {
+		    outstr (op) = ch
+		    op = op + 1
+		}
+		if (op > maxchars)
+		    call error (0, "qsmkmacro: Output buffer overflow")
+	    }
+	}
+end
+
+
+# QUADSECT -- ??
+
+bool procedure quadsect (qg, section, option, amp, x1, x2, y1, y2)
+
+pointer	qg			#I Pointer to initialised quadgeom structure.
+char	section[SZ_LINE]	#I Default section specification.
+int	option			#I Type of section required.
+int	amp			#I Amplifier for which section is required.
+int	x1, x2, y1, y2		#O Corners of specified section.
+bool	overlap			#O true if part of section read through amp.
+
+int	xskip, xsize, yskip, ysize
+int	sx1, sx2, sy1, sy2, sxs, sys
+int	dx1, dx2, dy1, dy2
+int	tx1, tx2, ty1, ty2
+int	bx1, bx2, by1, by2
+
+begin
+
+	# Decode input section
+	x1  = 1
+	x2  = QG_NX(qg, 0)
+	sxs = 1
+	y1  = 1
+	y2  = QG_NY(qg, 0)
+	sys = 1
+	call ccd_section (section, x1, x2, sxs, y1, y2, sys)
+	sx1 = min (x1, x2)
+	sx2 = max (x1, x2)
+	sy1 = min (y1, y2)
+	sy2 = max (y1, y2)
+
+	# Set up null return (overlap) values in case no part of section was
+	# read with this amplifier.
+	overlap = false
+	x1 = 0
+	x2 = 0
+	y1 = 0
+	y2 = 0
+
+	# Calculate suplimentary quantitiies as required.
+	switch (option) {
+
+	case OPT_REFLECT, OPT_DUPLICATE:
+	    xskip = sx1 - QG_DX1(qg, 0)
+	    xsize = sx2 - sx1 + 1
+	    yskip = sy1 - QG_DY1(qg, 0)
+	    ysize = sy2 - sy1 + 1
+	}
+
+	# Determine the intersection of the specified section with the portion
+	# of the image read through the specified readout.
+	switch (option) {
+
+	case OPT_BIASSEC:
+	    bx1 = QG_AX1(qg, amp) + QG_BX1(qg, amp) - 1
+	    bx2 = QG_AX1(qg, amp) + QG_BX2(qg, amp) - 1
+	    by1 = QG_AY1(qg, amp) + QG_BY1(qg, amp) - 1
+	    by2 = QG_AY1(qg, amp) + QG_BY2(qg, amp) - 1
+
+	    if (sx1 > bx2)
+		return (overlap)
+	    if (sx2 < bx1)
+		return (overlap)
+	    if (sy1 > by2)
+		return (overlap)
+	    if (sy2 < by1)
+		return (overlap)
+
+	    x1 = max (sx1, bx1)
+	    x2 = min (sx2, bx2)
+	    y1 = max (sy1, by1)
+	    y2 = min (sy2, by2)
+
+	case OPT_DATASEC:
+	    dx1 = QG_AX1(qg, amp) + QG_DX1(qg, amp) - 1
+	    dx2 = QG_AX1(qg, amp) + QG_DX2(qg, amp) - 1
+	    dy1 = QG_AY1(qg, amp) + QG_DY1(qg, amp) - 1
+	    dy2 = QG_AY1(qg, amp) + QG_DY2(qg, amp) - 1
+
+	    if (sx1 > dx2)
+		return (overlap)
+	    if (sx2 < dx1)
+		return (overlap)
+	    if (sy1 > dy2)
+		return (overlap)
+	    if (sy2 < dy1)
+		return (overlap)
+
+	    x1 = max (sx1, dx1)
+	    x2 = min (sx2, dx2)
+	    y1 = max (sy1, dy1)
+	    y2 = min (sy2, dy2)
+
+	case OPT_TRIMSEC:
+	    tx1 = QG_AX1(qg, amp) + QG_TX1(qg, amp) - 1
+	    tx2 = QG_AX1(qg, amp) + QG_TX2(qg, amp) - 1
+	    ty1 = QG_AY1(qg, amp) + QG_TY1(qg, amp) - 1
+	    ty2 = QG_AY1(qg, amp) + QG_TY2(qg, amp) - 1
+
+	    if (sx1 > tx2)
+		return (overlap)
+	    if (sx2 < tx1)
+		return (overlap)
+	    if (sy1 > ty2)
+		return (overlap)
+	    if (sy2 < ty1)
+		return (overlap)
+
+	    x1 = max (sx1, tx1)
+	    x2 = min (sx2, tx2)
+	    y1 = max (sy1, ty1)
+	    y2 = min (sy2, ty2)
+
+	case OPT_REFLECT:
+	    dx1 = QG_AX1(qg, amp) + QG_DX1(qg, amp) - 1
+	    dx2 = QG_AX1(qg, amp) + QG_DX2(qg, amp) - 1
+	    dy1 = QG_AY1(qg, amp) + QG_DY1(qg, amp) - 1
+	    dy2 = QG_AY1(qg, amp) + QG_DY2(qg, amp) - 1
+
+	    switch (QG_AMPTYPE(qg, amp)) {
+	    case AMP11:
+		x1  = dx1 + xskip
+		x2  = x1  + xsize - 1
+		y1  = dy1 + yskip
+		y2  = y1  + ysize - 1
+
+	    case AMP12:
+		x2  = dx2 - xskip
+		x1  = x2  - xsize + 1
+		y1  = dy1 + yskip
+		y2  = y1  + ysize - 1
+
+	    case AMP21:
+		x1  = dx1 + xskip
+		x2  = x1  + xsize - 1
+		y2  = dy2 - yskip
+		y1  = y2  - ysize + 1
+
+	    case AMP22:
+		x2  = dx2 - xskip
+		x1  = x2  - xsize + 1
+		y2  = dy2 - yskip
+		y1  = y2  - ysize + 1
+	    }
+		
+	    if (x1 > dx2)
+		return (overlap)
+	    if (x2 < dx1)
+		return (overlap)
+	    if (y1 > dy2)
+		return (overlap)
+	    if (y2 < dy1)
+		return (overlap)
+
+	    x1 = max (x1, dx1)
+	    x2 = min (x2, dx2)
+	    y1 = max (y1, dy1)
+	    y2 = min (y2, dy2)
+
+	case OPT_DUPLICATE:
+	    dx1 = QG_AX1(qg, amp) + QG_DX1(qg, amp) - 1
+	    dx2 = QG_AX1(qg, amp) + QG_DX2(qg, amp) - 1
+	    dy1 = QG_AY1(qg, amp) + QG_DY1(qg, amp) - 1
+	    dy2 = QG_AY1(qg, amp) + QG_DY2(qg, amp) - 1
+
+	    x1  = dx1 + xskip
+	    x2  = x1  + xsize - 1
+	    y1  = dy1 + yskip
+	    y2  = y1  + ysize - 1
+
+	    if (x1 > dx2)
+		return (overlap)
+	    if (x2 < dx1)
+		return (overlap)
+	    if (y1 > dy2)
+		return (overlap)
+	    if (y2 < dy1)
+		return (overlap)
+
+	    x1 = max (x1, dx1)
+	    x2 = min (x2, dx2)
+	    y1 = max (y1, dy1)
+	    y2 = min (y2, dy2)
+
+	}
+
+	overlap = true
+	return (overlap)
+
+end