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