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