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" 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 177 template <
class SrcIterator,
class SrcAccessor,
class DestImage>
179 prepareSlantedEdgeInput(SrcIterator sul, SrcIterator slr, SrcAccessor src, DestImage & res,
182 unsigned int w = slr.x - sul.x;
183 unsigned int h = slr.y - sul.y;
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 unsigned int c = std::min(w, h);
193 for(
unsigned int 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");
232 edgeWidth = options.desired_edge_width,
242 if(y1 - y0 >= (
int)minimumNumberOfLines)
247 "slantedEdgeMTF(): Input edge is too slanted or image is too small");
249 y0 = std::max(y0, 0);
250 y1 = std::min(y1+1, (
int)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,
unsigned int edgeWidth)
270 using namespace functor;
277 unsigned int w = i.width(),
280 Matrix<double> m(3,3), r(3, 1), l(3, 1);
281 for(
unsigned int y = 0; y < h; ++y)
283 for(
unsigned int x = 0; x < edgeWidth; ++x)
299 for(
unsigned int y = 0; y < h; ++y)
301 for(
unsigned int x = 0; x < w; ++x)
308 template <
class Image,
class BackInsertable>
309 void slantedEdgeSubpixelShift(Image
const & img, BackInsertable & centers,
double & angle,
312 unsigned int w = img.width();
313 unsigned int h = img.height();
322 double sy = 0.0, sx = 0.0, syy = 0.0, sxy = 0.0;
323 for(
unsigned int y = 0; y < h; ++y)
327 for(
unsigned int x = 0; x < w; ++x)
338 if(c-desiredEdgeWidth < 0)
340 if(c + ew + 1 >= (
int)w)
342 for(
int 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(
unsigned int 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 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;
375 Image r(ww, 1), s(ww, 1);
376 for(
unsigned int y = 0; y < h; ++y)
378 int x0 = int(centers[y]) - w2;
380 for(; x1 < (int)ww; x1 += 4)
382 r(x1, 0) += img(x0, y);
388 for(
unsigned int x = 0; x < ww; ++x)
390 vigra_precondition(s(x,0) > 0.0,
391 "slantedEdgeMTF(): Input edge is not slanted enough");
395 result.resize(ww-1, 1);
396 for(
unsigned int x = 0; x < ww-1; ++x)
398 result(x,0) = r(x+1,0) - r(x,0);
402 template <
class Image,
class BackInsertable>
403 void slantedEdgeMTFImpl(Image
const & i, BackInsertable & mtf,
double angle,
406 unsigned int w = i.width();
407 unsigned int h = i.height();
414 if(w - 2*desiredEdgeWidth < 64)
419 magnitude.resize(w, h);
420 for(
unsigned int y = 0; y < h; ++y)
422 for(
unsigned int x = 0; x < w; ++x)
424 magnitude(x, y) =
norm(otf(x, y));
441 destImage(noise,
Point2D(w/2, 0)));
446 magnitude.resize(w, h);
447 for(
unsigned int y = 0; y < h; ++y)
449 for(
unsigned int x = 0; x < w; ++x)
461 unsigned int ww = w/4;
462 double maxVal = smooth(0,0),
464 for(
unsigned int k = 1; k < ww; ++k)
466 if(smooth(k,0) >= 0.0 && smooth(k,0) < minVal)
467 minVal = smooth(k,0);
469 double norm = maxVal-minVal;
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)
476 double xx = 4.0*k/w/slantCorrection;
477 if(n < 0.0 || xx > 1.0)
479 mtf.push_back(Result(xx, n));
621 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
623 slantedEdgeMTF(SrcIterator sul, SrcIterator slr, SrcAccessor src, BackInsertable & mtf,
627 unsigned int edgeWidth = detail::prepareSlantedEdgeInput(sul, slr, src, preparedInput, options);
628 detail::slantedEdgeShadingCorrection(preparedInput, edgeWidth);
632 detail::slantedEdgeSubpixelShift(preparedInput, centers, angle, options);
635 detail::slantedEdgeCreateOversampledLine(preparedInput, centers, oversampledLine);
637 detail::slantedEdgeMTFImpl(oversampledLine, mtf, angle, options);
640 template <
class SrcIterator,
class SrcAccessor,
class BackInsertable>
642 slantedEdgeMTF(triple<SrcIterator, SrcIterator, SrcAccessor> src, BackInsertable & mtf,
645 slantedEdgeMTF(src.first, src.second, src.third, mtf, options);
648 template <
class T1,
class S1,
class BackInsertable>
653 slantedEdgeMTF(srcImageRange(src), mtf, options);
706 template <
class Vector>
709 double minVal = mtf[0][1];
710 for(
unsigned int k = 1; k < mtf.size(); ++k)
712 if(mtf[k][1] < minVal)
717 for(
unsigned int k = 1; k < mtf.size(); ++k)
721 double x = mtf[k][0],
725 if(mtf[k][1] == minVal)
735 #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
Generic 1 dimensional convolution kernel.
Definition: separableconvolution.hxx:52
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:707
SlantedEdgeMTFOptions & mtfSmoothingScale(double scale)
Definition: slanted_edge_mtf.hxx:150
void initGaussian(double std_dev, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2251
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
const_iterator begin() const
Definition: array_vector.hxx:223
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
Two dimensional difference vector.
Definition: diff2d.hxx:185
SlantedEdgeMTFOptions & minimumEdgeWidth(unsigned int n)
Definition: slanted_edge_mtf.hxx:139
Definition: accessor.hxx:43
void separableConvolveX(...)
Performs a 1 dimensional convolution in x direction.
Two dimensional size object.
Definition: diff2d.hxx:482
Two dimensional point or position.
Definition: diff2d.hxx:592
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
void cannyEdgelList(...)
Simple implementation of Canny's edge detector.
Definition: array_vector.hxx:58
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
value_type strength
Definition: edgedetection.hxx:1359
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:1457
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
Fundamental class template for images.
Definition: basicimage.hxx:473
void initGaussianDerivative(double std_dev, int order, value_type norm, double windowRatio=0.0)
Definition: separableconvolution.hxx:2376
SlantedEdgeMTFOptions()
Definition: slanted_edge_mtf.hxx:100
void slantedEdgeMTF(...)
Determine the magnitude transfer function of the camera.
Two dimensional rectangle.
Definition: diff2d.hxx:872
Base class for, and view to, vigra::MultiArray.
Definition: multi_array.hxx:655
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
const_iterator end() const
Definition: array_vector.hxx:237
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
Definition: edgedetection.hxx:1341
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616