diff --git a/include/IECore/KDTree.h b/include/IECore/KDTree.h
index a9969cc371..87fc52d4b2 100644
--- a/include/IECore/KDTree.h
+++ b/include/IECore/KDTree.h
@@ -36,10 +36,12 @@
 #define IE_CORE_KDTREE_H
 
 #include "IECore/Export.h"
+#include "IECore/TypeTraits.h"
 #include "IECore/VectorTraits.h"
 
 IECORE_PUSH_DEFAULT_VISIBILITY
 #include "Imath/ImathVec.h"
+#include "Imath/ImathBox.h"
 IECORE_POP_DEFAULT_VISIBILITY
 
 #include <set>
@@ -48,6 +50,15 @@ IECORE_POP_DEFAULT_VISIBILITY
 namespace IECore
 {
 
+namespace Detail
+{
+	// \todo - ugly machinery needed until we deprecate the old signature of enclosedPoints
+	template <typename, typename = void>
+	struct IsIterator : std::false_type {};
+	template <typename T>
+	struct IsIterator<T, std::void_t< typename std::iterator_traits<T>::iterator_category >> : std::true_type {};
+}
+
 /// The KDTree class provides accelerated searching of pointsets. It is
 /// templated so that it can operate on a wide variety of datatypes, and uses
 /// the VectorTraits.h and VectorOps.h functionality to assist in this.
@@ -93,21 +104,28 @@ class KDTree
 		/// \threading May be called by multiple concurrent threads.
 		PointIterator nearestNeighbour( const Point &p, BaseType &distSquared ) const;
 
-		/// Populates the passed vector of iterators with the neighbours of point p which are closer than radius r. Returns the number of points found.
-		/// \todo There should be a form where nearNeighbours is an output iterator, to allow any container to be filled.
-		/// See enclosedPoints for an example of this form.
-		/// \threading May be called by multiple concurrent threads provided they are each using a different vector for the result.
-		unsigned int nearestNeighbours( const Point &p, BaseType r, std::vector<PointIterator> &nearNeighbours ) const;
+		/// Call a functor for each neighbour of point p which is closer than radius r.
+		/// The functor must take a PointIterator.
+		template<typename F>
+		void nearestNeighbours( const Point &p, BaseType r, F &&functor ) const;
 
 		class Neighbour;
 		/// Populates the passed vector with the N closest neighbours to p, sorted with the closest first. Returns the number found.
 		/// \threading May be called by multiple concurrent threads provided they are each using a different vector for the result.
 		unsigned int nearestNNeighbours( const Point &p, unsigned int numNeighbours, std::vector<Neighbour> &nearNeighbours ) const;
 
-		/// Finds all the points contained by the specified bound, outputting them to the specified iterator.
+		/// Finds all the points contained by the specified bound, outputting them to the specified functor,
+		/// which must take a PointIterator.
 		/// \threading May be called by multiple concurrent threads.
-		template<typename Box, typename OutputIterator>
-		void enclosedPoints( const Box &bound, OutputIterator it ) const;
+		template<typename Box, typename F>
+		void enclosedPoints( const Box &bound, F &&functor ) const;
+
+		// Finds all the points contained within a set of half-spaces, passing them to the given
+		// functor which must take a PointIterator.
+		// A half-space is specified with an origin and a plane normal ( the normal points towards the region that
+		// is included )
+		template<typename F>
+		void enclosedPoints( const std::vector<Point> &normals, const std::vector<Point> &origins, F &&functor ) const;
 
 		/// Returns the number of nodes in the tree.
 		inline NodeIndex numNodes() const;
@@ -132,15 +150,30 @@ class KDTree
 
 		class AxisSort;
 
+		// -- Utilities used when building the tree	--
+
+		// Compute min/max of a list of points
+		Imath::Box<Point> bound( PermutationConstIterator permFirst, PermutationConstIterator permLast );
+		// Return which axis of the bounding box is largest
 		unsigned char majorAxis( PermutationConstIterator permFirst, PermutationConstIterator permLast );
+		// Recursively build the tree
 		void build( NodeIndex nodeIndex, PermutationIterator permFirst, PermutationIterator permLast );
 
+
+		// -- Walk functions that implement the recursive searches --
+
 		void nearestNeighbourWalk( NodeIndex nodeIndex, const Point &p, PointIterator &closestPoint, BaseType &distSquared ) const;
 
-		void nearestNeighboursWalk( NodeIndex nodeIndex, const Point &p, BaseType r2, std::vector<PointIterator> &nearNeighbours ) const;
+		template<typename F>
+		void nearestNeighboursWalk( NodeIndex nodeIndex, const Point &p, BaseType r2, F &&functor ) const;
+
+		template<typename Box, typename F>
+		void enclosedPointsWalk( NodeIndex nodeIndex, const Box &bound, F &&functor  ) const;
+
+		struct HalfSpaceWorkingData;
 
-		template<typename Box, typename OutputIterator>
-		void enclosedPointsWalk( NodeIndex nodeIndex, const Box &bound, OutputIterator it ) const;
+		template<typename F>
+		void enclosedPointsHalfSpacesWalk( NodeIndex nodeIndex, std::vector<HalfSpaceWorkingData> &working, F &&functor ) const;
 
 		void nearestNNeighboursWalk( NodeIndex nodeIndex, const Point &p, unsigned int numNeighbours, std::vector<Neighbour> &nearNeighbours, BaseType &maxDistSquared ) const;
 
@@ -148,6 +181,7 @@ class KDTree
 		NodeVector m_nodes;
 		int m_maxLeafSize;
 		PointIterator m_lastPoint;
+		Imath::Box< Point > m_bound;
 
 };
 
