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

vigra/multi_math.hxx VIGRA

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, 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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.8.0 (20 Sep 2011)