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));