Orthonormal Basis Transforms
Grant Schindler, 2007

 Click numbered buttons to change active image. Click & drag slider to change # basis vectors used. Click tabs to change active basis.

 Description Any image can be represented as a combination of simpler "basis" images. Normally, these basis images each consist of a single active pixel -- this describes the "Identity" basis -- but other sets of basis images, including sine waves (Fourier Analysis) and eigenvectors (Principal Components Analysis or PCA), can be chosen to describe the same image. In the example above, each row of the Transformed Image tells exactly how to combine the rows of the Orthogonal Basis in order to form the rows of the Reconstructed Image. The above images are visual representations of matrices where black = -1.0 and white = 1.0, while gray values lie somewhere between. Watch the image on the left slowly take shape as you add basis vectors to the reconstruction: Drag the slider all the way left and then use the '<' and '>' keys to add or remove one basis vector at a time. Do the same for each basis: Identity, PCA, Fourier, Haar, and Random. Use keys 'i','p','f','h', and 'r' to switch modes. Use number keys to switch images. Motivation After reading the statement "As is the case with SVD, Fourier analysis involves expansion of the original data in an orthogonal basis." in the paper Singular Value Decomposition and Principal Component Analysis, I set out to demonstrate this fact to myself and to understand its broader implications for all orthonormal bases. In physics and mathematics, such orthonormal basis transforms are called integral transforms. Technical Details The most important idea here is that the reconstructed image on the left is literally the product of the two matrices on the right: the transformed image times the basis equals the reconstructed image. The math is the same no matter what set of basis vectors we choose. We compute each set of basis vectors as follows: PCA: The eigenvectors of the original image are computed via power iteration. Just 3 iterations per eigenvector suffices. PCA is most often carried out with a toolbox Singular Value Decomposition (SVD) implementation which, while efficient, masks the simplicity with which eigenvectors can be extracted from data using power iteration. (Note: In this context, PCA is perhaps more properly called the Karhunen-Loeve Transform.) Fourier: Computed as per the Discrete Cosine Transform (DCT II). Haar: Computed as per 1D Haar Wavelets. Random: A random orthonormal basis is built up incrementally. At each stage, we append to the basis the normalized difference between a random vector and that vector's reconstruction using the basis thus far. The logic is that any part of a random vector unreachable by the current basis must, by definition, be orthogonal to it. Why orthonormal basis transforms and not just orthogonal basis transforms? If an image is transformed by multiplying it with a matrix, then the transform can be undone by multiplying the result with the inverse of the matrix. In the case of an orthonormal basis (having vectors of unit length), the inverse is just the transpose of the matrix. Thus, inverting an orthonormal basis transform is a trivial operation. Source Code Source code is provided below for educational purposes. Complete source files (including images) for the project are here: transforms.zip (152KB). Built with Processing. References Michael E. Wall, Andreas Rechtsteiner, Luis M. Rocha. Singular Value Decomposition and Principal Component Analysis. 2003.

 `//Orthogonal Basis Transforms//Grant Schindler, 2007int N = 256, H = 192; //Image Dimensionsfloat[][] T = new float[N][N]; //Transform Matrixfloat[][] T_; //Inverse Transform Matrixfloat[][] X = new float[H][N]; //Transformed Imagefloat[][] X_ = new float[H][N]; //Transformed Image with Zeroed Elementsfloat[][] I; //Original Imagefloat[][] I_ = new float[H][N]; //Reconstructed ImagePImage imgI,imgI_,imgT,imgX; //Displayable Image for Each Matrix//-------------------------------------------------------------------------------------// Orthonormal Basis Transform: 2 Simple Functions ------------------------------------//-------------------------------------------------------------------------------------//Transform - Use T to Compute Coefficients for Each Image Rowvoid Transform(){ for (int y=0; y < H; y++) X[y] = MatrixVectorMultiply(T,I[y],N,N);}//InverseTransform - Reconstruct Original Image by Inverting Transformvoid InverseTransform(){ for (int y=0; y < H; y++) I_[y] = MatrixVectorMultiply(T_,X_[y],N,N);}//-------------------------------------------------------------------------------------// Discrete Cosine Transform (Fourier Variant) ----------------------------------------//-------------------------------------------------------------------------------------//Create Matrix T of Cosine Wavesvoid GetCosineMatrix(){ for (int k=0; k < N; k++){ for (int n=0; n < N; n++){ T[k][n] = cos(PI * k * (n+0.5)/N); T[k][n] *= sqrt(2.0/N); //Make T orthogonal (inverse=transpose) if (k==0) T[k][n] *= 1.0/sqrt(2.0); //Make T orthogonal (inverse=transpose) } } T_ = Transpose(T,N,N); //Create Transpose/Inverse to Invert Transform Later }//-------------------------------------------------------------------------------------// Haar Wavelets ----------------------------------------------------------------------//-------------------------------------------------------------------------------------//Create Matrix T of Haar Waveletsvoid GetHaarMatrix(){ for (int k=0; k < N; k++){ float m = 8.0 - floor(log(max(1,k))/log(2)); //Level float o = k - pow(2.0, 8.0 - m); //Offset for (int n=0; n < N; n++){ T[k][n] = Haar(pow(2.0,-m)*n - o); if (k==0) T[k][n] = 1.0; T[k][n] *= pow(2.0,-m/2); //Make T orthogonal (inverse=transpose) } } T_ = Transpose(T,N,N); //Create Transpose/Inverse to Invert Transform Later }float Haar(float x) { if (x >= 0.0 && x < 0.5) {return 1.0;} else if (x >= 0.5 && x < 1.0) {return -1.0;} else {return 0.0;}}//-------------------------------------------------------------------------------------// Identity Matrix as Basis -----------------------------------------------------------//-------------------------------------------------------------------------------------//Create Identity Matrixvoid GetIdentity(){ T = Identity(N); T_ = Transpose(T,N,N); //Create Transpose/Inverse to Invert Transform Later }//-------------------------------------------------------------------------------------// Simple PCA - Power Iteration -> Eigenvectors ---------------------------------------//-------------------------------------------------------------------------------------void GetEigenvectors(boolean truncate){ //Create Eigenvalue Matrix - Set to All Zeros float[][] V = Zeros(N,N); //Construct Covariance Matrix (I-Transpose x I) - Has Same Eigenvectors as I float[][] A = MatrixMultiply(Transpose(I,H,N),I,N,H); float[] b = new float[N]; float Av = 1.0, lambda = 1.0; //Compute Largest Eigenvector, Subtract, Iterate for (int e = 0; e < N; e++){ //Randomly Initialize Eigenvector b = VectorRandom(N); //Iteratively Multiply Matrix A by Vector x - Coverges to Largest Eigenvector of Columns of A for (int i = 0; i < 3; i++){ //Amazing that 3 iterations is enough b = MatrixVectorMultiply(A,b,N,N); Av = b[0]; //A*v = lambda*v -> lambda = (A*v)[0]/v[0]; b = VectorNormalize(b,N); } //Copy New Eigenvector into Matrix arraycopy(b,V[e]); //Project Each Column of A onto Eigenvector to Get Coefficients float[] c = MatrixVectorMultiply(Transpose(A,N,N),b,N,N); //Reconstruct Covariance Matrix Using Only the Current Eigenvector float[][] R = OuterProduct(b,c,N,N); //Update A by Subtracting Out Reconstructed Matrix (Eigenvector x Coefficients) A = MatrixSubtract(A,R,N,N); //Check Eigenvalue Size -- Terminate if Below Threshold (Just Modelling Noise) lambda = Av/b[0]; if (lambda < 0.1 && truncate) e = N+1; } //Set Transform Matrix to Eigenvalue Matrix (Ditto for Inverse-Transform), Update Image T = V; T_ = Transpose(T,N,N);}//-------------------------------------------------------------------------------------// Random Orthonormal Basis -----------------------------------------------------------//-------------------------------------------------------------------------------------void GetRandomOrthogonalBasis(){ float[][] V = Zeros(N,N); float[] b = new float[N]; //Choose Random Vector, Orthogonalize, Iterate for (int e = 0; e < N; e++){ //Randomly Initialize Vector b = VectorRandom(N); b = VectorNormalize(b,N); if (e > 0){ //Project Random Vector onto Basis-So-Far float[] c = MatrixVectorMultiply(V,b,e,N); //Reconstruct Random Vector From Basis-So-Far float[] d = MatrixVectorMultiply(Transpose(V,e,N),c,N,e); //Subtract Reconstruction From Random Vector: Residual is Orthogonal to Subspace Spanned by Basis b = VectorSubtract(b,d,N); b = VectorNormalize(b,N); } //Copy New Vector into Matrix arraycopy(b,V[e]); } //Set Transform Matrix to Eigenvalue Matrix (Ditto for Inverse-Transform), Update Image T = V; T_ = Transpose(T,N,N);}//-------------------------------------------------------------------------------------// Matrix Functions (Basic) -----------------------------------------------------------//-------------------------------------------------------------------------------------//Empty Matrixfloat[][] Zeros(int m, int n){ float[][] Z = new float[m][n]; for (int i = 0; i < m; i++){ for (int j=0; j < n; j++){ Z[i][j] = 0.0;}} return Z;}//Identity Matrixfloat[][] Identity(int n){ float[][] ID = Zeros(n,n); for (int i=0; i < n; i++){ ID[i][i] = 1.0;} return ID;}//Transpose of a Matrixfloat[][] Transpose(float [][] M, int m, int n){ float[][] MT = new float[n][m]; for (int i=0; i < m; i++){ for (int j=0; j < n; j++){ MT[j][i] = M[i][j];}} return MT;}//Muliply two (nx1) vectors v1 and v2float InnerProduct(float[] v1, float[] v2, int n){ float sum = 0.0; for (int i=0; i < n; i++){ sum += v1[i] * v2[i];} return sum;}//Multiply mx1, 1xn vectors to get mxn matrixfloat[][] OuterProduct(float[] v1, float[] v2, int m, int n){ float[][] C = new float[m][n]; for (int i=0; i < m; i++){ for (int j=0; j < n; j++){ C[i][j] = v1[i] * v2[j];}} return C;}//Multiply (mxn) matrix M with (nx1) vector vfloat[] MatrixVectorMultiply(float[][] M, float[] v, int m, int n){ float[] x = new float[m]; for (int i=0; i < m; i++){ x[i] = InnerProduct(M[i],v,n);} return x;}float[][] MatrixMultiply(float[][] A, float[][] B, int m, int n){ float[][] B_ = Transpose(B,n,m); float[][] C = new float[m][m]; for (int i=0; i < m; i++){ for (int j=0; j < m; j++){ C[i][j] = InnerProduct(A[i],B_[j],n);}} return C;}//Multiply a vector by a scalarfloat[] VectorScalarMultiply(float[] v, float s, int n){ float[] result = new float[n]; for (int i=0; i < n; i++){ result[i] = s * v[i];} return result;}//Subtract matrix B from matrix Afloat[][] MatrixSubtract(float [][] A, float [][] B, int m, int n){ float[][] C = new float[m][n]; for (int i=0; i < m; i++){ C[i] = VectorSubtract(A[i],B[i],n);} return C;}float[] VectorSubtract(float[] a, float[] b, int n){ float[] c = new float[n]; for (int j=0; j < n; j++){ c[j] = a[j] - b[j];} return c;}//Normalize a vector to unit lengthfloat[] VectorNormalize(float[] v, int n){ float L = InnerProduct(v,v,n); return VectorScalarMultiply(v, 1.0/(sqrt(L)), n);}float[] VectorRandom(int n){ float[] b = new float[n]; for (int i=0; i < n; i++){ b[i] = random(-1,1);} return b;}//-------------------------------------------------------------------------------------// Coefficient Matrix Editing --------------------------------------------------------//-------------------------------------------------------------------------------------void ZeroColumns(){ float k3 = map(sliderVal,slidersLeft,slidersRight,0,256); for (int y = 0; y < H; y++) {arraycopy(X[y],X_[y]);} for (int k = (int)k3; k < N; k++) {Zero(X_,H,N,-1,k);} }//Zero Out Specified Row or Column of a Matrixvoid Zero(float [][] M, int m, int n, int r, int c){ for (int i=0; i < m; i++) for (int j=0; j < n; j++) if (i == r || j == c) M[i][j] = 0.0;}//-------------------------------------------------------------------------------------// Matrix-Image Conversion ------------------------------------------------------------//-------------------------------------------------------------------------------------//Convert a matrix to an image which can be displayed or savedPImage MatrixToImage(float[][] M, int m, int n, float lo, float hi){ PImage img = createImage(n,m,RGB); img.loadPixels(); for (int x=0; x < n; x++) for (int y=0; y < m; y++) img.pixels[y*n+x] = color(map(M[y][x],lo,hi,0,255)); img.updatePixels(); return img;}//Convert image to a matrixfloat[][] MatrixOfImage(PImage img){ float[][] I = new float[img.height][img.width]; for (int y=0; y < H; y++) for (int x = 0; x < N; x++) I[y][x] = map(brightness(img.pixels[y*N+x]),0,255,-1,1); return I;}//-------------------------------------------------------------------------------------// Image Handling ---------------------------------------------------------------------//-------------------------------------------------------------------------------------float[] s = {1, 0.25, sqrt(2.0/N), 0.001, 0.25}; //Basis Image Contrast Scaling Factorsvoid initTransform(){ if (mode == 1) GetEigenvectors(true); //PCA if (mode == 2) GetCosineMatrix(); //Fourier if (mode == 0) GetIdentity(); //Identity if (mode == 4) GetRandomOrthogonalBasis(); //Random Orthogonal Basis if (mode == 3) GetHaarMatrix(); //Haar Wavelet Basis Transform(); ZeroColumns(); InverseTransform(); imgT = MatrixToImage(T,N,N,-s[mode],s[mode]); imgI_ = MatrixToImage(I_,H,N,-1,1); imgX = MatrixToImage(X,H,N,-1,1);}void LoadImage(int i){ //Load Image loaded = i; imgI = loadImage((i+1) + ".jpg"); imgI.filter(GRAY); imgI.loadPixels(); //Convert Current Image to Matrix I = MatrixOfImage(imgI);}void updateReconstructedImage(){ //Only Recompute Inverse Transform if New Columns Have Been Zeroed Out of the Coefficient Matrix if (valid < 3){ ZeroColumns(); InverseTransform(); imgI_ = MatrixToImage(I_,H,N,-1,1); imgX = MatrixToImage(X_,H,N,-1,1); }}//-------------------------------------------------------------------------------------// User Interaction -------------------------------------------------------------------//-------------------------------------------------------------------------------------int buttonsY = 240 , buttonsX[] = new int[7]; //UI Element Positionsint tabsY = N + 38, tabsX[] = new int[5];float slidersLeft = N + 35;float slidersRight= N + slidersLeft + 1;float sliderVal = slidersRight;int slidersY = 225;PFont font10, font12;int loaded = 0, mode = 1, valid = 0;boolean sliderActive = false;String[] tabText = {"Identity", "PCA", "Fourier", "Haar", "Random"};//Initializationvoid setup(){ noLoop(); size(5+(N+30)*3-25,N+30 + 30); frameRate(15); font12 = loadFont("Helvetica-Bold-12.vlw"); font10 = loadFont("Helvetica-Bold-10.vlw"); buttonsX[0] = 42; //Image Number Buttons for (int i=1; i < 7; i++) {buttonsX[i] = buttonsX[i-1] + 30;} tabsX[0] = 30 + (N+30)*2; //Basis Tabs for (int i=1; i < 5; i++) {tabsX[i] = tabsX[i-1] + 51;} LoadImage(0); //Load Image initTransform(); //Initialize and Perform PCA valid = 0; loop();}//Draw function only really executes when a change has invalidated the current image reconstructionvoid draw(){ if (valid < 3){ background(221,221,204); stroke(0); fill(255); //Draw Images with Frames rect(3 ,23,N+3,H+3); image(imgI_, 5 , 25); rect(3 + N+30 ,23,N+3,H+3); image(imgX , 5 + N+30 , 25); rect(3 + (N+30)*2,23,N+3,N+3); image(imgT , 5 + (N+30)*2, 25); //Draw Text textFont(font12); textAlign(CENTER); fill(0); text("=",276,110+20); text("x",563,110+20); text("Reconstructed Image", 133 , 17); text("Transformed Image" , 133+ N+30 , 17); text("Orthogonal Basis" , 133+(N+30)*2 , 17); //Draw Image Number Buttons for (int i=0; i < 7; i++){ fill((loaded==i) ? 200:255); rect(buttonsX[i]-9,buttonsY-12,16,16); fill(0); text(i+1,buttonsX[i],buttonsY); } //Draw Basis Buttons textFont(font10); textAlign(CENTER); for (int i=0; i < 5; i++){ fill((mode==i) ? 200:255); rect(tabsX[i]-25,tabsY-12,51,16); fill(0); text(tabText[i],tabsX[i],tabsY); } //Draw Sliders fill(255); rect(slidersLeft,slidersY+3,slidersRight-slidersLeft,4); fill(128); rect(sliderVal+5,slidersY+3,slidersRight-(sliderVal+5),4); fill(0); rect(sliderVal-5,slidersY,10,10); int nrComponents = (int) map(sliderVal,slidersLeft,slidersRight,0,256); text(nrComponents + " Basis Vectors Used", 133+ N+30 , N + 5); valid = max(0,valid+1); //Formerly Boolean, but Sometimes Didn't Draw at Start :( }}//-------------------------------------------------------------------------------------// Mouse and Keyboard Interaction -----------------------------------------------------//-------------------------------------------------------------------------------------char[] numKeys = {'1','2','3','4','5','6','7'};char[] tabKeys = {'i','p','f','h','r'};void mouseReleased() {sliderActive = false;}void keyPressed(){ if (key == ',') {sliderVal = max(sliderVal-1,slidersLeft+2); valid--;} if (key == '.') {sliderVal = min(sliderVal+1,slidersRight); valid--;} for (int i=0; i < 7; i++) {if (key == numKeys[i]) {LoadImage(i); initTransform(); valid--;}} for (int i=0; i < 5; i++) {if (key == tabKeys[i]) {mode = i; initTransform(); valid--;}} updateReconstructedImage();}void mouseDragged(){ if (sliderActive){sliderVal = max(min(mouseX,slidersRight),slidersLeft+2); valid--;} updateReconstructedImage();}void mousePressed(){ for (int i=0; i < 7; i++){ //Image Number Buttons if (sq(mouseX-buttonsX[i]) + sq(mouseY-buttonsY) < sq(10)){ LoadImage(i); initTransform(); valid--;}} for (int i=0; i < 5; i++){ //Basis Tabs if (sq(mouseX-tabsX[i]) < sq(25) && sq(mouseY-tabsY) < sq(10)){ mode = i; initTransform(); valid--;}} if (mouseX > (slidersLeft-10) && mouseX <= (slidersRight+10)){ sliderActive = true;} //Slider mouseDragged();}`