[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
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> <<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>><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 datais 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, y1; 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 unisotropic, 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 unisotropic 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 algrithm. 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> <<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>><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 squence 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> <<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>><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) |
html generated using doxygen and Python
|