Initial commit

1 files changed, 293 insertions, 0 deletions

diff --git a/noao/imred/vtel/unwrap.x b/noao/imred/vtel/unwrap.x
new file mode 100644
index 00000000..a753ddf4
--- /dev/null
+++ b/noao/imred/vtel/unwrap.x
@@ -0,0 +1,293 @@
+include <mach.h>
+include <imhdr.h>
+
+define	MAXBADLINES	20		# maximum number of bad lines
+define	BADTHRESH	1000		# threshold for bad lines
+define	FIXWIDTH	20		# Width of average for fixline
+
+# UNWRAP -- Filter an iraf image.  This filter checks for binary wraparound
+# in IRAF images.  The algorithm is described in detail in the help page.
+# The program accepts templates for both input and output image lists.
+
+procedure t_unwrap()
+
+char	image[SZ_FNAME]			# input image template
+char	outimage[SZ_FNAME]		# output image template
+int	threshold1			# threshold value for first unwrap
+int	threshold2			# threshold value for second unwrap
+int	wrapval1			# wrapvalue for first unwrap
+int	wrapval2			# wrapvalue for second unwrap
+int	cstart				# column to start on
+int	step				# number of steps to perform
+bool	verbose				# verbose flag
+
+int	i, j 
+int	listin, listout
+int	length, nlines
+int	badlines[MAXBADLINES]
+int	diff, nbad
+char	tempimage[SZ_FNAME]
+pointer	im, imout, lgp, lgp2, lpp, cck, sp
+
+bool	clgetb()
+int	imtopenp(), imtlen(), imtgetim(), clgeti()
+pointer	immap(), imgl2s(), impl2s()
+errchk	immap, imgl2s, impl2s
+
+begin
+	# Get parameters from the CL.
+	listin = imtopenp ("image")
+	listout = imtopenp ("outimage")
+	threshold1 = clgeti("threshold1")
+	wrapval1 = clgeti("wrapval1")
+	threshold2 = clgeti("threshold2")
+	wrapval2 = clgeti("wrapval2")
+	cstart = clgeti("cstart")
+	step = clgeti("step")
+	verbose = clgetb("verbose")
+
+	if (verbose) {
+	    call printf ("\n\nUNWRAP: ")
+	    call printf ("threshold1 = %d\n")
+	    call pargi (threshold1)
+	    call printf ("\twrapval1   = %d\n")
+	    call pargi (wrapval1)
+	    call printf ("\tthreshold2 = %d\n")
+	    call pargi (threshold2)
+	    call printf ("\twrapval2   = %d\n")
+	    call pargi (wrapval2)
+	    call printf ("\tcstart     = %d\n")
+	    call pargi (cstart)
+	    call printf ("\tstep       = %d\n\n")
+	    call pargi (step)
+	    call flush (STDOUT)
+	}
+
+	# Check the number of elements.
+	if (imtlen (listin) != imtlen (listout)) {
+	    call imtclose (listin)
+	    call imtclose (listout)
+	    call error (1, "Wrong number of elements in the operand lists")
+	}
+
+	# Get the next images from the lists.
+	while (imtgetim (listin, image, SZ_FNAME) != EOF) {
+	    if (imtgetim (listout, outimage, SZ_FNAME) != EOF) {
+
+	        if (verbose) {
+		    # Write out about the input file name and output file name.
+		    call printf ("\tUnwrapping %s into %s.  ")
+		        call pargstr (image)
+		        call pargstr (outimage)
+		    call flush (STDOUT)
+	        }
+
+	        # Open images.
+	        iferr {
+	            im = immap (image, READ_WRITE, 0)
+	        } then {
+		    call eprintf ("Cannot open image %s.\n")
+		        call pargstr (image)
+		    next
+	        }
+
+	        call xt_mkimtemp (image, outimage, tempimage, SZ_FNAME)
+
+	        iferr {
+	            imout = immap (outimage, NEW_COPY, im)
+	        } then {
+		    call eprintf ("Cannot open image %s, (already exists?).\n")
+		        call pargstr (outimage)
+		    next
+	        }
+
+	        length = IM_LEN(im,1)
+	        nlines = IM_LEN(im,2)
+
+	        # Set up the column check array, then unwrap line by line.
+	        call smark (sp)
+	        call salloc (cck, nlines, TY_INT)
+	        call amovks (0, Memi[cck], nlines)
+	        do i = 1, nlines {
+	            lgp = imgl2s (im, i)
+	            lpp = impl2s (imout, i)
+	            call unwrapline (Mems[lgp], Mems[lpp], cck, length,
+			threshold1, wrapval1, threshold2, wrapval2, cstart,
+			step, i)
+	        }
+
+	        # Step 5 is the final step. (fixline)
+	        if (step == 5) {
+	            # Analyze the column, check for wraps.
+	            nbad = 0
+	            do i = 2, nlines {
+	    	        diff = Memi[cck+i-1] - Memi[cck+i-2]
+	    	        if (abs(diff) > BADTHRESH) {
+	    	            # Mark this line bad.
+	    	            nbad = nbad + 1
+	    	            if (nbad > MAXBADLINES)
+	    		        break
+	    	            badlines[nbad] = i
+	    	        }
+	            }
+	        }
+
+	        # If number bad lines <= than MAXBADLINES, fix em, else, quit.
+	        if (nbad <= MAXBADLINES && nbad > 0) {
+	            do i = 1, nbad {
+
+	    	        # GET the lines above and below the bad line and PUT the
+	    	        # bad line.  Then average the above and below lines and
+			# save in the bad line.
+
+	    	        if (badlines[i] != 1 && badlines[i] != nlines) {
+	    	            if ((badlines[i+1] - badlines[i]) == 1) {
+	                        lgp = imgl2s (imout, badlines[i]-1)
+	                        lgp2 = imgl2s (imout, badlines[i]+1)
+	                        lpp = impl2s (imout, badlines[i])
+	    	                do j = 1, length {
+	    		            Mems[lpp+j-1] = int((real(Mems[lgp+j-1]) +
+	    		                real(Mems[lgp2+j-1]))/2. + .5)
+	    		        }
+	    	            } 
+	    	        }
+	            }
+	        }
+
+	        if (verbose) {
+	            call printf ("number of bad lines = %d\n")
+	                call pargi (nbad)
+	            do i = 1, nbad {
+	                call printf ("\tbadlines[%d] = %d\n")
+	    	        call pargi (i)
+	    	        call pargi (badlines[i])
+	            }
+		    call printf ("\n")
+	            call flush (STDOUT)
+	        }
+
+	        # Unmap images.
+	        call imunmap (im)
+	        call imunmap (imout)
+	        call xt_delimtemp (outimage, tempimage)
+	        call sfree (sp)
+
+	    } # End of if (not EOF)
+	}  # End of while loop on input images
+
+	call imtclose (listin)
+	call imtclose (listout)
+end
+
+
+# UNWRAPLINE -- Unwrap a line of the image.
+
+procedure unwrapline (line1, line2, cck, numpix, threshold1, wrapval1,
+	threshold2, wrapval2, cstart, step, whichline)
+
+short	line1[numpix]			# input line
+short	line2[numpix]			# output line
+pointer	cck				# pointer to array for column check
+int	numpix				# number of pixels per line
+int	threshold1			# unwrap threshold for first unwrap
+int	wrapval1			# unwrap value for first unwrap
+int	threshold2			# unwrap threshold for second unwrap
+int	wrapval2			# unwrap value for second unwrap
+int	cstart				# column to start on
+int	step				# steps to complete
+int	whichline			# which line this is we are unwrapping
+
+pointer	tl1, tl2, tl3			# pointers of temporary arrays
+pointer	sp				# stack pointer
+int	i, diff, sum
+short	wrap				# wrap number
+
+begin
+	# Mark the stack and allcoate the temporary arrays.
+	call smark (sp)
+	call salloc (tl1, numpix, TY_SHORT)
+	call salloc (tl2, numpix, TY_SHORT)
+	call salloc (tl3, numpix, TY_SHORT)
+
+	# Initialize wrap.
+	wrap = 0
+
+	# Copy the input line into the output line and the temporary arrays.
+	call amovs (line1, line2, numpix)
+	call amovs (line1, Mems[tl1], numpix)
+	call amovs (line1, Mems[tl2], numpix)
+	call amovs (line1, Mems[tl3], numpix)
+
+	# Check the image width, do various things if the image is too small.
+
+	# Too small for anything.
+	if (numpix <= 4) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Too small for step 5 (fixline).
+	if (numpix <= FIXWIDTH && step == 5)
+	    step = 4
+
+	# Unwrap1 (step 1).
+	Mems[tl1+cstart-1] = line1[cstart]
+	do i = cstart+1, numpix {
+	    diff = line1[i] - line1[i-1]
+	    if (diff < -threshold1)
+		wrap = wrap + 1
+	    if (diff > threshold1)
+		wrap = wrap - 1
+
+	    Mems[tl1+i-1] = line1[i] + wrap * wrapval1
+	}
+	if (step == 1) {
+	    call amovs (Mems[tl1], line2, numpix)
+	    call sfree (sp)
+	    return
+	}
+
+	# If the user wants it, step 2 (dif).
+	do i = cstart, numpix
+	    Mems[tl2+i-1] = Mems[tl1+i-1] - Mems[tl1+i-2]
+
+	if (step == 2) {
+	    call amovs (Mems[tl2], line2, numpix)
+	    call sfree (sp)
+	    return
+	}
+
+	# If the user wants it, step 3 (unwrap2).
+	wrap = 0
+	line2[cstart] = Mems[tl2+cstart-1]
+	do i = cstart+1, numpix {
+	    diff = Mems[tl2+i-1] - Mems[tl2+i-2]
+	    if (diff < -threshold2)
+	        wrap = wrap + 1
+	    if (diff > threshold2)
+		wrap = wrap - 1
+
+	    line2[i] = Mems[tl2+i-1] + wrap * wrapval2
+	}
+	if (step == 3) {
+	    call sfree (sp)
+	    return
+	}
+
+	# If the user wants it, step 4 (reconstruct).
+	do i = cstart, numpix
+	    line2[i] = line2[i-1] + line2[i]
+	
+	if (step == 4) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Again, if the user wants it, save data for step 5, (fixline).
+	sum = 0
+	do i = numpix-FIXWIDTH+1, numpix
+	    sum = sum + line2[i]
+	Memi[cck+whichline-1] = int(real(sum)/real(FIXWIDTH) + .5)
+
+	call sfree (sp)
+end