Seg3D  2.4
Seg3D is a free volume segmentation and processing tool developed by the NIH Center for Integrative Biomedical Computing at the University of Utah Scientific Computing and Imaging (SCI) Institute.
mosaic_layout_common.hxx
1 /*
2  For more information, please see: http://software.sci.utah.edu
3 
4  The MIT License
5 
6  Copyright (c) 2016 Scientific Computing and Imaging Institute,
7  University of Utah.
8 
9 
10  Permission is hereby granted, free of charge, to any person obtaining a
11  copy of this software and associated documentation files (the "Software"),
12  to deal in the Software without restriction, including without limitation
13  the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  and/or sell copies of the Software, and to permit persons to whom the
15  Software is furnished to do so, subject to the following conditions:
16 
17  The above copyright notice and this permission notice shall be included
18  in all copies or substantial portions of the Software.
19 
20  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  DEALINGS IN THE SOFTWARE.
27  */
28 
29 // File : mosaic_layout_common.hxx
30 // Author : Pavel A. Koshevoy
31 // Created : Thu Mar 22 13:09:33 MDT 2007
32 // Copyright : (C) 2004-2008 University of Utah
33 // Description : Helper functions for automatic mosaic layout.
34 
35 #ifndef MOSAIC_LAYOUT_COMMON_HXX_
36 #define MOSAIC_LAYOUT_COMMON_HXX_
37 
38 // the includes:
39 #include <Core/ITKCommon/common.hxx>
40 #include <Core/ITKCommon/FFT/fft_common.hxx>
41 #include <Core/ITKCommon/grid_common.hxx>
42 #include <Core/ITKCommon/Transform/itkLegendrePolynomialTransform.h>
43 #include <Core/ITKCommon/Optimizers/itkRegularStepGradientDescentOptimizer2.h>
44 #include <Core/ITKCommon/Optimizers/itkImageMosaicVarianceMetric.h>
45 #include <Core/ITKCommon/the_utils.hxx>
46 
47 // ITK includes:
48 #include <itkCommand.h>
49 #include <itkImageMaskSpatialObject.h>
50 #include <itkImageRegistrationMethod.h>
51 #include <itkMeanSquaresImageToImageMetric.h>
52 
53 // boost:
54 #include <boost/filesystem.hpp>
55 
56 namespace bfs=boost::filesystem;
57 
58 #include <sstream>
59 #include <string>
60 
61 // affine_transform_t
62 //
64 
65 enum StrategyType
66 {
67  DEFAULT,
68  TOP_LEFT_BOOK,
69  TOP_LEFT_SNAKE
70 };
71 
72 
73 //----------------------------------------------------------------
74 // reset
75 //
76 extern void
77 reset(const unsigned int num_images,
78  const unsigned int max_cascade_len,
79  array3d(translate_transform_t::Pointer) & path,
80  array3d(double) & cost);
81 
82 //----------------------------------------------------------------
83 // setup_transform
84 //
85 // TODO: clean this code up
86 template <typename TImage>
87 affine_transform_t::Pointer
88 setup_transform(const TImage * image,
89  const base_transform_t * t0,
90  const base_transform_t * t1 = NULL,
91  const unsigned int samples = 16)
92 {
93  // shortcuts:
94 #define inv_transform \
95  inverse_transform<TImage, std::vector<base_transform_t::Pointer> >
96 
97 #define fwd_transform \
98  forward_transform<TImage, std::vector<const base_transform_t *> >
99 
100  typedef typename TImage::IndexType index_t;
101  index_t i00;
102  i00[0] = 0;
103  i00[1] = 0;
104 
105  typename TImage::PointType origin;
106  image->TransformIndexToPhysicalPoint(i00, origin);
107 
108  index_t i11;
109  i11[0] = 1;
110  i11[1] = 1;
111 
112  typename TImage::PointType spacing;
113  image->TransformIndexToPhysicalPoint(i11, spacing);
114  spacing[0] -= origin[0];
115  spacing[1] -= origin[1];
116 
117  typename TImage::SizeType sz = image->GetLargestPossibleRegion().GetSize();
118 
119  typename TImage::PointType tile_min = origin;
120  typename TImage::PointType tile_max;
121  tile_max[0] = tile_min[0] + spacing[0] * static_cast<double>(sz[0]);
122  tile_max[1] = tile_min[1] + spacing[1] * static_cast<double>(sz[1]);
123 
124  double w = sz[0] * spacing[0];
125  double h = sz[1] * spacing[1];
126  double Umax = w / 2.0;
127  double Vmax = h / 2.0;
128 
129  affine_transform_t::Pointer affine = affine_transform_t::New();
130  affine->setup(tile_min[0],
131  tile_max[0],
132  tile_min[1],
133  tile_max[1],
134  Umax,
135  Vmax);
136 
137  std::vector<const base_transform_t *> fwd_cascade(2);
138  std::vector<base_transform_t::Pointer> inv_cascade(2);
139  fwd_cascade[0] = t0;
140  fwd_cascade[1] = t1;
141 
142  inv_cascade[0] = t0->GetInverseTransform();
143  inv_cascade[1] = (t1 == NULL) ? NULL : t1->GetInverseTransform();
144 
145  const unsigned int cascade_len = (t1 == NULL) ? 1 : 2;
146 
147  // calculate the shift:
148  pnt2d_t center = pnt2d(tile_min[0] + (tile_max[0] - tile_min[0]) / 2.0,
149  tile_min[1] + (tile_max[1] - tile_min[1]) / 2.0);
150  pnt2d_t transformed_center = inv_transform(inv_cascade, cascade_len, center);
151  vec2d_t shift = center - transformed_center;
152  affine->setup_translation(shift[0], shift[1]);
153 
154  // initialize the local bounding box:
155  typename TImage::PointType mosaic_min;
156  typename TImage::PointType mosaic_max;
157  mosaic_min[0] = std::numeric_limits<double>::max();
158  mosaic_min[1] = mosaic_min[0];
159  mosaic_max[0] = -mosaic_min[0];
160  mosaic_max[1] = -mosaic_min[0];
161 
162  // find out where this image is mapped to in the mosaic
163  // by the given transform:
164  {
165  // a handy shortcut:
166  typename TImage::PointType tile_point;
167 #define UPDATE_BBOX( u, v ) \
168  tile_point[0] = tile_min[0] + u * (tile_max[0] - tile_min[0]); \
169  tile_point[1] = tile_min[1] + v * (tile_max[1] - tile_min[1]); \
170  update_bbox(mosaic_min, mosaic_max, \
171  inverse_transform<TImage, \
172  std::vector<base_transform_t::Pointer> > \
173  (inv_cascade, cascade_len, tile_point))
174 
175  // image corners:
176  UPDATE_BBOX(0.0, 0.0);
177  UPDATE_BBOX(1.0, 0.0);
178  UPDATE_BBOX(0.0, 1.0);
179  UPDATE_BBOX(1.0, 1.0);
180 
181  // points along the image edges:
182  for (unsigned int x = 0; x < samples; x++)
183  {
184  const double t = static_cast<double>(x + 1) / static_cast<double>(samples + 1);
185 
186  UPDATE_BBOX(t, 0.0);
187  UPDATE_BBOX(t, 1.0);
188  UPDATE_BBOX(0.0, t);
189  UPDATE_BBOX(1.0, t);
190  }
191 #undef UPDATE_BBOX
192  }
193 
194  const double W = mosaic_max[0] - mosaic_min[0];
195  const double H = mosaic_max[1] - mosaic_min[1];
196 
197  std::vector<pnt2d_t> uv(samples * samples);
198  std::vector<pnt2d_t> xy(samples * samples);
199 
200  // explore the transform with jittered sampling:
201  for (unsigned int i = 0; i < samples; i++)
202  {
203  for (unsigned int j = 0; j < samples; j++)
204  {
205  const unsigned int k = i * samples + j;
206 
207  uv[k][0] =
208  mosaic_min[0] + W * (static_cast<double>(i) + drand()) / static_cast<double>(samples);
209  uv[k][1] =
210  mosaic_min[1] + H * (static_cast<double>(j) + drand()) / static_cast<double>(samples);
211 
212  xy[k] = fwd_transform(fwd_cascade, cascade_len, uv[k]);
213  }
214  }
215 
216 #undef inv_transform
217 #undef fwd_transform
218 
219  affine->solve_for_parameters(0, affine_transform_t::Degree + 1, uv, xy);
220  return affine;
221 }
222 
223 //----------------------------------------------------------------
224 // assemble_cascades
225 //
226 extern void
227 assemble_cascades(const unsigned int num_images,
228  const unsigned int max_cascade_len,
229  array3d(translate_transform_t::Pointer) & path,
230  array3d(double) & cost,
231  const bool & cumulative_cost = false);
232 
233 //----------------------------------------------------------------
234 // establish_mappings
235 //
236 extern void
237 establish_mappings(const unsigned int num_images,
238  const unsigned int max_cascade_len,
239  const array3d(translate_transform_t::Pointer) & path,
240  const array3d(double) & cost,
241  array2d(translate_transform_t::Pointer) & mapping,
242  array2d(double) & mapping_cost);
243 
244 
245 
246 //----------------------------------------------------------------
247 // refine_one_pair
248 //
249 template <typename TImage, typename TMask>
250 double
251 refine_one_pair(const TImage * i0,
252  const TImage * i1,
253  const TMask * m0,
254  const TMask * m1,
255  translate_transform_t::Pointer & t01,
256 
257  const unsigned int iterations,
258  const double & min_step,
259  const double & max_step,
260 
261  const std::string & fn_prefix)
262 {
263  // setup the registration object:
264  typedef itk::ImageRegistrationMethod<TImage, TImage> registration_t;
265  typename registration_t::Pointer registration = registration_t::New();
266 
267  // setup the image interpolator:
268  typedef itk::LinearInterpolateImageFunction<TImage, double> interpolator_t;
269  registration->SetInterpolator(interpolator_t::New());
270 
271  if (! fn_prefix.empty())
272  {
273  save_rgb<TImage>(fn_prefix + "a.tif", i0, i1, t01, m0, m1);
274  }
275 
276  // setup the optimizer:
277  optimizer_t::Pointer optimizer = optimizer_t::New();
278  registration->SetOptimizer(optimizer);
279  optimizer_observer_t<optimizer_t>::Pointer observer =
281  optimizer->AddObserver(itk::IterationEvent(), observer);
282  optimizer->SetMinimize(true);
283  optimizer->SetNumberOfIterations(iterations);
284  optimizer->SetMinimumStepLength(min_step);
285  optimizer->SetMaximumStepLength(max_step);
286  optimizer->SetGradientMagnitudeTolerance(1e-6);
287  optimizer->SetRelaxationFactor(5e-1);
288  optimizer->SetPickUpPaceSteps(5);
289  optimizer->SetBackTracking(true);
290 
291  // FIXME: this is probably unnecessary:
292  typedef optimizer_t::ScalesType optimizer_scales_t;
293  optimizer_scales_t translation_scales(t01->GetNumberOfParameters());
294  translation_scales.Fill(1.0);
295 
296  try { optimizer->SetScales(translation_scales); }
297  catch (itk::ExceptionObject & exception) {}
298 
299  // setup the image-to-image metric:
300  typedef itk::MeanSquaresImageToImageMetric<TImage, TImage> metric_t;
301 
302  typename metric_t::Pointer metric = metric_t::New();
303  registration->SetMetric(metric);
304 
305  registration->SetTransform(t01);
306  registration->SetInitialTransformParameters(t01->GetParameters());
307 
308  // setup the masks:
309  typedef itk::ImageMaskSpatialObject<2> mask_so_t;
310  typename TMask::ConstPointer fi_mask = m0;
311  if (m0 != NULL)
312  {
313  mask_so_t::Pointer fi_mask_so = mask_so_t::New();
314  fi_mask_so->SetImage(fi_mask);
315  metric->SetFixedImageMask(fi_mask_so);
316  }
317 
318  typename TMask::ConstPointer mi_mask = m1;
319  if (m1 != NULL)
320  {
321  mask_so_t::Pointer mi_mask_so = mask_so_t::New();
322  mi_mask_so->SetImage(mi_mask);
323  metric->SetMovingImageMask(mi_mask_so);
324  }
325 
326  // setup the fixed and moving image:
327  typename TImage::ConstPointer fi = i0;
328  typename TImage::ConstPointer mi = i1;
329 
330  registration->SetFixedImageRegion(fi->GetLargestPossibleRegion());
331  registration->SetFixedImage(fi);
332  registration->SetMovingImage(mi);
333 
334  // evaluate the metric before the registration:
335  double metric_before =
336  eval_metric<metric_t, interpolator_t>(t01, fi, mi, fi_mask, mi_mask);
337  assert(metric_before != std::numeric_limits<double>::max());
338 
339  translate_transform_t::ParametersType params_before = t01->GetParameters();
340 
341  std::ostringstream oss;
342  // perform the registration:
343  try
344  {
345  registration->StartRegistration();
346  }
347  catch (itk::ExceptionObject & exception)
348  {
349  oss << std::endl << "image registration threw an exception: "
350  << std::endl << exception.what() << std::endl;
351  }
352  t01->SetParameters(optimizer->GetBestParams());
353 
354  // evaluate the metric after the registration:
355  double metric_after =
356  eval_metric<metric_t, interpolator_t>(t01, fi, mi, fi_mask, mi_mask);
357 
358  translate_transform_t::ParametersType params_after =
359  t01->GetParameters();
360 
361  oss << std::endl << "BEFORE: " << metric_before << std::endl
362  << "AFTER: " << metric_after << std::endl;
363 
364  if (metric_before <= metric_after || metric_after != metric_after)
365  {
366  oss << "NOTE: minimization failed, ignoring registration results..."
367  << std::endl;
368  t01->SetParameters(params_before);
369  metric_after = metric_before;
370  }
371 
372  if (! fn_prefix.empty())
373  {
374  save_rgb<TImage>(fn_prefix + "b.tif", i0, i1, t01, m0, m1);
375  }
376 
377  CORE_LOG_MESSAGE(oss.str());
378  return metric_after;
379 }
380 
381 
382 //----------------------------------------------------------------
383 // refine_one_pair
384 //
385 template <class TImage, class TMask>
386 double
387 refine_one_pair(const array2d(typename TImage::Pointer) & tile_pyramid,
388  const array2d(typename TMask::ConstPointer) & mask_pyramid,
389  const unsigned int & ia,
390  const unsigned int & ib,
391  translate_transform_t::Pointer & t01,
392 
393  const unsigned int iterations,
394  const double & min_step,
395  const double & max_step,
396 
397  const bfs::path & prefix)
398 {
399  const unsigned int num_levels = tile_pyramid.size();
400 
401  double metric = std::numeric_limits<double>::max();
402  for (unsigned int i = 0; i < num_levels; i++)
403  {
404  std::ostringstream fn_prefix;
405  if (! prefix.empty())
406  {
407  fn_prefix << prefix << "level-" << the_text_t::number(i) << "-";
408  }
409 
410  metric = refine_one_pair<TImage, TMask>(tile_pyramid[i][ia],
411  tile_pyramid[i][ib],
412  mask_pyramid[i][ia],
413  mask_pyramid[i][ib],
414  t01,
415  iterations,
416  min_step,
417  max_step /* * integer_power<double>(2.0, num_levels - i - 1)*/,
418  fn_prefix.str());
419  }
420 
421  return metric;
422 }
423 
424 
425 //----------------------------------------------------------------
426 // refine_pairs
427 //
428 template <class TImage, class TMask>
429 void
430 refine_pairs(const std::vector<typename TImage::Pointer> & image,
431  const std::vector<typename TMask::ConstPointer> & mask,
432  const double & overlap_min,
433  const double & overlap_max,
434  const unsigned int & iterations,
435  const double & min_step,
436  const double & max_step,
437 
438  const array2d(translate_transform_t::Pointer) & mapping,
439  array2d(translate_transform_t::Pointer) & path,
440  array2d(double) & cost,
441 
442  const bfs::path & prefix)
443 {
444  static unsigned int pass = 0;
445 
446  std::ostringstream oss;
447  const unsigned int num_images = image.size();
448  for (unsigned int i = 0; i < num_images - 1; i++)
449  {
450  for (unsigned int j = i + 1; j < num_images; j++)
451  {
452  // check whether the two images overlap:
453  double overlap = overlap_ratio<TImage>(image[i], image[j], mapping[i][j]);
454  if (overlap < overlap_min || overlap > overlap_max) continue;
455 
456  translate_transform_t::Pointer t_ij = translate_transform_t::New();
457  t_ij->SetParameters(mapping[i][j]->GetParameters());
458 
459  oss << "refining the mapping: " << std::setw(2) << j << " -> "
460  << std::setw(2) << i << std::endl
461  << "overlap: " << static_cast<int>(overlap * 100.0) << " percent" << std::endl;
462 
463  std::ostringstream fn_prefix;
464  if (! prefix.empty())
465  {
466  fn_prefix << prefix << "pair-" << the_text_t::number(i, 2, '0') << "-" <<
467  the_text_t::number(j, 2, '0') << "-pass-" << the_text_t::number(pass, 2, '0') << "-";
468  }
469 
470  cost[i][j] = refine_one_pair<TImage, TMask>(image[i],
471  image[j],
472  mask[i],
473  mask[j],
474  t_ij,
475  iterations,
476  min_step,
477  max_step,
478  fn_prefix.str());
479 
480  overlap = overlap_ratio<TImage>(image[i], image[j], t_ij);
481 
482  if (overlap < overlap_min || overlap > overlap_max)
483  {
484  path[i][j] = NULL;
485  cost[i][j] = std::numeric_limits<double>::max();
486  }
487  else
488  {
489  // TODO: encourage larger overlap:
490  // cost[i][j] = 1.0 - overlap; // encourage edge overlap
491  // cost[i][j] = overlap; // encourage corner overlap, works well !!!
492  path[i][j] = t_ij;
493  }
494 
495  path[j][i] = inverse(t_ij);
496  cost[j][i] = cost[i][j];
497  }
498  oss << std::endl;
499  }
500 
501  // int prev_precision = oss.precision(2);
502  oss << "\t";
503  for (unsigned int i = 0; i < num_images; i++) oss << std::setw(9) << i;
504  oss << std::endl;
505  oss << "\t+-";
506  for (unsigned int i = 0; i < num_images; i++) oss << "---------";
507  oss << std::endl;
508  for (unsigned int i = 0; i < num_images; i++)
509  {
510  oss << i << "\t| ";
511  for (unsigned int j = 0; j < num_images; j++)
512  {
513  if (path[i][j].GetPointer() != NULL)
514  {
515  oss << ' ' << cost[i][j];
516  }
517  else
518  {
519  oss << " ";
520  }
521  }
522  oss << std::endl;
523  }
524  oss << std::endl;
525  // oss.precision(prev_precision);
526  CORE_LOG_MESSAGE(oss.str());
527 }
528 
529 //----------------------------------------------------------------
530 // refine_pairs
531 //
532 template <class TImage, class TMask>
533 void
534 refine_pairs(const array2d(typename TImage::Pointer) & tile_pyramid,
535  const array2d(typename TMask::ConstPointer) & mask_pyramid,
536  const double & overlap_min,
537  const double & overlap_max,
538  const unsigned int & iterations,
539  const double & min_step,
540  const double & max_step,
541 
542  const array2d(translate_transform_t::Pointer) & mapping,
543  array2d(translate_transform_t::Pointer) & path,
544  array2d(double) & cost,
545 
546  const bfs::path & prefix)
547 {
548  static unsigned int pass = 0;
549  pass++;
550 
551  unsigned int pyramid_levels = tile_pyramid.size();
552  unsigned int high_res_level = pyramid_levels - 1;
553 
554  const unsigned int num_images = tile_pyramid[high_res_level].size();
555  std::ostringstream oss;
556  for (unsigned int i = 0; i < num_images - 1; i++)
557  {
558  for (unsigned int j = i + 1; j < num_images; j++)
559  {
560  // check whether the two images overlap:
561  double overlap =
562  overlap_ratio<TImage>(tile_pyramid[high_res_level][i],
563  tile_pyramid[high_res_level][j],
564  mapping[i][j]);
565  if (overlap < overlap_min || overlap > overlap_max) continue;
566 
567  translate_transform_t::Pointer t_ij = translate_transform_t::New();
568  t_ij->SetParameters(mapping[i][j]->GetParameters());
569 
570  oss << "refining the mapping: " << std::setw(2) << j << " -> "
571  << std::setw(2) << i << std::endl
572  << "overlap: " << static_cast<int>(overlap * 100.0) << " percent" << std::endl;
573 
574  std::ostringstream fn_prefix;
575  if (! prefix.empty())
576  {
577  fn_prefix << prefix << "pair-" << the_text_t::number(i, 2, '0') << "-" <<
578  the_text_t::number(j, 2, '0') << "-pass-" << the_text_t::number(pass, 2, '0') << "-";
579  }
580 
581  cost[i][j] = refine_one_pair<TImage, TMask>(tile_pyramid,
582  mask_pyramid,
583  i,
584  j,
585  t_ij,
586  iterations,
587  min_step,
588  max_step,
589  fn_prefix.str());
590 
591  overlap = overlap_ratio<TImage>(tile_pyramid[high_res_level][i],
592  tile_pyramid[high_res_level][j],
593  t_ij);
594 
595  if (overlap < overlap_min || overlap > overlap_max)
596  {
597  path[i][j] = NULL;
598  cost[i][j] = std::numeric_limits<double>::max();
599  }
600  else
601  {
602  // TODO: encourage larger overlap:
603  // cost[i][j] = 1.0 - overlap; // encourage edge overlap
604  // cost[i][j] = overlap; // encourage corner overlap, works well !!!
605  path[i][j] = t_ij;
606  }
607 
608  path[j][i] = inverse(t_ij);
609  cost[j][i] = cost[i][j];
610  }
611  oss << std::endl;
612  }
613 
614  // int prev_precision = oss.precision(2);
615  oss << "\t";
616  for (unsigned int i = 0; i < num_images; i++) oss << std::setw(9) << i;
617  oss << std::endl;
618  oss << "\t+-";
619  for (unsigned int i = 0; i < num_images; i++) oss << "---------";
620  oss << std::endl;
621  for (unsigned int i = 0; i < num_images; i++)
622  {
623  oss << i << "\t| ";
624  for (unsigned int j = 0; j < num_images; j++)
625  {
626  if (path[i][j].GetPointer() != NULL)
627  {
628  oss << ' ' << cost[i][j];
629  }
630  else
631  {
632  oss << " ";
633  }
634  }
635  oss << std::endl;
636  }
637  oss << std::endl;
638  // oss.precision(prev_precision);
639  CORE_LOG_MESSAGE(oss.str());
640 }
641 
642 
643 //----------------------------------------------------------------
644 // brute_force_one_pair
645 //
646 template <typename TImage, typename TMask>
647 double
648 brute_force_one_pair(const TImage * i0,
649  const TImage * i1,
650  const TMask * m0,
651  const TMask * m1,
652  translate_transform_t::Pointer & t01,
653 
654  const int & dx,
655  const int & dy,
656 
657  const std::string & fn_prefix)
658 {
659  translate_transform_t::ParametersType init_params = t01->GetParameters();
660  translate_transform_t::ParametersType best_params = init_params;
661  translate_transform_t::ParametersType curr_params = init_params;
662 
663  translate_transform_t::Pointer t = translate_transform_t::New();
664  double best_metric = std::numeric_limits<double>::max();
665  typename TImage::SpacingType sp = i0->GetSpacing();
666 
667  typedef typename itk::NearestNeighborInterpolateImageFunction
668  <TImage, double> nn_interpolator_t;
669  typename nn_interpolator_t::Pointer nn = nn_interpolator_t::New();
670  nn->SetInputImage(i1);
671 
672  if (! fn_prefix.empty())
673  {
674  save_rgb<TImage>(fn_prefix + "a.tif", i0, i1, t01, m0, m1);
675  }
676 
677  for (int x = -dx; x <= dx; x++)
678  {
679  for (int y = -dy; y <= dy; y++)
680  {
681  // the window:
682  curr_params[0] = init_params[0] + sp[0] * static_cast<double>(x);
683  curr_params[1] = init_params[1] + sp[1] * static_cast<double>(y);
684 
685  t->SetParameters(curr_params);
686 
687  double overlap_area = 0;
688  double metric =
689  my_metric<TImage, nn_interpolator_t>
690  (overlap_area, i0, i1, t, m0, m1, nn);
691 
692  if (metric < best_metric)
693  {
694  best_metric = metric;
695  best_params = curr_params;
696  }
697  }
698  }
699 
700  t01->SetParameters(best_params);
701 
702  if (! fn_prefix.empty())
703  {
704  save_rgb<TImage>(fn_prefix + "b.tif", i0, i1, t01, m0, m1);
705  }
706 
707  return best_metric;
708 }
709 
710 
711 //----------------------------------------------------------------
712 // brute_force_one_pair
713 //
714 template <typename TImage, typename TMask>
715 double
716 brute_force_one_pair(const array2d(typename TImage::Pointer) & tile_pyramid,
717  const array2d(typename TMask::ConstPointer) & mask_pyramid,
718  const unsigned int & ia,
719  const unsigned int & ib,
720  translate_transform_t::Pointer & t01,
721 
722  const int & dx,
723  const int & dy,
724  const unsigned int & coarse_to_fine_levels,
725 
726  const bfs::path & prefix)
727 {
728  unsigned int num_levels = std::min(coarse_to_fine_levels,
729  static_cast<unsigned int>(tile_pyramid.size()));
730  unsigned int high_res_level = num_levels - 1;
731 
732  for (unsigned int i = 0; i < num_levels; i++)
733  {
734  std::ostringstream fn_prefix;
735  if (! prefix.empty())
736  {
737  fn_prefix << prefix << "level-" << the_text_t::number(i) << "-";
738  }
739 
740  brute_force_one_pair<TImage, TMask>(tile_pyramid[i][ia].GetPointer(),
741  tile_pyramid[i][ib].GetPointer(),
742  mask_pyramid[i][ia].GetPointer(),
743  mask_pyramid[i][ib].GetPointer(),
744  t01,
745  dx,
746  dy,
747  fn_prefix.str());
748  }
749 
750  // setup the image-to-image metric:
751  typedef itk::MeanSquaresImageToImageMetric<TImage, TImage> metric_t;
752 
753  typedef itk::LinearInterpolateImageFunction<TImage, double> interpolator_t;
754  double metric =
755  eval_metric<metric_t, interpolator_t>(t01,
756  tile_pyramid[high_res_level][ia],
757  tile_pyramid[high_res_level][ib],
758  mask_pyramid[high_res_level][ia],
759  mask_pyramid[high_res_level][ib]);
760  return metric;
761 }
762 
763 
764 //----------------------------------------------------------------
765 // brute_force_pairs
766 //
767 template <typename TImage, typename TMask>
768 void
769 brute_force_pairs(const std::vector<typename TImage::Pointer> & image,
770  const std::vector<typename TMask::ConstPointer> & mask,
771  const double & overlap_min,
772  const double & overlap_max,
773  const int & dx,
774  const int & dy,
775 
776  const array2d(translate_transform_t::Pointer) & mapping,
777  array2d(translate_transform_t::Pointer) & path,
778  array2d(double) & cost,
779 
780  const bfs::path & prefix)
781 {
782  static unsigned int pass = 0;
783  pass++;
784 
785  const unsigned int num_images = image.size();
786  std::ostringstream oss;
787  for (unsigned int i = 0; i < num_images - 1; i++)
788  {
789  for (unsigned int j = i + 1; j < num_images; j++)
790  {
791  // check whether the two images overlap:
792  double overlap =
793  overlap_ratio<TImage>(image[i], image[j], mapping[i][j]);
794  if (overlap < overlap_min || overlap > overlap_max) continue;
795 
796  translate_transform_t::Pointer t_ij = translate_transform_t::New();
797  t_ij->SetParameters(mapping[i][j]->GetParameters());
798 
799  oss << "refining the mapping: " << std::setw(2) << j << " -> "
800  << std::setw(2) << i << std::endl
801  << "overlap: " << static_cast<int>(overlap * 100.0) << " percent" << std::endl;
802 
803  std::ostringstream fn_prefix;
804  if (! prefix.empty())
805  {
806  fn_prefix << prefix << "pair-" << the_text_t::number(i, 2, '0') << "-" <<
807  the_text_t::number(j, 2, '0') << "-pass-" << the_text_t::number(pass, 2, '0') << "-";
808  }
809 
810  cost[i][j] = brute_force_one_pair<TImage, TMask>(image[i],
811  image[j],
812  mask[i],
813  mask[j],
814  t_ij,
815  dx,
816  dy,
817  fn_prefix.str());
818 
819  overlap = overlap_ratio<TImage>(image[i], image[j], t_ij);
820 
821  if (overlap < overlap_min || overlap > overlap_max)
822  {
823  path[i][j] = NULL;
824  cost[i][j] = std::numeric_limits<double>::max();
825  }
826  else
827  {
828  // FIXME: encourage larger overlap:
829  // cost[i][j] = 1.0 - overlap; // encourage edge overlap
830  // cost[i][j] = overlap; // encourage corner overlap, works well !!!
831  path[i][j] = t_ij;
832  }
833 
834  path[j][i] = inverse(t_ij);
835  cost[j][i] = cost[i][j];
836  }
837  oss << std::endl;
838  }
839 
840  // int prev_precision = oss.precision(2);
841  oss << "\t";
842  for (unsigned int i = 0; i < num_images; i++) oss << std::setw(9) << i;
843  oss << std::endl;
844  oss << "\t+-";
845  for (unsigned int i = 0; i < num_images; i++) oss << "---------";
846  oss << std::endl;
847  for (unsigned int i = 0; i < num_images; i++)
848  {
849  oss << i << "\t| ";
850  for (unsigned int j = 0; j < num_images; j++)
851  {
852  if (path[i][j].GetPointer() != NULL)
853  {
854  oss << ' ' << cost[i][j];
855  }
856  else
857  {
858  oss << " ";
859  }
860  }
861  oss << std::endl;
862  }
863  oss << std::endl;
864  // oss.precision(prev_precision);
865  CORE_LOG_MESSAGE(oss.str());
866 }
867 
868 //----------------------------------------------------------------
869 // brute_force_pairs
870 //
871 template <typename TImage, typename TMask>
872 void
873 brute_force_pairs(const array2d(typename TImage::Pointer) & tile_pyramid,
874  const array2d(typename TMask::ConstPointer) & mask_pyramid,
875  const double & overlap_min,
876  const double & overlap_max,
877  const int & dx,
878  const int & dy,
879 
880  // number of levels to try, starting from the lowest:
881  const unsigned int & coarse_to_fine_levels,
882 
883  const array2d(translate_transform_t::Pointer) & mapping,
884  array2d(translate_transform_t::Pointer) & path,
885  array2d(double) & cost,
886 
887  const bfs::path & prefix)
888 {
889  static unsigned int pass = 0;
890  pass++;
891 
892  unsigned int pyramid_levels = tile_pyramid.size();
893  unsigned int high_res_level = pyramid_levels - 1;
894 
895  const unsigned int num_images = tile_pyramid[high_res_level].size();
896  std::ostringstream oss;
897  for (unsigned int i = 0; i < num_images - 1; i++)
898  {
899  for (unsigned int j = i + 1; j < num_images; j++)
900  {
901  // check whether the two images overlap:
902  double overlap =
903  overlap_ratio<TImage>(tile_pyramid[high_res_level][i],
904  tile_pyramid[high_res_level][j],
905  mapping[i][j]);
906  if (overlap < overlap_min || overlap > overlap_max) continue;
907 
908  translate_transform_t::Pointer t_ij = translate_transform_t::New();
909  t_ij->SetParameters(mapping[i][j]->GetParameters());
910 
911  oss << "refining the mapping: " << std::setw(2) << j << " -> "
912  << std::setw(2) << i << std::endl
913  << "overlap: " << static_cast<int>(overlap * 100.0) << " percent" << std::endl;
914 
915  std::ostringstream fn_prefix;
916  if (! prefix.empty())
917  {
918  fn_prefix << prefix << "pair-" << the_text_t::number(i, 2, '0') << "-" <<
919  the_text_t::number(j, 2, '0') << "-pass-" << the_text_t::number(pass, 2, '0') << "-";
920  }
921 
922  cost[i][j] = brute_force_one_pair<TImage, TMask>(tile_pyramid,
923  mask_pyramid,
924  i,
925  j,
926  t_ij,
927  dx,
928  dy,
929  coarse_to_fine_levels,
930  fn_prefix.str());
931 
932  overlap = overlap_ratio<TImage>(tile_pyramid[high_res_level][i],
933  tile_pyramid[high_res_level][j],
934  t_ij);
935 
936  if (overlap < overlap_min || overlap > overlap_max)
937  {
938  path[i][j] = NULL;
939  cost[i][j] = std::numeric_limits<double>::max();
940  }
941  else
942  {
943  // FIXME: encourage larger overlap:
944  // cost[i][j] = 1.0 - overlap; // encourage edge overlap
945  // cost[i][j] = overlap; // encourage corner overlap, works well !!!
946  path[i][j] = t_ij;
947  }
948 
949  path[j][i] = inverse(t_ij);
950  cost[j][i] = cost[i][j];
951  }
952  oss << std::endl;
953  }
954 
955  // int prev_precision = oss.precision(2);
956  oss << "\t";
957  for (unsigned int i = 0; i < num_images; i++) oss << std::setw(9) << i;
958  oss << std::endl;
959  oss << "\t+-";
960  for (unsigned int i = 0; i < num_images; i++) oss << "---------";
961  oss << std::endl;
962  for (unsigned int i = 0; i < num_images; i++)
963  {
964  oss << i << "\t| ";
965  for (unsigned int j = 0; j < num_images; j++)
966  {
967  if (path[i][j].GetPointer() != NULL)
968  {
969  oss << ' ' << cost[i][j];
970  }
971  else
972  {
973  oss << " ";
974  }
975  }
976  oss << std::endl;
977  }
978  oss << std::endl;
979  // oss.precision(prev_precision);
980  CORE_LOG_MESSAGE(oss.str());
981 }
982 
983 
984 //----------------------------------------------------------------
985 // dump_neighbors
986 //
987 template <typename TImage, typename TMask>
988 void
989 dump_neighbors(const std::vector<typename TImage::Pointer> & image,
990  const std::vector<typename TMask::ConstPointer> & mask,
991  const array2d(translate_transform_t::Pointer) & path,
992  const bfs::path & prefix)
993 {
994  static unsigned int pass = 0;
995  pass++;
996 
997  const unsigned int num_images = image.size();
998  std::ostringstream oss;
999  for (unsigned int i = 0; i < num_images; i++)
1000  {
1001  std::list<typename TImage::Pointer> image_list;
1002  std::list<typename TMask::ConstPointer> mask_list;
1003  std::list<base_transform_t::ConstPointer> transform_list;
1004 
1005  image_list.push_back(image[i]);
1006  mask_list.push_back(mask[i]);
1007  transform_list.push_back((identity_transform_t::New()).GetPointer());
1008 
1009  oss << std::endl;
1010  for (unsigned int j = 0; j < num_images; j++)
1011  {
1012  if (i == j) continue;
1013  if (path[i][j].GetPointer() == NULL) continue;
1014 
1015  oss << "found a mapping: " << std::setw(2) << j << " -> "
1016  << std::setw(2) << i << std::endl;
1017 
1018  image_list.push_back(image[j]);
1019  mask_list.push_back(mask[j]);
1020  transform_list.push_back(path[i][j].GetPointer());
1021  }
1022 
1023  CORE_LOG_MESSAGE(oss.str());
1024 
1025  std::vector<typename TImage::Pointer>
1026  images(image_list.begin(), image_list.end());
1027 
1028  std::vector<typename TMask::ConstPointer>
1029  masks(mask_list.begin(), mask_list.end());
1030 
1031  std::vector<base_transform_t::ConstPointer>
1032  transforms(transform_list.begin(), transform_list.end());
1033 
1034  typename TImage::Pointer mosaic[3];
1035  make_mosaic_rgb(mosaic,
1036  std::vector<bool>(images.size(), false),
1037  transforms,
1038  images,
1039  masks,
1040  FEATHER_NONE_E,
1041  255,
1042  true);
1043 
1044  std::ostringstream fn_save;
1045  fn_save << prefix << "neighbors-" << the_text_t::number(i, 2, '0') <<
1046  "-pass-" << the_text_t::number(pass) << ".png";
1047  save_rgb<typename TImage::Pointer>(mosaic, fn_save.str(), true);
1048  }
1049 }
1050 
1051 
1052 //----------------------------------------------------------------
1053 // calc_area_and_dist
1054 //
1055 template <class TImage, class transform_ptr_t>
1056 void
1057 calc_area_and_dist(const std::vector<typename TImage::Pointer> & image,
1058  const array2d(transform_ptr_t) & mapping,
1059  array2d(double) & overlap,
1060  array2d(double) & distance)
1061 {
1062  const unsigned int num_images = image.size();
1063  resize(overlap, num_images, num_images);
1064  resize(distance, num_images, num_images);
1065 
1066  for (unsigned int i = 0; i < num_images; i++)
1067  {
1068  if (image[i].GetPointer() == NULL) continue;
1069 
1070  // calculate the tile center:
1071  pnt2d_t ci;
1072  {
1073  typename TImage::SizeType sz =
1074  image[i]->GetLargestPossibleRegion().GetSize();
1075 
1076  typename TImage::SpacingType sp =
1077  image[i]->GetSpacing();
1078 
1079  pnt2d_t tile_min = image[i]->GetOrigin();
1080  pnt2d_t tile_max;
1081  tile_max[0] = tile_min[0] + sp[0] * static_cast<double>(sz[0]);
1082  tile_max[1] = tile_min[1] + sp[1] * static_cast<double>(sz[1]);
1083 
1084  ci = pnt2d(tile_min[0] + (tile_max[0] - tile_min[0]) / 2.0,
1085  tile_min[1] + (tile_max[1] - tile_min[1]) / 2.0);
1086  }
1087 
1088  overlap[i][i] = 1.0;
1089  distance[i][i] = 0.0;
1090 
1091  for (unsigned int j = i + 1; j < num_images; j++)
1092  {
1093  if (mapping[i][j].GetPointer() == NULL) continue;
1094  overlap[i][j] = overlap_ratio<TImage>(image[i], image[j],
1095  mapping[i][j]);
1096  overlap[j][i] = overlap[i][j];
1097 
1098  pnt2d_t cj;
1099  {
1100  typename TImage::SizeType sz =
1101  image[j]->GetLargestPossibleRegion().GetSize();
1102 
1103  typename TImage::SpacingType sp =
1104  image[j]->GetSpacing();
1105 
1106  pnt2d_t tile_min = image[j]->GetOrigin();
1107  pnt2d_t tile_max;
1108  tile_max[0] = tile_min[0] + sp[0] * static_cast<double>(sz[0]);
1109  tile_max[1] = tile_min[1] + sp[1] * static_cast<double>(sz[1]);
1110 
1111  pnt2d_t center =
1112  pnt2d(tile_min[0] + (tile_max[0] - tile_min[0]) / 2.0,
1113  tile_min[1] + (tile_max[1] - tile_min[1]) / 2.0);
1114 
1115  cj = mapping[j][i]->TransformPoint(center);
1116  }
1117 
1118  distance[i][j] = (ci - cj).GetNorm();
1119  distance[j][i] = distance[i][j];
1120  }
1121  }
1122 }
1123 
1124 
1125 //----------------------------------------------------------------
1126 // local_mapping_t
1127 //
1129 {
1130  public:
1131  local_mapping_t(const unsigned int & from = UINT_MAX,
1132  const unsigned int & to = UINT_MAX,
1133  const double & dist = std::numeric_limits<double>::max()):
1134  from_(from),
1135  to_(to),
1136  distance_(dist)
1137  {}
1138 
1139  inline bool operator < (const local_mapping_t & m) const
1140  {
1141  if (to_ == m.to_) return distance_ < m.distance_;
1142  return to_ < m.to_;
1143  }
1144 
1145  inline bool operator == (const local_mapping_t & m) const
1146  { return to_ == m.to_; }
1147 
1148  unsigned int from_;
1149  unsigned int to_;
1150  double distance_;
1151 };
1152 
1153 
1154 //----------------------------------------------------------------
1155 // neighbor_t
1156 //
1158 {
1159  public:
1160  neighbor_t():
1161  id_(UINT_MAX),
1162  metric_(std::numeric_limits<double>::max())
1163  {}
1164 
1165  neighbor_t(const unsigned int & id,
1166  const double & metric,
1167  translate_transform_t::Pointer & t):
1168  id_(id),
1169  metric_(metric),
1170  t_(t)
1171  {}
1172 
1173  inline bool operator == (const neighbor_t & d) const
1174  { return id_ == d.id_; }
1175 
1176  inline bool operator < (const neighbor_t & d) const
1177  { return metric_ < d.metric_; }
1178 
1179  unsigned int id_;
1180  double metric_;
1181  translate_transform_t::Pointer t_;
1182 };
1183 
1184 //----------------------------------------------------------------
1187 template <class TImage, class TMask>
1188 bool
1189 match_pairs(std::vector<typename TImage::Pointer> & image,
1190  std::vector<typename TMask::ConstPointer> & mask,
1191  const bool & images_were_resampled,
1192  const bool & use_std_mask,
1193 
1194  std::vector<bool> & matched,
1195  array2d(translate_transform_t::Pointer) & path,
1196  array2d(double) & cost,
1197 
1198  const double & overlap_min,
1199  const double & overlap_max,
1200 
1201  image_t::PointType offset_min,
1202  image_t::PointType offset_max,
1203 
1204  double min_peak = 0.1,
1205  double peak_threshold = 0.3)
1206 {
1207  // shortcut:
1208  const unsigned int num_images = image.size();
1209  if (num_images < 2) return false;
1210 
1211  // initialize the mappings and mapping costs:
1212  for (unsigned int i = 0; i < num_images; i++)
1213  {
1214  matched[i] = false;
1215  for (unsigned int j = 0; j < num_images; j++)
1216  {
1217  path[i][j] = NULL;
1218  cost[i][j] = std::numeric_limits<double>::max();
1219  }
1220  }
1221 
1222  // this map is used to match the mismatched images:
1223  // index = image id
1224  // first = neighbor image id
1225  // second = list of peaks in the corresponding displacement PDF
1226  std::vector<std::pair<unsigned int, std::list<local_max_t> > >
1227  peaks(num_images);
1228 
1229  std::list<unsigned int> perimeter;
1230  bool perimeter_seeded = false;
1231  std::ostringstream oss;
1232  for (unsigned int seed = 0;
1233  seed < num_images && (!perimeter_seeded || !perimeter.empty());
1234  seed++)
1235  {
1236  // bootstrap the perimeter:
1237  if (!perimeter_seeded) push_back_unique(perimeter, seed);
1238 
1239  // try to match the perimeter to the unmatched images; as new matches
1240  // are discovered they are added to the perimeter:
1241  while (!perimeter.empty())
1242  {
1243  oss << "- - - - - - - - - - - - - - - - - - - - - - - - - - -" << std::endl;
1244 
1245  // consider a node from the perimeter:
1246  unsigned int node_id = remove_head(perimeter);
1247  matched[node_id] = true;
1248 
1249  // match it with the remaining unmatched nodes:
1250  std::list<neighbor_t> md;
1251  unsigned int md_size = 0;
1252  for (unsigned int i = 0; i < matched.size(); i++)
1253  {
1254  if (matched[i]) continue;
1255 
1256  oss << std::endl << "matching image " << node_id << " to " << i << ", "
1257  << std::flush;
1258 
1259  const TImage * fi = image[node_id];
1260  const TImage * mi = image[i];
1261  const TMask * fi_mask = mask[node_id];
1262  const TMask * mi_mask = mask[i];
1263 
1264  neighbor_t best_match;
1265  best_match.id_ = i;
1266  best_match.metric_ =
1267  match_one_pair<TImage, TMask>(images_were_resampled,
1268  use_std_mask,
1269  fi,
1270  mi,
1271  fi_mask,
1272  mi_mask,
1273  overlap_min,
1274  overlap_max,
1275  offset_min,
1276  offset_max,
1277  node_id,
1278  best_match.t_,
1279  peaks[i]);
1280 
1281  if (best_match.metric_ != std::numeric_limits<double>::max())
1282  {
1283  md.push_back(best_match);
1284  md_size++;
1285  }
1286  }
1287 
1288  // pick the top few matches from the sea of possible matches:
1289  md.sort();
1290 
1291  std::list<unsigned int> next_perimeter;
1292 
1293  if (md_size != 0)
1294  {
1295  oss << std::endl << "METRICS: " << std::flush;
1296  }
1297  else
1298  {
1299  oss << std::endl << "WARNING: " << node_id
1300  << " did not match any other tile" << std::flush;
1301  }
1302 
1303  bool separator = false;
1304  for (std::list<neighbor_t>::iterator i = md.begin();
1305  i != md.end(); ++i)
1306  {
1307  const neighbor_t & best = *(md.begin());
1308  const neighbor_t & data = *i;
1309 
1310  double decay = (1.0 + data.metric_) / (1.0 + best.metric_);
1311  if (decay > 4.0)
1312  {
1313  if (!separator)
1314  {
1315  oss << ".... ";
1316  separator = true;
1317  }
1318  }
1319  else
1320  {
1321  next_perimeter.push_front(data.id_);
1322 
1323  path[node_id][data.id_] = data.t_;
1324  cost[node_id][data.id_] = data.metric_;
1325  path[data.id_][node_id] = inverse(data.t_);
1326  cost[data.id_][node_id] = data.metric_;
1327  }
1328 
1329  oss << data.metric_ << ' ' << std::flush;
1330  }
1331  oss << std::endl;
1332 
1333  // expand the perimeter (elements must remain unique):
1334  expand(perimeter, next_perimeter, true);
1335 
1336  // clear the perimeter seeding flag:
1337  perimeter_seeded = perimeter_seeded || (next_perimeter.size() != 0);
1338  }
1339 
1340  // don't try to match the unmatched images
1341  // until we have a decent perimeter:
1342  if (!perimeter_seeded) continue;
1343 
1344  // try to match the unmatched images:
1345  for (unsigned int i = 0; i < num_images; i++)
1346  {
1347  if (matched[i]) continue;
1348 
1349  unsigned int node_id = peaks[i].first;
1350  const TImage * fi = image[node_id];
1351  const TImage * mi = image[i];
1352  const TMask * fi_mask = mask[node_id];
1353  const TMask * mi_mask = mask[i];
1354 
1355  // choose the best peak:
1356  neighbor_t best_match;
1357  for (std::list<local_max_t>::iterator j = peaks[i].second.begin();
1358  j != peaks[i].second.end(); ++j)
1359  {
1360  oss << "evaluating permutations..." << std::endl;
1361  const local_max_t & lm = *j;
1362 
1363  image_t::PointType offset_min;
1364  offset_min[0] = -std::numeric_limits<double>::max();
1365  offset_min[1] = -std::numeric_limits<double>::max();
1366  image_t::PointType offset_max;
1367  offset_max[0] = std::numeric_limits<double>::max();
1368  offset_max[1] = std::numeric_limits<double>::max();
1369 
1370  translate_transform_t::Pointer ti = translate_transform_t::New();
1371  double metric = estimate_displacement<TImage>(fi,
1372  mi,
1373  lm,
1374  ti,
1375  offset_min,
1376  offset_max,
1377  overlap_min,
1378  overlap_max,
1379  fi_mask,
1380  mi_mask);
1381  if (metric > best_match.metric_) continue;
1382 
1383  best_match = neighbor_t(i, metric, ti);
1384  }
1385  oss << std::endl;
1386  if (best_match.id_ == UINT_MAX) continue;
1387  oss << "METRIC: " << best_match.metric_ << std::endl;
1388 
1389  perimeter.push_front(best_match.id_);
1390 
1391  path[node_id][best_match.id_] = best_match.t_;
1392  cost[node_id][best_match.id_] = best_match.metric_;
1393  path[best_match.id_][node_id] = inverse(best_match.t_);
1394  cost[best_match.id_][node_id] = best_match.metric_;
1395  }
1396  }
1397 
1398  CORE_LOG_MESSAGE(oss.str());
1399  return perimeter_seeded;
1400 }
1401 
1402 //----------------------------------------------------------------
1403 // match_pairs (strategy version)
1404 //
1405 template <class TImage, class TMask>
1406 bool
1407 match_pairs_strategy(std::vector<typename TImage::Pointer> & image,
1408  std::vector<typename TMask::ConstPointer> & mask,
1409  const bool & images_were_resampled,
1410  const bool & use_std_mask,
1411 
1412  std::vector<bool> & matched,
1413  array2d(translate_transform_t::Pointer) & path,
1414  array2d(double) & cost,
1415 
1416  const double & overlap_min,
1417  const double & overlap_max,
1418 
1419  image_t::PointType offset_min,
1420  image_t::PointType offset_max,
1421 
1422  size_t shrink,
1423  int edgeLen,
1424  const StrategyType strategy,
1425 
1426  double min_peak = 0.1,
1427  double peak_threshold = 0.3)
1428 {
1429 #ifndef NDEBUG
1430 if (strategy == DEFAULT)
1431  CORE_LOG_DEBUG("match_pairs_strategy: DEFAULT");
1432 else if (strategy == TOP_LEFT_BOOK)
1433  CORE_LOG_DEBUG("match_pairs_strategy: TOP_LEFT_BOOK");
1434 else if (strategy == TOP_LEFT_SNAKE)
1435  CORE_LOG_DEBUG("match_pairs_strategy: TOP_LEFT_SNAKE");
1436 else
1437  CORE_LOG_DEBUG("match_pairs_strategy: UNKNOWN");
1438 #endif
1439 
1440  // shortcut:
1441  const unsigned int num_images = image.size();
1442  if (num_images < 2) return false;
1443 
1444  typename TImage::IndexValueType numRows =
1445  (typename TImage::IndexValueType)
1446  (image[0]->GetLargestPossibleRegion().GetSize()[0]);
1447  typename TImage::IndexValueType numCols =
1448  (typename TImage::IndexValueType)
1449  (image[0]->GetLargestPossibleRegion().GetSize()[1]);
1450  std::vector<std::pair<unsigned int, std::list<local_max_t> > >
1451  peaks(num_images);
1452 
1453  std::ostringstream oss;
1454  for (int row = 0; row < edgeLen; row++) {
1455  for (int column = 0; column < edgeLen; column++) {
1456  oss << "- - - - - - - - - - - - - - - - - - - - - -" << std::endl;
1457  // match it with the surrounding nodes:
1458  int seed = column + row * edgeLen;
1459  if (seed > num_images) break;
1460  matched[seed] = true;
1461 
1462  for (int i = -1; i <= 1; i++) {
1463  for (int j = -1; j <= 1; j++) {
1464  //determine the shift based off of offsets and rows.
1465  if ((i == 0 && j == 0) ||
1466  i + row < 0 || i + row >= edgeLen ||
1467  j + column < 0 || j + column >= edgeLen)
1468  continue;
1469 
1470  //forget the diagonal neighbors
1471  if (i && j) continue;
1472 
1473  //get the neighbor
1474  int neighbor = seed + edgeLen * i + j;
1475  if (strategy == TOP_LEFT_SNAKE)
1476  {
1477  int col = neighbor % edgeLen;
1478  int rw = ((double)neighbor) / ((double)edgeLen);
1479  if ((rw % 2) == 1)
1480  col = edgeLen - col - 1;
1481  neighbor = col + rw * edgeLen;
1482  }
1483 
1484  if (neighbor > num_images || neighbor < 0)
1485  continue;
1486 
1487  std::cout << "matching image " <<
1488  seed << " to " << neighbor << ". " << std::endl;
1489  //oss << "matching image " <<
1490  // seed << " to " << neighbor << ". " << std::endl;
1491 
1492  const TImage * fi = image[seed];
1493  const TImage * mi = image[neighbor];
1494  const TMask * fi_mask = mask[seed];
1495  const TMask * mi_mask = mask[neighbor];
1496 
1497  neighbor_t best_match;
1498  best_match.id_ = neighbor;
1499  best_match.metric_ = match_one_pair<TImage, TMask>(
1500  images_were_resampled,
1501  use_std_mask,
1502  fi,
1503  mi,
1504  fi_mask,
1505  mi_mask,
1506  overlap_min,
1507  overlap_max,
1508  offset_min,
1509  offset_max,
1510  seed,
1511  best_match.t_,
1512  peaks[neighbor],
1513  shrink,
1514  min_peak,
1515  peak_threshold,
1516  static_cast<size_t>((i+1) + ((j+1)*4)));
1517 
1518  if (best_match.metric_ ==
1519  std::numeric_limits<double>::max()) {
1520  //set default translations to minimum with low priority
1521  translate_transform_t::Pointer t = translate_transform_t::New();
1522  typedef itk::TranslationTransform<double,2>::OutputVectorType VectorType;
1523  VectorType v;
1524  double rowLen = static_cast<double>(numRows) * shrink;
1525  double colLen = static_cast<double>(numCols) * shrink;
1526  v[0] = rowLen * overlap_min * ((column<edgeLen/2)?-1.:1.);
1527  v[1] = colLen * overlap_min * ((row <edgeLen/2)?-1.:1.);
1528  if (j == -1)
1529  v[0] = rowLen * (1. - overlap_min);
1530  else if (j == 1)
1531  v[0] = rowLen * (overlap_min - 1.);
1532  if (i == -1)
1533  v[1] = colLen * (1. - overlap_min);
1534  else if (i == 1)
1535  v[1] = colLen * (overlap_min - 1.);
1536  // Force the same overlaps as adjacent tiles
1537  // This only works if were at row > 0 and column > 0
1538  if (column > 0 ) {
1539  if (column + j > 0) {
1540  v[1] = path[seed - 1][neighbor - 1].
1541  GetPointer()->GetOffset()[1];
1542  // std::cout << "aligning y with: " << (seed-1)
1543  // << " and " << (neighbor-1) << ", ";
1544  } else if (seed - 1 != neighbor) {
1545  v[1] = path[seed - 1][neighbor].
1546  GetPointer()->GetOffset()[1];
1547  // std::cout << "aligning y with: " << (seed-1)
1548  // << " and " << neighbor << std::endl;
1549  //} else {
1550  // std::cout << "Going with default min y, " ;
1551  }
1552  //} else {
1553  // std::cout << "Going with default min y, ";
1554  }
1555  if (row > 0) {
1556  if (row + i > 0) {
1557  v[0] = path[seed - edgeLen][neighbor - edgeLen].
1558  GetPointer()->GetOffset()[0];
1559  // std::cout << "aligning x with: " << (seed-edgeLen)
1560  // << " and " << (neighbor-edgeLen) << std::endl;
1561  } else if (seed - edgeLen != neighbor) {
1562  v[0] = path[seed - edgeLen][neighbor].
1563  GetPointer()->GetOffset()[0];
1564  // std::cout << "aligning x with: " << (seed-edgeLen)
1565  // << " and " << (neighbor) << std::endl;
1566  //} else {
1567  // std::cout << "Going with default min x." << std::endl;
1568  }
1569  //} else {
1570  // std::cout << "Going with default min x." << std::endl;
1571  }
1572  t.GetPointer()->SetOffset(v);
1573  local_max_t lm(0.1,v[0],v[1],1);
1574  best_match.metric_ = estimate_displacement<TImage>(
1575  fi,
1576  mi,
1577  lm,
1578  t,
1579  offset_min,
1580  offset_max,
1581  overlap_min,
1582  overlap_max,
1583  fi_mask,
1584  mi_mask);;
1585  best_match.t_ = t;
1586  //} else {
1587  // std::cout << "Successful normal match." << std::endl;
1588  }
1589  path[seed][neighbor] = best_match.t_;
1590  cost[seed][neighbor] = best_match.metric_;
1591  //std::cout << "with x: " << path[seed][neighbor].GetPointer()->
1592  // GetOffset()[0] << " and y: " << path[seed][neighbor].
1593  // GetPointer()->GetOffset()[1] << std::endl;
1594  if (cost[neighbor][seed] > best_match.metric_) {
1595  path[neighbor][seed] = inverse(best_match.t_);
1596  cost[neighbor][seed] = best_match.metric_;
1597  }
1598  }
1599  }
1600  }
1601  }
1602  CORE_LOG_MESSAGE(oss.str());
1603  return true;
1604 }
1605 
1606 
1607 //----------------------------------------------------------------
1608 // layout_mosaic
1609 //
1610 template <class TImage, class TMask>
1611 void
1612 layout_mosaic(// multi-resolution image tiles and tile masks,
1613  // the pyramids are sorted low-res first to high-res last:
1614  const unsigned int & num_tiles,
1615  array2d(typename TImage::Pointer) & tile_pyramid,
1616  array2d(typename TMask::ConstPointer) & mask_pyramid,
1617 
1618  // maximum number of intermediate tiles in cascaded mappings:
1619  const unsigned int & max_cascade_len,
1620 
1621  // the final transforms associated with each tile:
1622  std::vector<affine_transform_t::Pointer> & affine,
1623 
1624  // a vector identifying which tiles have been successfully matched:
1625  std::vector<bool> & matched,
1626 
1627  // matched tile layout order:
1628  std::list<unsigned int> & tile_order,
1629 
1630  // cost[num_intermediate_tiles][i][j] contains the maximum error metric
1631  // due to the image mapping between each image in the cascading chain:
1632  array3d(double) & cost,
1633 
1634  // path[num_intermediate_tiles][i][j] contains the transforms
1635  // associated with the specified number of intermediate images
1636  // that maps image i to image j:
1637  array3d(translate_transform_t::Pointer) & path,
1638 
1639  // image matching parameters:
1640  bool images_were_resampled,
1641  bool use_std_mask,
1642 
1643  // min/max tile overlap parameters:
1644  double overlap_min,
1645  double overlap_max,
1646 
1647  image_t::PointType offset_min,
1648  image_t::PointType offset_max,
1649 
1650  double min_peak,
1651  double peak_threshold,
1652 
1653  size_t shrink,
1654  int & edgeLen,
1655  const StrategyType strategy = DEFAULT,
1656 
1657  bool try_refining = false,
1658  const bfs::path & prefix = "")
1659 {
1660  typedef itk::LinearInterpolateImageFunction<TImage, double> interpolator_t;
1661 
1662  static const unsigned int iterations_per_level = 8;
1663  unsigned int pyramid_levels = tile_pyramid.size();
1664  unsigned int high_res_level = pyramid_levels - 1;
1665 
1666  bool perimeter_seeded = false;
1667 
1668  if (strategy == DEFAULT)
1669  perimeter_seeded =
1670  match_pairs<TImage, TMask>(tile_pyramid[high_res_level],
1671  mask_pyramid[high_res_level],
1672  images_were_resampled,
1673  use_std_mask,
1674  matched,
1675  path[0],
1676  cost[0],
1677  overlap_min,
1678  overlap_max,
1679  offset_min,
1680  offset_max,
1681  0.,
1682  peak_threshold);
1683  else
1684  perimeter_seeded =
1685  match_pairs_strategy<TImage, TMask>(tile_pyramid[high_res_level],
1686  mask_pyramid[high_res_level],
1687  images_were_resampled,
1688  use_std_mask,
1689  matched,
1690  path[0],
1691  cost[0],
1692  overlap_min,
1693  overlap_max,
1694  offset_min,
1695  offset_max,
1696  shrink,
1697  edgeLen,
1698  strategy,
1699  min_peak,
1700  peak_threshold);
1701 
1702  std::ostringstream oss;
1703  if (! perimeter_seeded)
1704  {
1705  oss << "WARNING: none of the tiles matched..." << std::endl;
1706  std::cout << "WARNING: none of the tiles matched..." << std::endl;
1707  CORE_LOG_MESSAGE(oss.str());
1708  return;
1709  }
1710 
1711  // try to improve the mappings as much as possible:
1712  // assemble cascaded transform chains with a prescribed maximum number of
1713  // intermediate nodes, try to choose the chain that is the most stable:
1714  if (strategy == DEFAULT)
1715  assemble_cascades(num_tiles, max_cascade_len, path, cost, false);
1716 
1717  // collapse the graph into a mapping:
1718  array2d(translate_transform_t::Pointer) mapping;
1719  resize(mapping, num_tiles, num_tiles);
1720 
1721  array2d(double) mapping_cost;
1722  resize(mapping_cost, num_tiles, num_tiles);
1723 
1724  establish_mappings(num_tiles,
1725  max_cascade_len,
1726  path,
1727  cost,
1728  mapping,
1729  mapping_cost);
1730 
1731  if (strategy == DEFAULT)
1732  {
1733  reset(num_tiles, max_cascade_len, path, cost);
1734 
1735  brute_force_pairs<TImage, TMask>(tile_pyramid,
1736  mask_pyramid,
1737  overlap_min,
1738  overlap_max,
1739  7, // dx
1740  7, // dy
1741  2, // pyramid_levels, // refinement levels
1742  mapping,
1743  path[0],
1744  cost[0],
1745  prefix);
1746  }
1747 
1748  if (! prefix.empty())
1749  {
1750  dump_neighbors<TImage, TMask>(tile_pyramid[high_res_level],
1751  mask_pyramid[high_res_level],
1752  path[0],
1753  prefix);
1754  }
1755 
1756  // find local overlaps and distances between overlapping tiles:
1757  array2d(double) local_overlap;
1758  array2d(double) local_distance;
1759  calc_area_and_dist<TImage, translate_transform_t::Pointer>
1760  (tile_pyramid[high_res_level], path[0], local_overlap, local_distance);
1761 
1762  unsigned int best_target = ~0;
1763 
1764  //#if 1
1765  double best_overlap = 0.0;
1766  for (unsigned int i = 0; i < num_tiles; i++)
1767  {
1768  double overlap_sum = 0.0;
1769  for (unsigned int j = 0; j < num_tiles; j++)
1770  {
1771  if (j == i) continue;
1772 
1773  overlap_sum += local_overlap[i][j];
1774  }
1775 
1776  if (overlap_sum > best_overlap)
1777  {
1778  best_target = i;
1779  best_overlap = overlap_sum;
1780  }
1781  }
1782  // sort the tiles according to their distance to the target tile:
1783  std::list<local_mapping_t> perimeter;
1784  perimeter.push_back(local_mapping_t(best_target, best_target, 0.0));
1785 
1786  // setup a list of remaining tiles:
1787  std::list<unsigned int> tiles;
1788  for (unsigned int i = 0; i < num_tiles; i++)
1789  {
1790  if (i == best_target) continue;
1791  if (tile_pyramid[high_res_level][i].GetPointer() == NULL) continue;
1792 
1793  tiles.push_back(i);
1794  }
1795 
1796  // setup the initial transform:
1797  affine.resize(num_tiles);
1798  {
1799  identity_transform_t::Pointer identity = identity_transform_t::New();
1800  affine[best_target] =
1801  setup_transform<TImage>(tile_pyramid[high_res_level][best_target],
1802  identity, 0);
1803  }
1804 
1805  // setup the metric transform parameter masks:
1806  std::vector<bool> param_shared;
1807  std::vector<bool> param_active;
1808  affine_transform_t::setup_shared_params_mask(false, param_shared);
1809  param_active.assign(param_shared.size(), true);
1810 
1811  // layout and refine the mosaic incrementally as a flood fill:
1812  tile_order.clear();
1813  while (true)
1814  {
1815  std::list<local_mapping_t> candidates;
1816 
1817  // look at the current perimeter and generate the next perimeter:
1818  while (!perimeter.empty())
1819  {
1820  const unsigned int i = remove_head(perimeter).to_;
1821  tile_order.push_back(i);
1822 
1823  // look at the remaining tiles as the candidates for the next perimeter:
1824  for (std::list<unsigned int>::const_iterator jter = tiles.begin();
1825  jter != tiles.end(); ++jter)
1826  {
1827  const unsigned int j = *jter;
1828  if (path[0][i][j].GetPointer() == NULL) continue;
1829 
1830  local_mapping_t m(i, j, cost[0][i][j]);
1831  candidates.push_back(m);
1832  }
1833  }
1834 
1835  if (candidates.empty())
1836  {
1837  // done:
1838  break;
1839  }
1840 
1841  // remove duplicate mappings to the same tile, assemble the next perimeter:
1842  candidates.sort();
1843  perimeter.push_back(remove_head(candidates));
1844  tiles.remove(perimeter.back().to_);
1845 
1846  while (!candidates.empty())
1847  {
1848  const local_mapping_t & prev = perimeter.back();
1849  local_mapping_t next = remove_head(candidates);
1850  if (prev.to_ == next.to_) continue;
1851 
1852  perimeter.push_back(next);
1853  tiles.remove(next.to_);
1854  }
1855 
1856  // number of tiles that will be optimized in the next step:
1857  const unsigned int num_optimize = tile_order.size() + perimeter.size();
1858 
1859  // try incremental refinement of the mosaic:
1861  mosaic_metric_t;
1862 
1863  typename mosaic_metric_t::Pointer mosaic_metric = mosaic_metric_t::New();
1864  mosaic_metric->image_.resize(num_optimize);
1865  mosaic_metric->mask_.resize(num_optimize);
1866  mosaic_metric->transform_.resize(num_optimize);
1867 
1868  std::vector<unsigned int> tile_id(num_optimize);
1869 
1870  unsigned int i = 0;
1871  for (std::list<unsigned int>::const_iterator iter = tile_order.begin();
1872  iter != tile_order.end(); ++iter, i++)
1873  {
1874  tile_id[i] = *iter;
1875  mosaic_metric->image_[i] = tile_pyramid[high_res_level][tile_id[i]];
1876  mosaic_metric->mask_[i] = mask_pyramid[high_res_level][tile_id[i]];
1877  mosaic_metric->transform_[i] = affine[tile_id[i]].GetPointer();
1878  }
1879 
1880  for (std::list<local_mapping_t>::const_iterator iter = perimeter.begin();
1881  iter != perimeter.end(); ++iter, i++)
1882  {
1883  const local_mapping_t & m = *iter;
1884  tile_id[i] = m.to_;
1885  affine[tile_id[i]] =
1886  setup_transform<TImage>(tile_pyramid[high_res_level][tile_id[i]],
1887  affine[m.from_],
1888  path[0][m.from_][tile_id[i]]);
1889 
1890  mosaic_metric->image_[i] = tile_pyramid[high_res_level][tile_id[i]];
1891  mosaic_metric->mask_[i] = mask_pyramid[high_res_level][tile_id[i]];
1892  mosaic_metric->transform_[i] = affine[tile_id[i]].GetPointer();
1893  }
1894 
1895  if (!try_refining)
1896  {
1897  continue;
1898  }
1899 
1900  mosaic_metric->setup_param_map(param_shared, param_active);
1901  mosaic_metric->Initialize();
1902 
1903  // setup the optimizer scales:
1904  typename mosaic_metric_t::params_t parameter_scales =
1905  mosaic_metric->GetTransformParameters();
1906  parameter_scales.Fill(1.0);
1907 
1908  // discourage skewing and rotation:
1909  for (unsigned int k = 0; k < num_optimize; k++)
1910  {
1911  const unsigned int * addr = &(mosaic_metric->address_[k][0]);
1912 
1913  parameter_scales[addr[affine_transform_t::index_a(0, 1)]] = 1e+6;
1914  parameter_scales[addr[affine_transform_t::index_b(1, 0)]] = 1e+6;
1915  }
1916 
1917  for (unsigned int i = 0; i < num_optimize; i++)
1918  {
1919  mosaic_metric->image_[i] = tile_pyramid[high_res_level][tile_id[i]];
1920  }
1921 
1922  typename mosaic_metric_t::measure_t metric_before =
1923  mosaic_metric->GetValue(mosaic_metric->GetTransformParameters());
1924 
1925  // run several iterations of the optimizer:
1926  for (unsigned int k = 0; k < 1; k++)
1927  {
1928  typename mosaic_metric_t::params_t params_before =
1929  mosaic_metric->GetTransformParameters();
1930 
1931  typename mosaic_metric_t::measure_t metric_after =
1932  std::numeric_limits<typename mosaic_metric_t::measure_t>::max();
1933 
1934  // use global refinement:
1935  optimizer_t::Pointer optimizer = optimizer_t::New();
1936  optimizer_observer_t<optimizer_t>::Pointer observer =
1938  optimizer->AddObserver(itk::IterationEvent(), observer);
1939  optimizer->SetMinimize(true);
1940  optimizer->SetNumberOfIterations(iterations_per_level);
1941  optimizer->SetMinimumStepLength(1e-12); // min_step
1942  optimizer->SetMaximumStepLength(1e-5); // max_step
1943  optimizer->SetGradientMagnitudeTolerance(1e-6);
1944  optimizer->SetRelaxationFactor(5e-1);
1945  optimizer->SetCostFunction(mosaic_metric);
1946  optimizer->SetInitialPosition(params_before);
1947  optimizer->SetScales(parameter_scales);
1948  optimizer->SetPickUpPaceSteps(5);
1949  optimizer->SetBackTracking(true);
1950 
1951  // refine the mosaic:
1952  try
1953  {
1954  oss << "\n" << k << ": refining affine transforms" << std::endl;
1955  optimizer->StartOptimization();
1956  }
1957  catch (itk::ExceptionObject & exception)
1958  {
1959  // oops:
1960  oss << "optimizer threw an exception:" << std::endl << exception.what() << std::endl;
1961  }
1962  mosaic_metric->SetTransformParameters(optimizer->GetBestParams());
1963  metric_after = optimizer->GetBestValue();
1964 
1965  typename mosaic_metric_t::params_t params_after =
1966  mosaic_metric->GetTransformParameters();
1967 
1968  oss << "before: METRIC = " << metric_before
1969  << ", PARAMS = " << params_before << std::endl
1970  << "after: METRIC = " << metric_after
1971  << ", PARAMS = " << params_after << std::endl;
1972 
1973  // quantify the improvement:
1974  double improvement = 1.0 - metric_after / metric_before;
1975  bool failed_to_improve = (metric_after - metric_before) >= 0.0;
1976  bool negligible_improvement = !failed_to_improve && (improvement < 1e-3);
1977 
1978  if (!failed_to_improve)
1979  {
1980  oss << "IMPROVEMENT: " << std::setw(3) << static_cast<int>(100.0 * improvement) << "%" << std::endl;
1981  }
1982 
1983  if (failed_to_improve)
1984  {
1985  oss << "NOTE: minimization failed, ignoring registration results..." << std::endl;
1986 
1987  // previous transform was better:
1988  mosaic_metric->SetTransformParameters(params_before);
1989  break;
1990  }
1991  else if (negligible_improvement)
1992  {
1993  oss << "NOTE: improvement is negligible..." << std::endl;
1994  break;
1995  }
1996 
1997  // avoid recalculating the same metric:
1998  metric_before = metric_after;
1999  }
2000  }
2001  CORE_LOG_MESSAGE(oss.str());
2002 }
2003 
2004 
2005 #endif // MOSAIC_LAYOUT_COMMON_HXX_
Computes mean pixel variance within the overlapping regions of a mosaic.
Definition: itkImageMosaicVarianceMetric.h:69
Definition: mosaic_layout_common.hxx:1128
Definition: optimizer_observer.hxx:44
Definition: itkLegendrePolynomialTransform.h:81
Definition: fft_common.hxx:77
Definition: mosaic_layout_common.hxx:1157