mirror of
https://bitbucket.org/cosmicvoids/vide_public.git
synced 2025-07-04 23:31:12 +00:00
Added openmp-ization of jozov2_watershed
This commit is contained in:
Guilhem Lavaux
parent
9b0e7e115e
commit
b555533fa3
1 changed files with 190 additions and 165 deletions
|
|
@ -1,3 +1,7 @@
|
||||||
|
#ifdef OPENMP
|
||||||
|
#include <omp.h>
|
||||||
|
#endif
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
#include <fstream>
|
#include <fstream>
|
||||||
#include <boost/format.hpp>
|
#include <boost/format.hpp>
|
||||||
|
|
@ -18,16 +22,6 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
|
||||||
float maxdenscontrast = 0;
|
float maxdenscontrast = 0;
|
||||||
bool *done_zones;
|
bool *done_zones;
|
||||||
|
|
||||||
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);
|
|
||||||
|
|
||||||
double *sorter = new double[numZones+1];
|
double *sorter = new double[numZones+1];
|
||||||
/* Assign sorter by probability (could use volume instead) */
|
/* Assign sorter by probability (could use volume instead) */
|
||||||
for (int h = 0; h < numZones; h++)
|
for (int h = 0; h < numZones; h++)
|
||||||
|
|
@ -42,181 +36,212 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
|
||||||
findrtop(sorter, numZones, iord, numZones);
|
findrtop(sorter, numZones, iord, numZones);
|
||||||
delete[] sorter;
|
delete[] sorter;
|
||||||
|
|
||||||
nhl = 0;
|
#pragma omp parallel
|
||||||
for (int ii = 0; ii < numZones; ii++)
|
{
|
||||||
{
|
inyet = new char[numZones];
|
||||||
int nhlcount = 0;
|
inyet2 = new char[numZones];
|
||||||
int h = iord[ii];
|
zonelist = new int[numZones];
|
||||||
float lowvol;
|
zonelist2 = new int[numZones];
|
||||||
bool beaten;
|
done_zones = new bool[numZones];
|
||||||
|
|
||||||
for (int hl = 0; hl < nhl; hl++)
|
|
||||||
inyet[zonelist[hl]] = 0;
|
|
||||||
|
|
||||||
zonelist[0] = h;
|
fill(inyet, inyet + numZones, 0);
|
||||||
inyet[h] = 1;
|
fill(inyet2, inyet2 + numZones, 0);
|
||||||
nhl = 1;
|
fill(done_zones, done_zones + numZones, false);
|
||||||
z[h].npjoin = z[h].np;
|
|
||||||
do {
|
|
||||||
/* Find the lowest-volume (highest-density) adjacency */
|
|
||||||
int nl = 0;
|
|
||||||
|
|
||||||
beaten = false;
|
|
||||||
lowvol = BIGFLT;
|
|
||||||
|
|
||||||
|
nhl = 0;
|
||||||
|
#pragma omp for schedule(dynamic,100)
|
||||||
|
for (int ii = 0; ii < numZones; ii++)
|
||||||
|
{
|
||||||
|
int nhlcount = 0;
|
||||||
|
int h = iord[ii];
|
||||||
|
float lowvol;
|
||||||
|
bool beaten;
|
||||||
|
|
||||||
for (int hl = 0; hl < nhl; hl++)
|
for (int hl = 0; hl < nhl; hl++)
|
||||||
{
|
inyet[zonelist[hl]] = 0;
|
||||||
int h2 = zonelist[hl];
|
|
||||||
if (inyet[h2] == 1) { /* If it's not already identified as
|
zonelist[0] = h;
|
||||||
an interior zone, with inyet=2 */
|
inyet[h] = 1;
|
||||||
bool interior = true;
|
nhl = 1;
|
||||||
for (int za = 0; za < z[h2].nadj; za++)
|
z[h].npjoin = z[h].np;
|
||||||
|
do {
|
||||||
|
/* Find the lowest-volume (highest-density) adjacency */
|
||||||
|
int nl = 0;
|
||||||
|
|
||||||
|
beaten = false;
|
||||||
|
lowvol = BIGFLT;
|
||||||
|
|
||||||
|
for (int hl = 0; hl < nhl; hl++)
|
||||||
|
{
|
||||||
|
int h2 = zonelist[hl];
|
||||||
|
if (inyet[h2] == 1) { /* If it's not already identified as
|
||||||
|
an interior zone, with inyet=2 */
|
||||||
|
bool interior = true;
|
||||||
|
for (int za = 0; za < z[h2].nadj; za++)
|
||||||
|
{
|
||||||
|
if (inyet[z[h2].adj[za]] == 0)
|
||||||
|
{
|
||||||
|
interior = false;
|
||||||
|
if (z[h2].slv[za] == lowvol)
|
||||||
|
{
|
||||||
|
links[nl] = z[h2].adj[za];
|
||||||
|
nl ++;
|
||||||
|
if (nl == NLINKS)
|
||||||
|
{
|
||||||
|
printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl);
|
||||||
|
exit(0);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (z[h2].slv[za] < lowvol)
|
||||||
|
{
|
||||||
|
lowvol = z[h2].slv[za];
|
||||||
|
links[0] = z[h2].adj[za];
|
||||||
|
nl = 1;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (interior)
|
||||||
|
inyet[h2] = 2; /* No bordering exter. zones */
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (nl == 0)
|
||||||
|
{
|
||||||
|
beaten = true;
|
||||||
|
z[h].leak = maxvol;
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
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. */
|
||||||
|
beaten = true;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (beaten)
|
||||||
|
{
|
||||||
|
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)
|
||||||
{
|
{
|
||||||
if (inyet[z[h2].adj[za]] == 0)
|
zonelist2[nhl2] = links[l];
|
||||||
|
inyet2[links[l]] = 1;
|
||||||
|
nhl2 ++;
|
||||||
|
bool added = true;
|
||||||
|
while (added && !beaten)
|
||||||
{
|
{
|
||||||
interior = false;
|
added = false;
|
||||||
if (z[h2].slv[za] == lowvol)
|
for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
|
||||||
{
|
{
|
||||||
links[nl] = z[h2].adj[za];
|
int h2 = zonelist2[hl];
|
||||||
nl ++;
|
|
||||||
if (nl == NLINKS)
|
if (inyet2[h2] == 1) {
|
||||||
{
|
bool interior = true; /* Guilty until proven innocent */
|
||||||
printf("Too many links with the same rho_sl! Increase NLINKS from %d\n",nl);
|
for (int za = 0; za < z[h2].nadj; za ++) {
|
||||||
exit(0);
|
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)
|
||||||
|
beaten = true;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
zonelist2[nhl2] = link2;
|
||||||
|
inyet2[link2] = 1;
|
||||||
|
nhl2++;
|
||||||
|
added = true;
|
||||||
|
}
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
if (interior)
|
||||||
if (z[h2].slv[za] < lowvol)
|
inyet2[h2] = 2;
|
||||||
{
|
}
|
||||||
lowvol = z[h2].slv[za];
|
|
||||||
links[0] = z[h2].adj[za];
|
|
||||||
nl = 1;
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
if (interior)
|
}
|
||||||
inyet[h2] = 2; /* No bordering exter. zones */
|
for (int hl = 0; hl < nhl2; hl++)
|
||||||
|
inyet2[zonelist2[hl]] = 0;
|
||||||
|
|
||||||
|
/* See if there's a beater */
|
||||||
|
if (beaten) {
|
||||||
|
z[h].leak = lowvol;
|
||||||
|
} else {
|
||||||
|
for (int h2 = 0; h2 < nhl2; h2++) {
|
||||||
|
zonelist[nhl] = zonelist2[h2];
|
||||||
|
inyet[zonelist2[h2]] = 1;
|
||||||
|
nhl++;
|
||||||
|
z[h].npjoin += z[zonelist2[h2]].np;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
if (nhl/10000 > nhlcount) {
|
||||||
|
nhlcount = nhl/10000;
|
||||||
|
printf(" %d",nhl); FF;
|
||||||
|
}
|
||||||
|
} while((lowvol < BIGFLT) && (!beaten));
|
||||||
|
|
||||||
if (nl == 0)
|
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
|
||||||
{
|
if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
|
||||||
beaten = true;
|
|
||||||
z[h].leak = maxvol;
|
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
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. */
|
|
||||||
beaten = true;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (beaten)
|
|
||||||
{
|
|
||||||
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)
|
|
||||||
beaten = true;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
zonelist2[nhl2] = link2;
|
|
||||||
inyet2[link2] = 1;
|
|
||||||
nhl2++;
|
|
||||||
added = true;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (interior)
|
|
||||||
inyet2[h2] = 2;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
for (int hl = 0; hl < nhl2; hl++)
|
|
||||||
inyet2[zonelist2[hl]] = 0;
|
|
||||||
|
|
||||||
/* See if there's a beater */
|
/* Don't sort; want the core zone to be first */
|
||||||
if (beaten) {
|
|
||||||
z[h].leak = lowvol;
|
if (nhlcount > 0) { /* Outputs the number of zones in large voids */
|
||||||
} else {
|
printf(" h%d:%d\n",h,nhl);
|
||||||
for (int h2 = 0; h2 < nhl2; h2++) {
|
FF;
|
||||||
zonelist[nhl] = zonelist2[h2];
|
|
||||||
inyet[zonelist2[h2]] = 1;
|
|
||||||
nhl++;
|
|
||||||
z[h].npjoin += z[zonelist2[h2]].np;
|
|
||||||
}
|
|
||||||
}
|
}
|
||||||
if (nhl/10000 > nhlcount) {
|
/* Calculate volume */
|
||||||
nhlcount = nhl/10000;
|
z[h].voljoin = 0.;
|
||||||
printf(" %d",nhl); FF;
|
z[h].zonelist = new int[nhl];
|
||||||
|
z[h].numzones = nhl;
|
||||||
|
for (int q = 0; q < nhl; q++) {
|
||||||
|
z[h].voljoin += z[zonelist[q]].vol;
|
||||||
|
z[h].zonelist[q] = zonelist[q];
|
||||||
}
|
}
|
||||||
} while((lowvol < BIGFLT) && (!beaten));
|
|
||||||
|
done_zones[h] = true;
|
||||||
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
|
z[h].nhl = nhl;
|
||||||
if (z[h].denscontrast < 1.) z[h].denscontrast = 1.;
|
|
||||||
|
|
||||||
/* find biggest denscontrast */
|
|
||||||
if (z[h].denscontrast > maxdenscontrast) {
|
|
||||||
maxdenscontrast = (double)z[h].denscontrast;
|
|
||||||
}
|
}
|
||||||
|
delete[] zonelist;
|
||||||
/* Don't sort; want the core zone to be first */
|
delete[] zonelist2;
|
||||||
|
delete[] links;
|
||||||
if (nhlcount > 0) { /* Outputs the number of zones in large voids */
|
delete[] iord;
|
||||||
printf(" h%d:%d\n",h,nhl);
|
|
||||||
FF;
|
}
|
||||||
|
|
||||||
|
maxdenscontrast = 0;
|
||||||
|
#pragma omp parallel shared(maxdenscontrast)
|
||||||
|
{
|
||||||
|
double maxdenscontrast_local = 0;
|
||||||
|
|
||||||
|
#pragma omp for schedule(static)
|
||||||
|
for (int h = 0; h < numZones; h++)
|
||||||
|
{
|
||||||
|
/* find biggest denscontrast */
|
||||||
|
if (z[h].denscontrast > maxdenscontrast_local) {
|
||||||
|
maxdenscontrast_local = (double)z[h].denscontrast;
|
||||||
|
}
|
||||||
}
|
}
|
||||||
/* Calculate volume */
|
|
||||||
z[h].voljoin = 0.;
|
#pragma omp critical
|
||||||
z[h].zonelist = new int[nhl];
|
{
|
||||||
z[h].numzones = nhl;
|
if (maxdenscontrast_local > maxdenscontrast)
|
||||||
for (int q = 0; q < nhl; q++) {
|
maxdenscontrast = maxdenscontrast_local;
|
||||||
z[h].voljoin += z[zonelist[q]].vol;
|
|
||||||
z[h].zonelist[q] = zonelist[q];
|
|
||||||
}
|
|
||||||
|
|
||||||
done_zones[h] = true;
|
|
||||||
z[h].nhl = nhl;
|
|
||||||
}
|
}
|
||||||
delete[] zonelist;
|
}
|
||||||
delete[] zonelist2;
|
|
||||||
delete[] links;
|
|
||||||
delete[] iord;
|
|
||||||
|
|
||||||
|
|
||||||
cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl;
|
cout << format("Maxdenscontrast = %f.") % maxdenscontrast << endl;
|
||||||
}
|
}
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue