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

vigra/noise_normalization.hxx

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

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

html generated using doxygen and Python
VIGRA 1.6.0 (13 Aug 2008)