LSST Applications g070148d5b3+33e5256705,g0d53e28543+25c8b88941,g0da5cf3356+2dd1178308,g1081da9e2a+62d12e78cb,g17e5ecfddb+7e422d6136,g1c76d35bf8+ede3a706f7,g295839609d+225697d880,g2e2c1a68ba+cc1f6f037e,g2ffcdf413f+853cd4dcde,g38293774b4+62d12e78cb,g3b44f30a73+d953f1ac34,g48ccf36440+885b902d19,g4b2f1765b6+7dedbde6d2,g5320a0a9f6+0c5d6105b6,g56b687f8c9+ede3a706f7,g5c4744a4d9+ef6ac23297,g5ffd174ac0+0c5d6105b6,g6075d09f38+66af417445,g667d525e37+2ced63db88,g670421136f+2ced63db88,g71f27ac40c+2ced63db88,g774830318a+463cbe8d1f,g7876bc68e5+1d137996f1,g7985c39107+62d12e78cb,g7fdac2220c+0fd8241c05,g96f01af41f+368e6903a7,g9ca82378b8+2ced63db88,g9d27549199+ef6ac23297,gabe93b2c52+e3573e3735,gb065e2a02a+3dfbe639da,gbc3249ced9+0c5d6105b6,gbec6a3398f+0c5d6105b6,gc9534b9d65+35b9f25267,gd01420fc67+0c5d6105b6,geee7ff78d7+a14128c129,gf63283c776+ede3a706f7,gfed783d017+0c5d6105b6,w.2022.47
LSST Data Management Base Package
Loading...
Searching...
No Matches
Interp.cc
Go to the documentation of this file.
1// -*- LSST-C++ -*-
2
3/*
4 * LSST Data Management System
5 * Copyright 2008-2015 LSST Corporation.
6 *
7 * This product includes software developed by the
8 * LSST Project (http://www.lsst.org/).
9 *
10 * This program is free software: you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation, either version 3 of the License, or
13 * (at your option) any later version.
14 *
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
19 *
20 * You should have received a copy of the LSST License Statement and
21 * the GNU General Public License along with this program. If not,
22 * see <http://www.lsstcorp.org/LegalNotices/>.
23 */
24
32#include <stdexcept>
33#include <algorithm>
34#include <cassert>
35#include <string>
36#include <typeinfo>
37#include <limits>
38#include "boost/format.hpp"
39
40#include "lsst/geom.h"
41#include "lsst/pex/exceptions.h"
42#include "lsst/afw/image/MaskedImage.h"
43#include "lsst/meas/algorithms/Interp.h"
44
45namespace lsst {
46namespace meas {
47namespace algorithms {
48
49typedef std::vector<Defect::Ptr>::const_iterator DefectCIter;
50
51/************************************************************************************************************/
52/*
53 * Classify an vector of Defect::Ptr for the given row, returning a vector of 1-D
54 * Defects (i.e. y0 == y1). In general we can merge in saturated pixels at
55 * this step, although we don't currently do so.
56 *
57 * See comment above do_defects for a description of how to interpret DefectType
58 */
59static std::vector<Defect::Ptr> classify_defects(
60 std::vector<Defect::Ptr> const &badList, // list of bad things
61 int const y, // the row to process
62 int const ncol, // number of columns in image
63 int = 0 // number of rows in image
64) {
65 std::vector<Defect::Ptr> badList1D;
66
67 for (DefectCIter begin = badList.begin(), end = badList.end(), bri = begin; bri != end; ++bri) {
68 Defect::Ptr defect = *bri;
69
70 if (y < defect->getY0() || y > defect->getY1() || ncol < defect->getX0()) {
71 continue;
72 }
73
74 int const x0 = defect->getX0();
75 int x1 = defect->getX1();
76 //
77 // Look for other defects that touch this one, and push them onto badList1d
78 //
79 for (++bri; bri != end; ++bri) {
80 defect = *bri;
81
82 if (y < defect->getY0() || y > defect->getY1()) { // this defect doesn't concern this row
83 continue;
84 }
85 if (x1 < defect->getX0() - 1) { // no further defects can touch this one
86 --bri;
87 break;
88 }
89 if (defect->getX1() > x1) {
90 x1 = defect->getX1();
91 }
92 }
93
94 int const nbad = x1 - x0 + 1;
95 assert(nbad >= 1);
96
97 defect = Defect::Ptr(new Defect(geom::BoxI(geom::Point2I(x0, y), geom::Extent2I(nbad, 1))));
98 badList1D.push_back(defect);
99
100 if (bri == end) {
101 break;
102 }
103 }
104 //
105 // Now process our new list
106 //
107 for (DefectCIter begin = badList1D.begin(), end = badList1D.end(), bri = begin; bri != end; ++bri) {
108 Defect::Ptr defect = *bri;
109
110 int const nbad = defect->getX1() - defect->getX0() + 1;
111 assert(nbad >= 1);
112
113 if (defect->getX0() == 0) {
114 if (nbad >= Defect::WIDE_DEFECT) {
115 defect->classify(Defect::WIDE_LEFT, 03);
116 } else {
117 defect->classify(Defect::LEFT, 03 << nbad);
118 }
119 } else if (defect->getX0() == 1) { /* only second column is usable */
120 if (nbad >= Defect::WIDE_DEFECT) {
121 defect->classify(Defect::WIDE_NEAR_LEFT, (01 << 2) | 03);
122 } else {
123 defect->classify(Defect::NEAR_LEFT, (01 << (nbad + 2)) | 03);
124 }
125 } else if (defect->getX1() == ncol - 2) { /* use only penultimate column */
126 if (nbad >= Defect::WIDE_DEFECT) {
127 defect->classify(Defect::WIDE_NEAR_RIGHT, (03 << 2) | 02);
128 } else {
129 defect->classify(Defect::NEAR_RIGHT, (03 << (nbad + 2)) | 02);
130 }
131 } else if (defect->getX1() == ncol - 1) {
132 if (nbad >= Defect::WIDE_DEFECT) {
133 defect->classify(Defect::WIDE_RIGHT, 03);
134 } else {
135 defect->classify(Defect::RIGHT, 03 << nbad);
136 }
137 } else if (nbad >= Defect::WIDE_DEFECT) {
138 defect->classify(Defect::WIDE, (03 << 2) | 03);
139 } else {
140 defect->classify(Defect::MIDDLE, (03 << (nbad + 2)) | 03);
141 }
142 /*
143 * look for bad columns in regions that we'll get `good' values from.
144 *
145 * We know that no two Defects are adjacent.
146 */
147 int nshift = 0; // number of bits to shift to get to left edge of defect pattern
148 switch (defect->getPos()) {
149 case Defect::WIDE: // no bits
150 case Defect::WIDE_NEAR_LEFT: // are used to encode
151 case Defect::WIDE_NEAR_RIGHT: // the bad section of data
152 nshift = 0;
153 break;
154 default:
155 nshift = nbad;
156 break;
157 }
158
159 if (bri != begin) {
160 Defect::Ptr const defect_m = *(bri - 1);
161 assert(defect_m->getX1() < defect->getX0());
162
163 if (defect_m->getX1() == defect->getX0() - 2) {
164 defect->classify(defect->getPos(), (defect->getType() & ~(02 << (nshift + 2))));
165 }
166 }
167
168 if (bri + 1 != end) {
169 Defect::Ptr const defect_p = *(bri + 1);
170
171 if (defect->getX1() == defect_p->getX0() - 2) {
172 Defect::DefectPosition defectPos = defect->getPos();
173 if (defectPos == Defect::LEFT || defectPos == Defect::NEAR_LEFT) {
174 defect->classify(defect->getPos(), (defect->getType() & ~(02 << nshift)));
175 } else {
176 defect->classify(defectPos, (defect->getType() & ~01));
177 }
178 }
179 }
180 }
181
182 return badList1D;
183}
184
185/*****************************************************************************/
186/*
187 * Interpolate over the defects in a given line of data. In the comments,
188 * a bad pixel is written as ., a good one as #, and unknown but non-interpolated pixels as ?.
189 *
190 * This may be mapped to an int by replacing # with 1 and . or ? with 0. So "##..##" would mean, "I have two
191 * adjacent bad pixels with 2 good neighbours to both the left and right", and have a defectType of 110011 or
192 * 063. This defect is in the middle of the chip, so has DefectPosition MIDDLE.
193 *
194 * The other options are LEFT, NEAR_LEFT, WIDE_NEAR_LEFT, and the corresponding RIGHT positions, and WIDE.
195 * LEFT are defects that touch the left side, and NEAR_LEFT ones that come within a pixel. WIDE are encoded
196 * omitting the .. so "##..................##" would be 1111 == 017 and WIDE.
197
198 * The LEFT ones are actually a bit tricky as they'd have leading 0s, so they are inverted ("....##" is
199 * written as 110000 not 000011).
200 */
201template <typename ImageT>
202static void do_defects(std::vector<Defect::Ptr> const &badList, // list of bad things
203 int const y, // Row that we should fix
204 ImageT &data, // data to fix
205 typename ImageT::Pixel min, // minimum acceptable value
206 double fallbackValue, // Value to fallback to if all else fails
207 bool useFallbackValueAtEdge, // use fallbackValue at edge of chip?
208 int nUseInterp // no. of pixels to interpolate towards edge
209) {
210 typedef typename ImageT::Pixel ImagePixel;
211 ImagePixel out1_2, out1_1, out2_1, out2_2; // == out[badX1-2], ..., out[bad_x2+2]
212 ImagePixel val; // unpack a pixel value
213 //
214 // Get pointer to this row of data
215 //
216 int const ncol = data.getWidth();
217 typename ImageT::x_iterator out = data.row_begin(y);
218
219 for (DefectCIter ptr = badList.begin(), end = badList.end(); ptr != end; ++ptr) {
220 Defect::Ptr const defect = *ptr;
221
222 if (y < defect->getY0() || y > defect->getY1()) {
223 continue;
224 }
225
226 int badX0 = defect->getX0();
227 int badX1 = defect->getX1();
228
229 Defect::DefectPosition defectPos = defect->getPos();
230 unsigned int defectType = defect->getType();
231
232 int nbad = badX1 - badX0 + 1;
233
234 if (nbad > nUseInterp && useFallbackValueAtEdge) {
235 switch (defectPos) {
236 case Defect::LEFT:
237 case Defect::WIDE_LEFT:
238 assert(badX0 == 0);
239
240 if (badX1 == ncol - 1) { // also RIGHT --- spans the entire image
241 for (int i = 0; i != ncol; ++i) {
242 out[i] = fallbackValue;
243 }
244 continue;
245 }
246
247 for (; badX0 <= badX1 - nUseInterp; ++badX0) {
248 out[badX0] = fallbackValue;
249 }
250
251 if (defectPos == Defect::LEFT) {
252 defectType >>= nbad; // we just want the last 2 bits
253 switch (defectType) {
254 case 01:
255 defectType = 02;
256 break;
257 case 03:
258 defectType = 03;
259 break;
260 default:
261 throw std::runtime_error(str(
262 boost::format("Impossible value of defectType: 0%o") % defectType));
263 }
264 }
265 nbad = badX1 - badX0 + 1;
266 defectType = (03 << (nbad + 2)) | defectType;
267 defectPos = (badX0 > 1) ? ((badX1 < ncol - 2) ? Defect::MIDDLE : Defect::NEAR_RIGHT)
268 : Defect::NEAR_LEFT;
269 break;
270 case Defect::RIGHT:
271 case Defect::WIDE_RIGHT:
272 assert(badX1 == ncol - 1);
273 for (; badX1 >= badX0 + nUseInterp; --badX1) {
274 out[badX1] = fallbackValue;
275 }
276 nbad = badX1 - badX0 + 1;
277 defectType = (03 << (nbad + 2)) | 03;
278 defectPos = (badX1 < ncol - 2) ? Defect::MIDDLE : Defect::NEAR_RIGHT;
279 break;
280 default:
281 break;
282 }
283 }
284
285 switch (defectPos) {
286 case Defect::LEFT:
287 assert(badX0 >= 0 && badX1 + 2 < ncol);
288
289 out2_1 = out[badX1 + 1];
290 out2_2 = out[badX1 + 2];
291
292 switch (defectType) {
293 case 02: /* .#?, <noise^2> = 0 */
294 val = 1.0000 * out2_1;
295 out[badX1] = (val < min) ? out2_1 : val;
296
297 break;
298 case 06: /* .##, <noise^2> = 0 */
299 val = 1.4288 * out2_1 - 0.4288 * out2_2;
300 out[badX1] = (val < min) ? out2_1 : val;
301
302 break;
303 case 014: /* ..##, <noise^2> = 0 */
304 val = 1.0933 * out2_1 - 0.0933 * out2_2;
305 out[badX0] = (val < min) ? out2_1 : val;
306
307 val = 1.4288 * out2_1 - 0.4288 * out2_2;
308 out[badX1] = (val < min) ? out2_1 : val;
309
310 break;
311 case 04: /* ..#?, <noise^2> = 0 */
312 val = 1.000 * out2_1;
313 out[badX0] = (val < min) ? out2_1 : val;
314 out[badX1] = (val < min) ? out2_1 : val;
315
316 break;
317 case 030: /* ...##, <noise^2> = 0 */
318 val = 0.6968 * out2_1 + 0.3032 * out2_2;
319 out[badX0] = (val < min) ? out2_1 : val;
320
321 val = 1.0933 * out2_1 - 0.0933 * out2_2;
322 out[badX1 - 1] = (val < min) ? out2_1 : val;
323
324 val = 1.4288 * out2_1 - 0.4288 * out2_2;
325 out[badX1] = (val < min) ? out2_1 : val;
326
327 break;
328 case 010: /* ...#?, <noise^2> = 0 */
329 val = 1.000 * out2_1;
330
331 out[badX0] = (val < min) ? out2_1 : val;
332 out[badX1 - 1] = (val < min) ? out2_1 : val;
333 out[badX1] = (val < min) ? out2_1 : val;
334
335 break;
336 case 060: /* ....##, <noise^2> = 0 */
337 val = 0.5370 * out2_1 + 0.4630 * out2_2;
338 out[badX0] = (val < min) ? out2_1 : val;
339
340 val = 0.6968 * out2_1 + 0.3032 * out2_2;
341 out[badX0 + 1] = (val < min) ? out2_1 : val;
342
343 val = 1.0933 * out2_1 - 0.0933 * out2_2;
344 out[badX1 - 1] = (val < min) ? out2_1 : val;
345
346 val = 1.4288 * out2_1 - 0.4288 * out2_2;
347 out[badX1] = (val < min) ? out2_1 : val;
348
349 break;
350 case 020: /* ....#?, <noise^2> = 0 */
351 val = 1.0000 * out2_1;
352
353 out[badX0] = (val < min) ? out2_1 : val;
354 out[badX0 + 1] = (val < min) ? out2_1 : val;
355 out[badX1 - 1] = (val < min) ? out2_1 : val;
356 out[badX1] = (val < min) ? out2_1 : val;
357
358 break;
359 case 0140: /* .....##, <noise^2> = 0 */
360 val = 0.5041 * out2_1 + 0.4959 * out2_2;
361 out[badX0] = (val < min) ? out2_1 : val;
362
363 val = 0.5370 * out2_1 + 0.4630 * out2_2;
364 out[badX0 + 1] = (val < min) ? out2_1 : val;
365
366 val = 0.6968 * out2_1 + 0.3032 * out2_2;
367 out[badX1 - 2] = (val < min) ? out2_1 : val;
368
369 val = 1.0933 * out2_1 - 0.0933 * out2_2;
370 out[badX1 - 1] = (val < min) ? out2_1 : val;
371
372 val = 1.4288 * out2_1 - 0.4288 * out2_2;
373 out[badX1] = (val < min) ? out2_1 : val;
374
375 break;
376 case 040: /* .....#?, <noise^2> = 0 */
377 val = 1.0000 * out2_1;
378 out[badX0] = (val < min) ? out2_1 : val;
379 out[badX0 + 1] = (val < min) ? out2_1 : val;
380 out[badX1 - 2] = (val < min) ? out2_1 : val;
381 out[badX1 - 1] = (val < min) ? out2_1 : val;
382 out[badX1] = (val < min) ? out2_1 : val;
383
384 break;
385 case 0300: /* ......##, <noise^2> = 0 */
386 val = 0.5003 * out2_1 + 0.4997 * out2_2;
387 out[badX0] = (val < min) ? out2_1 : val;
388
389 val = 0.5041 * out2_1 + 0.4959 * out2_2;
390 out[badX0 + 1] = (val < min) ? out2_1 : val;
391
392 val = 0.5370 * out2_1 + 0.4630 * out2_2;
393 out[badX0 + 2] = (val < min) ? out2_1 : val;
394
395 val = 0.6968 * out2_1 + 0.3032 * out2_2;
396 out[badX1 - 2] = (val < min) ? out2_1 : val;
397
398 val = 1.0933 * out2_1 - 0.0933 * out2_2;
399 out[badX1 - 1] = (val < min) ? out2_1 : val;
400
401 val = 1.4288 * out2_1 - 0.4288 * out2_2;
402 out[badX1] = (val < min) ? out2_1 : val;
403
404 break;
405 case 0100: /* ......#?, <noise^2> = 0 */
406 val = 1.0000 * out2_1;
407
408 out[badX0] = (val < min) ? out2_1 : val;
409 out[badX0 + 1] = (val < min) ? out2_1 : val;
410 out[badX0 + 2] = (val < min) ? out2_1 : val;
411 out[badX1 - 2] = (val < min) ? out2_1 : val;
412 out[badX1 - 1] = (val < min) ? out2_1 : val;
413 out[badX1] = (val < min) ? out2_1 : val;
414
415 break;
416 case 0600: /* .......##, <noise^2> = 0 */
417 val = 0.5000 * out2_1 + 0.5000 * out2_2;
418 out[badX0] = (val < min) ? out2_1 : val;
419
420 val = 0.5003 * out2_1 + 0.4997 * out2_2;
421 out[badX0 + 1] = (val < min) ? out2_1 : val;
422
423 val = 0.5041 * out2_1 + 0.4959 * out2_2;
424 out[badX0 + 2] = (val < min) ? out2_1 : val;
425
426 val = 0.5370 * out2_1 + 0.4630 * out2_2;
427 out[badX1 - 3] = (val < min) ? out2_1 : val;
428
429 val = 0.6968 * out2_1 + 0.3032 * out2_2;
430 out[badX1 - 2] = (val < min) ? out2_1 : val;
431
432 val = 1.0933 * out2_1 - 0.0933 * out2_2;
433 out[badX1 - 1] = (val < min) ? out2_1 : val;
434
435 val = 1.4288 * out2_1 - 0.4288 * out2_2;
436 out[badX1] = (val < min) ? out2_1 : val;
437
438 break;
439 case 0200: /* .......#?, <noise^2> = 0 */
440 val = 1.0000 * out2_1;
441 out[badX0] = (val < min) ? out2_1 : val;
442 out[badX0 + 1] = (val < min) ? out2_1 : val;
443 out[badX0 + 2] = (val < min) ? out2_1 : val;
444 out[badX1 - 3] = (val < min) ? out2_1 : val;
445 out[badX1 - 2] = (val < min) ? out2_1 : val;
446 out[badX1 - 1] = (val < min) ? out2_1 : val;
447 out[badX1] = (val < min) ? out2_1 : val;
448
449 break;
450 case 01400: /* ........##, <noise^2> = 0 */
451 val = 0.5000 * out2_1 + 0.5000 * out2_2;
452 out[badX0] = (val < min) ? out2_1 : val;
453
454 val = 0.5000 * out2_1 + 0.5000 * out2_2;
455 out[badX0 + 1] = (val < min) ? out2_1 : val;
456
457 val = 0.5003 * out2_1 + 0.4997 * out2_2;
458 out[badX0 + 2] = (val < min) ? out2_1 : val;
459
460 val = 0.5041 * out2_1 + 0.4959 * out2_2;
461 out[badX0 + 3] = (val < min) ? out2_1 : val;
462
463 val = 0.5370 * out2_1 + 0.4630 * out2_2;
464 out[badX1 - 3] = (val < min) ? out2_1 : val;
465
466 val = 0.6968 * out2_1 + 0.3032 * out2_2;
467 out[badX1 - 2] = (val < min) ? out2_1 : val;
468
469 val = 1.0933 * out2_1 - 0.0933 * out2_2;
470 out[badX1 - 1] = (val < min) ? out2_1 : val;
471
472 val = 1.4288 * out2_1 - 0.4288 * out2_2;
473 out[badX1] = (val < min) ? out2_1 : val;
474
475 break;
476 case 0400: /* ........#?, <noise^2> = 0 */
477 val = 1.0000 * out2_1;
478 out[badX0] = (val < min) ? out2_1 : val;
479 out[badX0 + 1] = (val < min) ? out2_1 : val;
480 out[badX0 + 2] = (val < min) ? out2_1 : val;
481 out[badX0 + 3] = (val < min) ? out2_1 : val;
482 out[badX1 - 3] = (val < min) ? out2_1 : val;
483 out[badX1 - 2] = (val < min) ? out2_1 : val;
484 out[badX1 - 1] = (val < min) ? out2_1 : val;
485 out[badX1] = (val < min) ? out2_1 : val;
486
487 break;
488 case 03000: /* .........##, <noise^2> = 0 */
489 val = 0.5000 * out2_1 + 0.5000 * out2_2;
490 out[badX0] = (val < min) ? out2_1 : val;
491
492 val = 0.5000 * out2_1 + 0.5000 * out2_2;
493 out[badX0 + 1] = (val < min) ? out2_1 : val;
494
495 val = 0.5000 * out2_1 + 0.5000 * out2_2;
496 out[badX0 + 2] = (val < min) ? out2_1 : val;
497
498 val = 0.5003 * out2_1 + 0.4997 * out2_2;
499 out[badX0 + 3] = (val < min) ? out2_1 : val;
500
501 val = 0.5041 * out2_1 + 0.4959 * out2_2;
502 out[badX1 - 4] = (val < min) ? out2_1 : val;
503
504 val = 0.5370 * out2_1 + 0.4630 * out2_2;
505 out[badX1 - 3] = (val < min) ? out2_1 : val;
506
507 val = 0.6968 * out2_1 + 0.3032 * out2_2;
508 out[badX1 - 2] = (val < min) ? out2_1 : val;
509
510 val = 1.0933 * out2_1 - 0.0933 * out2_2;
511 out[badX1 - 1] = (val < min) ? out2_1 : val;
512
513 val = 1.4288 * out2_1 - 0.4288 * out2_2;
514 out[badX1] = (val < min) ? out2_1 : val;
515
516 break;
517 case 01000: /* .........#?, <noise^2> = 0 */
518 val = 1.0000 * out2_1;
519 out[badX0] = (val < min) ? out2_1 : val;
520 out[badX0 + 1] = (val < min) ? out2_1 : val;
521 out[badX0 + 2] = (val < min) ? out2_1 : val;
522 out[badX0 + 3] = (val < min) ? out2_1 : val;
523 out[badX1 - 4] = (val < min) ? out2_1 : val;
524 out[badX1 - 3] = (val < min) ? out2_1 : val;
525 out[badX1 - 2] = (val < min) ? out2_1 : val;
526 out[badX1 - 1] = (val < min) ? out2_1 : val;
527 out[badX1] = (val < min) ? out2_1 : val;
528
529 break;
530 case 06000: /* ..........##, <noise^2> = 0 */
531 val = 0.5000 * out2_1 + 0.5000 * out2_2;
532 out[badX0] = (val < min) ? out2_1 : val;
533
534 val = 0.5000 * out2_1 + 0.5000 * out2_2;
535 out[badX0 + 1] = (val < min) ? out2_1 : val;
536
537 val = 0.5000 * out2_1 + 0.5000 * out2_2;
538 out[badX0 + 2] = (val < min) ? out2_1 : val;
539
540 val = 0.5000 * out2_1 + 0.5000 * out2_2;
541 out[badX0 + 3] = (val < min) ? out2_1 : val;
542
543 val = 0.5003 * out2_1 + 0.4997 * out2_2;
544 out[badX0 + 4] = (val < min) ? out2_1 : val;
545
546 val = 0.5041 * out2_1 + 0.4959 * out2_2;
547 out[badX1 - 4] = (val < min) ? out2_1 : val;
548
549 val = 0.5370 * out2_1 + 0.4630 * out2_2;
550 out[badX1 - 3] = (val < min) ? out2_1 : val;
551
552 val = 0.6968 * out2_1 + 0.3032 * out2_2;
553 out[badX1 - 2] = (val < min) ? out2_1 : val;
554
555 val = 1.0933 * out2_1 - 0.0933 * out2_2;
556 out[badX1 - 1] = (val < min) ? out2_1 : val;
557
558 val = 1.4288 * out2_1 - 0.4288 * out2_2;
559 out[badX1] = (val < min) ? out2_1 : val;
560
561 break;
562 case 02000: /* ..........#?, <noise^2> = 0 */
563 val = 1.0000 * out2_1;
564
565 out[badX0] = (val < min) ? out2_1 : val;
566 out[badX0 + 1] = (val < min) ? out2_1 : val;
567 out[badX0 + 2] = (val < min) ? out2_1 : val;
568 out[badX0 + 3] = (val < min) ? out2_1 : val;
569 out[badX0 + 4] = (val < min) ? out2_1 : val;
570 out[badX1 - 4] = (val < min) ? out2_1 : val;
571 out[badX1 - 3] = (val < min) ? out2_1 : val;
572 out[badX1 - 2] = (val < min) ? out2_1 : val;
573 out[badX1 - 1] = (val < min) ? out2_1 : val;
574 out[badX1] = (val < min) ? out2_1 : val;
575
576 break;
577 default:
578 // shFatal("Unsupported defect type: LEFT 0%o", defectType);
579 break; /* NOTREACHED */
580 }
581 break;
582 case Defect::WIDE_LEFT:
583 assert(badX0 >= 0);
584 if (badX1 + 2 >= ncol) { /* left defect extends near
585 right edge of data! */
586 if (badX1 == ncol - 2) { /* one column remains */
587 val = out[ncol - 1];
588 } else {
589 val = fallbackValue; /* there is no information */
590 }
591 for (int j = badX0; j <= badX1; j++) {
592 out[j] = val;
593 }
594 break;
595 }
596 out2_1 = out[badX1 + 1];
597 out2_2 = out[badX1 + 2];
598
599 switch (defectType) {
600 case 02: /* ?#., <noise^2> = 0 */
601 val = 1.0000 * out2_1;
602 val = (val < min) ? out2_1 : val;
603
604 for (int j = badX0; j <= badX1; j++) {
605 out[j] = val;
606 }
607 break;
608 case 03: /* ?##, <noise^2> = 0 */
609 val = 0.5000 * out2_1 + 0.5000 * out2_2;
610 if (val < min) {
611 val = out2_1;
612 }
613
614 for (int j = badX0; j < badX1 - 5; j++) {
615 out[j] = val;
616 }
617
618 val = 0.5003 * out2_1 + 0.4997 * out2_2;
619 out[badX1 - 5] = (val < min) ? out2_1 : val;
620
621 val = 0.5041 * out2_1 + 0.4959 * out2_2;
622 out[badX1 - 4] = (val < min) ? out2_1 : val;
623
624 val = 0.5370 * out2_1 + 0.4630 * out2_2;
625 out[badX1 - 3] = (val < min) ? out2_1 : val;
626
627 val = 0.6968 * out2_1 + 0.3032 * out2_2;
628 out[badX1 - 2] = (val < min) ? out2_1 : val;
629
630 val = 1.0933 * out2_1 - 0.0933 * out2_2;
631 out[badX1 - 1] = (val < min) ? out2_1 : val;
632
633 val = 1.4288 * out2_1 - 0.4288 * out2_2;
634 out[badX1] = (val < min) ? out2_1 : val;
635
636 break;
637 default:
638 // shFatal("Unsupported defect type: WIDE_LEFT 0%o",defect[i].type);
639 break; /* NOTREACHED */
640 }
641
642 break;
643 case Defect::RIGHT:
644 assert(badX0 >= 2 && badX1 < ncol);
645
646 out1_2 = out[badX0 - 2];
647 out1_1 = out[badX0 - 1];
648
649 switch (defectType) {
650 case 06: /* ##., <noise^2> = 0 */
651 val = -0.4288 * out1_2 + 1.4288 * out1_1;
652 out[badX1] = (val < min) ? out1_1 : val;
653
654 break;
655 case 014: /* ##.., <noise^2> = 0 */
656 val = -0.4288 * out1_2 + 1.4288 * out1_1;
657 out[badX0] = (val < min) ? out1_1 : val;
658
659 val = -0.0933 * out1_2 + 1.0933 * out1_1;
660 out[badX1] = (val < min) ? out1_1 : val;
661
662 break;
663 case 030: /* ##..., <noise^2> = 0 */
664 val = -0.4288 * out1_2 + 1.4288 * out1_1;
665 out[badX0] = (val < min) ? out1_1 : val;
666
667 val = -0.0933 * out1_2 + 1.0933 * out1_1;
668 out[badX1 - 1] = (val < min) ? out1_1 : val;
669
670 val = 0.3032 * out1_2 + 0.6968 * out1_1;
671 out[badX1] = (val < min) ? out1_1 : val;
672
673 break;
674 case 060: /* ##...., <noise^2> = 0 */
675 val = -0.4288 * out1_2 + 1.4288 * out1_1;
676 out[badX0] = (val < min) ? out1_1 : val;
677
678 val = -0.0933 * out1_2 + 1.0933 * out1_1;
679 out[badX0 + 1] = (val < min) ? out1_1 : val;
680
681 val = 0.3032 * out1_2 + 0.6968 * out1_1;
682 out[badX1 - 1] = (val < min) ? out1_1 : val;
683
684 val = 0.4630 * out1_2 + 0.5370 * out1_1;
685 out[badX1] = (val < min) ? out1_1 : val;
686
687 break;
688 case 0140: /* ##....., <noise^2> = 0 */
689 val = -0.4288 * out1_2 + 1.4288 * out1_1;
690 out[badX0] = (val < min) ? out1_1 : val;
691
692 val = -0.0933 * out1_2 + 1.0933 * out1_1;
693 out[badX0 + 1] = (val < min) ? out1_1 : val;
694
695 val = 0.3032 * out1_2 + 0.6968 * out1_1;
696 out[badX1 - 2] = (val < min) ? out1_1 : val;
697
698 val = 0.4630 * out1_2 + 0.5370 * out1_1;
699 out[badX1 - 1] = (val < min) ? out1_1 : val;
700
701 val = 0.4959 * out1_2 + 0.5041 * out1_1;
702 out[badX1] = (val < min) ? out1_1 : val;
703
704 break;
705 case 0300: /* ##......, <noise^2> = 0 */
706 val = -0.4288 * out1_2 + 1.4288 * out1_1;
707 out[badX0] = (val < min) ? out1_1 : val;
708
709 val = -0.0933 * out1_2 + 1.0933 * out1_1;
710 out[badX0 + 1] = (val < min) ? out1_1 : val;
711
712 val = 0.3032 * out1_2 + 0.6968 * out1_1;
713 out[badX0 + 2] = (val < min) ? out1_1 : val;
714
715 val = 0.4630 * out1_2 + 0.5370 * out1_1;
716 out[badX1 - 2] = (val < min) ? out1_1 : val;
717
718 val = 0.4959 * out1_2 + 0.5041 * out1_1;
719 out[badX1 - 1] = (val < min) ? out1_1 : val;
720
721 val = 0.4997 * out1_2 + 0.5003 * out1_1;
722 out[badX1] = (val < min) ? out1_1 : val;
723
724 break;
725 case 0600: /* ##......., <noise^2> = 0 */
726 val = -0.4288 * out1_2 + 1.4288 * out1_1;
727 out[badX0] = (val < min) ? out1_1 : val;
728
729 val = -0.0933 * out1_2 + 1.0933 * out1_1;
730 out[badX0 + 1] = (val < min) ? out1_1 : val;
731
732 val = 0.3032 * out1_2 + 0.6968 * out1_1;
733 out[badX0 + 2] = (val < min) ? out1_1 : val;
734
735 val = 0.4630 * out1_2 + 0.5370 * out1_1;
736 out[badX1 - 3] = (val < min) ? out1_1 : val;
737
738 val = 0.4959 * out1_2 + 0.5041 * out1_1;
739 out[badX1 - 2] = (val < min) ? out1_1 : val;
740
741 val = 0.4997 * out1_2 + 0.5003 * out1_1;
742 out[badX1 - 1] = (val < min) ? out1_1 : val;
743
744 val = 0.5000 * out1_2 + 0.5000 * out1_1;
745 out[badX1] = (val < min) ? out1_1 : val;
746
747 break;
748 case 01400: /* ##........, <noise^2> = 0 */
749 val = -0.4288 * out1_2 + 1.4288 * out1_1;
750 out[badX0] = (val < min) ? out1_1 : val;
751
752 val = -0.0933 * out1_2 + 1.0933 * out1_1;
753 out[badX0 + 1] = (val < min) ? out1_1 : val;
754
755 val = 0.3032 * out1_2 + 0.6968 * out1_1;
756 out[badX0 + 2] = (val < min) ? out1_1 : val;
757
758 val = 0.4630 * out1_2 + 0.5370 * out1_1;
759 out[badX0 + 3] = (val < min) ? out1_1 : val;
760
761 val = 0.4959 * out1_2 + 0.5041 * out1_1;
762 out[badX1 - 3] = (val < min) ? out1_1 : val;
763
764 val = 0.4997 * out1_2 + 0.5003 * out1_1;
765 out[badX1 - 2] = (val < min) ? out1_1 : val;
766
767 val = 0.5000 * out1_2 + 0.5000 * out1_1;
768 out[badX1 - 1] = (val < min) ? out1_1 : val;
769
770 val = 0.5000 * out1_2 + 0.5000 * out1_1;
771 out[badX1] = (val < min) ? out1_1 : val;
772
773 break;
774 case 03000: /* ##........., <noise^2> = 0 */
775 val = -0.4288 * out1_2 + 1.4288 * out1_1;
776 out[badX0] = (val < min) ? out1_1 : val;
777
778 val = -0.0933 * out1_2 + 1.0933 * out1_1;
779 out[badX0 + 1] = (val < min) ? out1_1 : val;
780
781 val = 0.3032 * out1_2 + 0.6968 * out1_1;
782 out[badX0 + 2] = (val < min) ? out1_1 : val;
783
784 val = 0.4630 * out1_2 + 0.5370 * out1_1;
785 out[badX0 + 3] = (val < min) ? out1_1 : val;
786
787 val = 0.4959 * out1_2 + 0.5041 * out1_1;
788 out[badX1 - 4] = (val < min) ? out1_1 : val;
789
790 val = 0.4997 * out1_2 + 0.5003 * out1_1;
791 out[badX1 - 3] = (val < min) ? out1_1 : val;
792
793 val = 0.5000 * out1_2 + 0.5000 * out1_1;
794 out[badX1 - 2] = (val < min) ? out1_1 : val;
795
796 val = 0.5000 * out1_2 + 0.5000 * out1_1;
797 out[badX1 - 1] = (val < min) ? out1_1 : val;
798
799 val = 0.5000 * out1_2 + 0.5000 * out1_1;
800 out[badX1] = (val < min) ? out1_1 : val;
801
802 break;
803 case 06000: /* ##.........., <noise^2> = 0 */
804 val = -0.4288 * out1_2 + 1.4288 * out1_1;
805 out[badX0] = (val < min) ? out1_1 : val;
806
807 val = -0.0933 * out1_2 + 1.0933 * out1_1;
808 out[badX0 + 1] = (val < min) ? out1_1 : val;
809
810 val = 0.3032 * out1_2 + 0.6968 * out1_1;
811 out[badX0 + 2] = (val < min) ? out1_1 : val;
812
813 val = 0.4630 * out1_2 + 0.5370 * out1_1;
814 out[badX0 + 3] = (val < min) ? out1_1 : val;
815
816 val = 0.4959 * out1_2 + 0.5041 * out1_1;
817 out[badX0 + 4] = (val < min) ? out1_1 : val;
818
819 val = 0.4997 * out1_2 + 0.5003 * out1_1;
820 out[badX1 - 4] = (val < min) ? out1_1 : val;
821
822 val = 0.5000 * out1_2 + 0.5000 * out1_1;
823 out[badX1 - 3] = (val < min) ? out1_1 : val;
824
825 val = 0.5000 * out1_2 + 0.5000 * out1_1;
826 out[badX1 - 2] = (val < min) ? out1_1 : val;
827
828 val = 0.5000 * out1_2 + 0.5000 * out1_1;
829 out[badX1 - 1] = (val < min) ? out1_1 : val;
830
831 val = 0.5000 * out1_2 + 0.5000 * out1_1;
832 out[badX1] = (val < min) ? out1_1 : val;
833
834 break;
835 default:
836 // shFatal("Unsupported defect type: RIGHT 0%o",defect[i].type);
837 break; /* NOTREACHED */
838 }
839 break;
840 case Defect::WIDE_RIGHT:
841 assert(badX1 < ncol);
842
843 if (badX0 < 2) { /* right defect extends near
844 left edge of data! */
845 if (badX0 == 1) { /* one column remains */
846 val = out[0];
847 } else {
848 val = fallbackValue; /* there is no information */
849 }
850 for (int j = badX0; j <= badX1; j++) {
851 out[j] = val;
852 }
853 break;
854 }
855
856 out1_2 = out[badX0 - 2];
857 out1_1 = out[badX0 - 1];
858
859 switch (defectType) {
860 case 03: /* ##?, S/N = infty */
861 val = -0.4288 * out1_2 + 1.4288 * out1_1;
862 out[badX0] = (val < min) ? out1_1 : val;
863
864 val = -0.0933 * out1_2 + 1.0933 * out1_1;
865 out[badX0 + 1] = (val < min) ? out1_1 : val;
866
867 val = 0.3032 * out1_2 + 0.6968 * out1_1;
868 out[badX0 + 2] = (val < min) ? out1_1 : val;
869
870 val = 0.4630 * out1_2 + 0.5370 * out1_1;
871 out[badX0 + 3] = (val < min) ? out1_1 : val;
872
873 val = 0.4959 * out1_2 + 0.5041 * out1_1;
874 out[badX0 + 4] = (val < min) ? out1_1 : val;
875
876 val = 0.4997 * out1_2 + 0.5003 * out1_1;
877 out[badX0 + 5] = (val < min) ? out1_1 : val;
878
879 val = 0.5000 * out1_2 + 0.5000 * out1_1;
880 val = (val < min) ? out1_1 : val;
881
882 for (int j = badX0 + 6; j <= badX1; j++) {
883 out[j] = val;
884 }
885 break;
886 default:
887 // shFatal("Unsupported defect type: WIDE_RIGHT 0%o",defect[i].type);
888 break; /* NOTREACHED */
889 }
890 break;
891 case Defect::MIDDLE:
892 case Defect::NEAR_LEFT:
893 case Defect::NEAR_RIGHT:
894 if (defectPos == Defect::MIDDLE) {
895 assert(badX0 >= 2 && badX1 + 2 < ncol);
896 out1_2 = out[badX0 - 2];
897 out2_2 = out[badX1 + 2];
898 } else if (defectPos == Defect::NEAR_LEFT) {
899 assert(badX0 >= 1 && badX1 + 2 < ncol);
900 out1_2 = -1; /* NOTUSED */
901 out2_2 = out[badX1 + 2];
902 } else if (defectPos == Defect::NEAR_RIGHT) {
903 assert(badX0 >= 2 && badX1 + 1 < ncol);
904 out1_2 = out[badX0 - 2];
905 out2_2 = -1; /* NOTUSED */
906 } else {
907 // shFatal("Unknown defect classification %d (%s:%d)",defectPos, __FILE__,__LINE__);
908 out1_2 = out2_2 = -1; /* NOTUSED */
909 }
910 out1_1 = out[badX0 - 1];
911 out2_1 = out[badX1 + 1];
912
913 switch (defectType) {
914 case 012: /* #.#., <noise^2> = 0, sigma = 1 */
915 val = 0.5000 * out1_1 + 0.5000 * out2_1;
916 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
917
918 break;
919 case 013: /* #.##, <noise^2> = 0 */
920 val = 0.4875 * out1_1 + 0.8959 * out2_1 - 0.3834 * out2_2;
921 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
922
923 break;
924 case 022: /* #..#., <noise^2> = 0, sigma = 1 */
925 val = 0.7297 * out1_1 + 0.2703 * out2_1;
926 out[badX0] = (val < 0) ? 0 : val;
927
928 val = 0.2703 * out1_1 + 0.7297 * out2_1;
929 out[badX1] = (val < 0) ? 0 : val;
930
931 break;
932 case 023: /* #..##, <noise^2> = 0 */
933 val = 0.7538 * out1_1 + 0.5680 * out2_1 - 0.3218 * out2_2;
934 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
935
936 val = 0.3095 * out1_1 + 1.2132 * out2_1 - 0.5227 * out2_2;
937 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
938
939 break;
940 case 032: /* ##.#., <noise^2> = 0 */
941 val = -0.3834 * out1_2 + 0.8959 * out1_1 + 0.4875 * out2_1;
942 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
943
944 break;
945 case 033: /* ##.##, <noise^2> = 0 */
946 /* These coefficients are also available as
947 interp::interp_1_c1 and interp::interp_1_c2 */
948 val = -0.2737 * out1_2 + 0.7737 * out1_1 + 0.7737 * out2_1 - 0.2737 * out2_2;
949 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
950
951 break;
952 case 042: /* #...#., <noise^2> = 0, sigma = 1 */
953 val = 0.8430 * out1_1 + 0.1570 * out2_1;
954 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
955
956 val = 0.5000 * out1_1 + 0.5000 * out2_1;
957 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
958
959 val = 0.1570 * out1_1 + 0.8430 * out2_1;
960 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
961
962 break;
963 case 043: /* #...##, <noise^2> = 0 */
964 val = 0.8525 * out1_1 + 0.2390 * out2_1 - 0.0915 * out2_2;
965 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
966
967 val = 0.5356 * out1_1 + 0.8057 * out2_1 - 0.3413 * out2_2;
968 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
969
970 val = 0.2120 * out1_1 + 1.3150 * out2_1 - 0.5270 * out2_2;
971 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
972
973 break;
974 case 062: /* ##..#., <noise^2> = 0 */
975 val = -0.5227 * out1_2 + 1.2132 * out1_1 + 0.3095 * out2_1;
976 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
977
978 val = -0.3218 * out1_2 + 0.5680 * out1_1 + 0.7538 * out2_1;
979 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
980
981 break;
982 case 063: /* ##..##, <noise^2> = 0 */
983 val = -0.4793 * out1_2 + 1.1904 * out1_1 + 0.5212 * out2_1 - 0.2323 * out2_2;
984 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
985
986 val = -0.2323 * out1_2 + 0.5212 * out1_1 + 1.1904 * out2_1 - 0.4793 * out2_2;
987 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
988
989 break;
990 case 0102: /* #....#., <noise^2> = 0, sigma = 1 */
991 val = 0.8810 * out1_1 + 0.1190 * out2_1;
992 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
993
994 val = 0.6315 * out1_1 + 0.3685 * out2_1;
995 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
996
997 val = 0.3685 * out1_1 + 0.6315 * out2_1;
998 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
999
1000 val = 0.1190 * out1_1 + 0.8810 * out2_1;
1001 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1002
1003 break;
1004 case 0103: /* #....##, <noise^2> = 0 */
1005 val = 0.8779 * out1_1 + 0.0945 * out2_1 + 0.0276 * out2_2;
1006 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1007
1008 val = 0.6327 * out1_1 + 0.3779 * out2_1 - 0.0106 * out2_2;
1009 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1010
1011 val = 0.4006 * out1_1 + 0.8914 * out2_1 - 0.2920 * out2_2;
1012 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1013
1014 val = 0.1757 * out1_1 + 1.3403 * out2_1 - 0.5160 * out2_2;
1015 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1016
1017 break;
1018 case 0142: /* ##...#., <noise^2> = 0 */
1019 val = -0.5270 * out1_2 + 1.3150 * out1_1 + 0.2120 * out2_1;
1020 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1021
1022 val = -0.3413 * out1_2 + 0.8057 * out1_1 + 0.5356 * out2_1;
1023 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1024
1025 val = -0.0915 * out1_2 + 0.2390 * out1_1 + 0.8525 * out2_1;
1026 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1027
1028 break;
1029 case 0143: /* ##...##, <noise^2> = 0 */
1030 val = -0.5230 * out1_2 + 1.3163 * out1_1 + 0.2536 * out2_1 - 0.0469 * out2_2;
1031 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1032
1033 val = -0.3144 * out1_2 + 0.8144 * out1_1 + 0.8144 * out2_1 - 0.3144 * out2_2;
1034 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1035
1036 val = -0.0469 * out1_2 + 0.2536 * out1_1 + 1.3163 * out2_1 - 0.5230 * out2_2;
1037 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1038
1039 break;
1040 case 0202: /* #.....#., <noise^2> = 0, sigma = 1 */
1041 val = 0.8885 * out1_1 + 0.1115 * out2_1;
1042 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1043
1044 val = 0.6748 * out1_1 + 0.3252 * out2_1;
1045 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1046
1047 val = 0.5000 * out1_1 + 0.5000 * out2_1;
1048 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1049
1050 val = 0.3252 * out1_1 + 0.6748 * out2_1;
1051 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1052
1053 val = 0.1115 * out1_1 + 0.8885 * out2_1;
1054 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1055
1056 break;
1057 case 0203: /* #.....##, <noise^2> = 0 */
1058 val = 0.8824 * out1_1 + 0.0626 * out2_1 + 0.0549 * out2_2;
1059 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1060
1061 val = 0.6601 * out1_1 + 0.2068 * out2_1 + 0.1331 * out2_2;
1062 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1063
1064 val = 0.4938 * out1_1 + 0.4498 * out2_1 + 0.0564 * out2_2;
1065 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1066
1067 val = 0.3551 * out1_1 + 0.9157 * out2_1 - 0.2708 * out2_2;
1068 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1069
1070 val = 0.1682 * out1_1 + 1.3447 * out2_1 - 0.5129 * out2_2;
1071 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1072
1073 break;
1074 case 0302: /* ##....#., <noise^2> = 0 */
1075 val = -0.5160 * out1_2 + 1.3403 * out1_1 + 0.1757 * out2_1;
1076 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1077
1078 val = -0.2920 * out1_2 + 0.8914 * out1_1 + 0.4006 * out2_1;
1079 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1080
1081 val = -0.0106 * out1_2 + 0.3779 * out1_1 + 0.6327 * out2_1;
1082 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1083
1084 val = 0.0276 * out1_2 + 0.0945 * out1_1 + 0.8779 * out2_1;
1085 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1086
1087 break;
1088 case 0303: /* ##....##, <noise^2> = 0 */
1089 val = -0.5197 * out1_2 + 1.3370 * out1_1 + 0.1231 * out2_1 + 0.0596 * out2_2;
1090 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1091
1092 val = -0.2924 * out1_2 + 0.8910 * out1_1 + 0.3940 * out2_1 + 0.0074 * out2_2;
1093 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1094
1095 val = 0.0074 * out1_2 + 0.3940 * out1_1 + 0.8910 * out2_1 - 0.2924 * out2_2;
1096 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1097
1098 val = 0.0596 * out1_2 + 0.1231 * out1_1 + 1.3370 * out2_1 - 0.5197 * out2_2;
1099 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1100
1101 break;
1102 case 0402: /* #......#., <noise^2> = 0, sigma = 1 */
1103 val = 0.8893 * out1_1 + 0.1107 * out2_1;
1104 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1105
1106 val = 0.6830 * out1_1 + 0.3170 * out2_1;
1107 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1108
1109 val = 0.5435 * out1_1 + 0.4565 * out2_1;
1110 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1111
1112 val = 0.4565 * out1_1 + 0.5435 * out2_1;
1113 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1114
1115 val = 0.3170 * out1_1 + 0.6830 * out2_1;
1116 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1117
1118 val = 0.1107 * out1_1 + 0.8893 * out2_1;
1119 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1120
1121 break;
1122 case 0403: /* #......##, <noise^2> = 0 */
1123 val = 0.8829 * out1_1 + 0.0588 * out2_1 + 0.0583 * out2_2;
1124 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1125
1126 val = 0.6649 * out1_1 + 0.1716 * out2_1 + 0.1635 * out2_2;
1127 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1128
1129 val = 0.5212 * out1_1 + 0.2765 * out2_1 + 0.2024 * out2_2;
1130 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1131
1132 val = 0.4477 * out1_1 + 0.4730 * out2_1 + 0.0793 * out2_2;
1133 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1134
1135 val = 0.3465 * out1_1 + 0.9201 * out2_1 - 0.2666 * out2_2;
1136 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1137
1138 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1139 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1140
1141 break;
1142 case 0602: /* ##.....#., <noise^2> = 0 */
1143 val = -0.5129 * out1_2 + 1.3447 * out1_1 + 0.1682 * out2_1;
1144 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1145
1146 val = -0.2708 * out1_2 + 0.9157 * out1_1 + 0.3551 * out2_1;
1147 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1148
1149 val = 0.0564 * out1_2 + 0.4498 * out1_1 + 0.4938 * out2_1;
1150 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1151
1152 val = 0.1331 * out1_2 + 0.2068 * out1_1 + 0.6601 * out2_1;
1153 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1154
1155 val = 0.0549 * out1_2 + 0.0626 * out1_1 + 0.8824 * out2_1;
1156 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1157
1158 break;
1159 case 0603: /* ##.....##, <noise^2> = 0 */
1160 val = -0.5179 * out1_2 + 1.3397 * out1_1 + 0.0928 * out2_1 + 0.0854 * out2_2;
1161 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1162
1163 val = -0.2796 * out1_2 + 0.9069 * out1_1 + 0.2231 * out2_1 + 0.1495 * out2_2;
1164 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1165
1166 val = 0.0533 * out1_2 + 0.4467 * out1_1 + 0.4467 * out2_1 + 0.0533 * out2_2;
1167 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1168
1169 val = 0.1495 * out1_2 + 0.2231 * out1_1 + 0.9069 * out2_1 - 0.2796 * out2_2;
1170 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1171
1172 val = 0.0854 * out1_2 + 0.0928 * out1_1 + 1.3397 * out2_1 - 0.5179 * out2_2;
1173 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1174
1175 break;
1176 case 01002: /* #.......#., <noise^2> = 0, sigma = 1 */
1177 val = 0.8894 * out1_1 + 0.1106 * out2_1;
1178 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1179
1180 val = 0.6839 * out1_1 + 0.3161 * out2_1;
1181 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1182
1183 val = 0.5517 * out1_1 + 0.4483 * out2_1;
1184 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1185
1186 val = 0.5000 * out1_1 + 0.5000 * out2_1;
1187 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1188
1189 val = 0.4483 * out1_1 + 0.5517 * out2_1;
1190 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1191
1192 val = 0.3161 * out1_1 + 0.6839 * out2_1;
1193 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1194
1195 val = 0.1106 * out1_1 + 0.8894 * out2_1;
1196 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1197
1198 break;
1199 case 01003: /* #.......##, <noise^2> = 0 */
1200 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1201 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1202
1203 val = 0.6654 * out1_1 + 0.1676 * out2_1 + 0.1670 * out2_2;
1204 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1205
1206 val = 0.5260 * out1_1 + 0.2411 * out2_1 + 0.2329 * out2_2;
1207 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1208
1209 val = 0.4751 * out1_1 + 0.2995 * out2_1 + 0.2254 * out2_2;
1210 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1211
1212 val = 0.4390 * out1_1 + 0.4773 * out2_1 + 0.0836 * out2_2;
1213 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1214
1215 val = 0.3456 * out1_1 + 0.9205 * out2_1 - 0.2661 * out2_2;
1216 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1217
1218 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1219 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1220
1221 break;
1222 case 01402: /* ##......#., <noise^2> = 0 */
1223 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1224 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1225
1226 val = -0.2666 * out1_2 + 0.9201 * out1_1 + 0.3465 * out2_1;
1227 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1228
1229 val = 0.0793 * out1_2 + 0.4730 * out1_1 + 0.4477 * out2_1;
1230 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1231
1232 val = 0.2024 * out1_2 + 0.2765 * out1_1 + 0.5212 * out2_1;
1233 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1234
1235 val = 0.1635 * out1_2 + 0.1716 * out1_1 + 0.6649 * out2_1;
1236 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1237
1238 val = 0.0583 * out1_2 + 0.0588 * out1_1 + 0.8829 * out2_1;
1239 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1240
1241 break;
1242 case 01403: /* ##......##, <noise^2> = 0 */
1243 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0891 * out2_1 + 0.0886 * out2_2;
1244 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1245
1246 val = -0.2771 * out1_2 + 0.9095 * out1_1 + 0.1878 * out2_1 + 0.1797 * out2_2;
1247 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1248
1249 val = 0.0677 * out1_2 + 0.4614 * out1_1 + 0.2725 * out2_1 + 0.1984 * out2_2;
1250 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1251
1252 val = 0.1984 * out1_2 + 0.2725 * out1_1 + 0.4614 * out2_1 + 0.0677 * out2_2;
1253 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1254
1255 val = 0.1797 * out1_2 + 0.1878 * out1_1 + 0.9095 * out2_1 - 0.2771 * out2_2;
1256 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1257
1258 val = 0.0886 * out1_2 + 0.0891 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1259 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1260
1261 break;
1262 case 02002: /* #........#., <noise^2> = 0, sigma = 1 */
1263 val = 0.8894 * out1_1 + 0.1106 * out2_1;
1264 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1265
1266 val = 0.6839 * out1_1 + 0.3161 * out2_1;
1267 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1268
1269 val = 0.5526 * out1_1 + 0.4474 * out2_1;
1270 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1271
1272 val = 0.5082 * out1_1 + 0.4918 * out2_1;
1273 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1274
1275 val = 0.4918 * out1_1 + 0.5082 * out2_1;
1276 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1277
1278 val = 0.4474 * out1_1 + 0.5526 * out2_1;
1279 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1280
1281 val = 0.3161 * out1_1 + 0.6839 * out2_1;
1282 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1283
1284 val = 0.1106 * out1_1 + 0.8894 * out2_1;
1285 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1286
1287 break;
1288 case 02003: /* #........##, <noise^2> = 0 */
1289 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1290 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1291
1292 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1293 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1294
1295 val = 0.5265 * out1_1 + 0.2370 * out2_1 + 0.2365 * out2_2;
1296 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1297
1298 val = 0.4799 * out1_1 + 0.2641 * out2_1 + 0.2560 * out2_2;
1299 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1300
1301 val = 0.4664 * out1_1 + 0.3038 * out2_1 + 0.2298 * out2_2;
1302 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1303
1304 val = 0.4381 * out1_1 + 0.4778 * out2_1 + 0.0841 * out2_2;
1305 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1306
1307 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1308 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1309
1310 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1311 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1312
1313 break;
1314 case 03002: /* ##.......#., <noise^2> = 0 */
1315 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1316 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1317
1318 val = -0.2661 * out1_2 + 0.9205 * out1_1 + 0.3456 * out2_1;
1319 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1320
1321 val = 0.0836 * out1_2 + 0.4773 * out1_1 + 0.4390 * out2_1;
1322 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1323
1324 val = 0.2254 * out1_2 + 0.2995 * out1_1 + 0.4751 * out2_1;
1325 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1326
1327 val = 0.2329 * out1_2 + 0.2411 * out1_1 + 0.5260 * out2_1;
1328 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1329
1330 val = 0.1670 * out1_2 + 0.1676 * out1_1 + 0.6654 * out2_1;
1331 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1332
1333 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1334 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1335
1336 break;
1337 case 03003: /* ##.......##, <noise^2> = 0 */
1338 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0889 * out2_1 + 0.0888 * out2_2;
1339 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1340
1341 val = -0.2768 * out1_2 + 0.9098 * out1_1 + 0.1838 * out2_1 + 0.1832 * out2_2;
1342 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1343
1344 val = 0.0703 * out1_2 + 0.4639 * out1_1 + 0.2370 * out2_1 + 0.2288 * out2_2;
1345 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1346
1347 val = 0.2130 * out1_2 + 0.2870 * out1_1 + 0.2870 * out2_1 + 0.2130 * out2_2;
1348 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1349
1350 val = 0.2288 * out1_2 + 0.2370 * out1_1 + 0.4639 * out2_1 + 0.0703 * out2_2;
1351 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1352
1353 val = 0.1832 * out1_2 + 0.1838 * out1_1 + 0.9098 * out2_1 - 0.2768 * out2_2;
1354 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1355
1356 val = 0.0888 * out1_2 + 0.0889 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1357 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1358
1359 break;
1360 case 04002: /* #.........#., <noise^2> = 0 */
1361 val = 0.8894 * out1_1 + 0.1106 * out2_1;
1362 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1363
1364 val = 0.6839 * out1_1 + 0.3161 * out2_1;
1365 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1366
1367 val = 0.5527 * out1_1 + 0.4473 * out2_1;
1368 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1369
1370 val = 0.5091 * out1_1 + 0.4909 * out2_1;
1371 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1372
1373 val = 0.5000 * out1_1 + 0.5000 * out2_1;
1374 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1375
1376 val = 0.4909 * out1_1 + 0.5091 * out2_1;
1377 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1378
1379 val = 0.4473 * out1_1 + 0.5527 * out2_1;
1380 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1381
1382 val = 0.3161 * out1_1 + 0.6839 * out2_1;
1383 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1384
1385 val = 0.1106 * out1_1 + 0.8894 * out2_1;
1386 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1387
1388 break;
1389 case 04003: /* #.........##, <noise^2> = 0 */
1390 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1391 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1392
1393 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1394 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1395
1396 val = 0.5265 * out1_1 + 0.2368 * out2_1 + 0.2367 * out2_2;
1397 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1398
1399 val = 0.4804 * out1_1 + 0.2601 * out2_1 + 0.2595 * out2_2;
1400 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1401
1402 val = 0.4712 * out1_1 + 0.2685 * out2_1 + 0.2603 * out2_2;
1403 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1404
1405 val = 0.4654 * out1_1 + 0.3043 * out2_1 + 0.2302 * out2_2;
1406 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1407
1408 val = 0.4380 * out1_1 + 0.4778 * out2_1 + 0.0842 * out2_2;
1409 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1410
1411 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1412 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1413
1414 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1415 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1416
1417 break;
1418 case 06002: /* ##........#., <noise^2> = 0 */
1419 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1420 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1421
1422 val = -0.2661 * out1_2 + 0.9206 * out1_1 + 0.3455 * out2_1;
1423 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1424
1425 val = 0.0841 * out1_2 + 0.4778 * out1_1 + 0.4381 * out2_1;
1426 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1427
1428 val = 0.2298 * out1_2 + 0.3038 * out1_1 + 0.4664 * out2_1;
1429 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1430
1431 val = 0.2560 * out1_2 + 0.2641 * out1_1 + 0.4799 * out2_1;
1432 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1433
1434 val = 0.2365 * out1_2 + 0.2370 * out1_1 + 0.5265 * out2_1;
1435 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1436
1437 val = 0.1673 * out1_2 + 0.1673 * out1_1 + 0.6654 * out2_1;
1438 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1439
1440 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1441 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1442
1443 break;
1444 case 06003: /* ##........##, <noise^2> = 0 */
1445 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0888 * out2_1 + 0.0888 * out2_2;
1446 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1447
1448 val = -0.2768 * out1_2 + 0.9098 * out1_1 + 0.1835 * out2_1 + 0.1835 * out2_2;
1449 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1450
1451 val = 0.0705 * out1_2 + 0.4642 * out1_1 + 0.2329 * out2_1 + 0.2324 * out2_2;
1452 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1453
1454 val = 0.2155 * out1_2 + 0.2896 * out1_1 + 0.2515 * out2_1 + 0.2434 * out2_2;
1455 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1456
1457 val = 0.2434 * out1_2 + 0.2515 * out1_1 + 0.2896 * out2_1 + 0.2155 * out2_2;
1458 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1459
1460 val = 0.2324 * out1_2 + 0.2329 * out1_1 + 0.4642 * out2_1 + 0.0705 * out2_2;
1461 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1462
1463 val = 0.1835 * out1_2 + 0.1835 * out1_1 + 0.9098 * out2_1 - 0.2768 * out2_2;
1464 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1465
1466 val = 0.0888 * out1_2 + 0.0888 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1467 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1468
1469 break;
1470 case 010002: /* #..........#., <noise^2> = 0, sigma = 1 */
1471 val = 0.8894 * out1_1 + 0.1106 * out2_1;
1472 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1473
1474 val = 0.6839 * out1_1 + 0.3161 * out2_1;
1475 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1476
1477 val = 0.5527 * out1_1 + 0.4473 * out2_1;
1478 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1479
1480 val = 0.5092 * out1_1 + 0.4908 * out2_1;
1481 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1482
1483 val = 0.5009 * out1_1 + 0.4991 * out2_1;
1484 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1485
1486 val = 0.4991 * out1_1 + 0.5009 * out2_1;
1487 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1488
1489 val = 0.4908 * out1_1 + 0.5092 * out2_1;
1490 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1491
1492 val = 0.4473 * out1_1 + 0.5527 * out2_1;
1493 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1494
1495 val = 0.3161 * out1_1 + 0.6839 * out2_1;
1496 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1497
1498 val = 0.1106 * out1_1 + 0.8894 * out2_1;
1499 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1500
1501 break;
1502 case 010003: /* #..........##, <noise^2> = 0 */
1503 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1504 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1505
1506 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1507 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1508
1509 val = 0.5265 * out1_1 + 0.2367 * out2_1 + 0.2367 * out2_2;
1510 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1511
1512 val = 0.4804 * out1_1 + 0.2598 * out2_1 + 0.2598 * out2_2;
1513 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1514
1515 val = 0.4717 * out1_1 + 0.2644 * out2_1 + 0.2639 * out2_2;
1516 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1517
1518 val = 0.4703 * out1_1 + 0.2690 * out2_1 + 0.2608 * out2_2;
1519 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1520
1521 val = 0.4654 * out1_1 + 0.3043 * out2_1 + 0.2303 * out2_2;
1522 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1523
1524 val = 0.4380 * out1_1 + 0.4778 * out2_1 + 0.0842 * out2_2;
1525 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1526
1527 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1528 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1529
1530 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1531 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1532
1533 break;
1534 case 014002: /* ##.........#., <noise^2> = 0 */
1535 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1536 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1537
1538 val = -0.2661 * out1_2 + 0.9206 * out1_1 + 0.3455 * out2_1;
1539 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1540
1541 val = 0.0842 * out1_2 + 0.4778 * out1_1 + 0.4380 * out2_1;
1542 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1543
1544 val = 0.2302 * out1_2 + 0.3043 * out1_1 + 0.4654 * out2_1;
1545 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1546
1547 val = 0.2603 * out1_2 + 0.2685 * out1_1 + 0.4712 * out2_1;
1548 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1549
1550 val = 0.2595 * out1_2 + 0.2601 * out1_1 + 0.4804 * out2_1;
1551 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1552
1553 val = 0.2367 * out1_2 + 0.2368 * out1_1 + 0.5265 * out2_1;
1554 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1555
1556 val = 0.1673 * out1_2 + 0.1673 * out1_1 + 0.6654 * out2_1;
1557 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1558
1559 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1560 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1561
1562 break;
1563 case 014003: /* ##.........##, <noise^2> = 0 */
1564 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0888 * out2_1 + 0.0888 * out2_2;
1565 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1566
1567 val = -0.2768 * out1_2 + 0.9098 * out1_1 + 0.1835 * out2_1 + 0.1835 * out2_2;
1568 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1569
1570 val = 0.0705 * out1_2 + 0.4642 * out1_1 + 0.2326 * out2_1 + 0.2326 * out2_2;
1571 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1572
1573 val = 0.2158 * out1_2 + 0.2899 * out1_1 + 0.2474 * out2_1 + 0.2469 * out2_2;
1574 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1575
1576 val = 0.2459 * out1_2 + 0.2541 * out1_1 + 0.2541 * out2_1 + 0.2459 * out2_2;
1577 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1578
1579 val = 0.2469 * out1_2 + 0.2474 * out1_1 + 0.2899 * out2_1 + 0.2158 * out2_2;
1580 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1581
1582 val = 0.2326 * out1_2 + 0.2326 * out1_1 + 0.4642 * out2_1 + 0.0705 * out2_2;
1583 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1584
1585 val = 0.1835 * out1_2 + 0.1835 * out1_1 + 0.9098 * out2_1 - 0.2768 * out2_2;
1586 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1587
1588 val = 0.0888 * out1_2 + 0.0888 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1589 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1590
1591 break;
1592 case 020003: /* #...........##, <noise^2> = 0 */
1593 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1594 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1595
1596 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1597 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1598
1599 val = 0.5265 * out1_1 + 0.2367 * out2_1 + 0.2367 * out2_2;
1600 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1601
1602 val = 0.4804 * out1_1 + 0.2598 * out2_1 + 0.2598 * out2_2;
1603 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1604
1605 val = 0.4718 * out1_1 + 0.2641 * out2_1 + 0.2641 * out2_2;
1606 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1607
1608 val = 0.4708 * out1_1 + 0.2649 * out2_1 + 0.2644 * out2_2;
1609 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1610
1611 val = 0.4702 * out1_1 + 0.2690 * out2_1 + 0.2608 * out2_2;
1612 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1613
1614 val = 0.4654 * out1_1 + 0.3044 * out2_1 + 0.2303 * out2_2;
1615 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1616
1617 val = 0.4380 * out1_1 + 0.4778 * out2_1 + 0.0842 * out2_2;
1618 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1619
1620 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1621 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1622
1623 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1624 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1625
1626 break;
1627 case 030002: /* ##..........#., <noise^2> = 0 */
1628 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1629 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1630
1631 val = -0.2661 * out1_2 + 0.9206 * out1_1 + 0.3455 * out2_1;
1632 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1633
1634 val = 0.0842 * out1_2 + 0.4778 * out1_1 + 0.4380 * out2_1;
1635 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1636
1637 val = 0.2303 * out1_2 + 0.3043 * out1_1 + 0.4654 * out2_1;
1638 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1639
1640 val = 0.2608 * out1_2 + 0.2690 * out1_1 + 0.4703 * out2_1;
1641 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1642
1643 val = 0.2639 * out1_2 + 0.2644 * out1_1 + 0.4717 * out2_1;
1644 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1645
1646 val = 0.2598 * out1_2 + 0.2598 * out1_1 + 0.4804 * out2_1;
1647 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1648
1649 val = 0.2367 * out1_2 + 0.2367 * out1_1 + 0.5265 * out2_1;
1650 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1651
1652 val = 0.1673 * out1_2 + 0.1673 * out1_1 + 0.6654 * out2_1;
1653 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1654
1655 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1656 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1657
1658 break;
1659 case 030003: /* ##..........##, <noise^2> = 0 */
1660 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0888 * out2_1 + 0.0888 * out2_2;
1661 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1662
1663 val = -0.2768 * out1_2 + 0.9098 * out1_1 + 0.1835 * out2_1 + 0.1835 * out2_2;
1664 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1665
1666 val = 0.0705 * out1_2 + 0.4642 * out1_1 + 0.2326 * out2_1 + 0.2326 * out2_2;
1667 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1668
1669 val = 0.2158 * out1_2 + 0.2899 * out1_1 + 0.2472 * out2_1 + 0.2471 * out2_2;
1670 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1671
1672 val = 0.2462 * out1_2 + 0.2544 * out1_1 + 0.2500 * out2_1 + 0.2495 * out2_2;
1673 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1674
1675 val = 0.2495 * out1_2 + 0.2500 * out1_1 + 0.2544 * out2_1 + 0.2462 * out2_2;
1676 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1677
1678 val = 0.2471 * out1_2 + 0.2472 * out1_1 + 0.2899 * out2_1 + 0.2158 * out2_2;
1679 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1680
1681 val = 0.2326 * out1_2 + 0.2326 * out1_1 + 0.4642 * out2_1 + 0.0705 * out2_2;
1682 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1683
1684 val = 0.1835 * out1_2 + 0.1835 * out1_1 + 0.9098 * out2_1 - 0.2768 * out2_2;
1685 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1686
1687 val = 0.0888 * out1_2 + 0.0888 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1688 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1689
1690 break;
1691 case 040003: /* #............##, <noise^2> = 0 */
1692 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1693 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1694
1695 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1696 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1697
1698 val = 0.5265 * out1_1 + 0.2367 * out2_1 + 0.2367 * out2_2;
1699 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1700
1701 val = 0.4804 * out1_1 + 0.2598 * out2_1 + 0.2598 * out2_2;
1702 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1703
1704 val = 0.4718 * out1_1 + 0.2641 * out2_1 + 0.2641 * out2_2;
1705 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1706
1707 val = 0.4708 * out1_1 + 0.2646 * out2_1 + 0.2646 * out2_2;
1708 out[badX0 + 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1709
1710 val = 0.4707 * out1_1 + 0.2649 * out2_1 + 0.2644 * out2_2;
1711 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1712
1713 val = 0.4702 * out1_1 + 0.2690 * out2_1 + 0.2608 * out2_2;
1714 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1715
1716 val = 0.4654 * out1_1 + 0.3044 * out2_1 + 0.2303 * out2_2;
1717 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1718
1719 val = 0.4380 * out1_1 + 0.4778 * out2_1 + 0.0842 * out2_2;
1720 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1721
1722 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1723 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1724
1725 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1726 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1727
1728 break;
1729 case 060002: /* ##...........#., <noise^2> = 0 */
1730 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1731 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1732
1733 val = -0.2661 * out1_2 + 0.9206 * out1_1 + 0.3455 * out2_1;
1734 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1735
1736 val = 0.0842 * out1_2 + 0.4778 * out1_1 + 0.4380 * out2_1;
1737 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1738
1739 val = 0.2303 * out1_2 + 0.3044 * out1_1 + 0.4654 * out2_1;
1740 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1741
1742 val = 0.2608 * out1_2 + 0.2690 * out1_1 + 0.4702 * out2_1;
1743 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1744
1745 val = 0.2644 * out1_2 + 0.2649 * out1_1 + 0.4708 * out2_1;
1746 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1747
1748 val = 0.2641 * out1_2 + 0.2641 * out1_1 + 0.4718 * out2_1;
1749 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1750
1751 val = 0.2598 * out1_2 + 0.2598 * out1_1 + 0.4804 * out2_1;
1752 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1753
1754 val = 0.2367 * out1_2 + 0.2367 * out1_1 + 0.5265 * out2_1;
1755 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1756
1757 val = 0.1673 * out1_2 + 0.1673 * out1_1 + 0.6654 * out2_1;
1758 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1759
1760 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1761 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1762
1763 break;
1764 default:
1765 // shFatal("Unsupported defect type: MIDDLE 0%o",defect[i].type);
1766 break; /* NOTREACHED */
1767 }
1768 break;
1769 case Defect::WIDE:
1770 case Defect::WIDE_NEAR_LEFT:
1771 case Defect::WIDE_NEAR_RIGHT:
1772 if (defectPos == Defect::WIDE_NEAR_LEFT) {
1773 assert(badX0 >= 1);
1774
1775 if (badX1 + 2 >= ncol) { /* left defect extends near
1776 right edge of data! */
1777 if (badX1 == ncol - 2) { /* one column remains */
1778 val = out[ncol - 1];
1779 } else {
1780 val = fallbackValue; /* there is no information */
1781 }
1782 for (int j = badX0; j <= badX1; j++) {
1783 out[j] = val;
1784 }
1785 break;
1786 }
1787 out1_2 = -1; /* NOTUSED */
1788 out2_2 = out[badX1 + 2];
1789 } else if (defectPos == Defect::WIDE) {
1790 assert(badX0 >= 2 && badX1 + 2 < ncol);
1791 out1_2 = out[badX0 - 2];
1792 out2_2 = out[badX1 + 2];
1793 } else if (defectPos == Defect::WIDE_NEAR_RIGHT) {
1794 assert(badX1 + 1 < ncol);
1795
1796 if (badX0 < 2) { /* right defect extends near
1797 left edge of data! */
1798 if (badX0 == 1) { /* one column remains */
1799 val = out[0];
1800 } else {
1801 val = fallbackValue; /* there is no information */
1802 }
1803 for (int j = badX0; j <= badX1; j++) {
1804 out[j] = val;
1805 }
1806 break;
1807 }
1808 out1_2 = out[badX0 - 2];
1809 out2_2 = -1; /* NOTUSED */
1810 } else {
1811 // shFatal("Unknown defect classification %d (%s:%d)",defectPos, __FILE__,__LINE__);
1812 out1_2 = out2_2 = -1; /* NOTUSED */
1813 }
1814
1815 out1_1 = out[badX0 - 1];
1816 out2_1 = out[badX1 + 1];
1817
1818 switch (defectType) {
1819 case 06: /* #?#., <noise^2> = 0 */
1820 val = 0.8894 * out1_1 + 0.1106 * out2_1;
1821 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1822
1823 val = 0.6839 * out1_1 + 0.3161 * out2_1;
1824 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1825
1826 val = 0.5527 * out1_1 + 0.4473 * out2_1;
1827 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1828
1829 val = 0.5092 * out1_1 + 0.4908 * out2_1;
1830 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1831
1832 val = 0.5010 * out1_1 + 0.4990 * out2_1;
1833 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1834
1835 val = 0.5001 * out1_1 + 0.4999 * out2_1;
1836 out[badX0 + 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1837
1838 val = 0.5000 * out1_1 + 0.5000 * out2_1;
1839
1840 for (int j = badX0 + 6; j < badX1 - 5; j++) {
1841 out[j] = val;
1842 }
1843
1844 val = 0.4999 * out1_1 + 0.5001 * out2_1;
1845 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1846
1847 val = 0.4990 * out1_1 + 0.5010 * out2_1;
1848 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1849
1850 val = 0.4908 * out1_1 + 0.5092 * out2_1;
1851 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1852
1853 val = 0.4473 * out1_1 + 0.5527 * out2_1;
1854 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1855
1856 val = 0.3161 * out1_1 + 0.6839 * out2_1;
1857 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1858
1859 val = 0.1106 * out1_1 + 0.8894 * out2_1;
1860 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1861
1862 break;
1863 case 07: /* #?##, <noise^2> = 0 */
1864 val = 0.8829 * out1_1 + 0.0585 * out2_1 + 0.0585 * out2_2;
1865 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1866
1867 val = 0.6654 * out1_1 + 0.1673 * out2_1 + 0.1673 * out2_2;
1868 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1869
1870 val = 0.5265 * out1_1 + 0.2367 * out2_1 + 0.2367 * out2_2;
1871 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1872
1873 val = 0.4804 * out1_1 + 0.2598 * out2_1 + 0.2598 * out2_2;
1874 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1875
1876 val = 0.4718 * out1_1 + 0.2641 * out2_1 + 0.2641 * out2_2;
1877 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1878
1879 val = 0.4708 * out1_1 + 0.2646 * out2_1 + 0.2646 * out2_2;
1880 out[badX0 + 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1881
1882 val = 0.4707 * out[badX0 - 1] + 0.2646 * out[badX1 + 1] + 0.2646 * out[badX1 + 2];
1883
1884 for (int j = badX0 + 6; j < badX1 - 5; j++) {
1885 out[j] = val;
1886 }
1887
1888 val = 0.4707 * out1_1 + 0.2649 * out2_1 + 0.2644 * out2_2;
1889 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1890
1891 val = 0.4702 * out1_1 + 0.2690 * out2_1 + 0.2608 * out2_2;
1892 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1893
1894 val = 0.4654 * out1_1 + 0.3044 * out2_1 + 0.2303 * out2_2;
1895 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1896
1897 val = 0.4380 * out1_1 + 0.4778 * out2_1 + 0.0842 * out2_2;
1898 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1899
1900 val = 0.3455 * out1_1 + 0.9206 * out2_1 - 0.2661 * out2_2;
1901 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1902
1903 val = 0.1673 * out1_1 + 1.3452 * out2_1 - 0.5125 * out2_2;
1904 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1905
1906 break;
1907 case 016: /* ##?#., <noise^2> = 0 */
1908 val = -0.5125 * out1_2 + 1.3452 * out1_1 + 0.1673 * out2_1;
1909 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1910
1911 val = -0.2661 * out1_2 + 0.9206 * out1_1 + 0.3455 * out2_1;
1912 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1913
1914 val = 0.0842 * out1_2 + 0.4778 * out1_1 + 0.4380 * out2_1;
1915 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1916
1917 val = 0.2303 * out1_2 + 0.3044 * out1_1 + 0.4654 * out2_1;
1918 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1919
1920 val = 0.2608 * out1_2 + 0.2690 * out1_1 + 0.4702 * out2_1;
1921 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1922
1923 val = 0.2644 * out1_2 + 0.2649 * out1_1 + 0.4707 * out2_1;
1924 out[badX0 + 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1925
1926 val = 0.2646 * out1_2 + 0.2646 * out1_1 + 0.4707 * out2_1;
1927 val = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1928
1929 for (int j = badX0 + 6; j < badX1 - 5; j++) {
1930 out[j] = val;
1931 }
1932
1933 val = 0.2646 * out1_2 + 0.2646 * out1_1 + 0.4708 * out2_1;
1934 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1935
1936 val = 0.2641 * out1_2 + 0.2641 * out1_1 + 0.4718 * out2_1;
1937 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1938
1939 val = 0.2598 * out1_2 + 0.2598 * out1_1 + 0.4804 * out2_1;
1940 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1941
1942 val = 0.2367 * out1_2 + 0.2367 * out1_1 + 0.5265 * out2_1;
1943 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1944
1945 val = 0.1673 * out1_2 + 0.1673 * out1_1 + 0.6654 * out2_1;
1946 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1947
1948 val = 0.0585 * out1_2 + 0.0585 * out1_1 + 0.8829 * out2_1;
1949 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1950
1951 break;
1952 case 017: /* ##?##, S/N = infty */
1953 val = -0.5177 * out1_2 + 1.3400 * out1_1 + 0.0888 * out2_1 + 0.0888 * out2_2;
1954 out[badX0] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1955
1956 val = -0.2768 * out1_2 + 0.9098 * out1_1 + 0.1835 * out2_1 + 0.1835 * out2_2;
1957 out[badX0 + 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1958
1959 val = 0.0705 * out1_2 + 0.4642 * out1_1 + 0.2326 * out2_1 + 0.2326 * out2_2;
1960 out[badX0 + 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1961
1962 val = 0.2158 * out1_2 + 0.2899 * out1_1 + 0.2472 * out2_1 + 0.2472 * out2_2;
1963 out[badX0 + 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1964
1965 val = 0.2462 * out1_2 + 0.2544 * out1_1 + 0.2497 * out2_1 + 0.2497 * out2_2;
1966 out[badX0 + 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1967
1968 val = 0.2497 * out1_2 + 0.2503 * out1_1 + 0.2500 * out2_1 + 0.2500 * out2_2;
1969 out[badX0 + 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1970
1971 val = 0.2500 * out1_2 + 0.2500 * out1_1 + 0.2500 * out2_1 + 0.2500 * out2_2;
1972 val = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1973
1974 for (int j = badX0 + 6; j < badX1 - 5; j++) {
1975 out[j] = val;
1976 }
1977
1978 val = 0.2500 * out1_2 + 0.2500 * out1_1 + 0.2503 * out2_1 + 0.2497 * out2_2;
1979 out[badX1 - 5] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1980
1981 val = 0.2497 * out1_2 + 0.2497 * out1_1 + 0.2544 * out2_1 + 0.2462 * out2_2;
1982 out[badX1 - 4] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1983
1984 val = 0.2472 * out1_2 + 0.2472 * out1_1 + 0.2899 * out2_1 + 0.2158 * out2_2;
1985 out[badX1 - 3] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1986
1987 val = 0.2326 * out1_2 + 0.2326 * out1_1 + 0.4642 * out2_1 + 0.0705 * out2_2;
1988 out[badX1 - 2] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1989
1990 val = 0.1835 * out1_2 + 0.1835 * out1_1 + 0.9098 * out2_1 - 0.2768 * out2_2;
1991 out[badX1 - 1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1992
1993 val = 0.0888 * out1_2 + 0.0888 * out1_1 + 1.3400 * out2_1 - 0.5177 * out2_2;
1994 out[badX1] = (val < min) ? 0.5 * (out1_1 + out2_1) : val;
1995
1996 break;
1997 default:
1998 // shFatal("Unsupported defect type: WIDE 0%o",defect[i].type);
1999 break; /* NOTREACHED */
2000 }
2001 break;
2002 }
2003 }
2004}
2005
2006template <typename MaskT>
2007static void do_defects(std::vector<Defect::Ptr> const &badList, // list of bad things
2008 int const y, // Row that we should fix
2009 MaskT &mask, // mask to set
2010 typename MaskT::Pixel const interpBit, // bit to set for bad pixels
2011 bool useFallbackValueAtEdge, // use fallbackValue at edge of chip?
2012 int nUseInterp // no. of pixels to interpolate towards edge
2013) {
2014 typename MaskT::x_iterator mask_row = mask.row_begin(y); // pointer to this row of mask
2015
2016 for (DefectCIter ptr = badList.begin(), end = badList.end(); ptr != end; ++ptr) {
2017 Defect::Ptr const defect = *ptr;
2018
2019 if (y < defect->getY0() || y > defect->getY1()) {
2020 continue;
2021 }
2022
2023 int const badX0 = defect->getX0();
2024 int const badX1 = defect->getX1();
2025
2026 for (int c = badX0; c <= badX1; ++c) {
2027 mask_row[c] |= interpBit;
2028 }
2029 }
2030}
2031
2032/************************************************************************************************************/
2033
2034namespace {
2035template <typename T>
2036struct Sort_ByX0 {
2037 bool operator()(std::shared_ptr<T> const a, std::shared_ptr<T> const b) const {
2038 return a->getX0() < b->getX0();
2039 }
2040};
2041} // namespace
2042
2046template <typename MaskedImageT>
2047void interpolateOverDefects(MaskedImageT &mimage,
2048 afw::detection::Psf const &,
2049 std::vector<Defect::Ptr> &_badList,
2050 double fallbackValue,
2051 bool useFallbackValueAtEdge
2052) {
2053 /*
2054 * Allow for image's origin
2055 */
2056 int const width = mimage.getWidth();
2057 int const height = mimage.getHeight();
2058
2059 std::vector<Defect::Ptr> badList;
2060 badList.reserve(_badList.size());
2061 for (std::vector<Defect::Ptr>::iterator ptr = _badList.begin(), end = _badList.end(); ptr != end; ++ptr) {
2062 geom::BoxI bbox = (*ptr)->getBBox();
2063 bbox.shift(geom::ExtentI(-mimage.getX0(), -mimage.getY0())); // allow for image's origin
2064 geom::PointI min = bbox.getMin(), max = bbox.getMax();
2065 if (min.getX() >= width) {
2066 continue;
2067 } else if (min.getX() < 0) {
2068 if (max.getX() < 0) {
2069 continue;
2070 } else {
2071 min.setX(0);
2072 }
2073 }
2074
2075 if (max.getX() < 0) {
2076 continue;
2077 } else if (max.getX() >= width) {
2078 max.setX(width - 1);
2079 }
2080
2081 bbox = geom::BoxI(min, max);
2082 Defect::Ptr ndefect(new Defect(bbox));
2083 ndefect->classify((*ptr)->getPos(), (*ptr)->getType());
2084 badList.push_back(ndefect);
2085 }
2086
2087 sort(badList.begin(), badList.end(), Sort_ByX0<Defect>());
2088 /*
2089 * Go through the frame looking at each pixel (except the edge ones which we ignore)
2090 */
2091 typename MaskedImageT::Mask::Pixel const interpBit =
2092 mimage.getMask()->getPlaneBitMask("INTRP"); // interp'd pixels
2093
2094 constexpr int nUseInterp = 6; // no. of pixels to interpolate towards edge
2095 static_assert(nUseInterp < Defect::WIDE_DEFECT,
2096 "make sure that we can handle these defects using"
2097 "the full interpolation not edge code");
2098
2099 for (int y = 0; y != height; y++) {
2100 std::vector<Defect::Ptr> badList1D = classify_defects(badList, y, width);
2101
2102 do_defects(badList1D, y, *mimage.getImage(),
2103 -std::numeric_limits<typename MaskedImageT::Image::Pixel>::max(), fallbackValue,
2104 useFallbackValueAtEdge, nUseInterp);
2105
2106 do_defects(badList1D, y, *mimage.getMask(), interpBit, useFallbackValueAtEdge, nUseInterp);
2107
2108 do_defects(badList1D, y, *mimage.getVariance(),
2109 -std::numeric_limits<typename MaskedImageT::Image::Pixel>::max(), fallbackValue,
2110 useFallbackValueAtEdge, nUseInterp);
2111 }
2112}
2113
2114/*****************************************************************************/
2124template <typename MaskedImageT>
2125std::pair<bool, typename MaskedImageT::Image::Pixel> interp::singlePixel(
2126 int,
2127 int,
2128 MaskedImageT const &,
2129 bool,
2130 double
2131) {
2132#if defined(SDSS)
2133 BADCOLUMN defect; /* describe a bad column */
2134 PIX *data; /* temp array to interpolate in */
2135 int i;
2136 int i0, i1; /* data corresponds to range of
2137 {row,col} == [i0,i1] */
2138 int ndata; /* dimension of data */
2139 static int ndatamax = 40; /* largest allowable defect. XXX */
2140 int nrow, ncol; /* == reg->n{row,col} */
2141 PIX *val; /* pointer to pixel (rowc, colc) */
2142 int z1, z2; /* range of bad {row,columns} */
2143
2144 shAssert(badmask != NULL && badmask->type == shTypeGetFromName("OBJMASK"));
2145 shAssert(reg != NULL && reg->type == TYPE_PIX);
2146 nrow = reg->nrow;
2147 ncol = reg->ncol;
2148
2149 if (horizontal) {
2150 for (z1 = colc - 1; z1 >= 0; z1--) {
2151 if (!phPixIntersectMask(badmask, z1, rowc)) {
2152 break;
2153 }
2154 }
2155 z1++;
2156
2157 for (z2 = colc + 1; z2 < ncol; z2++) {
2158 if (!phPixIntersectMask(badmask, z2, rowc)) {
2159 break;
2160 }
2161 }
2162 z2--;
2163
2164 i0 = (z1 > 2) ? z1 - 2 : 0; /* origin of available required data */
2165 i1 = (z2 < ncol - 2) ? z2 + 2 : ncol - 1; /* end of " " " " */
2166
2167 if (i0 < 2 || i1 >= ncol - 2) { /* interpolation will fail */
2168 return (-1); /* failure */
2169 }
2170
2171 ndata = (i1 - i0 + 1);
2172 if (ndata > ndatamax) {
2173 return (-1); /* failure */
2174 }
2175
2176 data = alloca(ndata * sizeof(PIX));
2177 for (i = i0; i <= i1; i++) {
2178 data[i - i0] = reg->ROWS[rowc][i];
2179 }
2180 val = &data[colc - i0];
2181 } else {
2182 for (z1 = rowc - 1; z1 >= 0; z1--) {
2183 if (!phPixIntersectMask(badmask, colc, z1)) {
2184 break;
2185 }
2186 }
2187 z1++;
2188
2189 for (z2 = rowc + 1; z2 < nrow; z2++) {
2190 if (!phPixIntersectMask(badmask, colc, z2)) {
2191 break;
2192 }
2193 }
2194 z2--;
2195
2196 i0 = (z1 > 2) ? z1 - 2 : 0; /* origin of available required data */
2197 i1 = (z2 < nrow - 2) ? z2 + 2 : nrow - 1; /* end of " " " " */
2198
2199 if (i0 < 2 || i1 >= ncol - 2) { /* interpolation will fail */
2200 return (-1); /* failure */
2201 }
2202
2203 ndata = (i1 - i0 + 1);
2204 if (ndata > ndatamax) {
2205 return (-1); /* failure */
2206 }
2207
2208 data = alloca(ndata * sizeof(PIX));
2209 for (i = i0; i <= i1; i++) {
2210 data[i - i0] = reg->ROWS[i][colc];
2211 }
2212 val = &data[rowc - i0];
2213 }
2214
2215 defect.x1 = z1 - i0;
2216 defect.x2 = z2 - i0;
2217 classify_defects(&defect, 1, ndata);
2218 do_defect(&defect, 1, data, ndata, minval);
2219
2220 return (*val);
2221#endif
2222
2223 return std::make_pair(false, std::numeric_limits<typename MaskedImageT::Image::Pixel>::min());
2224}
2225
2226/************************************************************************************************************/
2227//
2228// Explicit instantiations
2229//
2230// \cond
2231
2232typedef float ImagePixel;
2233
2234template void interpolateOverDefects(afw::image::MaskedImage<ImagePixel, afw::image::MaskPixel> &image,
2235 afw::detection::Psf const &, std::vector<Defect::Ptr> &badList, double,
2236 bool);
2237template std::pair<bool, ImagePixel> interp::singlePixel(
2238 int x, int y, afw::image::MaskedImage<ImagePixel, afw::image::MaskPixel> const &image,
2239 bool horizontal, double minval);
2240//
2241// Why do we need double images?
2242//
2243#if 1
2244template void interpolateOverDefects(afw::image::MaskedImage<double, afw::image::MaskPixel> &image,
2245 afw::detection::Psf const &, std::vector<Defect::Ptr> &badList, double,
2246 bool);
2247
2248template std::pair<bool, double> interp::singlePixel(
2249 int x, int y, afw::image::MaskedImage<double, afw::image::MaskPixel> const &image, bool horizontal,
2250 double minval);
2251
2252#endif
2253// \endcond
2254
2255} // namespace algorithms
2256} // namespace meas
2257} // namespace lsst
bbox
AmpInfoBoxKey bbox
Definition: Amplifier.cc:117
data
char * data
Definition: BaseRecord.cc:61
min
int min
Definition: BoundedField.cc:103
end
int end
Definition: BoundedField.cc:105
max
int max
Definition: BoundedField.cc:104
x
double x
Definition: ChebyshevBoundedField.cc:276
mask
afw::table::Key< afw::table::Array< MaskPixelT > > mask
Definition: HeavyFootprint.cc:213
ptr
uint64_t * ptr
Definition: RangeSet.cc:88
y
int y
Definition: SpanSet.cc:48
b
table::Key< int > b
Definition: TransmissionCurve.cc:466
a
table::Key< int > a
Definition: TransmissionCurve.cc:465
std::vector::begin
T begin(T... args)
lsst::afw::detection::Psf
A polymorphic base class for representing an image's Point Spread Function.
Definition: Psf.h:76
lsst::afw::image::MaskedImage
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
lsst::geom::Box2I
An integer coordinate rectangle.
Definition: Box.h:55
lsst::meas::algorithms::Defect
Encapsulate information about a bad portion of a detector.
Definition: Interp.h:72
lsst::meas::algorithms::Defect::Ptr
std::shared_ptr< Defect > Ptr
shared pointer to Defect
Definition: Interp.h:74
lsst::meas::algorithms::Defect::WIDE
@ WIDE
defect is in middle, and wide
Definition: Interp.h:82
lsst::meas::algorithms::Defect::NEAR_RIGHT
@ NEAR_RIGHT
defect is near right boundary
Definition: Interp.h:84
lsst::meas::algorithms::Defect::NEAR_LEFT
@ NEAR_LEFT
defect is near left boundary
Definition: Interp.h:78
lsst::meas::algorithms::Defect::WIDE_LEFT
@ WIDE_LEFT
defect is wide at left boundary
Definition: Interp.h:79
lsst::meas::algorithms::Defect::WIDE_RIGHT
@ WIDE_RIGHT
defect is wide at right boundary
Definition: Interp.h:85
lsst::meas::algorithms::Defect::MIDDLE
@ MIDDLE
defect is in middle of frame
Definition: Interp.h:80
lsst::meas::algorithms::Defect::WIDE_NEAR_LEFT
@ WIDE_NEAR_LEFT
defect is near left, and wide
Definition: Interp.h:81
lsst::meas::algorithms::Defect::RIGHT
@ RIGHT
defect is at right boundary
Definition: Interp.h:86
lsst::meas::algorithms::Defect::LEFT
@ LEFT
defect is at left boundary
Definition: Interp.h:77
lsst::meas::algorithms::Defect::WIDE_NEAR_RIGHT
@ WIDE_NEAR_RIGHT
defect is near right, and wide
Definition: Interp.h:83
std::vector::end
T end(T... args)
std::make_pair
T make_pair(T... args)
lsst::geom::BoxI
Box2I BoxI
Definition: Box.h:780
lsst::meas::algorithms::interp::singlePixel
std::pair< bool, typename MaskedImageT::Image::Pixel > singlePixel(int x, int y, MaskedImageT const &image, bool horizontal, double minval)
Return a boolean status (true: interpolation is OK) and the interpolated value for a pixel,...
Definition: Interp.cc:2125
lsst::meas::algorithms::DefectCIter
std::vector< Defect::Ptr >::const_iterator DefectCIter
Definition: Interp.cc:49
lsst::meas::algorithms::interpolateOverDefects
void interpolateOverDefects(MaskedImageT &image, afw::detection::Psf const &psf, std::vector< Defect::Ptr > &badList, double fallbackValue=0.0, bool useFallbackValueAtEdge=false)
Process a set of known bad pixels in an image.
Definition: Interp.cc:2047
exceptions.h
std::vector::push_back
T push_back(T... args)
std::vector::reserve
T reserve(T... args)
std::vector::size
T size(T... args)
val
ImageT val
Definition: CR.cc:146
Generated on Thu Nov 24 2022 08:01:24 for LSST Applications by doxygen 1.9.5