diff --git a/include/IECore/KDTree.inl b/include/IECore/KDTree.inl
index 73b3b71f6c..4a63a69284 100644
--- a/include/IECore/KDTree.inl
+++ b/include/IECore/KDTree.inl
@@ -33,9 +33,11 @@
 //////////////////////////////////////////////////////////////////////////
 
 #include "IECore/BoxOps.h"
+#include "IECore/Exception.h"
 #include "IECore/VectorOps.h"
 
 #include <algorithm>
+#include "boost/container/small_vector.hpp"
 
 namespace IECore
 {
@@ -108,6 +110,14 @@ class KDTree<PointIterator>::AxisSort
 		const unsigned int m_axis;
 };
 
+template<class PointIterator>
+struct KDTree<PointIterator>::HalfSpaceWorkingData
+{
+	Point normal;
+	Point currentInnermost;
+	BaseType threshold;
+};
+
 // initialisation
 
 template<class PointIterator>
@@ -122,7 +132,7 @@ KDTree<PointIterator>::KDTree( PointIterator first, PointIterator last, int maxL
 }
 
 template<class PointIterator>
-void KDTree<PointIterator>::init( PointIterator first, PointIterator last, int maxLeafSize  )
+void KDTree<PointIterator>::init( PointIterator first, PointIterator last, int maxLeafSize )
 {
 	m_maxLeafSize = maxLeafSize;
 	m_lastPoint = last;
@@ -135,32 +145,44 @@ void KDTree<PointIterator>::init( PointIterator first, PointIterator last, int m
 
 	/// \todo Can we reserve() enough space for m_nodes before doing this?
 	build( rootIndex(), m_perm.begin(), m_perm.end() );
+
+	// The cut planes we store only limit the size of each Node within the interior of the KDTree.
+	// If we need accurate sizes for Nodes on the exterior of the tree ( rather than treating them as
+	// infinite ), we need to include the bound as well ( this can be particularly important when the
+	// data is an axis-aligned plane, where every node on the "exterior" in the Z axis ).
+	m_bound = bound( m_perm.begin(), m_perm.end() );
 }
 
 template<class PointIterator>
-unsigned char KDTree<PointIterator>::majorAxis( PermutationConstIterator permFirst, PermutationConstIterator permLast )
+Imath::Box<typename KDTree<PointIterator>::Point> KDTree<PointIterator>::bound( PermutationConstIterator permFirst, PermutationConstIterator permLast )
 {
-	Point min, max;
+	Imath::Box<Point> result;
 	for( unsigned char i=0; i<VectorTraits<Point>::dimensions(); i++ ) {
-		min[i] = std::numeric_limits<BaseType>::max();
-		max[i] = std::numeric_limits<BaseType>::lowest();
+		result.min[i] = std::numeric_limits<BaseType>::max();
+		result.max[i] = std::numeric_limits<BaseType>::lowest();
 	}
 	for( PermutationConstIterator it=permFirst; it!=permLast; it++ )
 	{
 		for( unsigned char i=0; i<VectorTraits<Point>::dimensions(); i++ )
 		{
-			if( (**it)[i] < min[i] )
+			if( (**it)[i] < result.min[i] )
 			{
-				min[i] = (**it)[i];
+				result.min[i] = (**it)[i];
 			}
-			if( (**it)[i] > max[i] )
+			if( (**it)[i] > result.max[i] )
 			{
-				max[i] = (**it)[i];
+				result.max[i] = (**it)[i];
 			}
 		}
 	}
+	return result;
+}
+
+template<class PointIterator>
+unsigned char KDTree<PointIterator>::majorAxis( PermutationConstIterator permFirst, PermutationConstIterator permLast )
+{
 	unsigned char major = 0;
-	Point size = max - min;
+	Point size = bound( permFirst, permLast ).size();
 	for( unsigned char i=1; i<VectorTraits<Point>::dimensions(); i++ )
 	{
 		if( size[i] > size[major] )
@@ -183,7 +205,7 @@ void KDTree<PointIterator>::build( NodeIndex nodeIndex, PermutationIterator perm
 	if( permLast - permFirst > m_maxLeafSize )
 	{
 		unsigned int cutAxis = majorAxis( permFirst, permLast );
-		PermutationIterator permMid = permFirst  + (permLast - permFirst)/2;
+		PermutationIterator permMid = permFirst + (permLast - permFirst)/2;
 		std::nth_element( permFirst, permMid, permLast, AxisSort( cutAxis ) );
 		BaseType cutValue = (**permMid)[cutAxis];
 		// insert node
@@ -219,20 +241,45 @@ PointIterator KDTree<PointIterator>::nearestNeighbour( const Point &p, BaseType
 }
 
 template<class PointIterator>
-unsigned int KDTree<PointIterator>::nearestNeighbours( const Point &p, BaseType r, std::vector<PointIterator> &nearNeighbours ) const
+template<typename F>
+void KDTree<PointIterator>::nearestNeighbours( const Point &p, BaseType r, F &&functor ) const
 {
-	nearNeighbours.clear();
-
-	nearestNeighboursWalk(rootIndex(), p, r*r, nearNeighbours );
+	nearestNeighboursWalk(rootIndex(), p, r*r, functor );
+}
 
-	return nearNeighbours.size();
+template<class PointIterator>
+template<typename Box, typename F >
+void KDTree<PointIterator>::enclosedPoints( const Box &bound, F &&functor ) const
+{
+	enclosedPointsWalk( rootIndex(), bound, functor );
 }
 
 template<class PointIterator>
-template<typename Box, typename OutputIterator>
-void KDTree<PointIterator>::enclosedPoints( const Box &bound, OutputIterator it ) const
+template<typename F>
+void KDTree<PointIterator>::enclosedPoints(
+	const std::vector<Point> &normals, const std::vector<Point> &origins, F &&functor
+) const
 {
-	enclosedPointsWalk( rootIndex(), bound, it );
+	std::vector<HalfSpaceWorkingData> workingData;
+	if( normals.size() != origins.size() )
+	{
+		throw IECore::Exception( "Mismatched normals and origins passed to enclosedPoints" );
+	}
+
+	workingData.resize( normals.size() );
+	for( size_t i = 0; i < normals.size(); i++ )
+	{
+		workingData[i].normal = normals[i];
+		workingData[i].threshold = normals[i].dot( origins[i] );
+
+		
+		for( unsigned char j=0; j<VectorTraits<Point>::dimensions(); j++ )
+		{
+			workingData[i].currentInnermost[j] = std::max( normals[i][j] * m_bound.min[j], normals[i][j] * m_bound.max[j] );
+		}
+	}
+
+	enclosedPointsHalfSpacesWalk( rootIndex(), workingData, functor );
 }
 
 template<class PointIterator>
@@ -294,7 +341,8 @@ void KDTree<PointIterator>::nearestNeighbourWalk( NodeIndex nodeIndex, const Poi
 }
 
 template<class PointIterator>
-void KDTree<PointIterator>::nearestNeighboursWalk( NodeIndex nodeIndex, const Point &p, BaseType r2, std::vector<PointIterator> &nearNeighbours ) const
+template<typename F>
+void KDTree<PointIterator>::nearestNeighboursWalk( NodeIndex nodeIndex, const Point &p, BaseType r2, F &&functor ) const
 {
 	const Node &node = m_nodes[nodeIndex];
 	if( node.isLeaf() )
@@ -307,7 +355,7 @@ void KDTree<PointIterator>::nearestNeighboursWalk( NodeIndex nodeIndex, const Po
 
 			if (dist2 < r2 )
 			{
-				nearNeighbours.push_back( *perm );
+				functor( *perm );
 			}
 		}
 	}
@@ -327,10 +375,10 @@ void KDTree<PointIterator>::nearestNeighboursWalk( NodeIndex nodeIndex, const Po
 			secondChild = highChildIndex( nodeIndex );
 		}
 
-		nearestNeighboursWalk( firstChild, p, r2, nearNeighbours );
+		nearestNeighboursWalk( firstChild, p, r2, functor );
 		if( d*d < r2 )
 		{
-			nearestNeighboursWalk( secondChild, p, r2, nearNeighbours );
+			nearestNeighboursWalk( secondChild, p, r2, functor );
 		}
 	}
 }
@@ -398,8 +446,8 @@ void KDTree<PointIterator>::nearestNNeighboursWalk( NodeIndex nodeIndex, const P
 }
 
 template<class PointIterator>
-template<typename Box, typename OutputIterator>
-void KDTree<PointIterator>::enclosedPointsWalk( NodeIndex nodeIndex, const Box &bound, OutputIterator it ) const
+template<typename Box, typename F>
+void KDTree<PointIterator>::enclosedPointsWalk( NodeIndex nodeIndex, const Box &bound, F &&functor ) const
 {
 	const Node &node = m_nodes[nodeIndex];
 
@@ -411,7 +459,7 @@ void KDTree<PointIterator>::enclosedPointsWalk( NodeIndex nodeIndex, const Box &
 			const Point &pp = **perm;
 			if( boxIntersects( bound, pp ) )
 			{
-				*it++ = *perm;
+				functor( *perm );
 			}
 		}
 	}
@@ -419,11 +467,97 @@ void KDTree<PointIterator>::enclosedPointsWalk( NodeIndex nodeIndex, const Box &
 	{
 		if( vecGet( BoxTraits<Box>::min( bound ), node.cutAxis() ) <= node.cutValue() )
 		{
-			enclosedPointsWalk( lowChildIndex( nodeIndex ), bound, it );
+			enclosedPointsWalk( lowChildIndex( nodeIndex ), bound, functor );
 		}
 		if( vecGet( BoxTraits<Box>::max( bound ), node.cutAxis() ) >= node.cutValue() )
 		{
-			enclosedPointsWalk( highChildIndex( nodeIndex ), bound, it );
+			enclosedPointsWalk( highChildIndex( nodeIndex ), bound, functor );
+		}
+	}
+}
+
+template<class PointIterator>
+template<typename F>
+void KDTree<PointIterator>::enclosedPointsHalfSpacesWalk( NodeIndex nodeIndex, std::vector<HalfSpaceWorkingData> &working, F &&functor ) const
+{
+	// TODO - normalize somewhere before this?
+	const Node &node = m_nodes[nodeIndex];
+
+	if( node.isLeaf() )
+	{
+		PointIterator *permLast = node.permLast();
+		for( PointIterator *perm = node.permFirst(); perm!=permLast; perm++ )
+		{
+			const Point &pp = **perm;
+			bool reject = false;
+			for( HalfSpaceWorkingData &halfSpace : working )
+			{
+				reject |= pp.dot( halfSpace.normal ) < halfSpace.threshold;
+			}
+
+			if( !reject )
+			{
+				functor( *perm );
+			}
+		}
+	}
+	else
+	{
+		unsigned char cutAxis = node.cutAxis();
+		BaseType cutValue = node.cutValue();
+
+		boost::container::small_vector<BaseType, 6> restoreInnermost;
+		restoreInnermost.reserve( working.size() );
+
+		for( HalfSpaceWorkingData &halfSpace : working )
+		{
+			restoreInnermost.push_back( halfSpace.currentInnermost[ cutAxis ] );
+		}
+
+		bool rejectLow = false;
+		for( HalfSpaceWorkingData &halfSpace : working )
+		{
+			if( halfSpace.normal[ cutAxis ] > BaseType( 0 ) )
+			{
+				halfSpace.currentInnermost[ cutAxis ] = halfSpace.normal[ cutAxis ] * cutValue;
+				if( vecSumElements( halfSpace.currentInnermost ) < halfSpace.threshold )
+				{
+					rejectLow = true;
+				}
+			}
+		}
+
+		if( !rejectLow )
+		{
+			enclosedPointsHalfSpacesWalk( lowChildIndex( nodeIndex ), working, functor );
+		}
+
+		bool rejectHigh = false;
+		for( size_t i = 0; i < working.size(); i++ )
+		{
+			HalfSpaceWorkingData &halfSpace = working[i];
+			if( halfSpace.normal[ cutAxis ] < BaseType( 0 ) )
+			{
+				halfSpace.currentInnermost[ cutAxis ] = halfSpace.normal[ cutAxis ] * cutValue;
+				if( vecSumElements( halfSpace.currentInnermost ) < halfSpace.threshold )
+				{
+					rejectHigh = true;
+				}
+			}
+			else
+			{
+				halfSpace.currentInnermost[ cutAxis ] = restoreInnermost[i];
+			}
+		}
+
+		if( !rejectHigh )
+		{
+			enclosedPointsHalfSpacesWalk( highChildIndex( nodeIndex ), working, functor );
+		}
+
+		for( size_t i = 0; i < working.size(); i++ )
+		{
+			working[i].currentInnermost[ cutAxis ] = restoreInnermost[i];
 		}
 	}
 }
diff --git a/include/IECore/VectorOps.h b/include/IECore/VectorOps.h
index 4e47e40165..2f8ccb98e1 100644
--- a/include/IECore/VectorOps.h
+++ b/include/IECore/VectorOps.h
@@ -164,6 +164,9 @@ inline T vecConstruct( const typename VectorTraits<T>::BaseType *components );
 template<typename T>
 inline T vecCross( const T &v1, const T &v2 );
 
+template<typename T>
+inline typename VectorTraits<T>::BaseType vecSumElements( const T &v );
+
 } // namespace IECore
 
 #include "IECore/VectorOps.inl"
diff --git a/include/IECore/VectorOps.inl b/include/IECore/VectorOps.inl
index f8f3639951..91b0320c29 100644
--- a/include/IECore/VectorOps.inl
+++ b/include/IECore/VectorOps.inl
@@ -262,6 +262,21 @@ inline T vecCross( const T &v1, const T &v2 )
 
 }
 
+template<typename T>
+inline typename VectorTraits<T>::BaseType vecSumElements( const T &v )
+{
+	// TODO - the conventions for VectorTraits for 1 dimensional vectors are weird ...
+	// maybe instead of worrying about that, maybe I should just make this private?
+	if constexpr( T::dimensions() == 2 )
+	{
+		return v[0] + v[1];
+	}
+	else if constexpr( T::dimensions() == 3 )
+	{
+		return v[0] + v[1] + v[2];
+	}
+}
+
 } // namespace IECore
 
 
