include include include define LEN_UA 20000 # Maximum user header define LEN_COMMENT 70 # Maximum comment length # Cosmic ray data structure define LEN_MKO 4 define MKO_X Memi[$1] # X position define MKO_Y Memi[$1+1] # Y position define MKO_Z Memi[$1+2] # Flux define MKO_SORT Memi[$1+3] # Sort index # T_MKNOISE -- Add cosmic rays and possion and readout noise to images. # New images may be created or noise added to existing images. # The noise is not completely random for reasons of speed. procedure t_mknoise () int ilist # Input image list int olist # Output image list int objects # List of cosmic ray files 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? long seed # Random number seed int nobjects # Number of random cosmic rays real energy # Maximum random energy (electrons) bool cmmts # Add comments? bool new, fcmmts long seed1 real x, y, z, dmin, dmax int i, j, k, l, nx, ny, nlines, c1, c2, c3, c4, l1, l2, l3, l4, irbuf, ipbuf pointer sp, input, output, fname, 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() int clgeti(), access(), nowhite(), open(), fscan(), nscan() real clgetr(), urand() pointer immap(), imgl2r(), impl2r() pointer mkt_star() errchk open, immap, imgl2r, impl2r, malloc, realloc, mkt_gstar 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 (comment, max (SZ_FNAME,LEN_COMMENT), TY_CHAR) # Get parameters which apply to all images. ilist = imtopenp ("input") olist = imtopenp ("output") objects = imtopenp ("cosrays") 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) seed = clgetl ("seed") if (IS_INDEFL(seed)) seed1 = seed1 + clktime (long (0)) else seed1 = seed cmmts = clgetb ("comments") if (imtlen (ilist) == 0) call error (1, "No input image list") # Loop through input, output, and cosmic ray lists. # Missing output images take the input image name. # The cosmic ray 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) i = imtgetim (objects, Memc[fname], SZ_FNAME) # Map images. Check for new, existing, and in-place images. if (streq (Memc[input], Memc[output])) { ifnoerr (in = immap (Memc[input], READ_WRITE, LEN_UA)) { out = in new = false } else { iferr (out = immap (Memc[output], NEW_IMAGE, LEN_UA)) { call erract (EA_WARN) next } in = out call clgstr ("header", Memc[comment], SZ_FNAME) iferr (call mkh_header (out, Memc[comment], true, false)) call erract (EA_WARN) IM_NDIM(in) = 2 IM_LEN(in,1) = clgeti ("ncols") IM_LEN(in,2) = clgeti ("nlines") IM_PIXTYPE(in) = TY_REAL call clgstr ("title", IM_TITLE(out), SZ_IMTITLE) 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 call imaddr (out, "gain", gain) call imaddr (out, "rdnoise", rdnoise * gain) # Read the object list. call malloc (mko, LEN_MKO, TY_STRUCT) call mkt_init () # Set cosmic ray templates. mkt = mkt_star ("gaussian") # Read the object list. If none or a nonexistent list is given # and a number of random events is specified then generate them. # If a nonexistent object list is given then write the random # events out. fcmmts = false energy = INDEF nobjects = 0 i = nowhite (Memc[fname], Memc[fname], SZ_FNAME) if (access (Memc[fname], 0, 0) == YES) { i = open (Memc[fname], READ_ONLY, TEXT_FILE) while (fscan (i) != EOF) { call gargr (x) call gargr (y) call gargr (z) if (nscan() < 3) { fcmmts = true next } if (x < 1 || x > nc || y < 1 || y > nl) next if (nobjects == 0) { j = 100 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_SORT(mko), j, TY_INT) } else if (nobjects == j) { j = j + 100 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_SORT(mko), j, TY_INT) } Memr[MKO_X(mko)+nobjects] = x Memr[MKO_Y(mko)+nobjects] = y Memr[MKO_Z(mko)+nobjects] = z / gain Memi[MKO_SORT(mko)+nobjects] = nobjects nobjects = nobjects + 1 } call close (i) } else { nobjects = clgeti ("ncosrays") if (nobjects > 0) { energy = clgetr ("energy") / gain call malloc (MKO_X(mko), nobjects, TY_REAL) call malloc (MKO_Y(mko), nobjects, TY_REAL) call malloc (MKO_Z(mko), nobjects, TY_REAL) call malloc (MKO_SORT(mko), nobjects, TY_INT) do i = 0, nobjects-1 { Memr[MKO_X(mko)+i] = 1 + (nc-1) * urand (seed1) Memr[MKO_Y(mko)+i] = 1 + (nl-1) * urand (seed1) Memr[MKO_Z(mko)+i] = energy * urand (seed1) Memi[MKO_SORT(mko)+i] = i } if (Memc[fname] != EOS) { i = open (Memc[fname], NEW_FILE, TEXT_FILE) do j = 0, nobjects-1 { call fprintf (i, "%g %g %g\n") call pargr (Memr[MKO_X(mko)+j]) call pargr (Memr[MKO_Y(mko)+j]) call pargr (gain * Memr[MKO_Z(mko)+j]) } call close (i) } } } # If no objects are requested then do the image I/O # line by line to add requested background and noise # and then go on to the next image. irbuf = 0 ipbuf = 0 if (nobjects == 0) { call mkt_free () 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 { do i = 1, nl { obuf = impl2r (out, i) ptr1 = imgl2r (in, 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) } } # 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 mknoise") call mkh_comment1 (out, "background", '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 cosmic rays. # # 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 (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 } do i = 0, nobjects-1 { j = Memi[MKO_SORT(mko)+i] x = Memr[MKO_X(mko)+j] y = Memr[MKO_Y(mko)+j] z = Memr[MKO_Z(mko)+j] call mkt_gstar (mkt, obj, nx, ny, x, y, z) c1 = x - nx/2 + 0.5 c2 = c1 + nx - 1 c3 = max (1, c1) c4 = min (nc, c2) l1 = y - ny/2 + 0.5 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) ptr1 = ptr1 + nx } } # Flush out the line buffer. 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. if (nlines == nl) { do i = 1, nlines { j = mod (i, nlines) ptr2 = buf + j * nc l = Memi[lines+j] 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) } 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 mkt_free () call mfree (MKO_X(mko), TY_REAL) call mfree (MKO_Y(mko), TY_REAL) call mfree (MKO_Z(mko), TY_REAL) 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 mknoise") call mkh_comment1 (out, "background", '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 && Memc[fname] != EOS) { call mkh_comment1 (out, "cosrays", 's') 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) } call sprintf (Memc[comment], LEN_COMMENT, "%9tncosrays%24t%d") call pargi (nobjects) call mkh_comment (out, Memc[comment]) if (!IS_INDEF (energy)) call mkh_comment1 (out, "energy", 'r') call mkh_comment1 (out, "radius", 'r') call mkh_comment1 (out, "ar", 'r') call mkh_comment1 (out, "pa", 'r') } IM_LIMTIME(out) = IM_MTIME(out) + 1 if (in != out) call imunmap (in) call imunmap (out) } call imtclose (ilist) call imtclose (olist) call sfree (sp) end # MKRNOISE -- Make gaussian read noise. A buffer of saved noise values may be # used to greatly speed up the noise. In this case new noise values # are randomly chosen from the buffer. procedure mkrnoise (data, ndata, rdnoise, buf, nbuf, ibuf, seed) real data[ndata] # Output data int ndata # Number of data points real rdnoise # Read noise (in data units) pointer buf # Random value buffer int nbuf # Size of random value buffer (may be zero) int ibuf # Number of random numbers saved # ibuf < nbuf Save new values # ibuf = nbuf Use saved values # ibuf > nbuf Use new values long seed # Random number seed int i real val, urand(), gasdev() begin if (ibuf == nbuf) do i = 1, ndata data[i] = data[i] + Memr[buf+int(nbuf*urand (seed))] else if (ibuf > nbuf) do i = 1, ndata data[i] = data[i] + rdnoise * gasdev (seed) else { do i = 1, ndata { if (ibuf < nbuf) { val = rdnoise * gasdev (seed) Memr[buf+ibuf] = val ibuf = ibuf + 1 } else val = Memr[buf+int(nbuf*urand (seed))] data[i] = data[i] + val } } end # MKPNOISE -- Make poisson noise. For speed, values greater than 20 # use a gaussian approximation with the square root of the value as # the sigma. The normalized gaussian values may be saved and reused # by random selection to speed things up. procedure mkpnoise (in, data, ndata, b, g, buf, nbuf, ibuf, seed) real in[ndata] # Data to add noise real data[ndata] # Output data int ndata # Number of data points real b # Background (in data units) real g # Gain pointer buf # Random value buffer int nbuf # Size of random value buffer (may be zero) int ibuf # Number of random numbers saved # ibuf < nbuf Save new values # ibuf = nbuf Use saved values # ibuf > nbuf Use new values long seed # Random number seed int i real v, gv, urand(), poidev(), gasdev() begin if (ibuf == nbuf) do i = 1, ndata { v = g * (in[i] + b) if (v < 20.) data[i] = data[i] + (poidev (v, seed) - v) / g else data[i] = data[i] + sqrt (v) * Memr[buf+int(nbuf*urand(seed))] / g } else if (ibuf > nbuf) do i = 1, ndata { v = g * (in[i] + b) data[i] = data[i] + (poidev (v, seed) - v) / g } else { do i = 1, ndata { v = g * (in[i] + b) if (v < 20.) data[i] = data[i] + (poidev (v, seed) - v) / g else { if (ibuf < nbuf) { gv = gasdev (seed) Memr[buf+ibuf] = gv ibuf = ibuf + 1 } else gv = Memr[buf+int(nbuf*urand (seed))] data[i] = data[i] + sqrt (v) * gv / g } } } end