Initial commit

1 files changed, 619 insertions, 0 deletions

diff --git a/pkg/images/imutil/src/t_imtile.x b/pkg/images/imutil/src/t_imtile.x
new file mode 100644
index 00000000..92f5cce0
--- /dev/null
+++ b/pkg/images/imutil/src/t_imtile.x
@@ -0,0 +1,619 @@
+include <imhdr.h>
+include <fset.h>
+include "imtile.h"
+
+
+# T_IMTILE -- Combine a list of same-size subrasters into a single large
+# mosaiced image.
+
+procedure t_imtile ()
+
+int	nimages, nmissing, subtract, verbose
+pointer	it, sp, outimage, trimsection, medsection, nullinput, ranges
+pointer	str, index, c1, c2, l1, l2, isnull, median, imlist, outim
+
+bool	clgetb()
+char	clgetc()
+int	btoi(), clgwrd(), imtlen(), clgeti(), decode_ranges(), it_get_imtype()
+pointer	imtopenp(), it_setim()
+real	clgetr()
+
+begin
+	call fseti (STDOUT, F_FLUSHNL, YES)
+	call malloc (it, LEN_IRSTRUCT, TY_STRUCT)
+
+	# Allocate temporary working space.
+	call smark (sp)
+	call salloc (outimage, SZ_FNAME, TY_CHAR)
+	call salloc (trimsection, SZ_FNAME, TY_CHAR)
+	call salloc (medsection, SZ_FNAME, TY_CHAR)
+	call salloc (nullinput, SZ_FNAME, TY_CHAR)
+	call salloc (ranges, 3 * MAX_NRANGES + 1, TY_INT)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the input image list and the output image name.
+	imlist = imtopenp ("input")
+	call clgstr ("output", Memc[outimage], SZ_FNAME)
+	call clgstr ("trim_section", Memc[trimsection], SZ_FNAME)
+	call clgstr ("missing_input", Memc[nullinput], SZ_FNAME)
+	call clgstr ("median_section", Memc[medsection], SZ_FNAME)
+	if (Memc[medsection] == EOS)
+	    subtract = NO
+	else
+	    subtract = btoi (clgetb ("subtract"))
+	verbose = btoi (clgetb ("verbose"))
+
+	# Get the mosaicing parameters.
+	IT_NXSUB(it) = clgeti ("nctile")
+	IT_NYSUB(it) = clgeti ("nltile")
+	IT_CORNER(it) = clgwrd ("start_tile", Memc[str], SZ_FNAME,
+	    ",ll,lr,ul,ur,")
+	if (clgetb ("row_order"))
+	    IT_ORDER(it) = IT_ROW
+	else
+	    IT_ORDER(it) = IT_COLUMN
+	IT_RASTER(it) = btoi (clgetb ("raster_order"))
+	IT_NXOVERLAP(it) = clgeti ("ncoverlap")
+	IT_NYOVERLAP(it) = clgeti ("nloverlap")
+	IT_OVAL(it) = clgetr ("ovalue")
+
+	# Check that the number of observed and missing images matches
+	# the number of specified subrasters.
+	if (Memc[nullinput] == EOS) {
+	    nmissing = 0
+	    Memi[ranges] = 0
+	    Memi[ranges+1] = 0
+	    Memi[ranges+2] = 1
+	    Memi[ranges+3] = NULL
+	} else {
+	    if (decode_ranges (Memc[nullinput], Memi[ranges], MAX_NRANGES,
+	        nmissing) == ERR)
+	        call error (0, "Error decoding list of unobserved rasters.")
+	}
+	nimages = imtlen (imlist) + nmissing
+	if (nimages != (IT_NXSUB(it) * IT_NYSUB(it)))
+	    call error (0,
+	        "The number of input images is not equal to nxsub * nysub.")
+
+	# Compute the output image characteristics and open the output image.
+	outim = it_setim (it, imlist, Memc[trimsection], Memc[outimage],
+	    clgeti ("ncols"), clgeti ("nlines"), it_get_imtype (clgetc (
+	    "opixtype")))
+
+	# Allocate space for and setup the section descriptors.
+	call salloc (index, nimages, TY_INT)
+	call salloc (c1, nimages, TY_INT)
+	call salloc (c2, nimages, TY_INT)
+	call salloc (l1, nimages, TY_INT)
+	call salloc (l2, nimages, TY_INT)
+	call salloc (isnull, nimages, TY_INT)
+	call salloc (median, nimages, TY_REAL)
+
+	call it_setup (it, imlist, Memi[ranges], Memc[trimsection],
+	    Memc[medsection], outim, Memi[index], Memi[c1], Memi[c2],
+	    Memi[l1], Memi[l2], Memi[isnull], Memr[median])
+
+	# Make the output image.
+	call it_mkmosaic (imlist, Memc[trimsection], outim, Memi[index],
+	    Memi[c1], Memi[c2], Memi[l1], Memi[l2], Memi[isnull],
+	    Memr[median], IT_NXSUB(it), IT_NYSUB(it), IT_OVAL(it), subtract)
+
+        # Printe the results.
+	if (verbose == YES) {
+            call it_show (imlist, Memc[trimsection], Memc[outimage],
+		Memi[index], Memi[c1], Memi[c2], Memi[l1], Memi[l2],
+		Memi[isnull], Memr[median], IT_NXSUB(it)*IT_NYSUB(it), subtract)
+	}
+
+	# Close up files and free space.
+	call imunmap (outim)
+	call clpcls (imlist)
+	call sfree (sp)
+	call mfree (it, TY_STRUCT)
+end
+
+
+define	NTYPES	7
+
+# IT_GET_IMTYPE -- Procedure to get the image type.
+
+int procedure it_get_imtype (c)
+
+char	c	# character denoting the image type
+
+int	i, typecodes[NTYPES]
+int	stridx()
+string	types "usilrdx"
+data	typecodes /TY_USHORT, TY_SHORT, TY_INT, TY_LONG, TY_REAL, TY_DOUBLE,
+		   TY_COMPLEX/
+
+begin
+	i = stridx (c, types)
+	if (i == 0)
+	    return (ERR)
+	else
+	    return (typecodes[i])
+end
+
+
+# IT_SETUP -- Setup the data base parameters for the images.
+
+procedure it_setup (it, imlist, ranges, trimsection, medsection, outim,
+	index, c1, c2, l1, l2, isnull, median)
+
+pointer	it			# pointer to the imtil structure
+pointer	imlist			# pointer to the list of input images
+int	ranges[ARB]		# list of missing subrasters
+char	trimsection[ARB]	# input image section for output
+char	medsection[ARB]		# input image section for median computation
+pointer	outim			# pointer to the output image
+int	index[ARB]		# index array
+int	c1[ARB]			# array of beginning column limits
+int	c2[ARB]			# array of ending column limits
+int	l1[ARB]			# array of beginning line limits
+int	l2[ARB]			# array of ending line limits
+int	isnull[ARB]		# output input image order number
+real	median[ARB]		# output median of input image
+
+int	i, j, k, nimrows, nimcols, imcount, next_null
+pointer	sp, imname, im, buf
+int	get_next_number(), imtgetim()
+pointer	immap(), imgs2r()
+real	amedr()
+
+begin
+	nimcols = IM_LEN(outim,1)
+	nimrows = IM_LEN(outim,2)
+
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+
+	imcount = 1
+	next_null = 0
+	if (get_next_number (ranges, next_null) == EOF)
+	    next_null = IT_NXSUB(it) * IT_NYSUB(it) + 1
+
+	# Loop over the input images.
+	do i = 1, IT_NXSUB(it) * IT_NYSUB(it) {
+
+	    # Set the indices array.
+	    call it_indices (i, j, k, IT_NXSUB(it), IT_NYSUB(it),
+	        IT_CORNER(it), IT_RASTER(it), IT_ORDER(it))
+	    index[i] = i
+	    c1[i] = max (1, min (1 + (j - 1) * (IT_NCOLS(it) -
+	        IT_NXOVERLAP(it)), nimcols))
+	    c2[i] = min (nimcols, max (1, c1[i] + IT_NCOLS(it) - 1))
+	    l1[i] = max (1, min (1 + (k - 1) * (IT_NROWS(it) -
+	        IT_NYOVERLAP(it)), nimrows))
+	    l2[i] = min (nimrows, max (1, l1[i] + IT_NROWS(it) - 1))
+
+	    # Set the index of each image in the image template
+	    # and compute the median of the subraster.
+	    if (i < next_null) {
+		isnull[i] = imcount
+		if (medsection[1] != EOS) {
+		    if (imtgetim (imlist, Memc[imname], SZ_FNAME) == EOF)
+			call error (0, "Error reading input image list.")
+		    call strcat (medsection, Memc[imname], SZ_FNAME)
+		    im = immap (Memc[imname], READ_ONLY, TY_CHAR)
+		    buf = imgs2r (im, 1, int (IM_LEN(im,1)), 1, int (IM_LEN(im,
+		        2)))
+		    median[i] = amedr (Memr[buf], int (IM_LEN(im,1)) *
+			    int (IM_LEN(im,2)))
+		    call imunmap (im)
+		} else
+		    median[i] = INDEFR
+		imcount = imcount + 1
+	    } else {
+		isnull[i] = 0
+		if (medsection[1] == EOS)
+		    median[i] = INDEFR
+		else
+		    median[i] = IT_OVAL(it)
+		if (get_next_number (ranges, next_null) == EOF)
+	    	    next_null = IT_NXSUB(it) * IT_NYSUB(it) + 1
+	    }
+
+	}
+
+	call imtrew (imlist)
+	call sfree (sp)
+end
+
+
+# IT_SETIM -- Procedure to set up the output image characteristics.
+
+pointer procedure it_setim (it, list, trimsection, outimage, nimcols, nimrows,
+	opixtype)
+
+pointer	it		# pointer to the imtile structure
+pointer	list		# pointer to list of input images
+char	trimsection[ARB]# input image section
+char	outimage[ARB]	# name of the output image
+int	nimcols		# number of output image columns
+int	nimrows		# number of output image rows
+int	opixtype	# output image pixel type
+
+int	ijunk, nc, nr
+pointer	sp, imname, im, outim
+int	imtgetim()
+pointer	immap()
+
+begin
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+
+	# Get the size of the first subraster.
+	if (imtgetim (list, Memc[imname], SZ_FNAME) != EOF) {
+	    call strcat (trimsection, Memc[imname], SZ_FNAME)
+	    im = immap (Memc[imname], READ_ONLY, 0)
+	    IT_NCOLS(it) = IM_LEN(im,1)
+	    IT_NROWS(it) = IM_LEN(im,2)
+	    call imunmap (im)
+	    call imtrew (list)
+	} else
+	    call error (0, "Error reading first input image.\n")
+
+	# Compute the size of the output image.
+	ijunk = IT_NXSUB(it) * IT_NCOLS(it) - (IT_NXSUB(it) - 1) *
+	    IT_NXOVERLAP(it)
+	if (IS_INDEFI(nimcols))
+	    nc = ijunk
+	else
+	    nc = max (nimcols, ijunk)
+	ijunk = IT_NYSUB(it) * IT_NROWS(it) - (IT_NYSUB(it) - 1) *
+	    IT_NYOVERLAP(it)
+	if (IS_INDEFI(ijunk))
+	    nr = ijunk
+	else
+	    nr = max (nimrows, ijunk)
+
+	# Set the output pixel type.
+	if (opixtype == ERR)
+	    opixtype = TY_REAL
+
+	# Open output image and set the parameters.
+	outim = immap (outimage, NEW_IMAGE, 0)
+	IM_NDIM(outim) = 2
+	IM_LEN(outim,1) = nc
+	IM_LEN(outim,2) = nr
+	IM_PIXTYPE(outim) = opixtype
+
+	call sfree (sp)
+
+	return (outim)
+end
+
+
+# IT_MKMOSAIC -- Procedure to make the mosaiced image.
+
+procedure it_mkmosaic (imlist, trimsection, outim, index, c1, c2, l1, l2,
+	isnull, median, nxsub, nysub, oval, subtract)
+
+pointer	imlist		# pointer to input image list
+char	trimsection[ARB]# input image section
+pointer	outim		# pointer to the output image
+int	index[ARB]	# index array for sorting the images
+int	c1[ARB]		# array of column beginnings
+int	c2[ARB]		# array of column endings
+int	l1[ARB]		# array of line beginnings
+int	l2[ARB]		# array of line endings
+int	isnull[ARB]	# index of input image in the template
+real	median[ARB]	# array of input image median values
+int	nxsub		# number of subrasters per output image column
+int	nysub		# number of subrasters per output image row
+real	oval		# pixel value of undefined output image regions
+int	subtract	# subtract the median off each subraster
+
+int	i, j, noutcols, noutlines, olineptr, ll1, ll2
+pointer	sp, inimage, imptrs, buf
+int	imtrgetim()
+pointer	immap(), impl2r()
+
+begin
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (imptrs, nxsub, TY_POINTER)
+	call salloc (inimage, SZ_FNAME, TY_CHAR)
+
+	# Sort the subrasters on the yindex.
+	do i = 1, nxsub * nysub
+	    index[i] = i
+	call rg_qsorti (l1, index, index, nxsub * nysub)
+
+	noutcols = IM_LEN(outim,1)
+	noutlines = IM_LEN(outim,2)
+
+	# Loop over the input images.
+	olineptr = 1
+	do i = 1, nxsub * nysub, nxsub {
+
+	    # Compute the line and column limits.
+	    ll1 = l1[index[i]]
+	    ll2 = l2[index[i]]
+
+	    # Open the nxsub input images.
+	    do j = i, i + nxsub - 1 {
+		if (isnull[index[j]] <= 0) {
+		    Memc[inimage] = EOS
+		    Memi[imptrs+j-i] = NULL
+		} else {
+		    if (imtrgetim (imlist, isnull[index[j]], Memc[inimage],
+		        SZ_FNAME) == EOF)
+			Memi[imptrs+j-i] = NULL
+		    else {
+			call strcat (trimsection, Memc[inimage], SZ_FNAME)
+			Memi[imptrs+j-i] = immap (Memc[inimage], READ_ONLY, 0)
+		    }
+		}
+	    }
+
+	    # Write out the undefined lines.
+	    while (olineptr < ll1) {
+		buf = impl2r (outim, olineptr)
+		call amovkr (oval, Memr[buf], noutcols)
+		olineptr = olineptr + 1
+	    }
+
+	    # Write the output lines.
+	    call it_mklines (Memi[imptrs], outim, index, c1, c2, ll1, ll2,
+	        median, i, nxsub, oval, subtract)
+	    olineptr = ll2 + 1
+
+	    # Close up the images.
+	    # Open the nxsub input images.
+	    do j = i, i + nxsub - 1 {
+		if (Memi[imptrs+j-i] != NULL)
+		    call imunmap (Memi[imptrs+j-i])
+	    }
+
+	}
+
+	# Write out the remaining undefined lines.
+	while (olineptr < noutlines) {
+	    buf = impl2r (outim, olineptr)
+	    call amovkr (oval, Memr[buf], noutcols)
+	    olineptr = olineptr + 1
+	}
+
+	call sfree (sp)
+end
+
+
+# IT_MKLINES -- Construct and output image lines.
+
+procedure it_mklines (imptrs, outim, index, c1, c2, l1, l2, meds, init, nsub, 
+	oval, subtract)
+
+pointer	imptrs[ARB]		# array of input image pointers
+pointer	outim			# output imnage pointer
+int	index[ARB]		# array of indices
+int	c1[ARB]			# array of beginning columns
+int	c2[ARB]			# array of ending columns
+int	l1			# beginning line
+int	l2			# ending line
+real	meds[ARB]		# array of median values
+int	init			# first index
+int	nsub			# number of subrasters
+real	oval			# output value
+int	subtract		# subtract the median value
+
+int	i, j, jj, noutcols
+pointer	obuf, ibuf
+pointer	impl2r(), imgl2r()
+
+begin
+	noutcols = IM_LEN(outim, 1)
+	do i = l1, l2 {
+	    obuf = impl2r (outim, i)
+	    call amovkr (oval, Memr[obuf],  noutcols)
+	    do j = 1, nsub {
+		jj = index[j+init-1]
+		if (imptrs[j] != NULL) {
+		    ibuf = imgl2r (imptrs[j], i - l1 + 1)
+		    if (subtract == YES)
+		        call asubkr (Memr[ibuf], meds[jj], Memr[obuf+c1[jj]-1],
+			    c2[jj] - c1[jj] + 1)
+		    else
+		        call amovr (Memr[ibuf], Memr[obuf+c1[jj]-1], c2[jj] -
+			    c1[jj] + 1)
+		}
+	    }
+	}
+end
+
+
+# IT_INDICES -- Given the number in the list for a missing subraster and
+# information about how the subrasters were written return the i and j
+# indices of the specified subrasters.
+
+procedure it_indices (num, i, j, nxsub, nysub, corner, raster, order)
+
+int	num		# number of the subraster
+int	i,j		# indices of the subraster
+int	nxsub,nysub	# number of subrasters in x and y
+int	corner		# starting corner
+int	raster		# raster order
+int	order		# column or row order
+
+begin
+	switch (corner) {
+	case IT_LL:
+	    if (order == IT_ROW) {
+		if (mod (num, nxsub) == 0) {
+		    j = num / nxsub
+		    if (raster == YES && mod (j,2) == 0)
+			i = 1
+		    else
+		        i = nxsub
+		} else {
+		    j = num / nxsub + 1
+		    if (raster == YES && mod (j,2) == 0)
+			i = nxsub - mod (num, nxsub) + 1
+		    else
+		        i = mod (num, nxsub)
+		}
+	    } else if (order == IT_COLUMN) {
+		if (mod (num, nysub) == 0) {
+		    i = num / nysub
+		    if (raster == YES && mod (i,2) == 0)
+			j = 1
+		    else
+		        j = nysub
+		} else {
+		    i = num / nysub + 1
+		    if (raster == YES && mod (i,2) == 0)
+			j = nysub - mod (num, nysub) + 1
+		    else
+		        j = mod (num, nysub)
+		}
+	    }
+	case IT_LR:
+	    if (order == IT_ROW) {
+		if (mod (num, nxsub) == 0) {
+		    j = num / nxsub
+		    if (raster == YES && mod (j,2) == 0)
+			i = nxsub
+		    else
+			i = 1
+		} else {
+		    j = num / nxsub + 1
+		    if (raster == YES && mod (j,2) == 0)
+			i = mod (num, nxsub)
+		    else
+			i = nxsub - mod (num, nxsub) + 1
+		}
+	    } else if (order == IT_COLUMN) {
+		if (mod (num, nysub) == 0) {
+		    i = nxsub - num / nysub + 1
+		    if (raster == YES && mod (i,2) != 0)
+			j = 1
+		    else
+		        j = nysub
+		} else {
+		    i = nxsub - num / nysub
+		    if (raster == YES && mod (i,2) != 0)
+			j = nysub - mod (num, nysub) + 1
+		    else
+		        j = mod (num, nysub)
+		}
+	    }
+	case IT_UL:
+	    if (order == IT_ROW) {
+		if (mod (num, nxsub) == 0) {
+		    j = nysub - num / nxsub + 1
+		    if (raster == YES && mod (j,2) != 0)
+			i = 1
+		    else
+		        i = nxsub
+		} else {
+		    j = nysub - num / nxsub
+		    if (raster == YES && mod (j,2) != 0)
+			i = nxsub - mod (num, nxsub) + 1
+		    else
+		        i = mod (num, nxsub)
+		}
+	    } else if (order == IT_COLUMN) {
+		if (mod (num, nysub) == 0) {
+		    i = num / nysub
+		    if (raster == YES && mod (i,2) == 0)
+			j = nysub
+		    else
+			j = 1
+		} else {
+		    i = num / nysub + 1
+		    if (raster == YES && mod (i,2) == 0)
+			j = mod (num, nysub)
+		    else
+			j = nysub - mod (num, nysub) + 1
+		}
+	    }
+	case IT_UR:
+	    if (order == IT_ROW) {
+		if (mod (num, nxsub) == 0) {
+		    j = nysub - num / nxsub + 1
+		    if (raster == YES && mod (j,2) != 0)
+			i = nxsub
+		    else
+			i = 1
+		} else {
+		    j = nysub - num / nxsub
+		    if (raster == YES && mod (j,2) != 0)
+			i = mod (num, nxsub)
+		    else
+			i = nxsub - mod (num, nxsub) + 1
+		}
+	    } else if (order == IT_COLUMN) {
+		if (mod (num, nysub) == 0) {
+		    i = nxsub - num / nysub + 1
+		    if (raster == YES && mod (i,2) != 0)
+			j = nysub
+		    else
+			j = 1
+		} else {
+		    i = nxsub - num / nysub
+		    if (raster == YES && mod (i,2) != 0)
+			j = mod (num, nysub)
+		    else
+			j = nysub - mod (num, nysub) + 1
+		}
+	    }
+	}
+end
+
+
+# IT_SHOW -- List the results.
+
+procedure it_show (imlist, trimsection, outimage, index, c1, c2, l1,
+        l2, isnull, median, nsub, subtract)
+
+int     imlist          # input image list
+char    trimsection[ARB]# trim section of input image
+char    outimage[ARB]   # output image
+int     index[ARB]      # array of sorted indices (not used at present)
+int     c1[ARB]         # array of beginning column limits
+int     c2[ARB]         # array of ending column limits
+int     l1[ARB]         # array of beginning line limits
+int     l2[ARB]         # array of ending line limits
+int     isnull[ARB]     # image name index
+real    median[ARB]     # array of medians
+int     nsub            # number of subrasters
+int     subtract        # subtract the median from the subraster
+
+int     i
+pointer sp, imname
+int     imtrgetim()
+
+begin
+        call smark (sp)
+        call salloc (imname, SZ_FNAME, TY_CHAR)
+
+        do i = 1, nsub {
+
+            if (isnull[i] <= 0)
+                call strcpy ("nullimage", Memc[imname], SZ_FNAME)
+            else if (imtrgetim (imlist, isnull[i], Memc[imname],
+                SZ_FNAME) != EOF)
+                call strcat (trimsection, Memc[imname], SZ_FNAME)
+            else
+                Memc[imname] = EOS
+
+            call printf ("imcopy  %s  %s[%d:%d,%d:%d]  %g  %g\n")
+                call pargstr (Memc[imname])
+                call pargstr (outimage)
+                call pargi (c1[i])
+                call pargi (c2[i])
+                call pargi (l1[i])
+                call pargi (l2[i])
+                call pargr (median[i])
+            if (subtract == YES)
+                call pargr (-median[i])
+            else
+                call pargr (0.0)
+        }
+
+        call sfree (sp)
+end
+
+
+