This command can import spots from label image data contained in a channel of the Big Data Viewer. The image data in that channel is assumed to represent an instance segmentation (i.e. a label image) that has been processed outside Mastodon. The existing labels in that channel will be used as spot names.
\n" + + "The BDV channel that contains the labels has to be chosen.
\n" + + "The ellipsoid scaling factor can be used to increase (>1) or decrease (<1) the size of the resulting ellipsoid. 1 is equivalent of ellipsoids drawn at 2.2σ.
\n" + + "\n" + + "\n"; + + @SuppressWarnings( "unused" ) + @Parameter + private ProjectModel projectModel; + + @Parameter( label = "Instance segmentation source", initializer = "initImgSourceChoices" ) + public String imgSourceChoice = ""; + + @SuppressWarnings( "all" ) + @Parameter( label = "Ellipsoid scaling factor", min = "0", description = "Changes the size of the resulting ellipsoid in all dimensions. 1 means that the ellipsoid is drawn at 2.2σ, which is the default." ) + private double scaleFactor = 1; + + @SuppressWarnings( "all" ) + @Parameter( label = "Link spots having the same labels in consecutive frames", description = "This option assumes that labels from the input images are unique for one tracklet." ) + boolean linkSpotsWithSameLabels = false; + + @SuppressWarnings( "unused" ) + private void initImgSourceChoices() + { + List< String > choices = LabelImageUtils.getSourceNames( projectModel.getSharedBdvData() ); + getInfo().getMutableInput( "imgSourceChoice", String.class ).setChoices( choices ); + } + + @Override + public void run() + { + Optional< SourceAndConverter< ? > > sourceAndConverter = projectModel.getSharedBdvData().getSources().stream() + .filter( source -> source.getSpimSource().getName().equals( imgSourceChoice ) ).findFirst(); + if ( !sourceAndConverter.isPresent() ) + return; + LabelImageUtils.importSpotsFromBdvChannel( projectModel, sourceAndConverter.get().getSpimSource(), scaleFactor, + linkSpotsWithSameLabels ); + } +} diff --git a/src/main/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusView.java b/src/main/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusView.java new file mode 100644 index 000000000..ae6ed9d25 --- /dev/null +++ b/src/main/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusView.java @@ -0,0 +1,88 @@ +package org.mastodon.mamut.io.importer.labelimage.ui; + +import net.imagej.ImgPlus; +import org.mastodon.mamut.ProjectModel; +import org.mastodon.mamut.io.importer.labelimage.LabelImageUtils; +import org.scijava.ItemIO; +import org.scijava.ItemVisibility; +import org.scijava.command.Command; +import org.scijava.command.ContextCommand; +import org.scijava.plugin.Parameter; +import org.scijava.plugin.Plugin; +import java.util.Arrays; + +@Plugin( type = Command.class, label = "Import spots from ImageJ image" ) +public class ImportSpotsFromImgPlusView< T > extends ContextCommand +{ + + private static final int WIDTH = 15; + + private static final String TEMPLATE = "\n" + + "\n" + + "This command can import spots from the active image in ImageJ that contains an instance segmentation that has been processed outside Mastodon. The label image is assumed to match the existing big data viewer image in all dimensions (x,y,z and t). The existing labels will be used as spot names.
\n" + + "The ellipsoid scaling factor can be used to increase (>1) or decrease (<1) the size of the resulting ellipsoid. 1 is equivalent of ellipsoids drawn at 2.2σ.
\n" + + "The active image in ImageJ is: %s.
\n"
+ + "
It has the these dimensions: x=%s, y=%s, z=%s, t=%s.
\n" + + "The big data viewer image has these dimensions: x=%s, y=%s, z=%s, t=%s.
\n" + + "%s
\n" + + "\n" + + "\n"; + + @Parameter( type = ItemIO.INPUT, validater = "validateImageData" ) + private ImgPlus< T > imgPlus; + + @SuppressWarnings( "unused" ) + @Parameter( visibility = ItemVisibility.MESSAGE, required = false, persist = false, initializer = "initMessage" ) + private String documentation; + + @SuppressWarnings( "unused" ) + @Parameter + private ProjectModel projectModel; + + @SuppressWarnings( "all" ) + @Parameter( label = "Ellipsoid scaling factor", min = "0", description = "Changes the size of the resulting ellipsoid in all dimensions. 1 means that the ellipsoid is drawn at 2.2σ, which is the default." ) + private double scaleFactor = 1; + + @SuppressWarnings( "all" ) + @Parameter( label = "Link spots having the same labels in consecutive frames", description = "This option assumes that labels from the input images are unique for one tracklet." ) + boolean linkSpotsWithSameLabels = false; + + @SuppressWarnings( "unused" ) + private void validateImageData() + { + if ( !LabelImageUtils.dimensionsMatch( projectModel.getSharedBdvData(), imgPlus ) ) + { + String imgPlusDimensions = Arrays.toString( LabelImageUtils.getImgPlusDimensions( imgPlus ) ); + String bdvDimensions = Arrays.toString( LabelImageUtils.getBdvDimensions( projectModel.getSharedBdvData() ) ); + cancel( String.format( + "The dimensions of the ImageJ image (%s) do not match the dimensions of the big data viewer image (%s)." + + "\nThus no spots can be imported." + + "\nPlease make sure the dimensions match.", + imgPlusDimensions, bdvDimensions ) ); + } + } + + @Override + public void run() + { + if ( isCanceled() ) + return; + LabelImageUtils.importSpotsFromImgPlus( projectModel, imgPlus, scaleFactor, linkSpotsWithSameLabels ); + } + + @SuppressWarnings( "unused" ) + private void initMessage() + { + if ( imgPlus == null ) + return; + long[] bdvDimensions = LabelImageUtils.getBdvDimensions( projectModel.getSharedBdvData() ); + long[] imgPlusDimensions = LabelImageUtils.getImgPlusDimensions( imgPlus ); + boolean dimensionsMatch = Arrays.equals( bdvDimensions, imgPlusDimensions ); + String color = dimensionsMatch ? "green" : "red"; + String doNot = dimensionsMatch ? "" : " do not"; + String dimensionMatch = "The dimensions" + doNot + " match."; + documentation = String.format( TEMPLATE, imgPlus.getName(), imgPlusDimensions[ 0 ], imgPlusDimensions[ 1 ], imgPlusDimensions[ 2 ], + imgPlusDimensions[ 3 ], bdvDimensions[ 0 ], bdvDimensions[ 1 ], bdvDimensions[ 2 ], bdvDimensions[ 3 ], dimensionMatch ); + } +} diff --git a/src/main/java/org/mastodon/mamut/util/LineageTreeUtils.java b/src/main/java/org/mastodon/mamut/util/LineageTreeUtils.java index 2eadbd5d3..5e4a60277 100644 --- a/src/main/java/org/mastodon/mamut/util/LineageTreeUtils.java +++ b/src/main/java/org/mastodon/mamut/util/LineageTreeUtils.java @@ -48,9 +48,13 @@ import org.mastodon.mamut.model.branch.BranchLink; import org.mastodon.mamut.model.branch.BranchSpot; import org.mastodon.mamut.model.branch.ModelBranchGraph; +import org.mastodon.spatial.SpatialIndex; import org.mastodon.util.TreeUtils; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; import javax.annotation.Nullable; +import java.lang.invoke.MethodHandles; import java.util.NoSuchElementException; import java.util.Set; import java.util.function.BooleanSupplier; @@ -59,6 +63,8 @@ public class LineageTreeUtils { + private static final Logger logger = LoggerFactory.getLogger( MethodHandles.lookup().lookupClass() ); + private LineageTreeUtils() { // prevent from instantiation @@ -283,4 +289,45 @@ public static < V extends Vertex< E >, E extends Edge< V > > RefSet< E > getAllE } return children; } + + /** + * This method adds edges to the graph of the model between spots that have the same label and are in consecutive time points. + *+ * Spot(label=0,X=1,Y=1,tp=0) Spot(label=1,X=0,Y=1,tp=0 ) Spot(label=2,X=2,Y=1,tp=0 ) + * + * Spot(label=0,X=1,Y=2,tp=1) Spot(label=1,X=0,Y=0,tp=1 ) Spot(label=2,X=2,Y=0,tp=1 ) + * + * Spot(label=0,X=1,Y=3,tp=2 ) Spot(label=1,X=0,Y=-1,tp=2 ) Spot(label=3,X=2,Y=-1,tp=2 ) + * + * Spot(label=0,X=1,Y=4,tp=3 ) Spot(label=0,X=0,Y=-2,tp=3 ) + *+ * This method will add edges between the following spots: + *
+ * Spot(label=0,X=1,Y=1,tp=0) Spot(label=1,X=0,Y=1,tp=0 ) Spot(label=2,X=2,Y=1,tp=0 ) + * │ │ │ + * Spot(label=0,X=1,Y=2,tp=1) Spot(label=1,X=0,Y=0,tp=1 ) Spot(label=2,X=2,Y=0,tp=1 ) + * │ │ + * Spot(label=0,X=1,Y=3,tp=2 ) Spot(label=1,X=0,Y=-1,tp=2 ) Spot(label=3,X=2,Y=-1,tp=2 ) + * ┌──────┴────────────────────────────┐ + * Spot(label=0,X=1,Y=4,tp=3 ) Spot(label=0,X=0,Y=-2,tp=3 ) + *+ * @param model the model to link spots in. + */ + public static void linkSpotsWithSameLabel( final Model model ) + { + Link edgeRef = model.getGraph().edgeRef(); + for ( int timepoint = TreeUtils.getMinTimepoint( model ); timepoint < TreeUtils.getMaxTimepoint( model ); timepoint++ ) + { + SpatialIndex< Spot > currentTimepoint = model.getSpatioTemporalIndex().getSpatialIndex( timepoint ); + SpatialIndex< Spot > nextTimepoint = model.getSpatioTemporalIndex().getSpatialIndex( timepoint + 1 ); + currentTimepoint.forEach( spotA -> nextTimepoint.forEach( spotB -> { + if ( spotA.getLabel().equals( spotB.getLabel() ) ) + model.getGraph().addEdge( spotA, spotB, edgeRef ).init(); + } ) ); + } + model.getGraph().releaseRef( edgeRef ); + logger.debug( "Added {} edge(s) to the graph.", model.getGraph().edges().size() ); + } } diff --git a/src/test/java/org/mastodon/mamut/feature/branch/exampleGraph/ExampleGraph6.java b/src/test/java/org/mastodon/mamut/feature/branch/exampleGraph/ExampleGraph6.java new file mode 100644 index 000000000..2540c39d0 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/feature/branch/exampleGraph/ExampleGraph6.java @@ -0,0 +1,59 @@ +package org.mastodon.mamut.feature.branch.exampleGraph; + +import org.mastodon.mamut.model.ModelGraph; +import org.mastodon.mamut.model.Spot; + +/** + * Represents a {@link AbstractExampleGraph} with the following {@link ModelGraph}: + * + *
+ * Spot( 0, X=1, Y=1, tp=0 ) Spot( 1, X=0, Y=1, tp=0 ) Spot( 2, X=2, Y=1, tp=0 ) + * + * Spot( 0, X=1, Y=2, tp=1 ) Spot( 1, X=0, Y=0, tp=1 ) Spot( 2, X=2, Y=0, tp=1 ) + * + * Spot( 0, X=1, Y=3, tp=2 ) Spot( 1, X=0, Y=-1, tp=2 ) Spot( 3, X=2, Y=-1, tp=2 ) + * + * Spot( 0, X=1, Y=4, tp=3 ) Spot( 0, X=0, Y=-2, tp=3 ) + *+ */ +public class ExampleGraph6 extends AbstractExampleGraph +{ + + public final Spot spot0; + + public final Spot spot1; + + public final Spot spot2; + + public final Spot spot3; + + public final Spot spot4; + + public final Spot spot5; + + public final Spot spot6; + + public final Spot spot7; + + public final Spot spot8; + + public final Spot spot9; + + public final Spot spot10; + + public ExampleGraph6() + { + spot0 = addNode( "0", 0, new double[] { 1d, 1d, 0d } ); + spot1 = addNode( "0", 1, new double[] { 1d, 2d, 0d } ); + spot2 = addNode( "0", 2, new double[] { 1d, 3d, 0d } ); + spot3 = addNode( "1", 0, new double[] { 0d, 1d, 0d } ); + spot4 = addNode( "1", 1, new double[] { 0d, 0d, 0d } ); + spot5 = addNode( "1", 2, new double[] { 0d, -1d, 0d } ); + spot6 = addNode( "2", 0, new double[] { 2d, 1d, 0d } ); + spot7 = addNode( "2", 1, new double[] { 2d, 0d, 0d } ); + spot8 = addNode( "3", 2, new double[] { 2d, -1d, 0d } ); + spot9 = addNode( "0", 3, new double[] { 1d, 4d, 0d } ); + spot10 = addNode( "0", 3, new double[] { 0d, -2d, 0d } ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/exporter/labelimage/ExportLabelImageControllerTest.java b/src/test/java/org/mastodon/mamut/io/exporter/labelimage/ExportLabelImageControllerTest.java index 93d5bbd1e..d9be57cc0 100644 --- a/src/test/java/org/mastodon/mamut/io/exporter/labelimage/ExportLabelImageControllerTest.java +++ b/src/test/java/org/mastodon/mamut/io/exporter/labelimage/ExportLabelImageControllerTest.java @@ -88,41 +88,43 @@ public void setUp() @Test public void testSaveLabelImageToFile() throws IOException { - AbstractSource< IntType > source = createRandomSource(); - Context context = new Context(); - TimePoint timePoint = new TimePoint( timepoint ); - List< TimePoint > timePoints = Collections.singletonList( timePoint ); - VoxelDimensions voxelDimensions = new DefaultVoxelDimensions( 3 ); - voxelDimensions.dimensions( new double[] { 1, 1, 1 } ); - ExportLabelImageController exportLabelImageController = - new ExportLabelImageController( model, timePoints, Cast.unchecked( source ), context, voxelDimensions ); - File outputSpot = getTempFile( "resultSpot" ); - File outputBranchSpot = getTempFile( "resultBranchSpot" ); - File outputTrack = getTempFile( "resultTrack" ); - exportLabelImageController.saveLabelImageToFile( LabelOptions.SPOT_ID, outputSpot, false, 1 ); - exportLabelImageController.saveLabelImageToFile( LabelOptions.BRANCH_SPOT_ID, outputBranchSpot, false, 1 ); - exportLabelImageController.saveLabelImageToFile( LabelOptions.TRACK_ID, outputTrack, false, 1 ); - - ImgOpener imgOpener = new ImgOpener( context ); - SCIFIOImgPlus< IntType > imgSpot = getIntTypeSCIFIOImgPlus( imgOpener, outputSpot ); - SCIFIOImgPlus< IntType > imgBranchSpot = getIntTypeSCIFIOImgPlus( imgOpener, outputBranchSpot ); - SCIFIOImgPlus< IntType > imgTrack = getIntTypeSCIFIOImgPlus( imgOpener, outputTrack ); - - // check that the spot id / branchSpot id / track id is used as value in the center of the spot - assertNotNull( imgSpot ); - assertEquals( 3, imgSpot.dimensionsAsLongArray().length ); - assertEquals( 100, imgSpot.dimension( 0 ) ); - assertEquals( 100, imgSpot.dimension( 1 ) ); - assertEquals( 100, imgSpot.dimension( 2 ) ); - assertEquals( spot.getInternalPoolIndex() + ExportLabelImageController.LABEL_ID_OFFSET, imgSpot.getAt( center ).get() ); - assertEquals( - branchSpot.getInternalPoolIndex() + ExportLabelImageController.LABEL_ID_OFFSET, imgBranchSpot.getAt( center ).get() ); - assertEquals( ExportLabelImageController.LABEL_ID_OFFSET, imgTrack.getAt( center ).get() ); - // check that there is no value set outside the ellipsoid of the spot - long[] corner = new long[] { 0, 0, 0 }; - assertEquals( 0, imgSpot.getAt( corner ).get() ); - assertEquals( 0, imgBranchSpot.getAt( corner ).get() ); - assertEquals( 0, imgTrack.getAt( corner ).get() ); + try (final Context context = new Context()) + { + AbstractSource< IntType > source = createRandomSource(); + TimePoint timePoint = new TimePoint( timepoint ); + List< TimePoint > timePoints = Collections.singletonList( timePoint ); + VoxelDimensions voxelDimensions = new DefaultVoxelDimensions( 3 ); + voxelDimensions.dimensions( new double[] { 1, 1, 1 } ); + ExportLabelImageController exportLabelImageController = + new ExportLabelImageController( model, timePoints, Cast.unchecked( source ), context, voxelDimensions ); + File outputSpot = getTempFile( "resultSpot" ); + File outputBranchSpot = getTempFile( "resultBranchSpot" ); + File outputTrack = getTempFile( "resultTrack" ); + exportLabelImageController.saveLabelImageToFile( LabelOptions.SPOT_ID, outputSpot, false, 1 ); + exportLabelImageController.saveLabelImageToFile( LabelOptions.BRANCH_SPOT_ID, outputBranchSpot, false, 1 ); + exportLabelImageController.saveLabelImageToFile( LabelOptions.TRACK_ID, outputTrack, false, 1 ); + + ImgOpener imgOpener = new ImgOpener( context ); + SCIFIOImgPlus< IntType > imgSpot = getIntTypeSCIFIOImgPlus( imgOpener, outputSpot ); + SCIFIOImgPlus< IntType > imgBranchSpot = getIntTypeSCIFIOImgPlus( imgOpener, outputBranchSpot ); + SCIFIOImgPlus< IntType > imgTrack = getIntTypeSCIFIOImgPlus( imgOpener, outputTrack ); + + // check that the spot id / branchSpot id / track id is used as value in the center of the spot + assertNotNull( imgSpot ); + assertEquals( 3, imgSpot.dimensionsAsLongArray().length ); + assertEquals( 100, imgSpot.dimension( 0 ) ); + assertEquals( 100, imgSpot.dimension( 1 ) ); + assertEquals( 100, imgSpot.dimension( 2 ) ); + assertEquals( spot.getInternalPoolIndex() + ExportLabelImageController.LABEL_ID_OFFSET, imgSpot.getAt( center ).get() ); + assertEquals( + branchSpot.getInternalPoolIndex() + ExportLabelImageController.LABEL_ID_OFFSET, imgBranchSpot.getAt( center ).get() ); + assertEquals( ExportLabelImageController.LABEL_ID_OFFSET, imgTrack.getAt( center ).get() ); + // check that there is no value set outside the ellipsoid of the spot + long[] corner = new long[] { 0, 0, 0 }; + assertEquals( 0, imgSpot.getAt( corner ).get() ); + assertEquals( 0, imgBranchSpot.getAt( corner ).get() ); + assertEquals( 0, imgTrack.getAt( corner ).get() ); + } } private static SCIFIOImgPlus< IntType > getIntTypeSCIFIOImgPlus( ImgOpener imgOpener, File outputSpot ) @@ -142,15 +144,18 @@ private static File getTempFile( final String prefix ) throws IOException @Test public void testExceptions() throws IOException { - ExportLabelImageController controller = - new ExportLabelImageController( model, Collections.emptyList(), null, new Context(), null ); - File file = File.createTempFile( "foo", "foo" ); - file.deleteOnExit(); - assertThrows( - IllegalArgumentException.class, - () -> controller.saveLabelImageToFile( LabelOptions.SPOT_ID, null, false, 1 ) - ); - assertThrows( IllegalArgumentException.class, () -> controller.saveLabelImageToFile( null, file, false, 1 ) ); + try (final Context context = new Context()) + { + ExportLabelImageController controller = + new ExportLabelImageController( model, Collections.emptyList(), null, context, null ); + File file = File.createTempFile( "foo", "foo" ); + file.deleteOnExit(); + assertThrows( + IllegalArgumentException.class, + () -> controller.saveLabelImageToFile( LabelOptions.SPOT_ID, null, false, 1 ) + ); + assertThrows( IllegalArgumentException.class, () -> controller.saveLabelImageToFile( null, file, false, 1 ) ); + } } private static AbstractSource< IntType > createRandomSource() diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/LabelImageUtilsTest.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/LabelImageUtilsTest.java new file mode 100644 index 000000000..372d7576a --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/LabelImageUtilsTest.java @@ -0,0 +1,305 @@ +package org.mastodon.mamut.io.importer.labelimage; + +import bdv.spimdata.SequenceDescriptionMinimal; +import bdv.util.AbstractSource; +import bdv.util.RandomAccessibleIntervalSource; +import mpicbg.spim.data.generic.sequence.AbstractSequenceDescription; +import mpicbg.spim.data.generic.sequence.BasicViewSetup; +import mpicbg.spim.data.sequence.FinalVoxelDimensions; +import mpicbg.spim.data.sequence.TimePoint; +import mpicbg.spim.data.sequence.TimePoints; +import mpicbg.spim.data.sequence.VoxelDimensions; +import net.imagej.ImgPlus; +import net.imagej.axis.Axes; +import net.imagej.axis.CalibratedAxis; +import net.imagej.axis.DefaultLinearAxis; +import net.imagej.patcher.LegacyInjector; +import net.imglib2.FinalDimensions; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.loops.LoopBuilder; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.view.Views; +import org.junit.Before; +import org.junit.Test; +import org.mastodon.mamut.ProjectModel; +import org.mastodon.mamut.feature.EllipsoidIterable; +import org.mastodon.mamut.io.importer.labelimage.util.DemoUtils; +import org.mastodon.mamut.model.Model; +import org.mastodon.mamut.model.Spot; +import org.mastodon.views.bdv.overlay.util.JamaEigenvalueDecomposition; +import org.scijava.Context; +import org.slf4j.Logger; +import org.slf4j.LoggerFactory; + +import java.lang.invoke.MethodHandles; +import java.util.Arrays; +import java.util.Collections; +import java.util.Iterator; +import java.util.List; +import java.util.Map; +import java.util.concurrent.atomic.AtomicInteger; +import java.util.function.IntFunction; + +import static org.junit.Assert.assertArrayEquals; +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertFalse; +import static org.junit.Assert.assertNotNull; +import static org.junit.Assert.assertThrows; +import static org.junit.Assert.assertTrue; + +public class LabelImageUtilsTest +{ + private static final Logger logger = LoggerFactory.getLogger( MethodHandles.lookup().lookupClass() ); + + private Model model; + + private AbstractSequenceDescription< ?, ?, ? > sequenceDescription; + + @Before + public void setUp() + { + model = new Model(); + TimePoints timePoints = new TimePoints( Collections.singletonList( new TimePoint( 0 ) ) ); + VoxelDimensions voxelDimensions = new FinalVoxelDimensions( "um", 1, 1, 1 ); + Map< Integer, ? extends BasicViewSetup > setups = + Collections.singletonMap( 0, new BasicViewSetup( 0, "setup 0", new FinalDimensions( 10, 10, 10 ), voxelDimensions ) ); + sequenceDescription = new SequenceDescriptionMinimal( timePoints, setups, null, null ); + } + + @Test + public void testScaleExceptions() + { + VoxelDimensions voxelDimensions = new FinalVoxelDimensions( "um", 1, 1, 1 ); + assertThrows( IllegalArgumentException.class, () -> LabelImageUtils.scale( new double[ 2 ][ 2 ], 1, voxelDimensions ) ); + assertThrows( IllegalArgumentException.class, () -> LabelImageUtils.scale( new double[ 3 ][ 2 ], 1, voxelDimensions ) ); + } + + @Test + public void testCreateSpotFromLabelImageEmpty() + { + RandomAccessibleIntervalSource< FloatType > img = + new RandomAccessibleIntervalSource<>( createImageCubeCorners( 0 ), new FloatType(), new AffineTransform3D(), + "Segmentation" ); + + IntFunction< RandomAccessibleInterval< RealType< ? > > > imgProvider = frameId -> Cast.unchecked( img.getSource( frameId, 0 ) ); + LabelImageUtils.createSpotsFromLabelImage( imgProvider, model, 1, false, sequenceDescription, null ); + assertEquals( 0, model.getGraph().vertices().size() ); + } + + @Test + public void testCreateSpotFromNonLabelImage() + { + AbstractSource< FloatType > img = createNonLabelImage(); + + IntFunction< RandomAccessibleInterval< RealType< ? > > > imgProvider = frameId -> Cast.unchecked( img.getSource( frameId, 0 ) ); + LabelImageUtils.createSpotsFromLabelImage( imgProvider, model, 1, false, sequenceDescription, null ); + assertEquals( 0, model.getGraph().vertices().size() ); + } + + @Test + public void testCreateSpotFromWrongVoxelDimensions() + { + + RandomAccessibleIntervalSource< FloatType > img = + new RandomAccessibleIntervalSource<>( createImageCubeCorners( 1 ), new FloatType(), new AffineTransform3D(), + "Segmentation" ); + + VoxelDimensions wrongDimensions = new FinalVoxelDimensions( "um", 1, 1 ); + TimePoints timePoints = new TimePoints( Collections.singletonList( new TimePoint( 0 ) ) ); + Map< Integer, ? extends BasicViewSetup > setups = + Collections.singletonMap( 0, new BasicViewSetup( 0, "setup 0", new FinalDimensions( 10, 10, 10 ), wrongDimensions ) ); + AbstractSequenceDescription< ?, ?, ? > faultySequenceDescription = new SequenceDescriptionMinimal( timePoints, setups, null, null ); + IntFunction< RandomAccessibleInterval< RealType< ? > > > imgProvider = frameId -> Cast.unchecked( img.getSource( frameId, 0 ) ); + assertThrows( IllegalArgumentException.class, + () -> LabelImageUtils.createSpotsFromLabelImage( imgProvider, model, 1, false, faultySequenceDescription, null ) ); + + } + + @Test + public void testImportSpotsFromBdvChannel() + { + LegacyInjector.preinit(); + try (Context context = new Context()) + { + int pixelValue = 1; + Img< FloatType > img = createImageCubeCorners( pixelValue ); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( img, model, context ); + LabelImageUtils.importSpotsFromBdvChannel( projectModel, projectModel.getSharedBdvData().getSources().get( 0 ).getSpimSource(), + 1, false ); + + Iterator< Spot > iter = model.getGraph().vertices().iterator(); + Spot spot = iter.next(); + double[][] covarianceMatrix = new double[ 3 ][ 3 ]; + spot.getCovariance( covarianceMatrix ); + final JamaEigenvalueDecomposition eigenvalueDecomposition = new JamaEigenvalueDecomposition( 3 ); + eigenvalueDecomposition.decomposeSymmetric( covarianceMatrix ); + final double[] eigenValues = eigenvalueDecomposition.getRealEigenvalues(); + double axisA = Math.sqrt( eigenValues[ 0 ] ); + double axisB = Math.sqrt( eigenValues[ 1 ] ); + double axisC = Math.sqrt( eigenValues[ 2 ] ); + + assertNotNull( spot ); + assertEquals( 0, spot.getTimepoint() ); + assertEquals( 2, spot.getDoublePosition( 0 ), 0.01 ); + assertEquals( 2, spot.getDoublePosition( 1 ), 0.01 ); + assertEquals( 2, spot.getDoublePosition( 2 ), 0.01 ); + assertEquals( 0, spot.getInternalPoolIndex() ); + assertEquals( String.valueOf( pixelValue ), spot.getLabel() ); + assertEquals( 2.2, axisA, 0.2d ); + assertEquals( 2.2, axisB, 0.2d ); + assertEquals( 2.2, axisC, 0.2d ); + assertEquals( 5d, spot.getBoundingSphereRadiusSquared(), 1d ); + assertFalse( iter.hasNext() ); + } + } + + @Test + public void testImportSpotsFromImgPlus() + { + LegacyInjector.preinit(); + try (Context context = new Context()) + { + double[] center = { 18, 21, 22 }; + double[][] givenCovariance = { + { 33, 14, 0 }, + { 14, 32, 0 }, + { 0, 0, 95 } + }; + Spot spot = model.getGraph().addVertex().init( 0, center, givenCovariance ); + int pixelValue = 1; + Img< FloatType > image = createImageFromSpot( spot, pixelValue ); + ImgPlus< FloatType > imgPlus = createImgPlus( image, new FinalVoxelDimensions( "um", 1, 1, 1 ) ); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( image, model, context ); + LabelImageUtils.importSpotsFromImgPlus( projectModel, imgPlus, 1, false ); + + Iterator< Spot > iterator = model.getGraph().vertices().iterator(); + iterator.next(); + Spot createdSpot = iterator.next(); + double[] mean = createdSpot.positionAsDoubleArray(); + double[][] computedCovariance = new double[ 3 ][ 3 ]; + createdSpot.getCovariance( computedCovariance ); + + logger.debug( "Given center: {}", Arrays.toString( center ) ); + logger.debug( "Computed mean: {}", Arrays.toString( mean ) ); + logger.debug( "Given covariance: {}", Arrays.deepToString( givenCovariance ) ); + logger.debug( "Computed covariance: {}", Arrays.deepToString( computedCovariance ) ); + + assertArrayEquals( center, mean, 0.01d ); + assertArrayEquals( givenCovariance[ 0 ], computedCovariance[ 0 ], 10d ); + assertArrayEquals( givenCovariance[ 1 ], computedCovariance[ 1 ], 10d ); + assertArrayEquals( givenCovariance[ 2 ], computedCovariance[ 2 ], 10d ); + assertEquals( String.valueOf( pixelValue ), createdSpot.getLabel() ); + } + } + + @Test + public void testImportSpotsFromImgPlusAndLinkSameLabels() + { + LegacyInjector.preinit(); + try (Context context = new Context()) + { + Img< FloatType > twoFramesImage = DemoUtils.generateExampleImage(); + ImgPlus< FloatType > imgPlus = createImgPlus( twoFramesImage, new FinalVoxelDimensions( "um", 1, 1, 1 ) ); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( twoFramesImage, model, context ); + LabelImageUtils.importSpotsFromImgPlus( projectModel, imgPlus, 1, true ); + + assertEquals( 2, model.getGraph().vertices().size() ); + assertEquals( 1, model.getGraph().edges().size() ); + assertEquals( 1, model.getSpatioTemporalIndex().getSpatialIndex( 0 ).size() ); + assertEquals( 1, model.getSpatioTemporalIndex().getSpatialIndex( 1 ).size() ); + } + } + + @Test + public void testDimensionsMatch() + { + LegacyInjector.preinit(); + try (final Context context = new Context()) + { + Img< FloatType > image = ArrayImgs.floats( 10, 10, 10, 2 ); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( image, model, context ); + ImgPlus< FloatType > imgPlus = createImgPlus( image, new FinalVoxelDimensions( "um", 1, 1, 1 ) ); + assertTrue( LabelImageUtils.dimensionsMatch( projectModel.getSharedBdvData(), imgPlus ) ); + } + } + + @Test + public void testGetImgPlusDimensions() + { + Img< FloatType > image = ArrayImgs.floats( 100, 100, 100, 2 ); + ImgPlus< FloatType > imgPlus = createImgPlus( image, new FinalVoxelDimensions( "um", 1, 1, 1 ) ); + long[] dimensions = LabelImageUtils.getImgPlusDimensions( imgPlus ); + assertArrayEquals( new long[] { 100, 100, 100, 2 }, dimensions ); + } + + @Test + public void testGetSourceNames() + { + LegacyInjector.preinit(); + try (final Context context = new Context()) + { + Img< FloatType > image = ArrayImgs.floats( 100, 100, 100, 2 ); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( image, model, context ); + List< String > sourceNames = LabelImageUtils.getSourceNames( projectModel.getSharedBdvData() ); + assertEquals( 1, sourceNames.size() ); + assertEquals( "image channel 1", sourceNames.get( 0 ) ); + } + } + + private ImgPlus< FloatType > createImgPlus( final Img< FloatType > img, final VoxelDimensions voxelDimensions ) + { + final CalibratedAxis[] axes = { new DefaultLinearAxis( Axes.X, voxelDimensions.dimension( 0 ) ), + new DefaultLinearAxis( Axes.Y, voxelDimensions.dimension( 1 ) ), + new DefaultLinearAxis( Axes.Z, voxelDimensions.dimension( 2 ) ), new DefaultLinearAxis( Axes.TIME ) }; + return new ImgPlus<>( img, "Result", axes ); + } + + private static Img< FloatType > createImageCubeCorners( int pixelValue ) + { + Img< FloatType > img = new ArrayImgFactory<>( new FloatType() ).create( 4, 4, 4 ); + RandomAccess< FloatType > ra = img.randomAccess(); + // 8 corners of a cube + ra.setPositionAndGet( 1, 1, 1 ).set( pixelValue ); + ra.setPositionAndGet( 1, 3, 1 ).set( pixelValue ); + ra.setPositionAndGet( 3, 1, 1 ).set( pixelValue ); + ra.setPositionAndGet( 3, 3, 1 ).set( pixelValue ); + ra.setPositionAndGet( 1, 1, 3 ).set( pixelValue ); + ra.setPositionAndGet( 1, 3, 3 ).set( pixelValue ); + ra.setPositionAndGet( 3, 1, 3 ).set( pixelValue ); + ra.setPositionAndGet( 3, 3, 3 ).set( pixelValue ); + + return img; + } + + private static AbstractSource< FloatType > createNonLabelImage() + { + Img< FloatType > img = new ArrayImgFactory<>( new FloatType() ).create( 25, 25, 25 ); + AtomicInteger value = new AtomicInteger( 0 ); + LoopBuilder.setImages( img ).forEachPixel( floatType -> floatType.set( value.incrementAndGet() ) ); + + return new RandomAccessibleIntervalSource<>( img, new FloatType(), new AffineTransform3D(), + "Segmentation" ); + } + + private static Img< FloatType > createImageFromSpot( final Spot spot, int pixelValue ) + { + long[] dimensions = { 40, 40, 40, 1 }; + Img< FloatType > image = ArrayImgs.floats( dimensions ); + AffineTransform3D transform = new AffineTransform3D(); + AbstractSource< FloatType > frame = + new RandomAccessibleIntervalSource<>( Views.hyperSlice( image, 3, 0 ), new FloatType(), transform, "Ellipsoids" ); + + final EllipsoidIterable< FloatType > ellipsoidIterable = new EllipsoidIterable<>( frame ); + ellipsoidIterable.reset( spot ); + ellipsoidIterable.forEach( pixel -> pixel.set( pixelValue ) ); + return image; + } +} + diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceMatrixTest.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceMatrixTest.java new file mode 100644 index 000000000..a5dc570e3 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceMatrixTest.java @@ -0,0 +1,40 @@ +package org.mastodon.mamut.io.importer.labelimage.math; + +import org.junit.Test; + +import java.util.Arrays; + +import static org.junit.Assert.assertArrayEquals; +import static org.junit.Assert.assertThrows; + +public class CovarianceMatrixTest +{ + @Test + public void testGet() + { + int[][] dataInt = { { 1, 2 }, { 2, 3 }, { 3, 4 }, { 4, 5 }, { 5, 6 } }; + CovarianceMatrix matrix = new CovarianceMatrix( 2 ); + for ( int[] values : dataInt ) + matrix.addValues( values ); + double[][] actual = matrix.get(); + + assertArrayEquals( new double[] { 3d, 4d }, matrix.getMeans(), 0.0001d ); + assertArrayEquals( new double[] { 2.5d, 2.5d }, actual[ 0 ], 0.0001d ); + assertArrayEquals( new double[] { 2.5d, 2.5d }, actual[ 1 ], 0.0001d ); + } + + @Test + public void testException1() + { + CovarianceMatrix covarianceMatrix = new CovarianceMatrix( 2 ); + covarianceMatrix.addValues( new int[] { 1, 1 } ); + assertThrows( IllegalArgumentException.class, covarianceMatrix::get ); + } + + @Test + public void testException2() + { + CovarianceMatrix covarianceMatrix = new CovarianceMatrix( 2 ); + assertThrows( IllegalArgumentException.class, () -> covarianceMatrix.addValues( new int[] { 1 } ) ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceTest.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceTest.java new file mode 100644 index 000000000..dd6e98d37 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/CovarianceTest.java @@ -0,0 +1,33 @@ +package org.mastodon.mamut.io.importer.labelimage.math; + +import org.junit.Test; + +import static org.junit.Assert.assertEquals; +import static org.junit.Assert.assertThrows; + +public class CovarianceTest +{ + @Test + public void testGet() + { + double[] x = { 1, 2, 3, 4, 5 }; + double[] y = { 2, 3, 4, 5, 6 }; + org.mastodon.mamut.io.importer.labelimage.math.Covariance covariance = + new org.mastodon.mamut.io.importer.labelimage.math.Covariance(); + for ( int i = 0; i < x.length; i++ ) + covariance.addValues( x[ i ], y[ i ] ); + double actual = covariance.get(); + assertEquals( 2.5d, actual, 0.0001d ); + assertEquals( 3d, covariance.getMeanX(), 0.0001d ); + assertEquals( 4d, covariance.getMeanY(), 0.0001d ); + } + + @Test + public void testException() + { + org.mastodon.mamut.io.importer.labelimage.math.Covariance covariance = + new org.mastodon.mamut.io.importer.labelimage.math.Covariance(); + covariance.addValues( 1, 1 ); + assertThrows( IllegalArgumentException.class, covariance::get ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/MeansVectorTest.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/MeansVectorTest.java new file mode 100644 index 000000000..f74e9b186 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/math/MeansVectorTest.java @@ -0,0 +1,15 @@ +package org.mastodon.mamut.io.importer.labelimage.math; + +import org.junit.Test; + +import static org.junit.Assert.assertThrows; + +public class MeansVectorTest +{ + @Test + public void testAddValuesException() + { + MeansVector meansVector = new MeansVector( 3 ); + assertThrows( IllegalArgumentException.class, () -> meansVector.addValues( new int[] { 1, 2 } ) ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromBdvChannelViewDemo.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromBdvChannelViewDemo.java new file mode 100644 index 000000000..d436abb71 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromBdvChannelViewDemo.java @@ -0,0 +1,36 @@ +package org.mastodon.mamut.io.importer.labelimage.ui; + +import net.imglib2.img.Img; +import net.imglib2.type.numeric.real.FloatType; +import org.mastodon.mamut.MainWindow; +import org.mastodon.mamut.ProjectModel; +import org.mastodon.mamut.io.importer.labelimage.util.DemoUtils; +import org.mastodon.mamut.model.Model; +import org.mastodon.mamut.views.bdv.MamutViewBdv; +import org.scijava.Context; +import org.scijava.command.CommandService; +import org.scijava.ui.UIService; + +public class ImportSpotsFromBdvChannelViewDemo +{ + + public static void main( String[] args ) + { + @SuppressWarnings( "all" ) + Context context = new Context(); + UIService ui = context.service( UIService.class ); + CommandService cmd = context.service( CommandService.class ); + + Img< FloatType > image = DemoUtils.generateExampleImage(); + + // show ImageJ + ui.showUI(); + // open the image in Mastodon + Model model = new Model(); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( image, model, context ); + new MainWindow( projectModel ).setVisible( true ); + projectModel.getWindowManager().createView( MamutViewBdv.class ); + // run import spots command + cmd.run( ImportSpotsFromBdvChannelView.class, true, "projectModel", projectModel ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusViewDemo.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusViewDemo.java new file mode 100644 index 000000000..9d46816c9 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/ui/ImportSpotsFromImgPlusViewDemo.java @@ -0,0 +1,44 @@ +package org.mastodon.mamut.io.importer.labelimage.ui; + +import ij.ImagePlus; +import net.imglib2.img.Img; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.numeric.real.FloatType; + +import org.mastodon.mamut.MainWindow; +import org.mastodon.mamut.ProjectModel; +import org.mastodon.mamut.io.importer.labelimage.util.DemoUtils; +import org.mastodon.mamut.model.Model; +import org.mastodon.mamut.views.bdv.MamutViewBdv; +import org.scijava.Context; +import org.scijava.command.CommandService; +import org.scijava.ui.UIService; + +public class ImportSpotsFromImgPlusViewDemo +{ + + public static void main( String[] args ) + { + @SuppressWarnings( "all" ) + Context context = new Context(); + UIService ui = context.service( UIService.class ); + CommandService cmd = context.service( CommandService.class ); + + Img< FloatType > image = DemoUtils.generateExampleImage(); + + // show ImageJ + ui.showUI(); + // show image in ImageJ + ImagePlus imagePlus = ImageJFunctions.wrap( image, "label image" ); + imagePlus.setDimensions( 1, 100, 2 ); + imagePlus.setZ( 50 ); + imagePlus.show(); + // open the image in Mastodon + Model model = new Model(); + ProjectModel projectModel = DemoUtils.wrapAsAppModel( image, model, context ); + new MainWindow( projectModel ).setVisible( true ); + projectModel.getWindowManager().createView( MamutViewBdv.class ); + // run import spots command + cmd.run( ImportSpotsFromImgPlusView.class, true, "projectModel", projectModel ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/ComputeMeanAndVarianceDemo.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/ComputeMeanAndVarianceDemo.java new file mode 100644 index 000000000..7761908e2 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/ComputeMeanAndVarianceDemo.java @@ -0,0 +1,156 @@ +/*- + * #%L + * mastodon-ellipsoid-fitting + * %% + * Copyright (C) 2015 - 2023 Tobias Pietzsch, Jean-Yves Tinevez + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +package org.mastodon.mamut.io.importer.labelimage.util; + +import mpicbg.spim.data.sequence.DefaultVoxelDimensions; +import mpicbg.spim.data.sequence.VoxelDimensions; +import net.imagej.patcher.LegacyInjector; +import net.imglib2.Cursor; +import net.imglib2.img.Img; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Cast; +import net.imglib2.util.LinAlgHelpers; +import net.imglib2.util.Pair; +import net.imglib2.util.StopWatch; +import org.mastodon.mamut.io.importer.labelimage.LabelImageUtils; +import org.mastodon.mamut.model.Model; +import org.scijava.Context; + +import java.util.Arrays; + +/** + * Computing the mean position and covariance matrix for a given segmented + * region of an image is an easy way to get good ellipsoid parameters for + * that segment. + *
+ * Here is an example of how to do that. + */ +public class ComputeMeanAndVarianceDemo +{ + + public static void main( String[] args ) + { + LegacyInjector.preinit(); + try (final Context context = new Context()) + { + double[] center = { 40, 50, 60 }; + double[][] givenCovariance = { + { 400, 20, -10 }, + { 20, 200, 30 }, + { -10, 30, 100 } + }; + + long[] dimensions = { 100, 100, 100 }; + int background = 0; + int pixelValue = 1; + Img< FloatType > image = DemoUtils.generateExampleImage( center, givenCovariance, dimensions, background, pixelValue ); + + StopWatch stopWatchOnline = StopWatch.createAndStart(); + Pair< double[], double[][] > results = + DemoUtils.computeMeanCovarianceOnline( Cast.unchecked( image ), pixelValue, 2.1 ); + double[] onlineMean = results.getA(); + double[][] onlineCovariance = results.getB(); + stopWatchOnline.stop(); + + StopWatch stopWatchTwoPass = StopWatch.createAndStart(); + double[] mean = computeMean( image, pixelValue ); + double[][] computedCovariance = computeCovariance( image, mean, pixelValue ); + stopWatchTwoPass.stop(); + + System.out.println( "Given center: " + Arrays.toString( center ) ); + System.out.println( "Computed mean: " + Arrays.toString( mean ) ); + System.out.println( "Computed mean online: " + Arrays.toString( onlineMean ) ); + System.out.println( "Given covariance: " + Arrays.deepToString( givenCovariance ) ); + System.out.println( "Computed covariance: " + Arrays.deepToString( computedCovariance ) ); + System.out.println( "Computed covariance online: " + Arrays.deepToString( onlineCovariance ) ); + System.out.println( "Time to compute mean and covariance: \n" + stopWatchTwoPass.nanoTime() / 1e9 + " nano seconds" ); + System.out.println( "Time to compute mean and covariance online: \n" + stopWatchOnline.nanoTime() / 1e9 + " nano seconds" ); + + Model model = new Model(); + model.getGraph().addVertex().init( 0, onlineMean, onlineCovariance ); + DemoUtils.showBdvWindow( DemoUtils.wrapAsAppModel( image, model, context ) ); + } + } + + /** + * Computes the mean position of the pixels whose value equals the given {@code pixelValue}. + * + * @param image the image + * @param pixelValue the pixel value + * @return the mean position + */ + private static double[] computeMean( final Img< FloatType > image, final int pixelValue ) + { + Cursor< FloatType > cursor = image.cursor(); + double[] sum = new double[ 3 ]; + double[] position = new double[ 3 ]; + long counter = 0; + while ( cursor.hasNext() ) + if ( cursor.next().get() == pixelValue ) + { + cursor.localize( position ); + LinAlgHelpers.add( sum, position, sum ); + counter++; + } + LinAlgHelpers.scale( sum, 1. / counter, sum ); + return sum; + } + + /** + * Computes the covariance matrix of the pixels whose value equals the given {@code pixelValue}. + * Uses a two-pass algorithm to compute the covariance matrix, cf. Two-pass algorithm for covariance + * + * @param image the image + * @param mean the mean position + * @param pixelValue the pixel value + */ + private static double[][] computeCovariance( final Img< FloatType > image, final double[] mean, final int pixelValue ) + { + Cursor< FloatType > cursor = image.cursor(); + long counter = 0; + double[] position = new double[ 3 ]; + double[][] covariance = new double[ 3 ][ 3 ]; + cursor.reset(); + while ( cursor.hasNext() ) + if ( cursor.next().get() == pixelValue ) + { + cursor.localize( position ); + LinAlgHelpers.subtract( position, mean, position ); + for ( int i = 0; i < 3; i++ ) + for ( int j = 0; j < 3; j++ ) + covariance[ i ][ j ] += position[ i ] * position[ j ]; + counter++; + } + VoxelDimensions voxelDimensions = new DefaultVoxelDimensions( 3 ); + LabelImageUtils.scale( covariance, Math.sqrt( 1d / counter ), voxelDimensions ); + // I don't know why the factor 5 is needed. But it works. + LabelImageUtils.scale( covariance, Math.sqrt( 5d ), voxelDimensions ); + return covariance; + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/DemoUtils.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/DemoUtils.java new file mode 100644 index 000000000..bc92f0885 --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/DemoUtils.java @@ -0,0 +1,156 @@ +/*- + * #%L + * mastodon-ellipsoid-fitting + * %% + * Copyright (C) 2015 - 2023 Tobias Pietzsch, Jean-Yves Tinevez + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +package org.mastodon.mamut.io.importer.labelimage.util; + +import ij.ImagePlus; +import mpicbg.spim.data.sequence.DefaultVoxelDimensions; +import net.imagej.ImgPlus; +import net.imagej.axis.Axes; +import net.imagej.axis.AxisType; +import net.imglib2.Cursor; +import net.imglib2.img.Img; +import net.imglib2.img.ImgView; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.display.imagej.ImgToVirtualStack; +import net.imglib2.loops.LoopBuilder; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Pair; +import net.imglib2.util.ValuePair; +import net.imglib2.view.Views; +import org.mastodon.mamut.ProjectModel; +import org.mastodon.mamut.io.importer.labelimage.LabelImageUtils; +import org.mastodon.mamut.io.importer.labelimage.math.CovarianceMatrix; +import org.mastodon.mamut.io.importer.labelimage.math.MeansVector; +import org.mastodon.mamut.model.Model; +import org.mastodon.mamut.views.bdv.MamutViewBdv; +import org.mastodon.views.bdv.SharedBigDataViewerData; +import org.scijava.Context; + +import javax.annotation.Nonnull; +import java.util.Arrays; +import java.util.List; +import java.util.Objects; + +public class DemoUtils +{ + private DemoUtils() + { + // prevent from instantiation + } + + public static ProjectModel wrapAsAppModel( final Img< FloatType > image, final Model model, final Context context ) + { + final SharedBigDataViewerData sharedBigDataViewerData = asSharedBdvDataXyz( image ); + return ProjectModel.create( context, model, sharedBigDataViewerData, null ); + } + + public static SharedBigDataViewerData asSharedBdvDataXyz( final Img< FloatType > image1 ) + { + final ImagePlus image = + ImgToVirtualStack.wrap( new ImgPlus<>( image1, "image", new AxisType[] { Axes.X, Axes.Y, Axes.Z, Axes.TIME } ) ); + return Objects.requireNonNull( SharedBigDataViewerData.fromImagePlus( image ) ); + } + + public static void showBdvWindow( @Nonnull final ProjectModel appModel ) + { + appModel.getWindowManager().createView( MamutViewBdv.class ); + } + + /** + * Returns an example image with a single ellipsoid. + * + * @param center center of the ellipsoid + * @param cov covariance matrix of the ellipsoid + * @param dimensions dimensions of the image + * @param background value of the background + * @param pixelValue value of the ellipsoid + */ + public static Img< FloatType > generateExampleImage( + final double[] center, final double[][] cov, final long[] dimensions, final int background, final int pixelValue + ) + { + Img< FloatType > image = ArrayImgs.floats( dimensions ); + MultiVariateNormalDistributionRenderer.renderMultivariateNormalDistribution( center, cov, image ); + LoopBuilder.setImages( image ).forEachPixel( pixel -> { + if ( pixel.get() > 500 ) + pixel.set( pixelValue ); + else + pixel.set( background ); + } ); + return image; + } + + /** + * Returns an example image with a single ellipsoid and black background. + */ + public static Img< FloatType > generateExampleImage() + { + double[] center = { 40, 80, 50 }; + double[][] givenCovariance = { + { 400, 20, -10 }, + { 20, 200, 30 }, + { -10, 30, 100 } + }; + long[] dimensions = { 100, 100, 100 }; + int background = 0; + int pixelValue = 1; + Img< FloatType > frame = generateExampleImage( center, givenCovariance, dimensions, background, pixelValue ); + List< Img< FloatType > > twoIdenticalFrames = Arrays.asList( frame, frame ); + return ImgView.wrap( Views.stack( twoIdenticalFrames ) ); + } + + /** + * Computes the covariance matrix of the pixels whose value equals the given {@code pixelValue}. + * Uses an online algorithm to compute the covariance matrix, cf. Online algorithm for covariance + * + * @param image the image + * @param pixelValue the pixel value + */ + public static Pair< double[], double[][] > computeMeanCovarianceOnline( final Img< RealType< ? > > image, final int pixelValue, + double sigma ) + { + Cursor< RealType< ? > > cursor = image.cursor(); + int[] position = new int[ 3 ]; + cursor.reset(); + MeansVector mean = new MeansVector( 3 ); + CovarianceMatrix cov = new CovarianceMatrix( 3 ); + while ( cursor.hasNext() ) + if ( cursor.next().getRealDouble() == pixelValue ) + { + cursor.localize( position ); + mean.addValues( position ); + cov.addValues( position ); + } + double[] means = mean.get(); + double[][] covariances = cov.get(); + LabelImageUtils.scale( covariances, sigma, new DefaultVoxelDimensions( 3 ) ); + return new ValuePair<>( means, covariances ); + } +} diff --git a/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/MultiVariateNormalDistributionRenderer.java b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/MultiVariateNormalDistributionRenderer.java new file mode 100644 index 000000000..045eccbbf --- /dev/null +++ b/src/test/java/org/mastodon/mamut/io/importer/labelimage/util/MultiVariateNormalDistributionRenderer.java @@ -0,0 +1,77 @@ +/*- + * #%L + * mastodon-ellipsoid-fitting + * %% + * Copyright (C) 2015 - 2023 Tobias Pietzsch, Jean-Yves Tinevez + * %% + * Redistribution and use in source and binary forms, with or without + * modification, are permitted provided that the following conditions are met: + * + * 1. Redistributions of source code must retain the above copyright notice, + * this list of conditions and the following disclaimer. + * 2. Redistributions in binary form must reproduce the above copyright notice, + * this list of conditions and the following disclaimer in the documentation + * and/or other materials provided with the distribution. + * + * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" + * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE + * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE + * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE + * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR + * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF + * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS + * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN + * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) + * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE + * POSSIBILITY OF SUCH DAMAGE. + * #L% + */ +package org.mastodon.mamut.io.importer.labelimage.util; + +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.loops.LoopBuilder; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Intervals; +import net.imglib2.util.LinAlgHelpers; + +public class MultiVariateNormalDistributionRenderer +{ + + /** + * Renders the density function of a multivariate normal distribution into a given image. + * @see Wikipedia Multivariate normal distribution + * @param center center of the distribution + * @param cov covariance matrix of the distribution (must be symmetric and positive definite) + * @param image the image to render into (image is a cube) + * + */ + public static void renderMultivariateNormalDistribution( double[] center, double[][] cov, + RandomAccessibleInterval< FloatType > image ) + { + AffineTransform3D sigma = new AffineTransform3D(); + sigma.set( + cov[ 0 ][ 0 ], cov[ 0 ][ 1 ], cov[ 0 ][ 2 ], 0, + cov[ 1 ][ 0 ], cov[ 1 ][ 1 ], cov[ 1 ][ 2 ], 0, + cov[ 2 ][ 0 ], cov[ 2 ][ 1 ], cov[ 2 ][ 2 ], 0 + ); + double[] coord = new double[ 3 ]; + double[] out = new double[ 3 ]; + LoopBuilder.setImages( Intervals.positions( image ), image ).forEachPixel( ( position, pixel ) -> { + position.localize( coord ); + LinAlgHelpers.subtract( coord, center, coord ); + sigma.applyInverse( out, coord ); + // leave out the 1 / (sqrt( ( 2 * pi ) ^ 3 * det( cov )) factor to make the image more visible + double value = Math.exp( -0.5 * scalarProduct( coord, out ) ); + pixel.setReal( 1000 * value ); + } ); + } + + /** + * Computes the scalar product of two vectors. + */ + public static double scalarProduct( double[] a, double[] b ) + { + return a[ 0 ] * b[ 0 ] + a[ 1 ] * b[ 1 ] + a[ 2 ] * b[ 2 ]; + } +} diff --git a/src/test/java/org/mastodon/mamut/util/LineageTreeUtilsTest.java b/src/test/java/org/mastodon/mamut/util/LineageTreeUtilsTest.java index c2394a95e..a373395c4 100644 --- a/src/test/java/org/mastodon/mamut/util/LineageTreeUtilsTest.java +++ b/src/test/java/org/mastodon/mamut/util/LineageTreeUtilsTest.java @@ -34,6 +34,7 @@ import org.mastodon.collection.RefSet; import org.mastodon.mamut.feature.branch.exampleGraph.ExampleGraph2; import org.mastodon.mamut.feature.branch.exampleGraph.ExampleGraph4; +import org.mastodon.mamut.feature.branch.exampleGraph.ExampleGraph6; import org.mastodon.mamut.model.Link; import org.mastodon.mamut.model.Model; import org.mastodon.mamut.model.Spot; @@ -57,11 +58,14 @@ public class LineageTreeUtilsTest private ExampleGraph4 graph4; + private ExampleGraph6 graph6; + @Before public void setUp() { graph2 = new ExampleGraph2(); graph4 = new ExampleGraph4(); + graph6 = new ExampleGraph6(); } @Test @@ -134,4 +138,26 @@ public void testGetAllEdgeSuccessors() actual = LineageTreeUtils.getAllEdgeSuccessors( graph4.branchSpotD, graph4.getModel().getBranchGraph() ); assertEquals( expected, actual ); } + + @Test + public void testLinkSpotsWithSameLabel() + { + assertEquals( 0, graph6.getModel().getGraph().edges().size() ); + LineageTreeUtils.linkSpotsWithSameLabel( graph6.getModel() ); + assertEquals( 7, graph6.getModel().getGraph().edges().size() ); + assertSpotEquals( graph6.spot0.outgoingEdges().get( 0 ).getTarget(), graph6.spot1 ); + assertSpotEquals( graph6.spot1.outgoingEdges().get( 0 ).getTarget(), graph6.spot2 ); + assertSpotEquals( graph6.spot3.outgoingEdges().get( 0 ).getTarget(), graph6.spot4 ); + assertSpotEquals( graph6.spot4.outgoingEdges().get( 0 ).getTarget(), graph6.spot5 ); + assertSpotEquals( graph6.spot6.outgoingEdges().get( 0 ).getTarget(), graph6.spot7 ); + assertEquals( 2, graph6.spot2.outgoingEdges().size() ); + assertSpotEquals( graph6.spot2.outgoingEdges().get( 0 ).getTarget(), graph6.spot10 ); + assertSpotEquals( graph6.spot2.outgoingEdges().get( 1 ).getTarget(), graph6.spot9 ); + assertEquals( 0, graph6.spot7.outgoingEdges().size() ); + } + + private void assertSpotEquals( final Spot expected, final Spot actual ) + { + assertEquals( expected, actual ); + } }