Skip to content

Commit e38efd5

Browse files
maarztmaarzt
committed
Smoothen the edge of the truncated Gaussian kernel in Gauss3
1 parent 4de0236 commit e38efd5

File tree

2 files changed

+248
-10
lines changed

2 files changed

+248
-10
lines changed

‎src/main/java/net/imglib2/algorithm/gauss3/Gauss3.java

Lines changed: 70 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -226,26 +226,86 @@ public static int[] halfkernelsizes( final double[] sigma )
226226
return size;
227227
}
228228

229+
/**
230+
* Returns a gaussian half kernel with the given sigma and size.
231+
* <p>
232+
* The edges are smoothed by a second degree polynomial.
233+
* This improves the first derivative and the fourier spectrum
234+
* of the gaussian kernel.
235+
*/
229236
public static double[] halfkernel( final double sigma, final int size, final boolean normalize )
230237
{
231-
final double two_sq_sigma = 2 * sigma * sigma;
238+
final double two_sq_sigma = 2 * square( sigma );
232239
final double[] kernel = new double[ size ];
233240

234241
kernel[ 0 ] = 1;
235242
for ( int x = 1; x < size; ++x )
236-
kernel[ x ] = Math.exp( -( x * x ) / two_sq_sigma );
243+
kernel[ x ] = Math.exp( -square( x ) / two_sq_sigma );
244+
245+
smoothEdge( kernel );
237246

238247
if ( normalize )
239-
{
240-
double sum = 0.5;
241-
for ( int x = 1; x < size; ++x )
242-
sum += kernel[ x ];
243-
sum *= 2;
248+
normalizeHalfkernel( kernel );
244249

245-
for ( int x = 0; x < size; ++x )
246-
kernel[ x ] /= sum;
250+
return kernel;
251+
}
252+
253+
/**
254+
* This method smooths the truncated end of the gaussian kernel.
255+
* The code is taken from ImageJ1 "Gaussian Blur ...".
256+
* <p>
257+
* Detailed explanation:
258+
* <p>
259+
* The values kernel[x] for r < x < L are replaced by the values
260+
* of a polynomial p(x) = slope * (x - L) ^ 2.
261+
* Where L equals kernel.length. And the "slope" and "r" are chosen
262+
* such that there is a smooth transition between the kernel and
263+
* the polynomial. Thus their value and first derivative
264+
* match for x = r: p(r) = kernel[r] and p'(r) = kernel'[r].
265+
*/
266+
private static void smoothEdge( double[] kernel )
267+
{
268+
int L = kernel.length;
269+
double slope = Double.MAX_VALUE;
270+
int r = L;
271+
while ( r > L / 2 )
272+
{
273+
r--;
274+
double a = kernel[ r ] / square( L - r );
275+
if ( a < slope )
276+
slope = a;
277+
else
278+
{
279+
r++;
280+
break;
281+
}
247282
}
283+
for ( int x = r + 1; x < L; x++ )
284+
kernel[ x ] = slope * square( L - x );
285+
}
248286

249-
return kernel;
287+
/**
288+
* Normalizes a half kernel such that the values sum up to 1.
289+
*/
290+
private static void normalizeHalfkernel( double[] kernel )
291+
{
292+
double sum = 0.5 * kernel[ 0 ];
293+
for ( int x = 1; x < kernel.length; ++x )
294+
sum += kernel[ x ];
295+
sum *= 2;
296+
multiply( kernel, 1 / sum );
297+
}
298+
299+
// -- Helper methods --
300+
301+
private static double square( double x )
302+
{
303+
return x * x;
304+
}
305+
306+
private static void multiply( double[] values, double factor )
307+
{
308+
for ( int x = 0; x < values.length; ++x )
309+
values[ x ] *= factor;
250310
}
251311
}

‎src/test/java/net/imglib2/algorithm/gauss3/Gauss3Test.java

