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

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

html generated using doxygen and Python
vigra 1.11.1