QuOc

 

ChanVese.h

Go to the documentation of this file.
00001 /* Copyright (c) 2001-2010 AG Rumpf, INS, Universitaet Bonn                      */
00002 /*                                                                               */
00003 /* The contents of this file are subject to the terms of the Common Development  */
00004 /* and Distribution License Version 1.0 (the "License"); you may not use         */
00005 /* this file except in compliance with the License. You may obtain a copy of     */
00006 /* the License at http://www.sun.com/cddl/                                       */
00007 /*                                                                               */
00008 /* Software distributed under the License is distributed on an "AS IS" basis,    */
00009 /* WITHOUT WARRANTY OF ANY KIND, either express or implied.                      */
00010 /*                                                                               */
00011 #ifndef __CHANVESE_H
00012 #define __CHANVESE_H
00013 
00014 #include <FEOpInterface.h>
00015 #include <RudinOsherFatemi.h>
00016 #include <mcm.h>
00017 
00018 namespace aol {
00019 
00020 template <typename ConfiguratorType, typename HeavisideFunctionType>
00021 class HeavisideFunctionMassOp : public aol::FELinScalarWeightedMassInterface<ConfiguratorType, HeavisideFunctionMassOp<ConfiguratorType, HeavisideFunctionType> > {
00022 public:
00023   typedef typename ConfiguratorType::RealType RealType;
00024 protected:
00025   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00026   const HeavisideFunctionType &_heavisideFunction;
00027 public:
00028   HeavisideFunctionMassOp( const typename ConfiguratorType::InitType &Initializer, 
00029                            const aol::Vector<RealType> &LevelsetFunction,
00030                            const HeavisideFunctionType &HeavisideFunction,
00031                            aol::OperatorType OpType = aol::ONTHEFLY ) 
00032     : aol::FELinScalarWeightedMassInterface<ConfiguratorType, HeavisideFunctionMassOp<ConfiguratorType, HeavisideFunctionType> >( Initializer, OpType ),
00033       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00034       _heavisideFunction( HeavisideFunction )
00035   {
00036   }
00037   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& ) const { 
00038     return 1./(_heavisideFunction.evaluateDerivative(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint)));
00039   }
00040 };
00041 
00042 template <typename ConfiguratorType, typename HeavisideFunctionType>
00043 class HeavisideFunctionLumpedMassOp : public aol::LumpedMassOpInterface<ConfiguratorType, HeavisideFunctionLumpedMassOp<ConfiguratorType, HeavisideFunctionType> > {
00044 public:
00045   typedef typename ConfiguratorType::RealType RealType;
00046 protected:
00047   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00048   const HeavisideFunctionType &_heavisideFunction;
00049 public:
00050   HeavisideFunctionLumpedMassOp( const typename ConfiguratorType::InitType &Initializer,
00051                                  const aol::Vector<RealType> &LevelsetFunction,
00052                                  const HeavisideFunctionType &HeavisideFunction,
00053                                  const bool Invert = false )
00054     : aol::LumpedMassOpInterface<ConfiguratorType, HeavisideFunctionLumpedMassOp<ConfiguratorType, HeavisideFunctionType> >( Initializer, Invert ),
00055       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00056       _heavisideFunction( HeavisideFunction )
00057   {
00058   }
00059   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& ) const {
00060     return 1./(_heavisideFunction.evaluateDerivative(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint)));
00061   }
00062 };
00063 
00064 /**
00065  * \author Berkels
00066  * \warning The default value for \a Epsilon differs from the default value used in 
00067  *          VariationOfHeavisideLevelsetLengthEnergy. This will be corrected
00068  *          in the future, but for legacy reasons is not yet done.
00069  */
00070 template <typename ConfiguratorType, typename HeavisideFunctionType, int numberOfLevelsetFunctions = 1>
00071 class HeavisideLevelsetLengthEnergy
00072 : public aol::FENonlinIntegrationVectorInterface<ConfiguratorType,
00073                                                    HeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType, numberOfLevelsetFunctions>,
00074                                                    numberOfLevelsetFunctions > {
00075 public:
00076   typedef typename ConfiguratorType::RealType RealType;
00077 protected:
00078   const HeavisideFunctionType &_heavisideFunction;
00079   const RealType _epsilonSqr;
00080 public:
00081   HeavisideLevelsetLengthEnergy( const typename ConfiguratorType::InitType &Initializer,
00082                                  const HeavisideFunctionType &HeavisideFunction,
00083                                  const RealType Epsilon = 0. )
00084     : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00085                                                HeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType, numberOfLevelsetFunctions>,
00086                                                numberOfLevelsetFunctions > ( Initializer ),
00087       _heavisideFunction ( HeavisideFunction ),
00088       _epsilonSqr ( aol::Sqr( Epsilon ) ) {
00089   }
00090 
00091   RealType evaluateIntegrand( const aol::auto_container<numberOfLevelsetFunctions,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00092                               const typename ConfiguratorType::ElementType &El,
00093                               int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
00094     typename ConfiguratorType::VecType gradPhi;
00095     RealType integrand = 0.;
00096     for( int i = 0; i < numberOfLevelsetFunctions; i++ ){
00097       discrFuncs[i].evaluateGradientAtQuadPoint( El, QuadPoint, gradPhi );
00098       integrand += ( _heavisideFunction.evaluateDerivative(discrFuncs[i].evaluateAtQuadPoint( El, QuadPoint))
00099                      * sqrt ( gradPhi.normSqr() + _epsilonSqr ) );
00100     }
00101     return integrand;
00102   }
00103 };
00104 
00105 /**
00106  * Calulates the variation of HeavisideLevelsetLengthEnergy without the derivative
00107  * of the Heaviside function part. This has to be handled by the metric of the
00108  * gradient flow.
00109  *
00110  * \author Berkels
00111  */
00112 template <typename ConfiguratorType, int numberOfLevelsetFunctions = 1>
00113 class VariationOfHeavisideLevelsetLengthEnergy
00114 : public Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
00115 public:
00116   typedef typename ConfiguratorType::RealType RealType;
00117 protected:
00118   const typename ConfiguratorType::InitType &_grid;
00119   const RealType _epsilon;
00120   const RealType _scale;
00121 public:
00122   /**
00123    * \param Epsilon parameter to regularize the absolute value
00124    * \param Scale scales the whole result MultiVector
00125    */
00126   VariationOfHeavisideLevelsetLengthEnergy( const typename ConfiguratorType::InitType &Initializer,
00127                                             const RealType Epsilon = 0.1,
00128                                             const RealType Scale = 1.)
00129     : _grid( Initializer ),
00130       _epsilon( Epsilon ),
00131       _scale( Scale )
00132   {
00133   }
00134 
00135   void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::MultiVector<RealType> &MDest ) const {
00136     qc::MCMStiffOp<ConfiguratorType> mcmStiffOp( _grid, aol::ONTHEFLY, _epsilon, _scale );
00137     for( int i = 0; i < numberOfLevelsetFunctions; i++ ){
00138       mcmStiffOp.setImageReference( MArg[i] );
00139       mcmStiffOp.applyAdd( MArg[i], MDest[i] );
00140     }
00141   }
00142 };
00143 
00144 /**
00145  * Computes \f$ \int H''(\phi(x))) \|\nabla\phi(x)\|\varphi_j dx \f$, where arg[0]=phi and HeavisideFunction=H.
00146  *
00147  * \author Berkels
00148  */
00149 template <typename ConfiguratorType, typename HeavisideFunctionType>
00150 class FirstPartOfFullVariationOfHeavisideLevelsetLengthEnergy
00151   : public aol::FENonlinOpInterface< ConfiguratorType, FirstPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType> > {
00152   typedef typename ConfiguratorType::RealType RealType;
00153   const HeavisideFunctionType &_heavisideFunction;
00154   const RealType _epsilonSqr;
00155 public:
00156   FirstPartOfFullVariationOfHeavisideLevelsetLengthEnergy( const typename ConfiguratorType::InitType &Initializer,
00157                                                            const HeavisideFunctionType &HeavisideFunction,
00158                                                            const RealType Epsilon )
00159     : aol::FENonlinOpInterface< ConfiguratorType,
00160                                  FirstPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType> > ( Initializer ),
00161       _heavisideFunction( HeavisideFunction ),
00162       _epsilonSqr( aol::Sqr(Epsilon) ) {
00163   }
00164 
00165   void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc, 
00166                         const typename ConfiguratorType::ElementType &El, 
00167                         int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/,
00168                         typename ConfiguratorType::RealType &NL ) const {
00169     const RealType hPrimePrime = _heavisideFunction.evaluateSecondDerivative(DiscFunc.evaluateAtQuadPoint( El, QuadPoint));
00170 
00171     typename ConfiguratorType::VecType gradPhi;
00172     DiscFunc.evaluateGradientAtQuadPoint( El, QuadPoint, gradPhi );
00173 
00174     NL = hPrimePrime * sqrt( gradPhi.normSqr() + _epsilonSqr );
00175   }
00176 };
00177 
00178 /**
00179  * Computes \f$ \int H'(\phi(x))) \frac{\nabla\phi(x)}{\|\nabla\phi(x)\|_\epsilon}\cdot \nabla\varphi_j dx \f$, where arg[0]=phi and HeavisideFunction=H.
00180  *
00181  * \author Berkels
00182  */
00183 template <typename ConfiguratorType, typename HeavisideFunctionType>
00184 class SecondPartOfFullVariationOfHeavisideLevelsetLengthEnergy
00185 : public aol::FENonlinDiffOpInterface<ConfiguratorType, SecondPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType> > {
00186 public:
00187   typedef typename ConfiguratorType::RealType RealType;
00188 protected:
00189   const HeavisideFunctionType &_heavisideFunction;
00190   const RealType _epsilonSqr;
00191 public:
00192   SecondPartOfFullVariationOfHeavisideLevelsetLengthEnergy( const typename ConfiguratorType::InitType &Initializer,
00193                                                             const HeavisideFunctionType &HeavisideFunction,
00194                                                             const RealType Epsilon )
00195   : aol::FENonlinDiffOpInterface< ConfiguratorType,
00196                                    SecondPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType> > ( Initializer ),
00197       _heavisideFunction( HeavisideFunction ),
00198       _epsilonSqr( aol::Sqr(Epsilon) ) {
00199   }
00200 
00201   inline void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
00202                                const typename ConfiguratorType::ElementType &El,
00203                                int QuadPoint,
00204                                const typename ConfiguratorType::DomVecType &/*RefCoord*/,
00205                                typename ConfiguratorType::VecType &NL ) const {
00206     const RealType hPrime = _heavisideFunction.evaluateDerivative(DiscFunc.evaluateAtQuadPoint( El, QuadPoint));
00207     DiscFunc.evaluateGradientAtQuadPoint( El, QuadPoint, NL );
00208     NL *= hPrime / (sqrt(NL.normSqr()+_epsilonSqr));
00209   }
00210 };
00211 
00212 /**
00213  * Calulates the variation of HeavisideLevelsetLengthEnergy including the derivative
00214  * of the Heaviside function part.
00215  *
00216  * This involves the second derivate of the Heaviside function. Therefore the modelling
00217  * should try to avoid to use this. Currently it's only meant to be used for testing
00218  * purposes.
00219  *
00220  * Instead of using this class, let the metric of the gradient flow handle the Heaviside
00221  * function part of the variation and use VariationOfHeavisideLevelsetLengthEnergy.
00222  *
00223  * \author Berkels
00224  */
00225 template <typename ConfiguratorType, typename HeavisideFunctionType, int numberOfLevelsetFunctions = 1>
00226 class FullVariationOfHeavisideLevelsetLengthEnergy
00227 : public Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
00228 public:
00229   typedef typename ConfiguratorType::RealType RealType;
00230 protected:
00231   const typename ConfiguratorType::InitType &_grid;
00232   const HeavisideFunctionType &_heavisideFunction;
00233   const RealType _epsilon;
00234 public:
00235   /**
00236    * \param Epsilon parameter to regularize the absolute value
00237    */
00238   FullVariationOfHeavisideLevelsetLengthEnergy( const typename ConfiguratorType::InitType &Initializer,
00239                                                 const HeavisideFunctionType &HeavisideFunction,
00240                                                 const RealType Epsilon = 0.1 )
00241     : _grid( Initializer ),
00242       _heavisideFunction( HeavisideFunction ),
00243       _epsilon( Epsilon )
00244   {
00245   }
00246 
00247   void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::MultiVector<RealType> &MDest ) const {
00248 
00249     FirstPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType>
00250       part1( _grid, _heavisideFunction, _epsilon );
00251     SecondPartOfFullVariationOfHeavisideLevelsetLengthEnergy<ConfiguratorType, HeavisideFunctionType>
00252       part2( _grid, _heavisideFunction, _epsilon );
00253 
00254     for( int i = 0; i < numberOfLevelsetFunctions; i++ ){
00255       part1.applyAdd( MArg[i], MDest[i] );
00256       part2.applyAdd( MArg[i], MDest[i] );
00257     }
00258   }
00259 };
00260 
00261 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00262 class HeavisideFENonlinIntegrationVectorInterface
00263 : public aol::FENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>, TargetOpType::NumOfComponents > {
00264 public:
00265   typedef typename ConfiguratorType::RealType RealType;
00266 protected:
00267   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00268   const HeavisideFunctionType &_heavisideFunction;
00269   const TargetOpType &_targetOp;
00270 public:
00271   HeavisideFENonlinIntegrationVectorInterface( const typename ConfiguratorType::InitType &Initializer,
00272                                                  const aol::Vector<RealType> &LevelsetFunction,
00273                                                  const HeavisideFunctionType &HeavisideFunction,
00274                                                  const TargetOpType &TargetOp )
00275   : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00276                                                HeavisideFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>,
00277                                                TargetOpType::NumOfComponents > ( Initializer ),
00278       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00279       _heavisideFunction( HeavisideFunction ),
00280       _targetOp( TargetOp ){
00281   }
00282 
00283   RealType evaluateIntegrand( const aol::auto_container<TargetOpType::NumOfComponents,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00284                               const typename ConfiguratorType::ElementType &El,
00285                               int QuadPoint, const typename ConfiguratorType::VecType &RefCoord ) const {
00286     return ( _heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint))
00287              * _targetOp.evaluateIntegrand( discrFuncs, El, QuadPoint, RefCoord ) );
00288   }
00289 };
00290 
00291 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00292 class OneMinusHeavisideFENonlinIntegrationVectorInterface
00293 : public aol::FENonlinIntegrationVectorInterface<ConfiguratorType, OneMinusHeavisideFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>, TargetOpType::NumOfComponents > {
00294 public:
00295   typedef typename ConfiguratorType::RealType RealType;
00296 protected:
00297   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00298   const HeavisideFunctionType &_heavisideFunction;
00299   const TargetOpType &_targetOp;
00300 public:
00301   OneMinusHeavisideFENonlinIntegrationVectorInterface( const typename ConfiguratorType::InitType &Initializer,
00302                                                          const aol::Vector<RealType> &LevelsetFunction,
00303                                                          const HeavisideFunctionType &HeavisideFunction,
00304                                                          const TargetOpType &TargetOp )
00305   : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00306                                                OneMinusHeavisideFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>,
00307                                                TargetOpType::NumOfComponents > ( Initializer ),
00308       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00309       _heavisideFunction( HeavisideFunction ),
00310       _targetOp( TargetOp ){
00311   }
00312 
00313   RealType evaluateIntegrand( const aol::auto_container<TargetOpType::NumOfComponents,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00314                               const typename ConfiguratorType::ElementType &El,
00315                               int QuadPoint, const typename ConfiguratorType::VecType &RefCoord ) const {
00316     return ( (1. -_heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint)))
00317              * _targetOp.evaluateIntegrand( discrFuncs, El, QuadPoint, RefCoord ) );
00318   }
00319 };
00320 
00321 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00322 class HeavisideAndOneMinusHFENonlinIntegrationVectorInterface
00323 : public aol::FENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideAndOneMinusHFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>, TargetOpType::NumOfComponents > {
00324 public:
00325   typedef typename ConfiguratorType::RealType RealType;
00326 protected:
00327   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00328   const HeavisideFunctionType &_heavisideFunction;
00329   const TargetOpType &_hsTargetOp;
00330   const TargetOpType &_oneMinusHSTargetOp;
00331 public:
00332   HeavisideAndOneMinusHFENonlinIntegrationVectorInterface( const typename ConfiguratorType::InitType &Initializer,
00333                                                              const aol::Vector<RealType> &LevelsetFunction,
00334                                                              const HeavisideFunctionType &HeavisideFunction,
00335                                                              const TargetOpType &HSTargetOp,
00336                                                              const TargetOpType &OneMinusHSTargetOp )
00337   : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00338                                                HeavisideAndOneMinusHFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>,
00339                                                TargetOpType::NumOfComponents > ( Initializer ),
00340       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00341       _heavisideFunction( HeavisideFunction ),
00342       _hsTargetOp( HSTargetOp ),
00343       _oneMinusHSTargetOp( OneMinusHSTargetOp ){
00344   }
00345 
00346   RealType evaluateIntegrand( const aol::auto_container<TargetOpType::NumOfComponents,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00347                               const typename ConfiguratorType::ElementType &El,
00348                               int QuadPoint, const typename ConfiguratorType::VecType &RefCoord ) const {
00349     const RealType h = _heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint));
00350     const RealType hsOp = _hsTargetOp.evaluateIntegrand( discrFuncs, El, QuadPoint, RefCoord );
00351     const RealType oneMinusHSOp = _oneMinusHSTargetOp.evaluateIntegrand( discrFuncs, El, QuadPoint, RefCoord );
00352     return ( h * hsOp + ( 1. - h ) * oneMinusHSOp );
00353   }
00354 };
00355 
00356 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00357 class OneMinusHeavisideFENonlinOpInterface : public aol::FENonlinOpInterface< ConfiguratorType, OneMinusHeavisideFENonlinOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > {
00358   typedef typename ConfiguratorType::RealType RealType;
00359 protected:
00360   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00361   const HeavisideFunctionType &_heavisideFunction;
00362   const TargetOpType &_targetOp;
00363 public:
00364   OneMinusHeavisideFENonlinOpInterface( const typename ConfiguratorType::InitType &Initializer,
00365                                         const aol::Vector<RealType> &LevelsetFunction,
00366                                         const HeavisideFunctionType &HeavisideFunction,
00367                                         const TargetOpType &TargetOp )
00368     : aol::FENonlinOpInterface< ConfiguratorType,
00369                                  OneMinusHeavisideFENonlinOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > ( Initializer ),
00370       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00371       _heavisideFunction( HeavisideFunction ),
00372       _targetOp( TargetOp ){
00373   }
00374 
00375   void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc, 
00376                         const typename ConfiguratorType::ElementType &El, 
00377                         int QuadPoint, const typename ConfiguratorType::VecType &RefCoord,
00378                         typename ConfiguratorType::RealType &NL ) const {
00379     _targetOp.getNonlinearity( DiscFunc, El, QuadPoint, RefCoord, NL );
00380     NL *= (1. -_heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint)));
00381   }
00382 };
00383 
00384 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00385 class HeavySideOneMinusHFENonlinOpInterface : public aol::FENonlinOpInterface< ConfiguratorType, HeavySideOneMinusHFENonlinOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > {
00386   typedef typename ConfiguratorType::RealType RealType;
00387 protected:
00388   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00389   const HeavisideFunctionType &_heavisideFunction;
00390   const TargetOpType &_hsTargetOp;
00391   const TargetOpType &_oneMinusHSTargetOp;
00392 public:
00393   HeavySideOneMinusHFENonlinOpInterface( const typename ConfiguratorType::InitType &Initializer,
00394                                          const aol::Vector<RealType> &LevelsetFunction,
00395                                          const HeavisideFunctionType &HeavisideFunction,
00396                                          const TargetOpType &HSTargetOp,
00397                                          const TargetOpType &OneMinusHSTargetOp )
00398     : aol::FENonlinOpInterface< ConfiguratorType,
00399                                  HeavySideOneMinusHFENonlinOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > ( Initializer ),
00400       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00401       _heavisideFunction( HeavisideFunction ),
00402       _hsTargetOp( HSTargetOp ),
00403       _oneMinusHSTargetOp( OneMinusHSTargetOp ){
00404   }
00405 
00406   void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc, 
00407                         const typename ConfiguratorType::ElementType &El, 
00408                         int QuadPoint, const typename ConfiguratorType::VecType &RefCoord,
00409                         typename ConfiguratorType::RealType &NL ) const {
00410     const RealType h = _heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint));
00411     RealType hsOp = 0.;
00412     _hsTargetOp.getNonlinearity( DiscFunc, El, QuadPoint, RefCoord, hsOp );
00413     RealType oneMinusHSOp = 0.;
00414     _oneMinusHSTargetOp.getNonlinearity( DiscFunc, El, QuadPoint, RefCoord, oneMinusHSOp );
00415     NL = ( h * hsOp + ( 1. - h ) * oneMinusHSOp );
00416   }
00417 };
00418 
00419 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
00420 class HeavisideAndOneMinusHFENonlinDiffOpInterface
00421 : public aol::FENonlinDiffOpInterface<ConfiguratorType, HeavisideAndOneMinusHFENonlinDiffOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > {
00422 public:
00423   typedef typename ConfiguratorType::RealType RealType;
00424 protected:
00425   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00426   const HeavisideFunctionType &_heavisideFunction;
00427   const TargetOpType &_hsTargetOp;
00428   const TargetOpType &_oneMinusHSTargetOp;
00429 public:
00430   HeavisideAndOneMinusHFENonlinDiffOpInterface( const typename ConfiguratorType::InitType &Initializer,
00431                                                  const aol::Vector<RealType> &LevelsetFunction,
00432                                                  const HeavisideFunctionType &HeavisideFunction,
00433                                                  const TargetOpType &HSTargetOp,
00434                                                  const TargetOpType &OneMinusHSTargetOp )
00435   : aol::FENonlinDiffOpInterface< ConfiguratorType,
00436                                    HeavisideAndOneMinusHFENonlinDiffOpInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType> > ( Initializer ),
00437       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00438       _heavisideFunction( HeavisideFunction ),
00439       _hsTargetOp( HSTargetOp ),
00440       _oneMinusHSTargetOp( OneMinusHSTargetOp ){
00441   }
00442 
00443   inline void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
00444                                const typename ConfiguratorType::ElementType &El,
00445                                int QuadPoint,
00446                                const typename ConfiguratorType::VecType &RefCoord,
00447                                typename ConfiguratorType::VecType &NL ) const {
00448     const RealType h = _heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint));
00449     typename ConfiguratorType::VecType temp;
00450 
00451     _hsTargetOp.getNonlinearity( DiscFunc, El, QuadPoint, RefCoord, temp );
00452     temp *= h;
00453     NL = temp;
00454 
00455     _oneMinusHSTargetOp.getNonlinearity( DiscFunc, El, QuadPoint, RefCoord, temp );
00456     temp *= (1.-h);
00457     NL += temp;
00458   }
00459 };
00460 
00461 template <typename ConfiguratorType, int NumCompArg, int NumCompDest, typename HeavisideFunctionType, typename TargetOpType>
00462 class OneMinusHeavisideFENonlinVectorDiffOpInterface : public aol::FENonlinVectorDiffOpInterface< ConfiguratorType, NumCompArg, NumCompDest, OneMinusHeavisideFENonlinVectorDiffOpInterface<ConfiguratorType, NumCompArg, NumCompDest, HeavisideFunctionType, TargetOpType> > {
00463   typedef typename ConfiguratorType::RealType RealType;
00464 protected:
00465   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00466   const HeavisideFunctionType &_heavisideFunction;
00467   const TargetOpType &_targetOp;
00468 public:
00469   OneMinusHeavisideFENonlinVectorDiffOpInterface ( const typename ConfiguratorType::InitType &Initializer,
00470                                                    const aol::Vector<RealType> &LevelsetFunction,
00471                                                    const HeavisideFunctionType &HeavisideFunction,
00472                                                    const TargetOpType &TargetOp )
00473     : aol::FENonlinVectorDiffOpInterface< ConfiguratorType,
00474                                           NumCompArg,
00475                                           NumCompDest,
00476                                           OneMinusHeavisideFENonlinVectorDiffOpInterface<ConfiguratorType, NumCompArg, NumCompDest, HeavisideFunctionType, TargetOpType> > ( Initializer ),
00477       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00478       _heavisideFunction( HeavisideFunction ),
00479       _targetOp( TargetOp ){
00480   }
00481 
00482   void getNonlinearity ( const auto_container<NumCompArg, aol::DiscreteFunctionDefault<ConfiguratorType> > &DiscFuncs,
00483                          const typename ConfiguratorType::ElementType &El,
00484                          int QuadPoint, const typename ConfiguratorType::VecType &RefCoord,
00485                          aol::Mat<NumCompDest, NumCompDest, typename ConfiguratorType::RealType> &NL ) const {
00486     _targetOp.getNonlinearity( DiscFuncs, El, QuadPoint, RefCoord, NL );
00487     NL *= (1. -_heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint)));
00488   }
00489 };
00490 
00491 /**
00492  * Computes \f$ \int (1-H(\phi(x))) f(x) dx \f$, where arg[0]=f, LevelsetFunction=phi and HeavisideFunction=H.
00493  *
00494  * \author Berkels
00495  */
00496 template <typename ConfiguratorType, typename HeavisideFunctionType>
00497 class OneMinusHeavisideIntegrateFEFunction
00498   : public FENonlinIntegrationScalarInterface<ConfiguratorType, OneMinusHeavisideIntegrateFEFunction<ConfiguratorType, HeavisideFunctionType> > {
00499 
00500   typedef typename ConfiguratorType::RealType RealType;
00501   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00502   const HeavisideFunctionType &_heavisideFunction;
00503 
00504 public:
00505   OneMinusHeavisideIntegrateFEFunction ( const typename ConfiguratorType::InitType &Initializer,
00506                                          const aol::Vector<RealType> &LevelsetFunction,
00507                                          const HeavisideFunctionType &HeavisideFunction )
00508     : FENonlinIntegrationScalarInterface<ConfiguratorType, OneMinusHeavisideIntegrateFEFunction<ConfiguratorType, HeavisideFunctionType> > ( Initializer ),
00509       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00510       _heavisideFunction( HeavisideFunction ) { }
00511 
00512   RealType evaluateIntegrand ( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
00513                                const typename ConfiguratorType::ElementType &El,
00514                                int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
00515     const RealType h = _heavisideFunction.evaluate(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint));
00516     return (ZOTrait<RealType>::one-h) * DiscFunc.evaluateAtQuadPoint ( El, QuadPoint );
00517   }
00518 };
00519 
00520 /**
00521  * Computes \f$ - \int H'(\phi(x)) f(x) \varphi_j dx \f$, where arg[0]=f, LevelsetFunction=phi and HeavisideFunction=H.
00522  *
00523  * \author Berkels
00524  */
00525 template <typename ConfiguratorType, typename HeavisideFunctionType>
00526 class FullVariationOfOneMinusHeavisideIntegrateFEFunction
00527   : public aol::FENonlinOpInterface< ConfiguratorType, FullVariationOfOneMinusHeavisideIntegrateFEFunction<ConfiguratorType, HeavisideFunctionType> > {
00528 
00529   typedef typename ConfiguratorType::RealType RealType;
00530   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrLevelsetFunction;
00531   const HeavisideFunctionType &_heavisideFunction;
00532 
00533 public:
00534   FullVariationOfOneMinusHeavisideIntegrateFEFunction ( const typename ConfiguratorType::InitType &Initializer,
00535                                                         const aol::Vector<RealType> &LevelsetFunction,
00536                                                         const HeavisideFunctionType &HeavisideFunction )
00537     : aol::FENonlinOpInterface< ConfiguratorType, FullVariationOfOneMinusHeavisideIntegrateFEFunction<ConfiguratorType, HeavisideFunctionType> > ( Initializer ),
00538       _discrLevelsetFunction( Initializer, LevelsetFunction ),
00539       _heavisideFunction( HeavisideFunction ) { }
00540 
00541   void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc, 
00542                         const typename ConfiguratorType::ElementType &El, 
00543                         int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/,
00544                         typename ConfiguratorType::RealType &NL ) const {
00545     const RealType hPrime = _heavisideFunction.evaluateDerivative(_discrLevelsetFunction.evaluateAtQuadPoint( El, QuadPoint));
00546     NL = -1 * hPrime * DiscFunc.evaluateAtQuadPoint ( El, QuadPoint );
00547   }
00548 };
00549 
00550 /**
00551  * \author Berkels
00552  */
00553 template <typename RealType>
00554 class SquareFunction{
00555 public:
00556   static RealType evaluate( const RealType x ) {
00557     return aol::Sqr ( x );
00558   }
00559   static RealType evaluateDerivative( const RealType x ) {
00560     return 2 * x;
00561   }
00562 };
00563 
00564 /**
00565  * \brief Heaviside function replacement to convert the Chan Vese model to the one of Chan Esedoglu Nikolova with exact penalty.
00566  *
00567  * \author Berkels
00568  */
00569 template <typename RealType, typename _ScalingFunctionType = aol::SquareFunction<RealType> >
00570 class IdentityFunction{
00571 public:
00572   typedef _ScalingFunctionType ScalingFunctionType;
00573 
00574   IdentityFunction( const RealType /*Epsilon*/ = 0. ) {
00575   }
00576   static RealType evaluate( const RealType x ) {
00577     return x;
00578   }
00579   static RealType evaluateDerivative( const RealType /*x*/ ) {
00580     return aol::NumberTrait<RealType>::one;
00581   }
00582   //! In a proper Chan Vese type model, you should not need the second derivative of the Heaviside function.
00583   static RealType evaluateSecondDerivative( const RealType /*x*/ ) {
00584     return aol::NumberTrait<RealType>::zero;
00585   }
00586 };
00587 
00588 template <typename RealType>
00589 class ArcTanHeavisideFunction{
00590   const RealType _epsilon;
00591   const RealType _oneOverPi;
00592 public:
00593   typedef IdentityFunction<RealType> ScalingFunctionType;
00594 
00595   ArcTanHeavisideFunction( const RealType Epsilon )
00596    : _epsilon(Epsilon),
00597     _oneOverPi( 1./aol::NumberTrait<RealType>::getPi() )
00598   {
00599   }
00600   RealType evaluate( const RealType x ) const{
00601     return 0.5 + _oneOverPi*atan( x/_epsilon );
00602   }
00603   RealType evaluateDerivative( const RealType x ) const{
00604     return _oneOverPi/(_epsilon +  aol::Sqr(x)/_epsilon );
00605   }
00606   //! In a proper Chan Vese type model, you should not need the second derivative of the Heaviside function.
00607   RealType evaluateSecondDerivative( const RealType x ) const{
00608     return -2*x*_oneOverPi/(_epsilon*aol::Sqr(_epsilon +  aol::Sqr(x)/_epsilon ));
00609   }
00610 };
00611 
00612 template <typename RealType>
00613 class OneMinusArcTanHeavisideFunction{
00614   const RealType _epsilon;
00615   const RealType _oneOverPi;
00616 public:
00617   OneMinusArcTanHeavisideFunction( const RealType Epsilon )
00618    : _epsilon(Epsilon),
00619      _oneOverPi( 1./aol::NumberTrait<RealType>::getPi() )
00620   {
00621   }
00622   RealType evaluate( const RealType x ) const{
00623     return 0.5 - _oneOverPi*atan( x/_epsilon );
00624   }
00625   RealType evaluateDerivative( const RealType x ) const{
00626     return -1.*_oneOverPi/(_epsilon +  aol::Sqr(x)/_epsilon );
00627   }
00628 };
00629 
00630 template <typename RealType>
00631 class PolynomialHeavisideFunction{
00632   const RealType _epsilon;
00633 public:
00634   PolynomialHeavisideFunction( const RealType Epsilon )
00635    : _epsilon( Epsilon )
00636   {
00637   }
00638   RealType evaluate( const RealType x ) const{
00639     if( x > _epsilon*0.5){
00640       return 1.;
00641     }
00642     else{
00643       if( x < (-1.)*_epsilon*0.5 ){
00644         return 0.;
00645       }
00646       else{
00647         const RealType position = -x/_epsilon+0.5;
00648         return aol::Sqr(position)*(2.*position-3.)+1.;
00649       }
00650     }
00651   }
00652   RealType evaluateDerivative( const RealType x ) const{
00653     if( fabs(x) > _epsilon*0.5){
00654       return 0.;
00655     }
00656     else{
00657       const RealType position = -x/_epsilon+0.5;
00658       return -6.*position*(position-1.)/_epsilon;
00659     }
00660   }
00661 };
00662 
00663 //! Prints central differential quotioent and implemented derivative of the Heaviside Funtion.
00664 //! If they don't match, it's likely that there is a bug in the Heaviside Funtion.
00665 //! Also calculates the integral over the derivative of the Heaviside Function.
00666 //! This should be one, since the derivatives approximates the delta distribution.
00667 template <typename RealType, typename HeavisideFunctionType>
00668 void consistencyCheckOfHeavisideFuntion( HeavisideFunctionType &H ){
00669   const RealType epsilon = 0.0001;
00670   aol::MixedFormat format ( 4, 12 );
00671   for( int i = -100; i < 100; i++ ){
00672     const RealType x = static_cast<RealType>(i)/1000.;
00673     const RealType differentialQuotient = (H.evaluate( x + epsilon ) - H.evaluate( x - epsilon ))/ (2.*epsilon);
00674     const RealType derivative = H.evaluateDerivative( x );
00675     cerr << format(differentialQuotient) << " " << format(derivative) << " " << format(differentialQuotient-derivative) << endl;
00676   }
00677   const RealType width = 1000.;
00678   const int N = 1000000;
00679   RealType integral = 0.;
00680   const RealType h = 2.*width/static_cast<RealType>(N-1);
00681   for( int i = 0; i < N; i++ ){
00682     const RealType x = -width + static_cast<RealType>(i)*h;
00683     integral += h * H.evaluateDerivative( x );
00684   }
00685   cerr << "Integral of derivative = " << integral << endl;
00686 }
00687 
00688 
00689 /**
00690  * \brief Heaviside function replacement to convert the Chan Vese model to the one of Chan Esedoglu Nikolova.
00691  *
00692  * Attention: You can't use the Armijo rule in a gradient descent for an energy using this function:
00693  * evaluate is clamped, evaluateDerivative is not, in this sense the derivative doesn't match the energy
00694  * and therefore the Armijo rule will reject the negative gradient as descent direction.
00695  *
00696  * \author Berkels
00697  */
00698 template <typename RealType>
00699 class ClampedIdentityFunction{
00700 public:
00701   ClampedIdentityFunction( const RealType /*Epsilon*/ = 0. ) {
00702   }
00703   static RealType evaluate( const RealType x ) {
00704     return aol::Clamp<RealType>( x, 0, 1 );
00705   }
00706   static RealType evaluateDerivative( const RealType /*x*/ ) {
00707     return aol::NumberTrait<RealType>::one;
00708   }
00709   //! In a proper Chan Vese type model, you should not need the second derivative of the Heaviside function.
00710   static RealType evaluateSecondDerivative( const RealType /*x*/ ) {
00711     return aol::NumberTrait<RealType>::zero;
00712   }
00713 };
00714 
00715 /**
00716  * Regularized version of \f$ \nu(x)=max\{0,2|x-\frac{1}{2}|-1\}\f$.
00717  *
00718  * Used as penalty term in the Chan Esedoglu Nikolova model.
00719  *
00720  * \author Berkels
00721  */
00722 template <typename RealType>
00723 class ZeroOneIntervalPenaltyFunction {
00724   const RealType _oneOverFourEpsilon, _oneOverTwoEpsilon, _epsilonMinusOneHalf, _oneHalfMinusEpsilon, _oneHalfPlusEpsilon;
00725 public:
00726   ZeroOneIntervalPenaltyFunction ( const RealType Epsilon )
00727    : _oneOverFourEpsilon( aol::NumberTrait<RealType>::one/(4*Epsilon) ),
00728      _oneOverTwoEpsilon( aol::NumberTrait<RealType>::one/(2*Epsilon) ),
00729      _epsilonMinusOneHalf ( Epsilon - 0.5 ),
00730      _oneHalfMinusEpsilon ( 0.5 - Epsilon ),
00731      _oneHalfPlusEpsilon ( 0.5 + Epsilon )
00732   {
00733   }
00734   RealType evaluate( const RealType x ) const {
00735     const RealType y = aol::Abs( x - 0.5 );
00736     if ( y < _oneHalfMinusEpsilon )
00737       return aol::NumberTrait<RealType>::zero;
00738     else if ( y > _oneHalfPlusEpsilon )
00739       return ( y - 0.5 );
00740     else
00741       return _oneOverFourEpsilon * aol::Sqr( y + _epsilonMinusOneHalf );
00742   }
00743   RealType evaluateDerivative( const RealType x ) const{
00744     const RealType y = aol::Abs( x - 0.5 );
00745     if ( y < _oneHalfMinusEpsilon )
00746       return aol::NumberTrait<RealType>::zero;
00747     else if ( y > _oneHalfPlusEpsilon )
00748       return ( aol::signum( x - 0.5 ) * aol::NumberTrait<RealType>::one );
00749     else
00750       return aol::signum( x - 0.5 ) * _oneOverTwoEpsilon * ( y + _epsilonMinusOneHalf );
00751   }
00752 };
00753 
00754 /**
00755  * \f$ \nu(x)=max\{0,-x|x|, (x-1)|x-1|\} \f$.
00756  *
00757  * \author Berkels
00758  */
00759 template <typename RealType>
00760 class QuadraticZeroOneIntervalPenaltyFunction {
00761 public:
00762   QuadraticZeroOneIntervalPenaltyFunction ( const RealType /*Epsilon*/ ) {}
00763 
00764   RealType evaluate( const RealType x ) const {
00765     const RealType y = aol::Abs( x - 0.5 );
00766     if ( y < 0.5 )
00767       return aol::NumberTrait<RealType>::zero;
00768     else
00769       return ( aol::Sqr ( y - 0.5 ) );
00770   }
00771   RealType evaluateDerivative( const RealType x ) const{
00772     const RealType y = aol::Abs( x - 0.5 );
00773     if ( y < 0.5 )
00774       return aol::NumberTrait<RealType>::zero;
00775     else
00776       return aol::signum( x - 0.5 ) * 2 * ( y -0.5 );
00777   }
00778 };
00779 
00780 /**
00781  * \f$ \nu(x)=max\{0,-x^3, (x-1)^3\} \f$.
00782  *
00783  * \author Berkels
00784  */
00785 template <typename RealType>
00786 class CubicZeroOneIntervalPenaltyFunction {
00787 public:
00788   CubicZeroOneIntervalPenaltyFunction ( const RealType /*Epsilon*/ ) {}
00789   RealType evaluate( const RealType x ) const {
00790     const RealType y = aol::Abs( x - 0.5 );
00791     if ( y < 0.5 )
00792       return aol::NumberTrait<RealType>::zero;
00793     else
00794       return ( aol::Cub ( y - 0.5 ) );
00795   }
00796   RealType evaluateDerivative( const RealType x ) const{
00797     const RealType y = aol::Abs( x - 0.5 );
00798     if ( y < 0.5 )
00799       return aol::NumberTrait<RealType>::zero;
00800     else
00801       return aol::signum( x - 0.5 ) * 3 * aol::Sqr( y -0.5 );
00802   }
00803   RealType evaluateSecondDerivative( const RealType x ) const {
00804     const RealType y = aol::Abs( x - 0.5 );
00805     if ( y < 0.5 )
00806       return aol::NumberTrait<RealType>::zero;
00807     else
00808       return aol::signum( x - 0.5 ) * 6 * ( y -0.5 );
00809   }
00810 };
00811 
00812 /**
00813  * calculates \f$ \int_\Omega \nu(u(x)) dx \f$, where \f$ u \f$ is given by Arg and \f$ \nu \f$ by PenaltyFunction.
00814  *
00815  * \author Berkels
00816  */
00817 template <typename ConfiguratorType, typename PenaltyFunctionType>
00818 class PointwiseConstraintPenaltyEnergy
00819   : public FENonlinIntegrationScalarInterface<ConfiguratorType, PointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> > {
00820   typedef typename ConfiguratorType::RealType RealType;
00821   const PenaltyFunctionType &_penaltyFunction;
00822 public:
00823   PointwiseConstraintPenaltyEnergy ( const typename ConfiguratorType::InitType &Initializer,
00824                                      const PenaltyFunctionType &PenaltyFunction )
00825     : FENonlinIntegrationScalarInterface<ConfiguratorType, PointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> > ( Initializer ),
00826       _penaltyFunction( PenaltyFunction ) { }
00827 
00828   RealType evaluateIntegrand ( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc,
00829                                const typename ConfiguratorType::ElementType &El,
00830                                int QuadPoint, const typename ConfiguratorType::DomVecType &/*RefCoord*/ ) const {
00831     return _penaltyFunction.evaluate ( DiscFunc.evaluateAtQuadPoint( El, QuadPoint ) );
00832   }
00833 };
00834 
00835 /**
00836  * calculates \f$ \int_\Omega \nu^\prime(u(x)) \varphi_i dx \f$, where \f$ u \f$ is given by Arg and \f$ \nu \f$ by PenaltyFunction.
00837  *
00838  * \author Berkels
00839  */
00840 template <typename ConfiguratorType, typename PenaltyFunctionType>
00841 class VariationOfPointwiseConstraintPenaltyEnergy
00842   : public aol::FENonlinOpInterface< ConfiguratorType, VariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> > {
00843   typedef typename ConfiguratorType::RealType RealType;
00844   const PenaltyFunctionType &_penaltyFunction;
00845 public:
00846   VariationOfPointwiseConstraintPenaltyEnergy( const typename ConfiguratorType::InitType &Initializer,
00847                                                const PenaltyFunctionType &PenaltyFunction )
00848     : aol::FENonlinOpInterface< ConfiguratorType,
00849                                 VariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> > ( Initializer ),
00850       _penaltyFunction( PenaltyFunction ) {
00851   }
00852 
00853   void getNonlinearity( const aol::DiscreteFunctionDefault<ConfiguratorType> &DiscFunc, 
00854                         const typename ConfiguratorType::ElementType &El, 
00855                         int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/,
00856                         typename ConfiguratorType::RealType &NL ) const {
00857     NL = _penaltyFunction.evaluateDerivative ( DiscFunc.evaluateAtQuadPoint( El, QuadPoint ) );
00858   }
00859 };
00860 
00861 /**
00862  * calculates \f$ \int_\Omega \nu^{\prime\prime}(u(x)) \varphi_i(x) \varphi_j(x) dx \f$, where \f$ u \f$ is given by UDofs and \f$ \nu \f$ by PenaltyFunction.
00863  *
00864  * \author Berkels
00865  */
00866 template <typename ConfiguratorType, typename PenaltyFunctionType>
00867 class SecondVariationOfPointwiseConstraintPenaltyEnergy : public aol::FELinScalarWeightedMassInterface<ConfiguratorType, SecondVariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> > {
00868 public:
00869   typedef typename ConfiguratorType::RealType RealType;
00870 protected:
00871   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrU;
00872   const PenaltyFunctionType &_penaltyFunction;
00873 public:
00874   SecondVariationOfPointwiseConstraintPenaltyEnergy( const typename ConfiguratorType::InitType &Initializer, 
00875                                                      const aol::Vector<RealType> &UDofs,
00876                                                      const PenaltyFunctionType &PenaltyFunction,
00877                                                      aol::OperatorType OpType = aol::ONTHEFLY ) 
00878     : aol::FELinScalarWeightedMassInterface<ConfiguratorType, SecondVariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> >( Initializer, OpType ),
00879       _discrU( Initializer, UDofs ),
00880       _penaltyFunction( PenaltyFunction )
00881   {
00882   }
00883   inline RealType getCoeff( const qc::Element &El, int QuadPoint, const typename ConfiguratorType::VecType& ) const { 
00884     return _penaltyFunction.evaluateSecondDerivative ( _discrU.evaluateAtQuadPoint ( El, QuadPoint ) );
00885   }
00886 };
00887 
00888 /**
00889  * \author Berkels
00890  */
00891 template <typename ConfiguratorType, typename PenaltyFunctionType, typename MatrixType = qc::FastUniformGridMatrix<typename ConfiguratorType::RealType,ConfiguratorType::Dim> >
00892 class SecondVariationOfEsedogluSegmentationModelWithPenalty : public aol::Op<aol::Vector<typename ConfiguratorType::RealType>, MatrixType > {
00893   typedef typename ConfiguratorType::RealType RealType;
00894   const typename ConfiguratorType::InitType &_grid;
00895   const PenaltyFunctionType &_penaltyFunction;
00896   const RealType _lengthFactor, _penaltyFactor, _epsilon;
00897 public:
00898   SecondVariationOfEsedogluSegmentationModelWithPenalty ( const typename ConfiguratorType::InitType &Initializer,
00899                                                           const PenaltyFunctionType &PenaltyFunction,
00900                                                           const RealType LengthFactor,
00901                                                           const RealType PenaltyFactor,
00902                                                           const RealType Epsilon )
00903   : _grid ( Initializer ),
00904     _penaltyFunction ( PenaltyFunction ),
00905     _lengthFactor ( LengthFactor ),
00906     _penaltyFactor ( PenaltyFactor ),
00907     _epsilon ( Epsilon )
00908   {
00909   }
00910   virtual void apply( const aol::Vector<RealType> &Arg, MatrixType &Dest ) const {
00911     Dest.setZero();
00912     aol::SecondVariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> DDEPenalty ( _grid, Arg, _penaltyFunction );
00913     DDEPenalty.assembleAddMatrix ( Dest, _penaltyFactor );
00914 
00915     aol::IsotropicMatrixStiffOp<ConfiguratorType> isotropicMatrixStiffOp ( _grid, Arg, _epsilon );
00916     isotropicMatrixStiffOp.assembleAddMatrix( Dest,_lengthFactor ); 
00917   }
00918   virtual void applyAdd( const aol::Vector<RealType> &/*Arg*/, MatrixType &/*Dest*/ ) const {
00919     throw aol::UnimplementedCodeException( "Not implemented", __FILE__, __LINE__);
00920   }
00921 };
00922 
00923 /**
00924  * calculates \f$ \sum_j\int_\Omega \nu(u_j(x)) dx \f$, where \f$ u_j \f$ is given by MArg and \f$ \nu \f$ by PenaltyFunction.
00925  *
00926  * \author Berkels
00927  */
00928 template <typename ConfiguratorType, typename PenaltyFunctionType, int numberOfFunctions>
00929 class PointwiseConstraintPenaltyEnergyMulti
00930   : public aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00931                                                     PointwiseConstraintPenaltyEnergyMulti<ConfiguratorType, PenaltyFunctionType, numberOfFunctions>,
00932                                                     numberOfFunctions > {
00933   typedef typename ConfiguratorType::RealType RealType;
00934   const PenaltyFunctionType &_penaltyFunction;
00935 public:
00936   PointwiseConstraintPenaltyEnergyMulti ( const typename ConfiguratorType::InitType &Initializer,
00937                                           const PenaltyFunctionType &PenaltyFunction )
00938     : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
00939                                                PointwiseConstraintPenaltyEnergyMulti<ConfiguratorType, PenaltyFunctionType, numberOfFunctions>,
00940                                                numberOfFunctions > ( Initializer ),
00941       _penaltyFunction( PenaltyFunction ) {
00942   }
00943 
00944   RealType evaluateIntegrand( const aol::auto_container<numberOfFunctions,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
00945                               const typename ConfiguratorType::ElementType &El,
00946                               int QuadPoint, const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
00947     RealType integrand = 0.;
00948     for( int i = 0; i < numberOfFunctions; i++ ){
00949       integrand += _penaltyFunction.evaluate ( discrFuncs[i].evaluateAtQuadPoint( El, QuadPoint ) );
00950     }
00951     return integrand;
00952   }
00953 };
00954 
00955 /**
00956  * calculates \f$ \int_\Omega \nu^\prime(u_j(x)) \varphi_i dx \f$, where \f$ u_j \f$ is given by MArg and \f$ \nu \f$ by PenaltyFunction.
00957  *
00958  * \author Berkels
00959  */
00960 template <typename ConfiguratorType, typename PenaltyFunctionType, int numberOfFunctions>
00961 class VariationOfPointwiseConstraintPenaltyEnergyMulti
00962 : public Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
00963   typedef typename ConfiguratorType::RealType RealType;
00964   const VariationOfPointwiseConstraintPenaltyEnergy<ConfiguratorType, PenaltyFunctionType> _DEPenaltySingle;
00965 public:
00966   VariationOfPointwiseConstraintPenaltyEnergyMulti ( const typename ConfiguratorType::InitType &Initializer,
00967                                                      const PenaltyFunctionType &PenaltyFunction )
00968     : _DEPenaltySingle ( Initializer, PenaltyFunction )
00969   {
00970   }
00971 
00972   void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::MultiVector<RealType> &MDest ) const {
00973     for( int i = 0; i < numberOfFunctions; i++ ){
00974       _DEPenaltySingle.applyAdd( MArg[i], MDest[i] );
00975     }
00976   }
00977 };
00978 
00979 /**
00980  * \author Berkels
00981  */
00982 template <typename ConfiguratorType, typename PenaltyFunctionType, int numberOfFunctions, typename MatrixType = qc::FastUniformGridMatrix<typename ConfiguratorType::RealType,ConfiguratorType::Dim> >
00983 class SecondVariationOfEsedogluSegmentationModelWithPenaltyMulti
00984 : public aol::Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::SparseBlockMatrix<MatrixType> > {
00985   typedef typename ConfiguratorType::RealType RealType;
00986   const typename ConfiguratorType::InitType &_grid;
00987   const aol::SecondVariationOfEsedogluSegmentationModelWithPenalty<ConfiguratorType, PenaltyFunctionType, MatrixType> _DFSingle;
00988 public:
00989   SecondVariationOfEsedogluSegmentationModelWithPenaltyMulti( const typename ConfiguratorType::InitType &Initializer,
00990                                                               const PenaltyFunctionType &PenaltyFunction,
00991                                                               const RealType LengthFactor,
00992                                                               const RealType PenaltyFactor,
00993                                                               const RealType Epsilon )
00994     : _grid ( Initializer ),
00995       _DFSingle ( Initializer, PenaltyFunction, LengthFactor, PenaltyFactor, Epsilon )
00996   {
00997   }
00998 
00999   virtual void apply( const aol::MultiVector<RealType> &MArg,  aol::SparseBlockMatrix<MatrixType> &MDest ) const{
01000     MDest.deleteBlockEntries();
01001 
01002     for( int i = 0; i < numberOfFunctions; i++ ){
01003       MatrixType &diagBlockEntryI = MDest.allocateMatrix( i, i, _grid );
01004       _DFSingle.apply( MArg[i], diagBlockEntryI );
01005     }
01006   }
01007 
01008   virtual void applyAdd( const aol::MultiVector<RealType> &/*MArg*/, aol::SparseBlockMatrix<MatrixType> &/*MDest*/ ) const{
01009     throw aol::UnimplementedCodeException( "Not implemented", __FILE__, __LINE__);
01010   }
01011 };
01012 
01013 /**
01014  * \brief Iterates over all possible vectors of length NumberOfComponents with component values of zero or one.
01015  *
01016  * A vector of length l with component values of zero or one can be considered as an element of the power set
01017  * of a set with l elements, therefore the decision for the name of this class.
01018  *
01019  * \author Berkels
01020  */
01021 class PowerSetIterator{
01022   const int _numberOfComponents;
01023   const int _numberOfModifiers;
01024   int _currentPostion;
01025 public:
01026   PowerSetIterator( int NumberOfComponents )
01027     : _numberOfComponents( NumberOfComponents ),
01028       _numberOfModifiers( 1<<_numberOfComponents ),
01029       _currentPostion( 0 )
01030   {
01031   }
01032   //! Step to the next vector
01033   void increment(){
01034     if( _currentPostion < _numberOfModifiers )
01035       _currentPostion++;
01036   }
01037   //! Returns true, if the iterations over all possible vectors is finished.
01038   bool end() const{
01039     return ( _currentPostion == _numberOfModifiers );
01040   }
01041   //! returns the I-th component of the current vector
01042   int getComponent( int I ) const{
01043     // 2^I is equal to the bit shift 1<<I
01044     // Division of an integer by 2^I is equal to the bit shift >>I
01045     return ((_currentPostion & 1<<I )>>I);
01046   }
01047   int getCurrentPosition() const{
01048     return _currentPostion;
01049   }
01050   /**
01051    * Returns the position number given an element of the power set, i.e.
01052    * an int vector with component values of zero or one.
01053    *
01054    * This is needed for example, if you want to know which parameter is
01055    * active at a certain position in the computational domain from the
01056    * signs of the levelset functions.
01057    */
01058   static int getPositionNumberFromVector( const aol::Vector<int> &Vec ){
01059     int position = 0;
01060     for( int i = 0; i < Vec.size(); i++ )
01061       position += Vec[i]<<i;
01062     return position;
01063   }
01064   /**
01065    * Returns the position number given a global index and a set of
01066    * levelset functions.
01067    *
01068    * This is needed for example, if you want to know which parameter is
01069    * active at a certain position in the computational domain specified
01070    * by the global index based on the signs of the levelset functions.
01071    */
01072   template <typename RealType>
01073   static int getPositionNumberFromLevelsetFunctions( const aol::MultiVector<RealType> &LevelsetFunctions,
01074                                                      const int GlobalIndex,
01075                                                      const RealType Isovalue = aol::NumberTrait<RealType>::zero ){
01076     static aol::Vector<int> intVec;
01077     intVec.resize( LevelsetFunctions.numComponents() );
01078     for ( int i = 0; i < LevelsetFunctions.numComponents(); i++ )
01079       intVec[i] = (LevelsetFunctions[i].get(GlobalIndex)>= Isovalue) ? 1 : 0;
01080     return getPositionNumberFromVector(intVec);
01081   }
01082 };
01083 
01084 /**
01085  * \brief A helper class for Chan Vese with multiple levelset functions.
01086  *
01087  * Has to be constructed (or initialized with initValues) for each point where you want to evaluate
01088  * the product of the Heaviside funtions of the levelset functions.
01089  *
01090  * Note: In a proper Chan Vese model evaluateDerivative is not necessary,
01091  * but caching the necessary values for this function causes a noticeably
01092  * performance hit. Terefore all code for this is commented out at the
01093  * moment.
01094  *
01095  * If HeavisideFunctionType::ScalingFunctionType is equal to IdentityFunction,
01096  * this class fits to the classical Chan Vese model.
01097  *
01098  * \author Berkels
01099  */
01100 template <typename ConfiguratorType, typename HeavisideFunctionType>
01101 class HeavisideFunctionProduct {
01102 public:
01103   typedef typename ConfiguratorType::RealType RealType;
01104 protected:
01105   typedef typename HeavisideFunctionType::ScalingFunctionType ScalingFunctionType;
01106   const HeavisideFunctionType &_heavisideFunction;
01107   const int _numLevelsetFunctions;
01108   aol::Vector<RealType> _hOfPhiValues;
01109   //aol::Vector<RealType> _hPrimeOfPhiValues;
01110   vector< aol::DiscreteFunctionDefault<ConfiguratorType>* > _discrLevesetFunctionsVec;
01111 
01112   //! May only be called once! Calling it more than once will lead to memory leaks.
01113   void initDiscrLevesetFunctionsVec ( const typename ConfiguratorType::InitType &Initializer,
01114                                       const aol::MultiVector<RealType> &LevelsetFunctions ) {
01115     _discrLevesetFunctionsVec.resize( _numLevelsetFunctions );
01116     for( int i = 0; i < _numLevelsetFunctions; i++ ){
01117       _discrLevesetFunctionsVec[i] = new aol::DiscreteFunctionDefault<ConfiguratorType> ( Initializer, LevelsetFunctions[i] );
01118     }
01119   }
01120 public:
01121   void initValues ( const typename ConfiguratorType::ElementType &El,
01122                     const int QuadPoint ) {
01123     for( int i = 0; i < _numLevelsetFunctions; i++ ){
01124       const RealType phiValue = _discrLevesetFunctionsVec[i]->evaluateAtQuadPoint( El, QuadPoint);
01125       _hOfPhiValues[i] = _heavisideFunction.evaluate( phiValue );
01126       //_hPrimeOfPhiValues[i] = _heavisideFunction.evaluateDerivative( phiValue );
01127     }
01128   }
01129   HeavisideFunctionProduct( const typename ConfiguratorType::InitType &Initializer,
01130                             const HeavisideFunctionType &HeavisideFunction,
01131                             const aol::MultiVector<RealType> &LevelsetFunctions )
01132   : _heavisideFunction( HeavisideFunction ),
01133     _numLevelsetFunctions( LevelsetFunctions.numComponents() ),
01134     _hOfPhiValues( _numLevelsetFunctions )//,
01135     //_hPrimeOfPhiValues( _numLevelsetFunctions )
01136   {
01137     initDiscrLevesetFunctionsVec ( Initializer, LevelsetFunctions );
01138   }
01139   ~HeavisideFunctionProduct () {
01140     for( unsigned int i = 0; i < _discrLevesetFunctionsVec.size(); i++ ){
01141       delete _discrLevesetFunctionsVec[i];
01142     }
01143   }
01144   //! Note: Calling this constructor on every quad point will be slow!
01145   HeavisideFunctionProduct( const typename ConfiguratorType::InitType &Initializer,
01146                             const HeavisideFunctionType &HeavisideFunction,
01147                             const aol::MultiVector<RealType> &LevelsetFunctions,
01148                             const typename ConfiguratorType::ElementType &El,
01149                             const int QuadPoint )
01150   : _heavisideFunction( HeavisideFunction ),
01151     _numLevelsetFunctions( LevelsetFunctions.numComponents() ),
01152     _hOfPhiValues( _numLevelsetFunctions )//,
01153     //_hPrimeOfPhiValues( _numLevelsetFunctions )
01154   {
01155     initDiscrLevesetFunctionsVec ( Initializer, LevelsetFunctions );
01156     initValues ( El, QuadPoint );
01157   }
01158   //! Note: Calling this constructor on every quad point will be slow!
01159   HeavisideFunctionProduct( const typename ConfiguratorType::InitType &Initializer,
01160                             const HeavisideFunctionType &HeavisideFunction,
01161                             const aol::MultiVector<RealType> &LevelsetFunctions,
01162                             const typename ConfiguratorType::ElementType &El,
01163                             const typename ConfiguratorType::VecType& RefCoord )
01164   : _heavisideFunction( HeavisideFunction ),
01165     _numLevelsetFunctions( LevelsetFunctions.numComponents() ),
01166     _hOfPhiValues( _numLevelsetFunctions )//,
01167     //_hPrimeOfPhiValues( _numLevelsetFunctions )
01168   {
01169     for( int i = 0; i < _numLevelsetFunctions; i++ ){
01170       aol::DiscreteFunctionDefault<ConfiguratorType> discrFunc( Initializer, LevelsetFunctions[i] );
01171       const RealType phiValue = discrFunc.evaluate( El, RefCoord );
01172       _hOfPhiValues[i] = _heavisideFunction.evaluate( phiValue );
01173       //_hPrimeOfPhiValues[i] = _heavisideFunction.evaluateDerivative( phiValue );
01174     }
01175   }
01176   /**
01177    * Helper function to decide whether to use "x" or "1-x" without using an if statement.
01178    * f(x,1) = x
01179    * f(x,0) = 1-x
01180    */
01181   inline RealType f( const RealType x, const int b ) const { 
01182     return ( 2*(b-0.5f)*(x - 0.5f) + 0.5f );
01183   }
01184   /**
01185    * Helper function to decide whether to use "x" or "-x" without using an if statement.
01186    * g(x,1) = x
01187    * g(x,0) = -x
01188    */
01189   inline RealType g( const RealType x, const int b ) const { 
01190     return ( 2*(b-0.5f)*x );
01191   }
01192   RealType evaluate( const aol::PowerSetIterator &Iterator ) const {
01193     RealType value = aol::ZOTrait<RealType>::one;
01194     for( int i = 0; i < _numLevelsetFunctions; i++ ){
01195       value *= ScalingFunctionType::evaluate ( f( _hOfPhiValues[i], Iterator.getComponent(i) ) );
01196     }
01197     return value;
01198   }
01199   RealType evaluateSkippingOneComponent( const int ComponentToSkip, const aol::PowerSetIterator &Iterator ) const {
01200     RealType value = aol::ZOTrait<RealType>::one;
01201 
01202     for( int i = 0; i < ComponentToSkip; i++ ){
01203       value *= ScalingFunctionType::evaluate ( f( _hOfPhiValues[i], Iterator.getComponent(i) ) );
01204     }
01205     for( int i = ComponentToSkip+1; i < _numLevelsetFunctions; i++ ){
01206       value *= ScalingFunctionType::evaluate ( f( _hOfPhiValues[i], Iterator.getComponent(i) ) );
01207     }
01208     return value;
01209   }
01210   /**
01211    * ATTENTION: If ScalingFunctionType::evaluateDerivative doesn't always return one,
01212    * evaluateSkippingOneComponentWithSign doesn't exactly do what the name suggests.
01213    */
01214   RealType evaluateSkippingOneComponentWithSign( const int ComponentToSkip, const aol::PowerSetIterator &Iterator ) const {
01215     const RealType sign = 2*(Iterator.getComponent(ComponentToSkip)-0.5f);
01216     return (sign * ScalingFunctionType::evaluateDerivative ( f( _hOfPhiValues[ComponentToSkip], Iterator.getComponent(ComponentToSkip) ) ) * evaluateSkippingOneComponent( ComponentToSkip, Iterator ));
01217   }
01218   RealType evaluateSkippingTwoComponents( const int ComponentToSkipOne, const int ComponentToSkipTwo, const aol::PowerSetIterator &Iterator ) const {
01219     RealType value = aol::ZOTrait<RealType>::one;
01220 
01221     for( int i = 0; i < _numLevelsetFunctions; i++ ){
01222       if ( ( i != ComponentToSkipOne ) && ( i != ComponentToSkipTwo ) )
01223         value *= ScalingFunctionType::evaluate ( f( _hOfPhiValues[i], Iterator.getComponent(i) ) );
01224     }
01225     return value;
01226   }
01227   RealType evaluateMixedSecondDerivativeForCEMModel( const int ComponentOne, const int ComponentTwo, const aol::PowerSetIterator &Iterator ) const {
01228     const RealType sign1 = 2*(Iterator.getComponent(ComponentOne)-0.5f);
01229     const RealType sign2 = 2*(Iterator.getComponent(ComponentTwo)-0.5f);
01230     return ( sign1 * ScalingFunctionType::evaluateDerivative ( f( _hOfPhiValues[ComponentOne], Iterator.getComponent(ComponentOne) ) )
01231              * sign2 * ScalingFunctionType::evaluateDerivative ( f( _hOfPhiValues[ComponentTwo], Iterator.getComponent(ComponentTwo) ) )
01232              * evaluateSkippingTwoComponents( ComponentOne, ComponentTwo, Iterator ));
01233   }
01234   /*
01235   RealType evaluateDerivative( const int ComponentToDerive, const aol::PowerSetIterator &Iterator ) const {
01236     const RealType value = g( _hPrimeOfPhiValues[ComponentToDerive], Iterator.getComponent(ComponentToDerive) );
01237     return (value * evaluateSkippingOneComponent( ComponentToDerive, Iterator ));
01238   }
01239   */
01240 };
01241 
01242 /**
01243  * \brief Given n levelset functions and 2^n ops, this class calculates the sum 
01244  *        over these ops, where the integrands of the ops are multiplied with
01245  *        the proper combination of heaviside functions of the levelset
01246  *        functions as needed in the Chan Vese model for more than two segments.
01247  *        The ops have to be derived from FENonlinIntegrationVectorInterface.
01248  *
01249  * \author Berkels
01250  */
01251 template <typename ConfiguratorType, typename HeavisideFunctionType, typename TargetOpType>
01252 class HeavisideProductFENonlinIntegrationVectorInterface
01253 : public aol::FENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideProductFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>, TargetOpType::NumOfComponents > {
01254 public:
01255   typedef typename ConfiguratorType::RealType RealType;
01256 protected:
01257   const aol::MultiVector<RealType> &_levelsetFunctions;
01258   const HeavisideFunctionType &_heavisideFunction;
01259   const std::vector<const TargetOpType*> _targetOpVec;
01260 public:
01261   HeavisideProductFENonlinIntegrationVectorInterface( const typename ConfiguratorType::InitType &Initializer,
01262                                                         const aol::MultiVector<RealType> &LevelsetFunctions,
01263                                                         const HeavisideFunctionType &HeavisideFunction,
01264                                                         const std::vector<const TargetOpType*> &TargetOpVec )
01265     : aol::FENonlinIntegrationVectorInterface< ConfiguratorType,
01266                                                  HeavisideProductFENonlinIntegrationVectorInterface<ConfiguratorType, HeavisideFunctionType, TargetOpType>,
01267                                                  TargetOpType::NumOfComponents > ( Initializer ),
01268       _levelsetFunctions( LevelsetFunctions ),
01269       _heavisideFunction( HeavisideFunction ),
01270       _targetOpVec( TargetOpVec ){
01271     // Given n levelset functions, we have 2^n possible segments and also need that many target ops.
01272     if ( 1<<_levelsetFunctions.numComponents() != static_cast<int>(_targetOpVec.size()) )
01273       throw aol::Exception( "pow( 2, _levelsetFunctions.numComponents()) != _targetOpVec.length()", __FILE__, __LINE__ );
01274   }
01275 
01276   RealType evaluateIntegrand( const aol::auto_container<TargetOpType::NumOfComponents,aol::DiscreteFunctionDefault<ConfiguratorType> > &discrFuncs,
01277                               const typename ConfiguratorType::ElementType &El,
01278                               int QuadPoint, const typename ConfiguratorType::VecType &RefCoord ) const {
01279     // Construct HeavisideFunctionProduct at the point, where we want to evaluate it.
01280     HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType>
01281       heavisideFunctionProduct( this->_initializer, _heavisideFunction, _levelsetFunctions, El, QuadPoint );
01282 
01283     RealType integrand = aol::ZOTrait<RealType>::zero;
01284     // Do the summation over all combinations of Levelset function signs.
01285     aol::PowerSetIterator iterator(_levelsetFunctions.numComponents());
01286     do
01287     {
01288       integrand += ( _targetOpVec[iterator.getCurrentPosition()]->evaluateIntegrand( discrFuncs, El, QuadPoint, RefCoord )
01289                      * heavisideFunctionProduct.evaluate( iterator ) );
01290       iterator.increment();
01291     } while ( !iterator.end() );
01292 
01293     return integrand;
01294   }
01295 };
01296 
01297 /**
01298  * Checks if the MultiVector is of the correct structure to be an
01299  * argument of one of the multi levelset Chan Vese operators.
01300  * Throws an expection if the check fails.
01301  *
01302  * \author Berkels
01303  */
01304 template <typename RealType>
01305 void checkMultiChanVeseStructure( const aol::MultiVector<RealType> &LevelsetfunctionsAndParameters, const int NumberOfLevelsetFunctions, const int NumberOfAdditionalVectors = 0 ){
01306   if ( LevelsetfunctionsAndParameters.numComponents() != NumberOfLevelsetFunctions + (1<<NumberOfLevelsetFunctions) + NumberOfAdditionalVectors )
01307     throw aol::Exception( "LevelsetfunctionsAndParameters.numComponents() != NumberOfLevelsetFunctions + 1<<numberOfLevelsetFunctions", __FILE__, __LINE__ );
01308 }
01309 
01310 /**
01311  * General Interface for multi levelset Chan Vese type energies.
01312  * Only evaluateIndicator has to be provided. For optimizations
01313  * cacheIndicatorParameters can be overloaded.
01314  *
01315  * Apply expects MArg to contain the n levelset functions in the
01316  * first n components and the indicator paramters in the next
01317  * 2^n components.
01318  *
01319  * If MArg contains more than the above mentioned n+2^n components,
01320  * _numOfAdditionalVectorComponentsInArgument has to be set
01321  * accordingly.
01322  *
01323  * \author Berkels
01324  */
01325 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ParameterDimension, typename Imp>
01326 class ChanVeseEnergyInterface
01327   : public Op<aol::MultiVector<typename ConfiguratorType::RealType>, aol::Scalar<typename ConfiguratorType::RealType> > {
01328 public:
01329   typedef typename ConfiguratorType::RealType RealType;
01330 protected:
01331   mutable ConfiguratorType _config;
01332   const typename ConfiguratorType::InitType &_initializer;
01333   const HeavisideFunctionType &_heavisideFunction;
01334   int _numOfAdditionalVectorComponentsInArgument;
01335 public:
01336   ChanVeseEnergyInterface ( const typename ConfiguratorType::InitType &Initializer,
01337                             const HeavisideFunctionType &HeavisideFunction ) :
01338     _config ( Initializer ),
01339     _initializer ( Initializer ),
01340     _heavisideFunction( HeavisideFunction ),
01341     _numOfAdditionalVectorComponentsInArgument( 0 )
01342   {}
01343 
01344   virtual ~ChanVeseEnergyInterface() {}
01345 
01346   virtual void cacheIndicatorParameters( const aol::MultiVector<RealType> &/*Parameters*/ ) const {}
01347 
01348   void setNumOfAdditionalVectorComponentsInArgument ( const int NumOfAdditionalVectorComponentsInArgument ){
01349     _numOfAdditionalVectorComponentsInArgument = NumOfAdditionalVectorComponentsInArgument;
01350   }
01351   
01352   void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::Scalar<RealType> &Dest ) const {
01353 
01354     checkMultiChanVeseStructure<RealType>(MArg, NumberOfLevelsetFunctions, _numOfAdditionalVectorComponentsInArgument);
01355 
01356     // Group the levelset functions contained in MArg into a new MultiVector.
01357     aol::MultiVector<RealType> levelsetFunctions( 0, 0 );
01358     for( int i = 0; i < NumberOfLevelsetFunctions; i++ )
01359       levelsetFunctions.appendReference( MArg[i] );
01360 
01361     // Group the parameter vectors contained in MArg into a new MultiVector.
01362     aol::MultiVector<RealType> parameters( 0, 0 );
01363     for( int i = 0; i < (1<<NumberOfLevelsetFunctions); i++ )
01364       parameters.appendReference( MArg[i+NumberOfLevelsetFunctions] );
01365 
01366     // Since the indicator parameters are fixed during this calculation
01367     // it may be possible to speed up evaluateIndicatorDerivative, by
01368     // pre-calculating all dependencies on indicator parameters here.
01369     cacheIndicatorParameters( parameters );
01370 
01371     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01372 
01373     RealType res = aol::ZOTrait<RealType>::zero;
01374 
01375     HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType>
01376       heavisideFunctionProduct( _initializer, _heavisideFunction, levelsetFunctions );
01377 
01378 #ifdef _OPENMP
01379     std::vector<IteratorType> itVec(0);
01380     for ( IteratorType it = _config.begin(); it != _config.end(); ++it )
01381       itVec.push_back( it );
01382 
01383     const int numberOfElements = static_cast<int>(itVec.size());
01384 #ifdef _OPENMP
01385 #pragma omp parallel for reduction(+: res)
01386 #endif
01387     for ( int elementNumber = 0; elementNumber < numberOfElements; elementNumber++ ) {
01388       IteratorType it = itVec[elementNumber];
01389 #else
01390     for ( IteratorType it = _config.begin(); it != _config.end(); ++it ) {
01391 #endif
01392       typedef typename ConfiguratorType::QuadType QType;
01393       RealType a = aol::ZOTrait<RealType>::zero;
01394       for ( int q = 0; q < QType::numQuadPoints; q++ ) {
01395 
01396         // Init HeavisideFunctionProduct at the point, where we want to evaluate it.
01397         heavisideFunctionProduct.initValues ( *it, q );
01398 
01399         // Loop over all combinations of Levelset function signs, each represents one of the indicator parameters.
01400         aol::PowerSetIterator iterator( levelsetFunctions.numComponents() );
01401         const RealType weight = _config.getBaseFunctionSet ( *it ).getWeight ( q );
01402         const typename ConfiguratorType::VecType &refCoord = _config.getBaseFunctionSet ( *it ).getRefCoord ( q );
01403         do
01404         {
01405           const int i = iterator.getCurrentPosition();
01406           a += ( this->asImp().evaluateIndicator( MArg[NumberOfLevelsetFunctions+i], i, *it, q, refCoord )
01407                  * heavisideFunctionProduct.evaluate( iterator ) )
01408                  * weight;
01409           iterator.increment();
01410         } while ( !iterator.end() );
01411 
01412       }
01413 
01414       a *= _config.vol ( *it );
01415 
01416       res += a;
01417     }
01418     Dest += res;
01419   }
01420 
01421   //! interface function, has to be provided in derived classes.
01422   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01423                                const int IndicatorParameterIndex,
01424                                const typename ConfiguratorType::ElementType &El,
01425                                int QuadPoint,
01426                                const typename ConfiguratorType::VecType &RefCoord ) const {
01427     throw aol::Exception ( "called the interface function", __FILE__, __LINE__ );
01428     return this->asImp().evaluateIndicator ( IndicatorParameter, IndicatorParameterIndex, El, QuadPoint, RefCoord );
01429   }
01430 
01431 protected:
01432   // barton-nackman
01433   inline Imp& asImp() {
01434     return static_cast<Imp&> ( *this );
01435   }
01436   inline const Imp& asImp() const {
01437     return static_cast<const Imp&> ( *this );
01438   }
01439 
01440 };
01441 
01442 
01443 /**
01444  * General Interface for the derivate of multi levelset Chan Vese type
01445  * energies with respect to the indicator parameters. Only
01446  * evaluateIndicatorDerivative has to be provided. For optimizations
01447  * cacheIndicatorParameters can be overloaded.
01448  *
01449  * \author Berkels
01450  */
01451 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ParameterDimension, typename Imp>
01452 class ChanVeseEnergyParameterDerivativeInterface
01453       : public Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
01454 public:
01455   typedef typename ConfiguratorType::RealType RealType;
01456 protected:
01457   mutable ConfiguratorType _config;
01458   const typename ConfiguratorType::InitType &_initializer;
01459   const aol::MultiVector<RealType> &_levelsetFunctions;
01460   const HeavisideFunctionType &_heavisideFunction;
01461 public:
01462   ChanVeseEnergyParameterDerivativeInterface ( const typename ConfiguratorType::InitType &Initializer,
01463                                                const aol::MultiVector<RealType> &LevelsetFunctions,
01464                                                const HeavisideFunctionType &HeavisideFunction ) :
01465     _config ( Initializer ),
01466     _initializer ( Initializer ),
01467     _levelsetFunctions( LevelsetFunctions ),
01468     _heavisideFunction( HeavisideFunction )
01469   {}
01470 
01471   virtual ~ChanVeseEnergyParameterDerivativeInterface() {}
01472 
01473   virtual void cacheIndicatorParameters( const aol::MultiVector<RealType> &/*Parameters*/ ) const {}
01474 
01475   void applyAdd ( const aol::MultiVector<RealType> &Arg, aol::MultiVector<RealType> &Dest ) const {
01476 
01477     // Since the indicator parameters are fixed during this calculation
01478     // it may be possible to speed up evaluateIndicatorDerivative, by
01479     // pre-calculating all dependencies on indicator parameters here.
01480     cacheIndicatorParameters( Arg );
01481 
01482     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01483 
01484     const int NumberOfSegments = 1<<NumberOfLevelsetFunctions;
01485 
01486     aol::Mat<NumberOfSegments, ParameterDimension, RealType> res;
01487     aol::Mat<NumberOfSegments, ParameterDimension, RealType> a;
01488     aol::Vec<ParameterDimension, RealType> tempVec;
01489 
01490     HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType>
01491       heavisideFunctionProduct( _initializer, _heavisideFunction, _levelsetFunctions );
01492 
01493     for ( IteratorType it = _config.begin(); it != _config.end(); ++it ) {
01494       typedef typename ConfiguratorType::QuadType QType;
01495       a.setZero();
01496       for ( int q = 0; q < QType::numQuadPoints; q++ ) {
01497 
01498         // Init HeavisideFunctionProduct at the point, where we want to evaluate it.
01499         heavisideFunctionProduct.initValues( *it, q );
01500 
01501         // Loop over all combinations of Levelset function signs, each represents one of the indicator parameters.
01502         aol::PowerSetIterator iterator(_levelsetFunctions.numComponents());
01503         const RealType weight = _config.getBaseFunctionSet ( *it ).getWeight ( q ); 
01504         do
01505         {
01506           const int i = iterator.getCurrentPosition();
01507           this->asImp().evaluateIndicatorDerivative( Arg[i], i, *it, q, _config.getBaseFunctionSet ( *it ).getRefCoord ( q ), tempVec );
01508           tempVec *= heavisideFunctionProduct.evaluate( iterator ) * weight;
01509           a[i] += tempVec;
01510           iterator.increment();
01511         } while ( !iterator.end() );
01512 
01513       }
01514 
01515       a *= _config.vol ( *it );
01516 
01517       res += a;
01518     }
01519     for ( int i = 0; i < NumberOfSegments; i++ ) {
01520       for ( int j = 0; j < ParameterDimension; j++ ) {
01521         Dest[i][j] += res[i][j];
01522       }
01523     }
01524   }
01525 
01526   //! interface function, has to be provided in derived classes.
01527   void evaluateIndicatorDerivative ( const aol::Vector<RealType> &IndicatorParameter,
01528                                      const int IndicatorParameterIndex,
01529                                      const typename ConfiguratorType::ElementType &El,
01530                                      int QuadPoint,
01531                                      const typename ConfiguratorType::VecType &RefCoord,
01532                                      aol::Vec<ParameterDimension, RealType> &Derivative) const {
01533     throw aol::Exception ( "called the interface function", __FILE__, __LINE__ );
01534     this->asImp().evaluateIndicatorDerivative ( IndicatorParameter, IndicatorParameterIndex, El, QuadPoint, RefCoord, Derivative );
01535   }
01536 
01537 protected:
01538   // barton-nackman
01539   inline Imp& asImp() {
01540     return static_cast<Imp&> ( *this );
01541   }
01542   inline const Imp& asImp() const {
01543     return static_cast<const Imp&> ( *this );
01544   }
01545 
01546 };
01547 
01548 
01549 /**
01550  * General Interface for the derivate of multi levelset Chan Vese type
01551  * energies with respect to the levelset functions. Only
01552  * evaluateIndicator has to be provided. For optimizations
01553  * cacheIndicatorParameters can be overloaded.
01554  *
01555  * \author Berkels
01556  */
01557 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, typename Imp>
01558 class ChanVeseEnergyLevelsetDerivativeInterface
01559       : public Op<aol::MultiVector<typename ConfiguratorType::RealType> > {
01560 public:
01561   typedef typename ConfiguratorType::RealType RealType;
01562 protected:
01563   mutable ConfiguratorType _config;
01564   const typename ConfiguratorType::InitType &_initializer;
01565   const aol::MultiVector<RealType> &_parameters;
01566   const HeavisideFunctionType &_heavisideFunction;
01567 public:
01568   ChanVeseEnergyLevelsetDerivativeInterface ( const typename ConfiguratorType::InitType &Initializer,
01569                                               const aol::MultiVector<RealType> &Parameters,
01570                                               const HeavisideFunctionType &HeavisideFunction ) :
01571     _config ( Initializer ),
01572     _initializer ( Initializer ),
01573     _parameters( Parameters ),
01574     _heavisideFunction( HeavisideFunction )
01575   {}
01576 
01577   virtual ~ChanVeseEnergyLevelsetDerivativeInterface() {}
01578 
01579   void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::MultiVector<RealType> &MDest ) const {
01580 
01581     typedef typename ConfiguratorType::ElementIteratorType IteratorType;
01582     aol::Vec<NumberOfLevelsetFunctions, RealType> a;
01583     aol::FullMatrix<RealType> indicatorCache(_config.maxNumQuadPoints(), _parameters.numComponents() );
01584 
01585     std::vector<HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType>*> heavisideFunctionProductCache(_config.maxNumQuadPoints());
01586     for ( int q = 0; q < _config.maxNumQuadPoints(); q++ ) {
01587       heavisideFunctionProductCache[q] = new HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType>( _initializer, _heavisideFunction, MArg );
01588     }
01589 
01590 #ifdef _OPENMP
01591     std::vector<IteratorType> itVec(0);
01592     for ( IteratorType it = _config.begin(); it != _config.end(); ++it )
01593       itVec.push_back( it );
01594 
01595     const int numberOfElements = static_cast<int>(itVec.size());
01596 #ifdef _OPENMP
01597 #pragma omp parallel for
01598 #endif
01599     for ( int elementNumber = 0; elementNumber < numberOfElements; elementNumber++ ) {
01600       IteratorType it = itVec[elementNumber];
01601 #else
01602     for ( IteratorType it = _config.begin(); it != _config.end(); ++it ) {
01603 #endif
01604       typedef typename ConfiguratorType::QuadType QType;
01605       a.setZero();
01606       const int numLocalDofs = _config.getNumLocalDofs ( *it );
01607 
01608       const typename ConfiguratorType::BaseFuncSetType &bfs = _config.getBaseFunctionSet ( *it );
01609       const int numQuadPoints = bfs.numQuadPoints( );
01610 
01611       // Cache the indicator values at all quad points. We do this outside the "dof" loop
01612       // since the values don't depend on dof.
01613       // The extra scope deletes iterator after the iteration.
01614       {
01615         aol::PowerSetIterator iterator( NumberOfLevelsetFunctions );
01616         do
01617         {
01618           const int i = iterator.getCurrentPosition();
01619           for ( int q = 0; q < numQuadPoints; q++ ) {
01620             indicatorCache.set(q,i,this->asImp().evaluateIndicator( _parameters[i], i, *it, q, _config.getBaseFunctionSet ( *it ).getRefCoord ( q ) ));
01621           }
01622           iterator.increment();
01623         } while ( !iterator.end() );
01624       }
01625 
01626       // Init HeavisideFunctionProduct at each quad point to cache it's values there.
01627       // We do this outside the "dof" loop since the values don't depend on dof.
01628       for ( int q = 0; q < numQuadPoints; q++ ) {
01629         heavisideFunctionProductCache[q]->initValues ( *it, q );
01630       }
01631 
01632       for ( int dof = 0; dof < numLocalDofs; dof++ ) {
01633 
01634         for ( int q = 0; q < numQuadPoints; q++ ) {
01635 
01636           // Loop over all combinations of Levelset function signs, each represents one of the indicator parameters.
01637           aol::PowerSetIterator iterator( NumberOfLevelsetFunctions );
01638           const RealType weight = bfs.getWeight ( q ); 
01639           const RealType basefunctionValue = bfs.evaluate ( dof, q );
01640           do
01641           {
01642             const int i = iterator.getCurrentPosition();
01643             const RealType temp = indicatorCache.get(q,i)*weight*basefunctionValue;
01644             for( int j = 0; j < NumberOfLevelsetFunctions; j++ )
01645               a[j] += temp * heavisideFunctionProductCache[q]->evaluateSkippingOneComponentWithSign( j, iterator );
01646             // Perhaps evaluateSkippingOneComponentWithSign should be cached somehow?
01647             iterator.increment();
01648           } while ( !iterator.end() );
01649 
01650         }
01651 
01652         a *= _config.vol ( *it );
01653 
01654         const int globalIndex = _config.localToGlobal ( *it, dof );
01655         for( int j = 0; j < NumberOfLevelsetFunctions; j++ )
01656         {
01657 #ifdef _OPENMP
01658 #pragma omp atomic
01659 #endif
01660           MDest[j][globalIndex] += a[j];
01661         }
01662       }
01663     }
01664     for ( int q = 0; q < _config.maxNumQuadPoints(); q++ ) {
01665       delete heavisideFunctionProductCache[q];
01666     }
01667   }
01668 
01669   //! interface function, has to be provided in derived classes.
01670   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01671                                const int IndicatorParameterIndex,
01672                                const typename ConfiguratorType::ElementType &El,
01673                                int QuadPoint,
01674                                const typename ConfiguratorType::VecType &RefCoord ) const {
01675     throw aol::Exception ( "called the interface function", __FILE__, __LINE__ );
01676     return this->asImp().evaluateIndicator ( IndicatorParameter, IndicatorParameterIndex, El, QuadPoint, RefCoord );
01677   }
01678 
01679 protected:
01680   // barton-nackman
01681   inline Imp& asImp() {
01682     return static_cast<Imp&> ( *this );
01683   }
01684   inline const Imp& asImp() const {
01685     return static_cast<const Imp&> ( *this );
01686   }
01687 
01688 };
01689 
01690 
01691 /**
01692  * General Interface for the second derivate of multi levelset
01693  * quadratic Chan Esedoglu Nikolova type energies with respect
01694  * to one of the levelset functions. Only evaluateIndicator has
01695  * to be provided. For optimizations cacheIndicatorParameters
01696  * can be overloaded.
01697  *
01698  * Assumes that LevelsetFunctions and Parameters don't change
01699  * after this class is created, because cacheIndicatorParameters
01700  * is only called once in the constructor and nowhere else.
01701  *
01702  * Only meant to be used with
01703  * HeavisideFuncType == aol::IdentityFunction<typename ConfiguratorType::RealType>
01704  * Probably the template argument should be removed completely.
01705  *
01706  * \author Berkels
01707  */
01708 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, typename Imp>
01709 class QuadraticCEMEnergySecondLevelsetDerivativeInterface
01710   : public aol::FELinScalarWeightedMassInterface<ConfiguratorType, QuadraticCEMEnergySecondLevelsetDerivativeInterface<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, Imp> > {
01711 public:
01712   typedef typename ConfiguratorType::RealType RealType;
01713 protected:
01714   const aol::MultiVector<RealType> &_levelsetFunctions;
01715   const aol::MultiVector<RealType> &_parameters;
01716   const HeavisideFunctionType &_heavisideFunction;
01717   mutable HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType> _heavisideFunctionProduct;
01718   const int _levelsetFunctionNumber;
01719 public:
01720   QuadraticCEMEnergySecondLevelsetDerivativeInterface ( const typename ConfiguratorType::InitType &Initializer,
01721                                                         const aol::MultiVector<RealType> &LevelsetFunctions,
01722                                                         const aol::MultiVector<RealType> &Parameters,
01723                                                         const HeavisideFunctionType &HeavisideFunction,
01724                                                         const int LevelsetFunctionNumber,
01725                                                         aol::OperatorType OpType = ONTHEFLY )
01726   : aol::FELinScalarWeightedMassInterface<ConfiguratorType, QuadraticCEMEnergySecondLevelsetDerivativeInterface<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, Imp> > ( Initializer, OpType ),
01727     _levelsetFunctions ( LevelsetFunctions ),
01728     _parameters( Parameters ),
01729     _heavisideFunction( HeavisideFunction ),
01730     _heavisideFunctionProduct( Initializer, HeavisideFunction, LevelsetFunctions ),
01731     _levelsetFunctionNumber ( LevelsetFunctionNumber )
01732   {
01733     if ( HeavisideFunctionType::ScalingFunctionType::evaluate != aol::SquareFunction<RealType>::evaluate )
01734       throw aol::Exception( "QuadraticCEMEnergySecondLevelsetDerivativeInterface is only meant to be used with HeavisideFuncType == aol::IdentityFunction<typename ConfiguratorType::RealType> ", __FILE__, __LINE__ );
01735 
01736     // Since the indicator parameters are assumed to be fixed here
01737     // it may be possible to speed up evaluateIndicatorDerivative, by
01738     // pre-calculating all dependencies on indicator parameters here.
01739     cacheIndicatorParameters( Parameters );
01740   }
01741 
01742   virtual ~QuadraticCEMEnergySecondLevelsetDerivativeInterface() {}
01743 
01744   virtual void cacheIndicatorParameters( const aol::MultiVector<RealType> &/*Parameters*/ ) const {}
01745 
01746   inline RealType getCoeff ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType& RefCoord ) const {
01747 
01748     _heavisideFunctionProduct.initValues ( El, QuadPoint );
01749 
01750     RealType a = aol::ZOTrait<RealType>::zero;
01751     // Loop over all combinations of Levelset function signs, each represents one of the indicator parameters.
01752     aol::PowerSetIterator iterator( NumberOfLevelsetFunctions );
01753     do
01754     {
01755       const int i = iterator.getCurrentPosition();
01756       a += ( this->asImp().evaluateIndicator( _parameters[i], i, El, QuadPoint, RefCoord )
01757         * _heavisideFunctionProduct.evaluateSkippingOneComponent( _levelsetFunctionNumber, iterator ) )
01758         * 2;
01759       iterator.increment();
01760     } while ( !iterator.end() );
01761 
01762     return a;
01763   }
01764 
01765   //! interface function, has to be provided in derived classes.
01766   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01767                                const int IndicatorParameterIndex,
01768                                const typename ConfiguratorType::ElementType &El,
01769                                int QuadPoint,
01770                                const typename ConfiguratorType::VecType &RefCoord ) const {
01771     throw aol::Exception ( "called the interface function", __FILE__, __LINE__ );
01772     return this->asImp().evaluateIndicator ( IndicatorParameter, IndicatorParameterIndex, El, QuadPoint, RefCoord );
01773   }
01774 
01775 protected:
01776   // barton-nackman
01777   inline Imp& asImp() {
01778     return static_cast<Imp&> ( *this );
01779   }
01780   inline const Imp& asImp() const {
01781     return static_cast<const Imp&> ( *this );
01782   }
01783 
01784 };
01785 
01786 /**
01787  * General Interface for the mixed second derivates of multi levelset
01788  * quadratic Chan Esedoglu Nikolova type energies with respect
01789  * to two different of the levelset functions. Only evaluateIndicator has
01790  * to be provided. For optimizations cacheIndicatorParameters
01791  * can be overloaded.
01792  *
01793  * Assumes that LevelsetFunctions and Parameters don't change
01794  * after this class is created, because cacheIndicatorParameters
01795  * is only called once in the constructor and nowhere else.
01796  *
01797  * Only meant to be used with
01798  * HeavisideFuncType == aol::IdentityFunction<typename ConfiguratorType::RealType>
01799  * Probably the template argument should be removed completely.
01800  *
01801  * \author Berkels
01802  */
01803 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, typename Imp>
01804 class QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface
01805   : public aol::FELinScalarWeightedMassInterface<ConfiguratorType, QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, Imp> > {
01806 public:
01807   typedef typename ConfiguratorType::RealType RealType;
01808 protected:
01809   const aol::MultiVector<RealType> &_levelsetFunctions;
01810   const aol::MultiVector<RealType> &_parameters;
01811   const HeavisideFunctionType &_heavisideFunction;
01812   mutable HeavisideFunctionProduct<ConfiguratorType, HeavisideFunctionType> _heavisideFunctionProduct;
01813   const int _levelsetFunctionNumberOne, _levelsetFunctionNumberTwo;
01814 public:
01815   QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface ( const typename ConfiguratorType::InitType &Initializer,
01816                                                              const aol::MultiVector<RealType> &LevelsetFunctions,
01817                                                              const aol::MultiVector<RealType> &Parameters,
01818                                                              const HeavisideFunctionType &HeavisideFunction,
01819                                                              const int LevelsetFunctionNumberOne,
01820                                                              const int LevelsetFunctionNumberTwo,
01821                                                              aol::OperatorType OpType = ONTHEFLY )
01822   : aol::FELinScalarWeightedMassInterface<ConfiguratorType, QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, Imp> > ( Initializer, OpType ),
01823     _levelsetFunctions ( LevelsetFunctions ),
01824     _parameters( Parameters ),
01825     _heavisideFunction( HeavisideFunction ),
01826     _heavisideFunctionProduct( Initializer, HeavisideFunction, LevelsetFunctions ),
01827     _levelsetFunctionNumberOne ( LevelsetFunctionNumberOne ),
01828     _levelsetFunctionNumberTwo ( LevelsetFunctionNumberTwo )
01829   {
01830     if ( HeavisideFunctionType::ScalingFunctionType::evaluate != aol::SquareFunction<RealType>::evaluate )
01831       throw aol::Exception( "QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface is only meant to be used with HeavisideFuncType == aol::IdentityFunction<typename ConfiguratorType::RealType> ", __FILE__, __LINE__ );
01832 
01833     // Since the indicator parameters are assumed to be fixed here
01834     // it may be possible to speed up evaluateIndicatorDerivative, by
01835     // pre-calculating all dependencies on indicator parameters here.
01836     cacheIndicatorParameters( Parameters );
01837   }
01838 
01839   virtual ~QuadraticCEMEnergyMixedSecondLevelsetDerivativeInterface() {}
01840 
01841   virtual void cacheIndicatorParameters( const aol::MultiVector<RealType> &/*Parameters*/ ) const {}
01842 
01843   inline RealType getCoeff ( const typename ConfiguratorType::ElementType &El, int QuadPoint, const typename ConfiguratorType::DomVecType& RefCoord ) const {
01844 
01845     _heavisideFunctionProduct.initValues ( El, QuadPoint );
01846 
01847     RealType a = aol::ZOTrait<RealType>::zero;
01848     // Loop over all combinations of Levelset function signs, each represents one of the indicator parameters.
01849     aol::PowerSetIterator iterator( NumberOfLevelsetFunctions );
01850     do
01851     {
01852       const int i = iterator.getCurrentPosition();
01853       a += ( this->asImp().evaluateIndicator( _parameters[i], i, El, QuadPoint, RefCoord )
01854         * _heavisideFunctionProduct.evaluateMixedSecondDerivativeForCEMModel( _levelsetFunctionNumberOne, _levelsetFunctionNumberTwo, iterator ) );
01855       iterator.increment();
01856     } while ( !iterator.end() );
01857 
01858     return a;
01859   }
01860 
01861   //! interface function, has to be provided in derived classes.
01862   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01863                                const int IndicatorParameterIndex,
01864                                const typename ConfiguratorType::ElementType &El,
01865                                int QuadPoint,
01866                                const typename ConfiguratorType::VecType &RefCoord ) const {
01867     throw aol::Exception ( "called the interface function", __FILE__, __LINE__ );
01868     return this->asImp().evaluateIndicator ( IndicatorParameter, IndicatorParameterIndex, El, QuadPoint, RefCoord );
01869   }
01870 
01871 protected:
01872   // barton-nackman
01873   inline Imp& asImp() {
01874     return static_cast<Imp&> ( *this );
01875   }
01876   inline const Imp& asImp() const {
01877     return static_cast<const Imp&> ( *this );
01878   }
01879 
01880 };
01881 
01882 /**
01883  * Given a Chan Vese energy op and a parameter vector, this class
01884  * builds an op, which can be applied on MultiVectors of levelset
01885  * functions and calculates the energy for the levelset functions
01886  * with the given parameter vector.
01887  *
01888  * \author Berkels
01889  */
01890 template <typename RealType>
01891 class ChanVeseEnergyLevelsetPart : public aol::Op<aol::MultiVector<RealType>, aol::Scalar<RealType> > {
01892 protected:
01893   const aol::Op<aol::MultiVector<RealType>, aol::Scalar<RealType> > &_chanVeseEnergyOp;
01894   const aol::MultiVector<RealType> &_parameters;
01895 public:
01896   ChanVeseEnergyLevelsetPart ( const aol::Op<aol::MultiVector<RealType>, aol::Scalar<RealType> > &ChanVeseEnergyOp,
01897                                const aol::MultiVector<RealType> &Parameters )
01898     : _chanVeseEnergyOp ( ChanVeseEnergyOp ),
01899       _parameters ( Parameters ) {}
01900   virtual void apply ( const aol::MultiVector<RealType> &MArg, aol::Scalar<RealType> &Dest ) const {
01901     aol::MultiVector<RealType> mtemp( MArg, aol::FLAT_COPY );
01902     mtemp.appendReference( _parameters );
01903     _chanVeseEnergyOp.apply ( mtemp, Dest );
01904   }
01905   virtual void applyAdd ( const aol::MultiVector<RealType> &MArg, aol::Scalar<RealType> &Dest ) const {
01906     aol::MultiVector<RealType> mtemp( MArg, aol::FLAT_COPY );
01907     mtemp.appendReference( _parameters );
01908     _chanVeseEnergyOp.applyAdd ( mtemp, Dest );
01909   }
01910 };
01911 
01912 /**
01913  * Calulates the classical Chan Vese energy for multi levelset gray value segmentation
01914  *
01915  * \author Berkels
01916  */
01917 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions>
01918 class ClassicalChanVeseEnergyMulti
01919 : public aol::ChanVeseEnergyInterface< ConfiguratorType,
01920                                        HeavisideFunctionType,
01921                                        NumberOfLevelsetFunctions,
01922                                        1,
01923                                        ClassicalChanVeseEnergyMulti<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> > {
01924 public:
01925   typedef typename ConfiguratorType::RealType RealType;
01926 protected:
01927   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrU0;
01928 public:
01929   ClassicalChanVeseEnergyMulti( const typename ConfiguratorType::InitType &Initializer,
01930                                 const HeavisideFunctionType &HeavisideFunction,
01931                                 const aol::Vector<RealType> &U0 )
01932   : aol::ChanVeseEnergyInterface< ConfiguratorType,
01933                                   HeavisideFunctionType,
01934                                   NumberOfLevelsetFunctions,
01935                                   1,
01936                                   ClassicalChanVeseEnergyMulti<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> >
01937                                 ( Initializer, HeavisideFunction ),
01938       _discrU0( Initializer, U0 )
01939   {
01940   }
01941   ~ClassicalChanVeseEnergyMulti () {}
01942 
01943   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01944                                const int /*IndicatorParameterIndex*/,
01945                                const typename ConfiguratorType::ElementType &El,
01946                                int QuadPoint,
01947                                const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
01948     return aol::Sqr(_discrU0.evaluateAtQuadPoint( El, QuadPoint) - IndicatorParameter[0] );
01949   }
01950 };
01951 
01952 /**
01953  * Calulates the variation of ClassicalChanVeseEnergyMulti with respect to the levelset functions.
01954  *
01955  * \author Berkels
01956  */
01957 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions>
01958 class ClassicalChanVeseEnergyMultiPhiVariation
01959 : public aol::ChanVeseEnergyLevelsetDerivativeInterface< ConfiguratorType,
01960                                                          HeavisideFunctionType,
01961                                                          NumberOfLevelsetFunctions,
01962                                                          ClassicalChanVeseEnergyMultiPhiVariation<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> > {
01963 public:
01964   typedef typename ConfiguratorType::RealType RealType;
01965 protected:
01966   const aol::DiscreteFunctionDefault<ConfiguratorType> _discrU0;
01967 public:
01968   ClassicalChanVeseEnergyMultiPhiVariation( const typename ConfiguratorType::InitType &Initializer,
01969                                             const HeavisideFunctionType &HeavisideFunction,
01970                                             const aol::Vector<RealType> &U0,
01971                                             const aol::MultiVector<RealType> &Parameters )
01972   : aol::ChanVeseEnergyLevelsetDerivativeInterface< ConfiguratorType,
01973                                                     HeavisideFunctionType,
01974                                                     NumberOfLevelsetFunctions,
01975                                                     ClassicalChanVeseEnergyMultiPhiVariation<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> >
01976                                                    ( Initializer, Parameters, HeavisideFunction ),
01977       _discrU0( Initializer, U0 )
01978   {
01979   }
01980   ~ClassicalChanVeseEnergyMultiPhiVariation () {}
01981 
01982   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
01983                                const int /*IndicatorParameterIndex*/,
01984                                const typename ConfiguratorType::ElementType &El,
01985                                int QuadPoint,
01986                                const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
01987     return aol::Sqr(_discrU0.evaluateAtQuadPoint( El, QuadPoint) - IndicatorParameter[0] );
01988   }
01989 };
01990 
01991 /**
01992  * Calulates the classical Chan Vese energy of vector valued images for multi levelset gray value segmentation
01993  *
01994  * \author Berkels
01995  */
01996 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ImageDimension>
01997 class ClassicalChanVeseVectorEnergyMulti
01998 : public aol::ChanVeseEnergyInterface< ConfiguratorType,
01999                                        HeavisideFunctionType,
02000                                        NumberOfLevelsetFunctions,
02001                                        ImageDimension,
02002                                        ClassicalChanVeseVectorEnergyMulti<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> > {
02003 public:
02004   typedef typename ConfiguratorType::RealType RealType;
02005 protected:
02006   aol::auto_container< ImageDimension, const aol::DiscreteFunctionDefault<ConfiguratorType> > _discrU0;
02007 public:
02008   ClassicalChanVeseVectorEnergyMulti( const typename ConfiguratorType::InitType &Initializer,
02009                                       const HeavisideFunctionType &HeavisideFunction,
02010                                       const aol::MultiVector<RealType> &U0 )
02011   : aol::ChanVeseEnergyInterface< ConfiguratorType,
02012                                   HeavisideFunctionType,
02013                                   NumberOfLevelsetFunctions,
02014                                   ImageDimension,
02015                                   ClassicalChanVeseVectorEnergyMulti<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> >
02016                                 ( Initializer, HeavisideFunction )
02017   {
02018     for ( int c=0; c<ImageDimension; c++ ) {
02019       aol::DiscreteFunctionDefault<ConfiguratorType> temp( Initializer, U0[c] );
02020       _discrU0.set_copy( c, temp );
02021     }
02022   }
02023   ~ClassicalChanVeseVectorEnergyMulti () {}
02024 
02025   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
02026                                const int /*IndicatorParameterIndex*/,
02027                                const typename ConfiguratorType::ElementType &El,
02028                                int QuadPoint,
02029                                const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
02030     RealType indicator = 0.;
02031     for ( int i = 0; i < ImageDimension; i++ )
02032       indicator += aol::Sqr(_discrU0[i].evaluateAtQuadPoint( El, QuadPoint) - IndicatorParameter[i] );
02033     return indicator;
02034   }
02035 };
02036 
02037 /**
02038  * Calulates the variation of ClassicalChanVeseVectorEnergyMulti with respect to the levelset functions.
02039  *
02040  * \author Berkels
02041  */
02042 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ImageDimension>
02043 class ClassicalChanVeseVectorEnergyMultiPhiVariation
02044 : public aol::ChanVeseEnergyLevelsetDerivativeInterface< ConfiguratorType,
02045                                                          HeavisideFunctionType,
02046                                                          NumberOfLevelsetFunctions,
02047                                                          ClassicalChanVeseVectorEnergyMultiPhiVariation<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> > {
02048 public:
02049   typedef typename ConfiguratorType::RealType RealType;
02050 protected:
02051   aol::auto_container< ImageDimension, const aol::DiscreteFunctionDefault<ConfiguratorType> > _discrU0;
02052 public:
02053   ClassicalChanVeseVectorEnergyMultiPhiVariation( const typename ConfiguratorType::InitType &Initializer,
02054                                             const HeavisideFunctionType &HeavisideFunction,
02055                                             const aol::MultiVector<RealType> &U0,
02056                                             const aol::MultiVector<RealType> &Parameters )
02057   : aol::ChanVeseEnergyLevelsetDerivativeInterface< ConfiguratorType,
02058                                                     HeavisideFunctionType,
02059                                                     NumberOfLevelsetFunctions,
02060                                                     ClassicalChanVeseVectorEnergyMultiPhiVariation<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> >
02061                                                    ( Initializer, Parameters, HeavisideFunction )
02062   {
02063     for ( int c=0; c<ImageDimension; c++ ) {
02064       aol::DiscreteFunctionDefault<ConfiguratorType> temp( Initializer, U0[c] );
02065       _discrU0.set_copy( c, temp );
02066     }
02067   }
02068   ~ClassicalChanVeseVectorEnergyMultiPhiVariation () {}
02069 
02070   RealType evaluateIndicator ( const aol::Vector<RealType> &IndicatorParameter,
02071                                const int /*IndicatorParameterIndex*/,
02072                                const typename ConfiguratorType::ElementType &El,
02073                                int QuadPoint,
02074                                const typename ConfiguratorType::VecType &/*RefCoord*/ ) const {
02075     RealType indicator = 0.;
02076     for ( int i = 0; i < ImageDimension; i++ )
02077       indicator += aol::Sqr(_discrU0[i].evaluateAtQuadPoint( El, QuadPoint) - IndicatorParameter[i] );
02078     return indicator;
02079   }
02080 };
02081 
02082 /**
02083  * Calculates the volume of each segment given by the levelset functions.
02084  *
02085  * Argument in apply/applyAdd is ignored.
02086  *
02087  * Is derived from ChanVeseEnergyParameterDerivativeInterface even though 
02088  * this is not a derivative. This can be done nevertheless.
02089  *
02090  * \author Berkels
02091  */
02092 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions>
02093 class MultiLevelsetVolumes
02094 : public aol::ChanVeseEnergyParameterDerivativeInterface< ConfiguratorType,
02095                                                           HeavisideFunctionType,
02096                                                           NumberOfLevelsetFunctions,
02097                                                           1,
02098                                                           MultiLevelsetVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> > {
02099 public:
02100   typedef typename ConfiguratorType::RealType RealType;
02101   MultiLevelsetVolumes( const typename ConfiguratorType::InitType &Initializer,
02102                         const aol::MultiVector<RealType> &LevelsetFunctions,
02103                         const HeavisideFunctionType &HeavisideFunction )
02104   : aol::ChanVeseEnergyParameterDerivativeInterface< ConfiguratorType,
02105                                                      HeavisideFunctionType,
02106                                                      NumberOfLevelsetFunctions,
02107                                                      1,
02108                                                      MultiLevelsetVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions> >
02109                                                    ( Initializer, LevelsetFunctions, HeavisideFunction )
02110   {
02111   }
02112   void evaluateIndicatorDerivative ( const aol::Vector<RealType> &/*IndicatorParameter*/,
02113                                      const int /*IndicatorParameterIndex*/,
02114                                      const typename ConfiguratorType::ElementType &/*El*/,
02115                                      int /*QuadPoint*/,
02116                                      const typename ConfiguratorType::VecType &/*RefCoord*/,
02117                                      aol::Vec<1, RealType> &Derivative ) const {
02118     Derivative[0] = 1.;
02119   }
02120 };
02121 
02122 /**
02123  * Calculates the volume of each segment given by the levelset functions,
02124  * weighted by the gray value of the image U0. If you divide the results
02125  * of this class component-wise by the ones of MultiLevelsetVolumes, you
02126  * get the vector containing the average gray values of the image in the
02127  * different segments.
02128  *
02129  * This also works for vector valued images U0: In this case you have to
02130  * divide the vector valued components of the result by the appropriate
02131  * scalar values (vectors of length one) obtained from MultiLevelsetVolumes. 
02132  *
02133  * Argument in apply/applyAdd is ignored.
02134  *
02135  * Is derived from ChanVeseEnergyParameterDerivativeInterface even though 
02136  * this is not a derivative. This can be done nevertheless.
02137  *
02138  * \author Berkels
02139  */
02140 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ImageDimension>
02141 class MultiLevelsetWeightedVolumes
02142 : public aol::ChanVeseEnergyParameterDerivativeInterface< ConfiguratorType,
02143                                                           HeavisideFunctionType,
02144                                                           NumberOfLevelsetFunctions,
02145                                                           ImageDimension,
02146                                                           MultiLevelsetWeightedVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> > {
02147   aol::auto_container< ImageDimension, const aol::DiscreteFunctionDefault<ConfiguratorType> > _discrU0;
02148 public:
02149   typedef typename ConfiguratorType::RealType RealType;
02150   MultiLevelsetWeightedVolumes( const typename ConfiguratorType::InitType &Initializer,
02151                                 const aol::MultiVector<RealType> &LevelsetFunctions,
02152                                 const HeavisideFunctionType &HeavisideFunction,
02153                                 const aol::MultiVector<RealType> &U0 )
02154   : aol::ChanVeseEnergyParameterDerivativeInterface< ConfiguratorType,
02155                                                      HeavisideFunctionType,
02156                                                      NumberOfLevelsetFunctions,
02157                                                      ImageDimension,
02158                                                      MultiLevelsetWeightedVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension> >
02159                                                    ( Initializer, LevelsetFunctions, HeavisideFunction )
02160   {
02161     for ( int c=0; c<ImageDimension; c++ ) {
02162       aol::DiscreteFunctionDefault<ConfiguratorType> temp( Initializer, U0[c] );
02163       _discrU0.set_copy( c, temp );
02164     }
02165   }
02166   void evaluateIndicatorDerivative ( const aol::Vector<RealType> &/*IndicatorParameter*/,
02167                                      const int /*IndicatorParameterIndex*/,
02168                                      const typename ConfiguratorType::ElementType &El,
02169                                      int QuadPoint,
02170                                      const typename ConfiguratorType::VecType &/*RefCoord*/,
02171                                      aol::Vec<ImageDimension, RealType> &Derivative) const {
02172     for ( int i = 0; i < ImageDimension; i++ )
02173       Derivative[i] = _discrU0[i].evaluateAtQuadPoint( El, QuadPoint);
02174   }
02175 };
02176 
02177 
02178 /**
02179  * Calculates the average mean values (weighted according to LevelsetFunctions
02180  * and HeavisideFunction) in the segments for the classical piecewise constant
02181  * Chan Vese segmentation.
02182  *
02183  * \author Berkels
02184  */
02185 template <typename ConfiguratorType, typename HeavisideFunctionType, int NumberOfLevelsetFunctions, int ImageDimension>
02186 class ClassicalChanVeseMeanValueUpdater {
02187 private:
02188   typedef typename ConfiguratorType::RealType RealType;
02189   const typename ConfiguratorType::InitType &_grid;
02190   const HeavisideFunctionType &_heavisideFunction;
02191   const aol::MultiVector<RealType> &_imageMVec;
02192 public:
02193   ClassicalChanVeseMeanValueUpdater( const typename ConfiguratorType::InitType &Initializer,
02194                                      const HeavisideFunctionType &HeavisideFunction,
02195                                      const aol::MultiVector<RealType> &ImageMVec )
02196     : _grid ( Initializer ),
02197       _heavisideFunction ( HeavisideFunction ),
02198       _imageMVec( ImageMVec )
02199   {
02200   }
02201   void update( const aol::MultiVector<RealType> &LevelsetFunctions, aol::MultiVector<RealType> &MeanGrayValues ) const {
02202     // Calculate the new average values based on the current levelset functions.
02203     aol::MultiLevelsetVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions>
02204       volE( _grid, LevelsetFunctions, _heavisideFunction );
02205 
02206     aol::MultiLevelsetWeightedVolumes<ConfiguratorType, HeavisideFunctionType, NumberOfLevelsetFunctions, ImageDimension>
02207       volEweight( _grid, LevelsetFunctions, _heavisideFunction, _imageMVec );
02208 
02209     aol::MultiVector<RealType> tempMVec( MeanGrayValues.numComponents(), 1 );
02210     aol::MultiVector<RealType> tempGrayValues( MeanGrayValues, aol::STRUCT_COPY );
02211     volE.apply( tempMVec, tempMVec );
02212     // MultiLevelsetWeightedVolumes only uses the structure of the first argument, not the actual content.
02213     volEweight.apply( MeanGrayValues, tempGrayValues );
02214     for( int i = 0; i < MeanGrayValues.numComponents(); i++ ) {
02215       // We can only calculate a meaningful mean gray value for non-empty segments.
02216       if ( appeq ( tempMVec[i][0], aol::ZOTrait<RealType>::zero ) == false ) {
02217         MeanGrayValues[i] = tempGrayValues[i];
02218         MeanGrayValues[i] /= tempMVec[i][0];
02219       }
02220     }
02221   }
02222 };
02223 
02224 } // namespace aol
02225 
02226 #endif

Generated on Fri Sep 9 2011 21:07:38 for QuocMesh by doxygen 1.7.1