Connected components (#54)

* Add connected component analysis

 - `ConnectedComponentAlaysis` implementation of connected component analysis
 - `UnionFind` implementation of union find (required for cca)
 - `ConnectedComponentAnalaysisTest` tests for 2D and 3D

* Switch to imglib2 release version.

Thanks to @tpietzsch for a prompt release!

* Remove explicit imglib2 version

Version specified by pom-scijava is the one we should use

* Fix javadoc error

* Improve connected components
 - union find interface to allow for sparse labels
 - add interface for caller to specify id for each voxel (instead of using IntervalIndexer)

Avoid use of streams

Add IntArrayUnionFind

Add union find for sparse labels with TLongMap as store

Add test with offset

Expose mapping from set root to id to caller and add IdFromIntervalIndexerWithInterval

Restructure CCA

* Format connected component analaysis sources

Move test into appropriate package

* Fix JavaDoc

* Add note that mask and labeling should have same mind and max to JavaDoc

Author: hanslovsky

src/main/java/net/imglib2/algorithm/labeling/ConnectedComponentAnalysis.java

@@ -0,0 +1,319 @@
+/*
+ * #%L
+ * ImgLib2: a general-purpose, multidimensional image processing library.
+ * %%
+ * Copyright (C) 2009 - 2016 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.algorithm.labeling;
+
+import java.util.function.LongFunction;
+import java.util.function.LongUnaryOperator;
+import java.util.function.ToLongBiFunction;
+
+import gnu.trove.map.hash.TLongLongHashMap;
+import net.imglib2.Cursor;
+import net.imglib2.FinalInterval;
+import net.imglib2.Interval;
+import net.imglib2.Localizable;
+import net.imglib2.RandomAccessible;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.algorithm.neighborhood.DiamondShape;
+import net.imglib2.algorithm.neighborhood.Neighborhood;
+import net.imglib2.algorithm.neighborhood.RectangleNeighborhood;
+import net.imglib2.algorithm.neighborhood.Shape;
+import net.imglib2.algorithm.util.unionfind.IntArrayRankedUnionFind;
+import net.imglib2.algorithm.util.unionfind.UnionFind;
+import net.imglib2.type.BooleanType;
+import net.imglib2.type.numeric.IntegerType;
+import net.imglib2.util.IntervalIndexer;
+import net.imglib2.util.Intervals;
+import net.imglib2.view.Views;
+
+/**
+ *
+ * @author Philipp Hanslovsky
+ *
+ */
+public class ConnectedComponentAnalysis
+{
+
+	private static class StartAtOneIdForNextSet implements LongUnaryOperator
+	{
+
+		private final TLongLongHashMap setMappings = new TLongLongHashMap();
+
+		@Override
+		public long applyAsLong( final long root )
+		{
+
+			if ( !setMappings.containsKey( root ) )
+			{
+				setMappings.put( root, setMappings.size() + 1 );
+			}
+			return setMappings.get( root );
+		}
+
+	}
+
+	private static final class IdFromIntervalIndexerWithInterval< T > implements ToLongBiFunction< Localizable, T >
+	{
+
+		private final Interval interval;
+
+		private IdFromIntervalIndexerWithInterval( final Interval interval )
+		{
+			this.interval = interval;
+		}
+
+		@Override
+		public long applyAsLong( final Localizable l, final T t )
+		{
+			return IntervalIndexer.positionToIndexForInterval( l, interval );
+		}
+	}
+
+	public static < T > ToLongBiFunction< Localizable, T > idFromIntervalIndexer( final Interval interval )
+	{
+		return new IdFromIntervalIndexerWithInterval<>( interval );
+	}
+
+	/**
+	 *
+	 * Implementation of connected component analysis that uses
+	 * {@link IntArrayRankedUnionFind} to find sets of pixels that are connected
+	 * with respect to a 4-neighborhood ({@link DiamondShape}) or the
+	 * generalization for higher dimenions over a binary mask. {@code mask} and
+	 * {@code labeling} are expected to have equal min and max.
+	 *
+	 * @param mask
+	 *            Boolean mask to distinguish foreground ({@code true}) from
+	 *            background ({@code false}).
+	 * @param labeling
+	 *            Output parameter to store labeling: background pixels are
+	 *            labeled zero, foreground pixels are greater than zero: 1, 2,
+	 *            ..., N. Note that initially all pixels are expected to be zero
+	 *            as background values will not be written.
+	 */
+	public static < B extends BooleanType< B >, L extends IntegerType< L > > void connectedComponents(
+			final RandomAccessibleInterval< B > mask,
+			final RandomAccessibleInterval< L > labeling )
+	{
+		connectedComponents( mask, labeling, new DiamondShape( 1 ) );
+	}
+
+	/**
+	 *
+	 * Implementation of connected component analysis that uses
+	 * {@link IntArrayRankedUnionFind} to find sets of pixels that are connected
+	 * with respect to a neighborhood ({@code shape}) over a binary mask. {@code mask}
+	 * and {@code labeling} are expected to have equal min and max.
+	 *
+	 * @param mask
+	 *            Boolean mask to distinguish foreground ({@code true}) from
+	 *            background ({@code false}).
+	 * @param labeling
+	 *            Output parameter to store labeling: background pixels are
+	 *            labeled zero, foreground pixels are greater than zero: 1, 2,
+	 *            ..., N. Note that initially all pixels are expected to be zero
+	 *            as background values will not be written.
+	 * @param shape
+	 *            Connectivity of connected components, e.g. 4-neighborhood
+	 *            ({@link DiamondShape}), 8-neighborhood
+	 *            ({@link RectangleNeighborhood}) and their generalisations for
+	 *            higher dimensions.
+	 */
+	public static < B extends BooleanType< B >, L extends IntegerType< L > > void connectedComponents(
+			final RandomAccessibleInterval< B > mask,
+			final RandomAccessibleInterval< L > labeling,
+			final Shape shape )
+	{
+		assert Intervals.numElements( labeling ) < Integer.MAX_VALUE: "Cannot Image Using array union find.";
+		connectedComponents(
+				mask,
+				labeling,
+				shape,
+				n -> new IntArrayRankedUnionFind( ( int ) n ),
+				idFromIntervalIndexer( labeling ),
+				new StartAtOneIdForNextSet() );
+	}
+
+	/**
+	 *
+	 * Implementation of connected component analysis that uses
+	 * {@link UnionFind} to find sets of pixels that are connected with respect
+	 * to a neighborhood ({@code shape}) over a binary mask. {@code mask} and
+	 * {@code labeling} are expected to have equal min and max.
+	 *
+	 * @param mask
+	 *            Boolean mask to distinguish foreground ({@code true}) from
+	 *            background ({@code false}).
+	 * @param labeling
+	 *            Output parameter to store labeling: background pixels are
+	 *            labeled zero, foreground pixels are greater than zero: 1, 2,
+	 *            ..., N. Note that this is expected to be zero as background
+	 *            values will not be written.
+	 * @param shape
+	 *            Connectivity of connected components, e.g. 4-neighborhood
+	 *            ({@link DiamondShape}), 8-neighborhood
+	 *            ({@link RectangleNeighborhood}) and their generalisations for
+	 *            higher dimensions.
+	 * @param unionFindFactory
+	 *            Creates appropriate {@link UnionFind} data structure for size
+	 *            of {@code labeling}, e.g. {@link IntArrayRankedUnionFind} of
+	 *            appropriate size.
+	 * @param idForPixel
+	 *            Create id from pixel location and value. Multiple calls with
+	 *            the same argument should always return the same result.
+	 * @param idForSet
+	 *            Create id for a set from the root id of a set. Multiple calls
+	 *            with the same argument should always return the same result.
+	 */
+	public static < B extends BooleanType< B >, L extends IntegerType< L > > void connectedComponents(
+			final RandomAccessibleInterval< B > mask,
+			final RandomAccessibleInterval< L > labeling,
+			final Shape shape,
+			final LongFunction< UnionFind > unionFindFactory,
+			final ToLongBiFunction< Localizable, L > idForPixel,
+			final LongUnaryOperator idForSet )
+	{
+		assert Intervals.contains( mask, labeling ) && Intervals.contains( labeling, mask ): "Mask and labeling are not the same size.";
+		assert Intervals.numElements( labeling ) <= Integer.MAX_VALUE: "Too many pixels for integer based union find.";
+
+		final UnionFind uf = makeUnion( mask, labeling, shape, unionFindFactory, idForPixel );
+		UnionFind.relabel( mask, labeling, uf, idForPixel, idForSet );
+	}
+
+	private static < B extends BooleanType< B >, L extends IntegerType< L > > UnionFind makeUnion(
+			final RandomAccessibleInterval< B > mask,
+			final RandomAccessibleInterval< L > labeling,
+			final Shape shape,
+			final LongFunction< UnionFind > unionFindFactory,
+			final ToLongBiFunction< Localizable, L > idForPixel )
+	{
+		final UnionFind uf = unionFindFactory.apply( Intervals.numElements( labeling ) );
+		if ( shape instanceof DiamondShape && ( ( DiamondShape ) shape ).getRadius() == 1 )
+		{
+			connectedComponentsDiamondShape( mask, labeling, uf, idForPixel );
+		}
+		else
+		{
+			connectedComponentsGeneralShape( mask, labeling, shape, uf, idForPixel );
+		}
+		return uf;
+	}
+
+	private static < B extends BooleanType< B >, L extends IntegerType< L > > void connectedComponentsDiamondShape(
+			final RandomAccessible< B > mask,
+			final RandomAccessibleInterval< L > labeling,
+			final UnionFind uf,
+			final ToLongBiFunction< Localizable, L > id )
+	{
+		final int nDim = labeling.numDimensions();
+		final long[] min = Intervals.minAsLongArray( labeling );
+		final long[] max = Intervals.maxAsLongArray( labeling );
+
+		for ( int d = 0; d < nDim; ++d )
+		{
+			final long[] minPlusOne = min.clone();
+			final long[] maxMinusOne = max.clone();
+			minPlusOne[ d ] += 1;
+			maxMinusOne[ d ] -= 1;
+			final Cursor< B > lower = Views.interval( mask, new FinalInterval( min, maxMinusOne ) ).cursor();
+			final Cursor< B > upper = Views.interval( mask, new FinalInterval( minPlusOne, max ) ).cursor();
+			final Cursor< L > lowerL = Views.interval( labeling, new FinalInterval( min, maxMinusOne ) ).localizingCursor();
+			final Cursor< L > upperL = Views.interval( labeling, new FinalInterval( minPlusOne, max ) ).localizingCursor();
+
+			while ( lower.hasNext() )
+			{
+				final boolean l = lower.next().get();
+				final boolean u = upper.next().get();
+				lowerL.fwd();
+				upperL.fwd();
+				if ( l && u )
+				{
+					final long r1 = uf.findRoot( id.applyAsLong( lowerL, lowerL.get() ) );
+					final long r2 = uf.findRoot( id.applyAsLong( upperL, upperL.get() ) );
+					if ( r1 != r2 )
+					{
+						uf.join( r1, r2 );
+					}
+				}
+			}
+		}
+	}
+
+	private static < B extends BooleanType< B >, L extends IntegerType< L > > void connectedComponentsGeneralShape(
+			final RandomAccessibleInterval< B > mask,
+			final RandomAccessibleInterval< L > labeling,
+			final Shape shape,
+			final UnionFind uf,
+			final ToLongBiFunction< Localizable, L > id )
+	{
+
+		final RandomAccessible< Neighborhood< B > > maskNeighborhood = shape.neighborhoodsRandomAccessible( Views.extendZero( mask ) );
+		final RandomAccessible< Neighborhood< L > > labelingNeighborhood = shape.neighborhoodsRandomAccessible( Views.extendZero( labeling ) );
+
+		final Cursor< Neighborhood< B > > mnhCursor = Views.flatIterable( Views.interval( maskNeighborhood, labeling ) ).cursor();
+		final Cursor< Neighborhood< L > > lnhCursor = Views.flatIterable( Views.interval( labelingNeighborhood, labeling ) ).cursor();
+		final Cursor< B > maskCursor = Views.flatIterable( mask ).cursor();
+		final Cursor< L > labelCursor = Views.flatIterable( labeling ).localizingCursor();
+
+		while ( maskCursor.hasNext() )
+		{
+			final B center = maskCursor.next();
+			final L label = labelCursor.next();
+			final Neighborhood< B > mnh = mnhCursor.next();
+			final Neighborhood< L > lnh = lnhCursor.next();
+			if ( center.get() )
+			{
+				final long index = id.applyAsLong( labelCursor, label );
+				final Cursor< L > lnhc = lnh.localizingCursor();
+				final Cursor< B > mnhc = mnh.cursor();
+				while ( mnhc.hasNext() )
+				{
+					final L l = lnhc.next();
+					final B m = mnhc.next();
+					if ( m.get() )
+					{
+						final long r1 = uf.findRoot( index );
+						final long r2 = uf.findRoot( id.applyAsLong( lnhc, l ) );
+						if ( r1 != r2 )
+						{
+							uf.join( r1, r2 );
+						}
+					}
+				}
+			}
+		}
+
+	}
+
+}