Computed Tomography - Generation of the Sinogram Image in C#

• Download MakeSinogramExe.zip - 9.7 KB
• Download source - 27.7 KB

1. Introduction

In the field of medical imaging, Computed Tomography (CT), sometimes called Computerized Axial Tomography (CAT), is an important imaging modality. This is a part of structural imaging, as against functional imaging modalities like Nuclear Medicine.

This article describes the generation of the Sinogram image, given the image of a slice being imaged.

In Computed Tomography, a beam of X-rays are made to pass through an anatomical slice of interest, and the resulting X-rays are captured on a detector on the opposite side of the source. The source-detector combination is made to rotate in a full circle around the object being imaged, so as to capture a large number of projection images. This constitutes the Image Acquisition aspect of CT.

The Image Acquisition step is followed by the Reconstruction step, where special reconstruction algorithms (back-projection and filtered back-projection being a couple of them) are used to get back the anatomical details of the imaged slice.

The above-mentioned problem of CT can be seen in two ways:

1. Forward Problem: Given the anatomical details of a slice, the projection images are to be calculated. This is a relatively easy problem, as the standard rules of X-ray attenuation can be applied to generate the projection images. Any number of projection images can be generated given an image slice; these projection images differing in the angle of projection. The projection images are generated through a mathematical entity called the Radon Transform. This article is about the use of Radon Transforms to generate the Sinogram of an image. The Sinogram is an image, and we give the detailed computational steps to generate the Sinogram.
2. Inverse Problem: Given the projection images, details of the anatomical slice are to be computed. This is much harder than the Forward Problem, and involves the Inverse Radon Transform in its most basic form. As with all inverse problems, accurate reconstruction of the slice is not guaranteed to happen for all kinds of inputs, especially if there are ambiguities like missing projection values. Solution of the Inverse Problem is for a subsequent article from us here.

2. Brief Outline of the Mathematics of the Problem

The Radon Transform is related the projection function. Projections are a set of measurements of the integrated values of a parameter of an object, these integrations being done along straight lines through the object. Such integrals are called line integrals. In case of CT, the parameter which is integrated is the X-ray attenuation constant; this X-ray attenuation constant varies throughout the anatomical slice being imaged, and the projection, for a particular angle, is a function of the distribution of this attenuation constant throughout the slice.

For purpose of this article, we consider Parallel Beam CT, and not Fan Beam CT. This is shown in the following figure, taken Chapter Three of Avinash Kak and Malcolm Slaney, Principles of Computerized Tomographic Imaging, available for free download here. As can be seen from this figure, the projections for the two different angles θ₁ and θ₂ are different. We don't present the equations, but the abovementioned book contains all the necessary math.

The projection function being a set of intensities, these intensities are depicted as a function of angle of incidence; the resulting figure is called a Sinogram. In a Sinogram, typically, the projection intensities corresponding to an angle of 0 degrees is shown as the left-most vertical line of the image. The projection intensities corresponding to an angle of 1 degree is shown as a vertical line one pixel to the right of the 0 degree line. And so on, till 180 lines are accumulated one by the side of another, corresponding to a range of 0 to 180 degrees. Some people prefer to show all the lines till 360 degrees, but for Parallel Beam CT, the projection intensities for an angle of θ are a mirror image of those for an angle of 180 + θ. For purposes of this article, we use 180 projections. Therefore the width of the Sinogram image is always 180 in our case.

If the input slice is a rectangular, it becomes convenient to convert that to a square, by padding pixels to the sides as appropriate. This software program does such a padding of pixels for rectangular images, to convert them to square.

This software application allows you to load an image, and computes the Sinogram of the image; and allows you to save the Sinogram image. Additionally, it allows you to invert the Sinogram too.

3. Code for Generating the Sinogram

The code for generating the sinogram is organized into three classes as described below:

1. SquarePaddedImage.cs: This class has the responsibility of converting a rectangular image to a square image, by appropriately padding dark coloured pixels to the image. If the original image has its width greater than its height, the dark pixels are added at the top and bottom of the image. On the other hand, if the original image has its width lesser than its height, the dark pixels are added at the left and right sides. Logic for adding the pixels on top and bottom, or left and right, is embedded inside the function ComputePaddedImage(). As an example, the logic for padding pixels at the top and bottom is given below:

   if (originalWidth > originalHeight)
   {
   	// Width is greater than height; so we pad background pixels at the top and bottom
   	// of the image, so that the original image comes in the middle of the square block
   	PaddedWidth = originalWidth;
   	PaddedHeight = originalWidth;
   	InitializePaddedPixels();
   
   	int diff = (PaddedHeight - originalHeight) / 2;
   
   	for (int i = diff * PaddedWidth; i < diff * PaddedWidth + originalWidth * originalHeight; ++i)
   	{
   		Pixels8PaddedRed[i] = pixels8OriginalRed[i - diff * PaddedWidth];
   		Pixels8PaddedGreen[i] = pixels8OriginalGreen[i - diff * PaddedWidth];
   		Pixels8PaddedBlue[i] = pixels8OriginalBlue[i - diff * PaddedWidth];
   	}
   }
   

   Additionally this class also reads in the original image, and also has the capability to create a bitmap from the padded image.
