Add Hessian eigenvalues

pom.xml

@@ -229,6 +229,11 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
 			<artifactId>trove4j</artifactId>
 			<version>3.0.3</version>
 		</dependency>
+		<dependency>
+			<groupId>org.ojalgo</groupId>
+			<artifactId>ojalgo</artifactId>
+			<version>43.0</version>
+		</dependency>
 
 		<!-- Test dependencies -->
 		<dependency>

src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java

@@ -11,13 +11,13 @@
  * %%
  * 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
@@ -34,20 +34,34 @@
 
 package net.imglib2.algorithm.gradient;
 
+import java.util.ArrayList;
+import java.util.List;
+import java.util.concurrent.Callable;
+import java.util.concurrent.ExecutionException;
+import java.util.concurrent.ExecutorService;
+import java.util.concurrent.Future;
+
 import net.imglib2.Cursor;
+import net.imglib2.FinalInterval;
 import net.imglib2.RandomAccess;
 import net.imglib2.RandomAccessible;
 import net.imglib2.RandomAccessibleInterval;
 import net.imglib2.type.numeric.NumericType;
 import net.imglib2.util.Intervals;
+import net.imglib2.view.IntervalView;
 import net.imglib2.view.Views;
 
+/**
+ * @author Tobias Pietzsch
+ * @author Philipp Hanslovsky
+ *
+ */
 public class PartialDerivative
 {
 	// nice version...
 	/**
 	 * Compute the partial derivative of source in a particular dimension.
-	 * 
+	 *
 	 * @param source
 	 *            source image, has to provide valid data in the interval of the
 	 *            gradient image plus a one pixel border in dimension.
@@ -69,10 +83,82 @@ public static < T extends NumericType< T > > void gradientCentralDifference2( fi
 		}
 	}
 
+	// parallel version...
+	/**
+	 * Compute the partial derivative of source in a particular dimension.
+	 *
+	 * @param source
+	 *            source image, has to provide valid data in the interval of the
+	 *            gradient image plus a one pixel border in dimension.
+	 * @param gradient
+	 *            output image
+	 * @param dimension
+	 *            along which dimension the partial derivatives are computed
+	 * @param nTasks
+	 *            Number of tasks for gradient computation.
+	 * @param es
+	 *            {@link ExecutorService} providing workers for gradient
+	 *            computation. Service is managed (created, shutdown) by caller.
+	 */
+	public static < T extends NumericType< T > > void gradientCentralDifferenceParallel(
+			final RandomAccessible< T > source,
+			final RandomAccessibleInterval< T > gradient,
+			final int dimension,
+			final int nTasks,
+			final ExecutorService es ) throws InterruptedException, ExecutionException
+	{
+		final int nDim = source.numDimensions();
+		if ( nDim < 2 )
+		{
+			gradientCentralDifference( source, gradient, dimension );
+			return;
+		}
+
+		long dimensionMax = Long.MIN_VALUE;
+		int dimensionArgMax = -1;
+
+		for ( int d = 0; d < nDim; ++d )
+		{
+			final long size = gradient.dimension( d );
+			if ( d != dimension && size > dimensionMax )
+			{
+				dimensionMax = size;
+				dimensionArgMax = d;
+			}
+		}
+
+		final long stepSize = Math.max( dimensionMax / nTasks, 1 );
+		final long stepSizeMinusOne = stepSize - 1;
+		final long min = gradient.min( dimensionArgMax );
+		final long max = gradient.max( dimensionArgMax );
+
+		final ArrayList< Callable< Void > > tasks = new ArrayList<>();
+		for ( long currentMin = min, minZeroBase = 0; minZeroBase < dimensionMax; currentMin += stepSize, minZeroBase += stepSize )
+		{
+			final long currentMax = Math.min( currentMin + stepSizeMinusOne, max );
+			final long[] mins = new long[ nDim ];
+			final long[] maxs = new long[ nDim ];
+			gradient.min( mins );
+			gradient.max( maxs );
+			mins[ dimensionArgMax ] = currentMin;
+			maxs[ dimensionArgMax ] = currentMax;
+			final IntervalView< T > currentInterval = Views.interval( gradient, new FinalInterval( mins, maxs ) );
+			tasks.add( () -> {
+				gradientCentralDifference( source, currentInterval, dimension );
+				return null;
+			} );
+		}
+
+		final List< Future< Void > > futures = es.invokeAll( tasks );
+
+		for ( final Future< Void > f : futures )
+			f.get();
+	}
+
 	// fast version
 	/**
 	 * Compute the partial derivative of source in a particular dimension.
-	 * 
+	 *
 	 * @param source
 	 *            source image, has to provide valid data in the interval of the
 	 *            gradient image plus a one pixel border in dimension.

src/main/java/net/imglib2/algorithm/gradient/ScaleAsFunctionOfPosition.java

@@ -0,0 +1,110 @@
+/*
+ * #%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.gradient;
+
+import java.util.function.ToDoubleFunction;
+
+import net.imglib2.Interval;
+import net.imglib2.Localizable;
+import net.imglib2.RandomAccess;
+import net.imglib2.RandomAccessible;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.converter.AbstractConvertedRandomAccess;
+import net.imglib2.converter.AbstractConvertedRandomAccessible;
+import net.imglib2.type.Type;
+import net.imglib2.type.operators.MulFloatingPoint;
+
+/**
+ *
+ * @author Philipp Hanslovsky
+ *
+ *         Multiply the value of a {@link RandomAccessible} depending on the
+ *         position of the access.
+ *
+ */
+public class ScaleAsFunctionOfPosition< T extends Type< T > & MulFloatingPoint > extends AbstractConvertedRandomAccessible< T, T >
+{
+
+	private final ToDoubleFunction< Localizable > scalingFunction;
+
+	public ScaleAsFunctionOfPosition( final RandomAccessibleInterval< T > source, final ToDoubleFunction< Localizable > scalingFunction )
+	{
+		super( source );
+		this.scalingFunction = scalingFunction;
+	}
+
+	@Override
+	public ScaledRandomAccess< T > randomAccess()
+	{
+		return new ScaledRandomAccess<>( source.randomAccess(), scalingFunction );
+	}
+
+	@Override
+	public AbstractConvertedRandomAccess< T, T > randomAccess( final Interval interval )
+	{
+		return randomAccess();
+	}
+
+	public static class ScaledRandomAccess< T extends Type< T > & MulFloatingPoint > extends AbstractConvertedRandomAccess< T, T >
+	{
+
+		private final T t;
+
+		private final ToDoubleFunction< Localizable > scalingFunction;
+
+		public ScaledRandomAccess( final RandomAccess< T > source, final ToDoubleFunction< Localizable > scalingFunction )
+		{
+			super( source );
+			this.t = source.get().createVariable();
+			this.scalingFunction = scalingFunction;
+		}
+
+		@Override
+		public T get()
+		{
+			t.set( source.get() );
+			t.mul( scalingFunction.applyAsDouble( source ) );
+			return t;
+		}
+
+		@Override
+		public ScaledRandomAccess< T > copy()
+		{
+			return new ScaledRandomAccess<>( source.copyRandomAccess(), scalingFunction );
+		}
+
+	}
+
+}

src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues.java

@@ -0,0 +1,88 @@
+/*
+ * #%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.linalg.eigen;
+
+import net.imglib2.type.numeric.ComplexType;
+import net.imglib2.type.numeric.RealType;
+import net.imglib2.view.composite.Composite;
+
+/**
+ *
+ * Interface for handling different cases, e.g. square, symmetric, or 2
+ * dimensional tensors.
+ *
+ */
+public interface EigenValues< T extends RealType< T >, U extends ComplexType< U > >
+{
+	public void compute( final Composite< T > matrix, final Composite< U > evs );
+
+	public default EigenValues< T, U > copy()
+	{
+		return this;
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues1D< T, U > oneDimensional()
+	{
+		return new EigenValues1D<>();
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues2DSquare< T, U > square2D()
+	{
+		return new EigenValues2DSquare<>();
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues2DSymmetric< T, U > symmetric2D()
+	{
+		return new EigenValues2DSymmetric<>();
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValuesSquare< T, U > square( final int nDim )
+	{
+		return new EigenValuesSquare<>( nDim );
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValuesSymmetric< T, U > symmetric( final int nDim )
+	{
+		return new EigenValuesSymmetric<>( nDim );
+	}
+
+	public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues< T, U > invalid()
+	{
+		return ( m, evs ) -> {
+			throw new UnsupportedOperationException( "EigenValues not implemented yet!" );
+		};
+	}
+
+}

src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues1D.java

@@ -0,0 +1,48 @@
+/*
+ * #%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.linalg.eigen;
+
+import net.imglib2.type.numeric.ComplexType;
+import net.imglib2.type.numeric.RealType;
+import net.imglib2.view.composite.Composite;
+
+public class EigenValues1D< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U >
+{
+	@Override
+	public void compute( final Composite< T > tensor, final Composite< U > evs )
+	{
+		evs.get( 0 ).setReal( tensor.get( 0 ).getRealDouble() );
+	}
+}