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

.gitignore

@@ -1,10 +1,6 @@
 *~
 *.o
 *.prog
-voz1b1
-jozov
-vozinit
-voztie
 .mydepend
 ./pipeline/miscTools/*.out
 

c_tools/mock/loaders/sdf_loader.cpp

@@ -166,7 +166,7 @@ public:
       }
     if (load_flags & NEED_VELOCITY)
       {
-        const char *name_template[3] = { "x", "y", "z" };
+        const char *name_template[3] = { "vx", "vy", "vz" };
         for (int q = 0; q < 3; q++)
           {
             d->Vel[q] = new float[numPartToLoad];

c_tools/zobov2/voz1b1/voz.h

@@ -0,0 +1,29 @@
+#ifndef __VOZ_H
+#define __VOZ_H
+
+#define MAXVERVER 100000
+#define NGUARD 84 /*Actually, the number of SPACES between guard points
+##define NGUARD 42 /*Actually, the number of SPACES between guard points
+		    in each dim */
+
+typedef int pid_t;
+
+typedef struct Partadj {
+  pid_t nadj;
+  pid_t *adj;
+} PARTADJ;
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs);
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol);
+int posread(char *posfile, float ***p, float fact);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif

c_tools/zobov2/voz1b1/voz1b1.cpp

