Go to the documentation of this file.00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_SECOND_DERIVATVE_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_SECOND_DERIVATVE_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/core/containers/grid/util/grid_gradient.h>
00013 #include <OpenTissue/core/containers/grid/util/grid_hessian.h>
00014 #include <OpenTissue/core/math/math_vector3.h>
00015 #include <OpenTissue/core/math/math_matrix3x3.h>
00016
00017 namespace OpenTissue
00018 {
00019 namespace grid
00020 {
00021
00022 template<typename grid_type, typename real_type>
00023 inline void second_derivative(
00024 grid_type const & grid
00025 , size_t i
00026 , size_t j
00027 , size_t k
00028 , real_type & derivative
00029 )
00030 {
00031 using std::fabs;
00032
00033 typedef OpenTissue::math::Vector3<real_type> vector3_type;
00034 typedef OpenTissue::math::Matrix3x3<real_type> matrix3x3_type;
00035
00036 static vector3_type g;
00037 static matrix3x3_type H;
00038
00039
00040 real_type const too_small = 10e-7;
00041
00042 gradient(grid, i, j, k, g );
00043
00044 real_type tmp = g*g;
00045 if(fabs(tmp) < too_small)
00046 {
00047 derivative = 0.0;
00048 return;
00049 }
00050
00051 hessian(grid, i, j, k, H);
00052 derivative = (1.0)/(tmp)* (g * (H *g));
00053 }
00054
00055 template<typename grid_iterator,typename real_type>
00056 inline void second_derivative (grid_iterator const & iter, real_type & derivative)
00057 {
00058 typedef typename grid_iterator::grid_type grid_type;
00059
00060 size_t i = iter.i();
00061 size_t j = iter.j();
00062 size_t k = iter.k();
00063
00064 grid_type const & grid = iter.get_grid();
00065 second_derivative(grid, i, j, k, derivative);
00066 }
00067
00068 template<typename grid_iterator>
00069 inline typename grid_iterator::math_types::real_type
00070 second_derivative( grid_iterator const & iter )
00071 {
00072 typedef typename grid_iterator::math_types::real_type real_type;
00073
00074 real_type derivative;
00075
00076 second_derivative(iter, derivative);
00077
00078 return derivative;
00079 }
00080
00081 }
00082
00083 }
00084
00085
00086 #endif