Add Spherical curvature distortion transformation

src/main/java/net/imglib2/realtransform/distortion/SphericalCurvatureDistortionTransform.java

@@ -0,0 +1,155 @@
+package net.imglib2.realtransform.distortion;
+
+import net.imglib2.RealLocalizable;
+import net.imglib2.RealPositionable;
+import net.imglib2.realtransform.InvertibleRealTransform;
+
+/**
+ * Applies spherical field curvature distortion along the axial dimension.
+ * <p>
+ * 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.
+ * <p>
+ * 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 );
+	}
+
+}

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;

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 );
+	}
+
+}