Skip to content

Commit c51c6e9

Browse files
tpietzschtpietzsch
authored
Merge pull request #41 from hanslovsky/hessian-eigenvalues-rebase
Add HessianMatrix and TensorEigenValues
2 parents f5332f6 + caef15f commit c51c6e9

17 files changed

+2153
-4
lines changed

‎pom.xml‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -229,6 +229,11 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
229229
<artifactId>trove4j</artifactId>
230230
<version>3.0.3</version>
231231
</dependency>
232+
<dependency>
233+
<groupId>org.ojalgo</groupId>
234+
<artifactId>ojalgo</artifactId>
235+
<version>43.0</version>
236+
</dependency>
232237

233238
<!-- Test dependencies -->
234239
<dependency>

‎src/main/java/net/imglib2/algorithm/gradient/HessianMatrix.java‎

Lines changed: 389 additions & 0 deletions
Large diffs are not rendered by default.

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

Lines changed: 90 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -11,13 +11,13 @@
1111
* %%
1212
* Redistribution and use in source and binary forms, with or without
1313
* modification, are permitted provided that the following conditions are met:
14-
*
14+
*
1515
* 1. Redistributions of source code must retain the above copyright notice,
1616
* this list of conditions and the following disclaimer.
1717
* 2. Redistributions in binary form must reproduce the above copyright notice,
1818
* this list of conditions and the following disclaimer in the documentation
1919
* and/or other materials provided with the distribution.
20-
*
20+
*
2121
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
2222
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
2323
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
@@ -34,20 +34,34 @@
3434

3535
package net.imglib2.algorithm.gradient;
3636

37+
import java.util.ArrayList;
38+
import java.util.List;
39+
import java.util.concurrent.Callable;
40+
import java.util.concurrent.ExecutionException;
41+
import java.util.concurrent.ExecutorService;
42+
import java.util.concurrent.Future;
43+
3744
import net.imglib2.Cursor;
45+
import net.imglib2.FinalInterval;
3846
import net.imglib2.RandomAccess;
3947
import net.imglib2.RandomAccessible;
4048
import net.imglib2.RandomAccessibleInterval;
4149
import net.imglib2.type.numeric.NumericType;
4250
import net.imglib2.util.Intervals;
51+
import net.imglib2.view.IntervalView;
4352
import net.imglib2.view.Views;
4453

