/home/oan/prosjekt/gotools/segmentation/gpl_distro/lsseg_1.0_gpl/src/Filters.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: Filters.C                                                           
00037 //                                                                           
00038 // Created: Wed Feb 22 13:34:51 2006                                         
00039 //                                                                           
00040 // Author: Odd A. Andersen <Odd.Andersen@sintef.no>
00041 //                                                                           
00042 // Revision: $Id: Filters.C,v 1.12 2006/11/25 20:08:26 oan Exp $
00043 //                                                                           
00044 // Description:
00047 //                                                                           
00048 //===========================================================================
00049 
00050 #include <cmath>
00051 #include <stdexcept>
00052 #include <assert.h>
00053 #include "EigValComp3x3.h"
00054 #include "cimg_dependent.h"
00055 #include "Filters.h"
00056 #include "DiscreteApproximations.h"
00057 
00058 using namespace std;
00059 
00060 namespace {
00061 
00062 void diffusion_1D(const double* const u_old, 
00063                   const double* const g,
00064                   double * u_new, 
00065                   double dt,
00066                   int dim);
00067 void tridiag(const double* const a, 
00068              const double* const b, 
00069              const double* const c,
00070              const double* const r, 
00071              double* const sol, int n);
00072 
00073 void accumulate_scale(const lsseg::Image<double>& img1,
00074                       const lsseg::Image<double>& img2,
00075                       lsseg::Image<double>& inv_scale, 
00076                       lsseg::Image<double>& time_span,
00077                       double dt);
00078 
00079 inline double p_func(double u_squared, double p);
00080 
00081 }; // end anonymous namespace 
00082 
00083 
00084 namespace lsseg {
00085 
00086 //===========================================================================
00087 void compute_structure_tensor_2D(const Image<double>& img, 
00088                                  Image<double>& G, 
00089                                  bool square_root)
00090 //===========================================================================
00091 {
00092     assert(img.dimz() == 1);
00093     G.resize(img.dimx(), img.dimy(), 1, 3);
00094     G = 0;
00095     const int X = img.dimx();
00096     const int Y = img.dimy();
00097     for (int ch = 0; ch != img.numChannels(); ++ch) {
00098         for (int y = 0; y < Y; ++y) {
00099             const int yp = (y > 0) ? y-1 : y;
00100             const int yn = (y < Y-1) ? y+1 : y;
00101             const double dy_inv = (yn-yp == 2) ? 0.5 : 1;
00102 
00103             for (int x = 0; x < X; ++x) {
00104                 const int xp = (x > 0) ? x-1 : x;
00105                 const int xn = (x < X-1) ? x+1 : x;
00106                 const double dx_inv = (xn-xp == 2) ? 0.5 : 1;
00107 
00108                 const double Ipc = img(xp, y, 0, ch);
00109                 const double Inc = img(xn, y, 0, ch);
00110                 const double Icp = img(x, yp, 0, ch);
00111                 const double Icn = img(x, yn, 0, ch);
00112                 const double ix = dx_inv * (Inc - Ipc);
00113                 const double iy = dy_inv * (Icn - Icp);
00114                 
00115                 double norm_inv = 1;
00116                 if (square_root) {
00117                     const double norm = sqrt(ix * ix + iy * iy);
00118                     norm_inv = (norm > 0) ? double(1) / norm : 1;
00119                 }
00120                 G(x, y, 0, 0) += ix * ix * norm_inv;
00121                 G(x, y, 0, 1) += ix * iy * norm_inv;
00122                 G(x, y, 0, 2) += iy * iy * norm_inv;
00123             }
00124         }
00125     }
00126 }
00127 
00128 //===========================================================================
00129 void compute_structure_tensor_3D(const Image<double>& img,
00130                                  Image<double>& G,
00131                                  bool square_root)
00132 //===========================================================================
00133 {
00134     G.resize(img.dimx(), img.dimy(), img.dimz(), 6);
00135     G = 0;
00136     const int X = img.dimx();
00137     const int Y = img.dimy();
00138     const int Z = img.dimz();
00139     for (int ch = 0; ch != img.numChannels(); ++ch) {
00140         for (int z = 0; z < Z; ++z) {
00141             const int zp = (z > 0) ? z-1 : z;
00142             const int zn = (z < Z-1) ? z+1 : z;
00143             const double dz_inv = (zn - zp == 2) ? 0.5 : 1;
00144 
00145             for (int y = 0; y < Y; ++y) {
00146                 const int yp = (y > 0) ? y-1 : y;
00147                 const int yn = (y < Y-1) ? y+1 : y;
00148                 const double dy_inv = (yn - yp == 2) ? 0.5 : 1;
00149                 
00150                 for (int x = 0; x < X; ++x) {
00151                     const int xp = (x > 0) ? x-1 : x;
00152                     const int xn = (x < X-1) ? x+1 : x;
00153                     const double dx_inv = (xn-xp == 2) ? 0.5 : 1;
00154                     
00155                     const double Ipcc = img(xp, y, z, ch);
00156                     const double Incc = img(xn, y, z, ch);
00157                     const double Icpc = img(x, yp, z, ch);
00158                     const double Icnc = img(x, yn, z, ch);
00159                     const double Iccp = img(x, y, zp, ch);
00160                     const double Iccn = img(x, y, zn, ch);
00161                     const double ix = dx_inv * (Incc - Ipcc);
00162                     const double iy = dy_inv * (Icnc - Icpc);
00163                     const double iz = dz_inv * (Iccn - Iccp);
00164 
00165                     double norm_inv = 1;
00166                     if(square_root) {
00167                         const double norm = sqrt(ix * ix + iy * iy + iz * iz);
00168                         norm_inv = (norm > 0) ? double(1)/ norm : 1;
00169                     }
00170                     G(x, y, z, 0) = ix * ix * norm_inv;  // x x
00171                     G(x, y, z, 1) = ix * iy * norm_inv;  // x y
00172                     G(x, y, z, 2) = iy * iy * norm_inv;  // y y
00173                     G(x, y, z, 3) = ix * iz * norm_inv;  // x z
00174                     G(x, y, z, 4) = iy * iz * norm_inv;  // y z
00175                     G(x, y, z, 5) = iz * iz * norm_inv;  // z z
00176                 }
00177             }
00178         }
00179     }
00180 }
00181 
00182 
00183 //===========================================================================
00184 void compute_scale_factor_2D(const Image<double>& img,
00185                              Image<double>& scale_accum,
00186                              Image<double>& time_accum,
00187                              double dt,
00188                              double T)
00189 //===========================================================================
00190 {
00191     assert(img.numChannels() == 1);
00192     assert(img.dimz() == 1);
00193     // applying isotropic, nonlinear smoothing on image, and measuring the 
00194     // rate of change
00195     
00196     Image<double> tmp1(img); // copy content
00197     Image<double> tmp2(img, false); // do not copy content
00198     
00199     scale_accum.resize(img);
00200     scale_accum = 0;
00201     time_accum.resize(img);
00202     time_accum = 0;
00203     for (double t = 0; t < T; t+= dt) {
00204         nonlinear_gauss_filter_2D(tmp1, tmp2, dt, 0, 1);
00205         accumulate_scale(tmp1, tmp2, scale_accum, time_accum, dt);
00206         tmp1.swap(tmp2);
00207     }
00208 }
00209 
00210 //===========================================================================
00211 void compute_scale_factor_3D(const Image<double>& img,
00212                              Image<double>& scale_accum,
00213                              Image<double>& time_accum,
00214                              double dt,
00215                              double T)
00216 //===========================================================================
00217 {
00218     assert(img.numChannels() == 1);
00219     // applying isotropic, nonlinear smoothing on image, and measuring the 
00220     // rate of change
00221     
00222     Image<double> tmp1(img); // copy content
00223     Image<double> tmp2(img, false); // do not copy content
00224     
00225     scale_accum.resize(img);
00226     scale_accum = 0;
00227     time_accum.resize(img);
00228     time_accum = 0;
00229     for (double t = 0; t < T; t+= dt) {
00230         nonlinear_gauss_filter_3D(tmp1, tmp2, dt, 0, 1);
00231         accumulate_scale(tmp1, tmp2, scale_accum, time_accum, dt);
00232         tmp1.swap(tmp2);
00233     }
00234 }
00235 
00236 //===========================================================================
00237 void nonlinear_gauss_filter_2D(const Image<double>& img_old,  
00238                                Image<double>& img_new,
00239                                double dt,
00240                                double sigma,
00241                                double p)
00242 //===========================================================================
00243 {
00244     assert(img_old.size_compatible(img_new));
00245     
00246     Image<double> img_cpy(img_old);
00247 
00248     if (sigma > 0) {
00249         blur_image(img_cpy, sigma);
00250     }
00251 
00252     // computing the squared gradient matrix
00253     Image<double> g; 
00254     compute_squared_gradient_sum_2D(img_cpy, g);
00255     for (double* d = g.begin(); d != g.end(); ++d) {
00256         *d = p_func(*d, p);
00257     }
00258 
00259     // temporary result image
00260     Image<double> tmp_res;
00261 
00262     int PERMXY[] = {1, 0, 2, 3};
00263     for (int run = 0; run < 2; ++run) {
00264         int len = img_cpy.dimx();
00265         double* gp = g.begin();
00266         double* trg = (run == 0) ? img_new.begin() : tmp_res.begin();
00267         for (double* cp = img_cpy.begin(); 
00268              cp != img_cpy.end(); cp += len, gp += len, trg += len) {
00269             if (gp == g.end()) { // reset pointer to g data
00270                 gp = g.begin();
00271             }
00272             diffusion_1D(cp, gp, trg, dt, len);
00273         }
00274         img_cpy.permute(PERMXY); // permuting x, y and z component
00275         g.permute(PERMXY);
00276         img_new.permute(PERMXY);
00277         if (run == 0) {
00278             tmp_res.resize(img_new.dimx(), 
00279                            img_new.dimy(), 
00280                            img_new.dimz(), 
00281                            img_new.numChannels());
00282         } else {
00283             tmp_res.permute(PERMXY);
00284             img_new += tmp_res;
00285         }
00286     }
00287     img_new /= 2;
00288 }
00289 
00290 
00291 //===========================================================================
00292 void nonlinear_gauss_filter_3D(const Image<double>& img_old,
00293                                Image<double>& img_new,
00294                                double dt,
00295                                double sigma,
00296                                double p)
00297 //===========================================================================
00298 {
00299     const int PERM[] = {1, 2, 0, 3};
00300     assert(img_old.size_compatible(img_new));
00301     
00302     Image<double> img_cpy(img_old);
00303 
00304     if (sigma > 0) {
00305         blur_image(img_cpy, sigma);
00306     }
00307 
00308     // computing the squared gradient matrix
00309     Image<double> g; 
00310     compute_squared_gradient_sum_3D(img_cpy, g);
00311     for (double* d = g.begin(); d != g.end(); ++d) {
00312         *d = p_func(*d, p);
00313     }
00314 
00315     // temporary result image
00316     Image<double> tmp_res;
00317 
00318     // treating image lines in first run, columns in the second loop and pillars in
00319     // third run
00320     for (int run = 0; run < 3; ++run) {
00321         int len = img_cpy.dimx();
00322         double* gp = g.begin();
00323         double* trg = (run == 0) ? img_new.begin() : tmp_res.begin();
00324         for (double* cp = img_cpy.begin(); 
00325              cp != img_cpy.end(); cp += len, gp += len, trg += len) {
00326             if (gp == g.end()) { // reset pointer to g data
00327                 gp = g.begin();
00328             }
00329             diffusion_1D(cp, gp, trg, dt, len);
00330         }
00331         img_cpy.permute(PERM); // permuting x, y and z component
00332         g.permute(PERM);
00333         img_new.permute(PERM);
00334         if (run == 0) {
00335             tmp_res.resize(img_new.dimx(), 
00336                            img_new.dimy(), 
00337                            img_new.dimz(), 
00338                            img_new.numChannels());
00339         } else {
00340             tmp_res.permute(PERM);
00341             img_new += tmp_res;
00342         }
00343     }
00344 
00345     img_new /= 3;
00346 }
00347 
00348 //===========================================================================
00349 void anisotropic_smoothing(const Image<double>& img_old, 
00350                            Image<double>& img_new, 
00351                            double dt, 
00352                            double sigma,
00353                            double p)
00354 //===========================================================================
00355 {
00356     assert(p >= 0 && sigma >= 0 && dt >= 0);
00357     img_new.resize(img_old);
00358     if (img_old.dimz() > 1) {
00359         nonlinear_gauss_filter_3D(img_old, img_new, dt, sigma, p);
00360     } else {
00361         nonlinear_gauss_filter_2D(img_old, img_new, dt, sigma, p);
00362     }
00363 }
00364 
00365 //===========================================================================
00366 void compute_smoothing_geometry_3D(const Image<double>& G,
00367                                    double p1,
00368                                    double p2,
00369                                    double p3,
00370                                    Image<double>& T,
00371                                    bool take_square_root)
00372 //===========================================================================
00373 {
00374     const double EPS = 1.0e-8; // @@ only used in assert below
00375     assert(G.numChannels() == 6);
00376     T.resize(G);
00377     const int X = G.dimx(), Y = G.dimy(), Z = G.dimz();
00378     if (take_square_root) {
00379         p1 *= 0.5;
00380         p2 *= 0.5;
00381         p3 *= 0.5;
00382     }
00383     vector<double> v1(3), v2(3),v3(3);
00384     double lambda1, lambda2, lambda3;
00385     for (int z = 0; z < Z; ++z) {
00386         for (int y = 0; y < Y; ++y) {
00387             for (int x = 0; x < X; ++x) {
00388                 // finding eigenvalues and vectors of current structure tensor
00389                 const double alpha1 = G(x, y, z, 0);  // dx*dx
00390                 const double alpha2 = G(x, y, z, 2);  // dy*dy
00391                 const double alpha3 = G(x, y, z, 5);  // dz*dz
00392                 const double beta = G(x, y, z, 1);    // dx*dy
00393                 const double gamma = G(x, y, z, 3);   // dx*dz
00394                 const double delta = G(x, y, z, 4);   // dy*dz
00395 
00396                 numeric_eigsys(alpha1, alpha2, alpha3, 
00397                                beta, gamma, delta, 
00398                                lambda1, lambda2, lambda3, 
00399                                &v1[0], &v2[0], &v3[0]);
00400 
00401                 // sorting eigen-elements
00402 
00403                 if (fabs(lambda3) > fabs(lambda2)) {
00404                     swap(lambda3, lambda2);
00405                     swap(v3, v2);
00406                 }
00407                 if (fabs(lambda2) > fabs(lambda1)) {
00408                     swap(lambda2, lambda1);
00409                     swap(v2, v1);
00410                     if (fabs(lambda3) > fabs(lambda2)) {
00411                         swap(lambda3, lambda2);
00412                         swap(v3, v2);
00413                     }
00414                 }
00415                 assert(fabs(lambda1) >= fabs(lambda2) && fabs(lambda2) >= fabs(lambda3));
00416 
00417                 // constructing T
00418                 assert(lambda1 >= -EPS && lambda2 >= -EPS && lambda3 >= -EPS);
00419                 const double tmp = 1 + lambda1 + lambda2 + lambda3;
00420                 const double fac_3 = double(1) / pow(tmp, p1); // biggest @@
00421                 const double fac_2 = double(1) / pow(tmp, p2); 
00422                 const double fac_1 = double(1) / pow(tmp, p3); // smallepst@@
00423 
00424                 assert(fac_3 >= fac_2 && fac_2 >= fac_1);
00425 
00426                 T(x, y, z, 0) = 
00427                     (fac_3 * v3[0] * v3[0]) + 
00428                     (fac_2 * v2[0] * v2[0]) + 
00429                     (fac_1 * v1[0] * v1[0]);
00430                 T(x, y, z, 1) = 
00431                     (fac_3 * v3[0] * v3[1]) + 
00432                     (fac_2 * v2[0] * v2[1]) + 
00433                     (fac_1 * v1[0] * v1[1]);
00434                 T(x, y, z, 2) = 
00435                     (fac_3 * v3[1] * v3[1]) + 
00436                     (fac_2 * v2[1] * v2[1]) + 
00437                     (fac_1 * v1[1] * v1[1]);
00438                 T(x, y, z, 3) = 
00439                     (fac_3 * v3[0] * v3[2]) + 
00440                     (fac_2 * v2[0] * v2[2]) + 
00441                     (fac_1 * v1[0] * v1[2]);
00442                 T(x, y, z, 4) = 
00443                     (fac_3 * v3[1] * v3[2]) + 
00444                     (fac_2 * v2[1] * v2[2]) + 
00445                     (fac_1 * v1[1] * v1[2]);
00446                 T(x, y, z, 5) = 
00447                     (fac_3 * v3[2] * v3[2]) + 
00448                     (fac_2 * v2[2] * v2[2]) + 
00449                     (fac_1 * v1[2] * v1[2]);
00450             }
00451         }
00452     }
00453 }
00454 
00455 
00456 //===========================================================================
00457 void compute_smoothing_geometry_2D(const Image<double>& G,
00458                                    double p1,
00459                                    double p2,
00460                                    Image<double>& T,
00461                                    bool take_square_root)
00462 //===========================================================================
00463 {
00464     const double EPS = 1.0e-8;  // used in computing eigenvectors
00465     assert(G.numChannels() == 3 && G.dimz() == 1);
00466     T.resize(G);
00467     const int X = G.dimx();
00468     const int Y = G.dimy();
00469 
00470     if (take_square_root) {
00471         p1 *= 0.5;
00472         p2 *= 0.5;
00473     }
00474     double lm, lp, xm, xp, ym, yp; // lambda plus/minus, x and y components of theta plus/minus
00475     for (int y = 0; y < Y; ++y) {
00476         for (int x = 0; x < X; ++x) {
00477             const double a = G(x, y, 0, 0);
00478             const double b = G(x, y, 0, 1);
00479             const double c = G(x, y, 0, 2);
00480 
00481             // computing eigenvalues/vecs
00482             if (fabs(b) < EPS) {
00483                 // we consider the off-diagonal values to be zero
00484                 if (a > c) { 
00485                     lp = a;  lm = c;
00486                     xp = 1;  xm = 0;
00487                     yp = 0;  ym = 1;
00488                 } else {
00489                     lp = c;  lm = a;
00490                     xp = 0;  xm = 1;
00491                     yp = 1;  ym = 0;
00492                 }
00493             } else {
00494                 const double amc = a - c;
00495                 const double apc = a + c;
00496                 const double b2 = b * b;
00497                 const double discr = sqrt(amc*amc + 4 * b2); // root of discriminant
00498                 assert(discr > 0); // should only happen if b=0, which is taken care of above
00499 
00500                 lp = 0.5 * (apc + discr);  
00501                 const double amlp = a - lp;
00502                 double norm_fac_p = sqrt(b2 + amlp * amlp);
00503                 assert(norm_fac_p > 0); // same reason as above
00504                 norm_fac_p = double(1) / norm_fac_p;
00505                 xp = -b * norm_fac_p;
00506                 yp = amlp * norm_fac_p;
00507 
00508                 lm = 0.5 * (apc - discr);
00509                 const double amlm = a - lm;
00510                 double norm_fac_m = sqrt(b2 + amlm * amlm);
00511                 assert(norm_fac_m > 0); // same reason as above
00512                 norm_fac_m = double(1) / norm_fac_m;
00513                 xm = -b * norm_fac_m;
00514                 ym = amlm * norm_fac_m;
00515             }
00516 
00517             // constructing T
00518             const double tmp = 1 + lp + lm;
00519             const double fac_m = double(1) / pow(tmp, p1);
00520             const double fac_p = double(1) / pow(tmp, p2); 
00521             T(x, y, 0, 0) = (fac_p * xp * xp) + (fac_m * xm * xm); // upper left value
00522             T(x, y, 0, 1) = (fac_p * xp * yp) + (fac_m * xm * ym); // symmetric off-diagonal value
00523             T(x, y, 0, 2) = (fac_p * yp * yp) + (fac_m * ym * ym); // lower right value 
00524  
00525         }
00526     }
00527 }
00528 
00529 
00530 
00531 }; // end namespace lsseg
00532 
00533 namespace {
00534 
00535 //===========================================================================
00536 void accumulate_scale(const lsseg::Image<double>& img1, 
00537                       const lsseg::Image<double>& img2, 
00538                       lsseg::Image<double>& inv_scale, 
00539                       lsseg::Image<double>& time_span,
00540                       double dt)
00541 //===========================================================================
00542 {
00543     const double EPS = 1.0e-5;
00544     assert(img1.size_compatible(img2));
00545     assert(img1.size_compatible(inv_scale));
00546     assert(img1.size_compatible(time_span));
00547     assert(img1.numChannels() == 1);
00548 
00549     const double* img1_it = img1.begin();
00550     const double* img2_it = img2.begin();
00551     double* scale_it = inv_scale.begin();
00552     double* time_it = time_span.begin();
00553 
00554     for ( ; img1_it != img1.end(); ++img1_it, ++img2_it, ++scale_it, ++time_it) {
00555         double norm = fabs(*img1_it - *img2_it);
00556 
00557         if (norm > EPS) {
00558             *scale_it += norm * dt;
00559             *time_it += 4 * dt;
00560         }
00561     }
00562 }
00563 
00564 //==============================================================================
00565 double p_func(double u_squared, double p)
00566 //==============================================================================
00567 {
00568     static const double EPS_SQUARED = 1.0e-4; // regularising constant
00569     
00570     if (fabs(p-1) < EPS_SQUARED) {
00571         return double(1)/sqrt(u_squared + EPS_SQUARED);
00572     } else {
00573         return double(1) / pow(u_squared + EPS_SQUARED, 0.5 * p);
00574     }
00575 }
00576 
00577 //==============================================================================
00578 void tridiag(const double* const a, 
00579              const double* const b, 
00580              const double* const c,
00581              const double* const r, 
00582              double* const sol, int n)
00583 //==============================================================================
00584 {
00585     double bet; // temporary variable
00586     vector<double> gam(n); // workspace
00587     if (b[0] == 0.0) throw runtime_error("Error 1 in tridiag.");
00588     sol[0] = r[0] / (bet = b[0]);
00589     for (int j = 1; j < n; ++j) {
00590         gam[j] = c[j-1] / bet;
00591         bet = b[j] - a[j] * gam[j];
00592         if (bet == 0.0) throw runtime_error("Error 2 in tridiag.");
00593         sol[j] = (r[j] - a[j] * sol[j-1]) / bet;
00594     }
00595     for (int j = (n-2); j >= 0; --j) {
00596         sol[j] -= gam[j+1] * sol[j+1];
00597     }
00598 }
00599 
00600 //===========================================================================
00601 void diffusion_1D(const double* const u_old, 
00602                   const double* const g,
00603                   double * u_new, 
00604                   double dt,
00605                   int dim)
00606 //===========================================================================
00607 {
00608     if (dim == 1) {
00609         // no diffusion possible, just copy input data
00610         *u_new = *u_old;
00611         return;
00612     }
00613     static vector<double> a, b, c;
00614     a.resize(dim);
00615     b.resize(dim);
00616     c.resize(dim);
00617     
00618     // border cases
00619     a[0] = 0;
00620     c[0] = -dt * (g[0] + g[1]);
00621     b[0] = 1 - c[0];
00622 
00623     a[dim-1] = -dt * (g[dim-2] + g[dim-1]);
00624     c[dim-1] = 0;
00625     b[dim-1] = 1 - a[dim-1];
00626 
00627     // interior
00628     for (int i = 1; i < dim-1; ++i) {
00629         a[i] = c[i-1];
00630         c[i] = -dt * (g[i] + g[i+1]);// * 0.5;
00631         b[i] = 1 - (a[i] + c[i]);
00632     }
00633 
00634     // solving linear system
00635     tridiag(&a[0], &b[0], &c[0], u_old, u_new, dim);
00636 
00637 }
00638 
00639 
00640 }; // end anonymous namespace 
00641 
00642 
00643 
00644 
00645 // //===========================================================================
00646 // void compute_scale_factor(const Image& img,
00647 //                        Image& scale_accum,
00648 //                        Image& time_accum,
00649 //                        double dt,
00650 //                        double T)
00651 // //===========================================================================
00652 // {
00653 //     assert(img.dim == 1);
00654 //     // applying isotropic, nonlinear smoothing on image, and measuring the
00655 //     // rate of change
00656 
00657 //     Image tmp1(img); // copy content 
00658 //     Image tmp2(img, false); // do not copy content
00659 //     scale_accum.resize(img);
00660 //     scale_accum.fill(0);
00661 //     time_accum.resize(img);
00662 //     time_accum.fill(0);
00663 //     for (double t = 0; t < T; t+= dt) {
00664 //      cout << "Time: " << t << endl;
00665 //      nonlinear_gauss_filter(tmp1, tmp2, dt, 0, 1);
00666 //      accumulate_scale(tmp1, tmp2, scale_accum, time_accum, dt);
00667 //      tmp1.swap(tmp2);
00668 //     }
00669 // }
00670 
00671 
00672 
00673 // //===========================================================================
00674 // void compute_structure_tensor(const Image& img, Image& G, bool square_root)
00675 // //===========================================================================
00676 // {
00677 //     assert(img.dim == 1); // if not, offset will not work
00678 //     G.resize(img.width, img.height, 3);
00679 //     fill(G.data.begin(), G.data.end(), 0);
00680     
00681 //     double Ipc, Inc, Icp, Icn;
00682 //     if (square_root) {
00683 //      int N = img.width;
00684 //      int offset = 0;
00685 //      for (int y = 0; y < img.height; ++y) {
00686 //          for (int x = 0; x < img.width; ++x, ++offset) {
00687 //              Ipc = (x>0) ? img.data[offset-1] : 2 * img.data[offset] - img.data[offset+1];
00688 //              Inc = (x<img.width-1) ? img.data[offset+1] : 2*img.data[offset]-img.data[offset-1];
00689 //              Icp = (y>0) ? img.data[offset - N] : 2 * img.data[offset] - img.data[offset+N];
00690 //              Icn = (y<img.height-1) ? img.data[offset+N] : 2*img.data[offset]-img.data[offset-N];
00691                 
00692 //              double ix = 0.5 * (Inc - Ipc);
00693 //              double iy = 0.5 * (Icn - Icp);
00694 //              double norm = sqrt(ix * ix + iy * iy);
00695 //              if (norm > 0) {
00696 //                  double norm_inv = double(1) / norm;
00697 //                  G.data[offset + 0 * img.size()] = ix * ix * norm_inv;
00698 //                  G.data[offset + 1 * img.size()] = ix * iy * norm_inv;
00699 //                  G.data[offset + 2 * img.size()] = iy * iy * norm_inv;
00700 //              } else {
00701 //                  G.data[offset + 0 * img.size()] = 0;
00702 //                  G.data[offset + 1 * img.size()] = 0;
00703 //                  G.data[offset + 2 * img.size()] = 0;
00704 //              }
00705 //          }
00706 //      }
00707 //     } else {
00708 //      int N = img.width;
00709 //      int offset = 0;
00710 //      for (int y = 0; y < img.height; ++y) {
00711 //          for (int x = 0; x < img.width; ++x, ++offset) {
00712 //              Ipc = (x>0) ? img.data[offset-1] : 2 * img.data[offset] - img.data[offset+1];
00713 //              Inc = (x<img.width-1) ? img.data[offset+1] : 2*img.data[offset]-img.data[offset-1];
00714 //              Icp = (y>0) ? img.data[offset - N] : 2 * img.data[offset] - img.data[offset+N];
00715 //              Icn = (y<img.height-1) ? img.data[offset+N] : 2*img.data[offset]-img.data[offset-N];
00716                 
00717 //              double ix = 0.5 * (Inc - Ipc);
00718 //              double iy = 0.5 * (Icn - Icp);
00719 //              G.data[offset + 0 * img.size()] = ix * ix;
00720 //              G.data[offset + 1 * img.size()] = ix * iy;
00721 //              G.data[offset + 2 * img.size()] = iy * iy;
00722 //          }
00723 //      }
00724 //     }
00725 // }
00726 

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