vigra - vigra/multi_math.hxx Source File

[ 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>+ - * / ! ~ % && || == != &lt; &lt;= &gt; &gt;= &lt;&lt; &gt;&gt; & | ^ = += -= *= /=</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, vigra::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.template 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.template 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) Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany | html generated using doxygen and Python vigra 1.9.0 (Tue Nov 6 2012)