[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_tensorutilities.hxx
1/************************************************************************/
2/* */
3/* Copyright 2009-2010 by Ullrich Koethe */
4/* */
5/* This file is part of the VIGRA computer vision library. */
6/* The VIGRA Website is */
7/* http://hci.iwr.uni-heidelberg.de/vigra/ */
8/* Please direct questions, bug reports, and contributions to */
9/* ullrich.koethe@iwr.uni-heidelberg.de or */
10/* vigra@informatik.uni-hamburg.de */
11/* */
12/* Permission is hereby granted, free of charge, to any person */
13/* obtaining a copy of this software and associated documentation */
14/* files (the "Software"), to deal in the Software without */
15/* restriction, including without limitation the rights to use, */
16/* copy, modify, merge, publish, distribute, sublicense, and/or */
17/* sell copies of the Software, and to permit persons to whom the */
18/* Software is furnished to do so, subject to the following */
19/* conditions: */
20/* */
21/* The above copyright notice and this permission notice shall be */
22/* included in all copies or substantial portions of the */
23/* Software. */
24/* */
25/* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26/* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27/* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28/* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29/* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30/* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31/* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32/* OTHER DEALINGS IN THE SOFTWARE. */
33/* */
34/************************************************************************/
35
36#ifndef VIGRA_MULTI_TENSORUTILITIES_HXX
37#define VIGRA_MULTI_TENSORUTILITIES_HXX
38
39#include <cmath>
40#include "utilities.hxx"
41#include "mathutil.hxx"
42#include "metaprogramming.hxx"
43#include "multi_shape.hxx"
44#include "multi_pointoperators.hxx"
45
46namespace vigra {
47
48namespace detail {
49
50template <int N, class ArgumentVector, class ResultVector>
51class OuterProductFunctor
52{
53public:
54 typedef ArgumentVector argument_type;
55 typedef ResultVector result_type;
56 typedef typename ArgumentVector::value_type ValueType;
57
58 result_type operator()(argument_type const & in) const
59 {
60 result_type res;
61 for(int b=0, i=0; i<N; ++i)
62 {
63 for(int j=i; j<N; ++j, ++b)
64 {
65 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
66 }
67 }
68 return res;
69 }
70};
71
72template <int N, class ArgumentVector>
73class TensorTraceFunctor
74{
75public:
76
77 typedef ArgumentVector argument_type;
78 typedef typename ArgumentVector::value_type result_type;
79
80 result_type exec(argument_type const & v, MetaInt<1>) const
81 {
82 return v[0];
83 }
84
85 result_type exec(argument_type const & v, MetaInt<2>) const
86 {
87 return v[0] + v[2];
88 }
89
90 result_type exec(argument_type const & v, MetaInt<3>) const
91 {
92 return v[0] + v[3] + v[5];
93 }
94
95 template <int N2>
96 void exec(argument_type const &, result_type &, MetaInt<N2>) const
97 {
98 vigra_fail("tensorTraceMultiArray(): Sorry, can only handle dimensions up to 3.");
99 }
100
101 result_type operator()( const argument_type & a ) const
102 {
103 return exec(a, MetaInt<N>());
104 }
105};
106
107template <int N, class ArgumentVector, class ResultVector>
108class EigenvaluesFunctor
109{
110public:
111
112 typedef ArgumentVector argument_type;
113 typedef ResultVector result_type;
114
115 void exec(argument_type const & v, result_type & r, MetaInt<1>) const
116 {
117 symmetric2x2Eigenvalues(v[0], &r[0]);
118 }
119
120 void exec(argument_type const & v, result_type & r, MetaInt<2>) const
121 {
122 symmetric2x2Eigenvalues(v[0], v[1], v[2], &r[0], &r[1]);
123 }
124
125 void exec(argument_type const & v, result_type & r, MetaInt<3>) const
126 {
127 symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r[0], &r[1], &r[2]);
128 }
129
130 template <int N2>
131 void exec(argument_type const &, result_type &, MetaInt<N2>) const
132 {
133 vigra_fail("tensorEigenvaluesMultiArray(): Sorry, can only handle dimensions up to 3.");
134 }
135
136 result_type operator()( const argument_type & a ) const
137 {
138 result_type res;
139 exec(a, res, MetaInt<N>());
140 return res;
141 }
142};
143
144
145template <int N, class ArgumentVector>
146class DeterminantFunctor
147{
148public:
149
150 typedef ArgumentVector argument_type;
151 typedef typename ArgumentVector::value_type result_type;
152
153 result_type exec(argument_type const & v, MetaInt<1>) const
154 {
155 return v[0];
156 }
157
158 result_type exec(argument_type const & v, MetaInt<2>) const
159 {
160 return v[0]*v[2] - sq(v[1]);
161 }
162
163 result_type exec(argument_type const & v, MetaInt<3>) const
164 {
165 result_type r0, r1, r2;
166 symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r0, &r1, &r2);
167 return r0*r1*r2;
168 }
169
170 template <int N2>
171 void exec(argument_type const &, result_type &, MetaInt<N2>) const
172 {
173 vigra_fail("tensorDeterminantMultiArray(): Sorry, can only handle dimensions up to 3.");
174 }
175
176 result_type operator()( const argument_type & a ) const
177 {
178 return exec(a, MetaInt<N>());
179 }
180};
181
182} // namespace detail
183
184
185/** \addtogroup MultiPointoperators
186*/
187//@{
188
189/********************************************************/
190/* */
191/* vectorToTensorMultiArray */
192/* */
193/********************************************************/
194
195/** \brief Calculate the tensor (outer) product of a N-D vector with itself.
196
197 This function is useful to transform vector arrays into a tensor representation
198 that can be used as input to tensor based processing and analysis functions
199 (e.g. tensor smoothing). When the input array has N dimensions, the input value_type
200 must be a vector of length N, whereas the output value_type mus be vectors of length
201 N*(N-1)/2 which will represent the upper triangular part of the resulting (symmetric)
202 tensor. That is, for 2D arrays the output contains the elements
203 <tt>[t11, t12 == t21, t22]</tt> in this order, whereas it contains the elements
204 <tt>[t11, t12, t13, t22, t23, t33]</tt> for 3D arrays.
205
206 Currently, <tt>N <= 3</tt> is required.
207
208 <b> Declarations:</b>
209
210 pass arbitrary-dimensional array views:
211 \code
212 namespace vigra {
213 template <unsigned int N, class T1, class S1,
214 class T2, class S2>
215 void
216 vectorToTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
217 MultiArrayView<N, T2, S2> dest);
218 }
219 \endcode
220
221 \deprecatedAPI{vectorToTensorMultiArray}
222 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
223 \code
224 namespace vigra {
225 template <class SrcIterator, class SrcShape, class SrcAccessor,
226 class DestIterator, class DestAccessor>
227 void
228 vectorToTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
229 DestIterator di, DestAccessor dest);
230 }
231 \endcode
232 use argument objects in conjunction with \ref ArgumentObjectFactories :
233 \code
234 namespace vigra {
235 template <class SrcIterator, class SrcShape, class SrcAccessor,
236 class DestIterator, class DestAccessor>
237 void
238 vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
239 pair<DestIterator, DestAccessor> d);
240 }
241 \endcode
242 \deprecatedEnd
243
244 <b> Usage:</b>
245
246 <b>\#include</b> <vigra/multi_tensorutilities.hxx><br/>
247 Namespace: vigra
248
249 \code
250 MultiArray<3, float> vol(shape);
251 MultiArray<3, TinyVector<float, 3> > gradient(shape);
252 MultiArray<3, TinyVector<float, 6> > tensor(shape);
253
254 gaussianGradientMultiArray(vol, gradient, 2.0);
255 vectorToTensorMultiArray(gradient, tensor);
256 \endcode
257*/
258doxygen_overloaded_function(template <...> void vectorToTensorMultiArray)
259
260template <class SrcIterator, class SrcShape, class SrcAccessor,
261 class DestIterator, class DestAccessor>
262void
265{
266 static const int N = SrcShape::static_size;
267 static const int M = N*(N+1)/2;
268
269 typedef typename SrcAccessor::value_type SrcType;
270 typedef typename DestAccessor::value_type DestType;
271
272 for(int k=0; k<N; ++k)
273 if(shape[k] <=0)
274 return;
275
276 vigra_precondition(N == (int)src.size(si),
277 "vectorToTensorMultiArray(): Wrong number of channels in input array.");
278 vigra_precondition(M == (int)dest.size(di),
279 "vectorToTensorMultiArray(): Wrong number of channels in output array.");
280
281 transformMultiArray(si, shape, src, di, dest,
282 detail::OuterProductFunctor<N, SrcType, DestType>());
283}
284
285template <class SrcIterator, class SrcShape, class SrcAccessor,
286 class DestIterator, class DestAccessor>
287inline void
288vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
289 pair<DestIterator, DestAccessor> d)
290{
291 vectorToTensorMultiArray(s.first, s.second, s.third, d.first, d.second);
292}
293
294template <unsigned int N, class T1, class S1,
295 class T2, class S2>
296inline void
297vectorToTensorMultiArray(MultiArrayView<N, T1, S1> const & source,
298 MultiArrayView<N, T2, S2> dest)
299{
300 vigra_precondition(source.shape() == dest.shape(),
301 "vectorToTensorMultiArray(): shape mismatch between input and output.");
302 vectorToTensorMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
303}
304
305/********************************************************/
306/* */
307/* tensorTraceMultiArray */
308/* */
309/********************************************************/
310
311/** \brief Calculate the tensor trace for every element of a N-D tensor array.
312
313 This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
314 see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
315 symmetric tensor into a scalar array holding the tensor trace.
316
317 Currently, <tt>N <= 3</tt> is required.
318
319 <b> Declarations:</b>
320
321 pass arbitrary-dimensional array views:
322 \code
323 namespace vigra {
324 template <unsigned int N, class T1, class S1,
325 class T2, class S2>
326 void
327 tensorTraceMultiArray(MultiArrayView<N, T1, S1> const & source,
328 MultiArrayView<N, T2, S2> dest);
329 }
330 \endcode
331
332 \deprecatedAPI{tensorTraceMultiArray}
333 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
334 \code
335 namespace vigra {
336 template <class SrcIterator, class SrcShape, class SrcAccessor,
337 class DestIterator, class DestAccessor>
338 void
339 tensorTraceMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
340 DestIterator di, DestAccessor dest);
341 }
342 \endcode
343 use argument objects in conjunction with \ref ArgumentObjectFactories :
344 \code
345 namespace vigra {
346 template <class SrcIterator, class SrcShape, class SrcAccessor,
347 class DestIterator, class DestAccessor>
348 void
349 tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
350 pair<DestIterator, DestAccessor> d);
351 }
352 \endcode
353 \deprecatedEnd
354
355 <b> Usage:</b>
356
357 <b>\#include</b> <vigra/multi_tensorutilities.hxx><br/>
358 Namespace: vigra
359
360 \code
361 MultiArray<3, float> vol(shape);
362 MultiArray<3, TinyVector<float, 6> > hessian(shape);
363 MultiArray<3, float> trace(shape);
364
365 hessianOfGaussianMultiArray(vol, hessian, 2.0);
366 tensorTraceMultiArray(hessian, trace);
367 \endcode
368
369 <b> Preconditions:</b>
370
371 <tt>N == 2</tt> or <tt>N == 3</tt>
372*/
373doxygen_overloaded_function(template <...> void tensorTraceMultiArray)
374
375template <class SrcIterator, class SrcShape, class SrcAccessor,
376 class DestIterator, class DestAccessor>
377void
380{
381 static const int N = SrcShape::static_size;
382 typedef typename SrcAccessor::value_type SrcType;
383
384 transformMultiArray(si, shape, src, di, dest,
385 detail::TensorTraceFunctor<N, SrcType>());
386}
387
388template <class SrcIterator, class SrcShape, class SrcAccessor,
389 class DestIterator, class DestAccessor>
390inline void
391tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
392 pair<DestIterator, DestAccessor> d)
393{
394 tensorTraceMultiArray(s.first, s.second, s.third, d.first, d.second);
395}
396
397template <unsigned int N, class T1, class S1,
398 class T2, class S2>
399inline void
400tensorTraceMultiArray(MultiArrayView<N, T1, S1> const & source,
401 MultiArrayView<N, T2, S2> dest)
402{
403 vigra_precondition(source.shape() == dest.shape(),
404 "tensorTraceMultiArray(): shape mismatch between input and output.");
405 tensorTraceMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
406}
407
408
409/********************************************************/
410/* */
411/* tensorEigenvaluesMultiArray */
412/* */
413/********************************************************/
414
415/** \brief Calculate the tensor eigenvalues for every element of a N-D tensor array.
416
417 This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
418 see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
419 symmetric tensor into a vector-valued array holding the tensor eigenvalues (thus,
420 the destination value_type must be vectors of length N).
421
422 Currently, <tt>N <= 3</tt> is required.
423
424 <b> Declarations:</b>
425
426 pass arbitrary-dimensional array views:
427 \code
428 namespace vigra {
429 template <unsigned int N, class T1, class S1,
430 class T2, class S2>
431 void
432 tensorEigenvaluesMultiArray(MultiArrayView<N, T1, S1> const & source,
433 MultiArrayView<N, T2, S2> dest);
434 }
435 \endcode
436
437 \deprecatedAPI{tensorEigenvaluesMultiArray}
438 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
439 \code
440 namespace vigra {
441 template <class SrcIterator, class SrcShape, class SrcAccessor,
442 class DestIterator, class DestAccessor>
443 void
444 tensorEigenvaluesMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
445 DestIterator di, DestAccessor dest);
446 }
447 \endcode
448 use argument objects in conjunction with \ref ArgumentObjectFactories :
449 \code
450 namespace vigra {
451 template <class SrcIterator, class SrcShape, class SrcAccessor,
452 class DestIterator, class DestAccessor>
453 void
454 tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
455 pair<DestIterator, DestAccessor> d);
456 }
457 \endcode
458 \deprecatedEnd
459
460 <b> Usage (MultiArrayView API):</b>
461
462 <b>\#include</b> <vigra/multi_tensorutilities.hxx><br/>
463 Namespace: vigra
464
465 \code
466 MultiArray<3, float> vol(shape);
467 MultiArray<3, TinyVector<float, 6> > hessian(shape);
468 MultiArray<3, TinyVector<float, 3> > eigenvalues(shape);
469
470 hessianOfGaussianMultiArray(vol, hessian, 2.0);
471 tensorEigenvaluesMultiArray(hessian, eigenvalues);
472 \endcode
473
474 <b> Preconditions:</b>
475
476 <tt>N == 2</tt> or <tt>N == 3</tt>
477*/
478doxygen_overloaded_function(template <...> void tensorEigenvaluesMultiArray)
479
480template <class SrcIterator, class SrcShape, class SrcAccessor,
481 class DestIterator, class DestAccessor>
482void
485{
486 static const int N = SrcShape::static_size;
487 static const int M = N*(N+1)/2;
488
489 typedef typename SrcAccessor::value_type SrcType;
490 typedef typename DestAccessor::value_type DestType;
491
492 for(int k=0; k<N; ++k)
493 if(shape[k] <=0)
494 return;
495
496 vigra_precondition(M == (int)src.size(si),
497 "tensorEigenvaluesMultiArray(): Wrong number of channels in input array.");
498 vigra_precondition(N == (int)dest.size(di),
499 "tensorEigenvaluesMultiArray(): Wrong number of channels in output array.");
500
501 transformMultiArray(si, shape, src, di, dest,
502 detail::EigenvaluesFunctor<N, SrcType, DestType>());
503}
504
505template <class SrcIterator, class SrcShape, class SrcAccessor,
506 class DestIterator, class DestAccessor>
507inline void
508tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
509 pair<DestIterator, DestAccessor> d)
510{
511 tensorEigenvaluesMultiArray(s.first, s.second, s.third, d.first, d.second);
512}
513
514template <unsigned int N, class T1, class S1,
515 class T2, class S2>
516inline void
517tensorEigenvaluesMultiArray(MultiArrayView<N, T1, S1> const & source,
518 MultiArrayView<N, T2, S2> dest)
519{
520 vigra_precondition(source.shape() == dest.shape(),
521 "tensorEigenvaluesMultiArray(): shape mismatch between input and output.");
522 tensorEigenvaluesMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
523}
524
525/********************************************************/
526/* */
527/* tensorDeterminantMultiArray */
528/* */
529/********************************************************/
530
531/** \brief Calculate the tensor determinant for every element of a ND tensor array.
532
533 This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
534 see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
535 symmetric tensor into the a scalar array holding the tensor determinant.
536
537 Currently, <tt>N <= 3</tt> is required.
538
539 <b> Declarations:</b>
540
541 pass arbitrary-dimensional array views:
542 \code
543 namespace vigra {
544 template <unsigned int N, class T1, class S1,
545 class T2, class S2>
546 void
547 tensorDeterminantMultiArray(MultiArrayView<N, T1, S1> const & source,
548 MultiArrayView<N, T2, S2> dest);
549 }
550 \endcode
551
552 \deprecatedAPI{tensorDeterminantMultiArray}
553 pass \ref MultiIteratorPage "MultiIterators" and \ref DataAccessors :
554 \code
555 namespace vigra {
556 template <class SrcIterator, class SrcShape, class SrcAccessor,
557 class DestIterator, class DestAccessor>
558 void
559 tensorDeterminantMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
560 DestIterator di, DestAccessor dest);
561 }
562 \endcode
563 use argument objects in conjunction with \ref ArgumentObjectFactories :
564 \code
565 namespace vigra {
566 template <class SrcIterator, class SrcShape, class SrcAccessor,
567 class DestIterator, class DestAccessor>
568 void
569 tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
570 pair<DestIterator, DestAccessor> d);
571 }
572 \endcode
573 \deprecatedEnd
574
575 <b> Usage (MultiArrayView API):</b>
576
577 <b>\#include</b> <vigra/multi_tensorutilities.hxx><br/>
578 Namespace: vigra
579
580 \code
581 MultiArray<3, float> vol(shape);
582 MultiArray<3, TinyVector<float, 6> > hessian(shape);
583 MultiArray<3, float> determinant(shape);
584
585 hessianOfGaussianMultiArray(vol, hessian, 2.0);
586 tensorDeterminantMultiArray(hessian, determinant);
587 \endcode
588
589 <b> Preconditions:</b>
590
591 <tt>N == 2</tt> or <tt>N == 3</tt>
592*/
593doxygen_overloaded_function(template <...> void tensorDeterminantMultiArray)
594
595template <class SrcIterator, class SrcShape, class SrcAccessor,
596 class DestIterator, class DestAccessor>
597void
600{
601 typedef typename SrcAccessor::value_type SrcType;
602
603 static const int N = SrcShape::static_size;
604 static const int M = N*(N+1)/2;
605
606 for(int k=0; k<N; ++k)
607 if(shape[k] <=0)
608 return;
609
610 vigra_precondition(M == (int)src.size(si),
611 "tensorDeterminantMultiArray(): Wrong number of channels in output array.");
612
613 transformMultiArray(si, shape, src, di, dest,
614 detail::DeterminantFunctor<N, SrcType>());
615}
616
617template <class SrcIterator, class SrcShape, class SrcAccessor,
618 class DestIterator, class DestAccessor>
619inline void
620tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
621 pair<DestIterator, DestAccessor> d)
622{
623 tensorDeterminantMultiArray(s.first, s.second, s.third, d.first, d.second);
624}
625
626template <unsigned int N, class T1, class S1,
627 class T2, class S2>
628inline void
629tensorDeterminantMultiArray(MultiArrayView<N, T1, S1> const & source,
630 MultiArrayView<N, T2, S2> dest)
631{
632 vigra_precondition(source.shape() == dest.shape(),
633 "tensorDeterminantMultiArray(): shape mismatch between input and output.");
634 tensorDeterminantMultiArray(srcMultiArrayRange(source), destMultiArray(dest));
635}
636
637//@}
638
639} // namespace vigra
640
641#endif /* VIGRA_MULTI_TENSORUTILITIES_HXX */
Class for a single RGB value.
Definition rgbvalue.hxx:128
size_type size() const
Definition tinyvector.hxx:913
void tensorEigenvaluesMultiArray(...)
Calculate the tensor eigenvalues for every element of a N-D tensor array.
void tensorDeterminantMultiArray(...)
Calculate the tensor determinant for every element of a ND tensor array.
void symmetric3x3Eigenvalues(T a00, T a01, T a02, T a11, T a12, T a22, T *r0, T *r1, T *r2)
Compute the eigenvalues of a 3x3 real symmetric matrix.
Definition mathutil.hxx:754
void transformMultiArray(...)
Transform a multi-dimensional array with a unary function or functor.
void symmetric2x2Eigenvalues(T a00, T a01, T a11, T *r0, T *r1)
Compute the eigenvalues of a 2x2 real symmetric matrix.
Definition mathutil.hxx:734
void tensorTraceMultiArray(...)
Calculate the tensor trace for every element of a N-D tensor array.
void vectorToTensorMultiArray(...)
Calculate the tensor (outer) product of a N-D vector with itself.

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1