[ 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> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><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 quess
    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             std::string message("noiseVarianceListMedianCut(): internal error (");
    00587             message += std::string("k: ") + asString(k) + ", ";
    00588             message += std::string("k1: ") + asString(k1) + ", ";
    00589             message += std::string("k2: ") + asString(k2) + ", ";
    00590             message += std::string("noise.size(): ") + asString(noise.size()) + ", ";
    00591             message += std::string("clusters.size(): ") + asString(clusters.size()) + ").";
    00592             vigra_invariant(k1 >= 0 && k1 < (int)noise.size() && k2 >= 0 && k2 < (int)noise.size(), message.c_str());
    00593
    00594             double diff = noise[k2][0] - noise[k1][0];
    00595             if(diff > diffMax)
    00596             {
    00597                 diffMax = diff;
    00598                 kMax = k;
    00599             }
    00600         }
    00601
    00602         if(diffMax == 0.0)
    00603             return; // all clusters have only one value
    00604
    00605         unsigned int k1 = clusters[kMax][0],
    00606                      k2 = clusters[kMax][1];
    00607         unsigned int kSplit = k1 + (k2 - k1) / 2;
    00608         clusters[kMax][1] = kSplit;
    00609         clusters.push_back(Result(kSplit, k2));
    00610     }
    00611 }
    00612
    00613 struct SortNoiseByMean
    00614 {
    00615     template <class T>
    00616     bool operator()(T const & l, T const & r) const
    00617     {
    00618         return l[0] < r[0];
    00619     }
    00620 };
    00621
    00622 struct SortNoiseByVariance
    00623 {
    00624     template <class T>
    00625     bool operator()(T const & l, T const & r) const
    00626     {
    00627         return l[1] < r[1];
    00628     }
    00629 };
    00630
    00631 template <class Vector1, class Vector2, class Vector3>
    00632 void noiseVarianceClusterAveraging(Vector1 & noise, Vector2 & clusters,
    00633                                    Vector3 & result, double quantile)
    00634 {
    00635     typedef typename Vector1::iterator Iter;
    00636     typedef typename Vector3::value_type Result;
    00637
    00638     for(unsigned int k=0; k<clusters.size(); ++k)
    00639     {
    00640         Iter i1 = noise.begin() + clusters[k][0];
    00641         Iter i2 = noise.begin() + clusters[k][1];
    00642
    00643         std::sort(i1, i2, SortNoiseByVariance());
    00644
    00645         std::size_t size = static_cast<std::size_t>(VIGRA_CSTD::ceil(quantile*(i2 - i1)));
    00646         if(static_cast<std::size_t>(i2 - i1) < size)
    00647             size = i2 - i1;
    00648         if(size < 1)
    00649             size = 1;
    00650         i2 = i1 + size;
    00651
    00652         double mean = 0.0,
    00653                variance = 0.0;
    00654         for(; i1 < i2; ++i1)
    00655         {
    00656             mean += (*i1)[0];
    00657             variance += (*i1)[1];
    00658         }
    00659
    00660         result.push_back(Result(mean / size, variance / size));
    00661     }
    00662 }
    00663
    00664 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00665 void noiseVarianceEstimationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00666                            BackInsertable & result,
    00667                            NoiseNormalizationOptions const & options)
    00668 {
    00669     typedef typename BackInsertable::value_type ResultType;
    00670
    00671     unsigned int w = slr.x - sul.x;
    00672     unsigned int h = slr.y - sul.y;
    00673
    00674     typedef typename NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
    00675     typedef BasicImage<TmpType> TmpImage;
    00676
    00677     TmpImage gradient(w, h);
    00678     symmetricDifferenceSquaredMagnitude(sul, slr, src, gradient.upperLeft(), gradient.accessor());
    00679
    00680     BImage homogeneous(w, h);
    00681     findHomogeneousRegions(gradient.upperLeft(), gradient.lowerRight(), gradient.accessor(),
    00682                                    homogeneous.upperLeft(), homogeneous.accessor());
    00683
    00684     // Generate noise of each of the remaining pixels == centers of homogenous areas (border is not used)
    00685     unsigned int windowRadius = options.window_radius;
    00686     for(unsigned int y=windowRadius; y<h-windowRadius; ++y)
    00687     {
    00688         for(unsigned int x=windowRadius; x<w-windowRadius; ++x)
    00689         {
    00690             if (! homogeneous(x, y))
    00691                 continue;
    00692
    00693             Diff2D center(x, y);
    00694             double mean = 0.0, variance = options.noise_variance_initial_guess;
    00695
    00696             bool success;
    00697
    00698             if(options.use_gradient)
    00699             {
    00700                 success = iterativeNoiseEstimationChi2(sul + center, src,
    00701                               gradient.upperLeft() + center, mean, variance,
    00702                               options.noise_estimation_quantile, windowRadius);
    00703             }
    00704             else
    00705             {
    00706                 success = iterativeNoiseEstimationGauss(sul + center, src,
    00707                               gradient.upperLeft() + center, mean, variance,
    00708                               options.noise_estimation_quantile, windowRadius);
    00709             }
    00710             if (success)
    00711             {
    00712                 result.push_back(ResultType(mean, variance));
    00713             }
    00714         }
    00715     }
    00716 }
    00717
    00718 template <class Vector, class BackInsertable>
    00719 void noiseVarianceClusteringImpl(Vector & noise, BackInsertable & result,
    00720                            unsigned int clusterCount, double quantile)
    00721 {
    00722     std::sort(noise.begin(), noise.end(), detail::SortNoiseByMean());
    00723
    00724     ArrayVector<TinyVector<unsigned int, 2> > clusters;
    00725     detail::noiseVarianceListMedianCut(noise, clusters, clusterCount);
    00726
    00727     std::sort(clusters.begin(), clusters.end(), detail::SortNoiseByMean());
    00728
    00729     detail::noiseVarianceClusterAveraging(noise, clusters, result, quantile);
    00730 }
    00731
    00732 template <class Functor,
    00733           class SrcIterator, class SrcAccessor,
    00734           class DestIterator, class DestAccessor>
    00735 bool
    00736 noiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00737                        DestIterator dul, DestAccessor dest,
    00738                        NoiseNormalizationOptions const & options)
    00739 {
    00740     ArrayVector<TinyVector<double, 2> > noiseData;
    00741     noiseVarianceEstimationImpl(sul, slr, src, noiseData, options);
    00742
    00743     if(noiseData.size() < 10)
    00744         return false;
    00745
    00746     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00747
    00748     noiseVarianceClusteringImpl(noiseData, noiseClusters,
    00749                                   options.cluster_count, options.averaging_quantile);
    00750
    00751     transformImage(sul, slr, src, dul, dest, Functor(noiseClusters));
    00752
    00753     return true;
    00754 }
    00755
    00756 template <class SrcIterator, class SrcAccessor,
    00757           class DestIterator, class DestAccessor>
    00758 bool
    00759 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00760                                     DestIterator dul, DestAccessor dest,
    00761                                     NoiseNormalizationOptions const & options,
    00762                                     VigraTrueType /* isScalar */)
    00763 {
    00764     typedef typename SrcAccessor::value_type SrcType;
    00765     typedef typename DestAccessor::value_type DestType;
    00766     return noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
    00767                                                          (sul, slr, src, dul, dest, options);
    00768 }
    00769
    00770 template <class SrcIterator, class SrcAccessor,
    00771           class DestIterator, class DestAccessor>
    00772 bool
    00773 nonparametricNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00774                               DestIterator dul, DestAccessor dest,
    00775                               NoiseNormalizationOptions const & options,
    00776                               VigraFalseType /* isScalar */)
    00777 {
    00778     int bands = src.size(sul);
    00779     for(int b=0; b<bands; ++b)
    00780     {
    00781         VectorElementAccessor<SrcAccessor> sband(b, src);
    00782         VectorElementAccessor<DestAccessor> dband(b, dest);
    00783         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00784         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00785
    00786         if(!noiseNormalizationImpl<NonparametricNoiseNormalizationFunctor<SrcType, DestType> >
    00787                                                            (sul, slr, sband, dul, dband, options))
    00788             return false;
    00789     }
    00790     return true;
    00791 }
    00792
    00793 template <class SrcIterator, class SrcAccessor,
    00794           class DestIterator, class DestAccessor>
    00795 bool
    00796 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00797                                     DestIterator dul, DestAccessor dest,
    00798                                     NoiseNormalizationOptions const & options,
    00799                                     VigraTrueType /* isScalar */)
    00800 {
    00801     typedef typename SrcAccessor::value_type SrcType;
    00802     typedef typename DestAccessor::value_type DestType;
    00803     return noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
    00804                                                          (sul, slr, src, dul, dest, options);
    00805 }
    00806
    00807 template <class SrcIterator, class SrcAccessor,
    00808           class DestIterator, class DestAccessor>
    00809 bool
    00810 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00811                               DestIterator dul, DestAccessor dest,
    00812                               NoiseNormalizationOptions const & options,
    00813                               VigraFalseType /* isScalar */)
    00814 {
    00815     int bands = src.size(sul);
    00816     for(int b=0; b<bands; ++b)
    00817     {
    00818         VectorElementAccessor<SrcAccessor> sband(b, src);
    00819         VectorElementAccessor<DestAccessor> dband(b, dest);
    00820         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00821         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00822
    00823         if(!noiseNormalizationImpl<QuadraticNoiseNormalizationFunctor<SrcType, DestType> >
    00824                                                            (sul, slr, sband, dul, dband, options))
    00825             return false;
    00826     }
    00827     return true;
    00828 }
    00829
    00830 template <class SrcIterator, class SrcAccessor,
    00831           class DestIterator, class DestAccessor>
    00832 void
    00833 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00834                               DestIterator dul, DestAccessor dest,
    00835                               double a0, double a1, double a2,
    00836                               VigraTrueType /* isScalar */)
    00837 {
    00838     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00839     noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
    00840     noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1 + a2));
    00841     noiseClusters.push_back(TinyVector<double, 2>(2.0, a0 + 2.0*a1 + 4.0*a2));
    00842     transformImage(sul, slr, src, dul, dest,
    00843                    QuadraticNoiseNormalizationFunctor<typename SrcAccessor::value_type,
    00844                                                    typename DestAccessor::value_type>(noiseClusters));
    00845 }
    00846
    00847 template <class SrcIterator, class SrcAccessor,
    00848           class DestIterator, class DestAccessor>
    00849 void
    00850 quadraticNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00851                               DestIterator dul, DestAccessor dest,
    00852                               double a0, double a1, double a2,
    00853                               VigraFalseType /* isScalar */)
    00854 {
    00855     int bands = src.size(sul);
    00856     for(int b=0; b<bands; ++b)
    00857     {
    00858         VectorElementAccessor<SrcAccessor> sband(b, src);
    00859         VectorElementAccessor<DestAccessor> dband(b, dest);
    00860         quadraticNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, a2, VigraTrueType());
    00861     }
    00862 }
    00863
    00864 template <class SrcIterator, class SrcAccessor,
    00865           class DestIterator, class DestAccessor>
    00866 bool
    00867 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00868                                     DestIterator dul, DestAccessor dest,
    00869                                     NoiseNormalizationOptions const & options,
    00870                                     VigraTrueType /* isScalar */)
    00871 {
    00872     typedef typename SrcAccessor::value_type SrcType;
    00873     typedef typename DestAccessor::value_type DestType;
    00874     return noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
    00875                                                          (sul, slr, src, dul, dest, options);
    00876 }
    00877
    00878 template <class SrcIterator, class SrcAccessor,
    00879           class DestIterator, class DestAccessor>
    00880 bool
    00881 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00882                               DestIterator dul, DestAccessor dest,
    00883                               NoiseNormalizationOptions const & options,
    00884                               VigraFalseType /* isScalar */)
    00885 {
    00886     int bands = src.size(sul);
    00887     for(int b=0; b<bands; ++b)
    00888     {
    00889         VectorElementAccessor<SrcAccessor> sband(b, src);
    00890         VectorElementAccessor<DestAccessor> dband(b, dest);
    00891         typedef typename VectorElementAccessor<SrcAccessor>::value_type SrcType;
    00892         typedef typename VectorElementAccessor<DestAccessor>::value_type DestType;
    00893
    00894         if(!noiseNormalizationImpl<LinearNoiseNormalizationFunctor<SrcType, DestType> >
    00895                                                            (sul, slr, sband, dul, dband, options))
    00896             return false;
    00897     }
    00898     return true;
    00899 }
    00900
    00901 template <class SrcIterator, class SrcAccessor,
    00902           class DestIterator, class DestAccessor>
    00903 void
    00904 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00905                               DestIterator dul, DestAccessor dest,
    00906                               double a0, double a1,
    00907                               VigraTrueType /* isScalar */)
    00908 {
    00909     ArrayVector<TinyVector<double, 2> > noiseClusters;
    00910     noiseClusters.push_back(TinyVector<double, 2>(0.0, a0));
    00911     noiseClusters.push_back(TinyVector<double, 2>(1.0, a0 + a1));
    00912     transformImage(sul, slr, src, dul, dest,
    00913                    LinearNoiseNormalizationFunctor<typename SrcAccessor::value_type,
    00914                                                    typename DestAccessor::value_type>(noiseClusters));
    00915 }
    00916
    00917 template <class SrcIterator, class SrcAccessor,
    00918           class DestIterator, class DestAccessor>
    00919 void
    00920 linearNoiseNormalizationImpl(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00921                               DestIterator dul, DestAccessor dest,
    00922                               double a0, double a1,
    00923                               VigraFalseType /* isScalar */)
    00924 {
    00925     int bands = src.size(sul);
    00926     for(int b=0; b<bands; ++b)
    00927     {
    00928         VectorElementAccessor<SrcAccessor> sband(b, src);
    00929         VectorElementAccessor<DestAccessor> dband(b, dest);
    00930         linearNoiseNormalizationImpl(sul, slr, sband, dul, dband, a0, a1, VigraTrueType());
    00931     }
    00932 }
    00933
    00934 } // namespace detail
    00935
    00936 template <bool P>
    00937 struct noiseVarianceEstimation_can_only_work_on_scalar_images
    00938 : vigra::staticAssert::AssertBool<P>
    00939 {};
    00940
    00941 /** \addtogroup NoiseNormalization Noise Normalization
    00942     Estimate noise with intensity-dependent variance and transform it into additive Gaussian noise.
    00943 */
    00944 //@{
    00945
    00946 /********************************************************/
    00947 /*                                                      */
    00948 /*                noiseVarianceEstimation               */
    00949 /*                                                      */
    00950 /********************************************************/
    00951
    00952 /** \brief Determine the noise variance as a function of the image intensity.
    00953
    00954     This operator applies an algorithm described in
    00955
    00956     W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
    00957     Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
    00958     Lecture Notes in Earth Science, Berlin: Springer, 1999
    00959
    00960     in order to estimate the noise variance as a function of the image intensity in a robust way,
    00961     i.e. so that intensity changes due to edges do not bias the estimate. The source value type
    00962     (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
    00963     The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible
    00964     from two <tt>double</tt> values. The following options can be set via the \a options object
    00965     (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
    00966
    00967     <tt>useGradient</tt>, <tt>windowRadius</tt>, <tt>noiseEstimationQuantile</tt>, <tt>noiseVarianceInitialGuess</tt>
    00968
    00969     <b> Declarations:</b>
    00970
    00971     pass arguments explicitly:
    00972     \code
    00973     namespace vigra {
    00974         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00975         void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    00976                                      BackInsertable & result,
    00977                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    00978     }
    00979     \endcode
    00980
    00981     use argument objects in conjunction with \ref ArgumentObjectFactories :
    00982     \code
    00983     namespace vigra {
    00984         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    00985         void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    00986                                      BackInsertable & result,
    00987                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    00988     }
    00989     \endcode
    00990
    00991     <b> Usage:</b>
    00992
    00993         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    00994     Namespace: vigra
    00995
    00996     \code
    00997     vigra::BImage src(w,h);
    00998     std::vector<vigra::TinyVector<double, 2> > result;
    00999
    01000     ...
    01001     vigra::noiseVarianceEstimation(srcImageRange(src), result,
    01002                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0));
    01003
    01004     // print the intensity / variance pairs found
    01005     for(int k=0; k<result.size(); ++k)
    01006         std::cout << "Intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
    01007     \endcode
    01008
    01009     <b> Required Interface:</b>
    01010
    01011     \code
    01012     SrcIterator upperleft, lowerright;
    01013     SrcAccessor src;
    01014
    01015     typedef SrcAccessor::value_type SrcType;
    01016     typedef NumericTraits<SrcType>::isScalar isScalar;
    01017     assert(isScalar::asBool == true);
    01018
    01019     double value = src(uperleft);
    01020
    01021     BackInsertable result;
    01022     typedef BackInsertable::value_type ResultType;
    01023     double intensity, variance;
    01024     result.push_back(ResultType(intensity, variance));
    01025     \endcode
    01026 */
    01027 doxygen_overloaded_function(template <...> void noiseVarianceEstimation)
    01028
    01029 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01030 inline
    01031 void noiseVarianceEstimation(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01032                            BackInsertable & result,
    01033                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01034 {
    01035     typedef typename BackInsertable::value_type ResultType;
    01036     typedef typename SrcAccessor::value_type SrcType;
    01037     typedef typename NumericTraits<SrcType>::isScalar isScalar;
    01038
    01039     VIGRA_STATIC_ASSERT((
    01040         noiseVarianceEstimation_can_only_work_on_scalar_images<(isScalar::asBool)>));
    01041
    01042     detail::noiseVarianceEstimationImpl(sul, slr, src, result, options);
    01043 }
    01044
    01045 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01046 inline
    01047 void noiseVarianceEstimation(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01048                            BackInsertable & result,
    01049                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01050 {
    01051     noiseVarianceEstimation(src.first, src.second, src.third, result, options);
    01052 }
    01053
    01054 /********************************************************/
    01055 /*                                                      */
    01056 /*                noiseVarianceClustering               */
    01057 /*                                                      */
    01058 /********************************************************/
    01059
    01060 /** \brief Determine the noise variance as a function of the image intensity and cluster the results.
    01061
    01062     This operator first calls \ref noiseVarianceEstimation() to obtain a sequence of intensity/variance pairs,
    01063     which are then clustered using the median cut algorithm. Then the cluster centers (i.e. average variance vs.
    01064     average intensity) are determined and returned in the \a result sequence.
    01065
    01066     In addition to the options valid for \ref noiseVarianceEstimation(), the following options can be set via
    01067     the \a options object (see \ref vigra::NoiseNormalizationOptions for details):<br><br>
    01068
    01069     <tt>clusterCount</tt>, <tt>averagingQuantile</tt>
    01070
    01071     <b> Declarations:</b>
    01072
    01073     pass arguments explicitly:
    01074     \code
    01075     namespace vigra {
    01076         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01077         void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01078                                 BackInsertable & result,
    01079                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01080     }
    01081     \endcode
    01082
    01083     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01084     \code
    01085     namespace vigra {
    01086         template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01087         void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01088                                 BackInsertable & result,
    01089                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01090     }
    01091     \endcode
    01092
    01093     <b> Usage:</b>
    01094
    01095         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01096     Namespace: vigra
    01097
    01098     \code
    01099     vigra::BImage src(w,h);
    01100     std::vector<vigra::TinyVector<double, 2> > result;
    01101
    01102     ...
    01103     vigra::noiseVarianceClustering(srcImageRange(src), result,
    01104                                   vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01105                                   clusterCount(15));
    01106
    01107     // print the intensity / variance pairs representing the cluster centers
    01108     for(int k=0; k<result.size(); ++k)
    01109         std::cout << "Cluster: " << k << ", intensity: " << result[k][0] << ", estimated variance: " << result[k][1] << std::endl;
    01110     \endcode
    01111
    01112     <b> Required Interface:</b>
    01113
    01114     same as \ref noiseVarianceEstimation()
    01115 */
    01116 doxygen_overloaded_function(template <...> void noiseVarianceClustering)
    01117
    01118 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01119 inline
    01120 void noiseVarianceClustering(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01121                            BackInsertable & result,
    01122                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01123 {
    01124     ArrayVector<TinyVector<double, 2> > variance;
    01125     noiseVarianceEstimation(sul, slr, src, variance, options);
    01126     detail::noiseVarianceClusteringImpl(variance, result, options.cluster_count, options.averaging_quantile);
    01127 }
    01128
    01129 template <class SrcIterator, class SrcAccessor, class BackInsertable>
    01130 inline
    01131 void noiseVarianceClustering(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01132                            BackInsertable & result,
    01133                            NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01134 {
    01135     noiseVarianceClustering(src.first, src.second, src.third, result, options);
    01136 }
    01137
    01138 /********************************************************/
    01139 /*                                                      */
    01140 /*             nonparametricNoiseNormalization          */
    01141 /*                                                      */
    01142 /********************************************************/
    01143
    01144 /** \brief Noise normalization by means of an estimated non-parametric noise model.
    01145
    01146     The original image is assumed to be corrupted by noise whose variance depends on the intensity in an unknown way.
    01147     The present functions first calls \ref noiseVarianceClustering() to obtain a sequence of intensity/variance pairs
    01148     (cluster centers) which estimate this dependency. The cluster centers are connected into a piecewise linear
    01149     function which is the inverted according to the formula derived in
    01150
    01151     W. F&ouml;rstner: <i>"Image Preprocessing for Feature Extraction in Digital Intensity, Color and Range Images"</i>,
    01152     Proc. Summer School on Data Analysis and the Statistical Foundations of Geomatics,
    01153     Lecture Notes in Earth Science, Berlin: Springer, 1999
    01154
    01155     The inverted formula defines a pixel-wise intensity transformation whose application turns the original image
    01156     into one that is corrupted by additive Gaussian noise with unit variance. Most subsequent algorithms will be able
    01157     to handle this type of noise much better than the original noise.
    01158
    01159     RGB and other multiband images will be processed one band at a time. The function returns <tt>true</tt> on success.
    01160     Noise normalization will fail if the original image does not contain sufficiently homogeneous regions that
    01161     allow robust estimation of the noise variance.
    01162
    01163     The \a options object may use all options described in \ref vigra::NoiseNormalizationOptions.
    01164
    01165     <b> Declarations:</b>
    01166
    01167     pass arguments explicitly:
    01168     \code
    01169     namespace vigra {
    01170         template <class SrcIterator, class SrcAccessor,
    01171                   class DestIterator, class DestAccessor>
    01172         bool nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01173                                              DestIterator dul, DestAccessor dest,
    01174                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01175     }
    01176     \endcode
    01177
    01178     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01179     \code
    01180     namespace vigra {
    01181         template <class SrcIterator, class SrcAccessor,
    01182                   class DestIterator, class DestAccessor>
    01183         bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01184                                              pair<DestIterator, DestAccessor> dest,
    01185                                              NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01186     }
    01187     \endcode
    01188
    01189     <b> Usage:</b>
    01190
    01191         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01192     Namespace: vigra
    01193
    01194     \code
    01195     vigra::BRGBImage src(w,h), dest(w, h);
    01196
    01197     ...
    01198     vigra::nonparametricNoiseNormalization(srcImageRange(src), destImage(dest),
    01199                                            vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01200                                            clusterCount(15));
    01201     \endcode
    01202
    01203     <b> Required Interface:</b>
    01204
    01205     same as \ref noiseVarianceEstimation()
    01206 */
    01207 doxygen_overloaded_function(template <...> bool nonparametricNoiseNormalization)
    01208
    01209 template <class SrcIterator, class SrcAccessor,
    01210           class DestIterator, class DestAccessor>
    01211 inline bool
    01212 nonparametricNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01213                                 DestIterator dul, DestAccessor dest,
    01214                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01215 {
    01216     typedef typename SrcAccessor::value_type SrcType;
    01217
    01218     return detail::nonparametricNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01219                                          typename NumericTraits<SrcType>::isScalar());
    01220 }
    01221
    01222 template <class SrcIterator, class SrcAccessor,
    01223           class DestIterator, class DestAccessor>
    01224 inline
    01225 bool nonparametricNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01226                                      pair<DestIterator, DestAccessor> dest,
    01227                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01228 {
    01229     return nonparametricNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01230 }
    01231
    01232 /********************************************************/
    01233 /*                                                      */
    01234 /*               quadraticNoiseNormalization            */
    01235 /*                                                      */
    01236 /********************************************************/
    01237
    01238 /** \brief Noise normalization by means of an estimated quadratic noise model.
    01239
    01240     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01241     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01242     quadratic function rather than a piecewise linear function. If the quadratic model is applicable, it leads
    01243     to a somewhat smoother transformation.
    01244
    01245     <b> Declarations:</b>
    01246
    01247     pass arguments explicitly:
    01248     \code
    01249     namespace vigra {
    01250         template <class SrcIterator, class SrcAccessor,
    01251                   class DestIterator, class DestAccessor>
    01252         bool quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01253                                          DestIterator dul, DestAccessor dest,
    01254                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01255     }
    01256     \endcode
    01257
    01258     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01259     \code
    01260     namespace vigra {
    01261         template <class SrcIterator, class SrcAccessor,
    01262                   class DestIterator, class DestAccessor>
    01263         bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01264                                          pair<DestIterator, DestAccessor> dest,
    01265                                          NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01266     }
    01267     \endcode
    01268
    01269     <b> Usage:</b>
    01270
    01271         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01272     Namespace: vigra
    01273
    01274     \code
    01275     vigra::BRGBImage src(w,h), dest(w, h);
    01276
    01277     ...
    01278     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01279                                        vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01280                                        clusterCount(15));
    01281     \endcode
    01282
    01283     <b> Required Interface:</b>
    01284
    01285     same as \ref noiseVarianceEstimation()
    01286 */
    01287 doxygen_overloaded_function(template <...> bool quadraticNoiseNormalization)
    01288
    01289 template <class SrcIterator, class SrcAccessor,
    01290           class DestIterator, class DestAccessor>
    01291 inline bool
    01292 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01293                                 DestIterator dul, DestAccessor dest,
    01294                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01295 {
    01296     typedef typename SrcAccessor::value_type SrcType;
    01297
    01298     return detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01299                                          typename NumericTraits<SrcType>::isScalar());
    01300 }
    01301
    01302 template <class SrcIterator, class SrcAccessor,
    01303           class DestIterator, class DestAccessor>
    01304 inline
    01305 bool quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01306                                      pair<DestIterator, DestAccessor> dest,
    01307                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01308 {
    01309     return quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01310 }
    01311
    01312 /********************************************************/
    01313 /*                                                      */
    01314 /*               quadraticNoiseNormalization            */
    01315 /*                                                      */
    01316 /********************************************************/
    01317
    01318 /** \brief Noise normalization by means of a given quadratic noise model.
    01319
    01320     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01321     functional dependency of the noise variance from the intensity is given (by a quadratic function)
    01322     rather than estimated:
    01323
    01324     \code
    01325     variance = a0 + a1 * intensity + a2 * sq(intensity)
    01326     \endcode
    01327
    01328     <b> Declarations:</b>
    01329
    01330     pass arguments explicitly:
    01331     \code
    01332     namespace vigra {
    01333         template <class SrcIterator, class SrcAccessor,
    01334                   class DestIterator, class DestAccessor>
    01335         void quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01336                                          DestIterator dul, DestAccessor dest,
    01337                                          double a0, double a1, double a2);
    01338     }
    01339     \endcode
    01340
    01341     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01342     \code
    01343     namespace vigra {
    01344         template <class SrcIterator, class SrcAccessor,
    01345                   class DestIterator, class DestAccessor>
    01346         void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01347                                         pair<DestIterator, DestAccessor> dest,
    01348                                         double a0, double a1, double a2);
    01349     }
    01350     \endcode
    01351
    01352     <b> Usage:</b>
    01353
    01354         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01355     Namespace: vigra
    01356
    01357     \code
    01358     vigra::BRGBImage src(w,h), dest(w, h);
    01359
    01360     ...
    01361     vigra::quadraticNoiseNormalization(srcImageRange(src), destImage(dest),
    01362                                        100, 0.02, 1e-6);
    01363     \endcode
    01364
    01365     <b> Required Interface:</b>
    01366
    01367     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01368     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01369     or a vector whose elements are assignable from <tt>double</tt>.
    01370 */
    01371 doxygen_overloaded_function(template <...> void quadraticNoiseNormalization)
    01372
    01373 template <class SrcIterator, class SrcAccessor,
    01374           class DestIterator, class DestAccessor>
    01375 inline void
    01376 quadraticNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01377                             DestIterator dul, DestAccessor dest,
    01378                             double a0, double a1, double a2)
    01379 {
    01380     typedef typename SrcAccessor::value_type SrcType;
    01381
    01382     detail::quadraticNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1, a2,
    01383                                          typename NumericTraits<SrcType>::isScalar());
    01384 }
    01385
    01386 template <class SrcIterator, class SrcAccessor,
    01387           class DestIterator, class DestAccessor>
    01388 inline
    01389 void quadraticNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01390                                  pair<DestIterator, DestAccessor> dest,
    01391                                  double a0, double a1, double a2)
    01392 {
    01393     quadraticNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1, a2);
    01394 }
    01395
    01396 /********************************************************/
    01397 /*                                                      */
    01398 /*                linearNoiseNormalization              */
    01399 /*                                                      */
    01400 /********************************************************/
    01401
    01402 /** \brief Noise normalization by means of an estimated linear noise model.
    01403
    01404     This function works in the same way as \ref nonparametricNoiseNormalization() with the exception of the
    01405     model for the dependency between intensity and noise variance: it assumes that this dependency is a
    01406     linear function rather than a piecewise linear function. If the linear model is applicable, it leads
    01407     to a very simple transformation which is similar to the familiar gamma correction.
    01408
    01409     <b> Declarations:</b>
    01410
    01411     pass arguments explicitly:
    01412     \code
    01413     namespace vigra {
    01414         template <class SrcIterator, class SrcAccessor,
    01415                 class DestIterator, class DestAccessor>
    01416         bool linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01417                                       DestIterator dul, DestAccessor dest,
    01418                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01419     }
    01420     \endcode
    01421
    01422     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01423     \code
    01424     namespace vigra {
    01425         template <class SrcIterator, class SrcAccessor,
    01426                   class DestIterator, class DestAccessor>
    01427         bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01428                                       pair<DestIterator, DestAccessor> dest,
    01429                                       NoiseNormalizationOptions const & options = NoiseNormalizationOptions());
    01430     }
    01431     \endcode
    01432
    01433     <b> Usage:</b>
    01434
    01435         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01436     Namespace: vigra
    01437
    01438     \code
    01439     vigra::BRGBImage src(w,h), dest(w, h);
    01440
    01441     ...
    01442     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01443                                     vigra::NoiseNormalizationOptions().windowRadius(9).noiseVarianceInitialGuess(25.0).
    01444                                     clusterCount(15));
    01445     \endcode
    01446
    01447     <b> Required Interface:</b>
    01448
    01449     same as \ref noiseVarianceEstimation()
    01450 */
    01451 doxygen_overloaded_function(template <...> bool linearNoiseNormalization)
    01452
    01453 template <class SrcIterator, class SrcAccessor,
    01454           class DestIterator, class DestAccessor>
    01455 inline bool
    01456 linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01457                                 DestIterator dul, DestAccessor dest,
    01458                                 NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01459 {
    01460     typedef typename SrcAccessor::value_type SrcType;
    01461
    01462     return detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, options,
    01463                                          typename NumericTraits<SrcType>::isScalar());
    01464 }
    01465
    01466 template <class SrcIterator, class SrcAccessor,
    01467           class DestIterator, class DestAccessor>
    01468 inline
    01469 bool linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01470                                      pair<DestIterator, DestAccessor> dest,
    01471                                      NoiseNormalizationOptions const & options = NoiseNormalizationOptions())
    01472 {
    01473     return linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, options);
    01474 }
    01475
    01476 /********************************************************/
    01477 /*                                                      */
    01478 /*                 linearNoiseNormalization             */
    01479 /*                                                      */
    01480 /********************************************************/
    01481
    01482 /** \brief Noise normalization by means of a given linear noise model.
    01483
    01484     This function works similar to \ref nonparametricNoiseNormalization() with the exception that the
    01485     functional dependency of the noise variance from the intensity is given (as a linear function)
    01486     rather than estimated:
    01487
    01488     \code
    01489     variance = a0 + a1 * intensity
    01490     \endcode
    01491
    01492     <b> Declarations:</b>
    01493
    01494     pass arguments explicitly:
    01495     \code
    01496     namespace vigra {
    01497         template <class SrcIterator, class SrcAccessor,
    01498                   class DestIterator, class DestAccessor>
    01499         void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01500                                       DestIterator dul, DestAccessor dest,
    01501                                       double a0, double a1);
    01502     }
    01503     \endcode
    01504
    01505     use argument objects in conjunction with \ref ArgumentObjectFactories :
    01506     \code
    01507     namespace vigra {
    01508         template <class SrcIterator, class SrcAccessor,
    01509                   class DestIterator, class DestAccessor>
    01510         void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01511                                       pair<DestIterator, DestAccessor> dest,
    01512                                       double a0, double a1);
    01513     }
    01514     \endcode
    01515
    01516     <b> Usage:</b>
    01517
    01518         <b>\#include</b> <<a href="noise__normalization_8hxx-source.html">vigra/noise_normalization.hxx</a>><br>
    01519     Namespace: vigra
    01520
    01521     \code
    01522     vigra::BRGBImage src(w,h), dest(w, h);
    01523
    01524     ...
    01525     vigra::linearNoiseNormalization(srcImageRange(src), destImage(dest),
    01526                                     100, 0.02);
    01527     \endcode
    01528
    01529     <b> Required Interface:</b>
    01530
    01531     The source value type must be convertible to <tt>double</tt> or must be a vector whose elements
    01532     are convertible to <tt>double</tt>. Likewise, the destination type must be assignable from <tt>double</tt>
    01533     or a vector whose elements are assignable from <tt>double</tt>.
    01534 */
    01535 doxygen_overloaded_function(template <...> void linearNoiseNormalization)
    01536
    01537 template <class SrcIterator, class SrcAccessor,
    01538           class DestIterator, class DestAccessor>
    01539 inline
    01540 void linearNoiseNormalization(SrcIterator sul, SrcIterator slr, SrcAccessor src,
    01541                               DestIterator dul, DestAccessor dest,
    01542                               double a0, double a1)
    01543 {
    01544     typedef typename SrcAccessor::value_type SrcType;
    01545
    01546     detail::linearNoiseNormalizationImpl(sul, slr, src, dul, dest, a0, a1,
    01547                                          typename NumericTraits<SrcType>::isScalar());
    01548 }
    01549
    01550 template <class SrcIterator, class SrcAccessor,
    01551           class DestIterator, class DestAccessor>
    01552 inline
    01553 void linearNoiseNormalization(triple<SrcIterator, SrcIterator, SrcAccessor> src,
    01554                               pair<DestIterator, DestAccessor> dest,
    01555                               double a0, double a1)
    01556 {
    01557     linearNoiseNormalization(src.first, src.second, src.third, dest.first, dest.second, a0, a1);
    01558 }
    01559
    01560 //@}
    01561
    01562 } // namespace vigra
    01563
    01564 #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.7.1 (3 Dec 2010)