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