conebeamreconstruction.java

ImageJ Conebeam Filtered Back Projection Plug-in. To use this plug-in, a stack of projection images

import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.gui.NewImage;
import ij.plugin.filter.PlugInFilter;
import ij.process.ImageProcessor;

public class ConeBeamReconstruction implements PlugInFilter {

	ImagePlus imp;
	protected ImageStack projectionStack;

	public int setup(String arg, ImagePlus imp) {
		if (arg.equals("about")) {
			showAbout();
			return DONE;
		}

		this.imp = imp;
		this.projectionStack = imp.getStack();
		return DOES_32 + STACK_REQUIRED + NO_CHANGES;
	}

	private void showAbout() {
		IJ.showMessage("About Cone Beam filtered Back Projection Reconstruction",
				"Welcome");
	}

	public void run(ImageProcessor ip) {

		int projectionWidth = this.projectionStack.getWidth();
		int projectionHeight = this.projectionStack.getHeight();
		int projectionNum = this.projectionStack.getSize();

		int reconWidth = projectionWidth;
		int reconHeight = projectionWidth;
		int reconNum = projectionHeight;

		ImagePlus reconImage = NewImage.createFloatImage(
				"Reconstruction Result", reconWidth, reconHeight, reconNum,
				NewImage.FILL_BLACK);

		float[][] reconPixels = new float[reconNum][];

		for (int stack = 1; stack <= reconNum; stack++) {
			reconPixels[stack - 1] = (float[]) reconImage.getStack()
					.getProcessor(stack).getPixels();
		}

		float[][] projectionPixels = new float[projectionNum][];
		for (int stack = 1; stack <= projectionNum; stack++) {
			projectionPixels[stack - 1] = (float[]) projectionStack
					.getProcessor(stack).getPixels();
		}

		float[] sinValue = new float[projectionNum];
		float[] cosValue = new float[projectionNum];

		this.getSinCosTable(sinValue, cosValue);

		int reconHalfWidth = reconWidth / 2;
		int reconHalfHeight = reconHeight / 2;
		int reconHalfNum = reconNum / 2;

		int xStart = -(int) Math.floor(reconWidth / 2) + 1;
		int yStart = -(int) Math.floor(reconHeight / 2) + 1;
		int zStart = -(int) Math.floor(reconNum / 2) + 1;

		int xIndex, yIndex, zIndex;
		float t, s, p, c;
		float dso = (float) IJ.getNumber("DSO", 1168);// DSO

		float mul = 0;
		float percentStep = 1 / (float) projectionNum;
		float percentWork = 0;

		int shiftLeft = (int) IJ.getNumber("ShiftLeft", -8);
		// int sliceNumber = (int) IJ.getNumber("Slice Number", reconHalfNum);
		long startTime = System.currentTimeMillis();
		for (int reconIndex = 1, z = zStart; reconIndex <= reconNum; reconIndex++, z++) {
			percentWork = 0;


			for (int ang = 0; ang < projectionNum; ang++) {
				for (int x = xStart; x <= reconHalfWidth; x++) {
					for (int y = yStart; y <= reconHalfHeight; y++) {

						s = x * cosValue[ang] + y * sinValue[ang];
						t = y * cosValue[ang] - x * sinValue[ang];

						p = (dso * t) / (dso - s);
						p = Math.round(p + reconHalfWidth) - shiftLeft;// Shifting


						c = (dso * z) / (dso - s);
						c = Math.round(c + reconHalfNum);



						if (p > 0 && c > 0 && p < reconWidth && c < reconNum) {

							//							
							xIndex = x + reconHalfWidth - 1;
							yIndex = y + reconHalfHeight - 1;
							zIndex = reconIndex - 1;
							mul = (dso * dso) / ((dso - s) * (dso - s));
							try {

								reconPixels[zIndex][yIndex * reconWidth
										+ xIndex] += (mul * projectionPixels[ang][(int) (c
										* projectionWidth + p)]) / 2;
							} catch (Exception e) {
								// IJ.showMessage("information" + "," + zIndex
								// + "," + yIndex + "," + xIndex + "," + c
								// + "," + p + "," + yIndex * reconWidth
								// + xIndex + ","
								// + (int) (c * projectionWidth + p));
							}

						}
					}

				}
				percentWork += percentStep;
				IJ.showStatus("Woriking:" + reconIndex + "/" + reconNum);
				IJ.showProgress(percentWork);
			}

		}
		reconImage.setSlice(reconHalfNum);
		reconImage.updateAndDraw();
		reconImage.show();
		long endTime = System.currentTimeMillis();
		IJ.showMessage("Consuming Time in Millis = " + (endTime - startTime));
		
		IJ.runPlugIn("ij.plugin.frame.ContrastAdjuster", "wl");

	}

	public void getSinCosTable(float[] sinValue, float[] cosValue) {
		if (sinValue.length != cosValue.length) {
			System.out.println("Sin and Cos table size not the same");
		}
		int AngleNum = sinValue.length;

		float angleStep = (float) (2 * Math.PI / AngleNum);

		float angle = 0;
		for (int i = 0; i < AngleNum; i++) {
			sinValue[i] = (float) Math.sin(angle);
			cosValue[i] = (float) Math.cos(angle);
			angle += angleStep;
		}
	}

	public static void main(String[] args) {
		int projectionNum = 306;
		float[] sinValue = new float[projectionNum];
		float[] cosValue = new float[projectionNum];

		new ConeBeamReconstruction().getSinCosTable(sinValue, cosValue);

		for (int i = 0; i < cosValue.length; i++) {
			System.out.println("i :" + i + ", Sin :" + sinValue[i] + ", Cos :"
					+ cosValue[i]);
		}
	}
}