[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

noise_normalization.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2006 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 
37 #ifndef VIGRA_NOISE_NORMALIZATION_HXX
38 #define VIGRA_NOISE_NORMALIZATION_HXX
39 
40 #include "utilities.hxx"
41 #include "tinyvector.hxx"
42 #include "stdimage.hxx"
43 #include "transformimage.hxx"
44 #include "combineimages.hxx"
45 #include "localminmax.hxx"
46 #include "functorexpression.hxx"
47 #include "numerictraits.hxx"
48 #include "separableconvolution.hxx"
49 #include "linear_solve.hxx"
50 #include "array_vector.hxx"
51 #include "static_assert.hxx"
52 #include "multi_shape.hxx"
53 
54 #include <algorithm>
55 
56 namespace vigra {
57 
58 /** \addtogroup NoiseNormalization Noise Normalization
59  Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
60 */
61 //@{
62 
63 /********************************************************/
64 /* */
65 /* NoiseNormalizationOptions */
66 /* */
67 /********************************************************/
68 
69 /** \brief Pass options to one of the noise normalization functions.
70 
71  <tt>NoiseNormalizationOptions</tt> is an argument object that holds various optional
72  parameters used by the noise normalization functions. If a parameter is not explicitly
73  set, a suitable default will be used.
74 
75  <b> Usage:</b>
76 
77  <b>\#include</b> <vigra/noise_normalization.hxx><br>
78  Namespace: vigra
79 
80  \code
81  MultiArray<2, float> src(w,h);
82  std::vector<TinyVector<double, 2> > result;
83 
84  ...
85  noiseVarianceEstimation(src, result,
86  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
87  \endcode
88 */
90 {
91  public:
92 
93  /** Initialize all options with default values.
94  */
96  : window_radius(6),
97  cluster_count(10),
98  noise_estimation_quantile(1.5),
99  averaging_quantile(0.8),
100  noise_variance_initial_guess(10.0),
101  use_gradient(true)
102  {}
103 
104  /** Select the noise estimation algorithm.
105 
106  If \a r is <tt>true</tt>, use the gradient-based noise estimator according to F&ouml;rstner (default).
107  Otherwise, use an algorithm that uses the intensity values directly.
108  */
110  {
111  use_gradient = r;
112  return *this;
113  }
114 
115  /** Set the window radius for a single noise estimate.
116  Every window of the given size gives raise to one intensity/variance pair.<br>
117  Default: 6 pixels
118  */
120  {
121  vigra_precondition(r > 0,
122  "NoiseNormalizationOptions: window radius must be > 0.");
123  window_radius = r;
124  return *this;
125  }
126 
127  /** Set the number of clusters for non-parametric noise normalization.
128  The intensity/variance pairs found are grouped into clusters before the noise
129  normalization transform is computed.<br>
130  Default: 10 clusters
131  */
133  {
134  vigra_precondition(c > 0,
135  "NoiseNormalizationOptions: cluster count must be > 0.");
136  cluster_count = c;
137  return *this;
138  }
139 
140  /** Set the quantile for cluster averaging.
141  After clustering, the cluster center (i.e. average noise variance as a function of the average
142  intensity in the cluster) is computed using only the cluster members whose estimated variance
143  is below \a quantile times the maximum variance in the cluster.<br>
144  Default: 0.8<br>
145  Precondition: 0 < \a quantile <= 1.0
146  */
148  {
149  vigra_precondition(quantile > 0.0 && quantile <= 1.0,
150  "NoiseNormalizationOptions: averaging quantile must be between 0 and 1.");
151  averaging_quantile = quantile;
152  return *this;
153  }
154 
155  /** Set the operating range of the robust noise estimator.
156  Intensity changes that are larger than \a quantile times the current estimate of the noise variance
157  are ignored by the robust noise estimator.<br>
158  Default: 1.5<br>
159  Precondition: 0 < \a quantile
160  */
162  {
163  vigra_precondition(quantile > 0.0,
164  "NoiseNormalizationOptions: noise estimation quantile must be > 0.");
165  noise_estimation_quantile = quantile;
166  return *this;
167  }
168 
169  /** Set the initial estimate of the noise variance.
170  Robust noise variance estimation is an iterative procedure starting at the given value.<br>
171  Default: 10.0<br>
172  Precondition: 0 < \a guess
173  */
175  {
176  vigra_precondition(guess > 0.0,
177  "NoiseNormalizationOptions: noise variance initial guess must be > 0.");
178  noise_variance_initial_guess = guess;
179  return *this;
180  }
181 
182  unsigned int window_radius, cluster_count;
183  double noise_estimation_quantile, averaging_quantile, noise_variance_initial_guess;
184  bool use_gradient;
185 };
186 
187 //@}
188 
189 template <class ArgumentType, class ResultType>
190 class NonparametricNoiseNormalizationFunctor
191 {
192  struct Segment
193  {
194  double lower, a, b, shift;
195  };
196 
197  ArrayVector<Segment> segments_;
198 
199  template <class T>
200  double exec(unsigned int k, T t) const
201  {
202  if(segments_[k].a == 0.0)
203  {
204  return t / VIGRA_CSTD::sqrt(segments_[k].b);
205  }
206  else
207  {
208  return 2.0 / segments_[k].a * VIGRA_CSTD::sqrt(std::max(0.0, segments_[k].a * t + segments_[k].b));
209  }
210  }
211 
212  public:
213  typedef ArgumentType argument_type;
214  typedef ResultType result_type;
215 
216  template <class Vector>
217  NonparametricNoiseNormalizationFunctor(Vector const & clusters)
218  : segments_(clusters.size()-1)
219  {
220  for(unsigned int k = 0; k<segments_.size(); ++k)
221  {
222  segments_[k].lower = clusters[k][0];
223  segments_[k].a = (clusters[k+1][1] - clusters[k][1]) / (clusters[k+1][0] - clusters[k][0]);
224  segments_[k].b = clusters[k][1] - segments_[k].a * clusters[k][0];
225  // FIXME: set a to zero if it is very small
226  // - determine what 'very small' means
227  // - shouldn't the two formulas (for a == 0, a != 0) be equal in the limit a -> 0 ?
228 
229  if(k == 0)
230  {
231  segments_[k].shift = segments_[k].lower - exec(k, segments_[k].lower);
232  }
233  else
234  {
235  segments_[k].shift = exec(k-1, segments_[k].lower) - exec(k, segments_[k].lower) + segments_[k-1].shift;
236  }
237  }
238  }
239 
240  result_type operator()(argument_type t) const
241  {
242  // find the segment
243  unsigned int k = 0;
244  for(; k < segments_.size(); ++k)
245  if(t < segments_[k].lower)
246  break;
247  if(k > 0)
248  --k;
249  return detail::RequiresExplicitCast<ResultType>::cast(exec(k, t) + segments_[k].shift);
250  }
251 };
252 
253 template <class ArgumentType, class ResultType>
254 class QuadraticNoiseNormalizationFunctor
255 {
256  double a, b, c, d, f, o;
257 
258  void init(double ia, double ib, double ic, double xmin)
259  {
260  a = ia;
261  b = ib;
262  c = ic;
263  d = VIGRA_CSTD::sqrt(VIGRA_CSTD::fabs(c));
264  if(c > 0.0)
265  {
266  o = VIGRA_CSTD::log(VIGRA_CSTD::fabs((2.0*c*xmin + b)/d + 2*VIGRA_CSTD::sqrt(c*sq(xmin) +b*xmin + a)))/d;
267  f = 0.0;
268  }
269  else
270  {
271  f = VIGRA_CSTD::sqrt(b*b - 4.0*a*c);
272  o = -VIGRA_CSTD::asin((2.0*c*xmin+b)/f)/d;
273  }
274  }
275 
276  public:
277  typedef ArgumentType argument_type;
278  typedef ResultType result_type;
279 
280  template <class Vector>
281  QuadraticNoiseNormalizationFunctor(Vector const & clusters)
282  {
283  double xmin = NumericTraits<double>::max();
284  Matrix<double> m(3,3), r(3, 1), l(3, 1);
285  for(unsigned int k = 0; k<clusters.size(); ++k)
286  {
287  l(0,0) = 1.0;
288  l(1,0) = clusters[k][0];
289  l(2,0) = sq(clusters[k][0]);
290  m += outer(l);
291  r += clusters[k][1]*l;
292  if(clusters[k][0] < xmin)
293  xmin = clusters[k][0];
294  }
295 
296  linearSolve(m, r, l);
297  init(l(0,0), l(1,0), l(2,0), xmin);
298  }
299 
300  result_type operator()(argument_type t) const
301  {
302  double r;
303  if(c > 0.0)
304  r = VIGRA_CSTD::log(VIGRA_CSTD::fabs((2.0*c*t + b)/d + 2.0*VIGRA_CSTD::sqrt(c*t*t +b*t + a)))/d-o;
305  else
306  r = -VIGRA_CSTD::asin((2.0*c*t+b)/f)/d-o;
307  return detail::RequiresExplicitCast<ResultType>::cast(r);
308  }
309 };
310 
311 template <class ArgumentType, class ResultType>
312 class LinearNoiseNormalizationFunctor
313 {
314  double a, b, o;
315 
316  void init(double ia, double ib, double xmin)
317  {
318  a = ia;
319  b = ib;
320  if(b != 0.0)
321  {
322  o = xmin - 2.0 / b * VIGRA_CSTD::sqrt(a + b * xmin);
323  }
324  else
325  {
326  o = xmin - xmin / VIGRA_CSTD::sqrt(a);
327  }
328  }
329 
330  public:
331  typedef ArgumentType argument_type;
332  typedef ResultType result_type;
333 
334  template <class Vector>
335  LinearNoiseNormalizationFunctor(Vector const & clusters)
336  {
337  double xmin = NumericTraits<double>::max();
338  Matrix<double> m(2,2), r(2, 1), l(2, 1);
339  for(unsigned int k = 0; k<clusters.size(); ++k)
340  {
341  l(0,0) = 1.0;
342  l(1,0) = clusters[k][0];
343  m += outer(l);
344  r += clusters[k][1]*l;
345  if(clusters[k][0] < xmin)
346  xmin = clusters[k][0];
347  }
348 
349  linearSolve(m, r, l);
350  init(l(0,0), l(1,0), xmin);
351  }
352 
353  result_type operator()(argument_type t) const
354  {
355  double r;
356  if(b != 0.0)
357  r = 2.0 / b * VIGRA_CSTD::sqrt(a + b*t) + o;
358  else
359  r = t / VIGRA_CSTD::sqrt(a) + o;
360  return detail::RequiresExplicitCast<ResultType>::cast(r);
361  }
362 };
363 
364 #define VIGRA_NoiseNormalizationFunctor(name, type, size) \
365 template <class ResultType> \
366 class name<type, ResultType> \
367 { \
368  ResultType lut_[size]; \
369  \
370  public: \
371  typedef type argument_type; \
372  typedef ResultType result_type; \
373  \
374  template <class Vector> \
375  name(Vector const & clusters) \
376  { \
377  name<double, ResultType> f(clusters); \
378  \
379  for(unsigned int k = 0; k < size; ++k) \
380  { \
381  lut_[k] = f(k); \
382  } \
383  } \
384  \
385  result_type operator()(argument_type t) const \
386  { \
387  return lut_[t]; \
388  } \
389 };
390 
391 VIGRA_NoiseNormalizationFunctor(NonparametricNoiseNormalizationFunctor, UInt8, 256)
392 VIGRA_NoiseNormalizationFunctor(NonparametricNoiseNormalizationFunctor, UInt16, 65536)
393 VIGRA_NoiseNormalizationFunctor(QuadraticNoiseNormalizationFunctor, UInt8, 256)
394 VIGRA_NoiseNormalizationFunctor(QuadraticNoiseNormalizationFunctor, UInt16, 65536)
395 VIGRA_NoiseNormalizationFunctor(LinearNoiseNormalizationFunctor, UInt8, 256)
396 VIGRA_NoiseNormalizationFunctor(LinearNoiseNormalizationFunctor, UInt16, 65536)
397 
398 #undef VIGRA_NoiseNormalizationFunctor
399 
400 namespace detail {
401 
402 template <class SrcIterator, class SrcAcessor,
403  class GradIterator>
404 bool
405 iterativeNoiseEstimationChi2(SrcIterator s, SrcAcessor src, GradIterator g,
406  double & mean, double & variance,
407  double robustnessThreshold, int windowRadius)
408 {
409  double l2 = sq(robustnessThreshold);
410  double countThreshold = 1.0 - VIGRA_CSTD::exp(-l2);
411  double f = (1.0 - VIGRA_CSTD::exp(-l2)) / (1.0 - (1.0 + l2)*VIGRA_CSTD::exp(-l2));
412 
413  Diff2D ul(-windowRadius, -windowRadius);
414  int r2 = sq(windowRadius);
415 
416  for(int iter=0; iter<100 ; ++iter) // maximum iteration 100 only for terminating
417  // if something is wrong
418  {
419  double sum=0.0;
420  double gsum=0.0;
421  unsigned int count = 0;
422  unsigned int tcount = 0;
423 
424  SrcIterator siy = s + ul;
425  GradIterator giy = g + ul;
426  for(int y=-windowRadius; y <= windowRadius; y++, ++siy.y, ++giy.y)
427  {
428  typename SrcIterator::row_iterator six = siy.rowIterator();
429  typename GradIterator::row_iterator gix = giy.rowIterator();
430  for(int x=-windowRadius; x <= windowRadius; x++, ++six, ++gix)
431  {
432  if (sq(x) + sq(y) > r2)
433  continue;
434 
435  ++tcount;
436  if (*gix < l2*variance)
437  {
438  sum += src(six);
439  gsum += *gix;
440  ++count;
441  }
442  }
443  }
444  if (count==0) // not homogeneous enough
445  return false;
446 
447  double oldvariance = variance;
448  variance= f * gsum / count;
449  mean = sum / count;
450 
451  if ( closeAtTolerance(oldvariance - variance, 0.0, 1e-10))
452  return (count >= tcount * countThreshold / 2.0); // sufficiently many valid points
453  }
454  return false; // no convergence
455 }
456 
457 template <class SrcIterator, class SrcAcessor,
458  class GradIterator>
459 bool
460 iterativeNoiseEstimationGauss(SrcIterator s, SrcAcessor src, GradIterator,
461  double & mean, double & variance,
462  double robustnessThreshold, int windowRadius)
463 {
464  double l2 = sq(robustnessThreshold);
465  double countThreshold = erf(VIGRA_CSTD::sqrt(0.5 * l2));
466  double f = countThreshold / (countThreshold - VIGRA_CSTD::sqrt(2.0/M_PI*l2)*VIGRA_CSTD::exp(-l2/2.0));
467 
468  mean = src(s);
469 
470  Diff2D ul(-windowRadius, -windowRadius);
471  int r2 = sq(windowRadius);
472 
473  for(int iter=0; iter<100 ; ++iter) // maximum iteration 100 only for terminating
474  // if something is wrong
475  {
476  double sum = 0.0;
477  double sum2 = 0.0;
478  unsigned int count = 0;
479  unsigned int tcount = 0;
480 
481  SrcIterator siy = s + ul;
482  for(int y=-windowRadius; y <= windowRadius; y++, ++siy.y)
483  {
484  typename SrcIterator::row_iterator six = siy.rowIterator();
485  for(int x=-windowRadius; x <= windowRadius; x++, ++six)
486  {
487  if (sq(x) + sq(y) > r2)
488  continue;
489 
490  ++tcount;
491  if (sq(src(six) - mean) < l2*variance)
492  {
493  sum += src(six);
494  sum2 += sq(src(six));
495  ++count;
496  }
497  }
498  }
499  if (count==0) // not homogeneous enough
500  return false;
501 
502  double oldmean = mean;
503  double oldvariance = variance;
504  mean = sum / count;
505  variance= f * (sum2 / count - sq(mean));
506 
507  if ( closeAtTolerance(oldmean - mean, 0.0, 1e-10) &&
508  closeAtTolerance(oldvariance - variance, 0.0, 1e-10))
509  return (count >= tcount * countThreshold / 2.0); // sufficiently many valid points
510  }
511  return false; // no convergence
512 }
513 
514 
515 template <class SrcIterator, class SrcAccessor,
516  class DestIterator, class DestAccessor>
517 void
518 symmetricDifferenceSquaredMagnitude(
519  SrcIterator sul, SrcIterator slr, SrcAccessor src,
520  DestIterator dul, DestAccessor dest)
521 {
522  using namespace functor;
523  int w = slr.x - sul.x;
524  int h = slr.y - sul.y;
525 
526  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
527  typedef BasicImage<TmpType> TmpImage;
528 
529  Kernel1D<double> mask;
530  mask.initSymmetricGradient();
531  mask.setBorderTreatment(BORDER_TREATMENT_REFLECT);
532 
533  TmpImage dx(w, h), dy(w, h);
534  separableConvolveX(srcIterRange(sul, slr, src), destImage(dx), kernel1d(mask));
535  separableConvolveY(srcIterRange(sul, slr, src), destImage(dy), kernel1d(mask));
536  combineTwoImages(srcImageRange(dx), srcImage(dy), destIter(dul, dest), Arg1()*Arg1() + Arg2()*Arg2());
537 }
538 
539 template <class SrcIterator, class SrcAccessor,
540  class DestIterator, class DestAccessor>
541 void
542 findHomogeneousRegionsFoerstner(
543  SrcIterator sul, SrcIterator slr, SrcAccessor src,
544  DestIterator dul, DestAccessor dest,
545  unsigned int windowRadius = 6, double homogeneityThreshold = 40.0)
546 {
547  using namespace vigra::functor;
548  int w = slr.x - sul.x;
549  int h = slr.y - sul.y;
550 
551  BImage btmp(w, h);
552  transformImage(srcIterRange(sul, slr, src), destImage(btmp),
553  ifThenElse(Arg1() <= Param(homogeneityThreshold), Param(1), Param(0)));
554 
555  // Erosion
556  discErosion(srcImageRange(btmp), destIter(dul, dest), windowRadius);
557 }
558 
559 template <class SrcIterator, class SrcAccessor,
560  class DestIterator, class DestAccessor>
561 void
562 findHomogeneousRegions(
563  SrcIterator sul, SrcIterator slr, SrcAccessor src,
564  DestIterator dul, DestAccessor dest)
565 {
566  localMinima(sul, slr, src, dul, dest);
567 }
568 
569 template <class Vector1, class Vector2>
570 void noiseVarianceListMedianCut(Vector1 const & noise, Vector2 & clusters,
571  unsigned int maxClusterCount)
572 {
573  typedef typename Vector2::value_type Result;
574 
575  clusters.push_back(Result(0, noise.size()));
576 
577  while(clusters.size() <= maxClusterCount)
578  {
579  // find biggest cluster
580  unsigned int kMax = 0;
581  double diffMax = 0.0;
582  for(unsigned int k=0; k < clusters.size(); ++k)
583  {
584  int k1 = clusters[k][0], k2 = clusters[k][1]-1;
585 
586 #if 0 // turned the "internal error" in a postcondition message
587  // for the most likely case
588  std::string message("noiseVarianceListMedianCut(): internal error (");
589  message += std::string("k: ") + asString(k) + ", ";
590  message += std::string("k1: ") + asString(k1) + ", ";
591  message += std::string("k2: ") + asString(k2) + ", ";
592  message += std::string("noise.size(): ") + asString(noise.size()) + ", ";
593  message += std::string("clusters.size(): ") + asString(clusters.size()) + ").";
594  vigra_invariant(k1 >= 0 && k1 < (int)noise.size() && k2 >= 0 && k2 < (int)noise.size(), message.c_str());
595 #endif
596 
597  vigra_postcondition(k1 >= 0 && k1 < (int)noise.size() &&
598  k2 >= 0 && k2 < (int)noise.size(),
599  "noiseVarianceClustering(): Unable to find homogeneous regions.");
600 
601  double diff = noise[k2][0] - noise[k1][0];
602  if(diff > diffMax)
603  {
604  diffMax = diff;
605  kMax = k;
606  }
607  }
608 
609  if(diffMax == 0.0)
610  return; // all clusters have only one value
611 
612  unsigned int k1 = clusters[kMax][0],
613  k2 = clusters[kMax][1];
614  unsigned int kSplit = k1 + (k2 - k1) / 2;
615  clusters[kMax][1] = kSplit;
616  clusters.push_back(Result(kSplit, k2));
617  }
618 }
619 
620 struct SortNoiseByMean
621 {
622  template <class T>
623  bool operator()(T const & l, T const & r) const
624  {
625  return l[0] < r[0];
626  }
627 };
628 
629 struct SortNoiseByVariance
630 {
631  template <class T>
632  bool operator()(T const & l, T const & r) const
633  {
634  return l[1] < r[1];
635  }
636 };
637 
638 template <class Vector1, class Vector2, class Vector3>
639 void noiseVarianceClusterAveraging(Vector1 & noise, Vector2 & clusters,
640  Vector3 & result, double quantile)
641 {
642  typedef typename Vector1::iterator Iter;
643  typedef typename Vector3::value_type Result;
644 
645  for(unsigned int k=0; k<clusters.size(); ++k)
646  {
647  Iter i1 = noise.begin() + clusters[k][0];
648  Iter i2 = noise.begin() + clusters[k][1];
649 
650  std::sort(i1, i2, SortNoiseByVariance());
651 
652  std::size_t size = static_cast<std::size_t>(VIGRA_CSTD::ceil(quantile*(i2 - i1)));
653  if(static_cast<std::size_t>(i2 - i1) < size)
654  size = i2 - i1;
655  if(size < 1)
656  size = 1;
657  i2 = i1 + size;
658 
659  double mean = 0.0,
660  variance = 0.0;
661  for(; i1 < i2; ++i1)
662  {
663  mean += (*i1)[0];
664  variance += (*i1)[1];
665  }
666 
667  result.push_back(Result(mean / size, variance / size));
668  }
669 }
670 
671 template <class SrcIterator, class SrcAccessor, class BackInsertable>
672 void noiseVarianceEstimationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
673  BackInsertable & result,
674  NoiseNormalizationOptions const & options)
675 {
676  typedef typename BackInsertable::value_type ResultType;
677 
678  unsigned int w = slr.x - sul.x;
679  unsigned int h = slr.y - sul.y;
680 
681  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
682  typedef BasicImage<TmpType> TmpImage;
683 
684  TmpImage gradient(w, h);
685  symmetricDifferenceSquaredMagnitude(sul, slr, src, gradient.upperLeft(), gradient.accessor());
686 
687  BImage homogeneous(w, h);
688  findHomogeneousRegions(gradient.upperLeft(), gradient.lowerRight(), gradient.accessor(),
689  homogeneous.upperLeft(), homogeneous.accessor());
690 
691  // Generate noise of each of the remaining pixels == centers of homogeneous areas (border is not used)
692  unsigned int windowRadius = options.window_radius;
693  for(unsigned int y=windowRadius; y<h-windowRadius; ++y)
694  {
695  for(unsigned int x=windowRadius; x<w-windowRadius; ++x)
696  {
697  if (! homogeneous(x, y))
698  continue;
699 
700  Diff2D center(x, y);
701  double mean = 0.0, variance = options.noise_variance_initial_guess;
702 
703  bool success;
704 
705  if(options.use_gradient)
706  {
707  success = iterativeNoiseEstimationChi2(sul + center, src,
708  gradient.upperLeft() + center, mean, variance,
709  options.noise_estimation_quantile, windowRadius);
710  }
711  else
712  {
713  success = iterativeNoiseEstimationGauss(sul + center, src,
714  gradient.upperLeft() + center, mean, variance,
715  options.noise_estimation_quantile, windowRadius);
716  }
717  if (success)
718  {
719  result.push_back(ResultType(mean, variance));
720  }
721  }
722  }
723 }
724 
725 template <class Vector, class BackInsertable>
726 void noiseVarianceClusteringImpl(Vector & noise, BackInsertable & result,
727  unsigned int clusterCount, double quantile)
728 {
729  std::sort(noise.begin(), noise.end(), detail::SortNoiseByMean());
730 
731  ArrayVector<TinyVector<unsigned int, 2> > clusters;
732  detail::noiseVarianceListMedianCut(noise, clusters, clusterCount);
733 
734  std::sort(clusters.begin(), clusters.end(), detail::SortNoiseByMean());
735 
736  detail::noiseVarianceClusterAveraging(noise, clusters, result, quantile);
737 }
738 
739 template <class Functor,
740  class SrcIterator, class SrcAccessor,
741  class DestIterator, class DestAccessor>
742 bool
743 noiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
744  DestIterator dul, DestAccessor dest,
745  NoiseNormalizationOptions const & options)
746 {
747  ArrayVector<TinyVector<double, 2> > noiseData;
748  noiseVarianceEstimationImpl(sul, slr, src, noiseData, options);
749 
750  if(noiseData.size() < 10)
751  return false;
752 
753  ArrayVector<TinyVector<double, 2> > noiseClusters;
754 
755  noiseVarianceClusteringImpl(noiseData, noiseClusters,
756  options.cluster_count, options.averaging_quantile);
757 
758  transformImage(sul, slr, src, dul, dest, Functor(noiseClusters));
759 
760  return true;
761 }
762 
763 template <class SrcIterator, class SrcAccessor,
764  class DestIterator, class DestAccessor>
765 bool
766 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
767  DestIterator dul, DestAccessor dest,
768  NoiseNormalizationOptions const & options,
769  VigraTrueType /* isScalar */)
770 {
771  typedef typename SrcAccessor::value_type SrcType;
772  typedef typename DestAccessor::value_type DestType;
773  return noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
774  (sul, slr, src, dul, dest, options);
775 }
776 
777 template <class SrcIterator, class SrcAccessor,
778  class DestIterator, class DestAccessor>
779 bool
780 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
781  DestIterator dul, DestAccessor dest,
782  NoiseNormalizationOptions const & options,
783  VigraFalseType /* isScalar */)
784 {
785  int bands = src.size(sul);
786  for(int b=0; b<bands; ++b)
787  {
788  VectorElementAccessor<SrcAccessor> sband(b, src);
789  VectorElementAccessor<DestAccessor> dband(b, dest);
790  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
791  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
792 
793  if(!noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
794  (sul, slr, sband, dul, dband, options))
795  return false;
796  }
797  return true;
798 }
799 
800 template <class SrcIterator, class SrcAccessor,
801  class DestIterator, class DestAccessor>
802 bool
803 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
804  DestIterator dul, DestAccessor dest,
805  NoiseNormalizationOptions const & options,
806  VigraTrueType /* isScalar */)
807 {
808  typedef typename SrcAccessor::value_type SrcType;
809  typedef typename DestAccessor::value_type DestType;
810  return noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
811  (sul, slr, src, dul, dest, options);
812 }
813 
814 template <class SrcIterator, class SrcAccessor,
815  class DestIterator, class DestAccessor>
816 bool
817 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
818  DestIterator dul, DestAccessor dest,
819  NoiseNormalizationOptions const & options,
820  VigraFalseType /* isScalar */)
821 {
822  int bands = src.size(sul);
823  for(int b=0; b<bands; ++b)
824  {
825  VectorElementAccessor<SrcAccessor> sband(b, src);
826  VectorElementAccessor<DestAccessor> dband(b, dest);
827  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
828  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
829 
830  if(!noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
831  (sul, slr, sband, dul, dband, options))
832  return false;
833  }
834  return true;
835 }
836 
837 template <class SrcIterator, class SrcAccessor,
838  class DestIterator, class DestAccessor>
839 void
840 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
841  DestIterator dul, DestAccessor dest,
842  double a0, double a1, double a2,
843  VigraTrueType /* isScalar */)
844 {
845  ArrayVector<TinyVector<double, 2> > noiseClusters;
846  noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
847  noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1 + a2));
848  noiseClusters.push_back(TinyVector<double, 2>(2.0, a0 + 2.0*a1 + 4.0*a2));
849  transformImage(sul, slr, src, dul, dest,
850  QuadraticNoiseNormalizationFunctor<typename SrcAccessor::value_type,
851  typename DestAccessor::value_type>(noiseClusters));
852 }
853 
854 template <class SrcIterator, class SrcAccessor,
855  class DestIterator, class DestAccessor>
856 void
857 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
858  DestIterator dul, DestAccessor dest,
859  double a0, double a1, double a2,
860  VigraFalseType /* isScalar */)
861 {
862  int bands = src.size(sul);
863  for(int b=0; b<bands; ++b)
864  {
865  VectorElementAccessor<SrcAccessor> sband(b, src);
866  VectorElementAccessor<DestAccessor> dband(b, dest);
867  quadraticNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, a2, VigraTrueType());
868  }
869 }
870 
871 template <class SrcIterator, class SrcAccessor,
872  class DestIterator, class DestAccessor>
873 bool
874 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
875  DestIterator dul, DestAccessor dest,
876  NoiseNormalizationOptions const & options,
877  VigraTrueType /* isScalar */)
878 {
879  typedef typename SrcAccessor::value_type SrcType;
880  typedef typename DestAccessor::value_type DestType;
881  return noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
882  (sul, slr, src, dul, dest, options);
883 }
884 
885 template <class SrcIterator, class SrcAccessor,
886  class DestIterator, class DestAccessor>
887 bool
888 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
889  DestIterator dul, DestAccessor dest,
890  NoiseNormalizationOptions const & options,
891  VigraFalseType /* isScalar */)
892 {
893  int bands = src.size(sul);
894  for(int b=0; b<bands; ++b)
895  {
896  VectorElementAccessor<SrcAccessor> sband(b, src);
897  VectorElementAccessor<DestAccessor> dband(b, dest);
898  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
899  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
900 
901  if(!noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
902  (sul, slr, sband, dul, dband, options))
903  return false;
904  }
905  return true;
906 }
907 
908 template <class SrcIterator, class SrcAccessor,
909  class DestIterator, class DestAccessor>
910 void
911 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
912  DestIterator dul, DestAccessor dest,
913  double a0, double a1,
914  VigraTrueType /* isScalar */)
915 {
916  ArrayVector<TinyVector<double, 2> > noiseClusters;
917  noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
918  noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1));
919  transformImage(sul, slr, src, dul, dest,
920  LinearNoiseNormalizationFunctor<typename SrcAccessor::value_type,
921  typename DestAccessor::value_type>(noiseClusters));
922 }
923 
924 template <class SrcIterator, class SrcAccessor,
925  class DestIterator, class DestAccessor>
926 void
927 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
928  DestIterator dul, DestAccessor dest,
929  double a0, double a1,
930  VigraFalseType /* isScalar */)
931 {
932  int bands = src.size(sul);
933  for(int b=0; b<bands; ++b)
934  {
935  VectorElementAccessor<SrcAccessor> sband(b, src);
936  VectorElementAccessor<DestAccessor> dband(b, dest);
937  linearNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, VigraTrueType());
938  }
939 }
940 
941 } // namespace detail
942 
943 template <bool P>
944 struct noiseVarianceEstimation_can_only_work_on_scalar_images
945 : vigra::staticAssert::AssertBool<P>
946 {};
947 
948 /** \addtogroup NoiseNormalization Noise Normalization
949  Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
950 */
951 //@{
952 
953 /********************************************************/
954 /* */
955 /* noiseVarianceEstimation */
956 /* */
957 /********************************************************/
958 
959 /** \brief Determine the noise variance as a function of the image intensity.
960 
961  This operator applies an algorithm described in
962 
963  W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
964  Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
965  Lecture Notes in Earth Science, Berlin: Springer, 1999
966 
967  in order to estimate the noise variance as a function of the image intensity in a robust way,
968  i.e. so that intensity changes due to edges do not bias the estimate. The source value type
969  (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
970  The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible
971  from two <tt>double</tt> values. The following options can be set via the \a options object
972  (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
973 
974  <tt>useGradient</tt>, <tt>windowRadius</tt>, <tt>noiseEstimationQuantile</tt>, <tt>noiseVarianceInitialGuess</tt>
975 
976  <b> Declarations:</b>
977 
978  pass 2D array views:
979  \code
980  namespace vigra {
981  template <class T1, class S1, class BackInsertable>
982  void
983  noiseVarianceEstimation(MultiArrayView<2, T1, S1> const & src,
984  BackInsertable & result,
985  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
986  }
987  \endcode
988 
989  \deprecatedAPI{noiseVarianceEstimation}
990  pass \ref ImageIterators and \ref DataAccessors :
991  \code
992  namespace vigra {
993  template <class SrcIterator, class SrcAccessor, class BackInsertable>
994  void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
995  BackInsertable & result,
996  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
997  }
998  \endcode
999  use argument objects in conjunction with \ref ArgumentObjectFactories :
1000  \code
1001  namespace vigra {
1002  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1003  void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1004  BackInsertable & result,
1005  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1006  }
1007  \endcode
1008  \deprecatedEnd
1009 
1010  <b> Usage:</b>
1011 
1012  <b>\#include</b> <vigra/noise_normalization.hxx><br>
1013  Namespace: vigra
1014 
1015  \code
1016  MultiArray<2, float> src(w,h);
1017  std::vector<TinyVector<double, 2> > result;
1018 
1019  ...
1020  noiseVarianceEstimation(src, result,
1021  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
1022 
1023  // print the intensity / variance pairs found
1024  for(int k=0; k<result.size(); ++k)
1025  std::cout << "Intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
1026  \endcode
1027 
1028  \deprecatedUsage{noiseVarianceEstimation}
1029  \code
1030  vigra::BImage src(w,h);
1031  std::vector<vigra::TinyVector<double, 2> > result;
1032 
1033  ...
1034  vigra::noiseVarianceEstimation(srcImageRange(src), result,
1035  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
1036 
1037  // print the intensity / variance pairs found
1038  for(int k=0; k<result.size(); ++k)
1039  std::cout << "Intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
1040  \endcode
1041  <b> Required Interface:</b>
1042  \code
1043  SrcIterator upperleft, lowerright;
1044  SrcAccessor src;
1045 
1046  typedef SrcAccessor::value_type SrcType;
1047  typedef NumericTraits<SrcType>::isScalar isScalar;
1048  assert(isScalar::asBool == true);
1049 
1050  double value = src(uperleft);
1051 
1052  BackInsertable result;
1053  typedef BackInsertable::value_type ResultType;
1054  double intensity, variance;
1055  result.push_back(ResultType(intensity, variance));
1056  \endcode
1057  \deprecatedEnd
1058 */
1060 
1061 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1062 inline
1063 void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1064  BackInsertable & result,
1065  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1066 {
1067  typedef typename SrcAccessor::value_type SrcType;
1068  typedef typename NumericTraits<SrcType>::isScalar isScalar;
1069 
1070  VIGRA_STATIC_ASSERT((
1071  noiseVarianceEstimation_can_only_work_on_scalar_images<(isScalar::asBool)>));
1072 
1073  detail::noiseVarianceEstimationImpl(sul, slr, src, result, options);
1074 }
1075 
1076 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1077 inline void
1078 noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1079  BackInsertable & result,
1080  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1081 {
1082  noiseVarianceEstimation(src.first, src.second, src.third, result, options);
1083 }
1084 
1085 template <class T1, class S1, class BackInsertable>
1086 inline void
1087 noiseVarianceEstimation(MultiArrayView<2, T1, S1> const & src,
1088  BackInsertable & result,
1089  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1090 {
1091  noiseVarianceEstimation(srcImageRange(src), result, options);
1092 }
1093 
1094 /********************************************************/
1095 /* */
1096 /* noiseVarianceClustering */
1097 /* */
1098 /********************************************************/
1099 
1100 /** \brief Determine the noise variance as a function of the image intensity and cluster the results.
1101 
1102  This operator first calls \ref noiseVarianceEstimation() to obtain a sequence of intensity/variance pairs,
1103  which are then clustered using the median cut algorithm. Then the cluster centers (i.e. average variance vs.
1104  average intensity) are determined and returned in the \a result sequence.
1105 
1106  In addition to the options valid for \ref noiseVarianceEstimation(), the following options can be set via
1107  the \a options object (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
1108 
1109  <tt>clusterCount</tt>, <tt>averagingQuantile</tt>
1110 
1111  <b> Declarations:</b>
1112 
1113  pass 2D array views:
1114  \code
1115  namespace vigra {
1116  template <class T1, class S1, class BackInsertable>
1117  void
1118  noiseVarianceClustering(MultiArrayView<2, T1, S1> const & src,
1119  BackInsertable & result,
1120  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1121  }
1122  \endcode
1123 
1124  \deprecatedAPI{noiseVarianceClustering}
1125  pass \ref ImageIterators and \ref DataAccessors :
1126  \code
1127  namespace vigra {
1128  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1129  void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1130  BackInsertable & result,
1131  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1132  }
1133  \endcode
1134  use argument objects in conjunction with \ref ArgumentObjectFactories :
1135  \code
1136  namespace vigra {
1137  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1138  void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1139  BackInsertable & result,
1140  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1141  }
1142  \endcode
1143  \deprecatedEnd
1144 
1145  <b> Usage:</b>
1146 
1147  <b>\#include</b> <vigra/noise_normalization.hxx><br>
1148  Namespace: vigra
1149 
1150  \code
1151  MultiArray<2, float> src(w,h);
1152  std::vector<TinyVector<double, 2> > result;
1153 
1154  ...
1155  noiseVarianceClustering(src, result,
1156  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1157  clusterCount(15));
1158 
1159  // print the intensity / variance pairs representing the cluster centers
1160  for(int k=0; k<result.size(); ++k)
1161  std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
1162  \endcode
1163 
1164  \deprecatedUsage{noiseVarianceClustering}
1165  \code
1166  vigra::BImage src(w,h);
1167  std::vector<vigra::TinyVector<double, 2> > result;
1168 
1169  ...
1170  vigra::noiseVarianceClustering(srcImageRange(src), result,
1171  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1172  clusterCount(15));
1173 
1174  // print the intensity / variance pairs representing the cluster centers
1175  for(int k=0; k<result.size(); ++k)
1176  std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
1177  \endcode
1178  <b> Required Interface:</b>
1179  same as \ref noiseVarianceEstimation()
1180  \deprecatedEnd
1181 */
1183 
1184 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1185 inline
1186 void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1187  BackInsertable & result,
1188  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1189 {
1190  ArrayVector<TinyVector<double, 2> > variance;
1191  noiseVarianceEstimation(sul, slr, src, variance, options);
1192  detail::noiseVarianceClusteringImpl(variance, result, options.cluster_count, options.averaging_quantile);
1193 }
1194 
1195 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1196 inline void
1197 noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1198  BackInsertable & result,
1199  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1200 {
1201  noiseVarianceClustering(src.first, src.second, src.third, result, options);
1202 }
1203 
1204 template <class T1, class S1, class BackInsertable>
1205 inline void
1206 noiseVarianceClustering(MultiArrayView<2, T1, S1> const & src,
1207  BackInsertable & result,
1208  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1209 {
1210  noiseVarianceClustering(srcImageRange(src), result, options);
1211 }
1212 
1213 /********************************************************/
1214 /* */
1215 /* nonparametricNoiseNormalization */
1216 /* */
1217 /********************************************************/
1218 
1219 /** \brief Noise normalization by means of an estimated non-parametric noise model.
1220 
1221  The original image is assumed to be corrupted by noise whose variance depends on the intensity in an unknown way.
1222  The present functions first calls \ref noiseVarianceClustering() to obtain a sequence of intensity/variance pairs
1223  (cluster centers) which estimate this dependency. The cluster centers are connected into a piecewise linear
1224  function which is the inverted according to the formula derived in
1225 
1226  W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
1227  Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
1228  Lecture Notes in Earth Science, Berlin: Springer, 1999
1229 
1230  The inverted formula defines a pixel-wise intensity transformation whose application turns the original image
1231  into one that is corrupted by additive Gaussian noise with unit variance. Most subsequent algorithms will be able
1232  to handle this type of noise much better than the original noise.
1233 
1234  RGB and other multiband images will be processed one band at a time. The function returns <tt>true</tt> on success.
1235  Noise normalization will fail if the original image does not contain sufficiently homogeneous regions that
1236  allow robust estimation of the noise variance.
1237 
1238  The \a options object may use all options described in \ref vigra::NoiseNormalizationOptions.
1239 
1240  The function returns <tt>false</tt> if the noise estimation failed, so that no normalization could be performed.
1241 
1242  <b> Declarations:</b>
1243 
1244  pass 2D array views:
1245  \code
1246  namespace vigra {
1247  template <class T1, class S1,
1248  class T2, class S2>
1249  bool
1250  nonparametricNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1251  MultiArrayView<2, T2, S2> dest,
1252  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1253  }
1254  \endcode
1255 
1256  \deprecatedAPI{nonparametricNoiseNormalization}
1257  pass \ref ImageIterators and \ref DataAccessors :
1258  \code
1259  namespace vigra {
1260  template <class SrcIterator, class SrcAccessor,
1261  class DestIterator, class DestAccessor>
1262  bool nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1263  DestIterator dul, DestAccessor dest,
1264  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1265  }
1266  \endcode
1267  use argument objects in conjunction with \ref ArgumentObjectFactories :
1268  \code
1269  namespace vigra {
1270  template <class SrcIterator, class SrcAccessor,
1271  class DestIterator, class DestAccessor>
1272  bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1273  pair<DestIterator, DestAccessor> dest,
1274  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1275  }
1276  \endcode
1277  \deprecatedEnd
1278 
1279  <b> Usage:</b>
1280 
1281  <b>\#include</b> <vigra/noise_normalization.hxx><br>
1282  Namespace: vigra
1283 
1284  \code
1285  MultiArray<2, RGBValue<float> > src(w,h), dest(w, h);
1286 
1287  ...
1288  nonparametricNoiseNormalization(src, dest,
1289  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1290  clusterCount(15));
1291  \endcode
1292 
1293  \deprecatedUsage{nonparametricNoiseNormalization}
1294  \code
1295  vigra::BRGBImage src(w,h), dest(w, h);
1296 
1297  ...
1298  vigra::nonparametricNoiseNormalization(srcImageRange(src), destImage(dest),
1299  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1300  clusterCount(15));
1301  \endcode
1302  <b> Required Interface:</b>
1303  same as \ref noiseVarianceEstimation()
1304  \deprecatedEnd
1305 */
1307 
1308 template <class SrcIterator, class SrcAccessor,
1309  class DestIterator, class DestAccessor>
1310 inline bool
1311 nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1312  DestIterator dul, DestAccessor dest,
1313  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1314 {
1315  typedef typename SrcAccessor::value_type SrcType;
1316 
1317  return detail::nonparametricNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1318  typename NumericTraits<SrcType>::isScalar());
1319 }
1320 
1321 template <class SrcIterator, class SrcAccessor,
1322  class DestIterator, class DestAccessor>
1323 inline bool
1324 nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1325  pair<DestIterator, DestAccessor> dest,
1326  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1327 {
1328  return nonparametricNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1329 }
1330 
1331 template <class T1, class S1,
1332  class T2, class S2>
1333 inline bool
1334 nonparametricNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1335  MultiArrayView<2, T2, S2> dest,
1336  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1337 {
1338  vigra_precondition(src.shape() == dest.shape(),
1339  "nonparametricNoiseNormalization(): shape mismatch between input and output.");
1340  return nonparametricNoiseNormalization(srcImageRange(src), destImage(dest), options);
1341 }
1342 
1343 /********************************************************/
1344 /* */
1345 /* quadraticNoiseNormalization */
1346 /* */
1347 /********************************************************/
1348 
1349 /** \brief Noise normalization by means of an estimated or given quadratic noise model.
1350 
1351  This function works like \ref nonparametricNoiseNormalization() excapt that the model for the
1352  dependency between intensity and noise variance is assumed to be a
1353  quadratic function rather than a piecewise linear function. If the data conform to the quadratic model,
1354  this leads to a somewhat smoother transformation. The function returns <tt>false</tt> if the noise
1355  estimation failed, so that no normalization could be performed.
1356 
1357  In the second variant of the function, the parameters of the quadratic model are not estimated,
1358  but explicitly given according to:
1359  \code
1360  variance = a0 + a1 * intensity + a2 * sq(intensity)
1361  \endcode
1362 
1363  <b> Declarations:</b>
1364 
1365  pass 2D array views:
1366  \code
1367  namespace vigra {
1368  // estimate and apply a quadratic noise model
1369  template <class T1, class S1,
1370  class T2, class S2>
1371  bool
1372  quadraticNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1373  MultiArrayView<2, T2, S2> dest,
1374  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1375 
1376  // apply a given quadratic noise model
1377  template <class T1, class S1,
1378  class T2, class S2>
1379  void
1380  quadraticNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1381  MultiArrayView<2, T2, S2> dest,
1382  double a0, double a1, double a2);
1383  }
1384  \endcode
1385 
1386  \deprecatedAPI{quadraticNoiseNormalization}
1387  pass \ref ImageIterators and \ref DataAccessors :
1388  \code
1389  namespace vigra {
1390  // estimate and apply a quadratic noise model
1391  template <class SrcIterator, class SrcAccessor,
1392  class DestIterator, class DestAccessor>
1393  bool quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1394  DestIterator dul, DestAccessor dest,
1395  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1396 
1397  // apply a given quadratic noise model
1398  template <class SrcIterator, class SrcAccessor,
1399  class DestIterator, class DestAccessor>
1400  void quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1401  DestIterator dul, DestAccessor dest,
1402  double a0, double a1, double a2);
1403  }
1404  \endcode
1405  use argument objects in conjunction with \ref ArgumentObjectFactories :
1406  \code
1407  namespace vigra {
1408  // estimate and apply a quadratic noise model
1409  template <class SrcIterator, class SrcAccessor,
1410  class DestIterator, class DestAccessor>
1411  bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1412  pair<DestIterator, DestAccessor> dest,
1413  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1414 
1415  // apply a given quadratic noise model
1416  template <class SrcIterator, class SrcAccessor,
1417  class DestIterator, class DestAccessor>
1418  void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1419  pair<DestIterator, DestAccessor> dest,
1420  double a0, double a1, double a2);
1421  }
1422  \endcode
1423  \deprecatedEnd
1424 
1425  <b> Usage:</b>
1426 
1427  <b>\#include</b> <vigra/noise_normalization.hxx><br>
1428  Namespace: vigra
1429 
1430  \code
1431  MultiArray<2, RGBValue<float> > src(w,h), dest(w, h);
1432 
1433  ...
1434  // estimate the noise model and apply it to normalize the noise variance
1435  bool success = quadraticNoiseNormalization(src, dest,
1436  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1437  clusterCount(15));
1438  vigra_postcondition(success, "quadraticNoiseNormalization(): Unable to estimate noise model.");
1439 
1440  // apply a pre-computed noise model
1441  quadraticNoiseNormalization(src, dest, 100, 0.02, 1e-6);
1442  \endcode
1443 
1444  \deprecatedUsage{quadraticNoiseNormalization}
1445  \code
1446  vigra::BRGBImage src(w,h), dest(w, h);
1447 
1448  ...
1449  // estimate the noise model and apply it to normalize the noise variance
1450  vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
1451  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1452  clusterCount(15));
1453 
1454  // apply a pre-computed noise model
1455  vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
1456  100, 0.02, 1e-6);
1457  \endcode
1458  \deprecatedEnd
1459 
1460  <b> Required Interface:</b>
1461 
1462  The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
1463  are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
1464  or a vector whose elements are assignable from <tt>double</tt>.
1465 */
1467 
1468 template <class SrcIterator, class SrcAccessor,
1469  class DestIterator, class DestAccessor>
1470 inline bool
1471 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1472  DestIterator dul, DestAccessor dest,
1473  NoiseNormalizationOptions const & options)
1474 {
1475  typedef typename SrcAccessor::value_type SrcType;
1476 
1477  return detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1478  typename NumericTraits<SrcType>::isScalar());
1479 }
1480 
1481 template <class SrcIterator, class SrcAccessor,
1482  class DestIterator, class DestAccessor>
1483 inline bool
1484 quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1485  pair<DestIterator, DestAccessor> dest,
1486  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1487 {
1488  return quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1489 }
1490 
1491 template <class T1, class S1,
1492  class T2, class S2>
1493 inline bool
1494 quadraticNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1495  MultiArrayView<2, T2, S2> dest,
1496  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1497 {
1498  vigra_precondition(src.shape() == dest.shape(),
1499  "quadraticNoiseNormalization(): shape mismatch between input and output.");
1500  return quadraticNoiseNormalization(srcImageRange(src), destImage(dest), options);
1501 }
1502 
1503 /********************************************************/
1504 /* */
1505 /* quadraticNoiseNormalization */
1506 /* (variant) */
1507 /* */
1508 /********************************************************/
1509 
1510 template <class SrcIterator, class SrcAccessor,
1511  class DestIterator, class DestAccessor>
1512 inline void
1513 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1514  DestIterator dul, DestAccessor dest,
1515  double a0, double a1, double a2)
1516 {
1517  typedef typename SrcAccessor::value_type SrcType;
1518 
1519  detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1, a2,
1520  typename NumericTraits<SrcType>::isScalar());
1521 }
1522 
1523 template <class SrcIterator, class SrcAccessor,
1524  class DestIterator, class DestAccessor>
1525 inline void
1526 quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1527  pair<DestIterator, DestAccessor> dest,
1528  double a0, double a1, double a2)
1529 {
1530  quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1, a2);
1531 }
1532 
1533 template <class T1, class S1,
1534  class T2, class S2>
1535 inline void
1536 quadraticNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1537  MultiArrayView<2, T2, S2> dest,
1538  double a0, double a1, double a2)
1539 {
1540  vigra_precondition(src.shape() == dest.shape(),
1541  "quadraticNoiseNormalization(): shape mismatch between input and output.");
1542  quadraticNoiseNormalization(srcImageRange(src), destImage(dest), a0, a1, a2);
1543 }
1544 
1545 /********************************************************/
1546 /* */
1547 /* linearNoiseNormalization */
1548 /* */
1549 /********************************************************/
1550 
1551 /** \brief Noise normalization by means of an estimated or given linear noise model.
1552 
1553  This function works like \ref nonparametricNoiseNormalization() excapt that the model for the
1554  dependency between intensity and noise variance is assumed to be a
1555  linear function rather than a piecewise linear function. If the data conform to the linear model,
1556  this leads to a very simple transformation which is similar to the familiar gamma correction.
1557  The function returns <tt>false</tt> if the noise estimation failed, so that no
1558  normalization could be performed.
1559 
1560  In the second variant of the function, the parameters of the linear model are not estimated,
1561  but explicitly given according to:
1562  \code
1563  variance = a0 + a1 * intensity
1564  \endcode
1565 
1566  <b> Declarations:</b>
1567 
1568  pass 2D array views:
1569  \code
1570  namespace vigra {
1571  // estimate and apply a linear noise model
1572  template <class T1, class S1,
1573  class T2, class S2>
1574  bool
1575  linearNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1576  MultiArrayView<2, T2, S2> dest,
1577  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1578 
1579  // apply a given linear noise model
1580  template <class T1, class S1,
1581  class T2, class S2>
1582  void
1583  linearNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1584  MultiArrayView<2, T2, S2> dest,
1585  double a0, double a1);
1586  }
1587  \endcode
1588 
1589  \deprecatedAPI{linearNoiseNormalization}
1590  pass \ref ImageIterators and \ref DataAccessors :
1591  \code
1592  namespace vigra {
1593  // estimate and apply a linear noise model
1594  template <class SrcIterator, class SrcAccessor,
1595  class DestIterator, class DestAccessor>
1596  bool linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1597  DestIterator dul, DestAccessor dest,
1598  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1599 
1600  // apply a given linear noise model
1601  template <class SrcIterator, class SrcAccessor,
1602  class DestIterator, class DestAccessor>
1603  void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1604  DestIterator dul, DestAccessor dest,
1605  double a0, double a1);
1606  }
1607  \endcode
1608  use argument objects in conjunction with \ref ArgumentObjectFactories :
1609  \code
1610  namespace vigra {
1611  // estimate and apply a linear noise model
1612  template <class SrcIterator, class SrcAccessor,
1613  class DestIterator, class DestAccessor>
1614  bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1615  pair<DestIterator, DestAccessor> dest,
1616  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1617 
1618  // apply a given linear noise model
1619  template <class SrcIterator, class SrcAccessor,
1620  class DestIterator, class DestAccessor>
1621  void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1622  pair<DestIterator, DestAccessor> dest,
1623  double a0, double a1);
1624  }
1625  \endcode
1626  \deprecatedEnd
1627 
1628  <b> Usage:</b>
1629 
1630  <b>\#include</b> <vigra/noise_normalization.hxx><br>
1631  Namespace: vigra
1632 
1633  \code
1634  vigra::BRGBImage src(w,h), dest(w, h);
1635 
1636  ...
1637  // estimate the noise model and apply it to normalize the noise variance
1638  bool success = linearNoiseNormalization(src, dest,
1639  NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1640  clusterCount(15));
1641  vigra_postcondition(success, "linearNoiseNormalization(): Unable to estimate noise model.");
1642 
1643  // apply a pre-computed noise model
1644  linearNoiseNormalization(src, dest, 100, 0.02);
1645  \endcode
1646 
1647  \deprecatedUsage{linearNoiseNormalization}
1648  \code
1649  vigra::BRGBImage src(w,h), dest(w, h);
1650 
1651  ...
1652  // estimate the noise model and apply it to normalize the noise variance
1653  vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
1654  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1655  clusterCount(15));
1656 
1657  // apply a pre-computed noise model
1658  vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
1659  100, 0.02);
1660  \endcode
1661  \deprecatedEnd
1662 
1663  <b> Required Interface:</b>
1664 
1665  The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
1666  are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
1667  or a vector whose elements are assignable from <tt>double</tt>.
1668 */
1670 
1671 template <class SrcIterator, class SrcAccessor,
1672  class DestIterator, class DestAccessor>
1673 inline bool
1674 linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1675  DestIterator dul, DestAccessor dest,
1676  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1677 {
1678  typedef typename SrcAccessor::value_type SrcType;
1679 
1680  return detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1681  typename NumericTraits<SrcType>::isScalar());
1682 }
1683 
1684 template <class SrcIterator, class SrcAccessor,
1685  class DestIterator, class DestAccessor>
1686 inline bool
1687 linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1688  pair<DestIterator, DestAccessor> dest,
1689  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1690 {
1691  return linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1692 }
1693 
1694 template <class T1, class S1,
1695  class T2, class S2>
1696 inline bool
1697 linearNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1698  MultiArrayView<2, T2, S2> dest,
1699  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1700 {
1701  vigra_precondition(src.shape() == dest.shape(),
1702  "linearNoiseNormalization(): shape mismatch between input and output.");
1703  return linearNoiseNormalization(srcImageRange(src), destImage(dest), options);
1704 }
1705 
1706 /********************************************************/
1707 /* */
1708 /* linearNoiseNormalization */
1709 /* (variant) */
1710 /* */
1711 /********************************************************/
1712 
1713 template <class SrcIterator, class SrcAccessor,
1714  class DestIterator, class DestAccessor>
1715 inline
1716 void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1717  DestIterator dul, DestAccessor dest,
1718  double a0, double a1)
1719 {
1720  typedef typename SrcAccessor::value_type SrcType;
1721 
1722  detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1,
1723  typename NumericTraits<SrcType>::isScalar());
1724 }
1725 
1726 template <class SrcIterator, class SrcAccessor,
1727  class DestIterator, class DestAccessor>
1728 inline void
1729 linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1730  pair<DestIterator, DestAccessor> dest,
1731  double a0, double a1)
1732 {
1733  linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1);
1734 }
1735 
1736 template <class T1, class S1,
1737  class T2, class S2>
1738 inline void
1739 linearNoiseNormalization(MultiArrayView<2, T1, S1> const & src,
1740  MultiArrayView<2, T2, S2> dest,
1741  double a0, double a1)
1742 {
1743  vigra_precondition(src.shape() == dest.shape(),
1744  "linearNoiseNormalization(): shape mismatch between input and output.");
1745  linearNoiseNormalization(srcImageRange(src), destImage(dest), a0, a1);
1746 }
1747 
1748 //@}
1749 
1750 } // namespace vigra
1751 
1752 #endif // VIGRA_NOISE_NORMALIZATION_HXX
NoiseNormalizationOptions & noiseEstimationQuantile(double quantile)
Definition: noise_normalization.hxx:161
detail::SelectIntegerType< 8, detail::UnsignedIntTypes >::type UInt8
8-bit unsigned int
Definition: sized_int.hxx:179
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
void transformImage(...)
Apply unary point transformation to each pixel.
void localMinima(...)
Find local minima in an image or multi-dimensional array.
bool nonparametricNoiseNormalization(...)
Noise normalization by means of an estimated non-parametric noise model.
linalg::TemporaryMatrix< T > asin(MultiArrayView< 2, T, C > const &v)
NoiseNormalizationOptions & noiseVarianceInitialGuess(double guess)
Definition: noise_normalization.hxx:174
detail::SelectIntegerType< 16, detail::UnsignedIntTypes >::type UInt16
16-bit unsigned int
Definition: sized_int.hxx:181
bool linearNoiseNormalization(...)
Noise normalization by means of an estimated or given linear noise model.
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
NoiseNormalizationOptions()
Definition: noise_normalization.hxx:95
void discErosion(...)
Apply erosion (minimum) filter with disc of given radius to image.
BasicImage< UInt8 > BImage
Definition: stdimage.hxx:62
NoiseNormalizationOptions & averagingQuantile(double quantile)
Definition: noise_normalization.hxx:147
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
NoiseNormalizationOptions & clusterCount(unsigned int c)
Definition: noise_normalization.hxx:132
std::string asString(T t)(...)
void combineTwoImages(...)
Combine two source images into destination image.
NoiseNormalizationOptions & useGradient(bool r)
Definition: noise_normalization.hxx:109
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1459
Pass options to one of the noise normalization functions.
Definition: noise_normalization.hxx:89
bool closeAtTolerance(T1 l, T2 r, typename PromoteTraits< T1, T2 >::Promote epsilon)
Tolerance based floating-point equality.
Definition: mathutil.hxx:1638
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
NoiseNormalizationOptions & windowRadius(unsigned int r)
Definition: noise_normalization.hxx:119
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
bool quadraticNoiseNormalization(...)
Noise normalization by means of an estimated or given quadratic noise model.
void noiseVarianceEstimation(...)
Determine the noise variance as a function of the image intensity.
void noiseVarianceClustering(...)
Determine the noise variance as a function of the image intensity and cluster the results...
void separableConvolveY(...)
Performs a 1 dimensional convolution in y direction.
bool linearSolve(...)
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)