Merge pull request #46 from maarzt/seperable-convolution

Seperable Convolution with Asymmetric Kernels + Fast Gauss

Author: tpietzsch

pom.xml

@@ -200,7 +200,7 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
 		<!-- NB: Deploy releases to the ImageJ Maven repository. -->
 		<releaseProfiles>deploy-to-imagej</releaseProfiles>
 
-		<imglib2.version>5.1.0</imglib2.version>
+		<imglib2.version>5.6.0</imglib2.version>
 		<imglib2-realtransform.version>2.0.0-beta-39</imglib2-realtransform.version>
 		<jitk-tps.version>3.0.0</jitk-tps.version>
 	</properties>
@@ -248,6 +248,18 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
 			<artifactId>junit</artifactId>
 			<scope>test</scope>
 		</dependency>
+		<dependency>
+			<groupId>org.openjdk.jmh</groupId>
+			<artifactId>jmh-core</artifactId>
+			<version>1.19</version>
+			<scope>test</scope>
+		</dependency>
+		<dependency>
+			<groupId>org.openjdk.jmh</groupId>
+			<artifactId>jmh-generator-annprocess</artifactId>
+			<version>1.19</version>
+			<scope>test</scope>
+		</dependency>
 	</dependencies>
 
 	<build>

src/main/java/net/imglib2/algorithm/convolution/AbstractMultiThreadedConvolution.java

@@ -0,0 +1,70 @@
+package net.imglib2.algorithm.convolution;
+
+import java.util.concurrent.ExecutorService;
+import java.util.concurrent.Executors;
+import java.util.concurrent.ThreadPoolExecutor;
+
+import net.imglib2.RandomAccessible;
+import net.imglib2.RandomAccessibleInterval;
+
+/**
+ * Abstract class to help implementing a Convolution, that is multi threaded
+ * using an {@link ExecutorService}. This implements the method
+ * {@link Convolution#setExecutor(ExecutorService)}.
+ * <p>
+ * Classes that derive from {@link AbstractMultiThreadedConvolution} must
+ * override
+ * {@link AbstractMultiThreadedConvolution#process(RandomAccessible, RandomAccessibleInterval, ExecutorService, int)}
+ *
+ * @author Matthias Arzt
+ */
+public abstract class AbstractMultiThreadedConvolution< T > implements Convolution< T >
+{
+
+	private ExecutorService executor;
+
+	abstract protected void process( RandomAccessible< ? extends T > source,
+			RandomAccessibleInterval< ? extends T > target,
+			ExecutorService executorService,
+			int numThreads );
+
+	@Override
+	public void setExecutor( final ExecutorService executor )
+	{
+		this.executor = executor;
+	}
+
+	@Override
+	final public void process( final RandomAccessible< ? extends T > source, final RandomAccessibleInterval< ? extends T > target )
+	{
+		if ( executor == null )
+		{
+			final int numThreads = suggestNumThreads();
+			final ExecutorService executor = Executors.newFixedThreadPool( numThreads );
+			try
+			{
+				process( source, target, executor, numThreads );
+			}
+			finally
+			{
+				executor.shutdown();
+			}
+		}
+		else
+		{
+			process( source, target, executor, getNumThreads( executor ) );
+		}
+	}
+
+	private int getNumThreads( final ExecutorService executor )
+	{
+		if ( executor instanceof ThreadPoolExecutor )
+			return ( ( ThreadPoolExecutor ) executor ).getMaximumPoolSize();
+		return suggestNumThreads();
+	}
+
+	private int suggestNumThreads()
+	{
+		return Runtime.getRuntime().availableProcessors();
+	}
+}

src/main/java/net/imglib2/algorithm/convolution/Concatenation.java