diff --git a/src/IECorePython/KDTreeBinding.cpp b/src/IECorePython/KDTreeBinding.cpp
index 37f0a7d6ac..b387a32598 100644
--- a/src/IECorePython/KDTreeBinding.cpp
+++ b/src/IECorePython/KDTreeBinding.cpp
@@ -75,10 +75,10 @@ struct KDTreeWrapper
 
 	PointDataPtr m_points;
 
-	KDTreeWrapper(PointDataPtr points)
+	KDTreeWrapper(PointDataPtr points, int maxLeafSize = 4)
 	{
 		m_points = points->copy();
-		m_tree = new T(m_points->readable().begin(), m_points->readable().end());
+		m_tree = new T(m_points->readable().begin(), m_points->readable().end(), maxLeafSize);
 	}
 
 	virtual ~KDTreeWrapper()
@@ -100,22 +100,15 @@ struct KDTreeWrapper
 	{
 		assert(m_tree);
 
-		typedef std::vector<typename T::Iterator> PointArray;
+		IntVectorDataPtr indicesData = new IntVectorData();
 
-		PointArray points;
+		auto &indices = indicesData->writable();
+		m_tree->nearestNeighbours(
+			p, r,
+			[&indices, this]( const typename T::Iterator i ){ indices.push_back( std::distance( m_points->readable().begin(), i ) ); }
+		);
 
-		unsigned int num = m_tree->nearestNeighbours(p, r, points);
-
-		IntVectorDataPtr indices = new IntVectorData();
-
-		indices->writable().reserve( num );
-
-		for (typename PointArray::const_iterator it = points.begin(); it != points.end(); ++it)
-		{
-			indices->writable().push_back(  std::distance( m_points->readable().begin(), *it ) );
-		}
-
-		return indices;
+		return indicesData;
 
 	}
 
