Go to the documentation of this file.00001 #ifndef OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_STRICT_GRID_LOCAL_MINIMA_H
00002 #define OPENTISSUE_CORE_CONTAINERS_GRID_UTIL_STRICT_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
00020 namespace detail
00021 {
00022
00033 template < typename grid_type >
00034 inline bool is_strict_local_minima(
00035 size_t i
00036 , size_t j
00037 , size_t k
00038 , grid_type const & phi
00039 )
00040 {
00041 using std::min;
00042 static size_t idx[27];
00043 size_t I = phi.I();
00044 size_t J = phi.J();
00045 size_t K = phi.K();
00046 size_t im1 = ( i ) ? i - 1 : 0;
00047 size_t jm1 = ( j ) ? j - 1 : 0;
00048 size_t km1 = ( k ) ? k - 1 : 0;
00049 size_t ip1 = min( i + 1u, I - 1u );
00050 size_t jp1 = min( j + 1u, J - 1u );
00051 size_t kp1 = min( k + 1u, K - 1u );
00052 idx[0] = ( kp1 * J + jp1 ) * I + im1;
00053 idx[1] = ( k * J + jp1 ) * I + im1;
00054 idx[2] = ( km1 * J + jp1 ) * I + im1;
00055 idx[3] = ( kp1 * J + j ) * I + im1;
00056 idx[4] = ( k * J + j ) * I + im1;
00057 idx[5] = ( km1 * J + j ) * I + im1;
00058 idx[6] = ( kp1 * J + jm1 ) * I + im1;
00059 idx[7] = ( k * J + jm1 ) * I + im1;
00060 idx[8] = ( km1 * J + jm1 ) * I + im1;
00061 idx[9] = ( kp1 * J + jp1 ) * I + i;
00062 idx[10] = ( k * J + jp1 ) * I + i;
00063 idx[11] = ( km1 * J + jp1 ) * I + i;
00064 idx[12] = ( kp1 * J + j ) * I + i;
00065 idx[13] = ( k * J + j ) * I + i;
00066 idx[14] = ( km1 * J + j ) * I + i;
00067 idx[15] = ( kp1 * J + jm1 ) * I + i;
00068 idx[16] = ( k * J + jm1 ) * I + i;
00069 idx[17] = ( km1 * J + jm1 ) * I + i;
00070 idx[18] = ( kp1 * J + jp1 ) * I + ip1;
00071 idx[19] = ( k * J + jp1 ) * I + ip1;
00072 idx[20] = ( km1 * J + jp1 ) * I + ip1;
00073 idx[21] = ( kp1 * J + j ) * I + ip1;
00074 idx[22] = ( k * J + j ) * I + ip1;
00075 idx[23] = ( km1 * J + j ) * I + ip1;
00076 idx[24] = ( kp1 * J + jm1 ) * I + ip1;
00077 idx[25] = ( k * J + jm1 ) * I + ip1;
00078 idx[26] = ( km1 * J + jm1 ) * I + ip1;
00079 for(size_t i=0;i<27u;++i)
00080 if( phi(idx[13]) >= phi(idx[i]) )
00081 return false;
00082 return true;
00083 }
00084
00085 }
00086
00096 template < typename grid_type,typename point_container >
00097 inline void strict_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_strict_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 strict_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_strict_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