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

| vigra/noise_normalization.hxx |
--- | --- | ---

    00001 /************************************************************************/
    00002 /*                                                                      */
    00003 /*               Copyright 1998-2006 by Ullrich Koethe                  */
    00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
    00005 /*                                                                      */
    00006 /*    This file is part of the VIGRA computer vision library.           */
    00007 /*    ( Version 1.5.0, Dec 07 2006 )                                    */
    00008 /*    The VIGRA Website is                                              */
    00009 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
    00010 /*    Please direct questions, bug reports, and contributions to        */
    00011 /*        koethe@informatik.uni-hamburg.de          or                  */
    00012 /*        vigra@kogs1.informatik.uni-hamburg.de                         */
    00013 /*                                                                      */
    00014 /*    Permission is hereby granted, free of charge, to any person       */
    00015 /*    obtaining a copy of this software and associated documentation    */
    00016 /*    files (the "Software"), to deal in the Software without           */
    00017 /*    restriction, including without limitation the rights to use,      */
    00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
    00019 /*    sell copies of the Software, and to permit persons to whom the    */
    00020 /*    Software is furnished to do so, subject to the following          */
    00021 /*    conditions:                                                       */
    00022 /*                                                                      */
    00023 /*    The above copyright notice and this permission notice shall be    */
    00024 /*    included in all copies or substantial portions of the             */
    00025 /*    Software.                                                         */
    00026 /*                                                                      */
    00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
    00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
    00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
    00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
    00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
    00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
    00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
    00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
    00035 /*                                                                      */
    00036 /************************************************************************/
    00037
    00038
    00039 #ifndef VIGRA_NOISE_NORMALIZATION_HXX
    00040 #define VIGRA_NOISE_NORMALIZATION_HXX
    00041
    00042 #include "utilities.hxx"
    00043 #include "tinyvector.hxx"
    00044 #include "stdimage.hxx"
    00045 #include "transformimage.hxx"
    00046 #include "combineimages.hxx"
    00047 #include "localminmax.hxx"
    00048 #include "functorexpression.hxx"
    00049 #include "numerictraits.hxx"
    00050 #include "separableconvolution.hxx"
    00051 #include "linear_solve.hxx"
    00052 #include "array_vector.hxx"
    00053 #include "static_assert.hxx"
    00054 #include <algorithm>
    00055
    00056 namespace vigra {
    00057
    00058 /** \addtogroup NoiseNormalization Noise Normalization
    00059     Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
    00060 */
    00061 //@{
    00062
    00063 /********************************************************/
    00064 /*                                                      */
    00065 /*               NoiseNormalizationOptions              */
    00066 /*                                                      */
    00067 /********************************************************/
    00068
    00069 /** \brief Pass options to one of the noise normalization functions.
    00070
    00071     <tt>NoiseNormalizationOptions</tt>  is an argument object that holds various optional
    00072     parameters used by the noise normalization functions. If a parameter is not explicitly
    00073     set, a suitable default will be used.
    00074
    00075     <b> Usage:</b>
    00076
    00077         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    00078     Namespace: vigra
    00079
    00080     \code
    00081     vigra::BImage src(w,h);
    00082     std::vector<vigra::TinyVector<double, 2> > result;
    00083
    00084     ...
    00085     vigra::noiseVarianceEstimation(srcImageRange(src), result,
    00086                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
    00087     \endcode
    00088 */
    00089 class NoiseNormalizationOptions
    00090 {
    00091   public:
    00092
    00093         /** Initialize all options with default values.
    00094         */
    00095     NoiseNormalizationOptions()
    00096     : window_radius(6),
    00097       cluster_count(10),
    00098       noise_estimation_quantile(1.5),
    00099       averaging_quantile(0.8),
    00100       noise_variance_initial_guess(10.0),
    00101       use_gradient(true)
    00102     {}
    00103
    00104         /** Select the noise estimation algorithm.
    00105
    00106             If \a r is <tt>true</tt>, use the gradient-based noise estimator according to F&ouml;rstner (default).
    00107             Otherwise, use an algorithm that uses the intensity values directly.
    00108         */
    00109     NoiseNormalizationOptions & useGradient(bool r)
    00110     {
    00111         use_gradient = r;
    00112         return *this;
    00113     }
    00114
    00115         /** Set the window radius for a single noise estimate.
    00116             Every window of the given size gives raise to one intensity/variance pair.<br>
    00117             Default: 6 pixels
    00118         */
    00119     NoiseNormalizationOptions & windowRadius(unsigned int r)
    00120     {
    00121         vigra_precondition(r > 0,
    00122             "NoiseNormalizationOptions: window radius must be > 0.");
    00123         window_radius = r;
    00124         return *this;
    00125     }
    00126
    00127         /** Set the number of clusters for non-parametric noise normalization.
    00128             The intensity/variance pairs found are grouped into clusters before the noise
    00129             normalization transform is computed.<br>
    00130             Default: 10 clusters
    00131         */
    00132     NoiseNormalizationOptions & clusterCount(unsigned int c)
    00133     {
    00134         vigra_precondition(c > 0,
    00135             "NoiseNormalizationOptions: cluster count must be > 0.");
    00136         cluster_count = c;
    00137         return *this;
    00138     }
    00139
    00140         /** Set the quantile for cluster averaging.
    00141             After clustering, the cluster center (i.e. average noise variance as a function of the average
    00142             intensity in the cluster) is computed using only the cluster members whose estimated variance
    00143             is below \a quantile times the maximum variance in the cluster.<br>
    00144             Default: 0.8<br>
    00145             Precondition: 0 < \a quantile <= 1.0
    00146         */
    00147     NoiseNormalizationOptions & averagingQuantile(double quantile)
    00148     {
    00149         vigra_precondition(quantile > 0.0 && quantile <= 1.0,
    00150             "NoiseNormalizationOptions: averaging quantile must be between 0 and 1.");
    00151         averaging_quantile = quantile;
    00152         return *this;
    00153     }
    00154
    00155         /** Set the operating range of the robust noise estimator.
    00156             Intensity changes that are larger than \a quantile times the current estimate of the noise variance
    00157             are ignored by the robust noise estimator.<br>
    00158             Default: 1.5<br>
    00159             Precondition: 0 < \a quantile
    00160         */
    00161     NoiseNormalizationOptions & noiseEstimationQuantile(double quantile)
    00162     {
    00163         vigra_precondition(quantile > 0.0,
    00164             "NoiseNormalizationOptions: noise estimation quantile must be > 0.");
    00165         noise_estimation_quantile = quantile;
    00166         return *this;
    00167     }
    00168
    00169         /** Set the initial estimate of the noise variance.
    00170             Robust noise variance estimation is an iterative procedure starting at the given value.<br>
    00171             Default: 10.0<br>
    00172             Precondition: 0 < \a quess
    00173         */
    00174     NoiseNormalizationOptions & noiseVarianceInitialGuess(double guess)
    00175     {
    00176         vigra_precondition(guess > 0.0,
    00177             "NoiseNormalizationOptions: noise variance initial guess must be > 0.");
    00178         noise_variance_initial_guess = guess;
    00179         return *this;
    00180     }
    00181
    00182     unsigned int window_radius, cluster_count;
    00183     double noise_estimation_quantile, averaging_quantile, noise_variance_initial_guess;
    00184     bool use_gradient;
    00185 };
    00186
    00187 //@}
    00188
    00189 template <class ArgumentType, class ResultType>
    00190 class NonparametricNoiseNormalizationFunctor
    00191 {
    00192     struct Segment
    00193     {
    00194         double lower, a, b, shift;
    00195     };
    00196
    00197     ArrayVector<Segment> segments_;
    00198
    00199     template <class T>
    00200     double exec(unsigned int k, T t) const
    00201     {
    00202         if(segments_[k].a == 0.0)
    00203         {
    00204             return t / VIGRA_CSTD::sqrt(segments_[k].b);
    00205         }
    00206         else
    00207         {
    00208             return 2.0 / segments_[k].a * VIGRA_CSTD::sqrt(std::max(0.0, segments_[k].a * t + segments_[k].b));
    00209         }
    00210     }
    00211
    00212   public:
    00213     typedef ArgumentType argument_type;
    00214     typedef ResultType result_type;
    00215
    00216     template <class Vector>
    00217     NonparametricNoiseNormalizationFunctor(Vector const & clusters)
    00218     : segments_(clusters.size()-1)
    00219     {
    00220         for(unsigned int k = 0; k<segments_.size(); ++k)
    00221         {
    00222             segments_[k].lower = clusters[k][0];
    00223             segments_[k].a = (clusters[k+1][1] - clusters[k][1]) / (clusters[k+1][0] - clusters[k][0]);
    00224             segments_[k].b = clusters[k][1] - segments_[k].a * clusters[k][0];
    00225             // FIXME: set a to zero if it is very small
    00226             //          - determine what 'very small' means
    00227             //          - shouldn't the two formulas (for a == 0, a != 0) be equal in the limit a -> 0 ?
    00228
    00229             if(k == 0)
    00230             {
    00231                 segments_[k].shift = segments_[k].lower - exec(k, segments_[k].lower);
    00232             }
    00233             else
    00234             {
    00235                 segments_[k].shift = exec(k-1, segments_[k].lower) - exec(k, segments_[k].lower) + segments_[k-1].shift;
    00236             }
    00237         }
    00238     }
    00239
    00240     result_type operator()(argument_type t) const
    00241     {
    00242         // find the segment
    00243         unsigned int k = 0;
    00244         for(; k < segments_.size(); ++k)
    00245             if(t < segments_[k].lower)
    00246                 break;
    00247         if(k > 0)
    00248             --k;
    00249         return detail::RequiresExplicitCast<ResultType>::cast(exec(k, t) + segments_[k].shift);
    00250     }
    00251 };
    00252
    00253 template <class ArgumentType, class ResultType>
    00254 class QuadraticNoiseNormalizationFunctor
    00255 {
    00256     double a, b, c, d, f, o;
    00257
    00258     void init(double ia, double ib, double ic, double xmin)
    00259     {
    00260         a = ia;
    00261         b = ib;
    00262         c = ic;
    00263         d = VIGRA_CSTD::sqrt(VIGRA_CSTD::fabs(c));
    00264         if(c > 0.0)
    00265         {
    00266             o = VIGRA_CSTD::log(VIGRA_CSTD::fabs((2.0*c*xmin + b)/d + 2*VIGRA_CSTD::sqrt(c*sq(xmin) +b*xmin + a)))/d;
    00267             f = 0.0;
    00268         }
    00269         else
    00270         {
    00271             f = VIGRA_CSTD::sqrt(b*b - 4.0*a*c);
    00272             o = -VIGRA_CSTD::asin((2.0*c*xmin+b)/f)/d;
    00273         }
    00274     }
    00275
    00276   public:
    00277     typedef ArgumentType argument_type;
    00278     typedef ResultType result_type;
    00279
    00280     template <class Vector>
    00281     QuadraticNoiseNormalizationFunctor(Vector const & clusters)
    00282     {
    00283         double xmin = NumericTraits<double>::max();
    00284         Matrix<double> m(3,3), r(3, 1), l(3, 1);
    00285         for(unsigned int k = 0; k<clusters.size(); ++k)
    00286         {
    00287             l(0,0) = 1.0;
    00288             l(1,0) = clusters[k][0];
    00289             l(2,0) = sq(clusters[k][0]);
    00290             m += outer(l);
    00291             r += clusters[k][1]*l;
    00292             if(clusters[k][0] < xmin)
    00293                 xmin = clusters[k][0];
    00294         }
    00295
    00296         linearSolve(m, r, l);
    00297         init(l(0,0), l(1,0), l(2,0), xmin);
    00298     }
    00299
    00300     result_type operator()(argument_type t) const
    00301     {
    00302         double r;
    00303         if(c > 0.0)
    00304             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;
    00305         else
    00306             r = -VIGRA_CSTD::asin((2.0*c*t+b)/f)/d-o;
    00307         return detail::RequiresExplicitCast<ResultType>::cast(r);
    00308     }
    00309 };
    00310
    00311 template <class ArgumentType, class ResultType>
    00312 class LinearNoiseNormalizationFunctor
    00313 {
    00314     double a, b, o;
    00315
    00316     void init(double ia, double ib, double xmin)
    00317     {
    00318         a = ia;
    00319         b = ib;
    00320         if(b != 0.0)
    00321         {
    00322             o = xmin - 2.0 / b * VIGRA_CSTD::sqrt(a + b * xmin);
    00323         }
    00324         else
    00325         {
    00326             o = xmin - xmin / VIGRA_CSTD::sqrt(a);
    00327         }
    00328     }
    00329
    00330   public:
    00331     typedef ArgumentType argument_type;
    00332     typedef ResultType result_type;
    00333
    00334     template <class Vector>
    00335     LinearNoiseNormalizationFunctor(Vector const & clusters)
    00336     {
    00337         double xmin = NumericTraits<double>::max();
    00338         Matrix<double> m(2,2), r(2, 1), l(2, 1);
    00339         for(unsigned int k = 0; k<clusters.size(); ++k)
    00340         {
    00341             l(0,0) = 1.0;
    00342             l(1,0) = clusters[k][0];
    00343             m += outer(l);
    00344             r += clusters[k][1]*l;
    00345             if(clusters[k][0] < xmin)
    00346                 xmin = clusters[k][0];
    00347         }
    00348
    00349         linearSolve(m, r, l);
    00350         init(l(0,0), l(1,0), xmin);
    00351     }
    00352
    00353     result_type operator()(argument_type t) const
    00354     {
    00355         double r;
    00356         if(b != 0.0)
    00357             r = 2.0 / b * VIGRA_CSTD::sqrt(a + b*t) + o;
    00358         else
    00359             r =  t / VIGRA_CSTD::sqrt(a) + o;
    00360         return detail::RequiresExplicitCast<ResultType>::cast(r);
    00361     }
    00362 };
    00363
    00364 #define VIGRA_NoiseNormalizationFunctor(name, type, size) \
    00365 template <class ResultType> \
    00366 class name<type, ResultType> \
    00367 { \
    00368     ResultType lut_[size]; \
    00369     \
    00370   public: \
    00371     typedef type argument_type; \
    00372     typedef ResultType result_type; \
    00373      \
    00374     template <class Vector> \
    00375     name(Vector const & clusters) \
    00376     { \
    00377         name<double, ResultType> f(clusters); \
    00378          \
    00379         for(unsigned int k = 0; k < size; ++k) \
    00380         { \
    00381             lut_[k] = f(k); \
    00382         } \
    00383     } \
    00384      \
    00385     result_type operator()(argument_type t) const \
    00386     { \
    00387         return lut_[t]; \
    00388     } \
    00389 };
    00390
    00391 VIGRA_NoiseNormalizationFunctor(NonparametricNoiseNormalizationFunctor, UInt8, 256)
    00392 VIGRA_NoiseNormalizationFunctor(NonparametricNoiseNormalizationFunctor, UInt16, 65536)
    00393 VIGRA_NoiseNormalizationFunctor(QuadraticNoiseNormalizationFunctor, UInt8, 256)
    00394 VIGRA_NoiseNormalizationFunctor(QuadraticNoiseNormalizationFunctor, UInt16, 65536)
    00395 VIGRA_NoiseNormalizationFunctor(LinearNoiseNormalizationFunctor, UInt8, 256)
    00396 VIGRA_NoiseNormalizationFunctor(LinearNoiseNormalizationFunctor, UInt16, 65536)
    00397
    00398 #undef VIGRA_NoiseNormalizationFunctor
    00399
    00400 namespace detail {
    00401
    00402 template <class SrcIterator, class SrcAcessor,
    00403           class GradIterator>
    00404 bool
    00405 iterativeNoiseEstimationChi2(SrcIterator s, SrcAcessor src, GradIterator g,
    00406                          double & mean, double & variance,
    00407                          double robustnessThreshold, int windowRadius)
    00408 {
    00409     double l2 = sq(robustnessThreshold);
    00410     double countThreshold = 1.0 - VIGRA_CSTD::exp(-l2);
    00411     double f = (1.0 - VIGRA_CSTD::exp(-l2)) / (1.0 - (1.0 + l2)*VIGRA_CSTD::exp(-l2));
    00412
    00413     Diff2D ul(-windowRadius, -windowRadius);
    00414     int r2 = sq(windowRadius);
    00415
    00416     for(int iter=0; iter<100 ; ++iter) // maximum iteration 100 only for terminating
    00417                                        // if something is wrong
    00418     {
    00419         double sum=0.0;
    00420         double gsum=0.0;
    00421         unsigned int count = 0;
    00422         unsigned int tcount = 0;
    00423
    00424         SrcIterator siy = s + ul;
    00425         GradIterator giy = g + ul;
    00426         for(int y=-windowRadius; y <= windowRadius; y++, ++siy.y, ++giy.y)
    00427         {
    00428             typename SrcIterator::row_iterator six = siy.rowIterator();
    00429             typename GradIterator::row_iterator gix = giy.rowIterator();
    00430             for(int x=-windowRadius; x <= windowRadius; x++, ++six, ++gix)
    00431             {
    00432                 if (sq(x) + sq(y) > r2)
    00433                     continue;
    00434
    00435                 ++tcount;
    00436                 if (*gix < l2*variance)
    00437                 {
    00438                     sum += src(six);
    00439                     gsum += *gix;
    00440                     ++count;
    00441                 }
    00442             }
    00443         }
    00444         if (count==0) // not homogeneous enough
    00445             return false;
    00446
    00447         double oldvariance = variance;
    00448         variance= f * gsum / count;
    00449         mean = sum / count;
    00450
    00451         if ( closeAtTolerance(oldvariance - variance, 0.0, 1e-10))
    00452             return (count >= tcount * countThreshold / 2.0); // sufficiently many valid points
    00453     }
    00454     return false; // no convergence
    00455 }
    00456
    00457 template <class SrcIterator, class SrcAcessor,
    00458           class GradIterator>
    00459 bool
    00460 iterativeNoiseEstimationGauss(SrcIterator s, SrcAcessor src, GradIterator,
    00461                          double & mean, double & variance,
    00462                          double robustnessThreshold, int windowRadius)
    00463 {
    00464     double l2 = sq(robustnessThreshold);
    00465     double countThreshold = erf(VIGRA_CSTD::sqrt(0.5 * l2));
    00466     double f = countThreshold / (countThreshold - VIGRA_CSTD::sqrt(2.0/M_PI*l2)*VIGRA_CSTD::exp(-l2/2.0));
    00467
    00468     mean = src(s);
    00469
    00470     Diff2D ul(-windowRadius, -windowRadius);
    00471     int r2 = sq(windowRadius);
    00472
    00473     for(int iter=0; iter<100 ; ++iter) // maximum iteration 100 only for terminating
    00474                                        // if something is wrong
    00475     {
    00476         double sum = 0.0;
    00477         double sum2 = 0.0;
    00478         unsigned int count = 0;
    00479         unsigned int tcount = 0;
    00480
    00481         SrcIterator siy = s + ul;
    00482         for(int y=-windowRadius; y <= windowRadius; y++, ++siy.y)
    00483         {
    00484             typename SrcIterator::row_iterator six = siy.rowIterator();
    00485             for(int x=-windowRadius; x <= windowRadius; x++, ++six)
    00486             {
    00487                 if (sq(x) + sq(y) > r2)
    00488                     continue;
    00489
    00490                 ++tcount;
    00491                 if (sq(src(six) - mean) < l2*variance)
    00492                 {
    00493                     sum += src(six);
    00494                     sum2 += sq(src(six));
    00495                     ++count;
    00496                 }
    00497             }
    00498         }
    00499         if (count==0) // not homogeneous enough
    00500             return false;
    00501
    00502         double oldmean = mean;
    00503         double oldvariance = variance;
    00504         mean = sum / count;
    00505         variance= f * (sum2 / count - sq(mean));
    00506
    00507         if ( closeAtTolerance(oldmean - mean, 0.0, 1e-10) &&
    00508              closeAtTolerance(oldvariance - variance, 0.0, 1e-10))
    00509             return (count >= tcount * countThreshold / 2.0); // sufficiently many valid points
    00510     }
    00511     return false; // no convergence
    00512 }
    00513
    00514
    00515 template <class SrcIterator, class SrcAccessor,
    00516           class DestIterator, class DestAccessor>
    00517 void
    00518 symmetricDifferenceSquaredMagnitude(
    00519      SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00520      DestIterator dul, DestAccessor dest)
    00521 {
    00522     using namespace functor;
    00523     int w = slr.x - sul.x;
    00524     int h = slr.y - sul.y;
    00525
    00526     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00527     typedef BasicImage<TmpType> TmpImage;
    00528
    00529     Kernel1D<double> mask;
    00530     mask.initSymmetricGradient();
    00531     mask.setBorderTreatment(BORDER_TREATMENT_REFLECT);
    00532
    00533     TmpImage dx(w, h), dy(w, h);
    00534     separableConvolveX(srcIterRange(sul, slr, src), destImage(dx),  kernel1d(mask));
    00535     separableConvolveY(srcIterRange(sul, slr, src), destImage(dy),  kernel1d(mask));
    00536     combineTwoImages(srcImageRange(dx), srcImage(dy), destIter(dul, dest), Arg1()*Arg1() + Arg2()*Arg2());
    00537 }
    00538
    00539 template <class SrcIterator, class SrcAccessor,
    00540           class DestIterator, class DestAccessor>
    00541 void
    00542 findHomogeneousRegionsFoerstner(
    00543      SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00544      DestIterator dul, DestAccessor dest,
    00545      unsigned int windowRadius = 6, double homogeneityThreshold = 40.0)
    00546 {
    00547     using namespace vigra::functor;
    00548     int w = slr.x - sul.x;
    00549     int h = slr.y - sul.y;
    00550
    00551     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00552     typedef BasicImage<TmpType> TmpImage;
    00553
    00554     BImage btmp(w, h);
    00555     transformImage(srcIterRange(sul, slr, src), destImage(btmp),
    00556                     ifThenElse(Arg1() <= Param(homogeneityThreshold), Param(1), Param(0)));
    00557
    00558     // Erosion
    00559     discErosion(srcImageRange(btmp), destIter(dul, dest), windowRadius);
    00560 }
    00561
    00562 template <class SrcIterator, class SrcAccessor,
    00563           class DestIterator, class DestAccessor>
    00564 void
    00565 findHomogeneousRegions(
    00566      SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00567      DestIterator dul, DestAccessor dest)
    00568 {
    00569     localMinima(sul, slr, src, dul, dest);
    00570 }
    00571
    00572 template <class Vector1, class Vector2>
    00573 void noiseVarianceListMedianCut(Vector1 const & noise, Vector2 & clusters,
    00574                                 unsigned int maxClusterCount)
    00575 {
    00576     typedef typename Vector2::value_type Result;
    00577
    00578     clusters.push_back(Result(0, noise.size()));
    00579
    00580     while(clusters.size() <= maxClusterCount)
    00581     {
    00582         // find biggest cluster
    00583         unsigned int kMax;
    00584         double diffMax = 0.0;
    00585         for(unsigned int k=0; k < clusters.size(); ++k)
    00586         {
    00587             double diff = noise[clusters[k][1]-1][0] - noise[clusters[k][0]][0];
    00588             if(diff > diffMax)
    00589             {
    00590                 diffMax = diff;
    00591                 kMax = k;
    00592             }
    00593         }
    00594
    00595         if(diffMax == 0.0)
    00596             return; // all clusters have only one value
    00597
    00598         unsigned int k1 = clusters[kMax][0],
    00599                      k2 = clusters[kMax][1];
    00600         unsigned int kSplit = k1 + (k2 - k1) / 2;
    00601         clusters[kMax][1] = kSplit;
    00602         clusters.push_back(Result(kSplit, k2));
    00603     }
    00604 }
    00605
    00606 struct SortNoiseByMean
    00607 {
    00608     template <class T>
    00609     bool operator()(T const & l, T const & r) const
    00610     {
    00611         return l[0] < r[0];
    00612     }
    00613 };
    00614
    00615 struct SortNoiseByVariance
    00616 {
    00617     template <class T>
    00618     bool operator()(T const & l, T const & r) const
    00619     {
    00620         return l[1] < r[1];
    00621     }
    00622 };
    00623
    00624 template <class Vector1, class Vector2, class Vector3>
    00625 void noiseVarianceClusterAveraging(Vector1 & noise, Vector2 & clusters,
    00626                                    Vector3 & result, double quantile)
    00627 {
    00628     typedef typename Vector1::iterator Iter;
    00629     typedef typename Vector3::value_type Result;
    00630
    00631     for(unsigned int k=0; k<clusters.size(); ++k)
    00632     {
    00633         Iter i1 = noise.begin() + clusters[k][0];
    00634         Iter i2 = noise.begin() + clusters[k][1];
    00635
    00636         std::sort(i1, i2, SortNoiseByVariance());
    00637
    00638         std::size_t size = static_cast<std::size_t>(VIGRA_CSTD::ceil(quantile*(i2 - i1)));
    00639         if(static_cast<std::size_t>(i2 - i1) < size)
    00640             size = i2 - i1;
    00641         if(size < 1)
    00642             size = 1;
    00643         i2 = i1 + size;
    00644
    00645         double mean = 0.0,
    00646                variance = 0.0;
    00647         for(; i1 < i2; ++i1)
    00648         {
    00649             mean += (*i1)[0];
    00650             variance += (*i1)[1];
    00651         }
    00652
    00653         result.push_back(Result(mean / size, variance / size));
    00654     }
    00655 }
    00656
    00657 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00658 void noiseVarianceEstimationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00659                            BackInsertable & result,
    00660                            NoiseNormalizationOptions const & options)
    00661 {
    00662     typedef typename BackInsertable::value_type ResultType;
    00663
    00664     unsigned int w = slr.x - sul.x;
    00665     unsigned int h = slr.y - sul.y;
    00666
    00667     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00668     typedef BasicImage<TmpType> TmpImage;
    00669
    00670     TmpImage gradient(w, h);
    00671     symmetricDifferenceSquaredMagnitude(sul, slr, src, gradient.upperLeft(), gradient.accessor());
    00672
    00673     BImage homogeneous(w, h);
    00674     findHomogeneousRegions(gradient.upperLeft(), gradient.lowerRight(), gradient.accessor(),
    00675                                    homogeneous.upperLeft(), homogeneous.accessor());
    00676
    00677     // Generate noise of each of the remaining pixels == centers of homogenous areas (border is not used)
    00678     unsigned int windowRadius = options.window_radius;
    00679     for(unsigned int y=windowRadius; y<h-windowRadius; ++y)
    00680     {
    00681         for(unsigned int x=windowRadius; x<w-windowRadius; ++x)
    00682         {
    00683             if (! homogeneous(x, y))
    00684                 continue;
    00685
    00686             Diff2D center(x, y);
    00687             double mean = 0.0, variance = options.noise_variance_initial_guess;
    00688
    00689             bool success;
    00690
    00691             if(options.use_gradient)
    00692             {
    00693                 success = iterativeNoiseEstimationChi2(sul + center, src,
    00694                               gradient.upperLeft() + center, mean, variance,
    00695                               options.noise_estimation_quantile, windowRadius);
    00696             }
    00697             else
    00698             {
    00699                 success = iterativeNoiseEstimationGauss(sul + center, src,
    00700                               gradient.upperLeft() + center, mean, variance,
    00701                               options.noise_estimation_quantile, windowRadius);
    00702             }
    00703             if (success)
    00704             {
    00705                 result.push_back(ResultType(mean, variance));
    00706             }
    00707         }
    00708     }
    00709 }
    00710
    00711 template <class Vector, class BackInsertable>
    00712 void noiseVarianceClusteringImpl(Vector & noise, BackInsertable & result,
    00713                            unsigned int clusterCount, double quantile)
    00714 {
    00715     std::sort(noise.begin(), noise.end(), detail::SortNoiseByMean());
    00716
    00717     ArrayVector<TinyVector<unsigned int, 2> > clusters;
    00718     detail::noiseVarianceListMedianCut(noise, clusters, clusterCount);
    00719
    00720     std::sort(clusters.begin(), clusters.end(), detail::SortNoiseByMean());
    00721
    00722     detail::noiseVarianceClusterAveraging(noise, clusters, result, quantile);
    00723 }
    00724
    00725 template <class Functor,
    00726           class SrcIterator, class SrcAccessor,
    00727           class DestIterator, class DestAccessor>
    00728 bool
    00729 noiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00730                        DestIterator dul, DestAccessor dest,
    00731                        NoiseNormalizationOptions const & options)
    00732 {
    00733     ArrayVector<TinyVector<double, 2> > noiseData;
    00734     noiseVarianceEstimationImpl(sul, slr, src, noiseData, options);
    00735
    00736     if(noiseData.size() < 10)
    00737         return false;
    00738
    00739     std::sort(noiseData.begin(), noiseData.end(), SortNoiseByMean());
    00740
    00741     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00742
    00743     noiseVarianceClusteringImpl(noiseData, noiseClusters,
    00744                                   options.cluster_count, options.averaging_quantile);
    00745
    00746     transformImage(sul, slr, src, dul, dest, Functor(noiseClusters));
    00747
    00748     return true;
    00749 }
    00750
    00751 template <class SrcIterator, class SrcAccessor,
    00752           class DestIterator, class DestAccessor>
    00753 bool
    00754 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00755                                     DestIterator dul, DestAccessor dest,
    00756                                     NoiseNormalizationOptions const & options,
    00757                                     VigraTrueType /* isScalar */)
    00758 {
    00759     typedef typename SrcAccessor::value_type SrcType;
    00760     typedef typename DestAccessor::value_type DestType;
    00761     return noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
    00762                                                          (sul, slr, src, dul, dest, options);
    00763 }
    00764
    00765 template <class SrcIterator, class SrcAccessor,
    00766           class DestIterator, class DestAccessor>
    00767 bool
    00768 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00769                               DestIterator dul, DestAccessor dest,
    00770                               NoiseNormalizationOptions const & options,
    00771                               VigraFalseType /* isScalar */)
    00772 {
    00773     int bands = src.size(sul);
    00774     for(int b=0; b<bands; ++b)
    00775     {
    00776         VectorElementAccessor<SrcAccessor> sband(b, src);
    00777         VectorElementAccessor<DestAccessor> dband(b, dest);
    00778         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00779         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00780
    00781         if(!noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
    00782                                                            (sul, slr, sband, dul, dband, options))
    00783             return false;
    00784     }
    00785     return true;
    00786 }
    00787
    00788 template <class SrcIterator, class SrcAccessor,
    00789           class DestIterator, class DestAccessor>
    00790 bool
    00791 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00792                                     DestIterator dul, DestAccessor dest,
    00793                                     NoiseNormalizationOptions const & options,
    00794                                     VigraTrueType /* isScalar */)
    00795 {
    00796     typedef typename SrcAccessor::value_type SrcType;
    00797     typedef typename DestAccessor::value_type DestType;
    00798     return noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
    00799                                                          (sul, slr, src, dul, dest, options);
    00800 }
    00801
    00802 template <class SrcIterator, class SrcAccessor,
    00803           class DestIterator, class DestAccessor>
    00804 bool
    00805 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00806                               DestIterator dul, DestAccessor dest,
    00807                               NoiseNormalizationOptions const & options,
    00808                               VigraFalseType /* isScalar */)
    00809 {
    00810     int bands = src.size(sul);
    00811     for(int b=0; b<bands; ++b)
    00812     {
    00813         VectorElementAccessor<SrcAccessor> sband(b, src);
    00814         VectorElementAccessor<DestAccessor> dband(b, dest);
    00815         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00816         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00817
    00818         if(!noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
    00819                                                            (sul, slr, sband, dul, dband, options))
    00820             return false;
    00821     }
    00822     return true;
    00823 }
    00824
    00825 template <class SrcIterator, class SrcAccessor,
    00826           class DestIterator, class DestAccessor>
    00827 void
    00828 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00829                               DestIterator dul, DestAccessor dest,
    00830                               double a0, double a1, double a2,
    00831                               VigraTrueType /* isScalar */)
    00832 {
    00833     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00834     noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
    00835     noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1 + a2));
    00836     noiseClusters.push_back(TinyVector<double, 2>(2.0, a0 + 2.0*a1 + 4.0*a2));
    00837     transformImage(sul, slr, src, dul, dest,
    00838                    QuadraticNoiseNormalizationFunctor<typename SrcAccessor::value_type,
    00839                                                    typename DestAccessor::value_type>(noiseClusters));
    00840 }
    00841
    00842 template <class SrcIterator, class SrcAccessor,
    00843           class DestIterator, class DestAccessor>
    00844 void
    00845 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00846                               DestIterator dul, DestAccessor dest,
    00847                               double a0, double a1, double a2,
    00848                               VigraFalseType /* isScalar */)
    00849 {
    00850     int bands = src.size(sul);
    00851     for(int b=0; b<bands; ++b)
    00852     {
    00853         VectorElementAccessor<SrcAccessor> sband(b, src);
    00854         VectorElementAccessor<DestAccessor> dband(b, dest);
    00855         quadraticNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, a2, VigraTrueType());
    00856     }
    00857 }
    00858
    00859 template <class SrcIterator, class SrcAccessor,
    00860           class DestIterator, class DestAccessor>
    00861 bool
    00862 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00863                                     DestIterator dul, DestAccessor dest,
    00864                                     NoiseNormalizationOptions const & options,
    00865                                     VigraTrueType /* isScalar */)
    00866 {
    00867     typedef typename SrcAccessor::value_type SrcType;
    00868     typedef typename DestAccessor::value_type DestType;
    00869     return noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
    00870                                                          (sul, slr, src, dul, dest, options);
    00871 }
    00872
    00873 template <class SrcIterator, class SrcAccessor,
    00874           class DestIterator, class DestAccessor>
    00875 bool
    00876 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00877                               DestIterator dul, DestAccessor dest,
    00878                               NoiseNormalizationOptions const & options,
    00879                               VigraFalseType /* isScalar */)
    00880 {
    00881     int bands = src.size(sul);
    00882     for(int b=0; b<bands; ++b)
    00883     {
    00884         VectorElementAccessor<SrcAccessor> sband(b, src);
    00885         VectorElementAccessor<DestAccessor> dband(b, dest);
    00886         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00887         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00888
    00889         if(!noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
    00890                                                            (sul, slr, sband, dul, dband, options))
    00891             return false;
    00892     }
    00893     return true;
    00894 }
    00895
    00896 template <class SrcIterator, class SrcAccessor,
    00897           class DestIterator, class DestAccessor>
    00898 void
    00899 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00900                               DestIterator dul, DestAccessor dest,
    00901                               double a0, double a1,
    00902                               VigraTrueType /* isScalar */)
    00903 {
    00904     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00905     noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
    00906     noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1));
    00907     transformImage(sul, slr, src, dul, dest,
    00908                    LinearNoiseNormalizationFunctor<typename SrcAccessor::value_type,
    00909                                                    typename DestAccessor::value_type>(noiseClusters));
    00910 }
    00911
    00912 template <class SrcIterator, class SrcAccessor,
    00913           class DestIterator, class DestAccessor>
    00914 void
    00915 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00916                               DestIterator dul, DestAccessor dest,
    00917                               double a0, double a1,
    00918                               VigraFalseType /* isScalar */)
    00919 {
    00920     int bands = src.size(sul);
    00921     for(int b=0; b<bands; ++b)
    00922     {
    00923         VectorElementAccessor<SrcAccessor> sband(b, src);
    00924         VectorElementAccessor<DestAccessor> dband(b, dest);
    00925         linearNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, VigraTrueType());
    00926     }
    00927 }
    00928
    00929 } // namespace detail
    00930
    00931 template <bool P>
    00932 struct noiseVarianceEstimation_can_only_work_on_scalar_images
    00933 : vigra::staticAssert::AssertBool<P>
    00934 {};
    00935
    00936 /** \addtogroup NoiseNormalization Noise Normalization
    00937     Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
    00938 */
    00939 //@{
    00940
    00941 /********************************************************/
    00942 /*                                                      */
    00943 /*                noiseVarianceEstimation               */
    00944 /*                                                      */
    00945 /********************************************************/
    00946
    00947 /** \brief Determine the noise variance as a function of the image intensity.
    00948
    00949     This operator applies an algorithm described in
    00950
    00951     W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
    00952     Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
    00953     Lecture Notes in Earth Science, Berlin: Springer, 1999
    00954
    00955     in order to estimate the noise variance as a function of the image intensity in a robust way,
    00956     i.e. so that intensity changes due to edges do not bias the estimate. The source value type
    00957     (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
    00958     The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible
    00959     from two <tt>double</tt> values. The following options can be set via the \a options object
    00960     (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
    00961
    00962     <tt>useGradient</tt>, <tt>windowRadius</tt>, <tt>noiseEstimationQuantile</tt>, <tt>noiseVarianceInitialGuess</tt>
    00963
    00964     <b> Declarations:</b>
    00965
    00966     pass arguments explicitly:
    00967     \code
    00968     namespace vigra {
    00969         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00970         void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00971                                      BackInsertable & result,
    00972                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    00973     }
    00974     \endcode
    00975
    00976     use argument objects in conjunction with \ref ArgumentObjectFactories:
    00977     \code
    00978     namespace vigra {
    00979         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00980         void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00981                                      BackInsertable & result,
    00982                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    00983     }
    00984     \endcode
    00985
    00986     <b> Usage:</b>
    00987
    00988         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    00989     Namespace: vigra
    00990
    00991     \code
    00992     vigra::BImage src(w,h);
    00993     std::vector<vigra::TinyVector<double, 2> > result;
    00994
    00995     ...
    00996     vigra::noiseVarianceEstimation(srcImageRange(src), result,
    00997                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
    00998
    00999     // print the intensity / variance pairs found
    01000     for(int k=0; k<result.size(); ++k)
    01001         std::cout << "Intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
    01002     \endcode
    01003
    01004     <b> Required Interface:</b>
    01005
    01006     \code
    01007     SrcIterator upperleft, lowerright;
    01008     SrcAccessor src;
    01009
    01010     typedef SrcAccessor::value_type SrcType;
    01011     typedef NumericTraits<SrcType>::isScalar isScalar;
    01012     assert(isScalar::asBool == true);
    01013
    01014     double value = src(uperleft);
    01015
    01016     BackInsertable result;
    01017     typedef BackInsertable::value_type ResultType;
    01018     double intensity, variance;
    01019     result.push_back(ResultType(intensity, variance));
    01020     \endcode
    01021 */
    01022 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01023 inline
    01024 void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01025                            BackInsertable & result,
    01026                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01027 {
    01028     typedef typename BackInsertable::value_type ResultType;
    01029     typedef typename SrcAccessor::value_type SrcType;
    01030     typedef typename NumericTraits<SrcType>::isScalar isScalar;
    01031
    01032     VIGRA_STATIC_ASSERT((
    01033         noiseVarianceEstimation_can_only_work_on_scalar_images<(isScalar::asBool)>));
    01034
    01035     detail::noiseVarianceEstimationImpl(sul, slr, src, result, options);
    01036 }
    01037
    01038 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01039 inline
    01040 void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01041                            BackInsertable & result,
    01042                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01043 {
    01044     noiseVarianceEstimation(src.first, src.second, src.third, result, options);
    01045 }
    01046
    01047 /********************************************************/
    01048 /*                                                      */
    01049 /*                noiseVarianceClustering               */
    01050 /*                                                      */
    01051 /********************************************************/
    01052
    01053 /** \brief Determine the noise variance as a function of the image intensity and cluster the results.
    01054
    01055     This operator first calls \ref noiseVarianceEstimation() to obtain a sequence of intensity/variance pairs,
    01056     which are then clustered using the median cut algorithm. Then the cluster centers (i.e. average variance vs.
    01057     average intensity) are determined and returned in the \a result sequence.
    01058
    01059     In addition to the options valid for \ref noiseVarianceEstimation(), the following options can be set via
    01060     the \a options object (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
    01061
    01062     <tt>clusterCount</tt>, <tt>averagingQuantile</tt>
    01063
    01064     <b> Declarations:</b>
    01065
    01066     pass arguments explicitly:
    01067     \code
    01068     namespace vigra {
    01069         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01070         void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01071                                 BackInsertable & result,
    01072                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01073     }
    01074     \endcode
    01075
    01076     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01077     \code
    01078     namespace vigra {
    01079         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01080         void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01081                                 BackInsertable & result,
    01082                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01083     }
    01084     \endcode
    01085
    01086     <b> Usage:</b>
    01087
    01088         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01089     Namespace: vigra
    01090
    01091     \code
    01092     vigra::BImage src(w,h);
    01093     std::vector<vigra::TinyVector<double, 2> > result;
    01094
    01095     ...
    01096     vigra::noiseVarianceClustering(srcImageRange(src), result,
    01097                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01098                                   clusterCount(15));
    01099
    01100     // print the intensity / variance pairs representing the cluster centers
    01101     for(int k=0; k<result.size(); ++k)
    01102         std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
    01103     \endcode
    01104
    01105     <b> Required Interface:</b>
    01106
    01107     same as \ref noiseVarianceEstimation()
    01108 */
    01109 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01110 inline
    01111 void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01112                            BackInsertable & result,
    01113                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01114 {
    01115     ArrayVector<TinyVector<double, 2> > variance;
    01116     noiseVarianceEstimation(sul, slr, src, variance, options);
    01117     detail::noiseVarianceClusteringImpl(variance, result, options.cluster_count, options.averaging_quantile);
    01118 }
    01119
    01120 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01121 inline
    01122 void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01123                            BackInsertable & result,
    01124                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01125 {
    01126     noiseVarianceClustering(src.first, src.second, src.third, result, options);
    01127 }
    01128
    01129 /********************************************************/
    01130 /*                                                      */
    01131 /*             nonparametricNoiseNormalization          */
    01132 /*                                                      */
    01133 /********************************************************/
    01134
    01135 /** \brief Noise normalization by means of an estimated non-parametric noise model.
    01136
    01137     The original image is assumed to be corrupted by noise whose variance depends on the intensity in an unknown way.
    01138     The present functions first calls \ref noiseVarianceClustering() to obtain a sequence of intensity/variance pairs
    01139     (cluster centers) which estimate this dependency. The cluster centers are connected into a piecewise linear
    01140     function which is the inverted according to the formula derived in
    01141
    01142     W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
    01143     Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
    01144     Lecture Notes in Earth Science, Berlin: Springer, 1999
    01145
    01146     The inverted formula defines a pixel-wise intensity transformation whose application turns the original image
    01147     into one that is corrupted by additive Gaussian noise with unit variance. Most subsequent algorithms will be able
    01148     to handle this type of noise much better than the original noise.
    01149
    01150     RGB and other multiband images will be processed one band at a time. The function returns <tt>true</tt> on success.
    01151     Noise normalization will fail if the original image does not contain sufficiently homogeneous regions that
    01152     allow robust estimation of the noise variance.
    01153
    01154     The \a options object may use all options described in \ref vigra::NoiseNormalizationOptions.
    01155
    01156     <b> Declarations:</b>
    01157
    01158     pass arguments explicitly:
    01159     \code
    01160     namespace vigra {
    01161         template <class SrcIterator, class SrcAccessor,
    01162                   class DestIterator, class DestAccessor>
    01163         bool nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01164                                              DestIterator dul, DestAccessor dest,
    01165                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01166     }
    01167     \endcode
    01168
    01169     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01170     \code
    01171     namespace vigra {
    01172         template <class SrcIterator, class SrcAccessor,
    01173                   class DestIterator, class DestAccessor>
    01174         bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01175                                              pair<DestIterator, DestAccessor> dest,
    01176                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01177     }
    01178     \endcode
    01179
    01180     <b> Usage:</b>
    01181
    01182         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01183     Namespace: vigra
    01184
    01185     \code
    01186     vigra::BRGBImage src(w,h), dest(w, h);
    01187
    01188     ...
    01189     vigra::nonparametricNoiseNormalization(srcImageRange(src), destImage(dest),
    01190                                            vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01191                                            clusterCount(15));
    01192     \endcode
    01193
    01194     <b> Required Interface:</b>
    01195
    01196     same as \ref noiseVarianceEstimation()
    01197 */
    01198 template <class SrcIterator, class SrcAccessor,
    01199           class DestIterator, class DestAccessor>
    01200 inline bool
    01201 nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01202                                 DestIterator dul, DestAccessor dest,
    01203                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01204 {
    01205     typedef typename SrcAccessor::value_type SrcType;
    01206
    01207     return detail::nonparametricNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01208                                          typename NumericTraits<SrcType>::isScalar());
    01209 }
    01210
    01211 template <class SrcIterator, class SrcAccessor,
    01212           class DestIterator, class DestAccessor>
    01213 inline
    01214 bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01215                                      pair<DestIterator, DestAccessor> dest,
    01216                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01217 {
    01218     return nonparametricNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01219 }
    01220
    01221 /********************************************************/
    01222 /*                                                      */
    01223 /*               quadraticNoiseNormalization            */
    01224 /*                                                      */
    01225 /********************************************************/
    01226
    01227 /** \brief Noise normalization by means of an estimated quadratic noise model.
    01228
    01229     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01230     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01231     quadratic function rather than a piecewise linear function. If the quadratic model is applicable, it leads
    01232     to a somewhat smoother transformation.
    01233
    01234     <b> Declarations:</b>
    01235
    01236     pass arguments explicitly:
    01237     \code
    01238     namespace vigra {
    01239         template <class SrcIterator, class SrcAccessor,
    01240                   class DestIterator, class DestAccessor>
    01241         bool quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01242                                          DestIterator dul, DestAccessor dest,
    01243                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01244     }
    01245     \endcode
    01246
    01247     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01248     \code
    01249     namespace vigra {
    01250         template <class SrcIterator, class SrcAccessor,
    01251                   class DestIterator, class DestAccessor>
    01252         bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01253                                          pair<DestIterator, DestAccessor> dest,
    01254                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01255     }
    01256     \endcode
    01257
    01258     <b> Usage:</b>
    01259
    01260         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01261     Namespace: vigra
    01262
    01263     \code
    01264     vigra::BRGBImage src(w,h), dest(w, h);
    01265
    01266     ...
    01267     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01268                                        vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01269                                        clusterCount(15));
    01270     \endcode
    01271
    01272     <b> Required Interface:</b>
    01273
    01274     same as \ref noiseVarianceEstimation()
    01275 */
    01276 template <class SrcIterator, class SrcAccessor,
    01277           class DestIterator, class DestAccessor>
    01278 inline bool
    01279 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01280                                 DestIterator dul, DestAccessor dest,
    01281                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01282 {
    01283     typedef typename SrcAccessor::value_type SrcType;
    01284
    01285     return detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01286                                          typename NumericTraits<SrcType>::isScalar());
    01287 }
    01288
    01289 template <class SrcIterator, class SrcAccessor,
    01290           class DestIterator, class DestAccessor>
    01291 inline
    01292 bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01293                                      pair<DestIterator, DestAccessor> dest,
    01294                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01295 {
    01296     return quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01297 }
    01298
    01299 /********************************************************/
    01300 /*                                                      */
    01301 /*               quadraticNoiseNormalization            */
    01302 /*                                                      */
    01303 /********************************************************/
    01304
    01305 /** \brief Noise normalization by means of a given quadratic noise model.
    01306
    01307     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01308     functional dependency of the noise variance from the intensity is given (by a quadratic function)
    01309     rather than estimated:
    01310
    01311     \code
    01312     variance = a0 + a1 * intensity + a2 * sq(intensity)
    01313     \endcode
    01314
    01315     <b> Declarations:</b>
    01316
    01317     pass arguments explicitly:
    01318     \code
    01319     namespace vigra {
    01320         template <class SrcIterator, class SrcAccessor,
    01321                   class DestIterator, class DestAccessor>
    01322         void quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01323                                          DestIterator dul, DestAccessor dest,
    01324                                          double a0, double a1, double a2);
    01325     }
    01326     \endcode
    01327
    01328     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01329     \code
    01330     namespace vigra {
    01331         template <class SrcIterator, class SrcAccessor,
    01332                   class DestIterator, class DestAccessor>
    01333         void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01334                                         pair<DestIterator, DestAccessor> dest,
    01335                                         double a0, double a1, double a2);
    01336     }
    01337     \endcode
    01338
    01339     <b> Usage:</b>
    01340
    01341         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01342     Namespace: vigra
    01343
    01344     \code
    01345     vigra::BRGBImage src(w,h), dest(w, h);
    01346
    01347     ...
    01348     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01349                                        100, 0.02, 1e-6);
    01350     \endcode
    01351
    01352     <b> Required Interface:</b>
    01353
    01354     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01355     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01356     or a vector whose elements are assignable from <tt>double</tt>.
    01357 */
    01358 template <class SrcIterator, class SrcAccessor,
    01359           class DestIterator, class DestAccessor>
    01360 inline void
    01361 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01362                             DestIterator dul, DestAccessor dest,
    01363                             double a0, double a1, double a2)
    01364 {
    01365     typedef typename SrcAccessor::value_type SrcType;
    01366
    01367     detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1, a2,
    01368                                          typename NumericTraits<SrcType>::isScalar());
    01369 }
    01370
    01371 template <class SrcIterator, class SrcAccessor,
    01372           class DestIterator, class DestAccessor>
    01373 inline
    01374 void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01375                                  pair<DestIterator, DestAccessor> dest,
    01376                                  double a0, double a1, double a2)
    01377 {
    01378     quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1, a2);
    01379 }
    01380
    01381 /********************************************************/
    01382 /*                                                      */
    01383 /*                linearNoiseNormalization              */
    01384 /*                                                      */
    01385 /********************************************************/
    01386
    01387 /** \brief Noise normalization by means of an estimated linear noise model.
    01388
    01389     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01390     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01391     linear function rather than a piecewise linear function. If the linear model is applicable, it leads
    01392     to a very simple transformation which is similar to the familiar gamma correction.
    01393
    01394     <b> Declarations:</b>
    01395
    01396     pass arguments explicitly:
    01397     \code
    01398     namespace vigra {
    01399         template <class SrcIterator, class SrcAccessor,
    01400                 class DestIterator, class DestAccessor>
    01401         bool linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01402                                       DestIterator dul, DestAccessor dest,
    01403                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01404     }
    01405     \endcode
    01406
    01407     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01408     \code
    01409     namespace vigra {
    01410         template <class SrcIterator, class SrcAccessor,
    01411                   class DestIterator, class DestAccessor>
    01412         bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01413                                       pair<DestIterator, DestAccessor> dest,
    01414                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01415     }
    01416     \endcode
    01417
    01418     <b> Usage:</b>
    01419
    01420         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01421     Namespace: vigra
    01422
    01423     \code
    01424     vigra::BRGBImage src(w,h), dest(w, h);
    01425
    01426     ...
    01427     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01428                                     vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01429                                     clusterCount(15));
    01430     \endcode
    01431
    01432     <b> Required Interface:</b>
    01433
    01434     same as \ref noiseVarianceEstimation()
    01435 */
    01436 template <class SrcIterator, class SrcAccessor,
    01437           class DestIterator, class DestAccessor>
    01438 inline bool
    01439 linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01440                                 DestIterator dul, DestAccessor dest,
    01441                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01442 {
    01443     typedef typename SrcAccessor::value_type SrcType;
    01444
    01445     return detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01446                                          typename NumericTraits<SrcType>::isScalar());
    01447 }
    01448
    01449 template <class SrcIterator, class SrcAccessor,
    01450           class DestIterator, class DestAccessor>
    01451 inline
    01452 bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01453                                      pair<DestIterator, DestAccessor> dest,
    01454                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01455 {
    01456     return linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01457 }
    01458
    01459 /********************************************************/
    01460 /*                                                      */
    01461 /*                 linearNoiseNormalization             */
    01462 /*                                                      */
    01463 /********************************************************/
    01464
    01465 /** \brief Noise normalization by means of a given linear noise model.
    01466
    01467     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01468     functional dependency of the noise variance from the intensity is given (as a linear function)
    01469     rather than estimated:
    01470
    01471     \code
    01472     variance = a0 + a1 * intensity
    01473     \endcode
    01474
    01475     <b> Declarations:</b>
    01476
    01477     pass arguments explicitly:
    01478     \code
    01479     namespace vigra {
    01480         template <class SrcIterator, class SrcAccessor,
    01481                   class DestIterator, class DestAccessor>
    01482         void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01483                                       DestIterator dul, DestAccessor dest,
    01484                                       double a0, double a1);
    01485     }
    01486     \endcode
    01487
    01488     use argument objects in conjunction with \ref ArgumentObjectFactories:
    01489     \code
    01490     namespace vigra {
    01491         template <class SrcIterator, class SrcAccessor,
    01492                   class DestIterator, class DestAccessor>
    01493         void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01494                                       pair<DestIterator, DestAccessor> dest,
    01495                                       double a0, double a1);
    01496     }
    01497     \endcode
    01498
    01499     <b> Usage:</b>
    01500
    01501         <b>\#include</b> "<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>"<br>
    01502     Namespace: vigra
    01503
    01504     \code
    01505     vigra::BRGBImage src(w,h), dest(w, h);
    01506
    01507     ...
    01508     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01509                                     100, 0.02);
    01510     \endcode
    01511
    01512     <b> Required Interface:</b>
    01513
    01514     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01515     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01516     or a vector whose elements are assignable from <tt>double</tt>.
    01517 */
    01518 template <class SrcIterator, class SrcAccessor,
    01519           class DestIterator, class DestAccessor>
    01520 inline
    01521 void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01522                               DestIterator dul, DestAccessor dest,
    01523                               double a0, double a1)
    01524 {
    01525     typedef typename SrcAccessor::value_type SrcType;
    01526
    01527     detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1,
    01528                                          typename NumericTraits<SrcType>::isScalar());
    01529 }
    01530
    01531 template <class SrcIterator, class SrcAccessor,
    01532           class DestIterator, class DestAccessor>
    01533 inline
    01534 void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01535                               pair<DestIterator, DestAccessor> dest,
    01536                               double a0, double a1)
    01537 {
    01538     linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1);
    01539 }
    01540
    01541 //@}
    01542
    01543 } // namespace vigra
    01544
    01545 #endif // VIGRA_NOISE_NORMALIZATION_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.5.0 (7 Dec 2006)