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

multi_math.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2010-2011 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_MATH_HXX
37 #define VIGRA_MULTI_MATH_HXX
38 
39 #include "multi_array.hxx"
40 #include "tinyvector.hxx"
41 #include "rgbvalue.hxx"
42 #include "mathutil.hxx"
43 #include <complex>
44 
45 namespace vigra {
46 
47 /** \defgroup MultiMathModule vigra::multi_math
48 
49  Namespace <tt>vigra::multi_math</tt> holds VIGRA's support for efficient arithmetic and algebraic functions on multi-dimensional arrays (that is, \ref MultiArrayView and its subclasses). All <tt>multi_math</tt> functions operate element-wise. If you need matrix multiplication, use \ref LinearAlgebraModule instead.
50 
51  In order to avoid overload ambiguities, multi-array arithmetic must be explicitly activated by
52  \code
53  using namespace vigra::multi_math;
54  \endcode
55  (this should not be done globally, but only in the scope where the functionality is actually used).
56 
57  You can then use the standard operators in the expected way:
58  \code
59  MultiArray<2, float> i(Shape2(100, 100)), j(Shape2(100, 100));
60 
61  MultiArray<2, float> h = i + 4.0 * j;
62  h += (i.transpose() - j) / 2.0;
63  \endcode
64  etc. (supported operators are <tt>+ - * / ! ~ % && || == != &lt; &lt;= &gt; &gt;= &lt;&lt; &gt;&gt; & | ^ = += -= *= /=</tt>, with both scalar and array arguments).
65 
66  Algebraic functions are available as well:
67  \code
68  h = exp(-(sq(i) + sq(j)));
69  h *= atan2(-i, j);
70  \endcode
71  The following functions are implemented: <tt>abs, erf, even, odd, sign, signi, round, roundi, sqrt, sqrti, sq,
72  norm, squaredNorm, gamma, loggamma, exp, log, log10, sin, sin_pi, cos, cos_pi, asin, acos, tan, atan,
73  floor, ceil, conj, real, imag, arg, atan2, pow, fmod, min, max</tt>,
74  provided the array's element type supports the respective function.
75 
76  Supported element types currently include the built-in numeric types, \ref TinyVector, \ref RGBValue,
77  <tt>std::complex</tt>, and \ref FFTWComplex.
78 
79  In addition, <tt>multi_math</tt> supports a number of functions that reduce arrays to scalars:
80  \code
81  double s = sum<double>(i); // compute the sum of the elements, using 'double' as accumulator type
82  double p = product<double>(abs(i)); // compute the product of the elements' absolute values
83 
84  bool a = any(i < 0.0); // check if any element of i is negative
85  bool b = all(i > 0.0); // check if all elements of i are positive
86  \endcode
87 
88  Expressions are expanded so that no temporary arrays have to be created. To optimize cache locality,
89  loops are executed in the stride ordering of the left-hand-side array.
90 
91  <b>\#include</b> <vigra/multi_math.hxx>
92 
93  Namespace: vigra::multi_math
94 */
95 namespace multi_math {
96 
97 template <class ARG>
98 struct MultiMathOperand
99 {
100  typedef typename ARG::result_type result_type;
101 
102  static const int ndim = ARG::ndim;
103 
104  MultiMathOperand(ARG const & a)
105  : arg_(a)
106  {}
107 
108  // Check if all arrays involved in the expression have compatible shapes
109  // (including transparent expansion of singleton axes).
110  // 's' is the shape of the LHS array. If 's' is zero (i.e. the LHS is
111  // not yet initialized), it is set to the maximal RHS shape.
112  //
113  template <class SHAPE>
114  bool checkShape(SHAPE & s) const
115  {
116  return arg_.checkShape(s);
117  }
118 
119  // increment the pointer of all RHS arrays along the given 'axis'
120  void inc(unsigned int axis) const
121  {
122  arg_.inc(axis);
123  }
124 
125  // reset the pointer of all RHS arrays along the given 'axis'
126  void reset(unsigned int axis) const
127  {
128  arg_.reset(axis);
129  }
130 
131  // get the value of the expression at the current pointer location
132  result_type operator*() const
133  {
134  return *arg_;
135  }
136 
137  // get the value of the expression at an offset of the current pointer location
138  template <class SHAPE>
139  result_type operator[](SHAPE const & s) const
140  {
141  return arg_[s];
142  }
143 
144  ARG arg_;
145 };
146 
147 template <unsigned int N, class T, class C>
148 struct MultiMathOperand<MultiArrayView<N, T, C> >
149 {
150  typedef MultiMathOperand AllowOverload;
151  typedef typename MultiArrayShape<N>::type Shape;
152 
153  typedef T result_type;
154 
155  static const int ndim = (int)N;
156 
157  MultiMathOperand(MultiArrayView<N, T, C> const & a)
158  : p_(a.data()),
159  shape_(a.shape()),
160  strides_(a.stride())
161  {
162  // allow for transparent expansion of singleton dimensions
163  for(unsigned int k=0; k<N; ++k)
164  if(shape_[k] == 1)
165  strides_[k] = 0;
166  }
167 
168  bool checkShape(Shape & s) const
169  {
170  // support:
171  // * transparent expansion of singleton dimensions
172  // * determining LHS shape in a constructor
173  for(unsigned int k=0; k<N; ++k)
174  {
175  if(shape_[k] == 0)
176  {
177  return false;
178  }
179  else if(s[k] <= 1)
180  {
181  s[k] = shape_[k];
182  }
183  else if(shape_[k] > 1 && shape_[k] != s[k])
184  {
185  return false;
186  }
187  }
188  return true;
189  }
190 
191  T const & operator[](Shape const & s) const
192  {
193  return p_[dot(s, strides_)];
194  }
195 
196  void inc(unsigned int axis) const
197  {
198  p_ += strides_[axis];
199  }
200 
201  void reset(unsigned int axis) const
202  {
203  p_ -= shape_[axis]*strides_[axis];
204  }
205 
206  result_type operator*() const
207  {
208  return *p_;
209  }
210 
211  mutable T const * p_;
212  Shape shape_, strides_;
213 };
214 
215 template <unsigned int N, class T, class A>
216 struct MultiMathOperand<MultiArray<N, T, A> >
217 : public MultiMathOperand<MultiArrayView<N, T, UnstridedArrayTag> >
218 {
219  typedef MultiMathOperand AllowOverload;
220 
221  MultiMathOperand(MultiArray<N, T, A> const & a)
222  : MultiMathOperand<MultiArrayView<N, T, UnstridedArrayTag> >(a)
223  {}
224 };
225 
226 template <class T>
227 struct MultiMathScalarOperand
228 {
229  typedef MultiMathOperand<T> AllowOverload;
230  typedef T result_type;
231 
232  static const int ndim = 0;
233 
234  MultiMathScalarOperand(T const & v)
235  : v_(v)
236  {}
237 
238  template <class SHAPE>
239  bool checkShape(SHAPE const &) const
240  {
241  return true;
242  }
243 
244  template <class SHAPE>
245  T const & operator[](SHAPE const &) const
246  {
247  return v_;
248  }
249 
250  void inc(unsigned int /* axis */) const
251  {}
252 
253  void reset(unsigned int /* axis */) const
254  {}
255 
256  T const & operator*() const
257  {
258  return v_;
259  }
260 
261  T v_;
262 };
263 
264 #define VIGRA_CONSTANT_OPERAND(template_dcl, type) \
265 template template_dcl \
266 struct MultiMathOperand<type > \
267 : MultiMathScalarOperand<type > \
268 { \
269  MultiMathOperand(type const & v) \
270  : MultiMathScalarOperand<type >(v) \
271  {} \
272 };
273 
274 VIGRA_CONSTANT_OPERAND(<>, signed char)
275 VIGRA_CONSTANT_OPERAND(<>, signed short)
276 VIGRA_CONSTANT_OPERAND(<>, signed int)
277 VIGRA_CONSTANT_OPERAND(<>, signed long)
278 VIGRA_CONSTANT_OPERAND(<>, signed long long)
279 VIGRA_CONSTANT_OPERAND(<>, unsigned char)
280 VIGRA_CONSTANT_OPERAND(<>, unsigned short)
281 VIGRA_CONSTANT_OPERAND(<>, unsigned int)
282 VIGRA_CONSTANT_OPERAND(<>, unsigned long)
283 VIGRA_CONSTANT_OPERAND(<>, unsigned long long)
284 VIGRA_CONSTANT_OPERAND(<>, float)
285 VIGRA_CONSTANT_OPERAND(<>, double)
286 VIGRA_CONSTANT_OPERAND(<>, long double)
287 VIGRA_CONSTANT_OPERAND(<class T>, std::complex<T>)
288 
289 #define VIGRA_TINYVECTOR_ARGS <class T, int N>
290 #define VIGRA_TINYVECTOR_DECL TinyVector<T, N>
291 VIGRA_CONSTANT_OPERAND(VIGRA_TINYVECTOR_ARGS, VIGRA_TINYVECTOR_DECL)
292 #undef VIGRA_TINYVECTOR_ARGS
293 #undef VIGRA_TINYVECTOR_DECL
294 
295 #define VIGRA_RGBVALUE_ARGS <class V, unsigned int R, unsigned int G, unsigned int B>
296 #define VIGRA_RGBVALUE_DECL RGBValue<V, R, G, B>
297 VIGRA_CONSTANT_OPERAND(VIGRA_RGBVALUE_ARGS, VIGRA_RGBVALUE_DECL)
298 #undef VIGRA_RGBVALUE_ARGS
299 #undef VIGRA_RGBVALUE_DECL
300 
301 #undef VIGRA_CONSTANT_OPERAND
302 
303 template <class O, class F>
304 struct MultiMathUnaryOperator
305 {
306  typedef typename F::template Result<typename O::result_type>::type result_type;
307 
308  static const int ndim = O::ndim;
309 
310  MultiMathUnaryOperator(O const & o)
311  : o_(o)
312  {}
313 
314  template <class SHAPE>
315  bool checkShape(SHAPE & s) const
316  {
317  return o_.checkShape(s);
318  }
319 
320  //
321  void inc(unsigned int axis) const
322  {
323  o_.inc(axis);
324  }
325 
326  void reset(unsigned int axis) const
327  {
328  o_.reset(axis);
329  }
330 
331  template <class POINT>
332  result_type operator[](POINT const & p) const
333  {
334  return f_(o_[p]);
335  }
336 
337  result_type operator*() const
338  {
339  return f_(*o_);
340  }
341 
342  O o_;
343  F f_;
344 };
345 
346 #define VIGRA_MULTIMATH_UNARY_OPERATOR(NAME, FCT, OPNAME, RESTYPE) \
347 namespace math_detail { \
348 struct NAME \
349 { \
350  template <class T> \
351  struct Result \
352  { \
353  typedef RESTYPE type; \
354  }; \
355  \
356  template <class T> \
357  typename Result<T>::type \
358  operator()(T const & t) const \
359  { \
360  return FCT(t); \
361  } \
362 }; \
363 } \
364  \
365 template <unsigned int N, class T, class C> \
366 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<MultiArrayView<N, T, C> >, \
367  math_detail::NAME> > \
368 OPNAME(MultiArrayView<N, T, C> const & v) \
369 { \
370  typedef MultiMathOperand<MultiArrayView<N, T, C> > O; \
371  typedef MultiMathUnaryOperator<O, math_detail::NAME> OP; \
372  return MultiMathOperand<OP>(OP(v)); \
373 } \
374  \
375 template <unsigned int N, class T, class A> \
376 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<MultiArray<N, T, A> >, \
377  math_detail::NAME> > \
378 OPNAME(MultiArray<N, T, A> const & v) \
379 { \
380  typedef MultiMathOperand<MultiArray<N, T, A> > O; \
381  typedef MultiMathUnaryOperator<O, math_detail::NAME> OP; \
382  return MultiMathOperand<OP>(OP(v)); \
383 } \
384  \
385 template <class T> \
386 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<T>, \
387  math_detail::NAME> > \
388 OPNAME(MultiMathOperand<T> const & v) \
389 { \
390  typedef MultiMathOperand<T> O; \
391  typedef MultiMathUnaryOperator<O, math_detail::NAME> OP; \
392  return MultiMathOperand<OP>(OP(v)); \
393 }
394 
395 #define VIGRA_REALPROMOTE typename NumericTraits<T>::RealPromote
396 
397 #ifndef DOXYGEN // doxygen gets confused by these macros
398 
399 VIGRA_MULTIMATH_UNARY_OPERATOR(Negate, -, operator-, T)
400 VIGRA_MULTIMATH_UNARY_OPERATOR(Not, !, operator!, T)
401 VIGRA_MULTIMATH_UNARY_OPERATOR(BitwiseNot, ~, operator~, T)
402 
403 using vigra::abs;
404 VIGRA_MULTIMATH_UNARY_OPERATOR(Abs, vigra::abs, abs, typename NormTraits<T>::NormType)
405 
406 using vigra::erf;
407 VIGRA_MULTIMATH_UNARY_OPERATOR(Erf, vigra::erf, erf, VIGRA_REALPROMOTE)
408 VIGRA_MULTIMATH_UNARY_OPERATOR(Even, vigra::even, even, bool)
409 VIGRA_MULTIMATH_UNARY_OPERATOR(Odd, vigra::odd, odd, bool)
410 VIGRA_MULTIMATH_UNARY_OPERATOR(Sign, vigra::sign, sign, T)
411 VIGRA_MULTIMATH_UNARY_OPERATOR(Signi, vigra::signi, signi, int)
412 
413 using vigra::round;
414 VIGRA_MULTIMATH_UNARY_OPERATOR(Round, vigra::round, round, T)
415 
416 VIGRA_MULTIMATH_UNARY_OPERATOR(Roundi, vigra::roundi, roundi, int)
417 VIGRA_MULTIMATH_UNARY_OPERATOR(Sqrti, vigra::sqrti, sqrti, T)
418 VIGRA_MULTIMATH_UNARY_OPERATOR(Sq, vigra::sq, sq, typename NumericTraits<T>::Promote)
419 VIGRA_MULTIMATH_UNARY_OPERATOR(Norm, vigra::norm, norm, typename NormTraits<T>::NormType)
420 VIGRA_MULTIMATH_UNARY_OPERATOR(SquaredNorm, vigra::squaredNorm, squaredNorm,
421  typename NormTraits<T>::SquaredNormType)
422 VIGRA_MULTIMATH_UNARY_OPERATOR(Sin_pi, vigra::sin_pi, sin_pi, VIGRA_REALPROMOTE)
423 VIGRA_MULTIMATH_UNARY_OPERATOR(Cos_pi, vigra::cos_pi, cos_pi, VIGRA_REALPROMOTE)
424 
425 using vigra::gamma;
426 VIGRA_MULTIMATH_UNARY_OPERATOR(Gamma, vigra::gamma, gamma, VIGRA_REALPROMOTE)
427 
428 using vigra::loggamma;
429 VIGRA_MULTIMATH_UNARY_OPERATOR(Loggamma, vigra::loggamma, loggamma, VIGRA_REALPROMOTE)
430 
431 VIGRA_MULTIMATH_UNARY_OPERATOR(Sqrt, std::sqrt, sqrt, VIGRA_REALPROMOTE)
432 using vigra::exp;
433 VIGRA_MULTIMATH_UNARY_OPERATOR(Exp, vigra::exp, exp, VIGRA_REALPROMOTE)
434 VIGRA_MULTIMATH_UNARY_OPERATOR(Log, std::log, log, VIGRA_REALPROMOTE)
435 VIGRA_MULTIMATH_UNARY_OPERATOR(Log10, std::log10, log10, VIGRA_REALPROMOTE)
436 VIGRA_MULTIMATH_UNARY_OPERATOR(Sin, std::sin, sin, VIGRA_REALPROMOTE)
437 VIGRA_MULTIMATH_UNARY_OPERATOR(Asin, std::asin, asin, VIGRA_REALPROMOTE)
438 VIGRA_MULTIMATH_UNARY_OPERATOR(Cos, std::cos, cos, VIGRA_REALPROMOTE)
439 VIGRA_MULTIMATH_UNARY_OPERATOR(Acos, std::acos, acos, VIGRA_REALPROMOTE)
440 VIGRA_MULTIMATH_UNARY_OPERATOR(Tan, std::tan, tan, VIGRA_REALPROMOTE)
441 VIGRA_MULTIMATH_UNARY_OPERATOR(Atan, std::atan, atan, VIGRA_REALPROMOTE)
442 VIGRA_MULTIMATH_UNARY_OPERATOR(Floor, std::floor, floor, VIGRA_REALPROMOTE)
443 VIGRA_MULTIMATH_UNARY_OPERATOR(Ceil, std::ceil, ceil, VIGRA_REALPROMOTE)
444 
445 VIGRA_MULTIMATH_UNARY_OPERATOR(Conj, conj, conj, T)
446 VIGRA_MULTIMATH_UNARY_OPERATOR(Real, real, real, typename T::value_type)
447 VIGRA_MULTIMATH_UNARY_OPERATOR(Imag, imag, imag, typename T::value_type)
448 VIGRA_MULTIMATH_UNARY_OPERATOR(Arg, arg, arg, typename T::value_type)
449 
450 #endif //DOXYGEN
451 
452 #undef VIGRA_REALPROMOTE
453 #undef VIGRA_MULTIMATH_UNARY_OPERATOR
454 
455 template <class O1, class O2, class F>
456 struct MultiMathBinaryOperator
457 {
458  typedef typename F::template Result<typename O1::result_type,
459  typename O2::result_type>::type result_type;
460 
461  static const int ndim = O1::ndim > O2::ndim ? O1::ndim : O2::ndim;
462 
463  MultiMathBinaryOperator(O1 const & o1, O2 const & o2)
464  : o1_(o1),
465  o2_(o2)
466  {}
467 
468  template <class SHAPE>
469  bool checkShape(SHAPE & s) const
470  {
471  return o1_.checkShape(s) && o2_.checkShape(s);
472  }
473 
474  template <class POINT>
475  result_type operator[](POINT const & p) const
476  {
477  return f_(o1_[p], o2_[p]);
478  }
479 
480  void inc(unsigned int axis) const
481  {
482  o1_.inc(axis);
483  o2_.inc(axis);
484  }
485 
486  void reset(unsigned int axis) const
487  {
488  o1_.reset(axis);
489  o2_.reset(axis);
490  }
491 
492  result_type operator*() const
493  {
494  return f_(*o1_, *o2_);
495  }
496 
497  O1 o1_;
498  O2 o2_;
499  F f_;
500 };
501 
502 
503 // In the sequel, the nested type 'MultiMathOperand<T>::AllowOverload'
504 // ensures that template functions only participate in overload
505 // resolution when this type is defined, i.e. when T is a number
506 // or array type. It thus prevents 'ambiguous overload' errors.
507 //
508 #define VIGRA_MULTIMATH_BINARY_OPERATOR(NAME, FCT, OPNAME, SEP, RESTYPE) \
509 \
510 namespace math_detail { \
511 struct NAME \
512 { \
513  template <class T1, class T2> \
514  struct Result \
515  { \
516  typedef RESTYPE type; \
517  }; \
518  \
519  template <class T1, class T2> \
520  typename Result<T1, T2>::type \
521  operator()(T1 const & t1, T2 const & t2) const \
522  { \
523  return FCT(t1 SEP t2); \
524  } \
525 }; \
526 } \
527  \
528 template <unsigned int N, class T1, class A1, class T2, class A2> \
529 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1> >, \
530  MultiMathOperand<MultiArrayView<N, T2> >, \
531  math_detail::NAME> > \
532 OPNAME(MultiArray<N, T1, A1> const & v1, MultiArray<N, T2, A2> const & v2) \
533 { \
534  typedef MultiMathOperand<MultiArrayView<N, T1> > O1; \
535  typedef MultiMathOperand<MultiArrayView<N, T2> > O2; \
536  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
537  return MultiMathOperand<OP>(OP((MultiArrayView<N, T1> const &)v1, (MultiArrayView<N, T2> const &)v2)); \
538 } \
539 \
540 template <unsigned int N, class T1, class C1, class T2, class C2> \
541 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \
542  MultiMathOperand<MultiArrayView<N, T2, C2> >, \
543  math_detail::NAME> > \
544 OPNAME(MultiArrayView<N, T1, C1> const & v1, MultiArrayView<N, T2, C2> const & v2) \
545 { \
546  typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \
547  typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \
548  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
549  return MultiMathOperand<OP>(OP(v1, v2)); \
550 } \
551 \
552 template <unsigned int N, class T1, class T2, class C2> \
553 MultiMathOperand<MultiMathBinaryOperator<typename MultiMathOperand<T1>::AllowOverload, \
554  MultiMathOperand<MultiArrayView<N, T2, C2> >, \
555  math_detail::NAME> > \
556 OPNAME(T1 const & v1, MultiArrayView<N, T2, C2> const & v2) \
557 { \
558  typedef MultiMathOperand<T1> O1; \
559  typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \
560  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
561  return MultiMathOperand<OP>(OP(v1, v2)); \
562 } \
563  \
564 template <unsigned int N, class T1, class C1, class T2> \
565 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \
566  typename MultiMathOperand<T2>::AllowOverload, \
567  math_detail::NAME> > \
568 OPNAME(MultiArrayView<N, T1, C1> const & v1, T2 const & v2) \
569 { \
570  typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \
571  typedef MultiMathOperand<T2> O2; \
572  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
573  return MultiMathOperand<OP>(OP(v1, v2)); \
574 } \
575  \
576 template <unsigned int N, class T1, class T2, class C2> \
577 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \
578  MultiMathOperand<MultiArrayView<N, T2, C2> >, \
579  math_detail::NAME> > \
580 OPNAME(MultiMathOperand<T1> const & v1, MultiArrayView<N, T2, C2> const & v2) \
581 { \
582  typedef MultiMathOperand<T1> O1; \
583  typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \
584  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
585  return MultiMathOperand<OP>(OP(v1, v2)); \
586 } \
587  \
588 template <unsigned int N, class T1, class C1, class T2> \
589 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \
590  MultiMathOperand<T2>, \
591  math_detail::NAME> > \
592 OPNAME(MultiArrayView<N, T1, C1> const & v1, MultiMathOperand<T2> const & v2) \
593 { \
594  typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \
595  typedef MultiMathOperand<T2> O2; \
596  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
597  return MultiMathOperand<OP>(OP(v1, v2)); \
598 } \
599  \
600 template <class T1, class T2> \
601 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \
602  MultiMathOperand<T2>, \
603  math_detail::NAME> > \
604 OPNAME(MultiMathOperand<T1> const & v1, MultiMathOperand<T2> const & v2) \
605 { \
606  typedef MultiMathOperand<T1> O1; \
607  typedef MultiMathOperand<T2> O2; \
608  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
609  return MultiMathOperand<OP>(OP(v1, v2)); \
610 } \
611 \
612 template <class T1, class T2> \
613 MultiMathOperand<MultiMathBinaryOperator<typename MultiMathOperand<T1>::AllowOverload, \
614  MultiMathOperand<T2>, \
615  math_detail::NAME> > \
616 OPNAME(T1 const & v1, MultiMathOperand<T2> const & v2) \
617 { \
618  typedef MultiMathOperand<T1> O1; \
619  typedef MultiMathOperand<T2> O2; \
620  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
621  return MultiMathOperand<OP>(OP(v1, v2)); \
622 } \
623 \
624 template <class T1, class T2> \
625 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \
626  typename MultiMathOperand<T2>::AllowOverload, \
627  math_detail::NAME> > \
628 OPNAME(MultiMathOperand<T1> const & v1, T2 const & v2) \
629 { \
630  typedef MultiMathOperand<T1> O1; \
631  typedef MultiMathOperand<T2> O2; \
632  typedef MultiMathBinaryOperator<O1, O2, math_detail::NAME> OP; \
633  return MultiMathOperand<OP>(OP(v1, v2)); \
634 }
635 
636 #define VIGRA_NOTHING
637 #define VIGRA_COMMA ,
638 #define VIGRA_PROMOTE typename PromoteTraits<T1, T2>::Promote
639 #define VIGRA_REALPROMOTE typename PromoteTraits<typename NumericTraits<T1>::RealPromote, \
640  typename NumericTraits<T2>::RealPromote>::Promote
641 
642 VIGRA_MULTIMATH_BINARY_OPERATOR(Plus, VIGRA_NOTHING, operator+, +, VIGRA_PROMOTE)
643 VIGRA_MULTIMATH_BINARY_OPERATOR(Minus, VIGRA_NOTHING, operator-, -, VIGRA_PROMOTE)
644 VIGRA_MULTIMATH_BINARY_OPERATOR(Multiplies, VIGRA_NOTHING, operator*, *, VIGRA_PROMOTE)
645 VIGRA_MULTIMATH_BINARY_OPERATOR(Divides, VIGRA_NOTHING, operator/, /, VIGRA_PROMOTE)
646 VIGRA_MULTIMATH_BINARY_OPERATOR(Modulo, VIGRA_NOTHING, operator%, %, VIGRA_PROMOTE)
647 VIGRA_MULTIMATH_BINARY_OPERATOR(And, VIGRA_NOTHING, operator&&, &&, VIGRA_PROMOTE)
648 VIGRA_MULTIMATH_BINARY_OPERATOR(Or, VIGRA_NOTHING, operator||, ||, VIGRA_PROMOTE)
649 VIGRA_MULTIMATH_BINARY_OPERATOR(Equal, VIGRA_NOTHING, operator==, ==, bool)
650 VIGRA_MULTIMATH_BINARY_OPERATOR(NotEqual, VIGRA_NOTHING, operator!=, !=, bool)
651 VIGRA_MULTIMATH_BINARY_OPERATOR(Less, VIGRA_NOTHING, operator<, <, bool)
652 VIGRA_MULTIMATH_BINARY_OPERATOR(LessEqual, VIGRA_NOTHING, operator<=, <=, bool)
653 VIGRA_MULTIMATH_BINARY_OPERATOR(Greater, VIGRA_NOTHING, operator>, >, bool)
654 VIGRA_MULTIMATH_BINARY_OPERATOR(GreaterEqual, VIGRA_NOTHING, operator>=, >=, bool)
655 VIGRA_MULTIMATH_BINARY_OPERATOR(Leftshift, VIGRA_NOTHING, operator<<, <<, VIGRA_PROMOTE)
656 VIGRA_MULTIMATH_BINARY_OPERATOR(Rightshift, VIGRA_NOTHING, operator>>, >>, VIGRA_PROMOTE)
657 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseAnd, VIGRA_NOTHING, operator&, &, VIGRA_PROMOTE)
658 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseOr, VIGRA_NOTHING, operator|, |, VIGRA_PROMOTE)
659 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseXor, VIGRA_NOTHING, operator^, ^, VIGRA_PROMOTE)
660 
661 VIGRA_MULTIMATH_BINARY_OPERATOR(Atan2, std::atan2, atan2, VIGRA_COMMA, VIGRA_REALPROMOTE)
662 VIGRA_MULTIMATH_BINARY_OPERATOR(Pow, vigra::pow, pow, VIGRA_COMMA, VIGRA_REALPROMOTE)
663 VIGRA_MULTIMATH_BINARY_OPERATOR(Fmod, std::fmod, fmod, VIGRA_COMMA, VIGRA_REALPROMOTE)
664 VIGRA_MULTIMATH_BINARY_OPERATOR(Min, std::min, min, VIGRA_COMMA, VIGRA_PROMOTE)
665 VIGRA_MULTIMATH_BINARY_OPERATOR(Max, std::max, max, VIGRA_COMMA, VIGRA_PROMOTE)
666 VIGRA_MULTIMATH_BINARY_OPERATOR(Minimum, std::min, minimum, VIGRA_COMMA, VIGRA_PROMOTE)
667 VIGRA_MULTIMATH_BINARY_OPERATOR(Maximum, std::max, maximum, VIGRA_COMMA, VIGRA_PROMOTE)
668 
669 #undef VIGRA_NOTHING
670 #undef VIGRA_COMMA
671 #undef VIGRA_PROMOTE
672 #undef VIGRA_REALPROMOTE
673 #undef VIGRA_MULTIMATH_BINARY_OPERATOR
674 
675 namespace math_detail {
676 
677 // We pass 'strideOrder' to the recursion in order to make sure
678 // that the inner loop iterates over the output's major axis.
679 // Of course, this does not help when the RHS arrays are ordered
680 // differently -- maybe it is better to find the most common order
681 // among all arguments (both RHS and LHS)?
682 //
683 template <unsigned int N, class Assign>
684 struct MultiMathExec
685 {
686  enum { LEVEL = N-1 };
687 
688  template <class T, class Shape, class Expression>
689  static void exec(T * data, Shape const & shape, Shape const & strides,
690  Shape const & strideOrder, Expression const & e)
691  {
692  MultiArrayIndex axis = strideOrder[LEVEL];
693  for(MultiArrayIndex k=0; k<shape[axis]; ++k, data += strides[axis], e.inc(axis))
694  {
695  MultiMathExec<N-1, Assign>::exec(data, shape, strides, strideOrder, e);
696  }
697  e.reset(axis);
698  data -= shape[axis]*strides[axis];
699  }
700 };
701 
702 template <class Assign>
703 struct MultiMathExec<1, Assign>
704 {
705  enum { LEVEL = 0 };
706 
707  template <class T, class Shape, class Expression>
708  static void exec(T * data, Shape const & shape, Shape const & strides,
709  Shape const & strideOrder, Expression const & e)
710  {
711  MultiArrayIndex axis = strideOrder[LEVEL];
712  for(MultiArrayIndex k=0; k<shape[axis]; ++k, data += strides[axis], e.inc(axis))
713  {
714  Assign::assign(data, e);
715  }
716  e.reset(axis);
717  data -= shape[axis]*strides[axis];
718  }
719 };
720 
721 #define VIGRA_MULTIMATH_ASSIGN(NAME, OP) \
722 struct MultiMath##NAME \
723 { \
724  template <class T, class Expression> \
725  static void assign(T * data, Expression const & e) \
726  { \
727  *data OP vigra::detail::RequiresExplicitCast<T>::cast(*e); \
728  } \
729 }; \
730  \
731 template <unsigned int N, class T, class C, class Expression> \
732 void NAME(MultiArrayView<N, T, C> a, MultiMathOperand<Expression> const & e) \
733 { \
734  typename MultiArrayShape<N>::type shape(a.shape()); \
735  \
736  vigra_precondition(e.checkShape(shape), \
737  "multi_math: shape mismatch in expression."); \
738  \
739  MultiMathExec<N, MultiMath##NAME>::exec(a.data(), a.shape(), a.stride(), \
740  a.strideOrdering(), e); \
741 } \
742  \
743 template <unsigned int N, class T, class A, class Expression> \
744 void NAME##OrResize(MultiArray<N, T, A> & a, MultiMathOperand<Expression> const & e) \
745 { \
746  typename MultiArrayShape<N>::type shape(a.shape()); \
747  \
748  vigra_precondition(e.checkShape(shape), \
749  "multi_math: shape mismatch in expression."); \
750  \
751  if(a.size() == 0) \
752  a.reshape(shape); \
753  \
754  MultiMathExec<N, MultiMath##NAME>::exec(a.data(), a.shape(), a.stride(), \
755  a.strideOrdering(), e); \
756 }
757 
758 VIGRA_MULTIMATH_ASSIGN(assign, =)
759 VIGRA_MULTIMATH_ASSIGN(plusAssign, +=)
760 VIGRA_MULTIMATH_ASSIGN(minusAssign, -=)
761 VIGRA_MULTIMATH_ASSIGN(multiplyAssign, *=)
762 VIGRA_MULTIMATH_ASSIGN(divideAssign, /=)
763 
764 #undef VIGRA_MULTIMATH_ASSIGN
765 
766 template <unsigned int N, class Assign>
767 struct MultiMathReduce
768 {
769  enum { LEVEL = N-1 };
770 
771  template <class T, class Shape, class Expression>
772  static void exec(T & t, Shape const & shape, Expression const & e)
773  {
774  for(MultiArrayIndex k=0; k<shape[LEVEL]; ++k, e.inc(LEVEL))
775  {
776  MultiMathReduce<N-1, Assign>::exec(t, shape, e);
777  }
778  e.reset(LEVEL);
779  }
780 };
781 
782 template <class Assign>
783 struct MultiMathReduce<1, Assign>
784 {
785  enum { LEVEL = 0 };
786 
787  template <class T, class Shape, class Expression>
788  static void exec(T & t, Shape const & shape, Expression const & e)
789  {
790  for(MultiArrayIndex k=0; k<shape[0]; ++k, e.inc(0))
791  {
792  Assign::assign(&t, e);
793  }
794  e.reset(0);
795  }
796 };
797 
798 struct MultiMathReduceAll
799 {
800  template <class T, class Expression>
801  static void assign(T * data, Expression const & e)
802  {
803  *data = *data && (*e != NumericTraits<typename Expression::result_type>::zero());
804  }
805 };
806 
807 struct MultiMathReduceAny
808 {
809  template <class T, class Expression>
810  static void assign(T * data, Expression const & e)
811  {
812  *data = *data || (*e != NumericTraits<typename Expression::result_type>::zero());
813  }
814 };
815 
816 
817 } // namespace math_detail
818 
819 template <class U, class T>
820 U
821 sum(MultiMathOperand<T> const & v, U res = NumericTraits<U>::zero())
822 {
823  static const int ndim = MultiMathOperand<T>::ndim;
824  typename MultiArrayShape<ndim>::type shape;
825  v.checkShape(shape);
826  math_detail::MultiMathReduce<ndim, math_detail::MultiMathplusAssign>::exec(res, shape, v);
827  return res;
828 }
829 
830 template <class U, unsigned int N, class T, class S>
831 U
832 sum(MultiArrayView<N, T, S> const & v, U res = NumericTraits<U>::zero())
833 {
834  return v.template sum<U>() + res;
835 }
836 
837 template <class U, class T>
838 U
839 product(MultiMathOperand<T> const & v, U res = NumericTraits<U>::one())
840 {
841  static const int ndim = MultiMathOperand<T>::ndim;
842  typename MultiArrayShape<ndim>::type shape;
843  v.checkShape(shape);
844  math_detail::MultiMathReduce<ndim, math_detail::MultiMathmultiplyAssign>::exec(res, shape, v);
845  return res;
846 }
847 
848 template <class U, unsigned int N, class T, class S>
849 U
850 product(MultiArrayView<N, T, S> const & v, U res = NumericTraits<U>::one())
851 {
852  return v.template product<U>() * res;
853 }
854 
855 template <class T>
856 bool
857 all(MultiMathOperand<T> const & v)
858 {
859  static const int ndim = MultiMathOperand<T>::ndim;
860  typename MultiArrayShape<ndim>::type shape;
861  v.checkShape(shape);
862  bool res = true;
863  math_detail::MultiMathReduce<ndim, math_detail::MultiMathReduceAll>::exec(res, shape, v);
864  return res;
865 }
866 
867 template <class T>
868 bool
869 any(MultiMathOperand<T> const & v)
870 {
871  static const int ndim = MultiMathOperand<T>::ndim;
872  typename MultiArrayShape<ndim>::type shape;
873  v.checkShape(shape);
874  bool res = false;
875  math_detail::MultiMathReduce<ndim, math_detail::MultiMathReduceAny>::exec(res, shape, v);
876  return res;
877 }
878 
879 
880 }} // namespace vigra::multi_math
881 
882 #endif // VIGRA_MULTI_MATH_HXX
REAL sin_pi(REAL x)
sin(pi*x).
Definition: mathutil.hxx:1166
linalg::TemporaryMatrix< T > acos(MultiArrayView< 2, T, C > const &v)
FixedPoint16< 2, OverflowHandling > atan2(FixedPoint16< IntBits, OverflowHandling > y, FixedPoint16< IntBits, OverflowHandling > x)
Arctangent. Accuracy better than 1/3 degree (9 significant bits).
Definition: fixedpoint.hxx:1654
PromoteTraits< V1, V2 >::Promote dot(RGBValue< V1, RIDX1, GIDX1, BIDX1 > const &r1, RGBValue< V2, RIDX2, GIDX2, BIDX2 > const &r2)
dot product
Definition: rgbvalue.hxx:907
FFTWComplex< R > conj(const FFTWComplex< R > &a)
complex conjugate
Definition: fftw3.hxx:1030
R imag(const FFTWComplex< R > &a)
imaginary part
Definition: fftw3.hxx:1023
R arg(const FFTWComplex< R > &a)
pahse
Definition: fftw3.hxx:1009
Int32 roundi(FixedPoint16< IntBits, OverflowHandling > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:1775
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
linalg::TemporaryMatrix< T > exp(MultiArrayView< 2, T, C > const &v)
FFTWComplex< R >::SquaredNormType squaredNorm(const FFTWComplex< R > &a)
squared norm (= squared magnitude)
Definition: fftw3.hxx:1044
R real(const FFTWComplex< R > &a)
real part
Definition: fftw3.hxx:1016
int round(FixedPoint< IntBits, FracBits > v)
rounding to the nearest integer.
Definition: fixedpoint.hxx:683
linalg::TemporaryMatrix< T > asin(MultiArrayView< 2, T, C > const &v)
int signi(T t)
The integer sign function.
Definition: mathutil.hxx:570
std::ptrdiff_t MultiArrayIndex
Definition: multi_shape.hxx:55
bool odd(int t)
Check if an integer is odd.
REAL cos_pi(REAL x)
cos(pi*x).
Definition: mathutil.hxx:1204
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
NumericTraits< T >::Promote sq(T t)
The square function.
Definition: mathutil.hxx:344
linalg::TemporaryMatrix< T > log10(MultiArrayView< 2, T, C > const &v)
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector&#39;s elements
Definition: tinyvector.hxx:1871
bool even(int t)
Check if an integer is even.
double gamma(double x)
The gamma function.
Definition: mathutil.hxx:1549
linalg::TemporaryMatrix< T > log(MultiArrayView< 2, T, C > const &v)
double loggamma(double x)
The natural logarithm of the gamma function.
Definition: mathutil.hxx:1565
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
linalg::TemporaryMatrix< T > atan(MultiArrayView< 2, T, C > const &v)
linalg::TemporaryMatrix< T > tan(MultiArrayView< 2, T, C > const &v)
Int32 sqrti(Int32 v)
Signed integer square root.
Definition: mathutil.hxx:500
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)
T sign(T t)
The sign function.
Definition: mathutil.hxx:553
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667
SquareRootTraits< FixedPoint< IntBits, FracBits > >::SquareRootResult sqrt(FixedPoint< IntBits, FracBits > v)
square root.
Definition: fixedpoint.hxx:616

© 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.10.0 (Tue Dec 31 2013)