36 #ifndef VIGRA_SYMMETRY_HXX
37 #define VIGRA_SYMMETRY_HXX
40 #include "numerictraits.hxx"
41 #include "stdimage.hxx"
42 #include "convolution.hxx"
43 #include "multi_shape.hxx"
178 template <
class SrcIterator,
class SrcAccessor,
179 class DestIterator,
class DestAccessor>
182 DestIterator dul, DestAccessor ad,
185 vigra_precondition(scale > 0.0,
186 "radialSymmetryTransform(): Scale must be > 0");
188 int w = slr.x - sul.x;
189 int h = slr.y - sul.y;
191 if(w <= 0 || h <= 0)
return;
194 NumericTraits<typename SrcAccessor::value_type>::RealPromote TmpType;
196 typedef BasicImage<TmpType> TmpImage;
197 typedef typename TmpImage::Iterator TmpIterator;
201 IImage orientationCounter(w,h);
202 TmpImage magnitudeAccumulator(w,h);
205 destImage(gx), destImage(gy),
208 orientationCounter.init(0);
209 magnitudeAccumulator.init(NumericTraits<TmpType>::zero());
211 TmpIterator gxi = gx.upperLeft();
212 TmpIterator gyi = gy.upperLeft();
214 for(y=0; y<h; ++y, ++gxi.y, ++gyi.y)
216 typename TmpIterator::row_iterator gxr = gxi.rowIterator();
217 typename TmpIterator::row_iterator gyr = gyi.rowIterator();
219 for(
int x = 0; x<w; ++x, ++gxr, ++gyr)
224 if(magnitude < NumericTraits<TmpType>::epsilon()*10.0)
227 int dx = NumericTraits<int>::fromRealPromote(scale *
VIGRA_CSTD::cos(angle));
228 int dy = NumericTraits<int>::fromRealPromote(scale *
VIGRA_CSTD::sin(angle));
233 if(xx >= 0 && xx < w && yy >= 0 && yy < h)
235 orientationCounter(xx, yy) += 1;
236 magnitudeAccumulator(xx, yy) += detail::RequiresExplicitCast<TmpType>::cast(magnitude);
242 if(xx >= 0 && xx < w && yy >= 0 && yy < h)
244 orientationCounter(xx, yy) -= 1;
245 magnitudeAccumulator(xx, yy) -= detail::RequiresExplicitCast<TmpType>::cast(magnitude);
250 int maxOrientation = 0;
251 TmpType maxMagnitude = NumericTraits<TmpType>::zero();
255 for(
int x = 0; x<w; ++x)
259 if(o > maxOrientation)
271 for(
int x = 0; x<w; ++x)
273 double o = (double)orientationCounter(x, y) / maxOrientation;
274 magnitudeAccumulator(x, y) = detail::RequiresExplicitCast<TmpType>::cast(o * o * magnitudeAccumulator(x, y) / maxMagnitude);
278 gaussianSmoothing(srcImageRange(magnitudeAccumulator), destIter(dul, ad), 0.25*scale);
281 template <
class SrcIterator,
class SrcAccessor,
282 class DestIterator,
class DestAccessor>
285 pair<DestIterator, DestAccessor> dest,
289 dest.first, dest.second,
293 template <
class T1,
class S1,
297 MultiArrayView<2, T2, S2> dest,
300 vigra_precondition(src.shape() == dest.shape(),
301 "radialSymmetryTransform(): shape mismatch between input and output.");
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)
void radialSymmetryTransform(...)
Find centers of radial symmetry in an image.
doxygen_overloaded_function(template<...> void separableConvolveBlockwise) template< unsigned int N
Separated convolution on ChunkedArrays.
void gaussianGradient(...)
Calculate the gradient vector by means of a 1st derivatives of Gaussian filter.
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
void gaussianSmoothing(...)
Perform isotropic Gaussian convolution.
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616
BasicImage< Int32 > IImage
Definition: stdimage.hxx:116