1; # this is just to tell the Octave Interpreter that this is not a function file


##################################################################################
# FOURIER JPEG (see also http://amath.colorado.edu/courses/5720/2000Spr/Labs/FA/)
##################################################################################

function [x,y] = gen_1D_basis(N,n)
% [x,y] = gen_1D_basis(N,n)
%
% y is the nth Fourier basis vector for the space of complex valued
% functions on N points. x is the domain, useful for plotting
% with the command `plot(x,y)`
%
  omega = exp(-2*pi*I / N);
  x = [0:N-1];
  y = arrayfun(@(x) omega^(x*n),x);
endfunction

function X = gen_2D_basis(N,n,m)
% X = gen_2D_basis(N,n,m)
%
% X is the (n,m)th Fourier basis vector (an N x N matrix) for the space of complex valued functions
% on the N x N square; i.e. (omega_N^{kx+my})_{(x,y)\in N^2} where omega = e^{-2*pi*I / N} is the nth root of unity.
%
  X = zeros(N);
  omega = exp(-2*pi*I / N);
  for x = [0:N-1]
    for y = [0:N-1]
      X(x+1,y+1) = omega^(x*n+y*m);
    endfor
  endfor
endfunction

function [C,U] = thresholdp_fourier(X,p)
% [C,U] = thresholdp_fourier(X,p)
% 
% C is the 2D-Fourier transform of X with its smallest p*100% entries
% zeroed.  U is the Fourier inverse of X (i.e. the de-compressed image)
%
  C = fft2(X);
  aC = abs(C);
  aC1 = aC / (max(max(C)));
  N = prod(size(aC1));
  oC = sort(reshape(aC1,N,1));
  epsilon = oC(ceil(p*N));
  C = C .* (aC1 > epsilon);
  U = ifft2(C);
endfunction

function [C,U] = low_pass(X,w)
% [C,U] = low_pass(X,w)
%
% C is the 2D-Fourier transform of X with its central square of length
% N-2*w zeroed.  U is the Fourier inverse of X (i.e. the de-compressed
% image)
%
  C = fft2(X);
  N = size(X)(1);
  C(w+1:N-w,w+1:N-w) = zeros(N-2*w);
  U = ifft2(C);
endfunction

function a = pcz(C)
% a = pcz(C)
%
% a is the percentage of zeros in the matrix C
%
  n = prod(size(C));
  a = (n - nnz(C)) / n*100;
endfunction

########################################################
# FOURIER EXAMPLES
#######################################################

# [x,y] = gen_1D_basis(512,3);
# plot(x,real(y),x,imag(y));

# X = gen_2D_basis(16,0,1); 
# imshow((real(X)+1)/2);

# im = imread("spiral.jpg");
# hist(reshape(abs(im),prod(size(im)),1), 100); # the distribution of the intensity of coefficients in the standard basis
# hist(reshape(abs(fft2(im)),prod(size(im)),1), 100); # the distribution of the absolute value of the coefficients in the 2D Fourier basis
# [C,U] = thresholdp_fourier(im, 0.95); # compress image at 95% threshold
# imshow(uint8(real(U))-im); # the error -- sharp edges 
# imshow(uint8(real(U))); # the compressed image
# imshow(~~C); # the non-zeroed frequency components in the 2D Fourier coefficient matrix.

# [C,U] = low_pass(im, 20); # remove components with frequencies in both directions greater than 20/512 per pixel 
# pcz(C) # percentage of zeros in C







##################################################################################
# WAVELET JPEG (see also
# http://amath.colorado.edu/courses/5720/2000Spr/Labs/Haar/haar.html)

# NB: You will need to generate the OCT files fhaar.oct and ifhaar.oct.
# If you have installed octave from source, or have the octave-devel
# package, available in some Linux distributions,
# it should be sufficient to type `make` at the command line
##################################################################################

if ((exist("ifhaar") == 3) && (exist("fhaar") == 3))

  function [C,U] = thresholdp_wavelet(X,p)
% [C,U] = thresholdp_wavelet(X,p)
% 
% C is the 2D-Haar transform of X with its smallest p*100% entries
% zeroed.  U is the Haar inverse of X (i.e. the de-compressed image)
%
    n = log2(size(X)(1));
    m = log2(size(X)(2));
    C = fhaar(X, min(n,m));
    if p > 0
      aC = abs(C);
      aC1 = aC / (max(max(aC)));
      N = prod(size(aC1));
      oC = sort(reshape(aC1,N,1));
      epsilon = oC(ceil(p*N));
      C = C .* (aC1 > epsilon);
    endif
    U = ifhaar(C, min(n,m));
  endfunction

else

  function [C,U] = thresholdp_wavelet(X,p)
    error("fhaar.oct and ifhaar.oct not compiled -- reload this script \
when they have been compiled");
  endfunction

endif


########################################################
# WAVELET EXAMPLES
#######################################################

# im = imread("spiral.jpg");
# [C,U] = thresholdp_wavelet(double(im),0);
# imshow(log(abs(C)) / (max(max(log(abs(C)))))); # the coefficient matrix in the 2D Haar basis.
# [C,U] = thresholdp_wavelet(double(im),0.95); # compress at 95% threshold
# imshow(uint8(U));