@@ -144,22 +137,31 @@ struct KDTreeWrapper
 
 	IntVectorDataPtr enclosedPoints( const Box &bound )
 	{
-		typedef std::vector<typename T::Iterator> PointArray;
-
-		PointArray points;
+		IntVectorDataPtr indicesData = new IntVectorData();
 
-		m_tree->enclosedPoints( bound, std::back_insert_iterator<PointArray>( points ) );
+		auto &indices = indicesData->writable();
+		m_tree->enclosedPoints(
+			bound,
+			[&indices, this]( const typename T::Iterator i ){ indices.push_back( std::distance( m_points->readable().begin(), i ) ); }
+		);
 
-		IntVectorDataPtr indices = new IntVectorData();
+		return indicesData;
+	}
 
-		indices->writable().reserve( points.size() );
+	IntVectorDataPtr enclosedPointsWithHalfSpaces(
+		const typename TypedData< std::vector<typename T::Point> >::ConstPtr &normals,
+		const typename TypedData< std::vector<typename T::Point> >::ConstPtr &origins
+	)
+	{
+		IntVectorDataPtr indicesData = new IntVectorData();
 
-		for (typename PointArray::const_iterator it = points.begin(); it != points.end(); ++it)
-		{
-			indices->writable().push_back(  std::distance( m_points->readable().begin(), *it ) );
-		}
+		auto &indices = indicesData->writable();
+		m_tree->enclosedPoints(
+			normals->readable(), origins->readable(),
+			[&indices, this]( const typename T::Iterator i ){ indices.push_back( std::distance( m_points->readable().begin(), i ) ); }
+		);
 
-		return indices;
+		return indicesData;
 	}
 
 };
