TrakEM2_/src/main/java/org/janelia/intensity/MatchIntensities.java

   1 /**
   2  *
   3  */
   4 package org.janelia.intensity;
   5
   6 import ij.IJ;
   7 import ij.ImagePlus;
   8 import ij.ImageStack;
   9 import ij.gui.GenericDialog;
  10 import ij.gui.Roi;
  11 import ij.measure.Calibration;
  12 import ij.process.ColorProcessor;
  13 import ij.process.FloatProcessor;
  14 import ini.trakem2.Project;
  15 import ini.trakem2.display.Display;
  16 import ini.trakem2.display.Displayable;
  17 import ini.trakem2.display.Layer;
  18 import ini.trakem2.display.LayerSet;
  19 import ini.trakem2.display.Patch;
  20 import ini.trakem2.persistence.FSLoader;
  21 import ini.trakem2.plugin.TPlugIn;
  22 import ini.trakem2.utils.Utils;
  23
  24 import java.awt.Rectangle;
  25 import java.io.File;
  26 import java.util.ArrayList;
  27 import java.util.Collection;
  28 import java.util.HashMap;
  29 import java.util.HashSet;
  30 import java.util.List;
  31 import java.util.Map.Entry;
  32 import java.util.concurrent.ExecutionException;
  33
  34 import mpicbg.models.Affine1D;
  35 import mpicbg.models.AffineModel1D;
  36 import mpicbg.models.IdentityModel;
  37 import mpicbg.models.IllDefinedDataPointsException;
  38 import mpicbg.models.InterpolatedAffineModel1D;
  39 import mpicbg.models.Model;
  40 import mpicbg.models.NotEnoughDataPointsException;
  41 import mpicbg.models.Point;
  42 import mpicbg.models.PointMatch;
  43 import mpicbg.models.Tile;
  44 import mpicbg.models.TileConfiguration;
  45 import mpicbg.models.TranslationModel1D;
  46 import net.imglib2.img.list.ListImg;
  47 import net.imglib2.img.list.ListRandomAccess;
  48
  49 /**
  50  *
  51  * @author Stephan Saalfeld <saalfelds@janelia.hhmi.org>
  52  * @author Philipp Hanslovsky
  53  */
  54 public class MatchIntensities implements TPlugIn
  55 {
  56         protected LayerSet layerset = null;
  57
  58         static protected int numCoefficients = 8;
  59         static protected double lambda1 = 0.01;
  60         static protected double lambda2 = 0.01;
  61         static protected double neighborWeight = 0.1;
  62         static protected int radius = 5;
  63         static protected int iterations = 2000;
  64         static protected double scale = -1;
  65
  66         private Layer currentLayer( final Object... params )
  67         {
  68                 final Layer layer;
  69                 if ( params != null && params[ 0 ] != null )
  70                 {
  71                         final Object param = params[ 0 ];
  72                         if ( Layer.class.isInstance( param ) )
  73                                 layer = ( Layer ) param;
  74                         else if ( LayerSet.class.isInstance( param ) )
  75                                 layer = ( ( LayerSet ) param ).getLayer( 0 );
  76                         else if ( Displayable.class.isInstance( param ) )
  77                                 layer = ( ( Displayable ) param ).getLayer();
  78                         else
  79                                 layer = null;
  80                 }
  81                 else
  82                 {
  83                         final Display front = Display.getFront();
  84                         if ( front == null )
  85                                 layer = Project.getProjects().get( 0 ).getRootLayerSet().getLayer( 0 );
  86                         else
  87                                 layer = front.getLayer();
  88                 }
  89                 return layer;
  90         }
  91
  92         private static Rectangle getRoi( final LayerSet layerset )
  93         {
  94                 final Roi roi;
  95                 final Display front = Display.getFront();
  96                 if ( front == null )
  97                         roi = null;
  98                 else
  99                         roi = front.getRoi();
 100                 if ( roi == null )
 101                         return new Rectangle( 0, 0, ( int ) layerset.getLayerWidth(), ( int ) layerset.getLayerHeight() );
 102                 else
 103                         return roi.getBounds();
 104         }
 105
 106         private double suggestScale( final List< Layer > layers )
 107         {
 108                 if ( layers.size() < 2 )
 109                         return 0.0;
 110
 111                 final Layer layer1 = layers.get( 0 );
 112                 final Layer layer2 = layers.get( 1 );
 113                 final Calibration calib = layer1.getParent().getCalibration();
 114                 final double width = ( calib.pixelWidth + calib.pixelHeight ) * 0.5;
 115                 final double depth = calib.pixelDepth;
 116
 117                 return ( layer2.getZ() - layer1.getZ() ) * width / depth;
 118         }
 119
 120         @Override
 121         public boolean setup( final Object... params )
 122         {
 123                 if ( params != null && params[ 0 ] != null )
 124                 {
 125                         final Object param = params[ 0 ];
 126                         if ( LayerSet.class.isInstance( param ) )
 127                                 layerset = ( LayerSet ) param;
 128                         else if ( Displayable.class.isInstance( param ) )
 129                                 layerset = ( ( Displayable ) param ).getLayerSet();
 130                         else
 131                                 return false;
 132                 }
 133                 else
 134                 {
 135                         final Display front = Display.getFront();
 136                         if ( front == null )
 137                                 layerset = Project.getProjects().get( 0 ).getRootLayerSet();
 138                         else
 139                                 layerset = front.getLayerSet();
 140                 }
 141                 return true;
 142         }
 143
 144         final static protected < T extends Model< T > & Affine1D< T > > HashMap< Patch, ArrayList< Tile< T > > > generateCoefficientsTiles(
 145                         final Collection< Patch > patches,
 146                         final T template,
 147                         final int nCoefficients )
 148         {
 149                 final HashMap< Patch, ArrayList< Tile< T > > > map = new HashMap< Patch, ArrayList< Tile< T > > >();
 150                 for ( final Patch p : patches )
 151                 {
 152                         final ArrayList< Tile< T > > coefficientModels = new ArrayList< Tile< T > >();
 153                         for ( int i = 0; i < nCoefficients; ++i )
 154                                 coefficientModels.add( new Tile< T >( template.copy() ) );
 155
 156                         map.put( p, coefficientModels );
 157                 }
 158                 return map;
 159         }
 160
 161         final static protected void identityConnect( final Tile< ? > t1, final Tile< ? > t2, final double weight )
 162         {
 163                 final ArrayList< PointMatch > matches = new ArrayList< PointMatch >();
 164                 matches.add( new PointMatch( new Point( new double[] { 0 } ), new Point( new double[] { 0 } ) ) );
 165                 matches.add( new PointMatch( new Point( new double[] { 1 } ), new Point( new double[] { 1 } ) ) );
 166                 t1.connect( t2, matches );
 167         }
 168
 169         @Override
 170         public Object invoke( final Object... params )
 171         {
 172                 if ( !setup( params ) )
 173                         return null;
 174
 175                 final Layer layer = currentLayer( params );
 176                 final GenericDialog gd = new GenericDialog( "Match intensities" );
 177                 Utils.addLayerRangeChoices( layer, gd );
 178                 gd.addMessage( "Layer neighborhood range :" );
 179                 gd.addNumericField( "test_maximally :", radius, 0, 6, "layers" );
 180                 gd.addMessage( "Optimizer :" );
 181                 gd.addNumericField( "iterations :", iterations, 0, 6, "" );
 182                 gd.addNumericField( "scale_regularization :", lambda1, 2, 6, "" );
 183                 gd.addNumericField( "translation_regularization :", lambda2, 2, 6, "" );
 184                 gd.addNumericField( "smoothness_regularization :", neighborWeight, 2, 6, "" );
 185                 gd.showDialog();
 186
 187                 if ( gd.wasCanceled() )
 188                         return null;
 189
 190                 final List< Layer > layers = layerset.getLayers().subList( gd.getNextChoiceIndex(), gd.getNextChoiceIndex() + 1 );
 191                 radius = ( int ) gd.getNextNumber();
 192                 iterations = ( int )gd.getNextNumber();
 193                 lambda1 = gd.getNextNumber();
 194                 lambda2 = gd.getNextNumber();
 195                 neighborWeight = gd.getNextNumber();
 196
 197                 try
 198                 {
 199                         run( layers, radius, scale, numCoefficients, lambda1, lambda2, neighborWeight, getRoi( layerset ) );
 200                 }
 201                 catch ( final InterruptedException e )
 202                 {
 203                         Utils.log( "Match intensities interrupted." );
 204                 }
 205                 catch ( final ExecutionException e )
 206                 {
 207                         Utils.log( "Match intenities ExecutiuonException occurred:" );
 208                         e.printStackTrace( System.out );
 209                 }
 210
 211                 return null;
 212         }
 213
 214     @Override
 215     public boolean applies( final Object ob )
 216     {
 217         return true;
 218     }
 219
 220
 221     /**
 222      * TODO Test!  And then desperately multi-thread coefficient collection.
 223      *
 224      * @param layers
 225      * @param radius
 226      * @param scale
 227      * @param numCoefficients
 228      * @param lambda1
 229      * @param lambda2
 230      * @param neighborWeight
 231      * @param roi
 232      */
 233         public < M extends Model< M > & Affine1D< M > > void run(
 234                         final List< Layer > layers,
 235                         final int radius,
 236                         final double scale,
 237                         final int numCoefficients,
 238                         final double lambda1,
 239                         final double lambda2,
 240                         final double neighborWeight,
 241                         final Rectangle roi ) throws InterruptedException, ExecutionException
 242         {
 243                 // final PointMatchFilter filter = new RansacRegressionFilter();
 244                 final PointMatchFilter filter = new RansacRegressionReduceFilter();
 245
 246                 /* collect patches */
 247                 final ArrayList< Patch > patches = new ArrayList< Patch >();
 248                 for ( final Layer layer : layers )
 249                         patches.addAll( ( ArrayList )layer.getDisplayables( Patch.class, roi ) );
 250
 251                 /* generate coefficient tiles for all patches
 252                  * TODO consider offering alternative models */
 253                 final HashMap< Patch, ArrayList< Tile< ? extends M > > > coefficientsTiles =
 254                                 ( HashMap ) generateCoefficientsTiles(
 255                                                 patches,
 256                                                 new InterpolatedAffineModel1D< InterpolatedAffineModel1D< AffineModel1D, TranslationModel1D >, IdentityModel >(
 257                                                                 new InterpolatedAffineModel1D< AffineModel1D, TranslationModel1D >(
 258                                                                                 new AffineModel1D(), new TranslationModel1D(), lambda1 ),
 259                                                                 new IdentityModel(), lambda2 ),
 260                                                 numCoefficients * numCoefficients );
 261
 262                 /* completed patches */
 263                 final HashSet< Patch > completedPatches = new HashSet< Patch >();
 264
 265                 for ( final Patch p1 : patches )
 266                 {
 267                         completedPatches.add( p1 );
 268
 269                         final Rectangle box1 = p1.getBoundingBox().intersection( roi );
 270
 271                         final ArrayList< Patch > p2s = new ArrayList< Patch >();
 272
 273                         /* across adjacent layers */
 274                         final int layerIndex = layerset.getLayerIndex( p1.getLayer().getId() );
 275                         for ( int i = layerIndex - radius; i <= layerIndex + radius; ++i )
 276                         {
 277                                 final Layer layer = layerset.getLayer( i );
 278                                 if ( layer != null )
 279                                         p2s.addAll( ( Collection ) layer.getDisplayables( Patch.class, box1 ) );
 280                         }
 281
 282                         /* get the coefficient tiles */
 283                         final ArrayList< Tile< ? extends M > > p1CoefficientsTiles = coefficientsTiles.get( p1 );
 284
 285                         for ( final Patch p2 : p2s )
 286                         {
 287                                 /*
 288                                  * if this patch had been processed earlier, all matches are
 289                                  * already in
 290                                  */
 291                                 if ( completedPatches.contains( p2 ) )
 292                                         continue;
 293
 294                                 /* render intersection */
 295                                 final Rectangle box2 = p2.getBoundingBox();
 296                                 final Rectangle box = box1.intersection( box2 );
 297
 298                                 final int w = ( int ) ( box.width * scale + 0.5 );
 299                                 final int h = ( int ) ( box.height * scale + 0.5 );
 300                                 final int n = w * h;
 301
 302                                 final FloatProcessor pixels1 = new FloatProcessor( w, h );
 303                                 final FloatProcessor weights1 = new FloatProcessor( w, h );
 304                                 final ColorProcessor coefficients1 = new ColorProcessor( w, h );
 305                                 final FloatProcessor pixels2 = new FloatProcessor( w, h );
 306                                 final FloatProcessor weights2 = new FloatProcessor( w, h );
 307                                 final ColorProcessor coefficients2 = new ColorProcessor( w, h );
 308
 309                                 Render.render( p1, numCoefficients, numCoefficients, pixels1, weights1, coefficients1, box.x, box.y, scale );
 310                                 Render.render( p2, numCoefficients, numCoefficients, pixels2, weights2, coefficients2, box.x, box.y, scale );
 311
 312                                 /*
 313                                  * generate a matrix of all coefficients in p1 to all
 314                                  * coefficients in p2 to store matches
 315                                  */
 316                                 final ArrayList< ArrayList< PointMatch > > list = new ArrayList< ArrayList< PointMatch > >();
 317                                 for ( int i = 0; i < numCoefficients * numCoefficients * numCoefficients * numCoefficients; ++i )
 318                                         list.add( new ArrayList< PointMatch >() );
 319                                 final ListImg< ArrayList< PointMatch > > matrix = new ListImg< ArrayList< PointMatch > >( list, numCoefficients * numCoefficients, numCoefficients * numCoefficients );
 320                                 final ListRandomAccess< ArrayList< PointMatch > > ra = matrix.randomAccess();
 321
 322                                 /*
 323                                  * iterate over all pixels and feed matches into the match
 324                                  * matrix
 325                                  */
 326                                 for ( int i = 0; i < n; ++i )
 327                                 {
 328                                         final int c1 = coefficients1.get( i );
 329                                         if ( c1 > 0 )
 330                                         {
 331                                                 final int c2 = coefficients2.get( i );
 332                                                 if ( c2 > 0 )
 333                                                 {
 334                                                         final double w1 = weights1.getf( i );
 335                                                         if ( w1 > 0 )
 336                                                         {
 337                                                                 final double w2 = weights2.getf( i );
 338                                                                 if ( w2 > 0 )
 339                                                                 {
 340                                                                         final double p = pixels1.getf( i );
 341                                                                         final double q = pixels2.getf( i );
 342                                                                         final PointMatch pq = new PointMatch( new Point( new double[] { p } ), new Point( new double[] { q } ), w1 * w2 );
 343
 344                                                                         /* first label is 1 */
 345                                                                         ra.setPosition( c1 - 1, 0 );
 346                                                                         ra.setPosition( c2 - 1, 1 );
 347                                                                         ra.get().add( pq );
 348                                                                 }
 349                                                         }
 350                                                 }
 351                                         }
 352                                 }
 353
 354                                 /* filter matches */
 355                                 final ArrayList< PointMatch > inliers = new ArrayList< PointMatch >();
 356                                 for ( final ArrayList< PointMatch > candidates : matrix )
 357                                 {
 358                                         inliers.clear();
 359                                         filter.filter( candidates, inliers );
 360                                         candidates.clear();
 361                                         candidates.addAll( inliers );
 362                                 }
 363
 364                                 /* get the coefficient tiles of p2 */
 365                                 final ArrayList< Tile< ? extends M > > p2CoefficientsTiles = coefficientsTiles.get( p2 );
 366
 367                                 /* connect tiles across patches */
 368                                 for ( int i = 0; i < numCoefficients * numCoefficients; ++i )
 369                                 {
 370                                         final Tile< ? > t1 = p1CoefficientsTiles.get( i );
 371                                         ra.setPosition( i, 0 );
 372                                         for ( int j = 0; j < numCoefficients * numCoefficients; ++j )
 373                                         {
 374                                                 ra.setPosition( j, 1 );
 375                                                 final ArrayList< PointMatch > matches = ra.get();
 376                                                 if ( matches.size() > 0 )
 377                                                 {
 378                                                         final Tile< ? > t2 = p2CoefficientsTiles.get( j );
 379                                                         t1.connect( t2, ra.get() );
 380                                                         IJ.log( "Connected patch " + p1.getId() + ", coefficient " + i + "  +  patch " + p2.getId() + ", coefficient " + j + " by " + matches.size() + " samples." );
 381                                                 }
 382                                         }
 383                                 }
 384                         }
 385
 386                         /* connect tiles within patch */
 387                         for ( int y = 1; y < numCoefficients; ++y )
 388                         {
 389                                 final int yr = numCoefficients * y;
 390                                 final int yr1 = yr - numCoefficients;
 391                                 for ( int x = 0; x < numCoefficients; ++x )
 392                                 {
 393                                         identityConnect( p1CoefficientsTiles.get( yr1 + x ), p1CoefficientsTiles.get( yr + x ), neighborWeight );
 394                                 }
 395                         }
 396                         for ( int y = 0; y < numCoefficients; ++y )
 397                         {
 398                                 final int yr = numCoefficients * y;
 399                                 for ( int x = 1; x < numCoefficients; ++x )
 400                                 {
 401                                         final int yrx = yr + x;
 402                                         identityConnect( p1CoefficientsTiles.get( yrx ), p1CoefficientsTiles.get( yrx - 1 ), neighborWeight );
 403                                 }
 404                         }
 405                 }
 406
 407                 /* optimize */
 408                 final TileConfiguration tc = new TileConfiguration();
 409                 for ( final ArrayList< Tile< ? extends M > > coefficients : coefficientsTiles.values() )
 410                 {
 411                         // for ( final Tile< ? > t : coefficients )
 412                         // if ( t.getMatches().size() == 0 )
 413                         // IJ.log( "bang" );
 414                         tc.addTiles( coefficients );
 415                 }
 416
 417                 try
 418                 {
 419                         tc.optimize( 0.01f, iterations, iterations, 0.75f );
 420                 }
 421                 catch ( final NotEnoughDataPointsException e )
 422                 {
 423                         // TODO Auto-generated catch block
 424                         e.printStackTrace();
 425                 }
 426                 catch ( final IllDefinedDataPointsException e )
 427                 {
 428                         // TODO Auto-generated catch block
 429                         e.printStackTrace();
 430                 }
 431
 432                 /* save coefficients */
 433                 final double[] ab = new double[ 2 ];
 434                 final FSLoader loader = ( FSLoader ) layerset.getProject().getLoader();
 435                 final String itsDir = loader.getUNUIdFolder() + "trakem2.its/";
 436                 for ( final Entry< Patch, ArrayList< Tile< ? extends M > > > entry : coefficientsTiles.entrySet() )
 437                 {
 438                         final FloatProcessor as = new FloatProcessor( numCoefficients, numCoefficients );
 439                         final FloatProcessor bs = new FloatProcessor( numCoefficients, numCoefficients );
 440
 441                         final Patch p = entry.getKey();
 442
 443                         final double min = p.getMin();
 444                         final double max = p.getMax();
 445
 446                         final ArrayList< Tile< ? extends M > > tiles = entry.getValue();
 447                         for ( int i = 0; i < numCoefficients * numCoefficients; ++i )
 448                         {
 449                                 final Tile< ? extends M > t = tiles.get( i );
 450                                 final Affine1D< ? > affine = t.getModel();
 451                                 affine.toArray( ab );
 452
 453                                 /* coefficients mapping into existing [min, max] */
 454                                 as.setf( i, ( float ) ab[ 0 ] );
 455                                 bs.setf( i, ( float ) ( ( max - min ) * ab[ 1 ] + min - ab[ 0 ] * min ) );
 456                         }
 457                         final ImageStack coefficientsStack = new ImageStack( numCoefficients, numCoefficients );
 458                         coefficientsStack.addSlice( as );
 459                         coefficientsStack.addSlice( bs );
 460
 461                         final String itsPath = itsDir + FSLoader.createIdPath( Long.toString( p.getId() ), "it", ".tif" );
 462                         new File( itsPath ).getParentFile().mkdirs();
 463                         IJ.saveAs( new ImagePlus( "", coefficientsStack ), "tif", itsPath );
 464
 465                 }
 466         }
 467 }