TrakEM2_/src/main/java/mpicbg/trakem2/align/Align.java

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