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 "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"
00059
00060 namespace {
00061
00062 void region_compete(const std::vector<lsseg::LevelSetFunction>& cterms,
00063 std::vector<lsseg::LevelSetFunction>& fterms,
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 };
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
00106 assert(num_regs > 0 && num_iter >= 0);
00107 assert(regs[0].phi.dimz() == 1);
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));
00111 }
00112
00113
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
00146
00147
00148
00149 regs[r].phi.curvature2D(curv_terms[r], &(region_masks[r]));
00150 curv_terms[r] *= regs[r].mu;
00151
00152
00153
00154 regs[r].fgen->update(regs[r].phi);
00155 regs[r].fgen->force(force_terms[r], &(region_masks[r]));
00156 }
00157
00158
00159 region_compete(curv_terms, force_terms, region_masks, regs);
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169 for (int r = 0; r < num_regs; ++r) {
00170 double H1, H2;
00171
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
00180
00181
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(®ion_masks[r]);
00202
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
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));
00223 }
00224
00225
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
00261 regs[r].phi.curvature3D(curv_terms[r], &(region_masks[r]));
00262 curv_terms[r] *= regs[r].mu;
00263
00264
00265 regs[r].fgen->update(regs[r].phi);
00266 regs[r].fgen->force(force_terms[r], &(region_masks[r]));
00267 }
00268
00269
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
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
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(®ion_masks[r]);
00306 }
00307 }
00308 cout << "dt: " << dt << endl;
00309 }
00310 return time_accum;
00311 }
00312
00313
00314 };
00315
00316 namespace {
00317
00318
00319 void region_compete(const std::vector<lsseg::LevelSetFunction>& cterms,
00320 std::vector<lsseg::LevelSetFunction>& fterms,
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];
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
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
00374 int ix_h = -1;
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
00395 vals[i] -= nhsf;
00396 assert(vals[i] >= 0);
00397 }
00398 }
00399 }
00400
00401
00402
00403 };