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