KWStyle - itkBSplineDeformableTransform.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkBSplineDeformableTransform.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:33 $
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 _itkBSplineDeformableTransform_txx
18	DEF	#define _itkBSplineDeformableTransform_txx
19
20		#include "itkBSplineDeformableTransform.h"
21		#include "itkContinuousIndex.h"
22		#include "itkImageRegionIterator.h"
23		#include "itkImageRegionConstIterator.h"
24		#include "itkImageRegionConstIteratorWithIndex.h"
25		#include "itkIdentityTransform.h"
26
27		namespace itk
28		{
29
30		// Constructor with default arguments
31		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
32		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
33		::BSplineDeformableTransform():Superclass(SpaceDimension,0)
34		{
35
36		// Instantiate a weights function
37		m_WeightsFunction = WeightsFunctionType::New();
38		m_SupportSize = m_WeightsFunction->GetSupportSize();
39
40		// Instantiate an identity transform
41		typedef IdentityTransform<ScalarType,SpaceDimension> IdentityTransformType;
42		typename IdentityTransformType::Pointer id = IdentityTransformType::New();
43		m_BulkTransform = id;
44
45		// Default grid size is zero
46		typename RegionType::SizeType size;
47		typename RegionType::IndexType index;
48		size.Fill( 0 );
49		index.Fill( 0 );
50		m_GridRegion.SetSize( size );
51		m_GridRegion.SetIndex( index );
52
53		m_GridOrigin.Fill( 0.0 );  // default origin is all zeros
54		m_GridSpacing.Fill( 1.0 ); // default spacing is all ones
55
56		m_InputParametersPointer = NULL;
57		m_InternalParametersBuffer = ParametersType(0);
58
59		// Initialize coeffient images
60		for ( unsigned int j = 0; j < SpaceDimension; j++ )
61		{
62		m_WrappedImage[j] = ImageType::New();
63		m_WrappedImage[j]->SetRegions( m_GridRegion );
64		m_WrappedImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
65		m_WrappedImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
66		m_CoefficientImage[j] = NULL;
67		}
68
69		// Setup variables for computing interpolation
70		m_Offset = SplineOrder / 2;
71		if ( SplineOrder % 2 )
72		{
73		m_SplineOrderOdd = true;
74		}
75		else
76		{
77		m_SplineOrderOdd = false;
78		}
79		m_ValidRegion = m_GridRegion;
80
81		// Initialize jacobian images
82		for ( unsigned int j = 0; j < SpaceDimension; j++ )
83		{
84		m_JacobianImage[j] = ImageType::New();
85		m_JacobianImage[j]->SetRegions( m_GridRegion );
86		m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
87		m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
88		}
89
90		/**********************************************************
91	IND	****Fixed Parameters store the following information:
92	IND	********Grid Size
93	IND	********Grid Origin
94	IND	********Grid Spacing
95	IND	*****The size of these is equal to the  NInputDimensions
96	IND	************************************************************/
97		this->m_FixedParameters.SetSize ( NDimensions * 3 );
98		this->m_FixedParameters.Fill ( 0.0 );
99
100		m_LastJacobianIndex = m_ValidRegion.GetIndex();
101
102		}
103
104
105		// Destructor
106		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
107		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
108		::~BSplineDeformableTransform()
109		{
110
111		}
112
113
114		// Get the number of parameters
115		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
116		unsigned int
117		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
118		::GetNumberOfParameters(void) const
119		{
120
121		// The number of parameters equal SpaceDimension * number of
122		// of pixels in the grid region.
123		return ( static_cast<unsigned int>( SpaceDimension ) *
124	IND	***********static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
125
126		}
127
128
129		// Get the number of parameters per dimension
130		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
131		unsigned int
132		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
133		::GetNumberOfParametersPerDimension(void) const
134		{
135		// The number of parameters per dimension equal number of
136		// of pixels in the grid region.
137		return ( static_cast<unsigned int>( m_GridRegion.GetNumberOfPixels() ) );
138
139		}
140
141
142		// Set the grid region
143		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
144		void
145		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
146		::SetGridRegion( const RegionType& region )
147		{
148		if ( m_GridRegion != region )
149		{
150
151		m_GridRegion = region;
152
153		// set regions for each coefficient and jacobian image
154		for ( unsigned int j = 0; j < SpaceDimension; j++ )
155		{
156		m_WrappedImage[j]->SetRegions( m_GridRegion );
157		m_JacobianImage[j]->SetRegions( m_GridRegion );
158		}
159
160		// Set the valid region
161		// If the grid spans the interval [start,last].
162	LEN	// The valid interval for evaluation is [start+offset,last-offset] when spline order is even.
163	LEN	// The valid interval for evaluation is [start+offset,last-offset) when spline order is odd
164		// Where offset = floor( spline / 2 ).
165	LEN	// Note that the last pixel is not included in the valid region with odd spline orders.
166		typename RegionType::SizeType size = m_GridRegion.GetSize();
167		typename RegionType::IndexType index = m_GridRegion.GetIndex();
168		for ( unsigned int j = 0; j < SpaceDimension; j++ )
169		{
170		index[j] +=
171		static_cast< typename RegionType::IndexValueType >( m_Offset );
172		size[j] -=
173	IND	********static_cast< typename RegionType::SizeValueType> ( 2 * m_Offset );
174		m_ValidRegionLast[j] = index[j] +
175	IND	********static_cast< typename RegionType::IndexValueType >( size[j] ) - 1;
176		}
177		m_ValidRegion.SetSize( size );
178		m_ValidRegion.SetIndex( index );
179
180		this->Modified();
181		}
182
183		}
184
185
186		// Set the grid spacing
187		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
188		void
189		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
190		::SetGridSpacing( const SpacingType& spacing )
191		{
192		if ( m_GridSpacing != spacing )
193		{
194		m_GridSpacing = spacing;
195
196		// set spacing for each coefficient and jacobian image
197		for ( unsigned int j = 0; j < SpaceDimension; j++ )
198		{
199		m_WrappedImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
200		m_JacobianImage[j]->SetSpacing( m_GridSpacing.GetDataPointer() );
201		}
202
203		this->Modified();
204		}
205
206		}
207
208
209		// Set the grid origin
210		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
211		void
212		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
213		::SetGridOrigin( const OriginType& origin )
214		{
215		if ( m_GridOrigin != origin )
216		{
217		m_GridOrigin = origin;
218
219		// set spacing for each coefficient and jacobianimage
220		for ( unsigned int j = 0; j < SpaceDimension; j++ )
221		{
222		m_WrappedImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
223		m_JacobianImage[j]->SetOrigin( m_GridOrigin.GetDataPointer() );
224		}
225
226		this->Modified();
227		}
228
229		}
230
231
232		// Set the parameters
233		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
234		void
235		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
236		::SetIdentity()
237		{
238	IND	***if( m_InputParametersPointer )
239	IND	*****{
240	LEN,IND	*****ParametersType * parameters = const_cast<ParametersType *>( m_InputParametersPointer );
241	IND	*****parameters->Fill( 0.0 );
242	IND	*****this->Modified();
243	IND	*****}
244		}
245
246
247		// Set the parameters
248		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
249		void
250		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
251		::SetParameters( const ParametersType & parameters )
252		{
253
254		// check if the number of parameters match the
255		// expected number of parameters
256		if ( parameters.Size() != this->GetNumberOfParameters() )
257		{
258	LEN	itkExceptionMacro(<<"Mismatched between parameters size " << parameters.size()
259	LEN	<< " and region size " << m_GridRegion.GetNumberOfPixels() );
260		}
261
262		// Clean up buffered parameters
263		m_InternalParametersBuffer = ParametersType( 0 );
264
265		// Keep a reference to the input parameters
266		m_InputParametersPointer = ¶meters;
267
268		// Wrap flat array as images of coefficients
269		this->WrapAsImages();
270
271		// Modified is always called since we just have a pointer to the
272		// parameters and cannot know if the parameters have changed.
273		this->Modified();
274		}
275
276		// Set the Fixed Parameters
277		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
278		void
279		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
280		::SetFixedParameters( const ParametersType & parameters )
281		{
282
283		// check if the number of parameters match the
284		// expected number of parameters
285		if ( parameters.Size() != NDimensions*3 )
286		{
287	LEN	itkExceptionMacro(<<"Mismatched between parameters size " << parameters.size()
288		<< " and number of fixed parameters " << NDimensions*3 );
289		}
290
291		/**********************************************************
292		Fixed Parameters store the following information:
293		Grid Size
294		Grid Origin
295		Grid Spacing
296		The size of these is equal to the  NInputDimensions
297		**********************************************************/
298
299		/*** Set the Grid Parameters ***/
300		SizeType   gridSize;
301		for (unsigned int i=0;i<NDimensions;i++)
302		{
303		gridSize[i] = static_cast<int> (parameters[i]);
304		}
305		RegionType bsplineRegion;
306		bsplineRegion.SetSize( gridSize );
307
308		/*** Set the Origin Parameters ***/
309		OriginType origin;
310		for (unsigned int i=0;i<NDimensions;i++)
311		{
312		origin[i] = parameters[NDimensions+i];
313		}
314
315		/*** Set the Spacing Parameters ***/
316		SpacingType spacing;
317		for (unsigned int i=0;i<NDimensions;i++)
318		{
319		spacing[i] = parameters[2*NDimensions+i];
320		}
321
322
323		this->SetGridSpacing( spacing );
324		this->SetGridOrigin( origin );
325		this->SetGridRegion( bsplineRegion );
326
327		this->Modified();
328		}
329
330
331		// Wrap flat parameters as images
332		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
333		void
334		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
335		::WrapAsImages()
336		{
337
338		/**
339		* Wrap flat parameters array into SpaceDimension number of ITK images
340		* NOTE: For efficiency, parameters are not copied locally. The parameters
341		* are assumed to be maintained by the caller.
342		*/
343	LEN	PixelType * dataPointer = const_cast<PixelType *>(( m_InputParametersPointer->data_block() ));
344		unsigned int numberOfPixels = m_GridRegion.GetNumberOfPixels();
345
346		for ( unsigned int j = 0; j < SpaceDimension; j++ )
347		{
348		m_WrappedImage[j]->GetPixelContainer()->
349	IND	******SetImportPointer( dataPointer, numberOfPixels );
350		dataPointer += numberOfPixels;
351		m_CoefficientImage[j] = m_WrappedImage[j];
352		}
353
354		/**
355		* Allocate memory for Jacobian and wrap into SpaceDimension number
356		* of ITK images
357		*/
358		this->m_Jacobian.set_size( SpaceDimension, this->GetNumberOfParameters() );
359		this->m_Jacobian.Fill( NumericTraits<JacobianPixelType>::Zero );
360		m_LastJacobianIndex = m_ValidRegion.GetIndex();
361		JacobianPixelType * jacobianDataPointer = this->m_Jacobian.data_block();
362
363		for ( unsigned int j = 0; j < SpaceDimension; j++ )
364		{
365		m_JacobianImage[j]->GetPixelContainer()->
366	IND	******SetImportPointer( jacobianDataPointer, numberOfPixels );
367		jacobianDataPointer += this->GetNumberOfParameters() + numberOfPixels;
368		}
369		}
370
371
372		// Set the parameters by value
373		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
374		void
375		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
376		::SetParametersByValue( const ParametersType & parameters )
377		{
378
379		// check if the number of parameters match the
380		// expected number of parameters
381		if ( parameters.Size() != this->GetNumberOfParameters() )
382		{
383	LEN	itkExceptionMacro(<<"Mismatched between parameters size " << parameters.size()
384	LEN	<< " and region size " << m_GridRegion.GetNumberOfPixels() );
385		}
386
387		// copy it
388		m_InternalParametersBuffer = parameters;
389		m_InputParametersPointer = &m_InternalParametersBuffer;
390
391		// wrap flat array as images of coefficients
392		this->WrapAsImages();
393
394		// Modified is always called since we just have a pointer to the
395		// parameters and cannot know if the parameters have changed.
396		this->Modified();
397
398		}
399
400		// Get the parameters
401		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
402		const
403		typename BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
404		::ParametersType &
405		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
406		::GetParameters( void ) const
407		{
408		/** NOTE: For efficiency, this class does not keep a copy of the parameters -
409		* it just keeps pointer to input parameters.
410		*/
411		return (*m_InputParametersPointer);
412		}
413
414
415		// Get the parameters
416		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
417		const
418		typename BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
419		::ParametersType &
420		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
421		::GetFixedParameters( void ) const
422		{
423		RegionType resRegion = this->GetGridRegion(  );
424
425		for (unsigned int i=0;i<NDimensions;i++)
426		{
427		this->m_FixedParameters[i] = (resRegion.GetSize())[i];
428		}
429		for (unsigned int i=0;i<NDimensions;i++)
430		{
431		this->m_FixedParameters[NDimensions+i] = (this->GetGridOrigin())[i];
432		}
433		for (unsigned int i=0;i<NDimensions;i++)
434		{
435		this->m_FixedParameters[2*NDimensions+i] =  (this->GetGridSpacing())[i];
436		}
437
438		return (this->m_FixedParameters);
439		}
440
441
442
443		// Set the B-Spline coefficients using input images
444		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
445		void
446		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
447		::SetCoefficientImage( ImagePointer images[] )
448		{
449		if ( images[0] )
450		{
451		this->SetGridRegion( images[0]->GetBufferedRegion() );
452		this->SetGridSpacing( images[0]->GetSpacing() );
453		this->SetGridOrigin( images[0]->GetOrigin() );
454
455		for( unsigned int j = 0; j < SpaceDimension; j++ )
456		{
457		m_CoefficientImage[j] = images[j];
458		}
459
460		// Clean up buffered parameters
461		m_InternalParametersBuffer = ParametersType( 0 );
462		m_InputParametersPointer  = NULL;
463
464		}
465
466		}
467
468		// Print self
469		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
470		void
471		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
472		::PrintSelf(std::ostream &os, Indent indent) const
473		{
474
475		unsigned int j;
476
477		this->Superclass::PrintSelf(os,indent);
478
479		os << indent << "GridRegion: " << m_GridRegion << std::endl;
480		os << indent << "GridOrigin: " << m_GridOrigin << std::endl;
481		os << indent << "GridSpacing: " << m_GridSpacing << std::endl;
482
483		os << indent << "CoefficientImage: [ ";
484		for ( j = 0; j < SpaceDimension - 1; j++ )
485		{
486		os << m_CoefficientImage[j].GetPointer() << ", ";
487		}
488		os << m_CoefficientImage[j].GetPointer() << " ]" << std::endl;
489
490		os << indent << "WrappedImage: [ ";
491		for ( j = 0; j < SpaceDimension - 1; j++ )
492		{
493		os << m_WrappedImage[j].GetPointer() << ", ";
494		}
495		os << m_WrappedImage[j].GetPointer() << " ]" << std::endl;
496
497		os << indent << "InputParametersPointer: "
498		<< m_InputParametersPointer << std::endl;
499		os << indent << "ValidRegion: " << m_ValidRegion << std::endl;
500		os << indent << "LastJacobianIndex: " << m_LastJacobianIndex << std::endl;
501		os << indent << "BulkTransform: ";
502		os << m_BulkTransform.GetPointer() << std::endl;
503		os << indent << "WeightsFunction: ";
504		os << m_WeightsFunction.GetPointer() << std::endl;
505
506		if ( m_BulkTransform )
507		{
508		os << indent << "BulkTransformType: "
509		<< m_BulkTransform->GetNameOfClass() << std::endl;
510		}
511
512		}
513
514		// Transform a point
515		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
516		bool
517		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
518		::InsideValidRegion(
519		const ContinuousIndexType& index ) const
520		{
521		bool inside = true;
522
523		if ( !m_ValidRegion.IsInside( index ) )
524		{
525		inside = false;
526		}
527
528		if ( inside && m_SplineOrderOdd )
529		{
530		typedef typename ContinuousIndexType::ValueType ValueType;
531		for( unsigned int j = 0; j < SpaceDimension; j++ )
532		{
533		if ( index[j] >= static_cast<ValueType>( m_ValidRegionLast[j] ) )
534		{
535		inside = false;
536		break;
537		}
538		}
539		}
540
541		return inside;
542		}
543
544		// Transform a point
545		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
546		void
547		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
548		::TransformPoint(
549		const InputPointType & point,
550		OutputPointType & outputPoint,
551		WeightsType & weights,
552		ParameterIndexArrayType & indices,
553		bool& inside ) const
554		{
555
556		unsigned int j;
557		IndexType supportIndex;
558
559		InputPointType transformedPoint;
560		if ( m_BulkTransform )
561		{
562		transformedPoint = m_BulkTransform->TransformPoint( point );
563		}
564		else
565		{
566		transformedPoint = point;
567		}
568
569		if ( m_CoefficientImage[0] )
570		{
571
572		ContinuousIndexType index;
573		for ( j = 0; j < SpaceDimension; j++ )
574		{
575		index[j] = ( point[j] - m_GridOrigin[j] ) / m_GridSpacing[j];
576		}
577
578		// NOTE: if the support region does not lie totally within the grid
579		// we assume zero displacement and return the input point
580		inside = this->InsideValidRegion( index );
581		if ( !inside )
582		{
583		outputPoint = transformedPoint;
584		return;
585		}
586
587		// Compute interpolation weights
588		m_WeightsFunction->Evaluate( index, weights, supportIndex );
589
590		// For each dimension, correlate coefficient with weights
591		RegionType supportRegion;
592		supportRegion.SetSize( m_SupportSize );
593		supportRegion.SetIndex( supportIndex );
594
595		outputPoint.Fill( NumericTraits<ScalarType>::Zero );
596
597		typedef ImageRegionConstIterator<ImageType> IteratorType;
598		IteratorType m_Iterator[ SpaceDimension ];
599		unsigned long counter = 0;
600		const PixelType * basePointer = m_CoefficientImage[0]->GetBufferPointer();
601
602		for ( j = 0; j < SpaceDimension; j++ )
603		{
604		m_Iterator[j] = IteratorType( m_CoefficientImage[j], supportRegion );
605		}
606
607		while ( ! m_Iterator[0].IsAtEnd() )
608		{
609
610		// multiply weigth with coefficient
611		for ( j = 0; j < SpaceDimension; j++ )
612		{
613		outputPoint[j] += static_cast<ScalarType>(
614		weights[counter] * m_Iterator[j].Get());
615		}
616
617		// populate the indices array
618		indices[counter] = &(m_Iterator[0].Value()) - basePointer;
619
620		// go to next coefficient in the support region
621		++ counter;
622		for ( j = 0; j < SpaceDimension; j++ )
623		{
624		++( m_Iterator[j] );
625		}
626		}
627
628		// return results
629		for ( j = 0; j < SpaceDimension; j++ )
630		{
631		outputPoint[j] += transformedPoint[j];
632		}
633
634		}
635		else
636		{
637
638		itkWarningMacro( << "B-spline coefficients have not been set" );
639
640		for ( j = 0; j < SpaceDimension; j++ )
641		{
642		outputPoint[j] = transformedPoint[j];
643		}
644
645		}
646
647		}
648
649
650	EML
651		// Transform a point
652		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
653		typename BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
654		::OutputPointType
655		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
656		::TransformPoint(const InputPointType &point) const
657		{
658
659		WeightsType weights( m_WeightsFunction->GetNumberOfWeights() );
660		ParameterIndexArrayType indices( m_WeightsFunction->GetNumberOfWeights() );
661		OutputPointType outputPoint;
662		bool inside;
663
664		this->TransformPoint( point, outputPoint, weights, indices, inside );
665
666		return outputPoint;
667
668		}
669
670
671		// Compute the Jacobian in one position
672		template<class TScalarType, unsigned int NDimensions, unsigned int VSplineOrder>
673		const
674		typename BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
675		::JacobianType &
676		BSplineDeformableTransform<TScalarType, NDimensions,VSplineOrder>
677		::GetJacobian( const InputPointType & point ) const
678		{
679		// Can only compute Jacobian if parameters are set via
680		// SetParameters or SetParametersByValue
681		if( m_InputParametersPointer == NULL )
682		{
683		itkExceptionMacro( <<"Cannot compute Jacobian: parameters not set" );
684		}
685
686		// Zero all components of jacobian
687		// NOTE: for efficiency, we only need to zero out the coefficients
688		// that got fill last time this method was called.
689		RegionType supportRegion;
690		supportRegion.SetSize( m_SupportSize );
691		supportRegion.SetIndex( m_LastJacobianIndex );
692
693		typedef ImageRegionIterator<JacobianImageType> IteratorType;
694		IteratorType m_Iterator[ SpaceDimension ];
695		unsigned int j;
696
697		for ( j = 0; j < SpaceDimension; j++ )
698		{
699		m_Iterator[j] = IteratorType( m_JacobianImage[j], supportRegion );
700		}
701
702		while ( ! m_Iterator[0].IsAtEnd() )
703		{
704
705		// zero out jacobian elements
706		for ( j = 0; j < SpaceDimension; j++ )
707		{
708		m_Iterator[j].Set( NumericTraits<JacobianPixelType>::Zero );
709		}
710
711		for ( j = 0; j < SpaceDimension; j++ )
712		{
713		++( m_Iterator[j] );
714		}
715		}
716
717
718		ContinuousIndexType index;
719		for ( j = 0; j < SpaceDimension; j++ )
720		{
721		index[j] = ( point[j] - m_GridOrigin[j] ) / m_GridSpacing[j];
722		}
723
724		// NOTE: if the support region does not lie totally within the grid
725		// we assume zero displacement and return the input point
726		if ( !this->InsideValidRegion( index ) )
727		{
728		return this->m_Jacobian;
729		}
730
731		// Compute interpolation weights
732		WeightsType weights( m_WeightsFunction->GetNumberOfWeights() );
733		IndexType supportIndex;
734
735		m_WeightsFunction->Evaluate( index, weights, supportIndex );
736		m_LastJacobianIndex = supportIndex;
737
738		// For each dimension, copy the weight to the support region
739		supportRegion.SetIndex( supportIndex );
740		unsigned long counter = 0;
741
742		for ( j = 0; j < SpaceDimension; j++ )
743		{
744		m_Iterator[j] = IteratorType( m_JacobianImage[j], supportRegion );
745		}
746
747		while ( ! m_Iterator[0].IsAtEnd() )
748		{
749
750		// copy weight to jacobian image
751		for ( j = 0; j < SpaceDimension; j++ )
752		{
753		m_Iterator[j].Set( static_cast<JacobianPixelType>( weights[counter] ) );
754		}
755
756		// go to next coefficient in the support region
757		++ counter;
758		for ( j = 0; j < SpaceDimension; j++ )
759		{
760		++( m_Iterator[j] );
761		}
762		}
763
764
765		// Return the results
766		return this->m_Jacobian;
767
768		}
769
770
771
772		} // namespace
773
774		#endif
775

---

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

© Kitware Inc.