54+
/**
55+
* @author Tobias Pietzsch
56+
* @author Philipp Hanslovsky
57+
*
58+
*/
4559
public class PartialDerivative
4660
{
4761
// nice version...
4862
/**
4963
* Compute the partial derivative of source in a particular dimension.
50-
*
64+
*
5165
* @param source
5266
* source image, has to provide valid data in the interval of the
5367
* gradient image plus a one pixel border in dimension.
@@ -69,10 +83,82 @@ public static < T extends NumericType< T > > void gradientCentralDifference2( fi
6983
}
7084
}
7185

86+
// parallel version...
87+
/**
88+
* Compute the partial derivative of source in a particular dimension.
89+
*
90+
* @param source
91+
* source image, has to provide valid data in the interval of the
92+
* gradient image plus a one pixel border in dimension.
93+
* @param gradient
94+
* output image
95+
* @param dimension
96+
* along which dimension the partial derivatives are computed
97+
* @param nTasks
98+
* Number of tasks for gradient computation.
99+
* @param es
100+
* {@link ExecutorService} providing workers for gradient
101+
* computation. Service is managed (created, shutdown) by caller.
102+
*/
103+
public static < T extends NumericType< T > > void gradientCentralDifferenceParallel(
104+
final RandomAccessible< T > source,
105+
final RandomAccessibleInterval< T > gradient,
106+
final int dimension,
107+
final int nTasks,
108+
final ExecutorService es ) throws InterruptedException, ExecutionException
109+
{
110+
final int nDim = source.numDimensions();
111+
if ( nDim < 2 )
112+
{
113+
gradientCentralDifference( source, gradient, dimension );
114+
return;
115+
}
116+
117+
long dimensionMax = Long.MIN_VALUE;
118+
int dimensionArgMax = -1;
119+
120+
for ( int d = 0; d < nDim; ++d )
121+
{
122+
final long size = gradient.dimension( d );
123+
if ( d != dimension && size > dimensionMax )
124+
{
125+
dimensionMax = size;
126+
dimensionArgMax = d;
127+
}
128+
}
129+
130+
final long stepSize = Math.max( dimensionMax / nTasks, 1 );
131+
final long stepSizeMinusOne = stepSize - 1;
132+
final long min = gradient.min( dimensionArgMax );
133+
final long max = gradient.max( dimensionArgMax );
134+
135+
final ArrayList< Callable< Void > > tasks = new ArrayList<>();
136+
for ( long currentMin = min, minZeroBase = 0; minZeroBase < dimensionMax; currentMin += stepSize, minZeroBase += stepSize )
137+
{
138+
final long currentMax = Math.min( currentMin + stepSizeMinusOne, max );
139+
final long[] mins = new long[ nDim ];
140+
final long[] maxs = new long[ nDim ];
141+
gradient.min( mins );
142+
gradient.max( maxs );
143+
mins[ dimensionArgMax ] = currentMin;
144+
maxs[ dimensionArgMax ] = currentMax;
145+
final IntervalView< T > currentInterval = Views.interval( gradient, new FinalInterval( mins, maxs ) );
146+
tasks.add( () -> {
147+
gradientCentralDifference( source, currentInterval, dimension );
148+
return null;
149+
} );
150+
}
151+
152+
final List< Future< Void > > futures = es.invokeAll( tasks );
153+
154+
for ( final Future< Void > f : futures )
155+
f.get();
156+
}
157+
72158
// fast version
73159
/**
74160
* Compute the partial derivative of source in a particular dimension.
75-
*
161+
*
76162
* @param source
77163
* source image, has to provide valid data in the interval of the
78164
* gradient image plus a one pixel border in dimension.

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

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
/*
2+
* #%L
3+
* ImgLib2: a general-purpose, multidimensional image processing library.
4+
* %%
5+
* Copyright (C) 2009 - 2016 Tobias Pietzsch, Stephan Preibisch, Stephan Saalfeld,
6+
* John Bogovic, Albert Cardona, Barry DeZonia, Christian Dietz, Jan Funke,
7+
* Aivar Grislis, Jonathan Hale, Grant Harris, Stefan Helfrich, Mark Hiner,
8+
* Martin Horn, Steffen Jaensch, Lee Kamentsky, Larry Lindsey, Melissa Linkert,
9+
* Mark Longair, Brian Northan, Nick Perry, Curtis Rueden, Johannes Schindelin,
10+
* Jean-Yves Tinevez and Michael Zinsmaier.
11+
* %%
12+
* Redistribution and use in source and binary forms, with or without
13+
* modification, are permitted provided that the following conditions are met:
14+
*
15+
* 1. Redistributions of source code must retain the above copyright notice,
16+
* this list of conditions and the following disclaimer.
17+
* 2. Redistributions in binary form must reproduce the above copyright notice,
18+
* this list of conditions and the following disclaimer in the documentation
19+
* and/or other materials provided with the distribution.
20+
*
21+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
25+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31+
* POSSIBILITY OF SUCH DAMAGE.
32+
* #L%
33+
*/
34+
35+
package net.imglib2.algorithm.gradient;
36+
37+
import java.util.function.ToDoubleFunction;
38+
39+
import net.imglib2.Interval;
40+
import net.imglib2.Localizable;
41+
import net.imglib2.RandomAccess;
42+
import net.imglib2.RandomAccessible;
43+
import net.imglib2.RandomAccessibleInterval;
44+
import net.imglib2.converter.AbstractConvertedRandomAccess;
45+
import net.imglib2.converter.AbstractConvertedRandomAccessible;
46+
import net.imglib2.type.Type;
47+
import net.imglib2.type.operators.MulFloatingPoint;
48+
49+
/**
50+
*
51+
* @author Philipp Hanslovsky
52+
*
53+
* Multiply the value of a {@link RandomAccessible} depending on the
54+
* position of the access.
55+
*
56+
*/
57+
public class ScaleAsFunctionOfPosition< T extends Type< T > & MulFloatingPoint > extends AbstractConvertedRandomAccessible< T, T >
58+
{
59+
60+
private final ToDoubleFunction< Localizable > scalingFunction;
61+
62+
public ScaleAsFunctionOfPosition( final RandomAccessibleInterval< T > source, final ToDoubleFunction< Localizable > scalingFunction )
63+
{
64+
super( source );
65+
this.scalingFunction = scalingFunction;
66+
}
67+
68+
@Override
69+
public ScaledRandomAccess< T > randomAccess()
70+
{
71+
return new ScaledRandomAccess<>( source.randomAccess(), scalingFunction );
72+
}
73+
74+
@Override
75+
public AbstractConvertedRandomAccess< T, T > randomAccess( final Interval interval )
76+
{
77+
return randomAccess();
78+
}
79+
80+
public static class ScaledRandomAccess< T extends Type< T > & MulFloatingPoint > extends AbstractConvertedRandomAccess< T, T >
81+
{
82+
83+
private final T t;
84+
85+
private final ToDoubleFunction< Localizable > scalingFunction;
86+
87+
public ScaledRandomAccess( final RandomAccess< T > source, final ToDoubleFunction< Localizable > scalingFunction )
88+
{
89+
super( source );
90+
this.t = source.get().createVariable();
91+
this.scalingFunction = scalingFunction;
92+
}
93+
94+
@Override
95+
public T get()
96+
{
97+
t.set( source.get() );
98+
t.mul( scalingFunction.applyAsDouble( source ) );
99+
return t;
100+
}
101+
102+
@Override
103+
public ScaledRandomAccess< T > copy()
104+
{
105+
return new ScaledRandomAccess<>( source.copyRandomAccess(), scalingFunction );
106+
}
107+
108+
}
109+
110+
}

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

Lines changed: 88 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,88 @@
1+
/*
2+
* #%L
3+
* ImgLib2: a general-purpose, multidimensional image processing library.
4+
* %%
5+
* Copyright (C) 2009 - 2016 Tobias Pietzsch, Stephan Preibisch, Stephan Saalfeld,
6+
* John Bogovic, Albert Cardona, Barry DeZonia, Christian Dietz, Jan Funke,
7+
* Aivar Grislis, Jonathan Hale, Grant Harris, Stefan Helfrich, Mark Hiner,
8+
* Martin Horn, Steffen Jaensch, Lee Kamentsky, Larry Lindsey, Melissa Linkert,
9+
* Mark Longair, Brian Northan, Nick Perry, Curtis Rueden, Johannes Schindelin,
10+
* Jean-Yves Tinevez and Michael Zinsmaier.
11+
* %%
12+
* Redistribution and use in source and binary forms, with or without
13+
* modification, are permitted provided that the following conditions are met:
14+
*
15+
* 1. Redistributions of source code must retain the above copyright notice,
16+
* this list of conditions and the following disclaimer.
17+
* 2. Redistributions in binary form must reproduce the above copyright notice,
18+
* this list of conditions and the following disclaimer in the documentation
19+
* and/or other materials provided with the distribution.
20+
*
21+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
25+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31+
* POSSIBILITY OF SUCH DAMAGE.
32+
* #L%
33+
*/
34+
35+
package net.imglib2.algorithm.linalg.eigen;
36+
37+
import net.imglib2.type.numeric.ComplexType;
38+
import net.imglib2.type.numeric.RealType;
39+
import net.imglib2.view.composite.Composite;
40+
41+
/**
42+
*
43+
* Interface for handling different cases, e.g. square, symmetric, or 2
44+
* dimensional tensors.
45+
*
46+
*/
47+
public interface EigenValues< T extends RealType< T >, U extends ComplexType< U > >
48+
{
49+
public void compute( final Composite< T > matrix, final Composite< U > evs );
50+
51+
public default EigenValues< T, U > copy()
52+
{
53+
return this;
54+
}
55+
56+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues1D< T, U > oneDimensional()
57+
{
58+
return new EigenValues1D<>();
59+
}
60+
61+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues2DSquare< T, U > square2D()
62+
{
63+
return new EigenValues2DSquare<>();
64+
}
65+
66+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues2DSymmetric< T, U > symmetric2D()
67+
{
68+
return new EigenValues2DSymmetric<>();
69+
}
70+
71+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValuesSquare< T, U > square( final int nDim )
72+
{
73+
return new EigenValuesSquare<>( nDim );
74+
}
75+
76+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValuesSymmetric< T, U > symmetric( final int nDim )
77+
{
78+
return new EigenValuesSymmetric<>( nDim );
79+
}
80+
81+
public static < T extends RealType< T >, U extends ComplexType< U > > EigenValues< T, U > invalid()
82+
{
83+
return ( m, evs ) -> {
84+
throw new UnsupportedOperationException( "EigenValues not implemented yet!" );
85+
};
86+
}
87+
88+
}

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

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
/*
2+
* #%L
3+
* ImgLib2: a general-purpose, multidimensional image processing library.
4+
* %%
5+
* Copyright (C) 2009 - 2016 Tobias Pietzsch, Stephan Preibisch, Stephan Saalfeld,
6+
* John Bogovic, Albert Cardona, Barry DeZonia, Christian Dietz, Jan Funke,
7+
* Aivar Grislis, Jonathan Hale, Grant Harris, Stefan Helfrich, Mark Hiner,
8+
* Martin Horn, Steffen Jaensch, Lee Kamentsky, Larry Lindsey, Melissa Linkert,
9+
* Mark Longair, Brian Northan, Nick Perry, Curtis Rueden, Johannes Schindelin,
10+
* Jean-Yves Tinevez and Michael Zinsmaier.
11+
* %%
12+
* Redistribution and use in source and binary forms, with or without
13+
* modification, are permitted provided that the following conditions are met:
14+
*
15+
* 1. Redistributions of source code must retain the above copyright notice,
16+
* this list of conditions and the following disclaimer.
17+
* 2. Redistributions in binary form must reproduce the above copyright notice,
18+
* this list of conditions and the following disclaimer in the documentation
19+
* and/or other materials provided with the distribution.
20+
*
21+
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22+
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23+
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
24+
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
25+
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
26+
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
27+
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
28+
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
29+
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
30+
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
31+
* POSSIBILITY OF SUCH DAMAGE.
32+
* #L%
33+
*/
34+
35+
package net.imglib2.algorithm.linalg.eigen;
36+
37+
import net.imglib2.type.numeric.ComplexType;
38+
import net.imglib2.type.numeric.RealType;
39+
import net.imglib2.view.composite.Composite;
40+
41+
public class EigenValues1D< T extends RealType< T >, U extends ComplexType< U > > implements EigenValues< T, U >
42+
{
43+
@Override
44+
public void compute( final Composite< T > tensor, final Composite< U > evs )
45+
{
46+
evs.get( 0 ).setReal( tensor.get( 0 ).getRealDouble() );
47+
}
48+
}

0 commit comments

Comments
 (0)