KWStyle - itkGaussianBlurImageFunction.txx

Matrix View

Description

---

1		/*=========================================================================
2
3		Program:   Insight Segmentation & Registration Toolkit
4		Module:    $RCSfile: itkGaussianBlurImageFunction.txx.html,v $
5		Language:  C++
6		Date:      $Date: 2006/01/17 19:15:36 $
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 _itkGaussianBlurImageFunction_txx
18	DEF	#define _itkGaussianBlurImageFunction_txx
19
20		#include "itkGaussianBlurImageFunction.h"
21		#include "itkImageLinearIteratorWithIndex.h"
22
23		namespace itk
24		{
25
26		/** Set the Input Image */
27		template <class TInputImage,class TOutput>
28		GaussianBlurImageFunction<TInputImage,TOutput>
29		::GaussianBlurImageFunction()
30		{
31		typename GaussianFunctionType::ArrayType mean;
32		mean[0]=0.0f;
33		for(unsigned int i=0;i<itkGetStaticConstMacro(ImageDimension);i++)
34		{
35		m_Sigma[i] = 1.0f;
36		m_MaximumError[i] = 0.001f;
37		m_MaximumKernelWidth = 32;
38		m_Extent[i] = 1.0f;
39		}
40		m_UseImageSpacing = true;
41
42		m_GaussianFunction = GaussianFunctionType::New();
43		m_GaussianFunction->SetMean(mean);
44		m_GaussianFunction->SetNormalized(false); // faster
45		m_OperatorImageFunction = OperatorImageFunctionType::New();
46		m_InternalImage = InternalImageType::New();
47		this->RecomputeGaussianKernel();
48		m_Caster = CastImageFilterType::New();
49		}
50
51		/** Set the input image */
52		template <class TInputImage,class TOutput>
53		void
54		GaussianBlurImageFunction<TInputImage,TOutput>
55		::SetInputImage( const InputImageType * ptr )
56		{
57		Superclass::SetInputImage(ptr);
58		m_Caster->SetInput(ptr);
59		m_Caster->Update();
60		m_OperatorImageFunction->SetInputImage(m_Caster->GetOutput());
61		}
62
63
64		/** Print self method */
65		template <class TInputImage,class TOutput>
66		void
67		GaussianBlurImageFunction<TInputImage,TOutput>
68		::PrintSelf(std::ostream& os, Indent indent) const
69		{
70		this->Superclass::PrintSelf(os,indent);
71
72		for(unsigned int i=0;i<itkGetStaticConstMacro(ImageDimension);i++)
73		{
74		os << indent << "Sigma["<< i << "] : " <<  m_Sigma[i] << std::endl;
75	LEN	os << indent << "MaximumError["<< i << "] : " << m_MaximumError[i] << std::endl;
76		os << indent << "Extent["<< i << "] : " << m_Extent[i] << std::endl;
77		}
78		os << indent << "MaximumKernelWidth: " << m_MaximumKernelWidth << std::endl;
79		os << indent << "UseImageSpacing: " << m_UseImageSpacing << std::endl;
80
81		os << indent << "Internal Image : " << m_InternalImage << std::endl;
82		os << indent << "Image Caster : " << m_Caster << std::endl;
83
84		}
85
86		/** Set the variance of the gaussian in each direction */
87		template <class TInputImage,class TOutput>
88		void
89		GaussianBlurImageFunction<TInputImage,TOutput>
90		::SetSigma(const double* sigma)
91		{
92		unsigned int i;
93		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
94		{
95		if ( sigma[i] != m_Sigma[i] )
96		{
97		break;
98		}
99		}
100		if ( i < itkGetStaticConstMacro(ImageDimension) )
101		{
102		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
103		{
104		m_Sigma[i] = sigma[i];
105		}
106		this->RecomputeGaussianKernel();
107		}
108		}
109
110
111		/** Set the variance of the gaussian in each direction */
112		template <class TInputImage,class TOutput>
113		void
114		GaussianBlurImageFunction<TInputImage,TOutput>
115		::SetSigma( const double sigma)
116		{
117		unsigned int i;
118		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
119		{
120		if ( sigma != m_Sigma[i] )
121		{
122		break;
123		}
124		}
125		if ( i < itkGetStaticConstMacro(ImageDimension) )
126		{
127		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
128		{
129		m_Sigma[i] = sigma;
130		}
131		this->RecomputeGaussianKernel();
132		}
133		}
134
135		/** Set the extent of the gaussian in each direction */
136		template <class TInputImage,class TOutput>
137		void
138		GaussianBlurImageFunction<TInputImage,TOutput>
139		::SetExtent(const double* extent)
140		{
141		unsigned int i;
142		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
143		{
144		if ( extent[i] != m_Extent[i] )
145		{
146		break;
147		}
148		}
149		if ( i < itkGetStaticConstMacro(ImageDimension) )
150		{
151		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
152		{
153		m_Extent[i] = extent[i];
154		}
155		this->RecomputeGaussianKernel();
156		}
157		}
158
159
160		/** Set the extent of the gaussian in each direction */
161		template <class TInputImage,class TOutput>
162		void
163		GaussianBlurImageFunction<TInputImage,TOutput>
164		::SetExtent( const double extent)
165		{
166		unsigned int i;
167		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
168		{
169		if ( extent != m_Extent[i] )
170		{
171		break;
172		}
173		}
174		if ( i < itkGetStaticConstMacro(ImageDimension) )
175		{
176		for (i=0; i<itkGetStaticConstMacro(ImageDimension); i++)
177		{
178		m_Extent[i] = extent;
179		}
180		this->RecomputeGaussianKernel();
181		}
182		}
183
184		/** Recompute the gaussian kernel used to evaluate indexes
185		*  And allocate the internal image for processing depending on
186		*  the size of the operator */
187		template <class TInputImage,class TOutput>
188		void
189		GaussianBlurImageFunction<TInputImage,TOutput>
190		::RecomputeGaussianKernel()
191		{
192
193		typename InternalImageType::SizeType size;
194		// Compute the convolution of each kernel in each direction
195	LEN	for(unsigned int direction=0;direction<itkGetStaticConstMacro(ImageDimension);direction++)
196		{
197		GaussianOperatorType gaussianOperator;
198
199		gaussianOperator.SetDirection(direction);
200		gaussianOperator.SetMaximumError(m_MaximumError[direction]);
201		gaussianOperator.SetMaximumKernelWidth(m_MaximumKernelWidth);
202
203		if( (m_UseImageSpacing == true) && (this->GetInputImage()) )
204		{
205		if (this->GetInputImage()->GetSpacing()[direction] == 0.0)
206		{
207		itkExceptionMacro(<< "Pixel spacing cannot be zero");
208		}
209		else
210		{
211	LEN	gaussianOperator.SetVariance(m_Sigma[direction]*m_Sigma[direction]  / this->GetInputImage()->GetSpacing()[direction]);
212		}
213		}
214		else
215		{
216		gaussianOperator.SetVariance(m_Sigma[direction]*m_Sigma[direction]);
217		}
218
219		gaussianOperator.CreateDirectional();
220		m_OperatorArray[direction] = gaussianOperator;
221		size[direction] = gaussianOperator.GetSize()[direction];
222		}
223
224		// Allocate the internal image
225		m_InternalImage = InternalImageType::New();
226		typename InternalImageType::RegionType region;
227		region.SetSize(size);
228		m_InternalImage->SetRegions(region);
229		m_InternalImage->Allocate();
230		m_InternalImage->FillBuffer(0);
231		}
232
233		/** Evaluate the function at the specifed point */
234		template <class TInputImage,class TOutput>
235		TOutput
236		GaussianBlurImageFunction<TInputImage,TOutput>
237		::EvaluateAtIndex(const IndexType& index) const
238		{
239		// First time we use the complete image and fill the internal image
240		m_OperatorImageFunction->SetInputImage(m_Caster->GetOutput());
241		m_OperatorImageFunction->SetOperator(m_OperatorArray[0]);
242
243		// if 1D Image we return the result
244		if(itkGetStaticConstMacro(ImageDimension) == 1)
245		{
246		return m_OperatorImageFunction->EvaluateAtIndex(index);
247		}
248
249		// Compute the centered index fo the neighborhood
250		IndexType centerIndex;
251		for(unsigned int i=0;i<itkGetStaticConstMacro(ImageDimension);i++)
252		{
253	LEN	centerIndex[i] = (unsigned long)((float)m_InternalImage->GetBufferedRegion().GetSize()[i]/2.0);
254		}
255
256		// first direction
257		typename InternalImageType::IndexType ind;
258		ind = index;
259
260		//Define the region of the iterator
261		typename InternalImageType::RegionType region;
262	LEN	typename InternalImageType::SizeType size = m_InternalImage->GetBufferedRegion().GetSize();
263		size[0]=1;
264		region.SetSize(size);
265
266		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
267		{
268		if(i != 0)
269		{
270		ind[i] -= centerIndex[i];
271		}
272		}
273		region.SetIndex(ind);
274
275		typename InternalImageType::RegionType regionN;
276		regionN.SetSize(size);
277		ind = centerIndex;
278		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
279		{
280		if(i != 0)
281		{
282		ind[i] = 0;
283		}
284		}
285		regionN.SetIndex(ind);
286
287	LEN	itk::ImageLinearConstIteratorWithIndex<InternalImageType> it(m_Caster->GetOutput(),region);
288	LEN	itk::ImageLinearIteratorWithIndex<InternalImageType> itN(m_InternalImage,regionN);
289		it.SetDirection(1);
290		itN.SetDirection(1);
291		it.GoToBeginOfLine();
292		itN.GoToBeginOfLine();
293		while(!it.IsAtEnd())
294		{
295		while(!it.IsAtEndOfLine())
296		{
297	LEN	if(m_Caster->GetOutput()->GetLargestPossibleRegion().IsInside(it.GetIndex()))
298		{
299		itN.Set(m_OperatorImageFunction->EvaluateAtIndex(it.GetIndex()));
300		}
301		++it;
302		++itN;
303		}
304		it.NextLine();
305		itN.NextLine();
306		}
307
308		// Do the convolution in other directions
309	LEN	for(unsigned int direction=1;direction<itkGetStaticConstMacro(ImageDimension);direction++)
310		{
311
312		size[direction] = 1;
313		ind = centerIndex;
314		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
315		{
316		if(i > direction)
317		{
318		ind[i] = 0;
319		}
320		}
321		region.SetSize(size);
322		region.SetIndex(ind);
323
324
325		m_OperatorImageFunction->SetInputImage(m_InternalImage);
326		m_OperatorImageFunction->SetOperator(m_OperatorArray[direction]);
327
328	LEN	itk::ImageLinearIteratorWithIndex<InternalImageType> it(m_InternalImage,region);
329
330		unsigned int dir = direction +1;
331		if(dir == itkGetStaticConstMacro(ImageDimension))
332		{
333		dir = itkGetStaticConstMacro(ImageDimension)-1;
334		}
335
336		it.SetDirection(dir);
337		it.GoToBeginOfLine();
338		while(!it.IsAtEnd())
339		{
340		while(!it.IsAtEndOfLine())
341		{
342		it.Set(m_OperatorImageFunction->EvaluateAtIndex(it.GetIndex()));
343		++it;
344		}
345		it.NextLine();
346		}
347		}
348
349		return m_InternalImage->GetPixel(centerIndex);
350
351		}
352
353
354		/** Recompute the gaussian kernel used to evaluate indexes
355		*  The variance should be uniform */
356		template <class TInputImage,class TOutput>
357		void
358		GaussianBlurImageFunction<TInputImage,TOutput>
359		::RecomputeContinuousGaussianKernel(const double* offset) const
360		{
361	LEN	for(unsigned int direction=0;direction<itkGetStaticConstMacro(ImageDimension);direction++)
362		{
363		NeighborhoodType gaussianNeighborhood;
364		typename GaussianFunctionType::InputType pt;
365		typename NeighborhoodType::SizeType size;
366		size.Fill(0);
367		size[direction] = (unsigned long)(m_Sigma[direction]*m_Extent[direction]);
368
369		gaussianNeighborhood.SetRadius(size);
370
371		typename NeighborhoodType::Iterator it = gaussianNeighborhood.Begin();
372
373		itk::FixedArray<double,1> s;
374		s[0]=m_Sigma[direction];
375		m_GaussianFunction->SetSigma(s);
376
377		unsigned int i=0;
378		float sum = 0;
379		while(it != gaussianNeighborhood.End() )
380		{
381		pt[0]= gaussianNeighborhood.GetOffset(i)[direction]-offset[direction];
382		if( (m_UseImageSpacing == true) && (this->GetInputImage()) )
383		{
384		if (this->GetInputImage()->GetSpacing()[direction] == 0.0)
385		{
386		itkExceptionMacro(<< "Pixel spacing cannot be zero");
387		}
388		else
389		{
390		pt[0] *= this->GetInputImage()->GetSpacing()[direction];
391		}
392		}
393
394		(*it)= m_GaussianFunction->Evaluate(pt);
395		sum += (*it);
396		i++;
397		it++;
398		}
399
400		// Make the filter DC-Constant
401		it = gaussianNeighborhood.Begin();
402		while(it != gaussianNeighborhood.End() )
403		{
404		(*it) /= sum;
405		it++;
406		}
407		m_ContinuousOperatorArray[direction] = gaussianNeighborhood;
408		}
409		}
410
411		/** Evaluate the function at the specifed point */
412		template <class TInputImage,class TOutput>
413		TOutput
414		GaussianBlurImageFunction<TInputImage,TOutput>
415		::Evaluate(const PointType& point) const
416		{
417
418		IndexType index;
419
420		double offset[itkGetStaticConstMacro(ImageDimension)];
421		for(unsigned int i=0; i<itkGetStaticConstMacro(ImageDimension);i++)
422		{
423		index[i] = (unsigned long)point[i];
424		offset[i] = point[i]-index[i];
425		}
426
427		this->RecomputeContinuousGaussianKernel(offset);
428
429		// First time we use the complete image and fill the internal image
430		m_OperatorImageFunction->SetInputImage(m_Caster->GetOutput());
431		m_OperatorImageFunction->SetOperator(m_ContinuousOperatorArray[0]);
432
433		// if 1D Image we return the result
434		if(itkGetStaticConstMacro(ImageDimension) == 1)
435		{
436		return m_OperatorImageFunction->EvaluateAtIndex(index);
437		}
438
439		// Compute the centered index fo the neighborhood
440		IndexType centerIndex;
441		for(unsigned int i=0;i<itkGetStaticConstMacro(ImageDimension);i++)
442		{
443	LEN	centerIndex[i] = (unsigned long)((float)m_InternalImage->GetBufferedRegion().GetSize()[i]/2.0);
444		}
445
446		// first direction
447		typename InternalImageType::IndexType ind;
448		ind = index;
449
450		//Define the region of the iterator
451		typename InternalImageType::RegionType region;
452	LEN	typename InternalImageType::SizeType size = m_InternalImage->GetBufferedRegion().GetSize();
453		size[0]=1;
454		region.SetSize(size);
455
456		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
457		{
458		if(i != 0)
459		{
460		ind[i] -= centerIndex[i];
461		}
462		}
463		region.SetIndex(ind);
464
465		typename InternalImageType::RegionType regionN;
466		regionN.SetSize(size);
467		ind = centerIndex;
468		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
469		{
470		if(i != 0)
471		{
472		ind[i] = 0;
473		}
474		}
475		regionN.SetIndex(ind);
476
477	LEN	itk::ImageLinearConstIteratorWithIndex<InternalImageType> it(m_Caster->GetOutput(),region);
478	LEN	itk::ImageLinearIteratorWithIndex<InternalImageType> itN(m_InternalImage,regionN);
479		it.SetDirection(1);
480		itN.SetDirection(1);
481		it.GoToBeginOfLine();
482		itN.GoToBeginOfLine();
483		while(!it.IsAtEnd())
484		{
485		while(!it.IsAtEndOfLine())
486		{
487	LEN	if(m_Caster->GetOutput()->GetLargestPossibleRegion().IsInside(it.GetIndex()))
488		{
489		itN.Set(m_OperatorImageFunction->EvaluateAtIndex(it.GetIndex()));
490		}
491		++it;
492		++itN;
493		}
494		it.NextLine();
495		itN.NextLine();
496		}
497
498		// Do the convolution in other directions
499	LEN	for(unsigned int direction=1;direction<itkGetStaticConstMacro(ImageDimension);direction++)
500		{
501
502		size[direction] = 1;
503		ind = centerIndex;
504		for(unsigned int i = 0;i<itkGetStaticConstMacro(ImageDimension);i++)
505		{
506		if(i > direction)
507		{
508		ind[i] = 0;
509		}
510		}
511		region.SetSize(size);
512		region.SetIndex(ind);
513
514
515		m_OperatorImageFunction->SetInputImage(m_InternalImage);
516		m_OperatorImageFunction->SetOperator(m_ContinuousOperatorArray[direction]);
517
518	LEN	itk::ImageLinearIteratorWithIndex<InternalImageType> it(m_InternalImage,region);
519
520		unsigned int dir = direction +1;
521		if(dir == itkGetStaticConstMacro(ImageDimension))
522		{
523		dir = itkGetStaticConstMacro(ImageDimension)-1;
524		}
525
526		it.SetDirection(dir);
527		it.GoToBeginOfLine();
528		while(!it.IsAtEnd())
529		{
530		while(!it.IsAtEndOfLine())
531		{
532		it.Set(m_OperatorImageFunction->EvaluateAtIndex(it.GetIndex()));
533		++it;
534		}
535		it.NextLine();
536		}
537		}
538
539		return m_InternalImage->GetPixel(centerIndex);
540		}
541
542		/** Evaluate the function at specified ContinousIndex position.*/
543		template <class TInputImage,class TOutput>
544		TOutput
545		GaussianBlurImageFunction<TInputImage,TOutput>
546		::EvaluateAtContinuousIndex( const ContinuousIndexType & index ) const
547		{
548		PointType point;
549		for(unsigned int i=0; i<itkGetStaticConstMacro(ImageDimension);i++)
550		{
551		point[i] = index[i];
552		}
553		return this->Evaluate(point);
554
555		}
556
557		} // namespace itk
558
559		#endif
560

---

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

© Kitware Inc.