Imported the simple3DFilter to the sample directory

sample/CMakeLists.txt

@@ -69,3 +69,6 @@ endif (ENABLE_SHARP AND SHARP_LIBRARY AND SHARP_INCLUDE_PATH AND EIGEN3_FOUND)
 add_executable(test_cosmopower test_cosmopower.cpp)
 target_link_libraries(test_cosmopower ${tolink})
 
+
+add_executable(simple3DFilter simple3DFilter.cpp)
+target_link_libraries(simple3DFilter ${tolink})

sample/simple3DFilter.cpp

@@ -0,0 +1,193 @@
+#include <cassert>
+#include <CosmoTool/yorick.hpp>
+#include <CosmoTool/sphSmooth.hpp>
+#include <CosmoTool/mykdtree.hpp>
+#include <CosmoTool/miniargs.hpp>
+
+using namespace std;
+using namespace CosmoTool;
+
+#define N_SPH 16
+
+struct VCoord{
+  float v[3];
+};
+
+template<int i>
+ComputePrecision getVelocity(const VCoord& v)
+{
+  return v.v[i];
+}
+
+ComputePrecision getUnity(const VCoord& v)
+{
+  return 1.0;
+}
+
+typedef SPHSmooth<VCoord> MySmooth;
+typedef MySmooth::SPHTree MyTree;
+typedef MyTree::Cell MyCell;
+
+int main(int argc, char **argv)
+{
+  char *fname1, *fname2;
+  double rLimit, boxsize, rLimit2, cx, cy, cz;
+  int Nres;
+
+  MiniArgDesc args[] = {
+    { "INPUT DATA1", &fname1, MINIARG_STRING },
+    { "RADIUS LIMIT", &rLimit, MINIARG_DOUBLE },
+    { "BOXSIZE", &boxsize, MINIARG_DOUBLE },
+    { "RESOLUTION", &Nres, MINIARG_INT },
+    { "CX", &cx, MINIARG_DOUBLE },
+    { "CY", &cy, MINIARG_DOUBLE },
+    { "CZ", &cz, MINIARG_DOUBLE },
+    { 0, 0, MINIARG_NULL }
+  };
+
+  if (!parseMiniArgs(argc, argv, args))
+    return 1;
+
+  float *v1_data;
+  uint32_t *dimList;
+  uint32_t rank;
+  uint32_t N1_points, N2_points;
+  int *bins = new int[Nres*Nres*Nres];
+
+  rLimit2 = rLimit*rLimit;
+
+  loadArray(fname1, v1_data, dimList, rank);
+  assert(rank == 2);
+  assert(dimList[1] == 6);
+
+  N1_points = dimList[0];
+  delete[] dimList;
+  
+  cout << "Got " << N1_points << " in the first file." << endl;
+
+  MyCell *allCells_1 = new MyCell[N1_points];
+  
+  for (long i = 0; i < Nres*Nres*Nres; i++)
+    bins[i] = 0;
+
+  cout << "Shuffling data in cells..." << endl;
+  for (int i = 0 ; i < N1_points; i++)
+    {
+      for (int j = 0; j < 3; j++)
+	allCells_1[i].coord[j] = v1_data[i*6 + j];
+      for (int k = 0; k < 3; k++)
+	allCells_1[i].val.pValue.v[k] = v1_data[i*6 + 3 + k];
+      allCells_1[i].active = true;
+      allCells_1[i].val.weight = 0.0;
+
+      long rx = floor((allCells_1[i].coord[0]+cx)*Nres/boxsize+0.5);
+      long ry = floor((allCells_1[i].coord[1]+cy)*Nres/boxsize+0.5);
+      long rz = floor((allCells_1[i].coord[2]+cz)*Nres/boxsize+0.5);
+      
+      if (rx < 0 || rx >= Nres || ry < 0 || ry >= Nres || rz < 0 || rz >= Nres)
+	  continue;
+	  
+      bins[rx + ry*Nres + rz*Nres*Nres]++;
+    }
+  delete[] v1_data;
+   uint32_t dims[3] = { Nres, Nres, Nres } ;
+   saveArray("num_in_cell.nc", bins, dims, 3);
+
+  cout << "Building trees..." << endl;
+  MyTree tree1(allCells_1, N1_points);
+
+  cout << "Creating smoothing filter..." << endl;
+  MySmooth smooth1(&tree1, N_SPH);
+
+  uint32_t outDimList[3] = { Nres, Nres, Nres };
+  uint32_t outDimList2[4] = { 3, Nres, Nres, Nres };
+  ProgressiveOutput<float> out_den_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("density.nc", outDimList, 3);
+  ProgressiveOutput<float> out_vel_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("v3d.nc", outDimList2, 4);
+  ProgressiveOutput<float> out_rad_1 = 
+    ProgressiveOutput<float>::saveArrayProgressive("rad.nc", outDimList, 3);
+ 
+  cout << "Weighing..." << endl;
+  for (int rz = 0; rz < Nres; rz++)
+   {
+      double pz = (rz)*boxsize/Nres-cz;
+
+      cout << rz << " / " << Nres << endl;
+      for (int ry = 0; ry < Nres; ry++)
+	{
+	  double py = (ry)*boxsize/Nres-cy;
+	  for (int rx = 0; rx < Nres; rx++)
+	    {
+	      double px = (rx)*boxsize/Nres-cx;
+	      
+	      MyTree::coords c = { px, py, pz };
+
+	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+	      if (r2 > rLimit2)
+		{
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + Nres*ry + Nres*Nres*rz];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+	      smooth1.addGridSite(c);
+	    }
+	}
+    }
+
+  
+  
+
+  cout << "Interpolating..." << endl;
+  for (int rz = 0; rz < Nres; rz++)
+    {
+      double pz = (rz)*boxsize/Nres-cz;
+
+      cout << rz << " / " << Nres << endl;
+      for (int ry = 0; ry < Nres; ry++)
+	{
+	  double py = (ry)*boxsize/Nres-cy;
+	  for (int rx = 0; rx < Nres; rx++)
+	    {
+	      double px = (rx)*boxsize/Nres-cx;
+	      
+	      MyTree::coords c = { px, py, pz };
+
+	      double r2 = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
+	      if (r2 > rLimit2)
+		{
+		  out_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_vel_1.put(0);
+		  out_den_1.put(0);
+		  out_rad_1.put(0);
+		  continue;
+		}
+
+	      uint32_t numInCell = bins[rx + ry*Nres + rz*Nres*Nres];
+	      if (numInCell > N_SPH)
+		smooth1.fetchNeighbours(c, numInCell);
+	      else
+		smooth1.fetchNeighbours(c);
+
+	      float val;
+
+	      out_rad_1.put(smooth1.getSmoothingLen());
+	      val = smooth1.computeSmoothedValue(c, getVelocity<0>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<1>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getVelocity<2>);
+	      out_vel_1.put(val);
+	      val = smooth1.computeSmoothedValue(c, getUnity);
+	      out_den_1.put(val);
+	    }
+	}
+    }
+
+  return 0;
+};