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

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

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.7.1 (Thu Jun 14 2012)