[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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.6.0, Aug 13 2008 ) */ 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 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00012 /* vigra@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 doxygen_overloaded_function(template <...> void slantedEdgeMTF) 00597 00598 template <class SrcIterator, class SrcAccessor, class BackInsertable> 00599 void 00600 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf, 00601 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions()) 00602 { 00603 DImage preparedInput; 00604 unsigned int edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options); 00605 detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth); 00606 00607 ArrayVector<double> centers; 00608 double angle = 0.0; 00609 detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options); 00610 00611 DImage oversampledLine; 00612 detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine); 00613 00614 detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options); 00615 } 00616 00617 template <class SrcIterator, class SrcAccessor, class BackInsertable> 00618 inline void 00619 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf, 00620 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions()) 00621 { 00622 slantedEdgeMTF(src.first, src.second, src.third, mtf, options); 00623 } 00624 00625 /********************************************************/ 00626 /* */ 00627 /* mtfFitGaussian */ 00628 /* */ 00629 /********************************************************/ 00630 00631 /** \brief Fit a Gaussian function to a given MTF. 00632 00633 This function expects a squence of frequency / attenuation pairs as produced by \ref slantedEdgeMTF() 00634 and finds the best fitting Gaussian point spread function (Gaussian functions are good approximations 00635 of the PSF of many real cameras). It returns the standard deviation (scale) of this function. The algorithm 00636 computes the standard deviation by means of a linear least square on the logarithm of the MTF, i.e. 00637 an algebraic fit rather than a Euclidean fit - thus, the resulting Gaussian may not be the one that 00638 intuitively fits the data optimally. 00639 00640 <b> Declaration:</b> 00641 00642 \code 00643 namespace vigra { 00644 template <class Vector> 00645 double mtfFitGaussian(Vector const & mtf); 00646 } 00647 \endcode 00648 00649 <b> Usage:</b> 00650 00651 <b>\#include</b> <<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>><br> 00652 Namespace: vigra 00653 00654 \code 00655 vigra::BImage src(w,h); 00656 std::vector<vigra::TinyVector<double, 2> > mtf; 00657 00658 ... 00659 vigra::slantedEdgeMTF(srcImageRange(src), mtf); 00660 double scale = vigra::mtfFitGaussian(mtf) 00661 00662 std::cout << "The camera PSF is approximately a Gaussian at scale " << scale << std::endl; 00663 \endcode 00664 00665 <b> Required Interface:</b> 00666 00667 \code 00668 Vector mtf; 00669 int numberOfMeasurements = mtf.size() 00670 00671 double frequency = mtf[0][0]; 00672 double attenuation = mtf[0][1]; 00673 \endcode 00674 */ 00675 template <class Vector> 00676 double mtfFitGaussian(Vector const & mtf) 00677 { 00678 double minVal = mtf[0][1]; 00679 for(unsigned int k = 1; k < mtf.size(); ++k) 00680 { 00681 if(mtf[k][1] < minVal) 00682 minVal = mtf[k][1]; 00683 } 00684 double x2 = 0.0, 00685 xy = 0.0; 00686 for(unsigned int k = 1; k < mtf.size(); ++k) 00687 { 00688 if(mtf[k][1] <= 0.0) 00689 break; 00690 double x = mtf[k][0], 00691 y = VIGRA_CSTD::sqrt(-VIGRA_CSTD::log(mtf[k][1])/2.0)/M_PI; 00692 x2 += vigra::sq(x); 00693 xy += x*y; 00694 if(mtf[k][1] == minVal) 00695 break; 00696 } 00697 return xy / x2; 00698 } 00699 00700 //@} 00701 00702 } // namespace vigra 00703 00704 #endif // VIGRA_SLANTED_EDGE_MTF_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|