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

vigra/slanted_edge_mtf.hxx

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

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

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)