Initial commit

1 files changed, 552 insertions, 0 deletions

diff --git a/noao/artdata/t_mkobjects.x b/noao/artdata/t_mkobjects.x
new file mode 100644
index 00000000..c9ad0ce3
--- /dev/null
+++ b/noao/artdata/t_mkobjects.x
@@ -0,0 +1,552 @@
+include	<error.h>
+include	<imhdr.h>
+include	<math.h>
+include	<mach.h>
+
+define	LEN_UA		20000			# Maximum user header
+define	LEN_COMMENT	70			# Maximum comment length
+
+# Object data structure
+define	LEN_MKO		9
+define	MKO_MKT		Memi[$1]		# Template
+define	MKO_X		Memi[$1+1]		# X position
+define	MKO_Y		Memi[$1+2]		# Y position
+define	MKO_Z		Memi[$1+3]		# Flux
+define	MKO_R		Memi[$1+4]		# Scale size
+define	MKO_AR		Memi[$1+5]		# Axial ratio
+define	MKO_PA		Memi[$1+6]		# Position angle
+define	MKO_SAVE	Memi[$1+7]		# Save template?
+define	MKO_SORT	Memi[$1+8]		# Sort index
+
+
+# T_MKOBJECTS -- Add stars and galaxies to images.
+# New images may be created with a background and noise.
+
+procedure t_mkobjects ()
+
+int	ilist				# Input image list
+int	olist				# Output image list
+int	objects				# List of model files
+real	xo				# X offset
+real	yo				# Y offset
+int	nl				# Number of lines
+int	nc				# Number of columns
+real	background			# Background level
+real	gain				# Gain (electrons/DN)
+int	ranbuf				# Random number buffer size
+real	rdnoise				# Read noise (in electrons)
+bool	poisson				# Add Poisson noise?
+real	exptime				# Exposure time
+real	distance			# Relative distance
+real	m0				# Magnitude zero point
+long	seed				# Random number seed
+
+int	nobjects, save
+real	x, y, z, r, ar, pa, dmin, dmax
+
+bool	new, bsave, cmmts, fcmmts
+int	i, j, k, l, nx, ny, nlines, c1, c2, c3, c4, l1, l2, l3, l4, irbuf, ipbuf
+long	seed1
+pointer	sp, input, output, fname, type, star, comment, rbuf, pbuf
+pointer	in, out, buf, obuf, lines, newlines, obj, ptr1, ptr2
+pointer	mko, mkt
+
+long	clgetl(), clktime()
+bool	clgetb(), streq()
+int	imtopenp(), imtlen(), imtgetim(), btoi()
+int	clgeti(), open(), fscan(), nscan()
+real	clgetr()
+pointer	immap(), imgl2r(), impl2r()
+pointer	mkt_star(), mkt_object()
+errchk	open, immap, imgl2r, impl2r, malloc, realloc, mkt_star, mkt_object
+
+int	mko_compare()
+extern	mko_compare
+pointer	mko_sort
+common	/mko_qsort/ mko_sort
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+	call salloc (type, SZ_FNAME, TY_CHAR)
+	call salloc (star, SZ_FNAME, TY_CHAR)
+	call salloc (comment, LEN_COMMENT, TY_CHAR)
+
+	# Get parameters which apply to all images.
+	ilist = imtopenp ("input")
+	olist = imtopenp ("output")
+	objects = imtopenp ("objects")
+	xo = clgetr ("xoffset")
+	yo = clgetr ("yoffset")
+	call clgstr ("star", Memc[star], SZ_FNAME)
+	distance = clgetr ("distance")
+	background = clgetr ("background")
+	gain = clgetr ("gain")
+	ranbuf = clgeti ("ranbuf")
+	if (ranbuf == 0)
+	    ranbuf = -1
+	rdnoise = clgetr ("rdnoise") / gain
+	if (rdnoise > 0. && ranbuf > 0)
+	    call salloc (rbuf, ranbuf, TY_REAL)
+	poisson = clgetb ("poisson")
+	if (poisson && ranbuf > 0)
+	    call salloc (pbuf, ranbuf, TY_REAL)
+	exptime = clgetr ("exptime")
+	m0 = clgetr ("magzero")
+	seed = clgetl ("seed")
+	if (IS_INDEFL(seed))
+	    seed1 = seed1 + clktime (long (0))
+	else
+	    seed1 = seed
+	cmmts = clgetb ("comments")
+
+	background = exptime * background
+
+	if (imtlen (objects) != imtlen (ilist))
+	    call error (1, "Input and objects lists don't match")
+
+	# Initialize the template library.
+	call mkt_init ()
+
+	# Loop through input, output, and object lists.
+	# Missing output images take the input image name.
+	# The object list will repeat if shorter than input list.
+
+	Memc[fname] = EOS
+	while (imtgetim (ilist, Memc[input], SZ_FNAME) != EOF) {
+	    if (imtgetim (olist, Memc[output], SZ_FNAME) == EOF)
+		call strcpy (Memc[input], Memc[output], SZ_FNAME)
+
+	    # Get and check object list.
+	    i = imtgetim (objects, Memc[fname], SZ_FNAME)
+	    iferr (i = open (Memc[fname], READ_ONLY, TEXT_FILE)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Map images.  Check for new, existing, and in-place images.
+	    if (streq (Memc[input], Memc[output])) {
+		ifnoerr (out = immap (Memc[output], READ_WRITE, LEN_UA)) {
+		    in = out
+		    new = false
+		} else {
+		    iferr (out = immap (Memc[output], NEW_IMAGE, LEN_UA)) {
+			call erract (EA_WARN)
+			next
+		    }
+
+		    call clgstr ("header", Memc[comment], LEN_COMMENT)
+		    iferr (call mkh_header (out, Memc[comment], false, false))
+			call erract (EA_WARN)
+
+		    IM_NDIM(out) = 2
+		    IM_LEN(out,1) = clgeti ("ncols")
+		    IM_LEN(out,2) = clgeti ("nlines")
+		    IM_PIXTYPE(out) = TY_REAL
+		    call clgstr ("title", IM_TITLE(out), SZ_IMTITLE)
+		    call imaddr (out, "exptime", exptime)
+		    call imaddr (out, "gain", gain)
+		    call imaddr (out, "rdnoise", rdnoise * gain)
+
+		    in = out
+		    new = true
+		}
+	    } else {
+		iferr (in = immap (Memc[input], READ_ONLY, LEN_UA)) {
+		    call erract (EA_WARN)
+		    next
+		}
+		iferr (out = immap (Memc[output], NEW_COPY, in)) {
+		    call erract (EA_WARN)
+		    call imunmap (in)
+		    next
+		}
+		new = false
+	    }
+	    nc = IM_LEN(in,1)
+	    nl = IM_LEN(in,2)
+	    IM_MIN(out) = MAX_REAL
+	    IM_MAX(out) = -MAX_REAL
+
+	    # Set star and seeing templates.
+	    mkt = mkt_star (Memc[star])
+
+	    # Read the object list.
+	    call malloc (mko, LEN_MKO, TY_STRUCT)
+	    fcmmts = false
+	    nobjects = 0
+	    while (fscan (i) != EOF) {
+		call gargr (x)
+		call gargr (y)
+		call gargr (z)
+		if (nscan() < 3) {
+		    fcmmts = true
+		    next
+		}
+		call gargwrd (Memc[type], SZ_FNAME)
+		call gargr (r)
+		call gargr (ar)
+		call gargr (pa)
+		call gargb (bsave)
+		x = xo + x / distance
+		y = yo + y / distance
+		if (x < 1 || x > nc || y < 1 || y > nl)
+		    next
+		if (nobjects == 0) {
+		    j = 100
+		    call malloc (MKO_MKT(mko), j, TY_POINTER)
+		    call malloc (MKO_X(mko), j, TY_REAL)
+		    call malloc (MKO_Y(mko), j, TY_REAL)
+		    call malloc (MKO_Z(mko), j, TY_REAL)
+		    call malloc (MKO_R(mko), j, TY_REAL)
+		    call malloc (MKO_AR(mko), j, TY_REAL)
+		    call malloc (MKO_PA(mko), j, TY_REAL)
+		    call malloc (MKO_SAVE(mko), j, TY_INT)
+		    call malloc (MKO_SORT(mko), j, TY_INT)
+		} else if (nobjects == j) {
+		    j = j + 100
+		    call realloc (MKO_MKT(mko), j, TY_POINTER)
+		    call realloc (MKO_X(mko), j, TY_REAL)
+		    call realloc (MKO_Y(mko), j, TY_REAL)
+		    call realloc (MKO_Z(mko), j, TY_REAL)
+		    call realloc (MKO_R(mko), j, TY_REAL)
+		    call realloc (MKO_AR(mko), j, TY_REAL)
+		    call realloc (MKO_PA(mko), j, TY_REAL)
+		    call realloc (MKO_SAVE(mko), j, TY_INT)
+		    call realloc (MKO_SORT(mko), j, TY_INT)
+		}
+
+		Memr[MKO_X(mko)+nobjects] = x
+		Memr[MKO_Y(mko)+nobjects] = y
+		Memr[MKO_Z(mko)+nobjects] =
+		    exptime / (distance * distance) * 10. ** (-0.4*(z-m0))
+		if (nscan() < 7)
+		    Memi[MKO_MKT(mko)+nobjects] = mkt_star (Memc[star])
+		else {
+		    Memi[MKO_MKT(mko)+nobjects] = mkt_object (Memc[type])
+		    Memr[MKO_R(mko)+nobjects] = r / distance
+		    Memr[MKO_AR(mko)+nobjects] = ar
+		    Memr[MKO_PA(mko)+nobjects] = DEGTORAD (pa)
+		    if (nscan() == 8)
+			Memi[MKO_SAVE(mko)+nobjects] = btoi (bsave)
+		    else
+			Memi[MKO_SAVE(mko)+nobjects] = NO
+		}
+		Memi[MKO_SORT(mko)+nobjects] = nobjects
+		nobjects = nobjects + 1
+	    }
+	    call close (i)
+
+	    # If no objects are requested then do the image I/O
+	    # line by line.  Add noise if creating a new image or
+	    # copy the input image if a new output image is desired.
+	    # Then go on to the next image.
+
+	    irbuf = 0
+	    ipbuf = 0
+	    if (nobjects == 0) {
+		call mfree (mko, TY_STRUCT)
+
+		if (new) {
+		    do i = 1, nl {
+			obuf = impl2r (out, i)
+			if (background == 0.)
+			    call aclrr (Memr[obuf], nc)
+			else
+			    call amovkr (background, Memr[obuf], nc)
+			if (poisson)
+			    call mkpnoise (Memr[obuf], Memr[obuf], nc, 0., gain,
+				pbuf, ranbuf, ipbuf, seed1)
+			if (rdnoise > 0.)
+			    call mkrnoise (Memr[obuf], nc, rdnoise,
+				rbuf, ranbuf, irbuf, seed1)
+			call alimr (Memr[obuf], nc, dmin, dmax)
+			IM_MIN(out) = min (IM_MIN(out), dmin)
+			IM_MAX(out) = max (IM_MAX(out), dmax)
+		    }
+		} else if (in != out) {
+		    do i = 1, nl {
+			obuf = impl2r (out, i)
+			call amovr (Memr[imgl2r(in,i)], Memr[obuf],
+			    IM_LEN(in,1))
+			call alimr (Memr[obuf], nc, dmin, dmax)
+			IM_MIN(out) = min (IM_MIN(out), dmin)
+			IM_MAX(out) = max (IM_MAX(out), dmax)
+		    }
+		}
+
+		# Add comment history of task parameters.
+		if (cmmts) {
+		    call strcpy ("# ", Memc[comment], LEN_COMMENT)
+		    call cnvtime (clktime (0), Memc[comment+2], LEN_COMMENT-2)
+		    call mkh_comment (out, Memc[comment])
+		    call mkh_comment (out, "begin        mkobjects")
+		    call mkh_comment (out, Memc[comment])
+		    call mkh_comment1 (out, "background", 'r')
+		    call mkh_comment1 (out, "xoffset", 'r')
+		    call mkh_comment1 (out, "yoffset", 'r')
+		    call mkh_comment1 (out, "distance", 'r')
+		    call mkh_comment1 (out, "exptime", 'r')
+		    call mkh_comment1 (out, "magzero", 'r')
+		    call mkh_comment1 (out, "gain", 'r')
+		    call mkh_comment1 (out, "rdnoise", 'r')
+		    call mkh_comment1 (out, "poisson", 'b')
+		    call mkh_comment1 (out, "seed", 'i')
+		}
+
+		IM_LIMTIME(out) = IM_MTIME(out) + 1
+		if (in != out)
+		    call imunmap (in)
+		call imunmap (out)
+		next
+	    }
+
+	    # Add the objects.
+	    #
+	    # The object list is first sorted in Y for efficiency.
+	    # Get buffer of as many lines as possible to minimize random
+	    # access and speed up adding the objects.  Ideally the whole
+	    # image should be in memory but if not we scroll a buffer
+	    # using the fact that the objects are ordered in Y.
+	    # Use error checking to determine how much memory is available.
+
+	    mko_sort = MKO_Y(mko)
+	    call qsort (Memi[MKO_SORT(mko)], nobjects, mko_compare)
+
+	    for (nlines=nl;; nlines = 0.8 * nlines)
+		ifnoerr (call malloc (buf, nlines * nc, TY_REAL))
+		    break
+	    call malloc (lines, nlines, TY_INT)
+	    call malloc (newlines, nl, TY_INT)
+	    call amovki (YES, Memi[newlines], nl)
+
+	    # Fill the line buffer.
+	    do l = 1, nlines {
+		j = mod (l, nlines) 
+		ptr2 = buf + j * nc
+		Memi[lines+j] = l
+		if (new)
+		    call aclrr (Memr[ptr2], nc)
+		else
+		    call amovr (Memr[imgl2r(in,l)], Memr[ptr2], nc) 
+		Memi[newlines+l-1] = NO
+	    }
+		
+	    # Generate the object subrasters, add noise if needed, and
+	    # add the data to the line buffer.  Check for parts of the
+	    # object off the image (the object center is guarenteed to
+	    # be on the image).  Do image I/O if needed.
+
+	    do i = 0, nobjects-1 {
+		j = Memi[MKO_SORT(mko)+i]
+		mkt = Memi[MKO_MKT(mko)+j]
+		if (mkt == NULL)
+		    next
+		x = Memr[MKO_X(mko)+j]
+		y = Memr[MKO_Y(mko)+j]
+		z = Memr[MKO_Z(mko)+j]
+		r = Memr[MKO_R(mko)+j]
+		ar = Memr[MKO_AR(mko)+j]
+		pa = Memr[MKO_PA(mko)+j]
+		save = Memi[MKO_SAVE(mko)+j]
+
+		call mkt_gobject (mkt, obj, nx, ny, x, y, z, r, ar, pa, save)
+
+		c1 = nint (x) - nx/2
+		c2 = c1 + nx - 1
+		c3 = max (1, c1)
+		c4 = min (nc, c2)
+		l1 = nint (y) - ny/2
+		l2 = l1 + ny - 1
+		l3 = max (1, l1)
+		l4 = min (nl, l2)
+		k = c4 - c3 + 1 
+		ptr1 = obj + (l3 - l1) * nx + c3 - c1
+		c3 = c3 - 1
+		do l = l3, l4 {
+		    j = mod (l, nlines)
+		    if (l != Memi[lines+j]) {
+			ptr2 = buf + j * nc
+			obuf = impl2r (out, Memi[lines+j])
+			call amovr (Memr[ptr2], Memr[obuf], nc)
+			call alimr (Memr[obuf], nc, dmin, dmax)
+			IM_MIN(out) = min (IM_MIN(out), dmin)
+			IM_MAX(out) = max (IM_MAX(out), dmax)
+			Memi[lines+j] = l
+			if (Memi[newlines+l-1] == NO)
+			    call amovr (Memr[imgl2r(out,l)], Memr[ptr2], nc)
+			else if (new)
+			    call aclrr (Memr[ptr2], nc)
+			else
+			    call amovr (Memr[imgl2r(in,l)], Memr[ptr2], nc)
+			Memi[newlines+l-1] = NO
+		    }
+		    ptr2 = buf + j * nc + c3
+		    call aaddr (Memr[ptr1], Memr[ptr2], Memr[ptr2], k)
+		    if (!new) {
+		        if (poisson)
+			    call mkpnoise (Memr[ptr1], Memr[ptr2], k,
+				background, gain, pbuf, ranbuf, ipbuf, seed1)
+		        if (rdnoise > 0.)
+			    call mkrnoise (Memr[ptr2], k,
+				rdnoise, rbuf, ranbuf, irbuf, seed1)
+		    }
+		    ptr1 = ptr1 + nx
+		}
+	    }
+
+	    # Flush out the line buffer.  A new image requires addition of
+	    # background and noise.  If the whole image is in memory then
+	    # we can add the background and noise before flushing the data.
+	    # Otherwise, we need a second pass reading the image in line
+	    # by line and adding the background and noise.  Note that if
+	    # the image was not new then noise was added only to the
+	    # objects.
+	    
+	    if (nlines == nl) {
+		do i = 1, nlines {
+		    j = mod (i, nlines) 
+		    ptr2 = buf + j * nc
+		    l = Memi[lines+j]
+		    if (new) {
+			if (background != 0.)
+			    call aaddkr (Memr[ptr2], background, Memr[ptr2], nc)
+			if (poisson)
+			    call mkpnoise (Memr[ptr2], Memr[ptr2], nc, 0., gain,
+				pbuf, ranbuf, ipbuf, seed1)
+			if (rdnoise > 0.)
+			    call mkrnoise (Memr[ptr2], nc, rdnoise,
+				rbuf, ranbuf, irbuf, seed1)
+		    }
+		    obuf = impl2r (out, l)
+		    call amovr (Memr[ptr2], Memr[obuf], nc) 
+		    call alimr (Memr[obuf], nc, dmin, dmax)
+		    IM_MIN(out) = min (IM_MIN(out), dmin)
+		    IM_MAX(out) = max (IM_MAX(out), dmax)
+		}
+	    } else {
+		do i = 1, nlines {
+		    j = mod (i, nlines) 
+		    ptr2 = buf + j * nc
+		    l = Memi[lines+j]
+		    obuf = impl2r (out, l)
+		    call amovr (Memr[ptr2], Memr[obuf], nc) 
+		    call alimr (Memr[obuf], nc, dmin, dmax)
+		    IM_MIN(out) = min (IM_MIN(out), dmin)
+		    IM_MAX(out) = max (IM_MAX(out), dmax)
+		}
+
+		if (new) {
+		    call imflush (out)
+		    do i = 1, nl {
+			obuf = impl2r (out, i)
+			ptr1 = imgl2r (out, i)
+			if (background == 0.)
+			    call amovr (Memr[ptr1], Memr[obuf], nc)
+			else
+			    call aaddkr (Memr[ptr1], background, Memr[obuf], nc)
+			if (poisson)
+			    call mkpnoise (Memr[obuf], Memr[obuf], nc, 0., gain,
+				pbuf, ranbuf, ipbuf, seed1)
+			if (rdnoise > 0.)
+			    call mkrnoise (Memr[obuf], nc, rdnoise,
+				rbuf, ranbuf, irbuf, seed1)
+			call alimr (Memr[obuf], nc, dmin, dmax)
+			IM_MIN(out) = min (IM_MIN(out), dmin)
+			IM_MAX(out) = max (IM_MAX(out), dmax)
+		    }
+		}
+	    }
+
+	    # Since each image is different and the object lists may be
+	    # different we free most of the memory within the image list
+	    # loop.
+
+	    call mfree (buf, TY_REAL)
+	    call mfree (lines, TY_INT)
+	    call mfree (newlines, TY_INT)
+	    call mfree (MKO_MKT(mko), TY_POINTER)
+	    call mfree (MKO_X(mko), TY_REAL)
+	    call mfree (MKO_Y(mko), TY_REAL)
+	    call mfree (MKO_Z(mko), TY_REAL)
+	    call mfree (MKO_R(mko), TY_REAL)
+	    call mfree (MKO_AR(mko), TY_REAL)
+	    call mfree (MKO_PA(mko), TY_REAL)
+	    call mfree (MKO_SAVE(mko), TY_INT)
+	    call mfree (MKO_SORT(mko), TY_INT)
+	    call mfree (mko, TY_STRUCT)
+
+	    # Add comment history of task parameters.
+	    if (cmmts) {
+		call strcpy ("# ", Memc[comment], LEN_COMMENT)
+		call cnvtime (clktime (0), Memc[comment+2], LEN_COMMENT-2)
+		call mkh_comment (out, Memc[comment])
+		call mkh_comment (out, "begin        mkobjects")
+		call sprintf (Memc[comment], LEN_COMMENT, "%9t%s%24t%s")
+		    call pargstr ("objects")
+		    call pargstr (Memc[fname])
+		call mkh_comment (out, Memc[comment])
+		call mkh_comment1 (out, "background", 'r')
+		call mkh_comment1 (out, "xoffset", 'r')
+		call mkh_comment1 (out, "yoffset", 'r')
+		call mkh_comment1 (out, "star", 's')
+		call mkh_comment1 (out, "radius", 'r')
+		call mkh_comment1 (out, "beta", 'r')
+		call mkh_comment1 (out, "ar", 'r')
+		call mkh_comment1 (out, "pa", 'r')
+		call mkh_comment1 (out, "distance", 'r')
+		call mkh_comment1 (out, "exptime", 'r')
+		call mkh_comment1 (out, "magzero", 'r')
+		call mkh_comment1 (out, "gain", 'r')
+		call mkh_comment1 (out, "rdnoise", 'r')
+		call mkh_comment1 (out, "poisson", 'b')
+		call mkh_comment1 (out, "seed", 'i')
+		if (fcmmts) {
+		    i = open (Memc[fname], READ_ONLY, TEXT_FILE)
+		    while (fscan (i) != EOF) {
+			call gargr (x)
+			call gargr (y)
+			call gargr (z)
+			if (nscan() < 3) {
+			    call reset_scan ()
+			    call gargstr (Memc[comment], LEN_COMMENT)
+			    call mkh_comment (out, Memc[comment])
+			}
+		    }
+		    call close (i)
+		}
+	    }
+
+	    IM_LIMTIME(out) = IM_MTIME(out) + 1
+	    if (in != out)
+		call imunmap (in)
+	    call imunmap (out)
+
+	}
+
+	call mkt_free ()
+	call imtclose (ilist)
+	call imtclose (olist)
+	call sfree (sp)
+end
+
+
+# MKO_COMPARE -- Compare two values in the mko_sort array.
+
+int procedure mko_compare (i, j)
+
+int	i, j		# Array indices to be compared.
+
+pointer	mko_sort
+common	/mko_qsort/ mko_sort
+
+begin
+	if (Memr[mko_sort+i] < Memr[mko_sort+j])
+	    return (-1)
+	else if (Memr[mko_sort+i] > Memr[mko_sort+j])
+	    return (1)
+	else
+	    return (0)
+end