[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | 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.6.0, Aug 13 2008 )                                    */
    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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00012 /*        vigra@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 doxygen_overloaded_function(template <...> void noiseVarianceEstimation)
    01023
    01024 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01025 inline
    01026 void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01027                            BackInsertable & result,
    01028                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01029 {
    01030     typedef typename BackInsertable::value_type ResultType;
    01031     typedef typename SrcAccessor::value_type SrcType;
    01032     typedef typename NumericTraits<SrcType>::isScalar isScalar;
    01033
    01034     VIGRA_STATIC_ASSERT((
    01035         noiseVarianceEstimation_can_only_work_on_scalar_images<(isScalar::asBool)>));
    01036
    01037     detail::noiseVarianceEstimationImpl(sul, slr, src, result, options);
    01038 }
    01039
    01040 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01041 inline
    01042 void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01043                            BackInsertable & result,
    01044                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01045 {
    01046     noiseVarianceEstimation(src.first, src.second, src.third, result, options);
    01047 }
    01048
    01049 /********************************************************/
    01050 /*                                                      */
    01051 /*                noiseVarianceClustering               */
    01052 /*                                                      */
    01053 /********************************************************/
    01054
    01055 /** \brief Determine the noise variance as a function of the image intensity and cluster the results.
    01056
    01057     This operator first calls \ref noiseVarianceEstimation() to obtain a sequence of intensity/variance pairs,
    01058     which are then clustered using the median cut algorithm. Then the cluster centers (i.e. average variance vs.
    01059     average intensity) are determined and returned in the \a result sequence.
    01060
    01061     In addition to the options valid for \ref noiseVarianceEstimation(), the following options can be set via
    01062     the \a options object (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
    01063
    01064     <tt>clusterCount</tt>, <tt>averagingQuantile</tt>
    01065
    01066     <b> Declarations:</b>
    01067
    01068     pass arguments explicitly:
    01069     \code
    01070     namespace vigra {
    01071         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01072         void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01073                                 BackInsertable & result,
    01074                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01075     }
    01076     \endcode
    01077
    01078     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01079     \code
    01080     namespace vigra {
    01081         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01082         void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01083                                 BackInsertable & result,
    01084                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01085     }
    01086     \endcode
    01087
    01088     <b> Usage:</b>
    01089
    01090         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01091     Namespace: vigra
    01092
    01093     \code
    01094     vigra::BImage src(w,h);
    01095     std::vector<vigra::TinyVector<double, 2> > result;
    01096
    01097     ...
    01098     vigra::noiseVarianceClustering(srcImageRange(src), result,
    01099                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01100                                   clusterCount(15));
    01101
    01102     // print the intensity / variance pairs representing the cluster centers
    01103     for(int k=0; k<result.size(); ++k)
    01104         std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
    01105     \endcode
    01106
    01107     <b> Required Interface:</b>
    01108
    01109     same as \ref noiseVarianceEstimation()
    01110 */
    01111 doxygen_overloaded_function(template <...> void noiseVarianceClustering)
    01112
    01113 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01114 inline
    01115 void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01116                            BackInsertable & result,
    01117                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01118 {
    01119     ArrayVector<TinyVector<double, 2> > variance;
    01120     noiseVarianceEstimation(sul, slr, src, variance, options);
    01121     detail::noiseVarianceClusteringImpl(variance, result, options.cluster_count, options.averaging_quantile);
    01122 }
    01123
    01124 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01125 inline
    01126 void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01127                            BackInsertable & result,
    01128                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01129 {
    01130     noiseVarianceClustering(src.first, src.second, src.third, result, options);
    01131 }
    01132
    01133 /********************************************************/
    01134 /*                                                      */
    01135 /*             nonparametricNoiseNormalization          */
    01136 /*                                                      */
    01137 /********************************************************/
    01138
    01139 /** \brief Noise normalization by means of an estimated non-parametric noise model.
    01140
    01141     The original image is assumed to be corrupted by noise whose variance depends on the intensity in an unknown way.
    01142     The present functions first calls \ref noiseVarianceClustering() to obtain a sequence of intensity/variance pairs
    01143     (cluster centers) which estimate this dependency. The cluster centers are connected into a piecewise linear
    01144     function which is the inverted according to the formula derived in
    01145
    01146     W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
    01147     Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
    01148     Lecture Notes in Earth Science, Berlin: Springer, 1999
    01149
    01150     The inverted formula defines a pixel-wise intensity transformation whose application turns the original image
    01151     into one that is corrupted by additive Gaussian noise with unit variance. Most subsequent algorithms will be able
    01152     to handle this type of noise much better than the original noise.
    01153
    01154     RGB and other multiband images will be processed one band at a time. The function returns <tt>true</tt> on success.
    01155     Noise normalization will fail if the original image does not contain sufficiently homogeneous regions that
    01156     allow robust estimation of the noise variance.
    01157
    01158     The \a options object may use all options described in \ref vigra::NoiseNormalizationOptions.
    01159
    01160     <b> Declarations:</b>
    01161
    01162     pass arguments explicitly:
    01163     \code
    01164     namespace vigra {
    01165         template <class SrcIterator, class SrcAccessor,
    01166                   class DestIterator, class DestAccessor>
    01167         bool nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01168                                              DestIterator dul, DestAccessor dest,
    01169                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01170     }
    01171     \endcode
    01172
    01173     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01174     \code
    01175     namespace vigra {
    01176         template <class SrcIterator, class SrcAccessor,
    01177                   class DestIterator, class DestAccessor>
    01178         bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01179                                              pair<DestIterator, DestAccessor> dest,
    01180                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01181     }
    01182     \endcode
    01183
    01184     <b> Usage:</b>
    01185
    01186         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01187     Namespace: vigra
    01188
    01189     \code
    01190     vigra::BRGBImage src(w,h), dest(w, h);
    01191
    01192     ...
    01193     vigra::nonparametricNoiseNormalization(srcImageRange(src), destImage(dest),
    01194                                            vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01195                                            clusterCount(15));
    01196     \endcode
    01197
    01198     <b> Required Interface:</b>
    01199
    01200     same as \ref noiseVarianceEstimation()
    01201 */
    01202 doxygen_overloaded_function(template <...> bool nonparametricNoiseNormalization)
    01203
    01204 template <class SrcIterator, class SrcAccessor,
    01205           class DestIterator, class DestAccessor>
    01206 inline bool
    01207 nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01208                                 DestIterator dul, DestAccessor dest,
    01209                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01210 {
    01211     typedef typename SrcAccessor::value_type SrcType;
    01212
    01213     return detail::nonparametricNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01214                                          typename NumericTraits<SrcType>::isScalar());
    01215 }
    01216
    01217 template <class SrcIterator, class SrcAccessor,
    01218           class DestIterator, class DestAccessor>
    01219 inline
    01220 bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01221                                      pair<DestIterator, DestAccessor> dest,
    01222                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01223 {
    01224     return nonparametricNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01225 }
    01226
    01227 /********************************************************/
    01228 /*                                                      */
    01229 /*               quadraticNoiseNormalization            */
    01230 /*                                                      */
    01231 /********************************************************/
    01232
    01233 /** \brief Noise normalization by means of an estimated quadratic noise model.
    01234
    01235     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01236     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01237     quadratic function rather than a piecewise linear function. If the quadratic model is applicable, it leads
    01238     to a somewhat smoother transformation.
    01239
    01240     <b> Declarations:</b>
    01241
    01242     pass arguments explicitly:
    01243     \code
    01244     namespace vigra {
    01245         template <class SrcIterator, class SrcAccessor,
    01246                   class DestIterator, class DestAccessor>
    01247         bool quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01248                                          DestIterator dul, DestAccessor dest,
    01249                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01250     }
    01251     \endcode
    01252
    01253     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01254     \code
    01255     namespace vigra {
    01256         template <class SrcIterator, class SrcAccessor,
    01257                   class DestIterator, class DestAccessor>
    01258         bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01259                                          pair<DestIterator, DestAccessor> dest,
    01260                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01261     }
    01262     \endcode
    01263
    01264     <b> Usage:</b>
    01265
    01266         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01267     Namespace: vigra
    01268
    01269     \code
    01270     vigra::BRGBImage src(w,h), dest(w, h);
    01271
    01272     ...
    01273     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01274                                        vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01275                                        clusterCount(15));
    01276     \endcode
    01277
    01278     <b> Required Interface:</b>
    01279
    01280     same as \ref noiseVarianceEstimation()
    01281 */
    01282 doxygen_overloaded_function(template <...> bool quadraticNoiseNormalization)
    01283
    01284 template <class SrcIterator, class SrcAccessor,
    01285           class DestIterator, class DestAccessor>
    01286 inline bool
    01287 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01288                                 DestIterator dul, DestAccessor dest,
    01289                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01290 {
    01291     typedef typename SrcAccessor::value_type SrcType;
    01292
    01293     return detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01294                                          typename NumericTraits<SrcType>::isScalar());
    01295 }
    01296
    01297 template <class SrcIterator, class SrcAccessor,
    01298           class DestIterator, class DestAccessor>
    01299 inline
    01300 bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01301                                      pair<DestIterator, DestAccessor> dest,
    01302                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01303 {
    01304     return quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01305 }
    01306
    01307 /********************************************************/
    01308 /*                                                      */
    01309 /*               quadraticNoiseNormalization            */
    01310 /*                                                      */
    01311 /********************************************************/
    01312
    01313 /** \brief Noise normalization by means of a given quadratic noise model.
    01314
    01315     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01316     functional dependency of the noise variance from the intensity is given (by a quadratic function)
    01317     rather than estimated:
    01318
    01319     \code
    01320     variance = a0 + a1 * intensity + a2 * sq(intensity)
    01321     \endcode
    01322
    01323     <b> Declarations:</b>
    01324
    01325     pass arguments explicitly:
    01326     \code
    01327     namespace vigra {
    01328         template <class SrcIterator, class SrcAccessor,
    01329                   class DestIterator, class DestAccessor>
    01330         void quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01331                                          DestIterator dul, DestAccessor dest,
    01332                                          double a0, double a1, double a2);
    01333     }
    01334     \endcode
    01335
    01336     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01337     \code
    01338     namespace vigra {
    01339         template <class SrcIterator, class SrcAccessor,
    01340                   class DestIterator, class DestAccessor>
    01341         void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01342                                         pair<DestIterator, DestAccessor> dest,
    01343                                         double a0, double a1, double a2);
    01344     }
    01345     \endcode
    01346
    01347     <b> Usage:</b>
    01348
    01349         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01350     Namespace: vigra
    01351
    01352     \code
    01353     vigra::BRGBImage src(w,h), dest(w, h);
    01354
    01355     ...
    01356     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01357                                        100, 0.02, 1e-6);
    01358     \endcode
    01359
    01360     <b> Required Interface:</b>
    01361
    01362     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01363     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01364     or a vector whose elements are assignable from <tt>double</tt>.
    01365 */
    01366 doxygen_overloaded_function(template <...> void quadraticNoiseNormalization)
    01367
    01368 template <class SrcIterator, class SrcAccessor,
    01369           class DestIterator, class DestAccessor>
    01370 inline void
    01371 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01372                             DestIterator dul, DestAccessor dest,
    01373                             double a0, double a1, double a2)
    01374 {
    01375     typedef typename SrcAccessor::value_type SrcType;
    01376
    01377     detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1, a2,
    01378                                          typename NumericTraits<SrcType>::isScalar());
    01379 }
    01380
    01381 template <class SrcIterator, class SrcAccessor,
    01382           class DestIterator, class DestAccessor>
    01383 inline
    01384 void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01385                                  pair<DestIterator, DestAccessor> dest,
    01386                                  double a0, double a1, double a2)
    01387 {
    01388     quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1, a2);
    01389 }
    01390
    01391 /********************************************************/
    01392 /*                                                      */
    01393 /*                linearNoiseNormalization              */
    01394 /*                                                      */
    01395 /********************************************************/
    01396
    01397 /** \brief Noise normalization by means of an estimated linear noise model.
    01398
    01399     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01400     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01401     linear function rather than a piecewise linear function. If the linear model is applicable, it leads
    01402     to a very simple transformation which is similar to the familiar gamma correction.
    01403
    01404     <b> Declarations:</b>
    01405
    01406     pass arguments explicitly:
    01407     \code
    01408     namespace vigra {
    01409         template <class SrcIterator, class SrcAccessor,
    01410                 class DestIterator, class DestAccessor>
    01411         bool linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01412                                       DestIterator dul, DestAccessor dest,
    01413                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01414     }
    01415     \endcode
    01416
    01417     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01418     \code
    01419     namespace vigra {
    01420         template <class SrcIterator, class SrcAccessor,
    01421                   class DestIterator, class DestAccessor>
    01422         bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01423                                       pair<DestIterator, DestAccessor> dest,
    01424                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01425     }
    01426     \endcode
    01427
    01428     <b> Usage:</b>
    01429
    01430         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01431     Namespace: vigra
    01432
    01433     \code
    01434     vigra::BRGBImage src(w,h), dest(w, h);
    01435
    01436     ...
    01437     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01438                                     vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01439                                     clusterCount(15));
    01440     \endcode
    01441
    01442     <b> Required Interface:</b>
    01443
    01444     same as \ref noiseVarianceEstimation()
    01445 */
    01446 doxygen_overloaded_function(template <...> bool linearNoiseNormalization)
    01447
    01448 template <class SrcIterator, class SrcAccessor,
    01449           class DestIterator, class DestAccessor>
    01450 inline bool
    01451 linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01452                                 DestIterator dul, DestAccessor dest,
    01453                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01454 {
    01455     typedef typename SrcAccessor::value_type SrcType;
    01456
    01457     return detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01458                                          typename NumericTraits<SrcType>::isScalar());
    01459 }
    01460
    01461 template <class SrcIterator, class SrcAccessor,
    01462           class DestIterator, class DestAccessor>
    01463 inline
    01464 bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01465                                      pair<DestIterator, DestAccessor> dest,
    01466                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01467 {
    01468     return linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01469 }
    01470
    01471 /********************************************************/
    01472 /*                                                      */
    01473 /*                 linearNoiseNormalization             */
    01474 /*                                                      */
    01475 /********************************************************/
    01476
    01477 /** \brief Noise normalization by means of a given linear noise model.
    01478
    01479     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01480     functional dependency of the noise variance from the intensity is given (as a linear function)
    01481     rather than estimated:
    01482
    01483     \code
    01484     variance = a0 + a1 * intensity
    01485     \endcode
    01486
    01487     <b> Declarations:</b>
    01488
    01489     pass arguments explicitly:
    01490     \code
    01491     namespace vigra {
    01492         template <class SrcIterator, class SrcAccessor,
    01493                   class DestIterator, class DestAccessor>
    01494         void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01495                                       DestIterator dul, DestAccessor dest,
    01496                                       double a0, double a1);
    01497     }
    01498     \endcode
    01499
    01500     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01501     \code
    01502     namespace vigra {
    01503         template <class SrcIterator, class SrcAccessor,
    01504                   class DestIterator, class DestAccessor>
    01505         void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01506                                       pair<DestIterator, DestAccessor> dest,
    01507                                       double a0, double a1);
    01508     }
    01509     \endcode
    01510
    01511     <b> Usage:</b>
    01512
    01513         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01514     Namespace: vigra
    01515
    01516     \code
    01517     vigra::BRGBImage src(w,h), dest(w, h);
    01518
    01519     ...
    01520     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01521                                     100, 0.02);
    01522     \endcode
    01523
    01524     <b> Required Interface:</b>
    01525
    01526     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01527     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01528     or a vector whose elements are assignable from <tt>double</tt>.
    01529 */
    01530 doxygen_overloaded_function(template <...> void linearNoiseNormalization)
    01531
    01532 template <class SrcIterator, class SrcAccessor,
    01533           class DestIterator, class DestAccessor>
    01534 inline
    01535 void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01536                               DestIterator dul, DestAccessor dest,
    01537                               double a0, double a1)
    01538 {
    01539     typedef typename SrcAccessor::value_type SrcType;
    01540
    01541     detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1,
    01542                                          typename NumericTraits<SrcType>::isScalar());
    01543 }
    01544
    01545 template <class SrcIterator, class SrcAccessor,
    01546           class DestIterator, class DestAccessor>
    01547 inline
    01548 void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01549                               pair<DestIterator, DestAccessor> dest,
    01550                               double a0, double a1)
    01551 {
    01552     linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1);
    01553 }
    01554
    01555 //@}
    01556
    01557 } // namespace vigra
    01558
    01559 #endif // VIGRA_NOISE_NORMALIZATION_HXX

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python VIGRA 1.6.0 (13 Aug 2008)