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

details vigra/slanted_edge_mtf.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 1998-2006 by Ullrich Koethe                  */
00004 /*       Cognitive Systems Group, University of Hamburg, Germany        */
00005 /*                                                                      */
00006 /*    This file is part of the VIGRA computer vision library.           */
00007 /*    ( Version 1.5.0, Dec 07 2006 )                                    */
00008 /*    The VIGRA Website is                                              */
00009 /*        http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/      */
00010 /*    Please direct questions, bug reports, and contributions to        */
00011 /*        koethe@informatik.uni-hamburg.de          or                  */
00012 /*        vigra@kogs1.informatik.uni-hamburg.de                         */
00013 /*                                                                      */
00014 /*    Permission is hereby granted, free of charge, to any person       */
00015 /*    obtaining a copy of this software and associated documentation    */
00016 /*    files (the "Software"), to deal in the Software without           */
00017 /*    restriction, including without limitation the rights to use,      */
00018 /*    copy, modify, merge, publish, distribute, sublicense, and/or      */
00019 /*    sell copies of the Software, and to permit persons to whom the    */
00020 /*    Software is furnished to do so, subject to the following          */
00021 /*    conditions:                                                       */
00022 /*                                                                      */
00023 /*    The above copyright notice and this permission notice shall be    */
00024 /*    included in all copies or substantial portions of the             */
00025 /*    Software.                                                         */
00026 /*                                                                      */
00027 /*    THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND    */
00028 /*    EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   */
00029 /*    OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          */
00030 /*    NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT       */
00031 /*    HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,      */
00032 /*    WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING      */
00033 /*    FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR     */
00034 /*    OTHER DEALINGS IN THE SOFTWARE.                                   */
00035 /*                                                                      */
00036 /************************************************************************/
00037 
00038 
00039 #ifndef VIGRA_SLANTED_EDGE_MTF_HXX
00040 #define VIGRA_SLANTED_EDGE_MTF_HXX
00041 
00042 #include <algorithm>
00043 #include "array_vector.hxx"
00044 #include "basicgeometry.hxx"
00045 #include "edgedetection.hxx"
00046 #include "fftw3.hxx"
00047 #include "functorexpression.hxx"
00048 #include "linear_solve.hxx"
00049 #include "mathutil.hxx"
00050 #include "numerictraits.hxx"
00051 #include "separableconvolution.hxx"
00052 #include "static_assert.hxx"
00053 #include "stdimage.hxx"
00054 #include "transformimage.hxx"
00055 #include "utilities.hxx"
00056 
00057 namespace vigra {
00058 
00059 /** \addtogroup SlantedEdgeMTF Camera MTF Estimation
00060     Determine the magnitude transfer function (MTF) of a camera using the slanted edge method.
00061 */
00062 //@{ 
00063                                     
00064 /********************************************************/
00065 /*                                                      */
00066 /*                  SlantedEdgeMTFOptions               */
00067 /*                                                      */
00068 /********************************************************/
00069 
00070 /** \brief Pass options to one of the \ref slantedEdgeMTF() functions.
00071 
00072     <tt>SlantedEdgeMTFOptions</tt>  is an argument objects that holds various optional
00073     parameters used by the \ref slantedEdgeMTF() functions. If a parameter is not explicitly
00074     set, a suitable default will be used. Changing the defaults is only necessary if you can't 
00075     obtain good input data, but absolutely need an MTF estimate.
00076     
00077     <b> Usage:</b>
00078     
00079         <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
00080     Namespace: vigra
00081     
00082     \code
00083     vigra::BImage src(w,h);
00084     std::vector<vigra::TinyVector<double, 2> > mtf;
00085     
00086     ...
00087     vigra::slantedEdgeMTF(srcImageRange(src), mtf,
00088                           vigra::SlantedEdgeMTFOptions().mtfSmoothingScale(1.0));
00089     
00090     // print the frequency / attenuation pairs found
00091     for(int k=0; k<result.size(); ++k)
00092         std::cout << "frequency: " << mtf[k][0] << ", estimated attenuation: " << mtf[k][1] << std::endl;
00093     \endcode
00094 */
00095 
00096 class SlantedEdgeMTFOptions
00097 {
00098   public:
00099         /** Initialize all options with default values.
00100         */
00101     SlantedEdgeMTFOptions()
00102     : minimum_number_of_lines(20),
00103       desired_edge_width(10),
00104       minimum_edge_width(5),
00105       mtf_smoothing_scale(2.0)
00106     {}
00107 
00108         /** Minimum number of pixels the edge must cross.
00109         
00110             The longer the edge the more accurate the resulting MTF estimate. If you don't have good
00111             data, but absolutely have to compute an MTF, you may force a lower value here.<br>
00112             Default: 20
00113         */
00114     SlantedEdgeMTFOptions & minimumNumberOfLines(unsigned int n)
00115     {
00116         minimum_number_of_lines = n;
00117         return *this;
00118     }
00119 
00120         /** Desired number of pixels perpendicular to the edge.
00121         
00122             The larger the regions to either side of the edge, 
00123             the more accurate the resulting MTF estimate. If you don't have good
00124             data, but absolutely have to compute an MTF, you may force a lower value here.<br>
00125             Default: 10
00126         */
00127     SlantedEdgeMTFOptions & desiredEdgeWidth(unsigned int n)
00128     {
00129         desired_edge_width = n;
00130         return *this;
00131     }
00132 
00133         /** Minimum acceptable number of pixels perpendicular to the edge.
00134         
00135             The larger the regions to either side of the edge, 
00136             the more accurate the resulting MTF estimate. If you don't have good
00137             data, but absolutely have to compute an MTF, you may force a lower value here.<br>
00138             Default: 5
00139         */
00140     SlantedEdgeMTFOptions & minimumEdgeWidth(unsigned int n)
00141     {
00142         minimum_edge_width = n;
00143         return *this;
00144     }
00145 
00146         /** Amount of smoothing of the computed MTF.
00147         
00148             If the datais noisy, so will be the MTF. Thus, some smoothing is useful.<br>
00149             Default: 2.0
00150         */
00151     SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
00152     {
00153         vigra_precondition(scale >= 0.0,
00154             "SlantedEdgeMTFOptions: MTF smoothing scale must not be < 0.");
00155         mtf_smoothing_scale = scale;
00156         return *this;
00157     }
00158 
00159     unsigned int minimum_number_of_lines, desired_edge_width, minimum_edge_width;
00160     double mtf_smoothing_scale;
00161 };
00162 
00163 //@}
00164 
00165 namespace detail {
00166 
00167 struct SortEdgelsByStrength
00168 {
00169     bool operator()(Edgel const & l, Edgel const & r) const
00170     {
00171         return l.strength > r.strength;
00172     }
00173 };
00174 
00175     /* Make sure that the edge runs vertically, intersects the top and bottom border
00176        of the image, and white is on the left.
00177     */
00178 template <class SrcIterator, class SrcAccessor, class DestImage>
00179 unsigned int
00180 prepareSlantedEdgeInput(SrcIterator sul, SrcIterator slr, SrcAccessor src, DestImage & res,
00181                         SlantedEdgeMTFOptions const & options)
00182 {
00183     unsigned int w = slr.x - sul.x;
00184     unsigned int h = slr.y - sul.y;
00185 
00186     // find rough estimate of the edge
00187     ArrayVector<Edgel> edgels;
00188     cannyEdgelList(sul, slr, src, edgels, 2.0);
00189     std::sort(edgels.begin(), edgels.end(), SortEdgelsByStrength());
00190 
00191     double x = 0.0, y = 0.0, x2 = 0.0, y2 = 0.0, xy = 0.0;
00192     unsigned int c = std::min(w, h);
00193 
00194     for(unsigned int k = 0; k < c; ++k)
00195     {
00196         x += edgels[k].x;
00197         y += edgels[k].y;
00198         x2 += sq(edgels[k].x);
00199         xy += edgels[k].x*edgels[k].y;
00200         y2 += sq(edgels[k].y);
00201     }
00202     double xc = x / c, yc = y / c;
00203     x2 = x2 / c - sq(xc);
00204     xy = xy / c - xc*yc;
00205     y2 = y2 / c - sq(yc);
00206     double angle = 0.5*VIGRA_CSTD::atan2(2*xy, x2 - y2);
00207 
00208     DestImage tmp;
00209     // rotate image when slope is less than +-45 degrees
00210     if(VIGRA_CSTD::fabs(angle) < M_PI / 4.0)
00211     {
00212         xc = yc;
00213         yc = w - xc - 1;
00214         std::swap(w, h);
00215         tmp.resize(w, h);
00216         rotateImage(srcIterRange(sul, slr, src), destImage(tmp), 90);
00217         angle += M_PI / 2.0;
00218     }
00219     else
00220     {
00221         tmp.resize(w, h);
00222         copyImage(srcIterRange(sul, slr, src), destImage(tmp));
00223         if(angle < 0.0)
00224             angle += M_PI;
00225     }
00226     // angle is now between pi/4 and 3*pi/4
00227     double slope = VIGRA_CSTD::tan(M_PI/2.0 - angle);
00228     vigra_precondition(slope != 0.0,
00229           "slantedEdgeMTF(): Input edge is not slanted");
00230 
00231     // trim image so that the edge only intersects the top and bottom border
00232     unsigned int minimumNumberOfLines = options.minimum_number_of_lines, //20,
00233                  edgeWidth = options.desired_edge_width, // 10
00234                  minimumEdgeWidth = options.minimum_edge_width; // 5
00235 
00236     int y0, y1;
00237     for(; edgeWidth >= minimumEdgeWidth; --edgeWidth)
00238     {
00239         y0 = int(VIGRA_CSTD::floor((edgeWidth - xc) / slope + yc + 0.5));
00240         y1 = int(VIGRA_CSTD::floor((w - edgeWidth - 1 - xc) / slope + yc + 0.5));
00241         if(slope < 0.0)
00242             std::swap(y0, y1);
00243         if(y1 - y0 >= (int)minimumNumberOfLines)
00244             break;
00245     }
00246 
00247     vigra_precondition(edgeWidth >= minimumEdgeWidth,
00248         "slantedEdgeMTF(): Input edge is too slanted or image is too small");
00249 
00250     y0 = std::max(y0, 0);
00251     y1 = std::min(y1+1, (int)h);
00252 
00253     res.resize(w, y1-y0);
00254 
00255     // ensure that white is on the left
00256     if(tmp(0,0) < tmp(w-1, h-1))
00257     {
00258         rotateImage(srcIterRange(tmp.upperLeft() + Diff2D(0, y0), tmp.upperLeft() + Diff2D(w, y1), tmp.accessor()),
00259                     destImage(res), 180);
00260     }
00261     else
00262     {
00263         copyImage(srcImageRange(tmp), destImage(res));
00264     }
00265     return edgeWidth;
00266 }
00267 
00268 template <class Image>
00269 void slantedEdgeShadingCorrection(Image & i, unsigned int edgeWidth)
00270 {
00271     using namespace functor;
00272 
00273     // after prepareSlantedEdgeInput(), the white region is on the left
00274     // find a plane that approximates the logarithm of the white ROI
00275 
00276     transformImage(srcImageRange(i), destImage(i), log(Arg1() + Param(1.0)));
00277 
00278     unsigned int w = i.width(),
00279                  h = i.height(),
00280                  s = edgeWidth*h;
00281 
00282     Matrix<double> m(3,3), r(3, 1), l(3, 1);
00283     for(unsigned int y = 0; y < h; ++y)
00284     {
00285         for(unsigned int x = 0; x < edgeWidth; ++x)
00286         {
00287             l(0,0) = x;
00288             l(1,0) = y;
00289             l(2,0) = 1.0;
00290             m += outer(l);
00291             r += i(x,y)*l;
00292         }
00293     }
00294 
00295     linearSolve(m, r, l);
00296     double a = l(0,0),
00297            b = l(1,0),
00298            c = l(2,0);
00299 
00300     // subtract the plane and go back to the non-logarithmic representation
00301     for(unsigned int y = 0; y < h; ++y)
00302     {
00303         for(unsigned int x = 0; x < w; ++x)
00304         {
00305             i(x, y) = VIGRA_CSTD::exp(i(x,y) - a*x - b*y - c);
00306         }
00307     }
00308 }
00309 
00310 template <class Image, class BackInsertable>
00311 void slantedEdgeSubpixelShift(Image const & img, BackInsertable & centers, double & angle,
00312                               SlantedEdgeMTFOptions const & options)
00313 {
00314     unsigned int w = img.width();
00315     unsigned int h = img.height();
00316 
00317     Image grad(w, h);
00318     Kernel1D<double> kgrad;
00319     kgrad.initGaussianDerivative(1.0, 1);
00320 
00321     separableConvolveX(srcImageRange(img), destImage(grad), kernel1d(kgrad));
00322 
00323     int desiredEdgeWidth = (int)options.desired_edge_width;
00324     double sy = 0.0, sx = 0.0, syy = 0.0, sxy = 0.0;
00325     for(unsigned int y = 0; y < h; ++y)
00326     {
00327         double a = 0.0,
00328                b = 0.0;
00329         for(unsigned int x = 0; x < w; ++x)
00330         {
00331             a += x*grad(x,y);
00332             b += grad(x,y);
00333         }
00334         int c = int(a / b);
00335         // c is biased because the numbers of black and white pixels differ
00336         // repeat the analysis with a symmetric window around the edge
00337         a = 0.0;
00338         b = 0.0;
00339         int ew = desiredEdgeWidth;
00340         if(c-desiredEdgeWidth < 0)
00341             ew = c;
00342         if(c + ew + 1 >= (int)w)
00343             ew = w - c - 1;
00344         for(int x = c-ew; x < c+ew+1; ++x)
00345         {
00346             a += x*grad(x,y);
00347             b += grad(x,y);
00348         }
00349         double sc = a / b;
00350         sy += y;
00351         sx += sc;
00352         syy += sq(y);
00353         sxy += sc*y;
00354     }
00355     // fit a line to the subpixel locations
00356     double a = (h * sxy - sx * sy) / (h * syy - sq(sy));
00357     double b = (sx * syy - sxy * sy) / (h * syy - sq(sy));
00358 
00359     // compute the regularized subpixel values of the edge location
00360     angle = VIGRA_CSTD::atan(a);
00361     for(unsigned int y = 0; y < h; ++y)
00362     {
00363         centers.push_back(a * y + b);
00364     }
00365 }
00366 
00367 template <class Image, class Vector>
00368 void slantedEdgeCreateOversampledLine(Image const & img, Vector const & centers,
00369                                       Image & result)
00370 {
00371     unsigned int w = img.width();
00372     unsigned int h = img.height();
00373     unsigned int w2 = std::min(std::min(int(centers[0]), int(centers[h-1])),
00374                                std::min(int(w - centers[0]) - 1, int(w - centers[h-1]) - 1));
00375     unsigned int ww = 8*w2;
00376 
00377     Image r(ww, 1), s(ww, 1);
00378     for(unsigned int y = 0; y < h; ++y)
00379     {
00380         int x0 = int(centers[y]) - w2;
00381         int x1 = int((VIGRA_CSTD::ceil(centers[y]) - centers[y])*4);
00382         for(; x1 < (int)ww; x1 += 4)
00383         {
00384             r(x1, 0) += img(x0, y);
00385             ++s(x1, 0);
00386             ++x0;
00387         }
00388     }
00389 
00390     for(unsigned int x = 0; x < ww; ++x)
00391     {
00392         vigra_precondition(s(x,0) > 0.0,
00393            "slantedEdgeMTF(): Input edge is not slanted enough");
00394         r(x,0) /= s(x,0);
00395     }
00396 
00397     result.resize(ww-1, 1);
00398     for(unsigned int x = 0; x < ww-1; ++x)
00399     {
00400         result(x,0) = r(x+1,0) - r(x,0);
00401     }
00402 }
00403 
00404 template <class Image, class BackInsertable>
00405 void slantedEdgeMTFImpl(Image const & i, BackInsertable & mtf, double angle,
00406                         SlantedEdgeMTFOptions const & options)
00407 {
00408     unsigned int w = i.width();
00409     unsigned int h = i.height();
00410 
00411     double slantCorrection = VIGRA_CSTD::cos(angle);
00412     int desiredEdgeWidth = 4*options.desired_edge_width;
00413 
00414     Image magnitude;
00415 
00416     if(w - 2*desiredEdgeWidth < 64)
00417     {
00418         FFTWComplexImage otf(w, h);
00419         fourierTransform(srcImageRange(i), destImage(otf));
00420 
00421         magnitude.resize(w, h);
00422         for(unsigned int y = 0; y < h; ++y)
00423         {
00424             for(unsigned int x = 0; x < w; ++x)
00425             {
00426                 magnitude(x, y) = norm(otf(x, y));
00427             }
00428         }
00429     }
00430     else
00431     {
00432         w -= 2*desiredEdgeWidth;
00433         FFTWComplexImage otf(w, h);
00434         fourierTransform(srcImageRange(i, Rect2D(Point2D(desiredEdgeWidth, 0), Size2D(w, h))),
00435                          destImage(otf));
00436 
00437         // create an image where the edge is skipped - presumably it only contains the
00438         // noise which can then be subtracted
00439         Image noise(w,h);
00440         copyImage(srcImageRange(i, Rect2D(Point2D(0,0), Size2D(w/2, h))),
00441                   destImage(noise));
00442         copyImage(srcImageRange(i, Rect2D(Point2D(i.width()-w/2, 0), Size2D(w/2, h))),
00443                   destImage(noise, Point2D(w/2, 0)));
00444         FFTWComplexImage fnoise(w, h);
00445         fourierTransform(srcImageRange(noise), destImage(fnoise));
00446 
00447         // subtract the noise power spectrum from the total power spectrum
00448         magnitude.resize(w, h);
00449         for(unsigned int y = 0; y < h; ++y)
00450         {
00451             for(unsigned int x = 0; x < w; ++x)
00452             {
00453                 magnitude(x, y) = VIGRA_CSTD::sqrt(std::max(0.0, squaredNorm(otf(x, y))-squaredNorm(fnoise(x, y))));
00454             }
00455         }
00456     }
00457 
00458     Kernel1D<double> gauss;
00459     gauss.initGaussian(options.mtf_smoothing_scale);
00460     Image smooth(w,h);
00461     separableConvolveX(srcImageRange(magnitude), destImage(smooth), kernel1d(gauss));
00462 
00463     unsigned int ww = w/4;
00464     double maxVal = smooth(0,0),
00465            minVal = maxVal;
00466     for(unsigned int k = 1; k < ww; ++k)
00467     {
00468         if(smooth(k,0) >= 0.0 && smooth(k,0) < minVal)
00469             minVal = smooth(k,0);
00470     }
00471     double norm = maxVal-minVal;
00472 
00473     typedef typename BackInsertable::value_type Result;
00474     mtf.push_back(Result(0.0, 1.0));
00475     for(unsigned int k = 1; k < ww; ++k)
00476     {
00477         double n = (smooth(k,0) - minVal)/norm*sq(M_PI*k/w/VIGRA_CSTD::sin(M_PI*k/w));
00478         double xx = 4.0*k/w/slantCorrection;
00479         if(n < 0.0 || xx > 1.0)
00480             break;
00481         mtf.push_back(Result(xx, n));
00482     }
00483 }
00484 
00485 } // namespace detail
00486 
00487 /** \addtogroup SlantedEdgeMTF Camera MTF Estimation
00488     Determine the magnitude transfer function (MTF) of a camera using the slanted edge method.
00489 */
00490 //@{ 
00491                                     
00492 /********************************************************/
00493 /*                                                      */
00494 /*                     slantedEdgeMTF                   */
00495 /*                                                      */
00496 /********************************************************/
00497 
00498 /** \brief Determine the magnitude transfer function of the camera.
00499 
00500     This operator estimates the magnitude transfer function (MTF) of a camera by means of the 
00501     slanted edge method described in:
00502     
00503     ISO Standard No. 12233: <i>"Photography - Electronic still picture cameras - Resolution measurements"</i>, 2000
00504     
00505     The input must be an image that contains a single step edge with bright pixels on one side and dark pixels on 
00506     the other. However, the intensity values must be neither saturated nor zero. The algorithms computes the MTF
00507     from the Fourier transform of the edge's derivative. Thus, if the actual MTF is unisotropic, the estimated 
00508     MTF does actually only apply in the direction perpendicular to the edge - several edges at different 
00509     orientations are required to estimate an unisotropic MTF.
00510     
00511     The algorithm returns a sequence of frequency / attenuation pairs. The frequency axis is normalized so that the
00512     Nyquist frequency of the original image is 0.5. Since the edge's derivative is computed with subpixel accuracy,
00513     the attenuation can usually be computed for frequencies significantly above the Nyquist frequency as well. The 
00514     MTF estimate ends at either the first zero crossing of the MTF or at frequency 1, whichever comes earlier.
00515     
00516     The present implementation improves the original slanted edge algorithm according to ISO 12233 in a number of
00517     ways:
00518     
00519     <ul>
00520     <li> The edge is not required to run nearly vertically or horizontally (i.e. with a slant of approximately 5 degrees).
00521          The algorithm will automatically compute the edge's actual angle and adjust estimates accordingly. 
00522          However, it is still necessary for the edge to be somewhat slanted, because subpixel-accurate estimation 
00523          of the derivative is impossible otherwise (i.e. the edge position perpendicular to the edge direction must 
00524          differ by at least 1 pixel between the two ends of the edge). 
00525          
00526     <li> Our implementation uses a more accurate subpixel derivative algrithm. In addition, we first perform a shading 
00527          correction in order to reduce possible derivative bias due to nonuniform illumination.
00528 
00529     <li> If the input image is large enough (i.e. there are at least 20 pixels on either side of the edge over
00530          the edge's entire length), our algorithm attempts to subtract the estimated noise power spectrum
00531          from the estimated MTF.
00532     </ul>
00533     
00534     The source value type (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
00535     The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible 
00536     from two <tt>double</tt> values. Algorithm options can be set via the \a options object 
00537     (see \ref vigra::NoiseNormalizationOptions for details).
00538     
00539     <b> Declarations:</b>
00540     
00541     pass arguments explicitly:
00542     \code
00543     namespace vigra {
00544         template <class SrcIterator, class SrcAccessor, class BackInsertable>
00545         void
00546         slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
00547                     SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions());
00548     }
00549     \endcode
00550     
00551     use argument objects in conjunction with \ref ArgumentObjectFactories:
00552     \code
00553     namespace vigra {
00554         template <class SrcIterator, class SrcAccessor, class BackInsertable>
00555         void
00556         slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
00557                        SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
00558     }
00559     \endcode
00560     
00561     <b> Usage:</b>
00562     
00563         <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
00564     Namespace: vigra
00565     
00566     \code
00567     vigra::BImage src(w,h);
00568     std::vector<vigra::TinyVector<double, 2> > mtf;
00569     
00570     ...
00571     vigra::slantedEdgeMTF(srcImageRange(src), mtf);
00572     
00573     // print the frequency / attenuation pairs found
00574     for(int k=0; k<result.size(); ++k)
00575         std::cout << "frequency: " << mtf[k][0] << ", estimated attenuation: " << mtf[k][1] << std::endl;
00576     \endcode
00577 
00578     <b> Required Interface:</b>
00579     
00580     \code
00581     SrcIterator upperleft, lowerright;
00582     SrcAccessor src;
00583     
00584     typedef SrcAccessor::value_type SrcType;
00585     typedef NumericTraits<SrcType>::isScalar isScalar;
00586     assert(isScalar::asBool == true);
00587     
00588     double value = src(uperleft);
00589     
00590     BackInsertable result;
00591     typedef BackInsertable::value_type ResultType;    
00592     double intensity, variance;
00593     result.push_back(ResultType(intensity, variance));
00594     \endcode
00595 */
00596 template <class SrcIterator, class SrcAccessor, class BackInsertable>
00597 void
00598 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
00599                SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
00600 {
00601     DImage preparedInput;
00602     unsigned int edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options);
00603     detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth);
00604 
00605     ArrayVector<double> centers;
00606     double angle = 0.0;
00607     detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options);
00608 
00609     DImage oversampledLine;
00610     detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine);
00611 
00612     detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options);
00613 }
00614 
00615 template <class SrcIterator, class SrcAccessor, class BackInsertable>
00616 inline void
00617 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
00618                SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
00619 {
00620     slantedEdgeMTF(src.first, src.second, src.third, mtf, options);
00621 }
00622 
00623 /********************************************************/
00624 /*                                                      */
00625 /*                     mtfFitGaussian                   */
00626 /*                                                      */
00627 /********************************************************/
00628 
00629 /** \brief Fit a Gaussian function to a given MTF.
00630 
00631     This function expects a squence of frequency / attenuation pairs as produced by \ref slantedEdgeMTF()
00632     and finds the best fitting Gaussian point spread function (Gaussian functions are good approximations 
00633     of the PSF of many real cameras). It returns the standard deviation (scale) of this function. The algorithm
00634     computes the standard deviation by means of a linear least square on the logarithm of the MTF, i.e.
00635     an algebraic fit rather than a Euclidean fit - thus, the resulting Gaussian may not be the one that 
00636     intuitively fits the data optimally.
00637     
00638     <b> Declaration:</b>
00639     
00640     \code
00641     namespace vigra {
00642         template <class Vector>
00643         double mtfFitGaussian(Vector const & mtf);
00644     }
00645     \endcode
00646     
00647     <b> Usage:</b>
00648     
00649         <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
00650     Namespace: vigra
00651     
00652     \code
00653     vigra::BImage src(w,h);
00654     std::vector<vigra::TinyVector<double, 2> > mtf;
00655     
00656     ...
00657     vigra::slantedEdgeMTF(srcImageRange(src), mtf);
00658     double scale = vigra::mtfFitGaussian(mtf)
00659     
00660     std::cout << "The camera PSF is approximately a Gaussian at scale " << scale << std::endl;
00661     \endcode
00662 
00663     <b> Required Interface:</b>
00664     
00665     \code
00666     Vector mtf;
00667     int numberOfMeasurements = mtf.size()
00668     
00669     double frequency = mtf[0][0];
00670     double attenuation = mtf[0][1];
00671     \endcode
00672 */
00673 template <class Vector>
00674 double mtfFitGaussian(Vector const & mtf)
00675 {
00676     double minVal = mtf[0][1];
00677     for(unsigned int k = 1; k < mtf.size(); ++k)
00678     {
00679         if(mtf[k][1] < minVal)
00680             minVal = mtf[k][1];
00681     }
00682     double x2 = 0.0,
00683            xy = 0.0;
00684     for(unsigned int k = 1; k < mtf.size(); ++k)
00685     {
00686         if(mtf[k][1] <= 0.0)
00687             break;
00688         double x = mtf[k][0],
00689                y = VIGRA_CSTD::sqrt(-VIGRA_CSTD::log(mtf[k][1])/2.0)/M_PI;
00690         x2 += vigra::sq(x);
00691         xy += x*y;
00692         if(mtf[k][1] == minVal)
00693             break;
00694     }
00695     return xy / x2;
00696 }
00697 
00698 //@}
00699 
00700 } // namespace vigra
00701 
00702 #endif // VIGRA_SLANTED_EDGE_MTF_HXX

© Ullrich Köthe (koethe@informatik.uni-hamburg.de)
Cognitive Systems Group, University of Hamburg, Germany

html generated using doxygen and Python
VIGRA 1.5.0 (7 Dec 2006)