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