aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/ccdred/ccdtest/t_mkimage.x
blob: ff0d5f2680d7e4cd5abeec3aa40868d4315886ef (plain) (blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
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