aboutsummaryrefslogtreecommitdiff
path: root/noao/imred/quadred/src/ccdproc/setoverscan.x
blob: 2fef378adefeb64832f4b69fbad42b7f3ed82491 (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
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
include	<imhdr.h>
include	<imset.h>
include <pkg/gtools.h>
include	<pkg/xtanswer.h>
include	"ccdred.h"


# SET_OVERSCAN -- Set the overscan vector.
#
#   1.  Return immediately if the overscan correction is not requested or
#	if the image has been previously corrected.
#   2.	Determine the overscan columns or lines.  This may be specifed
#	directly or indirectly through the image header or symbol table.
#   3.	Average the overscan columns or lines.
#   4.	Fit a function with the ICFIT routines to smooth the overscan vector.
#   5.  Set the processing flag.
#   6.  Log the operation (to user, logfile, and output image header).

procedure set_overscan (ccd)

pointer	ccd			# CCD structure pointer

int	i, j, nsec, navg, npts, first, last
int	nc, nl, c1, c2, l1, l2
pointer	sp, str, errstr, buf, overscan, x, y, z

real	asumr()
bool	clgetb(), ccdflag()
pointer	imgl2r()
errchk	imgl2r, fit_overscan

begin
	# Check if the user wants this operation or if it has been done.
	if (!clgetb ("overscan") || ccdflag (IN_IM(ccd), "overscan"))
	    return

	call smark (sp)
	call salloc (str, SZ_LINE, TY_CHAR)
	call salloc (errstr, SZ_LINE, TY_CHAR)
	call imstats (IN_IM(ccd), IM_IMAGENAME, Memc[str], SZ_LINE)

	# Check bias section.
	nc = IM_LEN(IN_IM(ccd),1)
	nl = IM_LEN(IN_IM(ccd),2)

	c1 = BIAS_C1(ccd)
	c2 = BIAS_C2(ccd)
	l1 = BIAS_L1(ccd)
	l2 = BIAS_L2(ccd)
	navg = c2 - c1 + 1
	npts = CCD_L2(ccd) - CCD_L1(ccd) + 1

	nsec = max (1, IN_NSEC(ccd))
	do i = 1, nsec {
	    if (BIAS_SEC(ccd) != NULL) {
		c1 = BIAS_SC1(ccd,i)
		c2 = BIAS_SC2(ccd,i)
		l1 = BIAS_SL1(ccd,i)
		l2 = BIAS_SL2(ccd,i)
	    }
	    if ((c1 < 1) || (c2 > nc) || (l1 < 1) || (l2 > nl)) {
		call sprintf (Memc[errstr], SZ_LINE,
		"Error in bias section: image=%s[%d,%d], biassec=[%d:%d,%d:%d]")
		    call pargstr (Memc[str])
		    call pargi (nc)
		    call pargi (nl)
		    call pargi (c1)
		    call pargi (c2)
		    call pargi (l1)
		    call pargi (l2)
		call error (0, Memc[errstr])
	    }
	    if ((c1 == 1) && (c2 == nc) && (l1 == 1) && (l2 == nl)) {
		call error (0,
		    "Bias section not specified or given as full image")
	    }

	    # If no processing is desired then print overscan strip and return.
	    if (clgetb ("noproc")) {
		call eprintf (
		    "  [TO BE DONE] Overscan section is [%d:%d,%d:%d].\n")
		    call pargi (c1)
		    call pargi (c2)
		    call pargi (l1)
		    call pargi (l2)
		    call sfree (sp)
		    return
	    }
	}

	# Determine the overscan section parameters. The readout axis
	# determines the type of overscan.  The step sizes are ignored.
	# The limits in the long dimension are replaced by the trim limits.

	if (READAXIS(ccd) == 1) {
	    call salloc (buf, nsec*nl, TY_REAL)
	    z = buf
	    do i = 1, nl {
		y = imgl2r (IN_IM(ccd), i)
		do j = 1, nsec {
		    if (BIAS_SEC(ccd) != NULL) {
			l1 = BIAS_SL1(ccd,j)
			l2 = BIAS_SL2(ccd,j)
			if (i < l1 || i > l2)
			    next
			c1 = BIAS_SC1(ccd,j)
			c2 = BIAS_SC2(ccd,j)
			navg = c2 - c1 + 1
			z = buf + (j - 1) * nl
		    }
		    Memr[z+i-1] = asumr (Memr[y+c1-1], navg)
		}
	    }

	    # Trim the overscan vector and set the pixel coordinate.
	    call salloc (x, nl, TY_REAL)
	    call malloc (overscan, nsec*nl, TY_REAL)
	    y = overscan
	    do i = 1, nsec {
		if (BIAS_SEC(ccd) != NULL) {
		    c1 = BIAS_SC1(ccd,i)
		    c2 = BIAS_SC2(ccd,i)
		    l1 = BIAS_SL1(ccd,i)
		    l2 = BIAS_SL2(ccd,i)
		    navg = c2 - c1 + 1
		    npts = l2 - l1 + 1
		    y = overscan + (i - 1) * nl
		    z = buf + (i - 1) * nl
		}
		if (navg > 1)
		    call adivkr (Memr[z+l1-1], real (navg), Memr[z+l1-1],
			npts)
		call trim_overscan (Memr[z], npts, l1, Memr[x], Memr[y])
		call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
		    Memr[y], npts)
	    }

	} else {
	    first = l1
	    last = l2
	    navg = last - first + 1
	    npts = nc
	    call salloc (buf, npts, TY_REAL)
	    call aclrr (Memr[buf], npts)
	    do i = first, last
	        call aaddr (Memr[imgl2r(IN_IM(ccd),i)], Memr[buf], Memr[buf],
		    npts)
	    if (navg > 1)
	        call adivkr (Memr[buf], real (navg), Memr[buf], npts)

	    # Trim the overscan vector and set the pixel coordinate.
	    npts = CCD_C2(ccd) - CCD_C1(ccd) + 1
	    call malloc (overscan, npts, TY_REAL)
	    call salloc (x, npts, TY_REAL)
	    call trim_overscan (Memr[buf], npts, IN_C1(ccd), Memr[x],
		Memr[overscan])

	    call fit_overscan (Memc[str], c1, c2, l1, l2, Memr[x],
		Memr[overscan], npts)
	}
	
	# Set the CCD structure overscan parameters.
	CORS(ccd, OVERSCAN) = O
	COR(ccd) = YES
	OVERSCAN_VEC(ccd) = overscan

	# Log the operation.
	call strcpy ("overscan", Memc[errstr], SZ_LINE)
	y = overscan
	do i = 1, nsec {
	    if (BIAS_SEC(ccd) != NULL) {
		c1 = BIAS_SC1(ccd,i)
		c2 = BIAS_SC2(ccd,i)
		l1 = BIAS_SL1(ccd,i)
		l2 = BIAS_SL2(ccd,i)
		y = overscan + (i - 1) * nl
		npts = c2 - c1 + 1
		if (i > 1) {
		    call sprintf (Memc[errstr], SZ_LINE, "ovrscn%d")
			call pargi (i)
		}
	    }
	    call sprintf (Memc[str], SZ_LINE,
		"Overscan section is [%d:%d,%d:%d] with mean=%g")
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
		call pargr (asumr (Memr[y], npts) / npts)
	    call timelog (Memc[str], SZ_LINE)
	    call ccdlog (IN_IM(ccd), Memc[str])
	    call hdmpstr (OUT_IM(ccd), Memc[errstr], Memc[str])
	}

	call sfree (sp)
