+ * The axial offset follows the formula: z' = z + (R - √(R² - r²)), where R is the radius of curvature + * and r is the distance from the optical axis. + *
+ * Note: This transform is only invertible for points within radius R from the optical axis. + * Points at radius r greater than R will produce NaN values. + */ +public class SphericalCurvatureDistortionTransform implements InvertibleRealTransform +{ + private final int numDimensions; + + private final int axialDimension; + + private final double R; + + private final double sign; + + /** + * Creates a spherical curvature distortion transform. + * + * @param numDimensions number of dimensions + * @param axialDimension the dimension index of the optical axis where curvature is applied + * @param R radius of curvature + */ + public SphericalCurvatureDistortionTransform( final int numDimensions, final int axialDimension, double R ) + { + this(numDimensions, axialDimension, R, 1.0); + } + + private SphericalCurvatureDistortionTransform( final int numDimensions, final int axialDimension, double R, double sign) + { + this.numDimensions = numDimensions; + this.axialDimension = axialDimension; + this.R = R; + this.sign = sign; + } + + @Override + public int numSourceDimensions() + { + return numDimensions; + } + + @Override + public int numTargetDimensions() + { + return numDimensions; + } + + @Override + public void apply( double[] source, double[] target ) + { + for ( int i = 0; i < numDimensions; i++ ) + { + if ( i == axialDimension ) + target[ i ] = source[ i ] + sign * axialOffset( squaredRadius( source ) ); + else + target[ i ] = source[ i ]; + } + } + + private void apply( double s, double[] source, double[] target ) + { + for ( int i = 0; i < numDimensions; i++ ) + { + if ( i == axialDimension ) + target[ i ] = source[ i ] + s * sign * axialOffset( squaredRadius( source ) ); + else + target[ i ] = source[ i ]; + } + } + + @Override + public void apply( RealLocalizable source, RealPositionable target ) + { + + for ( int i = 0; i < numDimensions; i++ ) + { + if ( i == axialDimension ) + target.setPosition( source.getDoublePosition( i ) + sign * axialOffset( squaredRadius( source ) ), i ); + else + target.setPosition( source.getDoublePosition( i ), i ); + } + } + + private void apply( double s, RealLocalizable source, RealPositionable target ) + { + + for ( int i = 0; i < numDimensions; i++ ) + { + if ( i == axialDimension ) + target.setPosition( source.getDoublePosition( i ) + s * sign * axialOffset( squaredRadius( source ) ), i ); + else + target.setPosition( source.getDoublePosition( i ), i ); + } + } + + @Override + public void applyInverse( double[] source, double[] target ) + { + apply( -1, target, source ); + } + + @Override + public void applyInverse( RealPositionable source, RealLocalizable target ) + { + apply( -1, target, source ); + } + + @Override + public SphericalCurvatureDistortionTransform copy() + { + return new SphericalCurvatureDistortionTransform( numDimensions, axialDimension, R, sign ); + } + + @Override + public InvertibleRealTransform inverse() + { + return new SphericalCurvatureDistortionTransform( numDimensions, axialDimension, R, -1 * sign ); + } + + private double squaredRadius( double[] position ) + { + double r = 0; + for ( int i = 0; i < position.length; i++ ) + if ( i != axialDimension ) + r += position[ i ] * position[ i ]; + + return r; + } + + private double squaredRadius( RealLocalizable position ) + { + double r = 0; + for ( int i = 0; i < position.numDimensions(); i++ ) + if ( i != axialDimension ) + r += position.getDoublePosition( i ) * position.getDoublePosition( i ); + + return r; + } + + private double axialOffset( double r2 ) + { + return R - Math.sqrt( ( R * R ) - r2 ); + } + +} diff --git a/src/test/java/net/imglib2/realtransform/IterableInverseTests.java b/src/test/java/net/imglib2/realtransform/IterableInverseTests.java index 731fc77..be8cf7b 100644 --- a/src/test/java/net/imglib2/realtransform/IterableInverseTests.java +++ b/src/test/java/net/imglib2/realtransform/IterableInverseTests.java @@ -33,10 +33,16 @@ */ package net.imglib2.realtransform; +import static org.junit.Assert.assertArrayEquals; + +import java.util.Arrays; + import org.junit.Assert; import org.junit.Test; +import net.imglib2.RealInterval; import net.imglib2.RealPoint; +import net.imglib2.iterator.LocalizingRealIntervalIterator; import net.imglib2.position.FunctionRealRandomAccessible; import net.imglib2.realtransform.inverse.DifferentiableRealTransform; import net.imglib2.realtransform.inverse.InverseRealTransformGradientDescent; @@ -294,6 +300,39 @@ public void testAffineInverse() Assert.assertArrayEquals( "rotation matrix inverse", p, q, EPS ); } + /** + * Tests that the provided {@link InvertibleRealTransform} is within epsilon + * of being invertible over the given {@link RealInterval}. Points to test + * are sampled at a given spacing. + * + * @param tform + * the transformation to tests + * @param interval + * over which to test + * @param sampleSpacing + * the spacing between samples + * @param epsilon + * the allowable error + */ + public static void inverseTransformTestHelper( InvertibleRealTransform tform, final RealInterval interval, double[] sampleSpacing, final double epsilon ) + { + int nd = interval.numDimensions(); + final LocalizingRealIntervalIterator it = new LocalizingRealIntervalIterator( interval, sampleSpacing ); + final RealPoint p = new RealPoint( nd ); + final RealPoint q = new RealPoint( nd ); + while ( it.hasNext() ) + { + it.fwd(); + tform.apply( it, p ); + tform.applyInverse( q, p ); + assertArrayEquals( + String.format( "not invertible at point %s", Arrays.toString( it.positionAsDoubleArray() ) ), + it.positionAsDoubleArray(), + q.positionAsDoubleArray(), + epsilon ); + } + } + private class IterativeAffineInverse extends AffineTransform implements DifferentiableRealTransform { final AffineTransform jacobian; diff --git a/src/test/java/net/imglib2/realtransform/distortion/SphericalCurvatureDistortionTest.java b/src/test/java/net/imglib2/realtransform/distortion/SphericalCurvatureDistortionTest.java new file mode 100644 index 0000000..289bb0e --- /dev/null +++ b/src/test/java/net/imglib2/realtransform/distortion/SphericalCurvatureDistortionTest.java @@ -0,0 +1,23 @@ +package net.imglib2.realtransform.distortion; + +import org.junit.Test; + +import net.imglib2.FinalRealInterval; +import net.imglib2.realtransform.IterableInverseTests; +import net.imglib2.util.Intervals; + +public class SphericalCurvatureDistortionTest +{ + private static final double EPS = 1E-6; + + @Test + public void inverseTest() throws Exception + { + final SphericalCurvatureDistortionTransform tform = new SphericalCurvatureDistortionTransform( 3, 2, 1000.0 ); + final FinalRealInterval interval = Intervals.createMinMaxReal( -1.0, -1.0, -1.0, 2.01, 2.01, 2.01 ); + final double[] step = new double[] { 0.1, 0.1, 0.1 }; + + IterableInverseTests.inverseTransformTestHelper( tform, interval, step, EPS ); + } + +} \ No newline at end of file