00001
00009 #ifndef MUMFORDSHAHREGIONMODULE_H
00010 #define MUMFORDSHAHREGIONMODULE_H
00011
00012 #include <SegmentationModule/RegionModule.h>
00013
00014 namespace segmentation {
00015
00023 template <unsigned DIM, typename Pixel = float, typename PixelImage = float>
00024 class MumfordShahRegionModule : public RegionModule<DIM, Pixel, PixelImage>
00025 {
00026 public:
00027
00028 typedef RegionModule<DIM, Pixel, PixelImage> base;
00029 typedef typename base::LevelSet LevelSet;
00030 typedef typename base::Index Index;
00031 typedef typename base::Histogram Histogram;
00032 typedef typename base::ImageType ImageType;
00033 typedef typename base::Schemes Schemes;
00034
00037
00039 MumfordShahRegionModule(ImageType * _data,
00040 float _weight = 1.f,
00041 float _mu = 1.f,
00042 float _dt = 0.05f,
00043 int _iter = 100,
00044 bool _adaptive = true
00045 );
00046
00048 virtual ~MumfordShahRegionModule();
00049
00051
00053
00055 float& mu();
00057 int& iter();
00059 float& dt();
00060
00062
00064
00066 float value(const LevelSet & phi, const Index & ind);
00068 float energy(const LevelSet & phi, const Index & ind);
00070 void init(const LevelSet & phi);
00072 void update(const LevelSet & phi);
00074 void save_result(const LevelSet & phi, std::string evodir, int iter);
00075
00077 void update_mumfordshah_params(const LevelSet & phi);
00078
00080
00081 protected:
00082
00083
00084
00086 float poisson_mumford_shah_pos(const LevelSet & phi, const Index & ind);
00088 float poisson_mumford_shah_neg(const LevelSet & phi, const Index & ind);
00089
00092
00093 float m_mu;
00094 float m_dt;
00095 int m_iter;
00096 ImageType m_u_pos;
00097 ImageType m_u_neg;
00098
00100
00101 };
00102
00103
00105
00106
00107
00108 template <unsigned DIM, typename Pixel, typename PixelImage>
00109 MumfordShahRegionModule<DIM, Pixel, PixelImage>::MumfordShahRegionModule(ImageType * _data, float _weight, float _mu, float _dt, int _iter, bool _adaptive)
00110 :RegionModule<DIM, Pixel, PixelImage>(_data, _weight, _adaptive), m_mu(_mu), m_dt(_dt), m_iter(_iter), m_u_pos(_data->shape()), m_u_neg(_data->shape())
00111 { }
00112
00113 template <unsigned DIM, typename Pixel, typename PixelImage>
00114 MumfordShahRegionModule<DIM, Pixel, PixelImage>::~MumfordShahRegionModule()
00115 { }
00116
00117
00119
00120
00121
00122 template <unsigned DIM, typename Pixel, typename PixelImage>
00123 float& MumfordShahRegionModule<DIM, Pixel, PixelImage>::mu()
00124 {
00125 return m_mu;
00126 }
00127
00128 template <unsigned DIM, typename Pixel, typename PixelImage>
00129 int& MumfordShahRegionModule<DIM, Pixel, PixelImage>::iter()
00130 {
00131 return m_iter;
00132 }
00133
00134 template <unsigned DIM, typename Pixel, typename PixelImage>
00135 float& MumfordShahRegionModule<DIM, Pixel, PixelImage>::dt()
00136 {
00137 return m_dt;
00138 }
00139
00140
00142
00143
00144
00145 template <unsigned DIM, typename Pixel, typename PixelImage>
00146 float MumfordShahRegionModule<DIM, Pixel, PixelImage>::value(const LevelSet & phi, const Index & ind)
00147 {
00148 PixelImage val = (* base::m_data)(ind);
00149
00150 float ll_pos = Maths::sqr(val - m_u_pos(ind)) + m_mu*Schemes::sq_norm_gradient(m_u_pos,ind);
00151 float ll_neg = Maths::sqr(val - m_u_neg(ind)) + m_mu*Schemes::sq_norm_gradient(m_u_neg,ind);
00152
00153 return Globals::delta(phi(ind)) * (ll_pos - ll_neg);
00154 }
00155
00156 template <unsigned DIM, typename Pixel, typename PixelImage>
00157 float MumfordShahRegionModule<DIM, Pixel, PixelImage>::energy(const LevelSet & phi, const Index & ind)
00158 {
00159
00160 return 0.0f;
00161 }
00162
00163 template <unsigned DIM, typename Pixel, typename PixelImage>
00164 void MumfordShahRegionModule<DIM, Pixel, PixelImage>::init(const LevelSet & phi)
00165 {
00166 base::init(phi);
00167 update_mumfordshah_params(phi);
00168 }
00169
00170
00171 template <unsigned DIM, typename Pixel, typename PixelImage>
00172 void MumfordShahRegionModule<DIM, Pixel, PixelImage>::update(const LevelSet & phi)
00173 {
00174 update_mumfordshah_params(phi);
00175 }
00176
00177 template <unsigned DIM, typename Pixel, typename PixelImage>
00178 void MumfordShahRegionModule<DIM, Pixel, PixelImage>::update_mumfordshah_params(const LevelSet & phi)
00179 {
00180 m_u_pos = 0.0f;
00181 m_u_neg = 0.0f;
00182
00183 BaseImage<DIM, Pixel> temp_u_pos(base::m_data->shape());
00184 BaseImage<DIM, Pixel> temp_u_neg(base::m_data->shape());
00185
00186 #ifdef PRINT_TIME
00187 timeval _time;
00188 Globals::initChrono(_time);
00189 #endif
00190
00192
00193 Index tmp[phi.size()];
00194 int i = 0;
00195 for (typename LevelSet::template iterator<domain> it=phi.begin() ; it!=phi.end() ; ++it )
00196 tmp[i++] = it.position();
00197
00198 for (int t=0 ; t<m_iter ; t++) {
00199 #pragma omp parallel for
00200 for (int i = 0 ; i < phi.size() ; i++){
00201 Index & ind = tmp[i];
00202 if (phi(ind) > 0) temp_u_pos(ind) = poisson_mumford_shah_pos(phi, ind);
00203 else temp_u_neg(ind) = poisson_mumford_shah_neg(phi, ind);
00204 }
00205 #pragma omp parallel for
00206 for (int i = 0 ; i < phi.size() ; i++){
00207 Index & ind = tmp[i];
00208 if (phi(ind) > 0) m_u_pos(ind) += m_dt*temp_u_pos(ind);
00209 else m_u_neg(ind) += m_dt*temp_u_neg(ind);
00210 }
00211
00212 #if 0
00213
00214 for (int t=0 ; t<m_iter ; t++) {
00215 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00216 if (phi(i) > 0) temp_u_pos(i) = poisson_mumford_shah(m_u_pos, phi, i.position());
00217 else temp_u_neg(i) = poisson_mumford_shah(m_u_neg, phi, i.position());
00218 }
00219
00220 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00221 if (phi(i) > 0) m_u_pos(i) += m_dt*temp_u_pos(i);
00222 else m_u_neg(i) += m_dt*temp_u_neg(i);
00223 }
00224 #endif
00225
00226 #if 0
00227 ImageType tmp1(phi.shape());
00228 ImageType tmp2(phi.shape());
00229 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00230 if (phi(i) > 0) { tmp1(i) = temp_u_pos(i); tmp2(i) = m_u_pos(i);}
00231 else { tmp1(i) = temp_u_neg(i); tmp2(i) = m_u_neg(i);}
00232 }
00233
00234 std::ostringstream name;
00235 std::ofstream file1;
00236 std::ofstream file2;
00237
00238 name.str("");
00239 name << "residus_" << std::setw(6) << std::setfill('0') << t << ".inr";
00240 file1.open(name.str().c_str());
00241 file1 << format("inrimage") << tmp1 ;
00242 file1.close();
00243 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00244
00245 name.str("");
00246 name << "poisson_" << std::setw(6) << std::setfill('0') << t << ".inr";
00247 file2.open(name.str().c_str());
00248 file2 << format("inrimage") << tmp2 ;
00249 file2.close();
00250 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00251 #endif
00252 }
00253
00254 #ifdef PRINT_TIME
00255 Globals::printChrono(" - poisson PDE : ", _time);
00256 Globals::initChrono(_time);
00257 #endif
00258
00259
00260 phi.extend_data(m_u_pos,-1);
00261 phi.extend_data(m_u_neg,+1);
00262
00263 #ifdef PRINT_TIME
00264 Globals::printChrono(" - data extension in band : ", _time);
00265 #endif
00266 }
00267
00268
00270
00271
00272
00273 template <unsigned DIM, typename Pixel, typename PixelImage>
00274 void MumfordShahRegionModule<DIM, Pixel, PixelImage>::save_result(const LevelSet & phi, std::string evodir, int iter)
00275 {
00276 base::save_result(phi, evodir, iter);
00277
00278 ImageType tmp(phi.shape());
00279 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00280 if (phi(i) > 0) tmp(i) = m_u_pos(i);
00281 else tmp(i) = m_u_neg(i);
00282 }
00283
00284 std::ostringstream name;
00285 std::ofstream file;
00286
00287 name.str("");
00288 name << evodir << "MumfordShah_result_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00289 file.open(name.str().c_str());
00290 file << format("inrimage") << tmp ;
00291 file.close();
00292 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00293
00294 }
00295
00296
00297
00298 template <unsigned DIM, typename Pixel, typename PixelImage>
00299 float MumfordShahRegionModule<DIM, Pixel, PixelImage>::poisson_mumford_shah_pos(const LevelSet & phi, const Index & ind)
00300 {
00301 float result = DIM * m_u_pos(ind);
00302 result *= -2;
00303
00304 int inf_coords[DIM];
00305 int sup_coords[DIM];
00306 Schemes::find_neighbours(m_u_pos, ind, inf_coords, sup_coords);
00307
00308 Index sup;
00309 Index inf;
00310 for (int i=0 ; i<DIM ; i++){
00311 sup = ind;
00312 inf = ind;
00313 sup(i+1) = sup_coords[i];
00314 inf(i+1) = inf_coords[i];
00315 result += (phi(sup) <= 0) ? m_u_pos(ind) : m_u_pos(sup);
00316 result += (phi(inf) <= 0) ? m_u_pos(ind) : m_u_pos(inf);
00317 }
00318
00319 return m_mu*result + (* base::m_data)(ind) - m_u_pos(ind);
00320 }
00321
00322 template <unsigned DIM, typename Pixel, typename PixelImage>
00323 float MumfordShahRegionModule<DIM, Pixel, PixelImage>::poisson_mumford_shah_neg(const LevelSet & phi, const Index & ind)
00324 {
00325 float result = DIM * m_u_neg(ind);
00326 result *= -2;
00327
00328 int inf_coords[DIM];
00329 int sup_coords[DIM];
00330 Schemes::find_neighbours(m_u_neg, ind, inf_coords, sup_coords);
00331
00332 Index sup;
00333 Index inf;
00334 for (int i=0 ; i<DIM ; i++){
00335 sup = ind;
00336 inf = ind;
00337 sup(i+1) = sup_coords[i];
00338 inf(i+1) = inf_coords[i];
00339 result += (phi(sup) > 0) ? m_u_neg(ind) : m_u_neg(sup);
00340 result += (phi(inf) > 0) ? m_u_neg(ind) : m_u_neg(inf);
00341 }
00342
00343 return m_mu*result + (* base::m_data)(ind) - m_u_neg(ind);
00344 }
00345 }
00346
00347 #endif // MUMFORDSHAHREGIONMODULE_H