1 files changed, 132 insertions, 0 deletions
diff --git a/pkg/images/imutil/src/t_imdivide.x b/pkg/images/imutil/src/t_imdivide.x
new file mode 100644
index 00000000..510e49e5
--- /dev/null
+++ b/pkg/images/imutil/src/t_imdivide.x
@@ -0,0 +1,132 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	<imhdr.h>
+
+# T_IMDIVIDE -- Image division with rescaling.
+
+# Options for rescaling.
+define	NORESC	1		# Do not scale resultant image
+define	MEAN	2		# Scale resultant mean to given value
+define	NUMER	3		# Scale resultant mean to mean of numerator
+
+procedure t_imdivide ()
+
+char	image1[SZ_FNAME]			# Numerator image
+char	image2[SZ_FNAME]			# Denominator image
+char	image3[SZ_FNAME]			# Resultant image
+char	title[SZ_IMTITLE]			# Resultant image title
+int	rescale					# Option for rescaling
+real	constant				# Replacement for zero divide
+bool	verbose					# Verbose output?
+
+char	str[SZ_LINE]
+int	i, npix, ntotal
+real	sum1, sum2, sum3, scale
+long	line1[IM_MAXDIM], line2[IM_MAXDIM], line3[IM_MAXDIM]
+pointer	im1, im2, im3, data1, data2, data3
+
+int	clgwrd(), imgnlr(), impnlr()
+bool	clgetb(), strne()
+real	clgetr(), asumr(), ima_efncr()
+pointer	immap()
+extern	ima_efncr
+
+common	/imadcomr/ constant
+
+begin
+	# Access images and set parameters.
+	call clgstr ("numerator", image1, SZ_FNAME)
+	im1 = immap (image1, READ_ONLY, 0)
+	call clgstr ("denominator", image2, SZ_FNAME)
+	im2 = immap (image2, READ_ONLY, 0)
+	call clgstr ("resultant", image3, SZ_FNAME)
+	im3 = immap (image3, NEW_COPY, im1)
+
+	if (IM_NDIM (im1) != IM_NDIM (im2))
+	    call error (0, "Input images have different dimensions")
+	do i = 1, IM_NDIM (im1)
+	    if (IM_LEN (im1, i) != IM_LEN (im2, i))
+		call error (0, "Input images have different sizes")
+
+	call clgstr ("title", title, SZ_IMTITLE)
+	if (strne (title, "*"))
+	    call strcpy (title, IM_TITLE(im3), SZ_IMTITLE) 
+	IM_PIXTYPE(im3) = TY_REAL
+
+	constant = clgetr ("constant")
+	verbose = clgetb ("verbose")
+
+	# Initialize.
+	npix = IM_LEN(im1, 1)
+	ntotal = 0
+	sum1 = 0.
+	sum2 = 0.
+	sum3 = 0.
+	call amovkl (long(1), line1, IM_MAXDIM)
+	call amovkl (long(1), line2, IM_MAXDIM)
+	call amovkl (long(1), line3, IM_MAXDIM)
+
+	# Loop through the images doing the division.
+	# Accumulate the sums for mean values.
+	while (impnlr (im3, data3, line3) != EOF) {
+	    i = imgnlr (im1, data1, line1)
+	    i = imgnlr (im2, data2, line2)
+	    call advzr (Memr[data1], Memr[data2], Memr[data3], npix, ima_efncr)
+	    sum1 = sum1 + asumr (Memr[data1], npix)
+	    sum2 = sum2 + asumr (Memr[data2], npix)
+	    sum3 = sum3 + asumr (Memr[data3], npix)
+	    ntotal = ntotal + npix
+	}
+	sum1 = sum1 / ntotal
+	sum2 = sum2 / ntotal
+	sum3 = sum3 / ntotal
+
+	# Close the images.
+	call imunmap (im1)
+	call imunmap (im2)
+	call imunmap (im3)
+
+	# Print image means if verbose.
+	if (verbose) {
+	    call printf ("Task imdivide:\n")
+	    call printf ("    %s: Mean = %g\n")
+		call pargstr (image1)
+		call pargr (sum1)
+	    call printf ("    %s: Mean = %g\n")
+		call pargstr (image2)
+		call pargr (sum2)
+	    call printf ("    %s: Mean = %g\n")
+		call pargstr (image3)
+		call pargr (sum3)
+	}
+
+	# Determine resultant image rescaling.
+	rescale = clgwrd ("rescale", str, SZ_LINE, ",norescale,mean,numerator,")
+	switch (rescale) {
+	case NORESC:
+	    return
+	case MEAN:
+	    scale = clgetr ("mean") / sum3
+	case NUMER:
+	    scale = sum1 / sum3
+	}
+
+	if(verbose) {
+	    call printf ("    %s: Scale = %g\n")
+		call pargstr (image3)
+		call pargr (scale)
+	}
+
+	# Open image read_write and initialize line counters.
+	im1 = immap (image3, READ_WRITE, 0)
+	call amovkl (long(1), line1, IM_MAXDIM)
+	call amovkl (long(1), line2, IM_MAXDIM)
+
+	# Loop through the image rescaling the image lines.
+	while (imgnlr (im1, data1, line1) != EOF) {
+	    i = impnlr (im1, data2, line2)
+	    call amulkr (Memr[data1], scale, Memr[data2], npix)
+	}
+
+	call imunmap (im1)
+end