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