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

vigra/slanted_edge_mtf.hxx VIGRA

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)