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
|
include <error.h>
include <imhdr.h>
# T_CRGROW -- Grow cosmic ray mask identifications.
procedure t_crgrow ()
int input # Input masks
int output # Output masks
real radius # Radius
int inval # Input mask value to grow
int outval # Output grown mask value
pointer sp, inmask, outmask, temp1, temp2
pointer im, ptr
int imtopenp(), imtlen(), imtgetim()
bool strne()
int clgeti()
real clgetr()
pointer immap()
errchk immap, crgrow
begin
call smark (sp)
call salloc (inmask, SZ_FNAME, TY_CHAR)
call salloc (outmask, SZ_FNAME, TY_CHAR)
call salloc (temp1, SZ_FNAME, TY_CHAR)
call salloc (temp2, SZ_FNAME, TY_CHAR)
# Task parameters.
input = imtopenp ("input")
output = imtopenp ("output")
radius = max (0., clgetr ("radius"))
inval = clgeti ("inval")
outval = clgeti ("outval")
if (imtlen (output) != imtlen (input))
call error (1, "Input and output lists do not match")
# Grow the cosmic ray masks.
while (imtgetim (input, Memc[inmask], SZ_FNAME) != EOF) {
call strcpy (Memc[inmask], Memc[outmask], SZ_FNAME)
if (imtgetim (output, Memc[outmask], SZ_FNAME) == EOF)
call error (1, "Output list ended prematurely")
if (strne (Memc[inmask], Memc[outmask])) {
call imgcluster (Memc[inmask], Memc[temp1], SZ_FNAME)
call imgcluster (Memc[outmask], Memc[temp2], SZ_FNAME)
iferr (call imcopy (Memc[temp1], Memc[temp2])) {
call erract (EA_WARN)
next
}
im = immap (Memc[inmask], READ_ONLY, TY_CHAR)
iferr (call imgstr (im, "extname", Memc[temp1], SZ_FNAME))
call strcpy ("pl", Memc[temp1], SZ_FNAME)
call imunmap (im)
}
call xt_maskname (Memc[outmask], Memc[temp1], 0, Memc[outmask],
SZ_FNAME)
if (radius < 1.)
next
iferr {
im = NULL
ptr = immap (Memc[outmask], READ_WRITE, 0); im = ptr
call crgrow (im, radius, inval, outval)
} then {
call erract (EA_WARN)
if (strne (Memc[inmask], Memc[outmask])) {
if (im != NULL) {
call imunmap (im)
iferr (call imdelete (Memc[outmask]))
call erract (EA_WARN)
}
}
}
if (im != NULL)
call imunmap (im)
}
call imtclose (output)
call imtclose (input)
call sfree (sp)
end
# CRGROW -- Grow cosmic rays.
procedure crgrow (im, grow, inval, outval)
pointer im # Mask pointer (Read/Write)
real grow # Radius (pixels)
int inval # Input mask value for pixels to grow
int outval # Output mask value for grown pixels
int i, j, k, l, nc, nl, ngrow, nbufs, val1, val2
long v1[2], v2[2]
real grow2, y2
pointer , buf, buf1, buf2, ptr
int imgnli(), impnli()
errchk calloc, imgnli, impnli
begin
if (grow < 1. || inval == 0)
return
grow2 = grow * grow
ngrow = int (grow)
buf = NULL
iferr {
if (IM_NDIM(im) > 2)
call error (1,
"Only one or two dimensional masks are allowed")
nc = IM_LEN(im, 1)
if (IM_NDIM(im) > 1)
nl = IM_LEN(im,2)
else
nl = 1
# Initialize buffering.
nbufs = min (1 + 2 * ngrow, nl)
call calloc (buf, nc*nbufs, TY_INT)
call amovkl (long(1), v1, IM_NDIM(im))
call amovkl (long(1), v2, IM_NDIM(im))
while (imgnli (im, buf1, v1) != EOF) {
j = v1[2] - 1
buf2 = buf + mod (j, nbufs) * nc
do i = 1, nc {
val1 = Memi[buf1]
val2 = Memi[buf2]
if ((IS_INDEFI(inval) && val1 != 0) || val1 == inval) {
do k = max(1,j-ngrow), min (nl,j+ngrow) {
ptr = buf + mod (k, nbufs) * nc - 1
y2 = (k - j) ** 2
do l = max(1,i-ngrow), min (nc,i+ngrow) {
if ((l-i)**2 + y2 > grow2)
next
Memi[ptr+l] = -val1
}
}
} else {
if (val2 >= 0)
Memi[buf2] = val1
}
buf1 = buf1 + 1
buf2 = buf2 + 1
}
if (j > ngrow) {
while (impnli (im, buf2, v2) != EOF) {
k = v2[2] - 1
buf1 = buf + mod (k, nbufs) * nc
do i = 1, nc {
val1 = Memi[buf1]
if (val1 < 0) {
if (IS_INDEFI(outval))
Memi[buf2] = -val1
else
Memi[buf2] = outval
} else
Memi[buf2] = val1
Memi[buf1] = 0.
buf1 = buf1 + 1
buf2 = buf2 + 1
}
if (j != nl)
break
}
}
}
} then
call erract (EA_ERROR)
if (buf != NULL)
call mfree (buf, TY_INT)
end
|
