Go to the documentation of this file.00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_LOCAL_MINIMA_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_GRID_LOCAL_MINIMA_H
00003
00004
00005
00006
00007
00008
00009
00010 #include <OpenTissue/configuration.h>
00011
00012 #include <OpenTissue/core/containers/grid/util/grid_idx2coord.h>
00013 #include <cmath>
00014
00015 namespace OpenTissue
00016 {
00017 namespace grid
00018 {
00019 namespace detail
00020 {
00031 template < typename grid_type >
00032 inline bool is_local_minima(
00033 size_t i
00034 , size_t j
00035 , size_t k
00036 , grid_type const & phi
00037 )
00038 {
00039 using std::min;
00040 static size_t idx[27];
00041 size_t I = phi.I();
00042 size_t J = phi.J();
00043 size_t K = phi.K();
00044 size_t im1 = ( i ) ? i - 1 : 0;
00045 size_t jm1 = ( j ) ? j - 1 : 0;
00046 size_t km1 = ( k ) ? k - 1 : 0;
00047 size_t ip1 = min( i + 1u, I - 1u );
00048 size_t jp1 = min( j + 1u, J - 1u );
00049 size_t kp1 = min( k + 1u, K - 1u );
00050 idx[0] = ( kp1 * J + jp1 ) * I + im1;
00051 idx[1] = ( k * J + jp1 ) * I + im1;
00052 idx[2] = ( km1 * J + jp1 ) * I + im1;
00053 idx[3] = ( kp1 * J + j ) * I + im1;
00054 idx[4] = ( k * J + j ) * I + im1;
00055 idx[5] = ( km1 * J + j ) * I + im1;
00056 idx[6] = ( kp1 * J + jm1 ) * I + im1;
00057 idx[7] = ( k * J + jm1 ) * I + im1;
00058 idx[8] = ( km1 * J + jm1 ) * I + im1;
00059 idx[9] = ( kp1 * J + jp1 ) * I + i;
00060 idx[10] = ( k * J + jp1 ) * I + i;
00061 idx[11] = ( km1 * J + jp1 ) * I + i;
00062 idx[12] = ( kp1 * J + j ) * I + i;
00063 idx[13] = ( k * J + j ) * I + i;
00064 idx[14] = ( km1 * J + j ) * I + i;
00065 idx[15] = ( kp1 * J + jm1 ) * I + i;
00066 idx[16] = ( k * J + jm1 ) * I + i;
00067 idx[17] = ( km1 * J + jm1 ) * I + i;
00068 idx[18] = ( kp1 * J + jp1 ) * I + ip1;
00069 idx[19] = ( k * J + jp1 ) * I + ip1;
00070 idx[20] = ( km1 * J + jp1 ) * I + ip1;
00071 idx[21] = ( kp1 * J + j ) * I + ip1;
00072 idx[22] = ( k * J + j ) * I + ip1;
00073 idx[23] = ( km1 * J + j ) * I + ip1;
00074 idx[24] = ( kp1 * J + jm1 ) * I + ip1;
00075 idx[25] = ( k * J + jm1 ) * I + ip1;
00076 idx[26] = ( km1 * J + jm1 ) * I + ip1;
00077
00078 for(size_t i=0;i<27u;++i)
00079 if( phi(idx[13]) > phi(idx[i]) )
00080 return false;
00081 return true;
00082 }
00083
00084 }
00085
00086
00096 template < typename grid_type,typename point_container >
00097 inline void local_minima_as_points(
00098 grid_type const & phi
00099 , point_container & points
00100 )
00101 {
00102 typedef typename point_container::value_type vector3_type;
00103 typedef typename grid_type::const_index_iterator const_index_iterator;
00104 typedef typename vector3_type::value_type real_type;
00105 typedef typename grid_type::value_type value_type;
00106
00107 const_index_iterator pend = phi.end();
00108 const_index_iterator p = phi.begin();
00109 for(;p!=pend;++p)
00110 {
00111 if(detail::is_local_minima(p.i(),p.j(),p.k(),phi))
00112 {
00113 vector3_type point;
00114 idx2coord(phi,p.i(),p.j(),p.k(),point);
00115 points.push_back(point);
00116 }
00117 }
00118 }
00119
00129 template < typename grid_type,typename point_container >
00130 inline void local_minima_as_points(
00131 grid_type const & phi
00132 , grid_type const & mask
00133 , point_container & points
00134 )
00135 {
00136 typedef typename point_container::value_type vector3_type;
00137 typedef typename grid_type::const_index_iterator const_index_iterator;
00138 typedef typename vector3_type::value_type real_type;
00139 typedef typename grid_type::value_type value_type;
00140 const_index_iterator end = phi.end();
00141 const_index_iterator p = phi.begin();
00142 const_index_iterator m = mask.begin();
00143 for(;p!=end;++p,++m)
00144 {
00145 if( (*m) <= 0 )
00146 continue;
00147 if(detail::is_local_minima(p.i(),p.j(),p.k(),phi))
00148 {
00149 vector3_type point;
00150 idx2coord(phi,p.i(),p.j(),p.k(),point);
00151 points.push_back(point);
00152 }
00153 }
00154 }
00155
00156 }
00157 }
00158
00159
00160 #endif