diff --git a/pom.xml b/pom.xml index 33e419052..510508547 100644 --- a/pom.xml +++ b/pom.xml @@ -229,6 +229,11 @@ Jean-Yves Tinevez and Michael Zinsmaier. trove4j 3.0.3 + + org.ojalgo + ojalgo + 43.0 + diff --git a/src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java b/src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java new file mode 100644 index 000000000..10f738dab --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java @@ -0,0 +1,389 @@ +/* + * #%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.concurrent.ExecutionException; +import java.util.concurrent.ExecutorService; +import java.util.stream.IntStream; + +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.View; +import net.imglib2.algorithm.gauss3.Gauss3; +import net.imglib2.exception.IncompatibleTypeException; +import net.imglib2.outofbounds.OutOfBoundsFactory; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Util; +import net.imglib2.view.IntervalView; +import net.imglib2.view.MixedTransformView; +import net.imglib2.view.Views; + +/** + * + * Compute entries of n-dimensional Hessian matrix. + * + * @author Philipp Hanslovsky + * + */ +public class HessianMatrix +{ + + /** + * @param source + * n-dimensional input {@link RandomAccessibleInterval} + * @param gaussian + * n-dimensional output parameter + * {@link RandomAccessibleInterval} + * @param gradient + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * the gradients along all axes of the smoothed source (size of + * last dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * @param outOfBounds + * {@link OutOfBoundsFactory} that specifies how out of bound + * pixels of intermediate results should be handled (necessary + * for gradient computation). + * @param sigma + * Scale for Gaussian smoothing. + * + * @return Hessian matrix that was passed as output parameter. + * @throws IncompatibleTypeException + */ + public static < T extends RealType< T >, U extends RealType< U > > RandomAccessibleInterval< U > calculateMatrix( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > gaussian, + final RandomAccessibleInterval< U > gradient, + final RandomAccessibleInterval< U > hessianMatrix, + final OutOfBoundsFactory< U, ? super RandomAccessibleInterval< U > > outOfBounds, + final double... sigma ) throws IncompatibleTypeException + { + + if ( sigma.length == 1 ) + Gauss3.gauss( sigma[ 0 ], source, gaussian ); + else + Gauss3.gauss( sigma, source, gaussian ); + return calculateMatrix( Views.extend( gaussian, outOfBounds ), gradient, hessianMatrix, outOfBounds ); + } + + /** + * @param source + * n-dimensional pre-smoothed {@link RandomAccessible}. It is the + * callers responsibility to smooth the input at the desired + * scales. + * @param gradient + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * the gradients along all axes of the smoothed source (size of + * last dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * + * @param outOfBounds + * {@link OutOfBoundsFactory} that specifies how out of bound + * pixels of intermediate results should be handled (necessary + * for gradient computation). + * + * @return Hessian matrix that was passed as output parameter. + */ + public static < T extends RealType< T > > RandomAccessibleInterval< T > calculateMatrix( + final RandomAccessible< T > source, + final RandomAccessibleInterval< T > gradient, + final RandomAccessibleInterval< T > hessianMatrix, + final OutOfBoundsFactory< T, ? super RandomAccessibleInterval< T > > outOfBounds ) + { + + final int nDim = gradient.numDimensions() - 1; + + for ( long d = 0; d < nDim; ++d ) + PartialDerivative.gradientCentralDifference( source, Views.hyperSlice( gradient, nDim, d ), ( int ) d ); + + return calculateMatrix( Views.extend( gradient, outOfBounds ), hessianMatrix ); + } + + /** + * + * @param gradient + * n+1-dimensional {@link RandomAccessible} containing the + * gradients along all axes of the smoothed source (size of last + * dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * + * @return Hessian matrix that was passed as output parameter. + */ + public static < T extends RealType< T > > RandomAccessibleInterval< T > calculateMatrix( + final RandomAccessible< T > gradients, + final RandomAccessibleInterval< T > hessianMatrix ) + { + + final int nDim = gradients.numDimensions() - 1; + + long count = 0; + for ( long d1 = 0; d1 < nDim; ++d1 ) + { + final MixedTransformView< T > hs1 = Views.hyperSlice( gradients, nDim, d1 ); + for ( long d2 = d1; d2 < nDim; ++d2 ) + { + final IntervalView< T > hs2 = Views.hyperSlice( hessianMatrix, nDim, count ); + PartialDerivative.gradientCentralDifference( hs1, hs2, ( int ) d2 ); + ++count; + } + } + return hessianMatrix; + } + + // parallel + + /** + * @param source + * n-dimensional input {@link RandomAccessibleInterval} + * @param gaussian + * n-dimensional output parameter + * {@link RandomAccessibleInterval} + * @param gradient + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * the gradients along all axes of the smoothed source (size of + * last dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * @param outOfBounds + * {@link OutOfBoundsFactory} that specifies how out of bound + * pixels of intermediate results should be handled (necessary + * for gradient computation). + * @param nTasks + * Number of tasks used for parallel computation of eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * @param sigma + * Scale for Gaussian smoothing. + * + * @return Hessian matrix that was passed as output parameter. + * @throws IncompatibleTypeException + * @throws ExecutionException + * @throws InterruptedException + */ + public static < T extends RealType< T >, U extends RealType< U > > RandomAccessibleInterval< U > calculateMatrix( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > gaussian, + final RandomAccessibleInterval< U > gradient, + final RandomAccessibleInterval< U > hessianMatrix, + final OutOfBoundsFactory< U, ? super RandomAccessibleInterval< U > > outOfBounds, + final int nTasks, + final ExecutorService es, + final double... sigma ) throws IncompatibleTypeException, InterruptedException, ExecutionException + { + + if ( sigma.length == 1 ) + Gauss3.gauss( IntStream.range( 0, source.numDimensions() ).mapToDouble( i -> sigma[ 0 ] ).toArray(), source, gaussian, es ); + else + Gauss3.gauss( sigma, source, gaussian, es ); + return calculateMatrix( Views.extend( gaussian, outOfBounds ), gradient, hessianMatrix, outOfBounds, nTasks, es ); + } + + /** + * @param source + * n-dimensional pre-smoothed {@link RandomAccessible}. It is the + * callers responsibility to smooth the input at the desired + * scales. + * @param gradient + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * the gradients along all axes of the smoothed source (size of + * last dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * @param outOfBounds + * {@link OutOfBoundsFactory} that specifies how out of bound + * pixels of intermediate results should be handled (necessary + * for gradient computation). + * @param nTasks + * Number of tasks used for parallel computation of eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * + * @return Hessian matrix that was passed as output parameter. + * @throws IncompatibleTypeException + * @throws ExecutionException + * @throws InterruptedException + */ + public static < T extends RealType< T > > RandomAccessibleInterval< T > calculateMatrix( + final RandomAccessible< T > source, + final RandomAccessibleInterval< T > gradient, + final RandomAccessibleInterval< T > hessianMatrix, + final OutOfBoundsFactory< T, ? super RandomAccessibleInterval< T > > outOfBounds, + final int nTasks, + final ExecutorService es ) throws IncompatibleTypeException, InterruptedException, ExecutionException + { + + final int nDim = gradient.numDimensions() - 1; + + for ( long d = 0; d < nDim; ++d ) + PartialDerivative.gradientCentralDifferenceParallel( source, Views.hyperSlice( gradient, nDim, d ), ( int ) d, nTasks, es ); + + return calculateMatrix( Views.extend( gradient, outOfBounds ), hessianMatrix, nTasks, es ); + } + + /** + * + * @param gradient + * n+1-dimensional {@link RandomAccessible} containing the + * gradients along all axes of the smoothed source (size of last + * dimension is n) + * @param hessianMatrix + * n+1-dimensional {@link RandomAccessibleInterval} for storing + * all second partial derivatives as a linear representation + * (size of last dimension is n * ( n + 1 ) / 2) of upper + * triangular Hessian matrix: For n-dimensional input, + * hessianMatrix (m) will be populated along the nth + * dimension like this: [m11, m12, ... , m1n, m22, m23, ... , + * mnn] + * @param nTasks + * Number of tasks used for parallel computation of eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * + * @return Hessian matrix that was passed as output parameter. + * @throws IncompatibleTypeException + * @throws ExecutionException + * @throws InterruptedException + */ + public static < T extends RealType< T > > RandomAccessibleInterval< T > calculateMatrix( + final RandomAccessible< T > gradient, + final RandomAccessibleInterval< T > hessianMatrix, + final int nTasks, + final ExecutorService es ) throws IncompatibleTypeException, InterruptedException, ExecutionException + { + + final int nDim = gradient.numDimensions() - 1; + + long count = 0; + for ( long d1 = 0; d1 < nDim; ++d1 ) + { + final MixedTransformView< T > hs1 = Views.hyperSlice( gradient, nDim, d1 ); + for ( long d2 = d1; d2 < nDim; ++d2 ) + { + final IntervalView< T > hs2 = Views.hyperSlice( hessianMatrix, nDim, count ); + PartialDerivative.gradientCentralDifferenceParallel( hs1, hs2, ( int ) d2, nTasks, es ); + ++count; + } + } + return hessianMatrix; + } + + + + /** + * + * The {@link HessianMatrix#calculateMatrix} methods assume that the voxels + * of the input data are isotropic. This is not always the case and the + * finite differences would need to be normalized accordingly. As the + * derivative is a linear transformation re-normalization (scaling) of the + * Hessian matrix produces the desired result. In general, the normalization + * is arbitrary and we choose to normalize such that the second derivative + * remains unchanged along the dimension for which the extent of a voxel is + * the smallest. + * + * Note that the returned object is a {@link View} on the Hessian matrix and + * values are re-normalized as they are queried. For repeated use of the + * re-normalized Hessian matrix, it might be beneficial to persist the + * re-normalized values by copying into a {@link RandomAccessibleInterval}. + * + * @param hessian + * Hessian matrix to be re-normalized (scaled). + * hessian is an n+1 dimensional + * {@link RandomAccessibleInterval} that stores a linear + * representation of the upper triangular n-dimensional Hessian + * matrix along the last dimension as specified in + * {@link HessianMatrix#calculateMatrix(RandomAccessible, RandomAccessibleInterval)}. + * @param sigma + * Specify the voxel size for each dimension. + * @return Scaled {@link View} of the hessian matrix stored in + * hessian + */ + public static < T extends RealType< T > > IntervalView< T > scaleHessianMatrix( final RandomAccessibleInterval< T > hessian, final double[] sigma ) + { + + assert sigma.length == hessian.numDimensions() - 1; + assert sigma.length * ( sigma.length + 1 ) / 2 == hessian.dimension( sigma.length ); + final int maxD = sigma.length; + + final double minSigma = Util.min( sigma ); + final double minSigmaSq = minSigma * minSigma; + final double[] sigmaSquared = new double[ sigma.length * ( sigma.length + 1 ) / 2 ]; + for ( int i1 = 0, k = 0; i1 < sigma.length; ++i1 ) + for ( int i2 = i1; i2 < sigma.length; ++i2, ++k ) + sigmaSquared[ k ] = sigma[ i1 ] * sigma[ i2 ] / minSigmaSq; + + final ScaleAsFunctionOfPosition< T > scaledMatrix = new ScaleAsFunctionOfPosition<>( hessian, l -> { + return sigmaSquared[ l.getIntPosition( maxD ) ]; + } ); + + return Views.interval( scaledMatrix, hessian ); + } + + +} diff --git a/src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java b/src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java index 31c5c6136..b281077bb 100644 --- a/src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java +++ b/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. diff --git a/src/main/java/net/imglib2/algorithm/gradient/ScaleAsFunctionOfPosition.java b/src/main/java/net/imglib2/algorithm/gradient/ScaleAsFunctionOfPosition.java new file mode 100644 index 000000000..0624a091c --- /dev/null +++ b/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 ); + } + + } + +} diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues.java new file mode 100644 index 000000000..84808a896 --- /dev/null +++ b/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!" ); + }; + } + +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues1D.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues1D.java new file mode 100644 index 000000000..ea6d3937e --- /dev/null +++ b/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() ); + } +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSquare.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSquare.java new file mode 100644 index 000000000..7cf853e06 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSquare.java @@ -0,0 +1,67 @@ +/* + * #%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 EigenValues2DSquare< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U > +{ + @Override + public void compute( final Composite< T > tensor, final Composite< U > evs ) + { + final double x11 = tensor.get( 0 ).getRealDouble(); + final double x12 = tensor.get( 1 ).getRealDouble(); + final double x21 = tensor.get( 2 ).getRealDouble(); + final double x22 = tensor.get( 3 ).getRealDouble(); + final double sum = x11 + x22; + final double diff = x11 - x22; + final double radicand = 4 * x12 * x21 + diff * diff; + if ( radicand < 0.0d ) + { + final double halfSqrt = 0.5 * Math.sqrt( -radicand ); + final double halfSum = 0.5 * sum; + evs.get( 0 ).setComplexNumber( halfSum, halfSqrt ); + evs.get( 1 ).setComplexNumber( halfSum, -halfSqrt ); + } + else + { + final double sqrt = Math.sqrt( radicand ); + evs.get( 0 ).setComplexNumber( 0.5 * ( sum + sqrt ), 0.0d ); + evs.get( 1 ).setComplexNumber( 0.5 * ( sum - sqrt ), 0.0d ); + } + } +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSymmetric.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSymmetric.java new file mode 100644 index 000000000..ea72c800d --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValues2DSymmetric.java @@ -0,0 +1,55 @@ +/* + * #%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 EigenValues2DSymmetric< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U > +{ + @Override + public void compute( final Composite< T > tensor, final Composite< U > evs ) + { + final double x11 = tensor.get( 0 ).getRealDouble(); + final double x12 = tensor.get( 1 ).getRealDouble(); + final double x22 = tensor.get( 2 ).getRealDouble(); + final double sum = x11 + x22; + final double diff = x11 - x22; + final double sqrt = Math.sqrt( 4 * x12 * x12 + diff * diff ); + evs.get( 0 ).setReal( 0.5 * ( sum + sqrt ) ); + evs.get( 1 ).setReal( 0.5 * ( sum - sqrt ) ); + } +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSquare.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSquare.java new file mode 100644 index 000000000..5953eca3e --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSquare.java @@ -0,0 +1,82 @@ +/* + * #%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 java.util.Optional; + +import org.ojalgo.matrix.decomposition.Eigenvalue; + +import net.imglib2.algorithm.linalg.matrix.RealCompositeSquareMatrix; +import net.imglib2.type.numeric.ComplexType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.composite.Composite; + +public class EigenValuesSquare< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U > +{ + private final int nDim; + + private final double[] evs; + + private final RealCompositeSquareMatrix< T > m; + + private final Eigenvalue< Double > ed; + + private final Optional< double[] > emptyOptional = Optional.empty(); + + public EigenValuesSquare( final int nDim ) + { + super(); + this.nDim = nDim; + this.evs = new double[ nDim ]; + this.m = new RealCompositeSquareMatrix<>( null, nDim ); + this.ed = Eigenvalue.PRIMITIVE.make( this.m, false ); + } + + @Override + public void compute( final Composite< T > tensor, final Composite< U > evs ) + { + m.setData( tensor ); + ed.computeValuesOnly( this.m ); + ed.getEigenvalues( this.evs, emptyOptional ); + for ( int d = 0; d < this.evs.length; ++d ) + evs.get( d ).setReal( this.evs[ d ] ); + } + + @Override + public EigenValuesSquare< T, U > copy() + { + return new EigenValuesSquare<>( nDim ); + } +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSymmetric.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSymmetric.java new file mode 100644 index 000000000..398bbe7a4 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/eigen/EigenValuesSymmetric.java @@ -0,0 +1,82 @@ +/* + * #%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 java.util.Optional; + +import org.ojalgo.matrix.decomposition.Eigenvalue; + +import net.imglib2.algorithm.linalg.matrix.RealCompositeSymmetricMatrix; +import net.imglib2.type.numeric.ComplexType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.composite.Composite; + +public class EigenValuesSymmetric< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U > +{ + private final int nDim; + + private final double[] evs; + + private final RealCompositeSymmetricMatrix< T > m; + + private final Eigenvalue< Double > ed; + + private final Optional< double[] > emptyOptional = Optional.empty(); + + public EigenValuesSymmetric( final int nDim ) + { + super(); + this.nDim = nDim; + this.evs = new double[ nDim ]; + this.m = new RealCompositeSymmetricMatrix<>( null, nDim ); + this.ed = Eigenvalue.PRIMITIVE.make( this.m, true ); + } + + @Override + public void compute( final Composite< T > tensor, final Composite< U > evs ) + { + m.setData( tensor ); + ed.computeValuesOnly( m ); + ed.getEigenvalues( this.evs, emptyOptional ); + for ( int d = 0; d < this.evs.length; ++d ) + evs.get( d ).setReal( this.evs[ d ] ); + } + + @Override + public EigenValuesSymmetric< T, U > copy() + { + return new EigenValuesSymmetric<>( nDim ); + } +} \ No newline at end of file diff --git a/src/main/java/net/imglib2/algorithm/linalg/eigen/TensorEigenValues.java b/src/main/java/net/imglib2/algorithm/linalg/eigen/TensorEigenValues.java new file mode 100644 index 000000000..fe49ca503 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/eigen/TensorEigenValues.java @@ -0,0 +1,393 @@ +/* + * #%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 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.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.ImgFactory; +import net.imglib2.type.numeric.ComplexType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; +import net.imglib2.view.composite.NumericComposite; +import net.imglib2.view.composite.RealComposite; + +/** + * + * @author Philipp Hanslovsky + * + * Compute eigenvalues of rank 2 tensors. + * + */ +public class TensorEigenValues +{ + + // static methods + + // symmetric + + /** + * + * @param tensor + * Input that holds linear representation of upper triangular + * tensor in last dimension, i.e. if tensor (t) has n+1 + * dimensions, the last dimension must be of size n * ( n + 1 ) / + * 2, and the entries in the last dimension are arranged like + * this: [t11, t12, ... , t1n, t22, t23, ... , tnn] + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * the last dimension must be n. + * + */ + public static < T extends RealType< T >, U extends RealType< U > > RandomAccessibleInterval< U > calculateEigenValuesSymmetric( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues ) + { + + final int nDim = tensor.numDimensions() - 1; + assert eigenvalues.dimension( nDim ) * ( eigenvalues.dimension( nDim ) + 1 ) / 2 == tensor.dimension( nDim ); + + final EigenValues< T, U > ev; + if ( nDim == 1 ) + ev = EigenValues.oneDimensional(); + else if ( nDim == 2 ) + ev = EigenValues.symmetric2D(); + else if ( nDim > 2 ) + ev = EigenValues.symmetric( nDim ); + else + ev = EigenValues.invalid(); + + return calculateEigenValuesImpl( tensor, eigenvalues, ev ); + } + + /** + * + * @param tensor + * Input that holds linear representation of upper triangular + * tensor in last dimension, i.e. if tensor (t) has n+1 + * dimensions, the last dimension must be of size n * ( n + 1 ) / + * 2, and the entries in the last dimension are arranged like + * this: [t11, t12, ... , t1n, t22, t23, ... , tnn] + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * the last dimension must be n. + * @param nTasks + * Number of tasks used for parallel computation of eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * + */ + public static < T extends RealType< T >, U extends RealType< U > > RandomAccessibleInterval< U > calculateEigenValuesSymmetric( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues, + final int nTasks, + final ExecutorService es ) + { + + final int nDim = tensor.numDimensions() - 1; + assert eigenvalues.dimension( nDim ) * ( eigenvalues.dimension( nDim ) + 1 ) / 2 == tensor.dimension( nDim ); + + final EigenValues< T, U > ev; + if ( nDim == 1 ) + ev = EigenValues.oneDimensional(); + else if ( nDim == 2 ) + ev = EigenValues.symmetric2D(); + else if ( nDim > 2 ) + ev = EigenValues.symmetric( nDim ); + else + ev = EigenValues.invalid(); + + return calculateEigenValues( tensor, eigenvalues, ev, nTasks, es ); + } + + // square + + /** + * + * @param tensor + * Input that holds linear representation of tensor in last + * dimension, i.e. if tensor (t) has n+1 dimensions, the last + * dimension must be of size n * n, and the entries in the last + * dimension are arranged like this: [t11, t12, ... , t1n, t21, + * t22, t23, ... , tn1, ... , tnn] + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * + */ + public static < T extends RealType< T >, U extends ComplexType< U > > RandomAccessibleInterval< U > calculateEigenValuesSquare( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues ) + { + + final int nDim = tensor.numDimensions() - 1; + assert eigenvalues.dimension( nDim ) * ( eigenvalues.dimension( nDim ) + 1 ) / 2 == tensor.dimension( nDim ); + + final EigenValues< T, U > ev; + if ( nDim == 1 ) + ev = EigenValues.oneDimensional(); + else if ( nDim == 2 ) + ev = EigenValues.square2D(); + else if ( nDim > 2 ) + ev = EigenValues.square( nDim ); + else + ev = EigenValues.invalid(); + + return calculateEigenValuesImpl( tensor, eigenvalues, ev ); + } + + /** + * + * @param tensor + * Input that holds linear representation of tensor in last + * dimension, i.e. if tensor (t) has n+1 dimensions, the last + * dimension must be of size n * n, and the entries in the last + * dimension are arranged like this: [t11, t12, ... , t1n, t21, + * t22, t23, ... , tn1, ... , tnn] + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * the last dimension must be n. + * @param nTasks + * Number of threads/workers used for parallel computation of + * eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * + */ + public static < T extends RealType< T >, U extends ComplexType< U > > RandomAccessibleInterval< U > calculateEigenValuesSquare( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues, + final int nTasks, + final ExecutorService es ) + { + final int nDim = tensor.numDimensions() - 1; + assert eigenvalues.dimension( nDim ) * eigenvalues.dimension( nDim ) == tensor.dimension( nDim ); + + final EigenValues< T, U > ev; + if ( nDim == 1 ) + ev = EigenValues.oneDimensional(); + else if ( nDim == 2 ) + ev = EigenValues.square2D(); + else if ( nDim > 2 ) + ev = EigenValues.square( nDim ); + else + ev = EigenValues.invalid(); + + return calculateEigenValues( tensor, eigenvalues, ev, nTasks, es ); + } + + // general + + /** + * + * @param tensor + * Input that holds linear representation of tensor in last + * dimension. Parameter ev specifies representation. + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * the last dimension must be n. + * @param ev + * Implementation that specifies how to calculate eigenvalues + * from last dimension of input. + * + */ + public static < T extends RealType< T >, U extends ComplexType< U > > RandomAccessibleInterval< U > calculateEigenValues( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues, + final EigenValues< T, U > ev ) + { + return calculateEigenValuesImpl( tensor, eigenvalues, ev ); + } + + /** + * + * @param tensor + * Input that holds linear representation of tensor in last + * dimension. Parameter ev specifies representation. + * + * @param eigenvalues + * Target {@link RandomAccessibleInterval} for storing the + * resulting tensor eigenvalues. Number of dimensions must be the + * same as for input. For an n+1 dimensional input, the size of + * the last dimension must be n. + * @param ev + * Implementation that specifies how to calculate eigenvalues + * from last dimension of input. + * @param nTasks + * Number of threads/workers used for parallel computation of + * eigenvalues. + * @param es + * {@link ExecutorService} providing workers for parallel + * computation. Service is managed (created, shutdown) by caller. + * + */ + public static < T extends RealType< T >, U extends ComplexType< U > > RandomAccessibleInterval< U > calculateEigenValues( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues, + final EigenValues< T, U > ev, + final int nTasks, + final ExecutorService es ) + { + + assert nTasks > 0: "Passed nTasks < 1"; + + final int tensorDims = tensor.numDimensions(); + + long dimensionMax = Long.MIN_VALUE; + int dimensionArgMax = -1; + + for ( int d = 0; d < tensorDims - 1; ++d ) + { + final long size = tensor.dimension( d ); + if ( size > dimensionMax ) + { + dimensionMax = size; + dimensionArgMax = d; + } + } + + final long stepSize = Math.max( dimensionMax / nTasks, 1 ); + final long stepSizeMinusOne = stepSize - 1; + final long max = dimensionMax - 1; + + final ArrayList< Callable< RandomAccessibleInterval< U > > > tasks = new ArrayList<>(); + for ( long currentMin = 0; currentMin < dimensionMax; currentMin += stepSize ) + { + final long currentMax = Math.min( currentMin + stepSizeMinusOne, max ); + final long[] minT = new long[ tensorDims ]; + final long[] maxT = new long[ tensorDims ]; + final long[] minE = new long[ tensorDims ]; + final long[] maxE = new long[ tensorDims ]; + tensor.min( minT ); + tensor.max( maxT ); + eigenvalues.min( minE ); + eigenvalues.max( maxE ); + minE[ dimensionArgMax ] = minT[ dimensionArgMax ] = currentMin; + maxE[ dimensionArgMax ] = maxT[ dimensionArgMax ] = currentMax; + final IntervalView< T > currentTensor = Views.interval( tensor, new FinalInterval( minT, maxT ) ); + final IntervalView< U > currentEigenvalues = Views.interval( eigenvalues, new FinalInterval( minE, maxE ) ); + tasks.add( () -> calculateEigenValuesImpl( currentTensor, currentEigenvalues, ev.copy() ) ); + } + + + try + { + final List< Future< RandomAccessibleInterval< U > > > futures = es.invokeAll( tasks ); + for ( final Future< RandomAccessibleInterval< U > > f : futures ) + try + { + f.get(); + } + catch ( final ExecutionException e ) + { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + catch ( final InterruptedException e ) + { + // TODO Auto-generated catch block + e.printStackTrace(); + } + + return eigenvalues; + + + + } + + private static < T extends RealType< T >, U extends ComplexType< U > > RandomAccessibleInterval< U > calculateEigenValuesImpl( + final RandomAccessibleInterval< T > tensor, + final RandomAccessibleInterval< U > eigenvalues, + final EigenValues< T, U > ev ) + { + final Cursor< RealComposite< T > > m = Views.iterable( Views.collapseReal( tensor ) ).cursor(); + final Cursor< NumericComposite< U > > e = Views.iterable( Views.collapseNumeric( eigenvalues ) ).cursor(); + while ( m.hasNext() ) + ev.compute( m.next(), e.next() ); + return eigenvalues; + } + + /** + * + * Create appropriately sized image for tensor input. + * + * @param tensor + * n+1 dimensional {@link RandomAccessibleInterval}. + * @param factory + * {@link ImgFactory} used for creating the result image. + * @param u + * Variable necessary for creating the result image. + * @return n+1 dimensional {@link Img} with size n in the last dimension. + */ + public static < T extends RealType< T >, U extends RealType< U > > Img< U > createAppropriateResultImg( + final RandomAccessibleInterval< T > tensor, + final ImgFactory< U > factory, + final U u ) + { + final int nDim = tensor.numDimensions(); + final long[] dimensions = new long[ nDim ]; + tensor.dimensions( dimensions ); + dimensions[ nDim - 1 ] = nDim - 1; + return factory.create( dimensions, u ); + } + + +} diff --git a/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeMatrix.java b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeMatrix.java new file mode 100644 index 000000000..bba0f74f0 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeMatrix.java @@ -0,0 +1,127 @@ +/* + * #%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.matrix; + +import org.ojalgo.access.Access2D; +import org.ojalgo.matrix.store.PhysicalStore; + +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.composite.Composite; + +/** + * + * @author Philipp Hanslovsky + * + * {@link RealMatrix} that reads data from {@link Composite} (non-copy). + * + * @param + */ +public class RealCompositeMatrix< T extends RealType< T > > implements Access2D< Double >, Access2D.Collectable< Double, PhysicalStore< Double > > +{ + + protected Composite< T > data; + + protected final int nRows; + + protected final int nCols; + + protected final int length; + + public RealCompositeMatrix( final Composite< T > data, final int nRows, final int nCols ) + { + this( data, nRows, nCols, nRows * nCols ); + } + + protected RealCompositeMatrix( final Composite< T > data, final int nRows, final int nCols, final int length ) + { + super(); + + this.data = data; + this.nRows = nRows; + this.nCols = nCols; + this.length = length; + } + + public void setData( final Composite< T > data ) + { + this.data = data; + } + + public long rowAndColumnToLinear( final long row, final long col ) + { + return row * nCols + col; + } + + @Override + public long countColumns() + { + return nCols; + } + + @Override + public long countRows() + { + return nRows; + } + + @Override + public void supplyTo( final PhysicalStore< Double > receiver ) + { + final long tmpLimRows = Math.min( this.countRows(), receiver.countRows() ); + final long tmpLimCols = Math.min( this.countColumns(), receiver.countColumns() ); + + for ( long j = 0L; j < tmpLimCols; j++ ) + for (long i = 0L; i < tmpLimRows; i++) + receiver.set(i, j, this.doubleValue(i, j)); + } + + @Override + public double doubleValue( final long row, final long col ) + { + assert row >= 0 && row < this.nRows; + assert col >= 0 && col < this.nCols; + + final double val = data.get( rowAndColumnToLinear( row, col ) ).getRealDouble(); + + return val; + } + + @Override + public Double get( final long row, final long col ) + { + return doubleValue( row, col ); + } + +} diff --git a/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSquareMatrix.java b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSquareMatrix.java new file mode 100644 index 000000000..292a65a98 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSquareMatrix.java @@ -0,0 +1,62 @@ +/* + * #%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.matrix; + +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.composite.Composite; + +/** + * + * @author Philipp Hanslovsky + * + * Square {@link RealMatrix} that reads data from {@link Composite} + * (non-copy). + * + * @param + */ +public class RealCompositeSquareMatrix< T extends RealType< T > > extends RealCompositeMatrix< T > +{ + + public RealCompositeSquareMatrix( final Composite< T > data, final int nRowsOrCols ) + { + this( data, nRowsOrCols, nRowsOrCols * nRowsOrCols ); + } + + protected RealCompositeSquareMatrix( final Composite< T > data, final int nRowsOrCols, final int length ) + { + super( data, nRowsOrCols, nRowsOrCols, length ); + } + +} diff --git a/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSymmetricMatrix.java b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSymmetricMatrix.java new file mode 100644 index 000000000..17b97c0b5 --- /dev/null +++ b/src/main/java/net/imglib2/algorithm/linalg/matrix/RealCompositeSymmetricMatrix.java @@ -0,0 +1,86 @@ +/* + * #%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.matrix; + +import net.imglib2.type.numeric.RealType; +import net.imglib2.view.composite.Composite; + +/** + * + * @author Philipp Hanslovsky + * + * Symmetric {@link RealMatrix} that reads data from {@link Composite} + * (non-copy). + * + * @param + */ +public class RealCompositeSymmetricMatrix< T extends RealType< T > > extends RealCompositeSquareMatrix< T > +{ + + public RealCompositeSymmetricMatrix( final Composite< T > data, final int nRowsOrCols ) + { + this( data, nRowsOrCols, nRowsOrCols * ( nRowsOrCols + 1 ) / 2 ); + } + + protected RealCompositeSymmetricMatrix( final Composite< T > data, final int nRowsOrCols, final int length ) + { + super( data, nRowsOrCols, length ); + } + + @Override + public long rowAndColumnToLinear( final long row, final long col ) + { + + // total number of elements: length = nRows * ( nRows + 1 ) / 2 + // row - 1 complete rows + // number elements in non-complete rows: n = ( nRows - ( row - 1 ) ) * ( + // nRows -row ) / 2 + // number of elements total: length - n + ( col - row ) + + if ( row < col ) + { + final long rowDiff = nRows - row; + final long n = rowDiff * ( rowDiff + 1 ) / 2; + return length - n + col - row; + } + else + { + final long rowDiff = nRows - col; + final long n = rowDiff * ( rowDiff + 1 ) / 2; + return length - n + row - col; + } + } + +} diff --git a/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixEigenValuesTest.java b/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixEigenValuesTest.java new file mode 100644 index 000000000..3c0060b23 --- /dev/null +++ b/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixEigenValuesTest.java @@ -0,0 +1,116 @@ +/* + * #%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.hessian; + +import org.junit.Assert; +import org.junit.Test; + +import net.imglib2.Cursor; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.gauss3.Gauss3; +import net.imglib2.algorithm.gradient.HessianMatrix; +import net.imglib2.algorithm.linalg.eigen.TensorEigenValues; +import net.imglib2.exception.IncompatibleTypeException; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.DoubleArray; +import net.imglib2.outofbounds.OutOfBoundsBorderFactory; +import net.imglib2.type.numeric.real.DoubleType; +import net.imglib2.view.Views; + +public class HessianMatrixEigenValuesTest +{ + + @Test + public void test2D() throws IncompatibleTypeException + { + final int size = 5; + + final double[] imgArray = new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 4.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + final double[] eigenvaluesRef1 = new double[] { + 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 1.0, 0.0, 1.0, 0.0, + 1.0, 0.0, -2.0, 0.0, 1.0, + 0.0, 1.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0 + }; + + final double[] eigenvaluesRef2 = new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, -1.0, 0.0, -1.0, 0.0, + 0.0, 0.0, -2.0, 0.0, 0.0, + 0.0, -1.0, 0.0, -1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + final ArrayImg< DoubleType, DoubleArray > img = ArrayImgs.doubles( imgArray, size, size ); + + final double sigma = 0.1; + + final ArrayImg< DoubleType, DoubleArray > gaussian = ArrayImgs.doubles( size, size ); + Gauss3.gauss( sigma, Views.extendBorder( img ), gaussian ); + + final RandomAccessibleInterval< DoubleType > hessian = + HessianMatrix.calculateMatrix( Views.extendBorder( gaussian ), ArrayImgs.doubles( size, size, 2 ), ArrayImgs.doubles( size, size, 3 ), new OutOfBoundsBorderFactory<>() ); + + final RandomAccessibleInterval< DoubleType > eigenvalues = TensorEigenValues.calculateEigenValuesSymmetric( + hessian, + TensorEigenValues.createAppropriateResultImg( hessian, new ArrayImgFactory<>(), new DoubleType() ) ); + + Assert.assertEquals( img.numDimensions() + 1, eigenvalues.numDimensions() ); + for ( int d = 0; d < img.numDimensions(); ++d ) + Assert.assertEquals( img.dimension( d ), eigenvalues.dimension( d ) ); + + Assert.assertEquals( img.numDimensions(), eigenvalues.dimension( img.numDimensions() ) ); + + + for ( Cursor< DoubleType > r = ArrayImgs.doubles( eigenvaluesRef1, size, size ).cursor(), c = Views.hyperSlice( eigenvalues, img.numDimensions(), 0 ).cursor(); r.hasNext(); ) + Assert.assertEquals( r.next().get(), c.next().get(), 1.0e-20 ); + + for ( Cursor< DoubleType > r = ArrayImgs.doubles( eigenvaluesRef2, size, size ).cursor(), c = Views.hyperSlice( eigenvalues, img.numDimensions(), 1 ).cursor(); r.hasNext(); ) + Assert.assertEquals( r.next().get(), c.next().get(), 1.0e-20 ); + + + } + +} diff --git a/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixTest.java b/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixTest.java new file mode 100644 index 000000000..145cc697f --- /dev/null +++ b/src/test/java/net/imglib2/algorithm/hessian/HessianMatrixTest.java @@ -0,0 +1,273 @@ +/* + * #%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.hessian; + +import java.util.Arrays; +import java.util.stream.IntStream; + +import org.junit.Assert; +import org.junit.Test; + +import net.imglib2.Cursor; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.gauss3.Gauss3; +import net.imglib2.algorithm.gradient.HessianMatrix; +import net.imglib2.exception.IncompatibleTypeException; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.DoubleArray; +import net.imglib2.outofbounds.OutOfBoundsBorderFactory; +import net.imglib2.type.numeric.real.DoubleType; +import net.imglib2.util.Pair; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +public class HessianMatrixTest +{ + + public void test( + final RandomAccessibleInterval< DoubleType > img, + final RandomAccessibleInterval< DoubleType > hessian, + final double[][] refs, + final int size ) + { + Assert.assertEquals( img.numDimensions() + 1, hessian.numDimensions() ); + for ( int d = 0; d < img.numDimensions(); ++d ) + Assert.assertEquals( img.dimension( d ), hessian.dimension( d ) ); + + Assert.assertEquals( img.numDimensions() * ( img.numDimensions() + 1 ) / 2, hessian.dimension( img.numDimensions() ) ); + + for ( int i = 0; i < hessian.dimension( img.numDimensions() ); ++i ) + for ( Cursor< DoubleType > r = ArrayImgs.doubles( refs[ i ], size, size ).cursor(), c = Views.hyperSlice( hessian, img.numDimensions(), i ).cursor(); r.hasNext(); ) + Assert.assertEquals( r.next().get(), c.next().get(), 1e-20 ); + } + + @Test + public void test2D() throws IncompatibleTypeException + { + final int size = 5; + + final double[] imgArray = new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 4.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + final double[] dxxRef = new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 1.0, 0.0, -2.0, 0.0, 1.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + final double[] dxyRef = new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0, -1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, -1.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + final double[] dyyRef = new double[] { + 0.0, 0.0, 1.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, -2.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0 + }; + + final double[][] refs = { + dxxRef, + dxyRef, + dyyRef + }; + + final ArrayImg< DoubleType, DoubleArray > img = ArrayImgs.doubles( imgArray, size, size ); + + final double sigma = 0.1; + + final ArrayImg< DoubleType, DoubleArray > gaussian = ArrayImgs.doubles( size, size ); + Gauss3.gauss( sigma, Views.extendBorder( img ), gaussian ); + + final RandomAccessibleInterval< DoubleType > hessian = + HessianMatrix.calculateMatrix( Views.extendBorder( gaussian ), ArrayImgs.doubles( size, size, 2 ), ArrayImgs.doubles( size, size, 3 ), new OutOfBoundsBorderFactory<>() ); + + test( img, hessian, refs, size ); + + } + + @Test + public void test3D() throws IncompatibleTypeException + { + final int size = 3; + + final double[] imgArray = new double[] { + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 4.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0 + }; + + final double[] dxxRef = new double[] { + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + -1.0, -2.0, -1.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0 + }; + + final double[] dxyRef = new double[] { + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 1.0, 0.0, -1.0, + 0.0, 0.0, 0.0, + -1.0, 0.0, 1.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0 + }; + + final double[] dxzRef = new double[] { + 0.0, 0.0, 0.0, + 1.0, 0.0, -1.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + -1.0, 0.0, 1.0, + 0.0, 0.0, 0.0 + }; + + final double[] dyyRef = new double[] { + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, -1.0, 0.0, + 0.0, -2.0, 0.0, + 0.0, -1.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0 + }; + + final double[] dyzRef = new double[] { + 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, -1.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, -1.0, 0.0, + 0.0, 0.0, 0.0, + 0.0, 1.0, 0.0 + }; + + final double[] dzzRef = new double[] { + 0.0, 0.0, 0.0, + 0.0, -1.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, -2.0, 0.0, + 0.0, 0.0, 0.0, + + 0.0, 0.0, 0.0, + 0.0, -1.0, 0.0, + 0.0, 0.0, 0.0 + }; + + final double[][] refs = { + dxxRef, + dxyRef, + dxzRef, + dyyRef, + dyzRef, + dzzRef + }; + + final ArrayImg< DoubleType, DoubleArray > img = ArrayImgs.doubles( imgArray, size, size, size ); + + final double sigma = 0.1; + + final ArrayImg< DoubleType, DoubleArray > gaussian = ArrayImgs.doubles( size, size, size ); + Gauss3.gauss( sigma, Views.extendBorder( img ), gaussian ); + + final RandomAccessibleInterval< DoubleType > hessian = + HessianMatrix.calculateMatrix( Views.extendBorder( gaussian ), ArrayImgs.doubles( size, size, size, 3 ), ArrayImgs.doubles( size, size, size, 6 ), new OutOfBoundsBorderFactory<>() ); + + test( img, hessian, refs, size ); + } + + @Test + public void testScaling() + { + final int nDim = 3; + final double[] data = Arrays.stream( new double[ nDim * ( nDim + 1 ) / 2 ] ).map( d -> 1.0d ).toArray(); + final double[] sigma = IntStream.range( 0, nDim ).mapToDouble( i -> i + 1 ).toArray(); + final ArrayImg< DoubleType, DoubleArray > hessian = ArrayImgs.doubles( data, 1, 1, 1, data.length ); + final ArrayImg< DoubleType, DoubleArray > ref = ArrayImgs.doubles( new double[] { 1.0, 2.0, 3.0, 4.0, 6.0, 9.0 }, 1, 1, 1, data.length ); + final IntervalView< DoubleType > scaled = HessianMatrix.scaleHessianMatrix( hessian, sigma ); + for ( final Pair< DoubleType, DoubleType > p : Views.interval( Views.pair( ref, scaled ), ref ) ) + Assert.assertEquals( p.getA().get(), p.getB().get(), 0.0d ); + + } + +} diff --git a/src/test/java/net/imglib2/algorithm/hessian/ScaleAsFunctionOfPositionTest.java b/src/test/java/net/imglib2/algorithm/hessian/ScaleAsFunctionOfPositionTest.java new file mode 100644 index 000000000..96e02654b --- /dev/null +++ b/src/test/java/net/imglib2/algorithm/hessian/ScaleAsFunctionOfPositionTest.java @@ -0,0 +1,80 @@ +/* + * #%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.hessian; + +import java.util.stream.IntStream; + +import org.junit.Assert; +import org.junit.Test; + +import net.imglib2.Cursor; +import net.imglib2.algorithm.gradient.ScaleAsFunctionOfPosition; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.IntArray; +import net.imglib2.type.numeric.integer.IntType; +import net.imglib2.util.Intervals; +import net.imglib2.view.Views; + +public class ScaleAsFunctionOfPositionTest +{ + + long[] dim = new long[] { 2, 3, 4 }; + + int[] data = IntStream.range( 0, ( int ) Intervals.numElements( dim ) ).map( i -> 1 ).toArray(); + + ArrayImg< IntType, IntArray > img = ArrayImgs.ints( data, dim ); + + @Test + public void test() + { + final ScaleAsFunctionOfPosition< IntType > scaled = new ScaleAsFunctionOfPosition<>( img, l -> { + double s = l.getDoublePosition( 0 ); + for ( int d = 1; d < l.numDimensions(); ++d ) + s = Math.max( l.getDoublePosition( d ), s ); + return s; + } ); + + for ( final Cursor< IntType > c = Views.interval( scaled, img ).cursor(); c.hasNext(); ) { + final int v = c.next().get(); + int m = c.getIntPosition( 0 ); + for ( int d = 1; d < c.numDimensions(); ++d ) + m = Math.max( c.getIntPosition( d ), m ); + Assert.assertEquals( m, v ); + } + + } + +}