end


# FIT_OVERSCAN -- Fit a function to smooth the overscan vector.
#   The fitting uses the ICFIT procedures which may be interactive.
#   Changes to these parameters are "learned".  The user is queried with a four
#   valued logical query (XT_ANSWER routine) which may be turned off when
#   multiple images are processed.

procedure fit_overscan (image, c1, c2, l1, l2, x, overscan, npts)

char	image[ARB]		# Image name for query and title
int	c1, c2, l1, l2		# Overscan strip
real	x[npts]			# Pixel coordinates of overscan
real	overscan[npts]		# Input overscan and output fitted overscan
int	npts			# Number of data points

int	interactive, fd
pointer	sp, str, w, ic, cv, gp, gt

int	clgeti(), ic_geti(), open()
real	clgetr(), ic_getr()
pointer	gopen(), gt_init()
errchk	gopen, open

begin
	call smark (sp)
	call salloc (str, SZ_LINE, TY_CHAR)
	call salloc (w, npts, TY_REAL)
	call amovkr (1., Memr[w], npts)

	# Open the ICFIT procedures, get the fitting parameters, and
	# set the fitting limits.

	call ic_open (ic)
	call clgstr ("function", Memc[str], SZ_LINE)
	call ic_pstr (ic, "function", Memc[str])
	call ic_puti (ic, "order", clgeti ("order"))
	call clgstr ("sample", Memc[str], SZ_LINE)
	call ic_pstr (ic, "sample", Memc[str])
	call ic_puti (ic, "naverage", clgeti ("naverage"))
	call ic_puti (ic, "niterate", clgeti ("niterate"))
	call ic_putr (ic, "low", clgetr ("low_reject"))
	call ic_putr (ic, "high", clgetr ("high_reject"))
	call ic_putr (ic, "grow", clgetr ("grow"))
	call ic_putr (ic, "xmin", min (x[1], x[npts]))
	call ic_putr (ic, "xmax", max (x[1], x[npts]))
	call ic_pstr (ic, "xlabel", "Pixel")
	call ic_pstr (ic, "ylabel", "Overscan")

	# If the fitting is done interactively set the GTOOLS and GIO
	# pointers.  Also "learn" the fitting parameters since they may
	# be changed when fitting interactively.

	call sprintf (Memc[str], SZ_LINE,
	    "Fit overscan vector for %s[%d:%d,%d:%d] interactively")
	    call pargstr (image)
	    call pargi (c1)
	    call pargi (c2)
	    call pargi (l1)
	    call pargi (l2)
	call set_interactive (Memc[str], interactive)
	if ((interactive == YES) || (interactive == ALWAYSYES)) {
	    gt = gt_init ()
	    call sprintf (Memc[str], SZ_LINE,
	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
	        call pargstr (image)
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
	    call gt_sets (gt, GTTITLE, Memc[str])
	    call gt_sets (gt, GTTYPE, "line")
	    call gt_setr (gt, GTXMIN, x[1])
	    call gt_setr (gt, GTXMAX, x[npts])
	    call clgstr ("graphics", Memc[str], SZ_FNAME)
	    gp = gopen (Memc[str], NEW_FILE, STDGRAPH)

	    call icg_fit (ic, gp, "cursor", gt, cv, x, overscan, Memr[w], npts)

	    call ic_gstr (ic, "function", Memc[str], SZ_LINE)
	    call clpstr ("function", Memc[str])
	    call clputi ("order", ic_geti (ic, "order"))
	    call ic_gstr (ic, "sample", Memc[str], SZ_LINE)
	    call clpstr ("sample", Memc[str])
	    call clputi ("naverage", ic_geti (ic, "naverage"))
	    call clputi ("niterate", ic_geti (ic, "niterate"))
	    call clputr ("low_reject", ic_getr (ic, "low"))
	    call clputr ("high_reject", ic_getr (ic, "high"))
	    call clputr ("grow", ic_getr (ic, "grow"))

	    call gclose (gp)
	    call gt_free (gt)
	} else
	    call ic_fit (ic, cv, x, overscan, Memr[w], npts, YES, YES, YES, YES)

	# Make a log of the fit in the plot file if given.
	call clgstr ("plotfile", Memc[str], SZ_LINE)
	call xt_stripwhite (Memc[str])
	if (Memc[str] != EOS) {
	    fd = open (Memc[str], APPEND, BINARY_FILE)
	    gp = gopen ("stdvdm", NEW_FILE, fd)
	    gt = gt_init ()
	    call sprintf (Memc[str], SZ_LINE,
	        "Overscan vector for %s from section [%d:%d,%d:%d]\n")
	        call pargstr (image)
		call pargi (c1)
		call pargi (c2)
		call pargi (l1)
		call pargi (l2)
	    call gt_sets (gt, GTTITLE, Memc[str])
	    call gt_sets (gt, GTTYPE, "line")
	    call gt_setr (gt, GTXMIN, 1.)
	    call gt_setr (gt, GTXMAX, real (npts))
	    call icg_graphr (ic, gp, gt, cv, x, overscan, Memr[w], npts)
	    call gclose (gp)
	    call close (fd)
	    call gt_free (gt)
	}

	# Replace the raw overscan vector with the smooth fit.
	call cvvector (cv, x, overscan, npts)

	# Finish up.
	call ic_closer (ic)
	call cvfree (cv)
	call sfree (sp)
end


# TRIM_OVERSCAN -- Trim the overscan vector.

procedure trim_overscan (data, npts, start, x, overscan)

real	data[ARB]		# Full overscan vector
int	npts			# Length of trimmed vector
int	start			# Trim start
real	x[npts]			# Trimmed pixel coordinates (returned)
real	overscan[npts]		# Trimmed overscan vector (returned)

int	i, j

begin
	do i = 1, npts {
	    j = start + i - 1
	    x[i] = j
	    overscan[i] = data[j]
	}
end