From 4ee05287c32b7e6b69df13716ff06983cb871aab Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Wed, 9 Nov 2022 21:55:22 +0100
Subject: [PATCH 01/13] Add static ListImg.wrap() constructor method that wraps
 an existing List

---
 .../java/net/imglib2/img/list/ListImg.java    | 20 ++++++++++++++++---
 1 file changed, 17 insertions(+), 3 deletions(-)

diff --git a/src/main/java/net/imglib2/img/list/ListImg.java b/src/main/java/net/imglib2/img/list/ListImg.java
index f3187ebdb..024e753bf 100644
--- a/src/main/java/net/imglib2/img/list/ListImg.java
+++ b/src/main/java/net/imglib2/img/list/ListImg.java
@@ -37,6 +37,7 @@
 import java.util.ArrayList;
 import java.util.Collection;
 
+import java.util.List;
 import net.imglib2.img.Img;
 import net.imglib2.type.Type;
 
@@ -61,7 +62,7 @@
  */
 public class ListImg< T > extends AbstractListImg< T >
 {
-	final private ArrayList< T > pixels;
+	final private List< T > pixels;
 
 	public ListImg( final long[] dim, final T type )
 	{
@@ -86,9 +87,9 @@ public ListImg( final long[] dim, final T type )
 	public ListImg( final Collection< T > collection, final long... dim )
 	{
 		super( dim );
-		
+
 		assert numPixels == collection.size() : "Dimensions do not match number of pixels.";
-		
+
 		pixels = new ArrayList< T >( ( int ) numPixels );
 		pixels.addAll( collection );
 	}
@@ -130,4 +131,17 @@ public ListImg< T > copy()
 		}
 		return new ListImg< T >( this.pixels, dimension );
 	}
+
+	private ListImg( final List< T > pixels, final boolean dummy, final long... dim )
+	{
+		super( dim );
+
+		assert numPixels == pixels.size() : "Dimensions do not match number of pixels.";
+		this.pixels = pixels;
+	}
+
+	public static < T > ListImg< T > wrap( List< T > pixels, final long... dim )
+	{
+		return new ListImg<>( pixels, false, dim );
+	}
 }

From 8657f304e469ed560f3cabfb2e61768c0430ea7f Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Sat, 19 Nov 2022 22:09:55 +0100
Subject: [PATCH 02/13] Add default getDistance() implementation for neighbor
 search interfaces

---
 .../net/imglib2/neighborsearch/KNearestNeighborSearch.java   | 5 ++++-
 .../net/imglib2/neighborsearch/NearestNeighborSearch.java    | 5 ++++-
 .../net/imglib2/neighborsearch/RadiusNeighborSearch.java     | 5 ++++-
 3 files changed, 12 insertions(+), 3 deletions(-)

diff --git a/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearch.java b/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearch.java
index 781253fb6..edb3f6c8e 100644
--- a/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearch.java
+++ b/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearch.java
@@ -89,7 +89,10 @@ public interface KNearestNeighborSearch< T > extends NearestNeighborSearch< T >
 	 * the last search and the <em>i</em><sup>th</sup> nearest neighbor, ordered
 	 * by square Euclidean distance.
 	 */
-	double getDistance( int i );
+	default double getDistance( final int i )
+	{
+		return Math.sqrt( getSquareDistance( i ) );
+	}
 
 	/**
 	 * Create a copy.
diff --git a/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearch.java b/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearch.java
index 43d02daea..e2888df56 100644
--- a/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearch.java
+++ b/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearch.java
@@ -81,7 +81,10 @@ public interface NearestNeighborSearch< T > extends EuclideanSpace
 	 * the last search and the nearest neighbor, ordered by square Euclidean
 	 * distance.
 	 */
-	double getDistance();
+	default double getDistance()
+	{
+		return Math.sqrt( getSquareDistance() );
+	}
 
 	/**
 	 * Create a copy.
diff --git a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
index d55689753..db7d6913d 100644
--- a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
+++ b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
@@ -106,5 +106,8 @@ public interface RadiusNeighborSearch< T > extends EuclideanSpace
 	 * Access the Euclidean distance between the reference location as used for
 	 * the last search and the <em>i</em><sup>th</sup> neighbor.
 	 */
-	double getDistance( int i );
+	default double getDistance( final int i )
+	{
+		return Math.sqrt( getSquareDistance( i ) );
+	}
 }

From cdfbb073468325100907f20def4570e721173f52 Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Wed, 18 Jan 2023 13:23:15 +0100
Subject: [PATCH 03/13] Add RadiusNeighborSearch.copy() method

for consistency, because the other NeighborSearch interfaces have copy()
methods
---
 .../net/imglib2/neighborsearch/RadiusNeighborSearch.java | 5 +++++
 .../neighborsearch/RadiusNeighborSearchOnKDTree.java     | 9 +++++++++
 2 files changed, 14 insertions(+)

diff --git a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
index db7d6913d..02aa94df8 100644
--- a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
+++ b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearch.java
@@ -110,4 +110,9 @@ default double getDistance( final int i )
 	{
 		return Math.sqrt( getSquareDistance( i ) );
 	}
+
+	/**
+	 * Create a copy.
+	 */
+	RadiusNeighborSearch< T > copy();
 }
diff --git a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
index 28e00b450..909b5402e 100644
--- a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
+++ b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
@@ -146,4 +146,13 @@ public double getDistance( final int i )
 	{
 		return Math.sqrt( resultPoints.get( i ).b );
 	}
+
+	@Override
+	public RadiusNeighborSearchOnKDTree< T > copy()
+	{
+		final RadiusNeighborSearchOnKDTree< T > copy = new RadiusNeighborSearchOnKDTree<>( tree );
+		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
+		copy.resultPoints.addAll( resultPoints );
+		return copy;
+	}
 }

From 0dbb9b6be677c005947f01c8b516d826a5b0fa20 Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Fri, 14 Jul 2023 14:24:36 +0200
Subject: [PATCH 04/13] Revise KDTree implementation

---
 src/main/java/net/imglib2/KDTree.java         | 741 ++++--------------
 src/main/java/net/imglib2/KDTreeNode.java     | 185 ++---
 .../java/net/imglib2/kdtree/KDTreeData.java   | 323 ++++++++
 .../java/net/imglib2/kdtree/KDTreeImpl.java   | 229 ++++++
 .../java/net/imglib2/kdtree/KDTreeUtils.java  | 573 ++++++++++++++
 .../kdtree/KNearestNeighborSearchImpl.java    | 125 +++
 .../kdtree/NearestNeighborSearchImpl.java     |  95 +++
 .../kdtree/RadiusNeighborSearchImpl.java      | 155 ++++
 .../KNearestNeighborSearchOnKDTree.java       | 102 +--
 .../NearestNeighborSearchOnKDTree.java        |  73 +-
 .../RadiusNeighborSearchOnKDTree.java         | 102 +--
 .../net/imglib2/kdtree/KDTreeBenchmark.java   | 186 +++++
 .../net/imglib2/kdtree/KDTreeImplTest.java    | 151 ++++
 13 files changed, 2144 insertions(+), 896 deletions(-)
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeData.java
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeImpl.java
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeUtils.java
 create mode 100644 src/main/java/net/imglib2/kdtree/KNearestNeighborSearchImpl.java
 create mode 100644 src/main/java/net/imglib2/kdtree/NearestNeighborSearchImpl.java
 create mode 100644 src/main/java/net/imglib2/kdtree/RadiusNeighborSearchImpl.java
 create mode 100644 src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
 create mode 100644 src/test/java/net/imglib2/kdtree/KDTreeImplTest.java

diff --git a/src/main/java/net/imglib2/KDTree.java b/src/main/java/net/imglib2/KDTree.java
index 46bd1da70..6a95e53b8 100644
--- a/src/main/java/net/imglib2/KDTree.java
+++ b/src/main/java/net/imglib2/KDTree.java
@@ -1,177 +1,35 @@
-/*
- * #%L
- * ImgLib2: a general-purpose, multidimensional image processing library.
- * %%
- * Copyright (C) 2009 - 2023 Tobias Pietzsch, Stephan Preibisch, Stephan Saalfeld,
- * John Bogovic, Albert Cardona, Barry DeZonia, Christian Dietz, Jan Funke,
- * Aivar Grislis, Jonathan Hale, Grant Harris, Stefan Helfrich, Mark Hiner,
- * Martin Horn, Steffen Jaensch, Lee Kamentsky, Larry Lindsey, Melissa Linkert,
- * Mark Longair, Brian Northan, Nick Perry, Curtis Rueden, Johannes Schindelin,
- * Jean-Yves Tinevez and Michael Zinsmaier.
- * %%
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 
- * 1. Redistributions of source code must retain the above copyright notice,
- *    this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- *    this list of conditions and the following disclaimer in the documentation
- *    and/or other materials provided with the distribution.
- * 
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- * #L%
- */
-
 package net.imglib2;
 
-import java.util.ArrayDeque;
-import java.util.ArrayList;
-import java.util.Comparator;
 import java.util.List;
-import java.util.ListIterator;
-
-import net.imglib2.util.KthElement;
-
-/**
- * KDTree to access values at RealLocalizable positions.
- *
- * @param <T>
- *            type of values stored in the tree.
- *
- * @author Tobias Pietzsch
- */
-public class KDTree< T > implements EuclideanSpace, IterableRealInterval< T >
-{
-	/**
-	 * the number of dimensions.
-	 */
-	final protected int n;
-
-	final protected KDTreeNode< T > root;
+import java.util.function.IntFunction;
+import java.util.function.Supplier;
+import net.imglib2.converter.AbstractConvertedIterableRealInterval;
+import net.imglib2.converter.AbstractConvertedRealCursor;
+import net.imglib2.kdtree.KDTreeData;
+import net.imglib2.kdtree.KDTreeImpl;
 
-	/**
-	 * the number of nodes in the tree.
-	 */
-	final protected long size;
+import static net.imglib2.kdtree.KDTreeData.PositionsLayout.FLAT;
 
-	/**
-	 * minimum of each dimension.
-	 */
-	final protected double[] min;
+public class KDTree< T > implements EuclideanSpace, IterableRealInterval< T >
+{
+	private final KDTreeData< T > treeData;
 
-	/**
-	 * maximum of each dimension.
-	 */
-	final protected double[] max;
+	final KDTreeImpl impl;
 
 	/**
-	 * A KDTreeNode that stores its value as a reference.
+	 * Access to underlying data for serialization.
 	 */
-	protected final static class ValueNode< T > extends KDTreeNode< T >
+	public KDTreeData< T > treeData()
 	{
-		protected final T value;
-
-		/**
-		 * @param value
-		 *            reference to the node's value
-		 * @param position
-		 *            coordinates of this node
-		 * @param dimension
-		 *            dimension along which this node divides the space
-		 * @param left
-		 *            left child node
-		 * @param right
-		 *            right child node
-		 */
-		public ValueNode( final T value, final RealLocalizable position, final int dimension, final ValueNode< T > left, final ValueNode< T > right )
-		{
-			super( position, dimension, left, right );
-			this.value = value;
-		}
-
-		protected ValueNode( final ValueNode< T > node )
-		{
-			super( node );
-			this.value = node.value;
-		}
-
-		@Override
-		public T get()
-		{
-			return value;
-		}
-
-		@Override
-		public ValueNode< T > copy()
-		{
-			return new ValueNode< T >( this );
-		}
-
-		@Override
-		public String toString()
-		{
-			return "node " + getSplitDimension() + " ? " + getSplitCoordinate() + " | " + value;
-		}
+		return treeData;
 	}
 
 	/**
-	 * A KDTreeNode that stores its value as a Sampler.
+	 * Access to pure coordinate kD Tree implementation.
 	 */
-	protected static final class SamplerNode< T > extends KDTreeNode< T >
+	public KDTreeImpl impl()
 	{
-		protected final Sampler< T > sampler;
-
-		/**
-		 * @param sampler
-		 *            a sampler providing the node's value
-		 * @param position
-		 *            coordinates of this node
-		 * @param dimension
-		 *            dimension along which this node divides the space
-		 * @param left
-		 *            left child node
-		 * @param right
-		 *            right child node
-		 */
-		public SamplerNode( final Sampler< T > sampler, final RealLocalizable position, final int dimension, final SamplerNode< T > left, final SamplerNode< T > right )
-		{
-			super( position, dimension, left, right );
-			this.sampler = sampler;
-		}
-
-		protected SamplerNode( final SamplerNode< T > node )
-		{
-			super( node );
-			this.sampler = node.sampler.copy();
-		}
-
-		@Override
-		public T get()
-		{
-			return sampler.get();
-		}
-
-		@Override
-		public SamplerNode< T > copy()
-		{
-			return new SamplerNode< T >( this );
-		}
-
-		@Override
-		public String toString()
-		{
-			return "node " + getSplitDimension() + " ? " + getSplitCoordinate() + " | " + sampler.get();
-		}
+		return impl;
 	}
 
 	/**
@@ -183,58 +41,22 @@ public String toString()
 	 * </p>
 	 *
 	 * @param values
-	 *            a list of values
+	 * 		a list of values
 	 * @param positions
-	 *            a list of positions corresponding to the values
+	 * 		a list of positions corresponding to the values
 	 */
 	public < L extends RealLocalizable > KDTree( final List< T > values, final List< L > positions )
 	{
-		assert values.size() == positions.size();
-
-		this.n = positions.get( 0 ).numDimensions();
-		this.size = positions.size();
-
-		// test that dimensionality is preserved
-		assert ( verifyDimensions( positions, n ) );
-
-		this.min = new double[ n ];
-		this.max = new double[ n ];
-		for ( int d = 0; d < n; ++d )
-		{
-			min[ d ] = Double.MAX_VALUE;
-			max[ d ] = -Double.MAX_VALUE;
-		}
-		for ( final L position : positions )
-		{
-			for ( int d = 0; d < n; ++d )
-			{
-				final double x = position.getDoublePosition( d );
-
-				if ( x < min[ d ] )
-					min[ d ] = x;
-				if ( x > max[ d ] )
-					max[ d ] = x;
-			}
-		}
+		this( verifySize(values, positions), values, positions );
+	}
 
-		if ( values == positions )
-		{
-			if ( positions instanceof java.util.RandomAccess )
-				root = makeNode( positions, 0, positions.size() - 1, 0 );
-			else
-				root = makeNode( positions.listIterator(), positions.listIterator( positions.size() ), 0 );
-		}
-		else
-		{
-			final int[] permutation = new int[ positions.size() ];
-			for ( int k = 0; k < permutation.length; ++k )
-				permutation[ k ] = k;
-
-			if ( positions instanceof java.util.RandomAccess )
-				root = makeNode( positions, 0, positions.size() - 1, 0, values, permutation );
-			else
-				root = makeNode( positions.listIterator(), positions.listIterator( positions.size() ), 0, values, permutation );
-		}
+	private static int verifySize( final List< ? > values, final List< ? > positions )
+	{
+		if ( values.size() != positions.size() )
+			throw new IllegalArgumentException( "The list of values and the list of positions provided to KDTree should have the same size." );
+		if ( positions.isEmpty() )
+			throw new IllegalArgumentException( "List of positions is empty. At least one point is requires to construct a KDTree." );
+		return values.size();
 	}
 
 	/**
@@ -246,466 +68,209 @@ public < L extends RealLocalizable > KDTree( final List< T > values, final List<
 	 */
 	public KDTree( final IterableRealInterval< T > interval )
 	{
-		this.n = interval.numDimensions();
-		this.size = interval.size();
-		this.min = new double[ n ];
-		interval.realMin( this.min );
-		this.max = new double[ n ];
-		interval.realMax( this.max );
-		final ArrayList< RealCursor< T > > values = new ArrayList< RealCursor< T > >( ( int ) interval.size() );
-		final RealCursor< T > cursor = interval.localizingCursor();
-		while ( cursor.hasNext() )
-		{
-			cursor.next();
-			values.add( cursor.copy() );
-		}
-		root = makeSamplerNode( values, 0, values.size() - 1, 0 );
+		this( verifySize( interval ), interval, positionsIterable( interval ) );
 	}
 
-	/**
-	 * Check whether all positions in the positions list have dimension n.
-	 *
-	 * @return true, if all positions have dimension n.
-	 */
-	protected static < L extends RealLocalizable > boolean verifyDimensions( final List< L > positions, final int n )
-	{
-		for ( final L position : positions )
-			if ( position.numDimensions() != n )
-				return false;
-		return true;
-	}
+	private static final int MAX_ARRAY_SIZE = Integer.MAX_VALUE - 8;
 
-	/**
-	 * Compare RealLocalizables by comparing their coordinates in dimension d.
-	 */
-	public static final class DimComparator< L extends RealLocalizable > implements Comparator< L >
+	private static int verifySize( final IterableRealInterval< ? > interval )
 	{
-		final int d;
-
-		public DimComparator( final int d )
-		{
-			this.d = d;
-		}
-
-		@Override
-		public int compare( final L o1, final L o2 )
-		{
-			final float diff = o1.getFloatPosition( d ) - o2.getFloatPosition( d );
-			return ( diff < 0 ) ? -1 : ( diff > 0 ? 1 : 0 );
-		}
+		final long size = interval.size();
+		if ( size > MAX_ARRAY_SIZE )
+			throw new IllegalArgumentException( "Interval contains too many points to store in KDTree" );
+		else if ( size <= 0 )
+			throw new IllegalArgumentException( "Interval is empty. At least one point is requires to construct a KDTree." );
+		return ( int ) size;
 	}
 
-	/**
-	 * Construct the tree by recursively adding nodes. The sublist of positions
-	 * between indices i and j (inclusive) is split at the median element with
-	 * respect to coordinates in the given dimension d. The median becomes the
-	 * new node which is returned. The left and right partitions of the sublist
-	 * are processed recursively and form the left and right subtrees of the
-	 * node.
-	 *
-	 * @param positions
-	 *            list of positions
-	 * @param i
-	 *            start index of sublist to process
-	 * @param j
-	 *            end index of sublist to process
-	 * @param d
-	 *            dimension along which to split the sublist
-	 * @param values
-	 *            list of values corresponding to permuted positions
-	 * @param permutation
-	 *            the index of the values element at index k is permutation[k]
-	 * @return a new node containing the subtree of the given sublist of
-	 *         positions.
-	 */
-	protected < L extends RealLocalizable > ValueNode< T > makeNode( final List< L > positions, final int i, final int j, final int d, final List< T > values, final int[] permutation )
+	private static < A > Iterable< RealLocalizable > positionsIterable( IterableRealInterval< A > sourceInterval )
 	{
-		if ( j > i )
-		{
-			final int k = i + ( j - i ) / 2;
-			KthElement.kthElement( i, j, k, positions, permutation, new DimComparator< L >( d ) );
-
-			final int dChild = ( d + 1 == n ) ? 0 : d + 1;
-			return new ValueNode< T >( values.get( permutation[ k ] ), positions.get( k ), d, makeNode( positions, i, k - 1, dChild, values, permutation ), makeNode( positions, k + 1, j, dChild, values, permutation ) );
-		}
-		else if ( j == i )
-		{
-			return new ValueNode< T >( values.get( permutation[ i ] ), positions.get( i ), d, null, null );
-		}
-		else
+		return new AbstractConvertedIterableRealInterval< A, RealLocalizable >( sourceInterval )
 		{
-			return null;
-		}
-	}
 
-	/**
-	 * Construct the tree by recursively adding nodes. The sublist of positions
-	 * between iterators first and last is split at the median element with
-	 * respect to coordinates in the given dimension d. The median becomes the
-	 * new node which is returned. The left and right partitions of the sublist
-	 * are processed recursively and form the left and right subtrees of the
-	 * node.
-	 *
-	 * @param first
-	 *            first element of the sublist of positions
-	 * @param last
-	 *            last element of the sublist of positions
-	 * @param d
-	 *            dimension along which to split the sublist
-	 * @param values
-	 *            list of values corresponding to permuted positions
-	 * @param permutation
-	 *            the index of the values element at index k is permutation[k]
-	 * @return a new node containing the subtree of the given sublist of
-	 *         positions.
-	 */
-	protected < L extends RealLocalizable > ValueNode< T > makeNode( final ListIterator< L > first, final ListIterator< L > last, final int d, final List< T > values, final int[] permutation )
-	{
-		final int i = first.nextIndex();
-		final int j = last.previousIndex();
-		if ( j > i )
-		{
-			final int k = i + ( j - i ) / 2;
-			KthElement.kthElement( first, last, k, permutation, new DimComparator< L >( d ) );
-			first.previous();
-			final L current = first.next();
-
-			final int dChild = ( d + 1 == n ) ? 0 : d + 1;
-
-			// Node< T > right = makeNode( elements, k + 1, j, dChild );
-			for ( int c = j - last.previousIndex(); c > 0; --c )
-				last.next();
-			final ValueNode< T > right = makeNode( first, last, dChild, values, permutation );
-
-			// Node< T > left = makeNode( elements, i, k - 1, dChild );
-			for ( int c = first.nextIndex() - i; c > 0; --c )
-				first.previous();
-			for ( int c = last.nextIndex() - k; c > 0; --c )
-				last.previous();
-			final ValueNode< T > left = makeNode( first, last, dChild, values, permutation );
-
-			return new ValueNode< T >( values.get( permutation[ k ] ), current, d, left, right );
-		}
-		else if ( j == i )
-		{
-			final L current = first.next();
-			return new ValueNode< T >( values.get( permutation[ i ] ), current, d, null, null );
-		}
-		else
-		{
-			return null;
-		}
-	}
+			class Cursor extends AbstractConvertedRealCursor< A, RealLocalizable >
+			{
+				Cursor( final RealCursor< A > source )
+				{
+					super( source );
+				}
+
+				@Override
+				public RealLocalizable get()
+				{
+					return source;
+				}
+
+				@Override
+				public Cursor copy()
+				{
+					return new Cursor( source.copy() );
+				}
+			}
 
-	/**
-	 * See {@link #makeNode(List, int, int, int, List, int[])}. Here, no values are
-	 * attached to the nodes (or rather the positions are the values).
-	 *
-	 * @param elements
-	 *            list of elements (positions and values at the same time)
-	 * @param i
-	 *            start index of sublist to process
-	 * @param j
-	 *            end index of sublist to process
-	 * @param d
-	 *            dimension along which to split the sublist
-	 */
-	@SuppressWarnings( "unchecked" )
-	protected < L extends RealLocalizable > ValueNode< T > makeNode( final List< L > elements, final int i, final int j, final int d )
-	{
-		if ( j > i )
-		{
-			final int k = i + ( j - i ) / 2;
-			KthElement.kthElement( i, j, k, elements, new DimComparator< L >( d ) );
+			@Override
+			public AbstractConvertedRealCursor< A, RealLocalizable > cursor()
+			{
+				return new Cursor( sourceInterval.cursor() );
+			}
 
-			final int dChild = ( d + 1 == n ) ? 0 : d + 1;
-			return new ValueNode< T >( ( T ) elements.get( k ), elements.get( k ), d, makeNode( elements, i, k - 1, dChild ), makeNode( elements, k + 1, j, dChild ) );
-		}
-		else if ( j == i )
-		{
-			return new ValueNode< T >( ( T ) elements.get( i ), elements.get( i ), d, null, null );
-		}
-		else
-		{
-			return null;
-		}
+			@Override
+			public AbstractConvertedRealCursor< A, RealLocalizable > localizingCursor()
+			{
+				return new Cursor( sourceInterval.localizingCursor() );
+			}
+		};
 	}
 
-	/**
-	 * See {@link #makeNode(ListIterator, ListIterator, int, List, int[])}. Here, no
-	 * values are attached to the nodes (or rather the positions are the
-	 * values).
-	 *
-	 * @param first
-	 *            first element of the sublist to process
-	 * @param last
-	 *            last element of the sublist to process
-	 * @param d
-	 *            dimension along which to split the sublist
-	 */
-	@SuppressWarnings( "unchecked" )
-	protected < L extends RealLocalizable > ValueNode< T > makeNode( final ListIterator< L > first, final ListIterator< L > last, final int d )
+	public < L extends RealLocalizable > KDTree( final int numPoints, final Iterable< T > values, final Iterable< L > positions )
 	{
-		final int i = first.nextIndex();
-		final int j = last.previousIndex();
-		if ( j > i )
-		{
-			final int k = i + ( j - i ) / 2;
-			KthElement.kthElement( first, last, k, new DimComparator< L >( d ) );
-			first.previous();
-			final L current = first.next();
-
-			final int dChild = ( d + 1 == n ) ? 0 : d + 1;
-
-			// Node< T > right = makeNode( elements, k + 1, j, dChild );
-			for ( int c = j - last.previousIndex(); c > 0; --c )
-				last.next();
-			final ValueNode< T > right = makeNode( first, last, dChild );
-
-			// Node< T > left = makeNode( elements, i, k - 1, dChild );
-			for ( int c = first.nextIndex() - i; c > 0; --c )
-				first.previous();
-			for ( int c = last.nextIndex() - k; c > 0; --c )
-				last.previous();
-			final ValueNode< T > left = makeNode( first, last, dChild );
-
-			return new ValueNode< T >( ( T ) current, current, d, left, right );
-		}
-		else if ( j == i )
-		{
-			final L current = first.next();
-			return new ValueNode< T >( ( T ) current, current, d, null, null );
-		}
-		else
-		{
-			return null;
-		}
+		// TODO make storeValuesAsNativeImg a parameter
+		this( KDTreeData.create( numPoints, values, positions, true ) );
 	}
 
-	/**
-	 * Construct the tree by recursively adding nodes. The sublist of elements
-	 * between indices i and j (inclusive) is split at the median element with
-	 * respect to coordinates in the given dimension d. The median becomes the
-	 * new node which is returned. The left and right partitions of the sublist
-	 * are processed recursively and form the left and right subtrees of the
-	 * node. (The elements of the list are RealCursors which provide coordinates
-	 * and values.)
-	 *
-	 * @param elements
-	 *            list of elements (positions and values at the same time)
-	 * @param i
-	 *            start index of sublist to process
-	 * @param j
-	 *            end index of sublist to process
-	 * @param d
-	 *            dimension along which to split the sublist
-	 * @return a new node containing the subtree of the given sublist of
-	 *         elements
-	 */
-	protected SamplerNode< T > makeSamplerNode( final List< RealCursor< T > > elements, final int i, final int j, final int d )
+	// construct with pre-built data, e.g., from deserialization
+	public KDTree( final KDTreeData< T > data )
 	{
-		if ( j > i )
-		{
-			final int k = i + ( j - i ) / 2;
-			KthElement.kthElement( i, j, k, elements, new DimComparator< RealCursor< T > >( d ) );
-
-			final int dChild = ( d + 1 == n ) ? 0 : d + 1;
-			return new SamplerNode< T >( elements.get( k ), elements.get( k ), d, makeSamplerNode( elements, i, k - 1, dChild ), makeSamplerNode( elements, k + 1, j, dChild ) );
-		}
-		else if ( j == i )
-		{
-			return new SamplerNode< T >( elements.get( i ), elements.get( i ), d, null, null );
-		}
-		else
-		{
-			return null;
-		}
+		treeData = data;
+		impl = ( data.layout() == FLAT )
+				? KDTreeImpl.create( data.flatPositions(), data.numDimensions() )
+				: KDTreeImpl.create( data.positions() );
 	}
 
 	/**
 	 * Get the root node.
 	 *
 	 * @return the root node.
+	 *
+	 * @deprecated
+	 * {@link KDTreeNode} is now a re-usable proxy (like {@code NativeType}).
+	 * To work with existing code, {@link KDTreeNode#left()}, {@link
+	 * KDTreeNode#right()}, {@link KDTree#getRoot()} etc create new objects in each
+	 * call, instead of re-using existing proxies.
+	 * Code using that should be rewritten to reuse proxies, if possible.
 	 */
+	@Deprecated
 	public KDTreeNode< T > getRoot()
 	{
-		return root;
+		return new KDTreeNode<>( this ).setNodeIndex( impl.root() );
 	}
 
 	@Override
 	public int numDimensions()
 	{
-		return n;
-	}
-
-	public String toString( final KDTreeNode< T > left, final String indent )
-	{
-		if ( left == null )
-			return "";
-
-		return indent + "- " + left.toString() + "\n" + toString( left.left, indent + "  " ) + toString( left.right, indent + "  " );
-	}
-
-	@Override
-	public String toString()
-	{
-		return toString( root, "" );
+		return impl.numDimensions();
 	}
 
 	@Override
 	public double realMin( final int d )
 	{
-		return min[ d ];
-	}
-
-	@Override
-	public void realMin( final double[] m )
-	{
-		for ( int d = 0; d < n; ++d )
-			m[ d ] = min[ d ];
-	}
-
-	@Override
-	public void realMin( final RealPositionable m )
-	{
-		m.setPosition( min );
+		return treeData.boundingBox().realMin( d );
 	}
 
 	@Override
 	public double realMax( final int d )
 	{
-		return max[ d ];
-	}
-
-	@Override
-	public void realMax( final double[] m )
-	{
-		for ( int d = 0; d < n; ++d )
-			m[ d ] = max[ d ];
-	}
-
-	@Override
-	public void realMax( final RealPositionable m )
-	{
-		m.setPosition( max );
+		return treeData.boundingBox().realMax( d );
 	}
 
 	@Override
-	public long size()
+	public KDTreeCursor cursor()
 	{
-		return size;
+		return new KDTreeCursor();
 	}
 
-	@Override
-	public Object iterationOrder()
+	public final class KDTreeCursor extends KDTreeNode< T > implements RealCursor< T >
 	{
-		return this; // iteration order is only compatible with ourselves
-	}
-
-	public final class KDTreeCursor implements RealCursor< T >
-	{
-		private final KDTree< T > tree;
-
-		private final ArrayDeque< KDTreeNode< T > > nodes;
-
-		private KDTreeNode< T > currentNode;
-
-		public KDTreeCursor( final KDTree< T > kdtree )
+		KDTreeCursor()
 		{
-			this.tree = kdtree;
-			final int capacity = 2 + ( int ) ( Math.log( kdtree.size() ) / Math.log( 2 ) );
-			this.nodes = new ArrayDeque< KDTreeNode< T > >( capacity );
+			super( KDTree.this );
 			reset();
 		}
 
-		public KDTreeCursor( final KDTreeCursor c )
-		{
-			this.tree = c.tree;
-			this.nodes = c.nodes.clone();
-			this.currentNode = c.currentNode;
-		}
-
 		@Override
-		public void localize( final float[] position )
-		{
-			currentNode.localize( position );
-		}
-
-		@Override
-		public void localize( final double[] position )
-		{
-			currentNode.localize( position );
-		}
-
-		@Override
-		public float getFloatPosition( final int d )
-		{
-			return currentNode.getFloatPosition( d );
-		}
-
-		@Override
-		public double getDoublePosition( final int d )
+		public void fwd()
 		{
-			return currentNode.getDoublePosition( d );
+			setNodeIndex( nodeIndex() + 1 );
 		}
 
 		@Override
-		public int numDimensions()
+		public void reset()
 		{
-			return n;
+			setNodeIndex( -1 );
 		}
 
 		@Override
-		public T get()
+		public boolean hasNext()
 		{
-			return currentNode.get();
+			return nodeIndex() < impl.size() - 1;
 		}
 
 		@Override
 		public KDTreeCursor copy()
 		{
-			return new KDTreeCursor( this );
+			final KDTreeCursor copy = new KDTreeCursor();
+			copy.setNodeIndex( nodeIndex() );
+			return copy;
 		}
+	}
 
-		@Override
-		public void fwd()
-		{
-			if ( nodes.isEmpty() )
-				currentNode = null;
-			else
-			{
-				currentNode = nodes.pop();
-				if ( currentNode.left != null )
-					nodes.push( currentNode.left );
-				if ( currentNode.right != null )
-					nodes.push( currentNode.right );
-			}
-		}
+	@Override
+	public KDTreeCursor localizingCursor()
+	{
+		return cursor();
+	}
 
-		@Override
-		public void reset()
-		{
-			nodes.clear();
-			nodes.push( tree.getRoot() );
-			currentNode = null;
-		}
+	@Override
+	public KDTreeCursor iterator()
+	{
+		return cursor();
+	}
 
-		@Override
-		public boolean hasNext()
-		{
-			return !nodes.isEmpty();
-		}
+	@Override
+	public long size()
+	{
+		return impl.size();
 	}
 
 	@Override
-	public KDTreeCursor cursor()
+	public Object iterationOrder()
 	{
-		return new KDTreeCursor( this );
+		return this; // iteration order is only compatible with ourselves
 	}
 
 	@Override
-	public KDTreeCursor localizingCursor()
+	public String toString()
+	{
+		return toString( impl.root(), "", createNode() );
+	}
+
+	private String toString( final int node, final String indent, final KDTreeNode< T > ref )
+	{
+		if ( node < 0 )
+			return "";
+		return indent + "- " + ref.setNodeIndex( node ).toString() + "\n"
+				+ toString( impl.left( node ), indent + "  ", ref )
+				+ toString( impl.right( node ), indent + "  ", ref );
+	}
+
+	/**
+	 * Create a re-usable {@link KDTreeNode} proxy linked to this tree.
+	 * {@link KDTreeNode#setNodeIndex(int)} can be used to point the proxy to a
+	 * particular node in the tree.
+	 */
+	public KDTreeNode< T > createNode() {
+		return new KDTreeNode<>( this );
+	}
+
+	KDTreeNode< T > left( final KDTreeNode< T > parent )
+	{
+		final int c = impl.left( parent.nodeIndex() );
+		return c < 0 ? null : createNode().setNodeIndex( c );
+	}
+
+	KDTreeNode< T > right( final KDTreeNode< T > parent )
 	{
-		return new KDTreeCursor( this );
+		final int c = impl.right( parent.nodeIndex() );
+		return c < 0 ? null : createNode().setNodeIndex( c );
 	}
 }
