00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_GRADIENT_AT_POINT_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_GRADIENT_AT_POINT_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/core/math/math_trillinear.h>
00013 #include <OpenTissue/core/containers/grid/util/grid_enclosing_indices.h>
00014 #include <OpenTissue/core/containers/grid/util/grid_gradient.h>
00015
00016 namespace OpenTissue
00017 {
00018 namespace grid
00019 {
00020
00021 template<typename grid_type,typename vector3_type>
00022 inline vector3_type gradient_at_point (grid_type const & grid, vector3_type const & point )
00023 {
00024 typedef typename grid_type::value_type value_type;
00025
00026 const static value_type infty = grid.infinity();
00027 const static value_type unused = grid.unused();
00028
00029 size_t i0, j0, k0, i1, j1, k1;
00030 enclosing_indices(grid, point, i0, j0, k0, i1, j1, k1 );
00031
00032 vector3_type g000,g001,g010,g011,g100,g101,g110,g111;
00033 gradient( grid, i0, j0, k0, g000 );
00034 gradient( grid, i1, j0, k0, g001 );
00035 gradient( grid, i0, j1, k0, g010 );
00036 gradient( grid, i1, j1, k0, g011 );
00037 gradient( grid, i0, j0, k1, g100 );
00038 gradient( grid, i1, j0, k1, g101 );
00039 gradient( grid, i0, j1, k1, g110 );
00040 gradient( grid, i1, j1, k1, g111 );
00041
00042 if ( g000( 0 ) == unused )
00043 return vector3_type( infty, infty, infty );
00044 if ( g001( 0 ) == unused )
00045 return vector3_type( infty, infty, infty );
00046 if ( g010( 0 ) == unused )
00047 return vector3_type( infty, infty, infty );
00048 if ( g011( 0 ) == unused )
00049 return vector3_type( infty, infty, infty );
00050 if ( g100( 0 ) == unused )
00051 return vector3_type( infty, infty, infty );
00052 if ( g101( 0 ) == unused )
00053 return vector3_type( infty, infty, infty );
00054 if ( g110( 0 ) == unused )
00055 return vector3_type( infty, infty, infty );
00056 if ( g111( 0 ) == unused )
00057 return vector3_type( infty, infty, infty );
00058
00059 typename vector3_type::value_type s = ( point(0) - ( i0*grid.dx() + grid.min_coord(0) ) ) / grid.dx();
00060 typename vector3_type::value_type t = ( point(1) - ( j0*grid.dy() + grid.min_coord(1) ) ) / grid.dy();
00061 typename vector3_type::value_type u = ( point(2) - ( k0*grid.dz() + grid.min_coord(2) ) ) / grid.dz();
00062
00063 return OpenTissue::math::trillinear(g000, g001, g010, g011, g100, g101, g110, g111, s, t, u ) ;
00064 }
00065
00066 }
00067 }
00068
00069
00070 #endif