• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • Examples
  • File List
  • File Members

/home/hauberg/Dokumenter/Capture/humim-tracker-0.1/src/OpenTissue/OpenTissue/core/containers/grid/util/grid_gradient.h

Go to the documentation of this file.
00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_GRADIENT_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_GRADIENT_H
00003 //
00004 // OpenTissue Template Library
00005 // - A generic toolbox for physics-based modeling and simulation.
00006 // Copyright (C) 2008 Department of Computer Science, University of Copenhagen.
00007 //
00008 // OTTL is licensed under zlib: http://opensource.org/licenses/zlib-license.php
00009 //
00010 #include <OpenTissue/configuration.h>
00011 
00012 namespace OpenTissue
00013 {
00014   namespace grid
00015   {
00016 
00017     template<typename grid_type, typename vector3_type>
00018     inline void gradient(
00019       grid_type const & grid
00020       , size_t i
00021       , size_t j
00022       , size_t k
00023       , vector3_type & gradient 
00024       )
00025     {
00026       using std::min;
00027 
00028       typedef typename grid_type::value_type            value_type;
00029       typedef typename vector3_type::value_type         real_type;
00030 
00031       static value_type const unused = grid.unused();
00032 
00033       size_t const & I = grid.I();
00034       size_t const & J = grid.J();
00035       size_t const & K = grid.K();
00036 
00037       size_t im1 = 0, jm1 = 0, km1 = 0;
00038       if ( i>0 )
00039         im1 = i - 1;
00040       if ( j>0 )
00041         jm1 = j - 1;
00042       if ( k>0 )
00043         km1 = k - 1;
00044 
00045       size_t ip1 = min( i + 1u, I - 1u );
00046       size_t jp1 = min( j + 1u, J - 1u );
00047       size_t kp1 = min( k + 1u, K - 1u );
00048 
00049       size_t idx     = ( k   * J + j )   * I + i;
00050       size_t idx_im1 = ( k   * J + j )   * I + im1;
00051       size_t idx_ip1 = ( k   * J + j )   * I + ip1;
00052       size_t idx_jm1 = ( k   * J + jm1 ) * I + i;
00053       size_t idx_jp1 = ( k   * J + jp1 ) * I + i;
00054       size_t idx_km1 = ( km1 * J + j )   * I + i;
00055       size_t idx_kp1 = ( kp1 * J + j )   * I + i;
00056 
00057       //--- Return Unused-vector if any of the values in the map is unused.
00058       gradient = vector3_type(unused,unused,unused);
00059 
00060       value_type s_idx = grid( idx );    
00061       if ( s_idx == unused )
00062         return;
00063 
00064       value_type s_im1 = grid( idx_im1 );
00065       if ( s_im1 == unused )
00066         return;
00067 
00068       value_type s_ip1 = grid( idx_ip1 );
00069       if ( s_ip1 == unused )
00070         return;
00071 
00072       value_type s_jm1 = grid( idx_jm1 );
00073       if ( s_jm1 == unused )
00074         return;
00075 
00076       value_type s_jp1 = grid( idx_jp1 );
00077       if ( s_jp1 == unused )
00078         return;
00079 
00080       value_type s_km1 = grid( idx_km1 );
00081       if ( s_km1 == unused )
00082         return;
00083 
00084       value_type s_kp1 = grid( idx_kp1 );
00085       if ( s_kp1 == unused )
00086         return;
00087 
00088       real_type r_idx = static_cast<real_type>( s_idx );
00089       real_type r_im1 = static_cast<real_type>( s_im1 );
00090       real_type r_ip1 = static_cast<real_type>( s_ip1 );
00091       real_type r_jm1 = static_cast<real_type>( s_jm1 );
00092       real_type r_jp1 = static_cast<real_type>( s_jp1 );
00093       real_type r_km1 = static_cast<real_type>( s_km1 );
00094       real_type r_kp1 = static_cast<real_type>( s_kp1 );
00095 
00096       if ( i == 0 )
00097       {
00098         gradient( 0 ) = ( r_ip1 - r_idx );
00099         gradient( 0 ) /= grid.dx();
00100       }
00101       else if ( i == ( grid.I() - 1 ) )
00102       {
00103         gradient( 0 ) = ( r_idx - r_im1 );
00104         gradient( 0 ) /= grid.dx();
00105       }
00106       else
00107       {
00108         gradient( 0 ) = ( r_ip1 - r_im1 );
00109         gradient( 0 ) /= ( 2 * grid.dx() );
00110       }
00111       if ( j == 0 )
00112       {
00113         gradient( 1 ) = ( r_jp1 - r_idx );
00114         gradient( 1 ) /= grid.dy();
00115       }
00116       else if ( j == ( grid.J() - 1 ) )
00117       {
00118         gradient( 1 ) = ( r_idx - r_jm1 );
00119         gradient( 1 ) /= grid.dy();
00120       }
00121       else
00122       {
00123         gradient( 1 ) = ( r_jp1 - r_jm1 );
00124         gradient( 1 ) /= ( 2 * grid.dy() );
00125       }
00126 
00127       if ( k == 0 )
00128       {
00129         gradient( 2 ) = ( r_kp1 - r_idx );
00130         gradient( 2 ) /= grid.dz();
00131       }
00132       else if ( k == ( grid.K() - 1 ) )
00133       {
00134         gradient( 2 ) = ( r_idx - r_km1 );
00135         gradient( 2 ) /= grid.dz();
00136       }
00137       else
00138       {
00139         gradient( 2 ) = ( r_kp1 - r_km1 );
00140         gradient( 2 ) /= ( 2 * grid.dz() );
00141       }
00142     }
00143 
00144     template<typename grid_iterator,typename vector3_type>
00145     inline void gradient (grid_iterator const & iter, vector3_type & gradient )
00146     {
00147       typedef typename grid_iterator::grid_type               grid_type;
00148 
00149       size_t       i   = iter.i();
00150       size_t       j   = iter.j();
00151       size_t       k   = iter.k();
00152       grid_type const & grid = iter.get_grid();
00153       gradient(grid, i, j, k, gradient);
00154     }
00155 
00156     template<typename grid_iterator>
00157     inline typename grid_iterator::math_types::vector3_type gradient (grid_iterator const & iter )
00158     {
00159       typedef typename grid_iterator::math_types::vector3_type  vector3_type;
00160 
00161       vector3_type gradient;
00162 
00163       gradient(iter, gradient);
00164       return gradient;
00165     }
00166 
00167   } // namespace grid
00168 } // namespace OpenTissue
00169 
00170 // OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_GRADIENT_H
00171 #endif

Generated on Thu Dec 1 2011 12:51:05 for HUMIM Tracker by  doxygen 1.7.1