guilhem_lavaux/vide_public
0
0
Fork 0
mirror of https://bitbucket.org/cosmicvoids/vide_public.git synced 2025-07-04 23:31:12 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

Merge branch 'master' of bitbucket.org:cosmicvoids/void_identification

Browse source
This commit is contained in:
P.M. Sutter 2013-04-04 17:11:43 -05:00
commit c999e9e235
3 changed files with 182 additions and 151 deletions
Download patch file Download diff file Expand all files Collapse all files

291
c_tools/zobov2/jozov2_watershed.cpp
View file

@ -2,6 +2,7 @@
#include <omp.h> #include <omp.h>
#endif #endif
#include <queue>
#include <set> #include <set>
#include <iostream> #include <iostream>
#include <fstream> #include <fstream>
@ -13,196 +14,195 @@
using namespace std; using namespace std;
using boost::format; using boost::format;
struct ZoneDensityPair
{
int h;
double density;
double core;
bool operator<(const ZoneDensityPair& p2) const
{
return (density > p2.density) || (density==p2.density && core > p2.core);
}
};
typedef priority_queue<ZoneDensityPair> ZoneQueue;
static void build_process_queue(ZoneQueue& q, ZONE *z, PARTICLE *p, 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];
zdp.core = p[z[zdp.h].core].dens;
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) void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, float voltol)
{ {
int *iord;
double *sorter = new double[numZones+1];
/* Assign sorter by probability (could use volume instead) */
for (int h = 0; h < numZones; h++)
sorter[h] = (double)z[h].core;
/* Text output file */ /* Text output file */
printf("about to sort (pre-watershed sort) ...\n");FF;
iord = new int[numZones];
findrtop(sorter, numZones, iord, numZones);
delete[] sorter;
#pragma omp parallel #pragma omp parallel
{ {
char *inyet, *inyet2; char *inyet, *inyet2;
int *zonelist, *zonelist2; int *zonelist, *zonelist2;
int nhl; int nhl;
int *links = new int[NLINKS];
bool *done_zones;
int prev_ii = -1; int prev_ii = -1;
inyet = new char[numZones]; inyet = new char[numZones];
inyet2 = new char[numZones]; inyet2 = new char[numZones];
zonelist = new int[numZones]; zonelist = new int[numZones];
zonelist2 = new int[numZones]; zonelist2 = new int[numZones];
done_zones = new bool[numZones];
fill(inyet, inyet + numZones, 0); fill(inyet, inyet + numZones, 0);
fill(inyet2, inyet2 + numZones, 0); fill(inyet2, inyet2 + numZones, 0);
fill(done_zones, done_zones + numZones, false);
nhl = 0; nhl = 0;
#pragma omp for schedule(dynamic,1) #pragma omp for schedule(dynamic,1)
for (int h = 0; h < numZones; h++) for (int h = 0; h < numZones; h++)
{ {
int nhlcount = 0; int nhlcount = 0;
float lowvol; float previous_lowvol = BIGFLT, lowvol, z_cur_core_dens;
bool beaten; bool beaten;
set<int> to_process; priority_queue<ZoneDensityPair> to_process;
int link0;
int save_nhl;
pid_t save_npjoin;
for (int hl = 0; hl < nhl; hl++) for (int hl = 0; hl < nhl; hl++)
inyet[zonelist[hl]] = 0; inyet[zonelist[hl]] = 0;
zonelist[0] = h; zonelist[0] = h;
inyet[h] = 1; inyet[h] = 1;
nhl = 1; save_nhl = nhl = 1;
to_process.insert(h); save_npjoin = z[h].npjoin = z[h].np;
z[h].npjoin = z[h].np; z_cur_core_dens = p[z[h].core].dens;
build_process_queue(to_process, z, p, inyet, h);
do { do {
/* Find the lowest-volume (highest-density) adjacency */ /* Find the lowest-volume (highest-density) adjacency */
int nl = 0;
beaten = false; beaten = false;
lowvol = BIGFLT;
set<int>::iterator iter = to_process.begin(); if (to_process.empty())
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)
{ {
beaten = true; beaten = true;
z[h].leak = maxvol; z[h].leak = maxvol;
save_npjoin = z[h].npjoin;
save_nhl = nhl;
continue; continue;
} }
if (lowvol > voltol) 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())
{
save_npjoin = z[h].npjoin;
save_nhl = nhl;
beaten = true;
z[h].leak = maxvol;
continue;
}
/* See if there's a beater */
if (previous_lowvol != lowvol)
{
save_npjoin = z[h].npjoin;
save_nhl = nhl;
previous_lowvol = lowvol;
}
if (lowvol > voltol)
{ {
beaten = true; beaten = true;
z[h].leak = lowvol; z[h].leak = lowvol;
continue; continue;
} }
for (int l = 0; l < nl; l++) if (p[z[link0].core].dens < z_cur_core_dens)
{
//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)
{ {
beaten = true;
z[h].leak = lowvol; z[h].leak = lowvol;
continue; continue;
} }
/* Add everything linked to the link(s) */ /* Add everything linked to the link(s) */
int nhl2 = 0; int nhl2 = 0;
for (int l = 0; l < nl; l++) zonelist2[0] = link0;
{ inyet2[link0] = 1;
if (inyet2[links[l]] == 0) nhl2=1;
{
zonelist2[nhl2] = links[l]; bool added = true;
inyet2[links[l]] = 1; while (added && !beaten)
nhl2 ++; {
bool added = true; added = false;
while (added && !beaten) for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
{ {
added = false; int h2 = zonelist2[hl];
for (int hl = 0; (hl < nhl2) && (!beaten); hl++)
{ if (inyet2[h2] == 1) {
int h2 = zonelist2[hl]; bool interior = true; /* Guilty until proven innocent */
if (inyet2[h2] == 1) { for (int za = 0; za < z[h2].nadj; za ++) {
bool interior = true; /* Guilty until proven innocent */ int link2 = z[h2].adj[za];
for (int za = 0; za < z[h2].nadj; za ++) {
int link2 = z[h2].adj[za]; if ((inyet[link2]+inyet2[link2]) == 0) {
if ((inyet[link2]+inyet2[link2]) == 0) { interior = false;
interior = false; if (z[h2].slv[za] <= lowvol) {
if (z[h2].slv[za] <= lowvol) { if (p[z[link2].core].dens < z_cur_core_dens) {
//if (!done_zones[link2]) { // Equivalent to p[z[link2].core].dens < p[z[h].core].dens) beaten = true;
if (p[z[link2].core].dens < p[z[h].core].dens) { break;
beaten = true; }
break; zonelist2[nhl2] = link2;
} inyet2[link2] = 1;
zonelist2[nhl2] = link2; nhl2++;
inyet2[link2] = 1; added = true;
nhl2++; }
added = true; }
} }
}
} if (interior)
if (interior) inyet2[h2] = 2;
inyet2[h2] = 2; }
} }
} }
}
} /* See if there's a beater */
} if (beaten) {
z[h].leak = lowvol;
/* See if there's a beater */ } else {
if (beaten) {
z[h].leak = lowvol;
} else {
for (int h2 = 0; h2 < nhl2; h2++) {
int new_h = zonelist2[h2];
for (int h2 = 0; h2 < nhl2; h2++) {
int new_h = zonelist2[h2];
zonelist[nhl] = new_h; zonelist[nhl] = new_h;
assert(inyet[new_h] == 0); assert(inyet[new_h] == 0);
z[h].npjoin += z[new_h].np;
inyet[new_h] = 1; inyet[new_h] = 1;
if (inyet2[new_h] != 2) if (inyet2[new_h] != 2)
to_process.insert(new_h); build_process_queue(to_process, z, p, inyet, new_h);
nhl++; nhl++;
z[h].npjoin += z[new_h].np;
} }
} }
for (int hl = 0; hl < nhl2; hl++) for (int hl = 0; hl < nhl2; hl++)
inyet2[zonelist2[hl]] = 0; inyet2[zonelist2[hl]] = 0;
if (nhl/10000 > nhlcount) { if (nhl/10000 > nhlcount) {
@ -212,10 +212,18 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
(cout << format(" %d [%d]") % nhl % to_process.size()).flush(); (cout << format(" %d [%d]") % nhl % to_process.size()).flush();
nhlcount = nhl/10000; nhlcount = nhl/10000;
} }
} while((lowvol < BIGFLT) && (!beaten)); }
while((lowvol < BIGFLT) && (!beaten));
z[h].denscontrast = z[h].leak/p[z[h].core].dens;
if (z[h].denscontrast < 1.) z[h].denscontrast = 1.; if (!beaten)
{
save_npjoin = z[h].npjoin;
save_nhl = nhl;
}
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 */ /* Don't sort; want the core zone to be first */
@ -225,27 +233,24 @@ void doWatershed(PARTICLE *p, pid_t np, ZONE *z, int numZones, float maxvol, flo
FF; FF;
} }
/* Calculate volume */ /* Calculate volume */
z[h].npjoin = save_npjoin;
z[h].voljoin = 0.; z[h].voljoin = 0.;
z[h].zonelist = new int[nhl]; z[h].zonelist = new int[save_nhl];
z[h].numzones = nhl; z[h].numzones = save_nhl;
for (int q = 0; q < nhl; q++) { for (int q = 0; q < save_nhl; q++) {
z[h].voljoin += z[zonelist[q]].vol; z[h].voljoin += z[zonelist[q]].vol;
z[h].zonelist[q] = zonelist[q]; z[h].zonelist[q] = zonelist[q];
} }
done_zones[h] = true; z[h].nhl = save_nhl;
z[h].nhl = nhl;
} }
delete[] zonelist; delete[] zonelist;
delete[] zonelist2; delete[] zonelist2;
delete[] links;
delete[] inyet; delete[] inyet;
delete[] inyet2; delete[] inyet2;
delete[] done_zones;
} }
delete[] iord;
double maxdenscontrast = 0; double maxdenscontrast = 0;
#pragma omp parallel shared(maxdenscontrast) #pragma omp parallel shared(maxdenscontrast)
{ {

7
c_tools/zobov2/jozov2_zones.cpp
View file

@ -39,10 +39,13 @@ void buildInitialZones(PARTICLE *p, pid_t np, pid_t* jumped,
} }
(cout << "Post-jump ..." << endl).flush(); (cout << "Post-jump ..." << endl).flush();
numZones = 0; pid_t loc_NumZones = 0;
#pragma omp parallel for schedule(static) reduction(+:loc_NumZones)
for (pid_t i = 0; i < np; i++) for (pid_t i = 0; i < np; i++)
if (numinh[i] > 0) if (numinh[i] > 0)
numZones++; loc_NumZones++;
numZones = loc_NumZones;
cout << format("%d initial zones found") % numZones << endl; cout << format("%d initial zones found") % numZones << endl;
delete[] jumper; delete[] jumper;

35
zobov/voztie.c
View file

@ -2,8 +2,21 @@
#include <math.h> #include <math.h>
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h>
#include "voz.h" #include "voz.h"
int compare_int(const void *a, const void *b)
{
const int *ia = (const int *)a;
const int *ib = (const int *)b;
if ((*ia) < (*ib))
return -1;
else if ((*ia) > (*ib))
return 1;
else
return 0;
}
int main(int argc, char *argv[]) { int main(int argc, char *argv[]) {
FILE *part, *adj, *vol; FILE *part, *adj, *vol;
@ -135,15 +148,25 @@ int main(int argc, char *argv[]) {
fread(&na,1,sizeof(int),part); fread(&na,1,sizeof(int),part);
pid_t pid = orig[p]; pid_t pid = orig[p];
if (na > 0) { if (na > 0) {
assert(adjs[pid].nadj == 0);// || adjs[pid].nadj == na); pid_t *previous_adj = adjs[pid].adj;
adjs[orig[p]].nadj = na; assert(adjs[pid].nadj == 0 || adjs[pid].nadj == na);
if (adjs[pid].adj == 0) adjs[pid].nadj = na;
adjs[orig[p]].adj = (pid_t *)malloc(na*sizeof(pid_t)); adjs[pid].adj = (pid_t *)malloc(na*sizeof(pid_t));
if (adjs[orig[p]].adj == 0) { if (adjs[pid].adj == 0) {
printf("Couldn't allocate adjs[orig[%d]].adj.\n",p); printf("Couldn't allocate adjs[orig[%d]].adj.\n",p);
exit(0); exit(0);
} }
fread(adjs[orig[p]].adj,na,sizeof(pid_t),part); fread(adjs[pid].adj,na,sizeof(pid_t),part);
if (previous_adj != 0)
{
qsort(previous_adj, na, sizeof(pid_t), compare_int);
qsort(adjs[pid].adj, na, sizeof(pid_t), compare_int);
if (memcmp(previous_adj, adjs[pid].adj, sizeof(pid_t)*na) != 0)
{
printf("Inconsistent adjacencies between two divisions.\n");
abort();
}
}
} else { } else {
printf("0"); fflush(stdout); printf("0"); fflush(stdout);
} }