diff --git a/src/main/java/net/imglib2/KDTreeNode.java b/src/main/java/net/imglib2/KDTreeNode.java
index 4692a804a..6a0d4d421 100644
--- a/src/main/java/net/imglib2/KDTreeNode.java
+++ b/src/main/java/net/imglib2/KDTreeNode.java
@@ -1,107 +1,80 @@
-/*
- * #%L
- * ImgLib2: a general-purpose, multidimensional image processing library.
- * %%
- * Copyright (C) 2009 - 2023 Tobias Pietzsch, Stephan Preibisch, Stephan Saalfeld,
- * John Bogovic, Albert Cardona, Barry DeZonia, Christian Dietz, Jan Funke,
- * Aivar Grislis, Jonathan Hale, Grant Harris, Stefan Helfrich, Mark Hiner,
- * Martin Horn, Steffen Jaensch, Lee Kamentsky, Larry Lindsey, Melissa Linkert,
- * Mark Longair, Brian Northan, Nick Perry, Curtis Rueden, Johannes Schindelin,
- * Jean-Yves Tinevez and Michael Zinsmaier.
- * %%
- * Redistribution and use in source and binary forms, with or without
- * modification, are permitted provided that the following conditions are met:
- * 
- * 1. Redistributions of source code must retain the above copyright notice,
- *    this list of conditions and the following disclaimer.
- * 2. Redistributions in binary form must reproduce the above copyright notice,
- *    this list of conditions and the following disclaimer in the documentation
- *    and/or other materials provided with the distribution.
- * 
- * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
- * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
- * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
- * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
- * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
- * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
- * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
- * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
- * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
- * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
- * POSSIBILITY OF SUCH DAMAGE.
- * #L%
- */
-
 package net.imglib2;
 
+import java.util.function.IntFunction;
 
 /**
- * Abstract base class for nodes in a KDTree. A KDTreeNode has coordinates and a
- * value. It provides the coordinates via the {@link RealLocalizable} interface.
- * It provides the value via {@link Sampler#get()}.
+ * Proxy for a node in a KDTree. A KDTreeNode has coordinates and a value. It
+ * provides the coordinates via the {@link RealLocalizable} interface. It
+ * provides the value via {@link Sampler#get()}.
  *
  * @param <T>
  *            value type.
  *
- *
  * @author Tobias Pietzsch
  */
