diff --git a/.classpath b/.classpath index 31f9bd0c..1e9e77a1 100644 --- a/.classpath +++ b/.classpath @@ -43,7 +43,6 @@ - @@ -69,5 +68,6 @@ + diff --git a/src/edu/stanford/rsl/tutorial/dmip/Proton.png b/data/mipda/Proton.png similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/Proton.png rename to data/mipda/Proton.png diff --git a/src/edu/stanford/rsl/tutorial/dmip/Sinogram0.tif b/data/mipda/Sinogram0.tif similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/Sinogram0.tif rename to data/mipda/Sinogram0.tif diff --git a/src/edu/stanford/rsl/tutorial/dmip/Sinogram1.tif b/data/mipda/Sinogram1.tif similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/Sinogram1.tif rename to data/mipda/Sinogram1.tif diff --git a/src/edu/stanford/rsl/tutorial/dmip/Sinogram2.tif b/data/mipda/Sinogram2.tif similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/Sinogram2.tif rename to data/mipda/Sinogram2.tif diff --git a/src/edu/stanford/rsl/tutorial/dmip/T1.png b/data/mipda/T1.png similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/T1.png rename to data/mipda/T1.png diff --git a/src/edu/stanford/rsl/tutorial/dmip/frame32.jpg b/data/mipda/frame32.jpg similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/frame32.jpg rename to data/mipda/frame32.jpg diff --git a/src/edu/stanford/rsl/tutorial/dmip/frame90.jpg b/data/mipda/frame90.jpg similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/frame90.jpg rename to data/mipda/frame90.jpg diff --git a/src/edu/stanford/rsl/tutorial/dmip/mask.bmp b/data/mipda/mask.bmp similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/mask.bmp rename to data/mipda/mask.bmp diff --git a/src/edu/stanford/rsl/tutorial/dmip/mr12.dcm b/data/mipda/mr12.dcm similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/mr12.dcm rename to data/mipda/mr12.dcm diff --git a/data/mipda/mr_head_angio.jpg b/data/mipda/mr_head_angio.jpg new file mode 100644 index 00000000..453f5786 Binary files /dev/null and b/data/mipda/mr_head_angio.jpg differ diff --git a/data/mipda/mr_head_dorsal.jpg b/data/mipda/mr_head_dorsal.jpg new file mode 100644 index 00000000..a3c181e5 Binary files /dev/null and b/data/mipda/mr_head_dorsal.jpg differ diff --git a/data/mipda/output/info.txt b/data/mipda/output/info.txt new file mode 100644 index 00000000..61319ec6 --- /dev/null +++ b/data/mipda/output/info.txt @@ -0,0 +1,2 @@ +This folder is used for output data from MIPDA exercises. +Do not put anything else here because it might get overwritten! \ No newline at end of file diff --git a/src/edu/stanford/rsl/tutorial/dmip/testimg.bmp b/data/mipda/testimg.bmp similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/testimg.bmp rename to data/mipda/testimg.bmp diff --git a/src/edu/stanford/rsl/tutorial/dmip/undistorted.jpg b/data/mipda/undistorted.jpg similarity index 100% rename from src/edu/stanford/rsl/tutorial/dmip/undistorted.jpg rename to data/mipda/undistorted.jpg diff --git a/src/edu/stanford/rsl/tutorial/dmip/DMIP_FanBeamBackProjector2D.java b/src/edu/stanford/rsl/tutorial/dmip/DMIP_FanBeamBackProjector2D.java deleted file mode 100644 index 17780bc8..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/DMIP_FanBeamBackProjector2D.java +++ /dev/null @@ -1,329 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import edu.stanford.rsl.conrad.data.numeric.Grid1D; -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators; -import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; -import edu.stanford.rsl.conrad.filtering.redundancy.ParkerWeightingTool; -import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND; -import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import edu.stanford.rsl.tutorial.fan.FanBeamBackprojector2D; -import edu.stanford.rsl.tutorial.fan.redundancy.ParkerWeights; -import edu.stanford.rsl.tutorial.filters.RamLakKernel; -import edu.stanford.rsl.tutorial.filters.SheppLoganKernel; -import ij.IJ; -import ij.ImageJ; - -/** - * Exercise 7 of Diagnostic Medical Image Processing (DMIP) - * @author Marco Boegel - * - */ -public class DMIP_FanBeamBackProjector2D { - - //source to detector distance - private double focalLength; - //increment of rotation angle beta between two projections in [rad] - private double betaIncrement; - //number of pixels on the detector - private int detectorPixels; - //number of acquired projections - private int numProjs; - //spacing - mm per pixel - private double detectorSpacing; - //detector length in [mm] - private double detectorLength; - //max rotation angle beta in [rad] - private double maxBeta; - - public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN}; - - public DMIP_FanBeamBackProjector2D( ) - { - - } - - /** - * Initialize all relevant parameters for the reconstruction - * @param sino the sinogram - * @param focalLength source to detector distance - * @param maxRot maximum rotation angle in [rad] at which the sinogram was acquired - */ - public void setSinogramParams(Grid2D sino, double focalLength, double maxRot) - { - this.focalLength = focalLength; - this.numProjs = sino.getHeight(); - this.detectorPixels = sino.getWidth(); - this.detectorSpacing = sino.getSpacing()[0]; - this.detectorLength = detectorSpacing*detectorPixels; - - - double halfFanAngle = 0;//TODO - System.out.println("Half fan angle: " + halfFanAngle*180.0/Math.PI); - //TODO - this.betaIncrement = maxBeta /(double) numProjs; - System.out.println("Short-scan range: " + maxBeta*180/Math.PI); - } - - /** - * A pixel driven backprojection algorithm. Cosine, Redundancy and Ramp filters need to be applied separately beforehand - * @param sino the filtered sinogram - * @param recoSize the dimension of the output image - * @param spacing the spacing of the output image - * @return the reconstruction - */ - public Grid2D backproject(Grid2D sino, int[] recoSize, double[] spacing) - { - Grid2D result = new Grid2D(recoSize[0], recoSize[1]); - result.setSpacing(spacing[0], spacing[1]); - - for(int p = 0; p < numProjs; p++) - { - //First, compute the rotation angle beta and pre-compute cos(beta), sin(beta) - float beta = (float) (betaIncrement * p); - float cosBeta = (float) Math.cos(beta); - float sinBeta = (float) Math.sin(beta); - - //Compute direction and normal of the detector at the current rotation angle - final PointND detBorder = new PointND();//TODO - final SimpleVector dirDet = new SimpleVector();//TODO - final StraightLine detLine = new StraightLine(detBorder, dirDet); - - //Compute rotated source point - final PointND source = new PointND(focalLength * cosBeta, focalLength*sinBeta, 0.d); - - //pick current projection - Grid1D currProj = sino.getSubGrid(p); - - //pixel driven BP: iterate over all output image pixels - for(int x = 0; x < recoSize[0]; x++) - { - //transform the image pixel coordinates to world coordinates - float wx =0; //TODO - - for(int y = 0; y < recoSize[1]; y++) - { - float wy = 0;//TODO - - final PointND reconstructionPointWorld = new PointND(wx, wy, 0.d); - - //intersect the projection ray with the detector - //TODO - final PointND detPixel = new PointND();//TODO - - float valueTemp; - - if(detPixel != null) - { - //calculate position of the projected point on the 1D detector - final SimpleVector pt = SimpleOperators.subtract(detPixel.getAbstractVector(), detBorder.getAbstractVector()); - double length = pt.normL2(); - //deal with truncation - if(pt.getElement(0)*dirDet.getElement(0)+pt.getElement(1)*dirDet.getElement(1)<0) - { - length -= length; - } - double t = (length - 0.5d)/detectorSpacing; - - //check if projection of this world point hit the detector - if(currProj.getSize()[0] <= t+1 || t< 0) - continue; - - float value = InterpolationOperators.interpolateLinear(currProj, t); - - //Apply distance weighting - //see Fig 1a) exercise sheet - //TODO - //TODO - float dWeight = 0;//TODO - valueTemp = (float) (value / (dWeight*dWeight)); - } - else - { - //projected pixel lies outside the detector - valueTemp = 0.f; - } - - result.addAtIndex(x, y, valueTemp); - } - } - } - - //adjust scale. required because of shortscan - float normalizationFactor = (float) (numProjs / Math.PI); - NumericPointwiseOperators.divideBy(result, normalizationFactor); - - - return result; - } - - /** - * Cosine weighting for 2D sinograms - * @param sino the sinogram - * @return the cosine weighted sinogram - */ - public Grid2D applyCosineWeights(Grid2D sino) - { - Grid2D result = new Grid2D(sino); - - //Create 1D kernel (identical for each projection) - Grid1D cosineKernel = new Grid1D(detectorPixels); - for(int i=0; i Math.PI *2.d) { - continue; - } - - // implement the conditions as described in Parker's paper - if (beta <= 2 * (delta - gamma)) { - float val = 0; //TODO - if (Double.isNaN(val)){ - continue; - } - parker.setAtIndex(t, b , val); - - } else if (beta < Math.PI - 2.d * gamma) { - parker.setAtIndex(t, b , 1); - } - else if (beta <= (Math.PI + 2.d * delta) + 1e-12) { - float val = 0;//TODO - if (Double.isNaN(val)){ - continue; - } - parker.setAtIndex(t, b , val); - } - } - } - - // Correct for scaling due to varying angle - NumericPointwiseOperators.multiplyBy(parker, (float)( maxBeta / (Math.PI))); - parker.show(); - - //apply the parker weights to the sinogram - NumericPointwiseOperators.multiplyBy(result, parker); - - return result; - } - - - public static void main(String arg[]) - { - ImageJ ij = new ImageJ(); - - double focalLength = 600.f; - //maximum rotation angle in [rad] - double maxRot = Math.PI; - //choose the ramp filter (none, ramlak, shepplogan) - RampFilterType filter = RampFilterType.RAMLAK; - - DMIP_FanBeamBackProjector2D fbp = new DMIP_FanBeamBackProjector2D(); - - //Load and visualize the projection image data - String filename = "D:/02_lectures/DMIP/exercises/2014/6/Sinogram0.tif"; - Grid2D sino = ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); - sino.show("Sinogram"); - - //getSubGrid() only yields rows. -> Transpose so that each row is a projection - sino = sino.getGridOperator().transpose(sino); - sino.show(); - - - //Initialize parameters - double detectorSpacing = 0.5; - sino.setSpacing(detectorSpacing, 0); - fbp.setSinogramParams(sino, focalLength, maxRot); - sino.setSpacing(fbp.detectorSpacing, fbp.betaIncrement); - - //apply cosine weights - Grid2D cosineWeighted = fbp.applyCosineWeights(sino); - cosineWeighted.show("After Cosine Weighting"); - - //apply parker redundancy weights - Grid2D parkerWeighted = fbp.applyParkerWeights(cosineWeighted); - parkerWeighted.show("After Parker Weighting"); - - //apply ramp filter - switch(filter){ - case RAMLAK: - RamLakKernel ramLak = new RamLakKernel(fbp.detectorPixels, fbp.detectorSpacing); - for(int i = 0; i < fbp.numProjs; i++) - { - ramLak.applyToGrid(parkerWeighted.getSubGrid(i)); - } - parkerWeighted.show("After RamLak Filter"); - break; - case SHEPPLOGAN: - SheppLoganKernel sheppLogan = new SheppLoganKernel(fbp.detectorPixels, fbp.detectorSpacing); - for(int i = 0; i < fbp.numProjs; i++) - { - sheppLogan.applyToGrid(parkerWeighted.getSubGrid(i)); - } - parkerWeighted.show("After Shepp-Logan Filter"); - break; - case NONE: - default: - - } - - - //setup reconstruction image size - int[] recoSize = new int[]{128, 128}; - double[] spacing = new double[]{1.0, 1.0}; - - Grid2D reconstruction = fbp.backproject(parkerWeighted, recoSize, spacing); - reconstruction.show("Reconstructed image"); - - - //Compare our backprojection algorithm to the tutorial code - FanBeamBackprojector2D fbp2 = new FanBeamBackprojector2D(focalLength, fbp.detectorSpacing, fbp.betaIncrement,recoSize[0], recoSize[1]); - fbp2.backprojectPixelDriven(parkerWeighted).show("Tutorial reconstruction"); - - } -} - - diff --git a/src/edu/stanford/rsl/tutorial/dmip/DMIP_ParallelBeam.java b/src/edu/stanford/rsl/tutorial/dmip/DMIP_ParallelBeam.java deleted file mode 100644 index 90aa6017..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/DMIP_ParallelBeam.java +++ /dev/null @@ -1,306 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import java.util.ArrayList; - -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.conrad.data.numeric.InterpolationOperators; -import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; -import edu.stanford.rsl.conrad.geometry.shapes.simple.Box; -import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND; -import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine; -import edu.stanford.rsl.conrad.geometry.transforms.Transform; -import edu.stanford.rsl.conrad.geometry.transforms.Translation; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import edu.stanford.rsl.tutorial.phantoms.SheppLogan; -import ij.ImageJ;; - -/** - * Exercise 5 of Diagnostic Medical Image Processing (DMIP) - * @author Bastian Bier - * - */ - -public class DMIP_ParallelBeam { - - - public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN}; - - - /** - * Forward projection of the phantom onto the detector - * Rule of thumb: Always sample in the domain where you expect the output! - * Thus, we sample at the detector pixel positions and sum up the informations along one ray - * - * @param grid the image - * @param maxTheta the angular range in radians - * @param deltaTheta the angular step size in radians - * @param maxS the detector size in [mm] - * @param deltaS the detector element size in [mm] - */ - - public Grid2D projectRayDriven(Grid2D grid, double maxTheta, double deltaTheta, double maxS, double deltaS) { - - int maxSIndex = (int) (maxS / deltaS + 1); - int maxThetaIndex = (int) (maxTheta / deltaTheta + 1); - - final double samplingRate = 3.d; // # of samples per pixel - Grid2D sino = new Grid2D(new float[maxThetaIndex*maxSIndex], maxSIndex, maxThetaIndex); - sino.setSpacing(deltaS, deltaTheta); - - // set up image bounding box in WC - Translation trans = new Translation( - -(grid.getSize()[0] * grid.getSpacing()[0])/2, -(grid.getSize()[1] * grid.getSpacing()[1])/2, -1); - Transform inverse = trans.inverse(); - - Box b = new Box((grid.getSize()[0] * grid.getSpacing()[0]), (grid.getSize()[1] * grid.getSpacing()[1]), 2); - b.applyTransform(trans); - - for(int e=0; e points = b.intersect(line); - - // only if we have intersections - if (2 != points.size()){ - if(points.size() == 0) { - line.getDirection().multiplyBy(-1.d); - points = b.intersect(line); - } - if(points.size() == 0) - continue; - } - - PointND start = points.get(0); // [mm] - PointND end = points.get(1); // [mm] - - // get the normalized increment - SimpleVector increment = new SimpleVector( - end.getAbstractVector()); - increment.subtract(start.getAbstractVector()); - double distance = increment.normL2(); - increment.divideBy(distance * samplingRate); - - double sum = .0; - start = inverse.transform(start); - - // compute the integral along the line. - for (double t = 0.0; t < distance * samplingRate; ++t) { - PointND current = new PointND(start); - current.getAbstractVector().add(increment.multipliedBy(t)); - - double x = current.get(0) / grid.getSpacing()[0], - y = current.get(1) / grid.getSpacing()[1]; - - if (grid.getSize()[0] <= x + 1 - || grid.getSize()[1] <= y + 1 - || x < 0 || y < 0) - continue; - - sum += InterpolationOperators.interpolateLinear(grid, x, y); - } - - // normalize by the number of interpolation points - sum /= samplingRate; - // write integral value into the sinogram. - sino.setAtIndex(i, e, (float)sum); - } - } - return sino; - } - - /** - * Sampling of projections is defined in the constructor. - * Backprojection of the projections/sinogram - * Rule of thumb: Always sample in the domain where you expect the output! - * Here, we want to reconstruct the volume, thus we sample in the reconstructed grid! - * The projections are created pixel driven! - * - * @param sino the sinogram - * @param imageSizeX - * @param imageSizeY - * @param pxSzXMM - * @param pxSzYMM - */ - - public Grid2D backprojectPixelDriven(Grid2D sino, int imageSizeX, int imageSizeY, float pxSzXMM, float pxSzYMM) { - - int maxThetaIndex = sino.getSize()[1]; - double deltaTheta = sino.getSpacing()[1]; - int maxSIndex = sino.getSize()[0]; - double deltaS = sino.getSpacing()[0]; - double maxS = (maxSIndex-1) * deltaS; - - Grid2D grid = new Grid2D(imageSizeX, imageSizeY); - grid.setSpacing(pxSzXMM, pxSzYMM); - grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2); - - // loop over the projection angles - for (int i = 0; i < maxThetaIndex; i++) { - // compute actual value for theta - double theta = deltaTheta * i; - // precompute sine and cosines for faster computation - double cosTheta = Math.cos(theta); - double sinTheta = Math.sin(theta); - // get detector direction vector - SimpleVector dirDetector = new SimpleVector(sinTheta,cosTheta); - // loops over the image grid - for (int x = 0; x < grid.getSize()[0]; x++) { - for (int y = 0; y < grid.getSize()[1]; y++) { - // compute world coordinate of current pixel - double[] w = grid.indexToPhysical(x, y); - // wrap into vector - SimpleVector pixel = new SimpleVector(w[0], w[1]); - // project pixel onto detector - double s = SimpleOperators.multiplyInnerProd(pixel, dirDetector); - // compute detector element index from world coordinates - s += maxS/2; // [mm] - s /= deltaS; // [GU] - // get detector grid - Grid1D subgrid = sino.getSubGrid(i); - // check detector bounds, continue if out of array - if (subgrid.getSize()[0] <= s + 1 - || s < 0) - continue; - // get interpolated value - float val = InterpolationOperators.interpolateLinear(subgrid, s); - // sum value to sinogram - grid.addAtIndex(x, y, val); - } - - } - } - // apply correct scaling - NumericPointwiseOperators.divideBy(grid, (float) (maxThetaIndex / Math.PI)); - return grid; - } - - - /** - * Filtering the sinogram with a high pass filter - * - * The ramp filters are defined in the spatial domain but - * they are applied in the frequency domain. - * Remember: a convolution in the spatial domain corresponds to a - * multiplication in the frequency domain. - * - * Both, the sinogram and the ramp filter are transformed into - * the frequency domain and multiplied there. - * - * @param sinogram a line of the sinogram - * - */ - public Grid1D rampFiltering(Grid1D sinogram, RampFilterType filter){ - - double deltaS = 1; - - // Initialize the ramp filter - // Define the filter in the spatial domain on the full padded size! - Grid1DComplex ramp = new Grid1DComplex(sinogram.getSize()[0]); - - int paddedSize = ramp.getSize()[0]; - - if(filter == RampFilterType.RAMLAK) - { - // TODO: implement the ram-lak filter in the spatial domain - - } - else if(filter == RampFilterType.SHEPPLOGAN) - { - // TODO: implement the Shepp-Logan filter in the spatial domain - - } - else - { - // if no filtering is used - return sinogram; - } - - // TODO: Transform ramp filter into frequency domain - - - Grid1DComplex sinogramF = new Grid1DComplex(sinogram,true); - // TODO: Transform the input sinogram signal into the frequency domain - - - // TODO: Multiply the ramp filter with the transformed sinogram - - - // TODO: Backtransformation - - - // Crop the image to its initial size - Grid1D ret = new Grid1D(sinogram); - ret = sinogramF.getRealSubGrid(0, sinogram.getSize()[0]); - - return ret; - } - - public static void main(String[] args) - { - ImageJ ij = new ImageJ(); - - DMIP_ParallelBeam parallel = new DMIP_ParallelBeam(); - - // 0. Parameters - - // size of the phantom - int phantomSize = 128; - // projection image range - double angularRange = Math.PI; - // number of projection images - int projectionNumber = 180; - // angle in between adjacent projections - double angularStepSize = angularRange / projectionNumber; - // detector size in [mm] - float detectorSize = 200; - // size of a detector Element [mm] - float detectorSpacing = 1.0f; - // filterType: NONE, RAMLAK, SHEPPLOGAN - RampFilterType filter = RampFilterType.NONE; - - // 1. Create the Shepp Logan Phantom - SheppLogan sheppLoganPhantom = new SheppLogan(phantomSize); - sheppLoganPhantom.show(); - - // 2. Acquire forward projection images with a parallel projector - Grid2D sinogram = parallel.projectRayDriven(sheppLoganPhantom, angularRange, angularStepSize, detectorSize, detectorSpacing); - sinogram.show("The Sinogram"); - Grid2D filteredSinogram = new Grid2D(sinogram); - - // 4. Ramp Filtering - for (int theta = 0; theta < sinogram.getSize()[1]; ++theta) - { - // Filter each line of the sinogram independently - Grid1D tmp = parallel.rampFiltering(sinogram.getSubGrid(theta), filter); - - for(int i = 0; i < tmp.getSize()[0]; i++) - { - filteredSinogram.putPixelValue(i, theta, tmp.getAtIndex(i)); - } - } - - filteredSinogram.show("Filtered Sinogram"); - - // 5. Reconstruct the object with the information in the sinogram - Grid2D reco = parallel.backprojectPixelDriven(filteredSinogram, phantomSize, phantomSize, detectorSpacing, detectorSpacing); - reco.show("Reconstruction"); - } - -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/DefectPixelInterpolation.java b/src/edu/stanford/rsl/tutorial/dmip/DefectPixelInterpolation.java deleted file mode 100644 index f9b72faf..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/DefectPixelInterpolation.java +++ /dev/null @@ -1,404 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import java.util.ArrayList; - -import edu.stanford.rsl.conrad.data.generic.datatypes.Complex; -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.data.numeric.Grid2DComplex; -import edu.stanford.rsl.conrad.data.numeric.Grid3D; -import edu.stanford.rsl.conrad.data.numeric.NumericGridOperator; -import edu.stanford.rsl.conrad.filtering.MedianFilteringTool; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import ij.IJ; -import ij.ImageJ; - -/** - * - * Exercise 4 of Diagnostic Medical Image Processing (DMIP) - * @author Marco Boegel - * - */ -public class DefectPixelInterpolation { - - public DefectPixelInterpolation() - { - - } - - /** - * Spectral deconvolution for defect pixel interpolation - * "Defect Interpolation in Digital Radiography - How Object-Oriented Transform Coding Helps", - * T. Aach, V. Metzler, SPIE Vol. 4322: Medical Imaging, February 2001 - * @param image corrupted image with defect pixels - * @param mask binary image mask of the defect pixels (0: defect, 1: fine) - * @param maxIter maximum number of iterations - * @param zeroPadding enables zero padding for the input images. Images are enlarged to size of next power of 2 and filled with zeros - * @return image, where defect pixels are interpolated - */ - public Grid2D interpolateSpectral(Grid2D image, Grid2D mask, int maxIter, boolean zeroPadding) - { - //padding - //TODO - //TODO - - //fourier transform - Grid2DComplex G = new Grid2DComplex(1,1);//TODO - Grid2DComplex W = new Grid2DComplex(1,1);//TODO - //TODO - //TODO - - int[] dim = G.getSize(); - - int[] halfDim = {dim[0]/2, dim[1]/2}; - - double maxDeltaE_G_Ratio = Double.POSITIVE_INFINITY; - double maxDeltaE_G_Thresh = 1.0e-6; - - Grid2DComplex FHat = new Grid2DComplex(dim[0], dim[1], false); - Grid2DComplex FHatNext = new Grid2DComplex(dim[0], dim[1], false); - - //Setup visualization - double lastEredVal = 0; - //Visualize every 100th iteration - Grid3D visualize = new Grid3D(image.getWidth(), image.getHeight(), maxIter/100); - - for(int i = 0; i < maxIter; i++) - { - //Check for convergence - if(maxDeltaE_G_Ratio <= maxDeltaE_G_Thresh) - { - System.out.println("maxDeltaE_G_Ratio = " + maxDeltaE_G_Ratio); - break; - } - - //In the i-th iteration select line pair s1,t1 - //which maximizes the energy reduction [Paragraph after Eq. (16) in the paper] - double maxDeltaE_G = Double.NEGATIVE_INFINITY; - //create arraylist to store lines (in case multiple maxima are found) - //TODO - ArrayList sj1 = null; //HINT - for(int x = 0; x < dim[0]; x++) - { - for(int y = 0; y < dim[1]; y++) - { - float val = G.getAtIndex(x, y) ; - if( val > maxDeltaE_G) - { - //TODO - //TODO - //TODO - }else if(val == maxDeltaE_G) { - //TODO - } - } - } - //if there were more indices than one with the same max_value, pick a random one of these - int idx = (int) Math.floor(Math.random() * sj1.size()); - int s1 = sj1.get(idx)[0]; - int t1 = sj1.get(idx)[1]; - - //Calculate the ratio of energy reduction in comparison to the last iteration - if(i > 0) - { - maxDeltaE_G_Ratio = Math.abs(maxDeltaE_G - lastEredVal/maxDeltaE_G); - } - - lastEredVal = maxDeltaE_G; - - //Compute the corresponding linepair s2, t2: - //mirror the positions at halfDim - int s2 = 0;//TODO - int t2 = 0;//TODO - - //[Paragraph after Eq. (17) in the paper] - int twice_s1 = (2*s1) % dim[0]; - int twice_t1 = (2*t1) % dim[1]; - - - //Estimate FHat - //4 special cases, where only a single line can be selected: - //(0,0), (0, halfHeight), (halfWidth,0), (halfWidth, halfHeight) - boolean specialCase = false; - if( (s1 == 0 && t1 == 0 ) || - (s1 == halfDim[0] && t1 == 0) || - (s1 == 0 && t1 == halfDim[1]) || - (s1 == halfDim[0] && t1 == halfDim[1])) - { - System.out.println("Special Case"); - specialCase = true; - //Eq. 15 - //FHat = N*(G(s,t)/W(0,0)) - //TODO compute FHatNext, use Complex class - //TODO - //TODO - //TODO - Complex res = new Complex();//HINT - FHatNext.setRealAtIndex(s1, t1, (float) res.getReal()); - FHatNext.setImagAtIndex(s1, t1, (float) res.getImag()); - } - else - { - //General case - //Compute FHatNext for the general case Eq.9 - //TODO - //TODO - //TODO - //TODO - //TODO - Complex res_s1t1 = new Complex();//HINT - Complex res_s2t2 = new Complex();//HINT - FHatNext.setRealAtIndex(s1, t1, (float) res_s1t1.getReal()); - FHatNext.setImagAtIndex(s1, t1, (float) res_s1t1.getImag()); - FHatNext.setRealAtIndex(s2, t2, (float) res_s2t2.getReal()); - FHatNext.setImagAtIndex(s2, t2, (float) res_s2t2.getImag()); - } - - //End iteration step by forming the new error spectrum - updateErrorSpectrum(G, FHatNext, FHat, W, s1, t1, specialCase); - - //Get rid of rounding errors - //G(t1,s1) and G(t2,s2) should be zero - G.setAtIndex(s1, t1, 0); - if(!specialCase) - { - G.setAtIndex(s2, t2, 0); - } - - FHat = new Grid2DComplex(FHatNext); - - if(i % 100 == 0) - { - //For visualization, apply IFFT to the estimation - Grid2DComplex FHatV = new Grid2DComplex(FHat); - FHatV.transformInverse(); - Grid2D vis = new Grid2D(image); - - //Fill in the defect mask pixels with current estimation and remove the zero padding - for(int x = 0; x < vis.getWidth(); x++) - { - for(int y = 0; y < vis.getHeight(); y++) - { - if(mask.getAtIndex(x, y) == 0) - { - vis.setAtIndex(x, y, FHatV.getRealAtIndex(x, y)); - } - } - } - visualize.setSubGrid(i/100, vis); - - } - - } - - visualize.show(); - - //Compute the inverse fourier transform of the estimated image - FHat.transformInverse(); - - //Fill in the defect mask pixels with the current estimation and remove the zero padding - Grid2D result = new Grid2D(image); - //TODO - //TODO - //TODO - //TODO - - return result; - } - - - /** - * - * Do the convolution of the m-times-n matrix F and W - * s,t is the position of the selected line pair, the convolution is simplified in the following way: - * G(k1,k2) = F(k1,k2) 'conv' W(k1,k2) - * = (F(s,t)W(k1-s,k2-t) + F*(s,t)W(k1+s,k2+t)) / (MN) - * where F* is the conjugate complex. - * - * @param G Fourier transformation of input image - * @param FHatNext currently estimated FT of the fixed image - * @param FHat previous estimated FT of the fixed image - * @param W Fourier transformation of the mask image - * @param s1 position of the selected line pair - * @param t1 position of the selected line pair - * @param specialCase - */ - private static void updateErrorSpectrum(Grid2DComplex G, Grid2DComplex FHatNext, Grid2DComplex FHat, Grid2DComplex W, int s1, int t1, boolean specialCase) { - int[] sz = FHatNext.getSize(); - - // Accumulation: Update pair (s1,t1),(s2,t2) - Complex F_st = new Complex(FHatNext.getRealAtIndex(s1, t1) - FHat.getRealAtIndex(s1, t1), FHatNext.getImagAtIndex(s1, t1) - FHat.getImagAtIndex(s1, t1)); - Complex F_st_conj = F_st.getConjugate(); - - int MN = sz[0] * sz[1]; - - // Compute the new error spectrum - for(int j = 0; j < sz[1]; j++) - { - for(int i = 0; i < sz[0]; i++) - { - Complex GVal; - if(specialCase) - { - int xneg = (i - s1) % sz[0]; - int yneg = (j - t1) % sz[1]; - - if(xneg < 0) - { - xneg = sz[0] + xneg; - } - if(yneg < 0) - { - yneg = sz[1] + yneg; - } - - GVal = new Complex(G.getRealAtIndex(i, j), G.getImagAtIndex(i, j)); - Complex WNeg = new Complex(W.getRealAtIndex(xneg, yneg), W.getImagAtIndex(xneg, yneg)); - GVal.sub( ( F_st.mul(WNeg) ).div(MN) ); - - } - else - { - int xpos = (i + s1) % sz[0]; - int ypos = (j + t1) % sz[1]; - int xneg = (i - s1) % sz[0]; - int yneg = (j - t1) % sz[1]; - - if(xneg < 0) - { - xneg = sz[0] + xneg; - } - if(yneg < 0) - { - yneg = sz[1] + yneg; - } - - Complex WPos = new Complex(W.getRealAtIndex(xpos, ypos), W.getImagAtIndex(xpos, ypos)); - Complex WNeg = new Complex(W.getRealAtIndex(xneg, yneg), W.getImagAtIndex(xneg, yneg)); - GVal = new Complex(G.getRealAtIndex(i, j), G.getImagAtIndex(i, j)); - GVal = GVal.sub( ( ( F_st.mul(WNeg) ).add( F_st_conj.mul(WPos) ) ).div(MN) ); - - - } - - G.setRealAtIndex(i, j, (float) GVal.getReal()); - G.setImagAtIndex(i, j, (float) GVal.getImag()); - } - } - } - - public Grid2D interpolateMedian(Grid2D image, Grid2D defects, int kernelWidth, int kernelHeight) - { - //Pad the image. Otherwise, the filter will ignore kernelWidth/2 at each side of the image - Grid2D paddedImage = new Grid2D(image.getWidth()+kernelWidth, image.getHeight()+kernelHeight); - - for(int i = 0; i (d-1): extrema - // d=0: constant (horizontal line with y-intercept a_0 -> f(x)=a_0) - // d=1: oblique line with y-intercept a_0 & slope a_1 -> f(x)=a_0 + a_1 x - // d=2: parabola - // d>=2: continuous non-linear curve - // E.g. d=5: 4 extrema - // d = 10 -> NumKoeff: 66 -> but only 64 lattice points are known - int degree = 5; //Polynomial's degree: 2,...,10 - - // Number of Coefficients - // TODO: - int numCoeff = 0; - - // Number of Correspondences - // TODO: - int numCorresp = 0; - - // Print out of the used parameters - System.out.println("Polynom of degree: " + degree); - System.out.println("Number of Coefficients: " + numCoeff); - System.out.println("Number of Correspondences: " + numCorresp); - - // 3.Create the matrix A - SimpleMatrix A = new SimpleMatrix(numCorresp, numCoeff); - A.zeros(); - - // Realign the grid matrix into a vector - // Easier access in the next step - SimpleVector Xu2_vec = new SimpleVector(numCorresp); - SimpleVector Yu2_vec = new SimpleVector(numCorresp); - SimpleVector Xd2_vec = new SimpleVector(numCorresp); - SimpleVector Yd2_vec = new SimpleVector(numCorresp); - - for(int i = 0; i < ny; i++) - { - for(int j = 0; j < nx; j++) - { - Xu2_vec.setElementValue(i * ny + j, Xu2.getElement(j, i)); - Yu2_vec.setElementValue(i * ny + j, Yu2.getElement(j, i)); - Xd2_vec.setElementValue(i * ny + j, Xd2.getElement(j, i)); - Yd2_vec.setElementValue(i * ny + j, Yd2.getElement(j, i)); - } - } - - // Compute matrix A - for(int r = 0; r < numCorresp; r++) - { - int cc = 0; - for(int i = 0; i <= degree; i++) - { - for(int j = 0; j <= (degree-i); j++) - { - // TODO: - - } - } - } - - // Compute the pseudo-inverse of A with the help of the SVD (class: DecompositionSVD) - // TODO - // TODO - - - // Compute the distortion coefficients - // TODO - // TODO - - - // 4. Compute the distorted grid points (xDist, yDist) which are used to sample the - // distorted image to get the undistorted image - // (x,y) is the position in the undistorted image and (XDist,YDist) the - // position in the distorted (observed) X-ray image. - - Grid2D xDist = new Grid2D(imSize, imSize); - Grid2D yDist = new Grid2D(imSize, imSize); - - for(int x = 0; x < imSize; x++) - { - for(int y = 0; y < imSize; y++) - { - int cc = 0; - for(int k = 0; k <= degree; k++) - { - for(int l = 0; l <= degree - k; l++) - { - // TODO - // TODO - // TODO - } - } - } - } - - Grid2D undistortedImage = new Grid2D(imSize, imSize); - - for(int i = 0; i < imSize; i++) - { - for(int j = 0; j < imSize; j++) - { - // TODO - // TODO - } - } - undistortedImage.show("Undistorted Image"); - - Grid2D differenceImage = (Grid2D) NumericPointwiseOperators.subtractedBy(quadraticImage, undistortedImage); - differenceImage.show("diffImage"); - } - - - - public static void main(String[] args) { - ImageJ ij = new ImageJ(); - ImageUndistortion iu = new ImageUndistortion(); - iu.doImageUndstortion(); - } -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/Intro.java b/src/edu/stanford/rsl/tutorial/dmip/Intro.java deleted file mode 100644 index 3ed8c5bb..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/Intro.java +++ /dev/null @@ -1,253 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.numerics.DecompositionSVD; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix.MatrixNormType; -import edu.stanford.rsl.conrad.numerics.SimpleVector.VectorNormType; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import edu.stanford.rsl.conrad.utils.VisualizationUtil; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import ij.IJ; -import ij.ImageJ; -import ij.plugin.filter.Convolver; -import ij.process.FloatProcessor; - - -/** - * Introduction to the CONRAD Framework - * Exercise 1 of Diagnostic Medical Image Processing (DMIP) - * @author Marco Boegel - * - */ - -public class Intro { - - - public static void gridIntro(){ - - //Define the image size - int imageSizeX = 256; - int imageSizeY = 256; - - //Define an image - //Hint: Import the package edu.stanford.rsl.conrad.data.numeric.Grid2D - //TODO - - //Draw a circle - int radius = 50; - //Set all pixels within the circle to 100 - int insideVal = 100; - - //TODO - //TODO - //TODO - - //Show ImageJ GUI - ImageJ ij = new ImageJ(); - //Display image - //TODO - - //Copy an image - //TODO - //copy.show("Copy of circle"); - - - //Load an image from file - String filename = "D:/02_lectures/DMIP/exercises/2014/matlab_intro/mr12.dcm"; - //TODO. Hint: Use IJ and ImageUtil - //mrImage.show(); - - //convolution - //TODO - //TODO - - //define the kernel. Try simple averaging 3x3 filter - int kw = 3; - int kh = 3; - float[] kernel = new float[kw*kh]; - for(int i = 0; i < kernel.length; i++) - { - kernel[i] = 1.f / (kw*kh); - } - - //TODO - - - //write an image to disk, check the supported output formats - String outFilename ="D:/02_lectures/DMIP/exercises/2014/matlab_intro/mr12out.tif"; - //TODO - } - - - public static void signalIntro() - { - //How can I plot a sine function sin(2*PI*x)? - double stepSize = 0.01; - int plotLength = 500; - - double[] y = new double[plotLength]; - - for(int i = 0; i < y.length; i++) - { - //TODO - - } - - VisualizationUtil.createPlot(y).show(); - double[] x = new double [plotLength]; - for(int i = 0; i < x.length; i++) - { - x[i] = (double) i * stepSize; - } - - VisualizationUtil.createPlot(x, y, "sin(x)", "x", "y").show(); - - } - - public static void basicIntro() - { - //Display text - System.out.println("Creating a vector: v1 = [1.0; 2.0; 3.0]"); - - //create column vector - //TODO - //System.out.println("v1 = " + v1.toString()); - - //create a randomly initialized vector - SimpleVector vRand = new SimpleVector(3); - //TODO - //System.out.println("vRand = " + vRand.toString()); - - //create matrix M 3x3 1 2 3; 4 5 6; 7 8 9 - SimpleMatrix M = new SimpleMatrix(); - //TODO - //System.out.println("M = " + M.toString()); - - //determinant of M - //System.out.println("Determinant of matrix m: " + TODO ); - - //transpose M - //TODO - //copy matrix - //TODO - //transpose M inplace - //TODO - - //get size - int numRows = 0; - int numCols = 0; - //TODO - - //access elements of M - System.out.println("M: "); - for(int i = 0 ; i < numRows; i++) - { - for(int j = 0; j < numCols; j++) - { - //TODO - //System.out.print(element + " "); - } - System.out.println(); - } - - //Create 3x3 Matrix of 1's - SimpleMatrix Mones = new SimpleMatrix(3,3); - //TODO - //Create a 3x3 Matrix of 0's - SimpleMatrix Mzeros = new SimpleMatrix(3,3); - //TODO - //Create a 3x3 Identity matrix - SimpleMatrix Midentity = new SimpleMatrix(3,3); - //TODO - - //Matrix multiplication - //TODO - //System.out.println("M^T * M = " + ResMat.toString()); - - - //Matrix vector multiplication - //TODO - //System.out.println("M * v1 = " + resVec.toString()); - - - //Extract the last column vector from matrix M - //SimpleVector colVector = M.getCol(2); - //Extract the 1x2 subvector from the last column of matrix M - //TODO - //System.out.println("[m(0)(2); m(1)(2)] = " + subVector); - - //Matrix elementwise multiplication - //TODO - //System.out.println("M squared Elements: " + MsquaredElem.toString()); - - //round vectors - SimpleVector vRandCopy = new SimpleVector(vRand); - System.out.println("vRand = " + vRandCopy.toString()); - - vRandCopy.floor(); - System.out.println("vRand.floor() = " + vRandCopy.toString()); - - vRand.ceil(); - System.out.println("vRand.ceil() = " + vRand.toString()); - - //min, max, mean - //double minV1 = v1.min(); - //double maxV1 = v1.max(); - //System.out.println("Min(v1) = " + minV1 + " Max(v1) = " + maxV1); - - //for matrices: iterate over row or column vectors - SimpleVector maxVec = new SimpleVector(M.getCols()); - for(int i = 0; i < M.getCols(); i++) - { - maxVec.setElementValue(i, M.getCol(i).max()); - } - double maxM = maxVec.max(); - System.out.println("Max(M) = " + maxM); - - - - //Norms - //TODO matrix L1 - //TODO vector L2 - //System.out.println("||M||_F = " + matrixNormL1); - //System.out.println("||colVec||_2 = " + vecNormL2); - - //get normalized vector - //TODO - //normalize vector in-place - //TODO - //System.out.println("Normalized colVector: " + colVector.toString()); - - - //SVD - SimpleMatrix A = new SimpleMatrix(3,3); - A.setRowValue(0, new SimpleVector(11, 10, 14)); - A.setRowValue(1, new SimpleVector(12, 11, -13)); - A.setRowValue(2, new SimpleVector(14, 13, -66)); - - System.out.println("A = " + A.toString()); - - //TODO SVD - - //print singular matrix - //System.out.println(svd.getS().toString()); - - //get condition number - //System.out.println("Condition number of A: " + TODO ); - - //Re-compute A = U * S * V^T - //SimpleMatrix temp = SimpleOperators.multiplyMatrixProd(svd.getU(), svd.getS()); - //SimpleMatrix A2 = SimpleOperators.multiplyMatrixProd(temp, svd.getV().transposed()); - //System.out.println("U * S * V^T: " + A2.toString()); - - } - - public static void main(String arg[]) - { - basicIntro(); - gridIntro(); - signalIntro(); - } -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/RANSAC.java b/src/edu/stanford/rsl/tutorial/dmip/RANSAC.java deleted file mode 100644 index 35006f56..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/RANSAC.java +++ /dev/null @@ -1,212 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import ij.gui.Plot; -import java.util.ArrayList; - -/** - * Exercise 4 of Diagnostic Medical Image Processing (DMIP) - * @author Bastian Bier - * - */ - -public class RANSAC { - - /** - * Function calculation the ransac model estimate - * - * @param points The point cloud where the line should be fit through - * @param minParam The number of parameters required for the model - * @param p_opt The probability of picking the right points for the models during the iterations - * @param p_out The probability of an outliner - * - * @return solution the parameters of the resulting RANSAC line in a SimpleVector: [m,c] - */ - - public SimpleVector commonRansac(SimpleMatrix points, int mn, double p_opt, double p_out){ - - // Calculating the amount of required iterations - int it = (int) (Math.log(1 - p_opt) / Math.log(1 - Math.pow(1 - p_out, mn)) + 1 ); - - // Error of the best fitted line - double error = Double.POSITIVE_INFINITY; - - // Solution vector - SimpleVector solution = new SimpleVector(mn); - - for(int i = 0; i < it; i++) - { - // Select mn random points - // Calculate the indexes of the points - ArrayList indexes = new ArrayList(); - int randIdx = (int) (Math.random() / (1.f / points.getRows())); - indexes.add(randIdx); - - for(int n = 1; n < mn; n++){ - randIdx = (int) (Math.random() / (1.f / points.getRows())); - - while(indexes.contains(randIdx)) - { - randIdx = (int) (Math.random() / (1.f / points.getRows())); - } - indexes.add(randIdx); - } - - // Calculate the parameters - SimpleMatrix a = new SimpleMatrix(mn,mn); - - for(int n = 0; n < mn; n++){ - a.setRowValue(n, points.getRow(indexes.get(n))); - } - - SimpleVector lineParams= new SimpleVector(mn); - // TODO: estimate the line parameters for the selected points - - - - - // Calculate the error of the estimated line - // update the error and the parameters, if the current line has a smaller error - double cur_err = 0.0; - // TODO: calculate the error of the current line - - - if(cur_err < error) - { - error = cur_err; - solution = lineParams; - } - - } - - return solution; - } - - /** - * Function calculating a line through a point cloud using the SVD - * - * @param points The point cloud where the line should be fit through - * 2 points result in a exact line - * >2 points result in a regression line - * - * @return x_result the parameters of the line in a SimpleVector: [m,c] - */ - public SimpleVector fitline(SimpleMatrix points){ - - // Build up the measurement matrix - SimpleMatrix m = new SimpleMatrix(points.getRows(),2); - SimpleVector b = new SimpleVector(points.getRows()); - m.fill(1); - - for(int i = 0; i < points.getRows(); i++){ - m.setElementValue(i, 0, points.getElement(i, 0)); - b.setElementValue(i, points.getElement(i, 1)); - } - - // Solution vector containing the estimated parameters m and c - SimpleVector x_result = new SimpleVector(2); - - // Calculate the parameters using the Pseudo-Inverse - // TODO: calculate the line parameters, write them in x_result - - return x_result; - } - - /** - * Calculate the error of a line - * - * @param line_params Parameters of the line - * @param points The point cloud where the line should be fit through. - * - * @return error the calculated error - */ - public double lineError(SimpleVector line_params, SimpleMatrix points){ - - // Threshold defining the allowed distance of a point to the line - double thresh = 0.2; - - // TODO: line parameters - - - // TODO: get some point on the line - - - // TODO: calculate normal vector of the line - - - // TODO: calculate distance line to origin - - - // TODO: calculate the distance for each point to the line - // TODO: check if the distance is higher than the threshold - - - // TODO: return the error - return 0; - } - - - public static void main(String[] args) { - - // - RANSAC ransac = new RANSAC(); - - // - // The point cloud is defined - // - - SimpleMatrix pts = new SimpleMatrix(7,2); - pts.setRowValue(0, new SimpleVector(0,0)); - pts.setRowValue(1, new SimpleVector(1,1)); - pts.setRowValue(2, new SimpleVector(2,2)); - pts.setRowValue(3, new SimpleVector(3,3)); - pts.setRowValue(4, new SimpleVector(3.2,1.9)); - pts.setRowValue(5, new SimpleVector(4,4)); - pts.setRowValue(6, new SimpleVector(10,1.8)); - - - // - // Regression Line - // - - // Create a scatter plot of the point cloud and fit a regression line - Plot scatterPlot = new Plot("Regression Line", "X", "Y", Plot.DEFAULT_FLAGS); - scatterPlot.setLimits(0, 11, 0, 5); - scatterPlot.addPoints(pts.getCol(0).copyAsDoubleArray(), pts.getCol(1).copyAsDoubleArray(), Plot.BOX); - scatterPlot.show(); - - // Calculate the regression line through the given point cloud - SimpleVector regressionLine = ransac.fitline(pts); - - // Add the regression line - double y11 = regressionLine.getElement(0) * 11 + regressionLine.getElement(1); - double y0 = regressionLine.getElement(0) * 0 + regressionLine.getElement(1); - scatterPlot.drawLine(0, y0, 11, y11); - - - // - // RANSAC - // - - // Parameters for RANSAC - double p_opt = 0.9999; // probability how likely it is to pick the right mn points - double p_out = 0.2; // probability of an outlier - int min_number = 2; // minimum number of datapoints required to build the model - - // Create a scatter plot of the point cloud and fit a RANSAC line - Plot ransac_plot = new Plot("Ransac Line", "X", "Y", Plot.DEFAULT_FLAGS); - ransac_plot.setLimits(0, 11, 0, 5); - ransac_plot.addPoints(pts.getCol(0).copyAsDoubleArray(), pts.getCol(1).copyAsDoubleArray(), Plot.BOX); - - - // Compute a line using the RANSAC algorithm and plot it - SimpleVector ransacLine = ransac.commonRansac(pts, min_number, p_opt, p_out); - double y1 = ransacLine.getElement(0) * 0 + ransacLine.getElement(1); - double y2 = ransacLine.getElement(0) * 11 + ransacLine.getElement(1); - ransac_plot.drawLine(0, y1, 11, y2); - ransac_plot.show(); - - } -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/Registration1.java b/src/edu/stanford/rsl/tutorial/dmip/Registration1.java deleted file mode 100644 index 00b33287..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/Registration1.java +++ /dev/null @@ -1,153 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import ij.gui.Plot; - - -/** - * Exercise 7 of Diagnostic Medical Image Processing (DMIP) - * @author Bastian Bier - * - */ -public class Registration1 { - - /** - * Registration method - * Point Correspondences are used to calculate the rotation and translation - * Transformation from p to q - * @param p reference point cloud - * @param q point cloud to be registered - * @return translation and rotation for the registration (phi, t1, t2) - */ - private SimpleVector registrationUsingPointCorrespondences(SimpleMatrix p, SimpleMatrix q){ - - int numPoints = p.getRows(); - - // Build up measurement matrix m - SimpleMatrix m = new SimpleMatrix(numPoints * 2, 4); - - for(int i = 0; i < numPoints * 2; i++) - { - if(i < numPoints) - { - // TODO - // TODO - // TODO - // TODO - } - if(i >= numPoints) - { - // TODO - // TODO - // TODO - // TODO - } - } - - // Build up solution vector b - SimpleVector b = new SimpleVector(2 * numPoints); - - for(int i = 0; i < 2 * numPoints; i++) - { - if(i < numPoints) - { - // TODO - } - - if(i >= numPoints) - { - // TODO - } - } - - // Calculate Pseudo Inverse of m - SimpleMatrix m_inv; - m_inv = m.inverse(InversionType.INVERT_SVD); - - // Calculate Parameters with the help of the Pseudo Inverse - SimpleVector x; - x = SimpleOperators.multiply(m_inv, b); - - double r1 = x.getElement(0); - double r2 = x.getElement(1); - double t1 = x.getElement(2); - double t2 = x.getElement(3); - - - // TODO: normalize r - - - double phi = 0; // TODO - - // Write the result for the translation and the rotation into the result vector - SimpleVector result = new SimpleVector(phi, t1, t2); - - return result; - } - - - private SimpleMatrix applyTransform(SimpleMatrix points, double phi, SimpleVector translation){ - - SimpleMatrix r = new SimpleMatrix(2,2); - // TODO: fill the rotation matrix - - // TODO - // TODO - // TODO - // TODO - - SimpleMatrix transformedPoints = new SimpleMatrix(points.getRows(), points.getCols()); - - for(int i = 0; i < transformedPoints.getRows(); i++) - { - // TODO: transform points - - } - - return transformedPoints; - } - - public static void main(String[] args){ - - // Define Point Cloud - SimpleMatrix p_k = new SimpleMatrix(4,2); - SimpleMatrix q_k = new SimpleMatrix(4,2); - - p_k.setRowValue(0, new SimpleVector(1,2)); - p_k.setRowValue(1, new SimpleVector(3,6)); - p_k.setRowValue(2, new SimpleVector(4,6)); - p_k.setRowValue(3, new SimpleVector(3,4)); - - q_k.setRowValue(0, new SimpleVector(-0.7, 2.1)); - q_k.setRowValue(1, new SimpleVector(-2.1, 6.4)); - q_k.setRowValue(2, new SimpleVector(-1.4, 7.1)); - q_k.setRowValue(3, new SimpleVector(-0.7, 4.9)); - - // Plot Point Cloud - Plot plot = new Plot("Regression Line", "X", "Y", Plot.DEFAULT_FLAGS); - plot.setLimits(-10, 10, -10, 10); - plot.addPoints(p_k.getCol(0).copyAsDoubleArray(), p_k.getCol(1).copyAsDoubleArray(), Plot.BOX); - plot.addPoints(q_k.getCol(0).copyAsDoubleArray(), q_k.getCol(1).copyAsDoubleArray(), Plot.CIRCLE); - - // Calculate registration parameter - Registration1 reg = new Registration1(); - SimpleVector parameter = reg.registrationUsingPointCorrespondences(p_k, q_k); - - // Rotation and translation - double phi = parameter.getElement(0); - SimpleVector translation = new SimpleVector(parameter.getElement(1), parameter.getElement(2)); - - // Transform points - SimpleMatrix transformedPoints = reg.applyTransform(q_k, phi, translation); - - // Add transformed point cloud to the plot - plot.addPoints(transformedPoints.getCol(0).copyAsDoubleArray(), transformedPoints.getCol(1).copyAsDoubleArray(), Plot.CROSS); - plot.show(); - - } - - -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/Registration2.java b/src/edu/stanford/rsl/tutorial/dmip/Registration2.java deleted file mode 100644 index fa9b1f7c..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/Registration2.java +++ /dev/null @@ -1,325 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import java.util.ArrayList; -import java.util.Iterator; -import java.util.TreeMap; -import java.util.Map.Entry; - -import edu.stanford.rsl.conrad.data.numeric.Grid1D; -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.data.numeric.Grid3D; -import edu.stanford.rsl.conrad.geometry.transforms.AffineTransform; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import edu.stanford.rsl.conrad.utils.VisualizationUtil; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import edu.stanford.rsl.jpop.FunctionOptimizer; -import edu.stanford.rsl.jpop.OptimizableFunction; -import edu.stanford.rsl.jpop.OptimizationOutputFunction; -import edu.stanford.rsl.jpop.FunctionOptimizer.OptimizationMode; - -import ij.IJ; -import ij.ImageJ; -import ij.gui.PlotWindow; - - -/** - * Exercise 7 of Diagnostic Medical Image Processing (DMIP) - * Using Mutual Information to solve the registration problem - * @author Bastian Bier - * - */ -public class Registration2 { - - public Grid2D reference = null; - public Grid2D image = null; - - public Grid3D movingStack = null; - - public int iteration = 0; - double[] saveParameters; - - class CostFunction implements OptimizableFunction, OptimizationOutputFunction{ - - private TreeMap resultVisualizer; - - private PlotWindow resultVisualizerPlot; - - @Override - public void setNumberOfProcessingBlocks(int number) { - } - - @Override - public int getNumberOfProcessingBlocks() { - return 1; - } - - @Override - public double evaluate(double[] x, int block) { - - int nrChanges = 0; - for (int i = 0; i < x.length; i++) { - if(saveParameters[i]!=x[i]){ - nrChanges++; - } - } - - // Define Rotation - SimpleMatrix r = new SimpleMatrix(2, 2); - double phi2 = x[0] * (2 * Math.PI) / 360; - r.setElementValue(0, 0, Math.cos(phi2)); - r.setElementValue(0, 1, -Math.sin(phi2)); - r.setElementValue(1, 0, Math.sin(phi2)); - r.setElementValue(1, 1, Math.cos(phi2)); - - // Define translation - double t_x = x[1]; - double t_y = x[2]; - - Grid2D im_tmp = new Grid2D(image); - - // Perform rotation/translation - SimpleVector t = new SimpleVector(t_x, t_y); - AffineTransform affine = new AffineTransform(r, t); - im_tmp.applyTransform(affine); - - if (nrChanges >= 3) { - movingStack.setSubGrid(iteration, im_tmp); - iteration++; - } - - // Calculate the cost function - double cost = calculateMutualInformation(reference, im_tmp); - System.arraycopy(x, 0, saveParameters, 0, x.length); - - return cost; - } - - @Override - public void optimizerCallbackFunction(int currIterationNumber, double[] x, double currFctVal, - double[] gradientAtX) { - // Visualization of cost function value over time - if (this.resultVisualizer == null) - resultVisualizer = new TreeMap(); - resultVisualizer.put(currIterationNumber, currFctVal); - if (resultVisualizerPlot != null) - resultVisualizerPlot.close(); - - Grid1D out = new Grid1D(resultVisualizer.size()); - Iterator> it = resultVisualizer.entrySet().iterator(); - while (it.hasNext()) { - Entry e = it.next(); - out.setAtIndex(e.getKey(), e.getValue().floatValue()); - } - resultVisualizerPlot = VisualizationUtil.createPlot(out.getBuffer()).show(); - } - } - - /** - * Method to calculate the Mutual Information - * @param ref reference image - * @param mov moving image - * @return negative mutual information - */ - private double calculateMutualInformation(Grid2D ref, Grid2D mov){ - - int histSize = 256; - - // Step 1: Calculate joint histogram - SimpleMatrix jointHistogram = calculateJointHistogram(ref, mov); - - // Step 2: Get histogram for a single image from the joint histogram - // a) for the first image - SimpleVector histo1 = new SimpleVector(histSize); - histo1 = getHistogramFromJointHistogram(jointHistogram); - - // b) for the second image - SimpleVector histo2 = new SimpleVector(histSize); - SimpleMatrix jh_t = jointHistogram.transposed(); - histo2 = getHistogramFromJointHistogram(jh_t); - - // Step 3: Calculate the marginal entropies and the joint entropy - double entropy_jointHisto = 0; - double entropy_histo1 = 0; - double entropy_histo2 = 0; - - for(int i = 0; i < histSize; i++) - { - if(histo1.getElement(i) != 0) - { - // TODO: calculate entropy for histogram 1 - } - - if(histo2.getElement(i) != 0) - { - // TODO: calculate entropy for histogram 2 - } - } - - for (int i = 0; i < histSize; i++) - { - for (int j = 0; j < histSize; j++) - { - if(jointHistogram.getElement(i, j) != 0) - { - // TODO: calculate entropy of the joint histogram - } - } - } - - // make sure to consider the - in from of the sum (Entropy formula) - // TODO - // TODO - // TODO - - // Step 4: Calculate the mutual information - // Note: The mutual information is high for a good match - // but we require a minimization problem --> the result is inverted to fit the optimizer - double mutual_information = 0; - // TODO: calculate the mutual information - - return mutual_information * 1000; - } - - /** - * Method to calculate the joint histogram of two images - * @param im1 image1 - * @param im2 image2 - * @return a SimpleMatrix corresponding to the joint histogram - */ - private SimpleMatrix calculateJointHistogram(Grid2D im1, Grid2D im2){ - - // Calculate joint histogram - int histSize = 256; - SimpleMatrix jH = new SimpleMatrix(histSize, histSize); - - for (int i = 0; i < histSize; i++) { - for (int j = 0; j < histSize; j++) { - // TODO - } - } - - // Divide by the number of elements in order to get probabilities - for (int i = 0; i < histSize; i++) { - for (int j = 0; j < histSize; j++) { - // TODO - } - } - - return jH; - } - - /** - * Method to calculate a histogram from a joint histogram - * @param jH The joint histogram - * @return a SimpleVector corresponding to the marginal histogram - */ - private SimpleVector getHistogramFromJointHistogram(SimpleMatrix jH){ - - // Calculate histogram from joint histogram - int histSize = 256; - SimpleVector hist = new SimpleVector(histSize); - hist.zeros(); - - for(int i = 0; i < histSize; i++) - { - for(int j = 0; j < histSize; j++) - { - // TODO: sum up over the columns - } - } - - return hist; - } - - private double[] performOptimization(){ - - FunctionOptimizer fo = new FunctionOptimizer(); - - fo.setDimension(3); - fo.setNdigit(6); - fo.setItnlim(50); - fo.setMsg(16); - fo.setInitialX(new double[]{0,0,0}); - fo.setMaxima(new double[]{50,50,50}); - fo.setMinima(new double[]{-50,-50,-50}); - fo.setOptimizationMode(OptimizationMode.Function); - - CostFunction cF = new CostFunction(); - - movingStack = new Grid3D(reference.getWidth(), reference.getHeight(), 1000, false); - iteration = 0; - saveParameters = new double[]{Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY}; - - // Optimization visualized - ArrayList visFcts = new ArrayList(); - visFcts.add(cF); - fo.setCallbackFunctions(visFcts); - - double result[] = fo.optimizeFunction(cF); - - return result; - } - - public static void main(String[] args){ - - ImageJ ij = new ImageJ(); - - // Load images - // TODO Adjust paths - String filename1 = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/T1.png"; - String filename2 = "C:/StanfordRepo/CONRAD/src/edu/stanford/rsl/tutorial/dmip/Proton.png"; - - Grid2D image1 = ImageUtil.wrapImagePlus(IJ.openImage(filename1)).getSubGrid(0); - Grid2D image2 = ImageUtil.wrapImagePlus(IJ.openImage(filename2)).getSubGrid(0); - - image1.show("Input Image 1"); - image2.show("Input Image 2"); - - // Set the Origin of the image in its center - // The default origin of an image is in its top left corner - // Default Origin: [0.0, 0.0] - int w = image1.getWidth(); - int h = image1.getHeight(); - - image1.setOrigin(-(w-1) / 2 , -(h-1)/2); - image2.setOrigin(-(w-1) / 2 , -(h-1)/2); - image1.setSpacing(1); - image2.setSpacing(1); - - // Blurred Images for the registration to avoid local minima during optimization - Grid2D image1_blurred = new Grid2D(image1); - Grid2D image2_blurred = new Grid2D(image2); - - IJ.run(ImageUtil.wrapGrid(image1_blurred,""),"Gaussian Blur...", "sigma=4"); - IJ.run(ImageUtil.wrapGrid(image2_blurred,""),"Gaussian Blur...", "sigma=4"); - - Registration2 reg2 = new Registration2(); - reg2.reference = image1_blurred; - reg2.image = image2_blurred; - - // Perform Optimization - double res[] = reg2.performOptimization(); - - // Stack for visualization purposes only - Grid3D optimizationStepsGrid = new Grid3D(reg2.reference.getWidth(), reg2.reference.getHeight(), reg2.iteration, false); - for (int i = 0; i < optimizationStepsGrid.getSize()[2]; i++) { - optimizationStepsGrid.setSubGrid(i, reg2.movingStack.getSubGrid(i)); - } - optimizationStepsGrid.show("Optimization Steps"); - - // Transform image back - SimpleMatrix r = new SimpleMatrix(2,2); - Grid2D registeredImage = new Grid2D(image2); - double phi = (2 * Math.PI) / 360 * res[0]; - r.setElementValue(0, 0, Math.cos(phi)); - r.setElementValue(0, 1, -Math.sin(phi)); - r.setElementValue(1, 0, Math.sin(phi)); - r.setElementValue(1, 1, Math.cos(phi)); - - SimpleVector t2 = new SimpleVector(res[1], res[2]); - AffineTransform affine2 = new AffineTransform(r, t2); - registeredImage.applyTransform(affine2); - registeredImage.show("Registered Image"); - } -} \ No newline at end of file diff --git a/src/edu/stanford/rsl/tutorial/dmip/Registration3.java b/src/edu/stanford/rsl/tutorial/dmip/Registration3.java deleted file mode 100644 index 1edda7f4..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/Registration3.java +++ /dev/null @@ -1,254 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import java.util.ArrayList; -import java.util.Iterator; -import java.util.TreeMap; -import java.util.Map.Entry; - -import edu.stanford.rsl.conrad.data.numeric.Grid1D; -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.data.numeric.Grid3D; -import edu.stanford.rsl.conrad.geometry.transforms.AffineTransform; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import edu.stanford.rsl.conrad.utils.VisualizationUtil; -import edu.stanford.rsl.jpop.FunctionOptimizer; -import edu.stanford.rsl.jpop.FunctionOptimizer.OptimizationMode; -import edu.stanford.rsl.jpop.OptimizableFunction; -import edu.stanford.rsl.jpop.OptimizationOutputFunction; -import edu.stanford.rsl.tutorial.phantoms.SheppLogan; -import ij.IJ; -import ij.ImageJ; -import ij.gui.PlotWindow; - -/** - * Exercise 7 of Diagnostic Medical Image Processing (DMIP) - * Solve the Registration Problem using Sum of Squared Difference - * @author Bastian Bier - * - */ -public class Registration3 { - - public Grid2D reference = null; - public Grid2D image = null; - - public Grid3D movingStack = null; - - public int iteration = 0; - double[] saveParameters; - - public class CostFunction implements OptimizableFunction, OptimizationOutputFunction{ - - - private TreeMap resultVisualizer; - - private PlotWindow resultVisualizerPlot; - - @Override - public void setNumberOfProcessingBlocks(int number) { - } - - @Override - public int getNumberOfProcessingBlocks() { - return 1; - } - - /* - * This function gets a parameter vector and returns the result of the cost function - */ - @Override - public double evaluate(double[] x, int block) { - int nrChanges = 0; - for (int i = 0; i < x.length; i++) { - if(saveParameters[i]!=x[i]){ - nrChanges++; - } - } - - // Define Rotation - SimpleMatrix r = new SimpleMatrix(2,2); - double phi2 = x[0] * (2*Math.PI)/360; - r.setElementValue(0, 0, Math.cos(phi2)); - r.setElementValue(0, 1, - Math.sin(phi2)); - r.setElementValue(1, 0, Math.sin(phi2)); - r.setElementValue(1, 1, Math.cos(phi2)); - - // Define translation - double t_x = x[1]; - double t_y = x[2]; - - Grid2D im_tmp = new Grid2D(image); - - // Perform rotation/translation - SimpleVector t = new SimpleVector(t_x,t_y); - AffineTransform affine = new AffineTransform(r, t); - im_tmp.applyTransform(affine); - - if(nrChanges>=3){ - movingStack.setSubGrid(iteration, im_tmp); - iteration++; - } - - // Calculate the cost function - double cost = SumOfSquaredDifferences(reference, im_tmp); - System.arraycopy(x, 0, saveParameters, 0, x.length); - - return cost; - } - - @Override - public void optimizerCallbackFunction(int currIterationNumber, double[] x, double currFctVal, - double[] gradientAtX) { - // Visualization of cost function value over time - if (this.resultVisualizer == null) - resultVisualizer = new TreeMap(); - resultVisualizer.put(currIterationNumber, currFctVal); - if (resultVisualizerPlot != null) - resultVisualizerPlot.close(); - - Grid1D out = new Grid1D(resultVisualizer.size()); - Iterator> it = resultVisualizer.entrySet().iterator(); - while (it.hasNext()) { - Entry e = it.next(); - out.setAtIndex(e.getKey(), e.getValue().floatValue()); - } - resultVisualizerPlot = VisualizationUtil.createPlot(out.getBuffer()).show(); - } - - } - - private double SumOfSquaredDifferences(Grid2D ref, Grid2D imageMoving){ - - double sum = 0.0; - - for(int i = 0; i < ref.getWidth(); i++) - { - for(int j = 0; j < ref.getHeight(); j++) - { - // TODO: calculate SSD - } - } - - return sum/ref.getNumberOfElements(); - } - - public double[] performOptimization(){ - - // Initialize optimization class - FunctionOptimizer fo = new FunctionOptimizer(); - fo.setDimension(3); - fo.setItnlim(50); - movingStack = new Grid3D(reference.getWidth(), reference.getHeight(), 1000, false); - iteration = 0; - fo.setOptimizationMode(OptimizationMode.Function); - fo.setNdigit(8); - fo.setMsg(16); - fo.setInitialX(new double[]{0,0,0}); - fo.setMaxima(new double[]{50,50,50}); - fo.setMinima(new double[]{-50,-50,-50}); - - // Initialize the Costfunction of the optimization - CostFunction cf = new CostFunction(); - saveParameters = new double[]{Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY}; - - // Optimization visualized - ArrayList visFcts = new ArrayList(); - visFcts.add(cf); - fo.setCallbackFunctions(visFcts); - - // Perform the optimization with the given cost function - double[] result = fo.optimizeFunction(cf); - - return result; - } - - - - - public static void main(String[] args){ - - ImageJ ij = new ImageJ(); - - /////////////////////////////////////////////////////// - // Part 1: Apply a transformation on a phantom image // - /////////////////////////////////////////////////////// - - // Create Phantom - Grid2D phantom = new SheppLogan(256); - - // Set the Origin of the image in its center - // The default origin of an image is in its top left corner - // Default Origin: [0.0, 0.0] - int w = phantom.getWidth(); - int h = phantom.getHeight(); - phantom.setOrigin(-(w-1) / 2 , -(h-1)/2); - - Grid2D phantom_blurred = new Grid2D(phantom); - IJ.run(ImageUtil.wrapGrid(phantom_blurred,""),"Gaussian Blur...", "sigma=3"); - phantom.show("Phantom"); - - // Rotate the phantom by 45° and translate it with t = [20, 1] - - // Define Rotation and translation - SimpleMatrix r = new SimpleMatrix(2,2); - - // TODO: set phi - double phi = 0; - - // TODO: fill the rotation matrix - // TODO - // TODO - // TODO - // TODO - - // TODO: define translation - SimpleVector t = new SimpleVector(0,0); - - // Initialize transformed phantom - Grid2D transformedPhantom = new Grid2D(phantom); - Grid2D transformedPhantom_blurred = new Grid2D(phantom_blurred); - - // Create the affine transformation - AffineTransform affine = new AffineTransform(r, t); - - // Apply the transformation - transformedPhantom.applyTransform(affine); - transformedPhantom_blurred.applyTransform(affine); - - transformedPhantom.show("Transformed Phantom"); - - - ///////////////////////////////////// - // Part 2: Find the transformation // - ///////////////////////////////////// - - // Registration of the transformed image to the initial phantom - Registration3 reg3 = new Registration3(); - reg3.reference = phantom_blurred; - reg3.image = transformedPhantom_blurred; - - // Optimization - double[] res = reg3.performOptimization(); - - // Stack for visualization purposes only - Grid3D optimizationStepsGrid = new Grid3D(reg3.reference.getWidth(),reg3.reference.getHeight(),reg3.iteration,false); - for (int i = 0; i < optimizationStepsGrid.getSize()[2]; i++) { - optimizationStepsGrid.setSubGrid(i, reg3.movingStack.getSubGrid(i)); - } - optimizationStepsGrid.show(); - - // Transform image back - Grid2D backtransformedImage = new Grid2D(transformedPhantom); - phi = (2*Math.PI)/360 * res[0]; - r.setElementValue(0, 0, Math.cos(phi)); - r.setElementValue(0, 1, - Math.sin(phi)); - r.setElementValue(1, 0, Math.sin(phi)); - r.setElementValue(1, 1, Math.cos(phi)); - - SimpleVector t2 = new SimpleVector(res[1],res[2]); - AffineTransform affine2 = new AffineTransform(r, t2); - backtransformedImage.applyTransform(affine2); - backtransformedImage.show(); - } -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/SVDandFT.java b/src/edu/stanford/rsl/tutorial/dmip/SVDandFT.java deleted file mode 100644 index 49dedd8c..00000000 --- a/src/edu/stanford/rsl/tutorial/dmip/SVDandFT.java +++ /dev/null @@ -1,431 +0,0 @@ -package edu.stanford.rsl.tutorial.dmip; - -import edu.stanford.rsl.conrad.data.numeric.Grid2D; -import edu.stanford.rsl.conrad.data.numeric.Grid2DComplex; -import edu.stanford.rsl.conrad.data.numeric.Grid3D; -import edu.stanford.rsl.conrad.fitting.LinearFunction; -import edu.stanford.rsl.conrad.numerics.DecompositionSVD; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType; -import edu.stanford.rsl.conrad.numerics.SimpleMatrix.MatrixNormType; -import edu.stanford.rsl.conrad.numerics.SimpleOperators; -import edu.stanford.rsl.conrad.numerics.SimpleVector; -import edu.stanford.rsl.conrad.utils.ImageUtil; -import edu.stanford.rsl.conrad.utils.VisualizationUtil; -import ij.IJ; -import ij.ImageJ; - -/** - * - * Exercise 2 of Diagnostic Medical Image Processing (DMIP) - * @author Marco Boegel - * - */ -public class SVDandFT { - /* - public static void invertSVD(SimpleMatrix A) - { - - - System.out.println("A = " + A.toString()); - - //Compute the inverse of A without using inverse() - //TODO - - //Check output: re-compute A = U * S * V^T - SimpleMatrix temp = SimpleOperators.multiplyMatrixProd(svd.getU(), svd.getS()); - SimpleMatrix A2 = SimpleOperators.multiplyMatrixProd(temp, svd.getV().transposed()); - System.out.println("U * S * V^T: " + A2.toString()); - - //Moore-Penrose Pseudoinverse defined as V * Sinv * U^T - //Compute the inverse of the singular matrix S - SimpleMatrix Sinv = new SimpleMatrix( svd.getS().getRows(), svd.getS().getCols()); - Sinv.zeros(); - - int size = Math.min(Sinv.getCols(), Sinv.getRows()); - SimpleVector SinvDiag = new SimpleVector( size); - - //TODO - //TODO - //TODO - - Sinv.setDiagValue(SinvDiag); - - //Compare our implementation to svd.getreciprocalS() - System.out.println("Sinv = " + Sinv.toString()); - System.out.println("Srec = " +svd.getreciprocalS().toString()); - - - SimpleMatrix tempInv = SimpleOperators.multiplyMatrixProd(svd.getV(), svd.getreciprocalS()); - SimpleMatrix Ainv = SimpleOperators.multiplyMatrixProd(tempInv, svd.getU().transposed()); - System.out.println("Ainv = " + Ainv.toString()); - System.out.println("A.inverse() = " + A.inverse(InversionType.INVERT_SVD)); - - //Condition number - //TODO - System.out.println("Cond(A) = " + cond); - - //introduce a rank deficiency - double eps = 10e-3; - SimpleMatrix Slowrank = new SimpleMatrix(svd.getS()); - int sInd = Math.min(svd.getS().getRows(), svd.getS().getCols()); - for(int i = 0; i < sInd; i++) - { - double val = svd.getS().getElement(i, i); - //TODO - //TODO - //TODO - } - - SimpleMatrix templowrank = SimpleOperators.multiplyMatrixProd(svd.getU(), Slowrank); - SimpleMatrix Alowrank = SimpleOperators.multiplyMatrixProd(templowrank, svd.getV().transposed()); - System.out.println("A rank deficient = " + Alowrank.toString()); - - - - //Show that a change in a vector b of only 0.1% can lead to 240% change in the result of A*x=b - //Consider A as already defined - //And vector b - SimpleVector b = new SimpleVector(1.001, 0.999, 1.001); - - //solve for x - SimpleVector x = SimpleOperators.multiply(Ainv, b); - System.out.println(x.toString()); - - //if we set vector b to the vector consisting of 1's, we imply a 0.1% change - SimpleVector br = new SimpleVector(1,1,1); - - - SimpleVector xr = SimpleOperators.multiply(Ainv, br); - System.out.println(xr.toString()); - //We want only the difference caused by the change - SimpleVector xn = SimpleOperators.subtract(xr, x); - - //Alternatively: - //SimpleVector bn = SimpleOperators.subtract(br, b); - //SimpleVector xn = SimpleOperators.multiply(Ainv, bn); - - System.out.println("Modification of x " + SimpleOperators.divideElementWise(xn, x).toString()); - } - - public static void optimizationProblem1(SimpleMatrix A, int rankDeficiency) - { - System.out.println("Optimization Problem 1"); - //%%%%%%%%%%%%%%%%%%%%%%%%%%%% - // Optimization problem I - //%%%%%%%%%%%%%%%%%%%%%%%%%%% - - // Let us consider the following problem that appears in many image - // processing and computer vision problems: - // We computed a matrix A out of sensor data like an image. By theory - // the matrix A must have the singular values s1, s2, . . . , sp, where - // p = min(m, n). Of course, in practice A does not fulfill this constraint. - // Problem: What is the matrix A0 that is closest to A (according to the - // Frobenius norm) and has the required singular values? - - - // Let us assume that by theoretical arguments the matrix A is required - // to have a rank deficiency of one, and the two non-zero singular values - // are identical. The matrix A0 that is closest to A according to the - // Frobenius norm and fulfills the above requirements is: - - // Build mean of first two singular values to make them identical - // Set third singular value to zero - - DecompositionSVD svd = new DecompositionSVD(A); - int rank = svd.rank(); - - int newRank = rank - rankDeficiency; - - - - //Compute mean of remaining singular values - double mean = 0; - //TODO - //TODO - //TODO - - //Create new Singular matrix - SimpleMatrix Slowrank = new SimpleMatrix(svd.getS().getRows(), svd.getS().getCols()); - Slowrank.zeros(); - //Fill in remaining singular values with the mean. - //TODO - //TODO - //TODO - - //compute A0 - SimpleMatrix temp = SimpleOperators.multiplyMatrixProd(svd.getU(), Slowrank); - SimpleMatrix A0 = SimpleOperators.multiplyMatrixProd(temp, svd.getV().transposed()); - - System.out.println("A0 = " + A0.toString()); - - double normA = A.norm(MatrixNormType.MAT_NORM_FROBENIUS); - double normA0 = A0.norm(MatrixNormType.MAT_NORM_FROBENIUS); - - System.out.println("||A||_F = " + normA); - System.out.println("||A0||_F = " + normA0); - - - } - - public static void optimizationProblem4(double[] xCoords, double[] yCoords) - { - System.out.println("Optimization Problem 4"); - //%%%%%%%%%%%%%%%%%%%%%%%%%%% - //Optimization problem IV - //%%%%%%%%%%%%%%%%%%%%%%%%%%% - - // The Moore-Penrose pseudo-inverse is required to find the - // solution to the following optimization problem: - // Compute the regression line through the following 2-D points: - // We have sample points (x_i,y_i) and want to fit a line model - // y = mx + t (linear approximation) through the points - // y = A*b - // Note: y = m*x +t*1 Thus: A contains (x_i 1) in each row, and y contains the y_i coordinates - // We need to solve for b = (m t) - - SimpleMatrix A = new SimpleMatrix(xCoords.length, 2); - SimpleVector aCol = new SimpleVector(xCoords.length); - SimpleVector y = new SimpleVector(xCoords.length); - for(int i = 0; i < xCoords.length; i++) - { - aCol.setElementValue(i, xCoords[i]); - y.setElementValue(i, yCoords[i]); - } - A.setColValue(0, aCol); - SimpleVector aColOnes = new SimpleVector(xCoords.length); - aColOnes.ones(); - A.setColValue(1, aColOnes); - - - - // Note: We need to use SVD matrix inversion because the matrix is not square - // more samples (7) than unknowns (2) -> overdetermined system -> least-square - // solution. - - //get solution for b - //TODO - //TODO - - - LinearFunction lFunc = new LinearFunction(); - lFunc.setM(b.getElement(0)); - lFunc.setT(b.getElement(1)); - VisualizationUtil.createScatterPlot(xCoords, yCoords, lFunc).show(); - - } - - public static void optimizationProblem2(SimpleMatrix b) - { - System.out.println("Optimization Problem 2"); - // Estimate the matrix A in R^{2,2} such that for the following vectors - // the optimization problem gets solved: - // sum_{i=1}^{4} b'_{i} A b_{i} and ||A||_{F} = 1 - // The objective function is linear in the components of A, thus the whole - // sum can be rewritten in matrix notation: - // Ma = 0 - // where the measurement matrix M is built from single elements of the - // sum. - // Given the four points we get four equations and therefore four rows in M: - - SimpleMatrix M = new SimpleMatrix( b.getCols(), 4); - - for(int i = 0; i < b.getCols(); i++) - { - // i-th column of b contains a vector. First entry of M is x^2 - //Setup measurement matrix M - //TODO - //TODO - //TODO - //TODO - } - - // TASK: estimate the matrix A - // HINT: Nullspace - DecompositionSVD svd = new DecompositionSVD(M); - //TODO - //TODO - - //We need to reshape the vector back to the desired 2x2 matrix form - SimpleMatrix A = new SimpleMatrix(2,2); - //TODO - //TODO - - //check if Frobenius norm is 1.0 - double normF = A.norm(MatrixNormType.MAT_NORM_FROBENIUS); - System.out.println("||A||_F = " + normF); - - //check solution - - SimpleVector temp = SimpleOperators.multiply(M, a); - double result = SimpleOperators.multiplyInnerProd(a, temp); - System.out.println("Minimized error: " + result); - } - - public static void optimizationProblem3(Grid2D image, int rank) - { - //% - //%%%%%%%%%%%%%%%%%%%%%%%%%%% - // Optimization problem III - //%%%%%%%%%%%%%%%%%%%%%%%%%%% - - // The SVD can be used to compute the image matrix of rank k that best - // approximates an image. Figure 1 shows an example of an image I - // and its rank-1-approximation I0 = u_{1} * s_{1} * v'_{1}. - - //In order to apply the svd, we first need to transfer our Grid2D image to a matrix - SimpleMatrix I = new SimpleMatrix(image.getHeight(), image.getWidth()); - - //NOTE: indices of matrix and image are reversed - for(int i = 0; i < image.getHeight(); i++) - { - for(int j = 0; j < image.getWidth(); j++) - { - double val = image.getAtIndex(j, i); - I.setElementValue(i, j, val); - } - } - - DecompositionSVD svd = new DecompositionSVD(I); - - //output images - Grid3D imageRanks = new Grid3D(image.getWidth(), image.getHeight(), rank); - - //Create Rank k approximations - for(int k = 0; k < rank; k++) - { - //TODO - //TODO - - - //Transfer back to grid - Grid2D imageRank = new Grid2D(image.getWidth(),image.getHeight()); - for(int i = 0; i < image.getHeight(); i++) - { - for(int j = 0; j < image.getWidth(); j++) - { - if(k == 0) - { - imageRank.setAtIndex(j, i, (float) Iapprox.getElement(i, j)); - - } - else - { - //TODO - } - - - } - } - imageRanks.setSubGrid(k, imageRank); - } - - imageRanks.show(); - - - //Direct estimation of rank K - //TODO - //TODO - //Transfer back to grid - Grid2D imageRankK = new Grid2D(image.getWidth(), image.getHeight()); - for(int i = 0; i < image.getHeight(); i++) - { - for(int j = 0; j < image.getWidth(); j++) - { - imageRankK.setAtIndex(j, i, (float) IapproxK.getElement(i, j)); - } - } - imageRankK.show(); - - - } - - public static void fourierExercise(Grid2D image) - { - //TODO complex image - // Important: Grid2DComplex enlarges the original image to the next power of 2 - imageC.show(); - - - //Apply 2-D discrete fourier transform - //Puts the DC component of the signal in the upper left corner of the FFT - //TODO - imageC.show("Shepp-Logan FFT"); - - //TODO - imageC.show("Shepp-Logan FFTShift"); - - //Visualize log transformed FFT log(1+|FFTshift(image)|) - Grid2D logFFT = new Grid2D(imageC.getMagnSubGrid(0, 0, imageC.getWidth(), imageC.getHeight())); - logFFT.getGridOperator().addBy(logFFT, 1.f); - logFFT.getGridOperator().log(logFFT); - logFFT.show("Shepp-Logan FFTShift log"); - - // Important: Grid2DComplex enlarges the original image to the next power of 2 - // When transforming back, make sure to prune the image size to your original image - //TODO inverse FFT - //TODO prune - imageTransf.show(); - } - - public static void main(String[] args) { - - ImageJ ij = new ImageJ(); - //Create matrix A - // 11 10 14 - // 12 11 -13 - // 14 13 -66 - SimpleMatrix A = new SimpleMatrix(3,3); - A.setRowValue(0, new SimpleVector(11, 10, 14)); - A.setRowValue(1, new SimpleVector(12, 11, -13)); - A.setRowValue(2, new SimpleVector(14, 13, -66)); - - invertSVD(A); - optimizationProblem1(A, 1); - - //Data for problem 4 from lecture slides - double[] xCoords = new double[]{3.f, 2.f, 1.f, 0.f, -1.f, -1.f, -2.f}; - double[] yCoords = new double[]{2.f, 1.f, 2.f, 0.f, 1.f, -1.f, -1.f}; - - optimizationProblem4(xCoords, yCoords); - - - - // Data for problem 2 from lecture slides - SimpleMatrix vectors = new SimpleMatrix(2,4); - vectors.setColValue(0, new SimpleVector(1.f, 1.f)); - vectors.setColValue(1, new SimpleVector(-1.f, 2.f)); - vectors.setColValue(2, new SimpleVector(1.f, -3.f)); - vectors.setColValue(3, new SimpleVector(-1.f, -4.f)); - - optimizationProblem2(vectors); - - //Load an image from file - String filename = "D:/04_lectures/DMIP/exercises/2014/1/yu_fill.jpg"; - Grid2D image = ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); - image.show(); - - int rank = 100; - optimizationProblem3(image, rank); - - - //Load Data for Fourier Transform exercise - //1. Start the ReconstructionPipelineFrame (src/apps/gui/) - //2. In the Pipeline window, go to edit configuration and press "center volume" and save - //3. In the ImageJ window, navigate to Plugins->CONRAD->Create Numerical Phantom - //4. Choose Metric Volume Phantom - //5. Choose Shepp Logan Phantom - //6. Save the resulting volume. In the ImageJ window, File-Save As->Tiff... - - String filenameShepp = "D:/04_lectures/DMIP/exercises/2014/1/shepplogan.tif"; - Grid3D sheppLoganVolume = ImageUtil.wrapImagePlus(IJ.openImage(filenameShepp)); - //To work with a 2-D image, select slice 160 - Grid2D sheppLoganImage = sheppLoganVolume.getSubGrid(160); - sheppLoganImage.show(); - fourierExercise(sheppLoganImage); - - - - } -*/ -} diff --git a/src/edu/stanford/rsl/tutorial/dmip/yu_fill.jpg b/src/edu/stanford/rsl/tutorial/dmip/yu_fill.jpg deleted file mode 100644 index 5c19b407..00000000 Binary files a/src/edu/stanford/rsl/tutorial/dmip/yu_fill.jpg and /dev/null differ diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseFB.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseFB.java new file mode 100644 index 00000000..64ff7b34 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseFB.java @@ -0,0 +1,461 @@ +package edu.stanford.rsl.tutorial.mipda; + +import edu.stanford.rsl.conrad.data.numeric.Grid1D; +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators; +import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; +import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND; +import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import edu.stanford.rsl.conrad.utils.ImageUtil; +import edu.stanford.rsl.tutorial.filters.RamLakKernel; +import edu.stanford.rsl.tutorial.filters.SheppLoganKernel; +import ij.IJ; +import ij.ImageJ; + +/** + * Short Scan for Fan Beam Reconstruction + * Programming exercise for module "Fan Beam Reconstruction" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Marco Boegel, Anna Gebhard, Mena Abdelmalek + * + */ + +public class ExerciseFB { + + public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN}; + + //TODO: choose the ramp filter (NONE, RAMLAK, SHEPPLOGAN) + final RampFilterType filter = RampFilterType.RAMLAK; + //TODO: choose the sinogram data (Sinogram0.tif, Sinogram1.tif, Sinogram2.tif) + final String filename = "Sinogram0.tif"; + + // system parameters + final float focalLength = 600.f; // source to detector distance in [mm] + // + final int detectorPixels;// number of pixels on the detector + final float detectorSpacing = 0.5f;// spacing on the detector in [mm] per pixel + float detectorLength;// detector length in [mm] + // + final int numProjs;// number of acquired projections + float betaIncrement;// increment of rotation angle beta between two projections in [rad] + float maxBeta;// max rotation angle beta in [rad] + float halfFanAngle;// half fan angle in [rad] + + // parameters for reconstruction + int[] recoSize = new int[]{128, 128}; + float[] spacing = new float[]{1.0f, 1.0f}; + + // data containers + Grid2D sino;// original sinogram + Grid2D sino_cosW;// sinogram + cosine filter + Grid2D sino_cosW_parkerW;// sinogram + cosine filter + Parker weights + Grid2D sino_cosW_parkerW_filt;// sinogram + cosine filter + Parker weights + NONE/RAMLAK/SHEPPLOGAN + Grid2D parker;// Parker weights + Grid2D reconstruction;// reconstructed image + + public static void main(String arg[]) { + + ImageJ ij = new ImageJ(); + + //TODO: Look for the method initParameters() and do the first TODO-s there! + ExerciseFB fbp = new ExerciseFB(); + + fbp.sino.show("Sinogram"); + + System.out.println("Half fan angle [deg]: " + fbp.get_halfFanAngle()*180.0/Math.PI); + System.out.println("Short-scan range [deg]: " + fbp.get_maxBeta()*180/Math.PI); + + fbp.doShortScanReconstruction(); //more TODO-s here! + } + + void initParameters() { + + detectorLength = detectorPixels*detectorSpacing; + + halfFanAngle = 1.0f; //TODO: determine the half fan angle such that the complete fan covers the detector exactly (hint: trigonometry) + maxBeta = 1.0f; //TODO: compute the short scan minimum range using the value of halfFanAngle + betaIncrement = maxBeta/(float) numProjs; + } + + void doShortScanReconstruction() { + + //assume a flat detector (equally-spaced case) + + //apply cosine weights + sino_cosW = applyCosineWeights(sino); //TODO: Go to this method and implement the missing steps! + sino_cosW.show("After Cosine Weighting"); + + //apply Parker redundancy weights + sino_cosW_parkerW = applyParkerWeights(sino_cosW); //TODO: Go to this method and implement the missing steps! + parker.show("Parker Weights"); + sino_cosW_parkerW.show("After Parker Weighting"); + + //apply ramp filter + sino_cosW_parkerW_filt = applyRampFiltering(sino_cosW_parkerW); + switch(filter) { + case RAMLAK: + sino_cosW_parkerW_filt.show("After RamLak Filter"); + break; + case SHEPPLOGAN: + sino_cosW_parkerW_filt.show("After Shepp-Logan Filter"); + break; + case NONE: // + default: // + } + + //perform backprojection to reconstruct the image + Grid2D reconstruction = backproject(sino_cosW_parkerW_filt,recoSize,spacing); //TODO: Go to this method and implement the missing steps! + reconstruction.show("Reconstructed Image"); + } + + /** + * Cosine weighting for 2D sinograms + * @param sino the sinogram + * @return the cosine weighted sinogram + */ + public Grid2D applyCosineWeights(Grid2D sino) { + + Grid2D result = new Grid2D(sino); + + //create 1D kernel (identical for each projection) + Grid1D cosineKernel = new Grid1D(detectorPixels); + for(int i=0; i < detectorPixels; ++i){ + + float t = 0.f; //TODO: compute center of the current detector pixel (zero is center of detector) + cosineKernel.setAtIndex(i, 0.0f); //TODO (hint: use t and focalLength) + } + + //apply cosine weights to each projection + for(int p = 0; p < numProjs; p++){ + NumericPointwiseOperators.multiplyBy(result.getSubGrid(p), cosineKernel); + } + + return result; + } + + /** + * Parker redundancy weights for a 2D sinogram + * @param sino the sinogram + * @return parker weighted sinogram + */ + public Grid2D applyParkerWeights(Grid2D sino) { + + Grid2D result = new Grid2D(sino); + parker = new Grid2D(sino.getWidth(), sino.getHeight()); + + float delta = halfFanAngle; + float beta, gamma; + + // iterate over the detector elements + for (int i = 0; i < detectorPixels; ++i) { + + // compute gamma of the current ray (detector element & flat detector) + gamma = (float) Math.atan((i*detectorSpacing - detectorLength/2.f + 0.5f*detectorSpacing)/focalLength); + + // iterate over the projection angles + for (int b = 0; b < numProjs; ++b) { + + beta = b*betaIncrement; + + if (beta > 2.f*Math.PI) { + beta -= 2.f*Math.PI; + } + + // implement the conditions as described in Parker's paper ( https://dx.doi.org/10.1118/1.595078 ) + if (beta <= 2*(delta - gamma)) { + + double tmp = 0.f; //TODO + float val = 1.f; //TODO + + if (Double.isNaN(val)){ + continue; + } + + parker.setAtIndex(i, b, val); + } + else if (beta <= Math.PI - 2.f*gamma) { + parker.setAtIndex(i, b, 1); + } + else if (beta <= (Math.PI + 2.f*delta) + 1e-12) { + + double tmp = Math.PI/4.f * (Math.PI + 2.f*delta - beta)/(delta + gamma); // (error in the original paper) + float val = 1.f; //TODO + + if (Double.isNaN(val)){ + continue; + } + + parker.setAtIndex(i, b, val); + } + } + } + + // correct for scaling due to varying angle + NumericPointwiseOperators.multiplyBy(parker, (float)(maxBeta/(Math.PI))); + + //apply the Parker weights to the sinogram + NumericPointwiseOperators.multiplyBy(result, parker); + + return result; + } + + /** + * A pixel driven backprojection algorithm + * Cosine, redundancy weighting and ramp filters need to be applied separately beforehand! + * Remark: This implementation requires a detector whose center coincides with the center of rotation. + * @param sino the filtered sinogram + * @param recoSize the dimension of the output image + * @param spacing the spacing of the output image + * @return the reconstruction + */ + public Grid2D backproject(Grid2D sino, int[] recoSize, float[] spacing) { + + Grid2D result = new Grid2D(recoSize[0],recoSize[1]); // x and y flipped for ImageJ + result.setSpacing(spacing[0],spacing[1]); + + //remark: float values are used here to make a transfer to a GPU implementation easier + for(int p = 0; p < numProjs; p++) { + + //first, compute the rotation angle beta and pre-compute cos(beta), sin(beta) + final float beta = (float) (p*betaIncrement); + final float cosBeta = (float) Math.cos(beta); + final float sinBeta = (float) Math.sin(beta); + + //compute rotated source point (here: center of rotation coincides with the origin, beta==0 <=> source in top position) + final PointND source = new PointND(-focalLength*sinBeta, focalLength*cosBeta, 0.f); + + //compute the world coordinate point of the left detector edge (here: center of detector coincides with the origin) + float offset = -detectorLength/2.f; + final PointND detBorder = new PointND(0.f, 0.f, 0.f); //TODO + + //compute direction vector along the detector array (towards detector center) + SimpleVector dirDet = null; //TODO + if (dirDet == null) + dirDet = new SimpleVector(3); + dirDet.normalizeL2(); // ensure length == 1 + + final StraightLine detLine = new StraightLine(detBorder, dirDet); + detLine.setDirection(detLine.getDirection().normalizedL2()); // bugfix (to be sure) + + //pick current projection + Grid1D currProj = sino.getSubGrid(p); + + float wx, wy; + //pixel driven BP: iterate over all output image pixels + for(int x = 0; x < recoSize[0]; x++) { + + //transform the image pixel coordinates to world coordinates + wx = getWorldCoordinate(x, recoSize[0], spacing[0]); // TODO-s here! + + for(int y = 0; y < recoSize[1]; y++) { + + //transform the image pixel coordinates to world coordinates + wy = getWorldCoordinate(y, recoSize[1], spacing[1]); // TODO-s here! (in case you missed it above) + + final PointND reconstructionPointWorld = new PointND(wx, wy, 0.f); + + //intersect the projection ray with the detector + final StraightLine projectionLine = null; //TODO + if (projectionLine != null) + projectionLine.setDirection(projectionLine.getDirection().normalizedL2()); // bugfix + final PointND detPixel = null; //TODO + + float valueTemp; + if(detPixel != null) { + + //calculate position of the projected point on the 1D detector + float t = (float) SimpleOperators.multiplyInnerProd(detPixel.getAbstractVector(), dirDet) - offset; + + //check if projection of this world point hits the detector + if(t < 0 || t > detectorLength) + continue; + + float value = InterpolationOperators.interpolateLinear(currProj, t/detectorSpacing - 0.5f); + + //apply distance weighting + float dWeight = distanceWeight(reconstructionPointWorld, beta); // TODO-s here! + valueTemp = (float) (value/(dWeight*dWeight)); + } + else { + + //projected pixel lies outside of the detector + valueTemp = 0.f; + } + + result.addAtIndex(x, recoSize[1]-y-1, valueTemp); // flip to show top-down in ImageJ + } + } + } + + //adjust scale, required because of short scan + float normalizationFactor = (float) (numProjs/Math.PI); + NumericPointwiseOperators.divideBy(result, normalizationFactor); + + return result; + } + + /** + * Transformation from image coordinates to world coordinates (single axis) + * @param imCoordinate image coordinate + * @param dim the dimension of the image + * @param spacing the spacing of the image + * @return world coordinate + */ + public float getWorldCoordinate(int imCoordinate, int dim, float spacing) { + + float wCoordinate = 0.f; //TODO + + return wCoordinate; + } + + /** + * Compute distance weight (refer to the course slides) + * @param reconstructionPointWorld reconstruction point in world coordinates + * @param beta rotation angle beta + * @return distance weight + */ + public float distanceWeight(PointND reconstructionPointWorld, float beta) { + + //Compute distance weight + float radius = 0.f; //TODO + float phi = 0.f; //TODO + float U = 0.f; //TODO + + return U; + } + + + /** + * + * end of the exercise + */ + + public ExerciseFB() { + + //Load and visualize the projection image data + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + + sino = readImageFile(imageDataLoc + filename); + + // transpose so that each row is a projection + sino = sino.getGridOperator().transpose(sino); // warning "Error in transpose" can be ignored (too strict tolerance settings) + + detectorPixels = sino.getWidth(); + numProjs = sino.getHeight(); + + initParameters(); + + sino.setSpacing(detectorSpacing, betaIncrement); + } + + public Grid2D readImageFile(final String filename){ + return ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); + } + + private Grid2D applyRampFiltering(Grid2D sino) { + + Grid2D result = new Grid2D(sino); + + switch(filter) { + + case RAMLAK: + + RamLakKernel ramLak = new RamLakKernel(detectorPixels, detectorSpacing); + for(int i = 0; i < numProjs; i++) { + ramLak.applyToGrid(result.getSubGrid(i)); + } + + break; + + case SHEPPLOGAN: + + SheppLoganKernel sheppLogan = new SheppLoganKernel(detectorPixels, detectorSpacing); + for(int i = 0; i < numProjs; i++) { + sheppLogan.applyToGrid(result.getSubGrid(i)); + } + + break; + + case NONE: // + default: // + } + + return result; + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public RampFilterType get_filter() { + return filter; + } + // + public String get_filename() { + return filename; + } + // + public float get_focalLength() { + return focalLength; + } + // + public int get_detectorPixels() { + return detectorPixels; + } + // + public float get_detectorSpacing() { + return detectorSpacing; + } + // + public float get_detectorLength() { + return detectorLength; + } + // + public int get_numProjs() { + return numProjs; + } + // + public float get_betaIncrement() { + return betaIncrement; + } + // + public float get_maxBeta() { + return maxBeta; + } + // + public float get_halfFanAngle() { + return halfFanAngle; + } + // + public int[] get_recoSize() { + return recoSize; + } + // + public float[] get_spacing() { + return spacing; + } + // + public Grid2D get_sino() { + return sino; + } + // + public Grid2D get_sino_cosW() { + return sino_cosW; + } + // + public Grid2D get_sino_cosW_parkerW() { + return sino_cosW_parkerW; + } + // + public Grid2D get_sino_cosW_parkerW_filt() { + return sino_cosW_parkerW_filt; + } + // + public Grid2D get_parker() { + return parker; + } + // + public Grid2D get_reconstruction() { + return reconstruction; + } +} diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseGeoU.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseGeoU.java new file mode 100644 index 00000000..e1cd32eb --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseGeoU.java @@ -0,0 +1,562 @@ +package edu.stanford.rsl.tutorial.mipda; + +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; +import edu.stanford.rsl.conrad.utils.ImageUtil; +//import edu.stanford.rsl.conrad.utils.VisualizationUtil; +import ij.IJ; +import ij.ImageJ; +import edu.stanford.rsl.conrad.data.numeric.InterpolationOperators; +//import edu.stanford.rsl.conrad.data.numeric.NumericGridOperator; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import edu.stanford.rsl.conrad.numerics.DecompositionSVD; + +/** + * Geometric Undistortion + * Programming exercise 1 for module "Defect Pixel Interpolation" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Ashwini Jadhav, Anna Gebhard, Mena Abdelmalek + * + */ + +public class ExerciseGeoU { + + final static int caseNo = 0; // choose 0, 1, or 2 for different test images + + final static float a = 0.6f; // a: spread among x-direction + final static float b = 0.3f; // b: spread among y-direction + final static float d = 5.0f; // d: d/2 value of distortion field at level set for (X/a)^2+(Y/b)^2=1 + + String filename; // do not edit, choose a number in the main method instead + Grid2D originalImage; // do not edit, read from filename when instantiating the object + Grid2D distortedImage; // do not edit, generated in the main method + Grid2D undistortedImage; // do not edit, member variable for the output image + + // number of lattice points + final int nx = 0;//TODO: define the number of lattice point: nx (Positive number and less than 20) + final int ny = 0;//TODO: define the number of lattice point: ny (Positive number and less than 20) + + float fx; + float fy; + SimpleMatrix Xu; + SimpleMatrix Yu; + SimpleMatrix Xd; + SimpleMatrix Yd; + int numCoeff; + int numCorresp; + SimpleMatrix A; + DecompositionSVD svd; + SimpleMatrix A_pseudoinverse; + SimpleVector u_vec; + SimpleVector v_vec; + SimpleVector XuVector; + SimpleVector YuVector; + SimpleVector XdVector; + SimpleVector YdVector; + Grid2D xDist; + Grid2D yDist; + Grid2D xprime; // x-coordinates of point correspondences in the undistorted image + Grid2D yprime; // y-coordinates of point correspondences in the undistorted image + Grid2D x; // x-coordinates of point correspondences in the distorted image + Grid2D y; // y-coordinates of point correspondences in the distorted image + + + /** + * main method + * Here you can choose the image, and set the distortion parameters + */ + public static void main(String[] args) { + + new ImageJ(); + + // generate distorted image + ExerciseGeoU exObj = new ExerciseGeoU(caseNo); + exObj.init(a,b,d); + + // TASK: go to the method doImageUndistortion(...) and complete it + exObj.undistortedImage = exObj.doImageUndistortion(exObj.distortedImage); + + exObj.originalImage.show("Original Image"); + exObj.distortedImage.show("Distorted Image"); + exObj.undistortedImage.show("Undistorted Image"); + + Grid2D differenceImage = (Grid2D) NumericPointwiseOperators.subtractedBy( + exObj.originalImage, exObj.undistortedImage); + + differenceImage.show("Difference Original vs. Undistorted"); + + Grid2D distortedVSundistorted = (Grid2D) NumericPointwiseOperators.subtractedBy( + exObj.undistortedImage, exObj.distortedImage); + + distortedVSundistorted.show("Difference Undistorted vs. Distorted"); + } + + public Grid2D doImageUndistortion(Grid2D distortedImage){ + + Grid2D grid = new Grid2D(distortedImage); + + // check if image distortion routine has been run + if (distortedImage == null || undistortedImage == null + || xprime == null || yprime == null + || x == null || y == null) { + + System.err.println("Error in doImageUndistortion(): called before initialization of members."); + return grid; + } + + + getLatticePoints(distortedImage);// There are TODOs here + + + int degree = 5; // polynomial's degree: 2,...,10 + calcPolynomial(degree, Xd);// There are TODOs here + + + + computeMatrixA(degree, numCorresp, numCoeff);// There are TODOs here + + + computeDistortionCoeff(A, XuVector, YuVector);// There are TODOs here + + + grid = computeDistortedGrid(distortedImage, grid, degree);// There are TODOs here + + return grid; + } + + + + /** 1. Number of lattice points (this only works for symmetric images): + * nx, ny feature points: usually provided by tracking point + * correspondences in the phantom during the calibration step. + * Here, we use the distorted and undistorted coordinates generated + * in the course of the distortion simulation. + */ + public void getLatticePoints (Grid2D distortedImage){ + + + int imageWidth = distortedImage.getWidth(); + int imageHeight = distortedImage.getHeight(); + + // step size + // calculate the step size of the lattice points: fx and fy + fx = 0; //TODO + fy = 0; //TODO + + // fill the distorted and undistorted lattice points + // with data from the given correspondences + // matrix = number of rows (y) x number of columns (x) + Xu = new SimpleMatrix(ny,nx); + Yu = new SimpleMatrix(ny,nx); + Xd = new SimpleMatrix(ny,nx); + Yd = new SimpleMatrix(ny,nx); + + for(int i = 0; i < nx; i++){ + for(int j = 0; j < ny; j++){ + + // sample the distorted and undistorted grid points at the lattice points + //TODO: fill matrix Xu + //TODO: fill matrix Yu + //TODO: fill matrix Xd + //TODO: fill matrix Yd + } + } + } + + /** 2. Polynomial of degree d + * Polynomial of degree d -> (d-1): extrema (e.g., d=5: 4 extrema) + * d=0: constant (horizontal line with y-intercept a_0 -> f(x)=a_0) + * d=1: oblique line with y-intercept a_0 & slope a_1 -> f(x)=a_0 + a_1 x + * d=2: parabola + * d>=2: continuous non-linear curve + * ... + * d = 10 -> NumCoeff: 66 -> but only 64 lattice points are known (nx, ny = 8) + */ + + public void calcPolynomial(int degree, SimpleMatrix Xd){ + + // number of coefficients: numCoeff + // (hint: this is NOT the total number of multiplications!) + numCoeff = 0; //TODO + + // number of correspondences: numCorresp + numCorresp = 0; //TODO + + // Printout of the used parameters + System.out.println("Polynom of degree: " + degree); + System.out.println("Number of Coefficients: " + numCoeff); + System.out.println("Number of Correspondences: " + numCorresp); + + } + + /** + * 3. Create the matrix A + */ + + public void computeMatrixA(int degree, int numCorresp, int numCoeff){ + + A = new SimpleMatrix(numCorresp, numCoeff); + if (numCorresp <= 0 && numCoeff <= 0) + return; + A.zeros(); + + // Realign the grid matrix into a vector for easier access in the next step + XuVector = new SimpleVector(numCorresp); + YuVector = new SimpleVector(numCorresp); + XdVector = new SimpleVector(numCorresp); + YdVector = new SimpleVector(numCorresp); + + for(int i = 0; i < nx; i++){ + for(int j = 0; j < ny; j++){ + + XuVector.setElementValue(i*ny + j, Xu.getElement(j, i)); + YuVector.setElementValue(i*ny + j, Yu.getElement(j, i)); + XdVector.setElementValue(i*ny + j, Xd.getElement(j, i)); + YdVector.setElementValue(i*ny + j, Yd.getElement(j, i)); + } + } + + // Compute matrix A (coordinates from distorted image) + for(int r = 0; r < numCorresp; r++){ + + int cc = 0; + + for(int k = 0; k <= degree; k++){ + for(int l = 0; l <= (degree-k); l++){ + + // TODO: fill matrix A + cc++; + + } + } + } + } + + public void computeDistortionCoeff(SimpleMatrix A, SimpleVector XuVector, SimpleVector YuVector){ + + // Compute the pseudo-inverse of A with the help of the SVD (class: DecompositionSVD) + svd = null; // TODO + A_pseudoinverse = null; // TODO + + // Compute the distortion coefficients (solve for known corresponding undistorted points) + u_vec = null;// TODO + v_vec = null;// TODO + } + + /** + * 4. Compute the distorted grid points (xDist, yDist) which are used to sample the + * distorted image to get the undistorted image + * (xprime,yprime) is the position in the undistorted image + * (xDist,yDist) is the position in the distorted (observed) X-ray image. + */ + public Grid2D computeDistortedGrid(Grid2D distortedImage, Grid2D grid_out, int degree){ + + int imageWidth = distortedImage.getWidth(); + int imageHeight = distortedImage.getHeight(); + + xDist = new Grid2D(imageWidth,imageHeight); + yDist = new Grid2D(imageWidth,imageHeight); + + for(int i = 0; i < imageWidth; i++){ + for(int j = 0; j < imageHeight; j++){ + + float val1, val2; //variables Val1 and Val2 are used for intermediate computation of xDist and yDist respectively + + int cc = 0; + + for(int k = 0; k <= degree; k++){ + for(int l = 0; l <= degree - k; l++){ + + val1 = 0;// TODO + val2 = 0;// TODO + + // TODO: fill xDist + // TODO: fill yDist + + cc++; + } + } + } + } + + float val; // variable for intermediate computations of undistorted image + + for(int i = 0; i < imageWidth; i++){ + for(int j = 0; j < imageHeight; j++){ + + // hint: consider the fact that the coordinate origin is in the center of the image + val = 0;//TODO + //TODO: fill grid_out + } + } + + return grid_out; + + } + + + + /** + * + * end of the exercise + */ + + + /** + * standard constructor + */ + public ExerciseGeoU(){ + + this(0); + } + + /** + * parameterized constructor for the ExerciseDPI object + * @param caseNo: select one of three images to test your undistortion method + */ + public ExerciseGeoU(int caseNo){ + + switch (caseNo) { + case 0: + filename = "frame32.jpg"; + break; + case 1: + filename = "undistorted.jpg"; + break; + case 2: + filename = "frame90.jpg"; + break; + default: + filename = "frame32.jpg"; + break; + } + + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + + originalImage = ImageUtil.wrapImagePlus(IJ.openImage(imageDataLoc + filename)).getSubGrid(0); + } + + public void init(float a, float b, float d) { + + // Normalize intensity values to [0,1] + float max = NumericPointwiseOperators.max(originalImage); + float min = NumericPointwiseOperators.min(originalImage); + + for(int i = 0; i < originalImage.getWidth(); i++){ + for(int j = 0; j < originalImage.getHeight(); j++){ + + float scaledValue = (originalImage.getAtIndex(i, j) - min) / (max-min); + originalImage.setAtIndex(i, j, scaledValue); + } + } + + undistortedImage = new Grid2D(originalImage.getWidth(),originalImage.getHeight()); // initialization with zeros + + distortedImage = generateDistortedImage(originalImage,a,b,d); + } + + /** + * method to generate a distorted image + * An artificial distortion field is generated. With this field, + * a distorted image is generated from the original image. + * Remark: coordinates of the latest distortion simulated are saved in class member variables + * @param inputImage: undistorted image on which distortion is simulated + * @param a: elliptic extent of the distortion in x-direction + * @param b: elliptic extent of the distortion in y-direction + * @param d: strength of the distortion + * @return: distorted image + */ + + public Grid2D generateDistortedImage(Grid2D inputImage, float a, float b, float d){ + + Grid2D grid = new Grid2D(inputImage); + int imageWidth = grid.getWidth(); + int imageHeight = grid.getHeight(); + + // generate grids to sample coordinates X and Y of the original image + Grid2D X = new Grid2D(imageWidth,imageHeight); // x-coordinates with respect to origin in image center, spacing = 1 [unit] + Grid2D Y = new Grid2D(imageWidth,imageHeight); // y-coordinates with respect to origin in image center, spacing = 1 [unit] + + float halfWidth = imageWidth / 2.0f; + float halfHeight = imageHeight / 2.0f; + + for(int i = 0; i < X.getWidth(); i++){ + for(int j = 0; j < X.getHeight(); j++){ + + X.setAtIndex(i, j, i - halfWidth + 0.5f); // assign pixel centers in x-direction + Y.setAtIndex(i, j, j - halfHeight + 0.5f); // assign pixel center in y-direction + } + } + + // create an artificial elliptical distortion field R (additive field) + Grid2D R = new Grid2D(imageWidth,imageHeight); + + for(int i = 0; i < R.getWidth(); i++){ + for(int j = 0; j < R.getHeight(); j++){ + + // distortion gets stronger with the distance from the image center + float r = d * (float)(Math.sqrt( + Math.pow(X.getAtIndex(i, j)/(a*halfWidth),2) + + Math.pow(Y.getAtIndex(i, j)/(b*halfHeight),2))); + R.setAtIndex(i, j, r); + } + } + + // shifting the positive range to half positive/half negative range + float maxR = NumericPointwiseOperators.max(R); + NumericPointwiseOperators.subtractBy(R, maxR * 0.5f); + + // create the distorted image coordinates: + // distorted image coordinates = undistorted image points + artificial distortion field + + // distorted image coordinates Xd, Yd + Grid2D Xd = new Grid2D(X); + Grid2D Yd = new Grid2D(Y); + + NumericPointwiseOperators.addBy(Xd, R); + NumericPointwiseOperators.addBy(Yd, R); + + // create the distorted image: + // re-sample the original input image at the distorted image points (Xd,Yd) + // to create artificial distortion + for(int i = 0; i < imageWidth; i++){ + for(int j = 0; j < imageHeight; j++){ + + float distortedValue = InterpolationOperators.interpolateLinear( + inputImage, + Xd.getAtIndex(i, j) + halfWidth - 0.5f, // remember: coordinate origin was set to image center + Yd.getAtIndex(i, j) + halfHeight - 0.5f); + grid.setAtIndex(i, j, distortedValue); + } + } + + xprime = X; + yprime = Y; + x = Xd; + y = Yd; + + return grid; + } + + //getters for members + // variables which are checked (DO NOT CHANGE!) + + public static int get_caseNo() { + return caseNo; + } + + public Grid2D get_originalImage() { + return originalImage; + } + + public Grid2D get_distortedImage() { + return distortedImage; + } + + public Grid2D get_undistortedImage() { + return undistortedImage; + } + + public int get_nx() { + return nx; + } + + public int get_ny() { + return ny; + } + + public float get_fx() { + return fx; + } + + public float get_fy() { + return fy; + } + + public SimpleMatrix get_Xu() { + return Xu; + } + + public SimpleMatrix get_Yu() { + return Yu; + } + + public SimpleMatrix get_Xd() { + return Xd; + } + + public SimpleMatrix get_Yd() { + return Yd; + } + + public int get_numCoeff() { + return numCoeff; + } + + public int get_numCorresp() { + return numCorresp; + } + + public SimpleMatrix get_A() { + return A; + } + + public DecompositionSVD get_svd() { + return svd; + } + + public SimpleMatrix get_A_pseudoinverse() { + return A_pseudoinverse; + } + + public SimpleVector get_u_vec() { + return u_vec; + } + + public SimpleVector get_v_vec() { + return v_vec; + } + + public SimpleVector get_XuVector() { + return XuVector; + } + + public SimpleVector get_YuVector() { + return YuVector; + } + + public SimpleVector get_XdVector() { + return XdVector; + } + + public SimpleVector get_YdVector() { + return YdVector; + } + + public Grid2D get_xDist() { + return xDist; + } + + public Grid2D get_yDist() { + return yDist; + } + + public Grid2D get_xprime() { + return xprime; + } + + public Grid2D get_yprime() { + return yprime; + } + + public Grid2D get_x() { + return x; + } + + public Grid2D get_y() { + return y; + } + +} diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseMI.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseMI.java new file mode 100644 index 00000000..95c3508d --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseMI.java @@ -0,0 +1,463 @@ +package edu.stanford.rsl.tutorial.mipda; + +import java.util.ArrayList; +import java.util.Iterator; +import java.util.TreeMap; +import java.util.Map.Entry; + +import edu.stanford.rsl.conrad.data.numeric.Grid1D; +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.data.numeric.Grid3D; +import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; +import edu.stanford.rsl.conrad.geometry.transforms.AffineTransform; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix; +import edu.stanford.rsl.conrad.utils.ImageUtil; +import edu.stanford.rsl.conrad.utils.VisualizationUtil; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import edu.stanford.rsl.jpop.FunctionOptimizer; +import edu.stanford.rsl.jpop.OptimizableFunction; +import edu.stanford.rsl.jpop.OptimizationOutputFunction; +import edu.stanford.rsl.jpop.FunctionOptimizer.OptimizationMode; + +import ij.IJ; +import ij.ImageJ; + + +/** + * Using mutual information to solve the registration problem + * Programming exercise for chapter "Rigid Registration" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Bastian Bier, Anna Gebhard, Mena Abdelmalek, Ashwini Jadhav + * + */ + +public class ExerciseMI { + + final double blurSigma = 4.0; + final int iterLimit = 50; + + final Grid2D image1; + final Grid2D image2; + + public Grid2D reference_img; + public Grid2D moving_img; + + public Grid3D movingStack; + + + public static void main(String[] args){ + + ImageJ ij = new ImageJ(); + + ExerciseMI exObj = new ExerciseMI(); + + exObj.image1.show("Input Image 1"); + exObj.image2.show("Input Image 2"); + + // both images are blurred to avoid the optimizer running into local minima + // registering the blurred versions supposedly result in the same transformation parameters + exObj.reference_img.show("Reference image"); + exObj.moving_img.show("Moving image"); + + // initialize cost function + CostFunction costFunction = exObj.new CostFunction(); + + /* Hint for TODO-s: + * The cost function computes the difference of a certain measure for both images. + * You can look up the details in the respective nested class below if you want to, + * otherwise we spare you the details and point you right to the tasks. + * + * TASK: You have to implement the methods that are required to compute the similarity measure, + * which is mutual information here (more details are in the course description). + * These methods are + * - double calculateMutualInformation(Grid2D, Grid2D){} [lines 106-166] + * - SimpleMatrix calculateJointHistogram(Grid2D, Grid2D, int){} [lines 175-211] + * - SimpleVector getHistogramFromJointHistogram(Grid2D,boolean){} [lines 219-246] + */ + + // perform registration + double res[] = exObj.performRegistration(costFunction); + + // visualization of cost function value over time + VisualizationUtil.createPlot(costFunction.getCostPerIterGrid().getBuffer()).show(); + + // stack for visualization purposes only + if (exObj.movingStack.getNumberOfElements() > 0) + exObj.movingStack.show("Optimization Steps"); + + // transform image back + double phi = Math.PI/180*res[0]; + SimpleVector t = new SimpleVector(res[1], res[2]); + RigidTransform rigidTransform = exObj.new RigidTransform(phi,t); + + // show that the found transform registers the high-resolution moving image image2 well + Grid2D registeredImage = new Grid2D(exObj.image2); + registeredImage.applyTransform(rigidTransform); + registeredImage.show("Registered Image"); + } + + /** + * Method to calculate the mutual information (MI) + * @param ref reference image + * @param mov moving image + * @return negative mutual information + */ + // TODO-s here! + public double calculateMutualInformation(Grid2D ref, Grid2D mov){ + + final int histSize = 256; + + // Step 1: Calculate joint histogram + SimpleMatrix jointHistogram = calculateJointHistogram(ref, mov, histSize); + + // Step 2: Get histogram for a single image from the joint histogram + // a) for the first image + SimpleVector histo1 = new SimpleVector(histSize); + histo1 = getHistogramFromJointHistogram(jointHistogram,false); + + // b) for the second image + SimpleVector histo2 = new SimpleVector(histSize); + histo2 = getHistogramFromJointHistogram(jointHistogram,true); + + // Step 3: + // a) Calculate the joint entropy + double entropy_jointHisto = 0; + + for (int i = 0; i < histSize; i++) { + + for (int j = 0; j < histSize; j++){ + + if(jointHistogram.getElement(i, j) != 0) { + + // calculate entropy of the joint histogram (hint: use logarithm base 2 and use the correct sign) + entropy_jointHisto = 1;//TODO + } + } + } + + // b) Calculate the marginal entropies + double entropy_histo1 = 0; + double entropy_histo2 = 0; + + for(int i = 0; i < histSize; i++){ + + if(histo1.getElement(i) != 0){ + + // calculate entropy for histogram 1 (hint: use logarithm base 2 and use the correct sign) + entropy_histo1 = 1;//TODO + } + + if(histo2.getElement(i) != 0){ + + // calculate entropy for histogram 2 (hint: use logarithm base 2 and use the correct sign) + entropy_histo2 = 1;//TODO + } + } + + // Note: Make sure that you considered the correct sign in the entropy formulas! + + // Step 4: Calculate the mutual information + double mutual_information = entropy_histo1 + entropy_histo2 - entropy_jointHisto; + + // Mutual information is high for a good match. + // For the optimizer we require a minimization problem. + // The factor 100 serves as a trick to stabilize the optimization process. + return mutual_information *= -100; + } + + /** + * Method to calculate the joint histogram of two images + * @param im1 image1 + * @param im2 image2 + * @return a SimpleMatrix corresponding to the joint histogram + */ + // TODO-s here! + public SimpleMatrix calculateJointHistogram(Grid2D im1, Grid2D im2, int histSize){ + + SimpleMatrix jH = new SimpleMatrix(histSize, histSize); + + int imWidth = im1.getWidth(); + int imHeight = im1.getHeight(); + + if (imWidth != im2.getWidth() && imHeight != im2.getHeight()) { + + System.err.println("Image inputs have to have same size for joint histogram evaluation."); + return jH; + } + + float min1 = NumericPointwiseOperators.min(im1); + float scaleFactor1 = (histSize - 1)/(NumericPointwiseOperators.max(im1) - min1); + float min2 = NumericPointwiseOperators.min(im2); + float scaleFactor2 = (histSize - 1)/(NumericPointwiseOperators.max(im2) - min2); + + // calculate joint histogram + for (int i = 0; i < imWidth; i++) { + + for (int j = 0; j < imHeight; j++) { + + int value_ref = (int) (scaleFactor1*(im1.getAtIndex(i, j) - min1)); + int value_mov = (int) (scaleFactor2*(im2.getAtIndex(i, j) - min2)); + + // jH(k,l) counts how often the intensity pair k in the reference image, and l in the moving image occurs (at the corresponding location) + // use the correct indices and set the value for jH properly + //TODO + } + } + + // divide by the number of elements in order to get probabilities + //TODO + + return jH; + } + + /** + * Method to calculate a histogram from a joint histogram + * @param jH the joint histogram + * @return a SimpleVector corresponding to the marginal histogram + */ + // TODO-s here! + public SimpleVector getHistogramFromJointHistogram(SimpleMatrix jH, boolean sumRows){ + + int numCols = jH.getCols(); + int numRows = jH.getRows(); + + int histSize = sumRows ? numCols : numRows; + + SimpleVector hist = new SimpleVector(histSize); + + for(int i = 0; i < numRows; i++) { + + for(int j = 0; j < numCols; j++) { + + if (sumRows) { + + // sum up over the rows + //TODO + } + else { + + // sum up over the columns + //TODO + } + } + } + + return hist; + } + + /** + * + * end of the exercise + */ + + public ExerciseMI() { + + // Load images + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + + String filename1 = imageDataLoc + "T1.png"; + String filename2 = imageDataLoc + "Proton.png"; + + image1 = readImageFile(filename1); + image2 = readImageFile(filename2); + + // set the origin of the image in its center (the default origin of an image is in its top left corner) + centerOrigin(image1,new double[]{1.0,1.0}); + centerOrigin(image2,new double[]{1.0,1.0}); + + // blurred images for the registration to avoid local minima during optimization + reference_img = blurImage(image1,blurSigma); + moving_img = blurImage(image2,blurSigma); + + movingStack = new Grid3D(reference_img.getWidth(), reference_img.getHeight(), iterLimit, false); + } + + public Grid2D readImageFile(String filename){ + return ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); + } + + void centerOrigin(Grid2D image, double[] spacing) { + + image.setOrigin(-(image.getWidth()-1)/2 , -(image.getHeight()-1)/2); + image.setSpacing(spacing); + } + + Grid2D blurImage(Grid2D image, double sigma) { + + Grid2D image_blurred = new Grid2D(image); + IJ.run(ImageUtil.wrapGrid(image_blurred,""),"Gaussian Blur...", "sigma="+sigma); + + return image_blurred; + } + + public double[] performRegistration(CostFunction cF){ + + FunctionOptimizer fo = new FunctionOptimizer(); + + fo.setDimension(3); + fo.setNdigit(6); + fo.setItnlim(iterLimit); + fo.setMsg(16); + fo.setInitialX(new double[]{0,0,0}); + fo.setMaxima(new double[]{50,50,50}); + fo.setMinima(new double[]{-50,-50,-50}); + fo.setOptimizationMode(OptimizationMode.Function); + + // enable reading intermediate results + ArrayList visFcts = new ArrayList(); + visFcts.add(cF); + fo.setCallbackFunctions(visFcts); + + double result[] = fo.optimizeFunction(cF); + + movingStack = cropMovingStack(movingStack, cF.getIteration()); + + return result; + } + + private Grid3D cropMovingStack(Grid3D stack, int bound) { + + int[] mStackSize = stack.getSize(); + + ArrayList buffer = stack.getBuffer(); + + Grid3D newStack = new Grid3D(mStackSize[0], mStackSize[1], bound); + for (int i=0; i costPerIter_map; + + Grid2D tmp_img; + int iteration = 0; + int subIteration = 0; // only for debugging + + @Override + public void setNumberOfProcessingBlocks(int number) { + } + + @Override + public int getNumberOfProcessingBlocks() { + return 1; + } + + @Override + public double evaluate(double[] x, int block) { + + // define rotation and translation + double phi = x[0]*Math.PI/180; + SimpleVector t = new SimpleVector(x[1], x[2]); + RigidTransform rigid = new RigidTransform(phi, t); + + // perform rotation/translation + tmp_img = new Grid2D(moving_img); + tmp_img.applyTransform(rigid); + + // calculate the cost function + double cost = calculateMutualInformation(reference_img, tmp_img); + + subIteration++; + return cost; + } + + @Override + public void optimizerCallbackFunction(int currIterationNumber, double[] x, + double currFctVal, double[] gradientAtX) { + + iteration = currIterationNumber; + + // fill movingStack with current transformed image + movingStack.setSubGrid(currIterationNumber, tmp_img); + + // generate cost plot data + if (costPerIter_map == null) + costPerIter_map = new TreeMap(); + + costPerIter_map.put(currIterationNumber, currFctVal); + } + + public Grid1D getCostPerIterGrid() { + + if (costPerIter_map == null) + return new Grid1D(1); + + Grid1D costPerIter_grid = new Grid1D(costPerIter_map.size()); + + Iterator> it = costPerIter_map.entrySet().iterator(); + while (it.hasNext()) { + + Entry e = it.next(); + costPerIter_grid.setAtIndex(e.getKey(), e.getValue().floatValue()); + } + + return costPerIter_grid; + } + + public int getIteration() { + return iteration; + } + } + + public class RotationMatrix2D extends SimpleMatrix { + + private static final long serialVersionUID = 6708400687838433556L; + + public RotationMatrix2D() { + this(0.0); + } + + public RotationMatrix2D(double phi) { + + super(2,2); + + this.setElementValue(0, 0, Math.cos(phi)); + this.setElementValue(0, 1, - Math.sin(phi)); + this.setElementValue(1, 0, Math.sin(phi)); + this.setElementValue(1, 1, Math.cos(phi)); + } + } + + public class RigidTransform extends AffineTransform { + + private static final long serialVersionUID = 3469069367283792094L; + + public RigidTransform(double rotationAngle, SimpleVector translationVector) { + super(new RotationMatrix2D(rotationAngle), translationVector); + } + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public double get_blurSigma() { + return blurSigma; + } + // + public int get_iterLimit() { + return iterLimit; + } + // + public Grid2D get_image1() { + return image1; + } + // + public Grid2D get_image2() { + return image2; + } + // + public Grid2D get_reference_img() { + return reference_img; + } + // + public Grid2D get_moving_img() { + return moving_img; + } + // + public Grid3D get_movingStack() { + return movingStack; + } +} \ No newline at end of file diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExercisePB.java b/src/edu/stanford/rsl/tutorial/mipda/ExercisePB.java new file mode 100644 index 00000000..ad76fa56 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExercisePB.java @@ -0,0 +1,487 @@ +package edu.stanford.rsl.tutorial.mipda; + +import java.util.ArrayList; + +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.conrad.data.numeric.InterpolationOperators; +import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; +import edu.stanford.rsl.conrad.geometry.shapes.simple.Box; +import edu.stanford.rsl.conrad.geometry.shapes.simple.PointND; +import edu.stanford.rsl.conrad.geometry.shapes.simple.StraightLine; +import edu.stanford.rsl.conrad.geometry.transforms.Transform; +import edu.stanford.rsl.conrad.geometry.transforms.Translation; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import edu.stanford.rsl.tutorial.phantoms.SheppLogan; +import ij.ImageJ; + +/** + * Ramp Filter for Parallel Beam Reconstruction + * Programming exercise for module "Parallel Beam Reconstruction" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Bastian Bier, Ashwini Jadhav, Anna Gebhard, Mena Abdelmalek + * + */ + +public class ExercisePB { + + public enum RampFilterType {NONE, RAMLAK, SHEPPLOGAN}; + boolean filterShownOnce = false; + + RampFilterType filterType = RampFilterType.NONE; //TODO: Select one of the following values: NONE, RAMLAK, SHEPPLOGAN + // (make this choice when you have finished the rest of the exercise) + + SheppLogan sheppLoganPhantom; + + // parameters for the phantom + final int phantomSize = 256; // size of the phantom + final float world_unit = 1.0f; // pixel dimension + + // parameters for projection + final int projectionNumber = 180;// number of projection images + final double angularRange = Math.PI; // projection image range + final double angularStepLength = angularRange / projectionNumber; // angle in between adjacent projections + // + final float detectorSize = 250; // detector size in pixel + final float detectorSpacing = 1.0f; // size of a detector Element [mm] + final float detectorLength = detectorSize*detectorSpacing; + + // parameters for filtering + int paddedSize = 0; // determined later by phantomSize and Grid1DComplex + + // parameters for reconstruction + final float recoPixelDim = world_unit; // isotropic image pixels + final int recoSize = (int) ((int) phantomSize*world_unit/recoPixelDim); + + // outputs + Grid2D sinogram; + Grid1DComplex ramp; + Grid2D filteredSinogram; + Grid2D recoImage; + + + public static void main(String[] args) { + + //ImageJ ij = new ImageJ(); // optional TODO: uncomment if you want to analyze the output images + + ExercisePB parallel = new ExercisePB(); // Step 1: creates the Shepp-Logan phantom as example image + parallel.get_sheppLoganPhantom().show("Original Image (Shepp-Logan)"); + + // Step 2: acquire forward projection images with a parallel projector + parallel.sinogram = parallel.parallelProjection(); + parallel.get_sinogram().show("The Sinogram"); + + // Step 3: ramp filtering + parallel.filteredSinogram = parallel.filterSinogram(); // All TODO-s here! + parallel.get_filteredSinogram().show("Filtered Sinogram"); + + // Step 4: reconstruct the object by backprojection of the information in the filtered sinogram + parallel.recoImage = parallel.backprojectFilteredSinogram(); + parallel.get_recoImage().show("Reconstructed Image"); + } + + /** + * Filtering the sinogram with a high pass filter + * + * The ramp filters are defined in the spatial domain but + * they are applied in the frequency domain. + * Remember: a convolution in the spatial domain corresponds to a + * multiplication in the frequency domain. + * + * Both, the sinogram and the ramp filter are transformed into + * the frequency domain and multiplied there. + * + * @param first input is a line of the sinogram + * + */ + public Grid1D rampFiltering(Grid1D projection, RampFilterType filterType){ + + // initialize the ramp filter + ramp = new Grid1DComplex(projection.getSize()[0]); + + // define the filter in the spatial domain on the fully padded size + paddedSize = ramp.getSize()[0]; + + double deltaS = projection.getSpacing()[0]; + + // RampFilterType: RAMLAK, SHEPPLOGAN, NONE + switch (filterType) { + + case RAMLAK: // RamLak filter + ramlak(ramp, paddedSize, deltaS); // TODO: go to this method and implement it + break; + + case SHEPPLOGAN: // SheppLogan filter + sheppLogan(ramp, paddedSize, deltaS); // TODO: go to this method and implement it + break; + + default: // NONE (or other) + return projection; + } + + if (!filterShownOnce) + ramp.show("Ramp Filter in Spatial Domain (rearranged for FFT-Shift)"); + + // // TODO: Transform the filter into frequency domain (look for an appropriate method of Grid1DComplex) + + if(!filterShownOnce) { + + ramp.show("Ramp Filter in Frequency Domain"); + filterShownOnce = true; + } + + Grid1DComplex projectionF = null;// TODO: Transform the input sinogram signal ... + // // ... into the frequency domain (hint: similar to the last TODO) + + if (projectionF != null) { + for(int p = 0; p < projectionF.getSize()[0]; p++){ + // // TODO: Multiply the transformed sinogram with the ramp filter (complex multiplication) + } + } + + // // TODO: transform back to get the filtered sinogram (i.e. invert the Fourier transform) + + // crop the image to its initial size + Grid1D grid = new Grid1D(projection); + if (projectionF != null) { + grid = projectionF.getRealSubGrid(0, projection.getSize()[0]); + } + + NumericPointwiseOperators.multiplyBy(grid, (float) deltaS); // integral discretization (filtering) + + return grid; + } + + // TODO: implement the ram-lak filter in the spatial domain + public void ramlak(Grid1DComplex filterGrid, int paddedSize, double deltaS) { + + final float constantFactor = -1.f / ((float) ( Math.PI * Math.PI * deltaS * deltaS)); + + // // TODO: set correct value in filterGrid for zero frequency + + for (int i = 1; i < paddedSize/2; ++i) { // the "positive wing" of the filter + + if (false) {// TODO: condition -> only odd indices are nonzero + // // TODO: use setAtIndex and the constant "constantFactor" + } + } + + // remark: the sorting of frequencies in the Fourier domain is 0,...,k,-k,...,1 + // to implement the "negative wing", a periodicity assumption is used and the values are added accordingly + for (int i = paddedSize / 2; i < paddedSize; ++i) { + + final int tmp = paddedSize - i; // now we go back from N/2 to 1 + if (false) {// TODO: condition -> only odd indices are nonzero + // // TODO: use setAtIndex and the constant "constantFactor" + } + } + } + + // TODO: implement the Shepp-Logan filter in the spatial domain + public void sheppLogan(Grid1DComplex filterGrid, int paddedSize, double deltaS) { + + final float constantFactor = - 2.f / ((float) ( Math.PI * Math.PI * deltaS * deltaS)); + + // // TODO: set correct value in filterGrid for zero frequency + + for (int i = 1; i < paddedSize/2; ++i){ // the "positive wing" of the filter + // // TODO: use setAtIndex and the constant "constantFactor" + } + + // remark: the sorting of frequencies in the Fourier domain is 0,...,k,-k,...,1 + // to implement the "negative wing", a periodicity assumption is used and the values are added accordingly + for (int i = paddedSize / 2; i < paddedSize; ++i) { + + final float tmp = paddedSize - i; // now we go back from N/2 to 1 + // // TODO: use setAtIndex and the constant "constantFactor" + } + } + + /** + * The following content is read-only. + * Both methods + * 1. projectRayDriven(...) + * 2. and backprojectPixelDriven(...) + * are used above to + * 1. generate the sinogram + * 2. and to reconstruct the phantom from it. + * You can try to understand it if you want to see how these operators can be implemented. + * Otherwise the exercise is finished here. + */ + + /** + * Forward projection of the phantom onto the detector + * Rule of thumb: Always sample in the domain where you expect the output! + * Thus, we sample at the detector pixel positions and sum up the informations along one ray + * + * @param grid the image + * @param maxTheta the angular range in radians + * @param deltaTheta the angular step size in radians + * @param maxS the detector size in [mm] + * @param deltaS the detector element size in [mm] + */ + public Grid2D projectRayDriven(Grid2D grid, double maxTheta, double deltaTheta, double maxS, double deltaS) { + + final double samplingRate = 3.d; // # of samples per pixel + + // prepare output (sinogram) + int maxThetaIndex = (int) (maxTheta / deltaTheta); + int maxSIndex = (int) (maxS / deltaS); + + Grid2D sino = new Grid2D(new float[maxThetaIndex*maxSIndex], maxSIndex, maxThetaIndex); // rows: angle, columns: detector pixel + sino.setSpacing(deltaS, deltaTheta); + + // set up image bounding box in WC + Translation trans = new Translation( + -(grid.getSize()[0] * grid.getSpacing()[0])/2, -(grid.getSize()[1] * grid.getSpacing()[1])/2, -1); + Transform inverse = trans.inverse(); + // + Box b = new Box((grid.getSize()[0] * grid.getSpacing()[0]), (grid.getSize()[1] * grid.getSpacing()[1]), 2); + b.applyTransform(trans); // box now centered at world origin + + for(int e=0; e points = b.intersect(line); + + // only if we have intersections + if (points.size() != 2){ + if(points.size() == 0) { + line.getDirection().multiplyBy(-1.d); + points = b.intersect(line); + } + if(points.size() == 0) // ray passes the image area + continue; + } + + PointND start = points.get(0); // [mm] + PointND end = points.get(1); // [mm] + + // get the normalized increment + SimpleVector increment = new SimpleVector(end.getAbstractVector()); + increment.subtract(start.getAbstractVector()); + double distance = increment.normL2(); + increment.divideBy(distance * samplingRate); + + double sum = .0; + + start = inverse.transform(start); // use image coordinates + + // compute the integral along the line. + for (double t = 0.0; t < distance * samplingRate; ++t) { + + PointND current = new PointND(start); + current.getAbstractVector().add(increment.multipliedBy(t)); + + double x = current.get(0) / grid.getSpacing()[0]; + double y = current.get(1) / grid.getSpacing()[1]; + + if (grid.getSize()[0] <= x + 1 // ray outside right of image + || grid.getSize()[1] <= y + 1 // ray above image + || x < 0 // ray outside left of image + || y < 0) // ray below image + continue; + + sum += InterpolationOperators.interpolateLinear(grid, x, y); + } + + // normalize by the number of interpolation points + sum /= samplingRate; + // write integral value into the sinogram. + sino.setAtIndex(i, e, (float)sum); + } + } + + return sino; + } + + /** + * Backprojection of the projections/sinogram + * The projections are created pixel driven. + * -> Rule of thumb: Always sample in the domain where you expect the output! + * Here, we want to reconstruct the volume, thus we sample in the reconstructed grid! + * + * @param sino: the sinogram + * @param imageSizeX: x-dimension of reconstructed image + * @param imageSizeY: y-dimension of reconstructed image + * @param pxSzXMM: granularity in x + * @param pxSzYMM: granularity in y + */ + + public Grid2D backprojectPixelDriven(Grid2D sino, int imageSizeX, int imageSizeY, float pxSzXMM, float pxSzYMM) { + + int maxSIndex = sino.getSize()[0]; + double deltaS = sino.getSpacing()[0]; + int maxThetaIndex = sino.getSize()[1]; + double deltaTheta = sino.getSpacing()[1]; + + double maxS = maxSIndex * deltaS; // detector length + + Grid2D grid = new Grid2D(imageSizeX, imageSizeY); + grid.setSpacing(pxSzXMM, pxSzYMM); + grid.setOrigin(-(grid.getSize()[0]*grid.getSpacing()[0])/2, -(grid.getSize()[1]*grid.getSpacing()[1])/2); + + // loop over the projection angles + for (int i = 0; i < maxThetaIndex; i++) { + + // compute actual value for theta + double theta = deltaTheta * i; + + // get detector direction vector + double cosTheta = Math.cos(theta); + double sinTheta = Math.sin(theta); + SimpleVector dirDetector = new SimpleVector(cosTheta*grid.getSpacing()[0],sinTheta*grid.getSpacing()[1]); + dirDetector.normalizeL2(); + + // loops over the image grid + for (int x = 0; x < grid.getSize()[0]; x++) { + for (int y = 0; y < grid.getSize()[1]; y++) { + + // compute world coordinate of current pixel + double[] w = grid.indexToPhysical(x, y); + + // wrap into vector + SimpleVector pixel = new SimpleVector(w[0], w[1]); + + // project pixel onto detector + double s = SimpleOperators.multiplyInnerProd(pixel, dirDetector); + // compute detector element index from world coordinates + s += maxS/2; // [mm] + s /= deltaS; // [GU] + + // get detector grid + Grid1D subgrid = sino.getSubGrid(i); + + // check detector bounds, continue if out of array + if (subgrid.getSize()[0] <= s + 1 || s < 0) + continue; + + // get interpolated value and add value to sinogram + float val = InterpolationOperators.interpolateLinear(subgrid, s); + grid.addAtIndex(x, y, val); + } + } + } + + // apply correct scaling + NumericPointwiseOperators.multiplyBy(grid, (float) deltaTheta); // integral discretization (theta) + + return grid; + } + + /** + * + * end of the exercise + */ + + public ExercisePB() { + + sheppLoganPhantom = new SheppLogan(phantomSize); + sheppLoganPhantom.setSpacing(world_unit,world_unit); + } + + public Grid2D parallelProjection(){ + + Grid2D grid = projectRayDriven(sheppLoganPhantom, angularRange, angularStepLength, detectorLength, detectorSpacing); + return grid; + } + + public Grid2D filterSinogram(){ + + Grid2D grid = new Grid2D(sinogram); + + for (int theta = 0; theta < sinogram.getSize()[1]; ++theta) { + + // Filter each line of the sinogram independently + Grid1D tmp = rampFiltering(sinogram.getSubGrid(theta), filterType); + + for(int i = 0; i < tmp.getSize()[0]; i++) { + grid.putPixelValue(i, theta, tmp.getAtIndex(i)); + } + } + + return grid; + } + + public Grid2D backprojectFilteredSinogram(){ + + Grid2D reco = backprojectPixelDriven(filteredSinogram, recoSize, recoSize, recoPixelDim, recoPixelDim); + + return reco; + + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public RampFilterType get_RampFilterType() { + return filterType; + } + public SheppLogan get_sheppLoganPhantom() { + return sheppLoganPhantom; + } + public int get_phantomSize() { + return phantomSize; + } + public float get_world_unit() { + return world_unit; + } + public int get_projectionNumber() { + return projectionNumber; + } + public double get_angularRange() { + return angularRange; + } + public double get_angularStepLength() { + return angularStepLength; + } + public float get_detectorSize() { + return detectorSize; + } + public float get_detectorSpacing() { + return detectorSpacing; + } + public float get_detectorLength() { + return detectorLength; + } + public int get_paddedSize() { + return paddedSize; + } + public float get_recoPixelDim() { + return recoPixelDim; + } + public int get_recoSize() { + return recoSize; + } + public Grid2D get_sinogram() { + return sinogram; + } + public Grid1DComplex get_ramp() { + return ramp; + } + public Grid2D get_filteredSinogram() { + return filteredSinogram; + } + public Grid2D get_recoImage() { + return recoImage; + } +} + diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseREG.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseREG.java new file mode 100644 index 00000000..dd45d633 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseREG.java @@ -0,0 +1,242 @@ +package edu.stanford.rsl.tutorial.mipda; + +import java.awt.Color; + +import edu.stanford.rsl.conrad.numerics.SimpleMatrix; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import ij.gui.Plot; + + +/** + * Registration of Point Correspondences + * Programming exercise for chapter "Rigid Registration" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Bastian Bier, Frank Schebesch, Anna Gebhard, Mena Abdelmalek, Ashwini Jadhav + * + */ + +public class ExerciseREG { + + SimpleMatrix p_k; + SimpleMatrix q_k; + Plot plot; + + public static void main(String[] args){ + + ExerciseREG exREGobj = new ExerciseREG(); + + // calculate registration parameters: rotation angle and translation vector + RigidParameters parameter = exREGobj.registerPoints(exREGobj.p_k, exREGobj.q_k); //TODO: Go to this method and implement the missing steps! + + // transform points + SimpleMatrix transformedPoints = exREGobj.applyTransform(exREGobj.q_k, parameter); //TODO: Go to this method and implement the missing steps! + + // add transformed point cloud into the plot + exREGobj.plot.setColor(Color.red); + exREGobj.plot.addPoints(transformedPoints.getCol(0).copyAsDoubleArray(), transformedPoints.getCol(1).copyAsDoubleArray(), Plot.CROSS); + exREGobj.plot.show(); + } + + /** + * Registration method + * Point correspondences are used to calculate the optimal rotational and translational + * transformations from p to q + * @param p reference point cloud + * @param q point cloud to be registered + * @return parameter object, contains rotation angle and translation vector (phi, [t1, t2]) + */ + public RigidParameters registerPoints(SimpleMatrix p, SimpleMatrix q){ + + int numPoints = p.getRows(); + + // build up measurement matrix m for the complex numbers problem + SimpleMatrix m = buildMeasurementMatrix(q);// TODO-s here! + + // build up right hand side vector b for the complex numbers problem + SimpleVector b = buildRHSVector(p);// TODO-s here! + + // just some check (ignore) + assert(m.getRows() == 2*numPoints); + + // solve the problem + SimpleMatrix m_inv = m.inverse(InversionType.INVERT_SVD);// calculate pseudo inverse of m + SimpleVector x = SimpleOperators.multiply(m_inv, b);// solve for parameters + + double r1 = x.getElement(0); + double r2 = x.getElement(1); + double t1 = x.getElement(2); + double t2 = x.getElement(3); + + //calculate the absolute value of r = r1+i*r2 + double abs_r = 0.d;//TODO + //TODO: normalize r + //TODO: normalize r + + //calculate the angle phi + double phi = 0.d;//TODO + + // return both translation and rotation in a RigidParameters object + return new RigidParameters(phi,new SimpleVector(t1,t2)); + } + + public SimpleMatrix buildMeasurementMatrix(SimpleMatrix q) { + + int numPoints = q.getRows(); + + SimpleMatrix m = new SimpleMatrix(2*numPoints, 4); + + // build up measurement matrix m for the complex numbers problem + for(int i = 0; i < numPoints; i++) { + + //real part of the problem (even rows) + //TODO + //TODO + //TODO + //TODO + //imaginary part of the problem (odd rows) + //TODO + //TODO + //TODO + //TODO + } + + return m; + } + + public SimpleVector buildRHSVector(SimpleMatrix p) { + + int numPoints = p.getRows(); + + SimpleVector b = new SimpleVector(2*numPoints);// right hand side + + for(int i = 0; i < numPoints; i++) { + + //TODO: real part of the problem + //TODO: imaginary part of the problem + } + + return b; + } + + /** + * Transformation method + * Computes the registered points by transforming the original points "points" using the parameters from "parameter" + * @param points point cloud to be registered + * @param parameter rigid transformation parameters + * @return transformedPoints registered point cloud + */ + public SimpleMatrix applyTransform(SimpleMatrix points, RigidParameters parameter){ + + SimpleMatrix r = new RotationMatrix2D(parameter.getAngle()); + SimpleMatrix transformedPoints = new SimpleMatrix(points.getRows(), points.getCols()); + + // transform points (rotate and translate) + for(int i = 0; i < transformedPoints.getRows(); i++){ + + //TODO: rotate (you can use SimpleOperators) + //TODO: translate (you can use SimpleOperators) + } + + return transformedPoints; + } + + /** + * + * end of the exercise + */ + + public ExerciseREG(){ + + // Define Point Cloud + p_k = new SimpleMatrix(4,2); + q_k = new SimpleMatrix(4,2); + + p_k.setRowValue(0, new SimpleVector(1,2)); + p_k.setRowValue(1, new SimpleVector(3,6)); + p_k.setRowValue(2, new SimpleVector(4.5,5.5)); + p_k.setRowValue(3, new SimpleVector(3,3.5)); + + q_k.setRowValue(0, new SimpleVector(-0.7, 2.1)); + q_k.setRowValue(1, new SimpleVector(-2.1, 6.4)); + q_k.setRowValue(2, new SimpleVector(-1.4, 7.1)); + q_k.setRowValue(3, new SimpleVector(-0.7, 4.9)); + + // Plot Point Cloud + plot = new Plot("Regression Line", "X", "Y", Plot.DEFAULT_FLAGS); + plot.setLimits(-5, 5, 0, 10); + plot.setSize(512, 512); + plot.setScale(2.0f); + plot.setLineWidth(1.5f); + plot.addLegend("cloud {p_k}\ncloud {q_k}\nregistered {q_k}\nd"); + plot.addPoints(p_k.getCol(0).copyAsDoubleArray(), p_k.getCol(1).copyAsDoubleArray(), Plot.BOX); + plot.setColor(Color.blue); + plot.addPoints(q_k.getCol(0).copyAsDoubleArray(), q_k.getCol(1).copyAsDoubleArray(), Plot.CIRCLE); + + } + + public class RotationMatrix2D extends SimpleMatrix { + + private static final long serialVersionUID = -6939956675976343158L; + + public RotationMatrix2D() { + this(0.0); + } + + public RotationMatrix2D(double phi) { + + super(2,2); + + this.setElementValue(0, 0, Math.cos(phi)); + this.setElementValue(0, 1, - Math.sin(phi)); + this.setElementValue(1, 0, Math.sin(phi)); + this.setElementValue(1, 1, Math.cos(phi)); + } + } + + public class RigidParameters { + + private double phi; + private SimpleVector translation; + + public RigidParameters(){ + + phi = 0.0d; + translation = new SimpleVector(2); + } + public RigidParameters(final double phi, final SimpleVector translation){ + + this.phi = phi; + this.translation = translation; + } + // + public void setAngle(final double phi) { + this.phi = phi; + } + public void setTranslation(final SimpleVector translation) { + this.translation = translation; + } + // + public double getAngle() { + return phi; + } + public SimpleVector getTranslation() { + return translation; + } + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public SimpleMatrix get_p_k() { + return p_k; + } + // + public SimpleMatrix get_q_k() { + return p_k; + } + // + public Plot get_plot() { + return plot; + } +} diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseSVD.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseSVD.java new file mode 100644 index 00000000..b8b83e17 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseSVD.java @@ -0,0 +1,504 @@ +package edu.stanford.rsl.tutorial.mipda; + +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.data.numeric.Grid3D; +import edu.stanford.rsl.conrad.data.numeric.NumericGridOperator; +import edu.stanford.rsl.conrad.fitting.LinearFunction; +import edu.stanford.rsl.conrad.numerics.DecompositionSVD; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix.InversionType; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix.MatrixNormType; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import edu.stanford.rsl.conrad.utils.ImageUtil; +import edu.stanford.rsl.conrad.utils.VisualizationUtil; +import ij.IJ; +import ij.ImageJ; + +/** + * Singular Value Decomposition + * Programming exercise for module "Mathematical Tools" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Marco Boegel, Ashwini Jadhav, Mena Abdelmalek, Anna Gebhard + * + */ + +public class ExerciseSVD { + + float error[]; + + public static void main(String[] args) { + + ImageJ ij = new ImageJ(); + + //Create matrix A + // 11 10 14 + // 12 11 -13 + // 14 13 -66 + SimpleMatrix A = new SimpleMatrix(3,3); + A.setRowValue(0, new SimpleVector(11, 10, 14)); + A.setRowValue(1, new SimpleVector(12, 11, -13)); + A.setRowValue(2, new SimpleVector(14, 13, -66)); + + // Data for optimization problem 2 from lecture slides + SimpleMatrix vectors = new SimpleMatrix(2,4); + vectors.setColValue(0, new SimpleVector(1.f, 1.f)); + vectors.setColValue(1, new SimpleVector(-1.f, 2.f)); + vectors.setColValue(2, new SimpleVector(1.f, -3.f)); + vectors.setColValue(3, new SimpleVector(-1.f, -4.f)); + + //Load the head angiography image image + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + String filename = imageDataLoc + "mr_head_angio.jpg"; + Grid2D image = ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); + image.show("mr_head_angio:original"); + int rank = 150; + + //Data for problem 4 + double[] xCoords = new double[]{-3.f, -2.f, -1.f, 0.f, 1.5f, 2.f, 3.1f, 5.9f, 7.3f}; + double[] yCoords = new double[]{7.f, 8.f, 9.f, 3.3f, 2.f, -3.f, 4.f, -0.1, -0.5}; + + ExerciseSVD exsvd = new ExerciseSVD(); + + //create the svd of A + DecompositionSVD svd = exsvd.createSVD(A); //There are TODO-s here! + + //compute the pseudo inverse of A + exsvd.pseudoInverse(A); + + //Show that a change in a vector b of only 0.1% can lead to 240% change + //in the result of A*x=b + SimpleVector b = new SimpleVector(1.001, 0.999, 1.001); + //diff implies a 0.1% change + SimpleVector diff = new SimpleVector(-0.001, 0.001, -0.001); + //show the difference between the results + exsvd.showDifference(A, b, diff); //There are TODO-s here! + + //compute the condition number of A + exsvd.conditionNumber(A); + + //calculate the minimal dimension of S of the svd of A + int minDim = 0; + if (svd != null) { + minDim = exsvd.minDim(svd.getS()); //There are TODO-s here! + } + + //create a low rank matrix + exsvd.lowRankMatrix(A, minDim); //There are TODO-s here! + + //compute matrix closest to A according to the Frobenius norm + //with rank deficiency one and two identical non-zero singular values + int svdNewRank = 0; + if (svd != null) { + svdNewRank = exsvd.newRank(svd.rank(), 1); //There are TODO-s here! + } + exsvd.optimizationProblem1(A, svdNewRank); //There are TODO-s here! + + //estimate the matrix A that minimizes sum_{i=1}^{4} b'_{i} A b_{i} + //under the constraint ||A||_{F} = 1 + exsvd.optimizationProblem2(vectors); + + //compute the image matrix of rank 1 that best approximates an image + exsvd.optimizationProblem3(image, rank); //There are TODO-s here! + + //compute the regression line through the given 2-D points + exsvd.optimizationProblem4(xCoords, yCoords); //There are TODO-s here! + + } + + public DecompositionSVD createSVD(SimpleMatrix A) { + System.out.println("A = " + A.toString()); + + //Compute the SVD of A + DecompositionSVD svd = null; //TODO + + //Check output: re-compute A = U * S * V^T + if (svd != null) { + SimpleMatrix temp = SimpleOperators.multiplyMatrixProd(svd.getU(), svd.getS()); + SimpleMatrix A2 = SimpleOperators.multiplyMatrixProd(temp, svd.getV().transposed()); + System.out.println("U * S * V^T: " + A2.toString()); + } + + return svd; + } + + public int minDim(SimpleMatrix A) { + return 0; //TODO + } + + public SimpleMatrix lowRankMatrix(SimpleMatrix A, int minDim){ + + //introduce a rank deficiency + DecompositionSVD svd = new DecompositionSVD(A); + double eps = 10e-3; + SimpleMatrix Slowrank = new SimpleMatrix(svd.getS()); + + for(int i = 0; i < minDim; i++) + { + double val = svd.getS().getElement(i, i); + if(val> eps) + { + Slowrank.setElementValue(i, i, val); + } + } + + SimpleMatrix templowrank = SimpleOperators.multiplyMatrixProd(svd.getU(), Slowrank); + SimpleMatrix Alowrank = null; //TODO + if (Alowrank != null) { + System.out.println("A rank deficient = " + Alowrank.toString()); + } + return Alowrank; + + } + + public SimpleVector showDifference(SimpleMatrix A, SimpleVector b, SimpleVector diff) { + + SimpleMatrix Ainv = pseudoInverse(A); + SimpleVector br = new SimpleVector(b); + br.add(diff); + //solve for x and xr + SimpleVector x = SimpleOperators.multiply(Ainv, b); + System.out.println(x.toString()); + + SimpleVector xr = SimpleOperators.multiply(Ainv, br); + System.out.println(xr.toString()); + + //We want only the difference caused by the change + SimpleVector xn = vecDiff(xr, x); //There is a TODO here! + + //Alternatively: + //SimpleVector bn = SimpleOperators.subtract(br, b); + //SimpleVector xn = SimpleOperators.multiply(Ainv, bn); + + // compute and show percentual change + SimpleVector xPercentage = null; //TODO + //TODO + return xPercentage; + + } + + public SimpleVector vecDiff(SimpleVector x1, SimpleVector x2) { + return null; //TODO + } + + public int newRank(int oldRank, int rankDeficiency) { + return 0; //TODO + } + + public SimpleMatrix optimizationProblem1(SimpleMatrix A, int svdNewRank) + { + System.out.println("Optimization Problem 1"); + //%%%%%%%%%%%%%%%%%%%%%%%%%%%% + // Optimization problem I + //%%%%%%%%%%%%%%%%%%%%%%%%%%% + + // Let us consider the following problem that appears in many image + // processing and computer vision problems: + // We computed a matrix A out of sensor data like an image. By theory + // the matrix A must have the singular values s1, s2, . . . , sp, where + // p = min(m, n). Of course, in practice A does not fulfill this constraint. + // Problem: What is the matrix A0 that is closest to A (according to the + // Frobenius norm) and has the required singular values? + + + // Let us assume that by theoretical arguments the matrix A is required + // to have a rank deficiency of one, and the two non-zero singular values + // are identical. The matrix A0 that is closest to A according to the + // Frobenius norm and fulfills the above requirements is: + + // Build mean of first two singular values to make them identical + // Set third singular value to zero + + DecompositionSVD svd = new DecompositionSVD(A); + + //Compute mean of remaining singular values + double mean = 0; + + for(int i = 0; i < svdNewRank; i++) + { + mean += svd.getS().getElement(i, i); + } + + mean /= svdNewRank; + + //Create new Singular matrix + SimpleMatrix Slowrank = new SimpleMatrix(svd.getS().getRows(), svd.getS().getCols()); + Slowrank.zeros(); + //Fill in remaining singular values with the mean. + for(int i = 0; i < svdNewRank; i++) + { + Slowrank.setElementValue(i, i, mean); + } + + //compute A0 + SimpleMatrix tempA0 = null; //TODO + SimpleMatrix A0 = null; //TODO + + if (A0 != null) { + System.out.println("A0 = " + A0.toString()); + + double normA = A.norm(MatrixNormType.MAT_NORM_FROBENIUS); + double normA0 = A0.norm(MatrixNormType.MAT_NORM_FROBENIUS); + + System.out.println("||A||_F = " + normA); + System.out.println("||A0||_F = " + normA0); + } + + return A0; + + } + + public void optimizationProblem3(Grid2D image, int rank) + { + System.out.println("Optimization Problem 3"); + //% + //%%%%%%%%%%%%%%%%%%%%%%%%%%% + // Optimization problem III + //%%%%%%%%%%%%%%%%%%%%%%%%%%% + + // The SVD can be used to compute the image matrix of rank 1 that best + // approximates an image. Figure 1 shows an example of an image I + // and its rank-1-approximation I0 = u_{1} * s_{1} * v'_{1}. + + //In order to apply the svd, we first need to transfer our Grid2D image to a matrix + SimpleMatrix I = new SimpleMatrix(image.getHeight(), image.getWidth()); + + for(int i = 0; i < image.getHeight(); i++) + { + for(int j = 0; j < image.getWidth(); j++) + { + double val = image.getAtIndex(j, i); + I.setElementValue(i, j, val); + } + } + + DecompositionSVD svd = new DecompositionSVD(I); + + Grid3D imageRanks = new Grid3D(image.getWidth(), image.getHeight(), rank); + + // track error + error = new float[rank]; + + //Create Rank k approximations + for(int k = 0; k < rank; k++) + { + SimpleVector us = svd.getU().getCol(k).multipliedBy(svd.getSingularValues()[k]); + SimpleMatrix Iapprox = SimpleOperators.multiplyOuterProd(us, svd.getV().getCol(k)); + + + //Transfer back to grid + Grid2D imageRank = new Grid2D(image.getWidth(),image.getHeight()); + for(int i = 0; i < image.getHeight(); i++) + { + for(int j = 0; j < image.getWidth(); j++) + { + if(k == 0) + { + imageRank.setAtIndex(j, i, (float) Iapprox.getElement(i, j)); + + } + else + { + imageRank.setAtIndex(j, i, imageRanks.getAtIndex(j, i, k-1) + (float) Iapprox.getElement(i, j)); + } + + + } + } + imageRanks.setSubGrid(k, imageRank); + + // evaluate RMSE for the k-th approximation + error[k] = calculateRMSE(image,imageRank); //There is a TODO here! + } + + imageRanks.show("Stack of approximations"); + VisualizationUtil.createPlot("Rank-Error-Plot", error).show(); + + //Direct estimation of rank K (using rank 1-matrices, cf. lecture) [alternative solution] + SimpleMatrix usK = SimpleOperators.multiplyMatrixProd(svd.getU().getSubMatrix(0, 0, svd.getU().getRows(), rank), svd.getS().getSubMatrix(0, 0, rank, rank)); + SimpleMatrix IapproxK = SimpleOperators.multiplyMatrixProd(usK, svd.getV().getSubMatrix(0, 0, svd.getV().getRows(), rank).transposed()); + + //Transfer back to grid + Grid2D imageRankK = new Grid2D(image.getWidth(), image.getHeight()); + + for(int i = 0; i < image.getHeight(); i++) + { + for(int j = 0; j < image.getWidth(); j++) + { + imageRankK.setAtIndex(j, i, (float) IapproxK.getElement(i, j)); + } + } + imageRankK.show("Direct output from rank 1-matrices"); + + } + + public float calculateRMSE(Grid2D image1, Grid2D image2) { + NumericGridOperator op = new NumericGridOperator(); + float rmse = 0; //TODO + return rmse; + } + + public SimpleVector[] optimizationProblem4(double[] xCoords, double[] yCoords) + { + System.out.println("Optimization Problem 4"); + //%%%%%%%%%%%%%%%%%%%%%%%%%%% + //Optimization problem IV + //%%%%%%%%%%%%%%%%%%%%%%%%%%% + + // The Moore-Penrose pseudo-inverse is required to find the + // solution to the following optimization problem: + // Compute the regression line through the following 2-D points: + // We have sample points (x_i,y_i) and want to fit a line model + // y = mx + t (linear approximation) through the points + // y = A*b + // Note: y = m*x +t*1 Thus: A contains (x_i 1) in each row, and y contains the y_i coordinates + // We need to solve for b = (m t) + + SimpleMatrix A4 = new SimpleMatrix(xCoords.length, 2); + SimpleVector aCol = new SimpleVector(xCoords.length); + SimpleVector y = new SimpleVector(xCoords.length); + + for(int i = 0; i < xCoords.length; i++)//TODO + { + aCol.setElementValue(i, 0);//TODO + y.setElementValue(i, 0);//TODO + } + + A4.setColValue(0, aCol); + SimpleVector aColOnes = new SimpleVector(xCoords.length); + aColOnes.ones(); + A4.setColValue(1, aColOnes); + + // Note: We need to use SVD matrix inversion because the matrix is not square + // more samples (7) than unknowns (2) -> overdetermined system -> least-square + // solution. + SimpleMatrix Ainv = A4.inverse(InversionType.INVERT_SVD); + SimpleVector b = SimpleOperators.multiply(Ainv, y); + + LinearFunction lFunc = new LinearFunction(); + lFunc.setM(b.getElement(0)); + lFunc.setT(b.getElement(1)); + VisualizationUtil.createScatterPlot(xCoords, yCoords, lFunc).show(); + + SimpleVector[] params = new SimpleVector[2]; + params[0] = new SimpleVector(aCol); + params[1] = new SimpleVector(y); + return params; + } + + /** + * The following content is read-only. + * Both methods + * 1. conditionNumber(...) + * 2. optimizationProblem2(...) + * 3. pseudoInverse(...) + * are used above to + * 1. calculate the condition number of a matrix + * 2. estimate the matrix A that minimizes sum_{i=1}^{4} b'_{i} A b_{i} + * under the constraint ||A||_{F} = 1 + * 3. creates the pseudo inverse of a Matrix + * You can try to understand it if you want to see how these operators can be implemented. + * Otherwise the exercise is finished here. + */ + + public void conditionNumber(SimpleMatrix A){ + + DecompositionSVD svd = new DecompositionSVD(A); + double cond = svd.cond(); + System.out.println("Cond(A) = " + cond); + + } + + public void optimizationProblem2(SimpleMatrix b){ + + System.out.println("Optimization Problem 2"); + // Estimate the matrix A in R^{2,2} such that for the following vectors + // the optimization problem gets solved: + // sum_{i=1}^{4} b'_{i} A b_{i} and ||A||_{F} = 1 + // The objective function is linear in the components of A, thus the whole + // sum can be rewritten in matrix notation: + // Ma = 0 + // where the measurement matrix M is built from single elements of the + // sum. + // Given the four points we get four equations and therefore four rows in M: + + SimpleMatrix M = new SimpleMatrix(b.getCols(), 4); + + for(int i = 0; i < b.getCols(); i++) + { + // i-th column of b contains a vector. First entry of M is x^2 + M.setElementValue(i, 0, Math.pow(b.getElement(0, i), 2.f)); + M.setElementValue(i, 1, b.getElement(0, i)*b.getElement(1, i)); + M.setElementValue(i, 2, M.getElement(i, 1)); + M.setElementValue(i, 3, Math.pow(b.getElement(1, i), 2.f)); + } + + // TASK: estimate the matrix A + // SOLUTION: Compute vector that spans the nullspace of A; we can find it in + // the last column vector of V. + + DecompositionSVD svd = new DecompositionSVD(M); + int lastCol = svd.getV().getCols() -1; + SimpleVector a = svd.getV().getCol(lastCol); + + //We need to reshape the vector back to the desired 2x2 matrix form + SimpleMatrix A2 = new SimpleMatrix(2,2); + A2.setColValue(0, a.getSubVec(0, 2)); + A2.setColValue(1, a.getSubVec(2, 2)); + + //check if Frobenius norm is 1.0 + double normF = A2.norm(MatrixNormType.MAT_NORM_FROBENIUS); + System.out.println("||A2||_F = " + normF); + + //check solution + + SimpleVector temp = SimpleOperators.multiply(M, a); + double result = SimpleOperators.multiplyInnerProd(a, temp); + System.out.println("Minimized error: " + result); + } + + public SimpleMatrix pseudoInverse(SimpleMatrix A){ + + //Moore-Penrose Pseudoinverse defined as V * Sinv * U^T + //Compute the inverse of the singular matrix S + + DecompositionSVD svd = new DecompositionSVD(A); + SimpleMatrix Sinv = new SimpleMatrix( svd.getS().getRows(), svd.getS().getCols()); + Sinv.zeros(); + + int size = Math.min(Sinv.getCols(), Sinv.getRows()); + SimpleVector SinvDiag = new SimpleVector( size); + + for(int i = 0; i < Sinv.getRows(); i++) + { + double val = 1.0 / svd.getS().getElement(i, i); + SinvDiag.setElementValue(i, val); + } + + Sinv.setDiagValue(SinvDiag); + + //Compare our implementation to svd.getreciprocalS() + System.out.println("Sinv = " + Sinv.toString()); + System.out.println("Srec = " +svd.getreciprocalS().toString()); + + SimpleMatrix tempInv = SimpleOperators.multiplyMatrixProd(svd.getV(), svd.getreciprocalS()); + SimpleMatrix Ainv = SimpleOperators.multiplyMatrixProd(tempInv, svd.getU().transposed()); + System.out.println("Ainv = " + Ainv.toString()); + System.out.println("A.inverse() = " + A.inverse(InversionType.INVERT_SVD)); + return Ainv; + } + + + /** + * + * end of the exercise + */ + + // getters for members + // variables which are checked (DO NOT CHANGE!) + // + public float[] get_error() { + return error; + } + // + } diff --git a/src/edu/stanford/rsl/tutorial/mipda/ExerciseUM.java b/src/edu/stanford/rsl/tutorial/mipda/ExerciseUM.java new file mode 100644 index 00000000..c64b1ab5 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/ExerciseUM.java @@ -0,0 +1,452 @@ +package edu.stanford.rsl.tutorial.mipda; + +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.data.numeric.Grid2DComplex; +import edu.stanford.rsl.conrad.data.numeric.NumericPointwiseOperators; +import edu.stanford.rsl.conrad.filtering.MeanFilteringTool; +import edu.stanford.rsl.conrad.utils.ImageUtil; +import ij.IJ; +import ij.ImageJ; +import ij.ImagePlus; +import ij.process.FloatProcessor; + +/** + * Unsharp Masking + * Programming exercise for module "MR Inhomogeneities" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Ashwini Jadhav, Anna Gebhard, Mena Abdelmalek + * + */ + +public class ExerciseUM { + + Grid2D originalImage; + Grid2D biasedImage; + Grid2D correctedHUM; + Grid2D correctedF; + Grid2D correctedCut; + public static final String example = "head"; // select "box" or "head" + + public static void main(String[] args) { + + /* instantiate exercise object */ + ExerciseUM exObject = new ExerciseUM(example); + + /* start ImageJ */ + // always start an instance of ImageJ for image analysis + new ImageJ(); + + Grid2D image = exObject.originalImage; + Grid2D biased = exObject.biasedImage; + + image.show("Original image"); + biased.show("Biased image with noise"); + + /** + * preprocessing + */ + + + ImagePlus imp = new ImagePlus("",ImageUtil.wrapGrid2D(biased)); + IJ.run(imp, "Gaussian Blur...", "sigma=1.5"); // imp works on image directly + + exObject.correctedHUM = extClone(biased); // remark: extClone() is implemented only for this exercise! + exObject.correctedF = extClone(biased); + exObject.correctedCut = extClone(biased); + + + /** + * implement homomorphic unsharp masking + */ + int kernelSize = 31; + + exObject.correctedHUM = exObject.homomorphicFiltering(biased,kernelSize); // TODO: Go to this method and complete it! + exObject.correctedHUM.show("Corrected image using homomorphic unsharp masking"); + + /** + * frequency based correction + */ + double beta = 1.0d; + double sigma = 0.01d; + + exObject.correctedF = exObject.frequencyCorrection(biased, beta, sigma); // TODO: Go to this method and complete it! + exObject.correctedF.show("Corrected image using frequency filtering"); + + /** + * for comparison: hard frequency cutoff + */ + double cutoff = 0.01d; + + exObject.correctedCut = exObject.frequencyCutoff(exObject.correctedCut,cutoff); // TODO: Go to this method and complete it! + if (exObject.correctedCut != null) + exObject.correctedCut.show("Corrected image with frequency cutoff"); + + IJ.run(imp, "Window/Level...", ""); // open tool to window the images properly + } + + public Grid2D homomorphicFiltering(Grid2D grid, int kernelSize){ + + /* + * homomorphic unsharp masking: + * - compute global and local means + * - look for the Conrad API filtering tool MeanFilteringTool and use it to compute the mean values + * - weight the image by mu/mu_ij + */ + + Grid2D corrected = (Grid2D) grid.clone(); + + Grid2D mu_ij = null; //TODO: clone grid + + MeanFilteringTool localMeanFilter = new MeanFilteringTool(); + try { + //TODO // use configure() to configure localMeanFilter + //TODO // use the filter and applyToolToImage() such that mu_ij contains the result image of grid being average-filtered + } catch (Exception e) { + e.printStackTrace(); + } + + float mu = this.meanPositives(grid); + + float eps = 0.01f; + float tempVal = 0; + + for (int i = 0; i < corrected.getHeight(); i++){ + for (int j = 0; j < corrected.getWidth(); j++){ + + if (mu_ij==null) break; + // this is a workaround to avoid the quotient to get too large or too small + if (mu_ij.getAtIndex(j, i) <= eps*mu){ + tempVal = 1.0f/eps; + } + else if (mu_ij.getAtIndex(j, i) >= mu/eps){ + tempVal = eps; + } + else{ + tempVal = 0; // TODO + } + + // TODO: apply the bias field correction for the multiplicative model + } + } + + return corrected; + } + + public Grid2D frequencyCorrection(Grid2D grid, double beta, double sigma){ + + Grid2D corrected = (Grid2D) grid.clone(); + Grid2DComplex fftimage = new Grid2DComplex(grid); // check the Grid size and look up zero-padding + + // transform into Fourier domain + // TODO: Fourier Forward Transform + // TODO: Shift origin point to the center + + for (int k = 0; k < fftimage.getHeight(); k++){ + for (int l = 0; l < fftimage.getWidth(); l++){ + + double fk = (k-fftimage.getHeight()/2.0d)/fftimage.getWidth()/grid.getSpacing()[0]; + double fl = (l-fftimage.getWidth()/2.0d)/fftimage.getHeight()/grid.getSpacing()[1]; + + /* Apply the reverse Gaussian filter from the lecture */ + float val = 0; // TODO + // TODO: setRealAtIndex for fftimage + // TODO: setImageAtIndex for fftimage + } + } + + // transform back into spatial domain + // TODO: inverse shift origin for fftimage + // TODO: inverse Fourier Transform for fftimage + + for (int i = 0; i < corrected.getHeight(); i++){ + for (int j = 0; j < corrected.getWidth(); j++){ + corrected.setAtIndex(j, i, fftimage.getRealAtIndex(j, i)); + } + } + return corrected; + } + + public Grid2D frequencyCutoff(Grid2D grid,double cutoff) { + + Grid2D returnGrid = new Grid2D(grid); // TODO: apply a hard frequency cut-off (high pass filter), ... + returnGrid = null; // ... look for the appropriate class method in this file + + return returnGrid; + } + + public Grid2D lowPass(Grid2D grid, double cutoff){ + return passFilter(grid,true,cutoff); + } + + public Grid2D highPass(Grid2D grid, double cutoff){ + return passFilter(grid,false,cutoff); + } + + + + /** + * + * end of the exercise + */ + + // further class methods (DO NOT CHANGE!) + public Grid2D readImageFile(String filename){ + +// return ImageUtil.wrapImagePlus(IJ.openImage(getClass().getResource(filename).getPath())).getSubGrid(0); + return ImageUtil.wrapImagePlus(IJ.openImage(filename)).getSubGrid(0); + } + + private void simulateBiasedImage(Grid2D image, Grid2D biasField){ + NumericPointwiseOperators.multiplyBy(image, biasField); + } + + public void addNoise(Grid2D image, double noiseLevel){ + + FloatProcessor fp; + + fp = ImageUtil.wrapGrid2D(image); + fp.noise(noiseLevel*(fp.getMax()-fp.getMin())); // direct change of input image + } + + private float meanPositives(Grid2D grid){ + + int numNonZeros = 0; + float val = 0; + + for (int i = 0; i < grid.getHeight(); i++){ + for (int j = 0; j < grid.getWidth(); j++){ + if (grid.getAtIndex(j, i) > 0){ + val += grid.getAtIndex(j, i); + numNonZeros++; + } + } + } + + val/=numNonZeros; + + return val; + } + + private void generateExerciseGainField(Grid2D gainField){ + + float upPerc = 0.2f; + float lowPerc = 5.0f; + + for (int i = 0; i < gainField.getHeight(); i++){ + for (int j = 0; j < gainField.getWidth(); j++){ + + double[] coords = gainField.indexToPhysical(j, i); + + double fx = 0.5d*coords[0]/gainField.getWidth()*Math.PI; + double fy = 0.5d*coords[1]/gainField.getHeight()*Math.PI; + + float val = (float) (0.5*fx*fy + 0.3*fx*fx*fx + 0.7*fy*fy + 0.1*fx + 0.4*fy); + + gainField.setAtIndex(j, i, val); + } + } + + scaleGrid(gainField, upPerc, lowPerc); + } + + private void placeBoxOnGrid(Grid2D grid,double magnitude,double[] mu,double[] bds){ + + for (int i = 0; i < grid.getHeight(); i++){ + for (int j = 0; j < grid.getWidth(); j++){ + + double[] coords = grid.indexToPhysical(j,i); // get actual position (not indices) + + if (coords[0]>-bds[0]+mu[0] && coords[0]-bds[1]+mu[1] && coords[1] low pass + // type=false => high pass + + Grid2DComplex fftGrid = new Grid2DComplex(grid); + Grid2D filteredGrid = extClone(grid); + + fftGrid.transformForward(); // does not know of spacing or origin + fftGrid.setSpacing(1.0f/fftGrid.getWidth()/grid.getSpacing()[0],1.0f/fftGrid.getHeight()/grid.getSpacing()[1]); + fftGrid.fftshift(); + +// fftGrid.getMagnSubGrid(0,0,fftGrid.getWidth(),fftGrid.getHeight()).show("Magnitude of FFT image"); +// fftGrid.getPhaseSubGrid(0,0,fftGrid.getWidth(),fftGrid.getHeight()).show("Phase of FFT image"); + + for (int i = 0; i < fftGrid.getHeight(); i++){ + for (int j = 0; j < fftGrid.getWidth(); j++){ + + double[] frequencies = + new double[]{(i-fftGrid.getHeight()/2.0d)*fftGrid.getSpacing()[0], + (j-fftGrid.getWidth()/2.0d)*fftGrid.getSpacing()[1]}; + if (type){ + if (Math.pow(frequencies[0],2.0d) + Math.pow(frequencies[1],2.0d) > Math.pow(cutoff,2.0d)){ + fftGrid.setRealAtIndex(j, i, 0.0f); + fftGrid.setImagAtIndex(j, i, 0.0f); + } + } + else{ + if (Math.pow(frequencies[0],2.0d) + Math.pow(frequencies[1],2.0d) < Math.pow(cutoff,2.0d)){ + fftGrid.setRealAtIndex(j, i, 0.0f); + fftGrid.setImagAtIndex(j, i, 0.0f); + } + } + fftGrid.setImagAtIndex(j, i, 0.0f); + } + } + +// fftGrid.getMagnSubGrid(0,0,fftGrid.getWidth(),fftGrid.getHeight()).show("Magnitude of filtered FFT image"); + + fftGrid.ifftshift(); + fftGrid.transformInverse(); + + NumericPointwiseOperators.copy(filteredGrid,fftGrid.getRealSubGrid(0,0,grid.getWidth(),grid.getHeight())); + + return filteredGrid; + } + + private void normalize(Grid2D grid){ + + // Normalize intensity values to [0,1] + scaleGrid(grid, 0.0f, 1.0f); + } + + private void scaleGrid(Grid2D grid, float newMin, float newMax){ + + float max = NumericPointwiseOperators.max(grid); + float min = NumericPointwiseOperators.min(grid); + + for(int i = 0; i < grid.getHeight(); i++){ + for(int j = 0; j < grid.getWidth(); j++){ + + float scaledValue = (grid.getAtIndex(j, i) - min) / (max-min) * (newMax-newMin) + newMin; + grid.setAtIndex(j, i, scaledValue); + } + } + } + + public ExerciseUM(String example){ + + /* image parameters */ + // define virtual spacing + float imageXspacing = 1.0f; + float imageYspacing = 1.0f; + + // noise (Gaussian) + double noiseLevel = 0.01d; // percentage + + /* image container */ + Grid2D gainField; + + /* specify image container properties */ + switch (example) { + + case "head": + + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + + originalImage = readImageFile(imageDataLoc + "mr_head_dorsal.jpg"); + + originalImage.setSpacing(imageXspacing,imageYspacing); + // move origin to center + originalImage.setOrigin(-originalImage.getWidth()/2.0d*originalImage.getSpacing()[0], + -originalImage.getHeight()/2.0d*originalImage.getSpacing()[1]); + + break; + + default: + originalImage = new Grid2D(256,386); + + originalImage.setSpacing(imageXspacing,imageYspacing); + // move origin to center + originalImage.setOrigin(-originalImage.getWidth()/2.0d*originalImage.getSpacing()[0], + -originalImage.getHeight()/2.0d*originalImage.getSpacing()[1]); + + placeBoxOnGrid(originalImage,1.0d,new double[]{20.0d,30.0d},new double[]{30.0d,50.0d}); + + break; + } + + // Normalize intensity values to [0,1] + normalize(originalImage); + + // copy settings + gainField = extClone(originalImage); + + /* generate biased noisy image */ + generateExerciseGainField(gainField); + + biasedImage = extClone(originalImage); + simulateBiasedImage(biasedImage, gainField); + + addNoise(biasedImage,noiseLevel); + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public Grid2D get_originalImage() { + return originalImage; + } + // + public Grid2D get_biasedImage() { + return biasedImage; + } + // + public Grid2D get_correctedHUM() { + return correctedHUM; + } + // + public Grid2D get_correctedF() { + return correctedF; + } + // + public Grid2D get_correctedCut() { + return correctedCut; + } +} diff --git a/src/edu/stanford/rsl/tutorial/mipda/Intro.java b/src/edu/stanford/rsl/tutorial/mipda/Intro.java new file mode 100644 index 00000000..b44c7314 --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/Intro.java @@ -0,0 +1,485 @@ +package edu.stanford.rsl.tutorial.mipda; + +import edu.stanford.rsl.conrad.data.numeric.Grid2D; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix; +import edu.stanford.rsl.conrad.numerics.SimpleMatrix.MatrixNormType; +import edu.stanford.rsl.conrad.numerics.SimpleVector.VectorNormType; +import edu.stanford.rsl.conrad.utils.ImageUtil; +import edu.stanford.rsl.conrad.utils.VisualizationUtil; +import edu.stanford.rsl.conrad.numerics.SimpleOperators; +import edu.stanford.rsl.conrad.numerics.SimpleVector; +import ij.IJ; +import ij.ImageJ; +import ij.plugin.filter.Convolver; +import ij.process.FloatProcessor; + + +/** + * Introduction to the CONRAD Framework + * Programming exercise for module "Course Introduction" + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch, Marco Boegel, Ashwini Jadhav + * + */ + +public class Intro { + + // variables which are checked (DO NOT CHANGE!) + SimpleVector v1; + SimpleVector vRand; + SimpleMatrix M; + double Mdeterminant; + SimpleMatrix Mtrans; + SimpleMatrix Mcopy; + int numRows; + int numCols; + String elementsOutput = ""; + SimpleMatrix Mones; + SimpleMatrix Mzeros; + SimpleMatrix Midentity; + SimpleMatrix ResMat; + SimpleVector resVec; + SimpleVector subVector; + SimpleMatrix MsquaredElem ; + double matrixNormL1; + double vecNormL2; + SimpleVector v2; + SimpleVector tempColVector; + int plotLength = 500; + double[] y = new double[plotLength]; + double[] x = new double [plotLength]; + int imageSizeX = 256, imageSizeY = 256; + Grid2D image; + Grid2D copy; + String imageDataLoc = System.getProperty("user.dir") + "/data/" + "/mipda/"; + String filename = imageDataLoc + "mr12.dcm"; + Grid2D mrImage; + Convolver conv; + boolean convolution; + Grid2D convolvedImage; + String imageFileName; + String outFileName; + + /** + * You can change this to a location OUTSIDE of the Conrad workspace, where you want to save the result to. + * If you leave it as it is, you find the result in your home directory. + */ + String outputDataLoc = System.getProperty("user.home") + "/mipda/output/"; + + + /** + * In the basicIntro routine you find out how matrix and vector computations + * can be performed. Follow the code and fill in missing code indicated by TODOs. + * There are either empty lines which need exactly one line of code, + * or you have to initialize declared variables appropriately. + */ + public void basicIntro() { + + //display text + System.out.println("Creating a vector: v1 = [3.0; 2.0; 1.0]"); + + //create column vector with the entries 3.0, 2.0 and 1.0 as floats + v1 = null; //TODO + if (v1 != null) + System.out.println("v1 = " + v1.toString()); + + //create a randomly initialized vector with values between 0 and 1 + vRand = new SimpleVector(3); + //TODO + System.out.println("vRand = " + vRand.toString()); + + //create a 3x3 matrix M with row vectors (1, 2, 3), (4, 5, 6) and (7, 8, 9) + //hint: have a look at the method setColValue(int, SimpleVector) + M = null;//TODO + //TODO + //TODO + //TODO + if (M != null) + System.out.println("M = " + M.toString()); + + //Mdeterminant: compute the determinant of M + Mdeterminant = 0.0; //TODO + + System.out.println("Determinant of matrix m: " + Mdeterminant); + + //transpose M + Mtrans = null;//TODO + //copy matrix using copy constructor + Mcopy = null;//TODO + //transpose Mcopy in-place + //TODO + //get size + numRows = 0;//TODO + numCols = 0;//TODO + + //access and print the elements of M (the original matrix from above) + System.out.println("M: "); + + for(int i = 0 ; i < numRows; i++) + { + for(int j = 0; j < numCols; j++) + { + double element = 0;//TODO + + elementsOutput = elementsOutput + element + " "; + + } + elementsOutput += "\n"; + } + System.out.print(elementsOutput); + + //create a 3x3 matrix Mones of ones (all entries are 1) + Mones = new SimpleMatrix(3,3); + //TODO + //create a 3x3 matrix Mzeros of zeros (all entries are 0) + Mzeros = new SimpleMatrix(3,3); + //TODO + //create a 3x3 identity matrix + Midentity = new SimpleMatrix(3,3); + //TODO + + //matrix multiplication + //compute the matrix product of Mtrans and M + //hint: have a look at the class edu.stanford.rsl.conrad.numerics.SimpleOperators + ResMat = null;//TODO + if (ResMat != null) + System.out.println("M^T * M = " + ResMat.toString()); + + //matrix vector multiplication + //compute the matrix vector product of M and v1 + // -> SimpleOperators + resVec = null;//TODO + if (resVec != null) + System.out.println("M * v1 = " + resVec.toString()); + + //extract the last column vector from matrix M + SimpleVector colVector = null; + if (M != null) { + if (M.getCols() > 2) { + colVector = M.getCol(2); + } + } + + //extract the top 1x2 subvector from the last column of matrix M + subVector = null;//TODO + if (subVector != null) + System.out.println("[m(0)(2); m(1)(2)] = " + subVector); + + //matrix elementwise multiplication + //compute a matrix which has the squared value of the elements of M in each component + MsquaredElem = null;//TODO + if (MsquaredElem != null) + System.out.println("M squared Elements: " + MsquaredElem.toString()); + + //examples how to round vectors + // copy of the random vector vRand + SimpleVector vRandCopy = new SimpleVector(vRand); + System.out.println("vRand = " + vRandCopy.toString()); + // bring down to a round figure + vRandCopy.floor(); + System.out.println("vRand.floor() = " + vRandCopy.toString()); + // round up + vRand.ceil(); + System.out.println("vRand.ceil() = " + vRand.toString()); + + //compute min and max values + // for vectors + if (v1 != null) { + double minV1 = v1.min(); + double maxV1 = v1.max(); + System.out.println("Min(v1) = " + minV1 + " Max(v1) = " + maxV1); + } + // for matrices: iterate over row or column vectors + if (M != null) { + SimpleVector maxVec = new SimpleVector(M.getCols()); + for(int i = 0; i < M.getCols(); i++) + { + maxVec.setElementValue(i, M.getCol(i).max()); + } + double maxM = 0; //TODO: compute total maximum + System.out.println("Max(M) = " + maxM); + } + + //norms + matrixNormL1 = 0; //TODO Frobenius norm of M + vecNormL2 = 0; //TODO L2 vector norm of colVector + System.out.println("||M||_F = " + matrixNormL1); + System.out.println("||colVec||_2 = " + vecNormL2); + + //get the normalized vector from colVector without overwriting + v2 = null;//TODO + + if (colVector != null) + tempColVector = new SimpleVector(colVector); //this copy is necessary for unit testing + else + tempColVector = new SimpleVector(2); + + //normalize tempColVector vector in-place (overwrite it with normalized values) + //TODO + if (v2 != null) + System.out.println("Normalized colVector: " + v2.toString()); + System.out.println("||colVec||_2 = " + tempColVector.norm(VectorNormType.VEC_NORM_L2)); // if done correctly this yields 1 + } + + + /** + * The signalIntro routine shows how a sine function can be plotted. + * Follow the code and fill in missing code indicated by TODOs. + * There are either empty lines which need exactly one line of code, + * or you have to initialize declared variables appropriately. + */ + public void signalIntro() + { + //How is the sine function sin(2*PI*x) plotted using this framework? + // use Math.sin and Math.PI + double stepSize = 0.01; + //compute the image of sin(2*PI*x) using the given step size stepSize and starting at the origin + for(int i = 0; i < y.length; i++) + { + y[i] = 0.0; //TODO + + } + VisualizationUtil.createPlot(y).show(); + + // now plot it with the specified x values + // (such that it does not show the loop index but the actual x-values on the x-axis) + // x is increased after every iteration by stepSize + for(int i = 0; i < x.length; i++) + { + x[i] = 0.0; //TODO + } + + VisualizationUtil.createPlot(x, y, "sin(x)", "x", "y").show(); + } + + /** + * The gridIntro routine contains code for image manipulation with ImageJ. + * First a simple image is created which shows a sphere. + * Next, we load an image from the disk, implement an average filter, and save the result. + * Follow the code and fill in missing code indicated by TODOs. + * There are either empty lines to be filled or you have to initialize declared variables appropriately. + */ + public void gridIntro(){ + + //define the image size + int imageSizeX = 256; + int imageSizeY = 256; + + //define an image + //hint: use the package edu.stanford.rsl.conrad.data.numeric.Grid2D + image = null;//TODO + + //draw a filled circle + int radius = 50; + //set all pixels within the circle to 100 + int insideVal = 100; + + // fill 'image' with data that shows a sphere with radius 'radius' and intensity 'insideVal' + //TODO (multiple code lines, hint: two for loops) + + //show ImageJ GUI + ImageJ ij = new ImageJ(); + //display image using the Grid2D class methods + //TODO + + //copy an image + copy = null; //TODO + if (copy != null) + copy.show("Copied image of a circle"); + + //load an image from file + // first use IJ to open the image + // then wrap it using the class ImageUtil (static) + // finally you need to extract the image with getSubGrid(0) + mrImage = null;//TODO + if (mrImage != null) + mrImage.show("MR image"); + + // learn how to compute the convolution of an image with a filter kernel + conv = new Convolver(); + FloatProcessor ip = null; + if (mrImage != null) + ip = ImageUtil.wrapGrid2D(mrImage); + + //this defines a simple averaging 3x3 filter kernel + int kw = 3; + int kh = 3; + float[] kernel = new float[kw*kh]; + + for(int i = 0; i < kernel.length; i++) + { + kernel[i] = 1.f / (kw*kh); + } + + // trigger the convolution using ip, kernel, kw, kh + convolution = false;//TODO + if (ip != null) { + convolvedImage = ImageUtil.wrapFloatProcessor(ip); + convolvedImage.show("Convolved Image"); + } + + + //write an image to disk & check the supported output formats + imageFileName = "mr12out.tiff"; // TODO: choose a valid filename (just file+ending, no path) + + outFileName = outputDataLoc + imageFileName; + if (mrImage != null) { + //TODO: save the image using IJ and ImageUtil + } + } + + /** + * + * end of the exercise + */ + + + public static void main(String[] args) { + + (new Intro()).basicIntro(); + (new Intro()).signalIntro(); + (new Intro()).gridIntro(); + } + + // getters for members + // variables which are checked (DO NOT CHANGE!) + public SimpleVector get_v1() { + return v1; + } + // + public SimpleVector get_vRand() { + return vRand; + } + // + public SimpleMatrix get_M() { + return M; + } + // + public double get_Mdeterminant() { + return Mdeterminant; + } + // + public SimpleMatrix get_Mtrans() { + return Mtrans; + } + // + public SimpleMatrix get_Mcopy() { + return Mcopy; + } + // + public int get_numRows() { + return numRows; + } + // + public int get_numCols() { + return numCols; + } + // + public String get_elementsOutput() { + return elementsOutput; + } + // + public SimpleMatrix get_Mones() { + return Mones; + } + // + public SimpleMatrix get_Mzeros() { + return Mzeros; + } + // + public SimpleMatrix get_Midentity() { + return Midentity; + } + // + public SimpleMatrix get_ResMat() { + return ResMat; + } + // + public SimpleVector get_resVec() { + return resVec; + } + // + public SimpleVector get_subVector() { + return subVector; + } + // + public SimpleMatrix get_MsquaredElem() { + return MsquaredElem; + } + // + public double get_matrixNormL1() { + return matrixNormL1; + } + // + public double get_vecNormL2() { + return vecNormL2; + } + // + public SimpleVector get_v2() { + return v2; + } + // + public SimpleVector get_tempColVector() { + return tempColVector; + } + // + public int get_plotLength() { + return plotLength; + } + // + public double[] get_y() { + return y; + } + // + public double[] get_x() { + return x; + } + // + public int get_imageSizeX(){ + return imageSizeX; + } + // + public int get_imageSizeY() { + return imageSizeY; + } + // + public Grid2D get_image() { + return image; + } + // + public Grid2D get_copy() { + return copy; + } + // + public String get_imageDataLoc() { + return imageDataLoc; + } + // + public String get_filename() { + return filename; + } + // + public Grid2D get_mrImage() { + return mrImage; + } + // + public Convolver get_conv() { + return conv; + } + // + public boolean get_convolution() { + return convolution; + } + // + public Grid2D get_convolvedImage() { + return convolvedImage; + } + // + public String get_imageFileName() { + return imageFileName; + } + // + public String get_outFileName() { + return outFileName; + } +} diff --git a/src/edu/stanford/rsl/tutorial/mipda/RunExerciseEvaluation.java b/src/edu/stanford/rsl/tutorial/mipda/RunExerciseEvaluation.java new file mode 100644 index 00000000..6ab5b24e --- /dev/null +++ b/src/edu/stanford/rsl/tutorial/mipda/RunExerciseEvaluation.java @@ -0,0 +1,137 @@ +package edu.stanford.rsl.tutorial.mipda; + +import java.awt.EventQueue; + +import javax.swing.DefaultComboBoxModel; +import javax.swing.JComboBox; +import javax.swing.JLabel; +import javax.swing.JOptionPane; +import javax.swing.JPanel; +import javax.swing.SwingUtilities; +import java.lang.InterruptedException; +import java.lang.reflect.InvocationTargetException; + +import org.junit.internal.TextListener; +import org.junit.runner.JUnitCore; + +import ModuleCourseIntroduction.IntroTestClass; +import ModuleMathematicalTools.SVDTestClass; +import ModulePreprocessing.GeoUTestClass; +import ModulePreprocessing.UMTestClass; +import ModuleReconstruction.PBTestClass; +import ModuleReconstruction.FBTestClass; +import ModuleRegistration.REGTestClass; +import ModuleRegistration.MITestClass; + +/** + * Testing tool for programming exercises + * of the course "Medical Image Processing for Diagnostic Applications (MIPDA)" + * @author Frank Schebesch + * + */ + +public class RunExerciseEvaluation { + + private static final String[] exerciseNames = {"Intro","SVD","GeoU","UM","PB","FB","REG","MI"}; + + public static void main(String[] args) { + + JUnitCore junitCore = new JUnitCore(); + junitCore.addListener(new TextListener(System.out)); + + String choice = ask(exerciseNames); + System.out.println("Running test " + choice); + + if (choice.equals("Intro")) { + junitCore.run(IntroTestClass.class); + //org.junit.runner.JUnitCore.main("ModuleCourseIntroduction.IntroTestClass"); + } + else if (choice.equals("SVD")) { + junitCore.run(SVDTestClass.class); + } + else if (choice.equals("GeoU")) { + junitCore.run(GeoUTestClass.class); + } + else if (choice.equals("UM")) { + junitCore.run(UMTestClass.class); + } + else if (choice.equals("PB")) { + junitCore.run(PBTestClass.class); + } + else if (choice.equals("FB")) { + junitCore.run(FBTestClass.class); + } + else if (choice.equals("REG")) { + junitCore.run(REGTestClass.class); + } + else if (choice.equals("MI")) { + junitCore.run(MITestClass.class); + } + else { + System.out.println("Aborted, no test selected."); +// System.out.println("No such test available."); + } + } + + // thanks to https://stackoverflow.com/questions/13408238/simple-java-gui-as-a-popup-window-and-drop-down-menu + // (adapted) + public static String ask(final String... values) { + + String result = null; + + if (EventQueue.isDispatchThread()) { + + JPanel panel = new JPanel(); + panel.add(new JLabel("Select the exercise test that you wish to run:")); + + DefaultComboBoxModel model = new DefaultComboBoxModel<>(); + for (String value : values) { + model.addElement(value); + } + JComboBox comboBox = new JComboBox<>(model); + panel.add(comboBox); + + int iResult = JOptionPane.showConfirmDialog(null, panel, "jUnit Exercise Test Selection", JOptionPane.OK_CANCEL_OPTION, JOptionPane.QUESTION_MESSAGE); + switch (iResult) { + case JOptionPane.OK_OPTION: + result = (String) comboBox.getSelectedItem(); + break; + case JOptionPane.OK_CANCEL_OPTION: + return ""; + } + + } else { + + Response response = new Response(values); + try { + SwingUtilities.invokeAndWait(response); + result = response.getResponse(); + } catch (InterruptedException | InvocationTargetException ex) { + ex.printStackTrace(); + } + + } + + return result; + + } + // thanks to https://stackoverflow.com/questions/13408238/simple-java-gui-as-a-popup-window-and-drop-down-menu + public static class Response implements Runnable { + + private String[] values; + private String response; + + public Response(String... values) { + this.values = values; + } + + @Override + public void run() { + response = ask(values); + } + + public String getResponse() { + return response; + } + } +}