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