KWStyle - itkBSplineInterpolationWeightFunction.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkBSplineInterpolationWeightFunction.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 _itkBSplineInterpolationWeightFunction_txx
18	DEF	#define _itkBSplineInterpolationWeightFunction_txx
19
20		#include "itkBSplineInterpolationWeightFunction.h"
21		#include "itkImage.h"
22		#include "itkMatrix.h"
23		#include "itkImageRegionConstIteratorWithIndex.h"
24
25		// anonymous namespace
26	NMS	namespace
27		{
28		//--------------------------------------------------------------------------
29		// The 'floor' function on x86 and mips is many times slower than these
30		// and is used a lot in this code, optimize for different CPU architectures
31		inline int BSplineFloor(double x)
32		{
33		#if defined mips || defined sparc || defined __ppc__
34		return (int)((unsigned int)(x + 2147483648.0) - 2147483648U);
35		#elif defined i386 || defined _M_IX86
36		union { unsigned int hilo[2]; double d; } u;
37		u.d = x + 103079215104.0;
38		return (int)((u.hilo[1]<<16)|(u.hilo[0]>>16));
39	IND	#else
40		return int(floor(x));
41		#endif
42		}
43
44		}
45
46		namespace itk
47		{
48
49		/**
50		* Constructor
51		*/
52	LEN	template <class TCoordRep, unsigned int VSpaceDimension, unsigned int VSplineOrder>
53		BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>
54		::BSplineInterpolationWeightFunction()
55		{
56
57		// Initialize the number of weights;
58		m_NumberOfWeights =
59	IND	****static_cast<unsigned long>( pow( static_cast<double>(SplineOrder + 1),
60	IND	*************************************static_cast<double>(SpaceDimension) ) );
61
62		// Initialize support region is a hypercube of length SplineOrder + 1
63		m_SupportSize.Fill( SplineOrder + 1 );
64
65		// Initialize offset to index lookup table
66		m_OffsetToIndexTable.set_size( m_NumberOfWeights, SpaceDimension );
67
68		typedef Image<char,SpaceDimension> CharImageType;
69		typename CharImageType::Pointer tempImage = CharImageType::New();
70		tempImage->SetRegions( m_SupportSize );
71		tempImage->Allocate();
72		tempImage->FillBuffer( 0 );
73
74
75		typedef ImageRegionConstIteratorWithIndex<CharImageType> IteratorType;
76		IteratorType iterator( tempImage, tempImage->GetBufferedRegion() );
77		unsigned long counter = 0;
78
79		while ( !iterator.IsAtEnd() )
80		{
81		for(unsigned  int j = 0; j < SpaceDimension; j++ )
82		{
83		m_OffsetToIndexTable[counter][j] = iterator.GetIndex()[j];
84		}
85		++counter;
86		++iterator;
87		}
88
89
90		// Initialize the interpolation kernel
91		m_Kernel = KernelType::New();
92
93		}
94
95
96		/**
97		* Standard "PrintSelf" method
98		*/
99	LEN	template <class TCoordRep, unsigned int VSpaceDimension, unsigned int VSplineOrder>
100		void
101		BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>
102		::PrintSelf(
103		std::ostream& os,
104		Indent indent) const
105		{
106		Superclass::PrintSelf( os, indent );
107
108		os << indent << "NumberOfWeights: " << m_NumberOfWeights << std::endl;
109		os << indent << "SupportSize: " << m_SupportSize << std::endl;
110		}
111
112
113	EML
114		/**
115		* Compute weights for interpolation at continous index position
116		*/
117	LEN	template <class TCoordRep, unsigned int VSpaceDimension, unsigned int VSplineOrder>
118	LEN	typename BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>
119		::WeightsType
120		BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>
121		::Evaluate(
122		const ContinuousIndexType& index ) const
123		{
124
125		WeightsType weights( m_NumberOfWeights );
126		IndexType startIndex;
127
128		this->Evaluate( index, weights, startIndex );
129
130		return weights;
131		}
132
133
134	EML
135		/**
136		* Compute weights for interpolation at continous index position
137		*/
138	LEN	template <class TCoordRep, unsigned int VSpaceDimension, unsigned int VSplineOrder>
139	LEN	void BSplineInterpolationWeightFunction<TCoordRep, VSpaceDimension, VSplineOrder>
140		::Evaluate(
141		const ContinuousIndexType & index,
142		WeightsType & weights,
143		IndexType & startIndex ) const
144		{
145
146		unsigned int j, k;
147
148		// Find the starting index of the support region
149		for ( j = 0; j < SpaceDimension; j++ )
150		{
151		startIndex[j] = static_cast<typename IndexType::IndexValueType>(
152		BSplineFloor( index[j] - static_cast<double>( SplineOrder / 2 ) ) );
153		}
154
155		// Compute the weights
156		Matrix<double,SpaceDimension,SplineOrder + 1> weights1D;
157		for ( j = 0; j < SpaceDimension; j++ )
158		{
159		double x = index[j] - static_cast<double>( startIndex[j] );
160
161		for( k = 0; k <= SplineOrder; k++ )
162		{
163		weights1D[j][k] = m_Kernel->Evaluate( x );
164		x -= 1.0;
165		}
166		}
167
168		for ( k = 0; k < m_NumberOfWeights; k++ )
169		{
170
171		weights[k] = 1.0;
172
173		for ( j = 0; j < SpaceDimension; j++ )
174		{
175		weights[k] *= weights1D[j][ m_OffsetToIndexTable[k][j] ];
176		}
177
178		}
179
180		}
181
182
183		} // end namespace itk
184
185		#endif
186
187	EOF

---

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

© Kitware Inc.