@@ -170,10 +172,12 @@ void bindKDTree(const char *bindName)
 {
 	class_<KDTreeWrapper<T>, boost::noncopyable>(bindName, no_init)
 		.def(init< typename KDTreeWrapper<T>::PointDataPtr >() )
+		.def(init< typename KDTreeWrapper<T>::PointDataPtr, int >() )
 		.def("nearestNeighbour", &KDTreeWrapper<T>::nearestNeighbour )
 		.def("nearestNeighbours", &KDTreeWrapper<T>::nearestNeighbours )
 		.def("nearestNNeighbours", &KDTreeWrapper<T>::nearestNNeighbours )
 		.def("enclosedPoints", &KDTreeWrapper<T>::enclosedPoints )
+		.def("enclosedPoints", &KDTreeWrapper<T>::enclosedPointsWithHalfSpaces )
 		;
 }
 
diff --git a/test/IECore/KDTree.py b/test/IECore/KDTree.py
index c32362568a..399ea32784 100644
--- a/test/IECore/KDTree.py
+++ b/test/IECore/KDTree.py
@@ -61,15 +61,14 @@ def doNearestNeighbours(self, numPoints):
 
 		for i in range(0, numPoints):
 			for r in self.radii:
-				pIdxArray = self.tree.nearestNeighbours( self.points[i], r )
+				result = self.tree.nearestNeighbours( self.points[i], r )
 
