Skip to content

Commit c3232dd

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

File tree

3 files changed

+239
-96
lines changed

3 files changed

+239
-96
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/convolution/fast_gauss/FastGaussTest.java‎

Lines changed: 0 additions & 86 deletions
This file was deleted.

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

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

0 commit comments

Comments
 (0)