();
- QApplication application(argc,argv);
- Viewer3D viewer;
+ Viewer3D<> viewer;
viewer.setWindowTitle("simple Volume Viewer");
viewer.show();
@@ -108,6 +111,6 @@ int main( int argc, char** argv )
viewer << *it;
}
}
- viewer << Viewer3D::updateDisplay;
+ viewer << Viewer3D<>::updateDisplay;
return application.exec();
}
diff --git a/deformations/BinaryPredicates.h b/deformations/BinaryPredicates.h
new file mode 100644
index 0000000..b4fe43b
--- /dev/null
+++ b/deformations/BinaryPredicates.h
@@ -0,0 +1,91 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file BinaryPredicates.h
+ *
+ * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ * @date 2012/02/02
+ *
+ * This files contains several basic classes representing binary predicates
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if defined(BinaryPredicates_RECURSES)
+#error Recursive header files inclusion detected in BinaryPredicates.h
+#else // defined(BinaryPredicates_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define BinaryPredicates_RECURSES
+
+#if !defined BinaryPredicates_h
+/** Prevents repeated inclusion of headers. */
+#define BinaryPredicates_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include "DGtal/base/Common.h"
+//////////////////////////////////////////////////////////////////////////////
+
+namespace DGtal
+{
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class ConstantBinaryPredicate
+ /**
+ * Description of template class 'ConstantPointPredicate'
+ * \brief Aim: The predicate that returns always the same value boolCst
+ *
+ * @tparam boolCst any boolean value
+ */
+ template
+ struct ConstantBinaryPredicate
+ {
+ /**
+ * @param t1 first argument
+ * @param t2 second argument
+ * @tparam T1 type of the first argument
+ * @tparam T2 type of the second argument
+ * @return the value of the predicate
+ */
+ template
+ bool operator()( const T1& t1, const T2& t2 ) const
+ {
+ return boolCst;
+ }
+
+ }; // end of class ConstantBinaryPredicate
+
+ typedef ConstantBinaryPredicate TrueBinaryPredicate;
+ typedef ConstantBinaryPredicate FalseBinaryPredicate;
+
+} // namespace DGtal
+
+
+///////////////////////////////////////////////////////////////////////////////
+// Includes inline functions.
+//#include "BinaryPredicates.ih"
+
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // !defined BinaryPredicates_h
+
+#undef BinaryPredicates_RECURSES
+#endif // else defined(BinaryPredicates_RECURSES)
diff --git a/deformations/CMakeLists.txt b/deformations/CMakeLists.txt
index 8c9fd0b..f5605d0 100644
--- a/deformations/CMakeLists.txt
+++ b/deformations/CMakeLists.txt
@@ -1,5 +1,7 @@
PROJECT(Deformations)
+#C++11
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11")
#DGtal
CMAKE_MINIMUM_REQUIRED(VERSION 2.6)
@@ -12,10 +14,15 @@ LINK_DIRECTORIES(${DGTAL_LIBRARY_DIRS})
#Inclusion de Boost
include(FindBoost)
-find_package(Boost 1.36.0 REQUIRED program_options)
-message(STATUS "Found Boost: ${Boost_LIBRARIES} ")
+find_package(Boost 1.49.0 COMPONENTS program_options)
+message(STATUS "Found Boost: ${Boost_LIBRARIES} in ${Boost_LIBRARY_DIRS} and ${Boost_INCLUDE_DIRS} ")
+message(STATUS "More precisely: ${Boost_PROGRAM_OPTIONS_LIBRARY} ")
link_directories(${Boost_LIBRARY_DIRS})
include_directories(${Boost_INCLUDE_DIRS})
+link_directories(${Boost_LIBRARY_DIRS})
+
+# Workaround for strange bug while trying to include boost_program_options
+SET(BOOST_LIBRARIES ${BOOST_LIBRARIES} boost_program_options)
#fftw
INCLUDE_DIRECTORIES(/usr/include/)
@@ -25,6 +32,18 @@ LINK_DIRECTORIES(/usr/lib/)
ADD_EXECUTABLE(3dVolViewer 3dVolViewer)
TARGET_LINK_LIBRARIES(3dVolViewer ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES})
+ADD_EXECUTABLE(3dVolProjector 3dVolProjector)
+TARGET_LINK_LIBRARIES(3dVolProjector ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES})
+
+ADD_EXECUTABLE(generate3d generate3d)
+TARGET_LINK_LIBRARIES(generate3d ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES})
+
+#ADD_EXECUTABLE(testLocalDeformation3d testLocalDeformation3d)
+#TARGET_LINK_LIBRARIES(testLocalDeformation3d ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES})
+
+#ADD_EXECUTABLE(testLocalDeformation2d testLocalDeformation2d)
+#TARGET_LINK_LIBRARIES(testLocalDeformation2d ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES})
+
ADD_EXECUTABLE(testDisk testDisk)
TARGET_LINK_LIBRARIES(testDisk ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES} fftw3)
@@ -35,4 +54,5 @@ ADD_EXECUTABLE(deformation3d deformation3d)
TARGET_LINK_LIBRARIES(deformation3d ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES} fftw3)
ADD_EXECUTABLE(imageBlurring imageBlurring)
-TARGET_LINK_LIBRARIES(imageBlurring ${DGTAL_LIBRARIES} fftw3)
+TARGET_LINK_LIBRARIES(imageBlurring ${DGTAL_LIBRARIES} ${BOOST_LIBRARIES} fftw3)
+
diff --git a/deformations/DifferentialOperators.h b/deformations/DifferentialOperators.h
new file mode 100644
index 0000000..77d3bfa
--- /dev/null
+++ b/deformations/DifferentialOperators.h
@@ -0,0 +1,999 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file DifferentialOperators.h
+ * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ *
+ * @date 201/12/19
+ *
+ * Header file for module DifferentialOperators.cpp
+ *
+ * This file is part of the DGtal library.
+ */
+
+#if defined(DifferentialOperators_RECURSES)
+#error Recursive header files inclusion detected in DifferentialOperators.h
+#else // defined(DifferentialOperators_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define DifferentialOperators_RECURSES
+
+#if !defined DifferentialOperators_h
+/** Prevents repeated inclusion of headers. */
+#define DifferentialOperators_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include
+#include "DGtal/images/CConstImage.h"
+
+//////////////////////////////////////////////////////////////////////////////
+
+namespace DGtal
+{
+
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class ForwardDifference
+ /**
+ * Description of template class 'ForwardDifference'
+ * \brief Aim: Computes the forward difference at a point.
+ *
+ *
+ * @tparam TFonctor model of CPointFunctor
+ * @tparam TPointPredicate model of CPointPredicate
+ * @tparam TOutputValue type of returned value (default TFunctor::Value)
+ */
+ template
+ class ForwardDifference
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ ForwardDifference( Image& aStartingImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~ForwardDifference() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Difference.
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return first derivative along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myU;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class BackwardDifference
+ /**
+ * Description of template class 'BackwardDifference'
+ * \brief Aim: Computes the backward difference at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TImage type of image
+ * @tparam TOutputValue type of returned value (default TImage::Value)
+ */
+ template
+ class BackwardDifference
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ BackwardDifference( Image& aStartingImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~BackwardDifference() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Difference.
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return first derivative along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myU;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class CentralDifference
+ /**
+ * Description of template class 'CentralDifference'
+ * \brief Aim: Computes the backward difference at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TImage type of image
+ * @tparam TOutputValue type of returned value (default TImage::Value)
+ */
+ template
+ class CentralDifference
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ CentralDifference( Image& aStartingImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~CentralDifference() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Difference.
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return first derivative along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myU;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ };
+
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class Difference2
+ /**
+ * Description of template class 'Difference2'
+ * \brief Aim: Computes the second difference at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TImage type of image
+ * @tparam TOutputValue type of returned value (default TImage::Value)
+ */
+ template
+ class Difference2
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ Difference2( Image& aStartingImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~Difference2() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Second forward/backward difference.
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return first derivative along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myU;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class NormalizedDifference2
+ /**
+ * Description of template class 'NormalizedDifference2'
+ * \brief Aim: Computes the second difference at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TImage type of image
+ * @tparam TOutputValue type of returned value (default TImage::Value)
+ */
+ template
+ class NormalizedDifference2
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ NormalizedDifference2( Image& aStartingImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~NormalizedDifference2() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Second forward/backward difference
+ * normalized by the inverse of the gradient modulus
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return first derivative along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myU;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ // ------------------------- Internals --------------------------------
+ private:
+
+ /**
+ * Return the harmonic average of the inverse of @a aV1 and @a aV2
+ *
+ * @param aV1 a first value
+ * @param aV2 a second value
+ * @return the average
+ */
+ double average ( const Value& aV1, const Value& aN2 ) const;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class WeightedDifference2
+ /**
+ * Description of template class 'WeightedDifference2'
+ * \brief Aim: Computes the second difference at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TImage type of image
+ * @tparam TOutputValue type of returned value (default TImage::Value)
+ */
+ template
+ class WeightedDifference2
+ {
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef TOutputValue OutputValue;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+
+ /**
+ * Constructor.
+ *
+ * @param aImage any image
+ * @param aWImage any image of weights
+ * @param aGridStep any length (=1 by default)
+ */
+ WeightedDifference2( Image& aImage, Image& aWImage, const OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~WeightedDifference2() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Second forward/backward differences on @a aImage
+ * weighted by @a aWImage and normalized
+ * by the inverse of the gradient modulus
+ *
+ * @param aPoint the point where the derivative is computed
+ * @param aDim the axis along which the derivative is computed
+ * @return second derivative of @a aImage along axis @a aDim at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Dimension& aDim ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Reference on an image
+ */
+ Image& myImage;
+ /**
+ * Reference on the weights image
+ */
+ Image& myWImage;
+
+ /**
+ * Grid step
+ */
+ OutputValue myH;
+
+ // ------------------------- Internals --------------------------------
+ private:
+
+ /**
+ * Return the harmonic average of @a aV1 and @a aV2,
+ * @a aV1 and @a aV2 being normalized by @a aN1 and @a aN2.
+ *
+ * @param aV1 a first value
+ * @param aN1 any value dividing @a aV1
+ * @param aV2 a second value
+ * @param aN2 any value dividing @a aV2
+ * @return the normalized harmonic average of @a aV1 and @a aV2
+ */
+ double average ( const Value& aV1, const double& aN1,
+ const Value& aV2, const double& aN2 ) const;
+
+ };
+
+ ///////////////////////////////////////////////////////////////////
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class Gradient
+ /**
+ * Description of template class 'Gradient'
+ * \brief Aim: Computes the gradient at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TDifference type of directionnal differential operator
+ */
+ template
+ class Gradient
+ {
+
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TDifference FiniteDifference;
+ typedef typename FiniteDifference::Dimension Dimension;
+ static const typename FiniteDifference::Dimension dimension = FiniteDifference::dimension;
+ typedef typename FiniteDifference::Image Image;
+ typedef typename FiniteDifference::Point Point;
+
+ typedef DGtal::PointVector OutputValue;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+ /**
+ * Constructor.
+ *
+ * @param aD a finite difference operator
+ */
+ Gradient( const FiniteDifference& aD): myD(aD) {}
+
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ Gradient( typename FiniteDifference::Image& aStartingImage,
+ const typename FiniteDifference::OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~Gradient() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Gradient
+ *
+ * @param aPoint the point where the gradient is computed
+ * @return gradient of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Finite difference
+ */
+ FiniteDifference myD;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class UpwindGradient
+ /**
+ * Description of template class 'UpwindGradient'
+ * \brief Aim: Computes the gradient at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TDifference type of directionnal differential operator
+ */
+ template
+ class UpwindGradient
+ {
+
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ typedef DGtal::PointVector OutputValue;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ UpwindGradient( Image& aStartingImage,
+ const TOutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~UpwindGradient() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Gradient
+ *
+ * @param aPoint the point where the gradient is computed
+ * @param aVector displacement vector of the interface
+ * @return gradient of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, const Vector& aVector ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Forward difference
+ */
+ ForwardDifference myF;
+ /**
+ * Backward difference
+ */
+ BackwardDifference myB;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class GodunovGradient
+ /**
+ * Description of template class 'GodunovGradient'
+ * \brief Aim: Computes the gradient at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TDifference type of directionnal differential operator
+ */
+ template
+ class GodunovGradient
+ {
+
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ BOOST_CONCEPT_ASSERT(( concepts::CConstImage ));
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TImage Image;
+ typedef typename Image::Value Value;
+
+ typedef typename Image::Point Point;
+ typedef typename Image::Vector Vector;
+ typedef typename Image::Domain Domain;
+
+ typedef typename Image::Dimension Dimension;
+ static const typename Image::Dimension dimension = Image::dimension;
+
+ typedef DGtal::PointVector OutputValue;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image
+ * @param isPositive flag equal to 'true' if the
+ * interface moves in the direction of the normal
+ * and 'false' if it moves in the opposite direction
+ * (='true' by default)
+ * @param aGridStep any length (=1 by default)
+ */
+ GodunovGradient( Image& aStartingImage, bool isPositive = true,
+ const TOutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~GodunovGradient() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Gradient
+ *
+ * @param aPoint the point where the gradient is computed
+ * @return gradient of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint ) const;
+
+ /**
+ * Gradient
+ *
+ * @param aPoint the point where the gradient is computed
+ * @param isPositive flag equal to 'true' if the
+ * interface moves in the direction of the normal
+ * and 'false' if it moves in the opposite direction
+ * (='true' by default)
+ * @return gradient of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint, bool isPositive ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Forward difference
+ */
+ ForwardDifference myF;
+ /**
+ * Backward difference
+ */
+ BackwardDifference myB;
+ /**
+ * Flag
+ */
+ bool myIsPositive;
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class Divergence
+ /**
+ * Description of template class 'Divergence'
+ * \brief Aim: Computes the divergence at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TDifference type of directionnal differential operator
+ */
+ template
+ class Divergence
+ {
+
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TDifference FiniteDifference;
+ typedef typename FiniteDifference::Dimension Dimension;
+ static const typename FiniteDifference::Dimension dimension = FiniteDifference::dimension;
+ typedef typename FiniteDifference::Image Image;
+ typedef typename FiniteDifference::Point Point;
+ typedef typename FiniteDifference::OutputValue OutputValue;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+ /**
+ * Constructor.
+ *
+ * @param aD a finite difference operator
+ */
+ Divergence( const FiniteDifference& aD ): myD(aD) {}
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ Divergence( typename FiniteDifference::Image& aStartingImage,
+ const typename FiniteDifference::OutputValue& aGridStep = 1);
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~Divergence() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Divergence
+ *
+ * @param aPoint the point where the divergence is computed
+ * @return gradient of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Finite difference
+ */
+ FiniteDifference myD;
+
+ };
+
+ /////////////////////////////////////////////////////////////////////////////
+ // template class GradientModulus
+ /**
+ * Description of template class 'GradientModulus'
+ * \brief Aim: Computes the gradient modulus at a point
+ * of an image.
+ *
+ * @code
+ * @endcode
+ *
+ * @tparam TGradient type of gradient
+ */
+ template
+ class GradientModulus
+ {
+
+
+ // ----------------------- Types ------------------------------
+ public:
+
+ typedef TGradient Gradient;
+ typedef typename Gradient::Dimension Dimension;
+ static const typename Gradient::Dimension dimension = Gradient::dimension;
+ typedef typename Gradient::Point Point;
+
+ typedef double OutputValue;
+
+ // ----------------------- Standard services ------------------------------
+ public:
+
+ /**
+ * Constructor.
+ *
+ * @param aG any gradient operator
+ */
+ GradientModulus( const Gradient& aG) : myG(aG) {}
+
+ /**
+ * Constructor.
+ *
+ * @param aStartingImage any image of signed values
+ * @param aGridStep any length (=1 by default)
+ */
+ GradientModulus( typename Gradient::Image& aStartingImage,
+ const typename Gradient::OutputValue::Component& aGridStep = 1);
+
+
+ /**
+ * Destructor. Does nothing.
+ */
+ ~GradientModulus() {}
+
+ /**
+ * Checks the validity/consistency of the object.
+ * @return 'true' if the object is valid, 'false' otherwise.
+ */
+ bool isValid() const {return true;}
+
+ /**
+ * Gradient modulus
+ *
+ * @param aPoint the point where the gradient modulus is computed
+ * @return gradient modulus of @a myU at @ aPoint
+ */
+ OutputValue operator() ( const Point& aPoint ) const;
+
+ // ------------------------- Private Datas --------------------------------
+ private:
+
+ /**
+ * Gradient
+ */
+ Gradient myG;
+
+ };
+
+
+ ///////////////////////////////////////////////////////////////////
+
+
+
+} // namespace DGtal
+
+
+///////////////////////////////////////////////////////////////////////////////
+// Includes inline functions.
+#include "DifferentialOperators.ih"
+
+// //
+///////////////////////////////////////////////////////////////////////////////
+
+#endif // !defined DifferentialOperators_h
+
+#undef DifferentialOperators_RECURSES
+#endif // else defined(DifferentialOperators_RECURSES)
diff --git a/deformations/DifferentialOperators.ih b/deformations/DifferentialOperators.ih
new file mode 100644
index 0000000..a8132b5
--- /dev/null
+++ b/deformations/DifferentialOperators.ih
@@ -0,0 +1,470 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+
+/**
+ * @file DifferentialOperators.ih
+ * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
+ * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
+ *
+ * @date 2011/12/19
+ *
+ * @brief Implementation of inline methods defined in DifferentialOperators.h
+ *
+ * This file is part of the DGtal library.
+ */
+
+
+//////////////////////////////////////////////////////////////////////////////
+#include
+#include
+//////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+// IMPLEMENTATION of inline methods.
+///////////////////////////////////////////////////////////////////////////////
+
+///////////////////////////////////////////////////////////////////////////////
+// --------------------One-dimensional differential operators -----------------
+
+template
+inline
+DGtal::ForwardDifference
+::ForwardDifference( Image& aStartingImage, const OutputValue& aGridStep )
+ : myU( aStartingImage ), myH( aGridStep )
+{
+}
+
+template
+inline
+typename DGtal::ForwardDifference::OutputValue
+DGtal::ForwardDifference::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Point next = aPoint;
+ next[aDim] += 1;
+ if ( myU.domain().isInside(next) )
+ return ( (myU(next) - myU(aPoint)) / myH );
+ else
+ {
+ Point prev = aPoint;
+ prev[aDim] -= 1;
+ if ( myU.domain().isInside(prev) ) {
+ return ( (myU(aPoint) - myU(prev)) / myH );
+ } else {
+ ASSERT( false && "too small image" ) ;
+ throw std::domain_error("too small image");
+ }
+ }
+}
+
+template
+inline
+DGtal::BackwardDifference
+::BackwardDifference( Image& aStartingImage, const OutputValue& aGridStep )
+ : myU( aStartingImage ), myH( aGridStep )
+{
+}
+
+template
+inline
+typename DGtal::BackwardDifference::OutputValue
+DGtal::BackwardDifference::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Point prev = aPoint;
+ prev[aDim] -= 1;
+
+ if ( myU.domain().isInside(prev) )
+ return ( (myU(aPoint) - myU(prev)) / myH );
+ else
+ {
+ Point next = aPoint;
+ next[aDim] += 1;
+ if ( myU.domain().isInside(next) ) {
+ return ( (myU(next) - myU(aPoint)) / myH );
+ } else {
+ ASSERT( false && "too small image" ) ;
+ throw std::domain_error("too small image");
+ }
+ }
+}
+
+template
+inline
+DGtal::CentralDifference
+::CentralDifference( Image& aStartingImage, const OutputValue& aGridStep )
+ : myU( aStartingImage ), myH( aGridStep )
+{
+}
+
+template
+inline
+typename DGtal::CentralDifference::OutputValue
+DGtal::CentralDifference::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Domain d = myU.domain();
+ Point next = aPoint;
+ next[aDim] += 1;
+ bool nextInside = d.isInside(next);
+ Point prev = aPoint;
+ prev[aDim] -= 1;
+ bool prevInside = d.isInside(prev);
+
+ if (nextInside && prevInside )
+ return ( (myU(next) - myU(prev)) / (2*myH) );
+ else if ( nextInside && (!prevInside) )
+ return ( (myU(next) - myU(aPoint)) / myH );
+ else if ( (!nextInside) && prevInside )
+ return ( (myU(aPoint) - myU(prev)) / myH );
+ else {
+ ASSERT( false && "too small image" ) ;
+ throw std::domain_error("too small image");
+ }
+}
+
+template
+inline
+DGtal::Difference2
+::Difference2( Image& aStartingImage, const OutputValue& aGridStep )
+ : myU( aStartingImage ), myH( aGridStep )
+{
+}
+
+template
+inline
+typename DGtal::Difference2::OutputValue
+DGtal::Difference2::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Domain d = myU.domain();
+ Point next = aPoint;
+ next[aDim] += 1;
+ bool nextInside = d.isInside(next);
+ Point previous = aPoint;
+ previous[aDim] -= 1;
+ bool previousInside = d.isInside(previous);
+
+ if (nextInside && previousInside )
+ return ( (myU(next) - 2*myU(aPoint) + myU(previous) ) / (myH*myH) );
+ else
+ return 0;
+}
+
+template
+inline
+DGtal::NormalizedDifference2
+::NormalizedDifference2( Image& aStartingImage, const OutputValue& aGridStep )
+ : myU( aStartingImage ), myH( aGridStep )
+{
+}
+
+template
+inline
+typename DGtal::NormalizedDifference2::OutputValue
+DGtal::NormalizedDifference2::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Domain d = myU.domain();
+ Point next = aPoint;
+ next[aDim] += 1;
+ bool nextInside = d.isInside(next);
+ Point previous = aPoint;
+ previous[aDim] -= 1;
+ bool previousInside = d.isInside(previous);
+
+ typedef Gradient > Gradient;
+ GradientModulus m( myU );
+ double mod = m(aPoint);
+
+ if (nextInside && previousInside )
+ {
+ double nMod = m(next);
+ double pMod = m(previous);
+ return ( ( average(nMod,mod)*(myU(next) - myU(aPoint)) )
+ - ( average(mod,pMod)*(myU(aPoint) - myU(previous)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else if ( (!nextInside) && previousInside )
+ {
+ double pMod = m(previous);
+ return ( ( ( 1 - average(pMod,mod) )
+ *(myU(aPoint) - myU(previous)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else if ( (!previousInside) && nextInside )
+ {
+ double nMod = m(next);
+ return ( ( ( average(nMod,mod) - 1 )
+ *(myU(next) - myU(aPoint)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else {
+ ASSERT( false && "too small image" ) ;
+ throw std::domain_error("too small image");
+ }
+}
+
+template
+inline
+double
+DGtal::NormalizedDifference2::
+average( const Value& aV1, const Value& aV2 ) const
+{
+ double som = NumberTraits::castToDouble(aV1 + aV2);
+ if (som == 0) return 0;
+ else return 2/som;
+}
+
+template
+inline
+DGtal::WeightedDifference2
+::WeightedDifference2( Image& aImage, Image& aWImage, const OutputValue& aGridStep )
+ : myImage( aImage ), myWImage( aWImage), myH( aGridStep )
+{
+ ASSERT(aWImage.domain().lowerBound() == aImage.domain().lowerBound());
+ ASSERT(aWImage.domain().upperBound() == aImage.domain().upperBound());
+}
+
+template
+inline
+typename DGtal::WeightedDifference2::OutputValue
+DGtal::WeightedDifference2::
+operator() (const Point& aPoint, const Dimension& aDim) const
+{
+ ASSERT( (aDim >= 0) && (aDim < dimension) );
+
+ Domain d = myImage.domain();
+ Point next = aPoint;
+ next[aDim] += 1;
+ bool nextInside = d.isInside(next);
+ Point previous = aPoint;
+ previous[aDim] -= 1;
+ bool previousInside = d.isInside(previous);
+
+ typedef Gradient > Gradient;
+ GradientModulus m( myImage );
+ double mod = m(aPoint);
+
+ if (nextInside && previousInside )
+ {
+ double nMod = m(next);
+ double pMod = m(previous);
+ return ( ( average(myWImage(next),nMod,myWImage(aPoint),mod)
+ *(myImage(next) - myImage(aPoint)) )
+ - ( average(myWImage(aPoint),mod,myWImage(previous),pMod)
+ *(myImage(aPoint) - myImage(previous)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else if ( (!nextInside) && previousInside )
+ {
+ double pMod = m(previous);
+ return ( ( ( myWImage(aPoint)
+ - average(myWImage(previous),pMod,myWImage(aPoint),mod) )
+ *(myImage(aPoint) - myImage(previous)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else if ( (!previousInside) && nextInside )
+ {
+ double nMod = m(next);
+ return ( ( ( average(myWImage(next),nMod,myWImage(aPoint),mod)
+ - myWImage(aPoint) )
+ *(myImage(next) - myImage(aPoint)) )
+ / NumberTraits::castToDouble(myH*myH) );
+ }
+ else {
+ ASSERT( false && "too small image" ) ;
+ throw std::domain_error("too small image");
+ }
+
+}
+
+template
+inline
+double
+DGtal::WeightedDifference2::
+average( const Value& aV1, const double& aN1,
+ const Value& aV2, const double& aN2 ) const
+{
+ ASSERT(aV1 != 0);
+ ASSERT(aV2 != 0);
+ double q1 = aN1 / NumberTraits::castToDouble(aV1);
+ double q2 = aN2 / NumberTraits::castToDouble(aV2);
+ double sum = q1 + q2;
+ if (sum == 0) return 0;
+ else return 2/sum;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// -------------------- Gradient operators -----------------
+
+
+template
+inline
+DGtal::Gradient
+::Gradient( typename FiniteDifference::Image& aStartingImage,
+ const typename FiniteDifference::OutputValue& aGridStep )
+ : myD( FiniteDifference( aStartingImage, aGridStep ) )
+{
+}
+
+template
+inline
+typename DGtal::Gradient::OutputValue
+DGtal::Gradient::
+operator() (const Point& aPoint) const
+{
+ OutputValue g;
+ for (Dimension k = 0; k < dimension; ++k )
+ {
+ g[k] = myD(aPoint,k);
+ }
+ return g;
+}
+
+
+template
+inline
+DGtal::UpwindGradient
+::UpwindGradient( Image& aStartingImage, const TOutputValue& aGridStep )
+ : myF( ForwardDifference( aStartingImage, aGridStep ) ),
+ myB( BackwardDifference( aStartingImage, aGridStep ) )
+{
+}
+
+template
+inline
+typename DGtal::UpwindGradient::OutputValue
+DGtal::UpwindGradient::
+ operator() (const Point& aPoint, const Vector& aVector) const
+{
+ OutputValue g;
+ for (Dimension k = 0; k < dimension; ++k )
+ {
+ if (aVector[k] > 0)
+ g[k] = myB(aPoint,k);
+ else if (aVector[k] < 0)
+ g[k] = myF(aPoint,k);
+ else //== 0
+ g[k] = 0;
+ }
+ return g;
+}
+
+template
+inline
+DGtal::GodunovGradient
+::GodunovGradient( Image& aStartingImage, bool isPositive, const TOutputValue& aGridStep )
+ : myF( ForwardDifference( aStartingImage, aGridStep ) ),
+ myB( BackwardDifference( aStartingImage, aGridStep ) ),
+ myIsPositive( isPositive )
+{
+}
+
+template
+inline
+typename DGtal::GodunovGradient::OutputValue
+DGtal::GodunovGradient::
+ operator() (const Point& aPoint) const
+{
+ return operator() (aPoint,myIsPositive);
+}
+
+template
+inline
+typename DGtal::GodunovGradient::OutputValue
+DGtal::GodunovGradient::
+operator() (const Point& aPoint, bool isPositive) const
+{
+ OutputValue g;
+ for (Dimension k = 0; k < dimension; ++k )
+ {
+ if (isPositive)
+ {
+ typename OutputValue::Coordinate a = myB(aPoint,k);
+ if (a < 0) a = 0;
+ typename OutputValue::Coordinate b = myF(aPoint,k);
+ if (b > 0) b = 0;
+ g[k] = std::sqrt( std::max(a*a,b*b) );
+ }
+ else
+ {
+ typename OutputValue::Coordinate a = myB(aPoint,k);
+ if (a > 0) a = 0;
+ typename OutputValue::Coordinate b = myF(aPoint,k);
+ if (b < 0) b = 0;
+ g[k] = std::sqrt( std::max(a*a,b*b) );
+ }
+ }
+ return g;
+}
+
+
+template
+inline
+DGtal::GradientModulus
+::GradientModulus( typename Gradient::Image& aStartingImage,
+ const typename Gradient::OutputValue::Component& aGridStep )
+ : myG( Gradient( aStartingImage, aGridStep ) )
+{
+}
+
+template
+inline
+typename DGtal::GradientModulus::OutputValue
+DGtal::GradientModulus::
+operator() (const Point& aPoint) const
+{
+ typename Gradient::OutputValue g = myG( aPoint );
+ return g.norm( Gradient::OutputValue::L_2 );
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// ----------------------- Divergence Operator ------------------------------
+
+
+template
+inline
+DGtal::Divergence
+::Divergence( typename FiniteDifference::Image& aStartingImage,
+ const typename FiniteDifference::OutputValue& aGridStep )
+ : myD( FiniteDifference( aStartingImage, aGridStep ) )
+{
+}
+
+template
+inline
+typename DGtal::Divergence::OutputValue
+DGtal::Divergence::
+operator() (const Point& aPoint) const
+{
+ OutputValue sum = 0;
+ for (Dimension k = 0; k < dimension; ++k )
+ {
+ sum += myD(aPoint,k);
+ }
+ return sum;
+}
diff --git a/deformations/ExactDiffusionEvolver.h b/deformations/ExactDiffusionEvolver.h
index 804afa4..72f5cad 100644
--- a/deformations/ExactDiffusionEvolver.h
+++ b/deformations/ExactDiffusionEvolver.h
@@ -70,7 +70,7 @@ namespace DGtal
{
//ASSERT
- BOOST_CONCEPT_ASSERT(( CImage ));
+ BOOST_CONCEPT_ASSERT(( concepts::CImage ));
// ----------------------- Types ------------------------------
public:
diff --git a/deformations/ExactDiffusionEvolver.ih b/deformations/ExactDiffusionEvolver.ih
index a42aece..13cdf20 100644
--- a/deformations/ExactDiffusionEvolver.ih
+++ b/deformations/ExactDiffusionEvolver.ih
@@ -81,7 +81,7 @@ DGtal::ExactDiffusionEvolver::update(Image& aF, const double& aT)
for (unsigned int k = 0; k < Point::dimension; ++k)
{
//normalization and centering (to be in [-0.5,0.5[^d)
- double coord = (double) p.at(k)/v.at(k);
+ double coord = (double) p[k]/v[k];
if (coord >= 0.5) coord -= 1.0;
norm2 += coord*coord;
}
diff --git a/deformations/ExactReactionEvolver.h b/deformations/ExactReactionEvolver.h
index b3ac392..306c05e 100644
--- a/deformations/ExactReactionEvolver.h
+++ b/deformations/ExactReactionEvolver.h
@@ -65,7 +65,7 @@ namespace DGtal
{
//ASSERT
- BOOST_CONCEPT_ASSERT(( CImage ));
+ BOOST_CONCEPT_ASSERT(( concepts::CImage ));
// ----------------------- Types ------------------------------
public:
diff --git a/deformations/ExplicitReactionEvolver.h b/deformations/ExplicitReactionEvolver.h
index 1542b08..7521261 100644
--- a/deformations/ExplicitReactionEvolver.h
+++ b/deformations/ExplicitReactionEvolver.h
@@ -45,7 +45,7 @@
#include "DGtal/base/Common.h"
#include "DGtal/images/CImage.h"
-#include "DGtal/images/DifferentialOperators.h"
+#include "DifferentialOperators.h"
//////////////////////////////////////////////////////////////////////////////
@@ -68,7 +68,7 @@ namespace DGtal
{
//ASSERT
- BOOST_CONCEPT_ASSERT(( CImage ));
+ BOOST_CONCEPT_ASSERT(( concepts::CImage ));
// ----------------------- Types ------------------------------
public:
diff --git a/deformations/ExplicitReactionEvolver.ih b/deformations/ExplicitReactionEvolver.ih
index 33ea367..b8690e9 100644
--- a/deformations/ExplicitReactionEvolver.ih
+++ b/deformations/ExplicitReactionEvolver.ih
@@ -63,6 +63,9 @@ DGtal::ExplicitReactionEvolver::update(Image& aF, const double
{
ASSERT(aT >= 0);
+ // Result of the update
+ Image next_aF(aF.domain());
+
double k = 0;
if (myWithVolumeConservation)
k = force(aF);
@@ -80,12 +83,13 @@ DGtal::ExplicitReactionEvolver::update(Image& aF, const double
GodunovGradient gradient( aF, (myG >= 0) );
GradientModulus > m( gradient );
//new value
- double tmp = derivative (v)
- - myEpsilon * ( myG + f*m(*it) + k) * std::sqrt( 2*function(v) );
+ double tmp = derivative (v) - myEpsilon * ( myG + f*m(*it) + k) * std::sqrt( 2*function(v) );
Value newValue = v - ( aT / (myEpsilon * myEpsilon ) ) * tmp;
- aF.setValue( *it, static_cast(newValue) );
+ next_aF.setValue( *it, static_cast(newValue) );
}
+ // Update
+ aF = next_aF;
}
@@ -167,4 +171,6 @@ DGtal::operator<< ( std::ostream & out,
// //
///////////////////////////////////////////////////////////////////////////////
+/* GNU coding style */
+/* vim: set ts=2 sw=2 expandtab cindent cinoptions=>4,n-2,{2,^-2,:2,=2,g0,h2,p5,t0,+2,(0,u0,w1,m1 : */
diff --git a/deformations/FFT.h b/deformations/FFT.h
index c31794a..9e547bf 100644
--- a/deformations/FFT.h
+++ b/deformations/FFT.h
@@ -6,6 +6,7 @@
#include
#include
+#include
namespace DGtal{
@@ -20,7 +21,7 @@ namespace DGtal{
{
//ASSERT
- BOOST_CONCEPT_ASSERT(( CImage ));
+ BOOST_CONCEPT_ASSERT(( concepts::CImage ));
// ----------------------- Types ------------------------------
public:
diff --git a/deformations/FFT.ih b/deformations/FFT.ih
index 4a13d55..f16aad2 100644
--- a/deformations/FFT.ih
+++ b/deformations/FFT.ih
@@ -29,22 +29,21 @@ DGtal::FFT::compute(ComplexImage& anImage)
//image size
Vector v = myImage.extent();
- int t[dimension];
+ int t[dimension];
int n = 1;
- for (unsigned int i= 0; i < dimension; ++i)
- {
- t[i] = v.at(i);
- n *= t[i];
+ for (int i = 0; i < dimension; ++i)
+ { //row-major order
+ t[dimension-1-i] = v[i];
+ n *= v[i];
}
-//using fftw library
+ //using fftw library
//input and output arrays
fftw_complex* spatial_repr;
fftw_complex* frequency_repr;
spatial_repr=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*n);
frequency_repr=(fftw_complex*) fftw_malloc(sizeof(fftw_complex)*n);
- //fill the arrays
Domain d = myImage.domain();
typename Domain::ConstIterator it = d.begin();
typename Domain::ConstIterator itEnd = d.end();
diff --git a/deformations/FMM.h b/deformations/FMM.h
new file mode 100644
index 0000000..0c92582
--- /dev/null
+++ b/deformations/FMM.h
@@ -0,0 +1,498 @@
+/**
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program. If not, see .
+ *
+ **/
+
+#pragma once
+
+/**
+ * @file FMM.h
+ *
+ * @author Tristan Roussillon (\c
+ * tristan.roussillon@liris.cnrs.fr ) Laboratoire d'InfoRmatique en
+ * Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS,
+ * France
+ *
+ * @date 2012/01/17
+ *
+ * @brief Fast Marching Method for incremental distance transform
+ *
+ * This file is part of the DGtal library.
+ *
+ */
+
+#if defined(FMM_RECURSES)
+#error Recursive header files inclusion detected in FMM.h
+#else // defined(FMM_RECURSES)
+/** Prevents recursive inclusion of headers. */
+#define FMM_RECURSES
+
+#if !defined FMM_h
+/** Prevents repeated inclusion of headers. */
+#define FMM_h
+
+//////////////////////////////////////////////////////////////////////////////
+// Inclusions
+#include
+#include
+#include