-public abstract class KDTreeNode< T > implements RealLocalizable, Sampler< T >
+public class KDTreeNode< T > implements RealLocalizable, Sampler< T >
 {
-	/**
-	 * number of dimensions of the space (that is, k).
-	 */
-	protected final int n;
+	private final KDTree< T > tree;
+
+	private int nodeIndex;
+
+	private IntFunction< T > values;
+
+	KDTreeNode( final KDTree< T > tree )
+	{
+		this.tree = tree;
+	}
 
 	/**
-	 * coordinates of the node.
+	 * Make this proxy refer to the given {@code nodeIndex} in the associated tree.
+	 *
+	 * @return {@code this}
 	 */
-	protected final double[] pos;
+	public KDTreeNode< T > setNodeIndex( final int nodeIndex )
+	{
+		this.nodeIndex = nodeIndex;
+		return this;
+	}
 
 	/**
-	 * dimension along which this node divides the space.
+	 * Get the {@code nodeIndex} which this proxy currently refers to.
 	 */
-	protected final int splitDimension;
+	public int nodeIndex() {
+		return nodeIndex;
+	}
 
 	/**
 	 * Left child of this node. All nodes x in the left subtree have
 	 * {@code x.pos[splitDimension] <= this.pos[splitDimension]}.
+	 *
+	 * @deprecated
+	 * {@link KDTreeNode} is now a re-usable proxy (like {@code NativeType}).
+	 * To work with existing code, {@link KDTreeNode#left()}, {@link
+	 * KDTreeNode#right()}, {@link KDTree#getRoot()} etc create new objects in each
+	 * call, instead of re-using existing proxies.
+	 * Code using that should be rewritten to reuse proxies, if possible.
 	 */
-	public final KDTreeNode< T > left;
+	@Deprecated
+	public KDTreeNode< T > left()
+	{
+		return tree.left( this );
+	}
 
 	/**
 	 * Right child of this node. All nodes x in the right subtree have
 	 * {@code x.pos[splitDimension] >= this.pos[splitDimension]}.
+	 *
+	 * @deprecated
+	 * {@link KDTreeNode} is now a re-usable proxy (like {@code NativeType}).
+	 * To work with existing code, {@link KDTreeNode#left()}, {@link
+	 * KDTreeNode#right()}, {@link KDTree#getRoot()} etc create new objects in each
+	 * call, instead of re-using existing proxies.
+	 * Code using that should be rewritten to reuse proxies, if possible.
 	 */
-	public final KDTreeNode< T > right;
-
-	/**
-	 * @param position
-	 *            coordinates of this node
-	 * @param dimension
-	 *            dimension along which this node divides the space
-	 * @param left
-	 *            left child node
-	 * @param right
-	 *            right child node
-	 */
-	public KDTreeNode( final RealLocalizable position, final int dimension, final KDTreeNode< T > left, final KDTreeNode< T > right )
+	@Deprecated
+	public KDTreeNode< T > right()
 	{
-		this.n = position.numDimensions();
-		this.pos = new double[ n ];
-		position.localize( pos );
-		this.splitDimension = dimension;
-		this.left = left;
-		this.right = right;
-	}
-
-	protected KDTreeNode( final KDTreeNode< T > node )
-	{
-		this.n = node.n;
-		this.pos = node.pos.clone();
-		this.splitDimension = node.splitDimension;
-		this.left = node.left;
-		this.right = node.right;
+		return tree.right( this );
 	}
 
 	/**
@@ -111,79 +84,62 @@ protected KDTreeNode( final KDTreeNode< T > node )
 	 */
 	public final int getSplitDimension()
 	{
-		return splitDimension;
+		return tree.impl.splitDimension( nodeIndex );
 	}
 
 	/**
-	 * Get the position along {@link KDTreeNode#getSplitDimension()} where this
+	 * Get the position along {@link net.imglib2.KDTreeNode#getSplitDimension()} where this
 	 * node divides the space.
 	 *
 	 * @return splitting position.
 	 */
 	public final double getSplitCoordinate()
 	{
-		return pos[ splitDimension ];
+		return getDoublePosition( getSplitDimension() );
 	}
 
 	@Override
-	public final int numDimensions()
+	public int numDimensions()
 	{
-		return n;
+		return tree.numDimensions();
 	}
 
 	@Override
-	public final void localize( final float[] position )
+	public double getDoublePosition( final int d )
 	{
-		for ( int d = 0; d < n; ++d )
-			position[ d ] = ( float ) pos[ d ];
+		return tree.impl.getDoublePosition( nodeIndex, d );
 	}
 
 	@Override
-	public final void localize( final double[] position )
+	public T get()
 	{
-		for ( int d = 0; d < n; ++d )
-			position[ d ] = pos[ d ];
+		if ( values == null )
+			values = tree.treeData().valuesSupplier().get();
+		return values.apply( nodeIndex );
 	}
 
 	@Override
-	public final float getFloatPosition( final int d )
+	public KDTreeNode< T > copy()
 	{
-		return ( float ) pos[ d ];
+		final KDTreeNode< T > copy = new KDTreeNode<>( tree );
+		copy.setNodeIndex( nodeIndex );
+		return copy;
 	}
 
-	@Override
-	public final double getDoublePosition( final int d )
-	{
-		return pos[ d ];
-	}
-
-	@Override
-	public abstract KDTreeNode< T > copy();
-
 	/**
 	 * Compute the squared distance from p to this node.
 	 */
 	public final float squDistanceTo( final float[] p )
 	{
-		float sum = 0;
-		for ( int d = 0; d < n; ++d )
-		{
-			sum += ( pos[ d ] - p[ d ] ) * ( pos[ d ] - p[ d ] );
-		}
-		return sum;
+		return tree.impl.squDistance( nodeIndex, p );
 	}
 
 	/**
 	 * Compute the squared distance from p to this node.
 	 */
-	public final double squDistanceTo( final double[] p )
+	public double squDistanceTo( final double[] p )
 	{
-		double sum = 0;
-		for ( int d = 0; d < n; ++d )
-		{
-			sum += ( pos[ d ] - p[ d ] ) * ( pos[ d ] - p[ d ] );
-		}
-		return sum;
+		return tree.impl.squDistance( nodeIndex, p );
 	}
 
 	/**
@@ -191,11 +147,6 @@ public final double squDistanceTo( final double[] p )
 	 */
 	public final double squDistanceTo( final RealLocalizable p )
 	{
-		double sum = 0;
-		for ( int d = 0; d < n; ++d )
-		{
-			sum += ( pos[ d ] - p.getDoublePosition( d ) ) * ( pos[ d ] - p.getDoublePosition( d ) );
-		}
-		return sum;
+		return tree.impl.squDistance( nodeIndex, p );
 	}
 }
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
new file mode 100644
index 000000000..761fea25a
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -0,0 +1,323 @@
+package net.imglib2.kdtree;
+
+import java.util.List;
+import java.util.function.IntFunction;
+import java.util.function.Supplier;
+import net.imglib2.FinalRealInterval;
+import net.imglib2.RandomAccess;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.RealInterval;
+import net.imglib2.RealLocalizable;
+import net.imglib2.img.Img;
+import net.imglib2.img.list.ListImg;
+import net.imglib2.type.NativeType;
+import net.imglib2.util.Util;
+
+import static net.imglib2.kdtree.KDTreeData.PositionsLayout.FLAT;
+
+/**
+ * Stores the KDTree data, that is, positions and values.
+ * <p>
+ * Positions are stored in either {@code FLAT} or {@code NESTED} {@link
+ * PositionsLayout layout}. With {@code NESTED} layout, positions are stored as
+ * a nested {@code double[][]} array where {@code positions[d][i]} is dimension
+ * {@code d} of the {@code i}-th point. With {@code FLAT} layout, positions are
+ * stored as a flat {@code double[]} array, where {@code positions[d + i*n]} is
+ * dimension {@code d} of the {@code i}-th point, with {@code n} the number of
+ * dimensions.
+ * <p>
+ * Values (of type {@code T}) are stored as either a 1D {@code
+ * RandomAccessibleInterval<T>}, or a {@code List<T>}. Individual values can be
+ * accessed by {@link #valuesSupplier()}{@code .get().apply(i)}. {@code
+ * valueSupplier().get()} returns a reusable {@code IntFunction<T>}. Here {@code
+ * T} maybe a proxy that is reused in subsequent {@code apply(i)}.
+ * <p>
+ * {@link #values()} returns all values as a 1D {@code
+ * RandomAccessibleInterval<T>}. (If data is stored as {@code List<T>}, it is
+ * wrapped into a {@code ListImg}.)
+ * <p>
+ * {@link #positions()} returns positions in nested {@code double[][]} (which is
+ * created if internal storage is {@code FLAT}). {@link #flatPositions()}
+ * returns flat {@code double[]} if internal storage is {@code FLAT}, otherwise
+ * {@code null}.
+ *
+ * @param <T>
+ * 		the type of values stored in the tree.
+ */
+public class KDTreeData< T >
+{
+	public enum PositionsLayout
+	{
+		FLAT,
+		NESTED
+	}
+
+	private final int numDimensions;
+	private final int numPoints;
+
+	private final PositionsLayout layout;
+	private final double[][] positions;
+	private final double[] flatPositions;
+
+	private final List< T > valuesList;
+	private final RandomAccessibleInterval< T > valuesImg;
+	private final Supplier< IntFunction< T > > valuesSupplier;
+
+	private final T type;
+
+	private volatile RealInterval boundingBox;
+
+	public KDTreeData( double[][] positions, List< T > values )
+	{
+		numPoints = values.size();
+		numDimensions = positions.length;
+
+		layout = PositionsLayout.NESTED;
+		this.positions = positions;
+		flatPositions = null;
+
+		valuesList = values;
+		valuesImg = null;
+		final IntFunction< T > v = values::get;
+		valuesSupplier = () -> v;
+
+		type = KDTreeUtils.getType( values );
+	}
+
+	public KDTreeData( double[][] positions, List< T > values, RealInterval boundingBox )
+	{
+		this( positions, values );
+		this.boundingBox = boundingBox;
+	}
+
+	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
+	{
+		numPoints = ( int ) values.dimension( 0 );
+		numDimensions = positions.length;
+
+		layout = PositionsLayout.NESTED;
+		this.positions = positions;
+		flatPositions = null;
+
+		valuesList = null;
+		valuesImg = values;
+		valuesSupplier = () -> {
+			final RandomAccess<T> ra = valuesImg.randomAccess();
+			return i -> ra.setPositionAndGet(i);
+		};
+
+		type = Util.getTypeFromInterval( values );
+	}
+
+	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
+	{
+		this( positions, values );
+		this.boundingBox = boundingBox;
+	}
+
+	public KDTreeData( double[] positions, List< T > values )
+	{
+		numPoints = values.size();
+		numDimensions = positions.length / numPoints;
+
+		layout = FLAT;
+		this.positions = null;
+		flatPositions = positions;
+
+		valuesList = values;
+		valuesImg = null;
+		final IntFunction< T > v = values::get;
+		valuesSupplier = () -> v;
+
+		type = KDTreeUtils.getType( values );
+	}
+
+	public KDTreeData( double[] positions, List< T > values, RealInterval boundingBox )
+	{
+		this( positions, values );
+		this.boundingBox = boundingBox;
+	}
+
+	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
+	{
+		numPoints = ( int ) values.dimension( 0 );
+		numDimensions = positions.length / numPoints;
+
+		layout = FLAT;
+		this.positions = null;
+		flatPositions = positions;
+
+		valuesList = null;
+		valuesImg = values;
+		valuesSupplier = () -> {
+			final RandomAccess<T> ra = valuesImg.randomAccess();
+			return i -> ra.setPositionAndGet(i);
+		};
+
+		type = Util.getTypeFromInterval( values );
+	}
+
+	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
+	{
+		this( positions, values );
+		this.boundingBox = boundingBox;
+	}
+
+	// TODO could also be Class<T> instead? What is more useful?
+	public T type()
+	{
+		return type;
+	}
+
+	/**
+	 * Get positions of points in the tree as a nested {@code double[][]} array
+	 * where {@code positions[d][i]} is dimension {@code d} of the {@code i}-th
+	 * point.
+	 * <p>
+	 * For serialisation and usage by the tree.
+	 * <p>
+	 * Internal storage may be flattened into single {@code double[]} array. In
+	 * this case, the nested {@code double[][]} array is created here.
+	 */
+	public double[][] positions()
+	{
+		return layout == FLAT ? KDTreeUtils.unflatten( flatPositions, numDimensions ) : positions;
+	}
+
+	/**
+	 * Get the values as a 1D {@code RandomAccessibleInterval}, for
+	 * serialization. (If the underlying storage is a {@code List<T>}, it will
+	 * be wrapped as a {@code ListImg}.)
+	 */
+	public RandomAccessibleInterval< T > values()
+	{
+		return valuesImg != null
+				? valuesImg
+				: ListImg.wrap( valuesList, size() );
+	}
+
+	/**
+	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
+	 * values for a given node indices.. If the returned {@code IntFunction<T>}
+	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
+	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
+	 * {@code IntFunction<T>}.
+	 */
+	public Supplier< IntFunction< T > > valuesSupplier()
+	{
+		return valuesSupplier;
+	}
+
+	/**
+	 * Get positions of points in the tree as a flat {@code double[]} where
+	 * @code positions[d + i*n]} is dimension {@code d} of the {@code i}-th
+	 * point, with {@code n} the number of dimensions.
+	 * <p>
+	 * For serialisation and usage by the tree.
+	 * <p>
+	 * Internal storage may be a {@code NESTED} {@code double[][]} array. In
+	 * this case, {@code flatPositions()} returns {@code null}.
+	 */
+	public double[] flatPositions()
+	{
+		return flatPositions;
+	}
+
+	/**
+	 * Get the internal layout of positions.
+	 * <p>
+	 * Positions are stored in either {@code FLAT} or {@code NESTED} {@link
+	 * PositionsLayout layout}. With {@code NESTED} layout, positions are stored
+	 * as a nested {@code double[][]} array where {@code positions[d][i]} is
+	 * dimension {@code d} of the {@code i}-th point. With {@code FLAT} layout,
+	 * positions are stored as a flat {@code double[]} array, where {@code
+	 * positions[d + i*n]} is dimension {@code d} of the {@code i}-th point,
+	 * with {@code n} the number of dimensions.
+	 */
+	public PositionsLayout layout()
+	{
+		return layout;
+	}
+
+	/**
+	 * @return dimensionality of points in the tree
+	 */
+	public int numDimensions()
+	{
+		return numDimensions;
+	}
+
+	/**
+	 * @return number of points in the tree
+	 */
+	public int size()
+	{
+		return numPoints;
+	}
+
+	public RealInterval boundingBox()
+	{
+		if ( boundingBox == null )
+			boundingBox = createBoundingBox();
+		return boundingBox;
+	}
+
+	private RealInterval createBoundingBox()
+	{
+		final double[] min = new double[ numDimensions ];
+		final double[] max = new double[ numDimensions ];
+		if ( layout == FLAT )
+			KDTreeUtils.computeMinMax( flatPositions, min, max );
+		else
+			KDTreeUtils.computeMinMax( positions, min, max );
+		return FinalRealInterval.wrap( min, max );
+	}
+
+	/**
+	 * Create {@link KDTreeData} from the given {@code values} and {@code positions}).
+	 * (copies {@code positions} and sorts into a KDTree structure).
+	 *
+	 * @param numPoints
+	 * 		number of points (number of elements in {@code values} and {@code positions}).
+	 * @param values
+	 * 		values associated with points
+	 * @param positions
+	 * 		points positions
+	 * @param storeValuesAsNativeImg
+	 * 		If {@code true} and {@code T} is a {@code NativeType},
+	 * 		store values into {@code NativeImg}.
+	 * 		Otherwise, store values as a {@code List<T>}.
+	 */
+	public static < L extends RealLocalizable, T > KDTreeData< T > create(
+			final int numPoints,
+			final Iterable< T > values,
+			final Iterable< L > positions,
+			final boolean storeValuesAsNativeImg )
+	{
+		if ( numPoints <= 0 )
+			throw new IllegalArgumentException( "At least one point is required to construct a KDTree." );
+		final int numDimensions = KDTreeUtils.getNumDimensions( positions );
+		final double[][] points = KDTreeUtils.initPositions( numDimensions, numPoints, positions );
+		final int[] tree = KDTreeUtils.makeTree( points );
+		final int[] invtree = KDTreeUtils.invert( tree );
+
+		final boolean useFlatLayout = ( long ) numDimensions * numPoints <= KDTreeUtils.MAX_ARRAY_SIZE;
+		if ( storeValuesAsNativeImg && KDTreeUtils.getType( values ) instanceof NativeType )
+		{
+			@SuppressWarnings( "unchecked" )
+			final Img< T > treeValues = ( Img< T > ) KDTreeUtils.orderValuesImg( invtree, ( Iterable ) values );
+			if ( useFlatLayout )
+				return new KDTreeData<>(KDTreeUtils.reorderToFlatLayout( points, tree ), treeValues);
+			else
+				return new KDTreeData<>(KDTreeUtils.reorder( points, tree ), treeValues);
+		}
+		else
+		{
+			final List< T > treeValues = KDTreeUtils.orderValuesList( invtree, values );
+			if ( useFlatLayout )
+				return new KDTreeData<>(KDTreeUtils.reorderToFlatLayout( points, tree ), treeValues);
+			else
+				return new KDTreeData<>(KDTreeUtils.reorder( points, tree ), treeValues);
+		}
+	}
+}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeImpl.java b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
new file mode 100644
index 000000000..2f4468e77
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
@@ -0,0 +1,229 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RealLocalizable;
+
+import static net.imglib2.kdtree.KDTreeUtils.leftChildIndex;
+import static net.imglib2.kdtree.KDTreeUtils.parentIndex;
+import static net.imglib2.kdtree.KDTreeUtils.rightChildIndex;
+
+/**
+ * Represents the tree structure, and provides access to node coordinates (but
+ * not values).
+ * <p>
+ * The nodes in the tree are arranged in Eytzinger layout (children of i are at
+ * 2i and 2i+1). Additionally, pivot indices are chosen such that "leaf layers"
+ * are filled from the left.
+ * <p>
+ * For example 10 nodes will always be arranged like this:
+ * <pre>
+ *            0
+ *         /     \
+ *       1         2
+ *     /   \     /   \
+ *    3     4   5     6
+ *   / \   /
+ *  7   8 9
+ * </pre>
+ *
+ * never like this:
+ * <pre>
+ *            0
+ *         /     \
+ *       1         2
+ *     /   \     /   \
+ *    3     4   5     6
+ *   /         /     /
+ *  7         8     9
+ * </pre>
+ *
+ * By choosing pivots in this way, the tree structure is fully
+ * determined. For every node index, the child indices can be calculated
+ * without dependent reads. And iff the calculated child index is less
+ * than the number of nodes, the child exists.
+ */
+public abstract class KDTreeImpl
+{
+	final int numDimensions;
+
+	private final int numPoints;
+
+	static class Nested extends KDTreeImpl
+	{
+		private final double[][] positions;
+
+		Nested( final double[][] positions )
+		{
+			super( positions.length, positions[ 0 ].length );
+			this.positions = positions;
+		}
+
+		@Override
+		public double getDoublePosition( final int i, final int d )
+		{
+			return positions[ d ][ i ];
+		}
+	}
+
+	static class Flat extends KDTreeImpl
+	{
+		private final double[] positions;
+
+		Flat( final double[] positions, final int numDimensions )
+		{
+			super( numDimensions, positions.length / numDimensions );
+			this.positions = positions;
+		}
+
+		@Override
+		public double getDoublePosition( final int i, final int d )
+		{
+			return positions[ numDimensions * i + d ];
+		}
+	}
+
+	KDTreeImpl( final int numDimensions, final int numPoints )
+	{
+		this.numDimensions = numDimensions;
+		this.numPoints = numPoints;
+	}
+
+	/**
+	 * Get the root node of the tree.
+	 *
+	 * @return index of the root node
+	 */
+	public int root()
+	{
+		return 0;
+	}
+
+	/**
+	 * Get the left child of node {@code i}.
+	 *
+	 * @param i
+	 * 		node index
+	 *
+	 * @return index of left child or {@code -1} if no left child exists
+	 */
+	public int left( final int i )
+	{
+		return ifExists( leftChildIndex( i ) );
+	}
+
+	/**
+	 * Get the right child of node {@code i}.
+	 *
+	 * @param i
+	 * 		node index
+	 *
+	 * @return index of right child or {@code -1} if no right child exists
+	 */
+	public int right( final int i )
+	{
+		return ifExists( rightChildIndex( i ) );
+	}
+
+	/**
+	 * Get the parent of node {@code i}.
+	 *
+	 * @param i
+	 * 		node index
+	 *
+	 * @return index of parent
+	 */
+	public int parent( final int i )
+	{
+		return i == root() ? -1 : parentIndex( i );
+	}
+
+	/**
+	 * If a node with index {@code i} exists, returns {@code i}.
+	 * Otherwise, returns {@code -1}.
+	 */
+	private int ifExists( final int i )
+	{
+		return i < numPoints ? i : -1;
+	}
+
+	/**
+	 * Get the dimension along which node {@code i} divides the space.
+	 *
+	 * @param i
+	 * 		node index
+	 *
+	 * @return splitting dimension.
+	 */
+	public int splitDimension( final int i )
+	{
+		return ( 31 - Integer.numberOfLeadingZeros( i + 1 ) ) % numDimensions;
+	}
+
+	public abstract double getDoublePosition( final int i, final int d );
+
+	/**
+	 * Compute the squared distance from node {@code i} to {@code pos}.
+	 */
+	public float squDistance( final int i, final float[] pos )
+	{
+		float sum = 0;
+		for ( int d = 0; d < numDimensions; ++d )
+		{
+			final float diff = pos[ d ] - ( float ) getDoublePosition( i, d );
+			sum += diff * diff;
+		}
+		return sum;
+	}
+
+	/**
+	 * Compute the squared distance from node {@code i} to {@code pos}.
+	 */
+	public double squDistance( final int i, final double[] pos )
+	{
+		double sum = 0;
+		for ( int d = 0; d < numDimensions; ++d )
+		{
+			final double diff = pos[ d ] - getDoublePosition( i, d );
+			sum += diff * diff;
+		}
+		return sum;
+	}
+
+	/**
+	 * Compute the squared distance from node {@code i} to {@code pos}.
+	 */
+	public double squDistance( final int i, final RealLocalizable pos )
+	{
+		double sum = 0;
+		for ( int d = 0; d < numDimensions; ++d )
+		{
+			final double diff = pos.getDoublePosition( d ) - getDoublePosition( i, d );
+			sum += diff * diff;
+		}
+		return sum;
+	}
+
+	public int numDimensions()
+	{
+		return numDimensions;
+	}
+
+	public int size()
+	{
+		return numPoints;
+	}
+
+	public int depth()
+	{
+		return 32 - Integer.numberOfLeadingZeros( numPoints );
+	}
+
+	public static KDTreeImpl create( final double[][] positions )
+	{
+		return new Nested( positions );
+	}
+
+	public static KDTreeImpl create( final double[] positions, final int numDimensions )
+	{
+		return new Flat( positions, numDimensions );
+	}
+}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeUtils.java b/src/main/java/net/imglib2/kdtree/KDTreeUtils.java
new file mode 100644
index 000000000..c972ff12b
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeUtils.java
@@ -0,0 +1,573 @@
+package net.imglib2.kdtree;
+
+import java.util.Arrays;
+import java.util.Iterator;
+import java.util.List;
+import net.imglib2.RandomAccess;
+import net.imglib2.RealLocalizable;
+import net.imglib2.img.Img;
+import net.imglib2.img.array.ArrayImgFactory;
+import net.imglib2.type.NativeType;
+
+final class KDTreeUtils
+{
+	static final int MAX_ARRAY_SIZE = Integer.MAX_VALUE - 8;
+
+	/**
+	 * If the tree is flattened into an array the left child of node at
+	 * index {@code i} has index {@code 2 * i + 1}.
+	 */
+	static int leftChildIndex( final int i )
+	{
+		return 2 * i + 1;
+	}
+
+	/**
+	 * If the tree is flattened into an array the right child of node at
+	 * index {@code i} has index {@code 2 * i + 2}.
+	 */
+	static int rightChildIndex( final int i )
+	{
+		return 2 * i + 2;
+	}
+
+	/**
+	 * If the tree is flattened into an array the parent of node at
+	 * index {@code i} has index {@code (i - 1) / 2} (except for the
+	 * root node {@code i==0}).
+	 */
+	static int parentIndex( final int i )
+	{
+		return ( i - 1 ) / 2;
+	}
+
+	/**
+	 * Copy the coordinates of the given {@code points} to a new {@code double[][] positions} array.
+	 * The coordinate in dimension {@code d} of the {@code i}th point is stored at {@code positions[d][i]}.
+	 * That is, {@code positions[0]} has all X coordinates, {@code positions[0]} has all Y coordinates, and so on.
+	 */
+	static double[][] initPositions(
+			final int numDimensions,
+			final int numPoints,
+			final Iterable< ? extends RealLocalizable > points )
+	{
+		final double[][] positions = new double[ numDimensions ][ numPoints ];
+		final Iterator< ? extends RealLocalizable > ipos =  points.iterator();
+		for ( int i = 0; i < numPoints; ++i )
+		{
+			if ( !ipos.hasNext() )
+				throw new IllegalArgumentException( "positions Iterable is empty" );
+			final RealLocalizable pos = ipos.next();
+			for ( int d = 0; d < numDimensions; d++ )
+				positions[ d ][ i ] = pos.getDoublePosition( d );
+		}
+		return positions;
+	}
+
+	/**
+	 * Sort the given points into a k-d tree.
+	 * <p>
+	 * The tree is given as a flat (heap-like) array of point indices, {@code int[] tree}.
+	 * The index of the point chosen as the root node is {@code tree[0]}.
+	 * The coordinates of the root node are at {@code positions[d][tree[0]]}.
+	 * <p>
+	 * The indices of the children (less-or-equal and greater-or-equal, respectively) of root are at {@code tree[1]} and {@code tree[2]}, and so on.
+	 *
+	 * @param positions
+	 * 		The coordinates for the {@code i}th point are stored at {@code positions[d][i]} where {@code d} is the dimension.
+	 * 	    See {@link #initPositions(int, int, Iterable)}.
+	 *
+	 * @return flattened tree of point indices
+	 */
+	static int[] makeTree( double[][] positions )
+	{
+		return new MakeTree( positions ).tree;
+	}
+
+	/**
+	 * Re-order the node {@code positions} to form a tree corresponding to the index array {@code tree'={0,1,2,...}}.
+	 *
+	 * @param positions
+	 * @param tree
+	 *
+	 * @return
+	 */
+	static double[][] reorder( double[][] positions, int[] tree )
+	{
+		final int numDimensions = positions.length;
+		final int numPoints = positions[ 0 ].length;
+		assert tree.length == numPoints;
+		final double[][] reordered = new double[ numDimensions ][];
+		Arrays.setAll( reordered, d -> reorder( positions[ d ], tree ) );
+		return reordered;
+	}
+
+	/**
+	 * Create a new {@code double[]} array that contains the elements of {@code
+	 * values}, ordered such that {@code values[order[i]]} is at index {@code i}.
+	 */
+	static double[] reorder( final double[] values, final int[] order )
+	{
+		final int size = order.length;
+		final double[] reordered = new double[ size ];
+		Arrays.setAll( reordered, i -> values[ order[ i ] ] );
+		return reordered;
+	}
+
+	/**
+	 * Create a new {@code int[]} array that contains the elements of {@code
+	 * values}, ordered such that {@code values[order[i]]} is at index {@code i}.
+	 */
+	static int[] reorder( final int[] values, final int[] order )
+	{
+		final int size = order.length;
+		final int[] reordered = new int[ size ];
+		Arrays.setAll( reordered, i -> values[ order[ i ] ] );
+		return reordered;
+	}
+
+	/**
+	 * Re-order the node {@code positions} to form a tree corresponding to the index array {@code tree={0,1,2,...}}.
+	 * Then flatten the result into a 1-D array, interleaving coordinates in all dimensions.
+	 *
+	 * @param positions
+	 * @param tree
+	 *
+	 * @return
+	 */
+	static double[] reorderToFlatLayout( final double[][] positions, final int[] tree )
+	{
+		final int numDimensions = positions.length;
+		final int numPoints = positions[ 0 ].length;
+		assert tree.length == numPoints;
+		if ( ( long ) numDimensions * numPoints > MAX_ARRAY_SIZE )
+			throw new IllegalArgumentException( "positions[][] is too large to be stored in a flat array" );
+		final double[] reordered = new double[ numDimensions * numPoints ];
+		for ( int i = 0; i < numPoints; ++i )
+			for ( int d = 0; d < numDimensions; ++d )
+				reordered[ numDimensions * i + d ] = positions[ d ][ tree[ i ] ];
+		return reordered;
+	}
+
+	/**
+	 * Flatten the nested {@code positions} array.
+	 *
+	 * @param positions positions in nested layout
+	 * @return positions in flattened layout
+	 */
+	static double[] flatten( double[][] positions )
+	{
+		final int numDimensions = positions.length;
+		final int numPoints = positions[ 0 ].length;
+		final double[] flattened = new double[ numDimensions * numPoints ];
+		for ( int i = 0; i < numPoints; ++i )
+			for ( int d = 0; d < numDimensions; ++d )
+				flattened[ numDimensions * i + d ] = positions[ d ][ i ];
+		return flattened;
+
+	}
+
+	/**
+	 * Transform flat {@code positions} array into a nested {@code
+	 * double[numDimensions][numPoints]} array.
+	 * <p>
+	 * With flat layout, positions are stored as a flat {@code double[]} array,
+	 * where {@code positions[d + i*n]} is dimension {@code d} of the {@code
+	 * i}-th point, with {@code n} the number of dimensions.
+	 * <p>
+	 * With nested layout, positions are stored as a nested {@code double[][]}
+	 * array where {@code positions[d][i]} is dimension {@code d} of the {@code
+	 * i}-th point.
+	 *
+	 * @param positions
+	 * 		positions in flattened layout
+	 * @param n
+	 * 		number of dimensions
+	 *
+	 * @return positions in nested layout
+	 */
+	static double[][] unflatten( double[] positions, final int n )
+	{
+		final int numPoints = positions.length / n;
+		final double[][] unflattened = new double[ n ][ numPoints ];
+		for (int i = 0; i < positions.length; ++i )
+		{
+			final int d = i % n;
+			unflattened[ d ][ i / n + d ] = positions[ i ];
+		}
+		return unflattened;
+
+	}
+
+	/**
+	 * Compute bounding box of positions in nested layout
+	 */
+	static void computeMinMax(final double[][] positions, final double[] min, final double[] max) {
+		final int n = min.length;
+		for (int d = 0; d < n; d++) {
+			double maxd = Double.NEGATIVE_INFINITY;
+			double mind = Double.POSITIVE_INFINITY;
+			for (double v : positions[d]) {
+				if (v < mind) {
+					mind = v;
+				}
+				if (v > maxd) {
+					maxd = v;
+				}
+			}
+			min[d] = mind;
+			max[d] = maxd;
+		}
+	}
+
+	/**
+	 * Compute bounding box of positions in flat layout
+	 */
+	static void computeMinMax(final double[] flatPositions, final double[] min, final double[] max) {
+		final int n = min.length;
+		Arrays.fill(max, Double.NEGATIVE_INFINITY);
+		Arrays.fill(min, Double.POSITIVE_INFINITY);
+		int d = 0;
+		for (double v : flatPositions) {
+			if (v < min[d]) {
+				min[d] = v;
+			}
+			if (v > max[d]) {
+				max[d] = v;
+			}
+			d = (d + 1) % n;
+		}
+	}
+
+	/**
+	 * Invert the given permutation {@code tree}.
+	 * <p>
+	 * For example, {@code tree = {3, 4, 1, 0, 5, 2}} indicates that coordinates
+	 * and value for the node at heap index {@code i} can be found at index
+	 * {@code tree[i]} in the respective input list.
+	 * <p>
+	 * The inverse, {@code inv = {3, 4, 1, 0, 5, 2}} indicates that coordinates
+	 * and value at index {@code i} in the respective input list belong to the
+	 * node at heap index {@code inv[i]}.
+	 *
+	 * @param tree a permutation
+	 * @return the inverse permutation
+	 */
+	static int[] invert( int[] tree )
+	{
+		// For example:
+		// i =          0  1  2  3  4  5
+		// tree =      {3, 4, 1, 0, 5, 2}
+		// output =    {3, 2, 5, 0, 1, 4}
+
+		final int[] inv = new int[ tree.length ];
+		for ( int i = 0; i < tree.length; i++ )
+			inv[tree[i]] = i;
+		return inv;
+	}
+
+	/**
+	 * Re-order the node {@code values} to form a tree corresponding to the index array {@code tree'={0,1,2,...}}.
+	 * The tree is given as an {@link #invert(int[]) inverted permutation}, so that we can iterate through the {@code values} in order, putting each at the right index in the returned {@code List}.
+	 *
+	 * @param invtree
+	 * @param values
+	 * @param <T>
+	 *
+	 * @return
+	 */
+	static < T > List< T > orderValuesList(
+			final int[] invtree,
+			final Iterable< T > values )
+	{
+		final int size = invtree.length;
+		@SuppressWarnings( "unchecked" )
+		final T[] orderedValues = ( T[] ) new Object[ size ];
+		final Iterator< T > ival = values.iterator();
+		for ( final int i : invtree )
+		{
+			if ( !ival.hasNext() )
+				throw new IllegalArgumentException( "provided values Iterable has fewer elements than required" );
+			orderedValues[ i ] = ival.next();
+		}
+		return Arrays.asList( orderedValues );
+	}
+
+	/**
+	 * Re-order the node {@code values} to form a tree corresponding to the index array {@code tree'={0,1,2,...}}.
+	 * The tree is given as an {@link #invert(int[]) inverted permutation}, so that we can iterate through the {@code values} in order, putting each at the right index in the returned 1D {@code Img}.
+	 *
+	 * @param invtree
+	 * @param values
+	 * @param <T>
+	 *
+	 * @return
+	 */
+	static < T extends NativeType< T > > Img< T > orderValuesImg(
+			final int[] invtree,
+			final Iterable< T > values )
+	{
+		final int size = invtree.length;
+		final Img< T > img = new ArrayImgFactory<>( getType( values ) ).create( size );
+		final RandomAccess< T > orderedValues = img.randomAccess();
+		final Iterator< T > ival = values.iterator();
+		for ( final int i : invtree )
+		{
+			if ( !ival.hasNext() )
+				throw new IllegalArgumentException( "provided values Iterable has fewer elements than required" );
+			orderedValues.setPositionAndGet( i ).set( ival.next() );
+		}
+		return img;
+	}
+
+	/**
+	 * Returns the first element of {@code values}.
+	 *
+	 * @throws IllegalArgumentException
+	 * 		if {@code values} has no elements.
+	 */
+	static < T > T getType( Iterable< T > values )
+	{
+		final Iterator< T > ival = values.iterator();
+		if ( !ival.hasNext() )
+			throw new IllegalArgumentException( "values Iterable is empty" );
+		return ival.next();
+	}
+
+	/**
+	 * Returns the number of dimensions of the first element of {@code positions}.
+	 * If {@code positions} has no elements, returns {@code 0}.
+	 *
+	 * @param positions
+	 * 		list of points
+	 *
+	 * @return number of dimensions of the first point
+	 */
+	static int getNumDimensions( Iterable< ? extends RealLocalizable > positions )
+	{
+		final Iterator< ? extends RealLocalizable > ipos = positions.iterator();
+		return ipos.hasNext() ? ipos.next().numDimensions() : 0;
+	}
+
+	/**
+	 * Swap {@code order[i]} and {@code order[j]}.
+	 */
+	private static void swap( final int i, final int j, final int[] order )
+	{
+		final int tmp = order[ i ];
+		order[ i ] = order[ j ];
+		order[ j ] = tmp;
+	}
+
+	/**
+	 * Partition a sublist. The list is given by an immutable array of {@code
+	 * values}, and an index array that represents the {@code order} of values
+	 * in the list. This method only rearranges the {@code order} array.
+	 * <p>
+	 * A pivot element is chosen by median-of-three method. Then {@code
+	 * [i,j]} is reordered, such that all elements before the pivot are
+	 * smaller-equal and all elements after the pivot are larger-equal the
+	 * pivot. The index of the pivot element is returned.
+	 *
+	 * @param i
+	 *            index of first element of the sublist
+	 * @param j
+	 *            index of last element of the sublist
+	 * @param values
+	 * 			  the array of values of list elements.
+	 * @param order
+	 *            order of list elements. E.g., {@code order[0]=3} means that {@code values[3]} is the first element of the list.
+	 * @return index of pivot element
+	 */
+	static int partition( int i, int j, final double[] values, final int[] order )
+	{
+		final int len = j - i + 1;
+		if ( len <= 2 )
+		{
+			if ( len <= 0 )
+				throw new IllegalArgumentException();
+			if ( values[ order[ i ] ] > values[ order[ j ] ] )
+				swap( i, j, order );
+			return i;
+		}
+		else
+		{
+			final int m = ( i + j ) / 2;
+			if ( values[ order[ i ] ] > values[ order[ m ] ] )
+				swap( i, m, order );
+			if ( values[ order[ i ] ] > values[ order[ j ] ] )
+				swap( i, j, order );
+			if ( values[ order[ m ] ] > values[ order[ j ] ] )
+				swap( m, j, order );
+			swap( m, i + 1, order );
+			final int p = ++i;
+			final double pivot = values[ order[ p ] ];
+			while ( true )
+			{
+				while ( values[ order[ ++i ] ] < pivot )
+					;
+				while ( values[ order[ --j ] ] > pivot )
+					;
+				if ( j < i )
+					break;
+				swap( i, j, order );
+			}
+			swap( p, j, order );
+			return j;
+		}
+	}
+
+	/**
+	 * Sort a sublist. The list is given by an immutable array of {@code
+	 * values}, and an index array that represents the {@code order} of values
+	 * in the list. This method only rearranges the {@code order} array.
+	 *
+	 * @param i
+	 *            index of first element of the sublist
+	 * @param j
+	 *            index of last element of the sublist
+	 * @param values
+	 * 			  the array of values of list elements.
+	 * @param order
+	 *            order of list elements. E.g., {@code order[0]=3} means that {@code values[3]} is the first element of the list.
+	 */
+	static void quicksort( final int i, final int j, final double[] values, final int[] order )
+	{
+		if ( 0 <= i && i < j )
+		{
+			final int p = partition( i, j, values, order );
+			quicksort( i, p - 1, values, order );
+			quicksort( p + 1, j, values, order );
+		}
+	}
+
+	private static final class MakeTree
+	{
+		private final int numDimensions;
+
+		private final int numPoints;
+
+		/**
+		 * The coordinates for the {@code i}th point are stored at {@code positions[d][k]} where {@code d} is the dimension.
+		 */
+		private final double[][] positions;
+
+		/**
+		 * Temporary array to keep track of elements.
+		 * Initialized to {@code 0, 1, ... } and then permuted when sorting the elements into a tree.
+		 */
+		private final int[] indices;
+
+		/**
+		 * Node indices in a flattened (heap-like) array.
+		 * For example: the index of the root node is {@code tree[0]}.
+		 * The coordinates of the root node are at {@code positions[d][tree[0]]} where {@code d} is the dimension.
+		 * The children of the root node are at {@code tree[1]} and {@code tree[2]}, and so on.
+		 */
+		private final int[] tree;
+
+		private MakeTree( final double[][] positions )
+		{
+			this.positions = positions;
+			numDimensions = positions.length;
+			numPoints = positions[ 0 ].length;
+			indices = new int[ numPoints ];
+			tree = new int[ numPoints ];
+			Arrays.setAll( indices, j -> j );
+			makeNode( 0, numPoints - 1, 0, 0 );
+		}
+
+		/**
+		 * Calculate pivot index such that the tree will be arranged in a way that
+		 * "leaf layers" are filled from the left.
+		 * For example 10 nodes will always be arranged like this:
+		 * <pre>
+		 *            0
+		 *         /     \
+		 *       1         2
+		 *     /   \     /   \
+		 *    3     4   5     6
+		 *   / \   /
+		 *  7   8 9
+		 * </pre>
+		 *
+		 * never like this:
+		 * <pre>
+		 *            0
+		 *         /     \
+		 *       1         2
+		 *     /   \     /   \
+		 *    3     4   5     6
+		 *   /         /     /
+		 *  7         8     9
+		 * </pre>
+		 *
+		 * By choosing pivots in this way, the tree structure is fully
+		 * determined. For every node index, the child indices can be calculated
+		 * without dependent reads. And iff the calculated child index is less
+		 * than the number of nodes, the child exists.
+		 */
+		private static int pivot( final int len )
+		{
+			final int h = Integer.highestOneBit( len );
+			final int h2 = h >> 1;
+			return ( len - h >= h2 )
+					? h - 1
+					: len - h2;
+		}
+
+		private void makeNode( final int i, final int j, final int d, final int nodeIndex )
+		{
+			if ( j > i )
+			{
+				final int k = i + pivot( j - i + 1 );
+				kthElement( i, j, k, d );
+				tree[ nodeIndex ] = indices[ k ];
+				final int dChild = ( d + 1 ) % numDimensions;
+				makeNode( i, k - 1, dChild, leftChildIndex( nodeIndex ) );
+				makeNode( k + 1, j, dChild, rightChildIndex( nodeIndex ) );
+			}
+			else if ( j == i )
+			{
+				tree[ nodeIndex ] = indices[ i ];
+			}
+		}
+
+		/**
+		 * Partition a sublist of Nodes by their coordinate in the specified
+		 * dimension, such that the k-th smallest value is at position {@code
+		 * k}, elements before the k-th are smaller or equal and elements after
+		 * the k-th are larger or equal.
+		 *
+		 * @param i
+		 *            index of first element of the sublist
+		 * @param j
+		 *            index of last element of the sublist
+		 * @param k
+		 *            index for k-th smallest value. {@code i <= k <= j}.
+		 * @param compare_d
+		 *            dimension by which to order the sublist
+		 */
+		private void kthElement( int i, int j, final int k, final int compare_d )
+		{
+			while ( true )
+			{
+				final int pivotpos = partition( i, j, positions[ compare_d ], indices );
+				if ( pivotpos > k )
+				{
+					// partition lower half
+					j = pivotpos - 1;
+				}
+				else if ( pivotpos < k )
+				{
+					// partition upper half
+					i = pivotpos + 1;
+				}
+				else
+					return;
+			}
+		}
+	}
+
+	private KDTreeUtils() {}
+}
diff --git a/src/main/java/net/imglib2/kdtree/KNearestNeighborSearchImpl.java b/src/main/java/net/imglib2/kdtree/KNearestNeighborSearchImpl.java
new file mode 100644
index 000000000..fccec79bb
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KNearestNeighborSearchImpl.java
@@ -0,0 +1,125 @@
+package net.imglib2.kdtree;
+
+import java.util.Arrays;
+import net.imglib2.RealLocalizable;
+
+/**
+ * <em>k</em>-nearest-neighbor search on {@link KDTreeImpl}.
+ * Results are node indices.
+ */
+public class KNearestNeighborSearchImpl
+{
+	private final KDTreeImpl tree;
+	private final int numDimensions;
+	private final int numPoints;
+	private final double[] pos;
+	private final int k;
+	private final double[] bestSquDistance;
+	private final int[] bestIndex;
+	private final double[] axisDiffs;
+	private final int[] awayChilds;
+
+	public KNearestNeighborSearchImpl( final KDTreeImpl tree, final int k )
+	{
+		this.tree = tree;
+		numDimensions = tree.numDimensions();
+		numPoints = tree.size();
+		this.k = k;
+		pos = new double[ numDimensions ];
+		bestSquDistance = new double[ k ];
+		bestIndex = new int[ k ];
+		final int depth = tree.depth();
+		axisDiffs = new double[ depth + 1 ];
+		awayChilds = new int[ depth + 1 ];
+	}
+
+	/**
+	 * Insert index into list of best nodes.
+	 * Also checks whether index will be inserted at all, that is,
+	 * whether squDistance < bestSquDistance[k-1]
+	 */
+	private void insert( final double squDistance, final int index )
+	{
+		// first check whether index will be inserted at all
+		if ( squDistance < bestSquDistance[ k - 1 ] )
+		{
+
+			// find insertion point, shifting existing elements to make room
+			int i;
+			for ( i = k - 1; i > 0 && squDistance < bestSquDistance[ i - 1 ]; --i )
+			{
+				bestSquDistance[ i ] = bestSquDistance[ i - 1 ];
+				bestIndex[ i ] = bestIndex[ i - 1 ];
+			}
+
+			// insert index at i,
+			bestSquDistance[ i ] = squDistance;
+			bestIndex[ i ] = index;
+		}
+	}
+
+	public void search( final RealLocalizable p )
+	{
+		p.localize( pos );
+		int current = tree.root();
+		int depth = 0;
+		Arrays.fill( bestSquDistance, Double.POSITIVE_INFINITY );
+		Arrays.fill( bestIndex, -1 );
+		while ( true )
+		{
+			insert( tree.squDistance( current, pos ), current );
+
+			final int d = depth % numDimensions;
+			final double axisDiff = pos[ d ] - tree.getDoublePosition( current, d );
+			final boolean leftIsNearBranch = axisDiff < 0;
+
+			// search the near branch
+			final int nearChild = ( 2 * current ) + ( leftIsNearBranch ? 1 : 2 );
+			final int awayChild = ( 2 * current ) + ( leftIsNearBranch ? 2 : 1 );
+			++depth;
+			awayChilds[ depth ] = awayChild;
+			axisDiffs[ depth ] = axisDiff * axisDiff;
+			if ( nearChild >= numPoints )
+			{
+				while ( awayChilds[ depth ] >= numPoints || axisDiffs[ depth ] > bestSquDistance[ k - 1 ] )
+				{
+					if ( --depth == 0 )
+					{
+						return;
+					}
+				}
+				current = awayChilds[ depth ];
+				awayChilds[ depth ] = numPoints;
+			}
+			else
+			{
+				current = nearChild;
+			}
+		}
+	}
+
+	public int k()
+	{
+		return k;
+	}
+
+	public int bestIndex( final int i )
+	{
+		return bestIndex[ i ];
+	}
+
+	public double bestSquDistance( final int i )
+	{
+		return bestSquDistance[ i ];
+	}
+
+	public KNearestNeighborSearchImpl copy()
+	{
+		final KNearestNeighborSearchImpl copy = new KNearestNeighborSearchImpl( tree, k );
+		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
+		System.arraycopy( bestIndex, 0, copy.bestIndex, 0, bestIndex.length );
+		System.arraycopy( bestSquDistance, 0, copy.bestSquDistance, 0, bestSquDistance.length );
+		return copy;
+	}
+
+}
diff --git a/src/main/java/net/imglib2/kdtree/NearestNeighborSearchImpl.java b/src/main/java/net/imglib2/kdtree/NearestNeighborSearchImpl.java
new file mode 100644
index 000000000..d887bea94
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/NearestNeighborSearchImpl.java
@@ -0,0 +1,95 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RealLocalizable;
+
+/**
+ * Nearest-neighbor search on {@link KDTreeImpl}.
+ * Results are node indices.
+ */
+public class NearestNeighborSearchImpl
+{
+	private final KDTreeImpl tree;
+	private final int numDimensions;
+	private final int numPoints;
+	private final double[] pos;
+	private int bestIndex;
+	private double bestSquDistance;
+	private final double[] axisDiffs;
+	private final int[] awayChilds;
+
+	public NearestNeighborSearchImpl( final KDTreeImpl tree )
+	{
+		this.tree = tree;
+		numDimensions = tree.numDimensions();
+		numPoints = tree.size();
+		pos = new double[ numDimensions ];
+		bestIndex = -1;
+		bestSquDistance = Double.POSITIVE_INFINITY;
+		final int depth = tree.depth();
+		axisDiffs = new double[ depth + 1 ];
+		awayChilds = new int[ depth + 1 ];
+	}
+
+	public void search( final RealLocalizable p )
+	{
+		p.localize( pos );
+		int current = tree.root();
+		int depth = 0;
+		bestSquDistance = ( bestIndex >= 0 ) ? tree.squDistance( bestIndex, pos ) : Double.POSITIVE_INFINITY;
+		while ( true )
+		{
+			final double squDistance = tree.squDistance( current, pos );
+			if ( squDistance < bestSquDistance )
+			{
+				bestSquDistance = squDistance;
+				bestIndex = current;
+			}
+
+			final int d = depth % numDimensions;
+			final double axisDiff = pos[ d ] - tree.getDoublePosition( current, d );
+			final boolean leftIsNearBranch = axisDiff < 0;
+
+			// search the near branch
+			final int nearChild = ( 2 * current ) + ( leftIsNearBranch ? 1 : 2 );
+			final int awayChild = ( 2 * current ) + ( leftIsNearBranch ? 2 : 1 );
+			++depth;
+			awayChilds[ depth ] = awayChild;
+			axisDiffs[ depth ] = axisDiff * axisDiff;
+			if ( nearChild >= numPoints )
+			{
+				while ( awayChilds[ depth ] >= numPoints || axisDiffs[ depth ] > bestSquDistance )
+				{
+					if ( --depth == 0 )
+					{
+						return;
+					}
+				}
+				current = awayChilds[ depth ];
+				awayChilds[ depth ] = numPoints;
+			}
+			else
+			{
+				current = nearChild;
+			}
+		}
+	}
+
+	public int bestIndex()
+	{
+		return bestIndex;
+	}
+
+	public double bestSquDistance()
+	{
+		return bestSquDistance;
+	}
+
+	public NearestNeighborSearchImpl copy()
+	{
+		final NearestNeighborSearchImpl copy = new NearestNeighborSearchImpl( tree );
+		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
+		copy.bestIndex = bestIndex;
+		copy.bestSquDistance = bestSquDistance;
+		return copy;
+	}
+}
diff --git a/src/main/java/net/imglib2/kdtree/RadiusNeighborSearchImpl.java b/src/main/java/net/imglib2/kdtree/RadiusNeighborSearchImpl.java
new file mode 100644
index 000000000..8d45f2264
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/RadiusNeighborSearchImpl.java
@@ -0,0 +1,155 @@
+package net.imglib2.kdtree;
+
+import java.util.Arrays;
+import net.imglib2.RealLocalizable;
+
+import static net.imglib2.kdtree.KDTreeUtils.quicksort;
+import static net.imglib2.kdtree.KDTreeUtils.reorder;
+
+/**
+ * Radius neighbor search on {@link KDTreeImpl}.
+ * Results are node indices.
+ */
+public class RadiusNeighborSearchImpl
+{
+	private final KDTreeImpl tree;
+	private final int numDimensions;
+	private final int numPoints;
+	private final double[] pos;
+	private final double[] axisDiffs;
+	private final int[] awayChilds;
+	private final Neighbors neighbors;
+
+	public RadiusNeighborSearchImpl( final KDTreeImpl tree )
+	{
+		this.tree = tree;
+		numDimensions = tree.numDimensions();
+		numPoints = tree.size();
+		pos = new double[ numDimensions ];
+		final int depth = tree.depth();
+		axisDiffs = new double[ depth + 1 ];
+		awayChilds = new int[ depth + 1 ];
+		neighbors = new Neighbors();
+	}
+
+	public void search( final RealLocalizable p, final double radius, final boolean sortResults )
+	{
+		assert radius >= 0;
+		final double squRadius = radius * radius;
+		p.localize( pos );
+		neighbors.clear();
+		int current = tree.root();
+		int depth = 0;
+		while ( true )
+		{
+			final double squDistance = tree.squDistance( current, pos );
+			if ( squDistance < squRadius )
+			{
+				neighbors.add( squDistance, current );
+			}
+
+			final int d = depth % numDimensions;
+			final double axisDiff = pos[ d ] - tree.getDoublePosition( current, d );
+			final boolean leftIsNearBranch = axisDiff < 0;
+
+			// search the near branch
+			final int nearChild = ( 2 * current ) + ( leftIsNearBranch ? 1 : 2 );
+			final int awayChild = ( 2 * current ) + ( leftIsNearBranch ? 2 : 1 );
+			++depth;
+			awayChilds[ depth ] = awayChild;
+			axisDiffs[ depth ] = axisDiff * axisDiff;
+			if ( nearChild >= numPoints )
+			{
+				while ( awayChilds[ depth ] >= numPoints || axisDiffs[ depth ] > squRadius )
+				{
+					if ( --depth == 0 )
+					{
+						if ( sortResults )
+						{
+							neighbors.sort();
+						}
+						return;
+					}
+				}
+				current = awayChilds[ depth ];
+				awayChilds[ depth ] = numPoints;
+			}
+			else
+			{
+				current = nearChild;
+			}
+		}
+	}
+
+	public int numNeighbors()
+	{
+		return neighbors.size;
+	}
+
+	public int bestIndex( final int i )
+	{
+		return neighbors.indices[ i ];
+	}
+
+	public double bestSquDistance( final int i )
+	{
+		return neighbors.distances[ i ];
+	}
+
+	public RadiusNeighborSearchImpl copy()
+	{
+		final RadiusNeighborSearchImpl copy = new RadiusNeighborSearchImpl( tree );
+		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
+		copy.neighbors.makeCopyOf( neighbors );
+		return copy;
+	}
+
+	static class Neighbors
+	{
+		double[] distances;
+		int[] indices;
+		int size;
+
+		Neighbors()
+		{
+			final int capacity = 10;
+			distances = new double[ capacity ];
+			indices = new int[ capacity ];
+		}
+
+		void clear()
+		{
+			size = 0;
+		}
+
+		void add( final double distance, final int index )
+		{
+			if ( distances.length <= size )
+			{
+				// reallocate
+				final int newLength = distances.length * 2;
+				distances = Arrays.copyOf( distances, newLength );
+				indices = Arrays.copyOf( indices, newLength );
+			}
+			distances[ size ] = distance;
+			indices[ size ] = index;
+			++size;
+		}
+
+		void sort()
+		{
+			final int[] order = new int[ size ];
+			Arrays.setAll( order, i -> i );
+			quicksort( 0, size - 1, distances, order );
+			System.arraycopy( reorder( distances, order ), 0, distances, 0, size );
+			System.arraycopy( reorder( indices, order ), 0, indices, 0, size );
+		}
+
+		void makeCopyOf( final Neighbors other )
+		{
+			distances = other.distances.clone();
+			indices = other.indices.clone();
+			size = other.size;
+		}
+	}
+}
diff --git a/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearchOnKDTree.java b/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearchOnKDTree.java
index 621fa6c79..4dcc426d2 100644
--- a/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearchOnKDTree.java
+++ b/src/main/java/net/imglib2/neighborsearch/KNearestNeighborSearchOnKDTree.java
@@ -34,10 +34,12 @@
 
 package net.imglib2.neighborsearch;
 
