38 #ifndef VIGRA_MULTI_CONVOLUTION_H
39 #define VIGRA_MULTI_CONVOLUTION_H
41 #include "separableconvolution.hxx"
42 #include "array_vector.hxx"
43 #include "multi_array.hxx"
44 #include "accessor.hxx"
45 #include "numerictraits.hxx"
46 #include "navigator.hxx"
47 #include "metaprogramming.hxx"
48 #include "multi_pointoperators.hxx"
49 #include "functorexpression.hxx"
64 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
65 class DestIterator,
class DestAccessor,
class KernelIterator>
67 internalSeparableConvolveMultiArrayTmp(
68 SrcIterator si, SrcShape
const & shape, SrcAccessor src,
69 DestIterator di, DestAccessor dest, KernelIterator kit)
71 enum { N = 1 + SrcIterator::level };
73 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
76 ArrayVector<TmpType> tmp( shape[0] );
78 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
79 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
82 SNavigator snav( si, shape, 0 );
83 DNavigator dnav( di, shape, 0 );
85 for( ; snav.hasMore(); snav++, dnav++ )
88 copyLine( snav.begin(), snav.end(), src,
89 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
92 typename AccessorTraits<TmpType>::default_const_accessor()),
93 destIter( dnav.begin(), dest ),
100 for(
int d = 1; d < N; ++d, ++kit )
102 DNavigator dnav( di, shape, d );
104 tmp.resize( shape[d] );
106 for( ; dnav.hasMore(); dnav++ )
109 copyLine( dnav.begin(), dnav.end(), dest,
110 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
113 typename AccessorTraits<TmpType>::default_const_accessor()),
114 destIter( dnav.begin(), dest ),
226 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
227 class DestIterator,
class DestAccessor,
class KernelIterator>
230 DestIterator d, DestAccessor dest, KernelIterator kernels )
232 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
234 if(!IsSameType<TmpType, typename DestAccessor::value_type>::boolResult)
237 MultiArray<SrcShape::static_size, TmpType> tmpArray(shape);
238 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src,
239 tmpArray.traverser_begin(),
typename AccessorTraits<TmpType>::default_accessor(), kernels );
245 detail::internalSeparableConvolveMultiArrayTmp( s, shape, src, d, dest, kernels );
249 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
250 class DestIterator,
class DestAccessor,
class KernelIterator>
253 triple<SrcIterator, SrcShape, SrcAccessor>
const & source,
254 pair<DestIterator, DestAccessor>
const & dest, KernelIterator kit )
257 dest.first, dest.second, kit );
260 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
261 class DestIterator,
class DestAccessor,
class T>
264 DestIterator d, DestAccessor dest,
265 Kernel1D<T>
const & kernel )
267 ArrayVector<Kernel1D<T> > kernels(shape.size(), kernel);
272 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
273 class DestIterator,
class DestAccessor,
class T>
276 pair<DestIterator, DestAccessor>
const & dest,
277 Kernel1D<T>
const & kernel )
279 ArrayVector<Kernel1D<T> > kernels(source.second.size(), kernel);
282 dest.first, dest.second, kernels.begin() );
347 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
348 class DestIterator,
class DestAccessor,
class T>
351 DestIterator d, DestAccessor dest,
354 enum { N = 1 + SrcIterator::level };
355 vigra_precondition( dim < N,
356 "convolveMultiArrayOneDimension(): The dimension number to convolve must be smaller "
357 "than the data dimensionality" );
359 typedef typename NumericTraits<typename DestAccessor::value_type>::RealPromote TmpType;
360 ArrayVector<TmpType> tmp( shape[dim] );
362 typedef MultiArrayNavigator<SrcIterator, N> SNavigator;
363 typedef MultiArrayNavigator<DestIterator, N> DNavigator;
365 SNavigator snav( s, shape, dim );
366 DNavigator dnav( d, shape, dim );
368 for( ; snav.hasMore(); snav++, dnav++ )
371 copyLine( snav.begin(), snav.end(), src,
372 tmp.begin(),
typename AccessorTraits<TmpType>::default_accessor() );
374 convolveLine( srcIterRange( tmp.begin(), tmp.end(),
typename AccessorTraits<TmpType>::default_const_accessor()),
375 destIter( dnav.begin(), dest ),
376 kernel1d( kernel ) );
380 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
381 class DestIterator,
class DestAccessor,
class T>
384 pair<DestIterator, DestAccessor>
const & dest,
388 dest.first, dest.second, dim, kernel );
453 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
454 class DestIterator,
class DestAccessor>
457 DestIterator d, DestAccessor dest,
double sigma )
459 Kernel1D<double> gauss;
460 gauss.initGaussian( sigma );
465 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
466 class DestIterator,
class DestAccessor>
469 pair<DestIterator, DestAccessor>
const & dest,
473 dest.first, dest.second, sigma );
548 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
549 class DestIterator,
class DestAccessor>
552 DestIterator di, DestAccessor dest,
double sigma )
554 typedef typename DestAccessor::value_type DestType;
555 typedef typename DestType::value_type DestValueType;
556 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
558 static const int N = SrcShape::static_size;
560 for(
int k=0; k<N; ++k)
564 vigra_precondition(N == dest.size(di),
565 "gaussianGradientMultiArray(): Wrong number of channels in output array.");
567 vigra_precondition(sigma > 0.0,
"gaussianGradientMultiArray(): Scale must be positive.");
569 Kernel1D<KernelType> gauss, derivative;
570 gauss.initGaussian(sigma);
572 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
575 for(
int d = 0; d < N; ++d )
577 ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
578 kernels[d].initGaussianDerivative(sigma, 1);
583 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
584 class DestIterator,
class DestAccessor>
587 pair<DestIterator, DestAccessor>
const & dest,
double sigma )
590 dest.first, dest.second, sigma );
659 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
660 class DestIterator,
class DestAccessor>
663 DestIterator di, DestAccessor dest)
665 typedef typename DestAccessor::value_type DestType;
666 typedef typename DestType::value_type DestValueType;
667 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
669 static const int N = SrcShape::static_size;
671 for(
int k=0; k<N; ++k)
675 vigra_precondition(N == dest.size(di),
676 "symmetricGradientMultiArray(): Wrong number of channels in output array.");
678 Kernel1D<KernelType> filter;
679 filter.initSymmetricGradient();
681 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
684 for(
int d = 0; d < N; ++d )
687 di, ElementAccessor(d, dest),
692 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
693 class DestIterator,
class DestAccessor>
696 pair<DestIterator, DestAccessor>
const & dest )
699 dest.first, dest.second);
769 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
770 class DestIterator,
class DestAccessor>
773 DestIterator di, DestAccessor dest,
double sigma )
775 using namespace functor;
777 typedef typename DestAccessor::value_type DestType;
778 typedef typename NumericTraits<DestType>::RealPromote KernelType;
779 typedef typename AccessorTraits<KernelType>::default_accessor DerivativeAccessor;
781 static const int N = SrcShape::static_size;
783 vigra_precondition(sigma > 0.0,
"laplacianOfGaussianMultiArray(): Scale must be positive.");
785 Kernel1D<KernelType> gauss;
786 gauss.initGaussian(sigma);
788 MultiArray<N, KernelType> derivative(shape);
791 for(
int d = 0; d < N; ++d )
793 ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
794 kernels[d].initGaussianDerivative(sigma, 2);
798 di, dest, kernels.begin());
803 derivative.traverser_begin(), DerivativeAccessor(),
806 di, dest, Arg1() + Arg2() );
811 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
812 class DestIterator,
class DestAccessor>
815 pair<DestIterator, DestAccessor>
const & dest,
double sigma )
818 dest.first, dest.second, sigma );
888 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
889 class DestIterator,
class DestAccessor>
892 DestIterator di, DestAccessor dest,
double sigma )
894 typedef typename DestAccessor::value_type DestType;
895 typedef typename DestType::value_type DestValueType;
896 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
898 static const int N = SrcShape::static_size;
899 static const int M = N*(N+1)/2;
901 for(
int k=0; k<N; ++k)
905 vigra_precondition(M == dest.size(di),
906 "hessianOfGaussianMultiArray(): Wrong number of channels in output array.");
908 vigra_precondition(sigma > 0.0,
"hessianOfGaussianMultiArray(): Scale must be positive.");
910 Kernel1D<KernelType> gauss;
911 gauss.initGaussian(sigma);
913 typedef VectorElementAccessor<DestAccessor> ElementAccessor;
916 for(
int b=0, i=0; i<N; ++i)
918 for(
int j=i; j<N; ++j, ++b)
920 ArrayVector<Kernel1D<KernelType> > kernels(N, gauss);
923 kernels[i].initGaussianDerivative(sigma, 2);
927 kernels[i].initGaussianDerivative(sigma, 1);
928 kernels[j].initGaussianDerivative(sigma, 1);
936 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
937 class DestIterator,
class DestAccessor>
940 pair<DestIterator, DestAccessor>
const & dest,
double sigma )
943 dest.first, dest.second, sigma );
948 template<
int N,
class VectorType>
949 struct StructurTensorFunctor
951 typedef VectorType result_type;
952 typedef typename VectorType::value_type ValueType;
955 VectorType operator()(T
const & in)
const
958 for(
int b=0, i=0; i<N; ++i)
960 for(
int j=i; j<N; ++j, ++b)
962 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
1041 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1042 class DestIterator,
class DestAccessor>
1045 DestIterator di, DestAccessor dest,
1046 double innerScale,
double outerScale)
1048 static const int N = SrcShape::static_size;
1049 static const int M = N*(N+1)/2;
1051 typedef typename DestAccessor::value_type DestType;
1052 typedef typename DestType::value_type DestValueType;
1053 typedef typename NumericTraits<DestValueType>::RealPromote KernelType;
1054 typedef TinyVector<KernelType, N> GradientVector;
1055 typedef typename AccessorTraits<GradientVector>::default_accessor GradientAccessor;
1057 for(
int k=0; k<N; ++k)
1061 vigra_precondition(M == dest.size(di),
1062 "structureTensorMultiArray(): Wrong number of channels in output array.");
1064 vigra_precondition(innerScale > 0.0 && outerScale >= 0.0,
1065 "structureTensorMultiArray(): Scale must be positive.");
1067 MultiArray<N, GradientVector> gradient(shape);
1069 gradient.traverser_begin(), GradientAccessor(),
1074 detail::StructurTensorFunctor<N, DestType>());
1079 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
1080 class DestIterator,
class DestAccessor>
1083 pair<DestIterator, DestAccessor>
const & dest,
1084 double innerScale,
double outerScale)
1087 dest.first, dest.second, innerScale, outerScale );
1095 #endif //-- VIGRA_MULTI_CONVOLUTION_H