package edu.stanford.rsl.tutorial.truncation;
import edu.stanford.rsl.conrad.data.numeric.Grid1D;
import edu.stanford.rsl.conrad.data.numeric.Grid1DComplex;
import edu.stanford.rsl.conrad.data.numeric.Grid2D;
import edu.stanford.rsl.tutorial.parallel.ParallelBackprojector2D;
import edu.stanford.rsl.tutorial.parallel.ParallelProjector2D;
import edu.stanford.rsl.tutorial.phantoms.SimpleGridsForTruncationCorrection;
import ij.ImageJ;
/** Jennifer Maier
*
* This class creates the sinogram of a simple grid, which is afterwards truncated.
* The truncated sinogram is then corrected by extrapolating with a polynomial for each row.
*
*/
public class PolynomialTruncationCorrectionExample {
public PolynomialTruncationCorrectionExample(){}
public static void main(String[] args) {
new ImageJ();
PolynomialTruncationCorrectionExample T = new PolynomialTruncationCorrectionExample();
// create Phantom
int imageSize = 200;
double maxS = 200.0;
double deltaS = 1.0;
int numTheta = 360;
int truncationSize = 40;
// third parameter:
// 0: circle, 1: ellipse, 2: sqaure
SimpleGridsForTruncationCorrection original = new SimpleGridsForTruncationCorrection(imageSize, imageSize, 2);
original.show("original image");
// get and show Sinogram of the cylinder
ParallelProjector2D projector2d = new ParallelProjector2D(Math.PI, Math.PI/(numTheta-1), maxS, deltaS);
Grid2D sinogram = projector2d.projectRayDrivenCL(original);
sinogram.show("original sinogram");
// apply ramlak filter to sinogram
Grid2D origSino = new Grid2D(sinogram);
T.ramlak(origSino);
// backproject original
ParallelBackprojector2D backprojector = new ParallelBackprojector2D(imageSize, imageSize, 1, 1);
Grid2D backproj = backprojector.backprojectRayDriven(origSino);
backproj.show("original backprojected (ramlak filtered)");
// truncate Sinogram
Grid2D truncatedSinogram = T.getTruncatedSinogram(sinogram, truncationSize);
truncatedSinogram.setSpacing(sinogram.getSpacing());
truncatedSinogram.show("truncated sinogram");
Grid2D truncatedSinogramCopy = new Grid2D(truncatedSinogram);
T.ramlak(truncatedSinogramCopy);
Grid2D truncatedBackprojected = backprojector.backprojectRayDriven(truncatedSinogramCopy);
truncatedBackprojected.show("truncated backprojected (ramlak filtered)");
// correct truncation using polynomials
Grid2D correctedSinogram = T.polynomialTruncationCorrection(truncatedSinogram);
correctedSinogram.setSpacing(sinogram.getSpacing());
// add black regions on each side to ensure correct reconstruction
Grid2D broaderCorrectedSinogram = T.addBlackOnSides(correctedSinogram, sinogram);
broaderCorrectedSinogram.setSpacing(sinogram.getSpacing());
broaderCorrectedSinogram.show("corrected sinogram");
// apply ramlak filter to corrected sinogram
T.ramlak(broaderCorrectedSinogram);
// backproject corrected sinogram
ParallelBackprojector2D backprojector2 = new ParallelBackprojector2D(imageSize, imageSize, 1, 1);
Grid2D correctedBackproj = backprojector2.backprojectRayDriven(broaderCorrectedSinogram);
correctedBackproj.show("corrected backprojected (ramlak filtered)");
}
/**
* correct truncation with a polynomial in each row
*
* @param truncatedSinogram truncated sinogram that has to be corrected
*/
private Grid2D polynomialTruncationCorrection(Grid2D truncatedSinogram) {
// get L for each angle
int[] l = getL(truncatedSinogram);
int maxL = max(l);
// get derivative of truncated Sinogram
Grid2D firstDerivativeTruncatedSinogram = derive(truncatedSinogram);
// get width of the truncated sinogram
int truncatedSinogramWidth = getSinogramWidth(truncatedSinogram);
int R = truncatedSinogramWidth/2;
Grid2D correctedSinogram = new Grid2D(2*maxL, truncatedSinogram.getHeight());
int posDerivative = (int) Math.floor(0.5f*(truncatedSinogram.getWidth()-truncatedSinogramWidth));
// correct Sinogram on right and left side
for (int sinoRow = 0; sinoRow < correctedSinogram.getHeight(); sinoRow++) {
//get L, p', p for computation of a, b, c
double L = l[sinoRow];
double pPrime_left = firstDerivativeTruncatedSinogram.getAtIndex(posDerivative, sinoRow);
double p_left = truncatedSinogram.getAtIndex(posDerivative, sinoRow);
double pPrime_right = firstDerivativeTruncatedSinogram.getAtIndex(firstDerivativeTruncatedSinogram.getWidth()-posDerivative-1, sinoRow);
double p_right = truncatedSinogram.getAtIndex(truncatedSinogram.getWidth()-posDerivative-1, sinoRow);
double a_right = (p_right*p_right - 2*p_right*pPrime_right*(R-L))/(3*R*R-L*L-2*R*L);
double b_right = 2*p_right*pPrime_right + 2*a_right*R;
double c_right = -a_right*L*L - b_right*L;
double a_left = (p_left*p_left - 2*p_left*pPrime_left*((-R)-(-L)))/(3*(-R)*(-R)-(-L)*(-L)-2*(-R)*(-L));
double b_left = 2*p_left*pPrime_left + 2*a_left*(-R);
double c_left = -a_left*(-L)*(-L) - b_left*(-L);
// add polynomial on left and right side
for (int col = 0; col < L-truncatedSinogramWidth/2; col++) {
correctedSinogram.setAtIndex(maxL-R-col, sinoRow, (float) Math.sqrt(a_left*(-R-col)*(-R-col)+b_left*(-R-col)+c_left));
correctedSinogram.setAtIndex(maxL+R+col, sinoRow, (float) Math.sqrt(a_right*(R+col)*(R+col)+b_right*(R+col)+c_right));
}
// take original values in between
for (int col = 0; col < truncatedSinogramWidth; col++) {
correctedSinogram.setAtIndex(maxL-truncatedSinogramWidth/2 + col, sinoRow, truncatedSinogram.getAtIndex(posDerivative+col, sinoRow));
}
}
return correctedSinogram;
}
/**
* truncate sinogram: delete values at left and right side
*
* @param sinogram sinogram that will be truncated
* @param truncationWidth number of pixels that are deleted in each row
* starting from left and right image boarder
*/
private Grid2D getTruncatedSinogram(Grid2D sinogram, int truncationWidth) {
Grid2D truncatedSinogram = new Grid2D(sinogram);
for (int row = 0; row < sinogram.getHeight(); row++) {
for (int col = 0; col < truncationWidth; col++) {
truncatedSinogram.setAtIndex(col, row, 0.0f);
truncatedSinogram.setAtIndex(truncatedSinogram.getWidth()-1-col, row, 0.0f);
}
}
return truncatedSinogram;
}
/**
* add black on the sides of the corrected sinogram (for reconstruction)
*
* @param correctedSinogram sinogram to which black will be added
* @param originalSinogram the resulting sinogram will have the same width as this Grid2D
*/
private Grid2D addBlackOnSides(Grid2D correctedSinogram, Grid2D originalSinogram) {
Grid2D result = new Grid2D(originalSinogram.getWidth(), correctedSinogram.getHeight());
for (int row = 0; row < result.getHeight(); row ++) {
for (int col = 0; col < result.getWidth(); col++) {
result.setAtIndex(col, row, 0.0f);
}
}
int begin = (result.getWidth()-correctedSinogram.getWidth())/2;
for (int row = 0; row < result.getHeight(); row++) {
for (int col = 0; col < correctedSinogram.getWidth(); col++) {
result.setAtIndex(col+begin, row, correctedSinogram.getAtIndex(col, row));
}
}
return result;
}
/**
* get actual width of sinogram (from first to last non-zero value)
*
* @param truncatedSinogram
*/
private int getSinogramWidth(Grid2D truncatedSinogram) {
int start = truncatedSinogram.getWidth();
int end = 0;
for (int row = 0; row < truncatedSinogram.getHeight(); row++) {
boolean rowStartFound = false;
boolean rowEndFound = false;
for (int col = 0; col < truncatedSinogram.getWidth(); col++) {
if (rowStartFound == false
&& truncatedSinogram.getAtIndex(col, row) != 0.0f
&& !Double.isNaN(truncatedSinogram.getAtIndex(col, row))
) {
if (col < start) {
start = col;
}
rowStartFound = true;
}
if (rowEndFound == false
&& truncatedSinogram.getAtIndex(truncatedSinogram.getWidth()-1-col, row) != 0.0f
&& !Double.isNaN(truncatedSinogram.getAtIndex(truncatedSinogram.getWidth()-1-col, row))
){
if (truncatedSinogram.getWidth()-col > end) {
end = truncatedSinogram.getWidth()-col;
}
rowEndFound = true;
}
if (rowStartFound && rowEndFound) break;
}
}
return end-start;
}
/**
* get derivative along rows
*
* @param sinogram
*/
private Grid2D derive(Grid2D sinogram) {
Grid2D firstDerivative = new Grid2D(sinogram.getWidth()-1, sinogram.getHeight());
for (int row = 0; row < sinogram.getHeight(); row++) {
for (int col = 0; col < firstDerivative.getWidth(); col++) {
firstDerivative.setAtIndex(col, row, sinogram.getAtIndex(col+1, row) - sinogram.getAtIndex(col, row));
}
}
return firstDerivative;
}
/**
* find maximum of an positive integer array
*
* @param l
*/
private int max(int[] l) {
int curr_max = 0;
for (int i = 0; i < l.length; i++) {
if (l[i] > curr_max) {
curr_max = l[i];
}
}
return curr_max;
}
/**
* get width L of polynomial for each row in the polynomial truncation correction
*
* @param sinogram
*/
private int[] getL(Grid2D sinogram) {
double deltaTheta = 180.0f/sinogram.getHeight();
int[] L = new int[sinogram.getHeight()];
double mu_average = 1.0f;
for (int row = 0; row < sinogram.getHeight(); row++) {
double theta = row * deltaTheta;
double thetaRad = theta * (Math.PI/180);
double max = 0.0f;
for (int rowPrime = (int) (row - 90/deltaTheta); rowPrime < (int) (row + 90/deltaTheta); rowPrime++) {
boolean changedThetaPrimeMin = false;
boolean changedThetaPrimePlus = false;
if (rowPrime >= sinogram.getHeight()) {
rowPrime -= sinogram.getHeight();
changedThetaPrimeMin = true;
} else if (rowPrime < 0) {
rowPrime += sinogram.getHeight();
changedThetaPrimePlus = true;
}
double thetaPrime = rowPrime * deltaTheta;
double thetaPrimeRad = thetaPrime * (Math.PI/180);
for (int s = 0; s < sinogram.getWidth(); s++) {
double curr_L = sinogram.getAtIndex(s, rowPrime) * Math.sin(thetaRad - thetaPrimeRad) / mu_average;
if (Math.abs(curr_L) > max) {
max = Math.abs(curr_L);
}
}
if (changedThetaPrimeMin) {
rowPrime += sinogram.getHeight();
} else if (changedThetaPrimePlus) {
rowPrime -= sinogram.getHeight();
}
}
L[row] = (int) max/2;
}
return L;
}
/**
* RamLak filtering of the input sinogram
*
* @param sinogram
*/
public void ramlak(Grid2D sinogram) {
Grid1DComplex[] sinogramTransform = new Grid1DComplex[(sinogram.getHeight())];
Grid1D[] filteredSinogram = new Grid1D[(sinogram.getHeight())];
Grid1DComplex ramlacFilter = new Grid1DComplex(sinogram.getWidth());
int middle = (int) ramlacFilter.getSize()[0]/2;
ramlacFilter.setRealAtIndex(middle, 1.0f/4);
for (int i = middle + 1; i < ramlacFilter.getSize()[0]; i+=2) {
ramlacFilter.setRealAtIndex(i, (float) (-1.0f/((i-middle)*(i-middle)*Math.PI*Math.PI)));
ramlacFilter.setImagAtIndex(i, 0);
}
for (int i = middle - 1; i >= 0; i-=2) {
ramlacFilter.setRealAtIndex(i, (float) (-1.0f/((i-middle)*(i-middle)*Math.PI*Math.PI)));
ramlacFilter.setImagAtIndex(i, 0);
}
ramlacFilter.transformForward();
for (int row = 0; row < sinogram.getHeight(); row++) {
sinogramTransform[row] = new Grid1DComplex(sinogram.getSubGrid(row));
sinogramTransform[row].transformForward();
for (int i = 0; i < ramlacFilter.getSize()[0]; i++) {
sinogramTransform[row].multiplyAtIndex(i, ramlacFilter.getAtIndex(i));
}
sinogramTransform[row].transformInverse();
filteredSinogram[row] = sinogramTransform[row].getRealSubGrid(0,sinogram.getWidth());
}
for (int row = 0; row < sinogram.getHeight(); row++) {
for (int col = 0; col < sinogram.getWidth(); col++) {
sinogram.setAtIndex(col, row, filteredSinogram[row].getAtIndex(col));
}
}
}
}
/*
* Copyright (C) 2015 Jennifer Maier, jennifer.maier@fau.de
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/