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

slanted_edge_mtf.hxx VIGRA

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