KWStyle - itkKernelTransform.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkKernelTransform.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	DEF	#ifndef _itkKernelTransform_txx
18	DEF	#define _itkKernelTransform_txx
19		#include "itkKernelTransform.h"
20
21		namespace itk
22		{
23
24
25		/**
26		*
27		*/
28		template <class TScalarType, unsigned int NDimensions>
29		KernelTransform<TScalarType, NDimensions>::
30		KernelTransform():Superclass(
31		NDimensions,
32		NDimensions )
33		// the second NDimensions is associated is provided as
34		// a tentative number for initializing the Jacobian.
35		// The matrix can be resized at run time so this number
36		// here is irrelevant. The correct size of the Jacobian
37		// will be NDimension X NDimension.NumberOfLandMarks.
38		{
39
40		m_I.set_identity();
41		m_SourceLandmarks = PointSetType::New();
42		m_TargetLandmarks = PointSetType::New();
43		m_Displacements   = VectorSetType::New();
44		m_WMatrixComputed = false;
45
46		m_Stiffness = 0.0;
47		}
48
49		/**
50		*
51		*/
52		template <class TScalarType, unsigned int NDimensions>
53		KernelTransform<TScalarType, NDimensions>::
54		~KernelTransform()
55		{
56		}
57
58
59		/**
60		*
61		*/
62		template <class TScalarType, unsigned int NDimensions>
63		void
64		KernelTransform<TScalarType, NDimensions>::
65		SetSourceLandmarks(PointSetType * landmarks)
66		{
67		itkDebugMacro("setting SourceLandmarks to " << landmarks );
68		if (this->m_SourceLandmarks != landmarks)
69		{
70		this->m_SourceLandmarks = landmarks;
71		this->UpdateParameters();
72		this->Modified();
73		}
74		}
75
76
77		/**
78		*
79		*/
80		template <class TScalarType, unsigned int NDimensions>
81		void
82		KernelTransform<TScalarType, NDimensions>::
83		SetTargetLandmarks(PointSetType * landmarks)
84		{
85		itkDebugMacro("setting TargetLandmarks to " << landmarks );
86		if (this->m_TargetLandmarks != landmarks)
87		{
88		this->m_TargetLandmarks = landmarks;
89		this->UpdateParameters();
90		this->Modified();
91		}
92		}
93
94
95	EML
96		/**
97		*
98		*/
99		template <class TScalarType, unsigned int NDimensions>
100		const typename KernelTransform<TScalarType, NDimensions>::GMatrixType &
101		KernelTransform<TScalarType, NDimensions>::
102		ComputeG( const InputVectorType & ) const
103		{
104		//
105		// Should an Exception be thrown here  ?
106		//
107		itkWarningMacro(<< "ComputeG() should be reimplemented in the subclass !!");
108		return m_GMatrix;
109		}
110
111		/**
112		*
113		*/
114		template <class TScalarType, unsigned int NDimensions>
115		const typename KernelTransform<TScalarType, NDimensions>::GMatrixType &
116		KernelTransform<TScalarType, NDimensions>::
117		ComputeReflexiveG( PointsIterator ) const
118		{
119		m_GMatrix.fill( NumericTraits< TScalarType >::Zero );
120		m_GMatrix.fill_diagonal( m_Stiffness );
121
122		return m_GMatrix;
123		}
124
125
126	EML
127	EML
128		/**
129	WCM	* Default implementation of the the method. This can be overloaded
130		* in transforms whose kernel produce diagonal G matrices.
131		*/
132		template <class TScalarType, unsigned int NDimensions>
133		void
134		KernelTransform<TScalarType, NDimensions>::
135		ComputeDeformationContribution( const InputPointType  & thisPoint,
136		OutputPointType & result     ) const
137		{
138
139		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
140
141		PointsIterator sp  = m_SourceLandmarks->GetPoints()->Begin();
142
143		for(unsigned int lnd=0; lnd < numberOfLandmarks; lnd++ )
144		{
145		const GMatrixType & Gmatrix = ComputeG( thisPoint - sp->Value() );
146		for(unsigned int dim=0; dim < NDimensions; dim++ )
147		{
148		for(unsigned int odim=0; odim < NDimensions; odim++ )
149		{
150		result[ odim ] += Gmatrix(dim, odim ) * m_DMatrix(dim,lnd);
151		}
152		}
153		++sp;
154		}
155
156		}
157
158
159		/**
160		*
161		*/
162		template <class TScalarType, unsigned int NDimensions>
163		void KernelTransform<TScalarType, NDimensions>
164		::ComputeD(void)
165		{
166		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
167
168		PointsIterator sp  = m_SourceLandmarks->GetPoints()->Begin();
169		PointsIterator tp  = m_TargetLandmarks->GetPoints()->Begin();
170		PointsIterator end = m_SourceLandmarks->GetPoints()->End();
171
172		m_Displacements->Reserve( numberOfLandmarks );
173		typename VectorSetType::Iterator vt = m_Displacements->Begin();
174
175		while( sp != end )
176		{
177		vt->Value() = tp->Value() - sp->Value();
178		vt++;
179		sp++;
180		tp++;
181		}
182		}
183
184		/**
185		*
186		*/
187		template <class TScalarType, unsigned int NDimensions>
188		void KernelTransform<TScalarType, NDimensions>
189		::ComputeWMatrix(void)
190		{
191
192		typedef vnl_svd<TScalarType>  SVDSolverType;
193
194		this->ComputeL();
195		this->ComputeY();
196		SVDSolverType svd( m_LMatrix, 1e-8 );
197		m_WMatrix = svd.solve( m_YMatrix );
198
199		this->ReorganizeW();
200
201		}
202
203		/**
204		*
205		*/
206		template <class TScalarType, unsigned int NDimensions>
207		void KernelTransform<TScalarType, NDimensions>::
208		ComputeL(void)
209		{
210		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
211		vnl_matrix<TScalarType> O2(NDimensions*(NDimensions+1),
212	IND	*****************************NDimensions*(NDimensions+1), 0);
213
214		this->ComputeP();
215		this->ComputeK();
216
217		m_LMatrix.set_size( NDimensions*(numberOfLandmarks+NDimensions+1),
218	IND	********************NDimensions*(numberOfLandmarks+NDimensions+1) );
219		m_LMatrix.fill( 0.0 );
220
221		m_LMatrix.update( m_KMatrix, 0, 0 );
222		m_LMatrix.update( m_PMatrix, 0, m_KMatrix.columns() );
223		m_LMatrix.update( m_PMatrix.transpose(), m_KMatrix.rows(), 0);
224		m_LMatrix.update( O2, m_KMatrix.rows(), m_KMatrix.columns());
225
226		}
227
228
229		/**
230		*
231		*/
232		template <class TScalarType, unsigned int NDimensions>
233		void KernelTransform<TScalarType, NDimensions>::
234		ComputeK(void)
235		{
236		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
237		GMatrixType G;
238
239		this->ComputeD();
240
241		m_KMatrix.set_size( NDimensions * numberOfLandmarks,
242		NDimensions * numberOfLandmarks );
243
244		m_KMatrix.fill( 0.0 );
245
246		PointsIterator p1  = m_SourceLandmarks->GetPoints()->Begin();
247		PointsIterator end = m_SourceLandmarks->GetPoints()->End();
248
249		// K matrix is symmetric, so only evaluate the upper triangle and
250		// store the values in bot the upper and lower triangle
251		unsigned int i = 0;
252		while( p1 != end )
253		{
254		PointsIterator p2 = p1; // start at the diagonal element
255		unsigned int j = i;
256
257		// Compute the block diagonal element, i.e. kernel for pi->pi
258		G = ComputeReflexiveG(p1);
259		m_KMatrix.update(G, i*NDimensions, i*NDimensions);
260		p2++;
261		j++;
262
263		// Compute the upper (and copy into lower) triangular part of K
264		while( p2 != end )
265		{
266		const InputVectorType s = p1.Value() - p2.Value();
267		G = ComputeG(s);
268		// write value in upper and lower triangle of matrix
269		m_KMatrix.update(G, i*NDimensions, j*NDimensions);
270		m_KMatrix.update(G, j*NDimensions, i*NDimensions);
271		p2++;
272		j++;
273		}
274		p1++;
275		i++;
276		}
277		}
278
279
280	EML
281		/**
282		*
283		*/
284		template <class TScalarType, unsigned int NDimensions>
285		void KernelTransform<TScalarType, NDimensions>::
286		ComputeP()
287		{
288		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
289		IMatrixType I;
290		IMatrixType temp;
291		InputPointType p;
292
293		I.set_identity();
294		m_PMatrix.set_size( NDimensions*numberOfLandmarks,
295		NDimensions*(NDimensions+1) );
296		m_PMatrix.fill( 0.0 );
297		for (unsigned int i = 0; i < numberOfLandmarks; i++)
298		{
299		m_SourceLandmarks->GetPoint(i, &p);
300		for (unsigned int j = 0; j < NDimensions; j++)
301		{
302		temp = I * p[j];
303		m_PMatrix.update(temp, i*NDimensions, j*NDimensions);
304		}
305		m_PMatrix.update(I, i*NDimensions, NDimensions*NDimensions);
306		}
307		}
308
309
310	EML
311		/**
312		*
313		*/
314		template <class TScalarType, unsigned int NDimensions>
315		void KernelTransform<TScalarType, NDimensions>::
316		ComputeY(void)
317		{
318		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
319
320		typename VectorSetType::ConstIterator displacement =
321	IND	****m_Displacements->Begin();
322
323		m_YMatrix.set_size( NDimensions*(numberOfLandmarks+NDimensions+1), 1);
324
325		m_YMatrix.fill( 0.0 );
326
327		for (unsigned int i = 0; i < numberOfLandmarks; i++)
328		{
329		for (unsigned int j = 0; j < NDimensions; j++)
330		{
331		m_YMatrix.put(i*NDimensions+j, 0, displacement.Value()[j]);
332		}
333		displacement++;
334		}
335
336		for (unsigned int i = 0; i < NDimensions*(NDimensions+1); i++)
337		{
338		m_YMatrix.put(numberOfLandmarks*NDimensions+i, 0, 0);
339		}
340		}
341
342
343		/**
344		*
345		*/
346		template <class TScalarType, unsigned int NDimensions>
347		void
348		KernelTransform<TScalarType, NDimensions>
349		::ReorganizeW(void)
350		{
351		unsigned long numberOfLandmarks = m_SourceLandmarks->GetNumberOfPoints();
352
353		// The deformable (non-affine) part of the registration goes here
354		m_DMatrix.set_size(NDimensions,numberOfLandmarks);
355		unsigned int ci = 0;
356		for(unsigned int lnd=0; lnd < numberOfLandmarks; lnd++ )
357		{
358		for(unsigned int dim=0; dim < NDimensions; dim++ )
359		{
360		m_DMatrix(dim,lnd) = m_WMatrix(ci++,0);
361		}
362		}
363
364		// This matrix holds the rotational part of the Affine component
365		for(unsigned int j=0; j < NDimensions; j++ )
366		{
367		for(unsigned int i=0; i < NDimensions; i++ )
368		{
369		m_AMatrix(i,j) = m_WMatrix(ci++,0);
370		}
371		}
372
373		// This vector holds the translational part of the Affine component
374		for(unsigned int k=0; k < NDimensions; k++ )
375		{
376		m_BVector(k) = m_WMatrix(ci++,0);
377		}
378
379		// release WMatrix memory by assigning a small one.
380		m_WMatrix = WMatrixType(1,1);
381
382		}
383
384
385	EML
386		/**
387		*
388		*/
389		template <class TScalarType, unsigned int NDimensions>
390		typename KernelTransform<TScalarType, NDimensions>::OutputPointType
391		KernelTransform<TScalarType, NDimensions>
392		::TransformPoint(const InputPointType& thisPoint) const
393		{
394
395		OutputPointType result;
396
397		typedef typename OutputPointType::ValueType ValueType;
398
399		result.Fill( NumericTraits< ValueType >::Zero );
400
401		this->ComputeDeformationContribution( thisPoint, result );
402
403		// Add the rotational part of the Affine component
404		for(unsigned int j=0; j < NDimensions; j++ )
405		{
406		for(unsigned int i=0; i < NDimensions; i++ )
407		{
408		result[i] += m_AMatrix(i,j) * thisPoint[j];
409		}
410		}
411
412
413
414		// This vector holds the translational part of the Affine component
415		for(unsigned int k=0; k < NDimensions; k++ )
416		{
417		result[k] += m_BVector(k) + thisPoint[k];
418		}
419
420		return result;
421
422		}
423
424
425	EML
426	EML
427		// Compute the Jacobian in one position
428		template <class TScalarType, unsigned int NDimensions>
429		const typename KernelTransform<TScalarType,NDimensions>::JacobianType &
430		KernelTransform< TScalarType,NDimensions>::
431		GetJacobian( const InputPointType & ) const
432		{
433
434
435		this->m_Jacobian.Fill( 0.0 );
436
437		// TODO
438		// The Jacobian should be computable in terms of the matrices
439		// used to Transform points...
440
441		return this->m_Jacobian;
442
443		}
444
445
446		// Set the parameters
447		// NOTE that in this transformation both the Source and Target
448		// landmarks could be considered as parameters. It is assumed
449		// here that the Target landmarks are provided by the user and
450		// are not changed during the optimization process required for
451		// registration.
452		template <class TScalarType, unsigned int NDimensions>
453		void
454		KernelTransform<TScalarType, NDimensions>::
455		SetParameters( const ParametersType & parameters )
456		{
457		typename PointsContainer::Pointer landmarks = PointsContainer::New();
458		const unsigned int numberOfLandmarks =  parameters.Size() / NDimensions;
459		landmarks->Reserve( numberOfLandmarks );
460
461		PointsIterator itr = landmarks->Begin();
462		PointsIterator end = landmarks->End();
463
464		InputPointType  landMark;
465
466		unsigned int pcounter = 0;
467		while( itr != end )
468		{
469		for(unsigned int dim=0; dim<NDimensions; dim++)
470		{
471		landMark[ dim ] = parameters[ pcounter ];
472		pcounter++;
473		}
474		itr.Value() = landMark;
475		itr++;
476		}
477
478		m_SourceLandmarks->SetPoints( landmarks );
479		}
480
481		// Set the fixed parameters
482		// Since the API of the SetParameters() function sets the
483		// source landmarks, this function was added to support the
484		// setting of the target landmarks, and allowing the Transform
485		// I/O mechanism to be supported.
486		template <class TScalarType, unsigned int NDimensions>
487		void
488		KernelTransform<TScalarType, NDimensions>::
489		SetFixedParameters( const ParametersType & parameters )
490		{
491		typename PointsContainer::Pointer landmarks = PointsContainer::New();
492		const unsigned int numberOfLandmarks =  parameters.Size() / NDimensions;
493
494		landmarks->Reserve( numberOfLandmarks );
495
496		PointsIterator itr = landmarks->Begin();
497		PointsIterator end = landmarks->End();
498
499		InputPointType  landMark;
500
501		unsigned int pcounter = 0;
502		while( itr != end )
503		{
504		for(unsigned int dim=0; dim<NDimensions; dim++)
505		{
506		landMark[ dim ] = parameters[ pcounter ];
507		pcounter++;
508		}
509		itr.Value() = landMark;
510		itr++;
511		}
512
513		m_TargetLandmarks->SetPoints( landmarks );
514		}
515
516
517		// Update parameters array
518		// They are the components of all the landmarks in the source space
519		template <class TScalarType, unsigned int NDimensions>
520		void
521		KernelTransform<TScalarType, NDimensions>::
522		UpdateParameters( void ) const
523		{
524	LEN	this->m_Parameters = ParametersType( m_SourceLandmarks->GetNumberOfPoints() * NDimensions );
525
526		PointsIterator itr = m_SourceLandmarks->GetPoints()->Begin();
527		PointsIterator end = m_SourceLandmarks->GetPoints()->End();
528
529		unsigned int pcounter = 0;
530		while( itr != end )
531		{
532		InputPointType  landmark = itr.Value();
533		for(unsigned int dim=0; dim<NDimensions; dim++)
534		{
535		this->m_Parameters[ pcounter ] = landmark[ dim ];
536		pcounter++;
537		}
538		itr++;
539		}
540		}
541
542
543	EML
544	EML
545		// Get the parameters
546		// They are the components of all the landmarks in the source space
547		template <class TScalarType, unsigned int NDimensions>
548		const typename KernelTransform<TScalarType, NDimensions>::ParametersType &
549		KernelTransform<TScalarType, NDimensions>::
550		GetParameters( void ) const
551		{
552		this->UpdateParameters();
553		return this->m_Parameters;
554
555		}
556
557
558		// Get the fixed parameters
559		// This returns the target landmark locations
560		// This was added to support the Transform Reader/Writer mechanism
561		template <class TScalarType, unsigned int NDimensions>
562		const typename KernelTransform<TScalarType, NDimensions>::ParametersType &
563		KernelTransform<TScalarType, NDimensions>::
564		GetFixedParameters( void ) const
565		{
566	LEN	this->m_FixedParameters = ParametersType( m_TargetLandmarks->GetNumberOfPoints() * NDimensions );
567
568		PointsIterator itr = m_TargetLandmarks->GetPoints()->Begin();
569		PointsIterator end = m_TargetLandmarks->GetPoints()->End();
570
571		unsigned int pcounter = 0;
572		while( itr != end )
573		{
574		InputPointType  landmark = itr.Value();
575		for(unsigned int dim=0; dim<NDimensions; dim++)
576		{
577		this->m_FixedParameters[ pcounter ] = landmark[ dim ];
578		pcounter++;
579		}
580		itr++;
581		}
582
583	IND	*return this->m_FixedParameters;
584
585		}
586
587
588	EML
589		template <class TScalarType, unsigned int NDimensions>
590		void
591		KernelTransform<TScalarType, NDimensions>::
592		PrintSelf(std::ostream& os, Indent indent) const
593		{
594		Superclass::PrintSelf(os,indent);
595		if (m_SourceLandmarks)
596		{
597		os << indent << "SourceLandmarks: " << std::endl;
598		m_SourceLandmarks->Print(os,indent.GetNextIndent());
599		}
600		if (m_TargetLandmarks)
601		{
602		os << indent << "TargetLandmarks: " << std::endl;
603		m_TargetLandmarks->Print(os,indent.GetNextIndent());
604		}
605		if (m_Displacements)
606		{
607		os << indent << "Displacements: " << std::endl;
608		m_Displacements->Print(os,indent.GetNextIndent());
609		}
610		os << indent << "Stiffness: " << m_Stiffness << std::endl;
611		}
612		} // namespace itk
613
614		#endif
615

---

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

© Kitware Inc.