[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/multi_math.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 2010-2011 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_MULTI_MATH_HXX 00037 #define VIGRA_MULTI_MATH_HXX 00038 00039 #include "multi_array.hxx" 00040 #include "tinyvector.hxx" 00041 #include "rgbvalue.hxx" 00042 #include "mathutil.hxx" 00043 #include <complex> 00044 00045 namespace vigra { 00046 00047 /** \defgroup MultiMathModule vigra::multi_math 00048 00049 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. 00050 00051 In order to avoid overload ambiguities, multi-array arithmetic must be explicitly activated by 00052 \code 00053 using namespace vigra::multi_math; 00054 \endcode 00055 (this should not be done globally, but only in the scope where the functionality is actually used). 00056 00057 You can then use the standard operators in the expected way: 00058 \code 00059 MultiArray<2, float> i(Shape2(100, 100)), j(Shape2(100, 100)); 00060 00061 MultiArray<2, float> h = i + 4.0 * j; 00062 h += (i.transpose() - j) / 2.0; 00063 \endcode 00064 etc. (supported operators are <tt>+ - * / ! ~ % && || == != < <= > >= << >> & | ^ = += -= *= /=</tt>, with both scalar and array arguments). 00065 00066 Algebraic functions are available as well: 00067 \code 00068 h = exp(-(sq(i) + sq(j))); 00069 h *= atan2(-i, j); 00070 \endcode 00071 The following functions are implemented: <tt>abs, erf, even, odd, sign, signi, round, roundi, sqrt, sqrti, sq, 00072 norm, squaredNorm, gamma, loggamma, exp, log, log10, sin, sin_pi, cos, cos_pi, asin, acos, tan, atan, 00073 floor, ceil, conj, real, imag, arg, atan2, pow, fmod, min, max</tt>, 00074 provided the array's element type supports the respective function. 00075 00076 Supported element types currently include the built-in numeric types, \ref TinyVector, \ref RGBValue, 00077 <tt>std::complex</tt>, and \ref FFTWComplex. 00078 00079 In addition, <tt>multi_math</tt> supports a number of functions that reduce arrays to scalars: 00080 \code 00081 double s = sum<double>(i); // compute the sum of the elements, using 'double' as accumulator type 00082 double p = product<double>(abs(i)); // compute the product of the elements' absolute values 00083 00084 bool a = any(i < 0.0); // check if any element of i is negative 00085 bool b = all(i > 0.0); // check if all elements of i are positive 00086 \endcode 00087 00088 Expressions are expanded so that no temporary arrays have to be created. To optimize cache locality, 00089 loops are executed in the stride ordering of the left-hand-side array. 00090 00091 <b>\#include</b> <vigra/multi_math.hxx> 00092 00093 Namespace: vigra::multi_math 00094 */ 00095 namespace multi_math { 00096 00097 template <class ARG> 00098 struct MultiMathOperand 00099 { 00100 typedef typename ARG::result_type result_type; 00101 00102 static const int ndim = ARG::ndim; 00103 00104 MultiMathOperand(ARG const & a) 00105 : arg_(a) 00106 {} 00107 00108 // Check if all arrays involved in the expression have compatible shapes 00109 // (including transparent expansion of singleton axes). 00110 // 's' is the shape of the LHS array. If 's' is zero (i.e. the LHS is 00111 // not yet initialized), it is set to the maximal RHS shape. 00112 // 00113 template <class SHAPE> 00114 bool checkShape(SHAPE & s) const 00115 { 00116 return arg_.checkShape(s); 00117 } 00118 00119 // increment the pointer of all RHS arrays along the given 'axis' 00120 void inc(unsigned int axis) const 00121 { 00122 arg_.inc(axis); 00123 } 00124 00125 // reset the pointer of all RHS arrays along the given 'axis' 00126 void reset(unsigned int axis) const 00127 { 00128 arg_.reset(axis); 00129 } 00130 00131 // get the value of the expression at the current pointer location 00132 result_type operator*() const 00133 { 00134 return *arg_; 00135 } 00136 00137 // get the value of the expression at an offset of the current pointer location 00138 template <class SHAPE> 00139 result_type operator[](SHAPE const & s) const 00140 { 00141 return arg_[s]; 00142 } 00143 00144 ARG arg_; 00145 }; 00146 00147 template <unsigned int N, class T, class C> 00148 struct MultiMathOperand<MultiArrayView<N, T, C> > 00149 { 00150 typedef MultiMathOperand AllowOverload; 00151 typedef typename MultiArrayShape<N>::type Shape; 00152 00153 typedef T result_type; 00154 00155 static const int ndim = (int)N; 00156 00157 MultiMathOperand(MultiArrayView<N, T, C> const & a) 00158 : p_(a.data()), 00159 shape_(a.shape()), 00160 strides_(a.stride()) 00161 { 00162 // allow for transparent expansion of singleton dimensions 00163 for(unsigned int k=0; k<N; ++k) 00164 if(shape_[k] == 1) 00165 strides_[k] = 0; 00166 } 00167 00168 bool checkShape(Shape & s) const 00169 { 00170 // support: 00171 // * transparent expansion of singleton dimensions 00172 // * determining LHS shape in a constructor 00173 for(unsigned int k=0; k<N; ++k) 00174 { 00175 if(shape_[k] == 0) 00176 { 00177 return false; 00178 } 00179 else if(s[k] <= 1) 00180 { 00181 s[k] = shape_[k]; 00182 } 00183 else if(shape_[k] > 1 && shape_[k] != s[k]) 00184 { 00185 return false; 00186 } 00187 } 00188 return true; 00189 } 00190 00191 T const & operator[](Shape const & s) const 00192 { 00193 return p_[dot(s, strides_)]; 00194 } 00195 00196 void inc(unsigned int axis) const 00197 { 00198 p_ += strides_[axis]; 00199 } 00200 00201 void reset(unsigned int axis) const 00202 { 00203 p_ -= shape_[axis]*strides_[axis]; 00204 } 00205 00206 result_type operator*() const 00207 { 00208 return *p_; 00209 } 00210 00211 mutable T const * p_; 00212 Shape shape_, strides_; 00213 }; 00214 00215 template <unsigned int N, class T, class A> 00216 struct MultiMathOperand<MultiArray<N, T, A> > 00217 : public MultiMathOperand<MultiArrayView<N, T, UnstridedArrayTag> > 00218 { 00219 typedef MultiMathOperand AllowOverload; 00220 00221 MultiMathOperand(MultiArray<N, T, A> const & a) 00222 : MultiMathOperand<MultiArrayView<N, T, UnstridedArrayTag> >(a) 00223 {} 00224 }; 00225 00226 template <class T> 00227 struct MultiMathScalarOperand 00228 { 00229 typedef MultiMathOperand<T> AllowOverload; 00230 typedef T result_type; 00231 00232 static const int ndim = 0; 00233 00234 MultiMathScalarOperand(T const & v) 00235 : v_(v) 00236 {} 00237 00238 template <class SHAPE> 00239 bool checkShape(SHAPE const &) const 00240 { 00241 return true; 00242 } 00243 00244 template <class SHAPE> 00245 T const & operator[](SHAPE const &) const 00246 { 00247 return v_; 00248 } 00249 00250 void inc(unsigned int /* axis */) const 00251 {} 00252 00253 void reset(unsigned int /* axis */) const 00254 {} 00255 00256 T const & operator*() const 00257 { 00258 return v_; 00259 } 00260 00261 T v_; 00262 }; 00263 00264 #define VIGRA_CONSTANT_OPERAND(template_dcl, type) \ 00265 template template_dcl \ 00266 struct MultiMathOperand<type > \ 00267 : MultiMathScalarOperand<type > \ 00268 { \ 00269 MultiMathOperand(type const & v) \ 00270 : MultiMathScalarOperand<type >(v) \ 00271 {} \ 00272 }; 00273 00274 VIGRA_CONSTANT_OPERAND(<>, signed char) 00275 VIGRA_CONSTANT_OPERAND(<>, signed short) 00276 VIGRA_CONSTANT_OPERAND(<>, signed int) 00277 VIGRA_CONSTANT_OPERAND(<>, signed long) 00278 VIGRA_CONSTANT_OPERAND(<>, signed long long) 00279 VIGRA_CONSTANT_OPERAND(<>, unsigned char) 00280 VIGRA_CONSTANT_OPERAND(<>, unsigned short) 00281 VIGRA_CONSTANT_OPERAND(<>, unsigned int) 00282 VIGRA_CONSTANT_OPERAND(<>, unsigned long) 00283 VIGRA_CONSTANT_OPERAND(<>, unsigned long long) 00284 VIGRA_CONSTANT_OPERAND(<>, float) 00285 VIGRA_CONSTANT_OPERAND(<>, double) 00286 VIGRA_CONSTANT_OPERAND(<>, long double) 00287 VIGRA_CONSTANT_OPERAND(<class T>, std::complex<T>) 00288 00289 #define VIGRA_TINYVECTOR_ARGS <class T, int N> 00290 #define VIGRA_TINYVECTOR_DECL TinyVector<T, N> 00291 VIGRA_CONSTANT_OPERAND(VIGRA_TINYVECTOR_ARGS, VIGRA_TINYVECTOR_DECL) 00292 #undef VIGRA_TINYVECTOR_ARGS 00293 #undef VIGRA_TINYVECTOR_DECL 00294 00295 #define VIGRA_RGBVALUE_ARGS <class V, unsigned int R, unsigned int G, unsigned int B> 00296 #define VIGRA_RGBVALUE_DECL RGBValue<V, R, G, B> 00297 VIGRA_CONSTANT_OPERAND(VIGRA_RGBVALUE_ARGS, VIGRA_RGBVALUE_DECL) 00298 #undef VIGRA_RGBVALUE_ARGS 00299 #undef VIGRA_RGBVALUE_DECL 00300 00301 #undef VIGRA_CONSTANT_OPERAND 00302 00303 template <class O, class F> 00304 struct MultiMathUnaryOperator 00305 { 00306 typedef typename F::template Result<typename O::result_type>::type result_type; 00307 00308 static const int ndim = O::ndim; 00309 00310 MultiMathUnaryOperator(O const & o) 00311 : o_(o) 00312 {} 00313 00314 template <class SHAPE> 00315 bool checkShape(SHAPE & s) const 00316 { 00317 return o_.checkShape(s); 00318 } 00319 00320 // 00321 void inc(unsigned int axis) const 00322 { 00323 o_.inc(axis); 00324 } 00325 00326 void reset(unsigned int axis) const 00327 { 00328 o_.reset(axis); 00329 } 00330 00331 template <class POINT> 00332 result_type operator[](POINT const & p) const 00333 { 00334 return f_(o_[p]); 00335 } 00336 00337 result_type operator*() const 00338 { 00339 return f_(*o_); 00340 } 00341 00342 O o_; 00343 F f_; 00344 }; 00345 00346 #define VIGRA_MULTIMATH_UNARY_OPERATOR(NAME, FCT, OPNAME, RESTYPE) \ 00347 namespace detail { \ 00348 struct NAME \ 00349 { \ 00350 template <class T> \ 00351 struct Result \ 00352 { \ 00353 typedef RESTYPE type; \ 00354 }; \ 00355 \ 00356 template <class T> \ 00357 typename Result<T>::type \ 00358 operator()(T const & t) const \ 00359 { \ 00360 return FCT(t); \ 00361 } \ 00362 }; \ 00363 } \ 00364 \ 00365 template <unsigned int N, class T, class C> \ 00366 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<MultiArrayView<N, T, C> >, \ 00367 detail::NAME> > \ 00368 OPNAME(MultiArrayView<N, T, C> const & v) \ 00369 { \ 00370 typedef MultiMathOperand<MultiArrayView<N, T, C> > O; \ 00371 typedef MultiMathUnaryOperator<O, detail::NAME> OP; \ 00372 return MultiMathOperand<OP>(OP(v)); \ 00373 } \ 00374 \ 00375 template <unsigned int N, class T, class A> \ 00376 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<MultiArray<N, T, A> >, \ 00377 detail::NAME> > \ 00378 OPNAME(MultiArray<N, T, A> const & v) \ 00379 { \ 00380 typedef MultiMathOperand<MultiArray<N, T, A> > O; \ 00381 typedef MultiMathUnaryOperator<O, detail::NAME> OP; \ 00382 return MultiMathOperand<OP>(OP(v)); \ 00383 } \ 00384 \ 00385 template <class T> \ 00386 MultiMathOperand<MultiMathUnaryOperator<MultiMathOperand<T>, \ 00387 detail::NAME> > \ 00388 OPNAME(MultiMathOperand<T> const & v) \ 00389 { \ 00390 typedef MultiMathOperand<T> O; \ 00391 typedef MultiMathUnaryOperator<O, detail::NAME> OP; \ 00392 return MultiMathOperand<OP>(OP(v)); \ 00393 } 00394 00395 #define VIGRA_REALPROMOTE typename NumericTraits<T>::RealPromote 00396 00397 #ifndef DOXYGEN // doxygen gets confused by these macros 00398 00399 VIGRA_MULTIMATH_UNARY_OPERATOR(Negate, -, operator-, T) 00400 VIGRA_MULTIMATH_UNARY_OPERATOR(Not, !, operator!, T) 00401 VIGRA_MULTIMATH_UNARY_OPERATOR(BitwiseNot, ~, operator~, T) 00402 00403 using vigra::abs; 00404 VIGRA_MULTIMATH_UNARY_OPERATOR(Abs, vigra::abs, abs, typename NormTraits<T>::NormType) 00405 00406 using vigra::erf; 00407 VIGRA_MULTIMATH_UNARY_OPERATOR(Erf, vigra::erf, erf, VIGRA_REALPROMOTE) 00408 VIGRA_MULTIMATH_UNARY_OPERATOR(Even, vigra::even, even, bool) 00409 VIGRA_MULTIMATH_UNARY_OPERATOR(Odd, vigra::odd, odd, bool) 00410 VIGRA_MULTIMATH_UNARY_OPERATOR(Sign, vigra::sign, sign, T) 00411 VIGRA_MULTIMATH_UNARY_OPERATOR(Signi, vigra::signi, signi, int) 00412 00413 using vigra::round; 00414 VIGRA_MULTIMATH_UNARY_OPERATOR(Round, vigra::round, round, T) 00415 00416 VIGRA_MULTIMATH_UNARY_OPERATOR(Roundi, vigra::roundi, roundi, int) 00417 VIGRA_MULTIMATH_UNARY_OPERATOR(Sqrti, vigra::sqrti, sqrti, T) 00418 VIGRA_MULTIMATH_UNARY_OPERATOR(Sq, vigra::sq, sq, typename NumericTraits<T>::Promote) 00419 VIGRA_MULTIMATH_UNARY_OPERATOR(Norm, vigra::norm, norm, typename NormTraits<T>::NormType) 00420 VIGRA_MULTIMATH_UNARY_OPERATOR(SquaredNorm, vigra::squaredNorm, squaredNorm, 00421 typename NormTraits<T>::SquaredNormType) 00422 VIGRA_MULTIMATH_UNARY_OPERATOR(Sin_pi, vigra::sin_pi, sin_pi, VIGRA_REALPROMOTE) 00423 VIGRA_MULTIMATH_UNARY_OPERATOR(Cos_pi, vigra::cos_pi, cos_pi, VIGRA_REALPROMOTE) 00424 00425 using vigra::gamma; 00426 VIGRA_MULTIMATH_UNARY_OPERATOR(Gamma, vigra::gamma, gamma, VIGRA_REALPROMOTE) 00427 00428 using vigra::loggamma; 00429 VIGRA_MULTIMATH_UNARY_OPERATOR(Loggamma, vigra::loggamma, loggamma, VIGRA_REALPROMOTE) 00430 00431 VIGRA_MULTIMATH_UNARY_OPERATOR(Sqrt, std::sqrt, sqrt, VIGRA_REALPROMOTE) 00432 using vigra::exp; 00433 VIGRA_MULTIMATH_UNARY_OPERATOR(Exp, vigra::exp, exp, VIGRA_REALPROMOTE) 00434 VIGRA_MULTIMATH_UNARY_OPERATOR(Log, std::log, log, VIGRA_REALPROMOTE) 00435 VIGRA_MULTIMATH_UNARY_OPERATOR(Log10, std::log10, log10, VIGRA_REALPROMOTE) 00436 VIGRA_MULTIMATH_UNARY_OPERATOR(Sin, std::sin, sin, VIGRA_REALPROMOTE) 00437 VIGRA_MULTIMATH_UNARY_OPERATOR(Asin, std::asin, asin, VIGRA_REALPROMOTE) 00438 VIGRA_MULTIMATH_UNARY_OPERATOR(Cos, std::cos, cos, VIGRA_REALPROMOTE) 00439 VIGRA_MULTIMATH_UNARY_OPERATOR(Acos, std::acos, acos, VIGRA_REALPROMOTE) 00440 VIGRA_MULTIMATH_UNARY_OPERATOR(Tan, std::tan, tan, VIGRA_REALPROMOTE) 00441 VIGRA_MULTIMATH_UNARY_OPERATOR(Atan, std::atan, atan, VIGRA_REALPROMOTE) 00442 VIGRA_MULTIMATH_UNARY_OPERATOR(Floor, std::floor, floor, VIGRA_REALPROMOTE) 00443 VIGRA_MULTIMATH_UNARY_OPERATOR(Ceil, std::ceil, ceil, VIGRA_REALPROMOTE) 00444 00445 VIGRA_MULTIMATH_UNARY_OPERATOR(Conj, conj, conj, T) 00446 VIGRA_MULTIMATH_UNARY_OPERATOR(Real, real, real, typename T::value_type) 00447 VIGRA_MULTIMATH_UNARY_OPERATOR(Imag, imag, imag, typename T::value_type) 00448 VIGRA_MULTIMATH_UNARY_OPERATOR(Arg, arg, arg, typename T::value_type) 00449 00450 #endif //DOXYGEN 00451 00452 #undef VIGRA_REALPROMOTE 00453 #undef VIGRA_MULTIMATH_UNARY_OPERATOR 00454 00455 template <class O1, class O2, class F> 00456 struct MultiMathBinaryOperator 00457 { 00458 typedef typename F::template Result<typename O1::result_type, 00459 typename O2::result_type>::type result_type; 00460 00461 static const int ndim = O1::ndim > O2::ndim ? O1::ndim : O2::ndim; 00462 00463 MultiMathBinaryOperator(O1 const & o1, O2 const & o2) 00464 : o1_(o1), 00465 o2_(o2) 00466 {} 00467 00468 template <class SHAPE> 00469 bool checkShape(SHAPE & s) const 00470 { 00471 return o1_.checkShape(s) && o2_.checkShape(s); 00472 } 00473 00474 template <class POINT> 00475 result_type operator[](POINT const & p) const 00476 { 00477 return f_(o1_[p], o2_[p]); 00478 } 00479 00480 void inc(unsigned int axis) const 00481 { 00482 o1_.inc(axis); 00483 o2_.inc(axis); 00484 } 00485 00486 void reset(unsigned int axis) const 00487 { 00488 o1_.reset(axis); 00489 o2_.reset(axis); 00490 } 00491 00492 result_type operator*() const 00493 { 00494 return f_(*o1_, *o2_); 00495 } 00496 00497 O1 o1_; 00498 O2 o2_; 00499 F f_; 00500 }; 00501 00502 00503 // In the sequel, the nested type 'MultiMathOperand<T>::AllowOverload' 00504 // ensures that template functions only participate in overload 00505 // resolution when this type is defined, i.e. when T is a number 00506 // or array type. It thus prevents 'ambiguous overload' errors. 00507 // 00508 #define VIGRA_MULTIMATH_BINARY_OPERATOR(NAME, FCT, OPNAME, SEP, RESTYPE) \ 00509 \ 00510 namespace detail { \ 00511 struct NAME \ 00512 { \ 00513 template <class T1, class T2> \ 00514 struct Result \ 00515 { \ 00516 typedef RESTYPE type; \ 00517 }; \ 00518 \ 00519 template <class T1, class T2> \ 00520 typename Result<T1, T2>::type \ 00521 operator()(T1 const & t1, T2 const & t2) const \ 00522 { \ 00523 return FCT(t1 SEP t2); \ 00524 } \ 00525 }; \ 00526 } \ 00527 \ 00528 template <unsigned int N, class T1, class A1, class T2, class A2> \ 00529 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1> >, \ 00530 MultiMathOperand<MultiArrayView<N, T2> >, \ 00531 detail::NAME> > \ 00532 OPNAME(MultiArray<N, T1, A1> const & v1, MultiArray<N, T2, A2> const & v2) \ 00533 { \ 00534 typedef MultiMathOperand<MultiArrayView<N, T1> > O1; \ 00535 typedef MultiMathOperand<MultiArrayView<N, T2> > O2; \ 00536 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00537 return MultiMathOperand<OP>(OP((MultiArrayView<N, T1> const &)v1, (MultiArrayView<N, T2> const &)v2)); \ 00538 } \ 00539 \ 00540 template <unsigned int N, class T1, class C1, class T2, class C2> \ 00541 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \ 00542 MultiMathOperand<MultiArrayView<N, T2, C2> >, \ 00543 detail::NAME> > \ 00544 OPNAME(MultiArrayView<N, T1, C1> const & v1, MultiArrayView<N, T2, C2> const & v2) \ 00545 { \ 00546 typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \ 00547 typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \ 00548 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00549 return MultiMathOperand<OP>(OP(v1, v2)); \ 00550 } \ 00551 \ 00552 template <unsigned int N, class T1, class T2, class C2> \ 00553 MultiMathOperand<MultiMathBinaryOperator<typename MultiMathOperand<T1>::AllowOverload, \ 00554 MultiMathOperand<MultiArrayView<N, T2, C2> >, \ 00555 detail::NAME> > \ 00556 OPNAME(T1 const & v1, MultiArrayView<N, T2, C2> const & v2) \ 00557 { \ 00558 typedef MultiMathOperand<T1> O1; \ 00559 typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \ 00560 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00561 return MultiMathOperand<OP>(OP(v1, v2)); \ 00562 } \ 00563 \ 00564 template <unsigned int N, class T1, class C1, class T2> \ 00565 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \ 00566 typename MultiMathOperand<T2>::AllowOverload, \ 00567 detail::NAME> > \ 00568 OPNAME(MultiArrayView<N, T1, C1> const & v1, T2 const & v2) \ 00569 { \ 00570 typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \ 00571 typedef MultiMathOperand<T2> O2; \ 00572 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00573 return MultiMathOperand<OP>(OP(v1, v2)); \ 00574 } \ 00575 \ 00576 template <unsigned int N, class T1, class T2, class C2> \ 00577 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \ 00578 MultiMathOperand<MultiArrayView<N, T2, C2> >, \ 00579 detail::NAME> > \ 00580 OPNAME(MultiMathOperand<T1> const & v1, MultiArrayView<N, T2, C2> const & v2) \ 00581 { \ 00582 typedef MultiMathOperand<T1> O1; \ 00583 typedef MultiMathOperand<MultiArrayView<N, T2, C2> > O2; \ 00584 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00585 return MultiMathOperand<OP>(OP(v1, v2)); \ 00586 } \ 00587 \ 00588 template <unsigned int N, class T1, class C1, class T2> \ 00589 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<MultiArrayView<N, T1, C1> >, \ 00590 MultiMathOperand<T2>, \ 00591 detail::NAME> > \ 00592 OPNAME(MultiArrayView<N, T1, C1> const & v1, MultiMathOperand<T2> const & v2) \ 00593 { \ 00594 typedef MultiMathOperand<MultiArrayView<N, T1, C1> > O1; \ 00595 typedef MultiMathOperand<T2> O2; \ 00596 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00597 return MultiMathOperand<OP>(OP(v1, v2)); \ 00598 } \ 00599 \ 00600 template <class T1, class T2> \ 00601 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \ 00602 MultiMathOperand<T2>, \ 00603 detail::NAME> > \ 00604 OPNAME(MultiMathOperand<T1> const & v1, MultiMathOperand<T2> const & v2) \ 00605 { \ 00606 typedef MultiMathOperand<T1> O1; \ 00607 typedef MultiMathOperand<T2> O2; \ 00608 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00609 return MultiMathOperand<OP>(OP(v1, v2)); \ 00610 } \ 00611 \ 00612 template <class T1, class T2> \ 00613 MultiMathOperand<MultiMathBinaryOperator<typename MultiMathOperand<T1>::AllowOverload, \ 00614 MultiMathOperand<T2>, \ 00615 detail::NAME> > \ 00616 OPNAME(T1 const & v1, MultiMathOperand<T2> const & v2) \ 00617 { \ 00618 typedef MultiMathOperand<T1> O1; \ 00619 typedef MultiMathOperand<T2> O2; \ 00620 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00621 return MultiMathOperand<OP>(OP(v1, v2)); \ 00622 } \ 00623 \ 00624 template <class T1, class T2> \ 00625 MultiMathOperand<MultiMathBinaryOperator<MultiMathOperand<T1>, \ 00626 typename MultiMathOperand<T2>::AllowOverload, \ 00627 detail::NAME> > \ 00628 OPNAME(MultiMathOperand<T1> const & v1, T2 const & v2) \ 00629 { \ 00630 typedef MultiMathOperand<T1> O1; \ 00631 typedef MultiMathOperand<T2> O2; \ 00632 typedef MultiMathBinaryOperator<O1, O2, detail::NAME> OP; \ 00633 return MultiMathOperand<OP>(OP(v1, v2)); \ 00634 } 00635 00636 #define VIGRA_NOTHING 00637 #define VIGRA_COMMA , 00638 #define VIGRA_PROMOTE typename PromoteTraits<T1, T2>::Promote 00639 #define VIGRA_REALPROMOTE typename PromoteTraits<typename NumericTraits<T1>::RealPromote, \ 00640 typename NumericTraits<T2>::RealPromote>::Promote 00641 00642 VIGRA_MULTIMATH_BINARY_OPERATOR(Plus, VIGRA_NOTHING, operator+, +, VIGRA_PROMOTE) 00643 VIGRA_MULTIMATH_BINARY_OPERATOR(Minus, VIGRA_NOTHING, operator-, -, VIGRA_PROMOTE) 00644 VIGRA_MULTIMATH_BINARY_OPERATOR(Multiplies, VIGRA_NOTHING, operator*, *, VIGRA_PROMOTE) 00645 VIGRA_MULTIMATH_BINARY_OPERATOR(Divides, VIGRA_NOTHING, operator/, /, VIGRA_PROMOTE) 00646 VIGRA_MULTIMATH_BINARY_OPERATOR(Modulo, VIGRA_NOTHING, operator%, %, VIGRA_PROMOTE) 00647 VIGRA_MULTIMATH_BINARY_OPERATOR(And, VIGRA_NOTHING, operator&&, &&, VIGRA_PROMOTE) 00648 VIGRA_MULTIMATH_BINARY_OPERATOR(Or, VIGRA_NOTHING, operator||, ||, VIGRA_PROMOTE) 00649 VIGRA_MULTIMATH_BINARY_OPERATOR(Equal, VIGRA_NOTHING, operator==, ==, bool) 00650 VIGRA_MULTIMATH_BINARY_OPERATOR(NotEqual, VIGRA_NOTHING, operator!=, !=, bool) 00651 VIGRA_MULTIMATH_BINARY_OPERATOR(Less, VIGRA_NOTHING, operator<, <, bool) 00652 VIGRA_MULTIMATH_BINARY_OPERATOR(LessEqual, VIGRA_NOTHING, operator<=, <=, bool) 00653 VIGRA_MULTIMATH_BINARY_OPERATOR(Greater, VIGRA_NOTHING, operator>, >, bool) 00654 VIGRA_MULTIMATH_BINARY_OPERATOR(GreaterEqual, VIGRA_NOTHING, operator>=, >=, bool) 00655 VIGRA_MULTIMATH_BINARY_OPERATOR(Leftshift, VIGRA_NOTHING, operator<<, <<, VIGRA_PROMOTE) 00656 VIGRA_MULTIMATH_BINARY_OPERATOR(Rightshift, VIGRA_NOTHING, operator>>, >>, VIGRA_PROMOTE) 00657 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseAnd, VIGRA_NOTHING, operator&, &, VIGRA_PROMOTE) 00658 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseOr, VIGRA_NOTHING, operator|, |, VIGRA_PROMOTE) 00659 VIGRA_MULTIMATH_BINARY_OPERATOR(BitwiseXor, VIGRA_NOTHING, operator^, ^, VIGRA_PROMOTE) 00660 00661 VIGRA_MULTIMATH_BINARY_OPERATOR(Atan2, std::atan2, atan2, VIGRA_COMMA, VIGRA_REALPROMOTE) 00662 VIGRA_MULTIMATH_BINARY_OPERATOR(Pow, std::pow, pow, VIGRA_COMMA, VIGRA_REALPROMOTE) 00663 VIGRA_MULTIMATH_BINARY_OPERATOR(Fmod, std::fmod, fmod, VIGRA_COMMA, VIGRA_REALPROMOTE) 00664 VIGRA_MULTIMATH_BINARY_OPERATOR(Min, std::min, min, VIGRA_COMMA, VIGRA_PROMOTE) 00665 VIGRA_MULTIMATH_BINARY_OPERATOR(Max, std::max, max, VIGRA_COMMA, VIGRA_PROMOTE) 00666 VIGRA_MULTIMATH_BINARY_OPERATOR(Minimum, std::min, minimum, VIGRA_COMMA, VIGRA_PROMOTE) 00667 VIGRA_MULTIMATH_BINARY_OPERATOR(Maximum, std::max, maximum, VIGRA_COMMA, VIGRA_PROMOTE) 00668 00669 #undef VIGRA_NOTHING 00670 #undef VIGRA_COMMA 00671 #undef VIGRA_PROMOTE 00672 #undef VIGRA_REALPROMOTE 00673 #undef VIGRA_MULTIMATH_BINARY_OPERATOR 00674 00675 namespace detail { 00676 00677 // We pass 'strideOrder' to the recursion in order to make sure 00678 // that the inner loop iterates over the output's major axis. 00679 // Of course, this does not help when the RHS arrays are ordered 00680 // differently -- maybe it is better to find the most common order 00681 // among all arguments (both RHS and LHS)? 00682 // 00683 template <unsigned int N, class Assign> 00684 struct MultiMathExec 00685 { 00686 enum { LEVEL = N-1 }; 00687 00688 template <class T, class Shape, class Expression> 00689 static void exec(T * data, Shape const & shape, Shape const & strides, 00690 Shape const & strideOrder, Expression const & e) 00691 { 00692 MultiArrayIndex axis = strideOrder[LEVEL]; 00693 for(MultiArrayIndex k=0; k<shape[axis]; ++k, data += strides[axis], e.inc(axis)) 00694 { 00695 MultiMathExec<N-1, Assign>::exec(data, shape, strides, strideOrder, e); 00696 } 00697 e.reset(axis); 00698 data -= shape[axis]*strides[axis]; 00699 } 00700 }; 00701 00702 template <class Assign> 00703 struct MultiMathExec<1, Assign> 00704 { 00705 enum { LEVEL = 0 }; 00706 00707 template <class T, class Shape, class Expression> 00708 static void exec(T * data, Shape const & shape, Shape const & strides, 00709 Shape const & strideOrder, Expression const & e) 00710 { 00711 MultiArrayIndex axis = strideOrder[LEVEL]; 00712 for(MultiArrayIndex k=0; k<shape[axis]; ++k, data += strides[axis], e.inc(axis)) 00713 { 00714 Assign::assign(data, e); 00715 } 00716 e.reset(axis); 00717 data -= shape[axis]*strides[axis]; 00718 } 00719 }; 00720 00721 #define VIGRA_MULTIMATH_ASSIGN(NAME, OP) \ 00722 struct MultiMath##NAME \ 00723 { \ 00724 template <class T, class Expression> \ 00725 static void assign(T * data, Expression const & e) \ 00726 { \ 00727 *data OP vigra::detail::RequiresExplicitCast<T>::cast(*e); \ 00728 } \ 00729 }; \ 00730 \ 00731 template <unsigned int N, class T, class C, class Expression> \ 00732 void NAME(MultiArrayView<N, T, C> a, MultiMathOperand<Expression> const & e) \ 00733 { \ 00734 typename MultiArrayShape<N>::type shape(a.shape()); \ 00735 \ 00736 vigra_precondition(e.checkShape(shape), \ 00737 "multi_math: shape mismatch in expression."); \ 00738 \ 00739 MultiMathExec<N, MultiMath##NAME>::exec(a.data(), a.shape(), a.stride(), \ 00740 a.strideOrdering(), e); \ 00741 } \ 00742 \ 00743 template <unsigned int N, class T, class A, class Expression> \ 00744 void NAME##OrResize(MultiArray<N, T, A> & a, MultiMathOperand<Expression> const & e) \ 00745 { \ 00746 typename MultiArrayShape<N>::type shape(a.shape()); \ 00747 \ 00748 vigra_precondition(e.checkShape(shape), \ 00749 "multi_math: shape mismatch in expression."); \ 00750 \ 00751 if(a.size() == 0) \ 00752 a.reshape(shape); \ 00753 \ 00754 MultiMathExec<N, MultiMath##NAME>::exec(a.data(), a.shape(), a.stride(), \ 00755 a.strideOrdering(), e); \ 00756 } 00757 00758 VIGRA_MULTIMATH_ASSIGN(assign, =) 00759 VIGRA_MULTIMATH_ASSIGN(plusAssign, +=) 00760 VIGRA_MULTIMATH_ASSIGN(minusAssign, -=) 00761 VIGRA_MULTIMATH_ASSIGN(multiplyAssign, *=) 00762 VIGRA_MULTIMATH_ASSIGN(divideAssign, /=) 00763 00764 #undef VIGRA_MULTIMATH_ASSIGN 00765 00766 template <unsigned int N, class Assign> 00767 struct MultiMathReduce 00768 { 00769 enum { LEVEL = N-1 }; 00770 00771 template <class T, class Shape, class Expression> 00772 static void exec(T & t, Shape const & shape, Expression const & e) 00773 { 00774 for(MultiArrayIndex k=0; k<shape[LEVEL]; ++k, e.inc(LEVEL)) 00775 { 00776 MultiMathReduce<N-1, Assign>::exec(t, shape, e); 00777 } 00778 e.reset(LEVEL); 00779 } 00780 }; 00781 00782 template <class Assign> 00783 struct MultiMathReduce<1, Assign> 00784 { 00785 enum { LEVEL = 0 }; 00786 00787 template <class T, class Shape, class Expression> 00788 static void exec(T & t, Shape const & shape, Expression const & e) 00789 { 00790 for(MultiArrayIndex k=0; k<shape[0]; ++k, e.inc(0)) 00791 { 00792 Assign::assign(&t, e); 00793 } 00794 e.reset(0); 00795 } 00796 }; 00797 00798 struct MultiMathReduceAll 00799 { 00800 template <class T, class Expression> 00801 static void assign(T * data, Expression const & e) 00802 { 00803 *data = *data && (*e != NumericTraits<typename Expression::result_type>::zero()); 00804 } 00805 }; 00806 00807 struct MultiMathReduceAny 00808 { 00809 template <class T, class Expression> 00810 static void assign(T * data, Expression const & e) 00811 { 00812 *data = *data || (*e != NumericTraits<typename Expression::result_type>::zero()); 00813 } 00814 }; 00815 00816 00817 } // namespace detail 00818 00819 template <class U, class T> 00820 U 00821 sum(MultiMathOperand<T> const & v, U res = NumericTraits<U>::zero()) 00822 { 00823 static const int ndim = MultiMathOperand<T>::ndim; 00824 typename MultiArrayShape<ndim>::type shape; 00825 v.checkShape(shape); 00826 detail::MultiMathReduce<ndim, detail::MultiMathplusAssign>::exec(res, shape, v); 00827 return res; 00828 } 00829 00830 template <class U, unsigned int N, class T, class S> 00831 U 00832 sum(MultiArrayView<N, T, S> const & v, U res = NumericTraits<U>::zero()) 00833 { 00834 return v.sum<U>() + res; 00835 } 00836 00837 template <class U, class T> 00838 U 00839 product(MultiMathOperand<T> const & v, U res = NumericTraits<U>::one()) 00840 { 00841 static const int ndim = MultiMathOperand<T>::ndim; 00842 typename MultiArrayShape<ndim>::type shape; 00843 v.checkShape(shape); 00844 detail::MultiMathReduce<ndim, detail::MultiMathmultiplyAssign>::exec(res, shape, v); 00845 return res; 00846 } 00847 00848 template <class U, unsigned int N, class T, class S> 00849 U 00850 product(MultiArrayView<N, T, S> const & v, U res = NumericTraits<U>::one()) 00851 { 00852 return v.product<U>() * res; 00853 } 00854 00855 template <class T> 00856 bool 00857 all(MultiMathOperand<T> const & v) 00858 { 00859 static const int ndim = MultiMathOperand<T>::ndim; 00860 typename MultiArrayShape<ndim>::type shape; 00861 v.checkShape(shape); 00862 bool res = true; 00863 detail::MultiMathReduce<ndim, detail::MultiMathReduceAll>::exec(res, shape, v); 00864 return res; 00865 } 00866 00867 template <class T> 00868 bool 00869 any(MultiMathOperand<T> const & v) 00870 { 00871 static const int ndim = MultiMathOperand<T>::ndim; 00872 typename MultiArrayShape<ndim>::type shape; 00873 v.checkShape(shape); 00874 bool res = false; 00875 detail::MultiMathReduce<ndim, detail::MultiMathReduceAny>::exec(res, shape, v); 00876 return res; 00877 } 00878 00879 00880 }} // namespace vigra::multi_math 00881 00882 #endif // VIGRA_MULTI_MATH_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|