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);
}
*/
}