KWStyle - itkWindowedSincInterpolateImageFunction.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkWindowedSincInterpolateImageFunction.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:49 $
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 _itkWindowedSincInterpolateImageFunction_txx
18	DEF	#define _itkWindowedSincInterpolateImageFunction_txx
19
20		#include "itkWindowedSincInterpolateImageFunction.h"
21
22		#include "vnl/vnl_math.h"
23
24		namespace itk
25		{
26
27		/* Constant definitions for functions */
28		namespace Function {
29
30		template<unsigned int VRadius, class TInput, class TOutput>
31		const double
32		CosineWindowFunction<VRadius, TInput, TOutput>
33		::m_Factor = vnl_math::pi / ( 2 * VRadius );
34
35		template<unsigned int VRadius, class TInput, class TOutput>
36		const double
37		HammingWindowFunction<VRadius, TInput, TOutput>
38		::m_Factor = vnl_math::pi / VRadius;
39
40		template<unsigned int VRadius, class TInput, class TOutput>
41		const double
42		WelchWindowFunction<VRadius, TInput, TOutput>
43		::m_Factor = 1.0 / ( VRadius * VRadius );
44
45		template<unsigned int VRadius, class TInput, class TOutput>
46		const double
47		LancosWindowFunction<VRadius, TInput, TOutput>
48		::m_Factor = vnl_math::pi / VRadius;
49
50		template<unsigned int VRadius, class TInput, class TOutput>
51		const double
52		BlackmanWindowFunction<VRadius, TInput, TOutput>
53		::m_Factor1 = vnl_math::pi / VRadius;
54
55		template<unsigned int VRadius, class TInput, class TOutput>
56		const double
57		BlackmanWindowFunction<VRadius, TInput, TOutput>
58		::m_Factor2 = 2.0 * vnl_math::pi / VRadius;
59
60
61	EML
62	EML
63		} // end namespace Function
64
65
66		/**
67		* Window size constant
68		*/
69
70		template<class TInputImage, unsigned int VRadius,
71		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
72		const unsigned int
73		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
74		TWindowFunction,TBoundaryCondition,TCoordRep>
75		::m_WindowSize = VRadius << 1;
76
77		/**
78		* Constructor
79		*/
80		template<class TInputImage, unsigned int VRadius,
81		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
82		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
83		TWindowFunction,TBoundaryCondition,TCoordRep>
84		::WindowedSincInterpolateImageFunction()
85		{
86		unsigned int dim;
87
88		// Compute the offset table size
89		m_OffsetTableSize = 1;
90		for(dim=0;dim<ImageDimension;dim++)
91		{
92		m_OffsetTableSize *= m_WindowSize;
93		}
94
95		// Allocate the offset table
96		m_OffsetTable = new unsigned int[m_OffsetTableSize];
97
98		// Allocate the weights tables
99		m_WeightOffsetTable = new unsigned int *[m_OffsetTableSize];
100		for(unsigned int i=0;i<m_OffsetTableSize;i++)
101		{
102		m_WeightOffsetTable[i] = new unsigned int[ImageDimension];
103		}
104		}
105
106		/**
107		* Destructor
108		*/
109		template<class TInputImage, unsigned int VRadius,
110		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
111		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
112		TWindowFunction,TBoundaryCondition,TCoordRep>
113		::~WindowedSincInterpolateImageFunction()
114		{
115		// Clear the offset table
116		delete [] m_OffsetTable;
117
118		// Clear the weights tables
119		for(unsigned int i=0; i < m_OffsetTableSize; i++)
120		{
121		delete [] m_WeightOffsetTable[i];
122		}
123		delete[] m_WeightOffsetTable;
124		}
125
126		template<class TInputImage, unsigned int VRadius,
127		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
128		void
129		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
130		TWindowFunction,TBoundaryCondition,TCoordRep>
131		::SetInputImage(const ImageType *image)
132		{
133		unsigned int dim;
134
135		// Call the parent implementation
136		Superclass::SetInputImage(image);
137
138		if( image == NULL )
139		{
140		return;
141		}
142
143		// Set the radius for the neighborhood
144		Size<ImageDimension> radius;
145		radius.Fill(VRadius);
146
147		// Initialize the neighborhood
148		IteratorType it = IteratorType(radius, image, image->GetBufferedRegion());
149
150		// Compute the offset tables (we ignore all the zero indices
151		// in the neighborhood)
152		unsigned int iOffset = 0;
153		int empty = VRadius;
154		for(unsigned int iPos = 0; iPos < it.Size(); iPos++)
155		{
156		// Get the offset (index)
157		typename IteratorType::OffsetType off = it.GetOffset(iPos);
158
159		// Check if the offset has zero weights
160		bool nonzero = true;
161		for(dim = 0; dim < ImageDimension; dim++)
162		{
163		if(off[dim] == -empty)
164		{
165		nonzero = false;
166		break;
167		}
168		}
169
170		// Only use offsets with non-zero indices
171		if(nonzero)
172		{
173		// Set the offset index
174		m_OffsetTable[iOffset] = iPos;
175
176		// Set the weight table indices
177		for(dim = 0; dim < ImageDimension; dim++)
178		{
179		m_WeightOffsetTable[iOffset][dim] = off[dim] + VRadius - 1;
180		}
181
182		// Increment the index
183		iOffset++;
184		}
185		}
186		}
187
188		/**
189		* PrintSelf
190		*/
191		template<class TInputImage, unsigned int VRadius,
192		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
193		void
194		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
195		TWindowFunction,TBoundaryCondition,TCoordRep>
196		::PrintSelf(std::ostream& os, Indent indent) const
197		{
198		this->Superclass::PrintSelf(os,indent);
199		}
200
201
202		/**
203		* Evaluate at image index position
204		*/
205		template<class TInputImage, unsigned int VRadius,
206		class TWindowFunction, class TBoundaryCondition, class TCoordRep>
207		typename WindowedSincInterpolateImageFunction<TInputImage,VRadius,
208		TWindowFunction,TBoundaryCondition,TCoordRep>
209		::OutputType
210		WindowedSincInterpolateImageFunction<TInputImage,VRadius,
211		TWindowFunction,TBoundaryCondition,TCoordRep>
212		::EvaluateAtContinuousIndex(
213		const ContinuousIndexType& index) const
214		{
215		unsigned int dim;
216		Index<ImageDimension> baseIndex;
217		double distance[ImageDimension];
218
219		// Compute the integer index based on the continuous one by
220		// 'flooring' the index
221		for( dim = 0; dim < ImageDimension; dim++ )
222		{
223		// The following "if" block is equivalent to the following line without
224		// having to call floor.
225		//    baseIndex[dim] = (long) floor( index[dim] );
226		if (index[dim] >= 0.0)
227		{
228		baseIndex[dim] = (long) index[dim];
229		}
230		else
231		{
232		long tIndex = (long) index[dim];
233		if (double(tIndex) != index[dim])
234		{
235		tIndex--;
236		}
237		baseIndex[dim] = tIndex;
238		}
239
240		distance[dim] = index[dim] - double( baseIndex[dim] );
241		}
242
243		// cout << "Sampling at index " << index << " discrete " << baseIndex << endl;
244
245		// Position the neighborhood at the index of interest
246		Size<ImageDimension> radius;
247		radius.Fill(VRadius);
248	LEN	IteratorType nit = IteratorType( radius, this->GetInputImage(), this->GetInputImage()->GetBufferedRegion());
249		nit.SetLocation( baseIndex );
250
251		// Compute the sinc function for each dimension
252		double xWeight[ImageDimension][2 * VRadius];
253		for( dim = 0; dim < ImageDimension; dim++ )
254		{
255		// x is the offset, hence the parameter of the kernel
256		double x = distance[dim] + VRadius;
257
258		// If distance is zero, i.e. the index falls precisely on the
259		// pixel boundary, the weights form a delta function.
260		if(distance[dim] == 0.0)
261		{
262		for( unsigned int i = 0; i < m_WindowSize; i++)
263		{
264		xWeight[dim][i] =
265		static_cast<int>(i) == VRadius - 1 ? 1 : 0;
266		}
267		}
268		else
269		{
270		// i is the relative offset in dimension dim.
271		for( unsigned int i = 0; i < m_WindowSize; i++)
272		{
273		// Increment the offset, taking it through the range
274		// (dist + rad - 1, ..., dist - rad), i.e. all x
275		// such that abs(x) <= rad
276		x -= 1.0;
277
278		// Compute the weight for this m
279		xWeight[dim][i] = m_WindowFunction(x) * Sinc(x);
280		}
281		}
282		}
283
284		// Iterate over the neighborhood, taking the correct set
285		// of weights in each dimension
286		double xPixelValue = 0.0;
287		for(unsigned int j = 0; j < m_OffsetTableSize; j++)
288		{
289		// Get the offset for this neighbor
290		unsigned int off = m_OffsetTable[j];
291
292		// Get the intensity value at the pixel
293		double xVal = nit.GetPixel(off);
294
295		// Multiply the intensity by each of the weights. Gotta hope
296		// that the compiler will unwrap this loop and pipeline this!
297		for(dim = 0; dim < ImageDimension; dim++)
298		{
299		xVal *= xWeight[ dim ][ m_WeightOffsetTable[j][dim] ];
300		}
301
302		// Increment the pixel value
303		xPixelValue += xVal;
304		}
305
306		// Return the interpolated value
307		return static_cast<OutputType>(xPixelValue);
308		}
309
310		} // namespace itk
311
312		#endif
313

---

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

© Kitware Inc.