@@ -0,0 +1,404 @@
+#include <cassert>
+#include <boost/format.hpp>
+#include <iostream>
+#include <fstream>
+#include <CosmoTool/miniargs.hpp>
+#include <cmath>
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+#include "voz_io.hpp"
+
+using CosmoTool::MINIARG_STRING;
+using CosmoTool::MINIARG_DOUBLE;
+using CosmoTool::MINIARG_INT;
+using CosmoTool::MINIARG_NULL;
+using boost::format;
+using namespace std;
+
+#define DL for (d=0;d<3;d++)
+
+bool checkParameters(int *numdiv, int *b)
+{
+  for (int i = 0; i < 3; i++)
+    {
+      if (numdiv[i] < 0 || b[i] < 0 || b[i] >= numdiv[i])
+        return false;
+    }
+  return true;
+}
+
+
+struct BoxData
+{
+  float width[3], totwidth[3];
+  float width2[3], totwidth2[3];
+  float bf, s[3], g[3], c[3];
+  coordT *parts;
+  pid_t *orig;
+  pid_t nvpall, nvp, nvpbuf;
+  bool guard_added;
+  double xyz_min[3], xyz_max[3];
+
+  void prepareBox(const PositionData& pdata, int *b_id, int *numdivs);
+
+  void checkParticle(double *xyz, bool& in_main, bool& in_buf);
+  void addGuardPoints();
+  void findBoundary();
+};
+
+void BoxData::checkParticle(double *xyz, bool& in_main, bool& in_buf)
+{
+  in_main = in_buf = true;
+  for (int d = 0; d < 3; d++)
+    {
+      xyz[d] -= (double)c[d];
+      if (xyz[d] > width2[d])
+        xyz[d] -= width[d];
+      if (xyz[d] < -width2[d])
+        xyz[d] += width[d];
+
+      in_buf = in_buf && (fabs(xyz[d]) < totwidth2[d]);
+      in_main = in_main && (fabs(xyz[d]) <= width2[d]);
+    }
+  assert(!in_main || in_buf);
+}
+
+void BoxData::prepareBox(const PositionData& pdata, int *b_id, int *numdivs)
+{
+  double BF = std::numeric_limits<double>::max();
+  guard_added = false;
+
+  for (int i = 0; i < 3; i++)
+    {
+      width[i] = (pdata.xyz_max[i] - pdata.xyz_min[i])/numdivs[i];
+      width2[i] = 0.5*width[i];
+
+      totwidth[i] = width[i]+2.*bf;
+      totwidth2[i] = width2[i] + bf;
+
+      s[i] = width[i]/(float)NGUARD;
+      if ((bf*bf - 2.*s[i]*s[i]) < 0.)
+        {
+          printf("bf = %f, s = %f.\n",bf,s[i]);
+          printf("Not enough guard points for given border.\nIncrease guards to >= %f\n.",
+                 sqrt(2.)*width[i]/bf);
+          exit(0);
+        }
+      g[i] = (bf / 2.)*(1. + sqrt(1 - 2.*s[i]*s[i]/(bf*bf)));
+      cout << format("s[%d] = %f, bf = %f, g[%d] = %f.") % i%  s[i] % bf % i % g[i] << endl;
+
+      c[i] = b_id[i] * width[i];
+    }
+  cout.flush();
+
+  cout << format("c: %f,%f,%f") % c[0] % c[1] % c[2] << endl;
+
+  /* Assign temporary array*/
+  nvpbuf = 0; /* Number of particles to tesselate, including buffer */
+  nvp = 0; /* Without the buffer */
+
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+      bool is_it_in_buf, is_it_in_main;
+
+      checkParticle(xyz, is_it_in_main, is_it_in_buf);
+
+      if (is_it_in_buf)
+        nvpbuf++;
+      if (is_it_in_main)
+        nvp++;
+      assert(nvp <= nvpbuf);
+    }
+
+  nvpbuf += 6*(NGUARD+1)*(NGUARD+1); /* number of guard points */
+  nvpall = nvpbuf;
+
+  parts = new coordT[3*nvpall];
+  orig = new pid_t[nvpall];
+
+  if (parts == 0)
+    {
+      cout << "Unable to allocate parts" << endl;
+      exit(0);
+    }
+  if (orig == 0)
+    {
+      cout << "Unable to allocate orig" << endl;
+      exit(0);
+    }
+
+  nvp = 0; /* nvp = number of particles without buffer */
+
+  for (int j = 0; j < 3; j++)
+    {
+      xyz_min[j] = BF;
+      xyz_max[j] = -BF;
+    }
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      bool is_it_in_main, is_it_in_buf;
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+
+      checkParticle(xyz, is_it_in_main, is_it_in_buf);
+
+      if (is_it_in_main) {
+        assert(nvp < nvpall);
+        for (int j = 0; j < 3; j++)
+          {
+            parts[3*nvp+j] = xyz[j];
+            xyz_min[j] = min(xyz_min[j], xyz[j]);
+            xyz_max[j] = max(xyz_max[j], xyz[j]);
+          }
+        orig[nvp] = i;
+        nvp++;
+      }
+  }
+  cout << format("nvp = %d") %nvp << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl;
+  nvpbuf = nvp;
+  for (pid_t i = 0; i < pdata.np; i++)
+    {
+      bool is_it_in_main, is_it_in_buf;
+      double xyz[3] = { pdata.xyz[0][i], pdata.xyz[1][i], pdata.xyz[2][i] };
+
+      checkParticle(xyz, is_it_in_main, is_it_in_buf);
+
+      if (is_it_in_buf && !is_it_in_main)
+        {
+          assert(nvpbuf < nvpall);
+          for (int j = 0; j < 3; j++)
+            {
+              parts[3*nvpbuf+j] = xyz[j];
+              xyz_min[j] = min(xyz_min[j], xyz[j]);
+              xyz_max[j] = max(xyz_max[j], xyz[j]);
+            }
+          orig[nvpbuf] = i;
+
+          nvpbuf++;
+        }
+    }
+  cout << format("nvpbuf = %d") % nvpbuf << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f\n") % xyz_min[0] % xyz_max[0] % xyz_min[1] % xyz_max[1] % xyz_min[2] % xyz_max[2] << endl;
+
+  double predict = 1;
+  for (int j = 0; j < 3; j++)
+    predict *= (width[j]+2*bf);
+  predict *= pdata.np;
+  cout << format("There should be ~ %g points; there are %d\n") % predict % nvpbuf << endl;
+}
+
+void BoxData::addGuardPoints()
+{
+  int rot_g[3][3]  = { {0, 0, 1}, {0,1,0}, {1,0,0} };
+  int rot_si[3][3] = { {1, 0, 0}, {1,0,0}, {0,1,0} };
+  int rot_sj[3][3] = { {0, 1, 0}, {0,0,1}, {0,0,1} };
+
+  pid_t nvpcur = nvpbuf;
+
+  for (int k = 0; k < 3; k++)
+    {
+
+      /* Add guard points */
+      for (int i=0; i <= NGUARD; i++)
+        {
+          for (int j=0; j <= NGUARD; j++)
+            {
+              assert(nvpcur < nvpall);
+              /* Bottom */
+              for (int a = 0; a < 3; a++)
+                parts[3*nvpcur+a] = -width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] - rot_g[k][a] * g[a];
+
+              nvpcur++;
+
+              assert(nvpcur < nvpall);
+              /* Top */
+              for (int a = 0; a < 3; a++)
+                parts[3*nvpcur+a] = (2*rot_g[k][a]-1)*width2[a] + (realT)i * s[a] * rot_si[k][a] + (realT)j * s[a] * rot_sj[k][a] + rot_g[k][a] * g[a];
+
+              nvpcur++;
+            }
+        }
+    }
+}
+
+void BoxData::findBoundary()
+{
+  double BF = std::numeric_limits<double>::max();
+
+  for (int j = 0; j < 3; j++)
+    {
+      xyz_min[j] = BF;
+      xyz_max[j] = -BF;
+    }
+
+  for (pid_t i = nvpbuf; i < nvpall; i++) {
+    for (int j = 0; j < 3; j++)
+      {
+        xyz_min[j] = std::min(xyz_min[j], parts[3*i+j]);
+        xyz_max[j] = std::max(xyz_max[j], parts[3*i+j]);
+      }
+  }
+}
+
+void saveTesselation(const string& outfile, PositionData& pdata, BoxData& boxdata, PARTADJ *adjs, float *vols)
+{
+  ofstream out(outfile.c_str());
+  if (!out)
+    {
+      cout << format("Unable to open %s") % outfile << endl;
+      exit(0);
+    }
+  out.write((char *)&pdata.np, sizeof(pid_t));
+  out.write((char *)&boxdata.nvp, sizeof(pid_t));
+  cout << format("nvp = %d") % boxdata.nvp << endl;
+
+  /* Tell us where the original particles were */
+  out.write((char *)boxdata.orig, sizeof(pid_t)*boxdata.nvp);
+  /* Volumes*/
+  out.write((char *)vols,sizeof(float)*boxdata.nvp);
+  /* Adjacencies */
+  for (pid_t i = 0; i < boxdata.nvp; i++)
+    {
+      out.write((char*)&adjs[i].nadj, sizeof(pid_t));
+      if (adjs[i].nadj > 0)
+        out.write((char *)adjs[i].adj, adjs[i].nadj*sizeof(pid_t));
+      else
+        (cout << "0").flush();
+    }
+  out.close();
+}
+
+int main(int argc, char *argv[]) {
+  PositionData pdata;
+  BoxData boxdata;
+
+  coordT deladjs[3*MAXVERVER], points[3*MAXVERVER];
+  pointT intpoints[3*MAXVERVER];
+  string outfile;
+  ofstream out;
+
+  char *suffix, *outDir, *posfile;
+  PARTADJ *adjs;
+  float *vols;
+  double border, boxsize[3];
+  int b[3], numdiv[3];
+  double totalvol;
+
+  CosmoTool::MiniArgDesc args[] = {
+    { "POSITION FILE", &posfile, MINIARG_STRING },
+    { "BORDER SIZE", &border, MINIARG_DOUBLE },
+    { "BOX_X", &boxsize[0], MINIARG_DOUBLE },
+    { "BOX_Y", &boxsize[1], MINIARG_DOUBLE },
+    { "BOX_Z", &boxsize[2], MINIARG_DOUBLE },
+    { "SUFFIX", &suffix, MINIARG_STRING },
+    { "NUM_DIVISION_X", &numdiv[0], MINIARG_INT },
+    { "NUM_DIVISION_Y", &numdiv[1], MINIARG_INT },
+    { "NUM_DIVISION_Z", &numdiv[2], MINIARG_INT },
+    { "B0", &b[0], MINIARG_INT },
+    { "B1", &b[1], MINIARG_INT },
+    { "B2", &b[2], MINIARG_INT },
+    { "OUTPUT DIRECTORY", &outDir, MINIARG_STRING },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!CosmoTool::parseMiniArgs(argc, argv, args))
+    return 1;
+
+  if (!checkParameters(numdiv, b))
+    return 2;
+
+  /* Boxsize should be the range in r, yielding a range 0-1 */
+  if (!pdata.readFrom(posfile))
+    return 3;
+
+  (cout << pdata.np << " particles" << endl).flush();
+
+  pdata.findExtrema();
+
+  (cout << boost::format("np: %d, x: %f,%f; y: %f,%f; z: %f,%f")
+    % pdata.np
+    % pdata.xyz_min[0] % pdata.xyz_max[0]
+    % pdata.xyz_min[1] % pdata.xyz_max[1]
+    % pdata.xyz_min[2] % pdata.xyz_max[2]).flush();
+
+
+  if (border > 0.)
+    boxdata.bf = border;
+  else
+    boxdata.bf = 0.1;
+
+  boxdata.prepareBox(pdata, b, numdiv);
+  pdata.destroy();
+  boxdata.addGuardPoints();
+
+  adjs = new PARTADJ[boxdata.nvpall];
+  if (adjs == 0)
+    {
+      cout << "Unable to allocate adjs" << endl;
+      return 0;
+    }
+
+  boxdata.findBoundary();
+
+  cout << format("Added guard points to total %d points (should be %d)")
+            % boxdata.nvpall % (boxdata.nvpbuf + 6*(NGUARD+1)*(NGUARD+1)) << endl;
+  cout << format("x: %f,%f; y: %f,%f; z:%f,%f") % boxdata.xyz_min[0] % boxdata.xyz_max[0] % boxdata.xyz_min[1] % boxdata.xyz_max[1] % boxdata.xyz_min[2] % boxdata.xyz_max[2] << endl;
+
+  /* Do tesselation*/
+  printf("File read.  Tessellating ...\n"); fflush(stdout);
+  int exitcode = delaunadj(boxdata.parts, boxdata.nvp, boxdata.nvpbuf, boxdata.nvpall, &adjs);
+  if (exitcode != 0)
+   {
+     printf("Error while tesselating. Stopping here."); fflush(stdout);
+     return 4;
+   }
+
+  /* Now calculate volumes*/
+  printf("Now finding volumes ...\n"); fflush(stdout);
+  vols = new float[boxdata.nvp];
+
+  for (pid_t i = 0; i < boxdata.nvp; i++)
+    { /* Just the original particles
+       * Assign adjacency coordinate array*/
+      /* Volumes */
+      for (int j = 0; j < adjs[i].nadj; j++)
+        {
+          for (int d = 0; d < 3; d++)
+            {
+              deladjs[3*j + d] = boxdata.parts[3*adjs[i].adj[j]+d] - boxdata.parts[3*i+d];
+
+              if (deladjs[3*j+d] < -boxdata.width2[d])
+                deladjs[3*j+d] += boxdata.width[d];
+              if (deladjs[3*j+d] > boxdata.width2[d])
+                deladjs[3*j+d] -= boxdata.width[d];
+            }
+        }
+
+      exitcode = vorvol(deladjs, points, intpoints, adjs[i].nadj, &(vols[i]));
+      vols[i] *= pdata.np/pdata.V0;
+      if ((i % 1000) == 0)
+        cout << format("%d: %d, %f") % i % adjs[i].nadj % vols[i] << endl;
+    }
+
+  /* Get the adjacencies back to their original values */
+  for (pid_t i=0; i<boxdata.nvp; i++)
+    for (int j=0; j < adjs[i].nadj; j++)
+      adjs[i].adj[j] = boxdata.orig[adjs[i].adj[j]];
+
+  totalvol = 0.;
+  for (pid_t i=0;i<boxdata.nvp; i++)
+    totalvol += vols[i];
+
+  cout << format("Average volume = %g") % (totalvol/boxdata.nvp) << endl;
+
+  /* Now the output!
+     First number of particles */
+  outfile = str(format("%s/part.%s.%02d.%02d.%02d") % outDir % suffix % b[0] % b[1] % b[2]);
+
+  cout << format("Output to %s") %outfile << endl << endl;
+  saveTesselation(outfile, pdata, boxdata, adjs, vols);
+  delete[] adjs;
+
+  return(0);
+}