-				for pIdx in pIdxArray:
-					self.assertTrue(pIdx >= 0)
-					self.assertTrue(pIdx < numPoints)
-					nearestPt = self.points[pIdx]
-					distToNearest = (self.points[pIdx] - self.points[i]).length()
+				expectedResult = []
+				for j in range( self.points.size() ) :
+					if ( self.points[j] - self.points[i] ).length() < r:
+						expectedResult.append( j )
 
-					self.assertTrue( distToNearest <= r )
+				self.assertEqual( sorted( result ), expectedResult )
 
 	def doNearestNNeighbours(self, numPoints):
 
@@ -129,6 +128,40 @@ def doEnclosedPoints( self, numPoints ) :
 				else :
 					self.assertFalse( i in s )
 
+	def doEnclosedPointsHalfSpaces( self, numPoints ) :
+
+		self.makeTree( numPoints )
+
+		for i in range( 0, 50 ) :
+			if i < 10:
+				# Test with a single random halfspace
+				normals = type( self.points )( [ (2 * self.randomVec() - 1) ] )
+				origins = type( self.points )( [ self.randomVec() ] )
+			elif i < 20:
+				# Test with 6 halfspaces that are guaranteed to contain some actual volume
+				# ( since they are chosen to be tangent to the surface of a sphere )
+				normals = type( self.points )( [ (2 * self.randomVec() - 1).normalized() for i in range( 6 ) ] )
+				origins = type( self.points )( [ -n * 0.25 + 0.5 for n in normals ] )
+			else:
+				# Test with three totally random halfspaces that may or may not intersect
+				# to contain anything
+				normals = type( self.points )( [ (2 * self.randomVec() - 1) for i in range( 3 ) ] )
+				origins = type( self.points )( [ self.randomVec() for i in range( 3 ) ] )
+
+			result = self.tree.enclosedPoints( normals, origins )
+
+			expectedResult = []
+			for i in range( self.points.size() ) :
+				reject = False
+				for n,o in zip( normals, origins ):
+					if ( self.points[i] - o ).dot( n ) < 0:
+						reject = True
+
+				if not reject:
+					expectedResult.append( i )
+
+			self.assertEqual( sorted( result ), expectedResult )
+
 
 class TestKDTreeV2f(unittest.TestCase, TestKDTree):
 
