?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
forgot to set isAligned flag in elastic layer alignment parameter
[trakem2.git] / TrakEM2_ / src / main / java / mpicbg / trakem2 / align / ElasticLayerAlignment.java
blob1694e91dbe1dc7bde35fd1d7eac731f34ba9172d
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 ini.trakem2.Project;
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.HashSet;
38 import java.util.Iterator;
39 import java.util.List;
40 import java.util.Set;
41 import java.util.concurrent.atomic.AtomicInteger;
42
43 import mpicbg.ij.blockmatching.BlockMatching;
44 import mpicbg.imagefeatures.Feature;
45 import mpicbg.imagefeatures.FloatArray2DSIFT;
46 import mpicbg.models.AbstractModel;
47 import mpicbg.models.AffineModel2D;
48 import mpicbg.models.ErrorStatistic;
49 import mpicbg.models.HomographyModel2D;
50 import mpicbg.models.InvertibleCoordinateTransform;
51 import mpicbg.models.NotEnoughDataPointsException;
52 import mpicbg.models.Point;
53 import mpicbg.models.PointMatch;
54 import mpicbg.models.RigidModel2D;
55 import mpicbg.models.SimilarityModel2D;
56 import mpicbg.models.Spring;
57 import mpicbg.models.SpringMesh;
58 import mpicbg.models.Tile;
59 import mpicbg.models.TileConfiguration;
60 import mpicbg.models.Transforms;
61 import mpicbg.models.TranslationModel2D;
62 import mpicbg.models.Vertex;
63 import mpicbg.trakem2.transform.MovingLeastSquaresTransform2;
64 import mpicbg.trakem2.util.Triple;
65
66 /**
67 * @author Stephan Saalfeld <saalfeld@mpi-cbg.de>
68 */
69 public class ElasticLayerAlignment
70 {
71 final static public class Param extends AbstractLayerAlignmentParam implements Serializable
72 {
73 private static final long serialVersionUID = 398966126705836033L;
74
75 public boolean isAligned = false;
76
77 public float layerScale = 0.1f;
78 public float minR = 0.6f;
79 public float maxCurvatureR = 10f;
80 public float rodR = 0.9f;
81 public int searchRadius = 200;
82 public int blockRadius = -1;
83
84 public boolean useLocalSmoothnessFilter = true;
85 public int localModelIndex = 1;
86 public float localRegionSigma = searchRadius;
87 public float maxLocalEpsilon = searchRadius / 2;
88 public float maxLocalTrust = 3;
89
90 public int resolutionSpringMesh = 16;
91 public float stiffnessSpringMesh = 0.1f;
92 public float dampSpringMesh = 0.6f;
93 public float maxStretchSpringMesh = 2000.0f;
94 public int maxIterationsSpringMesh = 1000;
95 public int maxPlateauwidthSpringMesh = 200;
96
97 public boolean setup( final Rectangle box )
98 {
99 /* Block Matching */
100 if ( blockRadius < 0 )
101 {
102 blockRadius = box.width / resolutionSpringMesh / 2;
103 }
104 final GenericDialog gdBlockMatching = new GenericDialog( "Elastically align layers: Block Matching parameters" );
105
106 gdBlockMatching.addMessage( "Block Matching:" );
107 /* TODO suggest isotropic resolution for this parameter */
108 gdBlockMatching.addNumericField( "layer_scale :", layerScale, 2 );
109 gdBlockMatching.addNumericField( "search_radius :", searchRadius, 0, 6, "px" );
110 gdBlockMatching.addNumericField( "block_radius :", blockRadius, 0, 6, "px" );
111 /* TODO suggest a resolution that matches searchRadius */
112 gdBlockMatching.addNumericField( "resolution :", resolutionSpringMesh, 0 );
113
114 gdBlockMatching.addMessage( "Correlation Filters:" );
115 gdBlockMatching.addNumericField( "minimal_PMCC_r :", minR, 2 );
116 gdBlockMatching.addNumericField( "maximal_curvature_ratio :", maxCurvatureR, 2 );
117 gdBlockMatching.addNumericField( "maximal_second_best_r/best_r :", rodR, 2 );
118
119 gdBlockMatching.addMessage( "Local Smoothness Filter:" );
120 gdBlockMatching.addCheckbox( "use_local_smoothness_filter", useLocalSmoothnessFilter );
121 gdBlockMatching.addChoice( "approximate_local_transformation :", Param.modelStrings, Param.modelStrings[ localModelIndex ] );
122 gdBlockMatching.addNumericField( "local_region_sigma:", localRegionSigma, 2, 6, "px" );
123 gdBlockMatching.addNumericField( "maximal_local_displacement (absolute):", maxLocalEpsilon, 2, 6, "px" );
124 gdBlockMatching.addNumericField( "maximal_local_displacement (relative):", maxLocalTrust, 2 );
125
126 gdBlockMatching.addMessage( "Miscellaneous:" );
127 gdBlockMatching.addCheckbox( "layers_are_pre-aligned", isAligned );
128 gdBlockMatching.addNumericField( "test_maximally :", maxNumNeighbors, 0, 6, "layers" );
129
130 gdBlockMatching.showDialog();
131
132 if ( gdBlockMatching.wasCanceled() )
133 return false;
134
135 layerScale = ( float )gdBlockMatching.getNextNumber();
136 searchRadius = ( int )gdBlockMatching.getNextNumber();
137 blockRadius = ( int )gdBlockMatching.getNextNumber();
138 resolutionSpringMesh = ( int )gdBlockMatching.getNextNumber();
139 minR = ( float )gdBlockMatching.getNextNumber();
140 maxCurvatureR = ( float )gdBlockMatching.getNextNumber();
141 rodR = ( float )gdBlockMatching.getNextNumber();
142 useLocalSmoothnessFilter = gdBlockMatching.getNextBoolean();
143 localModelIndex = gdBlockMatching.getNextChoiceIndex();
144 localRegionSigma = ( float )gdBlockMatching.getNextNumber();
145 maxLocalEpsilon = ( float )gdBlockMatching.getNextNumber();
146 maxLocalTrust = ( float )gdBlockMatching.getNextNumber();
147 isAligned = gdBlockMatching.getNextBoolean();
148 maxNumNeighbors = ( int )gdBlockMatching.getNextNumber();
149
150
151 if ( !isAligned )
152 {
153 if ( !setupSIFT( "Elastically align layers: " ) )
154 return false;
155
156 /* Geometric filters */
157
158 final GenericDialog gd = new GenericDialog( "Elastically align layers: Geometric filters" );
159
160 gd.addNumericField( "maximal_alignment_error :", maxEpsilon, 2, 6, "px" );
161 gd.addNumericField( "minimal_inlier_ratio :", minInlierRatio, 2 );
162 gd.addNumericField( "minimal_number_of_inliers :", minNumInliers, 0 );
163 gd.addChoice( "approximate_transformation :", Param.modelStrings, Param.modelStrings[ expectedModelIndex ] );
164 gd.addCheckbox( "ignore constant background", rejectIdentity );
165 gd.addNumericField( "tolerance :", identityTolerance, 2, 6, "px" );
166 gd.addNumericField( "give_up_after :", maxNumFailures, 0, 6, "failures" );
167
168 gd.showDialog();
169
170 if ( gd.wasCanceled() )
171 return false;
172
173 maxEpsilon = ( float )gd.getNextNumber();
174 minInlierRatio = ( float )gd.getNextNumber();
175 minNumInliers = ( int )gd.getNextNumber();
176 expectedModelIndex = gd.getNextChoiceIndex();
177 rejectIdentity = gd.getNextBoolean();
178 identityTolerance = ( float )gd.getNextNumber();
179 maxNumFailures = ( int )gd.getNextNumber();
180 }
181
182
183 /* Optimization */
184 final GenericDialog gdOptimize = new GenericDialog( "Elastically align layers: Optimization" );
185
186 gdOptimize.addMessage( "Approximate Optimizer:" );
187 gdOptimize.addChoice( "approximate_transformation :", Param.modelStrings, Param.modelStrings[ desiredModelIndex ] );
188 gdOptimize.addNumericField( "maximal_iterations :", maxIterationsOptimize, 0 );
189 gdOptimize.addNumericField( "maximal_plateauwidth :", maxPlateauwidthOptimize, 0 );
190
191 gdOptimize.addMessage( "Spring Mesh:" );
192 gdOptimize.addNumericField( "stiffness :", stiffnessSpringMesh, 2 );
193 gdOptimize.addNumericField( "maximal_stretch :", maxStretchSpringMesh, 2, 6, "px" );
194 gdOptimize.addNumericField( "maximal_iterations :", maxIterationsSpringMesh, 0 );
195 gdOptimize.addNumericField( "maximal_plateauwidth :", maxPlateauwidthSpringMesh, 0 );
196
197 gdOptimize.showDialog();
198
199 if ( gdOptimize.wasCanceled() )
200 return false;
201
202 desiredModelIndex = gdOptimize.getNextChoiceIndex();
203 maxIterationsOptimize = ( int )gdOptimize.getNextNumber();
204 maxPlateauwidthOptimize = ( int )gdOptimize.getNextNumber();
205
206 stiffnessSpringMesh = ( float )gdOptimize.getNextNumber();
207 maxStretchSpringMesh = ( float )gdOptimize.getNextNumber();
208 maxIterationsSpringMesh = ( int )gdOptimize.getNextNumber();
209 maxPlateauwidthSpringMesh = ( int )gdOptimize.getNextNumber();
210
211 return true;
212 }
213
214 public Param() {}
215
216 public Param(
217 final int SIFTfdBins,
218 final int SIFTfdSize,
219 final float SIFTinitialSigma,
220 final int SIFTmaxOctaveSize,
221 final int SIFTminOctaveSize,
222 final int SIFTsteps,
223
224 final boolean clearCache,
225 final int maxNumThreadsSift,
226 final float rod,
227
228 final int desiredModelIndex,
229 final int expectedModelIndex,
230 final float identityTolerance,
231 final boolean isAligned,
232 final float maxEpsilon,
233 final int maxIterationsOptimize,
234 final int maxNumFailures,
235 final int maxNumNeighbors,
236 final int maxNumThreads,
237 final int maxPlateauwidthOptimize,
238 final float minInlierRatio,
239 final int minNumInliers,
240 final boolean multipleHypotheses,
241 final boolean rejectIdentity,
242 final boolean visualize,
243
244 final int blockRadius,
245 final float dampSpringMesh,
246 final float layerScale,
247 final int localModelIndex,
248 final float localRegionSigma,
249 final float maxCurvatureR,
250 final int maxIterationsSpringMesh,
251 final float maxLocalEpsilon,
252 final float maxLocalTrust,
253 final int maxPlateauwidthSpringMesh,
254 final float maxStretchSpringMesh,
255 final float minR,
256 final int resolutionSpringMesh,
257 final float rodR,
258 final int searchRadius,
259 final float stiffnessSpringMesh,
260 final boolean useLocalSmoothnessFilter )
261 {
262 super(
263 SIFTfdBins,
264 SIFTfdSize,
265 SIFTinitialSigma,
266 SIFTmaxOctaveSize,
267 SIFTminOctaveSize,
268 SIFTsteps,
269 clearCache,
270 maxNumThreadsSift,
271 rod,
272 desiredModelIndex,
273 expectedModelIndex,
274 identityTolerance,
275 maxEpsilon,
276 maxIterationsOptimize,
277 maxNumFailures,
278 maxNumNeighbors,
279 maxNumThreads,
280 maxPlateauwidthOptimize,
281 minInlierRatio,
282 minNumInliers,
283 multipleHypotheses,
284 rejectIdentity,
285 visualize );
286
287 this.isAligned = isAligned;
288 this.blockRadius = blockRadius;
289 this.dampSpringMesh = dampSpringMesh;
290 this.layerScale = layerScale;
291 this.localModelIndex = localModelIndex;
292 this.localRegionSigma = localRegionSigma;
293 this.maxCurvatureR = maxCurvatureR;
294 this.maxIterationsSpringMesh = maxIterationsSpringMesh;
295 this.maxLocalEpsilon = maxLocalEpsilon;
296 this.maxLocalTrust = maxLocalTrust;
297 this.maxPlateauwidthSpringMesh = maxPlateauwidthSpringMesh;
298 this.maxStretchSpringMesh = maxStretchSpringMesh;
299 this.minR = minR;
300 this.resolutionSpringMesh = resolutionSpringMesh;
301 this.rodR = rodR;
302 this.searchRadius = searchRadius;
303 this.stiffnessSpringMesh = stiffnessSpringMesh;
304 this.useLocalSmoothnessFilter = useLocalSmoothnessFilter;
305 }
306
307 @Override
308 public Param clone()
309 {
310 return new Param(
311 ppm.sift.fdBins,
312 ppm.sift.fdSize,
313 ppm.sift.initialSigma,
314 ppm.sift.maxOctaveSize,
315 ppm.sift.minOctaveSize,
316 ppm.sift.steps,
317
318 ppm.clearCache,
319 ppm.maxNumThreadsSift,
320 ppm.rod,
321
322 desiredModelIndex,
323 expectedModelIndex,
324 identityTolerance,
325 isAligned,
326 maxEpsilon,
327 maxIterationsOptimize,
328 maxNumFailures,
329 maxNumNeighbors,
330 maxNumThreads,
331 maxPlateauwidthOptimize,
332 minInlierRatio,
333 minNumInliers,
334 multipleHypotheses,
335 rejectIdentity,
336 visualize,
337
338 blockRadius,
339 dampSpringMesh,
340 layerScale,
341 localModelIndex,
342 localRegionSigma,
343 maxCurvatureR,
344 maxIterationsSpringMesh,
345 maxLocalEpsilon,
346 maxLocalTrust,
347 maxPlateauwidthSpringMesh,
348 maxStretchSpringMesh,
349 minR,
350 resolutionSpringMesh,
351 rodR,
352 searchRadius,
353 stiffnessSpringMesh,
354 useLocalSmoothnessFilter );
355 }
356 }
357
358 final static Param p = new Param();
359
360
361 final static private String layerName( final Layer layer )
362 {
363 return new StringBuffer( "layer z=" )
364 .append( String.format( "%.3f", layer.getZ() ) )
365 .append( " `" )
366 .append( layer.getTitle() )
367 .append( "'" )
368 .toString();
369
370 }
371
372
373 /**
374 *
375 * @param param
376 * @param layerRange
377 * @param fixedLayers
378 * @param emptyLayers
379 * @param box
380 * @param propagateTransform
381 * @param fov
382 * @param filter
383 * @throws Exception
384 */
385 final public void exec(
386 final Param param,
387 final Project project,
388 final List< Layer > layerRange,
389 final Set< Layer > fixedLayers,
390 final Set< Layer > emptyLayers,
391 final Rectangle box,
392 final boolean propagateTransformBefore,
393 final boolean propagateTransformAfter,
394 final Filter< Patch > filter ) throws Exception
395 {
396 final double scale = Math.min( 1.0, Math.min( ( double )param.ppm.sift.maxOctaveSize / ( double )box.width, ( double )param.ppm.sift.maxOctaveSize / ( double )box.height ) );
397
398
399 /* create tiles and models for all layers */
400 final ArrayList< Tile< ? > > tiles = new ArrayList< Tile< ? > >();
401 for ( int i = 0; i < layerRange.size(); ++i )
402 {
403 switch ( param.desiredModelIndex )
404 {
405 case 0:
406 tiles.add( new Tile< TranslationModel2D >( new TranslationModel2D() ) );
407 break;
408 case 1:
409 tiles.add( new Tile< RigidModel2D >( new RigidModel2D() ) );
410 break;
411 case 2:
412 tiles.add( new Tile< SimilarityModel2D >( new SimilarityModel2D() ) );
413 break;
414 case 3:
415 tiles.add( new Tile< AffineModel2D >( new AffineModel2D() ) );
416 break;
417 case 4:
418 tiles.add( new Tile< HomographyModel2D >( new HomographyModel2D() ) );
419 break;
420 default:
421 return;
422 }
423 }
424
425 /* collect all pairs of slices for which a model could be found */
426 final ArrayList< Triple< Integer, Integer, AbstractModel< ? > > > pairs = new ArrayList< Triple< Integer, Integer, AbstractModel< ? > > >();
427
428
429 if ( !param.isAligned )
430 {
431
432 /* extract and save features, overwrite cached files if requested */
433 try
434 {
435 AlignmentUtils.extractAndSaveLayerFeatures( layerRange, box, scale, filter, param.ppm.sift, param.ppm.clearCache, param.ppm.maxNumThreadsSift );
436 }
437 catch ( final Exception e )
438 {
439 return;
440 }
441
442 /* match and filter feature correspondences */
443 int numFailures = 0;
444
445 final double pointMatchScale = param.layerScale / scale;
446
447 for ( int i = 0; i < layerRange.size(); ++i )
448 {
449 final ArrayList< Thread > threads = new ArrayList< Thread >( param.maxNumThreads );
450
451 final int sliceA = i;
452 final Layer layerA = layerRange.get( i );
453 final int range = Math.min( layerRange.size(), i + param.maxNumNeighbors + 1 );
454
455 final String layerNameA = layerName( layerA );
456
457 J: for ( int j = i + 1; j < range; )
458 {
459 final int numThreads = Math.min( param.maxNumThreads, range - j );
460 final ArrayList< Triple< Integer, Integer, AbstractModel< ? > > > models =
461 new ArrayList< Triple< Integer, Integer, AbstractModel< ? > > >( numThreads );
462
463 for ( int k = 0; k < numThreads; ++k )
464 models.add( null );
465
466 for ( int t = 0; t < numThreads && j < range; ++t, ++j )
467 {
468 final int ti = t;
469 final int sliceB = j;
470 final Layer layerB = layerRange.get( j );
471
472 final String layerNameB = layerName( layerB );
473
474 final Thread thread = new Thread()
475 {
476 @Override
477 public void run()
478 {
479 IJ.showProgress( sliceA, layerRange.size() - 1 );
480
481 Utils.log( "matching " + layerNameB + " -> " + layerNameA + "..." );
482
483 ArrayList< PointMatch > candidates = null;
484 if ( !param.ppm.clearCache )
485 candidates = mpicbg.trakem2.align.Util.deserializePointMatches(
486 project, param.ppm, "layer", layerB.getId(), layerA.getId() );
487
488 if ( null == candidates )
489 {
490 final ArrayList< Feature > fs1 = mpicbg.trakem2.align.Util.deserializeFeatures(
491 project, param.ppm.sift, "layer", layerA.getId() );
492 final ArrayList< Feature > fs2 = mpicbg.trakem2.align.Util.deserializeFeatures(
493 project, param.ppm.sift, "layer", layerB.getId() );
494 candidates = new ArrayList< PointMatch >( FloatArray2DSIFT.createMatches( fs2, fs1, param.ppm.rod ) );
495
496 /* scale the candidates */
497 for ( final PointMatch pm : candidates )
498 {
499 final Point p1 = pm.getP1();
500 final Point p2 = pm.getP2();
501 final float[] l1 = p1.getL();
502 final float[] w1 = p1.getW();
503 final float[] l2 = p2.getL();
504 final float[] w2 = p2.getW();
505
506 l1[ 0 ] *= pointMatchScale;
507 l1[ 1 ] *= pointMatchScale;
508 w1[ 0 ] *= pointMatchScale;
509 w1[ 1 ] *= pointMatchScale;
510 l2[ 0 ] *= pointMatchScale;
511 l2[ 1 ] *= pointMatchScale;
512 w2[ 0 ] *= pointMatchScale;
513 w2[ 1 ] *= pointMatchScale;
514
515 }
516
517 if ( !mpicbg.trakem2.align.Util.serializePointMatches(
518 project, param.ppm, "layer", layerB.getId(), layerA.getId(), candidates ) )
519 Utils.log( "Could not store point match candidates for layers " + layerNameB + " and " + layerNameA + "." );
520 }
521
522 AbstractModel< ? > model;
523 switch ( param.expectedModelIndex )
524 {
525 case 0:
526 model = new TranslationModel2D();
527 break;
528 case 1:
529 model = new RigidModel2D();
530 break;
531 case 2:
532 model = new SimilarityModel2D();
533 break;
534 case 3:
535 model = new AffineModel2D();
536 break;
537 case 4:
538 model = new HomographyModel2D();
539 break;
540 default:
541 return;
542 }
543
544 final ArrayList< PointMatch > inliers = new ArrayList< PointMatch >();
545
546 boolean modelFound;
547 boolean again = false;
548 try
549 {
550 do
551 {
552 again = false;
553 modelFound = model.filterRansac(
554 candidates,
555 inliers,
556 1000,
557 param.maxEpsilon * param.layerScale,
558 param.minInlierRatio,
559 param.minNumInliers,
560 3 );
561 if ( modelFound && param.rejectIdentity )
562 {
563 final ArrayList< Point > points = new ArrayList< Point >();
564 PointMatch.sourcePoints( inliers, points );
565 if ( Transforms.isIdentity( model, points, param.identityTolerance * param.layerScale ) )
566 {
567 IJ.log( "Identity transform for " + inliers.size() + " matches rejected." );
568 candidates.removeAll( inliers );
569 inliers.clear();
570 again = true;
571 }
572 }
573 }
574 while ( again );
575 }
576 catch ( final NotEnoughDataPointsException e )
577 {
578 modelFound = false;
579 }
580
581 if ( modelFound )
582 {
583 Utils.log( layerNameB + " -> " + layerNameA + ": " + inliers.size() + " corresponding features with an average displacement of " + ( PointMatch.meanDistance( inliers ) / param.layerScale ) + "px identified." );
584 Utils.log( "Estimated transformation model: " + model );
585 models.set( ti, new Triple< Integer, Integer, AbstractModel< ? > >( sliceA, sliceB, model ) );
586 }
587 else
588 {
589 Utils.log( layerNameB + " -> " + layerNameA + ": no correspondences found." );
590 return;
591 }
592 }
593 };
594 threads.add( thread );
595 thread.start();
596 }
597
598 try
599 {
600 for ( final Thread thread : threads )
601 thread.join();
602 }
603 catch ( final InterruptedException e )
604 {
605 Utils.log( "Establishing feature correspondences interrupted." );
606 for ( final Thread thread : threads )
607 thread.interrupt();
608 try
609 {
610 for ( final Thread thread : threads )
611 thread.join();
612 }
613 catch ( final InterruptedException f ) {}
614 return;
615 }
616
617 threads.clear();
618
619 /* collect successfully matches pairs and break the search on gaps */
620 for ( int t = 0; t < models.size(); ++t )
621 {
622 final Triple< Integer, Integer, AbstractModel< ? > > pair = models.get( t );
623 if ( pair == null )
624 {
625 if ( ++numFailures > param.maxNumFailures )
626 break J;
627 }
628 else
629 {
630 numFailures = 0;
631 pairs.add( pair );
632 }
633 }
634 }
635 }
636 }
637 else
638 {
639 for ( int i = 0; i < layerRange.size(); ++i )
640 {
641 final int range = Math.min( layerRange.size(), i + param.maxNumNeighbors + 1 );
642
643 for ( int j = i + 1; j < range; ++j )
644 {
645 pairs.add( new Triple< Integer, Integer, AbstractModel< ? > >( i, j, new TranslationModel2D() ) );
646 }
647 }
648 }
649
650 /* Elastic alignment */
651
652 /* Initialization */
653 final TileConfiguration initMeshes = new TileConfiguration();
654
655 final int meshWidth = ( int )Math.ceil( box.width * param.layerScale );
656 final int meshHeight = ( int )Math.ceil( box.height * param.layerScale );
657
658 final ArrayList< SpringMesh > meshes = new ArrayList< SpringMesh >( layerRange.size() );
659 for ( int i = 0; i < layerRange.size(); ++i )
660 meshes.add(
661 new SpringMesh(
662 param.resolutionSpringMesh,
663 meshWidth,
664 meshHeight,
665 param.stiffnessSpringMesh,
666 param.maxStretchSpringMesh * param.layerScale,
667 param.dampSpringMesh ) );
668
669 //final int blockRadius = Math.max( 32, meshWidth / p.resolutionSpringMesh / 2 );
670 final int blockRadius = Math.max( 16, mpicbg.util.Util.roundPos( param.layerScale * param.blockRadius ) );
671
672 Utils.log( "effective block radius = " + blockRadius );
673
674 /* scale pixel distances */
675 final int searchRadius = ( int )Math.round( param.layerScale * param.searchRadius );
676 final float localRegionSigma = param.layerScale * param.localRegionSigma;
677 final float maxLocalEpsilon = param.layerScale * param.maxLocalEpsilon;
678
679 final AbstractModel< ? > localSmoothnessFilterModel = Util.createModel( param.localModelIndex );
680
681
682 for ( final Triple< Integer, Integer, AbstractModel< ? > > pair : pairs )
683 {
684 /* free memory */
685 project.getLoader().releaseAll();
686
687 final SpringMesh m1 = meshes.get( pair.a );
688 final SpringMesh m2 = meshes.get( pair.b );
689
690 final ArrayList< PointMatch > pm12 = new ArrayList< PointMatch >();
691 final ArrayList< PointMatch > pm21 = new ArrayList< PointMatch >();
692
693 final Collection< Vertex > v1 = m1.getVertices();
694 final Collection< Vertex > v2 = m2.getVertices();
695
696 final Layer layer1 = layerRange.get( pair.a );
697 final Layer layer2 = layerRange.get( pair.b );
698
699 final boolean layer1Fixed = fixedLayers.contains( layer1 );
700 final boolean layer2Fixed = fixedLayers.contains( layer2 );
701
702 final Tile< ? > t1 = tiles.get( pair.a );
703 final Tile< ? > t2 = tiles.get( pair.b );
704
705 if ( !( layer1Fixed && layer2Fixed ) )
706 {
707 final Image img1 = project.getLoader().getFlatAWTImage(
708 layer1,
709 box,
710 param.layerScale,
711 0xffffffff,
712 ImagePlus.COLOR_RGB,
713 Patch.class,
714 AlignmentUtils.filterPatches( layer1, filter ),
715 true,
716 new Color( 0x00ffffff, true ) );
717
718 final Image img2 = project.getLoader().getFlatAWTImage(
719 layer2,
720 box,
721 param.layerScale,
722 0xffffffff,
723 ImagePlus.COLOR_RGB,
724 Patch.class,
725 AlignmentUtils.filterPatches( layer2, filter ),
726 true,
727 new Color( 0x00ffffff, true ) );
728
729 final int width = img1.getWidth( null );
730 final int height = img1.getHeight( null );
731
732 final FloatProcessor ip1 = new FloatProcessor( width, height );
733 final FloatProcessor ip2 = new FloatProcessor( width, height );
734 final FloatProcessor ip1Mask = new FloatProcessor( width, height );
735 final FloatProcessor ip2Mask = new FloatProcessor( width, height );
736
737 mpicbg.trakem2.align.Util.imageToFloatAndMask( img1, ip1, ip1Mask );
738 mpicbg.trakem2.align.Util.imageToFloatAndMask( img2, ip2, ip2Mask );
739
740 final float springConstant = 1.0f / ( pair.b - pair.a );
741
742 if ( layer1Fixed )
743 initMeshes.fixTile( t1 );
744 else
745 {
746 try
747 {
748 BlockMatching.matchByMaximalPMCC(
749 ip1,
750 ip2,
751 ip1Mask,
752 ip2Mask,
753 1.0f,
754 ( ( InvertibleCoordinateTransform )pair.c ).createInverse(),
755 blockRadius,
756 blockRadius,
757 searchRadius,
758 searchRadius,
759 param.minR,
760 param.rodR,
761 param.maxCurvatureR,
762 v1,
763 pm12,
764 new ErrorStatistic( 1 ) );
765 }
766 catch ( final InterruptedException e )
767 {
768 Utils.log( "Block matching interrupted." );
769 IJ.showProgress( 1.0 );
770 return;
771 }
772 if ( Thread.interrupted() )
773 {
774 Utils.log( "Block matching interrupted." );
775 IJ.showProgress( 1.0 );
776 return;
777 }
778
779 if ( param.useLocalSmoothnessFilter )
780 {
781 Utils.log( pair.a + " > " + pair.b + ": found " + pm12.size() + " correspondence candidates." );
782 localSmoothnessFilterModel.localSmoothnessFilter( pm12, pm12, localRegionSigma, maxLocalEpsilon, param.maxLocalTrust );
783 Utils.log( pair.a + " > " + pair.b + ": " + pm12.size() + " candidates passed local smoothness filter." );
784 }
785 else
786 {
787 Utils.log( pair.a + " > " + pair.b + ": found " + pm12.size() + " correspondences." );
788 }
789
790 /* <visualisation> */
791 // final List< Point > s1 = new ArrayList< Point >();
792 // PointMatch.sourcePoints( pm12, s1 );
793 // final ImagePlus imp1 = new ImagePlus( i + " >", ip1 );
794 // imp1.show();
795 // imp1.setOverlay( BlockMatching.illustrateMatches( pm12 ), Color.yellow, null );
796 // imp1.setRoi( Util.pointsToPointRoi( s1 ) );
797 // imp1.updateAndDraw();
798 /* </visualisation> */
799
800 for ( final PointMatch pm : pm12 )
801 {
802 final Vertex p1 = ( Vertex )pm.getP1();
803 final Vertex p2 = new Vertex( pm.getP2() );
804 p1.addSpring( p2, new Spring( 0, springConstant ) );
805 m2.addPassiveVertex( p2 );
806 }
807
808 /*
809 * adding Tiles to the initialing TileConfiguration, adding a Tile
810 * multiple times does not harm because the TileConfiguration is
811 * backed by a Set.
812 */
813 if ( pm12.size() > pair.c.getMinNumMatches() )
814 {
815 initMeshes.addTile( t1 );
816 initMeshes.addTile( t2 );
817 t1.connect( t2, pm12 );
818 }
819 }
820
821 if ( layer2Fixed )
822 initMeshes.fixTile( t2 );
823 else
824 {
825 try
826 {
827 BlockMatching.matchByMaximalPMCC(
828 ip2,
829 ip1,
830 ip2Mask,
831 ip1Mask,
832 1.0f,
833 pair.c,
834 blockRadius,
835 blockRadius,
836 searchRadius,
837 searchRadius,
838 param.minR,
839 param.rodR,
840 param.maxCurvatureR,
841 v2,
842 pm21,
843 new ErrorStatistic( 1 ) );
844 }
845 catch ( final InterruptedException e )
846 {
847 Utils.log( "Block matching interrupted." );
848 IJ.showProgress( 1.0 );
849 return;
850 }
851 if ( Thread.interrupted() )
852 {
853 Utils.log( "Block matching interrupted." );
854 IJ.showProgress( 1.0 );
855 return;
856 }
857
858 if ( param.useLocalSmoothnessFilter )
859 {
860 Utils.log( pair.a + " < " + pair.b + ": found " + pm21.size() + " correspondence candidates." );
861 localSmoothnessFilterModel.localSmoothnessFilter( pm21, pm21, localRegionSigma, maxLocalEpsilon, param.maxLocalTrust );
862 Utils.log( pair.a + " < " + pair.b + ": " + pm21.size() + " candidates passed local smoothness filter." );
863 }
864 else
865 {
866 Utils.log( pair.a + " < " + pair.b + ": found " + pm21.size() + " correspondences." );
867 }
868
869 /* <visualisation> */
870 // final List< Point > s2 = new ArrayList< Point >();
871 // PointMatch.sourcePoints( pm21, s2 );
872 // final ImagePlus imp2 = new ImagePlus( i + " <", ip2 );
873 // imp2.show();
874 // imp2.setOverlay( BlockMatching.illustrateMatches( pm21 ), Color.yellow, null );
875 // imp2.setRoi( Util.pointsToPointRoi( s2 ) );
876 // imp2.updateAndDraw();
877 /* </visualisation> */
878
879 for ( final PointMatch pm : pm21 )
880 {
881 final Vertex p1 = ( Vertex )pm.getP1();
882 final Vertex p2 = new Vertex( pm.getP2() );
883 p1.addSpring( p2, new Spring( 0, springConstant ) );
884 m1.addPassiveVertex( p2 );
885 }
886
887 /*
888 * adding Tiles to the initialing TileConfiguration, adding a Tile
889 * multiple times does not harm because the TileConfiguration is
890 * backed by a Set.
891 */
892 if ( pm21.size() > pair.c.getMinNumMatches() )
893 {
894 initMeshes.addTile( t1 );
895 initMeshes.addTile( t2 );
896 t2.connect( t1, pm21 );
897 }
898 }
899
900 Utils.log( pair.a + " <> " + pair.b + " spring constant = " + springConstant );
901 }
902 }
903
904 /* pre-align by optimizing a piecewise linear model */
905 initMeshes.optimize(
906 param.maxEpsilon * param.layerScale,
907 param.maxIterationsSpringMesh,
908 param.maxPlateauwidthSpringMesh );
909 for ( int i = 0; i < layerRange.size(); ++i )
910 meshes.get( i ).init( tiles.get( i ).getModel() );
911
912 /* optimize the meshes */
913 try
914 {
915 final long t0 = System.currentTimeMillis();
916 Utils.log("Optimizing spring meshes...");
917
918 SpringMesh.optimizeMeshes(
919 meshes,
920 param.maxEpsilon * param.layerScale,
921 param.maxIterationsSpringMesh,
922 param.maxPlateauwidthSpringMesh,
923 param.visualize );
924
925 Utils.log("Done optimizing spring meshes. Took " + (System.currentTimeMillis() - t0) + " ms");
926
927 }
928 catch ( final NotEnoughDataPointsException e )
929 {
930 Utils.log( "There were not enough data points to get the spring mesh optimizing." );
931 e.printStackTrace();
932 return;
933 }
934
935 /* translate relative to bounding box */
936 for ( final SpringMesh mesh : meshes )
937 {
938 for ( final PointMatch pm : mesh.getVA().keySet() )
939 {
940 final Point p1 = pm.getP1();
941 final Point p2 = pm.getP2();
942 final float[] l = p1.getL();
943 final float[] w = p2.getW();
944 l[ 0 ] = l[ 0 ] / param.layerScale + box.x;
945 l[ 1 ] = l[ 1 ] / param.layerScale + box.y;
946 w[ 0 ] = w[ 0 ] / param.layerScale + box.x;
947 w[ 1 ] = w[ 1 ] / param.layerScale + box.y;
948 }
949 }
950
951 /* free memory */
952 project.getLoader().releaseAll();
953
954 final Layer first = layerRange.get( 0 );
955 final List< Layer > layers = first.getParent().getLayers();
956
957 /* transfer layer transform into patch transforms and append to patches */
958 if ( propagateTransformBefore || propagateTransformAfter )
959 {
960 if ( propagateTransformBefore )
961 {
962 final MovingLeastSquaresTransform2 mlt = makeMLST2( meshes.get( 0 ).getVA().keySet() );
963 final int firstLayerIndex = first.getParent().getLayerIndex( first.getId() );
964 for ( int i = 0; i < firstLayerIndex; ++i )
965 applyTransformToLayer( layers.get( i ), mlt, filter );
966 }
967 if ( propagateTransformAfter )
968 {
969 final Layer last = layerRange.get( layerRange.size() - 1 );
970 final MovingLeastSquaresTransform2 mlt = makeMLST2( meshes.get( meshes.size() - 1 ).getVA().keySet() );
971 final int lastLayerIndex = last.getParent().getLayerIndex( last.getId() );
972 for ( int i = lastLayerIndex + 1; i < layers.size(); ++i )
973 applyTransformToLayer( layers.get( i ), mlt, filter );
974 }
975 }
976 for ( int l = 0; l < layerRange.size(); ++l )
977 {
978 IJ.showStatus( "Applying transformation to patches ..." );
979 IJ.showProgress( 0, layerRange.size() );
980
981 final Layer layer = layerRange.get( l );
982
983 final MovingLeastSquaresTransform2 mlt = new MovingLeastSquaresTransform2();
984 mlt.setModel( AffineModel2D.class );
985 mlt.setAlpha( 2.0f );
986 mlt.setMatches( meshes.get( l ).getVA().keySet() );
987
988 applyTransformToLayer( layer, mlt, filter );
989
990 if ( Thread.interrupted() )
991 {
992 Utils.log( "Interrupted during applying transformations to patches. No all patches have been updated. Re-generate mipmaps manually." );
993 }
994
995 IJ.showProgress( l + 1, layerRange.size() );
996 }
997
998 /* update patch mipmaps */
999 final int firstLayerIndex;
1000 final int lastLayerIndex;
1001
1002 if ( propagateTransformBefore )
1003 firstLayerIndex = 0;
1004 else
1005 {
1006 firstLayerIndex = first.getParent().getLayerIndex( first.getId() );
1007 }
1008 if ( propagateTransformAfter )
1009 lastLayerIndex = layers.size() - 1;
1010 else
1011 {
1012 final Layer last = layerRange.get( layerRange.size() - 1 );
1013 lastLayerIndex = last.getParent().getLayerIndex( last.getId() );
1014 }
1015
1016 for ( int i = firstLayerIndex; i <= lastLayerIndex; ++i )
1017 {
1018 final Layer layer = layers.get( i );
1019 if ( !( emptyLayers.contains( layer ) || fixedLayers.contains( layer ) ) )
1020 {
1021 for ( final Patch patch : AlignmentUtils.filterPatches( layer, filter ) )
1022 patch.updateMipMaps();
1023 }
1024 }
1025
1026 Utils.log( "Done." );
1027 }
1028
1029 final static protected MovingLeastSquaresTransform2 makeMLST2( final Set< PointMatch > matches ) throws Exception
1030 {
1031 final MovingLeastSquaresTransform2 mlt = new MovingLeastSquaresTransform2();
1032 mlt.setModel( AffineModel2D.class );
1033 mlt.setAlpha( 2.0f );
1034 mlt.setMatches( matches );
1035 return mlt;
1036 }
1037
1038 final static protected void applyTransformToLayer( final Layer layer, final MovingLeastSquaresTransform2 mlt, final Filter< Patch > filter ) throws InterruptedException
1039 {
1040 /*
1041 * Setting a transformation to a patch can take some time because
1042 * the new bounding box needs to be estimated which requires the
1043 * TransformMesh to be generated and all vertices iterated.
1044 *
1045 * Therefore multithreading.
1046 */
1047 final List< Patch > patches = AlignmentUtils.filterPatches( layer, filter );
1048
1049 final ArrayList< Thread > applyThreads = new ArrayList< Thread >( p.maxNumThreads );
1050 final AtomicInteger ai = new AtomicInteger( 0 );
1051 for ( int t = 0; t < p.maxNumThreads; ++t )
1052 {
1053 final Thread thread = new Thread(
1054 new Runnable()
1055 {
1056 @Override
1057 final public void run()
1058 {
1059 try
1060 {
1061 for ( int i = ai.getAndIncrement(); i < patches.size() && !Thread.interrupted(); i = ai.getAndIncrement() )
1062 mpicbg.trakem2.align.Util.applyLayerTransformToPatch( patches.get( i ), mlt.copy() );
1063 }
1064 catch ( final Exception e )
1065 {
1066 e.printStackTrace();
1067 }
1068 }
1069 } );
1070 applyThreads.add( thread );
1071 thread.start();
1072 }
1073
1074 for ( final Thread thread : applyThreads )
1075 thread.join();
1076 }
1077
1078
1079 /**
1080 * Stateful. Changing the parameters of this instance. Do not use in parallel.
1081 *
1082 * @param layerRange
1083 * @param fixedLayers
1084 * @param propagateTransformBefore
1085 * @param propagateTransformAfter
1086 * @param fov
1087 * @param filter
1088 * @throws Exception
1089 */
1090 final public void exec(
1091 final Project project,
1092 final List< Layer > layerRange,
1093 final Set< Layer > fixedLayers,
1094 final boolean propagateTransformBefore,
1095 final boolean propagateTransformAfter,
1096 final Rectangle fov,
1097 final Filter< Patch > filter ) throws Exception
1098 {
1099 Rectangle box = null;
1100 final HashSet< Layer > emptyLayers = new HashSet< Layer >();
1101 for ( final Iterator< Layer > it = layerRange.iterator(); it.hasNext(); )
1102 {
1103 /* remove empty layers */
1104 final Layer la = it.next();
1105 if ( !la.contains( Patch.class, true ) )
1106 {
1107 emptyLayers.add( la );
1108 // it.remove();
1109 }
1110 else
1111 {
1112 /* accumulate boxes */
1113 if ( null == box ) // The first layer:
1114 box = la.getMinimalBoundingBox( Patch.class, true );
1115 else
1116 box = box.union( la.getMinimalBoundingBox( Patch.class, true ) );
1117 }
1118 }
1119
1120 if ( box == null )
1121 box = new Rectangle();
1122
1123 if ( fov != null )
1124 box = box.intersection( fov );
1125
1126 if ( box.width <= 0 || box.height <= 0 )
1127 {
1128 Utils.log( "Bounding box empty." );
1129 return;
1130 }
1131
1132 if ( layerRange.size() == emptyLayers.size() )
1133 {
1134 Utils.log( "All layers in range are empty!" );
1135 return;
1136 }
1137
1138 /* do not work if there is only one layer selected */
1139 if ( layerRange.size() - emptyLayers.size() < 2 )
1140 {
1141 Utils.log( "All except one layer in range are empty!" );
1142 return;
1143 }
1144
1145 if ( !p.setup( box ) ) return;
1146
1147 exec( p.clone(), project, layerRange, fixedLayers, emptyLayers, box, propagateTransformBefore, propagateTransformAfter, filter );
1148 }
1149
1150
1151 /**
1152 * Stateful. Changing the parameters of this instance. Do not use in parallel.
1153 *
1154 * @param layerSet
1155 * @param firstIn
1156 * @param lastIn
1157 * @param propagateTransformBefore
1158 * @param propagateTransformAfter
1159 * @param fov
1160 * @param filter
1161 */
1162 final public void exec(
1163 final LayerSet layerSet,
1164 final int firstIn,
1165 final int lastIn,
1166 final int ref,
1167 final boolean propagateTransformBefore,
1168 final boolean propagateTransformAfter,
1169 final Rectangle fov,
1170 final Filter< Patch > filter ) throws Exception
1171 {
1172 final int first = Math.min( firstIn, lastIn );
1173 final int last = Math.max( firstIn, lastIn );
1174
1175 /* always first index first despite the method would return inverse order if last > first */
1176 final List< Layer > layerRange = layerSet.getLayers( first, last );
1177 final HashSet< Layer > fixedLayers = new HashSet< Layer >();
1178
1179 if ( ref - firstIn >= 0 )
1180 fixedLayers.add( layerRange.get( ref - firstIn ) );
1181
1182 Utils.log( layerRange.size() + "" );
1183
1184 exec( layerSet.getProject(), layerRange, fixedLayers, propagateTransformBefore, propagateTransformAfter, fov, filter );
1185 }
1186
1187
1188 /**
1189 * Stateful. Changing the parameters of this instance. Do not use in parallel.
1190 *
1191 * @param layerSet
1192 * @param firstIn
1193 * @param lastIn
1194 * @param ref
1195 * @param propagateTransform
1196 * @param fov
1197 * @param filter
1198 */
1199 final public void exec(
1200 final LayerSet layerSet,
1201 final int firstIn,
1202 final int lastIn,
1203 final int ref,
1204 final boolean propagateTransform,
1205 final Rectangle fov,
1206 final Filter< Patch > filter ) throws Exception
1207 {
1208 if ( firstIn < lastIn )
1209 exec( layerSet, firstIn, lastIn, ref, false, propagateTransform, fov, filter );
1210 else
1211 exec( layerSet, firstIn, lastIn, ref, propagateTransform, false, fov, filter );
1212 }
1213
1214 /**
1215 * Stateful. Changing the parameters of this instance. Do not use in parallel.
1216 *
1217 * @param layerSet
1218 * @param firstIn
1219 * @param lastIn
1220 * @param ref1
1221 * @param ref2
1222 * @param propagateTransformAfter
1223 * @param fov
1224 * @param filter
1225 */
1226 final public void exec(
1227 final LayerSet layerSet,
1228 final int firstIn,
1229 final int lastIn,
1230 final int ref1,
1231 final int ref2,
1232 final boolean propagateTransformBefore,
1233 final boolean propagateTransformAfter,
1234 final Rectangle fov,
1235 final Filter< Patch > filter ) throws Exception
1236 {
1237 final int first = Math.min( firstIn, lastIn );
1238 final int last = Math.max( firstIn, lastIn );
1239
1240 /* always first index first despite the method would return inverse order if last > first */
1241 final List< Layer > layerRange = layerSet.getLayers( first, last );
1242 final HashSet< Layer > fixedLayers = new HashSet< Layer >();
1243
1244 if ( ref1 - first >= 0 )
1245 fixedLayers.add( layerRange.get( ref1 - first ) );
1246 if ( ref2 - first >= 0 )
1247 fixedLayers.add( layerRange.get( ref2 - first ) );
1248
1249 Utils.log( layerRange.size() + "" );
1250
1251 exec( layerSet.getProject(), layerRange, fixedLayers, propagateTransformBefore, propagateTransformAfter, fov, filter );
1252 }
1253 }
register/edit/analyze/3D model TB of images