From 91a222a3c3d073552ce84e8a54ca3759d6537a7d Mon Sep 17 00:00:00 2001 From: Guilhem Lavaux Date: Mon, 1 Apr 2013 21:30:45 -0500 Subject: [PATCH] Re-Implemented jozov2 using priority_queues (far less complicated and maybe better performing) --- c_tools/zobov2/jozov2.cpp | 4 +- c_tools/zobov2/jozov2_watershed.cpp | 244 ++++++++++++++-------------- 2 files changed, 124 insertions(+), 124 deletions(-) diff --git a/c_tools/zobov2/jozov2.cpp b/c_tools/zobov2/jozov2.cpp index 718b8e9..b33a906 100644 --- a/c_tools/zobov2/jozov2.cpp +++ b/c_tools/zobov2/jozov2.cpp @@ -139,10 +139,10 @@ int main(int argc,char **argv) if (z[i].np == 1) continue; - txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e") + txt << format("%d %d %d %e %e %d %d %e %d %f %6.2e %f") % (h+1) % i % z[i].core % p[z[i].core].dens % z[i].vol % z[i].np % z[i].nhl % z[i].voljoin % z[i].npjoin - % z[i].denscontrast % prob << endl; + % z[i].denscontrast % prob % z[i].leak << endl; } /* h+1 to start from 1, not zero */ txt.close(); diff --git a/c_tools/zobov2/jozov2_watershed.cpp b/c_tools/zobov2/jozov2_watershed.cpp index d6ef31d..3c46e97 100644 --- a/c_tools/zobov2/jozov2_watershed.cpp +++ b/c_tools/zobov2/jozov2_watershed.cpp @@ -2,6 +2,7 @@ #include #endif +#include #include #include #include @@ -13,6 +14,40 @@ using namespace std; using boost::format; +struct ZoneDensityPair +{ + int h; + double density; + + bool operator<(const ZoneDensityPair& p2) const + { + return density > p2.density; + } +}; + +typedef priority_queue ZoneQueue; + +static void build_process_queue(ZoneQueue& q, ZONE *z, char *inyet, int h) +{ + ZoneDensityPair zdp; + ZONE& z_h = z[h]; + bool interior = true; + + assert(inyet[h] == 1); + for (int za = 0; za < z_h.nadj; za++) + { + zdp.h = z_h.adj[za]; + zdp.density = z_h.slv[za]; + if (inyet[zdp.h] == 0) + { + q.push(zdp); + interior = false; + } + } + if (interior) + inyet[h] = 2; +} + void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol) { int *iord; @@ -28,36 +63,33 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo iord = new int[numZones]; - findrtop(sorter, numZones, iord, numZones); +// findrtop(sorter, numZones, iord, numZones); delete[] sorter; -#pragma omp parallel + //#pragma omp parallel { char *inyet, *inyet2; int *zonelist, *zonelist2; int nhl; - int *links = new int[NLINKS]; - bool *done_zones; int prev_ii = -1; inyet = new char[numZones]; inyet2 = new char[numZones]; zonelist = new int[numZones]; zonelist2 = new int[numZones]; - done_zones = new bool[numZones]; fill(inyet, inyet + numZones, 0); fill(inyet2, inyet2 + numZones, 0); - fill(done_zones, done_zones + numZones, false); nhl = 0; -#pragma omp for schedule(dynamic,1) + //#pragma omp for schedule(dynamic,1) for (int h = 0; h < numZones; h++) { int nhlcount = 0; - float lowvol; + float lowvol, z_cur_core_dens; bool beaten; - set to_process; + priority_queue to_process; + int link0; for (int hl = 0; hl < nhl; hl++) inyet[zonelist[hl]] = 0; @@ -65,144 +97,113 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo zonelist[0] = h; inyet[h] = 1; nhl = 1; - to_process.insert(h); z[h].npjoin = z[h].np; + z_cur_core_dens = p[z[h].core].dens; + + build_process_queue(to_process, z, inyet, h); + do { /* Find the lowest-volume (highest-density) adjacency */ int nl = 0; beaten = false; - lowvol = BIGFLT; - set::iterator iter = to_process.begin(); - while (iter != to_process.end()) - { - int h2 = *iter; - ZONE& z_cur = z[h2]; - bool interior = true, touched = false; - assert(inyet[h2] == 1); - - for (int za = 0; za < z_cur.nadj; za++) - { - if (inyet[z_cur.adj[za]] == 0) - { - interior = false; - if (z_cur.slv[za] == lowvol) - { - links[nl] = z_cur.adj[za]; - touched = true; - nl ++; - if (nl == NLINKS) - { - printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl); - exit(0); - } - } - else - if (z_cur.slv[za] < lowvol) - { - lowvol = z_cur.slv[za]; - links[0] = z_cur.adj[za]; - nl = 1; - touched = true; - } - } - } - ++iter; - if (interior || !touched) - to_process.erase(h2); - if (interior) - inyet[h2] = 2; /* No bordering exter. zones */ - } - - if (nl == 0) + if (to_process.empty()) { beaten = true; z[h].leak = maxvol; continue; } + + do + { + lowvol = to_process.top().density; + link0 = to_process.top().h; + to_process.pop(); + } + while ((inyet[link0] != 0) && (!to_process.empty())); + + if (to_process.empty()) + { + beaten = true; + z[h].leak = maxvol; + continue; + } + + nl++; - if (lowvol > voltol) + if (lowvol > voltol) { beaten = true; z[h].leak = lowvol; continue; } - for (int l = 0; l < nl; l++) - { - //if (!done_zones[links[l]]) /* Equivalent to p[z[links[l]].core].dens < p[z[h].core].dens) as zones are sorted. */ - if (p[z[links[l]].core].dens < p[z[h].core].dens) - beaten = true; - } - - if (beaten) + if (p[z[link0].core].dens < z_cur_core_dens) { + beaten = true; z[h].leak = lowvol; continue; } /* Add everything linked to the link(s) */ int nhl2 = 0; - for (int l = 0; l < nl; l++) - { - if (inyet2[links[l]] == 0) - { - zonelist2[nhl2] = links[l]; - inyet2[links[l]] = 1; - nhl2 ++; - bool added = true; - while (added && !beaten) - { - added = false; - for (int hl = 0; (hl < nhl2) && (!beaten); hl++) - { - int h2 = zonelist2[hl]; - - if (inyet2[h2] == 1) { - bool interior = true; /* Guilty until proven innocent */ - for (int za = 0; za < z[h2].nadj; za ++) { - int link2 = z[h2].adj[za]; - if ((inyet[link2]+inyet2[link2]) == 0) { - interior = false; - if (z[h2].slv[za] <= lowvol) { - //if (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens) - if (p[z[link2].core].dens < p[z[h].core].dens) { - beaten = true; - break; - } - zonelist2[nhl2] = link2; - inyet2[link2] = 1; - nhl2++; - added = true; - } - } - } - if (interior) - inyet2[h2] = 2; - } - } - } - } - } - - /* See if there's a beater */ - if (beaten) { - z[h].leak = lowvol; - } else { - for (int h2 = 0; h2 < nhl2; h2++) { - int new_h = zonelist2[h2]; - + zonelist2[0] = link0; + inyet2[link0] = 1; + nhl2=1; + + bool added = true; + while (added && !beaten) + { + added = false; + for (int hl = 0; (hl < nhl2) && (!beaten); hl++) + { + int h2 = zonelist2[hl]; + + if (inyet2[h2] == 1) { + bool interior = true; /* Guilty until proven innocent */ + + for (int za = 0; za < z[h2].nadj; za ++) { + int link2 = z[h2].adj[za]; + + if ((inyet[link2]+inyet2[link2]) == 0) { + interior = false; + if (z[h2].slv[za] <= lowvol) { + if (p[z[link2].core].dens < z_cur_core_dens) { + beaten = true; + break; + } + zonelist2[nhl2] = link2; + inyet2[link2] = 1; + nhl2++; + added = true; + } + } + } + + if (interior) + inyet2[h2] = 2; + } + } + } + + /* See if there's a beater */ + if (beaten) { + z[h].leak = lowvol; + } else { + for (int h2 = 0; h2 < nhl2; h2++) { + int new_h = zonelist2[h2]; + zonelist[nhl] = new_h; assert(inyet[new_h] == 0); inyet[new_h] = 1; if (inyet2[new_h] != 2) - to_process.insert(new_h); + build_process_queue(to_process, z, inyet, new_h); nhl++; z[h].npjoin += z[new_h].np; } } - + for (int hl = 0; hl < nhl2; hl++) inyet2[zonelist2[hl]] = 0; if (nhl/10000 > nhlcount) { @@ -212,10 +213,12 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo (cout << format(" %d [%d]") % nhl % to_process.size()).flush(); nhlcount = nhl/10000; } - } while((lowvol < BIGFLT) && (!beaten)); - - z[h].denscontrast = z[h].leak/p[z[h].core].dens; - if (z[h].denscontrast < 1.) z[h].denscontrast = 1.; + } + while((lowvol < BIGFLT) && (!beaten)); + + z[h].denscontrast = z[h].leak/p[z[h].core].dens; + if (z[h].denscontrast < 1.) + z[h].denscontrast = 1.; /* Don't sort; want the core zone to be first */ @@ -233,19 +236,16 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo z[h].zonelist[q] = zonelist[q]; } - done_zones[h] = true; z[h].nhl = nhl; } delete[] zonelist; delete[] zonelist2; - delete[] links; delete[] inyet; delete[] inyet2; - delete[] done_zones; - + } delete[] iord; - + double maxdenscontrast = 0; #pragma omp parallel shared(maxdenscontrast) {