Initial commit

1 files changed, 478 insertions, 0 deletions

diff --git a/noao/imred/generic/flat1d.x b/noao/imred/generic/flat1d.x
new file mode 100644
index 00000000..5a0797fc
--- /dev/null
+++ b/noao/imred/generic/flat1d.x
@@ -0,0 +1,478 @@
+include	<error.h>
+include	<imhdr.h>
+include <pkg/gtools.h>
+
+# FLAT1D -- Fit a function to image lines or columns and output an image
+# consisting of the ratio.  Set a minimum value test to the denominator.
+# The fitting parameters may be set interactively using the icfit package.
+
+procedure flat1d ()
+
+int	listin				# Input image list
+int	listout				# Output image list
+int	axis				# Image axis to fit
+real	minflat				# Minimum fit value for ratio
+bool	interactive			# Interactive?
+
+char	sample[SZ_LINE]			# Sample ranges
+int	naverage			# Sample averaging size
+char	function[SZ_LINE]		# Curve fitting function
+int	order				# Order of curve fitting function
+real	low_reject,  high_reject	# Rejection thresholds
+int	niterate			# Number of rejection iterations
+real	grow				# Rejection growing radius
+
+char	input[SZ_LINE]			# Input image
+char	output[SZ_FNAME]		# Output image
+pointer	in, out				# IMIO pointers
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+
+int	imtopen(), imtgetim(), imtlen(), gt_init()
+int	clgeti()
+real	clgetr()
+bool	clgetb()
+
+begin
+	# Get input and output lists and check that the number of images
+	# are the same.
+
+	call clgstr ("input", input, SZ_LINE)
+	listin = imtopen (input)
+	call clgstr ("output", input, SZ_LINE)
+	listout = imtopen (input)
+	if (imtlen (listin) != imtlen (listout)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call error (0, "Input and output image lists do not match")
+	}
+
+	# Get task parameters.
+
+	axis = clgeti ("axis")
+	minflat = clgetr ("minflat")
+	interactive = clgetb ("interactive")
+
+	# Initialize the ICFIT package.
+	call clgstr ("sample", sample, SZ_LINE)
+	naverage = clgeti ("naverage")
+	call clgstr ("function", function, SZ_LINE)
+	order = clgeti ("order")
+	low_reject = clgetr ("low_reject")
+	high_reject = clgetr ("high_reject")
+	niterate = clgeti ("niterate")
+	grow = clgetr ("grow")
+
+	call ic_open (ic)
+	call ic_pstr (ic, "sample", sample)
+	call ic_puti (ic, "naverage", naverage)
+	call ic_pstr (ic, "function", function)
+	call ic_puti (ic, "order", order)
+	call ic_putr (ic, "low", low_reject)
+	call ic_putr (ic, "high", high_reject)
+	call ic_puti (ic, "niterate", niterate)
+	call ic_putr (ic, "grow", grow)
+	call ic_pstr (ic, "ylabel", "")
+
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "line")
+
+	# Fit each input image.
+
+	while ((imtgetim (listin, input, SZ_LINE) != EOF) &&
+	    (imtgetim (listout, output, SZ_FNAME) != EOF)) {
+
+	    iferr (call f1d_immap (input, output, in, out)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    call f1d_flat1d (in, out, ic, gt, input, axis, minflat, interactive)
+	    call imunmap (in)
+	    call imunmap (out)
+	}
+
+	call ic_closer (ic)
+	call gt_free (gt)
+	call imtclose (listin)
+	call imtclose (listout)
+end
+
+
+# F1D_FLAT1D -- Given the image descriptor determine the fitting function
+# for each line or column and create an output image.  If the interactive flag
+# is set then set the fitting parameters interactively.
+
+define	MAXBUF	512 * 100		# Maximum number of pixels per block
+
+procedure f1d_flat1d (in, out, ic, gt, title, axis, minflat, interactive)
+
+pointer	in				# IMIO pointer for input image
+pointer	out				# IMIO pointer for output image
+pointer	ic				# ICFIT pointer
+pointer	gt				# GTOOLS pointer
+char	title[ARB]			# Title
+int	axis				# Image axis to fit
+real	minflat				# Minimum value for flat division
+bool	interactive			# Interactive?
+
+char	graphics[SZ_FNAME]
+int	i, nx, new
+real	mindata, maxdata
+pointer	cv, gp, sp, x, wts, indata, outdata
+
+int	f1d_getline(), f1d_getdata(), strlen()
+pointer	gopen()
+
+begin
+	# Error check.
+
+	if (IM_NDIM (in) > 2)
+	    call error (0, "Image dimensions > 2 are not implemented")
+	if (axis > IM_NDIM (in))
+	    call error (0, "Axis exceeds image dimension")
+
+	# Allocate memory for curve fitting.
+
+	nx = IM_LEN (in, axis)
+	call smark (sp)
+	call salloc (x, nx, TY_REAL)
+	call salloc (wts, nx, TY_REAL)
+
+	do i = 1, nx
+	    Memr[x+i-1] = i
+	call amovkr (1., Memr[wts], nx)
+
+	call ic_putr (ic, "xmin", Memr[x])
+	call ic_putr (ic, "xmax", Memr[x+nx-1])
+
+	# If the interactive flag is set then use icg_fit to set the
+	# fitting parameters.  Get_fitline returns EOF when the user
+	# is done.  The weights are reset since the user may delete
+	# points.
+
+	if (interactive) {
+	    call clgstr ("graphics", graphics, SZ_FNAME)
+	    gp = gopen (graphics, NEW_FILE, STDGRAPH)
+	    i = strlen (title)
+	    indata = NULL
+	    while (f1d_getline (ic, gt, in, axis, title, indata) != EOF) {
+		title[i + 1] = EOS
+	        call icg_fit (ic, gp, "cursor", gt, cv, Memr[x], Memr[indata],
+		    Memr[wts], nx)
+		call amovkr (1., Memr[wts], nx)
+	    }
+	    call gclose (gp)
+	}
+
+	# Loop through the input image and create an output image.
+
+	new = YES
+
+	while (f1d_getdata (in, out, axis, MAXBUF, indata, outdata) != EOF) {
+
+	    call alimr (Memr[indata], nx, mindata, maxdata)
+	    if (maxdata >= minflat) {
+	        call ic_fit (ic, cv, Memr[x], Memr[indata], Memr[wts],
+		    nx, new, YES, new, new)
+	        new = NO
+
+	        call cvvector (cv, Memr[x], Memr[outdata], nx)
+	    }
+
+	    call f1d_flat (Memr[indata], Memr[outdata], Memr[outdata], nx,
+		minflat, mindata, maxdata)
+	}
+
+	call imaddr (out, "ccdmean", 1.)
+
+	call cvfree (cv)
+	call sfree (sp)
+end
+
+
+# F1D_IMMAP -- Map images for flat1d.
+
+procedure f1d_immap (input, output, in, out)
+
+char	input[ARB]		# Input image
+char	output[ARB]		# Output image
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+
+pointer	sp, root, sect, line, data
+
+int	access(), impnlr()
+pointer	immap()
+errchk	immap
+
+begin
+	# Get the root name and section of the input image.
+
+	call smark (sp)
+	call salloc (root, SZ_FNAME, TY_CHAR)
+	call salloc (sect, SZ_FNAME, TY_CHAR)
+
+	call get_root (input, Memc[root], SZ_FNAME)
+	call get_section (input, Memc[sect], SZ_FNAME)
+
+	# If the output image is not accessible then create it as a new copy
+	# of the full input image and initialize to unit response.
+
+	if (access (output, READ_WRITE, BINARY_FILE) == NO) {
+	    in = immap (Memc[root], READ_ONLY, 0)
+	    out = immap (output, NEW_COPY, in)
+	    IM_PIXTYPE(out) = TY_REAL
+
+	    call salloc (line, IM_MAXDIM, TY_LONG)
+	    call amovkl (long (1), Meml[line], IM_MAXDIM)
+	    while (impnlr (out, data, Meml[line]) != EOF)
+	        call amovkr (1., Memr[data], IM_LEN(out, 1))
+
+	    call imunmap (in)
+	    call imunmap (out)
+	}
+
+	# Map the input and output images.
+
+	in = immap (input, READ_ONLY, 0)
+
+	call sprintf (Memc[root], SZ_FNAME, "%s%s")
+	    call pargstr (output)
+	    call pargstr (Memc[sect])
+	out = immap (Memc[root], READ_WRITE, 0)
+
+	call sfree (sp)
+end
+
+
+# F1D_GETDATA -- Get a line of image data.
+
+int procedure f1d_getdata (in, out, axis, maxbuf, indata, outdata)
+
+pointer	in			# Input IMIO pointer
+pointer	out			# Output IMIO pointer
+int	axis			# Image axis
+int	maxbuf			# Maximum buffer size for column axis
+pointer	indata			# Input data pointer
+pointer	outdata			# Output data pointer
+
+int	i, index, last_index, col1, col2, nc, ncols, nlines, ncols_block
+pointer	inbuf, outbuf, ptr
+
+pointer	imgl1r(), impl1r(), imgl2r(), impl2r(), imgs2r(), imps2r()
+
+data	index/0/
+
+begin
+	# Increment to the next image vector.
+
+	index = index + 1
+
+	# Initialize for the first vector.
+
+	if (index == 1) {
+	    ncols = IM_LEN (in, 1)
+	    if (IM_NDIM (in) == 1)
+		nlines = 1
+	    else
+		nlines = IM_LEN (in, 2)
+
+	    switch (axis) {
+	    case 1:
+		last_index = nlines
+	    case 2:
+		last_index = ncols
+	        ncols_block = max (1, min (ncols, maxbuf / nlines))
+		col2 = 0
+
+	        call malloc (indata, nlines, TY_REAL)
+	        call malloc (outdata, nlines, TY_REAL)
+	    }
+	}
+
+	# Finish up if the last vector has been done.
+
+	if (index > last_index) {
+	    if (axis == 2) {
+	        ptr = outbuf + index - 1 - col1
+	        do i = 1, nlines {
+		    Memr[ptr] = Memr[outdata+i-1]
+		    ptr = ptr + nc
+	        }
+
+	        call mfree (indata, TY_REAL)
+	        call mfree (outdata, TY_REAL)
+	    }
+
+	    index = 0
+	    return (EOF)
+	}
+
+	# Get the next image vector.
+
+	switch (axis) {
+	case 1:
+	    if (IM_NDIM (in) == 1) {
+		indata = imgl1r (in)
+		outdata = impl1r (out)
+	    } else {
+		indata = imgl2r (in, index)
+		outdata = impl2r (out, index)
+	    }
+	case 2:
+	    if (index > 1) {
+		ptr = outbuf + index - 1 - col1
+		do i = 1, nlines {
+		    Memr[ptr] = Memr[outdata+i-1]
+		    ptr = ptr + nc
+		}
+	    }
+
+	    if (index > col2) {
+		col1 = col2 + 1
+		col2 = min (ncols, col1 + ncols_block - 1)
+		inbuf = imgs2r (in, col1, col2, 1, nlines)
+		outbuf = imps2r (out, col1, col2, 1, nlines)
+		nc = col2 - col1 + 1
+	    }
+
+	    ptr = inbuf + index - col1
+	    do i = 1, nlines {
+		Memr[indata+i-1] = Memr[ptr]
+		ptr = ptr + nc
+	    }
+	}
+	return (index)
+end
+
+# F1D_FLAT -- For the flat field values by ratioing the image data by the fit.
+# If the fit value is less than minflat then the ratio is set to 1.
+
+procedure f1d_flat (data, fit, flat, npts, minflat, mindata, maxdata)
+
+real	data[npts]		# Image data
+real	fit[npts]		# Fit to image data
+real	flat[npts]		# Ratio of image data to the fit
+int	npts			# Number of points
+real	minflat			# Minimum fit value for ratio
+real	mindata			# Minimum data value
+real	maxdata			# Maximum data value
+
+int	i
+
+begin
+	if (mindata >= minflat)
+	    call adivr (data, fit, flat, npts)
+
+	else if (maxdata < minflat)
+	    call amovkr (1., flat, npts)
+
+	else {
+	    do i = 1, npts {
+	        if (fit[i] < minflat)
+		    flat[i] = 1.
+	        else
+		    flat[i] = data[i] / fit[i]
+	    }
+	}
+end
+
+
+# F1D_GETLINE -- Get image data to be fit interactively.  Return EOF
+# when the user enters EOF or CR.  Default is 1 and the out of bounds
+# requests are silently limited to the nearest in edge.
+
+int procedure f1d_getline (ic, gt, im, axis, title, data)
+
+pointer	ic			# ICFIT pointer
+pointer	gt			# GTOOLS pointer
+pointer	im			# IMIO pointer
+int	axis			# Image axis
+char	title[ARB]		# Title
+pointer	data			# Image data
+
+char	line[SZ_LINE]
+int	i, j, stat, imlen
+pointer	x
+
+int	getline(), nscan()
+pointer	imgl1r()
+
+data	stat/EOF/
+
+begin
+	# If the image is one dimensional do not prompt.
+
+	if (IM_NDIM (im) == 1) {
+	    if (stat == EOF) {
+		call sprintf (title, SZ_LINE, "%s\n%s")
+	    	    call pargstr (title)
+	    	    call pargstr (IM_TITLE(im))
+		call gt_sets (gt, GTTITLE, title)
+		call mfree (data, TY_REAL)
+		call malloc (data, IM_LEN(im, 1), TY_REAL)
+		call amovr (Memr[imgl1r(im)], Memr[data], IM_LEN(im, 1))
+		stat = OK
+	    } else
+		stat = EOF
+
+	    return (stat)
+	}
+
+	# If the image is two dimensional prompt for the line or column.
+
+	switch (axis) {
+	case 1:
+	    imlen = IM_LEN (im, 2)
+	    call sprintf (title, SZ_LINE, "%s: Fit line =")
+	        call pargstr (title)
+	case 2:
+	    imlen = IM_LEN (im, 1)
+	    call sprintf (title, SZ_LINE, "%s: Fit column =")
+	        call pargstr (title)
+	}
+
+	call printf ("%s ")
+	    call pargstr (title)
+	call flush (STDOUT)
+
+	if (getline(STDIN, line) == EOF)
+	    return (EOF)
+
+	call sscan (line)
+	call gargi (i)
+	call gargi (j)
+
+	switch (nscan()) {
+	case 0:
+	    stat = EOF
+	    return (stat)
+	case 1:
+	    i = max (1, min (imlen, i))
+	    j = i
+	case 2:
+	    i = max (1, min (imlen, i))
+	    j = max (1, min (imlen, j))
+	}
+
+	call sprintf (title, SZ_LINE, "%s %d - %d\n%s")
+	    call pargstr (title)
+	    call pargi (i)
+	    call pargi (j)
+	    call pargstr (IM_TITLE(im))
+
+	call gt_sets (gt, GTTITLE, title)
+
+	switch (axis) {
+	case 1:
+	    call ic_pstr (ic, "xlabel", "Column")
+	    call xt_21imavg (im, axis, 1, IM_LEN(im, 1), i, j, x, data, imlen)
+	case 2:
+	    call ic_pstr (ic, "xlabel", "Line")
+	    call xt_21imavg (im, axis, i, j, 1, IM_LEN(im, 2), x, data, imlen)
+	}
+	call mfree (x, TY_REAL)
+
+	stat = OK
+	return (stat)
+end