-import net.imglib2.KDTree;
-import net.imglib2.KDTreeNode;
+import java.util.Arrays;
 import net.imglib2.RealLocalizable;
 import net.imglib2.Sampler;
+import net.imglib2.KDTree;
+import net.imglib2.KDTreeNode;
+import net.imglib2.kdtree.KNearestNeighborSearchImpl;
 
 /**
  * Implementation of {@link KNearestNeighborSearch} search for kd-trees.
@@ -46,35 +48,38 @@
  */
 public class KNearestNeighborSearchOnKDTree< T > implements KNearestNeighborSearch< T >
 {
-	protected KDTree< T > tree;
-
-	protected final int n;
+	private final KDTree< T > tree;
 
-	protected final double[] pos;
+	private final int k;
 
-	protected final int k;
+	private final KNearestNeighborSearchImpl impl;
 
-	protected KDTreeNode< T >[] bestPoints;
-
-	protected double[] bestSquDistances;
+	private final KDTreeNode< T >[] matches;
 
 	@SuppressWarnings( "unchecked" )
 	public KNearestNeighborSearchOnKDTree( final KDTree< T > tree, final int k )
 	{
 		this.tree = tree;
-		this.n = tree.numDimensions();
-		this.pos = new double[ n ];
 		this.k = k;
-		this.bestPoints = new KDTreeNode[ k ];
-		this.bestSquDistances = new double[ k ];
-		for ( int i = 0; i < k; ++i )
-			bestSquDistances[ i ] = Double.MAX_VALUE;
+		impl = new KNearestNeighborSearchImpl( tree.impl(), k );
+		matches = new KDTreeNode[ k ];
+		Arrays.setAll( matches, i -> tree.createNode() );
+	}
+
+	@SuppressWarnings( "unchecked" )
+	private KNearestNeighborSearchOnKDTree( final KNearestNeighborSearchOnKDTree< T > knn )
+	{
+		tree = knn.tree;
+		k = knn.k;
+		impl = knn.impl.copy();
+		matches = new KDTreeNode[ k ];
+		Arrays.setAll( matches, i -> tree.createNode().setNodeIndex( impl.bestIndex( i ) ) );
 	}
 
 	@Override
 	public int numDimensions()
 	{
-		return n;
+		return tree.numDimensions();
 	}
 
 	@Override
@@ -84,67 +89,35 @@ public int getK()
 	}
 
 	@Override
