Orthonormal Basis Transforms
Grant Schindler, 2007

This browser does not have a Java Plug-in.
Get the latest Java Plug-in here.


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, 2007

int N = 256, H = 192;            //Image Dimensions
float[][] T = new float[N][N];   //Transform Matrix
float[][] T_;                    //Inverse Transform Matrix
float[][] X  = new float[H][N];  //Transformed Image
float[][] X_ = new float[H][N];  //Transformed Image with Zeroed Elements
float[][] I;                     //Original Image
float[][] I_ = new float[H][N];  //Reconstructed Image
PImage imgI,imgI_,imgT,imgX;     //Displayable Image for Each Matrix

//-------------------------------------------------------------------------------------
// Orthonormal Basis Transform: 2 Simple Functions ------------------------------------
//-------------------------------------------------------------------------------------

//Transform - Use T to Compute Coefficients for Each Image Row
void Transform(){
  for (int y=0; y < H; y++)
    X[y] = MatrixVectorMultiply(T,I[y],N,N);
}

//InverseTransform - Reconstruct Original Image by Inverting Transform
void 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 Waves
void 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 Wavelets
void 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 Matrix
void 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 Matrix
float[][] 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 Matrix
float[][] 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 Matrix
float[][] 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 v2
float 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 matrix
float[][] 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 v
float[] 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 scalar
float[] 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 A
float[][] 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 length
float[] 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 Matrix
void 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 saved
PImage 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 matrix
float[][] 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 Factors

void 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 Positions
int 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"};

//Initialization
void 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 reconstruction
void 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();
}