KWStyle - itkLandmarkBasedTransformInitializer.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkLandmarkBasedTransformInitializer.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:40 $
7		Version:   $Revision: 1.4 $
8
9		Copyright (c) Insight Software Consortium. All rights reserved.
10		See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11
12		This software is distributed WITHOUT ANY WARRANTY; without even
13		the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14		PURPOSE.  See the above copyright notices for more information.
15
16		=========================================================================*/
17
18		#ifndef __itkLandmarkBasedTransformInitializer_txx
19		#define __itkLandmarkBasedTransformInitializer_txx
20
21		#include "itkLandmarkBasedTransformInitializer.h"
22		#include "itkMatrix.h"
23		#include "itkSymmetricEigenAnalysis.h"
24
25		namespace itk
26		{
27
28
29		template < class TTransform, class TFixedImage, class TMovingImage >
30		LandmarkBasedTransformInitializer< TTransform, TFixedImage, TMovingImage >
31		::LandmarkBasedTransformInitializer()
32		{
33		}
34
35
36
37		template < class TTransform, class TFixedImage, class TMovingImage >
38		void
39		LandmarkBasedTransformInitializer<TTransform, TFixedImage, TMovingImage >
40		::InitializeTransform()
41		{
42		// Sanity check
43		if( !m_Transform )
44		{
45		itkExceptionMacro( "Transform has not been set" );
46		return;
47		}
48		if( m_FixedLandmarks.size() != m_MovingLandmarks.size() )
49		{
50		itkExceptionMacro("Different number of fixed and moving landmarks");
51		return;
52		}
53
54
55
56		// We will do an explicit typeid check here (via dynamic_cast) to check
57		// the transform type. The initialization scheme will generally be different
58		// based on the transform type and the dimension. As more transforms are
59		// supported in future, an explicit typeid check is expected to be done here.
60	LEN	// Note that the typeid is done via dynamic_cast. This means that as more transforms
61		// are added in future, you will have to order your checks from the bottom
62		// of the transform hierarchy, upwards.
63		//
64		InputTransformType  transformType = Else;
65	LEN	VersorRigid3DTransformType *testPtr = dynamic_cast< VersorRigid3DTransformType *>(
66		this->m_Transform.GetPointer() );
67		if( testPtr )
68		{
69		transformType = VersorRigid3Dtransform;
70		}
71	LEN	else if( dynamic_cast< Rigid2DTransformType *>(this->m_Transform.GetPointer() ))
72		{
73		transformType = Rigid2Dtransfrom;
74		}
75
76
77		unsigned int numberOfLandmarks = m_FixedLandmarks.size();
78
79
80		// If images come from filters, then update those filters.
81		switch( transformType )
82		{
83		case VersorRigid3Dtransform:
84		{
85		// Sanity check
86		if( FixedImageType::ImageDimension != 3 )
87		{
88		itkExceptionMacro(
89	LEN	"Transform is VersorRigid3DTransform and Fixed image dimension is not 3");
90		return;
91		}
92		if( MovingImageType::ImageDimension != 3 )
93		{
94		itkExceptionMacro(
95	LEN	"Transform is VersorRigid3DTransform and Moving image dimension is not 3");
96		return;
97		}
98
99	LEN	// --- compute the necessary transform to match the two sets of landmarks ---
100		//
101		//
102		//    The solution is based on
103		//    Berthold K. P. Horn (1987),
104	LEN	//    "Closed-form solution of absolute orientation using unit quaternions,"
105		//    Journal of the Optical Society of America A, 4:629-642
106		//
107		//
108		//    Original python implementation by David G. Gobbi
109		//    Readpted from the code in VTK: Hybrid/vtkLandmarkTransform
110		//
111	LEN	//----------------------------------------------------------------------------
112
113	LEN	VersorRigid3DTransformType *transform = dynamic_cast< VersorRigid3DTransformType *>(
114	LEN	this->m_Transform.GetPointer() );
115
116		typedef typename VersorRigid3DTransformType::OutputVectorType VectorType;
117	TDA	typedef typename VersorRigid3DTransformType::OutputPointType PointType;
118	LEN,TDA	typedef typename VersorRigid3DTransformType::CenterType RotationCenterType;
119
120		RotationCenterType rotationCenter = transform->GetCenter();
121		VectorType fixedVector;
122		fixedVector.Fill( 0.0 );
123
124		// Compute the centroids
125		PointsContainerConstIterator fixedItr = m_FixedLandmarks.begin();
126		while( fixedItr != m_FixedLandmarks.end() )
127		{
128		fixedVector[0] += (*fixedItr)[0] - rotationCenter[0];
129		fixedVector[1] += (*fixedItr)[1] - rotationCenter[1];
130		fixedVector[2] += (*fixedItr)[2] - rotationCenter[2];
131		++fixedItr;
132		}
133
134		VectorType movingVector;
135		movingVector.Fill( 0.0 );
136
137		PointsContainerConstIterator movingItr = m_MovingLandmarks.begin();
138		while( movingItr != m_MovingLandmarks.end() )
139		{
140		movingVector[0] += (*movingItr)[0] - rotationCenter[0];
141		movingVector[1] += (*movingItr)[1] - rotationCenter[1];
142		movingVector[2] += (*movingItr)[2] - rotationCenter[2];
143		++movingItr;
144		}
145
146		VectorType fixedCentroidFromRotationCenter;
147		VectorType movingCentroidFromRotationCenter;
148
149		for(unsigned int ic=0; ic<ImageDimension; ic++)
150		{
151	LEN	fixedCentroidFromRotationCenter[ic]  = fixedVector[ic]  / m_FixedLandmarks.size();
152	LEN	movingCentroidFromRotationCenter[ic] = movingVector[ic] / m_MovingLandmarks.size();
153		}
154
155	LEN	PointType fixedCentroid  = rotationCenter + fixedCentroidFromRotationCenter;
156	LEN	PointType movingCentroid = rotationCenter + movingCentroidFromRotationCenter;
157
158		itkDebugMacro(<< "fixed centroid  = " <<  fixedCentroid);
159		itkDebugMacro(<< "moving centroid  = " << movingCentroid);
160
161		typedef typename VersorRigid3DTransformType::VersorType VersorType;
162
163		VersorType versor;
164
165		// If we have at least 3 landmarks, we can compute a rotation.
166		// Otherwise the versor will be an identity versor.
167		if( numberOfLandmarks >= ImageDimension )
168		{
169		itk::Matrix<double,ImageDimension,ImageDimension> M;
170
171		fixedItr  = m_FixedLandmarks.begin();
172		movingItr = m_MovingLandmarks.begin();
173
174		VectorType fixedCentered;
175		VectorType movingCentered;
176
177		fixedCentered.Fill( 0.0 );
178		movingCentered.Fill( 0.0 );
179
180		int ii=0;
181		// Computations are relative to the Center of Rotation.
182		while( movingItr != m_MovingLandmarks.end() )
183		{
184		for(unsigned int i=0; i<ImageDimension; i++)
185		{
186		fixedCentered[i]  = (*fixedItr)[i]  - fixedCentroid[i];
187		movingCentered[i] = (*movingItr)[i] - movingCentroid[i];
188		}
189
190		for(unsigned int i=0; i<ImageDimension; i++)
191		{
192		for(unsigned int j=0; j<ImageDimension; j++)
193		{
194		// mmm this indices i,j may have to be reverted...
195		M[i][j] += fixedCentered[i] * movingCentered[j];
196		}
197		}
198
199		++ii;
200		itkDebugMacro(<< "f_" << ii << " = " << fixedCentered );
201		itkDebugMacro(<< "m_" << ii << " = " << movingCentered );
202
203		++movingItr;
204		++fixedItr;
205		}
206
207		// -- build the 4x4 matrix N --
208
209		itk::Matrix<double,4,4> N;
210
211		// on-diagonal elements
212		N[0][0] =  M[0][0] +M[1][1] +M[2][2];
213		N[1][1] =  M[0][0] -M[1][1] -M[2][2];
214		N[2][2] = -M[0][0] +M[1][1] -M[2][2];
215		N[3][3] = -M[0][0] -M[1][1] +M[2][2];
216		// off-diagonal elements
217		N[0][1] = N[1][0] = M[1][2] -M[2][1];
218		N[0][2] = N[2][0] = M[2][0] -M[0][2];
219		N[0][3] = N[3][0] = M[0][1] -M[1][0];
220
221		N[1][2] = N[2][1] = M[0][1] +M[1][0];
222		N[1][3] = N[3][1] = M[2][0] +M[0][2];
223		N[2][3] = N[3][2] = M[1][2] +M[2][1];
224
225
226		itkDebugMacro( << "For Closed form solution: ");
227		itkDebugMacro(<< "M matrix " << M );
228		itkDebugMacro(<< "N matrix " << N );
229
230		vnl_matrix<double> eigenVectors(4,4);
231		vnl_vector<double> eigenValues(4);
232
233		typedef itk::SymmetricEigenAnalysis<
234		itk::Matrix< double,4,4 >,
235	IND	**************vnl_vector< double >,
236	IND	**************vnl_matrix< double > > SymmetricEigenAnalysisType;
237		SymmetricEigenAnalysisType symmetricEigenSystem(4);
238
239	LEN	symmetricEigenSystem.ComputeEigenValuesAndVectors( N, eigenValues, eigenVectors );
240
241		itkDebugMacro( << "EigenVectors " << eigenVectors);
242		itkDebugMacro( << "EigenValues " << eigenValues);
243
244	LEN	// By default eigen values are sorted in ascending order.  therefore the maximum
245	LEN	// eigen value is the one  in the fourth place = index 3. We need the eigen
246	LEN	// vector associated with the maximum eigenvalue, so we take the eigenvector
247		// from the last row, index=3.
248
249		versor.Set( eigenVectors[3][1],
250		eigenVectors[3][2],
251		eigenVectors[3][3],
252		eigenVectors[3][0]  );
253		itkDebugMacro(<< "Resulting versor" << versor);
254		}
255		else
256		{
257	LEN	itkWarningMacro(<< "Less than 3 landmarks available. Rotation is not computed");
258		}
259
260		transform->SetCenter(fixedCentroid);
261		transform->SetRotation( versor );
262
263		VectorType translation = transform->GetTranslation();
264		translation += movingCentroid - fixedCentroid;
265		transform->SetTranslation( translation );
266
267		break;
268		}
269
270
271		case Rigid2Dtransfrom:
272		{
273		// Sanity check
274		if( FixedImageType::ImageDimension != 2 )
275		{
276		itkExceptionMacro(
277		"Transform is Rigid2DTransfrom and Fixed image dimension is not 2");
278		return;
279		}
280		if( MovingImageType::ImageDimension != 2 )
281		{
282		itkExceptionMacro(
283		"Transform is Rigid2DTransform and Moving image dimension is not 2");
284		return;
285		}
286
287		Rigid2DTransformType *transform = dynamic_cast< Rigid2DTransformType *>(
288		this->m_Transform.GetPointer() );
289
290		typedef typename Rigid2DTransformType::OutputVectorType VectorType;
291	TDA	typedef typename Rigid2DTransformType::OutputPointType PointType;
292
293		//Initialize the transform to identity
294		transform->SetIdentity();
295		PointType rotationCenter = transform->GetCenter();
296
297		VectorType fixedVector;
298		fixedVector.Fill( 0.0 );
299
300		// Compute the centroids
301		PointsContainerConstIterator fixedItr = m_FixedLandmarks.begin();
302		while( fixedItr != m_FixedLandmarks.end() )
303		{
304		fixedVector[0] += (*fixedItr)[0] - rotationCenter[0];
305		fixedVector[1] += (*fixedItr)[1] - rotationCenter[1];
306		++fixedItr;
307		}
308
309		VectorType movingVector;
310		movingVector.Fill( 0.0 );
311
312		PointsContainerConstIterator movingItr = m_MovingLandmarks.begin();
313		while( movingItr != m_MovingLandmarks.end() )
314		{
315		movingVector[0] += (*movingItr)[0] - rotationCenter[0];
316		movingVector[1] += (*movingItr)[1] - rotationCenter[1];
317		++movingItr;
318		}
319
320		VectorType fixedCentroidFromRotationCenter;
321		VectorType movingCentroidFromRotationCenter;
322
323		for(unsigned int ic=0; ic<ImageDimension; ic++)
324		{
325	LEN	fixedCentroidFromRotationCenter[ic]  = fixedVector[ic]  / m_FixedLandmarks.size();
326	LEN	movingCentroidFromRotationCenter[ic] = movingVector[ic] / m_MovingLandmarks.size();
327		}
328
329	LEN	PointType fixedCentroid  = rotationCenter + fixedCentroidFromRotationCenter;
330	LEN	PointType movingCentroid = rotationCenter + movingCentroidFromRotationCenter;
331
332		itkDebugMacro(<< "fixed centroid  = " <<  fixedCentroid);
333		itkDebugMacro(<< "moving centroid  = " << movingCentroid);
334
335		double rotationAngle = 0.0;
336
337		// If we have at least 2 landmarks, we can compute a rotation.
338		// Otherwise the rotation matrix will be identity.
339		//
340		// For the Rigid2DTransform, the least squares error will be minimized
341		// by choosing the offset as the distance between the two centroids,
342		// fixed centroid (after having undergone the rotation transform, that
343		// we must compute) and the moving centroid.
344		// The rotation angle will be given by the cross and dot products of the
345		// fixed and moving landmark vectors, the vectors being computed relative
346		// to the fixed and moving centroids.
347		if( numberOfLandmarks >= 2 )
348		{
349		fixedItr  = m_FixedLandmarks.begin();
350		movingItr = m_MovingLandmarks.begin();
351
352		VectorType fixedCentered;
353		VectorType movingCentered;
354
355		fixedCentered.Fill( 0.0 );
356		movingCentered.Fill( 0.0 );
357
358		int ii=0;
359
360		double s_dot   = 0;
361		double s_cross = 0;
362		// Computations are relative to the Center of Rotation.
363		while( movingItr != m_MovingLandmarks.end() )
364		{
365		fixedCentered[0]  = (*fixedItr)[0]  - fixedCentroid[0];
366		movingCentered[0] = (*movingItr)[0] - movingCentroid[0];
367		fixedCentered[1]  = (*fixedItr)[1]  - fixedCentroid[1];
368		movingCentered[1] = (*movingItr)[1] - movingCentroid[1];
369
370		s_dot += (movingCentered[0] * fixedCentered[0]) +
371	IND	*******************(movingCentered[1] * fixedCentered[1]);
372		s_cross += (movingCentered[1] * fixedCentered[0]) -
373	IND	*********************(movingCentered[0] * fixedCentered[1]);
374
375		++ii;
376		itkDebugMacro(<< "f_" << ii << " = " << fixedCentered );
377		itkDebugMacro(<< "m_" << ii << " = " << movingCentered );
378
379		++movingItr;
380		++fixedItr;
381		}
382
383	LEN	itkDebugMacro(<< "Dot Product of landmarks: " << s_dot << " Cross Product: " << s_cross);
384		if( fabs(s_dot) > 0.00005 )
385		{
386		rotationAngle = atan2(s_cross, s_dot);
387		}
388		}
389		else
390		{
391	LEN	itkWarningMacro(<< "Less than 2 landmarks available. Rotation is not computed");
392		}
393
394		typename Rigid2DTransformType::Pointer t = Rigid2DTransformType::New();
395		t->SetIdentity();
396		t->SetAngle( rotationAngle );
397
398		transform->SetCenter( fixedCentroid );
399		transform->SetAngle( rotationAngle );
400
401		VectorType translation = transform->GetTranslation();
402		itkDebugMacro(<< "Initial transform translation: " << translation);
403		translation += movingCentroid - fixedCentroid;
404	LEN	itkDebugMacro(<< "translation computed as difference of centroids: " << translation);
405		transform->SetTranslation( translation );
406
407		break;
408		}
409
410
411
412		case Else:
413	LEN,IND	******itkWarningMacro(<< "Landmark initialization using the specified input transform not implemented");
414		m_Transform->SetIdentity();
415
416
417		default:
418	LEN,IND	******itkWarningMacro(<< "Landmark initialization using the specified input transform not implemented");
419		m_Transform->SetIdentity();
420
421		}
422
423
424		}
425
426
427
428
429		template < class TTransform, class TFixedImage, class TMovingImage >
430		void
431		LandmarkBasedTransformInitializer<TTransform, TFixedImage, TMovingImage >
432		::PrintSelf(std::ostream& os, Indent indent) const
433		{
434		Superclass::PrintSelf(os,indent);
435
436		os << indent << "Transform   = " << std::endl;
437		if (m_Transform)
438		{
439		os << indent << m_Transform  << std::endl;
440		}
441		else
442		{
443		os << indent << "None" << std::endl;
444		}
445
446		os << indent << "FixedImage   = " << std::endl;
447		if (m_FixedImage)
448		{
449		os << indent << m_FixedImage  << std::endl;
450		}
451		else
452		{
453		os << indent << "None" << std::endl;
454		}
455
456		os << indent << "MovingImage   = " << std::endl;
457		if (m_MovingImage)
458		{
459		os << indent << m_MovingImage  << std::endl;
460		}
461		else
462		{
463		os << indent << "None" << std::endl;
464		}
465
466		os << indent << "Fixed Landmarks: " << std::endl;
467		PointsContainerConstIterator fitr = m_FixedLandmarks.begin();
468		while( fitr != m_FixedLandmarks.end() )
469		{
470		os << indent << *fitr << std::endl;
471		++fitr;
472		}
473		os << indent << "Moving Landmarks: " << std::endl;
474		PointsContainerConstIterator mitr = m_MovingLandmarks.begin();
475		while( mitr != m_MovingLandmarks.end() )
476		{
477		os << indent << *mitr << std::endl;
478		++mitr;
479		}
480
481		}
482
483		}  // namespace itk
484
485		#endif /* __itkLandmarkBasedTransformInitializer_txx */
486

---

Generated by KWStyle 1.0b on Tuesday January,17 at 02:14:18PM

© Kitware Inc.