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