c_tools/zobov2/voz1b1/voz_io.cpp

@@ -0,0 +1,67 @@
+#include <limits>
+#include <iostream>
+#include <fstream>
+#include "voz_io.hpp"
+#include "CosmoTool/fortran.hpp"
+#include <cmath>
+
+using namespace CosmoTool;
+using namespace std;
+
+bool PositionData::readFrom(const string& fname)
+{
+  try
+    {
+      UnformattedRead f(fname.c_str());
+
+      f.beginCheckpoint();
+      np = f.readInt32();
+      f.endCheckpoint();
+
+      for (int j = 0; j < 3; j++)
+        {
+          xyz[j] = new float[np];
+
+          f.beginCheckpoint();
+          for (int p = 0; p < np; p++)
+            xyz[j][p] = f.readReal32();
+          f.endCheckpoint();
+        }
+    }
+  catch (const NoSuchFileException& e)
+    {
+      return false;
+    }
+  catch (const InvalidUnformattedAccess& e)
+    {
+      return false;
+    }
+  catch (const EndOfFileException& e)
+    {
+      return false;
+    }
+
+  return true;
+}
+
+void PositionData::findExtrema()
+{
+  const float BF = std::numeric_limits<float>::max();
+
+  for (int i = 0; i < 3; i++)
+    {
+      xyz_min[i] = BF;
+      xyz_max[i] = -BF;
+    }
+
+  V0=1;
+  for (int i = 0; i < 3; i++)
+    {
+      for (pid_t p = 0; p < np; p++)
+        {
+          xyz_min[i] = min(xyz_min[i], xyz[i][p]);
+          xyz_max[i] = max(xyz_max[i], xyz[i][p]);
+        }
+      V0 *= (xyz_max[i]-xyz_min[i]);
+    }
+}

