Skip to content

Commit d49d200

Browse files
tpietzschtpietzsch
committed
Add symetric convolution UnaryBlockOperators
1 parent b1ae0d9 commit d49d200

File tree

2 files changed

+594
-0
lines changed

2 files changed

+594
-0
lines changed

‎src/main/java/net/imglib2/algorithm/blocks/convolve/Convolve.java

+182
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,182 @@
1+
package net.imglib2.algorithm.blocks.convolve;
2+
3+
import static net.imglib2.algorithm.blocks.ClampType.NONE;
4+
import static net.imglib2.type.PrimitiveType.FLOAT;
5+
6+
import java.util.function.Function;
7+
8+
import net.imglib2.algorithm.blocks.BlockSupplier;
9+
import net.imglib2.algorithm.blocks.ClampType;
10+
import net.imglib2.algorithm.blocks.ComputationType;
11+
import net.imglib2.algorithm.blocks.DefaultUnaryBlockOperator;
12+
import net.imglib2.algorithm.blocks.UnaryBlockOperator;
13+
import net.imglib2.algorithm.blocks.convolve.ConvolveProcessors.ConvolveDouble;
14+
import net.imglib2.algorithm.blocks.convolve.ConvolveProcessors.ConvolveFloat;
15+
import net.imglib2.algorithm.convolution.kernel.Kernel1D;
16+
import net.imglib2.algorithm.gauss3.Gauss3;
17+
import net.imglib2.type.NativeType;
18+
import net.imglib2.type.PrimitiveType;
19+
import net.imglib2.type.numeric.real.DoubleType;
20+
import net.imglib2.type.numeric.real.FloatType;
21+
import net.imglib2.util.Util;
22+
23+
/**
24+
* Separable convolution.
25+
* <p>
26+
* Supported types are {@code UnsignedByteType}, {@code UnsignedShortType},
27+
* {@code UnsignedIntType}, {@code ByteType}, {@code ShortType}, {@code
28+
* IntType}, {@code LongType}, {@code FloatType}, {@code DoubleType}).
29+
* <p>
30+
* For {@code T} other than {@code DoubleType} or {@code FloatType}, the input
31+
* will be converted to float/double for computation and the result converted
32+
* back to {@code T}. To avoid unnecessary conversions, if you want the result
33+
* as {@code FloatType} then you should explicitly convert to {@code FloatType}
34+
* <em>before</em> applying the convolve operator.
35+
* This code:
36+
* <pre>{@code
37+
* RandomAccessible< UnsignedByteType > input;
38+
* BlockSupplier< FloatType > convolved = BlockSupplier.of( input )
39+
* .andThen( Convert.convert( new FloatType() ) )
40+
* .andThen( Convolve.convolve( kernels ) );
41+
* }</pre>
42+
* avoids loss of precision and is more efficient than
43+
* <pre>{@code
44+
* RandomAccessible< UnsignedByteType > input;
45+
* BlockSupplier< FloatType > convolved = BlockSupplier.of( input )
46+
* .andThen( Convolve.convolve( kernels ) )
47+
* .andThen( Convert.convert( new FloatType() ) );
48+
* }</pre>
49+
*
50+
*/
51+
public class Convolve
52+
{
53+
54+
/**
55+
* Convolve blocks of the standard ImgLib2 {@code RealType}s with a Gaussian kernel.
56+
* <p>
57+
* Precision for intermediate values is chosen as to represent the
58+
* input/output type without loss of precision. That is, {@code FLOAT} for
59+
* u8, i8, u16, i16, i32, f32, and otherwise {@code DOUBLE} for u32, i64,
60+
* f64.
61+
* <p>
62+
* The returned factory function creates an operator matching the
63+
* type and dimensionality of a given input {@code BlockSupplier<T>}.
64+
*
65+
* @param sigma
66+
* sigmas in each dimension. if the image has fewer or more dimensions
67+
* than values given, values will be truncated or the final value
68+
* repeated as necessary.
69+
* @param <T>
70+
* the input/output type
71+
*
72+
* @return factory for {@code UnaryBlockOperator} to convolve blocks of type {@code T}
73+
*/
74+
public static < T extends NativeType< T > >
75+
Function< BlockSupplier< T >, UnaryBlockOperator< T, T > > gauss( final double... sigma )
76+
{
77+
return gauss( ComputationType.AUTO, sigma );
78+
}
79+
80+
/**
81+
* Convolve blocks of the standard ImgLib2 {@code RealType}s with a Gaussian kernel.
82+
* <p>
83+
* The returned factory function creates an operator matching the
84+
* type and dimensionality of a given input {@code BlockSupplier<T>}.
85+
*
86+
* @param computationType
87+
* specifies in which precision intermediate values should be
88+
* computed. For {@code AUTO}, the type that can represent the
89+
* input/output type without loss of precision is picked. That is,
90+
* {@code FLOAT} for u8, i8, u16, i16, i32, f32, and otherwise {@code
91+
* DOUBLE} for u32, i64, f64.
92+
* @param sigma
93+
* sigmas in each dimension. if the image has fewer or more dimensions
94+
* than values given, values will be truncated or the final value
95+
* repeated as necessary.
96+
* @param <T>
97+
* the input/output type
98+
*
99+
* @return factory for {@code UnaryBlockOperator} to convolve blocks of type {@code T}
100+
*/
101+
public static < T extends NativeType< T > >
102+
Function< BlockSupplier< T >, UnaryBlockOperator< T, T > > gauss( final ComputationType computationType, final double... sigma )
103+
{
104+
return s -> {
105+
final T type = s.getType();
106+
final int n = s.numDimensions();
107+
return createOperator( type, computationType, ClampType.NONE, gaussKernels( Util.expandArray( sigma, n ) ) );
108+
};
109+
}
110+
111+
static Kernel1D[] gaussKernels( final double[] sigma )
112+
{
113+
final Kernel1D[] kernels = new Kernel1D[ sigma.length ];
114+
for ( int d = 0; d < sigma.length; d++ )
115+
if ( sigma[ d ] > 0 )
116+
kernels[ d ] = Kernel1D.symmetric( Gauss3.halfkernel( sigma[ d ] ) );
117+
return kernels;
118+
}
119+
120+
/**
121+
* Create a {@code UnaryBlockOperator} to convolve with the given {@code kernels}.
122+
* {@code kernels.length} must match the dimensionality of the input images.
123+
* If {@code kernels[d]==null}, the convolution for dimension {@code d} is
124+
* skipped (equivalent to convolution with the kernel {@code {1}}).
125+
* <p>
126+
* Supported types are {@code UnsignedByteType}, {@code UnsignedShortType},
127+
* {@code UnsignedIntType}, {@code ByteType}, {@code ShortType}, {@code
128+
* IntType}, {@code LongType}, {@code FloatType}, {@code DoubleType}).
129+
*
130+
* @param type
131+
* instance of the input type
132+
* @param computationType
133+
* specifies in which precision intermediate values should be
134+
* computed. For {@code AUTO}, the type that can represent the
135+
* input/output type without loss of precision is picked. That is,
136+
* {@code FLOAT} for u8, i8, u16, i16, i32, f32, and otherwise {@code
137+
* DOUBLE} for u32, i64, f64.
138+
* @param kernels
139+
* kernel to apply in each dimension
140+
* @param <T>
141+
* the input/output type
142+
*
143+
* @return {@code UnaryBlockOperator} to downsample blocks of type {@code T}
144+
*/
145+
public static < T extends NativeType< T > >
146+
UnaryBlockOperator< T, T > createOperator( final T type, final ComputationType computationType, final ClampType clampType, final Kernel1D[] kernels )
147+
{
148+
final boolean processAsFloat;
149+
switch ( computationType )
150+
{
151+
case FLOAT:
152+
processAsFloat = true;
153+
break;
154+
case DOUBLE:
155+
processAsFloat = false;
156+
break;
157+
default:
158+
case AUTO:
159+
final PrimitiveType pt = type.getNativeTypeFactory().getPrimitiveType();
160+
processAsFloat = pt.equals( FLOAT ) || pt.getByteCount() < FLOAT.getByteCount();
161+
break;
162+
}
163+
final UnaryBlockOperator< ?, ? > op = processAsFloat
164+
? convolveFloat( kernels )
165+
: convolveDouble( kernels );
166+
return op.adaptSourceType( type, NONE ).adaptTargetType( type, clampType );
167+
}
168+
169+
private static UnaryBlockOperator< FloatType, FloatType > convolveFloat( final Kernel1D[] kernels )
170+
{
171+
final FloatType type = new FloatType();
172+
final int n = kernels.length;
173+
return new DefaultUnaryBlockOperator<>( type, type, n, n, new ConvolveFloat( kernels ) );
174+
}
175+
176+
private static UnaryBlockOperator< DoubleType, DoubleType > convolveDouble( final Kernel1D[] kernels )
177+
{
178+
final DoubleType type = new DoubleType();
179+
final int n = kernels.length;
180+
return new DefaultUnaryBlockOperator<>( type, type, n, n, new ConvolveDouble( kernels ) );
181+
}
182+
}

0 commit comments

Comments
 (0)