@@ -148,6 +181,9 @@ def randomBox( self ) :
 		max = min + imath.V2f( random.random(), random.random() )
 		return imath.Box2f( min, max )
 
+	def randomVec( self ):
+		return imath.V2f( random.random(), random.random() )
+
 	def testConstructors(self):
 		"""Test KDTreeV2f constructors"""
 
@@ -178,6 +214,12 @@ def testEnclosedPoints(self):
 		for t in self.treeSizes:
 			self.doEnclosedPoints(t)
 
+	def testEnclosedPointsHalfSpaces(self):
+		"""Test KDTreeV2f enclosedPointsHalfSpaces"""
+
+		for t in self.treeSizes:
+			self.doEnclosedPointsHalfSpaces(t)
+
 class TestKDTreeV2d(unittest.TestCase, TestKDTree):
 
 	def makeTree(self, numPoints):
@@ -196,6 +238,9 @@ def randomBox( self ) :
 		max = min + imath.V2d( random.random(), random.random() )
 		return imath.Box2d( min, max )
 
+	def randomVec( self ):
+		return imath.V2d( random.random(), random.random() )
+
 	def testConstructors(self):
 		"""Test KDTreeV2d constructors"""
 
@@ -226,6 +271,12 @@ def testEnclosedPoints(self):
 		for t in self.treeSizes:
 			self.doEnclosedPoints(t)
 
