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

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

© 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.7.1 (Wed Mar 12 2014)