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"
53 #include "utilities.hxx"
100 : minimum_number_of_lines(20),
101 desired_edge_width(10),
102 minimum_edge_width(5),
103 mtf_smoothing_scale(2.0)
114 minimum_number_of_lines = n;
127 desired_edge_width = n;
140 minimum_edge_width = n;
151 vigra_precondition(scale >= 0.0,
152 "SlantedEdgeMTFOptions: MTF smoothing scale must not be < 0.");
153 mtf_smoothing_scale = scale;
157 unsigned int minimum_number_of_lines, desired_edge_width, minimum_edge_width;
158 double mtf_smoothing_scale;
165 struct SortEdgelsByStrength
167 bool operator()(Edgel
const & l, Edgel
const & r)
const
169 return l.strength > r.strength;
176 template <
class SrcIterator,
class SrcAccessor,
class DestImage>
178 prepareSlantedEdgeInput(SrcIterator sul, SrcIterator slr, SrcAccessor src, DestImage & res,
179 SlantedEdgeMTFOptions
const & options)
181 unsigned int w = slr.x - sul.x;
182 unsigned int h = slr.y - sul.y;
185 ArrayVector<Edgel> edgels;
187 std::sort(edgels.begin(), edgels.end(), SortEdgelsByStrength());
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);
192 for(
unsigned int k = 0; k < c; ++k)
196 x2 +=
sq(edgels[k].x);
197 xy += edgels[k].x*edgels[k].y;
198 y2 +=
sq(edgels[k].y);
200 double xc = x / c, yc = y / c;
201 x2 = x2 / c -
sq(xc);
203 y2 = y2 / c -
sq(yc);
208 if(VIGRA_CSTD::fabs(angle) < M_PI / 4.0)
214 rotateImage(srcIterRange(sul, slr, src), destImage(tmp), 90);
220 copyImage(srcIterRange(sul, slr, src), destImage(tmp));
226 vigra_precondition(slope != 0.0,
227 "slantedEdgeMTF(): Input edge is not slanted");
230 unsigned int minimumNumberOfLines = options.minimum_number_of_lines,
231 edgeWidth = options.desired_edge_width,
232 minimumEdgeWidth = options.minimum_edge_width;
235 for(; edgeWidth >= minimumEdgeWidth; --edgeWidth)
241 if(y1 - y0 >= (
int)minimumNumberOfLines)
245 vigra_precondition(edgeWidth >= minimumEdgeWidth,
246 "slantedEdgeMTF(): Input edge is too slanted or image is too small");
248 y0 = std::max(y0, 0);
249 y1 = std::min(y1+1, (
int)h);
251 res.resize(w, y1-y0);
254 if(tmp(0,0) < tmp(w-1, h-1))
256 rotateImage(srcIterRange(tmp.upperLeft() + Diff2D(0, y0), tmp.upperLeft() + Diff2D(w, y1), tmp.accessor()),
257 destImage(res), 180);
261 copyImage(srcImageRange(tmp), destImage(res));
266 template <
class Image>
267 void slantedEdgeShadingCorrection(Image & i,
unsigned int edgeWidth)
269 using namespace functor;
276 unsigned int w = i.width(),
279 Matrix<double> m(3,3), r(3, 1), l(3, 1);
280 for(
unsigned int y = 0; y < h; ++y)
282 for(
unsigned int x = 0; x < edgeWidth; ++x)
298 for(
unsigned int y = 0; y < h; ++y)
300 for(
unsigned int x = 0; x < w; ++x)
307 template <
class Image,
class BackInsertable>
308 void slantedEdgeSubpixelShift(Image
const & img, BackInsertable & centers,
double & angle,
309 SlantedEdgeMTFOptions
const & options)
311 unsigned int w = img.width();
312 unsigned int h = img.height();
315 Kernel1D<double> kgrad;
316 kgrad.initGaussianDerivative(1.0, 1);
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)
326 for(
unsigned int x = 0; x < w; ++x)
336 int ew = desiredEdgeWidth;
337 if(c-desiredEdgeWidth < 0)
339 if(c + ew + 1 >= (
int)w)
341 for(
int x = c-ew; x < c+ew+1; ++x)
353 double a = (h * sxy - sx * sy) / (h * syy -
sq(sy));
354 double b = (sx * syy - sxy * sy) / (h * syy -
sq(sy));
358 for(
unsigned int y = 0; y < h; ++y)
360 centers.push_back(a * y + b);
364 template <
class Image,
class Vector>
365 void slantedEdgeCreateOversampledLine(Image
const & img, Vector
const & centers,
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;
374 Image r(ww, 1), s(ww, 1);
375 for(
unsigned int y = 0; y < h; ++y)
377 int x0 = int(centers[y]) - w2;
379 for(; x1 < (int)ww; x1 += 4)
381 r(x1, 0) += img(x0, y);
387 for(
unsigned int x = 0; x < ww; ++x)
389 vigra_precondition(s(x,0) > 0.0,
390 "slantedEdgeMTF(): Input edge is not slanted enough");
394 result.resize(ww-1, 1);
395 for(
unsigned int x = 0; x < ww-1; ++x)
397 result(x,0) = r(x+1,0) - r(x,0);
401 template <
class Image,
class BackInsertable>
402 void slantedEdgeMTFImpl(Image
const & i, BackInsertable & mtf,
double angle,
403 SlantedEdgeMTFOptions
const & options)
405 unsigned int w = i.width();
406 unsigned int h = i.height();
409 int desiredEdgeWidth = 4*options.desired_edge_width;
413 if(w - 2*desiredEdgeWidth < 64)
418 magnitude.resize(w, h);
419 for(
unsigned int y = 0; y < h; ++y)
421 for(
unsigned int x = 0; x < w; ++x)
423 magnitude(x, y) =
norm(otf(x, y));
429 w -= 2*desiredEdgeWidth;
431 fourierTransform(srcImageRange(i, Rect2D(Point2D(desiredEdgeWidth, 0), Size2D(w, h))),
437 copyImage(srcImageRange(i, Rect2D(Point2D(0,0), Size2D(w/2, h))),
439 copyImage(srcImageRange(i, Rect2D(Point2D(i.width()-w/2, 0), Size2D(w/2, h))),
440 destImage(noise, Point2D(w/2, 0)));
445 magnitude.resize(w, h);
446 for(
unsigned int y = 0; y < h; ++y)
448 for(
unsigned int x = 0; x < w; ++x)
455 Kernel1D<double> gauss;
456 gauss.initGaussian(options.mtf_smoothing_scale);
460 unsigned int ww = w/4;
461 double maxVal = smooth(0,0),
463 for(
unsigned int k = 1; k < ww; ++k)
465 if(smooth(k,0) >= 0.0 && smooth(k,0) < minVal)
466 minVal = smooth(k,0);
468 double norm = maxVal-minVal;
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)
475 double xx = 4.0*k/w/slantCorrection;
476 if(n < 0.0 || xx > 1.0)
478 mtf.push_back(Result(xx, n));
595 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
597 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
598 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
601 unsigned int edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options);
602 detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth);
604 ArrayVector<double> centers;
606 detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options);
609 detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine);
611 detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options);
614 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
616 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
617 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
672 template <
class Vector>
675 double minVal = mtf[0][1];
676 for(
unsigned int k = 1; k < mtf.size(); ++k)
678 if(mtf[k][1] < minVal)
683 for(
unsigned int k = 1; k < mtf.size(); ++k)
687 double x = mtf[k][0],
691 if(mtf[k][1] == minVal)
701 #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:1640
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:673
SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
Definition: slanted_edge_mtf.hxx:149
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:138
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)
Definition: mathutil.hxx:341
void copyImage(...)
Copy source image into destination image.
SlantedEdgeMTFOptions & minimumNumberOfLines(unsigned int n)
Definition: slanted_edge_mtf.hxx:112
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition: matrix.hxx:1457
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
SlantedEdgeMTFOptions()
Definition: slanted_edge_mtf.hxx:99
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:125
Pass options to one of the slantedEdgeMTF() functions.
Definition: slanted_edge_mtf.hxx:94
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
bool linearSolve(const MultiArrayView< 2, T, C1 > &A, const MultiArrayView< 2, T, C2 > &b, MultiArrayView< 2, T, C3 > &res, std::string method="QR")
Definition: linear_solve.hxx:1173
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition: fftw3.hxx:1192