Jerome Piovano

Research in Medical Imaging and Computer Vision

/home/jpiovano/src/Odyssee++/trunk/Libs/LevelSet/src/SegmentationModule/SoftShapeModule.h

00001 
00009 #ifndef SOFTSHAPEMODULE_H
00010 #define SOFTSHAPEMODULE_H
00011 
00012 #include <SegmentationModule/SegmentationModule.h>
00013 //#include <image/RecFilters.H>
00014 //#include <image/IIR.H>
00015 #include <DericheFilter.h> 
00016 
00017 
00018 namespace segmentation {
00019 
00020 
00028         template <unsigned DIM, typename Pixel = float>
00029         class SoftShapeModule : public SegmentationModule<DIM, Pixel>
00030         {
00031         public:
00032 
00033                 typedef           SegmentationModule<DIM, Pixel>  base;
00034                 typedef typename  base::LevelSet                  LevelSet; 
00035                 typedef typename  base::Index                     Index;
00036                 typedef typename  base::Schemes                   Schemes;
00037                 
00038                 typedef           BaseImage<DIM, float>           ImageFloat;
00039 
00040                 typedef           DericheFilter<DIM, float>       DericheFilter;
00041                 //typedef           Images::Filters::IIR::Separable::GaussianFilter<4,double> DericheFilter;
00042 
00043                 typedef           Pixels::RGB<unsigned char>      RGBPixel;
00044                 typedef           BaseImage<DIM, RGBPixel>        ColorImageType;
00045 
00046 
00049 
00051                 SoftShapeModule(float       _weight = 1.f,      
00052                                 float       _sigma = 20.0f,     
00053                                 int         _p = 2              
00054                                 );
00055 
00057                 virtual ~SoftShapeModule();
00058 
00060 
00062 
00064                 void set_sigma(float _sigma);
00066                 int & p();
00067 
00069 
00071 
00073                 float value(const LevelSet & phi, const Index & ind);
00075                 float energy(const LevelSet & phi, const Index & ind);
00077                 void init(const LevelSet & phi);
00079                 void update(const LevelSet & phi);
00081                 void save_result(const LevelSet & phi, std::string evodir, int iter);
00082 
00084                 void update_params(const LevelSet & phi);
00085         
00087 
00088         protected:
00089 
00090                 inline float _g(float z, float d, float epsilon)
00091                 {
00092                         if (z < d-epsilon) return 1;
00093                         if (z > d+epsilon) return 0;
00094                         return 0.5*(1 - (z-d)/epsilon - (1/M_PI)*sin(M_PI*(z-d)/epsilon) );
00095                 }
00096 
00099 
00100                 float          m_sigma;         
00101                 DericheFilter  m_deriche;       
00102                 ImageFloat     m_e1;            
00103                 ImageFloat     m_e2;            
00104                 int            m_p;             
00105 
00107 
00108         };
00109 
00110 
00114 
00115         // Constructors - Destructor
00116 
00117 
00118         template <unsigned DIM, typename Pixel>
00119                 SoftShapeModule<DIM, Pixel>::SoftShapeModule(float _weight, float _sigma, int _p) 
00120                 : SegmentationModule<DIM, Pixel>(_weight, true), m_sigma(_sigma), m_deriche(_sigma, false) /*m_deriche(_sigma, 0)*/, m_e1(), m_e2(), m_p(_p)
00121         { }
00122 
00123         template <unsigned DIM, typename Pixel>
00124                 SoftShapeModule<DIM, Pixel>::~SoftShapeModule() 
00125         { }
00126 
00127 
00129         // Modify/Access member data`
00130 
00131         
00132         template <unsigned DIM, typename Pixel>
00133                 void SoftShapeModule<DIM, Pixel>::set_sigma(float _sigma)
00134         { 
00135                 m_sigma = _sigma;
00136                 //std::cerr << "cant do that !" << std::endl;
00137                 m_deriche.set_isotropic_sigma(_sigma);
00138         }
00139 
00140         template <unsigned DIM, typename Pixel>
00141                 int & SoftShapeModule<DIM, Pixel>::p()
00142         {
00143                 return m_p;
00144         }
00145 
00146 
00148         // Member functions
00149 
00150 
00151         //ok donc globalement : dirac coupe pas mal 
00152         //=> histogramme de level set :  0     distance    0
00153         //                               10000 distance +- 0.5
00154         //                               15000 distance +- 0.7
00155         // dirac a +-0.5 => 0.5 
00156         // dirac a +-0.7 => 0.2 
00157         // max courbure * dirac est donc << au moin de moitiee a 1/rayoncourbure.
00158 
00159         // de plus foutre des dirac de partout meme sur le ballon pas top (mais marche ...) 
00160         // attention !! si dirac sur ballon et pas sur softshape produit des resultats differents.
00161         // second dirac ici diminue la fonction heaviside regul => 0.7 ou 0.5 a la place de 1.0
00162 
00163         // essayer peut etre a la place d empecher terme soft module aux endroit forte courbure d activer ballon => non car aucune evol au debut.
00164 
00165         // semble accentuer le bruit :/ du coup a la fin courbure estimee assez deg.
00166         
00167 
00168         template <unsigned DIM, typename Pixel>
00169                 float SoftShapeModule<DIM, Pixel>::value(const LevelSet & phi, const Index & ind)
00170         { 
00171           return - m_e2(ind);
00172 
00173                 //return Globals::delta(phi(ind))*(- m_e2(ind)); 
00174 
00175                 // ne pas utiliser avec thicknesss 0 car courbure foire au niveau des singularites
00176 
00177                 // 2/sigma car sigma = diametre ou cherche rayon.
00178                 // MAIS mean curv / 2 a cause dirac et statistic histogramme sur fonction dist => redeviens 1/m_sigma
00179                 // return Globals::delta(phi(ind))*_g(-Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)), 1/m_sigma - 0.1, 0.1)*(- m_e2(ind));
00180         }       
00181 
00182         template <unsigned DIM, typename Pixel>
00183                 float SoftShapeModule<DIM, Pixel>::energy(const LevelSet & phi, const Index & ind)
00184         {
00185                 //TODO 
00186                 return 0;
00187         }
00188 
00189         template <unsigned DIM, typename Pixel>
00190                 void SoftShapeModule<DIM, Pixel>::init(const LevelSet & phi)
00191         {
00192                 base::init(phi);
00193                 m_e1.resize(phi.shape());
00194                 m_e2.resize(phi.shape());
00195 
00196                 update_params(phi);
00197         }
00198 
00199         template <unsigned DIM, typename Pixel>
00200                 void SoftShapeModule<DIM, Pixel>::update(const LevelSet & phi)
00201         { 
00202                 update_params(phi);
00203         }
00204 
00205         template <unsigned DIM, typename Pixel>
00206                  void SoftShapeModule<DIM, Pixel>::update_params(const LevelSet & phi) 
00207         {
00208                 // computation of m_e1
00209                 for (typename ImageFloat::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){
00210                         const Pixel      _phi = phi(ind);
00211                         if (_phi>=0) m_e1(ind) = 1.0f;
00212                         else         m_e1(ind) = 0.0f;
00213                 }
00214 
00215                 #ifdef PRINT_TIME
00216                 timeval _time;
00217                 Globals::initChrono(_time);
00218                 #endif
00219 
00220                 //m_e1 = Images::Filter<ImageFloat>(m_e1, m_deriche);
00221                 m_deriche.filter(m_e1);
00222 
00223                 #ifdef PRINT_TIME
00224                 Globals::printChrono("          - 1st filter : ", _time);
00225                 Globals::initChrono(_time);
00226                 #endif
00227 
00228                 // computation of m_e2
00229                 for (typename ImageFloat::template iterator<domain> ind=phi.begin() ; ind!=phi.end() ; ++ind ){
00230                         const Pixel      _phi = phi(ind);
00231                         if (_phi>=0)  m_e2(ind) = pow(m_e1(ind), m_p-1);
00232                         else          m_e2(ind) = 0.0f;
00233                 }
00234 
00235                 #ifdef PRINT_TIME
00236                 Globals::initChrono(_time);
00237                 #endif
00238 
00239                 //m_e2 = Images::Filter<ImageFloat>(m_e2, m_deriche);
00240                 m_deriche.filter(m_e2);
00241 
00242                 #ifdef PRINT_TIME
00243                 Globals::printChrono("          - 2nd filter : ", _time);
00244                 #endif
00245 
00246                 m_e2 *= m_p;
00247                 m_e2 += pow(m_e1, m_p);
00248         }
00249 
00250         template <unsigned DIM, typename Pixel>
00251                 void SoftShapeModule<DIM, Pixel>::save_result(const LevelSet & phi, std::string evodir, int iter)
00252         {
00253                 ImageFloat tmp(phi.shape());
00254                 std::ostringstream name;
00255 
00256                 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00257                         if(phi(i) > 0) tmp(i) = m_e1(i);
00258                         else           tmp(i) = 0;
00259                 }
00260 
00261                 name.str("");
00262                 name << evodir << "e1_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00263                 std::ofstream file(name.str().c_str());
00264                 file << format("inrimage") << tmp;
00265                 std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
00266 
00267                 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00268                         if(phi(i) > 0) tmp(i) = m_e2(i);
00269                         else           tmp(i) = 0;
00270                 }
00271 
00272                 name.str("");
00273                 name << evodir << "e2_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00274                 std::ofstream file2(name.str().c_str());
00275                 file2 << format("inrimage") << tmp;
00276                 std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
00277 
00278         #if 0 // test changement topologie (marche pas)
00279 
00280                 //ColorImageType tmp(phi.shape());
00281                 ImageFloat tmp(phi.shape());
00282                 tmp = 255;
00283                 
00284                 //#pragma omp parallel for
00285                 for (int i = 0 ; i < phi.nb_points() ; i++){
00286 
00287                         const Index& ind = phi.get_point(i);
00288                         const Pixel _phi = phi(ind);
00289 
00290                         int inf_coords[DIM];
00291                         int sup_coords[DIM];
00292                         Schemes::find_neighbours(phi, ind, inf_coords, sup_coords);
00293 
00294                         for (int i=0 ; i<DIM ; i++){
00295                                 Index inf = ind;
00296                                 Index sup = ind;
00297                                 inf(i+1) = inf_coords[i];
00298                                 sup(i+1) = sup_coords[i];
00299                                 if ((_phi * phi(inf) <= 0) || (_phi * phi(sup) <= 0)){
00300                                         //std::cout  << M_PI << " " << 2.0/m_sigma << " " << -Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)) <<" | " <<_g(-Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)), 2.0/m_sigma - 0.05, 0.05) << std::endl;
00301                                         //RGBPixel color((int) (255.0f * (1.0f - _g(-Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)), 2/m_sigma - 0.05, 0.05))),0,0);
00302                                         //RGBPixel color((int) (235.0f * (_g(-Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)), 2/m_sigma - 0.05, 0.05))),0,0);
00303 
00304                                         //tmp(ind) = color;
00305                                         tmp(ind) = (int) (235.0f * (_g(-Schemes::mean_curvature(phi,ind)*Globals::delta(phi(ind)), 1/m_sigma - 0.1, 0.1)));
00306                                 }
00307                         }
00308                 }
00309 
00310                 std::ostringstream name;
00311                 name.str("");
00312                 name << evodir << "soft_shape_result_" << std::setw(6) << std::setfill('0') << iter << ".inr"; //<< ".ppm";
00313                 std::ofstream file(name.str().c_str());
00314                 file << format("inrimage") << tmp;
00315                 std::cout <<  " - Saved " << name.str() << " ..."  << std::endl;
00316         #endif
00317 
00318         }
00319 }
00320 
00321 #endif // SoftShapeModule_H
For further information, please contact Jerome Piovano - Last update 2008-02-08