1. /***************************************************************************
  2.  *   Copyright (c) Johannes Perl                                           *
  3.  *   All rights reserved.                                                  *
  4.  *                                                                         *
  5.  *   This program is free software; you can redistribute it and/or modify  *
  6.  *   it under the terms of the GNU General Public License as published by  *
  7.  *   the Free Software Foundation; either version 2 of the License, or     *
  8.  *   (at your option) any later version.                                   *
  9.  ***************************************************************************/
  10.  
  11. #ifndef _SheppLoganFilter_txx
  12. #define _SheppLoganFilter_txx
  13.  
  14. #include "SheppLoganFilter.h"
  15.  
  16. template < typename TInputImage, typename TOutputImage >
  17. SheppLoganFilter< TInputImage, TOutputImage >
  18. ::SheppLoganFilter()
  19. {
  20. 	m_VRadius = 0;
  21. }
  22.  
  23. template < typename TInputImage, typename TOutputImage >
  24. void SheppLoganFilter< TInputImage, TOutputImage >
  25. ::ThreadedGenerateData(const RegionType& 
  26. 					 regionForThread, int threadId )
  27. {
  28. 	if(m_VRadius <= 0)
  29. 		itkExceptionMacro(<< "Member variable VRadius not properly initialized.");
  30.  
  31. 	//pointers to output and input image
  32. 	ImageConstPointer inputImage  = this->GetInput();
  33. 	ImagePointer outputImage = this->GetOutput();
  34.  
  35. 	RegionType inputRegion = inputImage->GetLargestPossibleRegion();
  36.  
  37. 	RegionType outputRegion = outputImage->GetLargestPossibleRegion();
  38. 	IteratorType out(outputImage, outputRegion);
  39.  
  40. 	NeighborhoodIteratorType::RadiusType radius;
  41. 	radius.Fill(0);
  42. 	radius[1] = m_VRadius;
  43. 	int vDim = m_VRadius*2+1;
  44. 	NeighborhoodIteratorType it(radius, inputImage, inputRegion);
  45.  
  46. 	//pointer to Shepp-Logan filter weights
  47. 	float *sl = new float[vDim];
  48. 	//vector with the OffsetTypes for the neighborhood iterator
  49. 	std::vector< NeighborhoodIteratorType::OffsetType > offset;
  50.  
  51. 	int actIndex;
  52. 	for(int i = -m_VRadius; i <= m_VRadius; ++i)
  53. 	{
  54. 		actIndex = i + m_VRadius;
  55. 		*(sl+actIndex) = -2.0/(PI*PI*(4*i*i-1));
  56. 		NeighborhoodIteratorType::OffsetType tmpO = {{0,i}};
  57. 		offset.push_back(tmpO);
  58. 	}
  59.  
  60. 	std::vector< int > tmpImg;
  61.  
  62. 	for (it.GoToBegin(), out.GoToBegin(); !it.IsAtEnd(); ++it, ++out)
  63. 	{
  64. 		float sum = 0;
  65. 		for(int i = 0; i < m_VRadius*2+1; ++i)
  66. 		{
  67. 			sum += *(sl+i) * it.GetPixel(offset[i]);
  68. 		}
  69.  
  70. 		tmpImg.push_back(sum);
  71. 	}
  72.  
  73. 	int i = 0;
  74. 	for (it.GoToBegin(), out.GoToBegin(); !it.IsAtEnd(); ++it, ++out)
  75. 	{
  76. 		out.Set(tmpImg[i++]);
  77. 	}
  78. }
  79.  
  80. template < typename TInputImage, typename TOutputImage>
  81. void SheppLoganFilter< TInputImage, TOutputImage >
  82. ::PrintSelf(std::ostream& os, itk::Indent indent) const
  83. {
  84. 	Superclass::PrintSelf(os,indent);
  85. 	os << "VRadius = " << m_VRadius << std::endl;
  86. }
  87.  
  88. #endif
  89.