Lines changed: 178 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,178 @@
1+
package net.imglib2.algorithm.gauss3;
2+
3+
import net.imglib2.IterableInterval;
4+
import net.imglib2.Localizable;
5+
import net.imglib2.RandomAccessibleInterval;
6+
import net.imglib2.algorithm.convolution.fast_gauss.FastGauss;
7+
import net.imglib2.algorithm.convolution.kernel.Kernel1D;
8+
import net.imglib2.algorithm.convolution.kernel.SeparableKernelConvolution;
9+
import net.imglib2.algorithm.gradient.PartialDerivative;
10+
import net.imglib2.converter.Converters;
11+
import net.imglib2.img.Img;
12+
import net.imglib2.img.array.ArrayImgFactory;
13+
import net.imglib2.loops.LoopBuilder;
14+
import net.imglib2.type.NativeType;
15+
import net.imglib2.type.numeric.RealType;
16+
import net.imglib2.type.numeric.real.DoubleType;
17+
import net.imglib2.util.Intervals;
18+
import net.imglib2.util.Localizables;
19+
import net.imglib2.util.Util;
20+
import net.imglib2.view.Views;
21+
import org.junit.Ignore;
22+
import org.junit.Test;
23+
24+
import java.util.function.BiFunction;
25+
26+
import static org.junit.Assert.fail;
27+
28+
public class Gauss3Test< T extends RealType< T > & NativeType< T > >
29+
{
30+
31+
private T type = ( T ) new DoubleType();
32+
33+
private double sigma = 4;
34+
35+
private long center = (long) ( 12 * sigma );
36+
37+
private long width = center * 2;
38+
39+
private RandomAccessibleInterval< T > input = scaleAndAddOffset( dirac() );
40+
41+
private RandomAccessibleInterval< T > expected = scaleAndAddOffset( idealGaussian( sigma ) );
42+
43+
@Test
44+
public void testGauss3()
45+
{
46+
RandomAccessibleInterval< T > result = createEmptyImage();
47+
Gauss3.gauss( sigma, Views.extendBorder( input ), result );
48+
assertImagesEqual( 40, subtractOffset( expected ), subtractOffset( result ) );
49+
assertImagesEqual( 35, deriveX( expected ), deriveX( result ) );
50+
assertImagesEqual( 24, secondDerivativeX( expected ), secondDerivativeX( result ) );
51+
}
52+
53+
@Ignore( "The FastGauss currently deals poorly with an offset in the image." )
54+
@Test
55+
public void testFastGauss()
56+
{
57+
RandomAccessibleInterval< T > result = createEmptyImage();
58+
FastGauss.convolve( sigma, Views.extendBorder( input ), result );
59+
assertImagesEqual( 50, subtractOffset( expected ), subtractOffset( result ) );
60+
assertImagesEqual( 45, deriveX( expected ), deriveX( result ) );
61+
assertImagesEqual( 35, secondDerivativeX( expected ), secondDerivativeX( result ) );
62+
}
63+
64+
@Test
65+
public void testFastGaussWithoutOffset()
66+
{
67+
// NB: This test shows, that FastGauss promises to have a higher acurazy than
68+
// Gauss3, if the problem demonstrated by {@link #testFastGauss} is fixed.
69+
RandomAccessibleInterval< T > result = createEmptyImage();
70+
FastGauss.convolve( sigma, Views.extendBorder( subtractOffset( input ) ), result );
71+
assertImagesEqual( 50, subtractOffset( expected ), result );
72+
assertImagesEqual( 45, deriveX( expected ), deriveX( result ) );
73+
assertImagesEqual( 35, secondDerivativeX( expected ), secondDerivativeX( result ) );
74+
}
75+
76+
// -- Helper methods --
77+
78+
private RandomAccessibleInterval< T > subtractOffset( RandomAccessibleInterval< ? extends RealType< ? > > image )
79+
{
80+
return Converters.convert( image, ( i, o ) -> o.setReal( i.getRealDouble() - 80 ), type );
81+
}
82+
83+
private RandomAccessibleInterval< T > scaleAndAddOffset( RandomAccessibleInterval< T > dirac )
84+
{
85+
LoopBuilder.setImages( dirac ).forEachPixel( pixel -> pixel.setReal( 20 * pixel.getRealDouble() + 80 ) );
86+
return dirac;
87+
}
88+
89+
private RandomAccessibleInterval< T > idealGaussian( double sigma )
90+
{
91+
return createImage( ( x, y ) -> gauss( sigma, x ) * gauss( sigma, y ) );
92+
}
93+
94+
private double gauss( double sigma, double x )
95+
{
96+
double a = 1. / Math.sqrt( 2 * Math.PI * Math.pow( sigma, 2 ) );
97+
double b = -0.5 / Math.pow( sigma, 2 );
98+
return a * Math.exp( b * Math.pow( x, 2 ) );
99+
}
100+
101+
private RandomAccessibleInterval< T > dirac()
102+
{
103+
return createImage( ( x, y ) -> ( x == 0 ) && ( y == 0 ) ? 1. : 0. );
104+
}
105+
106+
private RandomAccessibleInterval< T > createImage( BiFunction< Long, Long, Double > content )
107+
{
108+
Img< T > image = createEmptyImage();
109+
RandomAccessibleInterval< Localizable > positions = Views.interval( Localizables.randomAccessible( image.numDimensions() ), image );
110+
LoopBuilder.setImages( positions, image ).forEachPixel( ( p, pixel ) -> {
111+
long x = p.getLongPosition( 0 ) - center;
112+
long y = p.getLongPosition( 1 ) - center;
113+
pixel.setReal( content.apply( x, y ) );
114+
} );
115+
return image;
116+
}
117+
118+
private Img< T > createEmptyImage()
119+
{
120+
return new ArrayImgFactory<>( type ).create( width, width );
121+
}
122+
123+
private < T extends RealType< T > & NativeType< T > > RandomAccessibleInterval< T > deriveX( RandomAccessibleInterval< T > input )
124+
{
125+
Img< T > result = new ArrayImgFactory<>( Util.getTypeFromInterval( input ) ).create( Intervals.dimensionsAsLongArray( input ) );
126+
PartialDerivative.gradientCentralDifference( Views.extendBorder( input ), result, 0 );
127+
return result;
128+
}
129+
130+
private RandomAccessibleInterval< T > secondDerivativeX( RandomAccessibleInterval< ? extends RealType< ? > > input )
131+
{
132+
Img< T > result = createEmptyImage();
133+
SeparableKernelConvolution.convolution1d( Kernel1D.centralAsymmetric( 1, -2, 1 ), 0 )
134+
.process( Views.extendBorder( input ), result );
135+
return result;
136+
}
137+
138+
private void assertImagesEqual( int expectedSnr, RandomAccessibleInterval< ? extends RealType< ? > > a, RandomAccessibleInterval< T > b )
139+
{
140+
double actualSnr = snr( a, b );
141+
if ( expectedSnr > actualSnr )
142+
fail( "The SNR is lower than expected, expected: " + expectedSnr + " dB actual: " + actualSnr + " dB" );
143+
}
144+
145+
private static double snr( RandomAccessibleInterval< ? extends RealType< ? > > expected,
146+
RandomAccessibleInterval< ? extends RealType< ? > > actual )
147+
{
148+
double signal = meanSquaredSum( expected );
149+
double noise = meanSquaredSum( subtract( actual, expected ) );
150+
if ( signal == 0.0 )
151+
return Float.NEGATIVE_INFINITY;
152+
return 10 * ( Math.log10( signal / noise ) );
153+
}
154+
155+
private static RandomAccessibleInterval< DoubleType > subtract(
156+
RandomAccessibleInterval< ? extends RealType > a,
157+
RandomAccessibleInterval< ? extends RealType > b )
158+
{
159+
return Views.interval( Converters.convert(
160+
Views.pair( a, b ),
161+
( pair, out ) -> out.setReal( pair.getA().getRealDouble() - pair.getB().getRealDouble() ),
162+
new DoubleType() ), a );
163+
}
164+
165+
private static double meanSquaredSum( RandomAccessibleInterval< ? extends RealType< ? > > image )
166+
{
167+
double sum = 0;
168+
IterableInterval< ? extends RealType< ? > > iterable = Views.iterable( image );
169+
for ( RealType< ? > pixel : iterable )
170+
sum += square( pixel.getRealDouble() );
171+
return sum / iterable.size();
172+
}
173+
174+
private static double square( double value )
175+
{
176+
return value * value;
177+
}
178+
}

0 commit comments

Comments
 (0)