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;
166struct SortEdgelsByStrength
168 bool operator()(Edgel
const & l, Edgel
const & r)
const
170 return l.strength > r.strength;
177template <
class SrcIterator,
class SrcAccessor,
class DestImage>
179prepareSlantedEdgeInput(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);
205 double angle = 0.5*VIGRA_CSTD::atan2(2*xy, x2 - y2);
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));
226 double slope = VIGRA_CSTD::tan(M_PI/2.0 - angle);
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)
238 y0 = ptrdiff_t(VIGRA_CSTD::floor((edgeWidth - xc) / slope + yc + 0.5));
239 y1 = ptrdiff_t(VIGRA_CSTD::floor((w - edgeWidth - 1 - xc) / slope + yc + 0.5));
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));
267template <
class Image>
268void slantedEdgeShadingCorrection(Image & i, ptrdiff_t edgeWidth)
270 using namespace functor;
275 transformImage(srcImageRange(i), destImage(i), log(Arg1() + Param(1.0)));
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)
303 i(x, y) = VIGRA_CSTD::exp(i(x,y) - a*x - b*y - c);
308template <
class Image,
class BackInsertable>
309void 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));
358 angle = VIGRA_CSTD::atan(a);
359 for(ptrdiff_t y = 0; y < h; ++y)
361 centers.push_back(a * y + b);
365template <
class Image,
class Vector>
366void 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;
380 ptrdiff_t x1 = ptrdiff_t((VIGRA_CSTD::ceil(centers[y]) - centers[y])*4);
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);
403template <
class Image,
class BackInsertable>
404void slantedEdgeMTFImpl(Image
const & i, BackInsertable & mtf,
double angle,
405 SlantedEdgeMTFOptions
const & options)
407 ptrdiff_t w = i.width();
408 ptrdiff_t h = i.height();
410 double slantCorrection = VIGRA_CSTD::cos(angle);
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)
476 double n = (smooth(k,0) - minVal)/
norm*sq(M_PI*k/w/VIGRA_CSTD::sin(M_PI*k/w));
477 double xx = 4.0*k/w/slantCorrection;
478 if(n < 0.0 || xx > 1.0)
480 mtf.push_back(Result(xx, n));
622template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
641template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
643slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
644 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
649template <
class T1,
class S1,
class BackInsertable>
651slantedEdgeMTF(MultiArrayView<2, T1, S1>
const & src, BackInsertable & mtf,
652 SlantedEdgeMTFOptions
const & options = SlantedEdgeMTFOptions())
707template <
class Vector>
722 double x =
mtf[
k][0],
723 y = VIGRA_CSTD::sqrt(-VIGRA_CSTD::log(
mtf[
k][1])/2.0)/M_PI;
Class for a single RGB value.
Definition rgbvalue.hxx:128
Pass options to one of the slantedEdgeMTF() functions.
Definition slanted_edge_mtf.hxx:96
SlantedEdgeMTFOptions()
Definition slanted_edge_mtf.hxx:100
SlantedEdgeMTFOptions & desiredEdgeWidth(unsigned int n)
Definition slanted_edge_mtf.hxx:126
SlantedEdgeMTFOptions & minimumNumberOfLines(unsigned int n)
Definition slanted_edge_mtf.hxx:113
SlantedEdgeMTFOptions & minimumEdgeWidth(unsigned int n)
Definition slanted_edge_mtf.hxx:139
SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
Definition slanted_edge_mtf.hxx:150
size_type size() const
Definition tinyvector.hxx:913
void outer(const MultiArrayView< 2, T, C1 > &x, const MultiArrayView< 2, T, C2 > &y, MultiArrayView< 2, T, C3 > &r)
Definition matrix.hxx:1459
void fourierTransform(...)
Compute forward and inverse Fourier transforms.
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
double mtfFitGaussian(Vector const &mtf)
Fit a Gaussian function to a given MTF.
Definition slanted_edge_mtf.hxx:708
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037
BasicImage< FFTWComplex<> > FFTWComplexImage
Definition fftw3.hxx:1192
NumericTraits< T >::Promote sq(T t)
The square function.
Definition mathutil.hxx:382
void cannyEdgelList(...)
Simple implementation of Canny's edge detector.
void slantedEdgeMTF(...)
Determine the magnitude transfer function of the camera.
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition fftw3.hxx:1044
void copyImage(...)
Copy source image into destination image.