diff options
Diffstat (limited to 'noao/imred/ccdred/ccdtest/t_mkimage.x')
| -rw-r--r-- | noao/imred/ccdred/ccdtest/t_mkimage.x | 204 |
1 files changed, 204 insertions, 0 deletions
diff --git a/noao/imred/ccdred/ccdtest/t_mkimage.x b/noao/imred/ccdred/ccdtest/t_mkimage.x new file mode 100644 index 00000000..ff0d5f26 --- /dev/null +++ b/noao/imred/ccdred/ccdtest/t_mkimage.x @@ -0,0 +1,204 @@ +include <imhdr.h> + +define OPTIONS "|make|replace|add|multiply|" +define MAKE 1 # Create a new image +define REPLACE 2 # Replace pixels +define ADD 3 # Add to pixels +define MULTIPLY 4 # Multiply pixels + +# T_MKIMAGE -- Make or edit an image with simple values. +# An image may be created of a specified size, dimensionality, and pixel +# datatype. The image may also be edited to replace, add, or multiply +# by specified values. The values may be a combination of a sloped plane +# (repeated for dimensions greater than 2) and Gaussian noise. +# The editing may be confined to sections of the image by use of image +# sections in the input image. This task is a simple tool for +# specialized uses in test applications. +# +# The sloped plane is defined such that: +# +# pix[i,j] = value + slope * ((ncols + nlines) / 2 - 1) + slope * (i + j) +# +# The interpretation of value is that it is the mean of the plane. +# +# The Gaussian noise is only approximately random for purposes of speed! + +procedure t_mkimage () + +char image[SZ_FNAME] # Image to edit +char option[7] # Edit option +real value # Edit value +real slope # Slope +real sigma # Gaussian noise sigma +long seed # Random number seed + +int i, op, ncols, nlines +long vin[IM_MAXDIM], vout[IM_MAXDIM] +pointer sp, rannums, im, buf, bufin, bufout + +int clgwrd(), clgeti(), clscan(), nscan() imgnlr(), impnlr() +char clgetc() +real clgetr() +long clgetl() +pointer immap() + +data seed/1/ + +begin + call smark (sp) + call clgstr ("image", image, SZ_FNAME) + op = clgwrd ("option", option, 7, OPTIONS) + value = clgetr ("value") + slope = clgetr ("slope") + sigma = clgetr ("sigma") + if (clgetl ("seed") > 0) + seed = clgetl ("seed") + + call amovkl (long (1), vin, IM_MAXDIM) + call amovkl (long (1), vout, IM_MAXDIM) + switch (op) { + case MAKE: + im = immap (image, NEW_IMAGE, 0) + IM_NDIM(im) = clgeti ("ndim") + i = clscan ("dims") + do i = 1, IM_NDIM(im) + call gargi (IM_LEN(im, i)) + if (nscan() != IM_NDIM(im)) + call error (0, "Bad dimension string") + switch (clgetc ("pixtype")) { + case 's': + IM_PIXTYPE(im) = TY_SHORT + case 'i': + IM_PIXTYPE(im) = TY_INT + case 'l': + IM_PIXTYPE(im) = TY_LONG + case 'r': + IM_PIXTYPE(im) = TY_REAL + case 'd': + IM_PIXTYPE(im) = TY_DOUBLE + default: + call error (0, "Bad pixel type") + } + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + call salloc (rannums, 2 * ncols, TY_REAL) + call mksigma (sigma, seed, Memr[rannums], 2*ncols) + + while (impnlr (im, bufout, vout) != EOF) + call mkline (value, slope, sigma, seed, Memr[rannums], + Memr[bufout], vout[2] - 1, ncols, nlines) + case REPLACE: + im = immap (image, READ_WRITE, 0) + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + call salloc (rannums, 2 * ncols, TY_REAL) + call mksigma (sigma, seed, Memr[rannums], 2*ncols) + + while (impnlr (im, bufout, vout) != EOF) + call mkline (value, slope, sigma, seed, Memr[rannums], + Memr[bufout], vout[2] - 1, ncols, nlines) + case ADD: + im = immap (image, READ_WRITE, 0) + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + call salloc (buf, ncols, TY_REAL) + call salloc (rannums, 2 * ncols, TY_REAL) + call mksigma (sigma, seed, Memr[rannums], 2*ncols) + + while (imgnlr (im, bufin, vin) != EOF) { + i = impnlr (im, bufout, vout) + call mkline (value, slope, sigma, seed, Memr[rannums], + Memr[buf], vout[2] - 1, ncols, nlines) + call aaddr (Memr[bufin], Memr[buf], Memr[bufout], ncols) + } + case MULTIPLY: + im = immap (image, READ_WRITE, 0) + + ncols = IM_LEN(im,1) + nlines = IM_LEN(im,2) + call salloc (buf, ncols, TY_REAL) + call salloc (rannums, 2 * ncols, TY_REAL) + call mksigma (sigma, seed, Memr[rannums], 2*ncols) + + while (imgnlr (im, bufin, vin) != EOF) { + i = impnlr (im, bufout, vout) + call mkline (value, slope, sigma, seed, Memr[rannums], + Memr[buf], vout[2] - 1, ncols, nlines) + call amulr (Memr[bufin], Memr[buf], Memr[bufout], ncols) + } + } + + call imunmap (im) + call sfree (sp) +end + + +# MKLINE -- Make a line of data. A slope of zero is a special case. +# The Gaussian random numbers are taken from the sequence of stored +# values with starting point chosen randomly in the interval 1 to ncols. +# This is not very random but is much more efficient. + +procedure mkline (value, slope, sigma, seed, rannums, data, line, ncols, nlines) + +real value # Mean value +real slope # Slope in mean +real sigma # Sigma about mean +long seed # Random number seed +real rannums[ARB] # Random numbers +real data[ncols] # Data for line +int line # Line number +int ncols # Number of columns +int nlines # Number of lines + +int i +real a, urand() + +begin + if (slope == 0.) + call amovkr (value, data, ncols) + else { + a = value + slope * (line - (ncols + nlines) / 2. - 1) + do i = 1, ncols + data[i] = a + slope * i + } + if (sigma > 0.) { + i = (ncols - 1) * urand (seed) + 1 + call aaddr (rannums[i], data, data, ncols) + } +end + + +# MKSIGMA -- A sequence of random numbers of the specified sigma and +# starting seed is generated. The random number generator is modeled after +# that in Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. + +procedure mksigma (sigma, seed, rannums, nnums) + +real sigma # Sigma for random numbers +long seed # Seed for random numbers +real rannums[nnums] # Random numbers +int nnums # Number of random numbers + +int i +real v1, v2, r, fac, urand() + +begin + if (sigma > 0.) { + for (i=1; i<=nnums; i=i+1) { + repeat { + v1 = 2 * urand (seed) - 1. + v2 = 2 * urand (seed) - 1. + r = v1 ** 2 + v2 ** 2 + } until ((r > 0) && (r < 1)) + fac = sqrt (-2. * log (r) / r) * sigma + rannums[i] = v1 * fac + if (i == nnums) + break + i = i + 1 + rannums[i] = v2 * fac + } + } +end |