Differences (#44)

* Add forward and backward difference.

* PartialDerivative: add tests

* POM: bump parent to pom-imglib2 10.0.2

This makes the class LoopBuilder available.

* PartialDerivative: use LoopBuilder to get clean code

* Add formulae for finite diffierence methods

Addresses comments in #44
Also adds formulae for central difference

* Update author list of PartialDerivative.java

Added:
 - @tibuch
 - @maarzt

Author: tibuch

pom.xml

@@ -5,7 +5,7 @@
 	<parent>
 		<groupId>net.imglib2</groupId>
 		<artifactId>pom-imglib2</artifactId>
-		<version>10.0.1</version>
+		<version>10.0.2</version>
 		<relativePath />
 	</parent>
 

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

@@ -43,9 +43,9 @@
 
 import net.imglib2.Cursor;
 import net.imglib2.FinalInterval;
-import net.imglib2.RandomAccess;
 import net.imglib2.RandomAccessible;
 import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.loops.LoopBuilder;
 import net.imglib2.type.numeric.NumericType;
 import net.imglib2.util.Intervals;
 import net.imglib2.view.IntervalView;
@@ -54,13 +54,18 @@
 /**
  * @author Tobias Pietzsch
  * @author Philipp Hanslovsky
+ * @author Tim-Oliver Buchholz
+ * @author Matthias Arzt
  *
  */
 public class PartialDerivative
 {
 	// nice version...
 	/**
-	 * Compute the partial derivative of source in a particular dimension.
+	 * Compute the partial derivative (central difference approximation) of source
+	 * in a particular dimension:
+	 * {@code d_f( x ) = ( f( x + e ) - f( x - e ) ) / 2},
+	 * where {@code e} is the unit vector along that dimension.
 	 *
 	 * @param source
 	 *            source image, has to provide valid data in the interval of the
@@ -85,7 +90,10 @@ public static < T extends NumericType< T > > void gradientCentralDifference2( fi
 
 	// parallel version...
 	/**
-	 * Compute the partial derivative of source in a particular dimension.
+	 * Compute the partial derivative (central difference approximation) of source
+	 * in a particular dimension:
+	 * {@code d_f( x ) = ( f( x + e ) - f( x - e ) ) / 2},
+	 * where {@code e} is the unit vector along that dimension.
 	 *
 	 * @param source
 	 *            source image, has to provide valid data in the interval of the
@@ -157,81 +165,73 @@ public static < T extends NumericType< T > > void gradientCentralDifferenceParal
 
 	// fast version
 	/**
-	 * Compute the partial derivative of source in a particular dimension.
+	 * Compute the partial derivative (central difference approximation) of source
+	 * in a particular dimension:
+	 * {@code d_f( x ) = ( f( x + e ) - f( x - e ) ) / 2},
+	 * where {@code e} is the unit vector along that dimension.
 	 *
 	 * @param source
 	 *            source image, has to provide valid data in the interval of the
 	 *            gradient image plus a one pixel border in dimension.
-	 * @param gradient
+	 * @param result
 	 *            output image
 	 * @param dimension
 	 *            along which dimension the partial derivatives are computed
 	 */
-	public static < T extends NumericType< T > > void gradientCentralDifference( final RandomAccessible< T > source, final RandomAccessibleInterval< T > gradient, final int dimension )
+	public static < T extends NumericType< T > > void gradientCentralDifference( final RandomAccessible< T > source,
+			final RandomAccessibleInterval< T > result, final int dimension )
 	{
-		final int n = gradient.numDimensions();
-
-		final long[] min = new long[ n ];
-		gradient.min( min );
-		final long[] max = new long[ n ];
-		gradient.max( max );
-		final long[] shiftback = new long[ n ];
-		for ( int d = 0; d < n; ++d )
-			shiftback[ d ] = min[ d ] - max[ d ];
-
-		final RandomAccess< T > result = gradient.randomAccess();
-		final RandomAccess< T > back = source.randomAccess( Intervals.translate( gradient, -1, dimension ) );
-		final RandomAccess< T > front = source.randomAccess( Intervals.translate( gradient, 1, dimension ) );
-
-		result.setPosition( min );
-		back.setPosition( min );
-		back.bck( dimension );
-		front.setPosition( min );
-		front.fwd( dimension );
-
-		final long max0 = max[ 0 ];
-		while ( true )
-		{
-			// process pixel
-			final T t = result.get();
-			t.set( front.get() );
-			t.sub( back.get() );
-			t.mul( 0.5 );
+		final RandomAccessibleInterval< T > back = Views.interval( source, Intervals.translate( result, -1, dimension ) );
+		final RandomAccessibleInterval< T > front = Views.interval( source, Intervals.translate( result, 1, dimension ) );
+
+		LoopBuilder.setImages( result, back, front ).forEachPixel( ( r, b, f ) -> {
+			r.set( f );
+			r.sub( b );
+			r.mul( 0.5 );
+		} );
+	}
 
-			// move to next pixel
-			// check dimension 0 separately to avoid the loop over d in most
-			// iterations
-			if ( result.getLongPosition( 0 ) == max0 )
-			{
-				if ( n == 1 )
-					return;
-				result.move( shiftback[ 0 ], 0 );
-				back.move( shiftback[ 0 ], 0 );
-				front.move( shiftback[ 0 ], 0 );
-				// now check the remaining dimensions
-				for ( int d = 1; d < n; ++d )
-					if ( result.getLongPosition( d ) == max[ d ] )
-					{
-						result.move( shiftback[ d ], d );
-						back.move( shiftback[ d ], d );
-						front.move( shiftback[ d ], d );
-						if ( d == n - 1 )
-							return;
-					}
-					else
-					{
-						result.fwd( d );
-						back.fwd( d );
-						front.fwd( d );
-						break;
-					}
-			}
-			else
-			{
-				result.fwd( 0 );
-				back.fwd( 0 );
-				front.fwd( 0 );
-			}
-		}
+	/**
+	 * Compute the backward difference of source in a particular dimension:
+	 * {@code d_f( x ) = ( f( x ) - f( x - e ) )}
+	 * where {@code e} is the unit vector along that dimension
+	 *
+	 * @param source source image, has to provide valid data in the interval of
+	 *            the gradient image plus a one pixel border in dimension.
+	 * @param result output image
+	 * @param dimension along which dimension the partial derivatives are computed
+	 */
+	public static < T extends NumericType< T > > void gradientBackwardDifference( final RandomAccessible< T > source,
+			final RandomAccessibleInterval< T > result, final int dimension )
+	{
+		final RandomAccessibleInterval< T > back = Views.interval( source, Intervals.translate( result, -1, dimension ) );
+		final RandomAccessibleInterval< T > front = Views.interval( source, result );
+
+		LoopBuilder.setImages( result, back, front ).forEachPixel( ( r, b, f ) -> {
+			r.set( f );
+			r.sub( b );
+		} );
+	}
+
+	/**
+	 * Compute the forward difference of source in a particular dimension:
+	 * {@code d_f( x ) = ( f( x + e ) - f( x ) )}
+	 * where {@code e} is the unit vector along that dimension
+	 
+	 * @param source source image, has to provide valid data in the interval of
+	 *            the gradient image plus a one pixel border in dimension.
+	 * @param result output image
+	 * @param dimension along which dimension the partial derivatives are computed
+	 */
+	public static < T extends NumericType< T > > void gradientForwardDifference( final RandomAccessible< T > source,
+			final RandomAccessibleInterval< T > result, final int dimension )
+	{
+		final RandomAccessibleInterval< T > back = Views.interval( source, result );
+		final RandomAccessibleInterval< T > front = Views.interval( source, Intervals.translate( result, 1, dimension ) );
+
+		LoopBuilder.setImages( result, back, front ).forEachPixel( ( r, b, f ) -> {
+			r.set( f );
+			r.sub( b );
+		} );
 	}
 }

src/test/java/net/imglib2/algorithm/gradient/PartialDerivativeTest.java

@@ -0,0 +1,100 @@
+package net.imglib2.algorithm.gradient;
+
+import net.imglib2.Interval;
+import net.imglib2.RandomAccessibleInterval;
+import net.imglib2.img.array.ArrayImgs;
+import net.imglib2.type.numeric.real.DoubleType;
+import net.imglib2.util.Intervals;
+import net.imglib2.view.Views;
+import org.junit.Test;
+
+import java.util.ArrayList;
+import java.util.List;
+
+import static junit.framework.TestCase.assertTrue;
+import static org.junit.Assert.assertArrayEquals;
+
+/**
+ * Tests {@link PartialDerivative}.
+ *
+ * @author Matthias Arzt
+ */
+public class PartialDerivativeTest
+{
+
+	private final double DELTA = 0.0001;
+
+	private final RandomAccessibleInterval< DoubleType > testImage =
+			createImage2dMinSizeValues( /* min: */ 0, 0, /* size: */ 5, 3, /* values: */
+					0.0, 0.0, 0.0, 0.0, 0.0,
+					0.0, 2.0, 0.0, -1.0, 0.0,
+					0.0, 0.0, 5.0, 0.0, 0.0
+			);
+
+	private final RandomAccessibleInterval< DoubleType > centralDifferenceExpected =
+			createImage2dMinSizeValues( /* min: */ 1, 0, /* size: */ 3, 3, /* values: */
+					0.0, 0.0, 0.0,
+					0.0, -1.5, 0.0,
+					2.5, 0.0, -2.5
+			);
+
+	private final RandomAccessibleInterval< DoubleType > backwardDifferenceExpected =
+			createImage2dMinSizeValues( /* min: */ 1, 0, /* size: */ 4, 3, /* values: */
+					0.0, 0.0, 0.0, 0.0,
+					2.0, -2.0, -1.0, 1.0,
+					0.0, 5.0, -5.0, 0.0
+			);
+
+	public static RandomAccessibleInterval< DoubleType > createImage2dMinSizeValues( long minX, long minY, long sizeX, long sizeY, double... values )
+	{
+		Interval interval = Intervals.createMinSize( minX, minY, sizeX, sizeY );
+		return Views.translate( ArrayImgs.doubles( values, Intervals.dimensionsAsLongArray( interval ) ),
+				Intervals.minAsLongArray( interval ) );
+	}
+
+	@Test
+	public void testGradientCentralDifferenceX()
+	{
+		RandomAccessibleInterval< DoubleType > data = testImage;
+		RandomAccessibleInterval< DoubleType > result = emptyArrayImg( centralDifferenceExpected );
+		PartialDerivative.gradientCentralDifference( data, result, 0 );
+		assertImagesEqual( centralDifferenceExpected, result );
+	}
+
+	@Test
+	public void testGradientBackwardDifferenceX()
+	{
+		RandomAccessibleInterval< DoubleType > data = testImage;
+		RandomAccessibleInterval< DoubleType > result = emptyArrayImg( backwardDifferenceExpected );
+		PartialDerivative.gradientBackwardDifference( data, result, 0 );
+		assertImagesEqual( backwardDifferenceExpected, result );
+	}
+
+	@Test
+	public void testForwardDifferenceX()
+	{
+		RandomAccessibleInterval< DoubleType > data = testImage;
+		RandomAccessibleInterval< DoubleType > expected = Views.translate( backwardDifferenceExpected, -1, 0 );
+		RandomAccessibleInterval< DoubleType > result = emptyArrayImg( expected );
+		PartialDerivative.gradientForwardDifference( data, result, 0 );
+		assertImagesEqual( expected, result );
+	}
+
+	private void assertImagesEqual( RandomAccessibleInterval< DoubleType > expected, RandomAccessibleInterval< DoubleType > actual )
+	{
+		assertTrue( Intervals.equals( expected, actual ) );
+		assertArrayEquals( imageAsArray( expected ), imageAsArray( actual ), DELTA );
+	}
+
+	private double[] imageAsArray( RandomAccessibleInterval< DoubleType > image )
+	{
+		List< Double > values = new ArrayList<>();
+		Views.iterable( image ).forEach( x -> values.add( x.getRealDouble() ) );
+		return values.stream().mapToDouble( x -> x ).toArray();
+	}
+
+	private RandomAccessibleInterval< DoubleType > emptyArrayImg( Interval interval )
+	{
+		return Views.translate( ArrayImgs.doubles( Intervals.dimensionsAsLongArray( interval ) ), Intervals.minAsLongArray( interval ) );
+	}
+}