36 #ifndef VIGRA_MULTI_TENSORUTILITIES_HXX
37 #define VIGRA_MULTI_TENSORUTILITIES_HXX
40 #include "utilities.hxx"
41 #include "mathutil.hxx"
42 #include "metaprogramming.hxx"
43 #include "multi_pointoperators.hxx"
49 template <
int N,
class ArgumentVector,
class ResultVector>
50 class OuterProductFunctor
53 typedef ArgumentVector argument_type;
54 typedef ResultVector result_type;
55 typedef typename ArgumentVector::value_type ValueType;
57 result_type operator()(argument_type
const & in)
const
60 for(
int b=0, i=0; i<N; ++i)
62 for(
int j=i; j<N; ++j, ++b)
64 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
71 template <
int N,
class ArgumentVector>
72 class TensorTraceFunctor
76 typedef ArgumentVector argument_type;
77 typedef typename ArgumentVector::value_type result_type;
79 result_type exec(argument_type
const & v, MetaInt<1>)
const
84 result_type exec(argument_type
const & v, MetaInt<2>)
const
89 result_type exec(argument_type
const & v, MetaInt<3>)
const
91 return v[0] + v[3] + v[5];
95 void exec(argument_type
const & v, result_type & r, MetaInt<N2>)
const
97 vigra_fail(
"tensorTraceMultiArray(): Sorry, can only handle dimensions up to 3.");
100 result_type operator()(
const argument_type & a )
const
102 return exec(a, MetaInt<N>());
106 template <
int N,
class ArgumentVector,
class ResultVector>
107 class EigenvaluesFunctor
111 typedef ArgumentVector argument_type;
112 typedef ResultVector result_type;
114 void exec(argument_type
const & v, result_type & r, MetaInt<1>)
const
119 void exec(argument_type
const & v, result_type & r, MetaInt<2>)
const
124 void exec(argument_type
const & v, result_type & r, MetaInt<3>)
const
130 void exec(argument_type
const & v, result_type & r, MetaInt<N2>)
const
132 vigra_fail(
"tensorEigenvaluesMultiArray(): Sorry, can only handle dimensions up to 3.");
135 result_type operator()(
const argument_type & a )
const
138 exec(a, res, MetaInt<N>());
144 template <
int N,
class ArgumentVector>
145 class DeterminantFunctor
149 typedef ArgumentVector argument_type;
150 typedef typename ArgumentVector::value_type result_type;
152 result_type exec(argument_type
const & v, MetaInt<1>)
const
157 result_type exec(argument_type
const & v, MetaInt<2>)
const
159 return v[0]*v[2] -
sq(v[1]);
162 result_type exec(argument_type
const & v, MetaInt<3>)
const
164 result_type r0, r1, r2;
170 void exec(argument_type
const & v, result_type & r, MetaInt<N2>)
const
172 vigra_fail(
"tensorDeterminantMultiArray(): Sorry, can only handle dimensions up to 3.");
175 result_type operator()(
const argument_type & a )
const
177 return exec(a, MetaInt<N>());
247 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
248 class DestIterator,
class DestAccessor>
251 DestIterator di, DestAccessor dest)
253 static const int N = SrcShape::static_size;
254 static const int M = N*(N+1)/2;
256 typedef typename SrcAccessor::value_type SrcType;
257 typedef typename DestAccessor::value_type DestType;
259 for(
int k=0; k<N; ++k)
263 vigra_precondition(N == (
int)src.size(si),
264 "vectorToTensorMultiArray(): Wrong number of channels in input array.");
265 vigra_precondition(M == (
int)dest.size(di),
266 "vectorToTensorMultiArray(): Wrong number of channels in output array.");
269 detail::OuterProductFunctor<N, SrcType, DestType>());
272 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
273 class DestIterator,
class DestAccessor>
276 pair<DestIterator, DestAccessor> d)
336 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
337 class DestIterator,
class DestAccessor>
340 DestIterator di, DestAccessor dest)
342 static const int N = SrcShape::static_size;
343 typedef typename SrcAccessor::value_type SrcType;
346 detail::TensorTraceFunctor<N, SrcType>());
349 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
350 class DestIterator,
class DestAccessor>
353 pair<DestIterator, DestAccessor> d)
415 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
416 class DestIterator,
class DestAccessor>
419 DestIterator di, DestAccessor dest)
421 static const int N = SrcShape::static_size;
422 static const int M = N*(N+1)/2;
424 typedef typename SrcAccessor::value_type SrcType;
425 typedef typename DestAccessor::value_type DestType;
427 for(
int k=0; k<N; ++k)
431 vigra_precondition(M == (
int)src.size(si),
432 "tensorEigenvaluesMultiArray(): Wrong number of channels in input array.");
433 vigra_precondition(N == (
int)dest.size(di),
434 "tensorEigenvaluesMultiArray(): Wrong number of channels in output array.");
437 detail::EigenvaluesFunctor<N, SrcType, DestType>());
440 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
441 class DestIterator,
class DestAccessor>
444 pair<DestIterator, DestAccessor> d)
504 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
505 class DestIterator,
class DestAccessor>
508 DestIterator di, DestAccessor dest)
510 typedef typename SrcAccessor::value_type SrcType;
512 static const int N = SrcShape::static_size;
513 static const int M = N*(N+1)/2;
515 for(
int k=0; k<N; ++k)
519 vigra_precondition(M == (
int)src.size(si),
520 "tensorDeterminantMultiArray(): Wrong number of channels in output array.");
523 detail::DeterminantFunctor<N, SrcType>());
526 template <
class SrcIterator,
class SrcShape,
class SrcAccessor,
527 class DestIterator,
class DestAccessor>
530 pair<DestIterator, DestAccessor> d)