00001
00009 #ifndef LOCALSTATSREGIONMODULE_H
00010 #define LOCALSTATSREGIONMODULE_H
00011
00012
00013
00014
00015
00016
00017
00018 #include <SegmentationModule/RegionModule.h>
00019 #include <DericheFilter.h>
00020 #include <Schemes.h>
00021 #include <image/RecFilters.H>
00022 #include <image/IIR.H>
00023
00024 namespace segmentation {
00025
00026
00033 template <unsigned DIM, typename Pixel = float, typename PixelImage = float>
00034 class LocalGaussianRegionModule : public RegionModule<DIM, Pixel, PixelImage>
00035 {
00036 public:
00037
00038 typedef RegionModule<DIM, Pixel, PixelImage> base;
00039 typedef typename base::LevelSet LevelSet;
00040 typedef typename base::Index Index;
00041 typedef typename base::ImageType ImageType;
00042
00043 typedef BaseImage<DIM, double> ImageFloat;
00044
00045 #ifdef BLUR_ODYSSEE_LIB
00046 typedef Images::Filters::IIR::Separable::GaussianFilter<4,double> DericheFilter;
00047 #else
00048 typedef DericheFilter<DIM, double> DericheFilter;
00049 #endif
00050
00053
00055 LocalGaussianRegionModule(ImageType * _data,
00056 float _weight = 1.f,
00057 float _sigma = 20.0f,
00058 bool _var = false
00059 );
00060
00062 virtual ~LocalGaussianRegionModule();
00063
00065
00067
00069 void set_sigma(float _sigma);
00071 void use_variance(bool _var=true);
00072
00074
00076
00078 float value(const LevelSet & phi, const Index & ind);
00080 float energy(const LevelSet & phi, const Index & ind);
00082 void measures(const LevelSet & phi, const Index & ind, float meas[]);
00084 void init(const LevelSet & phi);
00086 void update(const LevelSet & phi);
00088 void save_result(const LevelSet & phi, std::string evodir, int iter);
00089
00091 void update_params(const LevelSet & phi);
00092
00094
00096
00097 protected:
00098
00099 float m_sigma;
00100 DericheFilter m_deriche;
00101
00102 ImageFloat m_chi_pos;
00103 ImageFloat m_mean_pos;
00104 ImageFloat m_mean_neg;
00105 ImageFloat m_var_pos;
00106 ImageFloat m_var_neg;
00107 float m_image_variance;
00108 bool m_var;
00109
00110 #ifdef SECOND_TERM
00111 ImageFloat m_term1_pos;
00112 ImageFloat m_term1_neg;
00113 ImageFloat m_term2_pos;
00114 ImageFloat m_term2_neg;
00115 #endif
00116
00117
00118 };
00119
00120
00124
00125
00126
00127
00128 template <unsigned DIM, typename Pixel, typename PixelImage>
00129 LocalGaussianRegionModule<DIM, Pixel, PixelImage>::LocalGaussianRegionModule(ImageType * _data, float _weight, float _sigma, bool _var)
00130 : RegionModule<DIM, Pixel, PixelImage>(_data, _weight, true), m_sigma(_sigma), m_chi_pos(_data->shape()), m_mean_pos(_data->shape()), m_mean_neg(_data->shape()), m_var_pos(_data->shape()), m_var_neg(_data->shape()), m_image_variance(-1.0f), m_var(_var)
00131 #ifdef SECOND_TERM
00132 ,m_term1_pos(_data->shape()),m_term2_pos(_data->shape()),m_term1_neg(_data->shape()),m_term2_neg(_data->shape())
00133 #endif
00134 #ifdef BLUR_ODYSSEE_LIB
00135 ,m_deriche(_sigma, 0)
00136 #else
00137 ,m_deriche(_sigma)
00138 #endif
00139 { }
00140
00141 template <unsigned DIM, typename Pixel, typename PixelImage>
00142 LocalGaussianRegionModule<DIM, Pixel, PixelImage>::~LocalGaussianRegionModule()
00143 { }
00144
00145
00147
00148
00149
00150 template <unsigned DIM, typename Pixel, typename PixelImage>
00151 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::set_sigma(float _sigma)
00152 {
00153 m_sigma = _sigma;
00154
00155 #ifdef BLUR_ODYSSEE_LIB
00156 std::cerr << "cant do that !" << std::endl;
00157 #else
00158 m_deriche.set_isotropic_sigma(_sigma);
00159 #endif
00160 }
00161
00162 template <unsigned DIM, typename Pixel, typename PixelImage>
00163 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::use_variance(bool _var)
00164 {
00165 m_var = _var;
00166 }
00167
00168
00170
00171
00172
00173 template <unsigned DIM, typename Pixel, typename PixelImage>
00174 float LocalGaussianRegionModule<DIM, Pixel, PixelImage>::value(const LevelSet & phi, const Index & ind)
00175 {
00176 PixelImage val = (* base::m_data)(ind);
00177
00178 float mean_pos = m_mean_pos(ind)/m_chi_pos(ind);
00179 float mean_neg = m_mean_neg(ind)/(1-m_chi_pos(ind));
00180
00181 float var_pos = m_var ? m_var_pos(ind)/m_chi_pos(ind) - Maths::sqr(mean_pos) : m_image_variance;
00182 float var_neg = m_var ? m_var_neg(ind)/(1-m_chi_pos(ind)) - Maths::sqr(mean_neg) : m_image_variance;
00183
00184 var_pos = (var_pos<m_image_variance/2)? m_image_variance/2 : var_pos;
00185 var_neg = (var_neg<m_image_variance/2)? m_image_variance/2 : var_neg;
00186
00187 float ll_pos = -0.5f*log(2*M_PI*var_pos) - 0.5f*Maths::sqr(val-mean_pos)/var_pos;
00188 float ll_neg = -0.5f*log(2*M_PI*var_neg) - 0.5f*Maths::sqr(val-mean_neg)/var_neg;
00189
00190 #ifdef SECOND_TERM
00191 ll_pos += val*m_term1_pos(ind)/var_pos - m_term2_pos(ind)/var_pos;
00192 ll_neg += val*m_term1_neg(ind)/var_neg - m_term2_neg(ind)/var_neg;
00193 #endif
00194
00195
00196 float result = ((mean_pos-mean_neg < 0)? 1 : -1) * (ll_pos - ll_neg)/3.0;
00197
00198 return result;
00199 }
00200
00201 template <unsigned DIM, typename Pixel, typename PixelImage>
00202 float LocalGaussianRegionModule<DIM, Pixel, PixelImage>::energy(const LevelSet & phi, const Index & ind)
00203 {
00204
00205 return 0;
00206 }
00207
00208 template <unsigned DIM, typename Pixel, typename PixelImage>
00209 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::measures(const LevelSet & phi, const Index & ind, float meas[])
00210 {
00211 float mean_pos = m_mean_pos(ind)/m_chi_pos(ind);
00212 float mean_neg = m_mean_neg(ind)/(1-m_chi_pos(ind));
00213 float var_pos = m_var ? m_var_pos(ind)/m_chi_pos(ind) - Maths::sqr(mean_pos) : m_image_variance;
00214 float var_neg = m_var ? m_var_neg(ind)/(1-m_chi_pos(ind)) - Maths::sqr(mean_neg) : m_image_variance;
00215
00216
00217 meas[0] = mean_pos;
00218 meas[1] = mean_neg;
00219 meas[2] = var_pos;
00220 meas[3] = var_neg;
00221 }
00222
00223 template <unsigned DIM, typename Pixel, typename PixelImage>
00224 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::init(const LevelSet & phi)
00225 {
00226 base::init(phi);
00227 update_params(phi);
00228 }
00229
00230 template <unsigned DIM, typename Pixel, typename PixelImage>
00231 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::update(const LevelSet & phi)
00232 {
00233 update_params(phi);
00234 }
00235
00236 template <unsigned DIM, typename Pixel, typename PixelImage>
00237 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::update_params(const LevelSet & phi)
00238 {
00239
00240 if(m_image_variance < 0.f) {
00241 PixelImage min;
00242 PixelImage max;
00243 float mean;
00244
00245 Globals::stats<DIM,PixelImage>(*base::m_data, min, max, mean, &m_image_variance);
00246 }
00247
00248
00249 for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
00250 const Pixel _phi = phi(ind);
00251 const PixelImage _val = (*base::m_data)(ind);
00252
00253 if (_phi>=0) {
00254 m_chi_pos(ind) = 1.0f;
00255 m_mean_pos(ind) = _val;
00256 m_mean_neg(ind) = 0.0f;
00257 m_var_pos(ind) = Maths::sqr(_val);
00258 m_var_neg(ind) = 0.0f;
00259 }
00260 else {
00261 m_chi_pos(ind) = 0.0f;
00262 m_mean_pos(ind) = 0.0f;
00263 m_mean_neg(ind) = _val;
00264 m_var_pos(ind) = 0.0f;
00265 m_var_neg(ind) = Maths::sqr(_val);
00266 }
00267 }
00268
00269 #pragma omp parallel sections
00270 {
00271 #pragma omp section
00272
00273 #ifdef BLUR_ODYSSEE_LIB
00274
00275
00276 m_chi_pos = Images::Filter<ImageFloat>(m_chi_pos, m_deriche);
00277 #else
00278 m_deriche.filter(m_chi_pos);
00279 #endif
00280 #pragma omp section
00281
00282 #ifdef BLUR_ODYSSEE_LIB
00283 m_mean_pos = Images::Filter<ImageFloat>(m_mean_pos, m_deriche);
00284 #else
00285 m_deriche.filter(m_mean_pos);
00286 #endif
00287 #pragma omp section
00288
00289 #ifdef BLUR_ODYSSEE_LIB
00290 m_mean_neg = Images::Filter<ImageFloat>(m_mean_neg, m_deriche);
00291 #else
00292 m_deriche.filter(m_mean_neg);
00293 #endif
00294 }
00295
00296 if(m_var) {
00297 #pragma omp parallel sections
00298 {
00299 #pragma omp section
00300
00301 #ifdef BLUR_ODYSSEE_LIB
00302 m_var_pos = Images::Filter<ImageFloat>(m_var_pos, m_deriche);
00303 #else
00304 m_deriche.filter(m_var_pos);
00305 #endif
00306 #pragma omp section
00307
00308 #ifdef BLUR_ODYSSEE_LIB
00309 m_var_neg = Images::Filter<ImageFloat>(m_var_neg, m_deriche);
00310 #else
00311 m_deriche.filter(m_var_neg);
00312 #endif
00313 }
00314 }
00315
00316 #ifdef SECOND_TERM
00317
00318 for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
00319 const Pixel _phi = phi(ind);
00320 const PixelImage _val = (*base::m_data)(ind);
00321
00322 if (_phi>=0) {
00323 m_term1_pos(ind) = ((*base::m_data)(ind) - (m_mean_pos(ind)/m_chi_pos(ind))) / m_chi_pos(ind);
00324 m_term2_pos(ind) = ((*base::m_data)(ind) - (m_mean_pos(ind)/m_chi_pos(ind))) * (m_mean_pos(ind)/m_chi_pos(ind)) / m_chi_pos(ind);
00325 m_term1_neg(ind) = 0.0f;
00326 m_term2_neg(ind) = 0.0f;
00327 }
00328 else {
00329 m_term1_pos(ind) = 0.0f;
00330 m_term2_pos(ind) = 0.0f;
00331 m_term1_neg(ind) = ((*base::m_data)(ind) - (m_mean_neg(ind)/(1 - m_chi_pos(ind)))) / (1 - m_chi_pos(ind));
00332 m_term2_neg(ind) = ((*base::m_data)(ind) - (m_mean_neg(ind)/(1 - m_chi_pos(ind)))) * (m_mean_neg(ind)/(1 - m_chi_pos(ind))) / (1 - m_chi_pos(ind));
00333 }
00334 }
00335 #pragma omp parallel sections
00336 {
00337 #pragma omp section
00338
00339 #ifdef BLUR_ODYSSEE_LIB
00340 m_term1_pos = Images::Filter<ImageFloat>(m_term1_pos, m_deriche);
00341 #else
00342 m_deriche.filter(m_term1_pos);
00343 #endif
00344 #pragma omp section
00345
00346 #ifdef BLUR_ODYSSEE_LIB
00347 m_term2_pos = Images::Filter<ImageFloat>(m_term2_pos, m_deriche);
00348 #else
00349 m_deriche.filter(m_term2_pos);
00350 #endif
00351 #pragma omp section
00352
00353 #ifdef BLUR_ODYSSEE_LIB
00354 m_term1_neg = Images::Filter<ImageFloat>(m_term1_neg, m_deriche);
00355 #else
00356 m_deriche.filter(m_term1_neg);
00357 #endif
00358 #pragma omp section
00359
00360 #ifdef BLUR_ODYSSEE_LIB
00361 m_term2_neg = Images::Filter<ImageFloat>(m_term2_neg, m_deriche);
00362 #else
00363 m_deriche.filter(m_term2_neg);
00364 #endif
00365 }
00366 #endif
00367 }
00368
00369 template <unsigned DIM, typename Pixel, typename PixelImage>
00370 void LocalGaussianRegionModule<DIM, Pixel, PixelImage>::save_result(const LevelSet & phi, std::string evodir, int iter)
00371 {
00372
00373 ImageFloat tmp(phi.shape());
00374 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00375 if(phi(i) > 0) tmp(i) = m_mean_pos(i)/m_chi_pos(i);
00376 else tmp(i) = m_mean_neg(i)/(1-m_chi_pos(i));
00377 }
00378
00379 std::ostringstream name;
00380 std::ofstream file;
00381
00382 #if 1 //
00383 name.str("");
00384 name << evodir << "LocalGaussian_" << std::setw(3) << std::setfill('0') << m_sigma << "_mean_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00385 file.open(name.str().c_str());
00386 file << format("inrimage") << tmp ;
00387 file.close();
00388 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00389 #endif
00390
00391 #if 1 //... separements + chi
00392 name.str("");
00393 name << evodir << "LocalGaussian_chi_" << std::setw(3) << std::setfill('0') << m_sigma << "_chi_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00394 std::ofstream file2;
00395 file2.open(name.str().c_str());
00396 file2 << format("inrimage") << m_chi_pos ;
00397 file2.close();
00398 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00399
00400 name.str("");
00401 name << evodir << "LocalGaussian_mean_pos_" << std::setw(3) << std::setfill('0') << m_sigma << "_mean_pos_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00402 std::ofstream file3;
00403 file3.open(name.str().c_str());
00404 file3 << format("inrimage") << m_mean_pos ;
00405 file3.close();
00406 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00407
00408 name.str("");
00409 name << evodir << "LocalGaussian_mean_neg_" << std::setw(3) << std::setfill('0') << m_sigma << "_mean_neg_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00410 std::ofstream file4;
00411 file4.open(name.str().c_str());
00412 file4 << format("inrimage") << m_mean_neg ;
00413 file4.close();
00414 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00415 #endif
00416
00417 #if 0 // ... que pos et neg sans blur
00418
00419 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00420 if(phi(i) > 0) tmp(i) = (* base::m_data)(i);
00421 else tmp(i) = 0;
00422 }
00423
00424 name.str("");
00425 name << evodir << "LocalGaussian_pos_" << std::setw(3) << std::setfill('0') << m_sigma << "_chi_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00426 std::ofstream file_pos2;
00427 file_pos2.open(name.str().c_str());
00428 file_pos2 << format("inrimage") << tmp;
00429 file_pos2.close();
00430 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00431
00432 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00433 if(phi(i) > 0) tmp(i) = 0;
00434 else tmp(i) = (* base::m_data)(i);
00435 }
00436
00437 name.str("");
00438 name << evodir << "LocalGaussian_neg_" << std::setw(3) << std::setfill('0') << m_sigma << "_chi_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00439 std::ofstream file_neg2;
00440 file_neg2.open(name.str().c_str());
00441 file_neg2 << format("inrimage") << tmp;
00442 file_neg2.close();
00443 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00444 #endif
00445
00446 #if 1 // ... image normale + contour + voisinage couleur bleu et jaune
00447
00448 typedef Pixels::RGB<unsigned char> RGBPixel;
00449 typedef BaseImage<DIM, RGBPixel> ColorImageType;
00450 ColorImageType tmp2(phi.shape());
00451 tmp2 = RGBPixel(0,0,0);
00452
00453
00454
00455
00456 int point_x = 90;
00457 int point_y = 100;
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509
00510
00511
00512
00513
00514
00515
00516
00517
00518
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534
00535
00536
00537
00538
00539 point_x = 90;
00540 point_y = 100;
00541 int sigmtmp = 9;
00542
00543
00544 RGBPixel neg(215,152,0);
00545 RGBPixel pos(255,252,0);
00546
00547
00548
00549
00550
00551
00552 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00553 const Pixel _phi = phi(i);
00554 const PixelImage _val = (*base::m_data)(i);
00555
00556 double coeff = floor(30 - sqrt(Maths::sqr(i(1)-point_x) + Maths::sqr(i(2)-point_y))) ;
00557 coeff = (coeff < 0)? 0 : coeff;
00558 coeff = ( ((int) coeff)%3 == 0)? coeff/30.0 : 0;
00559
00560 if(phi(i) > 0) {
00561 tmp2(i) += coeff*pos;
00562
00563 PixelImage _min = 0;
00564 _min = ((tmp2(i).green() + _min > 230)? 230 - tmp2(i).green() :_min);
00565 _min = ((tmp2(i).red() + _min > 230)? 230 - tmp2(i).red() :_min);
00566 _min = ((tmp2(i).blue() + _min > 230)? 230 - tmp2(i).blue() :_min);
00567 tmp2(i) += 0;
00568
00569 }
00570 else {
00571 tmp2(i) += coeff*neg;
00572
00573 PixelImage _min = 0;
00574 _min = ((tmp2(i).green() + _min > 230)? 230 - tmp2(i).green() :_min);
00575 _min = ((tmp2(i).red() + _min > 230)? 230 - tmp2(i).red() :_min);
00576 _min = ((tmp2(i).blue() + _min > 230)? 230 - tmp2(i).blue() :_min);
00577 tmp2(i) += 0;
00578 }
00579 }
00580
00581
00582
00583 point_x = 120 + 22;
00584 point_y = 210 + 40;
00585 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00586 const Pixel _phi = phi(i);
00587 const PixelImage _val = (*base::m_data)(i);
00588 double coeff = exp(-0.5f*(Maths::sqr(i(1)-point_x) + Maths::sqr(i(2)-point_y))/Maths::sqr( sigmtmp ))/2;
00589
00590 if(phi(i) > 0) {
00591 tmp2(i) += coeff*pos;
00592
00593 PixelImage _min = 0;
00594 _min = ((tmp2(i).green() + _min > 230)? 230 - tmp2(i).green() :_min);
00595 _min = ((tmp2(i).red() + _min > 230)? 230 - tmp2(i).red() :_min);
00596 _min = ((tmp2(i).blue() + _min > 230)? 230 - tmp2(i).blue() :_min);
00597 tmp2(i) += 0;
00598 }
00599 else {
00600 tmp2(i) += coeff*neg ;
00601
00602 PixelImage _min = 0;
00603 _min = ((tmp2(i).green() + _min > 230)? 230 - tmp2(i).green() :_min);
00604 _min = ((tmp2(i).red() + _min > 230)? 230 - tmp2(i).red() :_min);
00605 _min = ((tmp2(i).blue() + _min > 230)? 230 - tmp2(i).blue() :_min);
00606 tmp2(i) += 0;
00607 }
00608
00609 if(phi.is_in_band(i.position())){
00610 const Index& ind = i.position();
00611 const float _phi = phi(ind);
00612 int inf_coords[DIM];
00613 int sup_coords[DIM];
00614 Schemes<DIM,float>::find_neighbours(phi, ind, inf_coords, sup_coords);
00615
00616 for (int i=0 ; i<DIM ; i++){
00617 Index inf = ind;
00618 Index sup = ind;
00619 inf(i+1) = inf_coords[i];
00620 sup(i+1) = sup_coords[i];
00621 if ((_phi * phi(inf) <= 0) || (_phi * phi(sup) <= 0))
00622 tmp2(ind) = (_phi>0)? pos : neg ;
00623 }
00624 }
00625 }
00626
00627
00628
00629 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00630 const PixelImage _val = (* base::m_data)(i);
00631
00632 PixelImage _min = _val;
00633 _min = ((tmp2(i).green() + _min > 230)? 230 - tmp2(i).green() :_min);
00634 _min = ((tmp2(i).red() + _min > 230)? 230 - tmp2(i).red() :_min);
00635 _min = ((tmp2(i).blue() + _min > 230)? 230 - tmp2(i).blue() :_min);
00636 tmp2(i) += _min;
00637
00638 if(phi.is_in_band(i.position())){
00639 const Index& ind = i.position();
00640 const float _phi = phi(ind);
00641 int inf_coords[DIM];
00642 int sup_coords[DIM];
00643 Schemes<DIM,float>::find_neighbours(phi, ind, inf_coords, sup_coords);
00644
00645 for (int i=0 ; i<DIM ; i++){
00646 Index inf = ind;
00647 Index sup = ind;
00648 inf(i+1) = inf_coords[i];
00649 sup(i+1) = sup_coords[i];
00650 if ((_phi * phi(inf) <= 0) || (_phi * phi(sup) <= 0))
00651 tmp2(ind) = (_phi>0)? pos: neg;
00652 }
00653 }
00654 }
00655
00656 name.str("");
00657 name << evodir << "LocalGaussian_couleur_" << std::setw(6) << std::setfill('0') << iter << ".ppm";
00658 file.open(name.str().c_str());
00659 file << format("rawppm") << tmp2 ;
00660 file.close();
00661 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00662
00663 #endif
00664
00665 #if 0 // pour voir la valeur de la variance locale et ou elle est utilisee (si > m_image_variance)
00666 if (m_var){
00667 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00668 if(phi(i) > 0) tmp(i) = m_var_pos(i)/m_chi_pos(i) - Maths::sqr(m_mean_pos(i)/m_chi_pos(i));
00669 else tmp(i) = m_var_neg(i)/(1-m_chi_pos(i)) - Maths::sqr(m_mean_neg(i)/(1-m_chi_pos(i)));
00670 }
00671
00672 std::ofstream file2;
00673
00674 name.str("");
00675 name << evodir << "LocalGaussian_" << std::setw(3) << std::setfill('0') << m_sigma << "_var_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00676 file2.open(name.str().c_str());
00677 file2 << format("inrimage") << tmp ;
00678 file2.close();
00679 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00680
00681
00682 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00683 tmp(i) = (tmp(i) > m_image_variance)?1:0;
00684 }
00685
00686 std::ofstream file3;
00687
00688 name.str("");
00689 name << evodir << "LocalGaussian_" << std::setw(3) << std::setfill('0') << m_sigma << "_var_mask_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00690 file3.open(name.str().c_str());
00691 file3 << format("inrimage") << tmp ;
00692 file3.close();
00693 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00694 }
00695 else{
00696 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00697 tmp(i) = m_image_variance;
00698
00699 }
00700
00701 std::ofstream file2;
00702
00703 name.str("");
00704 name << evodir << "LocalGaussian_" << std::setw(3) << std::setfill('0') << m_sigma << "_var_" << std::setw(6) << std::setfill('0') << iter << ".inr";
00705 file2.open(name.str().c_str());
00706 file2 << format("inrimage") << tmp ;
00707 file2.close();
00708 std::cout << " - Saved " << name.str() << " ..." << std::endl;
00709 }
00710 #endif
00711
00712 #if 0 // DEBUG ajout ligne a fichier resultat
00713
00714 int point_x = 533;
00715 int point_y = 500;
00716
00717 if(iter == 0){
00718 name.str("");
00719 name << evodir << "chi_pos_1" ;
00720 std::ofstream chi_pos_1(name.str().c_str(), std::ios_base::app);
00721 name.str("");
00722 name << evodir << "mean_pos_no_norm_1" ;
00723 std::ofstream mean_pos_no_norm_1(name.str().c_str(), std::ios_base::app);
00724 name.str("");
00725 name << evodir << "mean_neg_no_norm_1" ;
00726 std::ofstream mean_neg_no_norm_1(name.str().c_str(), std::ios_base::app);
00727 name.str("");
00728 name << evodir << "mean_pos_1" ;
00729 std::ofstream mean_pos_1(name.str().c_str(), std::ios_base::app);
00730 name.str("");
00731 name << evodir << "mean_neg_1" ;
00732 std::ofstream mean_neg_1(name.str().c_str(), std::ios_base::app);
00733 name.str("");
00734 name << evodir << "var_pos_1" ;
00735 std::ofstream var_pos_1(name.str().c_str(), std::ios_base::app);
00736 name.str("");
00737 name << evodir << "var_neg_1" ;
00738 std::ofstream var_neg_1(name.str().c_str(), std::ios_base::app);
00739 name.str("");
00740 name << evodir << "chi_pos_2" ;
00741 std::ofstream chi_pos_2(name.str().c_str(), std::ios_base::app);
00742 name.str("");
00743 name << evodir << "chi_neg_2" ;
00744 std::ofstream chi_neg_2(name.str().c_str(), std::ios_base::app);
00745 name.str("");
00746 name << evodir << "mean_pos_no_norm_2" ;
00747 std::ofstream mean_pos_no_norm_2(name.str().c_str(), std::ios_base::app);
00748 name.str("");
00749 name << evodir << "mean_neg_no_norm_2" ;
00750 std::ofstream mean_neg_no_norm_2(name.str().c_str(), std::ios_base::app);
00751 name.str("");
00752 name << evodir << "mean_pos_2" ;
00753 std::ofstream mean_pos_2(name.str().c_str(), std::ios_base::app);
00754 name.str("");
00755 name << evodir << "mean_neg_2" ;
00756 std::ofstream mean_neg_2(name.str().c_str(), std::ios_base::app);
00757 name.str("");
00758 name << evodir << "var_pos_2" ;
00759 std::ofstream var_pos_2(name.str().c_str(), std::ios_base::app);
00760 name.str("");
00761 name << evodir << "var_neg_2" ;
00762 std::ofstream var_neg_2(name.str().c_str(), std::ios_base::app);
00763 name.str("");
00764 name << evodir << "speed_var" ;
00765 std::ofstream speed_var(name.str().c_str(), std::ios_base::app);
00766 name.str("");
00767 name << evodir << "speed_no_var" ;
00768 std::ofstream speed_no_var(name.str().c_str(), std::ios_base::app);
00769 name.str("");
00770
00771 for (typename LevelSet::template iterator<domain> i=phi.begin() ; i!=phi.end() ; ++i ){
00772 if (((*i)(1) == point_x) && ((*i)(2) == point_y)){
00773 chi_pos_1 << m_sigma << " " << m_chi_pos(i) << std::endl;
00774 mean_pos_no_norm_1 << m_sigma << " " << m_mean_pos(i) << std::endl;
00775 mean_neg_no_norm_1 << m_sigma << " " << m_mean_neg(i) << std::endl;
00776 mean_pos_1 << m_sigma << " " << m_mean_pos(i)/m_chi_pos(i) << std::endl;
00777 mean_neg_1 << m_sigma << " " << m_mean_neg(i)/(1-m_chi_pos(i)) << std::endl;
00778 var_pos_1 << m_sigma << " " << m_var_pos(i)/m_chi_pos(i) - Maths::sqr(m_mean_pos(i)/m_chi_pos(i)) << std::endl;
00779 var_neg_1 << m_sigma << " " << m_var_neg(i)/(1-m_chi_pos(i)) - Maths::sqr(m_mean_neg(i)/(1-m_chi_pos(i))) << std::endl;
00780 m_var = true;
00781 speed_var << m_sigma << " " << this->value(phi, i.position()) << std::endl;
00782 m_var = false;
00783 speed_no_var << m_sigma << " " << this->value(phi, i.position()) << std::endl;
00784 m_var = true;
00785 }
00786 }
00787
00788 double mean_pos = 0;
00789 double mean_neg = 0;
00790 double chi_pos = 0;
00791 double chi_neg = 0;
00792 for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
00793 const Pixel _phi = phi(ind);
00794 const PixelImage _val = (*base::m_data)(ind);
00795 double coeff = exp(-0.5f*(Maths::sqr(ind(1)-point_x) + Maths::sqr(ind(2)-point_y))/Maths::sqr(m_sigma))/(Maths::sqr(m_sigma)*2*M_PI);
00796
00797 if (_phi>=0) {
00798 chi_pos += coeff;
00799 mean_pos += _val*coeff;
00800 }
00801 else {
00802 chi_neg += coeff;
00803 mean_neg += _val*coeff;
00804 }
00805 }
00806
00807 chi_pos_2 << m_sigma << " " << chi_pos << std::endl;
00808 chi_neg_2 << m_sigma << " " << chi_neg << std::endl;
00809 mean_pos_no_norm_2 << m_sigma << " " << mean_pos << std::endl;
00810 mean_neg_no_norm_2 << m_sigma << " " << mean_neg << std::endl;
00811
00812 mean_pos = mean_pos/chi_pos;
00813 mean_neg = mean_neg/chi_neg;
00814
00815 float var_pos = 0;
00816 float var_neg = 0;
00817 for (typename ImageType::template iterator<domain> ind=base::m_data->begin() ; ind!=base::m_data->end() ; ++ind ){
00818 const Pixel _phi = phi(ind);
00819 const PixelImage _val = (*base::m_data)(ind);
00820 double coeff = exp(-0.5f*(Maths::sqr(ind(1)-point_x) + Maths::sqr(ind(2)-point_y))/Maths::sqr(m_sigma))/(Maths::sqr(m_sigma)*2*M_PI);
00821
00822 if (_phi>=0) var_pos += Maths::sqr(_val-mean_pos)*coeff;
00823 else var_neg += Maths::sqr(_val-mean_neg)*coeff;
00824 }
00825 var_pos = var_pos/chi_pos;
00826 var_neg = var_neg/chi_neg;
00827
00828 mean_pos_2 << m_sigma << " " << mean_pos << std::endl;
00829 mean_neg_2 << m_sigma << " " << mean_neg << std::endl;
00830 var_pos_2 << m_sigma << " " << var_pos << std::endl;
00831 var_neg_2 << m_sigma << " " << var_neg << std::endl;
00832
00833 chi_pos_1.close();
00834 mean_pos_no_norm_1.close();
00835 mean_neg_no_norm_1.close();
00836 mean_pos_1.close();
00837 mean_neg_1.close();
00838 var_pos_1.close();
00839 var_neg_1.close();
00840 chi_pos_2.close();
00841 chi_neg_2.close();
00842 mean_pos_no_norm_2.close();
00843 mean_neg_no_norm_2.close();
00844 mean_pos_2.close();
00845 mean_neg_2.close();
00846 var_pos_2.close();
00847 var_neg_2.close();
00848 }
00849 #endif
00850 }
00851 }
00852
00853 #endif // LOCALSTATSREGIONMODULE_H