diff --git a/core/src/main/java/tgmm/Ellipsoid.java b/core/src/main/java/tgmm/Ellipsoid.java new file mode 100644 index 0000000000000000000000000000000000000000..9995887ada5b1a8d6e0b408bcdf37a3b530b5c35 --- /dev/null +++ b/core/src/main/java/tgmm/Ellipsoid.java @@ -0,0 +1,59 @@ +package tgmm; + +import net.imglib2.RealPoint; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.util.LinAlgHelpers; +import tgmm.TgmmXmlReader.Gaussian; +import Jama.EigenvalueDecomposition; +import Jama.Matrix; + +public class Ellipsoid extends RealPoint +{ + final double[] m; + + final double[][] W; + + final double[][] S; + + final double boundingSphereRadiusSquared; + + public Ellipsoid( final Gaussian gaussian, final AffineTransform3D transform ) + { + super( 3 ); + m = position; + + final double[] wtmp = new double[9]; + LinAlgHelpers.scale( gaussian.getW(), gaussian.getNu(), wtmp ); + + final Matrix precMat = new Matrix( wtmp, 3 ); + final Matrix covMat = precMat.inverse(); + S = covMat.getArray(); + + transform.apply( gaussian.getM(), m ); + + final double[][] T = new double[3][3]; + for ( int r = 0; r < 3; ++r ) + for ( int c = 0; c < 3; ++c ) + T[r][c] = transform.get( r, c ); + + final double[][] TS = new double[3][3]; + LinAlgHelpers.mult( T, S, TS ); + LinAlgHelpers.multABT( TS, T, S ); + + W = covMat.inverse().getArray(); + + final EigenvalueDecomposition eig = covMat.eig(); + final double[] eigVals = eig.getRealEigenvalues(); + double max = 0; + for ( int k = 0; k < eigVals.length; k++ ) + max = Math.max( max, eigVals[ k ] ); + boundingSphereRadiusSquared = max * EllipsoidRealRandomAccessible.nSigmasSquared; + } + + public static double squaredMahalanobis( final double[] x, final double[] m, final double[][] W, final double[] tmp1, final double[] tmp2 ) + { + LinAlgHelpers.subtract( x, m, tmp1 ); + LinAlgHelpers.mult( W, tmp1, tmp2 ); + return LinAlgHelpers.dot( tmp1, tmp2 ); + } +} \ No newline at end of file diff --git a/core/src/main/java/tgmm/EllipsoidRealRandomAccessible.java b/core/src/main/java/tgmm/EllipsoidRealRandomAccessible.java new file mode 100644 index 0000000000000000000000000000000000000000..edafe58cdba34c6ae1734c0872f4b6338d089574 --- /dev/null +++ b/core/src/main/java/tgmm/EllipsoidRealRandomAccessible.java @@ -0,0 +1,112 @@ +package tgmm; + +import java.util.List; + +import net.imglib2.RealInterval; +import net.imglib2.RealPoint; +import net.imglib2.RealRandomAccess; +import net.imglib2.RealRandomAccessible; +import net.imglib2.collection.KDTree; +import net.imglib2.neighborsearch.RadiusNeighborSearchOnKDTree; +import net.imglib2.type.numeric.integer.UnsignedShortType; + +public class EllipsoidRealRandomAccessible implements RealRandomAccessible< UnsignedShortType > +{ + final static double nSigmas = 2; + + final static double nSigmasSquared = nSigmas * nSigmas; + + final KDTree< Ellipsoid > tree; + + final double maxRadius; + + public EllipsoidRealRandomAccessible( final List< Ellipsoid > ellipsoids ) + { + tree = new KDTree< Ellipsoid >( ellipsoids, ellipsoids ); + double max = 0; + for ( final Ellipsoid e : ellipsoids ) + max = Math.max( max, e.boundingSphereRadiusSquared ); + maxRadius = Math.sqrt( max ); + } + + @Override + public int numDimensions() + { + return 3; + } + + @Override + public RealRandomAccess< UnsignedShortType > realRandomAccess( final RealInterval interval ) + { + return realRandomAccess(); + } + + class Access extends RealPoint implements RealRandomAccess< UnsignedShortType > + { + private final RadiusNeighborSearchOnKDTree< Ellipsoid > search; + + private final UnsignedShortType t = new UnsignedShortType(); + + final double[] tmp1 = new double[ 3 ]; + + final double[] tmp2 = new double[ 3 ]; + + public Access() + { + super( 3 ); + search = new RadiusNeighborSearchOnKDTree< Ellipsoid >( tree ); + } + + private KDTree< Ellipsoid > tree() + { + return tree; + } + + protected Access( final Access a ) + { + super( a ); + search = new RadiusNeighborSearchOnKDTree< Ellipsoid >( a.tree() ); + } + + @Override + public UnsignedShortType get() + { +// boolean inBoundingSphere = false; + search.search( this, maxRadius, false ); + final int numNeighbors = search.numNeighbors(); + for ( int i = 0; i < numNeighbors; ++i ) + { + final Ellipsoid ell = search.getSampler( i ).get(); + if ( search.getSquareDistance( i ) < ell.boundingSphereRadiusSquared ) + { +// inBoundingSphere = true; + if ( Ellipsoid.squaredMahalanobis( this.position, ell.m, ell.W, tmp1, tmp2 ) < nSigmasSquared ) + { + t.set( 1000 ); + return t; + } + } + } + t.set( 0 ); + return t; + } + + @Override + public Access copy() + { + return new Access( this ); + } + + @Override + public Access copyRealRandomAccess() + { + return copy(); + } + } + + @Override + public RealRandomAccess< UnsignedShortType > realRandomAccess() + { + return new Access(); + } +} diff --git a/core/src/main/java/tgmm/ShowEllipsoids.java b/core/src/main/java/tgmm/ShowEllipsoids.java new file mode 100644 index 0000000000000000000000000000000000000000..8a2e88a0af651c3c81866504ce96a578edbaedfd --- /dev/null +++ b/core/src/main/java/tgmm/ShowEllipsoids.java @@ -0,0 +1,157 @@ +package tgmm; + +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; +import java.util.HashMap; + +import mpicbg.spim.data.SpimDataException; +import net.imglib2.Interval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealLocalizable; +import net.imglib2.RealPoint; +import net.imglib2.RealRandomAccessible; +import net.imglib2.display.RealARGBColorConverter; +import net.imglib2.display.RealARGBColorConverter.Imp1; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.ARGBType; +import net.imglib2.type.numeric.integer.UnsignedShortType; +import net.imglib2.util.Intervals; +import net.imglib2.view.Views; + +import org.jdom2.JDOMException; + +import tgmm.TgmmXmlReader.Gaussian; +import bdv.BigDataViewer; +import bdv.export.ProgressWriterConsole; +import bdv.tools.brightness.RealARGBColorConverterSetup; +import bdv.util.IntervalBoundingBox; +import bdv.util.ModifiableInterval; +import bdv.viewer.Interpolation; +import bdv.viewer.Source; +import bdv.viewer.SourceAndConverter; +import bdv.viewer.state.ViewerState; + +public class ShowEllipsoids +{ + public static class EllipsoidsSource implements Source< UnsignedShortType > + { + private final ModifiableInterval interval; + + private final HashMap< Integer, EllipsoidRealRandomAccessible > timepointIndexToEllipsoidAccessible; + + private final String name; + + private final UnsignedShortType type; + + protected final AffineTransform3D identity; + + public EllipsoidsSource( final HashMap< Integer, EllipsoidRealRandomAccessible > timepointIndexToEllipsoidAccessible, final String name ) + { + this.timepointIndexToEllipsoidAccessible = timepointIndexToEllipsoidAccessible; + this.name = name; + this.interval = new ModifiableInterval( Intervals.createMinMax( -100, -100, -100, 100, 100, 100 ) ); + this.type = new UnsignedShortType(); + this.identity = new AffineTransform3D(); + } + + public final ModifiableInterval getInterval( final int t, final int level ) + { + return interval; + } + + @Override + public boolean isPresent( final int t ) + { + return timepointIndexToEllipsoidAccessible.containsKey( t ); + } + + @Override + public RandomAccessibleInterval< UnsignedShortType > getSource( final int t, final int level ) + { + return Views.interval( Views.raster( timepointIndexToEllipsoidAccessible.get( t ) ), interval ); + } + + @Override + public RealRandomAccessible< UnsignedShortType > getInterpolatedSource( final int t, final int level, final Interpolation method ) + { + return timepointIndexToEllipsoidAccessible.get( t ); + } + + @Override + public AffineTransform3D getSourceTransform( final int t, final int level ) + { + return identity; + } + + @Override + public UnsignedShortType getType() + { + return type; + } + + @Override + public String getName() + { + return name; + } + + @Override + public int getNumMipmapLevels() + { + return 1; + } + } + + public static void main( final String[] args ) throws IOException, JDOMException, SpimDataException + { + final int timepointIndex = 245; + final String fn = "/Users/pietzsch/desktop/data/BDV130418A325/BDV130418A325_NoTempReg.xml"; + + System.setProperty( "apple.laf.useScreenMenuBar", "true" ); + final BigDataViewer bdv = new BigDataViewer( fn, new File( fn ).getName(), new ProgressWriterConsole() ); + bdv.getViewer().setTimepoint( timepointIndex ); + + final ViewerState viewerState = bdv.getViewer().getState(); + final int currentSourceId = viewerState.getCurrentSource(); + final Source< ? > imgSource = viewerState.getSources().get( currentSourceId ).getSpimSource(); + final Interval imgSourceInterval = imgSource.getSource( timepointIndex, 0 ); + + final HashMap< Integer, EllipsoidRealRandomAccessible > timepointIndexToEllipsoidAccessible = new HashMap< Integer, EllipsoidRealRandomAccessible >(); + for ( int tp = 240; tp < 248; ++tp ) + { + final int timepointId = tp + 1; + final String tgmmFilename = "/Users/pietzsch/Desktop/data/Fernando/extract/GMEMfinalResult_frame0" + timepointId + ".xml"; + System.out.println( tgmmFilename ); + final ArrayList< Gaussian > gaussians = TgmmXmlReader.read( tgmmFilename ); + final ArrayList< Ellipsoid > ellipsoids = new ArrayList< Ellipsoid >(); + final AffineTransform3D imgSourceTransform = imgSource.getSourceTransform( tp, 0 ); + for ( final Gaussian g : gaussians ) + ellipsoids.add( new Ellipsoid( g, imgSourceTransform ) ); + final EllipsoidRealRandomAccessible ellipsoidsAccessible = new EllipsoidRealRandomAccessible( ellipsoids ); + timepointIndexToEllipsoidAccessible.put( tp, ellipsoidsAccessible ); + } + + final EllipsoidsSource ellipsoidsSource = new EllipsoidsSource( timepointIndexToEllipsoidAccessible, "ellipsoids" ); + final Imp1< UnsignedShortType > converter = new RealARGBColorConverter.Imp1< UnsignedShortType >( 0, 3000 ); + converter.setColor( new ARGBType( ARGBType.rgba( 255, 0, 0, 0 ) ) ); + final SourceAndConverter< UnsignedShortType > ellipsoidsSourceAndConverter = new SourceAndConverter< UnsignedShortType >( ellipsoidsSource, converter ); + + final ArrayList< RealPoint > sourceCorners = new ArrayList< RealPoint >(); + final AffineTransform3D imgSourceTransform = imgSource.getSourceTransform( timepointIndex, 0 ); + for ( final RealLocalizable corner : IntervalBoundingBox.getCorners( imgSourceInterval ) ) + { + final RealPoint sourceCorner = new RealPoint( 3 ); + imgSourceTransform.apply( corner, sourceCorner ); + sourceCorners.add( sourceCorner ); + } + final Interval bb = Intervals.smallestContainingInterval( IntervalBoundingBox.getBoundingBox( sourceCorners ) ); + ellipsoidsSource.getInterval( timepointIndex, 0 ).set( bb ); + + final int ellipsoidSetupId = 12469; // some non-existing setup id + final RealARGBColorConverterSetup cs = new RealARGBColorConverterSetup( ellipsoidSetupId, converter ); + cs.setViewer( bdv.getViewer() ); + bdv.getSetupAssignments().addSetup( cs ); + bdv.getViewer().addSource( ellipsoidsSourceAndConverter ); + } +} diff --git a/core/src/main/java/tgmm/TgmmXmlReader.java b/core/src/main/java/tgmm/TgmmXmlReader.java new file mode 100644 index 0000000000000000000000000000000000000000..759b0b8df5e94bb5bf0322db2c2f52dd6d46ca4f --- /dev/null +++ b/core/src/main/java/tgmm/TgmmXmlReader.java @@ -0,0 +1,159 @@ +package tgmm; + +import java.io.IOException; +import java.util.ArrayList; +import java.util.List; + +import org.jdom2.Document; +import org.jdom2.Element; +import org.jdom2.JDOMException; +import org.jdom2.input.SAXBuilder; + +public class TgmmXmlReader +{ + public static void main( final String[] args ) + { + try + { + final ArrayList< Gaussian > gaussians = read( "/Users/pietzsch/Desktop/data/Fernando/extract/GMEMfinalResult_frame0249.xml" ); +// for ( final Gaussian g : gaussians ) +// System.out.println( g ); + } + catch ( final IOException e ) + { + e.printStackTrace(); + } + catch ( final JDOMException e ) + { + e.printStackTrace(); + } + } + + public static class Gaussian + { + private final double nu; + + private final double[] m; + + private final double[] W; + + private final int id; + + private final int lineage; + + private final int parent; + + public Gaussian( + final double nu, + final double[] m, + final double[] W, + final int id, + final int lineage, + final int parent ) + { + this.nu = nu; + this.m = m; + this.W = W; + this.id = id; + this.lineage = lineage; + this.parent = parent; + } + + public double getNu() + { + return nu; + } + + public double[] getM() + { + return m; + } + + public double[] getW() + { + return W; + } + + public int getId() + { + return id; + } + + public int getLineage() + { + return lineage; + } + + public int getParent() + { + return parent; + } + + @Override + public String toString() + { + return "Gaussian( " + + "nu = " + nu + + "m = " + net.imglib2.util.Util.printCoordinates( m ) + + "W = " + net.imglib2.util.Util.printCoordinates( W ) + + "id = " + id + + "lineage = " + lineage + + "parent = " + parent + + " )"; + } + } + + public static double getDoubleAttribute( final Element parent, final String name ) + { + return Double.parseDouble( parent.getAttributeValue( name ) ); + } + + public static double[] getDoubleArrayAttribute( final Element parent, final String name ) + { + try + { + final String text = parent.getAttributeValue( name ); + final String[] entries = text.split( "\\s+" ); + final double[] array = new double[ entries.length ]; + for ( int i = 0; i < entries.length; ++i ) + array[ i ] = Double.parseDouble( entries[ i ] ); + return array; + } + catch ( final Exception e ) + { + final String text = parent.getAttributeValue( name ); + System.out.println( text ); + final String[] entries = text.split( "\\s+" ); + for ( int i = 0; i < entries.length; ++i ) + System.out.println( entries[ i ] ); + return null; + } + } + + public static int getIntAttribute( final Element parent, final String name ) + { + return Integer.parseInt( parent.getAttributeValue( name ) ); + } + + public static ArrayList< Gaussian > read( final String xmlFilename ) throws IOException, JDOMException + { + final SAXBuilder sax = new SAXBuilder(); + final Document doc = sax.build( xmlFilename ); + final Element root = doc.getRootElement(); + + final List< Element > gaussianMixtureModels = root.getChildren( "GaussianMixtureModel" ); + final ArrayList< Gaussian > gaussians = new ArrayList< Gaussian >(); + for ( final Element elem : gaussianMixtureModels ) + { + final double nu = getDoubleAttribute( elem, "nu" ); + final double[] m = getDoubleArrayAttribute( elem, "m" ); + final double[] W = getDoubleArrayAttribute( elem, "W" ); + final int id = getIntAttribute( elem, "id" ); + final int lineage = getIntAttribute( elem, "lineage" ); + final int parent = getIntAttribute( elem, "parent" ); + + gaussians.add( new Gaussian( nu, m, W, id, lineage, parent ) ); + } + return gaussians; + } +}