Skip to content

Commit 14a708f

Browse files
bogovicjbogovicj
authored
bspline interpolation (#86)
* direct on-the-fly bspline interpolator (order 3 only) * bspline interpolation with precomputed coefficients (orders 0-5) * lazy bspline interpolation with precomputed coefficients in a lazy cell image (orders 0-5) TODO check if performance of precomputed coefficients for full image is similar to lazy cell image with block size of full image. If yes, then remove full image implementation because redundant. * add Lazy class for lazy evaluation of op-like BiConsumers (cell loaders) TODO this should be replaced when a vetted similar thing is available in ImgLib2
1 parent 2fe0ce8 commit 14a708f

26 files changed

+4252
-2
lines changed

‎README.md‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,2 +1,3 @@
11
[![](https://travis-ci.org/imglib/imglib2-algorithm.svg?branch=master)](https://travis-ci.org/imglib/imglib2-algorithm)
22

3+
* [Notes and recommendations for bspline interpolation](doc/bsplineNotes.md)

‎doc/bsplineNotes.md‎

Lines changed: 95 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,95 @@
1+
## BSpline notes
2+
3+
There are three kinds of BSpline interpolation.
4+
1. **Precomputed coefficients** (fastest, high memory use)
5+
2. **Lazily precomputed coefficients** (medium speed, medium memory use)
6+
3. **On the fly interpolation** (slowest, low memory use)
7+
8+
### Recommendations
9+
10+
* Generally use precomputed coefficients (1)
11+
* Use Lazily precomputed coefficients (2) for large images
12+
* Use the largest practical block size (128 cubed or 512 squared)
13+
* Smaller block size only for spare, small numbers of interpolations
14+
* Avoid blocks with less than 32 pixels per side
15+
* Use on-the-fly interpolation (3) sparingly
16+
17+
Below we describe each method with code examples.
18+
19+
20+
### 1 Precomputed coefficients
21+
22+
Recommended for small to medium sized images or when memory is not a concern.
23+
24+
The method precomputes BSpline coefficients and is substantially faster if get is called many times, but at the cost of memory - it allocates a RandomAccessibleInterval<T> the same size as the image to be interpolated. T defaults to DoubleType. Use it like this:
25+
26+
```
27+
int splineOrder = ...;
28+
boolean clipping = ...;
29+
factory= new BSplineCoefficientsInterpolatorFactory<>( img, splineOrder, clipping );
30+
realImg = Views.interpolate( img, factory);
31+
```
32+
33+
here the factory constructor computes the coefficients - which is why it needs to be passed. Unfortunately, this means that the img argument to Views.interpolate is meaningless. For example,
34+
35+
```
36+
realImgFoo = Views.interpolate( img, factory);
37+
realImgBar = Views.interpolate( someOtherImg, factory);
38+
```
39+
40+
both result in the same thing. We stuck with this to keep the api consistent. No one should be doing the bottom thing anyway.
41+
Generating the coefficients
42+
43+
`BSplineDecomposition` and `BSplineCoefficientsInterpolator` give more control for those who need it. For example:
44+
45+
```
46+
/*
47+
* Compute the coefficients over an arbitrary interval
48+
*/
49+
BSplineDecomposition<T,S> decomp = new BSplineDecomposition<>( splineOrder,, extendedImg );
50+
long[] min = Intervals.minAsLongArray( interval );
51+
Img<S> coefficientsBase = coefficientFactory.create( interval );
52+
coefficients = Views.translate( coefficientsBase, min );
53+
decomp.accept( coefficients );
54+
55+
/*
56+
* Get a RealRandomAccess that uses the spline coefficients computed above
57+
* and using a custom extension of the coefficients
58+
*/
59+
RealRandomAccess interp = BSplineCoefficientsInterpolator.build(
60+
splineOrder,
61+
Views.extendMirrorSingle( coefficients),
62+
new FloatType());
63+
```
64+
65+
### 2 Lazily precomputed coefficients
66+
67+
Recommended for large images or when memory is a concern.
68+
69+
The third method also precomputes BSpline coefficients, but does so over small blocks of a specified size on-demand, rather than computing all coefficients.
70+
71+
```
72+
int splineOrder = ...;
73+
int[] blockSize = ...
74+
BSplineLazyCoefficientsInterpolatorFactory lazyFactory =
75+
new BSplineLazyCoefficientsInterpolatorFactory<>( img, splineOrder, blockSize );
76+
RealRandomAccessible realImgLazy = Views.interpolate(
77+
Views.extendZero( img ), lazyFactory );
78+
```
79+
80+
We generally recommend using the largest block size practical ( e.g. 128-cubed in 3D or 512 squared in 2d, or larger if possible). Larger block sizes improve accuracy and speed (after precomputation). Small block sizes could be useful when only interpolating at a small, sparse set of points, and the goal is to compute as few coefficients as possible.
81+
82+
### 3 On the fly interpolation
83+
84+
We generally do not recommend this approach since lazily computing coefficients gives
85+
significant speed-ups with a controllable memory cost, for similar accuracy.
86+
The code is most similar to current interpolation methods:
87+
88+
```
89+
int splineOrder = ...;
90+
boolean clipping = ...;
91+
int radius = ...;
92+
realImg = Views.interpolate( img, new BSplineInterpolatorFactory<>( splineOrder, clipping, radius ));
93+
```
94+
95+
This computes the interpolation kernel on-the-fly at every get. The kernel size is defined by radius. clipping=true sets the interpolator to prevent over/underflow, just like the `ClampingNLinearInterpolator`.

‎pom.xml‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -240,6 +240,10 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
240240
<groupId>org.ojalgo</groupId>
241241
<artifactId>ojalgo</artifactId>
242242
</dependency>
243+
<dependency>
244+
<groupId>net.imglib2</groupId>
245+
<artifactId>imglib2-cache</artifactId>
246+
</dependency>
243247

244248
<!-- Test dependencies -->
245249
<dependency>
@@ -250,13 +254,11 @@ Jean-Yves Tinevez and Michael Zinsmaier.</license.copyrightOwners>
250254
<dependency>
251255
<groupId>org.openjdk.jmh</groupId>
252256
<artifactId>jmh-core</artifactId>
253-
<version>1.19</version>
254257
<scope>test</scope>
255258
</dependency>
256259
<dependency>
257260
<groupId>org.openjdk.jmh</groupId>
258261
<artifactId>jmh-generator-annprocess</artifactId>
259-
<version>1.19</version>
260262
<scope>test</scope>
261263
</dependency>
262264
</dependencies>

‎src/main/java/net/imglib2/algorithm/bspline/AbstractBsplineKernel.java‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
/*
2+
* #%L
3+
* ImgLib2: a general-purpose, multidimensional image processing library.
4+
* %%
5+
* Copyright (C) 2009 - 2018 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+
package net.imglib2.algorithm.bspline;
35+
36+
public abstract class AbstractBsplineKernel
37+
{
38+
protected static final double ONESIXTH = 1.0 / 6.0;
39+
protected static final double TWOTHIRDS = 2.0 / 3.0;
40+
41+
public abstract double c0();
42+
43+
public abstract double evaluate( final double x );
44+
45+
public double evaluateNorm( final double x )
46+
{
47+
return c0() * evaluate( x ); // make more efficient (save a multiply)?
48+
}
49+
50+
protected static double powIntPositive( final double base, final int pow )
51+
{
52+
double result = 1;
53+
for ( int i = 0; i < pow; ++i )
54+
{
55+
result *= base;
56+
}
57+
return result;
58+
}
59+
60+
}

‎src/main/java/net/imglib2/algorithm/bspline/BSplineCoefficientsInterpolator.java‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
/*
2+
* #%L
3+
* ImgLib2: a general-purpose, multidimensional image processing library.
4+
* %%
5+
* Copyright (C) 2009 - 2018 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.bspline;
36+
37+
import net.imglib2.Localizable;
38+
import net.imglib2.RandomAccessible;
39+
import net.imglib2.RandomAccessibleInterval;
40+
import net.imglib2.RealInterval;
41+
import net.imglib2.RealRandomAccess;
42+
import net.imglib2.algorithm.neighborhood.GeneralRectangleShape;
43+
import net.imglib2.algorithm.neighborhood.RectangleShape;
44+
import net.imglib2.interpolation.InterpolatorFactory;
45+
import net.imglib2.type.numeric.RealType;
46+
47+
48+
/**
49+
* Performs b-spline interpolation for order up to and including 5.
50+
*
51+
* @author John Bogovic
52+
* @author Stephan Saalfeld
53+
*/
54+
public interface BSplineCoefficientsInterpolator<T extends RealType<T>> extends RealRandomAccess< T >, InterpolatorFactory< T, RandomAccessibleInterval<T> >, Localizable
55+
{
56+
public static <S extends RealType<S>> BSplineCoefficientsInterpolator<S> build(
57+
final int order, final RandomAccessible< S > coefficients, final S type )
58+
{
59+
if( order % 2 == 0 )
60+
return new BSplineCoefficientsInterpolatorEven<S>( order, coefficients, type );
61+
else
62+
return new BSplineCoefficientsInterpolatorOdd<S>( order, coefficients, type );
63+
}
64+
65+
@Override
66+
public BSplineCoefficientsInterpolator<T> copy();
67+
68+
@Override
69+
default RealRandomAccess<T> create( RandomAccessibleInterval<T> f )
70+
{
71+
return copy();
72+
}
73+
74+
@Override
75+
default RealRandomAccess< T > create(RandomAccessibleInterval< T > f, RealInterval interval)
76+
{
77+
return copy();
78+
}
79+
80+
@Override
81+
default RealRandomAccess< T > copyRealRandomAccess()
82+
{
83+
return copy();
84+
}
85+
86+
/**
87+
* Return a RectangleShape
88+
*
89+
* This is necessary because rectangles always have an extents equal to an
90+
* odd number of pixels. Need to round up to the nearest odd number when
91+
* an even number of samples are needed.
92+
*
93+
* @return the rectangle shape
94+
*/
95+
public static RectangleShape shapeFromOrder( int bsplineOrder )
96+
{
97+
assert( bsplineOrder <= 5 && bsplineOrder >= 0 );
98+
99+
switch ( bsplineOrder ){
100+
case 0:
101+
return new RectangleShape( 0, false ); // need one sample - correct
102+
case 1:
103+
return new RectangleShape( 1, false ); // need two samples - round up
104+
case 2:
105+
return new RectangleShape( 1, false ); // need three samples - correct
106+
case 3:
107+
return new GeneralRectangleShape( 4, -1, false ); // need four samples - round up
108+
case 4:
109+
return new RectangleShape( 2, false ); // need five samples - round up
110+
case 5:
111+
return new GeneralRectangleShape( 6, -2, false ); // need four samples - round up
112+
default:
113+
return null;
114+
}
115+
}
116+
117+
}

0 commit comments

Comments
 (0)