[ 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  double diff = noise[clusters[k][1]-1][0] - noise[clusters[k][0]][0];
586  if(diff > diffMax)
587  {
588  diffMax = diff;
589  kMax = k;
590  }
591  }
592 
593  if(diffMax == 0.0)
594  return; // all clusters have only one value
595 
596  unsigned int k1 = clusters[kMax][0],
597  k2 = clusters[kMax][1];
598  unsigned int kSplit = k1 + (k2 - k1) / 2;
599  clusters[kMax][1] = kSplit;
600  clusters.push_back(Result(kSplit, k2));
601  }
602 }
603 
604 struct SortNoiseByMean
605 {
606  template <class T>
607  bool operator()(T const & l, T const & r) const
608  {
609  return l[0] < r[0];
610  }
611 };
612 
613 struct SortNoiseByVariance
614 {
615  template <class T>
616  bool operator()(T const & l, T const & r) const
617  {
618  return l[1] < r[1];
619  }
620 };
621 
622 template <class Vector1, class Vector2, class Vector3>
623 void noiseVarianceClusterAveraging(Vector1 & noise, Vector2 & clusters,
624  Vector3 & result, double quantile)
625 {
626  typedef typename Vector1::iterator Iter;
627  typedef typename Vector3::value_type Result;
628 
629  for(unsigned int k=0; k<clusters.size(); ++k)
630  {
631  Iter i1 = noise.begin() + clusters[k][0];
632  Iter i2 = noise.begin() + clusters[k][1];
633 
634  std::sort(i1, i2, SortNoiseByVariance());
635 
636  std::size_t size = static_cast<std::size_t>(VIGRA_CSTD::ceil(quantile*(i2 - i1)));
637  if(static_cast<std::size_t>(i2 - i1) < size)
638  size = i2 - i1;
639  if(size < 1)
640  size = 1;
641  i2 = i1 + size;
642 
643  double mean = 0.0,
644  variance = 0.0;
645  for(; i1 < i2; ++i1)
646  {
647  mean += (*i1)[0];
648  variance += (*i1)[1];
649  }
650 
651  result.push_back(Result(mean / size, variance / size));
652  }
653 }
654 
655 template <class SrcIterator, class SrcAccessor, class BackInsertable>
656 void noiseVarianceEstimationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
657  BackInsertable & result,
658  NoiseNormalizationOptions const & options)
659 {
660  typedef typename BackInsertable::value_type ResultType;
661 
662  unsigned int w = slr.x - sul.x;
663  unsigned int h = slr.y - sul.y;
664 
665  typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
666  typedef BasicImage<TmpType> TmpImage;
667 
668  TmpImage gradient(w, h);
669  symmetricDifferenceSquaredMagnitude(sul, slr, src, gradient.upperLeft(), gradient.accessor());
670 
671  BImage homogeneous(w, h);
672  findHomogeneousRegions(gradient.upperLeft(), gradient.lowerRight(), gradient.accessor(),
673  homogeneous.upperLeft(), homogeneous.accessor());
674 
675  // Generate noise of each of the remaining pixels == centers of homogenous areas (border is not used)
676  unsigned int windowRadius = options.window_radius;
677  for(unsigned int y=windowRadius; y<h-windowRadius; ++y)
678  {
679  for(unsigned int x=windowRadius; x<w-windowRadius; ++x)
680  {
681  if (! homogeneous(x, y))
682  continue;
683 
684  Diff2D center(x, y);
685  double mean = 0.0, variance = options.noise_variance_initial_guess;
686 
687  bool success;
688 
689  if(options.use_gradient)
690  {
691  success = iterativeNoiseEstimationChi2(sul + center, src,
692  gradient.upperLeft() + center, mean, variance,
693  options.noise_estimation_quantile, windowRadius);
694  }
695  else
696  {
697  success = iterativeNoiseEstimationGauss(sul + center, src,
698  gradient.upperLeft() + center, mean, variance,
699  options.noise_estimation_quantile, windowRadius);
700  }
701  if (success)
702  {
703  result.push_back(ResultType(mean, variance));
704  }
705  }
706  }
707 }
708 
709 template <class Vector, class BackInsertable>
710 void noiseVarianceClusteringImpl(Vector & noise, BackInsertable & result,
711  unsigned int clusterCount, double quantile)
712 {
713  std::sort(noise.begin(), noise.end(), detail::SortNoiseByMean());
714 
715  ArrayVector<TinyVector<unsigned int, 2> > clusters;
716  detail::noiseVarianceListMedianCut(noise, clusters, clusterCount);
717 
718  std::sort(clusters.begin(), clusters.end(), detail::SortNoiseByMean());
719 
720  detail::noiseVarianceClusterAveraging(noise, clusters, result, quantile);
721 }
722 
723 template <class Functor,
724  class SrcIterator, class SrcAccessor,
725  class DestIterator, class DestAccessor>
726 bool
727 noiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
728  DestIterator dul, DestAccessor dest,
729  NoiseNormalizationOptions const & options)
730 {
731  ArrayVector<TinyVector<double, 2> > noiseData;
732  noiseVarianceEstimationImpl(sul, slr, src, noiseData, options);
733 
734  if(noiseData.size() < 10)
735  return false;
736 
737  ArrayVector<TinyVector<double, 2> > noiseClusters;
738 
739  noiseVarianceClusteringImpl(noiseData, noiseClusters,
740  options.cluster_count, options.averaging_quantile);
741 
742  transformImage(sul, slr, src, dul, dest, Functor(noiseClusters));
743 
744  return true;
745 }
746 
747 template <class SrcIterator, class SrcAccessor,
748  class DestIterator, class DestAccessor>
749 bool
750 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
751  DestIterator dul, DestAccessor dest,
752  NoiseNormalizationOptions const & options,
753  VigraTrueType /* isScalar */)
754 {
755  typedef typename SrcAccessor::value_type SrcType;
756  typedef typename DestAccessor::value_type DestType;
757  return noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
758  (sul, slr, src, dul, dest, options);
759 }
760 
761 template <class SrcIterator, class SrcAccessor,
762  class DestIterator, class DestAccessor>
763 bool
764 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
765  DestIterator dul, DestAccessor dest,
766  NoiseNormalizationOptions const & options,
767  VigraFalseType /* isScalar */)
768 {
769  int bands = src.size(sul);
770  for(int b=0; b<bands; ++b)
771  {
772  VectorElementAccessor<SrcAccessor> sband(b, src);
773  VectorElementAccessor<DestAccessor> dband(b, dest);
774  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
775  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
776 
777  if(!noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
778  (sul, slr, sband, dul, dband, options))
779  return false;
780  }
781  return true;
782 }
783 
784 template <class SrcIterator, class SrcAccessor,
785  class DestIterator, class DestAccessor>
786 bool
787 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
788  DestIterator dul, DestAccessor dest,
789  NoiseNormalizationOptions const & options,
790  VigraTrueType /* isScalar */)
791 {
792  typedef typename SrcAccessor::value_type SrcType;
793  typedef typename DestAccessor::value_type DestType;
794  return noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
795  (sul, slr, src, dul, dest, options);
796 }
797 
798 template <class SrcIterator, class SrcAccessor,
799  class DestIterator, class DestAccessor>
800 bool
801 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
802  DestIterator dul, DestAccessor dest,
803  NoiseNormalizationOptions const & options,
804  VigraFalseType /* isScalar */)
805 {
806  int bands = src.size(sul);
807  for(int b=0; b<bands; ++b)
808  {
809  VectorElementAccessor<SrcAccessor> sband(b, src);
810  VectorElementAccessor<DestAccessor> dband(b, dest);
811  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
812  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
813 
814  if(!noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
815  (sul, slr, sband, dul, dband, options))
816  return false;
817  }
818  return true;
819 }
820 
821 template <class SrcIterator, class SrcAccessor,
822  class DestIterator, class DestAccessor>
823 void
824 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
825  DestIterator dul, DestAccessor dest,
826  double a0, double a1, double a2,
827  VigraTrueType /* isScalar */)
828 {
829  ArrayVector<TinyVector<double, 2> > noiseClusters;
830  noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
831  noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1 + a2));
832  noiseClusters.push_back(TinyVector<double, 2>(2.0, a0 + 2.0*a1 + 4.0*a2));
833  transformImage(sul, slr, src, dul, dest,
834  QuadraticNoiseNormalizationFunctor<typename SrcAccessor::value_type,
835  typename DestAccessor::value_type>(noiseClusters));
836 }
837 
838 template <class SrcIterator, class SrcAccessor,
839  class DestIterator, class DestAccessor>
840 void
841 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
842  DestIterator dul, DestAccessor dest,
843  double a0, double a1, double a2,
844  VigraFalseType /* isScalar */)
845 {
846  int bands = src.size(sul);
847  for(int b=0; b<bands; ++b)
848  {
849  VectorElementAccessor<SrcAccessor> sband(b, src);
850  VectorElementAccessor<DestAccessor> dband(b, dest);
851  quadraticNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, a2, VigraTrueType());
852  }
853 }
854 
855 template <class SrcIterator, class SrcAccessor,
856  class DestIterator, class DestAccessor>
857 bool
858 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
859  DestIterator dul, DestAccessor dest,
860  NoiseNormalizationOptions const & options,
861  VigraTrueType /* isScalar */)
862 {
863  typedef typename SrcAccessor::value_type SrcType;
864  typedef typename DestAccessor::value_type DestType;
865  return noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
866  (sul, slr, src, dul, dest, options);
867 }
868 
869 template <class SrcIterator, class SrcAccessor,
870  class DestIterator, class DestAccessor>
871 bool
872 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
873  DestIterator dul, DestAccessor dest,
874  NoiseNormalizationOptions const & options,
875  VigraFalseType /* isScalar */)
876 {
877  int bands = src.size(sul);
878  for(int b=0; b<bands; ++b)
879  {
880  VectorElementAccessor<SrcAccessor> sband(b, src);
881  VectorElementAccessor<DestAccessor> dband(b, dest);
882  typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
883  typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
884 
885  if(!noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
886  (sul, slr, sband, dul, dband, options))
887  return false;
888  }
889  return true;
890 }
891 
892 template <class SrcIterator, class SrcAccessor,
893  class DestIterator, class DestAccessor>
894 void
895 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
896  DestIterator dul, DestAccessor dest,
897  double a0, double a1,
898  VigraTrueType /* isScalar */)
899 {
900  ArrayVector<TinyVector<double, 2> > noiseClusters;
901  noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
902  noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1));
903  transformImage(sul, slr, src, dul, dest,
904  LinearNoiseNormalizationFunctor<typename SrcAccessor::value_type,
905  typename DestAccessor::value_type>(noiseClusters));
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  VigraFalseType /* isScalar */)
915 {
916  int bands = src.size(sul);
917  for(int b=0; b<bands; ++b)
918  {
919  VectorElementAccessor<SrcAccessor> sband(b, src);
920  VectorElementAccessor<DestAccessor> dband(b, dest);
921  linearNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, VigraTrueType());
922  }
923 }
924 
925 } // namespace detail
926 
927 template <bool P>
928 struct noiseVarianceEstimation_can_only_work_on_scalar_images
929 : vigra::staticAssert::AssertBool<P>
930 {};
931 
932 /** \addtogroup NoiseNormalization Noise Normalization
933  Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
934 */
935 //@{
936 
937 /********************************************************/
938 /* */
939 /* noiseVarianceEstimation */
940 /* */
941 /********************************************************/
942 
943 /** \brief Determine the noise variance as a function of the image intensity.
944 
945  This operator applies an algorithm described in
946 
947  W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
948  Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
949  Lecture Notes in Earth Science, Berlin: Springer, 1999
950 
951  in order to estimate the noise variance as a function of the image intensity in a robust way,
952  i.e. so that intensity changes due to edges do not bias the estimate. The source value type
953  (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
954  The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible
955  from two <tt>double</tt> values. The following options can be set via the \a options object
956  (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
957 
958  <tt>useGradient</tt>, <tt>windowRadius</tt>, <tt>noiseEstimationQuantile</tt>, <tt>noiseVarianceInitialGuess</tt>
959 
960  <b> Declarations:</b>
961 
962  pass arguments explicitly:
963  \code
964  namespace vigra {
965  template <class SrcIterator, class SrcAccessor, class BackInsertable>
966  void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
967  BackInsertable & result,
968  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
969  }
970  \endcode
971 
972  use argument objects in conjunction with \ref ArgumentObjectFactories :
973  \code
974  namespace vigra {
975  template <class SrcIterator, class SrcAccessor, class BackInsertable>
976  void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
977  BackInsertable & result,
978  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
979  }
980  \endcode
981 
982  <b> Usage:</b>
983 
984  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
985  Namespace: vigra
986 
987  \code
988  vigra::BImage src(w,h);
989  std::vector<vigra::TinyVector<double, 2> > result;
990 
991  ...
992  vigra::noiseVarianceEstimation(srcImageRange(src), result,
993  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
994 
995  // print the intensity / variance pairs found
996  for(int k=0; k<result.size(); ++k)
997  std::cout << "Intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
998  \endcode
999 
1000  <b> Required Interface:</b>
1001 
1002  \code
1003  SrcIterator upperleft, lowerright;
1004  SrcAccessor src;
1005 
1006  typedef SrcAccessor::value_type SrcType;
1007  typedef NumericTraits<SrcType>::isScalar isScalar;
1008  assert(isScalar::asBool == true);
1009 
1010  double value = src(uperleft);
1011 
1012  BackInsertable result;
1013  typedef BackInsertable::value_type ResultType;
1014  double intensity, variance;
1015  result.push_back(ResultType(intensity, variance));
1016  \endcode
1017 */
1019 
1020 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1021 inline
1022 void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1023  BackInsertable & result,
1024  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1025 {
1026  typedef typename BackInsertable::value_type ResultType;
1027  typedef typename SrcAccessor::value_type SrcType;
1028  typedef typename NumericTraits<SrcType>::isScalar isScalar;
1029 
1030  VIGRA_STATIC_ASSERT((
1031  noiseVarianceEstimation_can_only_work_on_scalar_images<(isScalar::asBool)>));
1032 
1033  detail::noiseVarianceEstimationImpl(sul, slr, src, result, options);
1034 }
1035 
1036 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1037 inline
1038 void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1039  BackInsertable & result,
1040  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1041 {
1042  noiseVarianceEstimation(src.first, src.second, src.third, result, options);
1043 }
1044 
1045 /********************************************************/
1046 /* */
1047 /* noiseVarianceClustering */
1048 /* */
1049 /********************************************************/
1050 
1051 /** \brief Determine the noise variance as a function of the image intensity and cluster the results.
1052 
1053  This operator first calls \ref noiseVarianceEstimation() to obtain a sequence of intensity/variance pairs,
1054  which are then clustered using the median cut algorithm. Then the cluster centers (i.e. average variance vs.
1055  average intensity) are determined and returned in the \a result sequence.
1056 
1057  In addition to the options valid for \ref noiseVarianceEstimation(), the following options can be set via
1058  the \a options object (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
1059 
1060  <tt>clusterCount</tt>, <tt>averagingQuantile</tt>
1061 
1062  <b> Declarations:</b>
1063 
1064  pass arguments explicitly:
1065  \code
1066  namespace vigra {
1067  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1068  void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1069  BackInsertable & result,
1070  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1071  }
1072  \endcode
1073 
1074  use argument objects in conjunction with \ref ArgumentObjectFactories :
1075  \code
1076  namespace vigra {
1077  template <class SrcIterator, class SrcAccessor, class BackInsertable>
1078  void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1079  BackInsertable & result,
1080  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1081  }
1082  \endcode
1083 
1084  <b> Usage:</b>
1085 
1086  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1087  Namespace: vigra
1088 
1089  \code
1090  vigra::BImage src(w,h);
1091  std::vector<vigra::TinyVector<double, 2> > result;
1092 
1093  ...
1094  vigra::noiseVarianceClustering(srcImageRange(src), result,
1095  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1096  clusterCount(15));
1097 
1098  // print the intensity / variance pairs representing the cluster centers
1099  for(int k=0; k<result.size(); ++k)
1100  std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
1101  \endcode
1102 
1103  <b> Required Interface:</b>
1104 
1105  same as \ref noiseVarianceEstimation()
1106 */
1108 
1109 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1110 inline
1111 void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1112  BackInsertable & result,
1113  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1114 {
1115  ArrayVector<TinyVector<double, 2> > variance;
1116  noiseVarianceEstimation(sul, slr, src, variance, options);
1117  detail::noiseVarianceClusteringImpl(variance, result, options.cluster_count, options.averaging_quantile);
1118 }
1119 
1120 template <class SrcIterator, class SrcAccessor, class BackInsertable>
1121 inline
1122 void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1123  BackInsertable & result,
1124  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1125 {
1126  noiseVarianceClustering(src.first, src.second, src.third, result, options);
1127 }
1128 
1129 /********************************************************/
1130 /* */
1131 /* nonparametricNoiseNormalization */
1132 /* */
1133 /********************************************************/
1134 
1135 /** \brief Noise normalization by means of an estimated non-parametric noise model.
1136 
1137  The original image is assumed to be corrupted by noise whose variance depends on the intensity in an unknown way.
1138  The present functions first calls \ref noiseVarianceClustering() to obtain a sequence of intensity/variance pairs
1139  (cluster centers) which estimate this dependency. The cluster centers are connected into a piecewise linear
1140  function which is the inverted according to the formula derived in
1141 
1142  W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
1143  Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
1144  Lecture Notes in Earth Science, Berlin: Springer, 1999
1145 
1146  The inverted formula defines a pixel-wise intensity transformation whose application turns the original image
1147  into one that is corrupted by additive Gaussian noise with unit variance. Most subsequent algorithms will be able
1148  to handle this type of noise much better than the original noise.
1149 
1150  RGB and other multiband images will be processed one band at a time. The function returns <tt>true</tt> on success.
1151  Noise normalization will fail if the original image does not contain sufficiently homogeneous regions that
1152  allow robust estimation of the noise variance.
1153 
1154  The \a options object may use all options described in \ref vigra::NoiseNormalizationOptions.
1155 
1156  <b> Declarations:</b>
1157 
1158  pass arguments explicitly:
1159  \code
1160  namespace vigra {
1161  template <class SrcIterator, class SrcAccessor,
1162  class DestIterator, class DestAccessor>
1163  bool nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1164  DestIterator dul, DestAccessor dest,
1165  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1166  }
1167  \endcode
1168 
1169  use argument objects in conjunction with \ref ArgumentObjectFactories :
1170  \code
1171  namespace vigra {
1172  template <class SrcIterator, class SrcAccessor,
1173  class DestIterator, class DestAccessor>
1174  bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1175  pair<DestIterator, DestAccessor> dest,
1176  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1177  }
1178  \endcode
1179 
1180  <b> Usage:</b>
1181 
1182  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1183  Namespace: vigra
1184 
1185  \code
1186  vigra::BRGBImage src(w,h), dest(w, h);
1187 
1188  ...
1189  vigra::nonparametricNoiseNormalization(srcImageRange(src), destImage(dest),
1190  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1191  clusterCount(15));
1192  \endcode
1193 
1194  <b> Required Interface:</b>
1195 
1196  same as \ref noiseVarianceEstimation()
1197 */
1199 
1200 template <class SrcIterator, class SrcAccessor,
1201  class DestIterator, class DestAccessor>
1202 inline bool
1203 nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1204  DestIterator dul, DestAccessor dest,
1205  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1206 {
1207  typedef typename SrcAccessor::value_type SrcType;
1208 
1209  return detail::nonparametricNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1210  typename NumericTraits<SrcType>::isScalar());
1211 }
1212 
1213 template <class SrcIterator, class SrcAccessor,
1214  class DestIterator, class DestAccessor>
1215 inline
1216 bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1217  pair<DestIterator, DestAccessor> dest,
1218  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1219 {
1220  return nonparametricNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1221 }
1222 
1223 /********************************************************/
1224 /* */
1225 /* quadraticNoiseNormalization */
1226 /* */
1227 /********************************************************/
1228 
1229 /** \brief Noise normalization by means of an estimated quadratic noise model.
1230 
1231  This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
1232  model for the dependency between intensity and noise variance: it assumes that this dependency is a
1233  quadratic function rather than a piecewise linear function. If the quadratic model is applicable, it leads
1234  to a somewhat smoother transformation.
1235 
1236  <b> Declarations:</b>
1237 
1238  pass arguments explicitly:
1239  \code
1240  namespace vigra {
1241  template <class SrcIterator, class SrcAccessor,
1242  class DestIterator, class DestAccessor>
1243  bool quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1244  DestIterator dul, DestAccessor dest,
1245  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1246  }
1247  \endcode
1248 
1249  use argument objects in conjunction with \ref ArgumentObjectFactories :
1250  \code
1251  namespace vigra {
1252  template <class SrcIterator, class SrcAccessor,
1253  class DestIterator, class DestAccessor>
1254  bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1255  pair<DestIterator, DestAccessor> dest,
1256  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1257  }
1258  \endcode
1259 
1260  <b> Usage:</b>
1261 
1262  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1263  Namespace: vigra
1264 
1265  \code
1266  vigra::BRGBImage src(w,h), dest(w, h);
1267 
1268  ...
1269  vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
1270  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1271  clusterCount(15));
1272  \endcode
1273 
1274  <b> Required Interface:</b>
1275 
1276  same as \ref noiseVarianceEstimation()
1277 */
1279 
1280 template <class SrcIterator, class SrcAccessor,
1281  class DestIterator, class DestAccessor>
1282 inline bool
1283 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1284  DestIterator dul, DestAccessor dest,
1285  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1286 {
1287  typedef typename SrcAccessor::value_type SrcType;
1288 
1289  return detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1290  typename NumericTraits<SrcType>::isScalar());
1291 }
1292 
1293 template <class SrcIterator, class SrcAccessor,
1294  class DestIterator, class DestAccessor>
1295 inline
1296 bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1297  pair<DestIterator, DestAccessor> dest,
1298  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1299 {
1300  return quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1301 }
1302 
1303 /********************************************************/
1304 /* */
1305 /* quadraticNoiseNormalization */
1306 /* */
1307 /********************************************************/
1308 
1309 /** \brief Noise normalization by means of a given quadratic noise model.
1310 
1311  This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
1312  functional dependency of the noise variance from the intensity is given (by a quadratic function)
1313  rather than estimated:
1314 
1315  \code
1316  variance = a0 + a1 * intensity + a2 * sq(intensity)
1317  \endcode
1318 
1319  <b> Declarations:</b>
1320 
1321  pass arguments explicitly:
1322  \code
1323  namespace vigra {
1324  template <class SrcIterator, class SrcAccessor,
1325  class DestIterator, class DestAccessor>
1326  void quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1327  DestIterator dul, DestAccessor dest,
1328  double a0, double a1, double a2);
1329  }
1330  \endcode
1331 
1332  use argument objects in conjunction with \ref ArgumentObjectFactories :
1333  \code
1334  namespace vigra {
1335  template <class SrcIterator, class SrcAccessor,
1336  class DestIterator, class DestAccessor>
1337  void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1338  pair<DestIterator, DestAccessor> dest,
1339  double a0, double a1, double a2);
1340  }
1341  \endcode
1342 
1343  <b> Usage:</b>
1344 
1345  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1346  Namespace: vigra
1347 
1348  \code
1349  vigra::BRGBImage src(w,h), dest(w, h);
1350 
1351  ...
1352  vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
1353  100, 0.02, 1e-6);
1354  \endcode
1355 
1356  <b> Required Interface:</b>
1357 
1358  The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
1359  are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
1360  or a vector whose elements are assignable from <tt>double</tt>.
1361 */
1363 
1364 template <class SrcIterator, class SrcAccessor,
1365  class DestIterator, class DestAccessor>
1366 inline void
1367 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1368  DestIterator dul, DestAccessor dest,
1369  double a0, double a1, double a2)
1370 {
1371  typedef typename SrcAccessor::value_type SrcType;
1372 
1373  detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1, a2,
1374  typename NumericTraits<SrcType>::isScalar());
1375 }
1376 
1377 template <class SrcIterator, class SrcAccessor,
1378  class DestIterator, class DestAccessor>
1379 inline
1380 void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1381  pair<DestIterator, DestAccessor> dest,
1382  double a0, double a1, double a2)
1383 {
1384  quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1, a2);
1385 }
1386 
1387 /********************************************************/
1388 /* */
1389 /* linearNoiseNormalization */
1390 /* */
1391 /********************************************************/
1392 
1393 /** \brief Noise normalization by means of an estimated linear noise model.
1394 
1395  This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
1396  model for the dependency between intensity and noise variance: it assumes that this dependency is a
1397  linear function rather than a piecewise linear function. If the linear model is applicable, it leads
1398  to a very simple transformation which is similar to the familiar gamma correction.
1399 
1400  <b> Declarations:</b>
1401 
1402  pass arguments explicitly:
1403  \code
1404  namespace vigra {
1405  template <class SrcIterator, class SrcAccessor,
1406  class DestIterator, class DestAccessor>
1407  bool linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1408  DestIterator dul, DestAccessor dest,
1409  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1410  }
1411  \endcode
1412 
1413  use argument objects in conjunction with \ref ArgumentObjectFactories :
1414  \code
1415  namespace vigra {
1416  template <class SrcIterator, class SrcAccessor,
1417  class DestIterator, class DestAccessor>
1418  bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1419  pair<DestIterator, DestAccessor> dest,
1420  NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
1421  }
1422  \endcode
1423 
1424  <b> Usage:</b>
1425 
1426  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1427  Namespace: vigra
1428 
1429  \code
1430  vigra::BRGBImage src(w,h), dest(w, h);
1431 
1432  ...
1433  vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
1434  vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
1435  clusterCount(15));
1436  \endcode
1437 
1438  <b> Required Interface:</b>
1439 
1440  same as \ref noiseVarianceEstimation()
1441 */
1443 
1444 template <class SrcIterator, class SrcAccessor,
1445  class DestIterator, class DestAccessor>
1446 inline bool
1447 linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1448  DestIterator dul, DestAccessor dest,
1449  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1450 {
1451  typedef typename SrcAccessor::value_type SrcType;
1452 
1453  return detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
1454  typename NumericTraits<SrcType>::isScalar());
1455 }
1456 
1457 template <class SrcIterator, class SrcAccessor,
1458  class DestIterator, class DestAccessor>
1459 inline
1460 bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1461  pair<DestIterator, DestAccessor> dest,
1462  NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
1463 {
1464  return linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
1465 }
1466 
1467 /********************************************************/
1468 /* */
1469 /* linearNoiseNormalization */
1470 /* */
1471 /********************************************************/
1472 
1473 /** \brief Noise normalization by means of a given linear noise model.
1474 
1475  This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
1476  functional dependency of the noise variance from the intensity is given (as a linear function)
1477  rather than estimated:
1478 
1479  \code
1480  variance = a0 + a1 * intensity
1481  \endcode
1482 
1483  <b> Declarations:</b>
1484 
1485  pass arguments explicitly:
1486  \code
1487  namespace vigra {
1488  template <class SrcIterator, class SrcAccessor,
1489  class DestIterator, class DestAccessor>
1490  void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1491  DestIterator dul, DestAccessor dest,
1492  double a0, double a1);
1493  }
1494  \endcode
1495 
1496  use argument objects in conjunction with \ref ArgumentObjectFactories :
1497  \code
1498  namespace vigra {
1499  template <class SrcIterator, class SrcAccessor,
1500  class DestIterator, class DestAccessor>
1501  void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1502  pair<DestIterator, DestAccessor> dest,
1503  double a0, double a1);
1504  }
1505  \endcode
1506 
1507  <b> Usage:</b>
1508 
1509  <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
1510  Namespace: vigra
1511 
1512  \code
1513  vigra::BRGBImage src(w,h), dest(w, h);
1514 
1515  ...
1516  vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
1517  100, 0.02);
1518  \endcode
1519 
1520  <b> Required Interface:</b>
1521 
1522  The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
1523  are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
1524  or a vector whose elements are assignable from <tt>double</tt>.
1525 */
1527 
1528 template <class SrcIterator, class SrcAccessor,
1529  class DestIterator, class DestAccessor>
1530 inline
1531 void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
1532  DestIterator dul, DestAccessor dest,
1533  double a0, double a1)
1534 {
1535  typedef typename SrcAccessor::value_type SrcType;
1536 
1537  detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1,
1538  typename NumericTraits<SrcType>::isScalar());
1539 }
1540 
1541 template <class SrcIterator, class SrcAccessor,
1542  class DestIterator, class DestAccessor>
1543 inline
1544 void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
1545  pair<DestIterator, DestAccessor> dest,
1546  double a0, double a1)
1547 {
1548  linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1);
1549 }
1550 
1551 //@}
1552 
1553 } // namespace vigra
1554 
1555 #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 (Thu Jun 14 2012)