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

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

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

html generated using doxygen and Python
vigra 1.7.0 (Thu Aug 25 2011)