Workshop 05

Workshop 05 - nD Images

← Home →

Implementing a CT slice viewer

The Computed Tomography (CT) DICOM files contain both data => Hounsfield Units (HU) matrices, and metadata => information on patients, CT scanner used, etc. The HU scale represents a linear transformation of the X-Rays attenuation coefficient μ. Thus, distilled water at standard pressure and temperature (STP) is defined as 0 HU, while the radiodensity of air at STP is defined as −1000 HU.

Substance	HU
Air	-1000
Fat	−120 to −90
Bone	300 to 1900
Blod	13 to 50
Water	0
Lung	-700 to -600
Kidney	20 to 45
Liver	60 ± 6
Muscle	35 to 55
White matter	20 to 30
Gray matter	37 to 45
Windowpane glass	500
Limestone (Calcar)	2800
Cooper	14000

One could use the previously created Matrix class for loading a CT slice stored as a csv file: CT Slice csv example file
For viewing CT slices, one needs to transform it into a grayscale Bitmap image. The first step is applying a linear transformation, as bitmaps are limited to 256 values.

		
 public class LinearTransformation
 {
     public double A { get; private set; }
     public double B { get; private set; }

     public LinearTransformation(double x1, double y1, double x2, double y2)
     {
         A = (y2 - y1) / (x2 - x1);
         B = y1 - A * x1;
     }

     public double F(double x)
     {
         return A * x + B;
     }
 }

The next step is generating a Bitmap (see the below ToBitmap extension):

	
public static Bitmap NewBitmap(int width, int height, Color background)
{
	var bitmap = new Bitmap(width, height, PixelFormat.Format24bppRgb);

	using (var graph = Graphics.FromImage(bitmap))
	using (var brush = new SolidBrush(background))
	{
		var rectangle = new Rectangle(0, 0, width, height);
		graph.FillRectangle(brush, rectangle);
	}
	return bitmap;
}

		
public static Bitmap ToBitmap(this Matrix image, int minValue, int maxValue)
{
   var bitmap = NewBitmap(image.Cols, image.Rows, Color.White);

   var transform = new LinearTransformation(minValue, 0, maxValue, 255);

   for (int x = 0; x < image.Cols; x++)
   {
	   for (int y = 0; y < image.Rows; y++)
	   {
		   var value = image[y, x];
		   var imgValue = (int)transform.F(value);

		   if (imgValue < 0) imgValue = 0;
		   if (imgValue > 255) imgValue = 255;

		   var color = Color.FromArgb(imgValue, imgValue, imgValue);

		   bitmap.SetPixel(x, y, color);
	   }
   }

   return bitmap;
}

For viewing the image, one could now assign the bitmap object to the Image property of a PictureBox Win Forms control:

		
var repo = new MatrixDataAccess();
var matrix = repo.Read("C:\\...\\slice.csv");
var bitmap = matrix.ToBitmap(0, 100);
pictureBox1.Width = matrix.Cols;
pictureBox1.Height = matrix.Rows;
pictureBox1.Image = bitmap;

Try generating the bitmap for various HU windows.
It is important to note the difficulties in viewing/understanding the whole information contained in one HU matrix.
For storing 3D data, CTs are composed of various slice series (native, arterial contrast, venous contrast).

Multi-dimensional images

For storing a set of multi-dimensional points (e.g. Hyperspectral Images) one could use a csv file with a different structure: one column for each nD coordinate. In a Fuzzy approach, an additional column could be considered for storing the Degree of Membership of each nD point.

Find below a nD vector implementation:

		
public class NVector
{
	private double[] _components;

	public int Dimension
	{
		get
		{
			return _components.Length;
		}
	}

	public NVector(int dimension)
	{
		_components = new double[dimension];
	}

	public NVector(double[] components)
	{
		_components = components;
	}

	public void SetAllValues(double value)
	{
		for (int index = 0; index < _components.Length; index++)
		{
			_components[index] = value;
		}
	}

	public double this[int index]
	{
		get
		{
			return _components[index];
		}
		set
		{
			_components[index] = value;
		}
	}

	public NVector Clone()
	{
		var vector = new NVector(Dimension);
		var index = 0;
		foreach (var component in _components)
		{
			vector[index] = component;
			index++;
		}
		return vector;
	}

	public string ToCsv()
	{
		return string.Join(Csv.CSV_VALUES_SEPARATOR, _components);
	}

	public override string ToString()
	{
		return ToCsv();
	}

	public override bool Equals(object obj)
	{
		if (obj == null) return false;
		if (!(obj is NVector)) return false;
		var vector = (NVector)obj;
		for (int index = 0; index < Dimension; index++)
		{
			if (_components[index] != vector[index])
			{
				return false;
			}
		}
		return true;
	}

	public override int GetHashCode()
	{
		var hashString = String.Join("_", _components);
		return hashString.GetHashCode();
	}