+	def testEnclosedPointsHalfSpaces(self):
+		"""Test KDTreeV2d enclosedPointsHalfSpaces"""
+
+		for t in self.treeSizes:
+			self.doEnclosedPointsHalfSpaces(t)
+
 class TestKDTreeV3f(unittest.TestCase, TestKDTree):
 
 	def makeTree(self, numPoints):
@@ -244,6 +295,9 @@ def randomBox( self ) :
 		max = min + imath.V3f( random.random(), random.random(), random.random() )
 		return imath.Box3f( min, max )
 
+	def randomVec( self ):
+		return imath.V3f( random.random(), random.random(), random.random() )
+
 	def testConstructors(self):
 		"""Test KDTreeV3f constructors"""
 
@@ -274,6 +328,12 @@ def testEnclosedPoints(self):
 		for t in self.treeSizes:
 			self.doEnclosedPoints(t)
 
+	def testEnclosedPointsHalfSpaces(self):
+		"""Test KDTreeV3f enclosedPointsHalfSpaces"""
+
+		for t in self.treeSizes:
+			self.doEnclosedPointsHalfSpaces(t)
+
 class TestKDTreeV3d(unittest.TestCase, TestKDTree):
 
 	def makeTree(self, numPoints):
@@ -292,6 +352,9 @@ def randomBox( self ) :
 		max = min + imath.V3d( random.random(), random.random(), random.random() )
 		return imath.Box3d( min, max )
 
+	def randomVec( self ):
+		return imath.V3d( random.random(), random.random(), random.random() )
+
 	def testConstructors(self):
 		"""Test KDTreeV3d constructors"""
 
@@ -322,6 +385,12 @@ def testEnclosedPoints(self):
 		for t in self.treeSizes:
 			self.doEnclosedPoints(t)
 
+	def testEnclosedPointsHalfSpaces(self):
+		"""Test KDTreeV3d enclosedPointsHalfSpaces"""
+
+		for t in self.treeSizes:
+			self.doEnclosedPointsHalfSpaces(t)
+
 
 if __name__ == "__main__":
 	unittest.main()
diff --git a/test/IECore/KDTreeTest.inl b/test/IECore/KDTreeTest.inl
index a1ae16eef6..4d0e186d4b 100644
--- a/test/IECore/KDTreeTest.inl
+++ b/test/IECore/KDTreeTest.inl
@@ -87,9 +87,12 @@ void KDTreeTest<T>::testNearestNeighbours()
 	for( typename Tree::Iterator it=m_points.begin(); it!=m_points.end(); it++ )
 	{
 		typename T::BaseType radius = 0.05;
-		unsigned int numNeighbours = m_tree->nearestNeighbours( *it, radius, nearNeighbours );
+		nearNeighbours.clear();
+		m_tree->nearestNeighbours( *it, radius, [&nearNeighbours]( const typename Tree::Iterator &i ){
+			nearNeighbours.push_back( i );
+		} );
 
-		BOOST_CHECK(numNeighbours <= m_numPoints);
+		BOOST_CHECK(nearNeighbours.size() <= m_numPoints);
 
 		typename IteratorVector::const_iterator nit = nearNeighbours.begin();
 		for (; nit != nearNeighbours.end(); ++nit)