vigra - vigra/slanted_edge_mtf.hxx Source File

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