Skip to content

Commit

Permalink
Differences (#44)
Browse files Browse the repository at this point in the history
* 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
  • Loading branch information
tibuch authored and hanslovsky committed Mar 27, 2018
1 parent 2b71c6d commit 0454000
Show file tree
Hide file tree
Showing 3 changed files with 170 additions and 70 deletions.
2 changes: 1 addition & 1 deletion pom.xml
Original file line number Diff line number Diff line change
Expand Up @@ -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>

Expand Down
138 changes: 69 additions & 69 deletions src/main/java/net/imglib2/algorithm/gradient/PartialDerivative.java
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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 );
} );
}
}
Original file line number Diff line number Diff line change
@@ -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 ) );
}
}

0 comments on commit 0454000

Please sign in to comment.