00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
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 };
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;
00171 G(x, y, z, 1) = ix * iy * norm_inv;
00172 G(x, y, z, 2) = iy * iy * norm_inv;
00173 G(x, y, z, 3) = ix * iz * norm_inv;
00174 G(x, y, z, 4) = iy * iz * norm_inv;
00175 G(x, y, z, 5) = iz * iz * norm_inv;
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
00194
00195
00196 Image<double> tmp1(img);
00197 Image<double> tmp2(img, false);
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
00220
00221
00222 Image<double> tmp1(img);
00223 Image<double> tmp2(img, false);
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
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
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()) {
00270 gp = g.begin();
00271 }
00272 diffusion_1D(cp, gp, trg, dt, len);
00273 }
00274 img_cpy.permute(PERMXY);
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
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
00316 Image<double> tmp_res;
00317
00318
00319
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()) {
00327 gp = g.begin();
00328 }
00329 diffusion_1D(cp, gp, trg, dt, len);
00330 }
00331 img_cpy.permute(PERM);
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;
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
00389 const double alpha1 = G(x, y, z, 0);
00390 const double alpha2 = G(x, y, z, 2);
00391 const double alpha3 = G(x, y, z, 5);
00392 const double beta = G(x, y, z, 1);
00393 const double gamma = G(x, y, z, 3);
00394 const double delta = G(x, y, z, 4);
00395
00396 numeric_eigsys(alpha1, alpha2, alpha3,
00397 beta, gamma, delta,
00398 lambda1, lambda2, lambda3,
00399 &v1[0], &v2[0], &v3[0]);
00400
00401
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
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);
00421 const double fac_2 = double(1) / pow(tmp, p2);
00422 const double fac_1 = double(1) / pow(tmp, p3);
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;
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;
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
00482 if (fabs(b) < EPS) {
00483
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);
00498 assert(discr > 0);
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);
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);
00512 norm_fac_m = double(1) / norm_fac_m;
00513 xm = -b * norm_fac_m;
00514 ym = amlm * norm_fac_m;
00515 }
00516
00517
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);
00522 T(x, y, 0, 1) = (fac_p * xp * yp) + (fac_m * xm * ym);
00523 T(x, y, 0, 2) = (fac_p * yp * yp) + (fac_m * ym * ym);
00524
00525 }
00526 }
00527 }
00528
00529
00530
00531 };
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;
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;
00586 vector<double> gam(n);
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
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
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
00628 for (int i = 1; i < dim-1; ++i) {
00629 a[i] = c[i-1];
00630 c[i] = -dt * (g[i] + g[i+1]);
00631 b[i] = 1 - (a[i] + c[i]);
00632 }
00633
00634
00635 tridiag(&a[0], &b[0], &c[0], u_old, u_new, dim);
00636
00637 }
00638
00639
00640 };
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665
00666
00667
00668
00669
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697
00698
00699
00700
00701
00702
00703
00704
00705
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726