1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2006 by Ullrich Koethe */
4 /* Cognitive Systems Group, University of Hamburg, Germany */
5 /* */
6 /* This file is part of the VIGRA computer vision library. */
7 /* ( Version 1.5.0, Dec 07 2006 ) */
8 /* The VIGRA Website is */
9 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */
10 /* Please direct questions, bug reports, and contributions to */
11 /* koethe@informatik.uni-hamburg.de or */
12 /* vigra@kogs1.informatik.uni-hamburg.de */
13 /* */
14 /* Permission is hereby granted, free of charge, to any person */
15 /* obtaining a copy of this software and associated documentation */
16 /* files (the "Software"), to deal in the Software without */
17 /* restriction, including without limitation the rights to use, */
18 /* copy, modify, merge, publish, distribute, sublicense, and/or */
19 /* sell copies of the Software, and to permit persons to whom the */
20 /* Software is furnished to do so, subject to the following */
21 /* conditions: */
22 /* */
23 /* The above copyright notice and this permission notice shall be */
24 /* included in all copies or substantial portions of the */
25 /* Software. */
26 /* */
27 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
28 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
29 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
30 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
31 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
32 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
33 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
34 /* OTHER DEALINGS IN THE SOFTWARE. */
35 /* */
36 /************************************************************************/
37
38
39 #ifndef VIGRA_SLANTED_EDGE_MTF_HXX
40 #define VIGRA_SLANTED_EDGE_MTF_HXX
41
42 #include <algorithm>
43 #include "array_vector.hxx"
44 #include "basicgeometry.hxx"
45 #include "edgedetection.hxx"
46 #include "fftw3.hxx"
47 #include "functorexpression.hxx"
48 #include "linear_solve.hxx"
49 #include "mathutil.hxx"
50 #include "numerictraits.hxx"
51 #include "separableconvolution.hxx"
52 #include "static_assert.hxx"
53 #include "stdimage.hxx"
54 #include "transformimage.hxx"
55 #include "utilities.hxx"
56
57 namespace vigra {
58
59 /** \addtogroup SlantedEdgeMTF Camera MTF Estimation
60 Determine the magnitude transfer function (MTF) of a camera using the slanted edge method.
61 */
62 //@{
63
64 /********************************************************/
65 /* */
66 /* SlantedEdgeMTFOptions */
67 /* */
68 /********************************************************/
69
70 /** \brief Pass options to one of the \ref slantedEdgeMTF() functions.
71
72 <tt>SlantedEdgeMTFOptions</tt> is an argument objects that holds various optional
73 parameters used by the \ref slantedEdgeMTF() functions. If a parameter is not explicitly
74 set, a suitable default will be used. Changing the defaults is only necessary if you can't
75 obtain good input data, but absolutely need an MTF estimate.
76
77 <b> Usage:</b>
78
79 <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
80 Namespace: vigra
81
82 \code
83 vigra::BImage src(w,h);
84 std::vector<vigra::TinyVector<double, 2> > mtf;
85
86 ...
87 vigra::slantedEdgeMTF(srcImageRange(src), mtf,
88 vigra::SlantedEdgeMTFOptions().mtfSmoothingScale(1.0));
89
90 // print the frequency / attenuation pairs found
91 for(int k=0; k<result.size(); ++k)
92 std::cout << "frequency: " << mtf[k][0] << ", estimated attenuation: " << mtf[k][1] << std::endl;
93 \endcode
94 */
95
96 class SlantedEdgeMTFOptions
97 {
98 public:
99 /** Initialize all options with default values.
100 */
SlantedEdgeMTFOptions()101 SlantedEdgeMTFOptions()
102 : minimum_number_of_lines(20),
103 desired_edge_width(10),
104 minimum_edge_width(5),
105 mtf_smoothing_scale(2.0)
106 {}
107
108 /** Minimum number of pixels the edge must cross.
109
110 The longer the edge the more accurate the resulting MTF estimate. If you don't have good
111 data, but absolutely have to compute an MTF, you may force a lower value here.<br>
112 Default: 20
113 */
minimumNumberOfLines(unsigned int n)114 SlantedEdgeMTFOptions & minimumNumberOfLines(unsigned int n)
115 {
116 minimum_number_of_lines = n;
117 return *this;
118 }
119
120 /** Desired number of pixels perpendicular to the edge.
121
122 The larger the regions to either side of the edge,
123 the more accurate the resulting MTF estimate. If you don't have good
124 data, but absolutely have to compute an MTF, you may force a lower value here.<br>
125 Default: 10
126 */
desiredEdgeWidth(unsigned int n)127 SlantedEdgeMTFOptions & desiredEdgeWidth(unsigned int n)
128 {
129 desired_edge_width = n;
130 return *this;
131 }
132
133 /** Minimum acceptable number of pixels perpendicular to the edge.
134
135 The larger the regions to either side of the edge,
136 the more accurate the resulting MTF estimate. If you don't have good
137 data, but absolutely have to compute an MTF, you may force a lower value here.<br>
138 Default: 5
139 */
minimumEdgeWidth(unsigned int n)140 SlantedEdgeMTFOptions & minimumEdgeWidth(unsigned int n)
141 {
142 minimum_edge_width = n;
143 return *this;
144 }
145
146 /** Amount of smoothing of the computed MTF.
147
148 If the datais noisy, so will be the MTF. Thus, some smoothing is useful.<br>
149 Default: 2.0
150 */
mtfSmoothingScale(double scale)151 SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
152 {
153 vigra_precondition(scale >= 0.0,
154 "SlantedEdgeMTFOptions: MTF smoothing scale must not be < 0.");
155 mtf_smoothing_scale = scale;
156 return *this;
157 }
158
159 unsigned int minimum_number_of_lines, desired_edge_width, minimum_edge_width;
160 double mtf_smoothing_scale;
161 };
162
163 //@}
164
165 namespace detail {
166
167 struct SortEdgelsByStrength
168 {
operator ()vigra::detail::SortEdgelsByStrength169 bool operator()(Edgel const & l, Edgel const & r) const
170 {
171 return l.strength > r.strength;
172 }
173 };
174
175 /* Make sure that the edge runs vertically, intersects the top and bottom border
176 of the image, and white is on the left.
177 */
178 template <class SrcIterator, class SrcAccessor, class DestImage>
179 unsigned int
prepareSlantedEdgeInput(SrcIterator sul,SrcIterator slr,SrcAccessor src,DestImage & res,SlantedEdgeMTFOptions const & options)180 prepareSlantedEdgeInput(SrcIterator sul, SrcIterator slr, SrcAccessor src, DestImage & res,
181 SlantedEdgeMTFOptions const & options)
182 {
183 unsigned int w = slr.x - sul.x;
184 unsigned int h = slr.y - sul.y;
185
186 // find rough estimate of the edge
187 ArrayVector<Edgel> edgels;
188 cannyEdgelList(sul, slr, src, edgels, 2.0);
189 std::sort(edgels.begin(), edgels.end(), SortEdgelsByStrength());
190
191 double x = 0.0, y = 0.0, x2 = 0.0, y2 = 0.0, xy = 0.0;
192 unsigned int c = std::min(w, h);
193
194 for(unsigned int k = 0; k < c; ++k)
195 {
196 x += edgels[k].x;
197 y += edgels[k].y;
198 x2 += sq(edgels[k].x);
199 xy += edgels[k].x*edgels[k].y;
200 y2 += sq(edgels[k].y);
201 }
202 double xc = x / c, yc = y / c;
203 x2 = x2 / c - sq(xc);
204 xy = xy / c - xc*yc;
205 y2 = y2 / c - sq(yc);
206 double angle = 0.5*VIGRA_CSTD::atan2(2*xy, x2 - y2);
207
208 DestImage tmp;
209 // rotate image when slope is less than +-45 degrees
210 if(VIGRA_CSTD::fabs(angle) < M_PI / 4.0)
211 {
212 xc = yc;
213 yc = w - xc - 1;
214 std::swap(w, h);
215 tmp.resize(w, h);
216 rotateImage(srcIterRange(sul, slr, src), destImage(tmp), 90);
217 angle += M_PI / 2.0;
218 }
219 else
220 {
221 tmp.resize(w, h);
222 copyImage(srcIterRange(sul, slr, src), destImage(tmp));
223 if(angle < 0.0)
224 angle += M_PI;
225 }
226 // angle is now between pi/4 and 3*pi/4
227 double slope = VIGRA_CSTD::tan(M_PI/2.0 - angle);
228 vigra_precondition(slope != 0.0,
229 "slantedEdgeMTF(): Input edge is not slanted");
230
231 // trim image so that the edge only intersects the top and bottom border
232 unsigned int minimumNumberOfLines = options.minimum_number_of_lines, //20,
233 edgeWidth = options.desired_edge_width, // 10
234 minimumEdgeWidth = options.minimum_edge_width; // 5
235
236 int y0, y1;
237 for(; edgeWidth >= minimumEdgeWidth; --edgeWidth)
238 {
239 y0 = int(VIGRA_CSTD::floor((edgeWidth - xc) / slope + yc + 0.5));
240 y1 = int(VIGRA_CSTD::floor((w - edgeWidth - 1 - xc) / slope + yc + 0.5));
241 if(slope < 0.0)
242 std::swap(y0, y1);
243 if(y1 - y0 >= (int)minimumNumberOfLines)
244 break;
245 }
246
247 vigra_precondition(edgeWidth >= minimumEdgeWidth,
248 "slantedEdgeMTF(): Input edge is too slanted or image is too small");
249
250 y0 = std::max(y0, 0);
251 y1 = std::min(y1+1, (int)h);
252
253 res.resize(w, y1-y0);
254
255 // ensure that white is on the left
256 if(tmp(0,0) < tmp(w-1, h-1))
257 {
258 rotateImage(srcIterRange(tmp.upperLeft() + Diff2D(0, y0), tmp.upperLeft() + Diff2D(w, y1), tmp.accessor()),
259 destImage(res), 180);
260 }
261 else
262 {
263 copyImage(srcImageRange(tmp), destImage(res));
264 }
265 return edgeWidth;
266 }
267
268 template <class Image>
slantedEdgeShadingCorrection(Image & i,unsigned int edgeWidth)269 void slantedEdgeShadingCorrection(Image & i, unsigned int edgeWidth)
270 {
271 using namespace functor;
272
273 // after prepareSlantedEdgeInput(), the white region is on the left
274 // find a plane that approximates the logarithm of the white ROI
275
276 transformImage(srcImageRange(i), destImage(i), log(Arg1() + Param(1.0)));
277
278 unsigned int w = i.width(),
279 h = i.height(),
280 s = edgeWidth*h;
281
282 Matrix<double> m(3,3), r(3, 1), l(3, 1);
283 for(unsigned int y = 0; y < h; ++y)
284 {
285 for(unsigned int x = 0; x < edgeWidth; ++x)
286 {
287 l(0,0) = x;
288 l(1,0) = y;
289 l(2,0) = 1.0;
290 m += outer(l);
291 r += i(x,y)*l;
292 }
293 }
294
295 linearSolve(m, r, l);
296 double a = l(0,0),
297 b = l(1,0),
298 c = l(2,0);
299
300 // subtract the plane and go back to the non-logarithmic representation
301 for(unsigned int y = 0; y < h; ++y)
302 {
303 for(unsigned int x = 0; x < w; ++x)
304 {
305 i(x, y) = VIGRA_CSTD::exp(i(x,y) - a*x - b*y - c);
306 }
307 }
308 }
309
310 template <class Image, class BackInsertable>
slantedEdgeSubpixelShift(Image const & img,BackInsertable & centers,double & angle,SlantedEdgeMTFOptions const & options)311 void slantedEdgeSubpixelShift(Image const & img, BackInsertable & centers, double & angle,
312 SlantedEdgeMTFOptions const & options)
313 {
314 unsigned int w = img.width();
315 unsigned int h = img.height();
316
317 Image grad(w, h);
318 Kernel1D<double> kgrad;
319 kgrad.initGaussianDerivative(1.0, 1);
320
321 separableConvolveX(srcImageRange(img), destImage(grad), kernel1d(kgrad));
322
323 int desiredEdgeWidth = (int)options.desired_edge_width;
324 double sy = 0.0, sx = 0.0, syy = 0.0, sxy = 0.0;
325 for(unsigned int y = 0; y < h; ++y)
326 {
327 double a = 0.0,
328 b = 0.0;
329 for(unsigned int x = 0; x < w; ++x)
330 {
331 a += x*grad(x,y);
332 b += grad(x,y);
333 }
334 int c = int(a / b);
335 // c is biased because the numbers of black and white pixels differ
336 // repeat the analysis with a symmetric window around the edge
337 a = 0.0;
338 b = 0.0;
339 int ew = desiredEdgeWidth;
340 if(c-desiredEdgeWidth < 0)
341 ew = c;
342 if(c + ew + 1 >= (int)w)
343 ew = w - c - 1;
344 for(int x = c-ew; x < c+ew+1; ++x)
345 {
346 a += x*grad(x,y);
347 b += grad(x,y);
348 }
349 double sc = a / b;
350 sy += y;
351 sx += sc;
352 syy += sq(y);
353 sxy += sc*y;
354 }
355 // fit a line to the subpixel locations
356 double a = (h * sxy - sx * sy) / (h * syy - sq(sy));
357 double b = (sx * syy - sxy * sy) / (h * syy - sq(sy));
358
359 // compute the regularized subpixel values of the edge location
360 angle = VIGRA_CSTD::atan(a);
361 for(unsigned int y = 0; y < h; ++y)
362 {
363 centers.push_back(a * y + b);
364 }
365 }
366
367 template <class Image, class Vector>
slantedEdgeCreateOversampledLine(Image const & img,Vector const & centers,Image & result)368 void slantedEdgeCreateOversampledLine(Image const & img, Vector const & centers,
369 Image & result)
370 {
371 unsigned int w = img.width();
372 unsigned int h = img.height();
373 unsigned int w2 = std::min(std::min(int(centers[0]), int(centers[h-1])),
374 std::min(int(w - centers[0]) - 1, int(w - centers[h-1]) - 1));
375 unsigned int ww = 8*w2;
376
377 Image r(ww, 1), s(ww, 1);
378 for(unsigned int y = 0; y < h; ++y)
379 {
380 int x0 = int(centers[y]) - w2;
381 int x1 = int((VIGRA_CSTD::ceil(centers[y]) - centers[y])*4);
382 for(; x1 < (int)ww; x1 += 4)
383 {
384 r(x1, 0) += img(x0, y);
385 ++s(x1, 0);
386 ++x0;
387 }
388 }
389
390 for(unsigned int x = 0; x < ww; ++x)
391 {
392 vigra_precondition(s(x,0) > 0.0,
393 "slantedEdgeMTF(): Input edge is not slanted enough");
394 r(x,0) /= s(x,0);
395 }
396
397 result.resize(ww-1, 1);
398 for(unsigned int x = 0; x < ww-1; ++x)
399 {
400 result(x,0) = r(x+1,0) - r(x,0);
401 }
402 }
403
404 template <class Image, class BackInsertable>
slantedEdgeMTFImpl(Image const & i,BackInsertable & mtf,double angle,SlantedEdgeMTFOptions const & options)405 void slantedEdgeMTFImpl(Image const & i, BackInsertable & mtf, double angle,
406 SlantedEdgeMTFOptions const & options)
407 {
408 unsigned int w = i.width();
409 unsigned int h = i.height();
410
411 double slantCorrection = VIGRA_CSTD::cos(angle);
412 int desiredEdgeWidth = 4*options.desired_edge_width;
413
414 Image magnitude;
415
416 if(w - 2*desiredEdgeWidth < 64)
417 {
418 FFTWComplexImage otf(w, h);
419 fourierTransform(srcImageRange(i), destImage(otf));
420
421 magnitude.resize(w, h);
422 for(unsigned int y = 0; y < h; ++y)
423 {
424 for(unsigned int x = 0; x < w; ++x)
425 {
426 magnitude(x, y) = norm(otf(x, y));
427 }
428 }
429 }
430 else
431 {
432 w -= 2*desiredEdgeWidth;
433 FFTWComplexImage otf(w, h);
434 fourierTransform(srcImageRange(i, Rect2D(Point2D(desiredEdgeWidth, 0), Size2D(w, h))),
435 destImage(otf));
436
437 // create an image where the edge is skipped - presumably it only contains the
438 // noise which can then be subtracted
439 Image noise(w,h);
440 copyImage(srcImageRange(i, Rect2D(Point2D(0,0), Size2D(w/2, h))),
441 destImage(noise));
442 copyImage(srcImageRange(i, Rect2D(Point2D(i.width()-w/2, 0), Size2D(w/2, h))),
443 destImage(noise, Point2D(w/2, 0)));
444 FFTWComplexImage fnoise(w, h);
445 fourierTransform(srcImageRange(noise), destImage(fnoise));
446
447 // subtract the noise power spectrum from the total power spectrum
448 magnitude.resize(w, h);
449 for(unsigned int y = 0; y < h; ++y)
450 {
451 for(unsigned int x = 0; x < w; ++x)
452 {
453 magnitude(x, y) = VIGRA_CSTD::sqrt(std::max(0.0, squaredNorm(otf(x, y))-squaredNorm(fnoise(x, y))));
454 }
455 }
456 }
457
458 Kernel1D<double> gauss;
459 gauss.initGaussian(options.mtf_smoothing_scale);
460 Image smooth(w,h);
461 separableConvolveX(srcImageRange(magnitude), destImage(smooth), kernel1d(gauss));
462
463 unsigned int ww = w/4;
464 double maxVal = smooth(0,0),
465 minVal = maxVal;
466 for(unsigned int k = 1; k < ww; ++k)
467 {
468 if(smooth(k,0) >= 0.0 && smooth(k,0) < minVal)
469 minVal = smooth(k,0);
470 }
471 double norm = maxVal-minVal;
472
473 typedef typename BackInsertable::value_type Result;
474 mtf.push_back(Result(0.0, 1.0));
475 for(unsigned int k = 1; k < ww; ++k)
476 {
477 double n = (smooth(k,0) - minVal)/norm*sq(M_PI*k/w/VIGRA_CSTD::sin(M_PI*k/w));
478 double xx = 4.0*k/w/slantCorrection;
479 if(n < 0.0 || xx > 1.0)
480 break;
481 mtf.push_back(Result(xx, n));
482 }
483 }
484
485 } // namespace detail
486
487 /** \addtogroup SlantedEdgeMTF Camera MTF Estimation
488 Determine the magnitude transfer function (MTF) of a camera using the slanted edge method.
489 */
490 //@{
491
492 /********************************************************/
493 /* */
494 /* slantedEdgeMTF */
495 /* */
496 /********************************************************/
497
498 /** \brief Determine the magnitude transfer function of the camera.
499
500 This operator estimates the magnitude transfer function (MTF) of a camera by means of the
501 slanted edge method described in:
502
503 ISO Standard No. 12233: <i>"Photography - Electronic still picture cameras - Resolution measurements"</i>, 2000
504
505 The input must be an image that contains a single step edge with bright pixels on one side and dark pixels on
506 the other. However, the intensity values must be neither saturated nor zero. The algorithms computes the MTF
507 from the Fourier transform of the edge's derivative. Thus, if the actual MTF is unisotropic, the estimated
508 MTF does actually only apply in the direction perpendicular to the edge - several edges at different
509 orientations are required to estimate an unisotropic MTF.
510
511 The algorithm returns a sequence of frequency / attenuation pairs. The frequency axis is normalized so that the
512 Nyquist frequency of the original image is 0.5. Since the edge's derivative is computed with subpixel accuracy,
513 the attenuation can usually be computed for frequencies significantly above the Nyquist frequency as well. The
514 MTF estimate ends at either the first zero crossing of the MTF or at frequency 1, whichever comes earlier.
515
516 The present implementation improves the original slanted edge algorithm according to ISO 12233 in a number of
517 ways:
518
519 <ul>
520 <li> The edge is not required to run nearly vertically or horizontally (i.e. with a slant of approximately 5 degrees).
521 The algorithm will automatically compute the edge's actual angle and adjust estimates accordingly.
522 However, it is still necessary for the edge to be somewhat slanted, because subpixel-accurate estimation
523 of the derivative is impossible otherwise (i.e. the edge position perpendicular to the edge direction must
524 differ by at least 1 pixel between the two ends of the edge).
525
526 <li> Our implementation uses a more accurate subpixel derivative algrithm. In addition, we first perform a shading
527 correction in order to reduce possible derivative bias due to nonuniform illumination.
528
529 <li> If the input image is large enough (i.e. there are at least 20 pixels on either side of the edge over
530 the edge's entire length), our algorithm attempts to subtract the estimated noise power spectrum
531 from the estimated MTF.
532 </ul>
533
534 The source value type (<TT>SrcAccessor::value_type</TT>) must be a scalar type which is convertible to <tt>double</tt>.
535 The result is written into the \a result sequence, whose <tt>value_type</tt> must be constructible
536 from two <tt>double</tt> values. Algorithm options can be set via the \a options object
537 (see \ref vigra::NoiseNormalizationOptions for details).
538
539 <b> Declarations:</b>
540
541 pass arguments explicitly:
542 \code
543 namespace vigra {
544 template <class SrcIterator, class SrcAccessor, class BackInsertable>
545 void
546 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
547 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions());
548 }
549 \endcode
550
551 use argument objects in conjunction with \ref ArgumentObjectFactories:
552 \code
553 namespace vigra {
554 template <class SrcIterator, class SrcAccessor, class BackInsertable>
555 void
556 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
557 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
558 }
559 \endcode
560
561 <b> Usage:</b>
562
563 <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
564 Namespace: vigra
565
566 \code
567 vigra::BImage src(w,h);
568 std::vector<vigra::TinyVector<double, 2> > mtf;
569
570 ...
571 vigra::slantedEdgeMTF(srcImageRange(src), mtf);
572
573 // print the frequency / attenuation pairs found
574 for(int k=0; k<result.size(); ++k)
575 std::cout << "frequency: " << mtf[k][0] << ", estimated attenuation: " << mtf[k][1] << std::endl;
576 \endcode
577
578 <b> Required Interface:</b>
579
580 \code
581 SrcIterator upperleft, lowerright;
582 SrcAccessor src;
583
584 typedef SrcAccessor::value_type SrcType;
585 typedef NumericTraits<SrcType>::isScalar isScalar;
586 assert(isScalar::asBool == true);
587
588 double value = src(uperleft);
589
590 BackInsertable result;
591 typedef BackInsertable::value_type ResultType;
592 double intensity, variance;
593 result.push_back(ResultType(intensity, variance));
594 \endcode
595 */
596 template <class SrcIterator, class SrcAccessor, class BackInsertable>
597 void
slantedEdgeMTF(SrcIterator sul,SrcIterator slr,SrcAccessor src,BackInsertable & mtf,SlantedEdgeMTFOptions const & options=SlantedEdgeMTFOptions ())598 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
599 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
600 {
601 DImage preparedInput;
602 unsigned int edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options);
603 detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth);
604
605 ArrayVector<double> centers;
606 double angle = 0.0;
607 detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options);
608
609 DImage oversampledLine;
610 detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine);
611
612 detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options);
613 }
614
615 template <class SrcIterator, class SrcAccessor, class BackInsertable>
616 inline void
slantedEdgeMTF(triple<SrcIterator,SrcIterator,SrcAccessor> src,BackInsertable & mtf,SlantedEdgeMTFOptions const & options=SlantedEdgeMTFOptions ())617 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
618 SlantedEdgeMTFOptions const & options = SlantedEdgeMTFOptions())
619 {
620 slantedEdgeMTF(src.first, src.second, src.third, mtf, options);
621 }
622
623 /********************************************************/
624 /* */
625 /* mtfFitGaussian */
626 /* */
627 /********************************************************/
628
629 /** \brief Fit a Gaussian function to a given MTF.
630
631 This function expects a squence of frequency / attenuation pairs as produced by \ref slantedEdgeMTF()
632 and finds the best fitting Gaussian point spread function (Gaussian functions are good approximations
633 of the PSF of many real cameras). It returns the standard deviation (scale) of this function. The algorithm
634 computes the standard deviation by means of a linear least square on the logarithm of the MTF, i.e.
635 an algebraic fit rather than a Euclidean fit - thus, the resulting Gaussian may not be the one that
636 intuitively fits the data optimally.
637
638 <b> Declaration:</b>
639
640 \code
641 namespace vigra {
642 template <class Vector>
643 double mtfFitGaussian(Vector const & mtf);
644 }
645 \endcode
646
647 <b> Usage:</b>
648
649 <b>\#include</b> "<a href="slanted__edge__mtf_8hxx-source.html">vigra/slanted_edge_mtf.hxx</a>"<br>
650 Namespace: vigra
651
652 \code
653 vigra::BImage src(w,h);
654 std::vector<vigra::TinyVector<double, 2> > mtf;
655
656 ...
657 vigra::slantedEdgeMTF(srcImageRange(src), mtf);
658 double scale = vigra::mtfFitGaussian(mtf)
659
660 std::cout << "The camera PSF is approximately a Gaussian at scale " << scale << std::endl;
661 \endcode
662
663 <b> Required Interface:</b>
664
665 \code
666 Vector mtf;
667 int numberOfMeasurements = mtf.size()
668
669 double frequency = mtf[0][0];
670 double attenuation = mtf[0][1];
671 \endcode
672 */
673 template <class Vector>
mtfFitGaussian(Vector const & mtf)674 double mtfFitGaussian(Vector const & mtf)
675 {
676 double minVal = mtf[0][1];
677 for(unsigned int k = 1; k < mtf.size(); ++k)
678 {
679 if(mtf[k][1] < minVal)
680 minVal = mtf[k][1];
681 }
682 double x2 = 0.0,
683 xy = 0.0;
684 for(unsigned int k = 1; k < mtf.size(); ++k)
685 {
686 if(mtf[k][1] <= 0.0)
687 break;
688 double x = mtf[k][0],
689 y = VIGRA_CSTD::sqrt(-VIGRA_CSTD::log(mtf[k][1])/2.0)/M_PI;
690 x2 += vigra::sq(x);
691 xy += x*y;
692 if(mtf[k][1] == minVal)
693 break;
694 }
695 return xy / x2;
696 }
697
698 //@}
699
700 } // namespace vigra
701
702 #endif // VIGRA_SLANTED_EDGE_MTF_HXX
703