loading and displaying TGMM ellipsoids. Cannot deal with #QNAN etc

core/src/main/java/tgmm/Ellipsoid.java

+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

core/src/main/java/tgmm/EllipsoidRealRandomAccessible.java

+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();
+	}
+}

core/src/main/java/tgmm/ShowEllipsoids.java

+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 );
+	}
+}

core/src/main/java/tgmm/TgmmXmlReader.java

+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;
+	}
+}