From bcf6b1a703a98a6310a927b5fafc791ea2c0d1e3 Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Wed, 30 Oct 2024 15:37:12 -0400 Subject: [PATCH] fix: jacobian determinant source for 2D transformations * minor cleanup and optimization --- .../JacobianDeterminantRandomAccess.java | 42 +++++++++++-------- .../source/JacobianDeterminantSource.java | 35 ++++++---------- 2 files changed, 36 insertions(+), 41 deletions(-) diff --git a/src/main/java/bigwarp/source/JacobianDeterminantRandomAccess.java b/src/main/java/bigwarp/source/JacobianDeterminantRandomAccess.java index bd11f416..43d7c38d 100644 --- a/src/main/java/bigwarp/source/JacobianDeterminantRandomAccess.java +++ b/src/main/java/bigwarp/source/JacobianDeterminantRandomAccess.java @@ -40,8 +40,9 @@ public class JacobianDeterminantRandomAccess< T extends RealType> extends Abs { protected DifferentiableRealTransform transform; - final private T value; - final double[] warpRes; + private final T value; + private final double[] x; + private final double[] x3; protected JacobianDeterminantRandomAccess( double[] dimensions ) { @@ -50,10 +51,11 @@ protected JacobianDeterminantRandomAccess( double[] dimensions ) protected JacobianDeterminantRandomAccess( final double[] dimensions, final T value, final DifferentiableRealTransform transform ) { - super( dimensions.length ); + super( 3 ); setTransform( transform ); - this.value = value; - warpRes = new double[ numDimensions() ]; + this.value = value.copy(); + x = new double[ dimensions.length ]; + x3 = new double[3]; } public void setTransform( final DifferentiableRealTransform transform ) @@ -67,29 +69,29 @@ public void setTransform( final DifferentiableRealTransform transform ) @Override public T get() { + // TODO check the below // copy value here so this is thread safe - T out = value.copy(); if( transform == null ) { - out.setZero(); - return out; + value.setOne(); + return value; } // compute the jacobian determinant at this point - double[] x = new double[ numDimensions() ]; - localize( x ); - AffineTransform jacobian = transform.jacobian( x ); - DMatrixRMaj jacMtx = new DMatrixRMaj(); + localize( x3 ); + System.arraycopy(x3, 0, x, 0, x.length); + + final AffineTransform jacobian = transform.jacobian( x ); + final DMatrixRMaj jacMtx = new DMatrixRMaj(); jacMtx.data = jacobian.getRowPackedCopy(); - out.setReal( CommonOps_DDRM.det( jacMtx ) ); + value.setReal( CommonOps_DDRM.det( jacMtx ) ); - return out; + return value; } public RealRandomAccess copy() { - return new JacobianDeterminantRandomAccess< T >( new double[ position.length ], value.copy(), - transform ); + return new JacobianDeterminantRandomAccess(new double[x.length], value.copy(), transform); } public RealRandomAccess copyRandomAccess() @@ -266,10 +268,14 @@ public static class JacobianDeterminantRandomAccessibleInterval ra; - public JacobianDeterminantRandomAccessibleInterval( Interval interval, T t, DifferentiableRealTransform warp ) + public JacobianDeterminantRandomAccessibleInterval( Interval interval, T t, DifferentiableRealTransform transform ) { super( interval ); - ra = new JacobianDeterminantRandomAccess< T >( new double[ interval.numDimensions() ], t, warp ); + int nd = interval.dimension(2) < 2 ? 2 : 3; + if( nd == 2) + ra = new JacobianDeterminantRandomAccess< T >( new double[ nd ], t, transform ); + else + ra = new JacobianDeterminantRandomAccess< T >( new double[ nd ], t, transform ); } @Override diff --git a/src/main/java/bigwarp/source/JacobianDeterminantSource.java b/src/main/java/bigwarp/source/JacobianDeterminantSource.java index e2b3fdaa..130f7646 100644 --- a/src/main/java/bigwarp/source/JacobianDeterminantSource.java +++ b/src/main/java/bigwarp/source/JacobianDeterminantSource.java @@ -21,8 +21,6 @@ */ package bigwarp.source; -import java.util.Arrays; - import bdv.viewer.Interpolation; import bdv.viewer.Source; import bigwarp.BigWarpData; @@ -32,7 +30,6 @@ import net.imglib2.Cursor; import net.imglib2.Interval; import net.imglib2.RandomAccessibleInterval; -import net.imglib2.RealRandomAccess; import net.imglib2.RealRandomAccessible; import net.imglib2.realtransform.AffineTransform3D; import net.imglib2.realtransform.BoundingBoxEstimation; @@ -58,20 +55,13 @@ public JacobianDeterminantSource( String name, BigWarpData data, T t ) { this.name = name; this.type = t; - sourceData = data; - //RandomAccessibleInterval fixedsrc = sourceData.sources.get( 1 ).getSpimSource().getSource( 0, 0 ); -// interval = sourceData.sources.get( sourceData.targetSourceIndices[ 0 ] ).getSpimSource().getSource( 0, 0 ); -// interval = sourceData.getTargetSource( 0 ).getSpimSource().getSource( 0, 0 ); - final BoundingBoxEstimation bbe = new BoundingBoxEstimation(); final AffineTransform3D affine = new AffineTransform3D(); data.getTargetSource( 0 ).getSpimSource().getSourceTransform( 0, 0, affine ); - interval = bbe.estimatePixelInterval( affine, data.getTargetSource( 0 ).getSpimSource().getSource( 0, 0 ) ); + interval = bbe.estimatePixelInterval( affine, data.getTargetSource( 0 ).getSpimSource().getSource( 0, 0 ) ); - -// VoxelDimensions srcVoxDims = sourceData.sources.get( sourceData.targetSourceIndices[ 0 ] ).getSpimSource().getVoxelDimensions(); VoxelDimensions srcVoxDims = sourceData.getTargetSource( 0 ).getSpimSource().getVoxelDimensions(); String unit = "pix"; if( srcVoxDims != null ) @@ -81,7 +71,7 @@ public JacobianDeterminantSource( String name, BigWarpData data, T t ) jacDetImg = new JacobianDeterminantRandomAccess.JacobianDeterminantRandomAccessibleInterval< T >( interval, t, null ); } - + public double getMax( LandmarkTableModel lm ) { double maxVal = 0.0; @@ -108,31 +98,30 @@ public void setTransform( final DifferentiableRealTransform transform ) jacDetImg.setTransform( transform ); } + @Deprecated public void debug( double[] pt ) { - RealRandomAccess rra = jacDetImg.realRandomAccess(); - - rra.setPosition( pt ); - - System.out.println("at : " + Arrays.toString( pt ) ); - System.out.println( "get val: " + rra.get()); - +// RealRandomAccess rra = jacDetImg.realRandomAccess(); +// +// rra.setPosition( pt ); +// +// System.out.println("at : " + Arrays.toString( pt ) ); +// System.out.println( "get val: " + rra.get()); +// // double[] warpRes = new double[ jacDetImg.ra.warp.numTargetDimensions() ]; // jacDetImg.ra.warp.apply( pt, warpRes ); // // System.out.println( "base res: " + baseRes[0] + " " + baseRes[1]); // System.out.println( "warp res: " + warpRes[0] + " " + warpRes[1]); - } public double[] minMax() { - double[] minmax = new double[ 2 ]; + final double[] minmax = new double[ 2 ]; minmax[ 0 ] = Double.MAX_VALUE; minmax[ 1 ] = Double.MIN_VALUE; - Cursor curs = Views.iterable( this.getSource( 0,0 ) ).cursor(); - + final Cursor curs = this.getSource( 0,0 ).cursor(); while( curs.hasNext() ) { double val = curs.next().getRealDouble();