-	public void search( final RealLocalizable reference )
+	public void search( final RealLocalizable p )
 	{
-		reference.localize( pos );
-		for ( int i = 0; i < k; ++i )
-			bestSquDistances[ i ] = Double.MAX_VALUE;
-		searchNode( tree.getRoot() );
-	}
-
-	protected void searchNode( final KDTreeNode< T > current )
-	{
-		// consider the current node
-		final double squDistance = current.squDistanceTo( pos );
-		if ( squDistance < bestSquDistances[ k - 1 ] )
-		{
-			int i = k - 1;
-			for ( int j = i - 1; i > 0 && squDistance < bestSquDistances[ j ]; --i, --j )
-			{
-				bestSquDistances[ i ] = bestSquDistances[ j ];
-				bestPoints[ i ] = bestPoints[ j ];
-			}
-			bestSquDistances[ i ] = squDistance;
-			bestPoints[ i ] = current;
-		}
-
-		final double axisDiff = pos[ current.getSplitDimension() ] - current.getSplitCoordinate();
-		final double axisSquDistance = axisDiff * axisDiff;
-		final boolean leftIsNearBranch = axisDiff < 0;
-
-		// search the near branch
-		final KDTreeNode< T > nearChild = leftIsNearBranch ? current.left : current.right;
-		final KDTreeNode< T > awayChild = leftIsNearBranch ? current.right : current.left;
-		if ( nearChild != null )
-			searchNode( nearChild );
-
-		// search the away branch - maybe
-		if ( ( axisSquDistance <= bestSquDistances[ k - 1 ] ) && ( awayChild != null ) )
-			searchNode( awayChild );
+		impl.search( p );
+		for ( int i = 0; i < k; i++ )
+			matches[ i ].setNodeIndex( impl.bestIndex( i ) );
 	}
 
 	@Override
 	public Sampler< T > getSampler( final int i )
 	{
-		return bestPoints[ i ];
+		return matches[ i ];
 	}
 
 	@Override
 	public RealLocalizable getPosition( final int i )
 	{
-		return bestPoints[ i ];
+		return matches[ i ];
 	}
 
 	@Override
 	public double getSquareDistance( final int i )
 	{
-		return bestSquDistances[ i ];
+		return impl.bestSquDistance( i );
 	}
 
 	@Override
-	public double getDistance( final int i )
+	public KNearestNeighborSearchOnKDTree< T > copy()
 	{
-		return Math.sqrt( bestSquDistances[ i ] );
+		return new KNearestNeighborSearchOnKDTree<>( this );
 	}
 
 	/* NearestNeighborSearch */
@@ -172,17 +145,4 @@ public double getDistance()
 	{
 		return getDistance( 0 );
 	}
-
-	@Override
-	public KNearestNeighborSearchOnKDTree< T > copy()
-	{
-		final KNearestNeighborSearchOnKDTree< T > copy = new KNearestNeighborSearchOnKDTree< T >( tree, k );
-		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
-		for ( int i = 0; i < k; ++i )
-		{
-			copy.bestPoints[ i ] = bestPoints[ i ];
-			copy.bestSquDistances[ i ] = bestSquDistances[ i ];
-		}
-		return copy;
-	}
 }
diff --git a/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearchOnKDTree.java b/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearchOnKDTree.java
index ce8980b1e..d94514d5b 100644
--- a/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearchOnKDTree.java
+++ b/src/main/java/net/imglib2/neighborsearch/NearestNeighborSearchOnKDTree.java
@@ -34,10 +34,11 @@
 
 package net.imglib2.neighborsearch;
 
-import net.imglib2.KDTree;
-import net.imglib2.KDTreeNode;
 import net.imglib2.RealLocalizable;
 import net.imglib2.Sampler;
+import net.imglib2.KDTree;
+import net.imglib2.KDTreeNode;
+import net.imglib2.kdtree.NearestNeighborSearchImpl;
 
 /**
  * Implementation of {@link NearestNeighborSearch} search for kd-trees.
@@ -47,60 +48,38 @@
  */
 public class NearestNeighborSearchOnKDTree< T > implements NearestNeighborSearch< T >
 {
-	protected KDTree< T > tree;
-
-	protected final int n;
+	private final KDTree< T > tree;
 
-	protected final double[] pos;
+	private final NearestNeighborSearchImpl impl;
 
-	protected KDTreeNode< T > bestPoint;
-
-	protected double bestSquDistance;
+	private final KDTreeNode< T > bestPoint;
 
 	public NearestNeighborSearchOnKDTree( final KDTree< T > tree )
 	{
-		n = tree.numDimensions();
-		pos = new double[ n ];
 		this.tree = tree;
+		impl = new NearestNeighborSearchImpl( tree.impl() );
+		bestPoint = tree.createNode();
 	}
 
-	@Override
-	public int numDimensions()
+	private NearestNeighborSearchOnKDTree( final NearestNeighborSearchOnKDTree< T > nn )
 	{
-		return n;
+		tree = nn.tree;
+		impl = nn.impl.copy();
+		bestPoint = tree.createNode();
+		bestPoint.setNodeIndex( nn.impl.bestIndex() );
 	}
 
 	@Override
-	public void search( final RealLocalizable p )
+	public int numDimensions()
 	{
-		p.localize( pos );
-		bestSquDistance = Double.MAX_VALUE;
-		searchNode( tree.getRoot() );
+		return tree.numDimensions();
 	}
 
-	protected void searchNode( final KDTreeNode< T > current )
+	@Override
+	public void search( final RealLocalizable p )
 	{
-		// consider the current node
-		final double distance = current.squDistanceTo( pos );
-		if ( distance < bestSquDistance )
-		{
-			bestSquDistance = distance;
-			bestPoint = current;
-		}
-
-		final double axisDiff = pos[ current.getSplitDimension() ] - current.getSplitCoordinate();
-		final double axisSquDistance = axisDiff * axisDiff;
-		final boolean leftIsNearBranch = axisDiff < 0;
-
-		// search the near branch
-		final KDTreeNode< T > nearChild = leftIsNearBranch ? current.left : current.right;
-		final KDTreeNode< T > awayChild = leftIsNearBranch ? current.right : current.left;
-		if ( nearChild != null )
-			searchNode( nearChild );
-
-		// search the away branch - maybe
-		if ( ( axisSquDistance <= bestSquDistance ) && ( awayChild != null ) )
-			searchNode( awayChild );
+		impl.search( p );
+		bestPoint.setNodeIndex( impl.bestIndex() );
 	}
 
 	@Override
@@ -118,22 +97,12 @@ public RealLocalizable getPosition()
 	@Override
 	public double getSquareDistance()
 	{
-		return bestSquDistance;
-	}
-
-	@Override
-	public double getDistance()
-	{
-		return Math.sqrt( bestSquDistance );
+		return impl.bestSquDistance();
 	}
 
 	@Override
 	public NearestNeighborSearchOnKDTree< T > copy()
 	{
-		final NearestNeighborSearchOnKDTree< T > copy = new NearestNeighborSearchOnKDTree< T >( tree );
-		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
-		copy.bestPoint = bestPoint;
-		copy.bestSquDistance = bestSquDistance;
-		return copy;
+		return new NearestNeighborSearchOnKDTree<>( this );
 	}
 }
diff --git a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
index 909b5402e..4fd670808 100644
--- a/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
+++ b/src/main/java/net/imglib2/neighborsearch/RadiusNeighborSearchOnKDTree.java
@@ -35,124 +35,90 @@
 package net.imglib2.neighborsearch;
 
 import java.util.ArrayList;
-import java.util.Collections;
-import java.util.Comparator;
-
-import net.imglib2.KDTree;
-import net.imglib2.KDTreeNode;
+import java.util.List;
 import net.imglib2.RealLocalizable;
 import net.imglib2.Sampler;
-import net.imglib2.util.ValuePair;
+import net.imglib2.KDTree;
+import net.imglib2.KDTreeNode;
+import net.imglib2.kdtree.RadiusNeighborSearchImpl;
 
 /**
  * Implementation of {@link RadiusNeighborSearch} search for kd-trees.
- * 
+ *
  * @author Tobias Pietzsch
  */
 public class RadiusNeighborSearchOnKDTree< T > implements RadiusNeighborSearch< T >
 {
-	protected KDTree< T > tree;
+	private final KDTree< T > tree;
 
-	protected final int n;
+	private final RadiusNeighborSearchImpl impl;
 
-	protected final double[] pos;
-
-	protected ArrayList< ValuePair< KDTreeNode< T >, Double > > resultPoints;
+	private final List< KDTreeNode< T > > matches;
 
 	public RadiusNeighborSearchOnKDTree( final KDTree< T > tree )
 	{
 		this.tree = tree;
-		this.n = tree.numDimensions();
-		this.pos = new double[ n ];
-		this.resultPoints = new ArrayList< ValuePair< KDTreeNode< T >, Double > >();
+		impl = new RadiusNeighborSearchImpl( tree.impl() );
+		matches = new ArrayList<>();
 	}
 
-	@Override
-	public void search( final RealLocalizable reference, final double radius, final boolean sortResults )
+	private RadiusNeighborSearchOnKDTree( final RadiusNeighborSearchOnKDTree< T > other )
 	{
-		assert radius >= 0;
-		reference.localize( pos );
-		resultPoints.clear();
-		searchNode( tree.getRoot(), radius * radius );
-		if ( sortResults )
-		{
-			Collections.sort( resultPoints, new Comparator< ValuePair< KDTreeNode< T >, Double > >()
-			{
-				@Override
-				public int compare( final ValuePair< KDTreeNode< T >, Double > o1, final ValuePair< KDTreeNode< T >, Double > o2 )
-				{
-					return Double.compare( o1.b, o2.b );
-				}
-			} );
-		}
+		tree = other.tree;
+		impl = other.impl.copy();
+		matches = new ArrayList<>();
+		for ( final KDTreeNode< T > match : other.matches )
+			matches.add( tree.createNode().setNodeIndex( match.nodeIndex() ) );
 	}
 
 	@Override
-	public int numDimensions()
+	public void search( final RealLocalizable reference, final double radius, final boolean sortResults )
 	{
-		return n;
+		impl.search( reference, radius, sortResults );
+		fillMatches( 0, impl.numNeighbors() - 1 );
 	}
 
-	protected void searchNode( final KDTreeNode< T > current, final double squRadius )
+	private void fillMatches( final int first, final int last )
 	{
-		// consider the current node
-		final double squDistance = current.squDistanceTo( pos );
-		if ( squDistance <= squRadius )
-		{
-			resultPoints.add( new ValuePair< KDTreeNode< T >, Double >( current, squDistance ) );
-		}
-
-		final double axisDiff = pos[ current.getSplitDimension() ] - current.getSplitCoordinate();
-		final double axisSquDistance = axisDiff * axisDiff;
-		final boolean leftIsNearBranch = axisDiff < 0;
-
-		// search the near branch
-		final KDTreeNode< T > nearChild = leftIsNearBranch ? current.left : current.right;
-		final KDTreeNode< T > awayChild = leftIsNearBranch ? current.right : current.left;
-		if ( nearChild != null )
-			searchNode( nearChild, squRadius );
+		while ( matches.size() < last + 1 )
+			matches.add( tree.createNode() );
+		for ( int i = 0; i <= last; ++i )
+			matches.get( i ).setNodeIndex( impl.bestIndex( i ) );
+	}
 
-		// search the away branch - maybe
-		if ( ( axisSquDistance <= squRadius ) && ( awayChild != null ) )
-			searchNode( awayChild, squRadius );
+	@Override
+	public int numDimensions()
+	{
+		return tree.numDimensions();
 	}
 
 	@Override
 	public int numNeighbors()
 	{
-		return resultPoints.size();
+		return impl.numNeighbors();
 	}
 
 	@Override
 	public Sampler< T > getSampler( final int i )
 	{
-		return resultPoints.get( i ).a;
+		return matches.get( i );
 	}
 
 	@Override
 	public RealLocalizable getPosition( final int i )
 	{
-		return resultPoints.get( i ).a;
+		return matches.get( i );
 	}
 
 	@Override
 	public double getSquareDistance( final int i )
 	{
-		return resultPoints.get( i ).b;
-	}
-
-	@Override
-	public double getDistance( final int i )
-	{
-		return Math.sqrt( resultPoints.get( i ).b );
+		return impl.bestSquDistance( i );
 	}
 
 	@Override
 	public RadiusNeighborSearchOnKDTree< T > copy()
 	{
-		final RadiusNeighborSearchOnKDTree< T > copy = new RadiusNeighborSearchOnKDTree<>( tree );
-		System.arraycopy( pos, 0, copy.pos, 0, pos.length );
-		copy.resultPoints.addAll( resultPoints );
-		return copy;
+		return new RadiusNeighborSearchOnKDTree<>( this );
 	}
 }
diff --git a/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java b/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
new file mode 100644
index 000000000..79c4e8fdd
--- /dev/null
+++ b/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
@@ -0,0 +1,186 @@
+package net.imglib2.kdtree;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Random;
+import java.util.concurrent.TimeUnit;
+import net.imglib2.KDTree;
+import net.imglib2.RealLocalizable;
+import net.imglib2.RealPoint;
+import net.imglib2.neighborsearch.NearestNeighborSearchOnKDTree;
+import net.imglib2.neighborsearch.KNearestNeighborSearchOnKDTree;
+import net.imglib2.neighborsearch.RadiusNeighborSearchOnKDTree;
+import org.openjdk.jmh.annotations.Benchmark;
+import org.openjdk.jmh.annotations.BenchmarkMode;
+import org.openjdk.jmh.annotations.Mode;
+import org.openjdk.jmh.annotations.OutputTimeUnit;
+import org.openjdk.jmh.annotations.Scope;
+import org.openjdk.jmh.annotations.State;
+import org.openjdk.jmh.runner.Runner;
+import org.openjdk.jmh.runner.RunnerException;
+import org.openjdk.jmh.runner.options.Options;
+import org.openjdk.jmh.runner.options.OptionsBuilder;
+import org.openjdk.jmh.runner.options.TimeValue;
+import org.openjdk.jmh.annotations.Setup;
+
+@State( Scope.Benchmark )
+public class KDTreeBenchmark
+{
+//	@Param({"3"})
+//	public int n;
+//
+//	@Param({"10000", "100000", "1000000"})
+//	public int numDataVertices;
+//
+//	@Param({"1000"})
+//	public int numTestVertices;
+//
+	public int n = 3;
+	public int k = 10;
+	public int radius = 1;
+	public int numDataVertices = 100000;
+	public int numTestVertices = 1000;
+	public double minCoordinateValue = -5;
+	public double maxCoordinateValue = 5;
+
+	List< RealPoint > dataVertices;
+	List< RealPoint > testVertices;
+
+	private KDTree< RealPoint > kdtree;
+
+	@Setup
+	public void setup()
+	{
+		createVertices();
+//		createVerticesSeqTest();
+		kdtree = new KDTree<>( dataVertices, dataVertices );
+//		spoil();
+	}
+
+	public void spoil() {
+		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
+		final int[] tree = KDTreeUtils.makeTree( points );
+		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
+		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final NearestNeighborSearchImpl search = new NearestNeighborSearchImpl( impl );
+		for ( RealPoint testVertex : testVertices )
+			search.search( testVertex );
+	}
+
+	@Benchmark
+	@BenchmarkMode( Mode.AverageTime )
+	@OutputTimeUnit( TimeUnit.MILLISECONDS )
+	public void createKDTree()
+	{
+		new KDTree<>( dataVertices, dataVertices );
+	}
+
+	@Benchmark
+	@BenchmarkMode( Mode.AverageTime )
+	@OutputTimeUnit( TimeUnit.MILLISECONDS )
+	public void nearestNeighborSearch()
+	{
+		final NearestNeighborSearchOnKDTree< RealPoint > kd = new NearestNeighborSearchOnKDTree<>( kdtree );
+		for ( final RealLocalizable t : testVertices )
+		{
+			kd.search( t );
+			kd.getSampler().get();
+		}
+	}
+
+	@Benchmark
+	@BenchmarkMode( Mode.AverageTime )
+	@OutputTimeUnit( TimeUnit.MILLISECONDS )
+	public void kNearestNeighborSearch()
+	{
+		final KNearestNeighborSearchOnKDTree< RealPoint > kd = new KNearestNeighborSearchOnKDTree<>( kdtree, k );
+		for ( final RealLocalizable t : testVertices )
+		{
+			kd.search( t );
+			kd.getSampler().get();
+//			for ( int i = 0; i < k; i++ )
+//			{
+//				kd.getSampler( i ).get();
+//			}
+		}
+	}
+
+	@Benchmark
+	@BenchmarkMode( Mode.AverageTime )
+	@OutputTimeUnit( TimeUnit.MILLISECONDS )
+	public void radiusNeighborSearch()
+	{
+		final RadiusNeighborSearchOnKDTree< RealPoint > kd = new RadiusNeighborSearchOnKDTree<>( kdtree );
+		for ( final RealLocalizable t : testVertices )
+		{
+			kd.search( t, radius, true );
+			for ( int i = 0; i < Math.min( kd.numNeighbors(), k ); i++ )
+			{
+				kd.getSampler( i ).get();
+			}
+		}
+	}
+
+	private void createVertices()
+	{
+		final double[] p = new double[ n ];
+		final double size = ( maxCoordinateValue - minCoordinateValue );
+		final Random rnd = new Random( 4379 );
+		dataVertices = new ArrayList<>();
+		for ( int i = 0; i < numDataVertices; ++i )
+		{
+			for ( int d = 0; d < n; ++d )
+				p[ d ] = rnd.nextDouble() * size + minCoordinateValue;
+			dataVertices.add( new RealPoint( p ) );
+		}
+		testVertices = new ArrayList<>();
+		for ( int i = 0; i < numTestVertices; ++i )
+		{
+			for ( int d = 0; d < n; ++d )
+				p[ d ] = rnd.nextDouble() * 2 * size + minCoordinateValue - size / 2;
+			testVertices.add( new RealPoint( p )  );
+		}
+	}
+
+	private void createVerticesSeqTest()
+	{
+		final double[] p = new double[ n ];
+		final double size = ( maxCoordinateValue - minCoordinateValue );
+		final Random rnd = new Random( 4379 );
+		dataVertices = new ArrayList<>();
+		for ( int i = 0; i < numDataVertices; ++i )
+		{
+			for ( int d = 0; d < n; ++d )
+				p[ d ] = rnd.nextDouble() * size + minCoordinateValue;
+			dataVertices.add( new RealPoint( p ) );
+		}
+		testVertices = new ArrayList<>();
+		for ( int i = 0; i < numTestVertices; ++i )
+		{
+			if ( rnd.nextDouble() < 0.8 )
+			{
+				int d = rnd.nextInt( n );
+				p[ d ] += size / 10000.0;
+			}
+			else
+			{
+				for ( int d = 0; d < n; ++d )
+					p[ d ] = rnd.nextDouble() * 2 * size + minCoordinateValue - size / 2;
+			}
+			testVertices.add( new RealPoint( p )  );
+		}
+	}
+
+	public static void main( final String... args ) throws RunnerException
+	{
+		final Options opt = new OptionsBuilder()
+				.include( KDTreeBenchmark.class.getSimpleName() )
+				.forks( 0 )
+				.warmupIterations( 4 )
+				.measurementIterations( 8 )
+				.warmupTime( TimeValue.milliseconds( 500 ) )
+				.measurementTime( TimeValue.milliseconds( 500 ) )
+				.build();
+		new Runner( opt ).run();
+	}
+}
diff --git a/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java b/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java
new file mode 100644
index 000000000..d3095f98f
--- /dev/null
+++ b/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java
@@ -0,0 +1,151 @@
+package net.imglib2.kdtree;
+
+import java.util.Arrays;
+import java.util.Comparator;
+import java.util.List;
+import net.imglib2.RealLocalizable;
+import net.imglib2.RealPoint;
+import net.imglib2.util.LinAlgHelpers;
+import org.junit.Assert;
+import org.junit.Before;
+
+import java.util.ArrayList;
+import java.util.Random;
+import org.junit.Test;
+
+public class KDTreeImplTest {
+
+	public int n = 3;
+	public int numDataVertices = 100;
+	public int numTestVertices = 10;
+	public double minCoordinateValue = -5;
+	public double maxCoordinateValue = 5;
+
+	public List< RealPoint > dataVertices;
+	public List< RealPoint > testVertices;
+
+    @Before
+    public void init()
+    {
+		final double[] p = new double[ n ];
+		final double size = ( maxCoordinateValue - minCoordinateValue );
+		final Random rnd = new Random( 4379 );
+		dataVertices = new ArrayList<>();
+		for ( int i = 0; i < numDataVertices; ++i )
+		{
+			for ( int d = 0; d < n; ++d )
+				p[ d ] = rnd.nextDouble() * size + minCoordinateValue;
+			dataVertices.add( new RealPoint( p ) );
+		}
+		testVertices = new ArrayList<>();
+		for ( int i = 0; i < numTestVertices; ++i )
+		{
+			for ( int d = 0; d < n; ++d )
+				p[ d ] = rnd.nextDouble() * 2 * size + minCoordinateValue - size / 2;
+			testVertices.add( new RealPoint( p ) );
+		}
+    }
+
+	@Test
+	public void testNearestNeighborSearch()
+	{
+		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
+		final int[] tree = KDTreeUtils.makeTree( points );
+		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
+		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final NearestNeighborSearchImpl search = new NearestNeighborSearchImpl( impl );
+
+		for ( RealPoint testVertex : testVertices )
+		{
+			final int expected = findNearestNeighborExhaustive( testVertex );
+
+			search.search( testVertex );
+			final int actual = tree[ search.bestIndex() ];
+
+//			System.out.println( "actual = " + actual + ", expected = " + expected );
+			Assert.assertEquals( expected, actual );
+		}
+	}
+
+	@Test
+	public void testKNearestNeighborSearch()
+	{
+		final int k = 10;
+
+		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
+		final int[] tree = KDTreeUtils.makeTree( points );
+		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
+		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final KNearestNeighborSearchImpl search = new KNearestNeighborSearchImpl( impl, k );
+
+		for ( RealPoint testVertex : testVertices )
+		{
+			final int[] expecteds = findNearestNeighborsExhaustive( testVertex, k );
+
+			search.search( testVertex );
+			final int[] actuals = new int[ k ];
+			Arrays.setAll( actuals, i -> tree[ search.bestIndex( i ) ] );
+
+//			System.out.println( "actual = " + actual + ", expected = " + expected );
+			Assert.assertArrayEquals( expecteds, actuals );
+		}
+	}
+
+	@Test
+	public void testRadiusNeighborSearch()
+	{
+		final double radius = 7;
+
+		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
+		final int[] tree = KDTreeUtils.makeTree( points );
+		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
+		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final RadiusNeighborSearchImpl search = new RadiusNeighborSearchImpl( impl );
+
+		for ( RealPoint testVertex : testVertices )
+		{
+			final int[] expecteds = findRadiusNeighborsExhaustive( testVertex, radius );
+
+			search.search( testVertex, radius, true );
+			final int[] actuals = new int[ search.numNeighbors() ];
+			Arrays.setAll( actuals, i -> tree[ search.bestIndex( i ) ] );
+
+			Assert.assertArrayEquals( expecteds, actuals );
+		}
+	}
+
+	private int findNearestNeighborExhaustive( final RealLocalizable point )
+	{
+		return findNearestNeighborsExhaustive( point, 1 )[ 0 ];
+	}
+
+	private int[] findNearestNeighborsExhaustive( final RealLocalizable point, final int k )
+	{
+		final List< RealPoint > sorted = new ArrayList<>();
+		sorted.addAll( dataVertices );
+		sorted.sort( Comparator.comparing( p -> distance( point, p ) ) );
+		final int[] neighbors = new int[ k ];
+		Arrays.setAll(neighbors, i -> dataVertices.indexOf( sorted.get( i ) ) );
+		return neighbors;
+	}
+
+	private static double distance( final RealLocalizable p1, final RealLocalizable p2 )
+	{
+		return LinAlgHelpers.distance( p2.positionAsDoubleArray(), p1.positionAsDoubleArray() );
+	}
+
+	private int[] findRadiusNeighborsExhaustive( final RealLocalizable point, final double radius )
+	{
+		final List< RealPoint > sorted = new ArrayList<>();
+		dataVertices.forEach( p -> {
+			if ( distance( point, p ) <= radius )
+			{
+				sorted.add( p );
+			}
+		} );
+		sorted.sort( Comparator.comparing( p -> distance( point, p ) ) );
+		final int[] neighbors = new int[ sorted.size() ];
+		Arrays.setAll(neighbors, i -> dataVertices.indexOf( sorted.get( i ) ) );
+		return neighbors;
+	}
+}

