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

multi_math.hxx |

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

© 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