	public double Module
	{
		get
		{
			var module = 0.0;
			for (int index = 0; index < Dimension; index++)
			{
				module += _components[index] * _components[index];
			}
			module = Math.Sqrt(module);
			return module;
		}
	}

	public static NVector operator +(NVector p, NVector q)
	{
		var vector = new NVector(p.Dimension);

		for (int index = 0; index < p.Dimension; index++)
		{
			vector[index] = p[index] + q[index];
		}

		return vector;
	}
	public static NVector operator -(NVector p, NVector q)
	{
		var vector = new NVector(p.Dimension);

		for (int index = 0; index < p.Dimension; index++)
		{
			vector[index] = p[index] - q[index];
		}

		return vector;
	}

	public static NVector operator *(double scalar, NVector p)
	{
		var vector = new NVector(p.Dimension);

		for (int index = 0; index < p.Dimension; index++)
		{
			vector[index] = scalar * p[index];
		}

		return vector;
	}

	public static NVector operator *(NVector p, double scalar)
	{
		return scalar * p;
	}

	public static NVector operator /(NVector p, double scalar)
	{
		var vector = new NVector(p.Dimension);

		for (int index = 0; index < p.Dimension; index++)
		{
			vector[index] = p[index] / scalar;
		}

		return vector;
	}
}

Thus, an nD image could be represented as a generic list of NVector elements.
Consider also the corresponding data access implementation:

		
public interface INVectorDataAccess
{
	void Write(IEnumerable<NVector> image, string file);
}

public class NVectorDataAccess : INVectorDataAccess
{
	public static StringBuilder ToCsv(IEnumerable<NVector> image)
	{
	  var sb = new StringBuilder();
	  foreach (var nVector in image)
	  {
		  sb.AppendLine(nVector.ToCsv()); 
	  }
	  return sb;
	}

	public void Write(IEnumerable<NVector> image, string file)
	{
	  var sb = ToCsv(image);
	  File.WriteAllText(file, sb.ToString());
	}
}

For viewing a 2D projection of the nD image, one could consider the following method:

	
public static Bitmap ToBitmap(this List<NVector> image, int width, int height, LinearTransformation transform, int xi = 0, int xj = 1)
{
	var bitmap = NewBitmap(width, height, Color.White);

	foreach (var nD in image)
	{
		var color = Color.FromArgb(0, 0, 0);

		bitmap.SetPixel((int)transform.F(nD[xi]), (int)transform.F(nD[xj]), color);
	}

	return bitmap;
}

nD Exercises:

Generate the coordinates of a Circle:

		
public static List<NVector> Circle(double radius, double step = 0.1)
{
    var img = new List<NVector>();
    for(var x = -radius; x <= radius; x += step)
    {
        var delta = radius * radius - x * x;

        var nV1 = new NVector(2);
        nV1[0] = x;
        nV1[1] = Math.Sqrt(delta);
        img.Add(nV1);
        
        var nV2 = new NVector(2);
        nV2[0] = x;
        nV2[1] = -nV1[1];
        img.Add(nV2);
    }
    return img;
}

Generate the coordinates of a Sphere and test the result by viewing various 2D projections:

		
public static List<NVector> Sphere(double radius, double step = 0.1)
{
	var img = new List<NVector>();
	for (var x = -radius; x <= radius; x += step)
	{
		for (var y = -radius; y <= radius; y += step)
		{
			var delta = radius * radius - (x * x + y * y);
			if (delta < 0) continue;

			var nV1 = new NVector(3);
			nV1[0] = x;
			nV1[1] = y;
			nV1[2] = Math.Sqrt(delta);
			img.Add(nV1);

			var nV2 = new NVector(3);
			nV2[0] = x;
			nV2[1] = y;
			nV2[2] = -nV1[2];
			img.Add(nV2);
		}
	}
	return img;
}

Generate the coordinates of a 4D HyperSphere and test the result by viewing various 2D projections:

		
public static List<NVector> HyperSphere(double radius, double step = 0.1)
{
	var img = new List<NVector>();
	for (var x = -radius; x <= radius; x += step)
	{
		for (var y = -radius; y <= radius; y += step)
		{
			for (var z = -radius; z <= radius; z += step)
			{
				var delta = radius * radius - (x * x + y * y + z * z);
				if (delta < 0) continue;

				var nV1 = new NVector(4);
				nV1[0] = x;
				nV1[1] = y;
				nV1[2] = z;
				nV1[3] = Math.Sqrt(delta);
				img.Add(nV1);

				var nV2 = new NVector(4);
				nV2[0] = x;
				nV2[1] = y;
				nV2[2] = z;
				nV2[3] = -nV1[3];
				img.Add(nV2);
			}
		}
	}
	return img;
}

References:

Application for estimating the fuzzy fractal dimension of hyperspectral satellite ocean color images
Optimization technique for increasing resolution in computed tomography imaging
.Net WPF module for RGB 3D reconstruction of medical three-channel images

← Home →