Initial commit

1 files changed, 159 insertions, 0 deletions

diff --git a/noao/imred/quadred/src/quad/quadscale.x b/noao/imred/quadred/src/quad/quadscale.x
new file mode 100644
index 00000000..5e594eb4
--- /dev/null
+++ b/noao/imred/quadred/src/quad/quadscale.x
@@ -0,0 +1,159 @@
+include <imhdr.h>
+include "quadgeom.h"
+
+define	OPPERATIONS	"|multiply|divide|"
+define	OPMULTIPLY	1
+define	OPDIVIDE	2
+define	TOL1		 1.0e-7
+define	TOL2		-1.0e-7
+
+procedure t_quadscale ()
+
+char	input[SZ_FNAME]		#TI Input image name.
+char	output[SZ_FNAME]	#TI Output image name.
+char    instrument[SZ_FNAME]    #TI Instrument translation file
+
+real	gain[QG_MAXAMPS]	#TI Gain factor for each quadrant
+int	op			#TI Multiply or divide by gain factors
+
+int	i
+pointer	in, out, qg
+char	buffer[SZ_LINE]
+
+real	clgetr()
+int	clgwrd(), hdmaccf()
+pointer	immap()
+
+begin
+        # Open instrument file
+        call clgstr ("instrument",  instrument,  SZ_FNAME)
+        call hdmopen (instrument)
+
+	# Open input image
+	call clgstr ("input",  input,  SZ_FNAME)
+	in = immap  (input, READ_ONLY, 0)
+
+	# Allocate quadgeom structure and initialise it from image header
+	call quadalloc (qg)
+
+
+	# Fill out quadgeom structure from header depending on revision level
+	if (hdmaccf (in, "HDR_REV") == NO) {
+
+	    # Check to see if the image has been processed or not
+	    if (hdmaccf (in, "ccdproc") == YES) {
+		call quadgeomred (in, qg)
+	    } else {
+		call quadgeom (in, qg, "", "")
+	    }
+
+	} else {
+	    call qghdr2 (in, qg)
+	}
+
+	# Open output image
+	call clgstr ("output",  output,  SZ_FNAME)
+	out = immap (output, NEW_COPY, in)
+	IM_PIXTYPE(out) = TY_REAL
+
+	# Get gain factors
+	gain[1] = clgetr ("gain11")
+	gain[2] = clgetr ("gain12")
+	gain[3] = clgetr ("gain21")
+	gain[4] = clgetr ("gain22")
+
+	# Get direction of opperation
+	op = clgwrd ("opperation", buffer, SZ_LINE, OPPERATIONS)
+
+	switch (op) {
+	case OPMULTIPLY:
+	    call quadmult (in, out, gain, qg)
+
+	case OPDIVIDE:
+	    # Check for zero gain --> divide by zero
+	    do i = 1, 4 {
+		if ((gain[i] < TOL1) && (gain[i] > TOL2)) {
+		    call error (0, "Attempt to divide by gain value of zero")
+		}
+	    }
+	    call quaddiv (in, out, gain, qg)
+
+	}
+
+	call imunmap (in)
+	call imunmap (out)
+	call quadfree (qg)
+	call hdmclose ()
+end
+
+procedure quadmult (in, out, gain, qg)
+
+pointer	in				#I imio pointer for input image.
+pointer	out				#I imio pointer for output image.
+real	gain[ARB]			#I Array of gain factors.
+pointer	qg				#I Pointer to quadgeom structure.
+
+pointer	inbuf, obuf
+int	junk, nx, x, y, line, amp, amp2, off
+long	invec[IM_MAXDIM], ovec[IM_MAXDIM]
+
+int	imgnlr(), impnlr()
+
+begin
+
+	# Setup start vector for sequential reads
+	call amovkl (long(1), invec, IM_MAXDIM)
+	call amovkl (long(1), ovec,  IM_MAXDIM)
+
+	do y = 1, QG_NAMPS(qg) {
+	    amp2 = QG_AMP(qg, 1, y)
+	    do line = 1, QG_NY(qg, amp2) {
+		junk = imgnlr (in,  inbuf, invec)
+		junk = impnlr (out, obuf,  ovec) 
+		off = 0
+		do x = 1, QG_NAMPSX(qg) {
+		    amp = QG_AMP(qg, x, y)
+		    nx = QG_NX(qg, amp)
+		    call amulkr (Memr[inbuf+off], gain[amp], Memr[obuf+off], nx)
+		    off = off + nx
+		}
+	    }
+	}
+
+end
+
+procedure quaddiv (in, out, gain, qg)
+
+pointer	in				#I imio pointer for input image.
+pointer	out				#I imio pointer for output image.
+real	gain[ARB]			#I Array of gain factors.
+pointer	qg				#I Pointer to quadgeom structure.
+
+pointer	inbuf, obuf
+int	junk, nx, x, y, line, amp, amp2, off
+long	invec[IM_MAXDIM], ovec[IM_MAXDIM]
+
+int	imgnlr(), impnlr()
+
+begin
+
+	# Setup start vector for sequential reads
+	call amovkl (long(1), invec, IM_MAXDIM)
+	call amovkl (long(1), ovec,  IM_MAXDIM)
+
+	do y = 1, QG_NAMPS(qg) {
+	    amp2 = QG_AMP(qg, 1, y)
+	    do line = 1, QG_NY(qg, amp2) {
+		junk = imgnlr (in,  inbuf, invec)
+		junk = impnlr (out, obuf,  ovec) 
+		off = 0
+		do x = 1, QG_NAMPSX(qg) {
+		    amp = QG_AMP(qg, x, y)
+		    nx = QG_NX(qg, amp)
+		    call adivkr (Memr[inbuf+off], gain[amp], Memr[obuf+off], nx)
+		    off = off + nx
+		}
+	    }
+	}
+
+end