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. /*!
  12.  * bp2d.cpp
  13.  * \author Johannes Perl
  14.  * \brief Performs a unfiltered backprojection on raw CT data.
  15.  */
  16.  
  17. #include "bp2d.h"
  18.  
  19. using namespace std;
  20.  
  21. int main(int argc, char* argv[])
  22. {
  23. 	string filename("images/radon180.png");
  24. 	float pStart = 0;
  25. 	float pEnd = 179;
  26. 	float pDist = 1;
  27. 	bp2d(filename, pStart, pEnd, pDist);
  28. }
  29.  
  30. int bp2d(string filename, float pStart, float pEnd, float pDist)
  31. {
  32. 	//defining image properties
  33. 	typedef unsigned int uint;
  34. 	typedef unsigned char uchar;
  35. 	typedef uchar PixelType;
  36. 	const int Dimension = 2;
  37. 	typedef itk::Image< PixelType, 2 > ImageType;
  38. 	typedef itk::ImageRegionIterator< ImageType > IteratorType;
  39.  
  40. 	//smartpointer to image to be read
  41. 	ImageType::Pointer image = ImageType::New();
  42.  
  43. 	//reading image from disk
  44. 	typedef itk::ImageFileReader< ImageType > ReaderType;
  45. 	ReaderType::Pointer reader = ReaderType::New();
  46. 	reader->SetFileName(filename);
  47. 	try
  48. 	{
  49. 		reader->Update();
  50. 	}
  51. 	catch( itk::ExceptionObject exp )
  52. 	{
  53. 		cerr << "Reading image " << filename << " caused a problem." << endl;
  54. 		cerr << exp << endl;
  55. 		return 1;
  56. 	}
  57.  
  58. 	image = reader->GetOutput();
  59.  
  60. 	ImageType::RegionType region = image->GetLargestPossibleRegion();
  61. 	ImageType::SizeType size = region.GetSize();
  62. 	unsigned int w = size[0];
  63. 	unsigned int h = size[1];
  64.  
  65. 	unsigned int bpImageDim = h;
  66.  
  67. 	cout << "Input image dimensions" << endl;
  68. 	cout << h << " " << w << endl;
  69. 	cout << "Output image dimensions" << endl;
  70. 	cout << h << " x " << h << endl;
  71.  
  72. 	ImageType::Pointer bpImage =  ImageType::New();
  73. 	ImageType::RegionType outputRegion;
  74. 	size[0] = size[1] = bpImageDim;
  75. 	outputRegion.SetSize(size);
  76. 	ImageType::IndexType start;
  77. 	start[0] = start[1] = 0;
  78. 	outputRegion.SetIndex(start);
  79. 	bpImage->SetRegions(outputRegion);
  80. 	bpImage->Allocate();
  81.  
  82. 	double middle = (bpImageDim+1) / 2;
  83.  
  84. 	//projection index
  85. 	ImageType::IndexType pIndex;
  86. 	pIndex[0] = 0; // x position
  87. 	pIndex[1] = 0; // y position
  88.  
  89. 	//vector for temporarily storing pixel values
  90. 	vector< vector < long > > tmpImg(bpImageDim, std::vector< long >(bpImageDim, 0));
  91.  
  92. 	int numP = 0;
  93.  
  94. 	//loop foreach projection
  95. 	for(float p = pStart; p <= pEnd; p += pDist)
  96. 	{
  97. 		++numP;
  98. 		pIndex[0] = p;
  99. 		double angle = p * PI / 180.0; //convert angle to radians
  100.  
  101. 		cout << "angle " << p << " degrees" << endl;
  102.  
  103. 		int pInd; //projection index
  104.  
  105. 		for(int x = 0; x < bpImageDim; ++x)
  106. 		{
  107. 			for(int y = 0; y < bpImageDim; ++y)
  108. 			{
  109. 				pInd = (middle + (x - middle) * cos(angle) - (y - middle) * sin(angle)) + 0.5;
  110. 				if(pInd > 0 && pInd < bpImageDim)
  111. 				{
  112. 					pIndex[1] = pInd;
  113. 					tmpImg[x][y] += image->GetPixel(pIndex);
  114. 				}
  115. 			}
  116. 		}
  117. 	}
  118.  
  119. 	//write acquired data into image normed by the number of projections
  120. 	IteratorType it(bpImage, bpImage->GetLargestPossibleRegion());
  121. 	for(it.GoToBegin(); !it.IsAtEnd(); ++it)
  122. 	{
  123. 		long x = it.GetIndex()[0];
  124. 		long y = it.GetIndex()[1];
  125. 		it.Set(tmpImg[x][y] / numP);
  126. 	}
  127.  
  128. 	//write image to file
  129. 	typedef itk::ImageFileWriter< ImageType > WriterType;
  130. 	WriterType::Pointer writer = WriterType::New();
  131. 	ostringstream fname;
  132. 	fname.fill('0');
  133. 	fname << "images/bpimage_" << setw(3) << static_cast<int>((pEnd - pStart) / pDist)+1 << ".png";
  134. 	writer->SetFileName(fname.str());
  135. 	writer->SetInput(bpImage);
  136.  
  137. 	try
  138. 	{
  139. 		writer->Update();
  140. 	}
  141. 	catch( itk::ExceptionObject exp )
  142. 	{
  143. 		cerr << "Exception caught !" << endl;
  144. 		cerr << exp << endl;
  145. 		return 1;
  146. 	}
  147.  
  148. 	return 0;
  149. }
  150.