From 1ae72331397c9953c6908d24212ef8ad09b7248f Mon Sep 17 00:00:00 2001
From: Guilhem Lavaux <lavaux@illinois.edu>
Date: Wed, 22 Jun 2011 13:57:06 -0500
Subject: [PATCH] Have a recursion that only count particles

---
 src/mykdtree.hpp |  2 ++
 src/mykdtree.tcc | 41 +++++++++++++++++++++++++++++++----------
 src/yorick.cpp   |  6 ++++--
 3 files changed, 37 insertions(+), 12 deletions(-)

diff --git a/src/mykdtree.hpp b/src/mykdtree.hpp
index b695ada..d661376 100644
--- a/src/mykdtree.hpp
+++ b/src/mykdtree.hpp
@@ -97,6 +97,7 @@ namespace CosmoTool {
 			     CoordType *distances,
 			     uint32_t numCells)
       throw (NotEnoughCells);
+    uint32_t countCells(const coords& x, CoordType r);
 
     Cell *getNearestNeighbour(const coords& x);
 
@@ -146,6 +147,7 @@ namespace CosmoTool {
 		    coords minBound,
 		    coords maxBound);
     
+    template<bool justCount>
     void recursiveIntersectionCells(RecursionInfoCells<N,ValType, CType>& info,
 				    Node *node,
 				    int level)
diff --git a/src/mykdtree.tcc b/src/mykdtree.tcc
index 2a233a4..9899d60 100644
--- a/src/mykdtree.tcc
+++ b/src/mykdtree.tcc
@@ -92,7 +92,7 @@ namespace CosmoTool {
     info.numCells = numCells;
     info.distances = 0;
 
-    recursiveIntersectionCells(info, root, 0);
+    recursiveIntersectionCells<false>(info, root, 0);
     return info.currentRank;
   }
 
@@ -113,11 +113,29 @@ namespace CosmoTool {
     info.numCells = numCells;
     info.distances = distances;
 
-    recursiveIntersectionCells(info, root, 0);
+    recursiveIntersectionCells<false>(info, root, 0);
     return info.currentRank;
   }
 
   template<int N, typename ValType, typename CType>
+  uint32_t KDTree<N,ValType,CType>::countCells(const coords& x, CoordType r)
+  {
+    RecursionInfoCells<N,ValType> info;
+
+    memcpy(info.x, x, sizeof(x));
+    info.r = r;
+    info.r2 = r*r;
+    info.cells = 0;
+    info.currentRank = 0;
+    info.numCells = 0;
+    info.distances = 0;
+
+    recursiveIntersectionCells<true>(info, root, 0);
+    return info.currentRank;
+  }
+
+  template<int N, typename ValType, typename CType>
+  template<bool justCount>
   void KDTree<N,ValType,CType>::recursiveIntersectionCells(RecursionInfoCells<N,ValType,CType>& info, 
 							   Node *node,
 							   int level)
@@ -135,11 +153,14 @@ namespace CosmoTool {
 	  }
 	if (d2 < info.r2)
 	  {
-	    if (info.currentRank == info.numCells)
-	      throw NotEnoughCells();
-	    info.cells[info.currentRank] = node->value;
-	    if (info.distances)
-	      info.distances[info.currentRank] = d2;
+	    if (!justCount)
+	      {
+		if (info.currentRank == info.numCells)
+		  throw NotEnoughCells();
+		info.cells[info.currentRank] = node->value;
+		if (info.distances)
+		  info.distances[info.currentRank] = d2;
+	      }
 	    info.currentRank++;
 	  }
       }
@@ -149,15 +170,15 @@ namespace CosmoTool {
 	((info.x[axis]-info.r) < node->value->coord[axis]))
       {
 	if (node->children[0] != 0)
-	  recursiveIntersectionCells(info, node->children[0],
+	  recursiveIntersectionCells<justCount>(info, node->children[0],
 				     level+1);
       }
     if (((info.x[axis]+info.r) > node->value->coord[axis]) &&
 	((info.x[axis]-info.r) < node->maxBound[axis]))
       {
 	if (node->children[1] != 0)
-	  recursiveIntersectionCells(info, node->children[1],
-				     level+1);
+	  recursiveIntersectionCells<justCount>(info, node->children[1],
+					       level+1);
       }
   }
 
diff --git a/src/yorick.cpp b/src/yorick.cpp
index 78c5faf..65ed7c5 100644
--- a/src/yorick.cpp
+++ b/src/yorick.cpp
@@ -236,8 +236,10 @@ namespace CosmoTool {
 	dimList[i] = edge[i];
 	fullSize *= edge[i];
       }
-    array = new T[fullSize];
-    v->get(array, edge);
+    if (fullSize != 0) {
+      array = new T[fullSize];
+      v->get(array, edge);
+    }
     delete[] edge;   
   }