36 #ifndef VIGRA_ORIENTEDTENSORFILTERS_HXX
37 #define VIGRA_ORIENTEDTENSORFILTERS_HXX
40 #include "utilities.hxx"
41 #include "initimage.hxx"
42 #include "stdconvolution.hxx"
128 template <
class SrcIterator,
class SrcAccessor,
129 class DestIterator,
class DestAccessor>
130 void hourGlassFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
131 DestIterator dul, DestAccessor dest,
132 double sigma,
double rho)
134 vigra_precondition(sigma >= 0.0 && rho >= 0.0,
135 "hourGlassFilter(): sigma and rho must be >= 0.0");
136 vigra_precondition(src.size(sul) == 3,
137 "hourGlassFilter(): input image must have 3 bands.");
138 vigra_precondition(dest.size(dul) == 3,
139 "hourGlassFilter(): output image must have 3 bands.");
143 int w = slr.x - sul.x;
144 int h = slr.y - sul.y;
147 double sigma2 = -0.5 / sigma / sigma;
148 double rho2 = -0.5 / rho / rho;
149 double norm = 1.0 / (2.0 * M_PI * sigma * sigma);
151 initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
153 for(
int y=0; y<h; ++y, ++sul.y, ++dul.y)
156 DestIterator d = dul;
157 for(
int x=0; x<w; ++x, ++s.x, ++d.x)
160 2.0*src.getComponent(s,1),
161 (double)src.getComponent(s,0) - src.getComponent(s,2));
165 double x0 = x - radius < 0 ? -x : -radius;
166 double y0 = y - radius < 0 ? -y : -radius;
167 double x1 = x + radius >= w ? w - x - 1 : radius;
168 double y1 = y + radius >= h ? h - y - 1 : radius;
170 DestIterator dwul = d + Diff2D((
int)x0, (
int)y0);
172 for(
double yy=y0; yy <= y1; ++yy, ++dwul.y)
174 typename DestIterator::row_iterator dw = dwul.rowIterator();
175 for(
double xx=x0; xx <= x1; ++xx, ++dw)
177 double r2 = xx*xx + yy*yy;
178 double p = u*xx - v*yy;
179 double q = v*xx + u*yy;
180 double kernel = (p == 0.0) ?
185 dest.set(dest(dw) + kernel*src(s), dw);
192 template <
class SrcIterator,
class SrcAccessor,
193 class DestIterator,
class DestAccessor>
196 pair<DestIterator, DestAccessor> d,
197 double sigma,
double rho)
199 hourGlassFilter(s.first, s.second, s.third, d.first, d.second, sigma, rho);
208 template <
class SrcIterator,
class SrcAccessor,
209 class DestIterator,
class DestAccessor>
210 void ellipticGaussian(SrcIterator sul, SrcIterator slr, SrcAccessor src,
211 DestIterator dul, DestAccessor dest,
212 double sigmax,
double sigmin)
214 vigra_precondition(sigmax >= sigmin && sigmin >= 0.0,
215 "ellipticGaussian(): "
216 "sigmax >= sigmin and sigmin >= 0.0 required");
218 int w = slr.x - sul.x;
219 int h = slr.y - sul.y;
222 double sigmin2 = -0.5 / sigmin / sigmin;
223 double norm = 1.0 / (2.0 * M_PI * sigmin * sigmax);
225 initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<typename DestAccessor::value_type>::zero());
227 for(
int y=0; y<h; ++y, ++sul.y, ++dul.y)
230 DestIterator d = dul;
231 for(
int x=0; x<w; ++x, ++s.x, ++d.x)
234 NumericTraits<typename SrcAccessor::component_type>::RealPromote TmpType;
235 TmpType d1 = src.getComponent(s,0) + src.getComponent(s,2);
236 TmpType d2 = src.getComponent(s,0) - src.getComponent(s,2);
237 TmpType d3 = 2.0 * src.getComponent(s,1);
239 TmpType excentricity = 1.0 - (d1 - d4) / (d1 + d4);
240 double sigmax2 = -0.5 /
sq((sigmax - sigmin)*excentricity + sigmin);
246 double x0 = x - radius < 0 ? -x : -radius;
247 double y0 = y - radius < 0 ? -y : -radius;
248 double x1 = x + radius >= w ? w - x - 1 : radius;
249 double y1 = y + radius >= h ? h - y - 1 : radius;
251 DestIterator dwul = d + Diff2D((
int)x0, (
int)y0);
253 for(
double yy=y0; yy <= y1; ++yy, ++dwul.y)
255 typename DestIterator::row_iterator dw = dwul.rowIterator();
256 for(
double xx=x0; xx <= x1; ++xx, ++dw)
258 double p = u*xx - v*yy;
259 double q = v*xx + u*yy;
261 dest.set(dest(dw) + kernel*src(s), dw);
268 template <
class SrcIterator,
class SrcAccessor,
269 class DestIterator,
class DestAccessor>
271 void ellipticGaussian(triple<SrcIterator, SrcIterator, SrcAccessor> s,
272 pair<DestIterator, DestAccessor> d,
273 double sigmax,
double sigmin)
275 ellipticGaussian(s.first, s.second, s.third, d.first, d.second, sigmax, sigmin);
284 class FoerstnerKernelBase
287 typedef double ResultType;
288 typedef double WeightType;
290 typedef VectorType::value_type ValueType;
292 FoerstnerKernelBase(
double scale,
bool ringShaped =
false)
293 : radius_((int)(3.0*scale+0.5)),
294 weights_(2*radius_+1, 2*radius_+1),
295 vectors_(2*radius_+1, 2*radius_+1)
297 double norm = 1.0 / (2.0 * M_PI * scale * scale);
298 double s2 = -0.5 / scale / scale;
300 for(
int y = -radius_; y <= radius_; ++y)
302 for(
int x = -radius_; x <= radius_; ++x)
304 double d2 = x*x + y*y;
306 vectors_(x+radius_,y+radius_) = d != 0.0 ?
307 VectorType(x/d, -y/d) :
308 VectorType(ValueType(0), ValueType(0));
309 weights_(x+radius_,y+radius_) = ringShaped ?
311 norm * VIGRA_CSTD::
exp(d2 * s2);
316 ResultType operator()(
int x,
int y, VectorType
const &)
const
319 return weights_(radius_, radius_);
327 class FoerstnerRingKernelBase
328 :
public FoerstnerKernelBase
331 FoerstnerRingKernelBase(
double scale)
332 : FoerstnerKernelBase(scale, true)
337 :
public FoerstnerKernelBase
340 typedef double ResultType;
341 typedef double WeightType;
344 Cos2RingKernel(
double scale)
345 : FoerstnerKernelBase(scale, true)
348 ResultType operator()(
int x,
int y, VectorType
const & v)
const
351 return weights_(radius_, radius_);
352 double d =
dot(vectors_(x+radius_, y+radius_), v);
353 return d * d * weights_(x+radius_, y+radius_);
358 :
public FoerstnerKernelBase
361 typedef double ResultType;
362 typedef double WeightType;
365 Cos2Kernel(
double scale)
366 : FoerstnerKernelBase(scale, false)
369 ResultType operator()(
int x,
int y, VectorType
const & v)
const
372 return weights_(radius_, radius_);
373 double d =
dot(vectors_(x+radius_, y+radius_), v);
374 return d * d * weights_(x+radius_, y+radius_);
379 :
public FoerstnerKernelBase
382 typedef double ResultType;
383 typedef double WeightType;
386 Sin2RingKernel(
double scale)
387 : FoerstnerKernelBase(scale, true)
390 ResultType operator()(
int x,
int y, VectorType
const & v)
const
393 return weights_(radius_, radius_);
394 double d =
dot(vectors_(x+radius_, y+radius_), v);
395 return (1.0 - d * d) * weights_(x+radius_, y+radius_);
400 :
public FoerstnerKernelBase
403 typedef double ResultType;
404 typedef double WeightType;
407 Sin2Kernel(
double scale)
408 : FoerstnerKernelBase(scale, false)
411 ResultType operator()(
int x,
int y, VectorType
const & v)
const
414 return weights_(radius_, radius_);
415 double d =
dot(vectors_(x+radius_, y+radius_), v);
416 return (1.0 - d * d) * weights_(x+radius_, y+radius_);
421 :
public FoerstnerKernelBase
424 typedef double ResultType;
425 typedef double WeightType;
428 Sin6RingKernel(
double scale)
429 : FoerstnerKernelBase(scale, true)
432 ResultType operator()(
int x,
int y, VectorType
const & v)
const
435 return weights_(radius_, radius_);
436 double d =
dot(vectors_(x+radius_, y+radius_), v);
437 return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
442 :
public FoerstnerKernelBase
445 typedef double ResultType;
446 typedef double WeightType;
449 Sin6Kernel(
double scale)
450 : FoerstnerKernelBase(scale, false)
453 ResultType operator()(
int x,
int y, VectorType
const & v)
const
456 return weights_(radius_, radius_);
457 double d =
dot(vectors_(x+radius_, y+radius_), v);
458 return VIGRA_CSTD::pow(1.0 - d * d, 3) * weights_(x+radius_, y+radius_);
463 :
public FoerstnerKernelBase
466 typedef double ResultType;
467 typedef double WeightType;
470 Cos6RingKernel(
double scale)
471 : FoerstnerKernelBase(scale, true)
474 ResultType operator()(
int x,
int y, VectorType
const & v)
const
477 return weights_(radius_, radius_);
478 double d =
dot(vectors_(x+radius_, y+radius_), v);
479 return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
484 :
public FoerstnerKernelBase
487 typedef double ResultType;
488 typedef double WeightType;
491 Cos6Kernel(
double scale)
492 : FoerstnerKernelBase(scale, false)
495 ResultType operator()(
int x,
int y, VectorType
const & v)
const
498 return weights_(radius_, radius_);
499 double d =
dot(vectors_(x+radius_, y+radius_), v);
500 return (1.0 - VIGRA_CSTD::pow(1.0 - d * d, 3)) * weights_(x+radius_, y+radius_);
510 template <
class SrcIterator,
class SrcAccessor,
511 class DestIterator,
class DestAccessor,
513 void orientedTrigonometricFilter(SrcIterator sul, SrcIterator slr, SrcAccessor src,
514 DestIterator dul, DestAccessor dest,
515 Kernel
const & kernel)
517 vigra_precondition(src.size(sul) == 2,
518 "orientedTrigonometricFilter(): input image must have 2 bands.");
519 vigra_precondition(dest.size(dul) == 3,
520 "orientedTrigonometricFilter(): output image must have 3 bands.");
522 int w = slr.x - sul.x;
523 int h = slr.y - sul.y;
524 int radius = kernel.radius_;
526 typedef typename SrcAccessor::value_type VectorType;
527 typedef typename DestAccessor::value_type TensorType;
529 initImage(dul, dul+Diff2D(w,h), dest, NumericTraits<TensorType>::zero());
531 for(
int y=0; y<h; ++y, ++sul.y, ++dul.y)
534 DestIterator d = dul;
535 for(
int x=0; x<w; ++x, ++s.x, ++d.x)
537 int x0 = x - radius < 0 ? -x : -radius;
538 int y0 = y - radius < 0 ? -y : -radius;
539 int x1 = x + radius >= w ? w - x - 1 : radius;
540 int y1 = y + radius >= h ? h - y - 1 : radius;
542 VectorType v(src(s));
543 TensorType t(
sq(v[0]), v[0]*v[1],
sq(v[1]));
544 double sqMag = t[0] + t[2];
551 for(dd.y = y0; dd.y <= y1; ++dd.y)
553 for(dd.x = x0; dd.x <= x1; ++dd.x)
555 dest.set(dest(d, dd) + kernel(dd.x, dd.y, v) * t, d, dd);
562 template <
class SrcIterator,
class SrcAccessor,
563 class DestIterator,
class DestAccessor,
566 void orientedTrigonometricFilter(triple<SrcIterator, SrcIterator, SrcAccessor> s,
567 pair<DestIterator, DestAccessor> d,
568 Kernel
const & kernel)
570 orientedTrigonometricFilter(s.first, s.second, s.third, d.first, d.second, kernel);