2. ImageRotation.cs: This class has the responsibility of rotating a square image by a given angle. For example, if the rotation angle is specified as 3 degrees, this class can perform such a rotation, and give the resulting pixel array as an output. Additionally, it can also output back the resulting bitmap.

    

   Rotation of a square image is achieved by means of the rotation formula, followed by the bilinear interpolation. The code for rotation along with bilinear interpolation is achieved by means of code adapted from something taken from the Internet.
3. Sinogram.cs: This is the most important class in this program, and does the function of computing the Sinogram image. First, the input image is padded to make it a square image. Then, in a loop, the image is rotated by an angle, and the pixel values are accumulated into a sum variables. These sum values are all placed appropriately into a sinogramValues buffer, which is then used to generate the Sinogram image. The relevant code extract to compute the sinogramValues buffer is given below:

   // Populate the sinogram buffer
   for (int k = 0; k < numberOfAngles; ++k)
   {
   	angleDegrees = -k;
   	// Just watch the console for large images. 
   	// It should go on until 180.
   	Console.Write(k + " "); 
   	ir.RotateImage(angleDegrees);
   	pixels8RotatedRed = ir.Pixels8RotatedRed;
   	pixels8RotatedGreen = ir.Pixels8RotatedGreen;
   	pixels8RotatedBlue = ir.Pixels8RotatedBlue;
   
   	for (int i = 0; i < SquareWidth; ++i)
   	{
   		sum = 0;
   		index1 = i * numberOfAngles + k;
   		for (int j = 0; j < SquareHeight; ++j)
   		{
   			index2 = j * SquareWidth + i;
   			red = pixels8RotatedRed[index2];
   			green = pixels8RotatedGreen[index2];
   			blue = pixels8RotatedBlue[index2];
   			gray = red * 0.3 + green * 0.59 + blue * 0.11;
   			gray /= maxsum;
   			sum += gray;
   		}
   		sinogramValues[index1] = Math.Exp(-sum);
   	} 
   }  

4. These three files, along with the Main Window XAML and CS files are packaged into a Visual Studio 2010 Solution.

4. Validation

For purposes of validation, I took the support of my students, who are quite familiar with the MATLAB software and have access to it. We have taken a set of representational images, and have run them through our program, and also through MATLAB, to check whether two aspects match well:

• Shapes of the sinograms, and
• Intensities of the sinograms.

The MATLAB code for creating a Sinogram is:

I1 = imread('Fig1.png');
I = rgb2gray(I1);
theta = 0:180;
[R,xp] = radon(I,theta);
imagesc(theta,xp,R);
colormap(gray);

A tabulation of the results is given below:

Original Image	Sinogram from this Program	Sinogram from MATLAB

The last of these is the well known Shepp-Logan phantom.

It is seen that the shapes of the sinograms match well with those obtained from MATLAB, for this set of representational images. This is barring the fact that the scales of the two sinograms are different, as it appears that MATLAB scales the images with its specific geometric scale factor.

As regards intensity, you may notice a discernible difference in intensity between the two sets of sinograms, for some of our representational images. This may be a result of differences in intensity-scaling between our program and MATLAB. In our case, we take the maximum intensity of the sinogram as mapping to 255, and the minimum intensity of the sinogram as mapping to 0. We add and subtract a tiny epsilon (0.00001) value from the maximum and minimum values respectively, just to make sure that the maximum pixel value does not exceed 255 and the minimum pixel value does not go below 0. This intensity scaling may be somewhat different in MATLAB, and this is explainable as the cause of differences in intensity.

5. Points of Interest

In this program, we have not focussed much on performance of the program, and you may find this computation a little slow. Our focus has mainly been on functionality. While running this program in Debug mode, you can see a counter of the angle values incrementing from 0 to 180 in steps of 1, in the debug output window.

Best results are obtained if the region of interest for which the Sinogram is to be computed is within the inscribed circle corresponding to the square image. There are different ways of handling the corner regions lying outside the inscribed circle, and we prefer not to do any special handling here.

The image of the Sinogram in the GUI is shown in Fill mode, so that it not seen as squished. This is irrespective of the fact that its width is always 180 pixels.

By changing the value of numberOfAngles in the file Sinogram.cs, the Sinogram can be computed for 360 degrees instead of 180 degrees.

We will focus on the Reconstruction aspect in a subsequent article.

6. Closure

In this article, a method to compute the Sinogram given an image was presented. The first step was to appropriately pad the image so as to convert any rectangular image to a square image. The Sinogram was created in the next two steps (i) To rotate the square image by a certain angle, and (ii) To compute the sum of the pixel values corresponding to the rotated image; these were performed in a loop. The resulting Sinogram was validated by comparing with MATLAB.

It is hoped that this article has helped in understanding of how a Sinogram is generated from a slice. Please experiment with generating Sinograms not only for anatomical images, but also any kind of general images and share your experiences through the Comments section below. If you feel anything amiss in the computation of the Sinogram, do write back. Even otherwise, do provide your valuable comments.

History

• Version 1.0: 19 July 2016.

Acknowledgements

I would like to thank my students Atheeth P, Jagruthi A, Kishore KS, Lakshmi SD, Mamatha K, Savitha VS, Sutapa B, for patiently listening to my lectures, and for helping me by running our representational images on MATLAB, and helping me to compare the results.