37 #ifndef VIGRA_SLANTED_EDGE_MTF_HXX
38 #define VIGRA_SLANTED_EDGE_MTF_HXX
41 #include "array_vector.hxx"
42 #include "basicgeometry.hxx"
43 #include "edgedetection.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"
54 #include "multi_shape.hxx"
101 : minimum_number_of_lines(20),
102 desired_edge_width(10),
103 minimum_edge_width(5),
104 mtf_smoothing_scale(2.0)
115 minimum_number_of_lines = n;
128 desired_edge_width = n;
141 minimum_edge_width = n;
152 vigra_precondition(scale >= 0.0,
153 "SlantedEdgeMTFOptions: MTF smoothing scale must not be < 0.");
154 mtf_smoothing_scale = scale;
158 unsigned int minimum_number_of_lines, desired_edge_width, minimum_edge_width;
159 double mtf_smoothing_scale;
166 struct SortEdgelsByStrength
168 bool operator()(Edgel
const & l, Edgel
const & r)
const
170 return l.strength > r.strength;
177 template <
class SrcIterator,
class SrcAccessor,
class DestImage>
179 prepareSlantedEdgeInput(SrcIterator sul, SrcIterator slr, SrcAccessor src, DestImage & res,
180 SlantedEdgeMTFOptions
const & options)
182 ptrdiff_t w = slr.x - sul.x;
183 ptrdiff_t h = slr.y - sul.y;
186 ArrayVector<Edgel> edgels;
188 std::sort(edgels.begin(), edgels.end(), SortEdgelsByStrength());
190 double x = 0.0, y = 0.0, x2 = 0.0, y2 = 0.0, xy = 0.0;
191 ptrdiff_t c = std::min(w, h);
193 for(ptrdiff_t k = 0; k < c; ++k)
197 x2 +=
sq(edgels[k].x);
198 xy += edgels[k].x*edgels[k].y;
199 y2 +=
sq(edgels[k].y);
201 double xc = x / c, yc = y / c;
202 x2 = x2 / c -
sq(xc);
204 y2 = y2 / c -
sq(yc);
209 if(VIGRA_CSTD::fabs(angle) < M_PI / 4.0)
215 rotateImage(srcIterRange(sul, slr, src), destImage(tmp), 90);
221 copyImage(srcIterRange(sul, slr, src), destImage(tmp));
227 vigra_precondition(slope != 0.0,
228 "slantedEdgeMTF(): Input edge is not slanted");
231 ptrdiff_t minimumNumberOfLines = options.minimum_number_of_lines,
232 edgeWidth = options.desired_edge_width,
233 minimumEdgeWidth = options.minimum_edge_width;
235 ptrdiff_t y0 = 0, y1 = h;
236 for(; edgeWidth >= minimumEdgeWidth; --edgeWidth)
242 if(y1 - y0 >= (ptrdiff_t)minimumNumberOfLines)
246 vigra_precondition(edgeWidth >= minimumEdgeWidth,
247 "slantedEdgeMTF(): Input edge is too slanted or image is too small");
249 y0 = std::max(y0, ptrdiff_t(0));
250 y1 = std::min(y1+1, h);
252 res.resize(w, y1-y0);
255 if(tmp(0,0) < tmp(w-1, h-1))
257 rotateImage(srcIterRange(tmp.upperLeft() + Diff2D(0, y0), tmp.upperLeft() + Diff2D(w, y1), tmp.accessor()),
258 destImage(res), 180);
262 copyImage(srcImageRange(tmp), destImage(res));
267 template <
class Image>
268 void slantedEdgeShadingCorrection(Image & i, ptrdiff_t edgeWidth)
270 using namespace functor;
277 ptrdiff_t w = i.width(),
280 Matrix<double> m(3,3), r(3, 1), l(3, 1);
281 for(ptrdiff_t y = 0; y < h; ++y)
283 for(ptrdiff_t x = 0; x < edgeWidth; ++x)
299 for(ptrdiff_t y = 0; y < h; ++y)
301 for(ptrdiff_t x = 0; x < w; ++x)
308 template <
class Image,
class BackInsertable>
309 void slantedEdgeSubpixelShift(Image
const & img, BackInsertable & centers,
double & angle,
310 SlantedEdgeMTFOptions
const & options)
312 ptrdiff_t w = img.width();
313 ptrdiff_t h = img.height();
316 Kernel1D<double> kgrad;
317 kgrad.initGaussianDerivative(1.0, 1);
321 ptrdiff_t desiredEdgeWidth = (ptrdiff_t)options.desired_edge_width;
322 double sy = 0.0, sx = 0.0, syy = 0.0, sxy = 0.0;
323 for(ptrdiff_t y = 0; y < h; ++y)
327 for(ptrdiff_t x = 0; x < w; ++x)
332 ptrdiff_t c = ptrdiff_t(a / b);
337 ptrdiff_t ew = desiredEdgeWidth;
338 if(c-desiredEdgeWidth < 0)
342 for(ptrdiff_t x = c-ew; x < c+ew+1; ++x)
354 double a = (h * sxy - sx * sy) / (h * syy -
sq(sy));
355 double b = (sx * syy - sxy * sy) / (h * syy -
sq(sy));
359 for(ptrdiff_t y = 0; y < h; ++y)
361 centers.push_back(a * y + b);
365 template <
class Image,
class Vector>
366 void slantedEdgeCreateOversampledLine(Image
const & img, Vector
const & centers,
369 ptrdiff_t w = img.width();
370 ptrdiff_t h = img.height();
371 ptrdiff_t w2 = std::min(std::min(ptrdiff_t(centers[0]), ptrdiff_t(centers[h-1])),
372 std::min(ptrdiff_t(w - centers[0]) - 1,
373 ptrdiff_t(w - centers[h-1]) - 1));
376 Image r(ww, 1), s(ww, 1);
377 for(ptrdiff_t y = 0; y < h; ++y)
379 ptrdiff_t x0 = ptrdiff_t(centers[y]) - w2;
381 for(; x1 < ww; x1 += 4)
383 r(x1, 0) += img(x0, y);
389 for(ptrdiff_t x = 0; x < ww; ++x)
391 vigra_precondition(s(x,0) > 0.0,
392 "slantedEdgeMTF(): Input edge is not slanted enough");
396 result.resize(ww-1, 1);
397 for(ptrdiff_t x = 0; x < ww-1; ++x)
399 result(x,0) = r(x+1,0) - r(x,0);
403 template <
class Image,
class BackInsertable>
404 void slantedEdgeMTFImpl(Image
const & i, BackInsertable & mtf,
double angle,
405 SlantedEdgeMTFOptions
const & options)
407 ptrdiff_t w = i.width();
408 ptrdiff_t h = i.height();
411 ptrdiff_t desiredEdgeWidth = 4*options.desired_edge_width;
415 if(w - 2*desiredEdgeWidth < 64)
420 magnitude.resize(w, h);
421 for(ptrdiff_t y = 0; y < h; ++y)
423 for(ptrdiff_t x = 0; x < w; ++x)
425 magnitude(x, y) =
norm(otf(x, y));
431 w -= 2*desiredEdgeWidth;
433 fourierTransform(srcImageRange(i, Rect2D(Point2D(desiredEdgeWidth, 0), Size2D(w, h))),
439 copyImage(srcImageRange(i, Rect2D(Point2D(0,0), Size2D(w/2, h))),
441 copyImage(srcImageRange(i, Rect2D(Point2D(i.width()-w/2, 0), Size2D(w/2, h))),
442 destImage(noise, Point2D(w/2, 0)));
447 magnitude.resize(w, h);
448 for(ptrdiff_t y = 0; y < h; ++y)
450 for(ptrdiff_t x = 0; x < w; ++x)
457 Kernel1D<double> gauss;
458 gauss.initGaussian(options.mtf_smoothing_scale);
463 double maxVal = smooth(0,0),
465 for(ptrdiff_t k = 1; k < ww; ++k)
467 if(smooth(k,0) >= 0.0 && smooth(k,0) < minVal)
468 minVal = smooth(k,0);
470 double norm = maxVal-minVal;
472 typedef typename BackInsertable::value_type Result;
473 mtf.push_back(Result(0.0, 1.0));
474 for(ptrdiff_t k = 1; k < ww; ++k)
477 double xx = 4.0*k/w/slantCorrection;
478 if(n < 0.0 || xx > 1.0)
480 mtf.push_back(Result(xx, n));
622 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
624 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
625 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
628 ptrdiff_t edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options);
629 detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth);
631 ArrayVector<double> centers;
633 detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options);
636 detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine);
638 detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options);
641 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
643 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
644 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
649 template <
class T1,
class S1,
class BackInsertable>
651 slantedEdgeMTF(MultiArrayView<2, T1, S1>
const & src, BackInsertable & mtf,
652 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
707 template <
class Vector>
710 double minVal = mtf[0][1];
711 for(ptrdiff_t k = 1; k < (ptrdiff_t)mtf.size(); ++k)
713 if(mtf[k][1] < minVal)
718 for(ptrdiff_t k = 1; k < (ptrdiff_t)mtf.size(); ++k)
722 double x = mtf[k][0],
726 if(mtf[k][1] == minVal)
736 #endif // VIGRA_SLANTED_EDGE_MTF_HXX
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
double mtfFitGaussian(Vector const &mtf)
Fit a Gaussian function to a given MTF.
Definition: slanted_edge_mtf.hxx:708
SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
Definition: slanted_edge_mtf.hxx:150
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
SlantedEdgeMTFOptions & minimumEdgeWidth(unsigned int n)
Definition: slanted_edge_mtf.hxx:139
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
void cannyEdgelList(...)
Simple implementation of Canny's edge detector.
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:382
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void copyImage(...)
Copy source image into destination image.
SlantedEdgeMTFOptions & minimumNumberOfLines(unsigned int n)
Definition: slanted_edge_mtf.hxx:113
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1459
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
SlantedEdgeMTFOptions()
Definition: slanted_edge_mtf.hxx:100
void slantedEdgeMTF(...)
Determine the magnitude transfer function of the camera.
linalg::TemporaryMatrix< T > atan(MultiArrayView< 2, T, C > const &v)
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
SlantedEdgeMTFOptions & desiredEdgeWidth(unsigned int n)
Definition: slanted_edge_mtf.hxx:126
Pass options to one of the slantedEdgeMTF() functions.
Definition: slanted_edge_mtf.hxx:95
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
BasicImage< double > DImage
Definition: stdimage.hxx:153
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition: fftw3.hxx:1192