c_tools/zobov2/voz1b1/voz_io.hpp

@@ -0,0 +1,24 @@
+#ifndef __VOZ_IO_HPP
+#define __VOZ_IO_HPP
+
+#include <string>
+
+struct PositionData
+{
+  float *xyz[3];
+  pid_t np;
+  float xyz_min[3], xyz_max[3];
+  float V0;
+
+  void destroy()
+  {
+    for (int j = 0; j < 3; j++)
+      delete[] xyz[j];
+  }
+
+  void findExtrema();
+
+  bool readFrom(const std::string& fname);
+};
+
+#endif

c_tools/zobov2/voz1b1/vozutil.c

@@ -0,0 +1,245 @@
+#include "libqhull/qhull_a.h"
+#include "voz.h"
+
+#define FOREACHvertex2_(vertices) FOREACHsetelement_(vertexT, vertices2,vertex2)
+
+/*char qh_version[] = "user_eg 3.1 2001/10/04";  */
+
+int compar(const void * n1, const void * n2) {
+  int i1,i2;
+
+  i1 = *(int *)n1;
+  i2 = *(int *)n2;
+  return 2*(i1 > i2) - 1 + (i1 == i2);
+}
+
+/* Finds Delaunay adjacencies of a set of points */
+int delaunadj (coordT *points, int nvp, int nvpbuf, int nvpall, PARTADJ **adjs) {
+  
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= stdout;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+  int i, ver, count;
+  int numfacets, numsimplicial, numridges, totneighbors, numneighbors,
+    numcoplanars, numtricoplanars;
+  setT *vertices, *vertices2, *vertex_points, *coplanar_points;
+  vertexT *vertex, **vertexp;
+  vertexT *vertex2, **vertex2p;
+  int vertex_i, vertex_n;
+  facetT *facet, **facetp, *neighbor, **neighborp;
+  pointT *point, **pointp;
+  int numdiv;
+  
+  PARTADJ adjst;
+  
+  int errorReported = 0;
+
+  adjst.adj = (int *)malloc(MAXVERVER*sizeof(int));
+  if (adjst.adj == NULL) {
+    printf("Unable to allocate adjst.adj\n");
+    exit(0);
+  }
+
+  /* Delaunay triangulation*/
+  sprintf (flags, "qhull s d Qt");
+  exitcode= qh_new_qhull (dim, nvpall, points, ismalloc,
+                      flags, outfile, errfile); 
+  
+
+  if (!exitcode) {                  /* if no error */
+    /* 'qh facet_list' contains the convex hull */
+    
+    /* From qh_printvneighbors */
+    qh_countfacets(qh facet_list, NULL, 0, &numfacets, &numsimplicial, 
+		   &totneighbors, &numridges, &numcoplanars, &numtricoplanars);
+    qh_vertexneighbors();
+    vertices= qh_facetvertices (qh facet_list, NULL, 0);
+    vertex_points= qh_settemp (nvpall);
+    coplanar_points= qh_settemp (nvpall);
+    qh_setzero (vertex_points, 0, nvpall);
+    qh_setzero (coplanar_points, 0, nvpall);
+    FOREACHvertex_(vertices)
+      qh_point_add (vertex_points, vertex->point, vertex);
+    FORALLfacet_(qh facet_list) {
+      FOREACHpoint_(facet->coplanarset)
+	qh_point_add (coplanar_points, point, facet);
+    }
+    ver = 0;
+    FOREACHvertex_i_(vertex_points) {
+      (*adjs)[ver].nadj = 0;
+      if (vertex) { 
+	/* Count the neighboring vertices, check that all are real
+	   neighbors */
+	adjst.nadj = 0;
+	FOREACHneighbor_(vertex) {
+	  if ((*adjs)[ver].nadj > -1) {
+	    if (neighbor->visitid) {
+	      vertices2 = neighbor->vertices;
+	      FOREACHvertex2_(vertices2) {
+		if (ver != qh_pointid(vertex2->point)) {
+		  adjst.adj[adjst.nadj] = (int)qh_pointid(vertex2->point);
+		  adjst.nadj ++;
+		  if (adjst.nadj > MAXVERVER) {
+		    printf("Increase MAXVERVER to at least %d!\n",adjst.nadj);
+		    exit(0);
+		  }
+		}
+	      }
+	    } else {
+	      printf(" %d",ver);
+	      (*adjs)[ver].nadj = -1; /* There are unreal vertices here */
+	    }
+	  }
+	}
+      } else (*adjs)[ver].nadj = -2;
+
+      /* Enumerate the unique adjacencies*/
+      if (adjst.nadj >= 4) {
+	qsort((void *)adjst.adj, adjst.nadj, sizeof(int), &compar);
+	count = 1;
+
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    if (adjst.adj[i] >= nvpbuf && !errorReported) {
+        errorReported = 1;
+	      printf("Guard point encountered.  Increase border and/or nguard.\n");
+	      printf("P:(%f,%f,%f), G: (%f,%f,%f)\n",points[3*ver],points[3*ver+1],points[3*ver+2],
+		     points[3*adjst.adj[i]],points[3*adjst.adj[i]+1],points[3*adjst.adj[i]+2]);
+	    }
+	    count++;
+	  }
+	(*adjs)[ver].adj = (int *)malloc(count*sizeof(int));
+	if ((*adjs)[ver].adj == NULL) {
+	  printf("Unable to allocate (*adjs)[ver].adj\n");
+	  exit(0);
+	}
+	(*adjs)[ver].adj[0] = adjst.adj[0];
+	count = 1;
+	for (i=1; i<adjst.nadj; i++)
+	  if (adjst.adj[i] != adjst.adj[i-1]) {
+	    (*adjs)[ver].adj[count] = adjst.adj[i];
+	    count++;
+	  }
+	(*adjs)[ver].nadj = count;
+      } else {
+	printf("Number of adjacencies %d < 4, particle %d -> %d\n",adjst.nadj,ver,ver);
+	exit(0);
+      }
+      ver++;
+      if (ver == nvp) break;
+    }
+    qh_settempfree (&coplanar_points);
+    qh_settempfree (&vertex_points);
+    qh_settempfree (&vertices);
+  }
+  qh_freeqhull(!qh_ALL);                 /* free long memory */
+  qh_memfreeshort (&curlong, &totlong);  /* free short memory and memory allocator */
+  if (curlong || totlong) 
+    fprintf (errfile, "qhull internal warning (delaunadj): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  free(adjst.adj);
+  return exitcode;
+}
+  
+/* Calculates the Voronoi volume from a set of Delaunay adjacencies */
+int vorvol (coordT *deladjs, coordT *points, pointT *intpoints, int numpoints, float *vol) {
+  int dim= 3;	            /* dimension of points */
+  boolT ismalloc= False;    /* True if qhull should free points in qh_freeqhull() or reallocation */
+  char flags[250];          /* option flags for qhull, see qh_opt.htm */
+  FILE *outfile= NULL;    /* output from qh_produce_output()
+			       use NULL to skip qh_produce_output() */
+  FILE *errfile= stderr;    /* error messages from qhull code */
+  int exitcode;             /* 0 if no error from qhull */
+  facetT *facet;	    /* set by FORALLfacets */
+  int curlong, totlong;	    /* memory remaining after qh_memfreeshort */
+
+  coordT *point, *normp, *coordp, *feasiblep, *deladj;
+  int i, j, k;
+  boolT zerodiv;
+  float runsum;
+  char region;
+  /*coordT *points;
+    pointT *intpoints;*/
+  
+  /* make point array from adjacency coordinates (add offset)*/
+  /*points = (coordT *)malloc(4*numpoints*sizeof(coordT));
+  if (points == NULL) {
+    printf("Unable to allocate points\n");
+    exit(0);
+    }*/
+  for (i=0; i<numpoints; i++) {
+    runsum = 0.;
+    deladj = deladjs + 3*i;
+    point = points + 4*i;
+    for (j=0;j<3;j++) {
+      runsum += deladj[j]*deladj[j];
+      point[j] = deladj[j];
+    }
+    point[3] = -0.5*runsum;
+  }
+  sprintf (flags, "qhull H0");
+
+  exitcode= qh_new_qhull (4, numpoints, points, ismalloc,
+			  flags, outfile, errfile);
+  
+  numpoints = 0;
+  if (!exitcode) {                  /* if no error */
+    FORALLfacets {
+      numpoints++;
+    }
+    /* Now we know how many points */
+    /*intpoints = (pointT *)malloc(dim*numpoints*sizeof(pointT));
+    if (intpoints == NULL) {
+    printf("Unable to allocate intpoints\n");
+    exit(0);
+    }*/
+
+    j = 0;
+    FORALLfacets {
+      if (!qh feasible_point) {
+	fprintf (stdout, "qhull input error (qh_printafacet): option 'Fp' needs qh feasible_point\n");
+	qh_errexit( qh_ERRinput, NULL, NULL);
+      }
+      point= coordp= intpoints + j*3;
+      j++;
+      normp= facet->normal;
+      feasiblep= qh feasible_point;
+      if (facet->offset < -qh MINdenom) {
+	for (k= qh hull_dim; k--; )
+	  *(coordp++)= (*(normp++) / - facet->offset) + *(feasiblep++);
+      }else {
+	for (k= qh hull_dim; k--; ) {
+	  *(coordp++)= qh_divzero (*(normp++), facet->offset, qh MINdenom_1,
+				   &zerodiv) + *(feasiblep++);
+	  if (zerodiv) {
+	    qh_memfree (point, qh normal_size);
+	    printf("LABELprintinfinite\n");
+	    exit(0);
+	  }
+	}
+      }
+    }
+  }
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+
+  /* Now we calculate the volume */
+  sprintf (flags, "qhull FA");
+  exitcode= qh_new_qhull (dim, numpoints, intpoints, ismalloc,
+			  flags, outfile, errfile);
+  
+  qh_getarea(qh facet_list);
+  *vol = qh totvol;
+  
+  qh_freeqhull (!qh_ALL);
+  qh_memfreeshort (&curlong, &totlong);
+  if (curlong || totlong) 
+    fprintf (errfile, "qhull internal warning (vorvol): did not free %d bytes of long memory (%d pieces)\n", totlong, curlong);
+  /*free(points); free(intpoints);*/
+  
+  return exitcode;
+}

python_tools/void_python_tools/backend/launchers.py

@@ -276,7 +276,7 @@ def launchZobov(sample, binPath, zobovDir=None, logDir=None, continueRun=None,
     cmd = "./%s >> %s 2>&1" % (vozScript, logFile)
     os.system(cmd)
 
-    cmd = "%s/../c_tools/zobov2/jozov2 %s %s %s %s %s %g %s >> %s 2>&1" % \
+    cmd = "%s/../c_tools/zobov2/jozov2/jozov2 %s %s %s %s %s %g %s >> %s 2>&1" % \
           (binPath, \
            zobovDir+"/adj_"+sampleName+".dat", \
            zobovDir+"/vol_"+sampleName+".dat", \