Fraxinus  2023.01.05-dev+develop.0da12
An IGT application
cxLandmarkTranslationRegistration.cpp
Go to the documentation of this file.
1 /*=========================================================================
2 This file is part of CustusX, an Image Guided Therapy Application.
3 
4 Copyright (c) SINTEF Department of Medical Technology.
5 All rights reserved.
6 
7 CustusX is released under a BSD 3-Clause license.
8 
9 See Lisence.txt (https://github.com/SINTEFMedtek/CustusX/blob/master/License.txt) for details.
10 =========================================================================*/
12 
13 /*=========================================================================
14 
15  Program: Insight Segmentation & Registration Toolkit
16  Module: $RCSfile: IterativeClosestPoint1.cxx,v $
17  Language: C++
18  Date: $Date: 2009-12-04 03:25:22 $
19  Version: $Revision: 1.7 $
20 
21  Copyright (c) Insight Software Consortium. All rights reserved.
22  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
23 
24  This software is distributed WITHOUT ANY WARRANTY; without even
25  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
26  PURPOSE. See the above copyright notices for more information.
27 
28  =========================================================================*/
29 #ifdef _MSC_VER
30 #pragma warning ( disable : 4786 )
31 #endif
32 
33 // Software Guide : BeginLatex
34 //
35 // This example illustrates how to perform Iterative Closest Point (ICP)
36 // registration in ITK. The main class featured in this section is the
37 // \doxygen{EuclideanDistancePointMetric}.
38 //
39 // Software Guide : EndLatex
40 
41 // Software Guide : BeginCodeSnippet
42 #include "itkTranslationTransform.h"
43 #include "itkEuclideanDistancePointMetric.h"
44 #include "itkLevenbergMarquardtOptimizer.h"
45 #include "itkPointSet.h"
46 #include "itkPointSetToPointSetRegistrationMethod.h"
47 // Software Guide : EndCodeSnippet
48 
49 #include <iostream>
50 #include <fstream>
51 
52 class CommandIterationUpdate: public itk::Command
53 {
54 public:
56  typedef itk::Command Superclass;
57  typedef itk::SmartPointer<Self> Pointer;
58  itkNewMacro( Self );
59 
60 protected:
62 
63 public:
64  typedef itk::LevenbergMarquardtOptimizer OptimizerType;
65  typedef const OptimizerType * OptimizerPointer;
66 
67  void Execute(itk::Object *caller, const itk::EventObject & event)
68  {
69  Execute((const itk::Object *) caller, event);
70  }
71 
72  void Execute(const itk::Object * object, const itk::EventObject & event)
73  {
74  OptimizerPointer optimizer = dynamic_cast<OptimizerPointer>(object);
75 
76  if (!itk::IterationEvent().CheckEvent(&event))
77  {
78  return;
79  }
80 
81  std::cout << "iteration Value = " << optimizer->GetCachedValue() << std::endl;
82  std::cout << "iteration Position = " << optimizer->GetCachedCurrentPosition();
83  std::cout << std::endl << std::endl;
84 
85  }
86 
87 };
88 
90 {
91 public:
92  static const unsigned int Dimension = 3;
93 
94  typedef itk::PointSet<float, Dimension> PointSetType;
95 
96 // PointSetType::Pointer fixedPointSet = PointSetType::New();
97 // PointSetType::Pointer movingPointSet = PointSetType::New();
98 
99  typedef PointSetType::PointType PointType;
100  typedef PointSetType::PointsContainer PointsContainer;
101 
102  //-----------------------------------------------------------
103  // Set up the Metric
104  //-----------------------------------------------------------
105  typedef itk::EuclideanDistancePointMetric<PointSetType, PointSetType> MetricType;
106 
107  typedef MetricType::TransformType TransformBaseType;
108  typedef TransformBaseType::ParametersType ParametersType;
109  typedef TransformBaseType::JacobianType JacobianType;
110 
111  //-----------------------------------------------------------
112  // Set up a Transform
113  //-----------------------------------------------------------
114 
115  typedef itk::TranslationTransform<double, Dimension> TransformType;
116 
117  // Optimizer Type
118  typedef itk::LevenbergMarquardtOptimizer OptimizerType;
119 
120  // Registration Method
121  typedef itk::PointSetToPointSetRegistrationMethod<PointSetType, PointSetType> RegistrationType;
122 
123  // methods:
125 
126  bool registerPoints(std::vector<cx::Vector3D> ref, std::vector<cx::Vector3D> target);
127  bool registerPoints(PointSetType::Pointer fixedPointSet, PointSetType::Pointer movingPointSet);
128  PointSetType::Pointer toItk(std::vector<cx::Vector3D> input);
129  int test();
130 
131 };
132 
133 LandmarkTranslation::PointSetType::Pointer LandmarkTranslation::toItk(std::vector<cx::Vector3D> input)
134 {
135  PointSetType::Pointer pointSet = PointSetType::New();
136  PointsContainer::Pointer pointContainer = PointsContainer::New();
137  PointType point;
138 
139  for (unsigned i = 0; i < input.size(); ++i)
140  {
141  point[0] = input[i][0];
142  point[1] = input[i][1];
143  point[2] = input[i][2];
144  pointContainer->InsertElement(i, point);
145  }
146  pointSet->SetPoints(pointContainer);
147  return pointSet;
148 }
149 
150 bool LandmarkTranslation::registerPoints(std::vector<cx::Vector3D> ref, std::vector<cx::Vector3D> target)
151 {
152  return registerPoints(toItk(ref), toItk(target));
153 }
154 
155 bool LandmarkTranslation::registerPoints(PointSetType::Pointer fixedPointSet, PointSetType::Pointer movingPointSet)
156 {
157 
158 //-----------------------------------------------------------
159 // Set up the Metric
160 //-----------------------------------------------------------
161  MetricType::Pointer metric = MetricType::New();
162 
163 //-----------------------------------------------------------
164 // Set up a Transform
165 //-----------------------------------------------------------
166 
167  TransformType::Pointer transform = TransformType::New();
168 
169  // Optimizer Type
170  OptimizerType::Pointer optimizer = OptimizerType::New();
171  optimizer->SetUseCostFunctionGradient(false);
172 
173  // Registration Method
174  RegistrationType::Pointer registration = RegistrationType::New();
175 
176  // Scale the translation components of the Transform in the Optimizer
177  OptimizerType::ScalesType scales(transform->GetNumberOfParameters());
178  std::cout << "transform->GetNumberOfParameters() " << transform->GetNumberOfParameters() << std::endl;
179  scales.Fill(0.01);
180 
181  unsigned long numberOfIterations = 100;
182  double gradientTolerance = 1e-5; // convergence criterion
183  double valueTolerance = 1e-5; // convergence criterion
184  double epsilonFunction = 1e-6; // convergence criterion
185 
186  optimizer->SetScales(scales);
187  optimizer->SetNumberOfIterations(numberOfIterations);
188  optimizer->SetValueTolerance(valueTolerance);
189  optimizer->SetGradientTolerance(gradientTolerance);
190  optimizer->SetEpsilonFunction(epsilonFunction);
191 
192  // Start from an Identity transform (in a normal case, the user
193  // can probably provide a better guess than the identity...
194  transform->SetIdentity();
195 
196  registration->SetInitialTransformParameters(transform->GetParameters());
197 
198  //------------------------------------------------------
199  // Connect all the components required for Registration
200  //------------------------------------------------------
201  registration->SetMetric(metric);
202  registration->SetOptimizer(optimizer);
203  registration->SetTransform(transform);
204  registration->SetFixedPointSet(fixedPointSet);
205  registration->SetMovingPointSet(movingPointSet);
206 
207 // // Connect an observer
208 // CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
209 // optimizer->AddObserver( itk::IterationEvent(), observer );
210 
211  try
212  {
213 // registration->StartRegistration();
214  registration->Update();
215  } catch (itk::ExceptionObject & e)
216  {
217  std::cout << e << std::endl;
218  return false;
219  }
220 
221  mResult = cx::Transform3D::Identity();
222  for (unsigned i = 0; i < transform->GetNumberOfParameters(); ++i)
223  mResult(i, 3) = transform->GetParameters()[i];
224 
225  //std::cout << "Solution = " << transform->GetParameters() << "\n" << mResult << std::endl;
226  //std::cout << "Solution = " << registration->GetLastTransformParameters() << std::endl;
227 
228  return true;
229 
230 }
231 
232 namespace cx
233 {
234 
242 Transform3D LandmarkTranslationRegistration::registerPoints(std::vector<Vector3D> ref,
243  std::vector<Vector3D> target, bool* ok)
244 {
245  if (ref.size() != target.size() || ref.empty())
246  {
247  std::cout << "Different sizes in ref and target: aborting registration." << std::endl;
248  *ok = false;
249  return Transform3D::Identity();
250  }
251 
252  // ad-hoc solution for one and two points: itk doesn't handle this for some reason.
253  if (ref.size() == 1)
254  {
255  Vector3D t = ref[0] - target[0];
256  *ok = true;
257  return createTransformTranslate(t);
258  }
259  if (ref.size() == 2)
260  {
261  Vector3D rr = (ref[0] + ref[1]) / 2.0;
262  Vector3D tt = (target[0] + target[1]) / 2.0;
263  Vector3D t = rr - tt;
264  *ok = true;
265  return createTransformTranslate(t);
266  }
267 
268  LandmarkTranslation registrator;
269  *ok = registrator.registerPoints(ref, target);
270  return registrator.mResult;
271 }
272 
273 } // namespace cx
itk::PointSet< float, Dimension > PointSetType
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
const unsigned int Dimension
void Execute(itk::Object *caller, const itk::EventObject &event)
void Execute(const itk::Object *object, const itk::EventObject &event)
itk::TranslationTransform< double, Dimension > TransformType
Transform3D createTransformTranslate(const Vector3D &translation)
MetricType::TransformType TransformBaseType
PointSetType::Pointer toItk(std::vector< cx::Vector3D > input)
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:42
TransformBaseType::ParametersType ParametersType
itk::PointSetToPointSetRegistrationMethod< PointSetType, PointSetType > RegistrationType
itk::LevenbergMarquardtOptimizer OptimizerType
PointSetType::PointsContainer PointsContainer
TransformBaseType::JacobianType JacobianType
cx::Transform3D mResult
transform of movingPointSet
bool registerPoints(std::vector< cx::Vector3D > ref, std::vector< cx::Vector3D > target)
itk::EuclideanDistancePointMetric< PointSetType, PointSetType > MetricType
float4 transform(float16 matrix, float4 voxel)
itk::LevenbergMarquardtOptimizer OptimizerType
Namespace for all CustusX production code.