Initial commit

1 files changed, 308 insertions, 0 deletions

diff --git a/noao/imred/bias/colbias.x b/noao/imred/bias/colbias.x
new file mode 100644
index 00000000..f7aac840
--- /dev/null
+++ b/noao/imred/bias/colbias.x
@@ -0,0 +1,308 @@
+include	<imhdr.h>
+include	<imio.h>
+include <pkg/gtools.h>
+include	<pkg/xtanswer.h>
+
+# COLBIAS -- Remove line by line bias from images.
+#
+# A one dimensional bias vector is extracted from the bias columns.
+# A function is fit to the bias vector and the function is subtracted
+# from the image lines.  A trim section may be specified to output
+# only a part of the bias subtracted image.
+
+# Control procedure for mapping the images.
+#
+# The input and output images are given by image templates.  The
+# number of output images must match the number of input images.  Image
+# sections are not allowed.  The output image may be the same as the input
+# image.
+
+procedure colbias ()
+
+int	listin				# List of input images
+int	listout				# List of output images
+int	logfiles			# List of log files
+char	biassec[SZ_FNAME]		# Bias section
+char	trimsec[SZ_FNAME]		# Trim section
+int	median				# Use median of bias section?
+int	interactive			# Interactive?
+
+char	function[SZ_LINE]		# Curve fitting function
+int	order				# Order of curve fitting function
+
+char	image[SZ_FNAME]
+char	input[SZ_FNAME]
+char	biasimage[SZ_FNAME]
+char	output[SZ_FNAME]
+char	logfile[SZ_FNAME]
+char	original[SZ_FNAME]
+char	title[SZ_LINE]
+
+int	logfd
+pointer	in, bias, out, ic, gt
+
+int	clgeti(), clpopnu(), clgfil(), open(), gt_init(), nowhite()
+int	imtopen(), imtlen(), imtgetim(), btoi()
+bool	clgetb()
+long	clktime()
+pointer	immap()
+real	clgetr()
+
+begin
+	# Get input and output lists and check that the number of images
+	# are the same.
+
+	call clgstr ("input", title, SZ_LINE)
+	listin = imtopen (title)
+	call clgstr ("output", title, SZ_LINE)
+	listout = imtopen (title)
+	if (imtlen (listin) != imtlen (listout)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call error (0, "Input and output image lists do not match")
+	}
+
+	# Get the bias and trim sections.
+
+	call clgstr ("bias", biassec, SZ_FNAME)
+	call clgstr ("trim", trimsec, SZ_FNAME)
+	if (nowhite (biassec, biassec, SZ_FNAME) == 0)
+	    ;
+	if (nowhite (trimsec, trimsec, SZ_FNAME) == 0)
+	    ;
+	median = btoi (clgetb ("median"))
+
+	# Determine if the task is interactive.  If not set the interactive
+	# flag to always no.
+
+	if (clgetb ("interactive"))
+	    interactive = YES
+	else
+	    interactive = ALWAYSNO
+
+	# Initialize the curve fitting package.
+
+	call ic_open (ic)
+	call clgstr ("function", function, SZ_LINE)
+	call ic_pstr (ic, "function", function)
+	order = clgeti ("order")
+	call ic_puti (ic, "order", order)
+	call ic_putr (ic, "low", clgetr ("low_reject"))
+	call ic_putr (ic, "high", clgetr ("high_reject"))
+	call ic_puti (ic, "niterate", clgeti ("niterate"))
+	call ic_pstr (ic, "xlabel", "Line")
+	call ic_pstr (ic, "ylabel", "Bias")
+
+	gt = gt_init ()
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Get the log files.
+
+	logfiles = clpopnu ("logfiles")
+
+	# For each input and output image map the bias image, the
+	# trimmed input image, and the output image.  Use a temporary
+	# image header for overwritting the input image.
+
+	while ((imtgetim (listin, image, SZ_FNAME) != EOF) &&
+	    (imtgetim (listout, output, SZ_FNAME) != EOF)) {
+
+	    call sprintf (biasimage, SZ_FNAME, "%s%s")
+		call pargstr (image)
+		call pargstr (biassec)
+	    call sprintf (input, SZ_FNAME, "%s%s")
+		call pargstr (image)
+		call pargstr (trimsec)
+
+	    in = immap (input, READ_ONLY, 0)
+	    bias = immap (biasimage, READ_ONLY, 0)
+	    call xt_mkimtemp (image, output, original, SZ_FNAME)
+	    out = immap (output, NEW_COPY, in)
+	    IM_PIXTYPE(out) = TY_REAL
+
+	    call sprintf (title, SZ_LINE, "colbias %s")
+		call pargstr (image)
+	    call xt_answer (title, interactive)
+	    call gt_sets (gt, GTTITLE, title)
+
+	    # Enter log header.
+
+	    while (clgfil (logfiles, logfile, SZ_FNAME) != EOF) {
+		logfd = open (logfile, APPEND, TEXT_FILE)
+	        call cnvtime (clktime (0), title, SZ_LINE)
+	        call fprintf (logfd, "\nCOLBIAS:  %s\n")
+		    call pargstr (title)
+		call fprintf (logfd, "input = %s\noutput = %s\nbias = %s\n")
+		    call pargstr (input)
+		    call pargstr (output)
+		    call pargstr (biasimage)
+		if (median == YES)
+		    call fprintf (logfd, "Median of bias section used.\n")
+		call close (logfd)
+	    }
+	    call clprew (logfiles)
+
+	    call cb_colbias (in, bias, out, ic, gt, median, logfiles,
+		interactive)
+
+	    call imunmap (in)
+	    call imunmap (bias)
+	    call imunmap (out)
+	    call xt_delimtemp (output, original)
+	}
+
+	call ic_closer (ic)
+	call gt_free (gt)
+	call clpcls (logfiles)
+	call imtclose (listin)
+	call imtclose (listout)
+end
+
+
+# CB_COLBIAS -- Get an average column bias vector from the bias image.
+# Fit a function to the bias vector and subtract it from the input image
+# to form the output image.  Line coordinates are in terms of the full
+# input image.
+
+procedure cb_colbias (in, bias, out, ic, gt, median, logfiles, interactive)
+
+pointer	in			# Input image pointer
+pointer	bias			# Bias image pointer
+pointer	out			# Output image pointer
+pointer	ic			# ICFIT pointer
+pointer	gt			# GTOOLS pointer
+int	median			# Median of bias section?
+int	logfiles		# List of log files
+int	interactive		# Interactive curve fitting?
+
+char	graphics[SZ_FNAME]	# Graphics output device
+char	logfile[SZ_FNAME]
+int	i, nbias, nx, ny, ydim, yoff, ystep, ylen
+real	y, z
+pointer	cv, gp, sp, ybias, zbias, wts
+
+int	clgfil()
+real	cveval()
+pointer	gopen(), imgl2r(), impl2r()
+
+begin
+	# The bias coordinates are in terms of the full input image because
+	# the input and bias images may have different sections.
+
+	nx = IM_LEN(in, 1)
+	ny = IM_LEN(in, 2)
+
+	ydim = IM_VMAP(in, 2)
+	yoff = IM_VOFF(in, ydim)
+	ystep = IM_VSTEP(in, ydim)
+	ylen = IM_SVLEN(in, ydim)
+
+	# Get the bias vector and set the weights.
+
+	call cb_getcolbias (bias, ybias, zbias, nbias, median)
+	call smark (sp)
+	call salloc (wts, nbias, TY_REAL)
+	call amovkr (1., Memr[wts], nbias)
+
+	# Do the curve fitting using the interactive curve fitting package.
+	# Free memory when the fit is complete.
+
+	call ic_putr (ic, "xmin", 1.)
+	call ic_putr (ic, "xmax", real (ylen))
+	if ((interactive == YES) || (interactive == ALWAYSYES)) {
+	    call clgstr ("graphics", graphics, SZ_FNAME)
+	    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+	    call gt_setr (gt, GTXMIN, 1.)
+	    call gt_setr (gt, GTXMAX, real (ylen))
+	    call icg_fit (ic, gp, "cursor", gt, cv, Memr[ybias], Memr[zbias],
+		Memr[wts], nbias)
+	    call gclose (gp)
+	} else {
+	    call ic_fit (ic, cv, Memr[ybias], Memr[zbias], Memr[wts], nbias,
+		YES, YES, YES, YES)
+	}
+
+	# Log the fitting information.
+
+	while (clgfil (logfiles, logfile, SZ_FNAME) != EOF) {
+	    call ic_show (ic, logfile, gt)
+	    call ic_errors (ic, logfile, cv, Memr[ybias], Memr[zbias],
+		Memr[wts], nbias)
+	}
+	call clprew (logfiles)
+
+	call mfree (ybias, TY_REAL)
+	call mfree (zbias, TY_REAL)
+	call sfree (sp)
+
+	# Subtract the bias function from the input image.
+
+	do i = 1, ny {
+	    y = yoff + i * ystep
+	    z = cveval (cv, y)
+	    call asubkr (Memr[imgl2r(in,i)], z, Memr[impl2r(out,i)], nx)
+	}
+
+	# Free curve fitting memory.
+
+	call cvfree (cv)
+end
+
+
+# CB_GETCOLBIAS -- Get the column bias vector.
+# The ybias line values are in terms of the full image.
+
+procedure cb_getcolbias (bias, ybias, zbias, nbias, median)
+
+pointer	bias			# Bias image pointer
+pointer	ybias, zbias		# Bias vector
+int	nbias			# Number of bias points
+int	median			# Median of bias section?
+
+int	i, nx, ny, ydim, yoff, ystep
+
+real	amedr(), asumr()
+pointer	imgl1r(), imgl2r()
+
+begin
+	# Check for a bias consisting of a single column which is turned
+	# into a 1D image by IMIO.
+	if (IM_NDIM(bias) == 1) {
+	    ny = IM_LEN(bias, 1)
+	    ydim = IM_VMAP(bias, 1)
+	    yoff = IM_VOFF(bias, ydim)
+	    ystep = IM_VSTEP(bias, ydim)
+
+	    nbias = ny
+	    call malloc (ybias, nbias, TY_REAL)
+	    call malloc (zbias, nbias, TY_REAL)
+
+	    do i = 1, nbias
+		Memr[ybias+i-1] = yoff + i * ystep
+	    call amovr (Memr[imgl1r(bias)], Memr[zbias], nbias)
+
+	    return
+	}
+		
+	nx = IM_LEN(bias, 1)
+	ny = IM_LEN(bias, 2)
+	ydim = IM_VMAP(bias, 2)
+	yoff = IM_VOFF(bias, ydim)
+	ystep = IM_VSTEP(bias, ydim)
+	
+	nbias = ny
+	call malloc (ybias, nbias, TY_REAL)
+	call malloc (zbias, nbias, TY_REAL)
+
+	if (median == YES) {
+	    do i = 1, ny {
+	        Memr[ybias+i-1] = yoff + i * ystep
+	        Memr[zbias+i-1] = amedr (Memr[imgl2r(bias,i)], nx)
+	    }
+	} else {
+	    do i = 1, ny {
+	        Memr[ybias+i-1] = yoff + i * ystep
+	        Memr[zbias+i-1] = asumr (Memr[imgl2r(bias,i)], nx) / nx
+	    }
+	}
+end