aboutsummaryrefslogtreecommitdiff
path: root/noao/artdata/t_mkobjects.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/artdata/t_mkobjects.x')
-rw-r--r--noao/artdata/t_mkobjects.x552
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
generated by cgit (git 2.34.1) at 2025-09-20 07:11:50 -0400