From b18786510ce8ac2252491b5290255691fd997975 Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Wed, 26 Feb 2025 15:00:54 -0500 Subject: [PATCH 1/7] chore: DistanceTransform use new methods * getType * addDimension with bounds --- .../morphology/distance/DistanceTransform.java | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index 24f4e2482..02e4c35aa 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -528,7 +528,7 @@ public static < T extends RealType< T >, U extends RealType< U >, V extends Real { transformAlongDimension( ( RandomAccessible< T > ) Views.addDimension( source ), - Views.interval( Views.addDimension( target ), new FinalInterval( target.dimension( 0 ), 1 ) ), + Views.addDimension( target, 0, 0 ), d, 0 ); } @@ -601,7 +601,7 @@ public static < T extends RealType< T >, U extends RealType< U >, V extends Real { transformAlongDimensionParallel( ( RandomAccessible< T > ) Views.addDimension( source ), - Views.interval( Views.addDimension( target ), new FinalInterval( target.dimension( 0 ), 1 ) ), + Views.addDimension( target, 0, 0 ), d, 0, es, @@ -734,8 +734,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final DISTANCE_TYPE distanceType, final double... weights ) { - - final U maxVal = Util.getTypeFromInterval( tmp ).createVariable(); + final U maxVal = tmp.getType().copy(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -789,7 +788,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final int nTasks, final double... weights ) throws InterruptedException, ExecutionException { - final U maxVal = Util.getTypeFromInterval( tmp ).createVariable(); + final U maxVal = tmp.getType().copy(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -916,7 +915,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final RandomAccessibleInterval< V > target, final Distance d ) { - final U maxVal = Util.getTypeFromInterval( tmp ).createVariable(); + final U maxVal = tmp.getType().copy(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -963,7 +962,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final ExecutorService es, final int nTasks ) throws InterruptedException, ExecutionException { - final U maxVal = Util.getTypeFromInterval( tmp ).createVariable(); + final U maxVal = tmp.getType().copy(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -1005,7 +1004,7 @@ private static < T extends RealType< T >, U extends RealType< U >, V extends Rea { transformL1AlongDimension( ( RandomAccessible< T > ) Views.addDimension( source ), - Views.interval( Views.addDimension( target ), new FinalInterval( target.dimension( 0 ), 1 ) ), + Views.addDimension( target, 0, 0 ), 0, weights[ 0 ] ); } @@ -1073,7 +1072,7 @@ private static < T extends RealType< T >, U extends RealType< U >, V extends Rea { transformL1AlongDimensionParallel( ( RandomAccessible< T > ) Views.addDimension( source ), - Views.interval( Views.addDimension( target ), new FinalInterval( target.dimension( 0 ), 1 ) ), + Views.addDimension( target, 0, 0 ), 0, weights[ 0 ], es, From 6452bcc10f8607f6e2da56bb7a042e219d080b31 Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Wed, 26 Feb 2025 16:17:37 -0500 Subject: [PATCH 2/7] wip/feat: toward propagating labels in DistanceTransform --- .../distance/DistanceTransform.java | 355 ++++++++++++++++++ 1 file changed, 355 insertions(+) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index 02e4c35aa..7f5756dcf 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -56,8 +56,10 @@ import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.img.cell.CellImg; import net.imglib2.img.cell.CellImgFactory; +import net.imglib2.loops.LoopBuilder; import net.imglib2.type.BooleanType; import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; import net.imglib2.type.numeric.integer.LongType; import net.imglib2.type.numeric.real.DoubleType; @@ -391,6 +393,130 @@ public static < T extends RealType< T > > void transform( transform( source, source, d ); } + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate and final results will be stored in + * {@code source} ({@link DoubleType} recommended). + *

+ * Also propagates a set of input labels so that it contains the nearest label + * everywhere. + * + * @param source + * Input function on which distance transform should be computed. + * @param labels + * Labels to be propagated. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input + * @param + * {@link IntegerType} the label type + */ + public static < T extends RealType< T >, L extends IntegerType< L > > void transformPropagateLabels( + final RandomAccessibleInterval< T > source, + final RandomAccessibleInterval< L > labels, + final DISTANCE_TYPE distanceType ) + { + // TODO generalize distance + transformPropagateLabels( source, source, labels, labels, new EuclidianDistanceIsotropic( 1 ) ); + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate and final results will be stored in + * {@code target} ({@link DoubleType} recommended). + * + * @param source + * Input function on which distance transform should be computed. + * @param target + * Final result of distance transform. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input distance type + * @param + * {@link RealType} result type + * @param + * {@link IntegerType} label type + * @param + * {@link IntegerType} label result type + */ + public static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > target, + final RandomAccessibleInterval< L > labels, + final RandomAccessibleInterval< M > labelsResult, + final Distance d ) + { + transformPropagateLabels( source, target, target, labels, labelsResult, d ); + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate results will be stored in {@code tmp} + * ({@link DoubleType} recommended). The output will be written into + * {@code target}. + * + * @param source + * Input function on which distance transform should be computed. + * @param tmp + * Storage for intermediate results. + * @param target + * Final result of distance transform. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input + * @param + * {@link RealType} intermediate results + * @param + * {@link RealType} output + */ + public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > tmp, + final RandomAccessibleInterval< V > target, + final RandomAccessible< L > labels, + final RandomAccessible< M > labelsResult, + final Distance d ) + { + assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; + final int nDim = source.numDimensions(); + final int lastDim = nDim - 1; + + if ( nDim == 1 ) + { + transformAlongDimensionPropagateLabels( + ( RandomAccessible< T > ) Views.addDimension( source ), + Views.addDimension( target, 0, 0 ), + Views.addDimension( labels ), + Views.addDimension( labelsResult ), + d, 0 ); + } + else + { + transformAlongDimensionPropagateLabels( source, tmp, labels, labelsResult, d, 0 ); + } + + for ( int dim = 1; dim < nDim; ++dim ) + { + if ( dim == lastDim ) + { + transformAlongDimensionPropagateLabels( tmp, target, labels, labelsResult, d, dim ); + } + else + { + transformAlongDimensionPropagateLabels( tmp, tmp, labels, labelsResult, d, dim ); + } + } + } + /** * Create * distance @@ -969,6 +1095,42 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R transform( converted, tmp, target, d, es, nTasks ); } + @SuppressWarnings( "unchecked" ) + public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( + final RandomAccessibleInterval< L > labels, + final DISTANCE_TYPE distanceType, + final double... weights ) + { + final RandomAccessibleInterval< DoubleType > distance = Util.getSuitableImgFactory( labels, new DoubleType()).create( labels ); + distance.getAt( 0 ).setReal( 0 ); + + labelTransform( labels, distance, distance, distanceType, weights ); + return ( RandomAccessibleInterval< T > ) distance; + } + + public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( + final RandomAccessible< L > source, + final RandomAccessibleInterval< T > distance, + final DISTANCE_TYPE distanceType, + final double... weights ) + { + labelTransform( source, distance, distance, distanceType, weights ); + } + + public static < L extends IntegerType< L >, T extends RealType< T >, U extends RealType< U > > void labelTransform( + final RandomAccessible< L > source, + final RandomAccessibleInterval< T > tmp, + final RandomAccessibleInterval< U > distance, + final DISTANCE_TYPE distanceType, + final double... weights ) + { + final U maxVal = distance.getType().copy(); + maxVal.setReal( maxVal.getMaxValue() ); + final Converter< L, U > converter = new LabelsToCost<>( 0, maxVal ); + final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); + transform( converted, tmp, distance, distanceType, weights ); + } + /** * Create binary distance transform on {@code source} using L1 distance. * Intermediate results will be stored in {@code tmp} ({@link DoubleType} @@ -1105,9 +1267,11 @@ private static < T extends RealType< T >, U extends RealType< U > > void transfo final int lastDim = target.numDimensions() - 1; final long size = target.dimension( dim ); final RealComposite< DoubleType > tmp = Views.collapseReal( createAppropriateOneDimensionalImage( size, new DoubleType() ) ).randomAccess().get(); + // do not permute if we already work on last dimension final Cursor< RealComposite< T > > s = Views.flatIterable( Views.collapseReal( dim == lastDim ? Views.interval( source, target ) : Views.permute( Views.interval( source, target ), dim, lastDim ) ) ).cursor(); final Cursor< RealComposite< U > > t = Views.flatIterable( Views.collapseReal( dim == lastDim ? target : Views.permute( target, dim, lastDim ) ) ).cursor(); + final RealComposite< LongType > lowerBoundDistanceIndex = Views.collapseReal( createAppropriateOneDimensionalImage( size, new LongType() ) ).randomAccess().get(); final RealComposite< DoubleType > envelopeIntersectLocation = Views.collapseReal( createAppropriateOneDimensionalImage( size + 1, new DoubleType() ) ).randomAccess().get(); @@ -1161,6 +1325,58 @@ private static < T extends RealType< T >, U extends RealType< U > > void transfo invokeAllAndWait( es, tasks ); } + + private static < T extends RealType< T >, U extends RealType< U > > void transformSingleColumn( + final RealComposite< T > source, + final RealComposite< U > target, + final RealComposite< LongType > lowerBoundDistanceIndex, + final RealComposite< DoubleType > envelopeIntersectLocation, + final RealComposite< LongType > nearestLabel, + final Distance d, + final int dim, + final long size ) + { + long k = 0; + + lowerBoundDistanceIndex.get( 0 ).set( 0 ); + envelopeIntersectLocation.get( 0 ).set( Double.NEGATIVE_INFINITY ); + envelopeIntersectLocation.get( 1 ).set( Double.POSITIVE_INFINITY ); + for ( long position = 1; position < size; ++position ) + { + long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + final double sourceAtPosition = source.get( position ).getRealDouble(); + + // the point at which these parabolas intersect + double s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + + for ( double envelopeValueAtK = envelopeIntersectLocation.get( k ).get(); s <= envelopeValueAtK; envelopeValueAtK = envelopeIntersectLocation.get( k ).get() ) + { + --k; + envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + } + ++k; + lowerBoundDistanceIndex.get( k ).set( position ); + envelopeIntersectLocation.get( k ).set( s ); + envelopeIntersectLocation.get( k + 1 ).set( Double.POSITIVE_INFINITY ); + } + + k = 0; + + for ( long position = 0; position < size; ++position ) + { + while ( envelopeIntersectLocation.get( k + 1 ).get() < position ) + { + ++k; + } + final long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + // copy necessary because of the following line, access to source + // after write to source -> source and target cannot be the same + target.get( position ).setReal( d.evaluate( position, envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), dim ) ); + } + + } + private static < T extends RealType< T >, U extends RealType< U > > void transformSingleColumn( final RealComposite< T > source, final RealComposite< U > target, @@ -1209,6 +1425,95 @@ private static < T extends RealType< T >, U extends RealType< U > > void transfo } + private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType< L >, M extends IntegerType< M > > void transformAlongDimensionPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > target, + final RandomAccessible< L > labelSource, + final RandomAccessible< M > labelTarget, + final Distance d, + final int dim ) + { + final int lastDim = target.numDimensions() - 1; + final long size = target.dimension( dim ); + final RealComposite< DoubleType > tmp = Views.collapseReal( createAppropriateOneDimensionalImage( size, new DoubleType() ) ).randomAccess().get(); + + // do not permute if we already work on last dimension + final Cursor< RealComposite< T > > s = Views.flatIterable( Views.collapseReal( dim == lastDim ? Views.interval( source, target ) : Views.permute( Views.interval( source, target ), dim, lastDim ) ) ).cursor(); + final Cursor< RealComposite< U > > t = Views.flatIterable( Views.collapseReal( dim == lastDim ? target : Views.permute( target, dim, lastDim ) ) ).cursor(); + + final Cursor< RealComposite< L > > ls = Views.flatIterable( + Views.collapseReal( dim == lastDim ? Views.interval( labelSource, target ) : Views.permute( Views.interval( labelSource, target ), dim, lastDim ) ) ).cursor(); + + final Cursor< RealComposite< M > > lt = Views.flatIterable( + Views.collapseReal( dim == lastDim ? Views.interval( labelTarget, target ) : Views.permute( Views.interval( labelTarget, target ), dim, lastDim ) ) ).cursor(); + + final RealComposite< LongType > lowerBoundDistanceIndex = Views.collapseReal( createAppropriateOneDimensionalImage( size, new LongType() ) ).randomAccess().get(); + final RealComposite< DoubleType > envelopeIntersectLocation = Views.collapseReal( createAppropriateOneDimensionalImage( size + 1, new DoubleType() ) ).randomAccess().get(); + + while ( s.hasNext() ) + { + final RealComposite< T > sourceComp = s.next(); + final RealComposite< U > targetComp = t.next(); + final RealComposite< L > labelComp = ls.next(); + final RealComposite< M > labelTargetComp = lt.next(); + for ( long i = 0; i < size; ++i ) + { + tmp.get( i ).set( sourceComp.get( i ).getRealDouble() ); + } + transformSingleColumnPropagateLabels( tmp, targetComp, labelComp, labelTargetComp, lowerBoundDistanceIndex, envelopeIntersectLocation, d, dim, size ); + } + } + + private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformSingleColumnPropagateLabels( + final RealComposite< T > source, + final RealComposite< U > target, + final RealComposite< L > labelsSource, + final RealComposite< M > labelsResult, + final RealComposite< LongType > lowerBoundDistanceIndex, + final RealComposite< DoubleType > envelopeIntersectLocation, + final Distance d, + final int dim, + final long size ) + { + long k = 0; + + lowerBoundDistanceIndex.get( 0 ).set( 0 ); + envelopeIntersectLocation.get( 0 ).set( Double.NEGATIVE_INFINITY ); + envelopeIntersectLocation.get( 1 ).set( Double.POSITIVE_INFINITY ); + for ( long position = 1; position < size; ++position ) + { + long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + final double sourceAtPosition = source.get( position ).getRealDouble(); + double s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + + for ( double envelopeValueAtK = envelopeIntersectLocation.get( k ).get(); s <= envelopeValueAtK; envelopeValueAtK = envelopeIntersectLocation.get( k ).get() ) + { + --k; + envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + } + ++k; + lowerBoundDistanceIndex.get( k ).set( position ); + envelopeIntersectLocation.get( k ).set( s ); + envelopeIntersectLocation.get( k + 1 ).set( Double.POSITIVE_INFINITY ); + } + + k = 0; + for ( long position = 0; position < size; ++position ) + { + while ( envelopeIntersectLocation.get( k + 1 ).get() < position ) + { + ++k; + } + final long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + // copy necessary because of the following line, access to source + // after write to source -> source and target cannot be the same + target.get( position ).setReal( d.evaluate( position, envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), dim ) ); + labelsResult.get( position ).setInteger( labelsSource.get( envelopeIndexAtK ).getIntegerLong() ); + } + + } + private static < T extends RealType< T >, U extends RealType< U > > void transformL1AlongDimension( final RandomAccessible< T > source, final RandomAccessibleInterval< U > target, @@ -1335,6 +1640,32 @@ public static int getLargestDimension( final Interval interval ) return IntStream.range( 0, interval.numDimensions() ).mapToObj( i -> new ValuePair<>( i, interval.dimension( i ) ) ).max( ( p1, p2 ) -> Long.compare( p1.getB(), p2.getB() ) ).get().getA(); } + private static , M extends IntegerType< M > > RandomAccessibleInterval< M > makeLabels( + final RandomAccessibleInterval< L > labels, + final M type) { + + final RandomAccessibleInterval< M > propagedLabels = Util.getSuitableImgFactory( labels, type ).create( labels ); + LoopBuilder.setImages( labels, propagedLabels ).forEachPixel( (s,d) -> { + d.setInteger( s.getIntegerLong() ); + }); + + return ( RandomAccessibleInterval< M > ) propagedLabels; + } + + private static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > makeDistances( + final long label, + final RandomAccessibleInterval< L > labels, + final T type) { + + final RandomAccessibleInterval< T > distances = Util.getSuitableImgFactory( labels, type ).create( labels ); + Views.pair( labels, distances ).view().interval( labels ).forEach( pair -> { + if( pair.getA().getIntegerLong() == label ) + pair.getB().setReal( Double.MAX_VALUE ); + }); + + return distances; + } + private static class BinaryMaskToCost< B extends BooleanType< B >, R extends RealType< R > > implements Converter< B, R > { @@ -1357,4 +1688,28 @@ public void convert( final B input, final R output ) } + private static class LabelsToCost< L extends IntegerType< L >, R extends RealType< R > > implements Converter< L, R > + { + private final R maxValForR; + + private final R zero; + + private final long label; + + public LabelsToCost( final long label, final R maxValForR ) + { + this.maxValForR = maxValForR; + this.label = label; + this.zero = maxValForR.createVariable(); + zero.setZero(); + } + + @Override + public void convert( final L input, final R output ) + { + output.set( input.getIntegerLong() == label ? zero : maxValForR ); + } + + } + } From 19cb4ec3d217af39308b5a4bd9c904f4d7439cef Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Thu, 27 Feb 2025 15:42:58 -0500 Subject: [PATCH 3/7] feat: add parallel methods for propagating labels * refactoring --- .../distance/DistanceTransform.java | 758 ++++++++++-------- 1 file changed, 406 insertions(+), 352 deletions(-) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index 7f5756dcf..1546afc08 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -56,7 +56,6 @@ import net.imglib2.img.array.ArrayImgFactory; import net.imglib2.img.cell.CellImg; import net.imglib2.img.cell.CellImgFactory; -import net.imglib2.loops.LoopBuilder; import net.imglib2.type.BooleanType; import net.imglib2.type.NativeType; import net.imglib2.type.numeric.IntegerType; @@ -83,6 +82,7 @@ *

* * @author Philipp Hanslovsky + * @author John Bogovic */ public class DistanceTransform { @@ -393,129 +393,7 @@ public static < T extends RealType< T > > void transform( transform( source, source, d ); } - /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate and final results will be stored in - * {@code source} ({@link DoubleType} recommended). - *

- * Also propagates a set of input labels so that it contains the nearest label - * everywhere. - * - * @param source - * Input function on which distance transform should be computed. - * @param labels - * Labels to be propagated. - * @param d - * {@link Distance} between two points. - * @param - * {@link RealType} input - * @param - * {@link IntegerType} the label type - */ - public static < T extends RealType< T >, L extends IntegerType< L > > void transformPropagateLabels( - final RandomAccessibleInterval< T > source, - final RandomAccessibleInterval< L > labels, - final DISTANCE_TYPE distanceType ) - { - // TODO generalize distance - transformPropagateLabels( source, source, labels, labels, new EuclidianDistanceIsotropic( 1 ) ); - } - - /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate and final results will be stored in - * {@code target} ({@link DoubleType} recommended). - * - * @param source - * Input function on which distance transform should be computed. - * @param target - * Final result of distance transform. - * @param d - * {@link Distance} between two points. - * @param - * {@link RealType} input distance type - * @param - * {@link RealType} result type - * @param - * {@link IntegerType} label type - * @param - * {@link IntegerType} label result type - */ - public static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > target, - final RandomAccessibleInterval< L > labels, - final RandomAccessibleInterval< M > labelsResult, - final Distance d ) - { - transformPropagateLabels( source, target, target, labels, labelsResult, d ); - } - - /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate results will be stored in {@code tmp} - * ({@link DoubleType} recommended). The output will be written into - * {@code target}. - * - * @param source - * Input function on which distance transform should be computed. - * @param tmp - * Storage for intermediate results. - * @param target - * Final result of distance transform. - * @param d - * {@link Distance} between two points. - * @param - * {@link RealType} input - * @param - * {@link RealType} intermediate results - * @param - * {@link RealType} output - */ - public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > tmp, - final RandomAccessibleInterval< V > target, - final RandomAccessible< L > labels, - final RandomAccessible< M > labelsResult, - final Distance d ) - { - assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; - final int nDim = source.numDimensions(); - final int lastDim = nDim - 1; - - if ( nDim == 1 ) - { - transformAlongDimensionPropagateLabels( - ( RandomAccessible< T > ) Views.addDimension( source ), - Views.addDimension( target, 0, 0 ), - Views.addDimension( labels ), - Views.addDimension( labelsResult ), - d, 0 ); - } - else - { - transformAlongDimensionPropagateLabels( source, tmp, labels, labelsResult, d, 0 ); - } - for ( int dim = 1; dim < nDim; ++dim ) - { - if ( dim == lastDim ) - { - transformAlongDimensionPropagateLabels( tmp, target, labels, labelsResult, d, dim ); - } - else - { - transformAlongDimensionPropagateLabels( tmp, tmp, labels, labelsResult, d, dim ); - } - } - } /** * Create @@ -1095,42 +973,6 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R transform( converted, tmp, target, d, es, nTasks ); } - @SuppressWarnings( "unchecked" ) - public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( - final RandomAccessibleInterval< L > labels, - final DISTANCE_TYPE distanceType, - final double... weights ) - { - final RandomAccessibleInterval< DoubleType > distance = Util.getSuitableImgFactory( labels, new DoubleType()).create( labels ); - distance.getAt( 0 ).setReal( 0 ); - - labelTransform( labels, distance, distance, distanceType, weights ); - return ( RandomAccessibleInterval< T > ) distance; - } - - public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( - final RandomAccessible< L > source, - final RandomAccessibleInterval< T > distance, - final DISTANCE_TYPE distanceType, - final double... weights ) - { - labelTransform( source, distance, distance, distanceType, weights ); - } - - public static < L extends IntegerType< L >, T extends RealType< T >, U extends RealType< U > > void labelTransform( - final RandomAccessible< L > source, - final RandomAccessibleInterval< T > tmp, - final RandomAccessibleInterval< U > distance, - final DISTANCE_TYPE distanceType, - final double... weights ) - { - final U maxVal = distance.getType().copy(); - maxVal.setReal( maxVal.getMaxValue() ); - final Converter< L, U > converter = new LabelsToCost<>( 0, maxVal ); - final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); - transform( converted, tmp, distance, distanceType, weights ); - } - /** * Create binary distance transform on {@code source} using L1 distance. * Intermediate results will be stored in {@code tmp} ({@link DoubleType} @@ -1325,58 +1167,6 @@ private static < T extends RealType< T >, U extends RealType< U > > void transfo invokeAllAndWait( es, tasks ); } - - private static < T extends RealType< T >, U extends RealType< U > > void transformSingleColumn( - final RealComposite< T > source, - final RealComposite< U > target, - final RealComposite< LongType > lowerBoundDistanceIndex, - final RealComposite< DoubleType > envelopeIntersectLocation, - final RealComposite< LongType > nearestLabel, - final Distance d, - final int dim, - final long size ) - { - long k = 0; - - lowerBoundDistanceIndex.get( 0 ).set( 0 ); - envelopeIntersectLocation.get( 0 ).set( Double.NEGATIVE_INFINITY ); - envelopeIntersectLocation.get( 1 ).set( Double.POSITIVE_INFINITY ); - for ( long position = 1; position < size; ++position ) - { - long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - final double sourceAtPosition = source.get( position ).getRealDouble(); - - // the point at which these parabolas intersect - double s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); - - for ( double envelopeValueAtK = envelopeIntersectLocation.get( k ).get(); s <= envelopeValueAtK; envelopeValueAtK = envelopeIntersectLocation.get( k ).get() ) - { - --k; - envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); - } - ++k; - lowerBoundDistanceIndex.get( k ).set( position ); - envelopeIntersectLocation.get( k ).set( s ); - envelopeIntersectLocation.get( k + 1 ).set( Double.POSITIVE_INFINITY ); - } - - k = 0; - - for ( long position = 0; position < size; ++position ) - { - while ( envelopeIntersectLocation.get( k + 1 ).get() < position ) - { - ++k; - } - final long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - // copy necessary because of the following line, access to source - // after write to source -> source and target cannot be the same - target.get( position ).setReal( d.evaluate( position, envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), dim ) ); - } - - } - private static < T extends RealType< T >, U extends RealType< U > > void transformSingleColumn( final RealComposite< T > source, final RealComposite< U > target, @@ -1425,136 +1215,47 @@ private static < T extends RealType< T >, U extends RealType< U > > void transfo } - private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType< L >, M extends IntegerType< M > > void transformAlongDimensionPropagateLabels( + private static < T extends RealType< T >, U extends RealType< U > > void transformL1AlongDimension( final RandomAccessible< T > source, final RandomAccessibleInterval< U > target, - final RandomAccessible< L > labelSource, - final RandomAccessible< M > labelTarget, - final Distance d, - final int dim ) + final int dim, + final double weight ) { final int lastDim = target.numDimensions() - 1; final long size = target.dimension( dim ); final RealComposite< DoubleType > tmp = Views.collapseReal( createAppropriateOneDimensionalImage( size, new DoubleType() ) ).randomAccess().get(); - // do not permute if we already work on last dimension final Cursor< RealComposite< T > > s = Views.flatIterable( Views.collapseReal( dim == lastDim ? Views.interval( source, target ) : Views.permute( Views.interval( source, target ), dim, lastDim ) ) ).cursor(); final Cursor< RealComposite< U > > t = Views.flatIterable( Views.collapseReal( dim == lastDim ? target : Views.permute( target, dim, lastDim ) ) ).cursor(); - final Cursor< RealComposite< L > > ls = Views.flatIterable( - Views.collapseReal( dim == lastDim ? Views.interval( labelSource, target ) : Views.permute( Views.interval( labelSource, target ), dim, lastDim ) ) ).cursor(); - - final Cursor< RealComposite< M > > lt = Views.flatIterable( - Views.collapseReal( dim == lastDim ? Views.interval( labelTarget, target ) : Views.permute( Views.interval( labelTarget, target ), dim, lastDim ) ) ).cursor(); - - final RealComposite< LongType > lowerBoundDistanceIndex = Views.collapseReal( createAppropriateOneDimensionalImage( size, new LongType() ) ).randomAccess().get(); - final RealComposite< DoubleType > envelopeIntersectLocation = Views.collapseReal( createAppropriateOneDimensionalImage( size + 1, new DoubleType() ) ).randomAccess().get(); - while ( s.hasNext() ) { final RealComposite< T > sourceComp = s.next(); final RealComposite< U > targetComp = t.next(); - final RealComposite< L > labelComp = ls.next(); - final RealComposite< M > labelTargetComp = lt.next(); for ( long i = 0; i < size; ++i ) { tmp.get( i ).set( sourceComp.get( i ).getRealDouble() ); } - transformSingleColumnPropagateLabels( tmp, targetComp, labelComp, labelTargetComp, lowerBoundDistanceIndex, envelopeIntersectLocation, d, dim, size ); + transformL1SingleColumn( tmp, targetComp, weight, size ); } } - private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformSingleColumnPropagateLabels( - final RealComposite< T > source, - final RealComposite< U > target, - final RealComposite< L > labelsSource, - final RealComposite< M > labelsResult, - final RealComposite< LongType > lowerBoundDistanceIndex, - final RealComposite< DoubleType > envelopeIntersectLocation, - final Distance d, + private static < T extends RealType< T >, U extends RealType< U > > void transformL1AlongDimensionParallel( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > target, final int dim, - final long size ) + final double weight, + final ExecutorService es, + final int nTasks ) throws InterruptedException, ExecutionException { - long k = 0; - - lowerBoundDistanceIndex.get( 0 ).set( 0 ); - envelopeIntersectLocation.get( 0 ).set( Double.NEGATIVE_INFINITY ); - envelopeIntersectLocation.get( 1 ).set( Double.POSITIVE_INFINITY ); - for ( long position = 1; position < size; ++position ) + int largestDim = getLargestDimension( Views.hyperSlice( target, dim, target.min( dim ) ) ); + // ignore dimension along which we calculate transform + if ( largestDim >= dim ) { - long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - final double sourceAtPosition = source.get( position ).getRealDouble(); - double s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); - - for ( double envelopeValueAtK = envelopeIntersectLocation.get( k ).get(); s <= envelopeValueAtK; envelopeValueAtK = envelopeIntersectLocation.get( k ).get() ) - { - --k; - envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); - } - ++k; - lowerBoundDistanceIndex.get( k ).set( position ); - envelopeIntersectLocation.get( k ).set( s ); - envelopeIntersectLocation.get( k + 1 ).set( Double.POSITIVE_INFINITY ); + largestDim += 1; } - - k = 0; - for ( long position = 0; position < size; ++position ) - { - while ( envelopeIntersectLocation.get( k + 1 ).get() < position ) - { - ++k; - } - final long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); - // copy necessary because of the following line, access to source - // after write to source -> source and target cannot be the same - target.get( position ).setReal( d.evaluate( position, envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), dim ) ); - labelsResult.get( position ).setInteger( labelsSource.get( envelopeIndexAtK ).getIntegerLong() ); - } - - } - - private static < T extends RealType< T >, U extends RealType< U > > void transformL1AlongDimension( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > target, - final int dim, - final double weight ) - { - final int lastDim = target.numDimensions() - 1; - final long size = target.dimension( dim ); - final RealComposite< DoubleType > tmp = Views.collapseReal( createAppropriateOneDimensionalImage( size, new DoubleType() ) ).randomAccess().get(); - // do not permute if we already work on last dimension - final Cursor< RealComposite< T > > s = Views.flatIterable( Views.collapseReal( dim == lastDim ? Views.interval( source, target ) : Views.permute( Views.interval( source, target ), dim, lastDim ) ) ).cursor(); - final Cursor< RealComposite< U > > t = Views.flatIterable( Views.collapseReal( dim == lastDim ? target : Views.permute( target, dim, lastDim ) ) ).cursor(); - - while ( s.hasNext() ) - { - final RealComposite< T > sourceComp = s.next(); - final RealComposite< U > targetComp = t.next(); - for ( long i = 0; i < size; ++i ) - { - tmp.get( i ).set( sourceComp.get( i ).getRealDouble() ); - } - transformL1SingleColumn( tmp, targetComp, weight, size ); - } - } - - private static < T extends RealType< T >, U extends RealType< U > > void transformL1AlongDimensionParallel( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > target, - final int dim, - final double weight, - final ExecutorService es, - final int nTasks ) throws InterruptedException, ExecutionException - { - int largestDim = getLargestDimension( Views.hyperSlice( target, dim, target.min( dim ) ) ); - // ignore dimension along which we calculate transform - if ( largestDim >= dim ) - { - largestDim += 1; - } - final long size = target.dimension( dim ); - final long stepPerChunk = Math.max( size / nTasks, 1 ); + final long size = target.dimension( dim ); + final long stepPerChunk = Math.max( size / nTasks, 1 ); final long[] min = Intervals.minAsLongArray( target ); final long[] max = Intervals.maxAsLongArray( target ); @@ -1626,6 +1327,395 @@ private static < T extends NativeType< T > & RealType< T > > Img< T > createAppr return size > Integer.MAX_VALUE ? new CellImgFactory<>( t, Integer.MAX_VALUE ).create( dim ) : new ArrayImgFactory<>( t ).create( dim ); } + @SuppressWarnings( "unchecked" ) + public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( + final RandomAccessibleInterval< L > labels, + final long backgroundLabel, + final double... weights ) + { + final RandomAccessibleInterval< DoubleType > distance = makeDistances( backgroundLabel, labels, new DoubleType() ); + labelTransform( labels, distance, weights ); + return ( RandomAccessibleInterval< T > ) distance; + } + + @SuppressWarnings( "unchecked" ) + public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( + final RandomAccessibleInterval< L > labels, + final long backgroundLabel, + final ExecutorService es, + final int nTasks, + final double... weights ) + { + final RandomAccessibleInterval< DoubleType > distance = makeDistances( backgroundLabel, labels, new DoubleType() ); + labelTransform( labels, distance, weights ); + return ( RandomAccessibleInterval< T > ) distance; + } + + /** + * Compute the distance transform, of the given distance image, storing the + * result in place. Simultaneously, modifies the labels image so that the + * value at every point is the closest label, where initial distances are + * given by the values in the distance argument. + * + * @param + * label type + * @param + * distance type + * @param labels + * the label iamge + * @param distance + * the distance image + * @param weights + * for distance computation per dimension + */ + public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( + final RandomAccessibleInterval< L > labels, + final RandomAccessibleInterval< T > distance, + final double... weights ) + { + final T maxVal = distance.getType().copy(); + maxVal.setReal( maxVal.getMaxValue() ); + final boolean isIsotropic = weights.length <= 1; + final double[] w = weights.length == labels.numDimensions() ? weights : DoubleStream.generate( () -> weights.length == 0 ? 1.0 : weights[ 0 ] ).limit( labels.numDimensions() ).toArray(); + + // TODO add L1 distance + final Distance distanceFun = isIsotropic ? new EuclidianDistanceIsotropic( w[ 0 ] ) : new EuclidianDistanceAnisotropic( w ); + transformPropagateLabels( distance, distance, distance, labels, labels, distanceFun ); + } + + public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( + final RandomAccessibleInterval< L > labels, + final RandomAccessibleInterval< T > distance, + final ExecutorService es, + final int nTasks, + final double... weights ) throws InterruptedException, ExecutionException + { + final T maxVal = distance.getType().copy(); + maxVal.setReal( maxVal.getMaxValue() ); + final boolean isIsotropic = weights.length <= 1; + final double[] w = weights.length == labels.numDimensions() ? weights : DoubleStream.generate( () -> weights.length == 0 ? 1.0 : weights[ 0 ] ).limit( labels.numDimensions() ).toArray(); + + // TODO add L1 distance + final Distance distanceFun = isIsotropic ? new EuclidianDistanceIsotropic( w[ 0 ] ) : new EuclidianDistanceAnisotropic( w ); + transformPropagateLabels( distance, distance, distance, labels, labels, distanceFun, es, nTasks ); + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate and final results will be stored in + * {@code source} ({@link DoubleType} recommended). + *

+ * Also propagates a set of input labels so that it contains the nearest label + * everywhere. + * + * @param source + * Input function on which distance transform should be computed. + * @param labels + * Labels to be propagated. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input + * @param + * {@link IntegerType} the label type + */ + public static < T extends RealType< T >, L extends IntegerType< L > > void transformPropagateLabels( + final RandomAccessibleInterval< T > distance, + final RandomAccessibleInterval< L > labels, + final DISTANCE_TYPE distanceType ) + { + // TODO generalize distance + transformPropagateLabels( distance, distance, labels, labels, new EuclidianDistanceIsotropic( 1 ) ); + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate and final results will be stored in + * {@code target} ({@link DoubleType} recommended). + * + * Simultaneously, modifies the labels image so that the label at every + * point is the closest label, where initial distances are given by the + * values in the distance argument. + */ + public static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( + final RandomAccessible< T > distance, + final RandomAccessibleInterval< U > distanceResult, + final RandomAccessible< L > labels, + final RandomAccessibleInterval< M > labelsResult, + final Distance d ) + { + transformPropagateLabels( distance, distanceResult, distanceResult, labels, labelsResult, d ); + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate results will be stored in {@code tmp} + * ({@link DoubleType} recommended). The output will be written into + * {@code target}. + * + * @param source + * Input function on which distance transform should be computed. + * @param tmp + * Storage for intermediate results. + * @param target + * Final result of distance transform. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input + * @param + * {@link RealType} intermediate results + * @param + * {@link RealType} output + */ + public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > tmp, + final RandomAccessibleInterval< V > target, + final RandomAccessible< L > labels, + final RandomAccessibleInterval< M > labelsResult, + final Distance d ) + { + assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; + final int nDim = source.numDimensions(); + final int lastDim = nDim - 1; + + if ( nDim == 1 ) + { + transformAlongDimensionPropagateLabels( + ( RandomAccessible< T > ) Views.addDimension( source ), + Views.addDimension( target, 0, 0 ), + Views.addDimension( labels ), + Views.addDimension( labelsResult ), + d, 0 ); + } + else + { + transformAlongDimensionPropagateLabels( source, tmp, labels, labelsResult, d, 0 ); + } + + for ( int dim = 1; dim < nDim; ++dim ) + { + if ( dim == lastDim ) + { + transformAlongDimensionPropagateLabels( tmp, target, labels, labelsResult, d, dim ); + } + else + { + transformAlongDimensionPropagateLabels( tmp, tmp, labels, labelsResult, d, dim ); + } + } + } + + /** + * Create + * distance + * transforms of sampled functions on {@code source} using arbitrary + * {@link Distance} d. Intermediate results will be stored in {@code tmp} + * ({@link DoubleType} recommended). The output will be written into + * {@code target}. + * + * @param source + * Input function on which distance transform should be computed. + * @param tmp + * Storage for intermediate results. + * @param target + * Final result of distance transform. + * @param d + * {@link Distance} between two points. + * @param + * {@link RealType} input + * @param + * {@link RealType} intermediate results + * @param + * {@link RealType} output + * @throws ExecutionException + * @throws InterruptedException + */ + public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > tmp, + final RandomAccessibleInterval< V > target, + final RandomAccessible< L > labels, + final RandomAccessibleInterval< M > labelsResult, + final Distance d, + final ExecutorService es, + final int nTasks) throws InterruptedException, ExecutionException + { + assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; + final int nDim = source.numDimensions(); + final int lastDim = nDim - 1; + + if ( nDim == 1 ) + { + transformAlongDimensionPropagateLabels( + ( RandomAccessible< T > ) Views.addDimension( source ), + Views.addDimension( target, 0, 0 ), + Views.addDimension( labels ), + Views.addDimension( labelsResult ), + d, 0 ); + } + else + { + transformAlongDimensionPropagateLabelsParallel( source, tmp, labels, labelsResult, d, 0, es, nTasks ); + } + + for ( int dim = 1; dim < nDim; ++dim ) + { + if ( dim == lastDim ) + { + transformAlongDimensionPropagateLabelsParallel( tmp, target, labels, labelsResult, d, dim, es, nTasks); + } + else + { + transformAlongDimensionPropagateLabelsParallel( tmp, tmp, labels, labelsResult, d, dim, es, nTasks); + } + } + } + + private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType< L >, M extends IntegerType< M > > void transformAlongDimensionPropagateLabels( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > target, + final RandomAccessible< L > labelSource, + final RandomAccessible< M > labelTarget, + final Distance d, + final int dim ) + { + final int lastDim = target.numDimensions() - 1; + final long size = target.dimension( dim ); + + final Img< DoubleType > tmpImg = createAppropriateOneDimensionalImage( size, new DoubleType() ); + final RealComposite< DoubleType > tmp = Views.collapseReal( tmpImg ).randomAccess().get(); + + final Img< L > tmpLabelImg = Util.getSuitableImgFactory( tmpImg, labelSource.getType() ).create( tmpImg ); + final RealComposite< L > tmpLabel = Views.collapseReal( tmpLabelImg ).randomAccess().get(); + + // do not permute if we already work on last dimension + final Cursor< RealComposite< T > > s = Views.flatIterable( Views.collapseReal( dim == lastDim ? Views.interval( source, target ) : Views.permute( Views.interval( source, target ), dim, lastDim ) ) ).cursor(); + final Cursor< RealComposite< U > > t = Views.flatIterable( Views.collapseReal( dim == lastDim ? target : Views.permute( target, dim, lastDim ) ) ).cursor(); + + final Cursor< RealComposite< L > > ls = Views.flatIterable( + Views.collapseReal( dim == lastDim ? Views.interval( labelSource, target ) : Views.permute( Views.interval( labelSource, target ), dim, lastDim ) ) ).cursor(); + + final Cursor< RealComposite< M > > lt = Views.flatIterable( + Views.collapseReal( dim == lastDim ? Views.interval( labelTarget, target ) : Views.permute( Views.interval( labelTarget, target ), dim, lastDim ) ) ).cursor(); + + final RealComposite< LongType > lowerBoundDistanceIndex = Views.collapseReal( createAppropriateOneDimensionalImage( size, new LongType() ) ).randomAccess().get(); + final RealComposite< DoubleType > envelopeIntersectLocation = Views.collapseReal( createAppropriateOneDimensionalImage( size + 1, new DoubleType() ) ).randomAccess().get(); + + while ( s.hasNext() ) + { + final RealComposite< T > sourceComp = s.next(); + final RealComposite< U > targetComp = t.next(); + final RealComposite< L > labelComp = ls.next(); + final RealComposite< M > labelTargetComp = lt.next(); + for ( long i = 0; i < size; ++i ) + { + tmp.get( i ).set( sourceComp.get( i ).getRealDouble() ); + tmpLabel.get( i ).setInteger( labelComp.get( i ).getIntegerLong() ); + } + transformSingleColumnPropagateLabels( tmp, targetComp, tmpLabel, labelTargetComp, lowerBoundDistanceIndex, envelopeIntersectLocation, d, dim, size ); + } + } + + private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformSingleColumnPropagateLabels( + final RealComposite< T > source, + final RealComposite< U > target, + final RealComposite< L > labelsSource, + final RealComposite< M > labelsResult, + final RealComposite< LongType > lowerBoundDistanceIndex, + final RealComposite< DoubleType > envelopeIntersectLocation, + final Distance d, + final int dim, + final long size ) + { + long k = 0; + + lowerBoundDistanceIndex.get( 0 ).set( 0 ); + envelopeIntersectLocation.get( 0 ).set( Double.NEGATIVE_INFINITY ); + envelopeIntersectLocation.get( 1 ).set( Double.POSITIVE_INFINITY ); + for ( long position = 1; position < size; ++position ) + { + long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + final double sourceAtPosition = source.get( position ).getRealDouble(); + double s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + + for ( double envelopeValueAtK = envelopeIntersectLocation.get( k ).get(); s <= envelopeValueAtK; envelopeValueAtK = envelopeIntersectLocation.get( k ).get() ) + { + --k; + envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + s = d.intersect( envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), position, sourceAtPosition, dim ); + } + ++k; + lowerBoundDistanceIndex.get( k ).set( position ); + envelopeIntersectLocation.get( k ).set( s ); + envelopeIntersectLocation.get( k + 1 ).set( Double.POSITIVE_INFINITY ); + } + + k = 0; + for ( long position = 0; position < size; ++position ) + { + while ( envelopeIntersectLocation.get( k + 1 ).get() < position ) + { + ++k; + } + final long envelopeIndexAtK = lowerBoundDistanceIndex.get( k ).get(); + // copy necessary because of the following line, access to source + // after write to source -> source and target cannot be the same + target.get( position ).setReal( d.evaluate( position, envelopeIndexAtK, source.get( envelopeIndexAtK ).getRealDouble(), dim ) ); + labelsResult.get( position ).setInteger( labelsSource.get( envelopeIndexAtK ).getIntegerLong() ); + } + + } + + private static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType< L >, M extends IntegerType< M > > void transformAlongDimensionPropagateLabelsParallel( + final RandomAccessible< T > source, + final RandomAccessibleInterval< U > target, + final RandomAccessible< L > labelSource, + final RandomAccessible< M > labelTarget, + final Distance d, + final int dim, + final ExecutorService es, + final int nTasks ) throws InterruptedException, ExecutionException + { + int largestDim = getLargestDimension( Views.hyperSlice( target, dim, target.min( dim ) ) ); + // ignore dimension along which we calculate transform + if ( largestDim >= dim ) + { + largestDim += 1; + } + final long size = target.dimension( dim ); + final long stepPerChunk = Math.max( size / nTasks, 1 ); + + final long[] min = Intervals.minAsLongArray( target ); + final long[] max = Intervals.maxAsLongArray( target ); + + final long largestDimMin = target.min( largestDim ); + final long largestDimMax = target.max( largestDim ); + + final ArrayList< Callable< Void > > tasks = new ArrayList<>(); + for ( long m = largestDimMin, M = largestDimMin + stepPerChunk - 1; m <= largestDimMax; m += stepPerChunk, M += stepPerChunk ) + { + min[ largestDim ] = m; + max[ largestDim ] = Math.min( M, largestDimMax ); + final Interval fi = new FinalInterval( min, max ); + tasks.add( () -> { + transformAlongDimensionPropagateLabels( source, Views.interval( target, fi ), + labelSource, Views.interval( labelTarget, fi ), + d, dim ); + return null; + } ); + } + + invokeAllAndWait( es, tasks ); + } + /** * Convenience method to find largest dimension of {@link Interval} * interval. @@ -1640,18 +1730,6 @@ public static int getLargestDimension( final Interval interval ) return IntStream.range( 0, interval.numDimensions() ).mapToObj( i -> new ValuePair<>( i, interval.dimension( i ) ) ).max( ( p1, p2 ) -> Long.compare( p1.getB(), p2.getB() ) ).get().getA(); } - private static , M extends IntegerType< M > > RandomAccessibleInterval< M > makeLabels( - final RandomAccessibleInterval< L > labels, - final M type) { - - final RandomAccessibleInterval< M > propagedLabels = Util.getSuitableImgFactory( labels, type ).create( labels ); - LoopBuilder.setImages( labels, propagedLabels ).forEachPixel( (s,d) -> { - d.setInteger( s.getIntegerLong() ); - }); - - return ( RandomAccessibleInterval< M > ) propagedLabels; - } - private static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > makeDistances( final long label, final RandomAccessibleInterval< L > labels, @@ -1688,28 +1766,4 @@ public void convert( final B input, final R output ) } - private static class LabelsToCost< L extends IntegerType< L >, R extends RealType< R > > implements Converter< L, R > - { - private final R maxValForR; - - private final R zero; - - private final long label; - - public LabelsToCost( final long label, final R maxValForR ) - { - this.maxValForR = maxValForR; - this.label = label; - this.zero = maxValForR.createVariable(); - zero.setZero(); - } - - @Override - public void convert( final L input, final R output ) - { - output.set( input.getIntegerLong() == label ? zero : maxValForR ); - } - - } - } From 1a9958626ca22a2d06d10bc3d87d5a22a19c0077 Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Thu, 27 Feb 2025 15:44:40 -0500 Subject: [PATCH 4/7] test: label propagation with DistanceTransform --- .../distance/DistanceTransformTest.java | 191 +++++++++++++++++- 1 file changed, 186 insertions(+), 5 deletions(-) diff --git a/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java b/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java index ff623beaf..f73668642 100644 --- a/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java +++ b/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java @@ -34,8 +34,14 @@ package net.imglib2.algorithm.morphology.distance; +import static org.junit.Assert.assertTrue; + +import java.util.ArrayList; import java.util.Arrays; +import java.util.HashSet; +import java.util.List; import java.util.Random; +import java.util.Set; import java.util.concurrent.ExecutionException; import java.util.concurrent.ExecutorService; import java.util.concurrent.Executors; @@ -46,6 +52,7 @@ import org.junit.Test; import net.imglib2.Cursor; +import net.imglib2.Interval; import net.imglib2.Localizable; import net.imglib2.Point; import net.imglib2.RandomAccess; @@ -58,16 +65,20 @@ import net.imglib2.img.basictypeaccess.array.DoubleArray; import net.imglib2.img.basictypeaccess.array.LongArray; import net.imglib2.type.logic.BitType; +import net.imglib2.type.numeric.IntegerType; import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.LongType; import net.imglib2.type.numeric.real.DoubleType; import net.imglib2.util.Intervals; +import net.imglib2.util.Localizables; import net.imglib2.util.Pair; +import net.imglib2.util.Util; import net.imglib2.view.Views; /** * * @author Philipp Hanslovsky - * + * @author John Bogovic */ public class DistanceTransformTest { @@ -159,9 +170,7 @@ private void testBinary( final DISTANCE_TYPE dt, final DistanceCalculator distan private static void compareRAIofRealType( final RandomAccessibleInterval< ? extends RealType< ? > > ref, final RandomAccessibleInterval< ? extends RealType< ? > > comp, final double tolerance ) { - Assert.assertArrayEquals( Intervals.dimensionsAsLongArray( ref ), Intervals.dimensionsAsLongArray( comp ) ); - Assert.assertArrayEquals( Intervals.minAsLongArray( ref ), Intervals.minAsLongArray( comp ) ); - Assert.assertArrayEquals( Intervals.maxAsLongArray( ref ), Intervals.maxAsLongArray( comp ) ); + assertTrue( Intervals.equals( ref, comp ) ); for ( final Pair< ? extends RealType< ? >, ? extends RealType< ? > > p : Views.flatIterable( Views.interval( Views.pair( ref, comp ), ref ) ) ) { Assert.assertEquals( p.getA().getRealDouble(), p.getB().getRealDouble(), tolerance ); @@ -440,7 +449,7 @@ private static < T extends RealType< T > > void checkDistance( final double[] weights, final DistanceCalculator distanceCalculator ) { - for ( final Cursor< T > c = Views.iterable( dist ).localizingCursor(); c.hasNext(); ) + for ( final Cursor< T > c = dist.localizingCursor(); c.hasNext(); ) { final double actual = c.next().getRealDouble(); final double expected = atSamePosition( foreground, c ) ? 0.0 : distanceCalculator.dist( foreground, c, weights ); @@ -448,4 +457,176 @@ private static < T extends RealType< T > > void checkDistance( } } + @Test + public void testLabelPropagation() + { + /* + * Iterate over numReplicates = [0..9] numDimensions = [2, 3] numLabels + * = [1..5] + */ + final int firstReplicate = 0; + final int lastReplicate = 9; + + final int firstNumDimensions = 2; + final int lastNumDimensions = 3; + + final int firstNumLabels = 2; + final int lastNumLabels = 5; + + final RandomAccessibleInterval< Localizable > parameters = Localizables.randomAccessibleInterval( + Intervals.createMinMax( + firstReplicate, firstNumDimensions, + firstNumLabels, lastReplicate, + lastNumDimensions, lastNumLabels ) ); + + parameters.forEach( params -> { + + @SuppressWarnings( "unused" ) + final int replicate = params.getIntPosition( 0 ); + final int numDimensions = params.getIntPosition( 1 ); + final int numLabels = params.getIntPosition( 2 ); + + testLabelPropagationHelper( numDimensions, numLabels ); + testLabelPropagationHelperParallel( numDimensions, numLabels ); + } ); + + } + + /** + * Creates an label and distances images with the requested number of dimensions (ndims), + * and places nLabels points with non-zero label. Checks that the propagated labels correctly + * reflect the nearest label (ties are allowed: any label equi-distant to a point passes). + * + * @param ndims number of dimensions + * @param nLabels number of labels + */ + private void testLabelPropagationHelper( int ndims, int nLabels ) + { + + final long[] imgDims = LongStream.iterate( dimensionSize, d -> d - 1 ).limit( ndims ).toArray(); + final ArrayImg< LongType, LongArray > labels = ArrayImgs.longs( imgDims ); + + final Set< PointAndLabel > points = initializeLabels( rng, nLabels, labels ); + DistanceTransform.labelTransform( labels, 0 ); + validateLabelsSet( "serial", points, labels ); + } + + /** + * Creates an label and distances images with the requested number of dimensions (ndims), + * and places nLabels points with non-zero label. Checks that the propagated labels correctly + * reflect the nearest label (ties are allowed: any label equi-distant to a point passes). + * + * @param ndims number of dimensions + * @param nLabels number of labels + */ + private void testLabelPropagationHelperParallel( int ndims, int nLabels ) + { + + final long[] imgDims = LongStream.iterate( dimensionSize, d -> d - 1 ).limit( ndims ).toArray(); + final ArrayImg< LongType, LongArray > labels = ArrayImgs.longs( imgDims ); + final Set< PointAndLabel > points = initializeLabels( rng, nLabels, labels ); + DistanceTransform.labelTransform( labels, 0, es, 3 * nThreads ); + validateLabelsSet( "parallel", points, labels ); + } + + private ArrayImg< LongType, LongArray > copyLongArrayImg( ArrayImg< LongType, LongArray > img ) + { + + final long[] dataOrig = img.getAccessType().getCurrentStorageArray(); + final long[] dataCopy = new long[ dataOrig.length ]; + System.arraycopy( dataOrig, 0, dataCopy, 0, dataOrig.length ); + return ArrayImgs.longs( dataCopy, img.dimensionsAsLongArray() ); + } + + private static Point randomPointInInterval( final Random rng, final Interval itvl ) + { + final int[] coords = IntStream.range( 0, itvl.numDimensions() ).map( i -> { + return rng.nextInt( ( int ) itvl.dimension( i ) ); + } ).toArray(); + return new Point( coords ); + } + + private static < T extends RealType< T >, L extends IntegerType< L > > Set< PointAndLabel > initializeLabels( Random random, int numLabels, RandomAccessibleInterval< L > labels ) + { + labels.forEach( p -> p.setZero() ); // Initialize all labels to 0 + Set< PointAndLabel > positions = new HashSet<>(); + + int currentLabel = 1; + // Set numLabels different random positions to a non-zero label + while ( positions.size() < numLabels ) + { + final Point pt = randomPointInInterval( random, labels ); + if ( !positions.contains( pt ) ) + { + + final PointAndLabel candidate = new PointAndLabel( currentLabel, pt.positionAsLongArray() ); + if ( !positions.contains( candidate ) ) + { + positions.add( candidate ); + labels.randomAccess().setPositionAndGet( pt ).setInteger( currentLabel ); + currentLabel++; + } + + } + } + return positions; + } + + /** + * Return the set of points within epsilon distance of the query point + * + * @param query point + * @param pointSet set of candidate points + * @param epsilon distance threshold + * @return the set of close points + */ + private static List< PointAndLabel > closestSet( Localizable query, Set< PointAndLabel > pointSet, final double epsilon ) + { + + final List< PointAndLabel > listOfEquidistant = new ArrayList<>(); + + double mindist = Double.MAX_VALUE; + for ( PointAndLabel pt : pointSet ) + { + double dist = Util.distance( query, pt ); + + if ( Math.abs( dist - mindist ) < epsilon ) + { + listOfEquidistant.add( pt ); + } + else if ( dist < mindist ) + { + mindist = dist; + listOfEquidistant.clear(); + listOfEquidistant.add( pt ); + } + } + + return listOfEquidistant; + } + + private static < T extends RealType< T >, L extends IntegerType< L > > void validateLabelsSet( final String prefix, final Set< PointAndLabel > points, final RandomAccessibleInterval< L > labels ) + { + final double EPS = 0.01; + final Cursor< L > c = labels.cursor(); + while ( c.hasNext() ) + { + c.fwd(); + final boolean labelIsClosest = closestSet( c, points, EPS ).stream().anyMatch( p -> p.label == c.get().getIntegerLong() ); + assertTrue( prefix + " point: " + Arrays.toString( c.positionAsLongArray() ), labelIsClosest ); + } + } + + private static class PointAndLabel extends Point + { + + long label; + + public PointAndLabel( long label, long[] position ) + { + super( position ); + this.label = label; + } + } + } From c2bab84d547d50dc120c31612c4c385d676fd57a Mon Sep 17 00:00:00 2001 From: John Bogovic Date: Fri, 28 Feb 2025 11:48:40 -0500 Subject: [PATCH 5/7] doc and rename labelTransform to voronoiDistanceTransform * better javadoc * create and use a shared createEuclideanDistance method * makeDistances sets maxVal of the given type * no L1 distance option for now --- .../distance/DistanceTransform.java | 364 +++++++++++++----- .../distance/DistanceTransformTest.java | 19 +- 2 files changed, 274 insertions(+), 109 deletions(-) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index 1546afc08..8ecdbf083 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -1327,28 +1327,102 @@ private static < T extends NativeType< T > & RealType< T > > Img< T > createAppr return size > Integer.MAX_VALUE ? new CellImgFactory<>( t, Integer.MAX_VALUE ).create( dim ) : new ArrayImgFactory<>( t ).create( dim ); } - @SuppressWarnings( "unchecked" ) - public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( + /** + * Compute the distance and nearest neighbors of a label image. + *

+ * Returns the joint distance transform of all non-background labels in the + * {@code labels} image. This distance transform will be the distance to the + * nearest non-background label at every point. + *

+ * Simultaneously, modifies the {@code labels} image so that the label at + * every point is the closest label, where initial distances are given by + * the values in the distance argument. + *

+ * This method uses 0 (zero) as the background label. + * + * @param + * the label type + * @param labels + * the label image + * @param weights + * for distance computation per dimension + * @return the distance map + */ + public static < L extends IntegerType< L > > RandomAccessibleInterval< DoubleType > voronoiDistanceTransform( + final RandomAccessibleInterval< L > labels, + final double... weights ) + { + return voronoiDistanceTransform( labels, 0, weights ); + } + + /** + * Compute the distance and nearest neighbors of a label image. + *

+ * Returns the joint distance transform of all non-background labels in the + * {@code labels} image. This distance transform will be the distance to the + * nearest non-background label at every point. + *

+ * Simultaneously, modifies the {@code labels} image so that the label at + * every point is the closest label, where initial distances are given by + * the values in the distance argument. + * + * @param + * the label type + * @param labels + * the label image + * @param backgroundLabel + * the background label + * @param weights + * for distance computation per dimension + * @return the distance map + */ + public static < L extends IntegerType< L > > RandomAccessibleInterval< DoubleType > voronoiDistanceTransform( final RandomAccessibleInterval< L > labels, final long backgroundLabel, final double... weights ) { final RandomAccessibleInterval< DoubleType > distance = makeDistances( backgroundLabel, labels, new DoubleType() ); - labelTransform( labels, distance, weights ); - return ( RandomAccessibleInterval< T > ) distance; + voronoiDistanceTransform( labels, distance, weights ); + return distance; } - @SuppressWarnings( "unchecked" ) - public static < L extends IntegerType< L >, T extends RealType< T > > RandomAccessibleInterval< T > labelTransform( + /** + * Compute the distance and nearest neighbors of a label image in parallel. + *

+ * Returns the joint distance transform of all non-background labels in the + * {@code labels} image. This distance transform will be the distance to the + * nearest non-background label at every point. + *

+ * Simultaneously, modifies the {@code labels} image so that the label at + * every point is the closest label, where initial distances are given by + * the values in the distance argument. + * + * @param + * the label type + * @param labels + * the label image + * @param backgroundLabel + * the background label + * @param es + * the ExecutorService + * @param nTasks + * the number of tasks in which to split the computation + * @param weights + * for distance computation per dimension + * @return the distance map + * @throws ExecutionException + * @throws InterruptedException + */ + public static < L extends IntegerType< L > > RandomAccessibleInterval< DoubleType > voronoiDistanceTransform( final RandomAccessibleInterval< L > labels, final long backgroundLabel, final ExecutorService es, final int nTasks, - final double... weights ) + final double... weights ) throws InterruptedException, ExecutionException { final RandomAccessibleInterval< DoubleType > distance = makeDistances( backgroundLabel, labels, new DoubleType() ); - labelTransform( labels, distance, weights ); - return ( RandomAccessibleInterval< T > ) distance; + labelTransform( labels, distance, es, nTasks, weights ); + return distance; } /** @@ -1362,27 +1436,48 @@ public static < L extends IntegerType< L >, T extends RealType< T > > RandomAcce * @param * distance type * @param labels - * the label iamge + * the label image * @param distance * the distance image * @param weights * for distance computation per dimension */ - public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( + public static < L extends IntegerType< L >, T extends RealType< T > > void voronoiDistanceTransform( final RandomAccessibleInterval< L > labels, final RandomAccessibleInterval< T > distance, final double... weights ) { - final T maxVal = distance.getType().copy(); - maxVal.setReal( maxVal.getMaxValue() ); - final boolean isIsotropic = weights.length <= 1; - final double[] w = weights.length == labels.numDimensions() ? weights : DoubleStream.generate( () -> weights.length == 0 ? 1.0 : weights[ 0 ] ).limit( labels.numDimensions() ).toArray(); - - // TODO add L1 distance - final Distance distanceFun = isIsotropic ? new EuclidianDistanceIsotropic( w[ 0 ] ) : new EuclidianDistanceAnisotropic( w ); + final Distance distanceFun = createEuclideanDistance(labels.numDimensions(), weights) ; transformPropagateLabels( distance, distance, distance, labels, labels, distanceFun ); } + /** + * Compute the distance and nearest neighbors of a label image in parallel. + *

+ * Returns the joint distance transform of all non-background labels in the + * {@code labels} image. This distance transform will be the distance to the + * nearest non-background label at every point. + *

+ * Simultaneously, modifies the {@code labels} image so that the label at + * every point is the closest label, where initial distances are given by + * the values in the distance argument. + * + * @param + * the label type + * @param labels + * the label image + * @param backgroundLabel + * the background label + * @param es + * the ExecutorService + * @param nTasks + * the number of tasks in which to split the computation + * @param weights + * for distance computation per dimension + * @return the distance map + * @throws ExecutionException + * @throws InterruptedException + */ public static < L extends IntegerType< L >, T extends RealType< T > > void labelTransform( final RandomAccessibleInterval< L > labels, final RandomAccessibleInterval< T > distance, @@ -1390,32 +1485,28 @@ public static < L extends IntegerType< L >, T extends RealType< T > > void label final int nTasks, final double... weights ) throws InterruptedException, ExecutionException { - final T maxVal = distance.getType().copy(); - maxVal.setReal( maxVal.getMaxValue() ); - final boolean isIsotropic = weights.length <= 1; - final double[] w = weights.length == labels.numDimensions() ? weights : DoubleStream.generate( () -> weights.length == 0 ? 1.0 : weights[ 0 ] ).limit( labels.numDimensions() ).toArray(); - - // TODO add L1 distance - final Distance distanceFun = isIsotropic ? new EuclidianDistanceIsotropic( w[ 0 ] ) : new EuclidianDistanceAnisotropic( w ); + final Distance distanceFun = createEuclideanDistance(labels.numDimensions(), weights) ; transformPropagateLabels( distance, distance, distance, labels, labels, distanceFun, es, nTasks ); } /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate and final results will be stored in - * {@code source} ({@link DoubleType} recommended). + * Computes the distance transform of the input distance map {@code distance} + * and simultaneously propagate a set of corresponding {@code labels}. + * Distance results are stored in the input {@code targetDistance} image. *

- * Also propagates a set of input labels so that it contains the nearest label - * everywhere. + * Simultaneously, propagates labels stored in the {@code labels} image so + * that the label at every point is the closest label, where initial + * distances are given by the values in the distance argument. Results of + * label propagation are stored in {@code labelsResult} + *

+ * This implementation operates in-place for both the {@code distance} and + * {@code labels} images. It uses an isotropic distance function with distance 1 + * between samples. * - * @param source + * @param distance * Input function on which distance transform should be computed. * @param labels * Labels to be propagated. - * @param d - * {@link Distance} between two points. * @param * {@link RealType} input * @param @@ -1423,158 +1514,211 @@ public static < L extends IntegerType< L >, T extends RealType< T > > void label */ public static < T extends RealType< T >, L extends IntegerType< L > > void transformPropagateLabels( final RandomAccessibleInterval< T > distance, - final RandomAccessibleInterval< L > labels, - final DISTANCE_TYPE distanceType ) + final RandomAccessibleInterval< L > labels) { - // TODO generalize distance transformPropagateLabels( distance, distance, labels, labels, new EuclidianDistanceIsotropic( 1 ) ); } /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate and final results will be stored in - * {@code target} ({@link DoubleType} recommended). + * Computes the distance transform of the input distance map {@code distance} + * and simultaneously propagate a set of corresponding {@code labels}. + * Distance results are stored in the input {@code targetDistance} image. + *

+ * Simultaneously, propagates labels stored in the {@code labels} image so + * that the label at every point is the closest label, where initial + * distances are given by the values in the distance argument. Results of + * label propagation are stored in {@code labelsResult} + *

+ * This implementation uses the {@targetDistance} for temporary storage. * - * Simultaneously, modifies the labels image so that the label at every - * point is the closest label, where initial distances are given by the - * values in the distance argument. + * @param + * input distance {@link RealType} + * @param + * output distance {@link RealType} + * @param + * input label {@link IntegerType} + * @param + * output label {@link IntegerType} + * @param distance + * Input distance function from which distance transform should + * be computed. + * @param tmpDistance + * Storage for intermediate distance results. + * @param targetDistance + * Final result of distance transform. May be the same instance + * as distance. + * @param labels + * the image of labels to be propagated + * @param labelsResult + * the image in which to store the result of label propagation. + * May be the same instance as labels + * @param d + * the {@link Distance} function */ public static < T extends RealType< T >, U extends RealType< U >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( final RandomAccessible< T > distance, - final RandomAccessibleInterval< U > distanceResult, + final RandomAccessibleInterval< U > targetDistance, final RandomAccessible< L > labels, final RandomAccessibleInterval< M > labelsResult, final Distance d ) { - transformPropagateLabels( distance, distanceResult, distanceResult, labels, labelsResult, d ); + transformPropagateLabels( distance, targetDistance, targetDistance, labels, labelsResult, d ); } /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate results will be stored in {@code tmp} - * ({@link DoubleType} recommended). The output will be written into - * {@code target}. + * Computes the distance transform of the input distance map {@code distance} + * and simultaneously propagate a set of corresponding {@code labels}. + * Distance results are stored in the input {@code targetDistance} image. + *

+ * Simultaneously, propagates labels stored in the {@code labels} image so + * that the label at every point is the closest label, where initial + * distances are given by the values in the distance argument. Results of + * label propagation are stored in {@code labelsResult} * - * @param source - * Input function on which distance transform should be computed. - * @param tmp - * Storage for intermediate results. - * @param target - * Final result of distance transform. - * @param d - * {@link Distance} between two points. * @param - * {@link RealType} input + * input distance {@link RealType} * @param - * {@link RealType} intermediate results + * intermediate distance result {@link RealType} * @param - * {@link RealType} output + * output distance {@link RealType} + * @param + * input label {@link IntegerType} + * @param + * output label {@link IntegerType} + * @param distance + * Input distance function from which distance transform should + * be computed. + * @param tmpDistance + * Storage for intermediate distance results. + * @param targetDistance + * Final result of distance transform. May be the same instance + * as distance. + * @param labels + * the image of labels to be propagated + * @param labelsResult + * the image in which to store the result of label propagation. + * May be the same instance as labels + * @param d + * {@link Distance} between two points. */ public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > tmp, - final RandomAccessibleInterval< V > target, + final RandomAccessible< T > distance, + final RandomAccessibleInterval< U > tmpDistance, + final RandomAccessibleInterval< V > targetDistance, final RandomAccessible< L > labels, final RandomAccessibleInterval< M > labelsResult, final Distance d ) { - assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; - final int nDim = source.numDimensions(); + assert distance.numDimensions() == targetDistance.numDimensions(): "Dimension mismatch"; + final int nDim = distance.numDimensions(); final int lastDim = nDim - 1; if ( nDim == 1 ) { transformAlongDimensionPropagateLabels( - ( RandomAccessible< T > ) Views.addDimension( source ), - Views.addDimension( target, 0, 0 ), + ( RandomAccessible< T > ) Views.addDimension( distance ), + Views.addDimension( targetDistance, 0, 0 ), Views.addDimension( labels ), Views.addDimension( labelsResult ), d, 0 ); } else { - transformAlongDimensionPropagateLabels( source, tmp, labels, labelsResult, d, 0 ); + transformAlongDimensionPropagateLabels( distance, tmpDistance, labels, labelsResult, d, 0 ); } for ( int dim = 1; dim < nDim; ++dim ) { if ( dim == lastDim ) { - transformAlongDimensionPropagateLabels( tmp, target, labels, labelsResult, d, dim ); + transformAlongDimensionPropagateLabels( tmpDistance, targetDistance, labels, labelsResult, d, dim ); } else { - transformAlongDimensionPropagateLabels( tmp, tmp, labels, labelsResult, d, dim ); + transformAlongDimensionPropagateLabels( tmpDistance, tmpDistance, labels, labelsResult, d, dim ); } } } /** - * Create - * distance - * transforms of sampled functions on {@code source} using arbitrary - * {@link Distance} d. Intermediate results will be stored in {@code tmp} - * ({@link DoubleType} recommended). The output will be written into - * {@code target}. + * In parallel, computes the distance transform of the input distance map {@code distance} + * and simultaneously propagate a set of corresponding {@code labels}. + * Distance results are stored in the input {@code targetDistance} image. + *

+ * Simultaneously, propagates labels stored in the {@code labels} image so + * that the label at every point is the closest label, where initial + * distances are given by the values in the distance argument. Results of + * label propagation are stored in {@code labelsResult} * - * @param source - * Input function on which distance transform should be computed. - * @param tmp - * Storage for intermediate results. - * @param target - * Final result of distance transform. - * @param d - * {@link Distance} between two points. * @param - * {@link RealType} input + * input distance {@link RealType} * @param - * {@link RealType} intermediate results + * intermediate distance result {@link RealType} * @param - * {@link RealType} output - * @throws ExecutionException - * @throws InterruptedException + * output distance {@link RealType} + * @param + * input label {@link IntegerType} + * @param + * output label {@link IntegerType} + * @param distance + * Input distance function from which distance transform should + * be computed. + * @param tmpDistance + * Storage for intermediate distance results. + * @param targetDistance + * Final result of distance transform. May be the same instance + * as distance. + * @param labels + * the image of labels to be propagated + * @param labelsResult + * the image in which to store the result of label propagation. + * May be the same instance as labels + * @param d + * {@link Distance} between two points. + * @param es + * the ExecutorService + * @param nTasks + * the number of tasks in which to split the computation + * @throws ExecutionException + * @throws InterruptedException */ public static < T extends RealType< T >, U extends RealType< U >, V extends RealType< V >, L extends IntegerType, M extends IntegerType > void transformPropagateLabels( - final RandomAccessible< T > source, - final RandomAccessibleInterval< U > tmp, - final RandomAccessibleInterval< V > target, + final RandomAccessible< T > distance, + final RandomAccessibleInterval< U > tmpDistance, + final RandomAccessibleInterval< V > targetDistance, final RandomAccessible< L > labels, final RandomAccessibleInterval< M > labelsResult, final Distance d, final ExecutorService es, final int nTasks) throws InterruptedException, ExecutionException { - assert source.numDimensions() == target.numDimensions(): "Dimension mismatch"; - final int nDim = source.numDimensions(); + assert distance.numDimensions() == targetDistance.numDimensions(): "Dimension mismatch"; + final int nDim = distance.numDimensions(); final int lastDim = nDim - 1; if ( nDim == 1 ) { transformAlongDimensionPropagateLabels( - ( RandomAccessible< T > ) Views.addDimension( source ), - Views.addDimension( target, 0, 0 ), + ( RandomAccessible< T > ) Views.addDimension( distance ), + Views.addDimension( targetDistance, 0, 0 ), Views.addDimension( labels ), Views.addDimension( labelsResult ), d, 0 ); } else { - transformAlongDimensionPropagateLabelsParallel( source, tmp, labels, labelsResult, d, 0, es, nTasks ); + transformAlongDimensionPropagateLabelsParallel( distance, tmpDistance, labels, labelsResult, d, 0, es, nTasks ); } for ( int dim = 1; dim < nDim; ++dim ) { if ( dim == lastDim ) { - transformAlongDimensionPropagateLabelsParallel( tmp, target, labels, labelsResult, d, dim, es, nTasks); + transformAlongDimensionPropagateLabelsParallel( tmpDistance, targetDistance, labels, labelsResult, d, dim, es, nTasks); } else { - transformAlongDimensionPropagateLabelsParallel( tmp, tmp, labels, labelsResult, d, dim, es, nTasks); + transformAlongDimensionPropagateLabelsParallel( tmpDistance, tmpDistance, labels, labelsResult, d, dim, es, nTasks); } } } @@ -1716,6 +1860,13 @@ private static < T extends RealType< T >, U extends RealType< U >, L extends Int invokeAllAndWait( es, tasks ); } + private static Distance createEuclideanDistance( int numDimensions, double... weights ) + { + final boolean isIsotropic = weights.length <= 1; + final double[] w = weights.length == numDimensions ? weights : DoubleStream.generate( () -> weights.length == 0 ? 1.0 : weights[ 0 ] ).limit( numDimensions ).toArray(); + return isIsotropic ? new EuclidianDistanceIsotropic( w[ 0 ] ) : new EuclidianDistanceAnisotropic( w ); + } + /** * Convenience method to find largest dimension of {@link Interval} * interval. @@ -1735,10 +1886,13 @@ private static < L extends IntegerType< L >, T extends RealType< T > > RandomAcc final RandomAccessibleInterval< L > labels, final T type) { + final T maxVal = type.copy(); + maxVal.setReal( maxVal.getMaxValue() ); + final RandomAccessibleInterval< T > distances = Util.getSuitableImgFactory( labels, type ).create( labels ); Views.pair( labels, distances ).view().interval( labels ).forEach( pair -> { if( pair.getA().getIntegerLong() == label ) - pair.getB().setReal( Double.MAX_VALUE ); + pair.getB().set( maxVal ); }); return distances; diff --git a/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java b/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java index f73668642..9129ad728 100644 --- a/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java +++ b/src/test/java/net/imglib2/algorithm/morphology/distance/DistanceTransformTest.java @@ -35,6 +35,7 @@ package net.imglib2.algorithm.morphology.distance; import static org.junit.Assert.assertTrue; +import static org.junit.Assert.fail; import java.util.ArrayList; import java.util.Arrays; @@ -487,7 +488,15 @@ public void testLabelPropagation() final int numLabels = params.getIntPosition( 2 ); testLabelPropagationHelper( numDimensions, numLabels ); - testLabelPropagationHelperParallel( numDimensions, numLabels ); + try + { + testLabelPropagationHelperParallel( numDimensions, numLabels ); + } + catch ( Exception e ) + { + e.printStackTrace(); + fail(); + } } ); } @@ -507,7 +516,7 @@ private void testLabelPropagationHelper( int ndims, int nLabels ) final ArrayImg< LongType, LongArray > labels = ArrayImgs.longs( imgDims ); final Set< PointAndLabel > points = initializeLabels( rng, nLabels, labels ); - DistanceTransform.labelTransform( labels, 0 ); + DistanceTransform.voronoiDistanceTransform( labels, 0 ); validateLabelsSet( "serial", points, labels ); } @@ -518,14 +527,16 @@ private void testLabelPropagationHelper( int ndims, int nLabels ) * * @param ndims number of dimensions * @param nLabels number of labels + * @throws ExecutionException + * @throws InterruptedException */ - private void testLabelPropagationHelperParallel( int ndims, int nLabels ) + private void testLabelPropagationHelperParallel( int ndims, int nLabels ) throws InterruptedException, ExecutionException { final long[] imgDims = LongStream.iterate( dimensionSize, d -> d - 1 ).limit( ndims ).toArray(); final ArrayImg< LongType, LongArray > labels = ArrayImgs.longs( imgDims ); final Set< PointAndLabel > points = initializeLabels( rng, nLabels, labels ); - DistanceTransform.labelTransform( labels, 0, es, 3 * nThreads ); + DistanceTransform.voronoiDistanceTransform( labels, 0, es, 3 * nThreads ); validateLabelsSet( "parallel", points, labels ); } From aa69c03cd0f0afb2364c4d6f1b0eee537c777613 Mon Sep 17 00:00:00 2001 From: Stephan Saalfeld Date: Fri, 28 Feb 2025 14:41:49 -0500 Subject: [PATCH 6/7] doc: fix javadoc errors --- .../algorithm/morphology/distance/DistanceTransform.java | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index 8ecdbf083..b7bddf787 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -1466,15 +1466,12 @@ public static < L extends IntegerType< L >, T extends RealType< T > > void voron * the label type * @param labels * the label image - * @param backgroundLabel - * the background label * @param es * the ExecutorService * @param nTasks * the number of tasks in which to split the computation * @param weights * for distance computation per dimension - * @return the distance map * @throws ExecutionException * @throws InterruptedException */ @@ -1529,7 +1526,7 @@ public static < T extends RealType< T >, L extends IntegerType< L > > void trans * distances are given by the values in the distance argument. Results of * label propagation are stored in {@code labelsResult} *

- * This implementation uses the {@targetDistance} for temporary storage. + * This implementation uses the {@code targetDistance} for temporary storage. * * @param * input distance {@link RealType} @@ -1542,8 +1539,6 @@ public static < T extends RealType< T >, L extends IntegerType< L > > void trans * @param distance * Input distance function from which distance transform should * be computed. - * @param tmpDistance - * Storage for intermediate distance results. * @param targetDistance * Final result of distance transform. May be the same instance * as distance. From eb1268046afed99e51090ac5a0e9ac113d8fd9af Mon Sep 17 00:00:00 2001 From: Stephan Saalfeld Date: Fri, 28 Feb 2025 14:53:21 -0500 Subject: [PATCH 7/7] fix: use createVariable() instead of copy() where value is not required --- .../algorithm/morphology/distance/DistanceTransform.java | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java index b7bddf787..48d9182c4 100644 --- a/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java +++ b/src/main/java/net/imglib2/algorithm/morphology/distance/DistanceTransform.java @@ -738,7 +738,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final DISTANCE_TYPE distanceType, final double... weights ) { - final U maxVal = tmp.getType().copy(); + final U maxVal = tmp.getType().createVariable(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -792,7 +792,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final int nTasks, final double... weights ) throws InterruptedException, ExecutionException { - final U maxVal = tmp.getType().copy(); + final U maxVal = tmp.getType().createVariable(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -919,7 +919,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final RandomAccessibleInterval< V > target, final Distance d ) { - final U maxVal = tmp.getType().copy(); + final U maxVal = tmp.getType().createVariable(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() ); @@ -966,7 +966,7 @@ public static < B extends BooleanType< B >, U extends RealType< U >, V extends R final ExecutorService es, final int nTasks ) throws InterruptedException, ExecutionException { - final U maxVal = tmp.getType().copy(); + final U maxVal = tmp.getType().createVariable(); maxVal.setReal( maxVal.getMaxValue() ); final Converter< B, U > converter = new BinaryMaskToCost<>( maxVal ); final RandomAccessible< U > converted = Converters.convert( source, converter, maxVal.createVariable() );