[ 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 /* 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) |
html generated using doxygen and Python
|