mpicbg/trakem2/align/Align.java

   1 /**
   2  *
   3  */
   4 package mpicbg.trakem2.align;
   5
   6 import java.awt.Color;
   7 import java.awt.Rectangle;
   8 import java.awt.geom.AffineTransform;
   9 import java.io.Serializable;
  10 import java.util.ArrayList;
  11 import java.util.Collection;
  12 import java.util.List;
  13 import java.util.Set;
  14 import java.util.concurrent.atomic.AtomicInteger;
  15
  16 import ij.IJ;
  17 import ij.ImagePlus;
  18 import ij.gui.GenericDialog;
  19
  20 import ini.trakem2.display.Display;
  21 import ini.trakem2.display.Displayable;
  22 import ini.trakem2.display.Layer;
  23 import ini.trakem2.display.Patch;
  24 import ini.trakem2.display.Selection;
  25 import ini.trakem2.persistence.Loader;
  26 import ini.trakem2.persistence.FSLoader;
  27 import ini.trakem2.utils.Utils;
  28
  29 import mpicbg.ij.FeatureTransform;
  30 import mpicbg.ij.SIFT;
  31 import mpicbg.imagefeatures.Feature;
  32 import mpicbg.imagefeatures.FloatArray2DSIFT;
  33 import mpicbg.models.AbstractAffineModel2D;
  34 import mpicbg.models.AffineModel2D;
  35 import mpicbg.models.NotEnoughDataPointsException;
  36 import mpicbg.models.PointMatch;
  37 import mpicbg.models.SimilarityModel2D;
  38 import mpicbg.models.Tile;
  39 import mpicbg.trakem2.transform.MovingLeastSquaresTransform;
  40 import mpicbg.trakem2.transform.RigidModel2D;
  41 import mpicbg.trakem2.transform.TranslationModel2D;
  42
  43 /**
  44  * A collection of methods regarding SIFT-based alignment
  45  *
  46  * TODO Bring the methods and tasks into a class for each method and clean up this mess.
  47  *
  48  * @author Stephan Saalfeld <saalfeld@mpi-cbg.de>
  49  */
  50 public class Align
  51 {
  52         static public class Param implements Serializable
  53         {
  54                 private static final long serialVersionUID = -2247163691721712461L;
  55
  56                 final public FloatArray2DSIFT.Param sift = new FloatArray2DSIFT.Param();
  57
  58                 /**
  59                  * Closest/next closest neighbour distance ratio
  60                  */
  61                 public float rod = 0.92f;
  62
  63                 /**
  64                  * Maximal allowed alignment error in px
  65                  */
  66                 public float maxEpsilon = 100.0f;
  67
  68                 /**
  69                  * Inlier/candidates ratio
  70                  */
  71                 public float minInlierRatio = 0.2f;
  72
  73                 /**
  74                  * Implemeted transformation models for choice
  75                  */
  76                 final static public String[] modelStrings = new String[]{ "Translation", "Rigid", "Similarity", "Affine" };
  77                 public int expectedModelIndex = 1;
  78                 public int desiredModelIndex = 1;
  79
  80                 public float correspondenceWeight = 1;
  81
  82                 public Param()
  83                 {
  84                         sift.maxOctaveSize = 600;
  85                         sift.fdSize = 8;
  86                 }
  87
  88                 public void addFields( final GenericDialog gd )
  89                 {
  90                         SIFT.addFields( gd, sift );
  91
  92                         gd.addNumericField( "closest/next_closest_ratio :", rod, 2 );
  93
  94                         gd.addMessage( "Geometric Consensus Filter:" );
  95                         gd.addNumericField( "maximal_alignment_error :", maxEpsilon, 2, 6, "px" );
  96                         gd.addNumericField( "inlier_ratio :", minInlierRatio, 2 );
  97                         gd.addChoice( "expected_transformation :", modelStrings, modelStrings[ expectedModelIndex ] );
  98
  99                         gd.addMessage( "Alignment:" );
 100                         gd.addChoice( "desired_transformation :", modelStrings, modelStrings[ desiredModelIndex ] );
 101                         gd.addNumericField( "correspondence weight :", correspondenceWeight, 2 );
 102                 }
 103
 104                 public boolean readFields( final GenericDialog gd )
 105                 {
 106                         SIFT.readFields( gd, sift );
 107
 108                         rod = ( float )gd.getNextNumber();
 109
 110                         maxEpsilon = ( float )gd.getNextNumber();
 111                         minInlierRatio = ( float )gd.getNextNumber();
 112                         expectedModelIndex = gd.getNextChoiceIndex();
 113                         desiredModelIndex = gd.getNextChoiceIndex();
 114
 115                         correspondenceWeight = ( float )gd.getNextNumber();
 116
 117                         return !gd.invalidNumber();
 118                 }
 119
 120                 public boolean setup( final String title )
 121                 {
 122                         final GenericDialog gd = new GenericDialog( title );
 123
 124                         addFields( gd );
 125
 126                         do
 127                         {
 128                                 gd.showDialog();
 129                                 if ( gd.wasCanceled() ) return false;
 130                         }
 131                         while ( !readFields( gd ) );
 132
 133                         return true;
 134                 }
 135
 136                 public Param clone()
 137                 {
 138                         Param p = new Param();
 139
 140                         p.sift.initialSigma = this.sift.initialSigma;
 141                         p.sift.steps = this.sift.steps;
 142                         p.sift.minOctaveSize = this.sift.minOctaveSize;
 143                         p.sift.maxOctaveSize = this.sift.maxOctaveSize;
 144                         p.sift.fdSize = this.sift.fdSize;
 145                         p.sift.fdBins = this.sift.fdBins;
 146
 147                         p.rod = rod;
 148                         p.maxEpsilon = maxEpsilon;
 149                         p.minInlierRatio = minInlierRatio;
 150                         p.expectedModelIndex = expectedModelIndex;
 151                         p.desiredModelIndex = desiredModelIndex;
 152
 153                         p.correspondenceWeight = correspondenceWeight;
 154
 155                         return p;
 156                 }
 157
 158                 /**
 159                  * Check if two parameter sets are equal.  So far, this method ignores
 160                  * the parameter {@link #desiredModelIndex} which defines the
 161                  * transformation class to be used for {@link Tile} alignment.  This
 162                  * makes sense for the current use in {@link PointMatch} serialization
 163                  * but might be misleading for other applications.
 164                  *
 165                  * TODO Think about this.
 166                  *
 167                  * @param p
 168                  * @return
 169                  */
 170                 public boolean equals( Param p )
 171                 {
 172                         return
 173                                 sift.equals( p.sift ) &&
 174                                 ( rod == p.rod ) &&
 175                                 ( maxEpsilon == p.maxEpsilon ) &&
 176                                 ( minInlierRatio == p.minInlierRatio ) &&
 177                                 ( expectedModelIndex == p.expectedModelIndex );
 178 //                      && ( desiredModelIndex == p.desiredModelIndex );
 179                 }
 180         }
 181
 182         final static public Param param = new Param();
 183
 184         static public class ParamOptimize extends Param
 185         {
 186                 private static final long serialVersionUID = 970673723211054580L;
 187
 188                 /**
 189                  * Maximal number of iteration allowed for the optimizer.
 190                  */
 191                 public int maxIterations = 2000;
 192
 193                 /**
 194                  * Maximal number of iterations allowed to not change the parameter to
 195                  * be optimized.
 196                  */
 197                 public int maxPlateauwidth = 200;
 198
 199                 @Override
 200                 public void addFields( final GenericDialog gd )
 201                 {
 202                         super.addFields( gd );
 203
 204                         gd.addNumericField( "maximal_iterations :", maxIterations, 0 );
 205                         gd.addNumericField( "maximal_plateauwidth :", maxPlateauwidth, 0 );
 206                 }
 207
 208                 @Override
 209                 public boolean readFields( final GenericDialog gd )
 210                 {
 211                         super.readFields( gd );
 212
 213                         maxIterations = ( int )gd.getNextNumber();
 214                         maxPlateauwidth = ( int )gd.getNextNumber();
 215
 216                         return !gd.invalidNumber();
 217                 }
 218
 219                 @Override
 220                 final public boolean setup( final String title )
 221                 {
 222                         final GenericDialog gd = new GenericDialog( title );
 223
 224                         addFields( gd );
 225
 226                         do
 227                         {
 228                                 gd.showDialog();
 229                                 if ( gd.wasCanceled() ) return false;
 230                         }
 231                         while ( !readFields( gd ) );
 232
 233                         return true;
 234                 }
 235
 236                 @Override
 237                 final public ParamOptimize clone()
 238                 {
 239                         ParamOptimize p = new ParamOptimize();
 240
 241                         p.sift.initialSigma = this.sift.initialSigma;
 242                         p.sift.steps = this.sift.steps;
 243                         p.sift.minOctaveSize = this.sift.minOctaveSize;
 244                         p.sift.maxOctaveSize = this.sift.maxOctaveSize;
 245                         p.sift.fdSize = this.sift.fdSize;
 246                         p.sift.fdBins = this.sift.fdBins;
 247
 248                         p.rod = rod;
 249                         p.maxEpsilon = maxEpsilon;
 250                         p.minInlierRatio = minInlierRatio;
 251                         p.expectedModelIndex = expectedModelIndex;
 252
 253                         p.desiredModelIndex = desiredModelIndex;
 254                         p.maxIterations = maxIterations;
 255                         p.maxPlateauwidth = maxPlateauwidth;
 256
 257                         return p;
 258                 }
 259
 260                 public boolean equals( ParamOptimize p )
 261                 {
 262                         return
 263                                 super.equals( p ) &&
 264                                 ( maxIterations == p.maxIterations ) &&
 265                                 ( maxPlateauwidth == p.maxPlateauwidth );
 266                 }
 267         }
 268
 269         final static public ParamOptimize paramOptimize = new ParamOptimize();
 270
 271         final static private class Features implements Serializable
 272         {
 273                 private static final long serialVersionUID = 2689219384710526198L;
 274
 275                 FloatArray2DSIFT.Param p;
 276                 ArrayList< Feature > features;
 277                 Features( final FloatArray2DSIFT.Param p, final ArrayList< Feature > features )
 278                 {
 279                         this.p = p;
 280                         this.features = features;
 281                 }
 282         }
 283
 284         final static private class PointMatches implements Serializable
 285         {
 286                 private static final long serialVersionUID = -2564147268101223484L;
 287
 288                 Param p;
 289                 ArrayList< PointMatch > pointMatches;
 290                 PointMatches( final Param p, final ArrayList< PointMatch > pointMatches )
 291                 {
 292                         this.p = p;
 293                         this.pointMatches = pointMatches;
 294                 }
 295         }
 296
 297         /**
 298          * Extracts {@link Feature SIFT-features} from a {@link List} of
 299          * {@link AbstractAffineTile2D Tiles} and saves them to disk.
 300          */
 301         final static protected class ExtractFeaturesThread extends Thread
 302         {
 303                 final protected Param p;
 304                 final protected List< AbstractAffineTile2D< ? > > tiles;
 305 //              final protected HashMap< AbstractAffineTile2D< ? >, Collection< Feature > > tileFeatures;
 306                 final protected AtomicInteger ai;
 307                 final protected AtomicInteger ap;
 308                 final protected int steps;
 309
 310                 public ExtractFeaturesThread(
 311                                 final Param p,
 312                                 final List< AbstractAffineTile2D< ? > > tiles,
 313 //                              final HashMap< AbstractAffineTile2D< ? >, Collection< Feature > > tileFeatures,
 314                                 final AtomicInteger ai,
 315                                 final AtomicInteger ap,
 316                                 final int steps )
 317                 {
 318                         this.p = p;
 319                         this.tiles = tiles;
 320                         this.ai = ai;
 321                         this.ap = ap;
 322                         this.steps = steps;
 323                 }
 324
 325                 @Override
 326                 final public void run()
 327                 {
 328                         FloatArray2DSIFT sift = new FloatArray2DSIFT( p.sift );
 329                         SIFT ijSIFT = new SIFT( sift );
 330
 331                         for ( int i = ai.getAndIncrement(); i < tiles.size() && !isInterrupted(); i = ai.getAndIncrement() )
 332                         {
 333                                 if (isInterrupted()) return;
 334                                 AbstractAffineTile2D< ? > tile = tiles.get( i );
 335                                 Collection< Feature > features = deserializeFeatures( p, tile );
 336                                 if ( features == null )
 337                                 {
 338                                         /* extract features and, in case there is not enough memory available, try to free it and do again */
 339                                         boolean memoryFlushed;
 340                                         do
 341                                         {
 342                                                 try
 343                                                 {
 344                                                         features = new ArrayList< Feature >();
 345                                                         long s = System.currentTimeMillis();
 346                                                         ijSIFT.extractFeatures( tile.createMaskedByteImage(), features );
 347                                                         IJ.log( features.size() + " features extracted in tile " + i + " \"" + tile.getPatch().getTitle() + "\" (took " + ( System.currentTimeMillis() - s ) + " ms)." );
 348                                                         if ( !serializeFeatures( p, tile, features ) )
 349                                                                 IJ.log( "Saving features failed for tile \"" + tile.getPatch() + "\"" );
 350                                                         memoryFlushed = false;
 351                                                 }
 352                                                 catch ( OutOfMemoryError e )
 353                                                 {
 354                                                         Utils.log2( "Flushing memory for feature extraction" );
 355                                                         Loader.releaseAllCaches();
 356                                                         memoryFlushed = true;
 357                                                 }
 358                                         }
 359                                         while ( memoryFlushed );
 360                                 }
 361                                 else
 362                                 {
 363                                         IJ.log( features.size() + " features loaded for tile " + i + " \"" + tile.getPatch().getTitle() + "\"." );
 364                                 }
 365                                 IJ.showProgress( ap.getAndIncrement(), steps );
 366                         }
 367                 }
 368         }
 369
 370
 371         final static protected class MatchFeaturesAndFindModelThread extends Thread
 372         {
 373                 final protected Param p;
 374                 final protected List< AbstractAffineTile2D< ? > > tiles;
 375                 //final protected HashMap< AbstractAffineTile2D< ? >, Collection< Feature > > tileFeatures;
 376                 final protected List< AbstractAffineTile2D< ? >[] > tilePairs;
 377                 final protected AtomicInteger ai;
 378                 final protected AtomicInteger ap;
 379                 final protected int steps;
 380
 381                 public MatchFeaturesAndFindModelThread(
 382                                 final Param p,
 383                                 final List< AbstractAffineTile2D< ? > > tiles,
 384                                 final List< AbstractAffineTile2D< ? >[] > tilePairs,
 385                                 final AtomicInteger ai,
 386                                 final AtomicInteger ap,
 387                                 final int steps )
 388                 {
 389                         this.p = p;
 390                         this.tiles = tiles;
 391                         this.tilePairs = tilePairs;
 392                         this.ai = ai;
 393                         this.ap = ap;
 394                         this.steps = steps;
 395                 }
 396
 397                 @Override
 398                 final public void run()
 399                 {
 400                         final List< PointMatch > candidates = new ArrayList< PointMatch >();
 401
 402                         for ( int i = ai.getAndIncrement(); i < tilePairs.size() && !isInterrupted(); i = ai.getAndIncrement() )
 403                         {
 404                                 if (isInterrupted()) return;
 405                                 candidates.clear();
 406                                 final AbstractAffineTile2D< ? >[] tilePair = tilePairs.get( i );
 407
 408                                 Collection< PointMatch > inliers = deserializePointMatches( p, tilePair[ 0 ], tilePair[ 1 ] );
 409
 410                                 if ( inliers == null )
 411                                 {
 412                                         inliers = new ArrayList< PointMatch >();
 413
 414                                         long s = System.currentTimeMillis();
 415
 416                                         FeatureTransform.matchFeatures(
 417                                                 fetchFeatures( p, tilePair[ 0 ] ),
 418                                                 fetchFeatures( p, tilePair[ 1 ] ),
 419                                                 candidates,
 420                                                 p.rod );
 421
 422                                         /* find the model */
 423                                         final AbstractAffineModel2D< ? > model;
 424                                         switch ( p.expectedModelIndex )
 425                                         {
 426                                         case 0:
 427                                                 model = new TranslationModel2D();
 428                                                 break;
 429                                         case 1:
 430                                                 model = new RigidModel2D();
 431                                                 break;
 432                                         case 2:
 433                                                 model = new SimilarityModel2D();
 434                                                 break;
 435                                         case 3:
 436                                                 model = new AffineModel2D();
 437                                                 break;
 438                                         default:
 439                                                 return;
 440                                         }
 441
 442                                         boolean modelFound;
 443                                         try
 444                                         {
 445                                                 modelFound = model.filterRansac(
 446                                                                 candidates,
 447                                                                 inliers,
 448                                                                 1000,
 449                                                                 p.maxEpsilon,
 450                                                                 p.minInlierRatio,
 451                                                                 Math.max( 7, 3 * model.getMinNumMatches() ),
 452                                                                 3 );
 453                                         }
 454                                         catch ( NotEnoughDataPointsException e )
 455                                         {
 456                                                 modelFound = false;
 457                                         }
 458                                         if ( modelFound )
 459                                                 IJ.log( "Model found for tiles \"" + tilePair[ 0 ].getPatch() + "\" and \"" + tilePair[ 1 ].getPatch() + "\":\n  correspondences  " + inliers.size() + " of " + candidates.size() + "\n  average residual error  " + model.getCost() + " px\n  took " + ( System.currentTimeMillis() - s ) + " ms" );
 460                                         else
 461                                                 IJ.log( "No model found for tiles \"" + tilePair[ 0 ].getPatch() + "\" and \"" + tilePair[ 1 ].getPatch() + "\":\n  correspondence candidates  " + candidates.size() + "\n  took " + ( System.currentTimeMillis() - s ) + " ms" );
 462
 463                                         if ( !serializePointMatches( p, tilePair[ 0 ], tilePair[ 1 ], inliers ) )
 464                                                 IJ.log( "Saving point matches failed for tiles \"" + tilePair[ 0 ].getPatch() + "\" and \"" + tilePair[ 1 ].getPatch() + "\"" );
 465
 466                                 }
 467                                 else
 468                                         IJ.log( "Point matches for tiles \"" + tilePair[ 0 ].getPatch().getTitle() + "\" and \"" + tilePair[ 1 ].getPatch().getTitle() + "\" fetched from disk cache" );
 469
 470                                 if ( inliers != null && inliers.size() > 0 )
 471                                 {
 472                                         /* weight the inliers */
 473                                         for ( final PointMatch pm : inliers )
 474                                                 pm.setWeights( new float[]{ p.correspondenceWeight } );
 475
 476                                         synchronized ( tilePair[ 0 ] )
 477                                         {
 478                                                 synchronized ( tilePair[ 1 ] ) { tilePair[ 0 ].connect( tilePair[ 1 ], inliers ); }
 479                                                 tilePair[ 0 ].clearVirtualMatches();
 480                                         }
 481                                         synchronized ( tilePair[ 1 ] ) { tilePair[ 1 ].clearVirtualMatches(); }
 482                                 }
 483
 484                                 IJ.showProgress( ap.getAndIncrement(), steps );
 485                         }
 486                 }
 487         }
 488
 489         final static protected boolean serializeFeatures( final Param p, AbstractAffineTile2D< ? > t, final Collection< Feature > f )
 490         {
 491                 final ArrayList< Feature > list = new ArrayList< Feature >();
 492                 list.addAll( f );
 493                 final Patch patch = t.getPatch();
 494                 final Loader loader = patch.getProject().getLoader();
 495                 final Features fe = new Features( p.sift, list );
 496                 return loader.serialize( fe, new StringBuilder( loader.getUNUIdFolder() ).append( "features.ser/" )
 497                         .append( FSLoader.createIdPath( Long.toString( patch.getId() ), "features", ".ser" ) ).toString() );
 498         }
 499
 500         /**
 501          * Retrieve the features only if saved with the exact same relevant SIFT parameters.
 502          */
 503         final static protected Collection< Feature > deserializeFeatures( final Param p, final AbstractAffineTile2D< ? > t )
 504         {
 505                 final Patch patch = t.getPatch();
 506                 final Loader loader = patch.getProject().getLoader();
 507
 508                 final Object ob = loader.deserialize( new StringBuilder( loader.getUNUIdFolder() ).append( "features.ser/" )
 509                         .append( FSLoader.createIdPath( Long.toString( patch.getId() ), "features", ".ser" ) ).toString() );
 510                 if ( null != ob )
 511                 {
 512                         try
 513                         {
 514                                 final Features fe = ( Features )ob;
 515                                 if ( p.sift.equals( fe.p ) && null != fe.p )
 516                                 {
 517                                         return fe.features;
 518                                 }
 519                         }
 520                         catch ( Exception e )
 521                         {
 522                                 e.printStackTrace();
 523                         }
 524                 }
 525                 return null;
 526         }
 527
 528
 529         final static protected Collection< Feature > fetchFeatures(
 530                         final Param p,
 531                         final AbstractAffineTile2D< ? > t )
 532         {
 533                 Collection< Feature > features = deserializeFeatures( p, t );
 534                 if ( features == null )
 535                 {
 536                         final FloatArray2DSIFT sift = new FloatArray2DSIFT( p.sift );
 537                         final SIFT ijSIFT = new SIFT( sift );
 538                         features = new ArrayList< Feature >();
 539                         long s = System.currentTimeMillis();
 540                         ijSIFT.extractFeatures( t.createMaskedByteImage(), features );
 541                         IJ.log( features.size() + " features extracted in tile \"" + t.getPatch().getTitle() + "\" (took " + ( System.currentTimeMillis() - s ) + " ms)." );
 542                         if ( !serializeFeatures( p, t, features ) )
 543                                 IJ.log( "Saving features failed for tile: " + t.getPatch() );
 544                 }
 545                 return features;
 546         }
 547
 548
 549         /**
 550          * Save a {@link Collection} of {@link PointMatch PointMatches} two-sided.
 551          * Creates two serialization files which is desperately required to clean
 552          * up properly invalid serializations on change of a {@link Patch}.
 553          *
 554          * @param p
 555          * @param t1
 556          * @param t2
 557          * @param m
 558          * @return
 559          */
 560         final static protected boolean serializePointMatches(
 561                         final Param p,
 562                         final AbstractAffineTile2D< ? > t1,
 563                         final AbstractAffineTile2D< ? > t2,
 564                         final Collection< PointMatch > m )
 565         {
 566                 final ArrayList< PointMatch > list = new ArrayList< PointMatch >();
 567                 list.addAll( m );
 568                 final ArrayList< PointMatch > tsil = new ArrayList< PointMatch >();
 569                 PointMatch.flip( m, tsil );
 570                 final Patch p1 = t1.getPatch();
 571                 final Patch p2 = t2.getPatch();
 572                 final Loader loader = p1.getProject().getLoader();
 573                 return
 574                         loader.serialize(
 575                                 new PointMatches( p, list ),
 576                                 new StringBuilder( loader.getUNUIdFolder() ).append( "pointmatches.ser/" ).append( FSLoader.createIdPath( Long.toString( p1.getId() ) + "_" + Long.toString( p2.getId() ), "pointmatches", ".ser" ) ).toString() ) &&
 577                         loader.serialize(
 578                                 new PointMatches( p, tsil ),
 579                                 new StringBuilder( loader.getUNUIdFolder() ).append( "pointmatches.ser/" ).append( FSLoader.createIdPath( Long.toString( p2.getId() ) + "_" + Long.toString( p1.getId() ), "pointmatches", ".ser" ) ).toString() );
 580         }
 581
 582
 583         final static protected Collection< PointMatch > deserializePointMatches(
 584                         final Param p,
 585                         final AbstractAffineTile2D< ? > t1,
 586                         final AbstractAffineTile2D< ? > t2 )
 587         {
 588                 final Patch p1 = t1.getPatch();
 589                 final Patch p2 = t2.getPatch();
 590                 final Loader loader = p1.getProject().getLoader();
 591
 592                 final Object ob = loader.deserialize( new StringBuilder( loader.getUNUIdFolder() ).append( "pointmatches.ser/" )
 593                                 .append( FSLoader.createIdPath( Long.toString( p1.getId() ) + "_" + Long.toString( p2.getId() ), "pointmatches", ".ser" ) ).toString() );
 594
 595                 if ( null != ob )
 596                 {
 597                         try
 598                         {
 599                                 final PointMatches pm = ( PointMatches )ob;
 600                                 if ( p.equals( pm.p ) && null != pm.p )
 601                                 {
 602                                         return pm.pointMatches;
 603                                 }
 604                         }
 605                         catch ( Exception e )
 606                         {
 607                                 e.printStackTrace();
 608                         }
 609                 }
 610                 return null;
 611         }
 612
 613
 614         /**
 615          * Fetch a {@link Collection} of corresponding
 616          * {@link Feature SIFT-features}.  Both {@link Feature SIFT-features} and
 617          * {@linkplain PointMatch corresponding points} are cached to disk.
 618          *
 619          * @param p
 620          * @param t1
 621          * @param t2
 622          * @return
 623          *   <dl>
 624          *     <dt>null</dt><dd>if matching failed for some reasons</dd>
 625          *     <dt>empty {@link Collection}</dt><dd>if there was no consistent set
 626          *       of {@link PointMatch matches}</dd>
 627          *     <dt>{@link Collection} of {@link PointMatch PointMatches}</dt>
 628          *       <dd>if there was a consistent set of {@link PointMatch
 629          *         PointMatches}</dd>
 630          *   </dl>
 631          */
 632         final static protected Collection< PointMatch > fetchPointMatches(
 633                         final Param p,
 634                         final AbstractAffineTile2D< ? > t1,
 635                         final AbstractAffineTile2D< ? > t2 )
 636         {
 637                 Collection< PointMatch > pointMatches = deserializePointMatches( p, t1, t2 );
 638                 if ( pointMatches == null )
 639                 {
 640                         final List< PointMatch > candidates = new ArrayList< PointMatch >();
 641                         final List< PointMatch > inliers = new ArrayList< PointMatch >();
 642
 643                         long s = System.currentTimeMillis();
 644                         FeatureTransform.matchFeatures(
 645                                         fetchFeatures( p, t1 ),
 646                                         fetchFeatures( p, t2 ),
 647                                         candidates,
 648                                         p.rod );
 649
 650                         final AbstractAffineModel2D< ? > model;
 651                         switch ( p.expectedModelIndex )
 652                         {
 653                         case 0:
 654                                 model = new TranslationModel2D();
 655                                 break;
 656                         case 1:
 657                                 model = new RigidModel2D();
 658                                 break;
 659                         case 2:
 660                                 model = new SimilarityModel2D();
 661                                 break;
 662                         case 3:
 663                                 model = new AffineModel2D();
 664                                 break;
 665                         default:
 666                                 return null;
 667                         }
 668
 669                         boolean modelFound;
 670                         try
 671                         {
 672                                 modelFound = model.filterRansac(
 673                                                 candidates,
 674                                                 inliers,
 675                                                 1000,
 676                                                 p.maxEpsilon,
 677                                                 p.minInlierRatio,
 678                                                 3 * model.getMinNumMatches(),
 679                                                 3 );
 680                         }
 681                         catch ( NotEnoughDataPointsException e )
 682                         {
 683                                 modelFound = false;
 684                         }
 685
 686                         if ( modelFound )
 687                         {
 688                                 IJ.log( "Model found for tiles \"" + t1.getPatch() + "\" and \"" + t2.getPatch().getTitle() + "\":\n  correspondences  " + inliers.size() + " of " + candidates.size() + "\n  average residual error  " + model.getCost() + " px\n  took " + ( System.currentTimeMillis() - s ) + " ms" );
 689                         }
 690                         else
 691                                 IJ.log( "No model found for tiles " + t1 + " and " + t2 + "." );
 692
 693                         if ( !serializePointMatches( p, t1, t2, pointMatches ) )
 694                                 IJ.log( "Saving point matches failed for tile \"" + t1.getPatch() + "\" and tile \"" + t2.getPatch() + "\"" );
 695                 }
 696                 return pointMatches;
 697         }
 698
 699
 700         /**
 701          * Align a set of overlapping {@link AbstractAffineTile2D tiles} using
 702          * the following procedure:
 703          *
 704          * <ol>
 705          * <li>Extract {@link Feature SIFT-features} from all
 706          * {@link AbstractAffineTile2D tiles}.</li>
 707          * <li>Establish {@link PointMatch point-correspondences} from
 708          * consistent sets of {@link Feature feature} matches among overlapping
 709          * tiles.</li>
 710          * <li>Globally align the tile configuration.</li>
 711          * </ol>
 712          *
 713          * Both
 714          * {@link SIFT#extractFeatures(ij.process.ImageProcessor, Collection) feature extraction}
 715          * and {@link FeatureTransform#matchFeatures(Collection, Collection, List, float) matching}
 716          * are executed in multiple {@link Thread Threads}, with the number of
 717          * {@link Thread Threads} being a parameter of the method.
 718          *
 719          * @param p
 720          * @param tiles
 721          * @param numThreads
 722          */
 723         final static public void alignTiles(
 724                         final ParamOptimize p,
 725                         final List< AbstractAffineTile2D< ? > > tiles,
 726                         final List< AbstractAffineTile2D< ? > > fixedTiles,
 727                         final int numThreads )
 728         {
 729                 final ArrayList< AbstractAffineTile2D< ? >[] > tilePairs = new ArrayList< AbstractAffineTile2D<?>[] >();
 730                 AbstractAffineTile2D.pairOverlappingTiles( tiles, tilePairs );
 731                 connectTilePairs( p, tiles, tilePairs, numThreads );
 732                 optimizeTileConfiguration( p, tiles, fixedTiles );
 733         }
 734
 735         /**
 736          * Align a set of {@link AbstractAffineTile2D tiles} that are
 737          * interconnected by {@link PointMatch point-correspondences}.
 738          */
 739         final static public void optimizeTileConfiguration(
 740                         final ParamOptimize p,
 741                         final List< AbstractAffineTile2D< ? > > tiles,
 742                         final List< AbstractAffineTile2D< ? > > fixedTiles )
 743         {
 744                 final TileConfiguration tc = new TileConfiguration();
 745                 tc.addTiles( tiles );
 746
 747                 ArrayList< Set< Tile< ? > > > graphs = Tile.identifyConnectedGraphs( tiles );
 748                 for ( Set< Tile< ? > > graph : graphs )
 749                 {
 750                         boolean pleaseFix = true;
 751                         if ( fixedTiles != null )
 752                                 for ( final Tile< ? > t : fixedTiles )
 753                                         if ( graph.contains( t ) )
 754                                         {
 755                                                 pleaseFix = false;
 756                                         }
 757                         if ( pleaseFix )
 758                                 tc.fixTile( graph.iterator().next() );
 759                 }
 760                 for ( final Tile< ? > t : fixedTiles )
 761                         tc.fixTile( t );
 762
 763                 try
 764                 {
 765                         tc.optimize( p.maxEpsilon, p.maxIterations, p.maxPlateauwidth );
 766                 }
 767                 catch ( Exception e ) { IJ.error( e.getMessage() + " " + e.getStackTrace() ); }
 768         }
 769
 770
 771         final static protected void pairwiseAlign(
 772                         AbstractAffineTile2D< ? > tile,
 773                         Set< AbstractAffineTile2D< ? > > visited )
 774         {
 775                 visited.add( tile );
 776                 for ( Tile< ? > t : tile.getConnectedTiles() )
 777                 {
 778                         if ( visited.contains( t ) ) continue;
 779                         pairwiseAlign( ( AbstractAffineTile2D< ? > )t, visited );
 780                         // TODO Actually do it ...
 781                 }
 782         }
 783
 784
 785         final static public void pairwiseAlignTileConfiguration(
 786                         final List< AbstractAffineTile2D< ? > > tiles )
 787         {
 788                 // TODO Implement it
 789         }
 790
 791
 792         /**
 793          * Connect a {@link List} of {@link AbstractAffineTile2D Tiles} by
 794          * geometrically consistent {@link Feature SIFT-feature} correspondences.
 795          *
 796          * @param p
 797          * @param tiles
 798          * @param numThreads
 799          */
 800         final static public void connectTilePairs(
 801                         final Param p,
 802                         final List< AbstractAffineTile2D< ? > > tiles,
 803                         final List< AbstractAffineTile2D< ? >[] > tilePairs,
 804                         final int numThreads )
 805         {
 806                 final AtomicInteger ai = new AtomicInteger( 0 );
 807                 final AtomicInteger ap = new AtomicInteger( 0 );
 808                 final int steps = tiles.size() + tilePairs.size();
 809                 final List< ExtractFeaturesThread > extractFeaturesThreads = new ArrayList< ExtractFeaturesThread >();
 810                 final List< MatchFeaturesAndFindModelThread > matchFeaturesAndFindModelThreads = new ArrayList< MatchFeaturesAndFindModelThread >();
 811
 812                 /** Extract and save Features */
 813                 for ( int i = 0; i < numThreads; ++i )
 814                 {
 815                         final ExtractFeaturesThread thread = new ExtractFeaturesThread( p.clone(), tiles, ai, ap, steps );
 816                         extractFeaturesThreads.add( thread );
 817                         thread.start();
 818                 }
 819                 try
 820                 {
 821                         for ( final ExtractFeaturesThread thread : extractFeaturesThreads )
 822                                 thread.join();
 823                 }
 824                 catch ( InterruptedException e )
 825                 {
 826                         IJ.log( "Feature extraction failed.\n" + e.getMessage() + "\n" + e.getStackTrace() );
 827                         return;
 828                 }
 829
 830                 /** Establish correspondences */
 831                 ai.set( 0 );
 832                 for ( int i = 0; i < numThreads; ++i )
 833                 {
 834                         MatchFeaturesAndFindModelThread thread = new MatchFeaturesAndFindModelThread( p.clone(), tiles, tilePairs, ai, ap, steps );
 835                         matchFeaturesAndFindModelThreads.add( thread );
 836                         thread.start();
 837                 }
 838                 try
 839                 {
 840                         for ( final MatchFeaturesAndFindModelThread thread : matchFeaturesAndFindModelThreads )
 841                                 thread.join();
 842                 }
 843                 catch ( InterruptedException e )
 844                 {
 845                         IJ.log( "Establishing feature correspondences failed.\n" + e.getMessage() + "\n" + e.getStackTrace() );
 846                         return;
 847                 }
 848         }
 849
 850
 851         /**
 852          *
 853          * @param p
 854          * @param patches
 855          * @param fixedPatches
 856          * @param tiles will contain the generated
 857          *   {@link AbstractAffineTile2D Tiles}
 858          * @param fixedTiles will contain the {@link AbstractAffineTile2D Tiles}
 859          *   corresponding to the {@link Patch Patches} in fixedPatches
 860          */
 861         final static public void tilesFromPatches(
 862                         final Param p,
 863                         final List< ? extends Patch > patches,
 864                         final List< ? extends Patch > fixedPatches,
 865                         final List< AbstractAffineTile2D< ? > > tiles,
 866                         final List< AbstractAffineTile2D< ? > > fixedTiles )
 867         {
 868                 for ( final Patch patch : patches )
 869                 {
 870                         final AbstractAffineTile2D< ? > t;
 871                         switch ( p.desiredModelIndex )
 872                         {
 873                         case 0:
 874                                 t = new TranslationTile2D( patch );
 875                                 break;
 876                         case 1:
 877                                 t = new RigidTile2D( patch );
 878                                 break;
 879                         case 2:
 880                                 t = new SimilarityTile2D( patch );
 881                                 break;
 882                         case 3:
 883                                 t = new AffineTile2D( patch );
 884                                 break;
 885                         default:
 886                                 return;
 887                         }
 888                         tiles.add( t );
 889                         if ( ( fixedPatches != null && fixedPatches.contains( patch ) ) || patch.isLocked() )
 890                                 fixedTiles.add( t );
 891                 }
 892         }
 893
 894
 895         /**
 896          * Align a selection of {@link Patch patches} in a Layer.
 897          *
 898          * @param layer
 899          */
 900         final static public void alignSelectedPatches( Selection selection, final int numThreads )
 901         {
 902                 final List< Patch > patches = new ArrayList< Patch >();
 903                 for ( final Displayable d : selection.getSelected() )
 904                         if ( d instanceof Patch ) patches.add(  ( Patch )d );
 905
 906                 if ( patches.size() < 2 ) return;
 907
 908                 if ( !paramOptimize.setup( "Align selected patches" ) ) return;
 909
 910                 List< AbstractAffineTile2D< ? > > tiles = new ArrayList< AbstractAffineTile2D< ? > >();
 911                 List< AbstractAffineTile2D< ? > > fixedTiles = new ArrayList< AbstractAffineTile2D< ? > >();
 912                 List< Patch > fixedPatches = new ArrayList< Patch >();
 913                 final Displayable active = selection.getActive();
 914                 if ( active != null && active instanceof Patch )
 915                         fixedPatches.add( ( Patch )active );
 916                 tilesFromPatches( paramOptimize, patches, fixedPatches, tiles, fixedTiles );
 917
 918                 alignTiles( paramOptimize, tiles, fixedTiles, numThreads );
 919
 920                 for ( AbstractAffineTile2D< ? > t : tiles )
 921                         t.getPatch().setAffineTransform( t.getModel().createAffine() );
 922         }
 923
 924
 925         /**
 926          * Align all {@link Patch patches} in a Layer.
 927          *
 928          * @param layer
 929          */
 930         final static public void alignLayer( final Layer layer, final int numThreads )
 931         {
 932                 if ( !paramOptimize.setup( "Align patches in layer" ) ) return;
 933
 934                 List< Displayable > displayables = layer.getDisplayables( Patch.class );
 935                 List< Patch > patches = new ArrayList< Patch >();
 936                 for ( Displayable d : displayables )
 937                         patches.add( ( Patch )d );
 938                 List< AbstractAffineTile2D< ? > > tiles = new ArrayList< AbstractAffineTile2D< ? > >();
 939                 List< AbstractAffineTile2D< ? > > fixedTiles = new ArrayList< AbstractAffineTile2D< ? > >();
 940                 tilesFromPatches( paramOptimize, patches, null, tiles, fixedTiles );
 941
 942                 alignTiles( paramOptimize, tiles, fixedTiles, numThreads );
 943
 944                 for ( AbstractAffineTile2D< ? > t : tiles )
 945                         t.getPatch().setAffineTransform( t.getModel().createAffine() );
 946         }
 947
 948
 949         final static public void alignLayersLinearly( final List< Layer > layers, final int numThreads )
 950         {
 951                 param.sift.maxOctaveSize = 1600;
 952
 953                 if ( !param.setup( "Align layers linearly" ) ) return;
 954
 955                 final Rectangle box = layers.get( 0 ).getParent().getMinimalBoundingBox( Patch.class );
 956                 final float scale = Math.min(  1.0f, Math.min( ( float )param.sift.maxOctaveSize / ( float )box.width, ( float )param.sift.maxOctaveSize / ( float )box.height ) );
 957                 final Param p = param.clone();
 958                 p.maxEpsilon *= scale;
 959
 960                 final FloatArray2DSIFT sift = new FloatArray2DSIFT( p.sift );
 961                 final SIFT ijSIFT = new SIFT( sift );
 962
 963                 Rectangle box1 = null;
 964                 Rectangle box2 = null;
 965                 final Collection< Feature > features1 = new ArrayList< Feature >();
 966                 final Collection< Feature > features2 = new ArrayList< Feature >();
 967                 List< PointMatch > candidates = new ArrayList< PointMatch >();
 968                 List< PointMatch > inliers = new ArrayList< PointMatch >();
 969
 970                 AffineTransform a = new AffineTransform();
 971
 972                 int i = 0;
 973                 for ( Layer l : layers )
 974                 {
 975                         long s = System.currentTimeMillis();
 976
 977                         features1.clear();
 978                         features1.addAll( features2 );
 979                         features2.clear();
 980
 981                         final Rectangle box3 = l.getMinimalBoundingBox( Patch.class );
 982
 983                         if ( box3 == null ) continue;
 984
 985                         box1 = box2;
 986                         box2 = box3;
 987
 988                         ijSIFT.extractFeatures(
 989                                         l.getProject().getLoader().getFlatImage( l, box2, scale, 0xffffffff, ImagePlus.GRAY8, Patch.class, true ).getProcessor(),
 990                                         features2 );
 991                         IJ.log( features2.size() + " features extracted in layer \"" + l.getTitle() + "\" (took " + ( System.currentTimeMillis() - s ) + " ms)." );
 992
 993                         if ( features1.size() > 0 )
 994                         {
 995                                 s = System.currentTimeMillis();
 996
 997                                 candidates.clear();
 998
 999                                 FeatureTransform.matchFeatures(
1000                                         features2,
1001                                         features1,
1002                                         candidates,
1003                                         p.rod );
1004
1005                                 final AbstractAffineModel2D< ? > model;
1006                                 switch ( p.expectedModelIndex )
1007                                 {
1008                                 case 0:
1009                                         model = new TranslationModel2D();
1010                                         break;
1011                                 case 1:
1012                                         model = new RigidModel2D();
1013                                         break;
1014                                 case 2:
1015                                         model = new SimilarityModel2D();
1016                                         break;
1017                                 case 3:
1018                                         model = new AffineModel2D();
1019                                         break;
1020                                 default:
1021                                         return;
1022                                 }
1023
1024                                 boolean modelFound;
1025                                 try
1026                                 {
1027                                         modelFound = model.filterRansac(
1028                                                         candidates,
1029                                                         inliers,
1030                                                         1000,
1031                                                         p.maxEpsilon,
1032                                                         p.minInlierRatio,
1033                                                         3 * model.getMinNumMatches(),
1034                                                         3 );
1035                                 }
1036                                 catch ( NotEnoughDataPointsException e )
1037                                 {
1038                                         modelFound = false;
1039                                 }
1040
1041                                 if ( modelFound )
1042                                 {
1043                                         IJ.log( "Model found for layer \"" + l.getTitle() + "\" and its predecessor:\n  correspondences  " + inliers.size() + " of " + candidates.size() + "\n  average residual error  " + ( model.getCost() / scale ) + " px\n  took " + ( System.currentTimeMillis() - s ) + " ms" );
1044                                         final AffineTransform b = new AffineTransform();
1045                                         b.translate( box1.x, box1.y );
1046                                         b.scale( 1.0f / scale, 1.0f / scale );
1047                                         b.concatenate( model.createAffine() );
1048                                         b.scale( scale, scale );
1049                                         b.translate( -box2.x, -box2.y);
1050
1051                                         a.concatenate( b );
1052                                         l.apply( Displayable.class, a );
1053                                         Display.repaint( l );
1054                                 }
1055                                 else
1056                                         IJ.log( "No model found for layer \"" + l.getTitle() + "\" and its predecessor." );
1057                         }
1058                         IJ.showProgress( ++i, layers.size() );
1059                 }
1060         }
1061
1062         /**
1063          * Temporary helper method that creates
1064          * @param matches
1065          * @param alpha
1066          * @return
1067          * @throws Exception
1068          */
1069         final static public MovingLeastSquaresTransform createMLST( final Collection< PointMatch > matches, final float alpha ) throws Exception
1070         {
1071                 final MovingLeastSquaresTransform mlst = new MovingLeastSquaresTransform();
1072                 mlst.setAlpha( 1.0f );
1073                 Class< ? extends AbstractAffineModel2D< ? > > c = AffineModel2D.class;
1074                 switch ( matches.size() )
1075                 {
1076                         case 1:
1077                                 c = TranslationModel2D.class;
1078                                 break;
1079                         case 2:
1080                                 c = SimilarityModel2D.class;
1081                                 break;
1082                         default:
1083                                 break;
1084                 }
1085                 mlst.setModel( c );
1086                 mlst.setMatches( matches );
1087
1088                 return mlst;
1089         }
1090
1091
1092         /**
1093          * Align two collections of tiles
1094          * @param p
1095          * @param a
1096          * @param b
1097          */
1098         final static public void alignTileCollections( final Param p, final Collection< AbstractAffineTile2D< ? > > a, final Collection< AbstractAffineTile2D< ? > > b )
1099         {
1100                 final ArrayList< Patch > pa = new ArrayList< Patch >();
1101                 final ArrayList< Patch > pb = new ArrayList< Patch >();
1102                 for ( final AbstractAffineTile2D< ? > t : a )
1103                         pa.add( t.getPatch() );
1104                 for ( final AbstractAffineTile2D< ? > t : b )
1105                         pb.add( t.getPatch() );
1106
1107                 final Layer la = pa.iterator().next().getLayer();
1108                 final Layer lb = pb.iterator().next().getLayer();
1109
1110                 final Rectangle boxA = Displayable.getBoundingBox( pa, null );
1111                 final Rectangle boxB = Displayable.getBoundingBox( pb, null );
1112
1113                 final float scale = Math.min(
1114                                 1.0f,
1115                                 Math.min(
1116                                         Math.min(
1117                                                         ( float )p.sift.maxOctaveSize / ( float )boxA.width,
1118                                                         ( float )p.sift.maxOctaveSize / ( float )boxA.height ),
1119                                         Math.min(
1120                                                         ( float )p.sift.maxOctaveSize / ( float )boxB.width,
1121                                                         ( float )p.sift.maxOctaveSize / ( float )boxB.height ) ) );
1122
1123                 final Param pp = p.clone();
1124                 pp.maxEpsilon *= scale;
1125
1126                 final FloatArray2DSIFT sift = new FloatArray2DSIFT( pp.sift );
1127                 final SIFT ijSIFT = new SIFT( sift );
1128
1129                 final Collection< Feature > featuresA = new ArrayList< Feature >();
1130                 final Collection< Feature > featuresB = new ArrayList< Feature >();
1131                 List< PointMatch > candidates = new ArrayList< PointMatch >();
1132                 List< PointMatch > inliers = new ArrayList< PointMatch >();
1133
1134                 long s = System.currentTimeMillis();
1135                 ijSIFT.extractFeatures(
1136                                 la.getProject().getLoader().getFlatImage( la, boxA, scale, 0xffffffff, ImagePlus.GRAY8, null, pa, true, Color.GRAY ).getProcessor(), featuresA );
1137                 Utils.log( featuresA.size() + " features extracted in graph A in layer \"" + la.getTitle() + "\" (took " + ( System.currentTimeMillis() - s ) + " ms)." );
1138
1139                 s = System.currentTimeMillis();
1140                 ijSIFT.extractFeatures(
1141                                 lb.getProject().getLoader().getFlatImage( lb, boxB, scale, 0xffffffff, ImagePlus.GRAY8, null, pb, true, Color.GRAY ).getProcessor(), featuresB );
1142                 Utils.log( featuresB.size() + " features extracted in graph B in layer \"" + lb.getTitle() + "\" (took " + ( System.currentTimeMillis() - s ) + " ms)." );
1143
1144                 if ( featuresA.size() > 0 && featuresB.size() > 0 )
1145                 {
1146                         s = System.currentTimeMillis();
1147                         FeatureTransform.matchFeatures(
1148                                         featuresA,
1149                                         featuresB,
1150                                         candidates,
1151                                         pp.rod );
1152
1153                         final AbstractAffineModel2D< ? > model;
1154                         switch ( p.expectedModelIndex )
1155                         {
1156                         case 0:
1157                                 model = new TranslationModel2D();
1158                                 break;
1159                         case 1:
1160                                 model = new RigidModel2D();
1161                                 break;
1162                         case 2:
1163                                 model = new SimilarityModel2D();
1164                                 break;
1165                         case 3:
1166                                 model = new AffineModel2D();
1167                                 break;
1168                         default:
1169                                 return;
1170                         }
1171
1172                         boolean modelFound;
1173                         try
1174                         {
1175                                 modelFound = model.filterRansac(
1176                                                 candidates,
1177                                                 inliers,
1178                                                 1000,
1179                                                 p.maxEpsilon,
1180                                                 p.minInlierRatio,
1181                                                 3 * model.getMinNumMatches(),
1182                                                 3 );
1183                         }
1184                         catch ( NotEnoughDataPointsException e )
1185                         {
1186                                 modelFound = false;
1187                         }
1188
1189                         if ( modelFound )
1190                         {
1191                                 IJ.log( "Model found for graph A and B in layers \"" + la.getTitle() + "\" and \"" + lb.getTitle() + "\":\n  correspondences  " + inliers.size() + " of " + candidates.size() + "\n  average residual error  " + ( model.getCost() / scale ) + " px\n  took " + ( System.currentTimeMillis() - s ) + " ms" );
1192                                 final AffineTransform at = new AffineTransform();
1193                                 at.translate( boxA.x, boxA.y );
1194                                 at.scale( 1.0f / scale, 1.0f / scale );
1195                                 at.concatenate( model.createAffine() );
1196                                 at.scale( scale, scale );
1197                                 at.translate( -boxB.x, -boxB.y);
1198
1199                                 for ( final Patch t : pa )
1200                                         t.preTransform( at, false );
1201                                 Display.repaint( la );
1202                         }
1203                         else
1204                                 IJ.log( "No model found for graph A and B in layers \"" + la.getTitle() + "\" and \"" + lb.getTitle() + "\"." );
1205                 }
1206         }
1207 }