From 119e510e90a5a1a1edba966a60cc289226de6d90 Mon Sep 17 00:00:00 2001 From: Stephan Saalfeld Date: Tue, 1 Sep 2015 20:59:16 -0400 Subject: [PATCH] included match intensities plugin, build with SNAPSHOTS for testing mvn -Dimagej.app.directory=$HOME/packages/Fiji.app/ -Denforcer.skip install --- TrakEM2_/pom.xml | 18 +- .../src/main/java/ini/trakem2/display/Display.java | 8 +- .../org/janelia/intensity/LinearIntensityMap.java | 268 ++++++++++++ .../org/janelia/intensity/MatchIntensities.java | 467 +++++++++++++++++++++ .../org/janelia/intensity/PointMatchFilter.java | 20 + .../janelia/intensity/RansacRegressionFilter.java | 48 +++ .../intensity/RansacRegressionReduceFilter.java | 83 ++++ .../main/java/org/janelia/intensity/Render.java | 302 +++++++++++++ pom.xml | 2 +- 9 files changed, 1208 insertions(+), 8 deletions(-) create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/LinearIntensityMap.java create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/MatchIntensities.java create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/PointMatchFilter.java create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionFilter.java create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionReduceFilter.java create mode 100644 TrakEM2_/src/main/java/org/janelia/intensity/Render.java diff --git a/TrakEM2_/pom.xml b/TrakEM2_/pom.xml index 48d79eb9..d90fd335 100644 --- a/TrakEM2_/pom.xml +++ b/TrakEM2_/pom.xml @@ -8,7 +8,7 @@ sc.fiji pom-trakem2 - 1.3.3-SNAPSHOT + 1.3.4-SNAPSHOT TrakEM2_ @@ -71,12 +71,10 @@ mpicbg mpicbg - 1.0.1 - + mpicbg mpicbg_ - 1.0.1 sc.fiji @@ -89,8 +87,7 @@ sc.fiji legacy-imglib1 - 1.1.3-DEPRECATED - + sc.fiji 3D_Viewer @@ -136,6 +133,15 @@ ome formats-bsd + + net.imglib2 + imglib2 + + + net.imglib2 + imglib2-realtransform + + diff --git a/TrakEM2_/src/main/java/ini/trakem2/display/Display.java b/TrakEM2_/src/main/java/ini/trakem2/display/Display.java index 31b68677..bfe24ed9 100644 --- a/TrakEM2_/src/main/java/ini/trakem2/display/Display.java +++ b/TrakEM2_/src/main/java/ini/trakem2/display/Display.java @@ -174,6 +174,8 @@ import mpicbg.trakem2.transform.AffineModel3D; import mpicbg.trakem2.transform.CoordinateTransform; import mpicbg.trakem2.transform.CoordinateTransformList; +import org.janelia.intensity.MatchIntensities; + /** A Display is a class to show a Layer and enable mouse and keyboard manipulation of all its components. */ public final class Display extends DBObject implements ActionListener, IJEventListener { @@ -3233,7 +3235,8 @@ public final class Display extends DBObject implements ActionListener, IJEventLi item = new JMenuItem("Blend (layer-wise)..."); item.addActionListener(this); adjust_menu.add(item); item = new JMenuItem("Blend (selected images)..."); item.addActionListener(this); adjust_menu.add(item); if (selection.isEmpty()) item.setEnabled(false); - popup.add(adjust_menu); + item = new JMenuItem("Match intensities..."); item.addActionListener(this); adjust_menu.add(item); + popup.add(adjust_menu); final JMenu script = new JMenu("Script"); final MenuScriptListener msl = new MenuScriptListener(); @@ -5541,6 +5544,9 @@ public final class Display extends DBObject implements ActionListener, IJEventLi return patch.getTitle().matches(regex); } }); + } else if (command.equals("Match intensities...")) { + final MatchIntensities matching = new MatchIntensities(); + matching.invoke(getActive()); } else if (command.equals("Montage")) { final Set affected = new HashSet(selection.getAffected()); // make an undo step! diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/LinearIntensityMap.java b/TrakEM2_/src/main/java/org/janelia/intensity/LinearIntensityMap.java new file mode 100644 index 00000000..ccd36f4d --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/LinearIntensityMap.java @@ -0,0 +1,268 @@ +/** + * License: GPL + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License 2 + * as published by the Free Software Foundation. + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +package org.janelia.intensity; + +import ij.ImageJ; +import ij.ImagePlus; +import net.imglib2.Cursor; +import net.imglib2.Dimensions; +import net.imglib2.FinalInterval; +import net.imglib2.IterableInterval; +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealRandomAccessible; +import net.imglib2.img.array.ArrayImg; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.basictypeaccess.array.ByteArray; +import net.imglib2.img.basictypeaccess.array.FloatArray; +import net.imglib2.img.basictypeaccess.array.IntArray; +import net.imglib2.img.basictypeaccess.array.ShortArray; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.interpolation.InterpolatorFactory; +import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory; +import net.imglib2.interpolation.randomaccess.NearestNeighborInterpolatorFactory; +import net.imglib2.realtransform.RealViews; +import net.imglib2.realtransform.Scale; +import net.imglib2.realtransform.Translation; +import net.imglib2.type.numeric.ARGBType; +import net.imglib2.type.numeric.NumericType; +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.DoubleType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; +import net.imglib2.view.composite.CompositeIntervalView; +import net.imglib2.view.composite.RealComposite; + +/** + * Transfer intensities by a linear function y=ax+b + * with the coefficients a and b being stored as a + * {@link RealComposite RealComposite<T>} 2d-vector + * (a, b) in a map. The + * map is a {@link RealRandomAccessible} and the defined interval as passed as + * an independent parameter. If the coefficients are passed as a raster, the + * interval is that of a raster. The map is scaled such that it applies to + * the interval of the input image. + * + * @author Stephan Saalfeld + */ +public class LinearIntensityMap< T extends RealType< T > > +{ + static public enum Interpolation{ NN, NL }; + + final static private < T extends RealType< T > >InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory( final Interpolation interpolation ) + { + switch ( interpolation ) + { + case NN: + return new NearestNeighborInterpolatorFactory< RealComposite< T > >(); + default: + return new NLinearInterpolatorFactory< RealComposite< T > >(); + } + } + + final protected Dimensions dimensions; + final protected Translation translation; + final protected RealRandomAccessible< RealComposite< T > > coefficients; + + final protected InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory; + + public LinearIntensityMap( final RandomAccessibleInterval< T > source, final InterpolatorFactory< RealComposite< T >, RandomAccessible< RealComposite< T > > > interpolatorFactory ) + { + this.interpolatorFactory = interpolatorFactory; + final CompositeIntervalView< T, RealComposite< T > > collapsedSource = Views.collapseReal( source ); + dimensions = new FinalInterval( collapsedSource ); + final double[] shift = new double[ dimensions.numDimensions() ]; + for ( int d = 0; d < shift.length; ++d ) + shift[ d ] = 0.5; + translation = new Translation( shift ); + + final RandomAccessible< RealComposite< T > > extendedCollapsedSource = Views.extendBorder( collapsedSource ); + coefficients = Views.interpolate( extendedCollapsedSource, interpolatorFactory ); + } + + public LinearIntensityMap( final RandomAccessibleInterval< T > source ) + { + this( source, new NLinearInterpolatorFactory< RealComposite< T > >() ); + } + + public LinearIntensityMap( final RandomAccessibleInterval< T > source, final Interpolation interpolation ) + { + this( source, LinearIntensityMap.< T >interpolatorFactory( interpolation ) ); + } + + @SuppressWarnings( { "rawtypes", "unchecked" } ) + public < S extends NumericType< S > > void run( final RandomAccessibleInterval< S > image ) + { + assert image.numDimensions() == dimensions.numDimensions() : "Number of dimensions do not match."; + + final double[] s = new double[ dimensions.numDimensions() ]; + for ( int d = 0; d < s.length; ++d ) + s[ d ] = image.dimension( d ) / dimensions.dimension( d ); + final Scale scale = new Scale( s ); + + System.out.println( "translation-n " + translation.numDimensions() ); + + final RandomAccessibleInterval< RealComposite< T > > stretchedCoefficients = + Views.offsetInterval( + Views.raster( + RealViews.transform( + RealViews.transform( + coefficients, + translation ), + scale ) ), + image ); + + /* decide on type which mapping to use */ + final S t = image.randomAccess().get(); + + if ( ARGBType.class.isInstance( t ) ) + mapARGB( Views.flatIterable( ( RandomAccessibleInterval< ARGBType > )image ), Views.flatIterable( stretchedCoefficients ) ); + else if ( RealComposite.class.isInstance( t ) ) + mapComposite( Views.flatIterable( ( RandomAccessibleInterval )image ), Views.flatIterable( stretchedCoefficients ) ); + else if ( RealType.class.isInstance( t ) ) + { + final RealType< ? > r = ( RealType )t; + if ( r.getMinValue() > -Double.MAX_VALUE || r.getMaxValue() < Double.MAX_VALUE ) +// TODO Bug in javac does not enable cast from RandomAccessibleInterval< S > to RandomAccessibleInterval< RealType >, remove when fixed + mapCrop( Views.flatIterable( ( RandomAccessibleInterval< RealType > )( Object )image ), Views.flatIterable( stretchedCoefficients ) ); + else +// TODO Bug in javac does not enable cast from RandomAccessibleInterval< S > to RandomAccessibleInterval< RealType >, remove when fixed + map( Views.flatIterable( ( RandomAccessibleInterval< RealType > )( Object )image ), Views.flatIterable( stretchedCoefficients ) ); + } + + } + + final static protected < S extends RealType< S >, T extends RealType< T > > void map( + final IterableInterval< S > image, + final IterableInterval< RealComposite< T > > coefficients ) + { + final Cursor< S > cs = image.cursor(); + final Cursor< RealComposite< T > > ct = coefficients.cursor(); + + while ( cs.hasNext() ) + { + final S s = cs.next(); + final RealComposite< T > t = ct.next(); + s.setReal( s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() ); + } + } + + final static protected < S extends RealType< S >, T extends RealType< T > > void mapCrop( + final IterableInterval< S > image, + final IterableInterval< RealComposite< T > > coefficients ) + { + final Cursor< S > cs = image.cursor(); + final Cursor< RealComposite< T > > ct = coefficients.cursor(); + final S firstValue = cs.next(); + final double minS = firstValue.getMinValue(); + final double maxS = firstValue.getMaxValue(); + + while ( cs.hasNext() ) + { + final S s = cs.next(); + final RealComposite< T > t = ct.next(); + + s.setReal( Math.max( minS, Math.min( maxS, s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() ) ) ); + } + } + + final static protected < S extends RealType< S >, T extends RealType< T > > void mapComposite( + final IterableInterval< RealComposite< S > > image, + final IterableInterval< RealComposite< T > > coefficients ) + { + final Cursor< RealComposite< S > > cs = image.cursor(); + final Cursor< RealComposite< T > > ct = coefficients.cursor(); + + while ( cs.hasNext() ) + { + final RealComposite< S > c = cs.next(); + final RealComposite< T > t = ct.next(); + + for ( final S s : c ) + s.setReal( s.getRealDouble() * t.get( 0 ).getRealDouble() + t.get( 1 ).getRealDouble() ); + } + } + + final static protected < T extends RealType< T > > void mapARGB( + final IterableInterval< ARGBType > image, + final IterableInterval< RealComposite< T > > coefficients ) + { + final Cursor< ARGBType > cs = image.cursor(); + final Cursor< RealComposite< T > > ct = coefficients.cursor(); + + while ( cs.hasNext() ) + { + final RealComposite< T > t = ct.next(); + final double alpha = t.get( 0 ).getRealDouble(); + final double beta = t.get( 1 ).getRealDouble(); + + final ARGBType s = cs.next(); + final int argb = s.get(); + final int a = ( ( argb >> 24 ) & 0xff ); + final double r = ( ( argb >> 16 ) & 0xff ) * alpha + beta; + final double g = ( ( argb >> 8 ) & 0xff ) * alpha + beta; + final double b = ( argb & 0xff ) * alpha + beta; + + s.set( + ( a << 24 ) | + ( ( r < 0 ? 0 : r > 255 ? 255 : ( int )( r + 0.5 ) ) << 16 ) | + ( ( g < 0 ? 0 : g > 255 ? 255 : ( int )( g + 0.5 ) ) << 8 ) | + ( b < 0 ? 0 : b > 255 ? 255 : ( int )( b + 0.5 ) ) ); + } + } + + public static void main( final String[] args ) + { + new ImageJ(); + + final double[] coefficients = new double[]{ + 0, 2, 4, 8, + 1, 1, 1, 1, + 1, 10, 5, 1, + 1, 1, 1, 1, + + 0, 10, 20, 30, + 40, 50, 60, 70, + 80, 90, 100, 110, + 120, 130, 140, 150 + }; + + final LinearIntensityMap< DoubleType > transform = new LinearIntensityMap< DoubleType >( ArrayImgs.doubles( coefficients, 4, 4, 2 ) ); + + //final ImagePlus imp = new ImagePlus( "http://upload.wikimedia.org/wikipedia/en/2/24/Lenna.png" ); + final ImagePlus imp1 = new ImagePlus( "http://fly.mpi-cbg.de/~saalfeld/Pictures/norway.jpg"); + + final ArrayImg< FloatType, FloatArray > image1 = ArrayImgs.floats( ( float[] )imp1.getProcessor().convertToFloatProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() ); + final ArrayImg< UnsignedByteType, ByteArray > image2 = ArrayImgs.unsignedBytes( ( byte[] )imp1.getProcessor().convertToByteProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() ); + final ArrayImg< UnsignedShortType, ShortArray > image3 = ArrayImgs.unsignedShorts( ( short[] )imp1.getProcessor().convertToShortProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() ); + final ArrayImg< ARGBType, IntArray > image4 = ArrayImgs.argbs( ( int[] )imp1.getProcessor().getPixels(), imp1.getWidth(), imp1.getHeight() ); + + ImageJFunctions.show( ArrayImgs.doubles( coefficients, 4, 4, 2 ) ); + + transform.run( image1 ); + transform.run( image2 ); + transform.run( image3 ); + transform.run( image4 ); + + ImageJFunctions.show( image1 ); + ImageJFunctions.show( image2 ); + ImageJFunctions.show( image3 ); + ImageJFunctions.show( image4 ); + } +} diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/MatchIntensities.java b/TrakEM2_/src/main/java/org/janelia/intensity/MatchIntensities.java new file mode 100644 index 00000000..8b4ca260 --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/MatchIntensities.java @@ -0,0 +1,467 @@ +/** + * + */ +package org.janelia.intensity; + +import ij.IJ; +import ij.ImagePlus; +import ij.ImageStack; +import ij.gui.GenericDialog; +import ij.gui.Roi; +import ij.measure.Calibration; +import ij.process.ColorProcessor; +import ij.process.FloatProcessor; +import ini.trakem2.Project; +import ini.trakem2.display.Display; +import ini.trakem2.display.Displayable; +import ini.trakem2.display.Layer; +import ini.trakem2.display.LayerSet; +import ini.trakem2.display.Patch; +import ini.trakem2.persistence.FSLoader; +import ini.trakem2.plugin.TPlugIn; +import ini.trakem2.utils.Utils; + +import java.awt.Rectangle; +import java.io.File; +import java.util.ArrayList; +import java.util.Collection; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map.Entry; +import java.util.concurrent.ExecutionException; + +import mpicbg.models.Affine1D; +import mpicbg.models.AffineModel1D; +import mpicbg.models.IdentityModel; +import mpicbg.models.IllDefinedDataPointsException; +import mpicbg.models.InterpolatedAffineModel1D; +import mpicbg.models.Model; +import mpicbg.models.NotEnoughDataPointsException; +import mpicbg.models.Point; +import mpicbg.models.PointMatch; +import mpicbg.models.Tile; +import mpicbg.models.TileConfiguration; +import mpicbg.models.TranslationModel1D; +import net.imglib2.img.list.ListImg; +import net.imglib2.img.list.ListRandomAccess; + +/** + * + * @author Stephan Saalfeld + * @author Philipp Hanslovsky + */ +public class MatchIntensities implements TPlugIn +{ + protected LayerSet layerset = null; + + static protected int numCoefficients = 8; + static protected double lambda1 = 0.01; + static protected double lambda2 = 0.01; + static protected double neighborWeight = 0.1; + static protected int radius = 5; + static protected int iterations = 2000; + static protected double scale = -1; + + private Layer currentLayer( final Object... params ) + { + final Layer layer; + if ( params != null && params[ 0 ] != null ) + { + final Object param = params[ 0 ]; + if ( Layer.class.isInstance( param ) ) + layer = ( Layer ) param; + else if ( LayerSet.class.isInstance( param ) ) + layer = ( ( LayerSet ) param ).getLayer( 0 ); + else if ( Displayable.class.isInstance( param ) ) + layer = ( ( Displayable ) param ).getLayer(); + else + layer = null; + } + else + { + final Display front = Display.getFront(); + if ( front == null ) + layer = Project.getProjects().get( 0 ).getRootLayerSet().getLayer( 0 ); + else + layer = front.getLayer(); + } + return layer; + } + + private static Rectangle getRoi( final LayerSet layerset ) + { + final Roi roi; + final Display front = Display.getFront(); + if ( front == null ) + roi = null; + else + roi = front.getRoi(); + if ( roi == null ) + return new Rectangle( 0, 0, ( int ) layerset.getLayerWidth(), ( int ) layerset.getLayerHeight() ); + else + return roi.getBounds(); + } + + private double suggestScale( final List< Layer > layers ) + { + if ( layers.size() < 2 ) + return 0.0; + + final Layer layer1 = layers.get( 0 ); + final Layer layer2 = layers.get( 1 ); + final Calibration calib = layer1.getParent().getCalibration(); + final double width = ( calib.pixelWidth + calib.pixelHeight ) * 0.5; + final double depth = calib.pixelDepth; + + return ( layer2.getZ() - layer1.getZ() ) * width / depth; + } + + @Override + public boolean setup( final Object... params ) + { + if ( params != null && params[ 0 ] != null ) + { + final Object param = params[ 0 ]; + if ( LayerSet.class.isInstance( param ) ) + layerset = ( LayerSet ) param; + else if ( Displayable.class.isInstance( param ) ) + layerset = ( ( Displayable ) param ).getLayerSet(); + else + return false; + } + else + { + final Display front = Display.getFront(); + if ( front == null ) + layerset = Project.getProjects().get( 0 ).getRootLayerSet(); + else + layerset = front.getLayerSet(); + } + return true; + } + + final static protected < T extends Model< T > & Affine1D< T > > HashMap< Patch, ArrayList< Tile< T > > > generateCoefficientsTiles( + final Collection< Patch > patches, + final T template, + final int nCoefficients ) + { + final HashMap< Patch, ArrayList< Tile< T > > > map = new HashMap< Patch, ArrayList< Tile< T > > >(); + for ( final Patch p : patches ) + { + final ArrayList< Tile< T > > coefficientModels = new ArrayList< Tile< T > >(); + for ( int i = 0; i < nCoefficients; ++i ) + coefficientModels.add( new Tile< T >( template.copy() ) ); + + map.put( p, coefficientModels ); + } + return map; + } + + final static protected void identityConnect( final Tile< ? > t1, final Tile< ? > t2, final double weight ) + { + final ArrayList< PointMatch > matches = new ArrayList< PointMatch >(); + matches.add( new PointMatch( new Point( new double[] { 0 } ), new Point( new double[] { 0 } ) ) ); + matches.add( new PointMatch( new Point( new double[] { 1 } ), new Point( new double[] { 1 } ) ) ); + t1.connect( t2, matches ); + } + + @Override + public Object invoke( final Object... params ) + { + if ( !setup( params ) ) + return null; + + final Layer layer = currentLayer( params ); + final GenericDialog gd = new GenericDialog( "Match intensities" ); + Utils.addLayerRangeChoices( layer, gd ); + gd.addMessage( "Layer neighborhood range :" ); + gd.addNumericField( "test_maximally :", radius, 0, 6, "layers" ); + gd.addMessage( "Optimizer :" ); + gd.addNumericField( "iterations :", iterations, 0, 6, "" ); + gd.addNumericField( "scale_regularization :", lambda1, 2, 6, "" ); + gd.addNumericField( "translation_regularization :", lambda2, 2, 6, "" ); + gd.addNumericField( "smoothness_regularization :", neighborWeight, 2, 6, "" ); + gd.showDialog(); + + if ( gd.wasCanceled() ) + return null; + + final List< Layer > layers = layerset.getLayers().subList( gd.getNextChoiceIndex(), gd.getNextChoiceIndex() + 1 ); + radius = ( int ) gd.getNextNumber(); + iterations = ( int )gd.getNextNumber(); + lambda1 = gd.getNextNumber(); + lambda2 = gd.getNextNumber(); + neighborWeight = gd.getNextNumber(); + + try + { + run( layers, radius, scale, numCoefficients, lambda1, lambda2, neighborWeight, getRoi( layerset ) ); + } + catch ( final InterruptedException e ) + { + Utils.log( "Match intensities interrupted." ); + } + catch ( final ExecutionException e ) + { + Utils.log( "Match intenities ExecutiuonException occurred:" ); + e.printStackTrace( System.out ); + } + + return null; + } + + @Override + public boolean applies( final Object ob ) + { + return true; + } + + + /** + * TODO Test! And then desperately multi-thread coefficient collection. + * + * @param layers + * @param radius + * @param scale + * @param numCoefficients + * @param lambda1 + * @param lambda2 + * @param neighborWeight + * @param roi + */ + public < M extends Model< M > & Affine1D< M > > void run( + final List< Layer > layers, + final int radius, + final double scale, + final int numCoefficients, + final double lambda1, + final double lambda2, + final double neighborWeight, + final Rectangle roi ) throws InterruptedException, ExecutionException + { + // final PointMatchFilter filter = new RansacRegressionFilter(); + final PointMatchFilter filter = new RansacRegressionReduceFilter(); + + /* collect patches */ + final ArrayList< Patch > patches = new ArrayList< Patch >(); + for ( final Layer layer : layers ) + patches.addAll( ( ArrayList )layer.getDisplayables( Patch.class, roi ) ); + + /* generate coefficient tiles for all patches + * TODO consider offering alternative models */ + final HashMap< Patch, ArrayList< Tile< ? extends M > > > coefficientsTiles = + ( HashMap ) generateCoefficientsTiles( + patches, + new InterpolatedAffineModel1D< InterpolatedAffineModel1D< AffineModel1D, TranslationModel1D >, IdentityModel >( + new InterpolatedAffineModel1D< AffineModel1D, TranslationModel1D >( + new AffineModel1D(), new TranslationModel1D(), lambda1 ), + new IdentityModel(), lambda2 ), + numCoefficients * numCoefficients ); + + /* completed patches */ + final HashSet< Patch > completedPatches = new HashSet< Patch >(); + + for ( final Patch p1 : patches ) + { + completedPatches.add( p1 ); + + final Rectangle box1 = p1.getBoundingBox().intersection( roi ); + + final ArrayList< Patch > p2s = new ArrayList< Patch >(); + + /* across adjacent layers */ + final int layerIndex = layerset.getLayerIndex( p1.getLayer().getId() ); + for ( int i = layerIndex - radius; i <= layerIndex + radius; ++i ) + { + final Layer layer = layerset.getLayer( i ); + if ( layer != null ) + p2s.addAll( ( Collection ) layer.getDisplayables( Patch.class, box1 ) ); + } + + /* get the coefficient tiles */ + final ArrayList< Tile< ? extends M > > p1CoefficientsTiles = coefficientsTiles.get( p1 ); + + for ( final Patch p2 : p2s ) + { + /* + * if this patch had been processed earlier, all matches are + * already in + */ + if ( completedPatches.contains( p2 ) ) + continue; + + /* render intersection */ + final Rectangle box2 = p2.getBoundingBox(); + final Rectangle box = box1.intersection( box2 ); + + final int w = ( int ) ( box.width * scale + 0.5 ); + final int h = ( int ) ( box.height * scale + 0.5 ); + final int n = w * h; + + final FloatProcessor pixels1 = new FloatProcessor( w, h ); + final FloatProcessor weights1 = new FloatProcessor( w, h ); + final ColorProcessor coefficients1 = new ColorProcessor( w, h ); + final FloatProcessor pixels2 = new FloatProcessor( w, h ); + final FloatProcessor weights2 = new FloatProcessor( w, h ); + final ColorProcessor coefficients2 = new ColorProcessor( w, h ); + + Render.render( p1, numCoefficients, numCoefficients, pixels1, weights1, coefficients1, box.x, box.y, scale ); + Render.render( p2, numCoefficients, numCoefficients, pixels2, weights2, coefficients2, box.x, box.y, scale ); + + /* + * generate a matrix of all coefficients in p1 to all + * coefficients in p2 to store matches + */ + final ArrayList< ArrayList< PointMatch > > list = new ArrayList< ArrayList< PointMatch > >(); + for ( int i = 0; i < numCoefficients * numCoefficients * numCoefficients * numCoefficients; ++i ) + list.add( new ArrayList< PointMatch >() ); + final ListImg< ArrayList< PointMatch > > matrix = new ListImg< ArrayList< PointMatch > >( list, numCoefficients * numCoefficients, numCoefficients * numCoefficients ); + final ListRandomAccess< ArrayList< PointMatch > > ra = matrix.randomAccess(); + + /* + * iterate over all pixels and feed matches into the match + * matrix + */ + for ( int i = 0; i < n; ++i ) + { + final int c1 = coefficients1.get( i ); + if ( c1 > 0 ) + { + final int c2 = coefficients2.get( i ); + if ( c2 > 0 ) + { + final double w1 = weights1.getf( i ); + if ( w1 > 0 ) + { + final double w2 = weights2.getf( i ); + if ( w2 > 0 ) + { + final double p = pixels1.getf( i ); + final double q = pixels2.getf( i ); + final PointMatch pq = new PointMatch( new Point( new double[] { p } ), new Point( new double[] { q } ), w1 * w2 ); + + /* first label is 1 */ + ra.setPosition( c1 - 1, 0 ); + ra.setPosition( c2 - 1, 1 ); + ra.get().add( pq ); + } + } + } + } + } + + /* filter matches */ + final ArrayList< PointMatch > inliers = new ArrayList< PointMatch >(); + for ( final ArrayList< PointMatch > candidates : matrix ) + { + inliers.clear(); + filter.filter( candidates, inliers ); + candidates.clear(); + candidates.addAll( inliers ); + } + + /* get the coefficient tiles of p2 */ + final ArrayList< Tile< ? extends M > > p2CoefficientsTiles = coefficientsTiles.get( p2 ); + + /* connect tiles across patches */ + for ( int i = 0; i < numCoefficients * numCoefficients; ++i ) + { + final Tile< ? > t1 = p1CoefficientsTiles.get( i ); + ra.setPosition( i, 0 ); + for ( int j = 0; j < numCoefficients * numCoefficients; ++j ) + { + ra.setPosition( j, 1 ); + final ArrayList< PointMatch > matches = ra.get(); + if ( matches.size() > 0 ) + { + final Tile< ? > t2 = p2CoefficientsTiles.get( j ); + t1.connect( t2, ra.get() ); + IJ.log( "Connected patch " + p1.getId() + ", coefficient " + i + " + patch " + p2.getId() + ", coefficient " + j + " by " + matches.size() + " samples." ); + } + } + } + } + + /* connect tiles within patch */ + for ( int y = 1; y < numCoefficients; ++y ) + { + final int yr = numCoefficients * y; + final int yr1 = yr - numCoefficients; + for ( int x = 0; x < numCoefficients; ++x ) + { + identityConnect( p1CoefficientsTiles.get( yr1 + x ), p1CoefficientsTiles.get( yr + x ), neighborWeight ); + } + } + for ( int y = 0; y < numCoefficients; ++y ) + { + final int yr = numCoefficients * y; + for ( int x = 1; x < numCoefficients; ++x ) + { + final int yrx = yr + x; + identityConnect( p1CoefficientsTiles.get( yrx ), p1CoefficientsTiles.get( yrx - 1 ), neighborWeight ); + } + } + } + + /* optimize */ + final TileConfiguration tc = new TileConfiguration(); + for ( final ArrayList< Tile< ? extends M > > coefficients : coefficientsTiles.values() ) + { + // for ( final Tile< ? > t : coefficients ) + // if ( t.getMatches().size() == 0 ) + // IJ.log( "bang" ); + tc.addTiles( coefficients ); + } + + try + { + tc.optimize( 0.01f, iterations, iterations, 0.75f ); + } + catch ( final NotEnoughDataPointsException e ) + { + // TODO Auto-generated catch block + e.printStackTrace(); + } + catch ( final IllDefinedDataPointsException e ) + { + // TODO Auto-generated catch block + e.printStackTrace(); + } + + /* save coefficients */ + final double[] ab = new double[ 2 ]; + final FSLoader loader = ( FSLoader ) layerset.getProject().getLoader(); + final String itsDir = loader.getUNUIdFolder() + "trakem2.its/"; + for ( final Entry< Patch, ArrayList< Tile< ? extends M > > > entry : coefficientsTiles.entrySet() ) + { + final FloatProcessor as = new FloatProcessor( numCoefficients, numCoefficients ); + final FloatProcessor bs = new FloatProcessor( numCoefficients, numCoefficients ); + + final Patch p = entry.getKey(); + + final double min = p.getMin(); + final double max = p.getMax(); + + final ArrayList< Tile< ? extends M > > tiles = entry.getValue(); + for ( int i = 0; i < numCoefficients * numCoefficients; ++i ) + { + final Tile< ? extends M > t = tiles.get( i ); + final Affine1D< ? > affine = t.getModel(); + affine.toArray( ab ); + + /* coefficients mapping into existing [min, max] */ + as.setf( i, ( float ) ab[ 0 ] ); + bs.setf( i, ( float ) ( ( max - min ) * ab[ 1 ] + min - ab[ 0 ] * min ) ); + } + final ImageStack coefficientsStack = new ImageStack( numCoefficients, numCoefficients ); + coefficientsStack.addSlice( as ); + coefficientsStack.addSlice( bs ); + + final String itsPath = itsDir + FSLoader.createIdPath( Long.toString( p.getId() ), "it", ".tif" ); + new File( itsPath ).getParentFile().mkdirs(); + IJ.saveAs( new ImagePlus( "", coefficientsStack ), "tif", itsPath ); + + } + } +} diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/PointMatchFilter.java b/TrakEM2_/src/main/java/org/janelia/intensity/PointMatchFilter.java new file mode 100644 index 00000000..f55ad2c5 --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/PointMatchFilter.java @@ -0,0 +1,20 @@ +/** + * + */ +package org.janelia.intensity; + +import java.util.Collection; +import java.util.List; + +import mpicbg.models.PointMatch; + +/** + * @author Stephan Saalfeld + * + */ +public interface PointMatchFilter +{ + public void filter( + final List< PointMatch > candidates, + final Collection< PointMatch > inliers ); +} diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionFilter.java b/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionFilter.java new file mode 100644 index 00000000..6572a8cc --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionFilter.java @@ -0,0 +1,48 @@ +/** + * + */ +package org.janelia.intensity; + +import java.util.Collection; +import java.util.List; + +import mpicbg.models.AffineModel1D; +import mpicbg.models.Model; +import mpicbg.models.PointMatch; + +/** + * @author Stephan Saalfeld + * + */ +public class RansacRegressionFilter implements PointMatchFilter +{ + final protected Model< ? > model = new AffineModel1D(); + final protected int iterations = 1000; + final protected float maxEpsilon = 0.1f; + final protected float minInlierRatio = 0.1f; + final protected int minNumInliers = 10; + final protected float maxTrust = 3.0f; + + @Override + public void filter( List< PointMatch > candidates, Collection< PointMatch > inliers ) + { + try + { + if ( + !model.filterRansac( + candidates, + inliers, + iterations, + maxEpsilon, + minInlierRatio, + minNumInliers, + maxTrust ) ) + inliers.clear(); + } + catch ( Exception e ) + { + inliers.clear(); + } + } + +} diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionReduceFilter.java b/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionReduceFilter.java new file mode 100644 index 00000000..21ec00b2 --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/RansacRegressionReduceFilter.java @@ -0,0 +1,83 @@ +/** + * + */ +package org.janelia.intensity; + +import java.util.Collection; +import java.util.Iterator; +import java.util.List; + +import mpicbg.models.AffineModel1D; +import mpicbg.models.Model; +import mpicbg.models.Point; +import mpicbg.models.PointMatch; + +/** + * @author Stephan Saalfeld + * + */ +public class RansacRegressionReduceFilter implements PointMatchFilter +{ + final protected Model< ? > model = new AffineModel1D(); + final protected int iterations = 1000; + final protected double maxEpsilon = 0.1; + final protected double minInlierRatio = 0.1; + final protected int minNumInliers = 10; + final protected double maxTrust = 3.0; + + final static protected double[] minMax( final Iterable< PointMatch > matches ) + { + final Iterator< PointMatch > iter = matches.iterator(); + PointMatch m = iter.next(); + double min = m.getP1().getL()[ 0 ], max = min; + while ( iter.hasNext() ) + { + m = iter.next(); + final double x = m.getP1().getL()[ 0 ]; + if ( x < min ) + min = x; + else if ( x > max ) + max = x; + } + return new double[]{ min, max }; + } + + @Override + public void filter( final List< PointMatch > candidates, final Collection< PointMatch > inliers ) + { + try + { + if ( + model.filterRansac( + candidates, + inliers, + iterations, + maxEpsilon, + minInlierRatio, + minNumInliers, + maxTrust ) ) + { + model.fit( inliers ); + + + final double[] minMax = minMax( inliers ); + + inliers.clear(); + + final Point p1 = new Point( new double[]{ minMax[ 0 ] } ); + final Point p2 = new Point( new double[]{ minMax[ 1 ] } ); + p1.apply( model ); + p2.apply( model ); + inliers.add( new PointMatch( p1, new Point( p1.getW().clone() ) ) ); + inliers.add( new PointMatch( p2, new Point( p2.getW().clone() ) ) ); + } + else + inliers.clear(); + } + catch ( final Exception e ) + { + inliers.clear(); + } + } + +} diff --git a/TrakEM2_/src/main/java/org/janelia/intensity/Render.java b/TrakEM2_/src/main/java/org/janelia/intensity/Render.java new file mode 100644 index 00000000..9958fae7 --- /dev/null +++ b/TrakEM2_/src/main/java/org/janelia/intensity/Render.java @@ -0,0 +1,302 @@ +/** + * License: GPL + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License 2 + * as published by the Free Software Foundation. + * + * 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, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ +package org.janelia.intensity; + +import ij.ImageJ; +import ij.ImagePlus; +import ij.ImageStack; +import ij.process.ByteProcessor; +import ij.process.ColorProcessor; +import ij.process.FloatProcessor; +import ij.process.ImageProcessor; +import ini.trakem2.Project; +import ini.trakem2.display.Layer; +import ini.trakem2.display.Patch; + +import java.awt.Rectangle; +import java.awt.image.BufferedImage; +import java.awt.image.WritableRaster; +import java.io.File; +import java.io.IOException; +import java.util.ArrayList; + +import javax.imageio.ImageIO; + +import mpicbg.ij.TransformMeshMapping; +import mpicbg.models.AffineModel2D; +import mpicbg.models.CoordinateTransform; +import mpicbg.models.CoordinateTransformList; +import mpicbg.models.CoordinateTransformMesh; +import mpicbg.models.NotEnoughDataPointsException; +import mpicbg.models.Point; +import mpicbg.models.PointMatch; +import mpicbg.models.SimilarityModel2D; +import mpicbg.models.TransformMesh; +import mpicbg.trakem2.transform.TransformMeshMappingWithMasks; +import mpicbg.trakem2.transform.TransformMeshMappingWithMasks.ImageProcessorWithMasks; +import mpicbg.trakem2.util.Downsampler; + +/** + * Render a patch. + * + * @author Stephan Saalfeld + */ +public class Render +{ + private Render() {} + + /** + * Create a {@link BufferedImage} from an existing pixel array. Make sure + * that pixels.length == width * height. + * + * @param pixels + * @param width + * @param height + * + * @return BufferedImage + */ + final static public BufferedImage createARGBImage( final int[] pixels, final int width, final int height ) + { + assert( pixels.length == width * height ) : "The number of pixels is not equal to width * height."; + + final BufferedImage image = new BufferedImage( width, height, BufferedImage.TYPE_INT_ARGB ); + final WritableRaster raster = image.getRaster(); + raster.setDataElements( 0, 0, width, height, pixels ); + return image; + } + + + final static void saveImage( final BufferedImage image, final String path, final String format ) throws IOException + { + ImageIO.write( image, format, new File( path ) ); + } + + /** + * Sample the average scaling of a given {@link CoordinateTransform} by transferring + * a set of point samples using the {@link CoordinateTransform} and then + * least-squares fitting a {@link SimilarityModel2D} to it. + * + * @param ct + * @param width of the samples set + * @param height of the samples set + * @param dx spacing between samples + * + * @return average scale factor + */ + final static protected double sampleAverageScale( final CoordinateTransform ct, final int width, final int height, final double dx ) + { + final ArrayList< PointMatch > samples = new ArrayList< PointMatch >(); + for ( double y = 0; y < height; y += dx ) + { + for ( double x = 0; x < width; x += dx ) + { + final Point p = new Point( new double[]{ x, y } ); + p.apply( ct ); + samples.add( new PointMatch( p, p ) ); + } + } + final SimilarityModel2D model = new SimilarityModel2D(); + try + { + model.fit( samples ); + } + catch ( final NotEnoughDataPointsException e ) + { + e.printStackTrace( System.err ); + return 1; + } + final double[] data = new double[ 6 ]; + model.toArray( data ); + return Math.sqrt( data[ 0 ] * data[ 0 ] + data[ 1 ] * data[ 1 ] ); + } + + + final static protected int bestMipmapLevel( final double scale ) + { + int invScale = ( int )( 1.0 / scale ); + int scaleLevel = 0; + while ( invScale > 1 ) + { + invScale >>= 1; + ++scaleLevel; + } + return scaleLevel; + } + + /** + * Create an affine transformation that compensates for both scale and + * pixel shift of a mipmap level that was generated by top-left pixel + * averaging. + * + * @param scaleLevel + * @return + */ + final static protected AffineModel2D createScaleLevelTransform( final int scaleLevel ) + { + final AffineModel2D a = new AffineModel2D(); + final int scale = 1 << scaleLevel; + final double t = ( scale - 1 ) * 0.5; + a.set( scale, 0, 0, scale, t, t ); + return a; + } + + /** + * Renders a patch, mapping its intensities [min, max] → [0, 1] + * + * @param patch the patch to be rendered + * @param targetImage target pixels, specifies the target box + * @param targetWeight target weight pixels, depending on alpha + * @param x target box offset in world coordinates + * @param y target box offset in world coordinates + * @param scale target scale + */ + final static public void render( + final Patch patch, + final int coefficientsWidth, + final int coefficientsHeight, + final FloatProcessor targetImage, + final FloatProcessor targetWeight, + final ColorProcessor targetCoefficients, + final double x, + final double y, + final double scale ) + { + /* assemble coordinate transformations and add bounding box offset */ + final CoordinateTransformList< CoordinateTransform > ctl = new CoordinateTransformList< CoordinateTransform >(); + ctl.add( patch.getFullCoordinateTransform() ); + final AffineModel2D affineScale = new AffineModel2D(); + affineScale.set( scale, 0, 0, scale, -x * scale, -y * scale ); + ctl.add( affineScale ); + + /* estimate average scale and generate downsampled source */ + final int width = patch.getOWidth(), height = patch.getOHeight(); + final double s = sampleAverageScale( ctl, width, height, width / patch.getMeshResolution() ); + final int mipmapLevel = bestMipmapLevel( s ); + final ImageProcessor ipMipmap = Downsampler.downsampleImageProcessor( patch.getImageProcessor(), mipmapLevel ); + + /* create a target */ + final ImageProcessor tp = ipMipmap.createProcessor( targetImage.getWidth(), targetImage.getHeight() ); + + /* prepare and downsample alpha mask if there is one */ + final ByteProcessor bpMaskMipmap; + final ByteProcessor bpMaskTarget; + + final ByteProcessor bpMask = patch.getAlphaMask(); + + if ( bpMask == null ) + { + bpMaskMipmap = null; + bpMaskTarget = null; + } + else + { + bpMaskMipmap = bpMask == null ? null : Downsampler.downsampleByteProcessor( bpMask, mipmapLevel ); + bpMaskTarget = new ByteProcessor( tp.getWidth(), tp.getHeight() ); + } + + /* create coefficients map */ + final ColorProcessor cp = new ColorProcessor( ipMipmap.getWidth(), ipMipmap.getHeight() ); + final int w = cp.getWidth(); + final int h = cp.getHeight(); + for ( int yi = 0; yi < h; ++yi ) + { + final int yc = yi * coefficientsHeight / h; + final int ic = yc * coefficientsWidth; + final int iyi = yi * w; + for ( int xi = 0; xi < w; ++xi ) + cp.set( iyi + xi, ic + ( xi * coefficientsWidth / w ) + 1 ); + } + + /* attach mipmap transformation */ + final CoordinateTransformList< CoordinateTransform > ctlMipmap = new CoordinateTransformList< CoordinateTransform >(); + ctlMipmap.add( createScaleLevelTransform( mipmapLevel ) ); + ctlMipmap.add( ctl ); + + /* create mesh */ + final CoordinateTransformMesh mesh = new CoordinateTransformMesh( ctlMipmap, patch.getMeshResolution(), ipMipmap.getWidth(), ipMipmap.getHeight() ); + + /* render */ + final ImageProcessorWithMasks source = new ImageProcessorWithMasks( ipMipmap, bpMaskMipmap, null ); + final ImageProcessorWithMasks target = new ImageProcessorWithMasks( tp, bpMaskTarget, null ); + final TransformMeshMappingWithMasks< TransformMesh > mapping = new TransformMeshMappingWithMasks< TransformMesh >( mesh ); + mapping.mapInterpolated( source, target ); + + final TransformMeshMapping< TransformMesh > coefficientsMapMapping = new TransformMeshMapping< TransformMesh >( mesh ); + coefficientsMapMapping.map( cp, targetCoefficients ); + + /* set alpha channel */ + final byte[] alphaPixels; + + if ( bpMaskTarget != null ) + alphaPixels = ( byte[] )bpMaskTarget.getPixels(); + else + alphaPixels = ( byte[] )target.outside.getPixels(); + + /* convert */ + final double min = patch.getMin(); + final double max = patch.getMax(); + final double a = 1.0 / ( max - min ); + final double b = 1.0 / 255.0; + + for ( int i = 0; i < alphaPixels.length; ++i ) + targetImage.setf( i, ( float )( ( tp.getf( i ) - min ) * a ) ); + + for ( int i = 0; i < alphaPixels.length; ++i ) + targetWeight.setf( i, ( float )( ( alphaPixels[ i ] & 0xff ) * b ) ); + } + + final static public void main( final String... args ) + { + new ImageJ(); + + final double scale = 0.05; +// final double scale = 1; + final int numCoefficients = 4; + + final Project project = Project.openFSProject( "/home/saalfeld/tmp/bock-lens-correction/subproject.xml", false ); + final Layer layer = project.getRootLayerSet().getLayer( 0 ); + final ArrayList< Patch > patches = ( ArrayList )layer.getDisplayables( Patch.class ); + + final Patch patch1 = patches.get( 0 ); + final Patch patch2 = patches.get( 2 ); + final Rectangle box = patch1.getBoundingBox().intersection( patch2.getBoundingBox() ); + //final Rectangle box = patch1.getBoundingBox(); + + final int w = ( int )( box.width * scale + 0.5 ); + final int h = ( int )( box.height * scale + 0.5 ); + + final FloatProcessor pixels1 = new FloatProcessor( w, h ); + final FloatProcessor weights1 = new FloatProcessor( w, h ); + final ColorProcessor coefficients1 = new ColorProcessor( w, h ); + final FloatProcessor pixels2 = new FloatProcessor( w, h ); + final FloatProcessor weights2 = new FloatProcessor( w, h ); + final ColorProcessor coefficients2 = new ColorProcessor( w, h ); + + render( patch1, numCoefficients, numCoefficients, pixels1, weights1, coefficients1, box.x, box.y, scale ); + render( patch2, numCoefficients, numCoefficients, pixels2, weights2, coefficients2, box.x, box.y, scale ); + + final ImageStack stack = new ImageStack( w, h ); + stack.addSlice( pixels1 ); + stack.addSlice( pixels2 ); + stack.addSlice( weights1 ); + stack.addSlice( weights2 ); + stack.addSlice( coefficients1.convertToFloatProcessor() ); + stack.addSlice( coefficients2.convertToFloatProcessor() ); + + new ImagePlus( "", stack ).show(); + } +} diff --git a/pom.xml b/pom.xml index 4d023ed6..79663369 100644 --- a/pom.xml +++ b/pom.xml @@ -9,7 +9,7 @@ sc.fiji pom-fiji - 8.2.0 + 15.0.0 pom-trakem2 -- 2.11.4.GIT