[ 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_pointoperators.hxx"
44 
45 namespace vigra {
46 
47 namespace detail {
48 
49 template <int N, class ArgumentVector, class ResultVector>
50 class OuterProductFunctor
51 {
52 public:
53  typedef ArgumentVector argument_type;
54  typedef ResultVector result_type;
55  typedef typename ArgumentVector::value_type ValueType;
56 
57  result_type operator()(argument_type const & in) const
58  {
59  result_type res;
60  for(int b=0, i=0; i<N; ++i)
61  {
62  for(int j=i; j<N; ++j, ++b)
63  {
64  res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
65  }
66  }
67  return res;
68  }
69 };
70 
71 template <int N, class ArgumentVector>
72 class TensorTraceFunctor
73 {
74 public:
75 
76  typedef ArgumentVector argument_type;
77  typedef typename ArgumentVector::value_type result_type;
78 
79  result_type exec(argument_type const & v, MetaInt<1>) const
80  {
81  return v[0];
82  }
83 
84  result_type exec(argument_type const & v, MetaInt<2>) const
85  {
86  return v[0] + v[2];
87  }
88 
89  result_type exec(argument_type const & v, MetaInt<3>) const
90  {
91  return v[0] + v[3] + v[5];
92  }
93 
94  template <int N2>
95  void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
96  {
97  vigra_fail("tensorTraceMultiArray(): Sorry, can only handle dimensions up to 3.");
98  }
99 
100  result_type operator()( const argument_type & a ) const
101  {
102  return exec(a, MetaInt<N>());
103  }
104 };
105 
106 template <int N, class ArgumentVector, class ResultVector>
107 class EigenvaluesFunctor
108 {
109 public:
110 
111  typedef ArgumentVector argument_type;
112  typedef ResultVector result_type;
113 
114  void exec(argument_type const & v, result_type & r, MetaInt<1>) const
115  {
116  symmetric2x2Eigenvalues(v[0], &r[0]);
117  }
118 
119  void exec(argument_type const & v, result_type & r, MetaInt<2>) const
120  {
121  symmetric2x2Eigenvalues(v[0], v[1], v[2], &r[0], &r[1]);
122  }
123 
124  void exec(argument_type const & v, result_type & r, MetaInt<3>) const
125  {
126  symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r[0], &r[1], &r[2]);
127  }
128 
129  template <int N2>
130  void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
131  {
132  vigra_fail("tensorEigenvaluesMultiArray(): Sorry, can only handle dimensions up to 3.");
133  }
134 
135  result_type operator()( const argument_type & a ) const
136  {
137  result_type res;
138  exec(a, res, MetaInt<N>());
139  return res;
140  }
141 };
142 
143 
144 template <int N, class ArgumentVector>
145 class DeterminantFunctor
146 {
147 public:
148 
149  typedef ArgumentVector argument_type;
150  typedef typename ArgumentVector::value_type result_type;
151 
152  result_type exec(argument_type const & v, MetaInt<1>) const
153  {
154  return v[0];
155  }
156 
157  result_type exec(argument_type const & v, MetaInt<2>) const
158  {
159  return v[0]*v[2] - sq(v[1]);
160  }
161 
162  result_type exec(argument_type const & v, MetaInt<3>) const
163  {
164  result_type r0, r1, r2;
165  symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r0, &r1, &r2);
166  return r0*r1*r2;
167  }
168 
169  template <int N2>
170  void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
171  {
172  vigra_fail("tensorDeterminantMultiArray(): Sorry, can only handle dimensions up to 3.");
173  }
174 
175  result_type operator()( const argument_type & a ) const
176  {
177  return exec(a, MetaInt<N>());
178  }
179 };
180 
181 } // namespace detail
182 
183 
184 /** \addtogroup MultiPointoperators
185 */
186 //@{
187 
188 /********************************************************/
189 /* */
190 /* vectorToTensorMultiArray */
191 /* */
192 /********************************************************/
193 
194 /** \brief Calculate the tensor (outer) product of a N-D vector with itself.
195 
196  This function is useful to transform vector arrays into a tensor representation
197  that can be used as input to tensor based processing and analysis functions
198  (e.g. tensor smoothing). When the input array has N dimensions, the input value_type
199  must be a vector of length N, whereas the output value_type mus be vectors of length
200  N*(N-1)/2 which will represent the upper triangular part of the resulting (symmetric)
201  tensor. That is, for 2D arrays the output contains the elements
202  <tt>[t11, t12 == t21, t22]</tt> in this order, whereas it contains the elements
203  <tt>[t11, t12, t13, t22, t23, t33]</tt> for 3D arrays.
204 
205  Currently, <tt>N <= 3</tt> is required.
206 
207  <b> Declarations:</b>
208 
209  pass arguments explicitly:
210  \code
211  namespace vigra {
212  template <class SrcIterator, class SrcShape, class SrcAccessor,
213  class DestIterator, class DestAccessor>
214  void
215  vectorToTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
216  DestIterator di, DestAccessor dest);
217  }
218  \endcode
219 
220  use argument objects in conjunction with \ref ArgumentObjectFactories :
221  \code
222  namespace vigra {
223  template <class SrcIterator, class SrcShape, class SrcAccessor,
224  class DestIterator, class DestAccessor>
225  void
226  vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
227  pair<DestIterator, DestAccessor> d);
228  }
229  \endcode
230 
231  <b> Usage:</b>
232 
233  <b>\#include</b> <vigra/multi_tensorutilities.hxx>
234 
235  \code
236  MultiArray<3, float> vol(shape);
237  MultiArray<3, TinyVector<float, 3> > gradient(shape);
238  MultiArray<3, TinyVector<float, 6> > tensor(shape);
239 
240  gaussianGradientMultiArray(srcMultiArrayRange(vol), destMultiArray(gradient), 2.0);
241  vectorToTensorMultiArray(srcMultiArrayRange(gradient), destMultiArray(tensor));
242  \endcode
243 
244 */
245 doxygen_overloaded_function(template <...> void vectorToTensorMultiArray)
246 
247 template <class SrcIterator, class SrcShape, class SrcAccessor,
248  class DestIterator, class DestAccessor>
249 void
250 vectorToTensorMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
251  DestIterator di, DestAccessor dest)
252 {
253  static const int N = SrcShape::static_size;
254  static const int M = N*(N+1)/2;
255 
256  typedef typename SrcAccessor::value_type SrcType;
257  typedef typename DestAccessor::value_type DestType;
258 
259  for(int k=0; k<N; ++k)
260  if(shape[k] <=0)
261  return;
262 
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.");
267 
268  transformMultiArray(si, shape, src, di, dest,
269  detail::OuterProductFunctor<N, SrcType, DestType>());
270 }
271 
272 template <class SrcIterator, class SrcShape, class SrcAccessor,
273  class DestIterator, class DestAccessor>
274 inline void
275 vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
276  pair<DestIterator, DestAccessor> d)
277 {
278  vectorToTensorMultiArray(s.first, s.second, s.third, d.first, d.second);
279 }
280 
281 /********************************************************/
282 /* */
283 /* tensorTraceMultiArray */
284 /* */
285 /********************************************************/
286 
287 /** \brief Calculate the tensor trace for every element of a N-D tensor array.
288 
289  This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
290  see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
291  symmetric tensor into a scalar array holding the tensor trace.
292 
293  Currently, <tt>N <= 3</tt> is required.
294 
295  <b> Declarations:</b>
296 
297  pass arguments explicitly:
298  \code
299  namespace vigra {
300  template <class SrcIterator, class SrcShape, class SrcAccessor,
301  class DestIterator, class DestAccessor>
302  void
303  tensorTraceMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
304  DestIterator di, DestAccessor dest);
305  }
306  \endcode
307 
308 
309  use argument objects in conjunction with \ref ArgumentObjectFactories :
310  \code
311  namespace vigra {
312  template <class SrcIterator, class SrcShape, class SrcAccessor,
313  class DestIterator, class DestAccessor>
314  void
315  tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
316  pair<DestIterator, DestAccessor> d);
317  }
318  \endcode
319 
320  <b> Usage:</b>
321 
322  <b>\#include</b> <vigra/multi_tensorutilities.hxx>
323 
324  \code
325  MultiArray<3, float> vol(shape);
326  MultiArray<3, TinyVector<float, 6> > hessian(shape);
327  MultiArray<3, float> trace(shape);
328 
329  hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
330  tensorTraceMultiArray(srcMultiArrayRange(hessian), destMultiArray(trace));
331  \endcode
332 
333 */
334 doxygen_overloaded_function(template <...> void tensorTraceMultiArray)
335 
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  static const int N = SrcShape::static_size;
343  typedef typename SrcAccessor::value_type SrcType;
344 
345  transformMultiArray(si, shape, src, di, dest,
346  detail::TensorTraceFunctor<N, SrcType>());
347 }
348 
349 template <class SrcIterator, class SrcShape, class SrcAccessor,
350  class DestIterator, class DestAccessor>
351 inline void
352 tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
353  pair<DestIterator, DestAccessor> d)
354 {
355  tensorTraceMultiArray(s.first, s.second, s.third, d.first, d.second);
356 }
357 
358 
359 /********************************************************/
360 /* */
361 /* tensorEigenvaluesMultiArray */
362 /* */
363 /********************************************************/
364 
365 /** \brief Calculate the tensor eigenvalues for every element of a N-D tensor array.
366 
367  This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
368  see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
369  symmetric tensor into a vector-valued array holding the tensor eigenvalues (thus,
370  the destination value_type must be vectors of length N).
371 
372  Currently, <tt>N <= 3</tt> is required.
373 
374  <b> Declarations:</b>
375 
376  pass arguments explicitly:
377  \code
378  namespace vigra {
379  template <class SrcIterator, class SrcShape, class SrcAccessor,
380  class DestIterator, class DestAccessor>
381  void
382  tensorEigenvaluesMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
383  DestIterator di, DestAccessor dest);
384  }
385  \endcode
386 
387 
388  use argument objects in conjunction with \ref ArgumentObjectFactories :
389  \code
390  namespace vigra {
391  template <class SrcIterator, class SrcShape, class SrcAccessor,
392  class DestIterator, class DestAccessor>
393  void
394  tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
395  pair<DestIterator, DestAccessor> d);
396  }
397  \endcode
398 
399  <b> Usage:</b>
400 
401  <b>\#include</b> <vigra/multi_tensorutilities.hxx>
402 
403  \code
404  MultiArray<3, float> vol(shape);
405  MultiArray<3, TinyVector<float, 6> > hessian(shape);
406  MultiArray<3, TinyVector<float, 3> > eigenvalues(shape);
407 
408  hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
409  tensorEigenvaluesMultiArray(srcMultiArrayRange(hessian), destMultiArray(eigenvalues));
410  \endcode
411 
412 */
413 doxygen_overloaded_function(template <...> void tensorEigenvaluesMultiArray)
414 
415 template <class SrcIterator, class SrcShape, class SrcAccessor,
416  class DestIterator, class DestAccessor>
417 void
418 tensorEigenvaluesMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
419  DestIterator di, DestAccessor dest)
420 {
421  static const int N = SrcShape::static_size;
422  static const int M = N*(N+1)/2;
423 
424  typedef typename SrcAccessor::value_type SrcType;
425  typedef typename DestAccessor::value_type DestType;
426 
427  for(int k=0; k<N; ++k)
428  if(shape[k] <=0)
429  return;
430 
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.");
435 
436  transformMultiArray(si, shape, src, di, dest,
437  detail::EigenvaluesFunctor<N, SrcType, DestType>());
438 }
439 
440 template <class SrcIterator, class SrcShape, class SrcAccessor,
441  class DestIterator, class DestAccessor>
442 inline void
443 tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
444  pair<DestIterator, DestAccessor> d)
445 {
446  tensorEigenvaluesMultiArray(s.first, s.second, s.third, d.first, d.second);
447 }
448 
449 /********************************************************/
450 /* */
451 /* tensorDeterminantMultiArray */
452 /* */
453 /********************************************************/
454 
455 /** \brief Calculate the tensor determinant for every element of a ND tensor array.
456 
457  This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2,
458  see \ref vectorToTensorMultiArray()) representing the upper triangular part of a
459  symmetric tensor into the a scalar array holding the tensor determinant.
460 
461  Currently, <tt>N <= 3</tt> is required.
462 
463  <b> Declarations:</b>
464 
465  pass arguments explicitly:
466  \code
467  namespace vigra {
468  template <class SrcIterator, class SrcShape, class SrcAccessor,
469  class DestIterator, class DestAccessor>
470  void
471  tensorDeterminantMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
472  DestIterator di, DestAccessor dest);
473  }
474  \endcode
475 
476 
477  use argument objects in conjunction with \ref ArgumentObjectFactories :
478  \code
479  namespace vigra {
480  template <class SrcIterator, class SrcShape, class SrcAccessor,
481  class DestIterator, class DestAccessor>
482  void
483  tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
484  pair<DestIterator, DestAccessor> d);
485  }
486  \endcode
487 
488  <b> Usage:</b>
489 
490  <b>\#include</b> <vigra/multi_tensorutilities.hxx>
491 
492  \code
493  MultiArray<3, float> vol(shape);
494  MultiArray<3, TinyVector<float, 6> > hessian(shape);
495  MultiArray<3, float> determinant(shape);
496 
497  hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
498  tensorDeterminantMultiArray(srcMultiArrayRange(hessian), destMultiArray(determinant));
499  \endcode
500 
501 */
502 doxygen_overloaded_function(template <...> void tensorDeterminantMultiArray)
503 
504 template <class SrcIterator, class SrcShape, class SrcAccessor,
505  class DestIterator, class DestAccessor>
506 void
507 tensorDeterminantMultiArray(SrcIterator si, SrcShape const & shape, SrcAccessor src,
508  DestIterator di, DestAccessor dest)
509 {
510  typedef typename SrcAccessor::value_type SrcType;
511 
512  static const int N = SrcShape::static_size;
513  static const int M = N*(N+1)/2;
514 
515  for(int k=0; k<N; ++k)
516  if(shape[k] <=0)
517  return;
518 
519  vigra_precondition(M == (int)src.size(si),
520  "tensorDeterminantMultiArray(): Wrong number of channels in output array.");
521 
522  transformMultiArray(si, shape, src, di, dest,
523  detail::DeterminantFunctor<N, SrcType>());
524 }
525 
526 template <class SrcIterator, class SrcShape, class SrcAccessor,
527  class DestIterator, class DestAccessor>
528 inline void
529 tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
530  pair<DestIterator, DestAccessor> d)
531 {
532  tensorDeterminantMultiArray(s.first, s.second, s.third, d.first, d.second);
533 }
534 
535 //@}
536 
537 } // namespace vigra
538 
539 #endif /* VIGRA_MULTI_TENSORUTILITIES_HXX */

© 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.9.0 (Wed Feb 27 2013)