?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
use MovingLeastSquaresTransform2 in ElasticLayerAlignment
[trakem2.git] / mpicbg / trakem2 / align / ElasticLayerAlignment.java
blob09bfed0fb542fcac768c08a9409b559a683b568d
1 /**
2 * License: GPL
3 *
4 * This program is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU General Public License 2
6 * as published by the Free Software Foundation.
7 *
8 * This program is distributed in the hope that it will be useful,
9 * but WITHOUT ANY WARRANTY; without even the implied warranty of
10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 * GNU General Public License for more details.
12 *
13 * You should have received a copy of the GNU General Public License
14 * along with this program; if not, write to the Free Software
15 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
16 */
17 package mpicbg.trakem2.align;
18
19
20 import ij.IJ;
21 import ij.ImagePlus;
22 import ij.gui.GenericDialog;
23 import ij.process.FloatProcessor;
24 import ij.process.ImageProcessor;
25 import ini.trakem2.display.Layer;
26 import ini.trakem2.display.LayerSet;
27 import ini.trakem2.display.Patch;
28 import ini.trakem2.utils.Filter;
29 import ini.trakem2.utils.Utils;
30
31 import java.awt.Color;
32 import java.awt.Image;
33 import java.awt.Rectangle;
34 import java.io.Serializable;
35 import java.util.ArrayList;
36 import java.util.Collection;
37 import java.util.Iterator;
38 import java.util.List;
39 import java.util.concurrent.Callable;
40 import java.util.concurrent.ExecutionException;
41 import java.util.concurrent.ExecutorService;
42 import java.util.concurrent.Executors;
43 import java.util.concurrent.Future;
44 import java.util.concurrent.atomic.AtomicInteger;
45
46 import mpicbg.ij.SIFT;
47 import mpicbg.ij.blockmatching.BlockMatching;
48 import mpicbg.imagefeatures.Feature;
49 import mpicbg.imagefeatures.FloatArray2DSIFT;
50 import mpicbg.models.AbstractModel;
51 import mpicbg.models.AffineModel2D;
52 import mpicbg.models.ErrorStatistic;
53 import mpicbg.models.HomographyModel2D;
54 import mpicbg.models.InvertibleCoordinateTransform;
55 import mpicbg.models.NotEnoughDataPointsException;
56 import mpicbg.models.Point;
57 import mpicbg.models.PointMatch;
58 import mpicbg.models.RigidModel2D;
59 import mpicbg.models.SimilarityModel2D;
60 import mpicbg.models.Spring;
61 import mpicbg.models.SpringMesh;
62 import mpicbg.models.Tile;
63 import mpicbg.models.TileConfiguration;
64 import mpicbg.models.Transforms;
65 import mpicbg.models.TranslationModel2D;
66 import mpicbg.models.Vertex;
67 import mpicbg.trakem2.transform.MovingLeastSquaresTransform;
68 import mpicbg.trakem2.transform.MovingLeastSquaresTransform2;
69
70 /**
71 * @author Stephan Saalfeld <saalfeld@mpi-cbg.de>
72 * @version 0.1a
73 */
74 public class ElasticLayerAlignment extends AbstractElasticAlignment
75 {
76 final static protected class Param implements Serializable
77 {
78 private static final long serialVersionUID = 875358335352119865L;
79
80 final public ParamPointMatch ppm = new ParamPointMatch();
81 {
82 ppm.sift.fdSize = 8;
83 }
84
85 /**
86 * Maximal accepted alignment error in px
87 */
88 public float maxEpsilon = 200.0f;
89
90 /**
91 * Inlier/candidates ratio
92 */
93 public float minInlierRatio = 0.0f;
94
95 /**
96 * Minimal absolute number of inliers
97 */
98 public int minNumInliers = 12;
99
100 /**
101 * Transformation models for choice
102 */
103 final static public String[] modelStrings = new String[]{ "Translation", "Rigid", "Similarity", "Affine", "Perspective" };
104 public int modelIndex = 3;
105
106 /**
107 * Ignore identity transform up to a given tolerance
108 */
109 public boolean rejectIdentity = true;
110 public float identityTolerance = 5.0f;
111
112 /**
113 * Maximal number of consecutive sections to be tested for an alignment model
114 */
115 public int maxNumNeighbors = 10;
116
117 /**
118 * Maximal number of consecutive slices for which no model could be found
119 */
120 public int maxNumFailures = 3;
121
122 public float layerScale = 0.1f;
123 public float minR = 0.8f;
124 public float maxCurvatureR = 3f;
125 public float rodR = 0.8f;
126
127 public int modelIndexOptimize = 1;
128 public int maxIterationsOptimize = 1000;
129 public int maxPlateauwidthOptimize = 200;
130
131 public int resolutionSpringMesh = 16;
132 public float stiffnessSpringMesh = 0.1f;
133 public float dampSpringMesh = 0.6f;
134 public float maxStretchSpringMesh = 2000.0f;
135 public int maxIterationsSpringMesh = 1000;
136 public int maxPlateauwidthSpringMesh = 200;
137
138 public boolean visualize = false;
139
140
141 public int maxNumThreads = Runtime.getRuntime().availableProcessors();
142
143 public boolean setup()
144 {
145 /* SIFT */
146 final GenericDialog gd = new GenericDialog( "Elastically align layers: SIFT parameters" );
147
148 SIFT.addFields( gd, ppm.sift );
149
150 gd.addNumericField( "closest/next_closest_ratio :", ppm.rod, 2 );
151
152 gd.addMessage( "Geometric Consensus Filter:" );
153 gd.addNumericField( "maximal_alignment_error :", maxEpsilon, 2, 6, "px" );
154 gd.addNumericField( "minimal_inlier_ratio :", minInlierRatio, 2 );
155 gd.addNumericField( "minimal_number_of_inliers :", minNumInliers, 0 );
156 gd.addChoice( "approximate_transformation :", Param.modelStrings, Param.modelStrings[ modelIndex ] );
157 gd.addCheckbox( "ignore constant background", rejectIdentity );
158 gd.addNumericField( "tolerance :", identityTolerance, 2, 6, "px" );
159
160 gd.addMessage( "Matching:" );
161 gd.addNumericField( "test_maximally :", maxNumNeighbors, 0, 6, "layers" );
162 gd.addNumericField( "give_up_after :", maxNumFailures, 0, 6, "failures" );
163
164 gd.addMessage( "Caching:" );
165 gd.addCheckbox( "clear_cache", ppm.clearCache );
166
167 gd.showDialog();
168
169 if ( gd.wasCanceled() )
170 return false;
171
172 SIFT.readFields( gd, ppm.sift );
173
174 ppm.rod = ( float )gd.getNextNumber();
175
176 maxEpsilon = ( float )gd.getNextNumber();
177 minInlierRatio = ( float )gd.getNextNumber();
178 minNumInliers = ( int )gd.getNextNumber();
179 modelIndex = gd.getNextChoiceIndex();
180 rejectIdentity = gd.getNextBoolean();
181 identityTolerance = ( float )gd.getNextNumber();
182 maxNumNeighbors = ( int )gd.getNextNumber();
183 maxNumFailures = ( int )gd.getNextNumber();
184
185 ppm.clearCache = gd.getNextBoolean();
186
187
188 /* Block Matching */
189 final GenericDialog gdBlockMatching = new GenericDialog( "Elastically align layers: Block Matching parameters" );
190 gdBlockMatching.addMessage( "Block Matching:" );
191
192 /* TODO suggest isotropic resolution for this parameter */
193 gdBlockMatching.addNumericField( "layer_scale :", layerScale, 2 );
194 gdBlockMatching.addNumericField( "minimal_PMCC_r :", minR, 2 );
195 gdBlockMatching.addNumericField( "maximal_curvature_ratio :", maxCurvatureR, 2 );
196 gdBlockMatching.addNumericField( "maximal_second_best_r/best_r :", rodR, 2 );
197
198 /* TODO suggest a resolution that matches maxEpsilon */
199 gdBlockMatching.addNumericField( "resolution :", resolutionSpringMesh, 0 );
200
201 gdBlockMatching.showDialog();
202
203 if ( gdBlockMatching.wasCanceled() )
204 return false;
205
206 layerScale = ( float )gdBlockMatching.getNextNumber();
207 minR = ( float )gdBlockMatching.getNextNumber();
208 maxCurvatureR = ( float )gdBlockMatching.getNextNumber();
209 rodR = ( float )gdBlockMatching.getNextNumber();
210 resolutionSpringMesh = ( int )gdBlockMatching.getNextNumber();
211
212 /* Optimization */
213 final GenericDialog gdOptimize = new GenericDialog( "Elastically align layers: Optimization" );
214
215 gdOptimize.addMessage( "Approximate Optimizer:" );
216 gdOptimize.addChoice( "approximate_transformation :", Param.modelStrings, Param.modelStrings[ modelIndexOptimize ] );
217 gdOptimize.addNumericField( "maximal_iterations :", maxIterationsOptimize, 0 );
218 gdOptimize.addNumericField( "maximal_plateauwidth :", maxPlateauwidthOptimize, 0 );
219
220 gdOptimize.addMessage( "Spring Mesh:" );
221 gdOptimize.addNumericField( "stiffness :", stiffnessSpringMesh, 2 );
222 gdOptimize.addNumericField( "maximal_stretch :", maxStretchSpringMesh, 2, 6, "px" );
223 gdOptimize.addNumericField( "maximal_iterations :", maxIterationsSpringMesh, 0 );
224 gdOptimize.addNumericField( "maximal_plateauwidth :", maxPlateauwidthSpringMesh, 0 );
225
226 gdOptimize.showDialog();
227
228 if ( gdOptimize.wasCanceled() )
229 return false;
230
231 modelIndexOptimize = gdOptimize.getNextChoiceIndex();
232 maxIterationsOptimize = ( int )gdOptimize.getNextNumber();
233 maxPlateauwidthOptimize = ( int )gdOptimize.getNextNumber();
234
235 stiffnessSpringMesh = ( float )gdOptimize.getNextNumber();
236 maxStretchSpringMesh = ( float )gdOptimize.getNextNumber();
237 maxIterationsSpringMesh = ( int )gdOptimize.getNextNumber();
238 maxPlateauwidthSpringMesh = ( int )gdOptimize.getNextNumber();
239
240 return true;
241 }
242 }
243
244 final static Param p = new Param();
245
246
247 final static private String layerName( final Layer layer )
248 {
249 return new StringBuffer( "layer z=" )
250 .append( String.format( "%.3f", layer.getZ() ) )
251 .append( " `" )
252 .append( layer.getTitle() )
253 .append( "'" )
254 .toString();
255
256 }
257
258 /**
259 * Extract SIFT features and save them into the project folder.
260 *
261 * @param layerSet the layerSet that contains all layers
262 * @param layerRange the list of layers to be aligned
263 * @param box a rectangular region of interest that will be used for alignment
264 * @param scale scale factor <= 1.0
265 * @param filter a name based filter for Patches (can be null)
266 * @param p SIFT extraction parameters
267 * @throws Exception
268 */
269 final static protected void extractAndSaveLayerFeatures(
270 final LayerSet layerSet,
271 final List< Layer > layerRange,
272 final Rectangle box,
273 final float scale,
274 final Filter< Patch > filter,
275 final FloatArray2DSIFT.Param siftParam,
276 final boolean clearCache ) throws ExecutionException, InterruptedException
277 {
278 final ExecutorService exec = Executors.newFixedThreadPool( p.maxNumThreads );
279
280 /* extract features for all slices and store them to disk */
281 final AtomicInteger counter = new AtomicInteger( 0 );
282 final ArrayList< Future< ArrayList< Feature > > > siftTasks = new ArrayList< Future< ArrayList< Feature > > >();
283
284 for ( int i = 0; i < layerRange.size(); ++i )
285 {
286 final int layerIndex = i;
287 final Rectangle finalBox = box;
288 siftTasks.add(
289 exec.submit( new Callable< ArrayList< Feature > >()
290 {
291 public ArrayList< Feature > call()
292 {
293 final Layer layer = layerRange.get( layerIndex );
294
295 final String layerName = layerName( layer );
296
297 IJ.showProgress( counter.getAndIncrement(), layerRange.size() - 1 );
298
299 final List< Patch > patches = filterPatches( layer, filter );
300
301 ArrayList< Feature > fs = null;
302 if ( !clearCache )
303 fs = mpicbg.trakem2.align.Util.deserializeFeatures( layerSet.getProject(), siftParam, "layer", layer.getId() );
304
305 if ( null == fs )
306 {
307 /* free memory */
308 layer.getProject().getLoader().releaseAll();
309
310 final FloatArray2DSIFT sift = new FloatArray2DSIFT( siftParam );
311 final SIFT ijSIFT = new SIFT( sift );
312 fs = new ArrayList< Feature >();
313 final ImageProcessor ip = layer.getProject().getLoader().getFlatImage( layer, finalBox, scale, 0xffffffff, ImagePlus.GRAY8, Patch.class, patches, true ).getProcessor();
314 ijSIFT.extractFeatures( ip, fs );
315 Utils.log( fs.size() + " features extracted for " + layerName );
316
317 if ( !mpicbg.trakem2.align.Util.serializeFeatures( layerSet.getProject(), siftParam, "layer", layer.getId(), fs ) )
318 Utils.log( "FAILED to store serialized features for " + layerName );
319 }
320 else
321 Utils.log( fs.size() + " features loaded for " + layerName );
322
323 return fs;
324 }
325 } ) );
326 }
327
328 /* join */
329 try
330 {
331 for ( Future< ArrayList< Feature > > fu : siftTasks )
332 fu.get();
333 }
334 catch ( InterruptedException e )
335 {
336 Utils.log( "Feature extraction interrupted." );
337 siftTasks.clear();
338 exec.shutdown();
339 throw e;
340 }
341 catch ( ExecutionException e )
342 {
343 Utils.log( "Execution exception during feature extraction." );
344 siftTasks.clear();
345 exec.shutdown();
346 throw e;
347 }
348
349 siftTasks.clear();
350 exec.shutdown();
351 }
352
353
354
355 /**
356 * Stateful. Changing the parameters of this instance. Do not use in parallel.
357 *
358 * @param layerSet
359 * @param first
360 * @param last
361 * @param propagateTransform
362 * @param fov
363 * @param filter
364 */
365 final public void exec(
366 final LayerSet layerSet,
367 final int first,
368 final int last,
369 final boolean propagateTransform,
370 final Rectangle fov,
371 final Filter< Patch > filter ) throws Exception
372 {
373 if ( !p.setup() ) return;
374
375 final List< Layer > layerRange = layerSet.getLayers( first, last );
376
377 Utils.log( layerRange.size() + "" );
378
379 Rectangle box = null;
380 for ( Iterator< Layer > it = layerRange.iterator(); it.hasNext(); )
381 {
382 /* remove empty layers */
383 final Layer la = it.next();
384 if ( !la.contains( Patch.class, true ) )
385 {
386 it.remove();
387 continue;
388 }
389
390 /* accumulate boxes */
391 if ( null == box ) // The first layer:
392 box = la.getMinimalBoundingBox( Patch.class, true );
393 else
394 box = box.union( la.getMinimalBoundingBox( Patch.class, true ) );
395 }
396
397 if ( fov != null )
398 box = box.intersection( fov );
399
400 if ( box.width <= 0 || box.height <= 0 )
401 {
402 Utils.log( "Bounding box empty." );
403 return;
404 }
405
406 if ( layerRange.size() == 0 )
407 {
408 Utils.log( "All layers in range are empty!" );
409 return;
410 }
411
412 /* do not work if there is only one layer selected */
413 if ( layerRange.size() < 2 )
414 {
415 Utils.log( "All except one layer in range are empty!" );
416 return;
417 }
418
419 final float scale = Math.min( 1.0f, Math.min( ( float )p.ppm.sift.maxOctaveSize / ( float )box.width, ( float )p.ppm.sift.maxOctaveSize / ( float )box.height ) );
420
421 /* extract and save features, overwrite cached files if requested */
422 try
423 {
424 extractAndSaveLayerFeatures( layerSet, layerRange, box, scale, filter, p.ppm.sift, p.ppm.clearCache );
425 }
426 catch ( Exception e )
427 {
428 return;
429 }
430
431 /* create tiles and models for all layers */
432 final ArrayList< Tile< ? > > tiles = new ArrayList< Tile< ? > >();
433 for ( int i = 0; i < layerRange.size(); ++i )
434 {
435 switch ( p.modelIndexOptimize )
436 {
437 case 0:
438 tiles.add( new Tile< TranslationModel2D >( new TranslationModel2D() ) );
439 break;
440 case 1:
441 tiles.add( new Tile< RigidModel2D >( new RigidModel2D() ) );
442 break;
443 case 2:
444 tiles.add( new Tile< SimilarityModel2D >( new SimilarityModel2D() ) );
445 break;
446 case 3:
447 tiles.add( new Tile< AffineModel2D >( new AffineModel2D() ) );
448 break;
449 case 4:
450 tiles.add( new Tile< HomographyModel2D >( new HomographyModel2D() ) );
451 break;
452 default:
453 return;
454 }
455 }
456
457 /* collect all pairs of slices for which a model could be found */
458 final ArrayList< Triple< Integer, Integer, AbstractModel< ? > > > pairs = new ArrayList< Triple< Integer, Integer, AbstractModel< ? > > >();
459
460 int numFailures = 0;
461
462 final float pointMatchScale = p.layerScale / scale;
463
464 for ( int i = 0; i < layerRange.size(); ++i )
465 {
466 final ArrayList< Thread > threads = new ArrayList< Thread >( p.maxNumThreads );
467
468 final int sliceA = i;
469 final Layer layerA = layerRange.get( i );
470 final int range = Math.min( layerRange.size(), i + p.maxNumNeighbors + 1 );
471
472 final String layerNameA = layerName( layerA );
473
474 J: for ( int j = i + 1; j < range; )
475 {
476 final int numThreads = Math.min( p.maxNumThreads, range - j );
477 final ArrayList< Triple< Integer, Integer, AbstractModel< ? > > > models =
478 new ArrayList< Triple< Integer, Integer, AbstractModel< ? > > >( numThreads );
479
480 for ( int k = 0; k < numThreads; ++k )
481 models.add( null );
482
483 for ( int t = 0; t < numThreads && j < range; ++t, ++j )
484 {
485 final int ti = t;
486 final int sliceB = j;
487 final Layer layerB = layerRange.get( j );
488
489 final String layerNameB = layerName( layerB );
490
491 final Thread thread = new Thread()
492 {
493 public void run()
494 {
495 IJ.showProgress( sliceA, layerRange.size() - 1 );
496
497 Utils.log( "matching " + layerNameB + " -> " + layerNameA + "..." );
498
499 ArrayList< PointMatch > candidates = null;
500 if ( !p.ppm.clearCache )
501 candidates = mpicbg.trakem2.align.Util.deserializePointMatches(
502 layerSet.getProject(), p.ppm, "layer", layerB.getId(), layerA.getId() );
503
504 if ( null == candidates )
505 {
506 ArrayList< Feature > fs1 = mpicbg.trakem2.align.Util.deserializeFeatures(
507 layerSet.getProject(), p.ppm.sift, "layer", layerA.getId() );
508 ArrayList< Feature > fs2 = mpicbg.trakem2.align.Util.deserializeFeatures(
509 layerSet.getProject(), p.ppm.sift, "layer", layerB.getId() );
510 candidates = new ArrayList< PointMatch >( FloatArray2DSIFT.createMatches( fs2, fs1, p.ppm.rod ) );
511
512 /* scale the candidates */
513 for ( final PointMatch pm : candidates )
514 {
515 final Point p1 = pm.getP1();
516 final Point p2 = pm.getP2();
517 final float[] l1 = p1.getL();
518 final float[] w1 = p1.getW();
519 final float[] l2 = p2.getL();
520 final float[] w2 = p2.getW();
521
522 l1[ 0 ] *= pointMatchScale;
523 l1[ 1 ] *= pointMatchScale;
524 w1[ 0 ] *= pointMatchScale;
525 w1[ 1 ] *= pointMatchScale;
526 l2[ 0 ] *= pointMatchScale;
527 l2[ 1 ] *= pointMatchScale;
528 w2[ 0 ] *= pointMatchScale;
529 w2[ 1 ] *= pointMatchScale;
530
531 }
532
533 if ( !mpicbg.trakem2.align.Util.serializePointMatches(
534 layerSet.getProject(), p.ppm, "layer", layerB.getId(), layerA.getId(), candidates ) )
535 Utils.log( "Could not store point match candidates for layers " + layerNameB + " and " + layerNameA + "." );
536 }
537
538 AbstractModel< ? > model;
539 switch ( p.modelIndex )
540 {
541 case 0:
542 model = new TranslationModel2D();
543 break;
544 case 1:
545 model = new RigidModel2D();
546 break;
547 case 2:
548 model = new SimilarityModel2D();
549 break;
550 case 3:
551 model = new AffineModel2D();
552 break;
553 case 4:
554 model = new HomographyModel2D();
555 break;
556 default:
557 return;
558 }
559
560 final ArrayList< PointMatch > inliers = new ArrayList< PointMatch >();
561
562 boolean modelFound;
563 boolean again = false;
564 try
565 {
566 do
567 {
568 again = false;
569 modelFound = model.filterRansac(
570 candidates,
571 inliers,
572 1000,
573 p.maxEpsilon * p.layerScale,
574 p.minInlierRatio,
575 p.minNumInliers,
576 3 );
577 if ( modelFound && p.rejectIdentity )
578 {
579 final ArrayList< Point > points = new ArrayList< Point >();
580 PointMatch.sourcePoints( inliers, points );
581 if ( Transforms.isIdentity( model, points, p.identityTolerance * p.layerScale ) )
582 {
583 IJ.log( "Identity transform for " + inliers.size() + " matches rejected." );
584 candidates.removeAll( inliers );
585 inliers.clear();
586 again = true;
587 }
588 }
589 }
590 while ( again );
591 }
592 catch ( NotEnoughDataPointsException e )
593 {
594 modelFound = false;
595 }
596
597 if ( modelFound )
598 {
599 Utils.log( layerNameB + " -> " + layerNameA + ": " + inliers.size() + " corresponding features with an average displacement of " + ( PointMatch.meanDistance( inliers ) / p.layerScale ) + "px identified." );
600 Utils.log( "Estimated transformation model: " + model );
601 models.set( ti, new Triple< Integer, Integer, AbstractModel< ? > >( sliceA, sliceB, model ) );
602 }
603 else
604 {
605 Utils.log( layerNameB + " -> " + layerNameA + ": no correspondences found." );
606 return;
607 }
608 }
609 };
610 threads.add( thread );
611 thread.start();
612 }
613
614 try
615 {
616 for ( final Thread thread : threads )
617 thread.join();
618 }
619 catch ( InterruptedException e )
620 {
621 Utils.log( "Establishing feature correspondences interrupted." );
622 for ( final Thread thread : threads )
623 thread.interrupt();
624 try
625 {
626 for ( final Thread thread : threads )
627 thread.join();
628 }
629 catch ( InterruptedException f ) {}
630 return;
631 }
632
633 threads.clear();
634
635 /* collect successfully matches pairs and break the search on gaps */
636 for ( int t = 0; t < models.size(); ++t )
637 {
638 final Triple< Integer, Integer, AbstractModel< ? > > pair = models.get( t );
639 if ( pair == null )
640 {
641 if ( ++numFailures > p.maxNumFailures )
642 break J;
643 }
644 else
645 {
646 numFailures = 0;
647 pairs.add( pair );
648 }
649 }
650 }
651 }
652
653 /* Elastic alignment */
654
655 /* Initialization */
656 final TileConfiguration initMeshes = new TileConfiguration();
657
658 final int meshWidth = ( int )Math.ceil( box.width * p.layerScale );
659 final int meshHeight = ( int )Math.ceil( box.height * p.layerScale );
660
661 final ArrayList< SpringMesh > meshes = new ArrayList< SpringMesh >( layerRange.size() );
662 for ( int i = 0; i < layerRange.size(); ++i )
663 meshes.add(
664 new SpringMesh(
665 p.resolutionSpringMesh,
666 meshWidth,
667 meshHeight,
668 p.stiffnessSpringMesh,
669 p.maxStretchSpringMesh * p.layerScale,
670 p.dampSpringMesh ) );
671
672 final int blockRadius = Math.max( 32, meshWidth / p.resolutionSpringMesh / 2 );
673
674 /** TODO set this something more than the largest error by the approximate model */
675 final int searchRadius = ( int )Math.round( p.layerScale * p.maxEpsilon );
676
677 for ( final Triple< Integer, Integer, AbstractModel< ? > > pair : pairs )
678 {
679 /* free memory */
680 layerSet.getProject().getLoader().releaseAll();
681
682 final SpringMesh m1 = meshes.get( pair.a );
683 final SpringMesh m2 = meshes.get( pair.b );
684
685 ArrayList< PointMatch > pm12 = new ArrayList< PointMatch >();
686 ArrayList< PointMatch > pm21 = new ArrayList< PointMatch >();
687
688 final Collection< Vertex > v1 = m1.getVertices();
689 final Collection< Vertex > v2 = m2.getVertices();
690
691 final Layer layer1 = layerRange.get( pair.a );
692 final Layer layer2 = layerRange.get( pair.b );
693
694 final Image img1 = layerSet.getProject().getLoader().getFlatAWTImage(
695 layer1,
696 box,
697 p.layerScale,
698 0xffffffff,
699 ImagePlus.COLOR_RGB,
700 Patch.class,
701 filterPatches( layer1, filter ),
702 true,
703 new Color( 0x00ffffff, true ) );
704
705 final Image img2 = layerSet.getProject().getLoader().getFlatAWTImage(
706 layer2,
707 box,
708 p.layerScale,
709 0xffffffff,
710 ImagePlus.COLOR_RGB,
711 Patch.class,
712 filterPatches( layer2, filter ),
713 true,
714 new Color( 0x00ffffff, true ) );
715
716 final int width = img1.getWidth( null );
717 final int height = img1.getHeight( null );
718
719 final FloatProcessor ip1 = new FloatProcessor( width, height );
720 final FloatProcessor ip2 = new FloatProcessor( width, height );
721 final FloatProcessor ip1Mask = new FloatProcessor( width, height );
722 final FloatProcessor ip2Mask = new FloatProcessor( width, height );
723
724 mpicbg.trakem2.align.Util.imageToFloatAndMask( img1, ip1, ip1Mask );
725 mpicbg.trakem2.align.Util.imageToFloatAndMask( img2, ip2, ip2Mask );
726
727 try
728 {
729 BlockMatching.matchByMaximalPMCC(
730 ip1,
731 ip2,
732 ip1Mask,
733 ip2Mask,
734 1.0f,
735 ( ( InvertibleCoordinateTransform )pair.c ).createInverse(),
736 blockRadius,
737 blockRadius,
738 searchRadius,
739 searchRadius,
740 p.minR,
741 p.rodR,
742 p.maxCurvatureR,
743 v1,
744 pm12,
745 new ErrorStatistic( 1 ) );
746 }
747 catch ( InterruptedException e )
748 {
749 Utils.log( "Block matching interrupted." );
750 IJ.showProgress( 1.0 );
751 return;
752 }
753 if ( Thread.interrupted() )
754 {
755 Utils.log( "Block matching interrupted." );
756 IJ.showProgress( 1.0 );
757 return;
758 }
759
760 Utils.log( pair.a + " > " + pair.b + ": found " + pm12.size() + " correspondences." );
761
762 /* <visualisation> */
763 // final List< Point > s1 = new ArrayList< Point >();
764 // PointMatch.sourcePoints( pm12, s1 );
765 // final ImagePlus imp1 = new ImagePlus( i + " >", ip1 );
766 // imp1.show();
767 // imp1.setOverlay( BlockMatching.illustrateMatches( pm12 ), Color.yellow, null );
768 // imp1.setRoi( Util.pointsToPointRoi( s1 ) );
769 // imp1.updateAndDraw();
770 /* </visualisation> */
771
772 try
773 {
774 BlockMatching.matchByMaximalPMCC(
775 ip2,
776 ip1,
777 ip2Mask,
778 ip1Mask,
779 1.0f,
780 pair.c,
781 blockRadius,
782 blockRadius,
783 searchRadius,
784 searchRadius,
785 p.minR,
786 p.rodR,
787 p.maxCurvatureR,
788 v2,
789 pm21,
790 new ErrorStatistic( 1 ) );
791 }
792 catch ( InterruptedException e )
793 {
794 Utils.log( "Block matching interrupted." );
795 IJ.showProgress( 1.0 );
796 return;
797 }
798 if ( Thread.interrupted() )
799 {
800 Utils.log( "Block matching interrupted." );
801 IJ.showProgress( 1.0 );
802 return;
803 }
804
805 IJ.log( pair.a + " < " + pair.b + ": found " + pm21.size() + " correspondences." );
806
807 /* <visualisation> */
808 // final List< Point > s2 = new ArrayList< Point >();
809 // PointMatch.sourcePoints( pm21, s2 );
810 // final ImagePlus imp2 = new ImagePlus( i + " <", ip2 );
811 // imp2.show();
812 // imp2.setOverlay( BlockMatching.illustrateMatches( pm21 ), Color.yellow, null );
813 // imp2.setRoi( Util.pointsToPointRoi( s2 ) );
814 // imp2.updateAndDraw();
815 /* </visualisation> */
816
817 final float springConstant = 1.0f / ( pair.b - pair.a );
818 IJ.log( pair.a + " <> " + pair.b + " spring constant = " + springConstant );
819
820 for ( final PointMatch pm : pm12 )
821 {
822 final Vertex p1 = ( Vertex )pm.getP1();
823 final Vertex p2 = new Vertex( pm.getP2() );
824 p1.addSpring( p2, new Spring( 0, springConstant ) );
825 m2.addPassiveVertex( p2 );
826 }
827
828 for ( final PointMatch pm : pm21 )
829 {
830 final Vertex p1 = ( Vertex )pm.getP1();
831 final Vertex p2 = new Vertex( pm.getP2() );
832 p1.addSpring( p2, new Spring( 0, springConstant ) );
833 m1.addPassiveVertex( p2 );
834 }
835
836 final Tile< ? > t1 = tiles.get( pair.a );
837 final Tile< ? > t2 = tiles.get( pair.b );
838
839 /*
840 * adding Tiles to the initialing TileConfiguration, adding a Tile
841 * multiple times does not harm because the TileConfiguration is
842 * backed by a Set.
843 */
844 if ( pm12.size() > pair.c.getMinNumMatches() )
845 {
846 initMeshes.addTile( t1 );
847 initMeshes.addTile( t2 );
848 t1.connect( t2, pm12 );
849 }
850 if ( pm21.size() > pair.c.getMinNumMatches() )
851 {
852 initMeshes.addTile( t1 );
853 initMeshes.addTile( t2 );
854 t2.connect( t1, pm21 );
855 }
856 }
857
858 /* pre-align by optimizing a piecewise linear model */
859 initMeshes.optimize(
860 p.maxEpsilon * p.layerScale,
861 p.maxIterationsSpringMesh,
862 p.maxPlateauwidthSpringMesh );
863 for ( int i = 0; i < layerRange.size(); ++i )
864 meshes.get( i ).init( tiles.get( i ).getModel() );
865
866 /* optimize the meshes */
867 try
868 {
869 long t0 = System.currentTimeMillis();
870 IJ.log("Optimizing spring meshes...");
871
872 SpringMesh.optimizeMeshes(
873 meshes,
874 p.maxEpsilon * p.layerScale,
875 p.maxIterationsSpringMesh,
876 p.maxPlateauwidthSpringMesh,
877 p.visualize );
878
879 IJ.log("Done optimizing spring meshes. Took " + (System.currentTimeMillis() - t0) + " ms");
880
881 }
882 catch ( NotEnoughDataPointsException e )
883 {
884 Utils.log( "There were not enough data points to get the spring mesh optimizing." );
885 e.printStackTrace();
886 return;
887 }
888
889 /* translate relative to bounding box */
890 for ( final SpringMesh mesh : meshes )
891 {
892 for ( final PointMatch pm : mesh.getVA().keySet() )
893 {
894 final Point p1 = pm.getP1();
895 final Point p2 = pm.getP2();
896 final float[] l = p1.getL();
897 final float[] w = p2.getW();
898 l[ 0 ] = l[ 0 ] / p.layerScale + box.x;
899 l[ 1 ] = l[ 1 ] / p.layerScale + box.y;
900 w[ 0 ] = w[ 0 ] / p.layerScale + box.x;
901 w[ 1 ] = w[ 1 ] / p.layerScale + box.y;
902 }
903 }
904
905 /* free memory */
906 layerSet.getProject().getLoader().releaseAll();
907
908 /* transfer layer transform into patch transforms and append to patches */
909 for ( int i = 0; i < layerRange.size(); ++i )
910 {
911 final Layer layer = layerRange.get( i );
912
913 final MovingLeastSquaresTransform2 mlt = new MovingLeastSquaresTransform2();
914 mlt.setModel( AffineModel2D.class );
915 mlt.setAlpha( 2.0f );
916 mlt.setMatches( meshes.get( i ).getVA().keySet() );
917
918 for ( final Patch patch : filterPatches( layer, filter ) )
919 mpicbg.trakem2.align.Util.applyLayerTransformToPatch( patch, mlt );
920 }
921
922 /* update patch mipmaps */
923 for ( final Layer layer : layerRange )
924 for ( final Patch patch : filterPatches( layer, filter ) )
925 patch.updateMipMaps();
926
927 Utils.log( "Done." );
928 }
929 }
register/edit/analyze/3D model TB of images