From a1d0f967a313036560e6d0a2b3f3910f35311fb7 Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Fri, 27 Sep 2024 10:50:33 +0200 Subject: [PATCH 1/7] bugfix --- .../fusion/transformed/FusedRandomAccessibleInterval.java | 1 + 1 file changed, 1 insertion(+) diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/transformed/FusedRandomAccessibleInterval.java b/src/main/java/net/preibisch/mvrecon/process/fusion/transformed/FusedRandomAccessibleInterval.java index 3f5f8c4d..353dfcae 100644 --- a/src/main/java/net/preibisch/mvrecon/process/fusion/transformed/FusedRandomAccessibleInterval.java +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/transformed/FusedRandomAccessibleInterval.java @@ -53,6 +53,7 @@ public FusedRandomAccessibleInterval( this.n = interval.numDimensions(); this.interval = interval; this.images = images; + this.fusion = fusion; if ( weights == null || weights.size() == 0 ) { From 18ba4ef96816acd8adf1337259588663933ee1aa Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Wed, 25 Sep 2024 12:19:47 +0200 Subject: [PATCH 2/7] Simplify FusionTools.cacheRandomAccessibleInterval --- .../mvrecon/process/fusion/FusionTools.java | 56 ++++++------------- 1 file changed, 18 insertions(+), 38 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/FusionTools.java b/src/main/java/net/preibisch/mvrecon/process/fusion/FusionTools.java index e82caf82..def8bbfa 100644 --- a/src/main/java/net/preibisch/mvrecon/process/fusion/FusionTools.java +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/FusionTools.java @@ -47,7 +47,6 @@ import mpicbg.spim.data.generic.sequence.AbstractSequenceDescription; import mpicbg.spim.data.generic.sequence.BasicImgLoader; import mpicbg.spim.data.generic.sequence.BasicViewDescription; -import mpicbg.spim.data.generic.sequence.BasicViewSetup; import mpicbg.spim.data.registration.ViewRegistration; import mpicbg.spim.data.sequence.Angle; import mpicbg.spim.data.sequence.Channel; @@ -64,15 +63,19 @@ import net.imglib2.IterableInterval; import net.imglib2.RandomAccess; import net.imglib2.RandomAccessibleInterval; -import net.imglib2.cache.img.CellLoader; +import net.imglib2.cache.CacheLoader; +import net.imglib2.cache.img.RandomAccessibleCacheLoader; import net.imglib2.cache.img.ReadOnlyCachedCellImgFactory; import net.imglib2.cache.img.ReadOnlyCachedCellImgOptions; -import net.imglib2.cache.img.SingleCellArrayImg; import net.imglib2.cache.img.optional.CacheOptions.CacheType; import net.imglib2.converter.Converters; import net.imglib2.converter.RealTypeConverters; import net.imglib2.img.Img; import net.imglib2.img.ImgFactory; +import net.imglib2.img.basictypeaccess.AccessFlags; +import net.imglib2.img.basictypeaccess.array.ArrayDataAccess; +import net.imglib2.img.cell.Cell; +import net.imglib2.img.cell.CellGrid; import net.imglib2.img.imageplus.ImagePlusImgFactory; import net.imglib2.realtransform.AffineTransform3D; import net.imglib2.type.NativeType; @@ -783,49 +786,26 @@ public static < T extends NativeType< T > > RandomAccessibleInterval< T > cacheR return cacheRandomAccessibleInterval( input, -1, type, cellDim ); } - public static < T extends NativeType< T > > RandomAccessibleInterval< T > cacheRandomAccessibleInterval( + public static < T extends NativeType< T >, A extends ArrayDataAccess< A > > RandomAccessibleInterval< T > cacheRandomAccessibleInterval( final RandomAccessibleInterval< T > input, final long maxCacheSize, final T type, final int... cellDim ) { - final RandomAccessibleInterval< T > in; - - if ( Views.isZeroMin( input ) ) - in = input; - else - in = Views.zeroMin( input ); - - final ReadOnlyCachedCellImgOptions options; - - if ( maxCacheSize > 0 ) - options = new ReadOnlyCachedCellImgOptions().cellDimensions( cellDim ).maxCacheSize( maxCacheSize ); - else - options = new ReadOnlyCachedCellImgOptions().cellDimensions( cellDim ).cacheType( CacheType.SOFTREF ); - + final ReadOnlyCachedCellImgOptions options = ReadOnlyCachedCellImgOptions.options() + .cellDimensions( cellDim ) + .cacheType( maxCacheSize > 0 ? CacheType.BOUNDED : CacheType.SOFTREF ) + .maxCacheSize( maxCacheSize ); final ReadOnlyCachedCellImgFactory factory = new ReadOnlyCachedCellImgFactory( options ); - final CellLoader< T > loader = new CellLoader< T >() - { - @Override - public void load( final SingleCellArrayImg< T, ? > cell ) throws Exception - { - final Cursor< T > cursor = cell.localizingCursor(); - final RandomAccess< T > ra = in.randomAccess(); - - while( cursor.hasNext() ) - { - cursor.fwd(); - ra.setPosition( cursor ); - cursor.get().set( ra.get() ); - } - } - }; - - final long[] dim = new long[ in.numDimensions() ]; - in.dimensions( dim ); + final long[] dim = input.dimensionsAsLongArray(); + final CacheLoader< Long, Cell< A > > loader = RandomAccessibleCacheLoader.get( + new CellGrid( dim, cellDim ), + input.view().zeroMin(), + AccessFlags.setOf( AccessFlags.VOLATILE ) ); + final RandomAccessibleInterval copy = factory.createWithCacheLoader( dim, type, loader ); - return translateIfNecessary( input, factory.create( dim, type, loader ) ); + return translateIfNecessary( input, copy ); } public static < T extends Type< T > > RandomAccessibleInterval< T > copyImg( final RandomAccessibleInterval< T > input, final ImgFactory< T > factory, final T type, final ExecutorService service ) From 0a76c55eb52eda8ed52dbd7de22c72822ab4e772 Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Wed, 25 Sep 2024 17:03:12 +0200 Subject: [PATCH 3/7] ExportN5Api: clean up --- .../java/net/preibisch/mvrecon/process/export/ExportN5Api.java | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java b/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java index 0290fafa..dbd947de 100644 --- a/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java +++ b/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java @@ -63,9 +63,6 @@ import net.imglib2.RandomAccessibleInterval; import net.imglib2.type.NativeType; import net.imglib2.type.numeric.RealType; -import net.imglib2.type.numeric.integer.UnsignedByteType; -import net.imglib2.type.numeric.integer.UnsignedShortType; -import net.imglib2.type.numeric.real.FloatType; import net.imglib2.util.Intervals; import net.imglib2.util.Util; import net.imglib2.view.Views; From 35ffe4fc066f9ee983504b10c98124b35b15185a Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Fri, 13 Sep 2024 22:04:58 +0200 Subject: [PATCH 4/7] ExportN5Api: downsamplings == null fix --- .../java/net/preibisch/mvrecon/process/export/ExportN5Api.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java b/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java index dbd947de..a99244ec 100644 --- a/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java +++ b/src/main/java/net/preibisch/mvrecon/process/export/ExportN5Api.java @@ -228,7 +228,7 @@ else if ( storageType == StorageFormat.HDF5 ) bb.dimensionsAsLongArray(), compression, blocksize(), - this.downsampling, + this.downsampling == null ? new int[][] {{1,1,1}} : this.downsampling, viewId, path, xmlOut, From af4912fc970d115b013741c39d741b167ecd1608 Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Thu, 12 Sep 2024 17:50:04 +0200 Subject: [PATCH 5/7] BlkAffineFusion for AVG, AVG_BLEND, MAX, FIRST and 3D falls back to LazyAffineFusion for 2D and AVG_BLEND_CONTENT, AVG_CONTENT --- .../mvrecon/process/fusion/blk/Blending.java | 318 ++++++++++++++++++ .../process/fusion/blk/BlkAffineFusion.java | 239 +++++++++++++ .../mvrecon/process/fusion/blk/FirstWins.java | 193 +++++++++++ .../process/fusion/blk/LinearRange.java | 54 +++ .../mvrecon/process/fusion/blk/Masking.java | 216 ++++++++++++ .../process/fusion/blk/MaxIntensity.java | 120 +++++++ .../mvrecon/process/fusion/blk/Overlap.java | 186 ++++++++++ .../process/fusion/blk/WeightedAverage.java | 118 +++++++ 8 files changed, 1444 insertions(+) create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/Blending.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/BlkAffineFusion.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/FirstWins.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/LinearRange.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/Masking.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/MaxIntensity.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/Overlap.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/fusion/blk/WeightedAverage.java diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Blending.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Blending.java new file mode 100644 index 00000000..562a5fa2 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Blending.java @@ -0,0 +1,318 @@ +/*- + * #%L + * Spark-based parallel BigStitcher project. + * %% + * Copyright (C) 2021 - 2024 Developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.fusion.blk; + +import java.util.Arrays; + +import net.imglib2.Interval; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.real.FloatType; + +class Blending +{ + /** + * Conceptually,the given {@code interval} is filled with blending weights, then transformed with {@code transform}. + *

+ * Blending weights are {@code 0 <= w <= 1}. + *

+ * Weights are {@code w=0} for the outermost {@code border} pixels of {@code interval} (and outside of {@code interval}). + * Then weights transition from {@code 0<=w<=1} over {@code blending} pixels. + * Weights are {@code w=1} inside {@code border+blending} from the {@code interval} bounds. + * + * @param interval + * @param border + * @param blending + * @param transform + */ + public static BlockSupplier< FloatType > create( + final Interval interval, + final float[] border, + final float[] blending, + final AffineTransform3D transform ) + { + return new BlendingBlockSupplier( interval, border, blending, transform ); + } + + private static class BlendingBlockSupplier implements BlockSupplier< FloatType > + { + private final AffineTransform3D t; + + /** + * constant partial differential vector of t in X. + */ + private final double[] d0; + + private final int n = 3; + + /** + * min border distance. + * for {@code x + * Blending weights are {@code 0 <= w <= 1}. + *

+ * Weights are {@code w=0} for the outermost {@code border} pixels of {@code interval} (and outside of {@code interval}). + * Then weights transition from {@code 0<=w<=1} over {@code blending} pixels. + * Weights are {@code w=1} inside {@code border+blending} from the {@code interval} bounds. + * + * @param interval + * @param border + * @param blending + * @param transform + */ + BlendingBlockSupplier( + final Interval interval, + final float[] border, + final float[] blending, + final AffineTransform3D transform ) + { + // concatenate shift-to-interval-min to transform + t = new AffineTransform3D(); + t.translate( interval.minAsDoubleArray() ); + t.preConcatenate( transform ); + + d0 = t.inverse().d( 0 ).positionAsDoubleArray(); + + for ( int d = 0; d < n; ++d ) + { + final int dim = ( int ) interval.dimension( d ); + b0[ d ] = border[ d ]; + b1[ d ] = border[ d ] + blending[ d ]; + b2[ d ] = dim - 1 - border[ d ] - blending[ d ]; + b3[ d ] = dim - 1 - border[ d ]; + + if ( b1[ d ] > b2[ d ] ) // there is no "inside region" where w=1 + { + b1[ d ] = ( b1[ d ] + b2[ d ] ) / 2; + b2[ d ] = b1[ d ]; + } + + // TODO handle the case where border is so big that w=0 everywhere + } + + this.blending = blending.clone(); + } + + /** + * + * @param srcPos + * min coordinate of the block to copy + * @param dest + * {@code float[]} array to copy into. + * @param size + * the size of the block to copy + */ + @Override + public void copy( final long[] srcPos, final Object dest, final int[] size ) + { + final float[] weights = ( float[] ) dest; + final long x0 = srcPos[ 0 ]; + final long y0 = srcPos[ 1 ]; + final long z0 = srcPos[ 2 ]; + final int sx = size[ 0 ]; + final int sy = size[ 1 ]; + final int sz = size[ 2 ]; + final double[] p = { x0, 0, 0 }; + for ( int z = 0; z < sz; ++z ) + { + p[ 2 ] = z + z0; + for ( int y = 0; y < sy; ++y ) + { + p[ 1 ] = y + y0; + final int offset = ( z * sy + y ) * sx; + fill_range( weights, offset, sx, p ); + } + } + } + + @Override + public BlockSupplier< FloatType > threadSafe() + { + return this; + } + + @Override + public BlockSupplier< FloatType > independentCopy() + { + return this; + } + + @Override + public int numDimensions() + { + // TODO: do we need 2D ???? + return n; + } + + private static final FloatType type = new FloatType(); + + @Override + public FloatType getType() + { + return type; + } + + private static final float EPSILON = 0.0001f; + + private void fill_range( + float[] weights, + final int offset, + final int length, + double[] transformed_start_pos ) + { + final double[] pos = new double[ n ]; + t.applyInverse( pos, transformed_start_pos ); + Arrays.fill( weights, offset, offset + length, 1 ); + int from = 0; + int to = length; + for ( int d = 0; d < 3; ++d ) + { + final float l0 = ( float ) pos[ d ]; + final float dd = ( float ) d0[ d ]; + + final float blend; + final float b0d; + final float b1d; + final float b2d; + final float b3d; + if ( dd > EPSILON ) + { + blend = blending[ d ] / dd; + b0d = ( b0[ d ] - l0 ) / dd; + b1d = ( b1[ d ] - l0 ) / dd; + b2d = ( b2[ d ] - l0 ) / dd; + b3d = ( b3[ d ] - l0 ) / dd; + } + else if ( dd < -EPSILON ) + { + blend = blending[ d ] / -dd; + b0d = ( b3[ d ] - l0 ) / dd; + b1d = ( b2[ d ] - l0 ) / dd; + b2d = ( b1[ d ] - l0 ) / dd; + b3d = ( b0[ d ] - l0 ) / dd; + } + else + { + final float const_weight = computeWeight( l0, blending[ d ], b0[ d ], b1[ d ], b2[ d ], b3[ d ] ); + for ( int x = 0; x < length; ++x ) + weights[ offset + x ] *= const_weight; + continue; + } + + final int b3di = Math.max( from, Math.min( to, 1 + ( int ) Math.floor( b3d ) ) ); + final int b2di = Math.max( from, Math.min( b3di, 1 + ( int ) Math.floor( b2d ) ) ); + final int b1di = Math.max( from, Math.min( b2di, 1 + ( int ) Math.floor( b1d ) ) ); + final int b0di = Math.max( from, Math.min( b1di, 1 + ( int ) Math.floor( b0d ) ) ); + + for ( int x = from; x < b0di; ++x ) + weights[ offset + x ] = 0; + for ( int x = b0di; x < b1di; ++x ) + weights[ offset + x ] *= Lookup.get( ( x - b0d ) / blend ); + for ( int x = b2di; x < b3di; ++x ) + weights[ offset + x ] *= Lookup.get( ( b3d - x ) / blend ); + for ( int x = b3di; x < to; ++x ) + weights[ offset + x ] = 0; + + from = b0di; + to = b3di; + } + } + + private static float computeWeight( + final float l, + final float blending, + final float b0, + final float b1, + final float b2, + final float b3 ) + { + if ( l < b0 ) + return 0; + else if ( l < b1 ) + return Lookup.get( ( l - b0 ) / blending ); + else if ( l < b2 ) + return 1; + else if ( l < b3 ) + return Lookup.get( ( b3 - l ) / blending ); + else + return 0; + } + + /** + * Lookup table for blending weight function + * {@code fn(x) = (Math.cos((1 - x) * Math.PI) + 1) / 2} + */ + private static final class Lookup + { + private static final int n = 30; + + // static lookup table for the blending function + // size of the array is n + 2 + private static final float[] lookUp = createLookup( n ); + + private static float[] createLookup( final int n ) + { + final float[] lookup = new float[ n + 2 ]; + for ( int i = 0; i <= n; i++ ) + { + final double d = ( double ) i / n; + lookup[ i ] = ( float ) ( ( Math.cos( ( 1 - d ) * Math.PI ) + 1 ) / 2 ); + } + lookup[ n + 1 ] = lookup[ n ]; + return lookup; + } + + static float get( final float d ) + { + final int i = ( int ) ( d * n ); + final float s = ( d * n ) - i; + return lookUp[ i ] * (1.0f - s) + lookUp[ i + 1 ] * s; + } + } + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/BlkAffineFusion.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/BlkAffineFusion.java new file mode 100644 index 00000000..b97c7f2a --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/BlkAffineFusion.java @@ -0,0 +1,239 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import static net.imglib2.algorithm.blocks.transform.Transform.Interpolation.NEARESTNEIGHBOR; +import static net.imglib2.algorithm.blocks.transform.Transform.Interpolation.NLINEAR; + +import java.util.ArrayList; +import java.util.Collection; +import java.util.Collections; +import java.util.HashMap; +import java.util.List; +import java.util.Map; + +import mpicbg.models.AffineModel1D; +import mpicbg.spim.data.generic.sequence.BasicImgLoader; +import mpicbg.spim.data.generic.sequence.BasicViewDescription; +import mpicbg.spim.data.generic.sequence.BasicViewSetup; +import mpicbg.spim.data.sequence.ViewId; +import net.imglib2.Dimensions; +import net.imglib2.Interval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.blocks.BlockAlgoUtils; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.algorithm.blocks.ClampType; +import net.imglib2.algorithm.blocks.convert.Convert; +import net.imglib2.algorithm.blocks.transform.Transform; +import net.imglib2.algorithm.blocks.transform.Transform.Interpolation; +import net.imglib2.converter.Converter; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.integer.UnsignedByteType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.util.Util; +import net.imglib2.view.Views; +import net.preibisch.legacy.io.IOFunctions; +import net.preibisch.mvrecon.fiji.plugin.fusion.FusionGUI.FusionType; +import net.preibisch.mvrecon.process.downsampling.DownsampleTools; +import net.preibisch.mvrecon.process.fusion.FusionTools; +import net.preibisch.mvrecon.process.fusion.lazy.LazyAffineFusion; +import net.preibisch.mvrecon.process.fusion.lazy.LazyFusionTools; +import net.preibisch.mvrecon.process.interestpointregistration.pairwise.constellation.grouping.Group; + +public class BlkAffineFusion +{ + public static < T extends RealType< T > & NativeType< T > > RandomAccessibleInterval< T > init( + final Converter< FloatType, T > converter, + final BasicImgLoader imgloader, + final Collection< ? extends ViewId > viewIds, + final Map< ViewId, ? extends AffineTransform3D > viewRegistrations, + final Map< ViewId, ? extends BasicViewDescription< ? > > viewDescriptions, + final FusionType fusionType, + final int interpolationMethod, + final Map< ViewId, AffineModel1D > intensityAdjustments, + final Interval fusionInterval, + final T type, + final int[] blockSize ) + { + // go through the views and check if they are all 2-dimensional + final boolean is2d = viewIds.stream() + .map( viewDescriptions::get ) + .map( BasicViewDescription::getViewSetup ) + .filter( BasicViewSetup::hasSize ) + .allMatch( vs -> vs.getSize().dimension( 2 ) == 1 ); + + if ( !supports( is2d, fusionType, intensityAdjustments ) ) + { + IOFunctions.println( "BlkAffineFusion: Fusion method not supported (yet). Falling back to LazyAffineFusion." ); + return LazyAffineFusion.init( converter, imgloader, viewIds, viewRegistrations, viewDescriptions, fusionType, interpolationMethod, intensityAdjustments, fusionInterval, type, blockSize ); + } + + final HashMap< ViewId, Dimensions > viewDimensions = LazyFusionTools.assembleDimensions( viewIds, viewDescriptions ); + final Interpolation interpolation = ( interpolationMethod == 1 ) ? NLINEAR : NEARESTNEIGHBOR; + + // to be able to use the "lowest ViewId" wins strategy + final List< ? extends ViewId > sortedViewIds = new ArrayList<>( viewIds ); + Collections.sort( sortedViewIds ); + + // Which views to process (use un-altered bounding box and registrations). + // Final filtering happens per Cell. + // Here we just pre-filter everything outside the fusionInterval. + final Overlap overlap = new Overlap( + sortedViewIds, + viewRegistrations, + viewDimensions, + LazyFusionTools.defaultAffineExpansion, + is2d ? 2 : 3 ) + .filter( fusionInterval ) + .offset( fusionInterval.minAsLongArray() ); + + final List< BlockSupplier< FloatType > > images = new ArrayList<>( overlap.numViews() ); + final List< BlockSupplier< FloatType > > weights = new ArrayList<>( overlap.numViews() ); + final List< BlockSupplier< UnsignedByteType > > masks = new ArrayList<>( overlap.numViews() ); + + for ( final ViewId viewId : overlap.getViewIds() ) + { + final AffineTransform3D model = viewRegistrations.get( viewId ).copy(); + + // this modifies the model so it maps from a smaller image to the global coordinate space, + // which applies for the image itself as well as the weights since they also use the smaller + // input image as reference + final double[] usedDownsampleFactors = new double[ 3 ]; + RandomAccessibleInterval inputImg = DownsampleTools.openDownsampled( imgloader, viewId, model, usedDownsampleFactors ); + + final AffineTransform3D transform = concatenateBoundingBoxOffset( model, fusionInterval ); + + final BlockSupplier< FloatType > viewBlocks = transformedBlocks( + Cast.unchecked( inputImg ), + transform, interpolation ); + images.add( viewBlocks ); + + // instantiate blending if necessary + final float[] blending = Util.getArrayFromValue( FusionTools.defaultBlendingRange, 3 ); + final float[] border = Util.getArrayFromValue( FusionTools.defaultBlendingBorder, 3 ); + + // adjust both for z-scaling (anisotropy), downsampling, and registrations itself + FusionTools.adjustBlending( viewDimensions.get( viewId ), Group.pvid( viewId ), blending, border, model ); + + switch ( fusionType ) + { + case AVG: + weights.add( Masking.create( inputImg, border, transform ).andThen( Convert.convert( new FloatType() ) ) ); + break; + case AVG_BLEND: + weights.add( Blending.create( inputImg, border, blending, transform ) ); + break; + case MAX: + case FIRST: + masks.add( Masking.create( inputImg, border, transform ) ); + break; + case AVG_CONTENT: + case AVG_BLEND_CONTENT: + default: + // should never happen + throw new IllegalStateException(); + } + } + + final BlockSupplier< FloatType > floatBlocks; + switch ( fusionType ) + { + case AVG: + case AVG_BLEND: + floatBlocks = WeightedAverage.of( images, weights, overlap ); + break; + case MAX: + floatBlocks = MaxIntensity.of( images, masks, overlap ); + break; + case FIRST: + floatBlocks = FirstWins.of( images, masks, overlap ); + break; + case AVG_CONTENT: + case AVG_BLEND_CONTENT: + default: + // should never happen + throw new IllegalStateException(); + } + + final BlockSupplier< T > blocks = convertToOutputType( + floatBlocks, + converter, type ) + .tile( 32 ); + return BlockAlgoUtils.cellImg( blocks, fusionInterval.dimensionsAsLongArray(), blockSize ); + } + + private static < T extends NativeType< T > > BlockSupplier< T > convertToOutputType( + final BlockSupplier< FloatType > floatBlocks, + final Converter< FloatType, T > converter, + final T type ) + { + if ( converter == null ) + { + return floatBlocks.andThen( Convert.convert( type, ClampType.CLAMP ) ); + } + else if ( converter instanceof net.imglib2.display.LinearRange ) + { + final net.imglib2.display.LinearRange range = ( net.imglib2.display.LinearRange ) converter; + return floatBlocks + .andThen( LinearRange.linearRange( range.getMin(), range.getMax() ) ) + .andThen( Convert.convert( type, ClampType.CLAMP ) ); + } + else + { + return floatBlocks.andThen( Convert.convert( type, () -> converter ) ); + } + } + + + private static < T extends NativeType< T > > BlockSupplier< FloatType > transformedBlocks( + final RandomAccessibleInterval< T > inputImg, + final AffineTransform3D transform, + final Interpolation interpolation ) + { + return BlockSupplier.of( Views.extendBorder( ( inputImg ) ) ) + .andThen( Convert.convert( new FloatType() ) ) + .andThen( Transform.affine( transform, interpolation ) ); + } + + + private static AffineTransform3D concatenateBoundingBoxOffset( + final AffineTransform3D transformFromSource, + final Interval boundingBoxInTarget ) + { + final AffineTransform3D t = new AffineTransform3D(); + t.setTranslation( + -boundingBoxInTarget.min( 0 ), + -boundingBoxInTarget.min( 1 ), + -boundingBoxInTarget.min( 2 ) ); + t.concatenate( transformFromSource ); + return t; + } + + + private static < T extends RealType< T > & NativeType< T > > boolean supports( + final boolean is2d, + final FusionType fusionType, + final Map< ViewId, AffineModel1D > intensityAdjustments ) + { + if ( is2d ) + return false; // TODO + + switch ( fusionType ) + { + case AVG_BLEND: + case FIRST: + case MAX: + case AVG: + break; + case AVG_CONTENT: + case AVG_BLEND_CONTENT: + return false; // TODO + } + + if ( intensityAdjustments != null ) + return false; // TODO + + return true; + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/FirstWins.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/FirstWins.java new file mode 100644 index 00000000..7cd5470e --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/FirstWins.java @@ -0,0 +1,193 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import static net.imglib2.type.PrimitiveType.BYTE; +import static net.imglib2.type.PrimitiveType.FLOAT; +import static net.imglib2.util.Util.safeInt; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +import net.imglib2.algorithm.blocks.AbstractBlockSupplier; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.blocks.TempArray; +import net.imglib2.type.numeric.integer.UnsignedByteType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.util.Intervals; + +class FirstWins +{ + public static BlockSupplier< FloatType > of( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< UnsignedByteType > > masks, + final Overlap overlap ) + { + return new FirstWinsBlockSupplier( images, masks, overlap ); + } + + private static class FirstWinsBlockSupplier extends AbstractBlockSupplier< FloatType > + { + private final int numDimensions; + + private final List< BlockSupplier< FloatType > > images; + + private final List< BlockSupplier< UnsignedByteType > > masks; + + private final Overlap overlap; + + private final TempArray< byte[] > tempArrayM; + + private final TempArray< byte[] > tempArrayAccM; + + private final TempArray< float[] > tempArrayI; + + FirstWinsBlockSupplier( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< UnsignedByteType > > masks, + final Overlap overlap ) + { + this.numDimensions = images.get( 0 ).numDimensions(); + this.images = images; + this.masks = masks; + this.overlap = overlap; + tempArrayM = TempArray.forPrimitiveType( BYTE ); + tempArrayAccM = TempArray.forPrimitiveType( BYTE ); + tempArrayI = TempArray.forPrimitiveType( FLOAT ); + } + + private FirstWinsBlockSupplier( final FirstWinsBlockSupplier s ) + { + numDimensions = s.numDimensions; + images = new ArrayList<>( s.images.size() ); + masks = new ArrayList<>( s.masks.size() ); + s.images.forEach( i -> images.add( i.independentCopy() ) ); + s.masks.forEach( i -> masks.add( i.independentCopy() ) ); + overlap = s.overlap; + tempArrayM = TempArray.forPrimitiveType( BYTE ); + tempArrayAccM = TempArray.forPrimitiveType( BYTE ); + tempArrayI = TempArray.forPrimitiveType( FLOAT ); + } + + @Override + public void copy( final long[] srcPos, final Object dest, final int[] size ) + { + final int len = safeInt( Intervals.numElements( size ) ); + + final byte[] tmpM = tempArrayM.get( len ); + final byte[] accM = tempArrayAccM.get( len ); + final float[] tmpI = tempArrayI.get( len ); + final float[] fdest = Cast.unchecked( dest ); + + // algorithm: + // create byte[] accM and fill with 0 + // for ( i : overlapping ) { + // get masks(i) into tmpM byte[] + // iterate tmpM[x] and accM[x] + // if tmpM[x] == 1 and accM[x] == 0: + // set accM[x] = i+1 + // if any accM was set: + // get images(i) into tmpI float[] + // iterate tmpM[x], tmpI[x] and fdest[x] + // if accM[x] == i+1: + // set fdest[i] = tmpI[x] + // } + // finally, where accM[x] == 0, set fdest[i] = 0 + + Arrays.fill( accM, ( byte ) 0 ); + + // final long[] srcMax = new long[ srcPos.length ]; + // Arrays.setAll( srcMax, d -> srcPos[ d ] + size[ d ] - 1 ); + // final int[] overlapping = overlap.getOverlappingViewIndices( srcPos, srcMax ); + // for ( int i : overlapping ) + // { + // masks.get( i ).copy( srcPos, tmpM, size ); + // final byte flag = ( byte ) ( i + 1 ); + // boolean anySet = false; + // for ( int x = 0; x < len; ++x ) + // { + // if ( tmpM[ x ] == 1 && accM[ x ] == 0 ) + // { + // accM[ x ] = flag; + // anySet = true; + // } + // } + // if ( anySet ) + // { + // images.get( i ).copy( srcPos, tmpI, size ); + // for ( int x = 0; x < len; ++x ) + // { + // if ( accM[ x ] == flag ) + // fdest[ x ] = tmpI[ x ]; + // } + // } + // } + // for ( int x = 0; x < len; ++x ) + // { + // if ( accM[ x ] == 0 ) + // fdest[ x ] = 0; + // } + + final long[] srcMax = new long[ srcPos.length ]; + Arrays.setAll( srcMax, d -> srcPos[ d ] + size[ d ] - 1 ); + final int[] overlapping = overlap.getOverlappingViewIndices( srcPos, srcMax ); + int remaining = len; + for ( int i : overlapping ) + { + masks.get( i ).copy( srcPos, tmpM, size ); + final byte flag = ( byte ) ( i + 1 ); + final int remainingBefore = remaining; + for ( int x = 0; x < len; ++x ) + { + if ( tmpM[ x ] == 1 && accM[ x ] == 0 ) + { + accM[ x ] = flag; + --remaining; + } + } + if ( remaining != remainingBefore ) + { + images.get( i ).copy( srcPos, tmpI, size ); + for ( int x = 0; x < len; ++x ) + { + if ( accM[ x ] == flag ) + fdest[ x ] = tmpI[ x ]; + } + } + if ( remaining == 0 ) + return; + } + for ( int x = 0; x < len; ++x ) + { + if ( accM[ x ] == 0 ) + fdest[ x ] = 0; + } + + // TODO + // alternative to try: + // count down remaining accM[x]==0 from len + // when remaining==0 stop + // if after all masks remaining>0, then where accM[x] == 0, set fdest[i] = 0 + } + + @Override + public BlockSupplier< FloatType > independentCopy() + { + return new FirstWinsBlockSupplier( this ); + } + + @Override + public int numDimensions() + { + return numDimensions; + } + + private static final FloatType type = new FloatType(); + + @Override + public FloatType getType() + { + return type; + } + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/LinearRange.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/LinearRange.java new file mode 100644 index 00000000..3d4538f4 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/LinearRange.java @@ -0,0 +1,54 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import net.imglib2.algorithm.blocks.AbstractDimensionlessBlockProcessor; +import net.imglib2.algorithm.blocks.BlockProcessor; +import net.imglib2.algorithm.blocks.DefaultUnaryBlockOperator; +import net.imglib2.algorithm.blocks.UnaryBlockOperator; +import net.imglib2.type.PrimitiveType; +import net.imglib2.type.numeric.real.FloatType; + +class LinearRange +{ + public static UnaryBlockOperator< FloatType, FloatType > linearRange( final double min, final double max ) + { + final float offset = ( float ) (-min / ( max - min )); + final float scale = ( float ) (1.0 / ( max - min )); + final FloatType type = new FloatType(); + return new DefaultUnaryBlockOperator<>( type, type, 0, 0, new LinearRangeBlockProcessor( scale, offset ) ); + } + + private static class LinearRangeBlockProcessor extends AbstractDimensionlessBlockProcessor< float[], float[] > + { + private final float scale; + + private final float offset; + + LinearRangeBlockProcessor( final float scale, final float offset ) + { + super( PrimitiveType.FLOAT ); + this.scale = scale; + this.offset = offset; + } + + protected LinearRangeBlockProcessor( LinearRangeBlockProcessor proc ) + { + super( proc ); + this.scale = proc.scale; + this.offset = proc.offset; + } + + @Override + public BlockProcessor< float[], float[] > independentCopy() + { + return new LinearRangeBlockProcessor( this ); + } + + @Override + public void compute( final float[] src, final float[] dest ) + { + final int len = sourceLength(); + for ( int i = 0; i < len; i++ ) + dest[ i ] = src[ i ] * scale + offset; + } + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Masking.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Masking.java new file mode 100644 index 00000000..18dffc93 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Masking.java @@ -0,0 +1,216 @@ +/*- + * #%L + * Spark-based parallel BigStitcher project. + * %% + * Copyright (C) 2021 - 2024 Developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.fusion.blk; + +import java.util.Arrays; + +import net.imglib2.Interval; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.integer.UnsignedByteType; + +class Masking +{ + /** + * Conceptually,the given {@code interval} is filled with masking weights, then transformed with {@code transform}. + *

+ * Weights are {@code w=0} for the outermost {@code border} pixels of {@code interval}. + * Weights are {@code w=1} inside {@code border} from the {@code interval} bounds. + * + * @param interval + * @param border + * @param transform + */ + public static BlockSupplier< UnsignedByteType > create( + final Interval interval, + final float[] border, + final AffineTransform3D transform ) + { + return new MaskingBlockSupplier( interval, border, transform ); + } + + private static class MaskingBlockSupplier implements BlockSupplier< UnsignedByteType > + { + private final AffineTransform3D t; + + /** + * constant partial differential vector of t in X. + */ + private final double[] d0; + + private final int n = 3; + + /** + * min border distance. + * for {@code x + * Weights are {@code w=0} for the outermost {@code border} pixels of {@code interval}. + * Weights are {@code w=1} inside {@code border} from the {@code interval} bounds. + * + * @param interval + * @param border + * @param transform + */ + MaskingBlockSupplier( + final Interval interval, + final float[] border, + final AffineTransform3D transform ) + { + // concatenate shift-to-interval-min to transform + t = new AffineTransform3D(); + t.translate( interval.minAsDoubleArray() ); + t.preConcatenate( transform ); + + d0 = t.inverse().d( 0 ).positionAsDoubleArray(); + + for ( int d = 0; d < n; ++d ) + { + final int dim = ( int ) interval.dimension( d ); + b0[ d ] = border[ d ]; + b3[ d ] = dim - 1 - border[ d ]; + + // TODO handle the case where border is so big that w=0 everywhere + } + } + + /** + * + * @param srcPos + * min coordinate of the block to copy + * @param dest + * {@code float[]} array to copy into. + * @param size + * the size of the block to copy + */ + @Override + public void copy( final long[] srcPos, final Object dest, final int[] size ) + { + final byte[] weights = ( byte[] ) dest; + final long x0 = srcPos[ 0 ]; + final long y0 = srcPos[ 1 ]; + final long z0 = srcPos[ 2 ]; + final int sx = size[ 0 ]; + final int sy = size[ 1 ]; + final int sz = size[ 2 ]; + final double[] p = { x0, 0, 0 }; + for ( int z = 0; z < sz; ++z ) + { + p[ 2 ] = z + z0; + for ( int y = 0; y < sy; ++y ) + { + p[ 1 ] = y + y0; + final int offset = ( z * sy + y ) * sx; + fill_range( weights, offset, sx, p ); + } + } + } + + @Override + public BlockSupplier< UnsignedByteType > threadSafe() + { + return this; + } + + @Override + public BlockSupplier< UnsignedByteType > independentCopy() + { + return this; + } + + @Override + public int numDimensions() + { + // TODO: do we need 2D ???? + return n; + } + + private static final UnsignedByteType type = new UnsignedByteType(); + + @Override + public UnsignedByteType getType() + { + return type; + } + + private static final float EPSILON = 0.0001f; + + private void fill_range( + byte[] weights, + final int offset, + final int length, + double[] transformed_start_pos ) + { + final double[] pos = new double[ n ]; + t.applyInverse( pos, transformed_start_pos ); + int b0di = 0; + int b3di = length; + for ( int d = 0; d < 3; ++d ) + { + final float l0 = ( float ) pos[ d ]; + final float dd = ( float ) d0[ d ]; + + final float b0d; + final float b3d; + if ( dd > EPSILON ) + { + b0d = ( b0[ d ] - l0 ) / dd; + b3d = ( b3[ d ] - l0 ) / dd; + } + else if ( dd < -EPSILON ) + { + b0d = ( b3[ d ] - l0 ) / dd; + b3d = ( b0[ d ] - l0 ) / dd; + } + else + { + // this either sets everything to 0, or nothing. + if ( l0 < b0[ d ] || l0 >= b3[ d ] ) + { + Arrays.fill( weights, offset, offset + length, ( byte ) 0 ); + return; + } + continue; + } + + b3di = Math.max( b0di, Math.min( b3di, 1 + ( int ) Math.floor( b3d ) ) ); + b0di = Math.max( b0di, Math.min( b3di, 1 + ( int ) Math.floor( b0d ) ) ); + } + + Arrays.fill( weights, offset, offset + b0di, ( byte ) 0 ); + Arrays.fill( weights, offset + b0di, offset + b3di, ( byte ) 1 ); + Arrays.fill( weights, offset + b3di, offset + length, ( byte ) 0 ); + } + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/MaxIntensity.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/MaxIntensity.java new file mode 100644 index 00000000..d94b302f --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/MaxIntensity.java @@ -0,0 +1,120 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import static net.imglib2.type.PrimitiveType.BYTE; +import static net.imglib2.type.PrimitiveType.FLOAT; +import static net.imglib2.util.Util.safeInt; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +import net.imglib2.algorithm.blocks.AbstractBlockSupplier; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.blocks.TempArray; +import net.imglib2.type.numeric.integer.UnsignedByteType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.util.Intervals; + +class MaxIntensity +{ + public static BlockSupplier< FloatType > of( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< UnsignedByteType > > masks, + final Overlap overlap ) + { + return new MaxIntensityBlockSupplier( images, masks, overlap ); + } + + private static class MaxIntensityBlockSupplier extends AbstractBlockSupplier< FloatType > + { + private final int numDimensions; + + private final List< BlockSupplier< FloatType > > images; + + private final List< BlockSupplier< UnsignedByteType > > masks; + + private final Overlap overlap; + + private final TempArray< byte[] > tempArrayM; + + private final TempArray< float[] > tempArrayI; + + MaxIntensityBlockSupplier( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< UnsignedByteType > > masks, + final Overlap overlap ) + { + this.numDimensions = images.get( 0 ).numDimensions(); + this.images = images; + this.masks = masks; + this.overlap = overlap; + tempArrayM = TempArray.forPrimitiveType( BYTE ); + tempArrayI = TempArray.forPrimitiveType( FLOAT ); + } + + private MaxIntensityBlockSupplier( final MaxIntensityBlockSupplier s ) + { + numDimensions = s.numDimensions; + images = new ArrayList<>( s.images.size() ); + masks = new ArrayList<>( s.masks.size() ); + s.images.forEach( i -> images.add( i.independentCopy() ) ); + s.masks.forEach( i -> masks.add( i.independentCopy() ) ); + overlap = s.overlap; + tempArrayM = TempArray.forPrimitiveType( BYTE ); + tempArrayI = TempArray.forPrimitiveType( FLOAT ); + } + + @Override + public void copy( final long[] srcPos, final Object dest, final int[] size ) + { + final int len = safeInt( Intervals.numElements( size ) ); + + final byte[] tmpM = tempArrayM.get( len ); + final float[] tmpI = tempArrayI.get( len ); + final float[] fdest = Cast.unchecked( dest ); + + final long[] srcMax = new long[ srcPos.length ]; + Arrays.setAll( srcMax, d -> srcPos[ d ] + size[ d ] - 1 ); + final int[] overlapping = overlap.getOverlappingViewIndices( srcPos, srcMax ); + boolean first = true; + for ( int i : overlapping ) + { + masks.get( i ).copy( srcPos, tmpM, size ); + images.get( i ).copy( srcPos, tmpI, size ); + if ( first ) + { + first = false; + for ( int x = 0; x < len; ++x ) + fdest[ x ] = tmpM[ x ] == 1 ? tmpI[ x ] : 0; + } + else + { + for ( int x = 0; x < len; ++x ) + if ( tmpM[ x ] == 1 ) + fdest[ x ] = Math.max( fdest[ x ], tmpI[ x ] ); + } + } + } + + @Override + public BlockSupplier< FloatType > independentCopy() + { + return new MaxIntensityBlockSupplier( this ); + } + + @Override + public int numDimensions() + { + return numDimensions; + } + + private static final FloatType type = new FloatType(); + + @Override + public FloatType getType() + { + return type; + } + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Overlap.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Overlap.java new file mode 100644 index 00000000..a23d63da --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/Overlap.java @@ -0,0 +1,186 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; +import java.util.Map; + +import mpicbg.spim.data.sequence.ViewId; +import net.imglib2.Dimensions; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.RealInterval; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.util.Intervals; + +/** + * Compute and cache the expanded bounding boxes of all transformed views for per-block overlap determination. + * Sort order of {@code viewIds} is maintained for per-block queries. (for FIRST-WINS strategy) + */ +class Overlap +{ + private final List< ? extends ViewId > viewIds; + + private final long[] bb; + + private final int numDimensions; + + private final int numViews; + + public Overlap( + final List< ? extends ViewId > viewIds, + final Map< ? extends ViewId, ? extends AffineTransform3D > viewRegistrations, + final Map< ? extends ViewId, ? extends Dimensions > viewDimensions, + final int expandOverlap, + final int numDimensions ) + { + this.viewIds = viewIds; + this.numDimensions = numDimensions; + final List< Interval > bounds = new ArrayList<>(); + for ( final ViewId viewId : viewIds ) + { + // expand to be conservative ... + final AffineTransform3D t = viewRegistrations.get( viewId ); + final Dimensions dim = viewDimensions.get( viewId ); + final RealInterval ri = t.estimateBounds( new FinalInterval( dim ) ); + final Interval boundingBoxLocal = Intervals.largestContainedInterval( ri ); + bounds.add( Intervals.expand( boundingBoxLocal, expandOverlap ) ); + } + + numViews = bounds.size(); + bb = new long[ numViews * numDimensions * 2 ]; + for ( int i = 0; i < numViews; ++i ) + setBounds( i, bounds.get( i ) ); + } + + private void setBounds( int i, Interval interval ) + { + final int o_min = numDimensions * ( 2 * i ); + final int o_max = numDimensions * ( 2 * i + 1 ); + for ( int d = 0; d < numDimensions; d++ ) + { + bb[ o_min + d ] = interval.min( d ); + bb[ o_max + d ] = interval.max( d ); + } + } + + private long boundsMin( int i, int d ) + { + final int o_min = numDimensions * ( 2 * i ); + return bb[ o_min + d ]; + } + + private long boundsMax( int i, int d ) + { + final int o_max = numDimensions * ( 2 * i + 1 ); + return bb[ o_max + d ]; + } + + int[] getOverlappingViewIndices( final Interval interval ) + { + final long[] min = interval.minAsLongArray(); + final long[] max = interval.maxAsLongArray(); + return getOverlappingViewIndices( min, max ); + } + + int[] getOverlappingViewIndices( final long[] min, final long[] max ) + { + final int[] indices = new int[ numViews ]; + int j = 0; + for ( int i = 0; i < numViews; ++i ) + if ( isOverlapping( i, min, max ) ) + indices[ j++ ] = i; + return Arrays.copyOf( indices, j ); + } + + List< ViewId > getOverlappingViewIds( final Interval interval ) + { + final int[] indices = getOverlappingViewIndices( interval ); + final List< ViewId > views = new ArrayList<>(); + for ( final int i : indices ) + views.add( viewIds.get( i ) ); + return views; + } + + private boolean isOverlapping( final int i, final long[] min, final long[] max ) + { + for ( int d = 0; d < numDimensions; ++d ) + if ( min[ d ] > boundsMax( i, d ) || max[ d ] < boundsMin( i, d ) ) + return false; + return true; + } + + /** + * The ViewIds that are checked. + *

+ * Elements of {@code int[]} array returned by {@link #getOverlappingViewIds(Interval)} + * correspond to indices into this list. + */ + public List< ? extends ViewId > getViewIds() + { + return viewIds; + } + + public int numViews() + { + return numViews; + } + + private Overlap( + final List< ? extends ViewId > viewIds, + final long[] bb, + final int numDimensions ) + { + this.viewIds = viewIds; + this.bb = bb; + this.numDimensions = numDimensions; + numViews = viewIds.size(); + } + + /** + * Create a new {@code Overlap}, containing only those ViewIds of this + * {@code Overlap} which overlap the given {@code boundingBox}. + * + * @param boundingBox + * interval to check for overlap + * + * @return a new {@code Overlap}, containing only those ViewIds that overlap the given {@code boundingBox}. + */ + public Overlap filter( final Interval boundingBox ) + { + final int[] indices = getOverlappingViewIndices( boundingBox ); + final int filteredNumViews = indices.length; + final List< ViewId > filteredViews = new ArrayList<>( filteredNumViews ); + final long[] filteredBB = new long[ filteredNumViews * numDimensions * 2 ]; + for ( final int i : indices ) + { + final int o = numDimensions * ( 2 * i ); + final int fo = numDimensions * ( 2 * filteredViews.size() ); + for ( int k = 0; k < 2 * numDimensions; ++k ) + filteredBB[ fo + k ] = bb[ o + k ]; + filteredViews.add( viewIds.get( i ) ); + } + return new Overlap( filteredViews, filteredBB, numDimensions ); + } + + /** + * Return a new {@code Overlap}, with all bounding boxes shifted by the + * given {@code offset}. + * + * @param offset + * offset to apply + * + * @return a new {@code Overlap}, with shifted bounding boxes + */ + public Overlap offset( final long[] offset ) + { + final long[] offsetBB = new long[ bb.length ]; + for ( int i = 0; i < 2 * numViews; ++i ) + { + final int o = numDimensions * i; + for ( int d = 0; d < numDimensions; ++d ) + offsetBB[ o + d ] = bb[ o + d ] - offset[ d ]; + } + return new Overlap( viewIds, offsetBB, numDimensions ); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/fusion/blk/WeightedAverage.java b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/WeightedAverage.java new file mode 100644 index 00000000..a0719321 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/fusion/blk/WeightedAverage.java @@ -0,0 +1,118 @@ +package net.preibisch.mvrecon.process.fusion.blk; + +import static net.imglib2.type.PrimitiveType.FLOAT; +import static net.imglib2.util.Util.safeInt; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; + +import net.imglib2.algorithm.blocks.AbstractBlockSupplier; +import net.imglib2.algorithm.blocks.BlockSupplier; +import net.imglib2.blocks.TempArray; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.util.Intervals; + +class WeightedAverage +{ + public static BlockSupplier< FloatType > of( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< FloatType > > weights, + final Overlap overlap ) + { + return new WeightedAverageBlockSupplier( images, weights, overlap ); + } + + private static class WeightedAverageBlockSupplier extends AbstractBlockSupplier< FloatType > + { + private final int numDimensions; + + private final List< BlockSupplier< FloatType > > images; + + private final List< BlockSupplier< FloatType > > weights; + + private final Overlap overlap; + + private final TempArray< float[] >[] tempArrays; + + WeightedAverageBlockSupplier( + final List< BlockSupplier< FloatType > > images, + final List< BlockSupplier< FloatType > > weights, + final Overlap overlap ) + { + this.numDimensions = images.get( 0 ).numDimensions(); + this.images = images; + this.weights = weights; + this.overlap = overlap; + tempArrays = Cast.unchecked( new TempArray[ 4 ] ); + Arrays.setAll( tempArrays, i -> TempArray.forPrimitiveType( FLOAT ) ); + } + + private WeightedAverageBlockSupplier( final WeightedAverageBlockSupplier s ) + { + numDimensions = s.numDimensions; + images = new ArrayList<>( s.images.size() ); + weights = new ArrayList<>( s.weights.size() ); + s.images.forEach( i -> images.add( i.independentCopy() ) ); + s.weights.forEach( i -> weights.add( i.independentCopy() ) ); + overlap = s.overlap; + tempArrays = Cast.unchecked( new TempArray[ 4 ] ); + Arrays.setAll( tempArrays, i -> TempArray.forPrimitiveType( FLOAT ) ); + } + + @Override + public void copy( final long[] srcPos, final Object dest, final int[] size ) + { + final int len = safeInt( Intervals.numElements( size ) ); + final float[] tmpI = tempArrays[ 0 ].get( len ); + final float[] tmpW = tempArrays[ 1 ].get( len ); + final float[] sumI = tempArrays[ 2 ].get( len ); + final float[] sumW = tempArrays[ 3 ].get( len ); + + Arrays.fill( sumI, 0 ); + Arrays.fill( sumW, 0 ); + + final long[] srcMax = new long[ srcPos.length ]; + Arrays.setAll( srcMax, d -> srcPos[ d ] + size[ d ] - 1 ); + final int[] overlapping = overlap.getOverlappingViewIndices( srcPos, srcMax ); + for ( int i : overlapping ) + { + images.get( i ).copy( srcPos, tmpI, size ); + weights.get( i ).copy( srcPos, tmpW, size ); + for ( int x = 0; x < len; ++x ) + { + sumI[ x ] += tmpW[ x ] * tmpI[ x ]; + sumW[ x ] += tmpW[ x ]; + } + } + + final float[] fdest = Cast.unchecked( dest ); + for ( int x = 0; x < len; ++x ) + { + final float w = sumW[ x ]; + fdest[ x ] = ( w > 0 ) ? sumI[ x ] / w : 0; + } + } + + @Override + public BlockSupplier< FloatType > independentCopy() + { + return new WeightedAverageBlockSupplier( this ); + } + + @Override + public int numDimensions() + { + return numDimensions; + } + + private static final FloatType type = new FloatType(); + + @Override + public FloatType getType() + { + return type; + } + } +} From ac01d152ca5aa2b1e05afec3ed32b03d27445a9f Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Fri, 27 Sep 2024 21:47:58 +0200 Subject: [PATCH 6/7] Use BlkAffineFusion instead of LazyAffineFusion --- .../mvrecon/fiji/plugin/Image_Fusion.java | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/plugin/Image_Fusion.java b/src/main/java/net/preibisch/mvrecon/fiji/plugin/Image_Fusion.java index f7dfdd02..a0a781c8 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/plugin/Image_Fusion.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/plugin/Image_Fusion.java @@ -61,6 +61,7 @@ import net.preibisch.mvrecon.process.export.Calibrateable; import net.preibisch.mvrecon.process.export.ImgExport; import net.preibisch.mvrecon.process.fusion.FusionTools; +import net.preibisch.mvrecon.process.fusion.blk.BlkAffineFusion; import net.preibisch.mvrecon.process.fusion.lazy.LazyAffineFusion; import net.preibisch.mvrecon.process.fusion.lazy.LazyNonRigidFusion; import net.preibisch.mvrecon.process.fusion.transformed.TransformVirtual; @@ -245,7 +246,9 @@ else if ( fusion.getPixelType() == 1 ) } else { - lazy = LazyAffineFusion.init( + System.out.println( "Image_Fusion.fuse" ); +// lazy = LazyAffineFusion.init( + lazy = BlkAffineFusion.init( conv, spimData.getSequenceDescription().getImgLoader(), group.getViews(), @@ -257,20 +260,6 @@ else if ( fusion.getPixelType() == 1 ) fusion.getBoundingBox(), (RealType & NativeType)type, blocksize ); - - // TODO: replace with LazyAffineFusion and varying blocksizes depending on the task - /* - virtual = FusionTools.fuseVirtual( - spimData.getSequenceDescription().getImgLoader(), - registrations, - spimData.getSequenceDescription().getViewDescriptions(), - group.getViews(), - fusion.useBlending(), - fusion.useContentBased(), - fusion.getInterpolation(), - fusion.getBoundingBox(), - fusion.adjustIntensities() ? spimData.getIntensityAdjustments().getIntensityAdjustments() : null ); - */ } final String title = getTitle( fusion.getSplittingType(), group ); From ec107feb9599918f5e76220a16100b93118cdbb9 Mon Sep 17 00:00:00 2001 From: tpietzsch Date: Fri, 27 Sep 2024 22:55:43 +0200 Subject: [PATCH 7/7] Fix javadoc error --- .../mvrecon/fiji/datasetmanager/MultiViewDatasetDefinition.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/datasetmanager/MultiViewDatasetDefinition.java b/src/main/java/net/preibisch/mvrecon/fiji/datasetmanager/MultiViewDatasetDefinition.java index 031a1e00..6b7804aa 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/datasetmanager/MultiViewDatasetDefinition.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/datasetmanager/MultiViewDatasetDefinition.java @@ -49,7 +49,7 @@ public interface MultiViewDatasetDefinition * save it as an XML file. * * @param xmlFileName - the filename of the XML that was selected (no path!) - * @return - the saved {@link SpimData} object + * @return - the saved {@link SpimData2} object */ public SpimData2 createDataset( final String xmlFileName );