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