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"
178template <
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;
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());
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)
221 double angle = VIGRA_CSTD::atan2(-*gyr, *gxr);
222 double magnitude = VIGRA_CSTD::sqrt(*gxr * *gxr + *gyr * *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)
257 int o = VIGRA_CSTD::abs(orientationCounter(x,y));
259 if(o > maxOrientation)
262 TmpType m = VIGRA_CSTD::abs(magnitudeAccumulator(x,y));
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);
281template <
class SrcIterator,
class SrcAccessor,
282 class DestIterator,
class DestAccessor>
285 pair<DestIterator, DestAccessor> dest,
289 dest.first, dest.second,
293template <
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.");
Fundamental class template for images.
Definition basicimage.hxx:476
BasicImage & init(value_type const &pixel)
Definition basicimage.hxx:1130
traverser upperLeft()
Definition basicimage.hxx:925
void gaussianSmoothing(...)
Perform isotropic Gaussian convolution.
void radialSymmetryTransform(...)
Find centers of radial symmetry in an image.
void gaussianGradient(...)
Calculate the gradient vector by means of a 1st derivatives of Gaussian filter.