@@ -0,0 +1,116 @@
+package net.imglib2.algorithm.convolution;
+
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.List;
+import java.util.concurrent.ExecutorService;
+
+import net.imglib2.Interval;
+import net.imglib2.RandomAccessible;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.img.Img;
+import net.imglib2.type.NativeType;
+import net.imglib2.util.Intervals;
+import net.imglib2.util.Pair;
+import net.imglib2.util.Util;
+import net.imglib2.util.ValuePair;
+import net.imglib2.view.Views;
+
+/**
+ * Helper to implement {@link Convolution#concat}.
+ *
+ * @author Matthias Arzt
+ */
+class Concatenation< T > implements Convolution< T >
+{
+
+	private final List< Convolution< T > > steps;
+
+	Concatenation( final List< ? extends Convolution< T > > steps )
+	{
+		this.steps = new ArrayList<>( steps );
+	}
+
+	@Override
+	public void setExecutor( final ExecutorService executor )
+	{
+		steps.forEach( step -> step.setExecutor( executor ) );
+	}
+
+	@Override
+	public Interval requiredSourceInterval( final Interval targetInterval )
+	{
+		Interval result = targetInterval;
+		for ( int i = steps.size() - 1; i >= 0; i-- )
+			result = steps.get( i ).requiredSourceInterval( result );
+		return result;
+	}
+
+	@Override
+	public T preferredSourceType( T targetType )
+	{
+		for ( int i = steps.size() - 1; i >= 0; i-- )
+			targetType = steps.get( i ).preferredSourceType( targetType );
+		return targetType;
+	}
+
+	@Override
+	public void process( final RandomAccessible< ? extends T > source, final RandomAccessibleInterval< ? extends T > target )
+	{
+		final List< Pair< T, Interval > > srcIntervals = tmpIntervals( Util.getTypeFromInterval( target ), target );
+		RandomAccessibleInterval< ? extends T > currentSource = Views.interval( source, srcIntervals.get( 0 ).getB() );
+		RandomAccessibleInterval< ? extends T > available = null;
+
+		for ( int i = 0; i < steps.size(); i++ )
+		{
+			final Convolution< T > step = steps.get( i );
+			final T targetType = srcIntervals.get( i + 1 ).getA();
+			final Interval targetInterval = srcIntervals.get( i + 1 ).getB();
+			RandomAccessibleInterval< ? extends T > currentTarget =
+					( i == steps.size() - 1 ) ? target : null;
+
+			if ( currentTarget == null && available != null &&
+					Intervals.contains( available, targetInterval ) &&
+					Util.getTypeFromInterval( available ).getClass().equals( targetType.getClass() ) )
+				currentTarget = Views.interval( available, targetInterval );
+
+			if ( currentTarget == null )
+				currentTarget = createImage( uncheckedCast( targetType ), targetInterval );
+
+			step.process( currentSource, currentTarget );
+
+			if ( i > 0 )
+				available = currentSource;
+			currentSource = currentTarget;
+		}
+	}
+
+	private static < T extends NativeType< T > > RandomAccessibleInterval< T > createImage( final T targetType, final Interval targetInterval )
+	{
+		final long[] dimensions = Intervals.dimensionsAsLongArray( targetInterval );
+		final Img< T > ts = Util.getArrayOrCellImgFactory( targetInterval, targetType ).create( dimensions );
+		return Views.translate( ts, Intervals.minAsLongArray( targetInterval ) );
+	}
+
+	private static < T > T uncheckedCast( final Object in )
+	{
+		@SuppressWarnings( "unchecked" )
+		final
+		T in1 = ( T ) in;
+		return in1;
+	}
+
+	private List< Pair< T, Interval > > tmpIntervals( T type, Interval interval )
+	{
+		final List< Pair< T, Interval > > result = new ArrayList<>( Collections.nCopies( steps.size() + 1, null ) );
+		result.set( steps.size(), new ValuePair<>( type, interval ) );
+		for ( int i = steps.size() - 1; i >= 0; i-- )
+		{
+			final Convolution< T > step = steps.get( i );
+			interval = step.requiredSourceInterval( interval );
+			type = step.preferredSourceType( type );
+			result.set( i, new ValuePair<>( type, interval ) );
+		}
+		return result;
+	}
+}

src/main/java/net/imglib2/algorithm/convolution/Convolution.java

@@ -0,0 +1,69 @@
+package net.imglib2.algorithm.convolution;
+
+import java.util.Arrays;
+import java.util.List;
+import java.util.concurrent.ExecutorService;
+
+import net.imglib2.Interval;
+import net.imglib2.RandomAccessible;
+import net.imglib2.RandomAccessibleInterval;
+
+/**
+ * This interface allows the client to perform a convolution. But also to query
+ * for the required input image size and preferred input image type. The
+ * ExectorService can be set, to allow multi or single threaded operation.
+ * <p>
+ * Very importantly, multiple {@link Convolution}s can be easily concatenated.
+ *
+ * @author Matthias Arzt
+ */
+public interface Convolution< T >
+{
+	/**
+	 * Returns the required size for source image, to calculate the given target
+	 * interval.
+	 */
+	Interval requiredSourceInterval( Interval targetInterval );
+
+	/**
+	 * What's the preferred type for the source image, when target should have
+	 * the specified type?
+	 */
+	T preferredSourceType( T targetType );
+
+	/**
+	 * Set the {@link ExecutorService} to be used for convolution.
+	 */
+	default void setExecutor( final ExecutorService executor )
+	{}
+
+	/**
+	 * Fills the target image, with the convolution result.
+	 *
+	 * @param source
+	 *            Source image. It must allow pixel access in the interval
+	 *            returned by {@code requiredSourceInterval(target)}
+	 * @param target
+	 *            Target image.
+	 */
+	void process( RandomAccessible< ? extends T > source, RandomAccessibleInterval< ? extends T > target );
+
+	/**
+	 * Concatenate multiple {@link Convolution}s to one convolution. (e.g.
+	 * Concatenation of a gauss convolution in X, and a gauss convolution in Y
+	 * will be a 2d gauss convolution).
+	 */
+	static < T > Convolution< T > concat( final Convolution< T >... steps )
+	{
+		return concat( Arrays.asList( steps ) );
+	}
+
+	static < T > Convolution< T > concat( final List< ? extends Convolution< T > > steps )
+	{
+		if ( steps.isEmpty() )
+			throw new IllegalArgumentException( "Concat requires at least one convolution operation." );
+		if ( steps.size() == 1 )
+			return steps.get( 0 );
+		return new Concatenation<>( steps );
+	}
+}