From 17a51f8db86f49d4920cfdf81d15aab84d96c761 Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Sat, 15 Jul 2023 00:07:35 +0200
Subject: [PATCH 05/13] Replace explicit type arguments with <>

---
 .../imglib2/nearestneighbor/KDTreeTest.java   | 28 +++++++++----------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java b/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java
index 5ff8cb5a5..c244adf78 100644
--- a/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java
+++ b/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java
@@ -58,7 +58,7 @@ public class KDTreeTest
 {
 	protected static boolean testNearestNeighbor( final int numDimensions, final int numPoints, final int numTests, final float min, final float max )
 	{
-		final ArrayList< RealPoint > points = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > points = new ArrayList<>();
 		final Random rnd = new Random( 435435435 );
 
 		final float[] p = new float[ numDimensions ];
@@ -75,13 +75,13 @@ protected static boolean testNearestNeighbor( final int numDimensions, final int
 		}
 
 		long start = System.currentTimeMillis();
-		final KDTree< RealPoint > kdTree = new KDTree< RealPoint >( points, points );
-		final NearestNeighborSearchOnKDTree< RealPoint > kd = new NearestNeighborSearchOnKDTree< RealPoint >( kdTree );
+		final KDTree< RealPoint > kdTree = new KDTree<>( points, points );
+		final NearestNeighborSearchOnKDTree< RealPoint > kd = new NearestNeighborSearchOnKDTree<>( kdTree );
 		final long kdSetupTime = System.currentTimeMillis() - start;
 		System.out.println( "kdtree setup took: " + ( kdSetupTime ) + " ms." );
 
 		start = System.currentTimeMillis();
-		final ArrayList< RealPoint > testpoints = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > testpoints = new ArrayList<>();
 		for ( int i = 0; i < numTests; ++i )
 		{
 			for ( int d = 0; d < numDimensions; ++d )
@@ -163,7 +163,7 @@ private static RealPoint findNearestNeighborExhaustive( final ArrayList< RealPoi
 
 	protected static boolean testKNearestNeighbor( final int neighbors, final int numDimensions, final int numPoints, final int numTests, final float min, final float max )
 	{
-		final ArrayList< RealPoint > points = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > points = new ArrayList<>();
 		final Random rnd = new Random( 435435435 );
 
 		final float[] p = new float[ numDimensions ];
@@ -180,13 +180,13 @@ protected static boolean testKNearestNeighbor( final int neighbors, final int nu
 		}
 
 		long start = System.currentTimeMillis();
-		final KDTree< RealPoint > kdTree = new KDTree< RealPoint >( points, points );
-		final KNearestNeighborSearchOnKDTree< RealPoint > kd = new KNearestNeighborSearchOnKDTree< RealPoint >( kdTree, neighbors );
+		final KDTree< RealPoint > kdTree = new KDTree<>( points, points );
+		final KNearestNeighborSearchOnKDTree< RealPoint > kd = new KNearestNeighborSearchOnKDTree<>( kdTree, neighbors );
 		final long kdSetupTime = System.currentTimeMillis() - start;
 		System.out.println( "kdtree setup took: " + ( kdSetupTime ) + " ms." );
 
 		start = System.currentTimeMillis();
-		final ArrayList< RealPoint > testpoints = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > testpoints = new ArrayList<>();
 		for ( int i = 0; i < numTests; ++i )
 		{
 			for ( int d = 0; d < numDimensions; ++d )
@@ -287,7 +287,7 @@ private static RealPoint[] findKNearestNeighborExhaustive( final ArrayList< Real
 
 	protected static boolean testRadiusNeighbor( final int numDimensions, final int numPoints, final int numTests, final float min, final float max )
 	{
-		final ArrayList< RealPoint > points = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > points = new ArrayList<>();
 		final Random rnd = new Random( 435435435 );
 
 		final float[] p = new float[ numDimensions ];
@@ -306,13 +306,13 @@ protected static boolean testRadiusNeighbor( final int numDimensions, final int
 		final double radius = rnd.nextDouble() * size / 10;
 
 		long start = System.currentTimeMillis();
-		final KDTree< RealPoint > kdTree = new KDTree< RealPoint >( points, points );
-		final RadiusNeighborSearchOnKDTree< RealPoint > kd = new RadiusNeighborSearchOnKDTree< RealPoint >( kdTree );
+		final KDTree< RealPoint > kdTree = new KDTree<>( points, points );
+		final RadiusNeighborSearchOnKDTree< RealPoint > kd = new RadiusNeighborSearchOnKDTree<>( kdTree );
 		final long kdSetupTime = System.currentTimeMillis() - start;
 		System.out.println( "kdtree setup took: " + ( kdSetupTime ) + " ms." );
 
 		start = System.currentTimeMillis();
-		final ArrayList< RealPoint > testpoints = new ArrayList< RealPoint >();
+		final ArrayList< RealPoint > testpoints = new ArrayList<>();
 		for ( int i = 0; i < numTests; ++i )
 		{
 			for ( int d = 0; d < numDimensions; ++d )
@@ -384,7 +384,7 @@ protected static boolean testRadiusNeighbor( final int numDimensions, final int
 
 	private static ArrayList< ValuePair< RealPoint, Double > > findNeighborsRadiusExhaustive( final ArrayList< RealPoint > points, final RealPoint t, final double radius, final boolean sortResults )
 	{
-		final ArrayList< ValuePair< RealPoint, Double > > withinRadius = new ArrayList< ValuePair< RealPoint, Double > >();
+		final ArrayList< ValuePair< RealPoint, Double > > withinRadius = new ArrayList<>();
 
 		final int n = t.numDimensions();
 		final float[] tpos = new float[ n ];
@@ -400,7 +400,7 @@ private static ArrayList< ValuePair< RealPoint, Double > > findNeighborsRadiusEx
 			dist = Math.sqrt( dist );
 
 			if ( dist <= radius )
-				withinRadius.add( new ValuePair< RealPoint, Double >( p, dist ) );
+				withinRadius.add( new ValuePair<>( p, dist ) );
 		}
 
 		if ( sortResults )

From 053d8023c9cbf5658fb0bf0bd84b6e514b703aa6 Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Fri, 14 Jul 2023 22:58:59 +0200
Subject: [PATCH 06/13] Move neighbor search related tests to
 net,imglib2.neighborsearch

---
 .../imglib2/{nearestneighbor => neighborsearch}/KDTreeTest.java | 2 +-
 .../NearestNeighborSearchOnIterableRealIntervalTest.java        | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename src/test/java/net/imglib2/{nearestneighbor => neighborsearch}/KDTreeTest.java (99%)
 rename src/test/java/net/imglib2/{nearestneighbor => neighborsearch}/NearestNeighborSearchOnIterableRealIntervalTest.java (99%)

diff --git a/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java b/src/test/java/net/imglib2/neighborsearch/KDTreeTest.java
similarity index 99%
rename from src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java
rename to src/test/java/net/imglib2/neighborsearch/KDTreeTest.java
index c244adf78..0d9fb8d08 100644
--- a/src/test/java/net/imglib2/nearestneighbor/KDTreeTest.java
+++ b/src/test/java/net/imglib2/neighborsearch/KDTreeTest.java
@@ -32,7 +32,7 @@
  * #L%
  */
 
-package net.imglib2.nearestneighbor;
+package net.imglib2.neighborsearch;
 
 import static org.junit.Assert.assertTrue;
 
diff --git a/src/test/java/net/imglib2/nearestneighbor/NearestNeighborSearchOnIterableRealIntervalTest.java b/src/test/java/net/imglib2/neighborsearch/NearestNeighborSearchOnIterableRealIntervalTest.java
similarity index 99%
rename from src/test/java/net/imglib2/nearestneighbor/NearestNeighborSearchOnIterableRealIntervalTest.java
rename to src/test/java/net/imglib2/neighborsearch/NearestNeighborSearchOnIterableRealIntervalTest.java
index ab05d1444..9e82adca7 100644
--- a/src/test/java/net/imglib2/nearestneighbor/NearestNeighborSearchOnIterableRealIntervalTest.java
+++ b/src/test/java/net/imglib2/neighborsearch/NearestNeighborSearchOnIterableRealIntervalTest.java
@@ -32,7 +32,7 @@
  * #L%
  */
 
-package net.imglib2.nearestneighbor;
+package net.imglib2.neighborsearch;
 
 import static org.junit.Assert.assertTrue;
 import net.imglib2.RealCursor;

From 72b7e9c8bac8fef4e9d200555162003d11511b9c Mon Sep 17 00:00:00 2001
From: tpietzsch <tobias.pietzsch@gmail.com>
Date: Sat, 15 Jul 2023 00:47:25 +0200
Subject: [PATCH 07/13] fix javadoc error

---
 src/main/java/net/imglib2/kdtree/KDTreeData.java | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index 761fea25a..ca65d42ac 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -210,7 +210,7 @@ public Supplier< IntFunction< T > > valuesSupplier()
 
 	/**
 	 * Get positions of points in the tree as a flat {@code double[]} where
-	 * @code positions[d + i*n]} is dimension {@code d} of the {@code i}-th
+	 * {@code positions[d + i*n]} is dimension {@code d} of the {@code i}-th
 	 * point, with {@code n} the number of dimensions.
 	 * <p>
 	 * For serialisation and usage by the tree.

From 720267b1a72c39bb9e38939590d2f7bb6bda5785 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Thu, 10 Aug 2023 16:54:34 -0400
Subject: [PATCH 08/13] Carve out values from KDTreeData

---
 .../java/net/imglib2/kdtree/KDTreeData.java   | 46 +++------------
 .../java/net/imglib2/kdtree/KDTreeValues.java | 53 ++++++++++++++++++
 .../net/imglib2/kdtree/KDTreeValuesImg.java   | 56 +++++++++++++++++++
 .../net/imglib2/kdtree/KDTreeValuesList.java  | 51 +++++++++++++++++
 4 files changed, 168 insertions(+), 38 deletions(-)
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeValues.java
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreeValuesList.java

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index ca65d42ac..7f8218577 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -59,11 +59,7 @@ public enum PositionsLayout
 	private final double[][] positions;
 	private final double[] flatPositions;
 
-	private final List< T > valuesList;
-	private final RandomAccessibleInterval< T > valuesImg;
-	private final Supplier< IntFunction< T > > valuesSupplier;
-
-	private final T type;
+	private final KDTreeValues<T> values;
 
 	private volatile RealInterval boundingBox;
 
@@ -76,12 +72,7 @@ public KDTreeData( double[][] positions, List< T > values )
 		this.positions = positions;
 		flatPositions = null;
 
-		valuesList = values;
-		valuesImg = null;
-		final IntFunction< T > v = values::get;
-		valuesSupplier = () -> v;
-
-		type = KDTreeUtils.getType( values );
+		this.values = new KDTreeValuesList<>(values);
 	}
 
 	public KDTreeData( double[][] positions, List< T > values, RealInterval boundingBox )
@@ -99,14 +90,7 @@ public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
 		this.positions = positions;
 		flatPositions = null;
 
-		valuesList = null;
-		valuesImg = values;
-		valuesSupplier = () -> {
-			final RandomAccess<T> ra = valuesImg.randomAccess();
-			return i -> ra.setPositionAndGet(i);
-		};
-
-		type = Util.getTypeFromInterval( values );
+		this.values = new KDTreeValuesImg<>(values);
 	}
 
 	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
@@ -124,12 +108,7 @@ public KDTreeData( double[] positions, List< T > values )
 		this.positions = null;
 		flatPositions = positions;
 
-		valuesList = values;
-		valuesImg = null;
-		final IntFunction< T > v = values::get;
-		valuesSupplier = () -> v;
-
-		type = KDTreeUtils.getType( values );
+		this.values = new KDTreeValuesList<>(values);
 	}
 
 	public KDTreeData( double[] positions, List< T > values, RealInterval boundingBox )
@@ -147,14 +126,7 @@ public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
 		this.positions = null;
 		flatPositions = positions;
 
-		valuesList = null;
-		valuesImg = values;
-		valuesSupplier = () -> {
-			final RandomAccess<T> ra = valuesImg.randomAccess();
-			return i -> ra.setPositionAndGet(i);
-		};
-
-		type = Util.getTypeFromInterval( values );
+		this.values = new KDTreeValuesImg<>(values);
 	}
 
 	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
@@ -166,7 +138,7 @@ public KDTreeData( double[] positions, RandomAccessibleInterval< T > values, Rea
 	// TODO could also be Class<T> instead? What is more useful?
 	public T type()
 	{
-		return type;
+		return values.type();
 	}
 
 	/**
@@ -191,9 +163,7 @@ public double[][] positions()
 	 */
 	public RandomAccessibleInterval< T > values()
 	{
-		return valuesImg != null
-				? valuesImg
-				: ListImg.wrap( valuesList, size() );
+		return values.values();
 	}
 
 	/**
@@ -205,7 +175,7 @@ public RandomAccessibleInterval< T > values()
 	 */
 	public Supplier< IntFunction< T > > valuesSupplier()
 	{
-		return valuesSupplier;
+		return values.valuesSupplier();
 	}
 
 	/**
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValues.java b/src/main/java/net/imglib2/kdtree/KDTreeValues.java
new file mode 100644
index 000000000..ab66189ac
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValues.java
@@ -0,0 +1,53 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RandomAccessibleInterval;
+
+import java.util.function.IntFunction;
+import java.util.function.Supplier;
+
+/**
+ * Stores the KDTree values.
+ * <p>
+ * Values (of type {@code T}) are stored as either a 1D {@code
+ * RandomAccessibleInterval<T>}, or a {@code List<T>}. Individual values can be
+ * accessed by {@link #valuesSupplier()}{@code .get().apply(i)}. {@code
+ * valueSupplier().get()} returns a reusable {@code IntFunction<T>}. Here {@code
+ * T} maybe a proxy that is reused in subsequent {@code apply(i)}.
+ * <p>
+ * {@link #values()} returns all values as a 1D {@code
+ * RandomAccessibleInterval<T>}. (If data is stored as {@code List<T>}, it is
+ * wrapped into a {@code ListImg}.)
+ *
+ * @param <T>
+ * 		the type of values stored in the tree.
+ */
+abstract public class KDTreeValues<T> {
+	protected final T type;
+
+	public KDTreeValues(T type) {
+		this.type = type;
+	}
+
+	public T type() {
+		return type;
+	}
+
+	/**
+	 * @return number of values in the tree
+	 */
+	abstract public int size();
+
+	/**
+	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
+	 */
+	abstract public RandomAccessibleInterval<T> values();
+
+	/**
+	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
+	 * values for a given node indices. If the returned {@code IntFunction<T>}
+	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
+	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
+	 * {@code IntFunction<T>}.
+	 */
+	abstract public Supplier<IntFunction<T>> valuesSupplier();
+}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
new file mode 100644
index 000000000..569d8b03e
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
@@ -0,0 +1,56 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RandomAccess;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.img.list.ListImg;
+import net.imglib2.util.Util;
+
+import java.util.List;
+import java.util.function.IntFunction;
+import java.util.function.Supplier;
+
+/**
+ * Stores KDTree values as 1D {@code RandomAccessibleInterval<T>}.
+ *
+ * @param <T>
+ * 		the type of values stored in the tree.
+ */
+public class KDTreeValuesImg<T> extends KDTreeValues<T> {
+
+	private final RandomAccessibleInterval<T> values;
+	private final Supplier<IntFunction<T>> valuesSupplier;
+
+	public KDTreeValuesImg(RandomAccessibleInterval<T> values) {
+		super(Util.getTypeFromInterval(values));
+		this.values = values;
+		valuesSupplier = () -> {
+			final RandomAccess<T> ra = values.randomAccess();
+			return ra::setPositionAndGet;
+		};
+	}
+
+	/**
+	 * @return number of values in the tree
+	 */
+	public int size() {
+		return (int) values.dimension(0);
+	}
+
+	/**
+	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
+	 */
+	public RandomAccessibleInterval<T> values() {
+		return values;
+	}
+
+	/**
+	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
+	 * values for a given node indices. If the returned {@code IntFunction<T>}
+	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
+	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
+	 * {@code IntFunction<T>}.
+	 */
+	public Supplier<IntFunction<T>> valuesSupplier() {
+		return valuesSupplier;
+	}
+}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
new file mode 100644
index 000000000..a76ab85a5
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
@@ -0,0 +1,51 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.img.list.ListImg;
+
+import java.util.List;
+import java.util.function.IntFunction;
+import java.util.function.Supplier;
+
+/**
+ * Stores KDTree values as List.
+ *
+ * @param <T>
+ * 		the type of values stored in the tree.
+ */
+public class KDTreeValuesList<T> extends KDTreeValues<T> {
+
+	private final List<T> values;
+	private final Supplier<IntFunction<T>> valuesSupplier;
+
+	public KDTreeValuesList(List<T> values) {
+		super(KDTreeUtils.getType(values));
+		this.values = values;
+		valuesSupplier = () -> values::get;
+	}
+
+	/**
+	 * @return number of values in the tree
+	 */
+	public int size() {
+		return values.size();
+	}
+
+	/**
+	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
+	 */
+	public RandomAccessibleInterval<T> values() {
+		return ListImg.wrap(values, values.size());
+	}
+
+	/**
+	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
+	 * values for a given node indices. If the returned {@code IntFunction<T>}
+	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
+	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
+	 * {@code IntFunction<T>}.
+	 */
+	public Supplier<IntFunction<T>> valuesSupplier() {
+		return valuesSupplier;
+	}
+}

From ea26d26ffba98af25079135d02361972d7ed1f97 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Fri, 11 Aug 2023 15:09:17 -0400
Subject: [PATCH 09/13] Separate actual positions from tree implementation

---
 .../java/net/imglib2/kdtree/KDTreeData.java   |  5 ++
 .../java/net/imglib2/kdtree/KDTreeImpl.java   | 33 +++-----
 .../net/imglib2/kdtree/KDTreePositions.java   | 80 +++++++++++++++++++
 .../net/imglib2/kdtree/KDTreeValuesImg.java   |  3 +
 .../net/imglib2/kdtree/KDTreeValuesList.java  |  3 +
 5 files changed, 102 insertions(+), 22 deletions(-)
 create mode 100644 src/main/java/net/imglib2/kdtree/KDTreePositions.java

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index 7f8218577..463ec804c 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -58,6 +58,7 @@ public enum PositionsLayout
 	private final PositionsLayout layout;
 	private final double[][] positions;
 	private final double[] flatPositions;
+	private final KDTreePositions newPositions;
 
 	private final KDTreeValues<T> values;
 
@@ -71,6 +72,7 @@ public KDTreeData( double[][] positions, List< T > values )
 		layout = PositionsLayout.NESTED;
 		this.positions = positions;
 		flatPositions = null;
+		newPositions = new KDTreePositions.Nested(positions);
 
 		this.values = new KDTreeValuesList<>(values);
 	}
@@ -89,6 +91,7 @@ public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
 		layout = PositionsLayout.NESTED;
 		this.positions = positions;
 		flatPositions = null;
+		newPositions = new KDTreePositions.Nested(positions);
 
 		this.values = new KDTreeValuesImg<>(values);
 	}
@@ -107,6 +110,7 @@ public KDTreeData( double[] positions, List< T > values )
 		layout = FLAT;
 		this.positions = null;
 		flatPositions = positions;
+		newPositions = new KDTreePositions.Flat(positions, numDimensions);
 
 		this.values = new KDTreeValuesList<>(values);
 	}
@@ -125,6 +129,7 @@ public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
 		layout = FLAT;
 		this.positions = null;
 		flatPositions = positions;
+		newPositions = new KDTreePositions.Flat(positions, numDimensions);
 
 		this.values = new KDTreeValuesImg<>(values);
 	}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeImpl.java b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
index 2f4468e77..9c33cb9e5 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
@@ -47,37 +47,24 @@ public abstract class KDTreeImpl
 
 	private final int numPoints;
 
+	// TODO: make this final
+	protected KDTreePositions positions;
+
 	static class Nested extends KDTreeImpl
 	{
-		private final double[][] positions;
-
 		Nested( final double[][] positions )
 		{
 			super( positions.length, positions[ 0 ].length );
-			this.positions = positions;
-		}
-
-		@Override
-		public double getDoublePosition( final int i, final int d )
-		{
-			return positions[ d ][ i ];
+			this.positions = new KDTreePositions.Nested(positions);
 		}
 	}
 
 	static class Flat extends KDTreeImpl
 	{
-		private final double[] positions;
-
 		Flat( final double[] positions, final int numDimensions )
 		{
 			super( numDimensions, positions.length / numDimensions );
-			this.positions = positions;
-		}
-
-		@Override
-		public double getDoublePosition( final int i, final int d )
-		{
-			return positions[ numDimensions * i + d ];
+			this.positions = new KDTreePositions.Flat(positions, numDimensions);
 		}
 	}
 
@@ -158,7 +145,9 @@ public int splitDimension( final int i )
 		return ( 31 - Integer.numberOfLeadingZeros( i + 1 ) ) % numDimensions;
 	}
 
-	public abstract double getDoublePosition( final int i, final int d );
+	public double getDoublePosition(final int i, final int d) {
+		return positions.get(i, d);
+	}
 
 	/**
 	 * Compute the squared distance from node {@code i} to {@code pos}.
@@ -168,7 +157,7 @@ public float squDistance( final int i, final float[] pos )
 		float sum = 0;
 		for ( int d = 0; d < numDimensions; ++d )
 		{
-			final float diff = pos[ d ] - ( float ) getDoublePosition( i, d );
+			final float diff = pos[ d ] - ( float ) positions.get(i, d);
 			sum += diff * diff;
 		}
 		return sum;
@@ -182,7 +171,7 @@ public double squDistance( final int i, final double[] pos )
 		double sum = 0;
 		for ( int d = 0; d < numDimensions; ++d )
 		{
-			final double diff = pos[ d ] - getDoublePosition( i, d );
+			final double diff = pos[ d ] - positions.get(i, d);
 			sum += diff * diff;
 		}
 		return sum;
@@ -196,7 +185,7 @@ public double squDistance( final int i, final RealLocalizable pos )
 		double sum = 0;
 		for ( int d = 0; d < numDimensions; ++d )
 		{
-			final double diff = pos.getDoublePosition( d ) - getDoublePosition( i, d );
+			final double diff = pos.getDoublePosition( d ) - positions.get(i, d);
 			sum += diff * diff;
 		}
 		return sum;
diff --git a/src/main/java/net/imglib2/kdtree/KDTreePositions.java b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
new file mode 100644
index 000000000..86f3f8ef8
--- /dev/null
+++ b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
@@ -0,0 +1,80 @@
+package net.imglib2.kdtree;
+
+import net.imglib2.RealLocalizable;
+
+import static net.imglib2.kdtree.KDTreeUtils.leftChildIndex;
+import static net.imglib2.kdtree.KDTreeUtils.parentIndex;
+import static net.imglib2.kdtree.KDTreeUtils.rightChildIndex;
+
+/**
+ * Stores the positions of the nodes in a KDTree and provides access to them.
+ * <p>
+ * Currently, there are two implementations:
+ * - {@link Nested} stores the positions as a {@code double[][]}, allowing a
+ *   total of 2^31-8 nodes but doesn't keep the positions contiguous in memory.
+ * - {@link Flat} stores the positions as a {@code double[]}, meaning that the
+ *   positions are contiguous in memory but the number of nodes is limited to
+ *   (2^31-8)/numDimensions.
+ */
+public abstract class KDTreePositions {
+
+	protected final int numDimensions;
+
+	protected final int numPoints;
+
+	static class Nested extends KDTreePositions {
+		private final double[][] positions;
+
+		Nested(final double[][] positions) {
+			super(positions.length, positions[0].length);
+			this.positions = positions;
+		}
+
+		@Override
+		public double get(final int i, final int d) {
+			return positions[d][i];
+		}
+	}
+
+	static class Flat extends KDTreePositions {
+		private final double[] positions;
+
+		Flat(final double[] positions, final int numDimensions) {
+			super(numDimensions, positions.length / numDimensions);
+			this.positions = positions;
+		}
+
+		@Override
+		public double get(final int i, final int d) {
+			return positions[numDimensions * i + d];
+		}
+	}
+
+	KDTreePositions(final int numDimensions, final int numPoints) {
+		this.numDimensions = numDimensions;
+		this.numPoints = numPoints;
+	}
+
+	/**
+	 * Get the coordinates of the node {@code i} in dimension {@code d}.
+	 *
+	 * @return the coordinate
+	 */
+	public abstract double get(final int i, final int d);
+
+	public int numDimensions() {
+		return numDimensions;
+	}
+
+	public int size() {
+		return numPoints;
+	}
+
+	public static KDTreePositions create(final double[][] positions) {
+		return new Nested( positions );
+	}
+
+	public static KDTreePositions create(final double[] positions, final int numDimensions) {
+		return new Flat(positions, numDimensions);
+	}
+}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
index 569d8b03e..dddb1be10 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
@@ -32,6 +32,7 @@ public KDTreeValuesImg(RandomAccessibleInterval<T> values) {
 	/**
 	 * @return number of values in the tree
 	 */
+	@Override
 	public int size() {
 		return (int) values.dimension(0);
 	}
@@ -39,6 +40,7 @@ public int size() {
 	/**
 	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
 	 */
+	@Override
 	public RandomAccessibleInterval<T> values() {
 		return values;
 	}
@@ -50,6 +52,7 @@ public RandomAccessibleInterval<T> values() {
 	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
 	 * {@code IntFunction<T>}.
 	 */
+	@Override
 	public Supplier<IntFunction<T>> valuesSupplier() {
 		return valuesSupplier;
 	}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
index a76ab85a5..cb14c6ee2 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
@@ -27,6 +27,7 @@ public KDTreeValuesList(List<T> values) {
 	/**
 	 * @return number of values in the tree
 	 */
+	@Override
 	public int size() {
 		return values.size();
 	}
@@ -34,6 +35,7 @@ public int size() {
 	/**
 	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
 	 */
+	@Override
 	public RandomAccessibleInterval<T> values() {
 		return ListImg.wrap(values, values.size());
 	}
@@ -45,6 +47,7 @@ public RandomAccessibleInterval<T> values() {
 	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
 	 * {@code IntFunction<T>}.
 	 */
+	@Override
 	public Supplier<IntFunction<T>> valuesSupplier() {
 		return valuesSupplier;
 	}

From 1232e56a8259d6671bc4728a42984578334d8987 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Fri, 11 Aug 2023 15:51:28 -0400
Subject: [PATCH 10/13] Store only KDTreePositions in KDTreeData

---
 src/main/java/net/imglib2/KDTree.java         |   9 +-
 .../java/net/imglib2/kdtree/KDTreeData.java   | 130 ++----------------
 .../java/net/imglib2/kdtree/KDTreeImpl.java   |  41 +-----
 .../net/imglib2/kdtree/KDTreePositions.java   |  97 +++++++++++--
 .../net/imglib2/kdtree/KDTreeBenchmark.java   |   2 +-
 .../net/imglib2/kdtree/KDTreeImplTest.java    |   6 +-
 6 files changed, 113 insertions(+), 172 deletions(-)

diff --git a/src/main/java/net/imglib2/KDTree.java b/src/main/java/net/imglib2/KDTree.java
index 6a95e53b8..8aae46b91 100644
--- a/src/main/java/net/imglib2/KDTree.java
+++ b/src/main/java/net/imglib2/KDTree.java
@@ -1,15 +1,12 @@
 package net.imglib2;
 
 import java.util.List;
-import java.util.function.IntFunction;
-import java.util.function.Supplier;
+
 import net.imglib2.converter.AbstractConvertedIterableRealInterval;
 import net.imglib2.converter.AbstractConvertedRealCursor;
 import net.imglib2.kdtree.KDTreeData;
 import net.imglib2.kdtree.KDTreeImpl;
 
-import static net.imglib2.kdtree.KDTreeData.PositionsLayout.FLAT;
-
 public class KDTree< T > implements EuclideanSpace, IterableRealInterval< T >
 {
 	private final KDTreeData< T > treeData;
@@ -132,9 +129,7 @@ public < L extends RealLocalizable > KDTree( final int numPoints, final Iterable
 	public KDTree( final KDTreeData< T > data )
 	{
 		treeData = data;
-		impl = ( data.layout() == FLAT )
-				? KDTreeImpl.create( data.flatPositions(), data.numDimensions() )
-				: KDTreeImpl.create( data.positions() );
+		impl = new KDTreeImpl(treeData.positions);
 	}
 
 	/**
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index 463ec804c..7635de4f9 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -3,62 +3,36 @@
 import java.util.List;
 import java.util.function.IntFunction;
 import java.util.function.Supplier;
-import net.imglib2.FinalRealInterval;
-import net.imglib2.RandomAccess;
 import net.imglib2.RandomAccessibleInterval;
 import net.imglib2.RealInterval;
 import net.imglib2.RealLocalizable;
 import net.imglib2.img.Img;
-import net.imglib2.img.list.ListImg;
 import net.imglib2.type.NativeType;
-import net.imglib2.util.Util;
 
-import static net.imglib2.kdtree.KDTreeData.PositionsLayout.FLAT;
 
 /**
  * Stores the KDTree data, that is, positions and values.
  * <p>
- * Positions are stored in either {@code FLAT} or {@code NESTED} {@link
- * PositionsLayout layout}. With {@code NESTED} layout, positions are stored as
- * a nested {@code double[][]} array where {@code positions[d][i]} is dimension
- * {@code d} of the {@code i}-th point. With {@code FLAT} layout, positions are
- * stored as a flat {@code double[]} array, where {@code positions[d + i*n]} is
- * dimension {@code d} of the {@code i}-th point, with {@code n} the number of
- * dimensions.
+ * Positions are stored as {@link KDTreePositions} as either
+ * {@link KDTreePositions.Flat} or {@link KDTreePositions.Nested}.
  * <p>
- * Values (of type {@code T}) are stored as either a 1D {@code
- * RandomAccessibleInterval<T>}, or a {@code List<T>}. Individual values can be
- * accessed by {@link #valuesSupplier()}{@code .get().apply(i)}. {@code
- * valueSupplier().get()} returns a reusable {@code IntFunction<T>}. Here {@code
- * T} maybe a proxy that is reused in subsequent {@code apply(i)}.
+ * Values (of type {@code T}) are stored as {@link KDTreeValues} either a
+ * 1D {@code RandomAccessibleInterval<T>} in {@link KDTreeValuesImg}, or
+ * a {@code List<T>} in {@link KDTreeValuesList}.
  * <p>
  * {@link #values()} returns all values as a 1D {@code
- * RandomAccessibleInterval<T>}. (If data is stored as {@code List<T>}, it is
- * wrapped into a {@code ListImg}.)
- * <p>
- * {@link #positions()} returns positions in nested {@code double[][]} (which is
- * created if internal storage is {@code FLAT}). {@link #flatPositions()}
- * returns flat {@code double[]} if internal storage is {@code FLAT}, otherwise
- * {@code null}.
+ * RandomAccessibleInterval<T>}. (If data is stored as {@link
+ * KDTreeValuesList}, it is wrapped into a {@code ListImg}.)
  *
  * @param <T>
  * 		the type of values stored in the tree.
  */
 public class KDTreeData< T >
 {
-	public enum PositionsLayout
-	{
-		FLAT,
-		NESTED
-	}
-
 	private final int numDimensions;
 	private final int numPoints;
 
-	private final PositionsLayout layout;
-	private final double[][] positions;
-	private final double[] flatPositions;
-	private final KDTreePositions newPositions;
+	public final KDTreePositions positions;
 
 	private final KDTreeValues<T> values;
 
@@ -68,12 +42,7 @@ public KDTreeData( double[][] positions, List< T > values )
 	{
 		numPoints = values.size();
 		numDimensions = positions.length;
-
-		layout = PositionsLayout.NESTED;
-		this.positions = positions;
-		flatPositions = null;
-		newPositions = new KDTreePositions.Nested(positions);
-
+		this.positions = new KDTreePositions.Nested(positions);
 		this.values = new KDTreeValuesList<>(values);
 	}
 
@@ -87,12 +56,7 @@ public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
 	{
 		numPoints = ( int ) values.dimension( 0 );
 		numDimensions = positions.length;
-
-		layout = PositionsLayout.NESTED;
-		this.positions = positions;
-		flatPositions = null;
-		newPositions = new KDTreePositions.Nested(positions);
-
+		this.positions = new KDTreePositions.Nested(positions);
 		this.values = new KDTreeValuesImg<>(values);
 	}
 
@@ -106,12 +70,7 @@ public KDTreeData( double[] positions, List< T > values )
 	{
 		numPoints = values.size();
 		numDimensions = positions.length / numPoints;
-
-		layout = FLAT;
-		this.positions = null;
-		flatPositions = positions;
-		newPositions = new KDTreePositions.Flat(positions, numDimensions);
-
+		this.positions = new KDTreePositions.Flat(positions, numDimensions);
 		this.values = new KDTreeValuesList<>(values);
 	}
 
@@ -125,12 +84,7 @@ public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
 	{
 		numPoints = ( int ) values.dimension( 0 );
 		numDimensions = positions.length / numPoints;
-
-		layout = FLAT;
-		this.positions = null;
-		flatPositions = positions;
-		newPositions = new KDTreePositions.Flat(positions, numDimensions);
-
+		this.positions = new KDTreePositions.Flat(positions, numDimensions);
 		this.values = new KDTreeValuesImg<>(values);
 	}
 
@@ -146,19 +100,9 @@ public T type()
 		return values.type();
 	}
 
-	/**
-	 * Get positions of points in the tree as a nested {@code double[][]} array
-	 * where {@code positions[d][i]} is dimension {@code d} of the {@code i}-th
-	 * point.
-	 * <p>
-	 * For serialisation and usage by the tree.
-	 * <p>
-	 * Internal storage may be flattened into single {@code double[]} array. In
-	 * this case, the nested {@code double[][]} array is created here.
-	 */
-	public double[][] positions()
+	public KDTreePositions positions()
 	{
-		return layout == FLAT ? KDTreeUtils.unflatten( flatPositions, numDimensions ) : positions;
+		return positions;
 	}
 
 	/**
@@ -183,37 +127,6 @@ public Supplier< IntFunction< T > > valuesSupplier()
 		return values.valuesSupplier();
 	}
 
-	/**
-	 * Get positions of points in the tree as a flat {@code double[]} where
-	 * {@code positions[d + i*n]} is dimension {@code d} of the {@code i}-th
-	 * point, with {@code n} the number of dimensions.
-	 * <p>
-	 * For serialisation and usage by the tree.
-	 * <p>
-	 * Internal storage may be a {@code NESTED} {@code double[][]} array. In
-	 * this case, {@code flatPositions()} returns {@code null}.
-	 */
-	public double[] flatPositions()
-	{
-		return flatPositions;
-	}
-
-	/**
-	 * Get the internal layout of positions.
-	 * <p>
-	 * Positions are stored in either {@code FLAT} or {@code NESTED} {@link
-	 * PositionsLayout layout}. With {@code NESTED} layout, positions are stored
-	 * as a nested {@code double[][]} array where {@code positions[d][i]} is
-	 * dimension {@code d} of the {@code i}-th point. With {@code FLAT} layout,
-	 * positions are stored as a flat {@code double[]} array, where {@code
-	 * positions[d + i*n]} is dimension {@code d} of the {@code i}-th point,
-	 * with {@code n} the number of dimensions.
-	 */
-	public PositionsLayout layout()
-	{
-		return layout;
-	}
-
 	/**
 	 * @return dimensionality of points in the tree
 	 */
@@ -232,20 +145,7 @@ public int size()
 
 	public RealInterval boundingBox()
 	{
-		if ( boundingBox == null )
-			boundingBox = createBoundingBox();
-		return boundingBox;
-	}
-
-	private RealInterval createBoundingBox()
-	{
-		final double[] min = new double[ numDimensions ];
-		final double[] max = new double[ numDimensions ];
-		if ( layout == FLAT )
-			KDTreeUtils.computeMinMax( flatPositions, min, max );
-		else
-			KDTreeUtils.computeMinMax( positions, min, max );
-		return FinalRealInterval.wrap( min, max );
+		return (boundingBox != null) ? boundingBox : positions.boundingBox();
 	}
 
 	/**
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeImpl.java b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
index 9c33cb9e5..5732a1703 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeImpl.java
@@ -41,37 +41,18 @@
  * without dependent reads. And iff the calculated child index is less
  * than the number of nodes, the child exists.
  */
-public abstract class KDTreeImpl
+public class KDTreeImpl
 {
 	final int numDimensions;
 
 	private final int numPoints;
 
-	// TODO: make this final
-	protected KDTreePositions positions;
+	protected final KDTreePositions positions;
 
-	static class Nested extends KDTreeImpl
-	{
-		Nested( final double[][] positions )
-		{
-			super( positions.length, positions[ 0 ].length );
-			this.positions = new KDTreePositions.Nested(positions);
-		}
-	}
-
-	static class Flat extends KDTreeImpl
-	{
-		Flat( final double[] positions, final int numDimensions )
-		{
-			super( numDimensions, positions.length / numDimensions );
-			this.positions = new KDTreePositions.Flat(positions, numDimensions);
-		}
-	}
-
-	KDTreeImpl( final int numDimensions, final int numPoints )
-	{
-		this.numDimensions = numDimensions;
-		this.numPoints = numPoints;
+	public KDTreeImpl(final KDTreePositions positions) {
+		this.numDimensions = positions.numDimensions;
+		this.numPoints = positions.numPoints;
+		this.positions = positions;
 	}
 
 	/**
@@ -205,14 +186,4 @@ public int depth()
 	{
 		return 32 - Integer.numberOfLeadingZeros( numPoints );
 	}
-
-	public static KDTreeImpl create( final double[][] positions )
-	{
-		return new Nested( positions );
-	}
-
-	public static KDTreeImpl create( final double[] positions, final int numDimensions )
-	{
-		return new Flat( positions, numDimensions );
-	}
 }
diff --git a/src/main/java/net/imglib2/kdtree/KDTreePositions.java b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
index 86f3f8ef8..1e9e3911c 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreePositions.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
@@ -1,26 +1,32 @@
 package net.imglib2.kdtree;
 
-import net.imglib2.RealLocalizable;
+import net.imglib2.FinalRealInterval;
+import net.imglib2.RealInterval;
 
-import static net.imglib2.kdtree.KDTreeUtils.leftChildIndex;
-import static net.imglib2.kdtree.KDTreeUtils.parentIndex;
-import static net.imglib2.kdtree.KDTreeUtils.rightChildIndex;
 
 /**
  * Stores the positions of the nodes in a KDTree and provides access to them.
  * <p>
  * Currently, there are two implementations:
- * - {@link Nested} stores the positions as a {@code double[][]}, allowing a
- *   total of 2^31-8 nodes but doesn't keep the positions contiguous in memory.
- * - {@link Flat} stores the positions as a {@code double[]}, meaning that the
- *   positions are contiguous in memory but the number of nodes is limited to
- *   (2^31-8)/numDimensions.
+ * - {@link Nested} stores the positions as a {@code double[][]} where {@code
+ *   positions[d][i]} is dimension {@code d} of the {@code i}-th point. This
+ *   allows for a total of 2^31-8 nodes but doesn't keep the positions
+ *   contiguous in memory.
+ * - {@link Flat} stores the positions as a {@code double[]}where {@code
+ *   positions[d + i*n]} is dimension {@code d} of the {@code i}-th point,
+ *   with {@code n} the number of dimensions. This means that the positions
+ *   are contiguous in memory but the number of nodes is limited to (2^31-8)/n.
+ *  <p>
+ * {@link #getNestedPositions()} returns positions in nested {@code double[][]}
+ * (which is created if class is {@link Flat}). {@link #getFlatPositions()}
+ * returns flat {@code double[]} if class is {@link Flat}, otherwise {@code null}.
  */
 public abstract class KDTreePositions {
 
 	protected final int numDimensions;
 
 	protected final int numPoints;
+	private volatile RealInterval boundingBox;
 
 	static class Nested extends KDTreePositions {
 		private final double[][] positions;
@@ -30,10 +36,28 @@ static class Nested extends KDTreePositions {
 			this.positions = positions;
 		}
 
-		@Override
-		public double get(final int i, final int d) {
+		@Override public double get(final int i, final int d) {
 			return positions[d][i];
 		}
+
+		@Override
+		public double[] getFlatPositions() {
+			// positions in this case might be too large to fit in a single array
+			return null;
+		}
+
+		@Override
+		public double[][] getNestedPositions() {
+			return positions;
+		}
+
+		@Override
+		protected RealInterval createBoundingBox() {
+			final double[] min = new double[numDimensions];
+			final double[] max = new double[numDimensions];
+			KDTreeUtils.computeMinMax(positions, min, max);
+			return FinalRealInterval.wrap(min, max);
+		}
 	}
 
 	static class Flat extends KDTreePositions {
@@ -48,6 +72,24 @@ static class Flat extends KDTreePositions {
 		public double get(final int i, final int d) {
 			return positions[numDimensions * i + d];
 		}
+
+		@Override
+		public double[] getFlatPositions() {
+			return positions;
+		}
+
+		@Override
+		public double[][] getNestedPositions() {
+			return KDTreeUtils.unflatten(positions, numDimensions);
+		}
+
+		@Override
+		protected RealInterval createBoundingBox() {
+			final double[] min = new double[numDimensions];
+			final double[] max = new double[numDimensions];
+			KDTreeUtils.computeMinMax(positions, min, max);
+			return FinalRealInterval.wrap(min, max);
+		}
 	}
 
 	KDTreePositions(final int numDimensions, final int numPoints) {
@@ -62,6 +104,33 @@ public double get(final int i, final int d) {
 	 */
 	public abstract double get(final int i, final int d);
 
+	/**
+	 * Get positions of points in the tree as a flat {@code double[]} array
+	 * where {@code positions[d + i*n]} is dimension {@code d} of the {@code i}-th
+	 * point.
+	 * <p>
+	 * For serialisation and usage by the tree.
+	 * <p>
+	 * Internal storage may be nested in a {@code double[][]} array. In
+	 * this case, there may be too many points to fit in a single {@code
+	 * double[]} array, so {@code null} is returned.
+	 */
+	public abstract double[] getFlatPositions();
+
+	/**
+	 * Get positions of points in the tree as a nested {@code double[][]} array
+	 * where {@code positions[d][i]} is dimension {@code d} of the {@code i}-th
+	 * point.
+	 * <p>
+	 * For serialisation and usage by the tree.
+	 * <p>
+	 * Internal storage may be flattened into single {@code double[]} array. In
+	 * this case, the nested {@code double[][]} array is created here.
+	 */
+	public abstract double[][] getNestedPositions();
+
+	protected abstract RealInterval createBoundingBox();
+
 	public int numDimensions() {
 		return numDimensions;
 	}
@@ -70,6 +139,12 @@ public int size() {
 		return numPoints;
 	}
 
+	public RealInterval boundingBox() {
+		if (boundingBox == null)
+			boundingBox = createBoundingBox();
+		return boundingBox;
+	}
+
 	public static KDTreePositions create(final double[][] positions) {
 		return new Nested( positions );
 	}
diff --git a/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java b/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
index 79c4e8fdd..d7eef7c88 100644
--- a/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
+++ b/src/test/java/net/imglib2/kdtree/KDTreeBenchmark.java
@@ -61,7 +61,7 @@ public void spoil() {
 		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
 		final int[] tree = KDTreeUtils.makeTree( points );
 		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
-		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final KDTreeImpl impl = new KDTreeImpl(new KDTreePositions.Nested(treePoints));
 		final NearestNeighborSearchImpl search = new NearestNeighborSearchImpl( impl );
 		for ( RealPoint testVertex : testVertices )
 			search.search( testVertex );
diff --git a/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java b/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java
index d3095f98f..e683320a8 100644
--- a/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java
+++ b/src/test/java/net/imglib2/kdtree/KDTreeImplTest.java
@@ -52,7 +52,7 @@ public void testNearestNeighborSearch()
 		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
 		final int[] tree = KDTreeUtils.makeTree( points );
 		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
-		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final KDTreeImpl impl = new KDTreeImpl(new KDTreePositions.Nested(treePoints));
 		final NearestNeighborSearchImpl search = new NearestNeighborSearchImpl( impl );
 
 		for ( RealPoint testVertex : testVertices )
@@ -75,7 +75,7 @@ public void testKNearestNeighborSearch()
 		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
 		final int[] tree = KDTreeUtils.makeTree( points );
 		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
-		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final KDTreeImpl impl = new KDTreeImpl(new KDTreePositions.Nested(treePoints));
 		final KNearestNeighborSearchImpl search = new KNearestNeighborSearchImpl( impl, k );
 
 		for ( RealPoint testVertex : testVertices )
@@ -99,7 +99,7 @@ public void testRadiusNeighborSearch()
 		final double[][] points = KDTreeUtils.initPositions( n, numDataVertices, dataVertices );
 		final int[] tree = KDTreeUtils.makeTree( points );
 		final double[][] treePoints = KDTreeUtils.reorder( points, tree );
-		final KDTreeImpl impl = new KDTreeImpl.Nested( treePoints );
+		final KDTreeImpl impl = new KDTreeImpl(new KDTreePositions.Nested(treePoints));
 		final RadiusNeighborSearchImpl search = new RadiusNeighborSearchImpl( impl );
 
 		for ( RealPoint testVertex : testVertices )

From ac6cf63f96103e6874921725a4415f2e0bae3830 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Fri, 11 Aug 2023 17:07:00 -0400
Subject: [PATCH 11/13] Add some convenience functions used in serialization

---
 src/main/java/net/imglib2/kdtree/KDTreeData.java | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index 7635de4f9..e874aa83d 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -105,6 +105,18 @@ public KDTreePositions positions()
 		return positions;
 	}
 
+	public double[] getFlatPositions() {
+		return positions.getFlatPositions();
+	}
+
+	public double[][] getNestedPositions() {
+		return positions.getNestedPositions();
+	}
+
+	public boolean positionsIsFlatArray() {
+		return (positions instanceof KDTreePositions.Flat);
+	}
+
 	/**
 	 * Get the values as a 1D {@code RandomAccessibleInterval}, for
 	 * serialization. (If the underlying storage is a {@code List<T>}, it will

From 2368e5c7a35fc0dfc5e7bd4d6f3c1edcdb445998 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Fri, 11 Aug 2023 20:19:09 -0400
Subject: [PATCH 12/13] Consolidate KDTreeValues

---
 .../java/net/imglib2/kdtree/KDTreeData.java   | 14 ++--
 .../java/net/imglib2/kdtree/KDTreeValues.java | 81 +++++++++++++++++++
 .../net/imglib2/kdtree/KDTreeValuesImg.java   | 59 --------------
 .../net/imglib2/kdtree/KDTreeValuesList.java  | 54 -------------
 4 files changed, 88 insertions(+), 120 deletions(-)
 delete mode 100644 src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
 delete mode 100644 src/main/java/net/imglib2/kdtree/KDTreeValuesList.java

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index e874aa83d..80bb52b99 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -17,12 +17,12 @@
  * {@link KDTreePositions.Flat} or {@link KDTreePositions.Nested}.
  * <p>
  * Values (of type {@code T}) are stored as {@link KDTreeValues} either a
- * 1D {@code RandomAccessibleInterval<T>} in {@link KDTreeValuesImg}, or
- * a {@code List<T>} in {@link KDTreeValuesList}.
+ * 1D {@code RandomAccessibleInterval<T>} in {@link KDTreeValues.ImgValues}, or
+ * a {@code List<T>} in {@link KDTreeValues.ListValues}.
  * <p>
  * {@link #values()} returns all values as a 1D {@code
  * RandomAccessibleInterval<T>}. (If data is stored as {@link
- * KDTreeValuesList}, it is wrapped into a {@code ListImg}.)
+ * KDTreeValues.ListValues}, it is wrapped into a {@code ListImg}.)
  *
  * @param <T>
  * 		the type of values stored in the tree.
@@ -43,7 +43,7 @@ public KDTreeData( double[][] positions, List< T > values )
 		numPoints = values.size();
 		numDimensions = positions.length;
 		this.positions = new KDTreePositions.Nested(positions);
-		this.values = new KDTreeValuesList<>(values);
+		this.values = new KDTreeValues.ListValues<>(values);
 	}
 
 	public KDTreeData( double[][] positions, List< T > values, RealInterval boundingBox )
@@ -57,7 +57,7 @@ public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
 		numPoints = ( int ) values.dimension( 0 );
 		numDimensions = positions.length;
 		this.positions = new KDTreePositions.Nested(positions);
-		this.values = new KDTreeValuesImg<>(values);
+		this.values = new KDTreeValues.ImgValues<>(values);
 	}
 
 	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
@@ -71,7 +71,7 @@ public KDTreeData( double[] positions, List< T > values )
 		numPoints = values.size();
 		numDimensions = positions.length / numPoints;
 		this.positions = new KDTreePositions.Flat(positions, numDimensions);
-		this.values = new KDTreeValuesList<>(values);
+		this.values = new KDTreeValues.ListValues<>(values);
 	}
 
 	public KDTreeData( double[] positions, List< T > values, RealInterval boundingBox )
@@ -85,7 +85,7 @@ public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
 		numPoints = ( int ) values.dimension( 0 );
 		numDimensions = positions.length / numPoints;
 		this.positions = new KDTreePositions.Flat(positions, numDimensions);
-		this.values = new KDTreeValuesImg<>(values);
+		this.values = new KDTreeValues.ImgValues<>(values);
 	}
 
 	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValues.java b/src/main/java/net/imglib2/kdtree/KDTreeValues.java
index ab66189ac..74c476657 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeValues.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeValues.java
@@ -1,7 +1,11 @@
 package net.imglib2.kdtree;
 
+import net.imglib2.RandomAccess;
 import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.img.list.ListImg;
+import net.imglib2.util.Util;
 
+import java.util.List;
 import java.util.function.IntFunction;
 import java.util.function.Supplier;
 
@@ -24,6 +28,75 @@
 abstract public class KDTreeValues<T> {
 	protected final T type;
 
+	/**
+	 * Stores KDTree values as 1D {@code RandomAccessibleInterval<T>}.
+	 *
+	 * @param <T>
+	 * 		the type of values stored in the tree.
+	 */
+	public static class ImgValues<T> extends KDTreeValues<T> {
+
+		private final RandomAccessibleInterval<T> values;
+		private final Supplier<IntFunction<T>> valuesSupplier;
+
+		public ImgValues(RandomAccessibleInterval<T> values) {
+			super(Util.getTypeFromInterval(values));
+			this.values = values;
+			valuesSupplier = () -> {
+				final RandomAccess<T> ra = values.randomAccess();
+				return ra::setPositionAndGet;
+			};
+		}
+
+		@Override
+		public int size() {
+			return (int) values.dimension(0);
+		}
+
+		@Override
+		public RandomAccessibleInterval<T> values() {
+			return values;
+		}
+
+		@Override
+		public Supplier<IntFunction<T>> valuesSupplier() {
+			return valuesSupplier;
+		}
+	}
+
+	/**
+	 * Stores KDTree values as List.
+	 *
+	 * @param <T>
+	 * 		the type of values stored in the tree.
+	 */
+	public static class ListValues<T> extends KDTreeValues<T> {
+
+		private final List<T> values;
+		private final Supplier<IntFunction<T>> valuesSupplier;
+
+		public ListValues(List<T> values) {
+			super(KDTreeUtils.getType(values));
+			this.values = values;
+			valuesSupplier = () -> values::get;
+		}
+
+		@Override
+		public int size() {
+			return values.size();
+		}
+
+		@Override
+		public RandomAccessibleInterval<T> values() {
+			return ListImg.wrap(values, values.size());
+		}
+
+		@Override
+		public Supplier<IntFunction<T>> valuesSupplier() {
+			return valuesSupplier;
+		}
+	}
+
 	public KDTreeValues(T type) {
 		this.type = type;
 	}
@@ -50,4 +123,12 @@ public T type() {
 	 * {@code IntFunction<T>}.
 	 */
 	abstract public Supplier<IntFunction<T>> valuesSupplier();
+
+	public static <S> KDTreeValues<S> create(final RandomAccessibleInterval<S> values) {
+		return new KDTreeValues.ImgValues<>(values);
+	}
+
+	public static <S> KDTreeValues<S> create(final List<S> values) {
+		return new KDTreeValues.ListValues<>(values);
+	}
 }
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
deleted file mode 100644
index dddb1be10..000000000
--- a/src/main/java/net/imglib2/kdtree/KDTreeValuesImg.java
+++ /dev/null
@@ -1,59 +0,0 @@
-package net.imglib2.kdtree;
-
-import net.imglib2.RandomAccess;
-import net.imglib2.RandomAccessibleInterval;
-import net.imglib2.img.list.ListImg;
-import net.imglib2.util.Util;
-
-import java.util.List;
-import java.util.function.IntFunction;
-import java.util.function.Supplier;
-
-/**
- * Stores KDTree values as 1D {@code RandomAccessibleInterval<T>}.
- *
- * @param <T>
- * 		the type of values stored in the tree.
- */
-public class KDTreeValuesImg<T> extends KDTreeValues<T> {
-
-	private final RandomAccessibleInterval<T> values;
-	private final Supplier<IntFunction<T>> valuesSupplier;
-
-	public KDTreeValuesImg(RandomAccessibleInterval<T> values) {
-		super(Util.getTypeFromInterval(values));
-		this.values = values;
-		valuesSupplier = () -> {
-			final RandomAccess<T> ra = values.randomAccess();
-			return ra::setPositionAndGet;
-		};
-	}
-
-	/**
-	 * @return number of values in the tree
-	 */
-	@Override
-	public int size() {
-		return (int) values.dimension(0);
-	}
-
-	/**
-	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
-	 */
-	@Override
-	public RandomAccessibleInterval<T> values() {
-		return values;
-	}
-
-	/**
-	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
-	 * values for a given node indices. If the returned {@code IntFunction<T>}
-	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
-	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
-	 * {@code IntFunction<T>}.
-	 */
-	@Override
-	public Supplier<IntFunction<T>> valuesSupplier() {
-		return valuesSupplier;
-	}
-}
diff --git a/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java b/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
deleted file mode 100644
index cb14c6ee2..000000000
--- a/src/main/java/net/imglib2/kdtree/KDTreeValuesList.java
+++ /dev/null
@@ -1,54 +0,0 @@
-package net.imglib2.kdtree;
-
-import net.imglib2.RandomAccessibleInterval;
-import net.imglib2.img.list.ListImg;
-
-import java.util.List;
-import java.util.function.IntFunction;
-import java.util.function.Supplier;
-
-/**
- * Stores KDTree values as List.
- *
- * @param <T>
- * 		the type of values stored in the tree.
- */
-public class KDTreeValuesList<T> extends KDTreeValues<T> {
-
-	private final List<T> values;
-	private final Supplier<IntFunction<T>> valuesSupplier;
-
-	public KDTreeValuesList(List<T> values) {
-		super(KDTreeUtils.getType(values));
-		this.values = values;
-		valuesSupplier = () -> values::get;
-	}
-
-	/**
-	 * @return number of values in the tree
-	 */
-	@Override
-	public int size() {
-		return values.size();
-	}
-
-	/**
-	 * Get the values as a 1D {@code RandomAccessibleInterval}, for serialization.
-	 */
-	@Override
-	public RandomAccessibleInterval<T> values() {
-		return ListImg.wrap(values, values.size());
-	}
-
-	/**
-	 * Get a {@code Supplier} that return {@code IntFunction<T>} to provide
-	 * values for a given node indices. If the returned {@code IntFunction<T>}
-	 * is stateful ({@code T} maybe a proxy that is reused in subsequent {@code
-	 * apply(i)}} every {@link Supplier#get()} creates a new instance of the
-	 * {@code IntFunction<T>}.
-	 */
-	@Override
-	public Supplier<IntFunction<T>> valuesSupplier() {
-		return valuesSupplier;
-	}
-}

From afb0ff934d47bb6f0c93d13fad293191f87bde43 Mon Sep 17 00:00:00 2001
From: Michael Innerberger <innerbergerm@hhmi.org>
Date: Sat, 12 Aug 2023 15:59:56 -0400
Subject: [PATCH 13/13] Decouple creation logic of positions and values

---
 .../java/net/imglib2/kdtree/KDTreeData.java   | 94 ++++---------------
 .../net/imglib2/kdtree/KDTreePositions.java   |  4 +
 2 files changed, 23 insertions(+), 75 deletions(-)

diff --git a/src/main/java/net/imglib2/kdtree/KDTreeData.java b/src/main/java/net/imglib2/kdtree/KDTreeData.java
index 80bb52b99..db84eccdd 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreeData.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreeData.java
@@ -1,6 +1,5 @@
 package net.imglib2.kdtree;
 
-import java.util.List;
 import java.util.function.IntFunction;
 import java.util.function.Supplier;
 import net.imglib2.RandomAccessibleInterval;
@@ -38,59 +37,15 @@ public class KDTreeData< T >
 
 	private volatile RealInterval boundingBox;
 
-	public KDTreeData( double[][] positions, List< T > values )
-	{
-		numPoints = values.size();
-		numDimensions = positions.length;
-		this.positions = new KDTreePositions.Nested(positions);
-		this.values = new KDTreeValues.ListValues<>(values);
-	}
-
-	public KDTreeData( double[][] positions, List< T > values, RealInterval boundingBox )
-	{
-		this( positions, values );
-		this.boundingBox = boundingBox;
-	}
-
-	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values )
-	{
-		numPoints = ( int ) values.dimension( 0 );
-		numDimensions = positions.length;
-		this.positions = new KDTreePositions.Nested(positions);
-		this.values = new KDTreeValues.ImgValues<>(values);
-	}
-
-	public KDTreeData( double[][] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
-	{
-		this( positions, values );
-		this.boundingBox = boundingBox;
-	}
-
-	public KDTreeData( double[] positions, List< T > values )
-	{
-		numPoints = values.size();
-		numDimensions = positions.length / numPoints;
-		this.positions = new KDTreePositions.Flat(positions, numDimensions);
-		this.values = new KDTreeValues.ListValues<>(values);
-	}
-
-	public KDTreeData( double[] positions, List< T > values, RealInterval boundingBox )
-	{
-		this( positions, values );
-		this.boundingBox = boundingBox;
-	}
-
-	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values )
-	{
-		numPoints = ( int ) values.dimension( 0 );
-		numDimensions = positions.length / numPoints;
-		this.positions = new KDTreePositions.Flat(positions, numDimensions);
-		this.values = new KDTreeValues.ImgValues<>(values);
+	public KDTreeData(KDTreePositions positions, KDTreeValues<T> values) {
+		numPoints = positions.numPoints();
+		numDimensions = positions.numDimensions();
+		this.positions = positions;
+		this.values = values;
 	}
 
-	public KDTreeData( double[] positions, RandomAccessibleInterval< T > values, RealInterval boundingBox )
-	{
-		this( positions, values );
+	public KDTreeData(KDTreePositions positions, KDTreeValues<T> values, RealInterval boundingBox) {
+		this(positions, values);
 		this.boundingBox = boundingBox;
 	}
 
@@ -100,11 +55,6 @@ public T type()
 		return values.type();
 	}
 
-	public KDTreePositions positions()
-	{
-		return positions;
-	}
-
 	public double[] getFlatPositions() {
 		return positions.getFlatPositions();
 	}
@@ -188,23 +138,17 @@ public static < L extends RealLocalizable, T > KDTreeData< T > create(
 		final int[] tree = KDTreeUtils.makeTree( points );
 		final int[] invtree = KDTreeUtils.invert( tree );
 
-		final boolean useFlatLayout = ( long ) numDimensions * numPoints <= KDTreeUtils.MAX_ARRAY_SIZE;
-		if ( storeValuesAsNativeImg && KDTreeUtils.getType( values ) instanceof NativeType )
-		{
-			@SuppressWarnings( "unchecked" )
-			final Img< T > treeValues = ( Img< T > ) KDTreeUtils.orderValuesImg( invtree, ( Iterable ) values );
-			if ( useFlatLayout )
-				return new KDTreeData<>(KDTreeUtils.reorderToFlatLayout( points, tree ), treeValues);
-			else
-				return new KDTreeData<>(KDTreeUtils.reorder( points, tree ), treeValues);
-		}
-		else
-		{
-			final List< T > treeValues = KDTreeUtils.orderValuesList( invtree, values );
-			if ( useFlatLayout )
-				return new KDTreeData<>(KDTreeUtils.reorderToFlatLayout( points, tree ), treeValues);
-			else
-				return new KDTreeData<>(KDTreeUtils.reorder( points, tree ), treeValues);
-		}
+		final boolean storeAsImg = (storeValuesAsNativeImg && KDTreeUtils.getType(values) instanceof NativeType);
+		@SuppressWarnings("unchecked")
+		final KDTreeValues<T> treeValues = (storeAsImg)
+				? new KDTreeValues.ImgValues<>((Img<T>)  KDTreeUtils.orderValuesImg(invtree, (Iterable) values))
+				: new KDTreeValues.ListValues<>(KDTreeUtils.orderValuesList(invtree, values));
+
+		final boolean useFlatLayout = (long) numDimensions * numPoints <= KDTreeUtils.MAX_ARRAY_SIZE;
+		final KDTreePositions treePositions = (useFlatLayout)
+				? new KDTreePositions.Flat(KDTreeUtils.reorderToFlatLayout(points, tree), numDimensions)
+				: new KDTreePositions.Nested(KDTreeUtils.reorder(points, tree));
+
+		return new KDTreeData<>(treePositions, treeValues);
 	}
 }
diff --git a/src/main/java/net/imglib2/kdtree/KDTreePositions.java b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
index 1e9e3911c..8d8acc6db 100644
--- a/src/main/java/net/imglib2/kdtree/KDTreePositions.java
+++ b/src/main/java/net/imglib2/kdtree/KDTreePositions.java
@@ -135,6 +135,10 @@ public int numDimensions() {
 		return numDimensions;
 	}
 
+	public int numPoints() {
+		return numPoints;
+	}
+
 	public int size() {
 		return numPoints;
 	}