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

vigra/multi_tensorutilities.hxx VIGRA

00001 /************************************************************************/
00002 /*                                                                      */
00003 /*               Copyright 2009-2010 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_TENSORUTILITIES_HXX
00037 #define VIGRA_MULTI_TENSORUTILITIES_HXX
00038 
00039 #include <cmath>
00040 #include "utilities.hxx"
00041 #include "mathutil.hxx"
00042 #include "metaprogramming.hxx"
00043 #include "multi_pointoperators.hxx"
00044 
00045 namespace vigra {
00046 
00047 namespace detail {
00048 
00049 template <int N, class ArgumentVector, class ResultVector>
00050 class OuterProductFunctor
00051 {
00052 public:
00053     typedef ArgumentVector argument_type;
00054     typedef ResultVector result_type;
00055     typedef typename ArgumentVector::value_type ValueType;
00056     
00057     result_type operator()(argument_type const & in) const
00058     {
00059         result_type res;
00060         for(int b=0, i=0; i<N; ++i)
00061         {
00062             for(int j=i; j<N; ++j, ++b)
00063             {
00064                 res[b] = detail::RequiresExplicitCast<ValueType>::cast(in[i]*in[j]);
00065             }
00066         }
00067         return res;
00068     }
00069 };
00070 
00071 template <int N, class ArgumentVector>
00072 class TensorTraceFunctor
00073 {
00074 public:
00075 
00076     typedef ArgumentVector argument_type;
00077     typedef typename ArgumentVector::value_type result_type;
00078     
00079     result_type exec(argument_type const & v, MetaInt<1>) const
00080     {
00081         return v[0];
00082     }
00083     
00084     result_type exec(argument_type const & v, MetaInt<2>) const
00085     {
00086         return v[0] + v[2];
00087     }
00088     
00089     result_type exec(argument_type const & v, MetaInt<3>) const
00090     {
00091         return v[0] + v[3] + v[5];
00092     }
00093     
00094     template <int N2>
00095     void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
00096     {
00097         vigra_fail("tensorTraceMultiArray(): Sorry, can only handle dimensions up to 3.");
00098     }
00099 
00100     result_type operator()( const argument_type & a ) const
00101     {
00102         return exec(a, MetaInt<N>());
00103     }
00104 };
00105 
00106 template <int N, class ArgumentVector, class ResultVector>
00107 class EigenvaluesFunctor
00108 {
00109 public:
00110 
00111     typedef ArgumentVector argument_type;
00112     typedef ResultVector result_type;
00113     
00114     void exec(argument_type const & v, result_type & r, MetaInt<1>) const
00115     {
00116         symmetric2x2Eigenvalues(v[0], &r[0]);
00117     }
00118     
00119     void exec(argument_type const & v, result_type & r, MetaInt<2>) const
00120     {
00121         symmetric2x2Eigenvalues(v[0], v[1], v[2], &r[0], &r[1]);
00122     }
00123     
00124     void exec(argument_type const & v, result_type & r, MetaInt<3>) const
00125     {
00126         symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r[0], &r[1], &r[2]);
00127     }
00128     
00129     template <int N2>
00130     void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
00131     {
00132         vigra_fail("tensorEigenvaluesMultiArray(): Sorry, can only handle dimensions up to 3.");
00133     }
00134 
00135     result_type operator()( const argument_type & a ) const
00136     {
00137         result_type res;
00138         exec(a, res, MetaInt<N>());
00139         return res;
00140     }
00141 };
00142 
00143 
00144 template <int N, class ArgumentVector>
00145 class DeterminantFunctor
00146 {
00147 public:
00148 
00149     typedef ArgumentVector argument_type;
00150     typedef typename ArgumentVector::value_type result_type;
00151     
00152     result_type exec(argument_type const & v, MetaInt<1>) const
00153     {
00154         return v[0];
00155     }
00156     
00157     result_type exec(argument_type const & v, MetaInt<2>) const
00158     {
00159         return v[0]*v[2] - sq(v[1]);
00160     }
00161     
00162     result_type exec(argument_type const & v, MetaInt<3>) const
00163     {
00164         result_type r0, r1, r2;
00165         symmetric3x3Eigenvalues(v[0], v[1], v[2], v[3], v[4], v[5], &r0, &r1, &r2);
00166         return r0*r1*r2;
00167     }
00168     
00169     template <int N2>
00170     void exec(argument_type const & v, result_type & r, MetaInt<N2>) const
00171     {
00172         vigra_fail("tensorDeterminantMultiArray(): Sorry, can only handle dimensions up to 3.");
00173     }
00174 
00175     result_type operator()( const argument_type & a ) const
00176     {
00177         return exec(a, MetaInt<N>());
00178     }
00179 };
00180 
00181 } // namespace detail
00182 
00183 
00184 /** \addtogroup MultiPointoperators
00185 */
00186 //@{
00187 
00188 /********************************************************/
00189 /*                                                      */
00190 /*                vectorToTensorMultiArray              */
00191 /*                                                      */
00192 /********************************************************/
00193 
00194 /** \brief Calculate the tensor (outer) product of a N-D vector with itself.
00195 
00196     This function is useful to transform vector arrays into a tensor representation 
00197     that can be used as input to tensor based processing and analysis functions
00198     (e.g. tensor smoothing). When the input array has N dimensions, the input value_type 
00199     must be a vector of length N, whereas the output value_type mus be vectors of length 
00200     N*(N-1)/2 which will represent the upper triangular part of the resulting (symmetric) 
00201     tensor. That is, for 2D arrays the output contains the elements 
00202     <tt>[t11, t12 == t21, t22]</tt> in this order, whereas it contains the elements
00203     <tt>[t11, t12, t13, t22, t23, t33]</tt> for 3D arrays.
00204     
00205     Currently, <tt>N <= 3</tt> is required.
00206     
00207     <b> Declarations:</b>
00208 
00209     pass arguments explicitly:
00210     \code
00211     namespace vigra {
00212         template <class SrcIterator, class SrcShape, class SrcAccessor,
00213                   class DestIterator, class DestAccessor>
00214         void 
00215         vectorToTensorMultiArray(SrcIterator  si, SrcShape const & shape, SrcAccessor src,
00216                                  DestIterator di, DestAccessor dest);
00217     }
00218     \endcode
00219 
00220     use argument objects in conjunction with \ref ArgumentObjectFactories :
00221     \code
00222     namespace vigra {
00223         template <class SrcIterator, class SrcShape, class SrcAccessor,
00224                   class DestIterator, class DestAccessor>
00225         void 
00226         vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00227                                  pair<DestIterator, DestAccessor> d);
00228     }
00229     \endcode
00230 
00231     <b> Usage:</b>
00232 
00233     <b>\#include</b> <vigra/multi_tensorutilities.hxx>
00234 
00235     \code
00236     MultiArray<3, float> vol(shape);
00237     MultiArray<3, TinyVector<float, 3> > gradient(shape);
00238     MultiArray<3, TinyVector<float, 6> > tensor(shape);
00239     
00240     gaussianGradientMultiArray(srcMultiArrayRange(vol), destMultiArray(gradient), 2.0);
00241     vectorToTensorMultiArray(srcMultiArrayRange(gradient), destMultiArray(tensor));
00242     \endcode
00243 
00244 */
00245 doxygen_overloaded_function(template <...> void vectorToTensorMultiArray)
00246 
00247 template <class SrcIterator, class SrcShape, class SrcAccessor,
00248           class DestIterator, class DestAccessor>
00249 void 
00250 vectorToTensorMultiArray(SrcIterator  si, SrcShape const & shape, SrcAccessor src,
00251                          DestIterator di, DestAccessor dest)
00252 {
00253     static const int N = SrcShape::static_size;
00254     static const int M = N*(N+1)/2;
00255     
00256     typedef typename SrcAccessor::value_type  SrcType;
00257     typedef typename DestAccessor::value_type DestType;
00258 
00259     for(int k=0; k<N; ++k)
00260         if(shape[k] <=0)
00261             return;
00262 
00263     vigra_precondition(N == (int)src.size(si),
00264         "vectorToTensorMultiArray(): Wrong number of channels in input array.");
00265     vigra_precondition(M == (int)dest.size(di),
00266         "vectorToTensorMultiArray(): Wrong number of channels in output array.");
00267 
00268     transformMultiArray(si, shape, src, di, dest, 
00269                         detail::OuterProductFunctor<N, SrcType, DestType>());
00270 }
00271 
00272 template <class SrcIterator, class SrcShape, class SrcAccessor,
00273           class DestIterator, class DestAccessor>
00274 inline void 
00275 vectorToTensorMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00276                          pair<DestIterator, DestAccessor> d)
00277 {
00278     vectorToTensorMultiArray(s.first, s.second, s.third, d.first, d.second);
00279 }
00280 
00281 /********************************************************/
00282 /*                                                      */
00283 /*                tensorTraceMultiArray                 */
00284 /*                                                      */
00285 /********************************************************/
00286 
00287 /** \brief Calculate the tensor trace for every element of a N-D tensor array.
00288 
00289     This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2, 
00290     see \ref vectorToTensorMultiArray()) representing the upper triangular part of a 
00291     symmetric tensor into a scalar array holding the tensor trace.
00292     
00293     Currently, <tt>N <= 3</tt> is required.
00294     
00295     <b> Declarations:</b>
00296 
00297     pass arguments explicitly:
00298     \code
00299     namespace vigra {
00300         template <class SrcIterator, class SrcShape, class SrcAccessor,
00301                   class DestIterator, class DestAccessor>
00302         void 
00303         tensorTraceMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00304                               DestIterator di, DestAccessor dest);
00305     }
00306     \endcode
00307 
00308 
00309     use argument objects in conjunction with \ref ArgumentObjectFactories :
00310     \code
00311     namespace vigra {
00312         template <class SrcIterator, class SrcShape, class SrcAccessor,
00313                   class DestIterator, class DestAccessor>
00314         void 
00315         tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00316                               pair<DestIterator, DestAccessor> d);
00317     }
00318     \endcode
00319 
00320     <b> Usage:</b>
00321 
00322     <b>\#include</b> <vigra/multi_tensorutilities.hxx>
00323 
00324     \code
00325     MultiArray<3, float> vol(shape);
00326     MultiArray<3, TinyVector<float, 6> > hessian(shape);
00327     MultiArray<3, float> trace(shape);
00328     
00329     hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
00330     tensorTraceMultiArray(srcMultiArrayRange(hessian), destMultiArray(trace));
00331     \endcode
00332 
00333 */
00334 doxygen_overloaded_function(template <...> void tensorTraceMultiArray)
00335 
00336 template <class SrcIterator, class SrcShape, class SrcAccessor,
00337           class DestIterator, class DestAccessor>
00338 void 
00339 tensorTraceMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00340                       DestIterator di, DestAccessor dest)
00341 {
00342     static const int N = SrcShape::static_size;
00343     typedef typename SrcAccessor::value_type  SrcType;
00344 
00345     transformMultiArray(si, shape, src, di, dest, 
00346                         detail::TensorTraceFunctor<N, SrcType>());
00347 }
00348 
00349 template <class SrcIterator, class SrcShape, class SrcAccessor,
00350           class DestIterator, class DestAccessor>
00351 inline void 
00352 tensorTraceMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00353                       pair<DestIterator, DestAccessor> d)
00354 {
00355     tensorTraceMultiArray(s.first, s.second, s.third, d.first, d.second);
00356 }
00357 
00358 
00359 /********************************************************/
00360 /*                                                      */
00361 /*             tensorEigenvaluesMultiArray              */
00362 /*                                                      */
00363 /********************************************************/
00364 
00365 /** \brief Calculate the tensor eigenvalues for every element of a N-D tensor array.
00366 
00367     This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2, 
00368     see \ref vectorToTensorMultiArray()) representing the upper triangular part of a 
00369     symmetric tensor into a vector-valued array holding the tensor eigenvalues (thus,
00370     the destination value_type must be vectors of length N).
00371     
00372     Currently, <tt>N <= 3</tt> is required.
00373     
00374     <b> Declarations:</b>
00375 
00376     pass arguments explicitly:
00377     \code
00378     namespace vigra {
00379         template <class SrcIterator, class SrcShape, class SrcAccessor,
00380                   class DestIterator, class DestAccessor>
00381         void 
00382         tensorEigenvaluesMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00383                                     DestIterator di, DestAccessor dest);
00384     }
00385     \endcode
00386 
00387 
00388     use argument objects in conjunction with \ref ArgumentObjectFactories :
00389     \code
00390     namespace vigra {
00391         template <class SrcIterator, class SrcShape, class SrcAccessor,
00392                   class DestIterator, class DestAccessor>
00393         void 
00394         tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00395                                     pair<DestIterator, DestAccessor> d);
00396     }
00397     \endcode
00398 
00399     <b> Usage:</b>
00400 
00401     <b>\#include</b> <vigra/multi_tensorutilities.hxx>
00402 
00403     \code
00404     MultiArray<3, float> vol(shape);
00405     MultiArray<3, TinyVector<float, 6> > hessian(shape);
00406     MultiArray<3, TinyVector<float, 3> > eigenvalues(shape);
00407     
00408     hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
00409     tensorEigenvaluesMultiArray(srcMultiArrayRange(hessian), destMultiArray(eigenvalues));
00410     \endcode
00411 
00412 */
00413 doxygen_overloaded_function(template <...> void tensorEigenvaluesMultiArray)
00414 
00415 template <class SrcIterator, class SrcShape, class SrcAccessor,
00416           class DestIterator, class DestAccessor>
00417 void 
00418 tensorEigenvaluesMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00419                             DestIterator di, DestAccessor dest)
00420 {
00421     static const int N = SrcShape::static_size;
00422     static const int M = N*(N+1)/2;
00423     
00424     typedef typename SrcAccessor::value_type  SrcType;
00425     typedef typename DestAccessor::value_type DestType;
00426 
00427     for(int k=0; k<N; ++k)
00428         if(shape[k] <=0)
00429             return;
00430 
00431     vigra_precondition(M == (int)src.size(si),
00432         "tensorEigenvaluesMultiArray(): Wrong number of channels in input array.");
00433     vigra_precondition(N == (int)dest.size(di),
00434         "tensorEigenvaluesMultiArray(): Wrong number of channels in output array.");
00435 
00436     transformMultiArray(si, shape, src, di, dest, 
00437                         detail::EigenvaluesFunctor<N, SrcType, DestType>());
00438 }
00439 
00440 template <class SrcIterator, class SrcShape, class SrcAccessor,
00441           class DestIterator, class DestAccessor>
00442 inline void 
00443 tensorEigenvaluesMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00444                             pair<DestIterator, DestAccessor> d)
00445 {
00446     tensorEigenvaluesMultiArray(s.first, s.second, s.third, d.first, d.second);
00447 }
00448 
00449 /********************************************************/
00450 /*                                                      */
00451 /*             tensorDeterminantMultiArray              */
00452 /*                                                      */
00453 /********************************************************/
00454 
00455 /** \brief Calculate the tensor determinant for every element of a ND tensor array.
00456 
00457     This function turns a N-D tensor (whose value_type is a vector of length N*(N+1)/2, 
00458     see \ref vectorToTensorMultiArray()) representing the upper triangular part of a 
00459     symmetric tensor into the a scalar array holding the tensor determinant.
00460     
00461     Currently, <tt>N <= 3</tt> is required.
00462     
00463     <b> Declarations:</b>
00464 
00465     pass arguments explicitly:
00466     \code
00467     namespace vigra {
00468         template <class SrcIterator, class SrcShape, class SrcAccessor,
00469                   class DestIterator, class DestAccessor>
00470         void 
00471         tensorDeterminantMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00472                                     DestIterator di, DestAccessor dest);
00473     }
00474     \endcode
00475 
00476 
00477     use argument objects in conjunction with \ref ArgumentObjectFactories :
00478     \code
00479     namespace vigra {
00480         template <class SrcIterator, class SrcShape, class SrcAccessor,
00481                   class DestIterator, class DestAccessor>
00482         void 
00483         tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00484                                     pair<DestIterator, DestAccessor> d);
00485     }
00486     \endcode
00487 
00488     <b> Usage:</b>
00489 
00490     <b>\#include</b> <vigra/multi_tensorutilities.hxx>
00491 
00492     \code
00493     MultiArray<3, float> vol(shape);
00494     MultiArray<3, TinyVector<float, 6> > hessian(shape);
00495     MultiArray<3, float> determinant(shape);
00496     
00497     hessianOfGaussianMultiArray(srcMultiArrayRange(vol), destMultiArray(hessian), 2.0);
00498     tensorDeterminantMultiArray(srcMultiArrayRange(hessian), destMultiArray(determinant));
00499     \endcode
00500 
00501 */
00502 doxygen_overloaded_function(template <...> void tensorDeterminantMultiArray)
00503 
00504 template <class SrcIterator, class SrcShape, class SrcAccessor,
00505           class DestIterator, class DestAccessor>
00506 void 
00507 tensorDeterminantMultiArray(SrcIterator si,  SrcShape const & shape, SrcAccessor src,
00508                             DestIterator di, DestAccessor dest)
00509 {
00510     typedef typename SrcAccessor::value_type  SrcType;
00511 
00512     static const int N = SrcShape::static_size;
00513     static const int M = N*(N+1)/2;
00514     
00515     for(int k=0; k<N; ++k)
00516         if(shape[k] <=0)
00517             return;
00518 
00519     vigra_precondition(M == (int)src.size(si),
00520         "tensorDeterminantMultiArray(): Wrong number of channels in output array.");
00521 
00522     transformMultiArray(si, shape, src, di, dest, 
00523                         detail::DeterminantFunctor<N, SrcType>());
00524 }
00525 
00526 template <class SrcIterator, class SrcShape, class SrcAccessor,
00527           class DestIterator, class DestAccessor>
00528 inline void 
00529 tensorDeterminantMultiArray(triple<SrcIterator, SrcShape, SrcAccessor> s,
00530                             pair<DestIterator, DestAccessor> d)
00531 {
00532     tensorDeterminantMultiArray(s.first, s.second, s.third, d.first, d.second);
00533 }
00534 
00535 //@}
00536 
00537 } // namespace vigra
00538 
00539 #endif /* VIGRA_MULTI_TENSORUTILITIES_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)