/home/oan/prosjekt/gotools/segmentation/gpl_distro/lsseg_1.0_gpl/src/MultiRegionAlgorithm.C

Go to the documentation of this file.
00001 //===========================================================================
00002 // The Level-Set Segmentation Library (LSSEG)
00003 //
00004 //
00005 // Copyright (C) 2000-2005 SINTEF ICT, Applied Mathematics, Norway.
00006 //
00007 // This program is free software; you can redistribute it and/or          
00008 // modify it under the terms of the GNU General Public License            
00009 // as published by the Free Software Foundation version 2 of the License. 
00010 //
00011 // This program is distributed in the hope that it will be useful,        
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of         
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          
00014 // GNU General Public License for more details.                           
00015 //
00016 // You should have received a copy of the GNU General Public License      
00017 // along with this program; if not, write to the Free Software            
00018 // Foundation, Inc.,                                                      
00019 // 59 Temple Place - Suite 330,                                           
00020 // Boston, MA  02111-1307, USA.                                           
00021 //
00022 // Contact information: e-mail: tor.dokken@sintef.no                      
00023 // SINTEF ICT, Department of Applied Mathematics,                         
00024 // P.O. Box 124 Blindern,                                                 
00025 // 0314 Oslo, Norway.                                                     
00026 // 
00027 //
00028 // Other licenses are also available for this software, notably licenses
00029 // for:
00030 // - Building commercial software.                                        
00031 // - Building software whose source code you wish to keep private.        
00032 //
00033 //===========================================================================
00034 //===========================================================================
00035 //                                                                           
00036 // File: MultiRegionAlgorithm.C                                              
00037 //                                                                           
00038 // Created: Fri Mar  3 11:13:18 2006                                         
00039 //                                                                           
00040 // Author: Odd A. Andersen <Odd.Andersen@sintef.no>
00041 //                                                                           
00042 // Revision: $Id: MultiRegionAlgorithm.C,v 1.15 2006/11/25 20:08:26 oan Exp $
00043 //                                                                           
00044 // Description:
00047 //                                                                           
00048 //===========================================================================
00049 
00050 #include "BorderCondition.h"
00051 #include <assert.h>
00052 #include <limits>
00053 #include <time.h>
00054 #include "level_set.h"
00055 #include "Mask.h"
00056 #include "simple_tools.h"
00057 #include "Region.h"
00058 #include "cimg_dependent.h" // for debug purposes
00059 
00060 namespace { // anonymous namespace
00061 
00062 void region_compete(const std::vector<lsseg::LevelSetFunction>& cterms,
00063                     std::vector<lsseg::LevelSetFunction>& fterms, // will be updated
00064                     const std::vector<lsseg::Mask>& masks,
00065                     const lsseg::Region* regs);
00066 
00067 void resolve_competing(double* vals, int num_vals);
00068 
00069 inline double smooth_dirac(double val)
00070 {
00071     const double FAC = 3;
00072     val *= FAC;
00073     return double(1) / (1 + val * val);
00074 }
00075 
00076 inline bool regional_overlap(size_t k, const lsseg::Region* regs, int num_regs)
00077 {
00078     bool one_region_found = false;
00079     for (int r = 0; r < num_regs; ++r) {
00080         if (regs[r].phi[k] <= 0) {
00081             if (one_region_found) {
00082                 return true;
00083             } else {
00084                 one_region_found = true;
00085             }
00086         }
00087     }
00088     return false;
00089 }
00090 
00091 }; // end anonymous namespace
00092 
00093 using namespace std;
00094 
00095 namespace lsseg {
00096 
00097 //===========================================================================
00098 double develop_multiregion_2D(Region* regs, 
00099                               int num_regs,
00100                               int num_iter, 
00101                               int reinit_modulo,
00102                               const Mask* geom_mask)
00103 //===========================================================================
00104 {
00105     // initial assertions
00106     assert(num_regs > 0 && num_iter >= 0);
00107     assert(regs[0].phi.dimz() == 1); // function written for the 2D case
00108     assert(!geom_mask || regs[0].phi.size_compatible(*geom_mask));
00109     for (int i = 1; i < num_regs; ++i) {
00110         assert(regs[i].phi.size_compatible(regs[0].phi)); // all regions should be of same size
00111     }
00112 
00113     // constants that are candidates to be user-defined
00114     const int MASK_WIDTH = 2;
00115     const BorderCondition bcond = NEUMANN;
00116     const int Y = regs[0].phi.dimy();
00117     const int X = regs[0].phi.dimx();
00118 
00119     std::vector<LevelSetFunction> curv_terms(num_regs);
00120     std::vector<LevelSetFunction> force_terms(num_regs); 
00121     std::vector<Mask> region_masks(num_regs);
00122     std::vector<Mask> old_region_masks(num_regs);
00123     std::vector<double> max_allowed_timestep(num_regs);
00124     double time_accum = 0;
00125     for (int r = 0; r < num_regs; ++r) {
00126         curv_terms[r].resize(X, Y);
00127         force_terms[r].resize(X, Y);
00128         old_region_masks[r].resize(X, Y);
00129         if (geom_mask) {
00130             region_masks[r] = *geom_mask;
00131         } else {
00132             region_masks[r].resize(X, Y);
00133             region_masks[r] = 1;
00134         }
00135     }
00136     
00137 
00138     for (int i = 0; i < num_iter; ++i) {
00139         for (int r = 0; r < num_regs; ++r) {
00140             region_masks[r].swap(old_region_masks[r]);
00141             make_border_mask(regs[r].phi, region_masks[r], MASK_WIDTH, &old_region_masks[r]);
00142             if (geom_mask) {
00143                 region_masks[r].intersect(*geom_mask);
00144             }
00145             //region_masks[r] = 1; //@@ KRULL
00146             //display_image(region_masks[r]);
00147 
00148             // compute curvature-driven force
00149             regs[r].phi.curvature2D(curv_terms[r], &(region_masks[r]));
00150             curv_terms[r] *= regs[r].mu;
00151             //display_image(curv_terms[r]);
00152 
00153             // compute individual normal force term
00154             regs[r].fgen->update(regs[r].phi);
00155             regs[r].fgen->force(force_terms[r], &(region_masks[r]));
00156         }    
00157 
00158         // let regions compete for influencence
00159         region_compete(curv_terms, force_terms, region_masks, regs); 
00160 
00161         //display_image(force_terms[0]);
00162         //display_image(force_terms[1]);
00163 //      display_image(force_terms[2]);
00164 //      Image<double> diff;
00165 //      diff = force_terms[0];
00166 //      diff -= force_terms[1];
00167 //      display_image(diff);
00168 
00169         for (int r = 0; r < num_regs; ++r) {
00170             double H1, H2;
00171             //display_image(force_terms[r]);
00172             normal_direction_flow_2D(regs[r].phi,
00173                                      force_terms[r],
00174                                      bcond,
00175                                      H1,
00176                                      H2,
00177                                      &(region_masks[r]),
00178                                      geom_mask); 
00179             // computing maximum allowed time step
00180             //cout << "H1 and H2: " << H1 << " " << H2 << endl;
00181             //cout << "Mu: " << regs[r].mu << endl;
00182             max_allowed_timestep[r] = double(1) / (4 * regs[r].mu + H1 + H2);
00183         }
00184 
00185         double dt = *min_element(max_allowed_timestep.begin(), max_allowed_timestep.end());
00186         time_accum += dt;
00187         LevelSetFunction grad;
00188         for (int r = 0; r < num_regs; ++r) {
00189             regs[r].phi.gradientNorm2D(grad, &(region_masks[r]));
00190             for (int y = 0; y < Y; ++y) {
00191                 for (int x = 0; x < X; ++x) {
00192                     if (region_masks[r](x, y) ) {
00193                         const double gradphi = grad(x, y);
00194                         const double uval = 
00195                             (gradphi * curv_terms[r](x, y)) + force_terms[r](x, y); 
00196                         regs[r].phi(x, y) += dt * uval;
00197                     }
00198                 }
00199             }
00200             if (reinit_modulo && ((i+1) % reinit_modulo == 0)) {
00201                 regs[r].phi.reinitialize2D(&region_masks[r]);
00202                 //regs[r].phi.reinitialize2D();
00203             }
00204         }
00205         cout << "dt: " << dt << endl;
00206     }
00207     return time_accum;
00208 }
00209 
00210 //===========================================================================
00211 double develop_multiregion_3D(Region* regs, 
00212                               int num_regs,
00213                               int num_iter, 
00214                               int reinit_modulo,
00215                               const Mask* geom_mask)
00216 //===========================================================================
00217 {
00218     // initial assertions
00219     assert(num_regs > 0 && num_iter >= 0);
00220     assert(!geom_mask || regs[0].phi.size_compatible(*geom_mask));
00221     for (int i = 1; i < num_regs; ++i) {
00222         assert(regs[i].phi.size_compatible(regs[0].phi)); // all regions should be of same size
00223     }
00224 
00225     // constants that are candidates to be user-defined
00226     const int MASK_WIDTH = 2;
00227     const BorderCondition bcond = NEUMANN;
00228     const int Z = regs[0].phi.dimz();
00229     const int Y = regs[0].phi.dimy();
00230     const int X = regs[0].phi.dimx();
00231 
00232     std::vector<LevelSetFunction> curv_terms(num_regs);
00233     std::vector<LevelSetFunction> force_terms(num_regs); 
00234     std::vector<Mask> region_masks(num_regs);
00235     std::vector<Mask> old_region_masks(num_regs);
00236     std::vector<double> max_allowed_timestep(num_regs);
00237     double time_accum = 0;
00238     for (int r = 0; r < num_regs; ++r) {
00239         curv_terms[r].resize(X, Y, Z);
00240         force_terms[r].resize(X, Y, Z);
00241         region_masks[r].resize(X, Y, Z);
00242         old_region_masks[r].resize(X, Y, Z);
00243         if (geom_mask) {
00244             region_masks[r] = *geom_mask;
00245         } else {
00246             region_masks[r].resize(X, Y, Z);
00247             region_masks[r] = 1;
00248         }
00249 
00250     }
00251 
00252     for (int i = 0; i < num_iter; ++i) {
00253         for (int r = 0; r < num_regs; ++r) {
00254             region_masks[r].swap(old_region_masks[r]);
00255             make_border_mask(regs[r].phi, region_masks[r], MASK_WIDTH, &old_region_masks[r]);
00256             if (geom_mask) {
00257                 region_masks[r].intersect(*geom_mask);
00258             }
00259 
00260             // compute curvature-driven force
00261             regs[r].phi.curvature3D(curv_terms[r], &(region_masks[r]));
00262             curv_terms[r] *= regs[r].mu;
00263 
00264             // compute individual normal force term
00265             regs[r].fgen->update(regs[r].phi);
00266             regs[r].fgen->force(force_terms[r], &(region_masks[r]));
00267         }    
00268 
00269         // let regions compete for influencence
00270         region_compete(curv_terms, force_terms, region_masks, regs); 
00271 
00272         for (int r = 0; r < num_regs; ++r) {
00273             double H1, H2, H3;
00274             //display_image(force_terms[r]);
00275             normal_direction_flow_3D(regs[r].phi,
00276                                      force_terms[r],
00277                                      bcond,
00278                                      H1,
00279                                      H2,
00280                                      H3,
00281                                      &(region_masks[r]),
00282                                      geom_mask); 
00283             // computing maximum allowed time step
00284             max_allowed_timestep[r] = double(1) / (4 * regs[r].mu + H1 + H2 + H3);
00285         }
00286 
00287         double dt = *min_element(max_allowed_timestep.begin(), max_allowed_timestep.end());
00288         time_accum += dt;
00289         LevelSetFunction grad;
00290         for (int r = 0; r < num_regs; ++r) {
00291             regs[r].phi.gradientNorm3D(grad, &(region_masks[r]));
00292             for (int z= 0; z < Z; ++z) {
00293                 for (int y = 0; y < Y; ++y) {
00294                     for (int x = 0; x < X; ++x) {
00295                         if (region_masks[r](x, y, z) ) {
00296                             const double gradphi = grad(x, y, z);
00297                             const double uval = 
00298                                 (gradphi * curv_terms[r](x, y, z)) + force_terms[r](x, y, z); 
00299                             regs[r].phi(x, y, z) += dt * uval;
00300                         }
00301                     }
00302                 }
00303             }
00304             if (reinit_modulo && ((i+1) % reinit_modulo == 0)) {
00305                 regs[r].phi.reinitialize3D(&region_masks[r]);
00306             }
00307         }
00308         cout << "dt: " << dt << endl;
00309     }
00310     return time_accum;
00311 }
00312 
00313 
00314 }; // end namespace lsseg
00315 
00316 namespace {
00317 
00318 //===========================================================================
00319 void region_compete(const std::vector<lsseg::LevelSetFunction>& cterms,
00320                     std::vector<lsseg::LevelSetFunction>& fterms, // will be updated
00321                     const std::vector<lsseg::Mask>& masks,
00322                     const lsseg::Region* regs)
00323 //===========================================================================
00324 {
00325     const double MIN_INFTY = -std::numeric_limits<double>::infinity();
00326     const size_t num_regs = cterms.size();
00327     assert(fterms.size() == num_regs && masks.size() == num_regs);
00328     const size_t REG_SIZE = fterms[0].size();
00329     const double OUT_CT = 1;
00330     std::vector<double> competing_vals(num_regs);
00331 
00332     for (size_t k = 0; k < REG_SIZE; ++k) {
00333 
00334         int overlaps = 0;
00335         for (size_t r = 0; r < num_regs; ++r) {
00336             if (masks[r][k]) {
00337                 ++overlaps;
00338             }
00339         }
00340 
00341         if (overlaps > 1) {
00342             for (size_t r = 0; r < num_regs; ++r) {
00343                 competing_vals[r] =
00344                     (masks[r][k]) ? fterms[r][k] - cterms[r][k] : MIN_INFTY;
00345             }
00346             resolve_competing(&competing_vals[0], num_regs);
00347             
00348             for (size_t r = 0; r < num_regs; ++r) {
00349                 if (masks[r][k]) {
00350                     fterms[r][k] = competing_vals[r] + cterms[r][k]; // @@ ?? is this correct?
00351                 }
00352             }
00353         } else {
00354             for (size_t r = 0; r < num_regs; ++r) {
00355                 if (masks[r][k] && !regional_overlap(k, regs, num_regs)) {
00356                     if (fterms[r][k] <= 0) {
00357                         // empty area will "suck" the function into it.
00358                         fterms[r][k] = OUT_CT; 
00359                     }
00360                 }
00361             }
00362         }
00363     }
00364 }
00365 
00366 //===========================================================================
00367 void resolve_competing(double* vals, int num_vals)
00368 //===========================================================================
00369 {
00370     const double MIN_VAL = -std::numeric_limits<double>::infinity();
00371     assert(num_vals >= 2);
00372 
00373     // finding highest and next highest value, as well as their indexes;
00374     int ix_h = -1; // index of highest value
00375     double hsf = MIN_VAL;
00376     double nhsf = hsf;
00377     for (int i = 0; i < num_vals; ++i) {
00378         if (vals[i] > hsf) {
00379             nhsf = hsf;
00380             hsf = vals[i];
00381             ix_h = i;
00382         } else if (vals[i] > nhsf) {
00383             nhsf = vals[i];
00384         }
00385     }
00386     assert(ix_h >= 0);
00387     assert(nhsf > MIN_VAL);
00388     assert(hsf > MIN_VAL);
00389     for (int i = 0; i < num_vals; ++i) {
00390         if (i != ix_h) {
00391             vals[i] -= hsf;
00392             assert(vals[i] <= 0);
00393         } else {
00394             // highest value
00395             vals[i] -= nhsf;
00396             assert(vals[i] >= 0);
00397         }
00398     }
00399 }
00400 
00401 
00402 
00403 }; // end anonymous namespace

Generated on Tue Nov 28 18:35:47 2006 for lsseg by  doxygen 1.4.7