Fraxinus  17.12-rc4
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) 2008-2014, SINTEF Department of Medical Technology
5 All rights reserved.
6 
7 Redistribution and use in source and binary forms, with or without
8 modification, are permitted provided that the following conditions are met:
9 
10 1. Redistributions of source code must retain the above copyright notice,
11  this list of conditions and the following disclaimer.
12 
13 2. Redistributions in binary form must reproduce the above copyright notice,
14  this list of conditions and the following disclaimer in the documentation
15  and/or other materials provided with the distribution.
16 
17 3. Neither the name of the copyright holder nor the names of its contributors
18  may be used to endorse or promote products derived from this software
19  without specific prior written permission.
20 
21 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
22 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
24 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
25 FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
26 DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
27 SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
28 CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
29 OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
30 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
31 =========================================================================*/
33 
34 /*=========================================================================
35 
36  Program: Insight Segmentation & Registration Toolkit
37  Module: $RCSfile: IterativeClosestPoint1.cxx,v $
38  Language: C++
39  Date: $Date: 2009-12-04 03:25:22 $
40  Version: $Revision: 1.7 $
41 
42  Copyright (c) Insight Software Consortium. All rights reserved.
43  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
44 
45  This software is distributed WITHOUT ANY WARRANTY; without even
46  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
47  PURPOSE. See the above copyright notices for more information.
48 
49  =========================================================================*/
50 #ifdef _MSC_VER
51 #pragma warning ( disable : 4786 )
52 #endif
53 
54 // Software Guide : BeginLatex
55 //
56 // This example illustrates how to perform Iterative Closest Point (ICP)
57 // registration in ITK. The main class featured in this section is the
58 // \doxygen{EuclideanDistancePointMetric}.
59 //
60 // Software Guide : EndLatex
61 
62 // Software Guide : BeginCodeSnippet
63 #include "itkTranslationTransform.h"
64 #include "itkEuclideanDistancePointMetric.h"
65 #include "itkLevenbergMarquardtOptimizer.h"
66 #include "itkPointSet.h"
67 #include "itkPointSetToPointSetRegistrationMethod.h"
68 // Software Guide : EndCodeSnippet
69 
70 #include <iostream>
71 #include <fstream>
72 
73 class CommandIterationUpdate: public itk::Command
74 {
75 public:
77  typedef itk::Command Superclass;
78  typedef itk::SmartPointer<Self> Pointer;
79  itkNewMacro( Self );
80 
81 protected:
83 
84 public:
85  typedef itk::LevenbergMarquardtOptimizer OptimizerType;
86  typedef const OptimizerType * OptimizerPointer;
87 
88  void Execute(itk::Object *caller, const itk::EventObject & event)
89  {
90  Execute((const itk::Object *) caller, event);
91  }
92 
93  void Execute(const itk::Object * object, const itk::EventObject & event)
94  {
95  OptimizerPointer optimizer = dynamic_cast<OptimizerPointer>(object);
96 
97  if (!itk::IterationEvent().CheckEvent(&event))
98  {
99  return;
100  }
101 
102  std::cout << "iteration Value = " << optimizer->GetCachedValue() << std::endl;
103  std::cout << "iteration Position = " << optimizer->GetCachedCurrentPosition();
104  std::cout << std::endl << std::endl;
105 
106  }
107 
108 };
109 
111 {
112 public:
113  static const unsigned int Dimension = 3;
114 
115  typedef itk::PointSet<float, Dimension> PointSetType;
116 
117 // PointSetType::Pointer fixedPointSet = PointSetType::New();
118 // PointSetType::Pointer movingPointSet = PointSetType::New();
119 
120  typedef PointSetType::PointType PointType;
121  typedef PointSetType::PointsContainer PointsContainer;
122 
123  //-----------------------------------------------------------
124  // Set up the Metric
125  //-----------------------------------------------------------
126  typedef itk::EuclideanDistancePointMetric<PointSetType, PointSetType> MetricType;
127 
128  typedef MetricType::TransformType TransformBaseType;
129  typedef TransformBaseType::ParametersType ParametersType;
130  typedef TransformBaseType::JacobianType JacobianType;
131 
132  //-----------------------------------------------------------
133  // Set up a Transform
134  //-----------------------------------------------------------
135 
136  typedef itk::TranslationTransform<double, Dimension> TransformType;
137 
138  // Optimizer Type
139  typedef itk::LevenbergMarquardtOptimizer OptimizerType;
140 
141  // Registration Method
142  typedef itk::PointSetToPointSetRegistrationMethod<PointSetType, PointSetType> RegistrationType;
143 
144  // methods:
146 
147  bool registerPoints(std::vector<cx::Vector3D> ref, std::vector<cx::Vector3D> target);
148  bool registerPoints(PointSetType::Pointer fixedPointSet, PointSetType::Pointer movingPointSet);
149  PointSetType::Pointer toItk(std::vector<cx::Vector3D> input);
150  int test();
151 
152 };
153 
154 LandmarkTranslation::PointSetType::Pointer LandmarkTranslation::toItk(std::vector<cx::Vector3D> input)
155 {
156  PointSetType::Pointer pointSet = PointSetType::New();
157  PointsContainer::Pointer pointContainer = PointsContainer::New();
158  PointType point;
159 
160  for (unsigned i = 0; i < input.size(); ++i)
161  {
162  point[0] = input[i][0];
163  point[1] = input[i][1];
164  point[2] = input[i][2];
165  pointContainer->InsertElement(i, point);
166  }
167  pointSet->SetPoints(pointContainer);
168  return pointSet;
169 }
170 
171 bool LandmarkTranslation::registerPoints(std::vector<cx::Vector3D> ref, std::vector<cx::Vector3D> target)
172 {
173  return registerPoints(toItk(ref), toItk(target));
174 }
175 
176 bool LandmarkTranslation::registerPoints(PointSetType::Pointer fixedPointSet, PointSetType::Pointer movingPointSet)
177 {
178 
179 //-----------------------------------------------------------
180 // Set up the Metric
181 //-----------------------------------------------------------
182  MetricType::Pointer metric = MetricType::New();
183 
184 //-----------------------------------------------------------
185 // Set up a Transform
186 //-----------------------------------------------------------
187 
188  TransformType::Pointer transform = TransformType::New();
189 
190  // Optimizer Type
191  OptimizerType::Pointer optimizer = OptimizerType::New();
192  optimizer->SetUseCostFunctionGradient(false);
193 
194  // Registration Method
195  RegistrationType::Pointer registration = RegistrationType::New();
196 
197  // Scale the translation components of the Transform in the Optimizer
198  OptimizerType::ScalesType scales(transform->GetNumberOfParameters());
199  std::cout << "transform->GetNumberOfParameters() " << transform->GetNumberOfParameters() << std::endl;
200  scales.Fill(0.01);
201 
202  unsigned long numberOfIterations = 100;
203  double gradientTolerance = 1e-5; // convergence criterion
204  double valueTolerance = 1e-5; // convergence criterion
205  double epsilonFunction = 1e-6; // convergence criterion
206 
207  optimizer->SetScales(scales);
208  optimizer->SetNumberOfIterations(numberOfIterations);
209  optimizer->SetValueTolerance(valueTolerance);
210  optimizer->SetGradientTolerance(gradientTolerance);
211  optimizer->SetEpsilonFunction(epsilonFunction);
212 
213  // Start from an Identity transform (in a normal case, the user
214  // can probably provide a better guess than the identity...
215  transform->SetIdentity();
216 
217  registration->SetInitialTransformParameters(transform->GetParameters());
218 
219  //------------------------------------------------------
220  // Connect all the components required for Registration
221  //------------------------------------------------------
222  registration->SetMetric(metric);
223  registration->SetOptimizer(optimizer);
224  registration->SetTransform(transform);
225  registration->SetFixedPointSet(fixedPointSet);
226  registration->SetMovingPointSet(movingPointSet);
227 
228 // // Connect an observer
229 // CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
230 // optimizer->AddObserver( itk::IterationEvent(), observer );
231 
232  try
233  {
234 // registration->StartRegistration();
235  registration->Update();
236  } catch (itk::ExceptionObject & e)
237  {
238  std::cout << e << std::endl;
239  return false;
240  }
241 
242  mResult = cx::Transform3D::Identity();
243  for (unsigned i = 0; i < transform->GetNumberOfParameters(); ++i)
244  mResult(i, 3) = transform->GetParameters()[i];
245 
246  //std::cout << "Solution = " << transform->GetParameters() << "\n" << mResult << std::endl;
247  //std::cout << "Solution = " << registration->GetLastTransformParameters() << std::endl;
248 
249  return true;
250 
251 }
252 
253 namespace cx
254 {
255 
263 Transform3D LandmarkTranslationRegistration::registerPoints(std::vector<Vector3D> ref,
264  std::vector<Vector3D> target, bool* ok)
265 {
266  if (ref.size() != target.size() || ref.empty())
267  {
268  std::cout << "Different sizes in ref and target: aborting registration." << std::endl;
269  *ok = false;
270  return Transform3D::Identity();
271  }
272 
273  // ad-hoc solution for one and two points: itk doesn't handle this for some reason.
274  if (ref.size() == 1)
275  {
276  Vector3D t = ref[0] - target[0];
277  *ok = true;
278  return createTransformTranslate(t);
279  }
280  if (ref.size() == 2)
281  {
282  Vector3D rr = (ref[0] + ref[1]) / 2.0;
283  Vector3D tt = (target[0] + target[1]) / 2.0;
284  Vector3D t = rr - tt;
285  *ok = true;
286  return createTransformTranslate(t);
287  }
288 
289  LandmarkTranslation registrator;
290  *ok = registrator.registerPoints(ref, target);
291  return registrator.mResult;
292 }
293 
294 } // 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:63
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.