pobieranie: Matlab Tips and Tricks (ang)
Pobierz dokument pdfPrzeglądaj wersję html pliku:
Matlab Tips and Tricks
Gabriel Peyr´ e peyre@cmapx.polytechnique.fr August 10, 2004
First keep in mind that this is not a Matlab tutorial. This is just a list of tricks I have found useful while writing my toolboxes available on the Matlab Central repository http://www.mathworks.com/matlabcentral/ You can e-mail me if you have corrections about these pieces of code, or if you would like to add your own tips to those described in this document.
Contents
1 General Programming Tips 2 General Mathematical Tips 3 Advanced Mathematical Tips 4 Signal and Image Processing Tips 5 Graph Theory Tips 1 4 4 6 10
1
General Programming Tips
• Suppress entries in a vector.
x( 3:5 ) = [];
• Reverse a vector.
x = x(end:-1:1);
• Compute the running time of a function call.
1
tic; fft(rand(500)); disp( [’it takes ’ num2str(toc) ’s.’]
);
• Make a array full of a
% guess which one is the fastest ? tic; NaN*ones(2000,2000); toc; tic; repmat(NaN,2000,2000); toc;
• Turn an nD array into a vector.
x = x(:);
• Compute the maximum value of an nD array.
m = max(x(:));
• Access a matrix from a list of entries. Here, we have I = [I1; I2] and y(i) = M( I1(i), I2(i) )
J = sub2ind(size(M), I(1,:),I(2,:) y = M(J); );
• Create a function that take optional arguments in a struct.
function y = f(x,options) % parse the struct if nargin<2 options.null = 0; % force creation of options end if isfield(options, ’a’) options.a = 1; % default value end a = options.a; if isfield(options, ’b’) options.b = 1; % default value end b = options.b; % Here the body of the function ...
• Create a graphical waitbar.
n = 100; h = waitbar(0,’Waiting ...’); for i=1:n waitbar(i/n); % here perform some stuff end close(h);
• How to duplicate a character n times.
2
str = char( zeros(n,1)+’*’ );
• Output a string without carriage return.
fprintf(’Some Text’);
• Assign value v in a nD array at a position ind (lenth-n vector).
ind = num2cell(ind); x( ind{:} ) = v;
• Save the current figure as an image in e.g. EPS file format.
saveas(gcf, str, ’png’);
• Remove the ticks from a drawing.
set(gca, ’XTick’, []); set(gca, ’YTick’, []);
• Saving and loading an image.
saveas(gcf, ’my image’, ’png’); % save M = double( imread( ’my image.png’ ) ); % load
• Saving and loading a matrix M in a binary file.
[n,p] = size(M); % saving str = ’my file’; % name of the file fid = fopen(str,’wb’); if fid<0 error([’error writing to file ’, str]); end fwrite(fid,M,’double’); fclose(fid); % loading fid = fopen(str,’rb’); if fid<0 error([’error reading file ’,str]); end [M, cnt] = fread(fid,[n,p],’double’); fclose(fid); if cnt =n*p error([’Error reading file ’, str]); end
• Find the angle that makes a 2D vector x with the vector [1,0]
% just the angle theta = atan2(x(2),x(1)); % if you want to compute the full polar decomposition [theta,r] = cart2pol(x);
3
2
General Mathematical Tips
• Rescale the entries of a vector x so that it spans [0, 1]
m = min(x(:)); M = max(x(:)); x = (b-a) * (x-m)/(M-m) + a;
• Generate n points evenly sampled.
x = 0:1/(n-1):1; % faster than linspace
• Compute the L2 squared norm of a vector or matrix x.
m = sum(x(:).ˆ2);
• Subsample a vector x or an image M by a factor 2.
x = x(1:2:end); % useful for wavelet transform M = M(1:2:end,1:2:end);
• Compute centered finite differences.
D1 = [x(2:end),x(end)]; D2 = [x(1),x(1:end-1)]; y = (D1-D2)/2;
• Compute the prime number just before n
n = 150; P = primes(n); n = P(end);
• Compute J, the reverse of a permutation I, i.e. an array which contains the number 1:n in arbitrary order.
J(I) = 1:length(I);
• Shuffle an array x.
y = x( randperm(length(x)) );
3
Advanced Mathematical Tips
• Generate n points x sampled uniformly at random on a sphere.
4
% tensor product gaussian is isotropic x = randn(3,n); d = sqrt( x(1,:).ˆ2+x(2,:).ˆ2+x(2,:).ˆ2 ); x(1,:) = x(1,:)./d; x(2,:) = x(2,:)./d; x(3,:) = x(3,:)./d;
• Construct a polygon x whose ith sidelength is s(i). Here x(i) is the complex affix of the ith vertex.
theta = [0;cumsum(s)]; theta = theta/theta(end); theta = theta(1:(end-1)); x = exp(2i*pi*theta); L = abs(x(1)-x(2)); x = x*s(1)/L; % rescale the result
• Compute y, the inverse of an integer x modulo a prime p.
% use Bezout thm [u,y,d] = gcd(x,p); y = mod(y,p);
• Compute the curvilinear abscise s of a curve c. Here, c(:,i) is the ith point of the curve.
D = c(:,2:end)-c(:,1:(end-1)); s = zeros(size(c,2),1); s(2:end) = sqrt( D(1,:).ˆ2 + D(2,:).ˆ2 ); s = cumsum(s);
• Compute the 3D rotation matrix M around an axis v
% v S M taken from the OpenGL red book = v/norm(v,’fro’); = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]; = v*transp(v) + cos(alpha)*(eye(3) - v*transp(v)) + sin(alpha)*S;
• Compute a VanderMonde matrix M i.e. M(i,j)=x(i)ˆj for j=0:d.
n = length(x); % first method [J,I] = meshgrid(0:d,1:n); A = x(I).ˆJ; % second method, less elegant but faster A = ones(n); for j = 2:n A(:,j) = x.*A(:,j-1); end
• Threshold (i.e. set to 0) the entries below T.
5
% x % I
first solution = (abs(x)>=T) .* x; second one : nearly 2 times slower = find(abs(x)<T); x(I) = 0;
• Keep only the n biggest coefficients of a signal x (set the others to 0).
[,I] = sort(abs(x(:))); x( I(1:end-n) ) = 0;
• Draw a 3D sphere.
p = 20; % precision t = 0:1/(p-1):1; [th,ph] = meshgrid( t*pi,t*2*pi ); x = cos(th); y = sin(th).*cos(ph); z = sin(th).*sin(ph); surf(x,y,z, z.*0); % some pretty rendering options shading interp; lighting gouraud; camlight infinite; axis square; axis off;
• Project 3D points on a 2D plane (best fit plane). P(:,k) is the kth point.
for i=1:3 % substract mean P(i,:) = P(i,:) - mean(P(i,:)); end C = P*P’; % covariance matrix % project on the two most important eigenvectors [V,D] = eigs(C); Q = V(:,1:2)’*P;
4
Signal and Image Processing Tips
• Compute circular convolution of x and y.
% use the Fourier convolution thm z = real( ifft( fft(x).*fft(y) ) );
• Display the result of an FFT with the 0 frequency in the middle.
x = peaks(256); imagesc( real( fftshift( fft2(x) ) ) );
• Resize an image M (new size is (p1,q1)).
[p,q] = size(M); % the original image [X,Y] = meshgrid( (0:p-1)/(p-1), (0:q-1)/(q-1) ); % new sampling location [XI,YI] = meshgrid( (0:p1-1)/(p1-1) , (0:q1-1)/(q1-1) ); M1 = interp2( X,Y, M, XI,YI ,’cubic’); % the new image
6
• Build a 1D gaussian filter of variance s.
x = -1/2:1/(n-1):1/2; f = exp( -(x.ˆ2)/(2*sˆ2) ); f = f / sum(sum(f));
• Build a 2D gaussian filter of variance s.
x = -1/2:1/(n-1):1/2; [Y,X] = meshgrid(x,x); f = exp( -(X.ˆ2+Y.ˆ2)/(2*sˆ2) ); f = f / sum(f(:));
• Perform a 1D convolution of signal f and filter h with symmetric boundary conditions. The center of the filter is 0 for odd length filter, and 1/2 otherwise
n = length(x); p = length(h); if mod(p,2)==1 d1 = (p-1)/2; d2 = (p-1)/2; else d1 = p/2-1; d2 = p/2; end xx = [ x(d1:-1:1); x; x(end:-1:end-d2+1) ]; y = conv(xx,h); y = y( (2*d1+1):(2*d1+n) );
• Same but for 2D signals
n = length(x); p = length(h); if mod(p,2)==1 d1 = (p-1)/2; d2 = (p-1)/2; else d1 = p/2-1; d2 = p/2; end xx = [ x(d1:-1:1,:); x; x(end:-1:end-d2+1,:) ]; xx = [ xx(:,d1:-1:1), xx, xx(:,end:-1:end-d2+1) ]; y = conv2(xx,h); y = y( (2*d1+1):(2*d1+n), (2*d1+1):(2*d1+n) );
• Extract all 0th level curves from an image M an put these curves into a cell array c list.
c = contourc(M,[0,0]); k = 0; p = 1; while p < size(c, 2) % parse the result lc = c(2,p); % length of the curve cc = c(:,(p+1):(p+lc)); p = p+lc+1; k = k+1; c list{k} = cc; end
• Quick computation of the integral y of an image M along a 2D curve c (the
7
curve is assumed in [0, 1]2 )
cs = c*(n-1) + 1; % scale to [1,n] I = round(cs); J = sub2ind(size(M), I(1,:),I(2,:) ); y = sum(M(J));
• Draw the image of a disk and a square.
n = 100; x = -1:2/(n-1):1; [Y,X] = meshgrid(x,x); c = [0,0]; r = 0.4; % center and radius of the disk D = (X-c(1)).ˆ2 + (Y-c(2)).ˆ2 < rˆ2; imagesc(D); % a disk C = max(abs(X-c(1)),abs(Y-c(2)))<r; imagesc(C); % a square
• Draw a 2D function whose value z is known only at scattered 2D points (x,y).
n = 400; x = rand(n,1); y = rand(n,1); % this is an example of surface z = cos(pi*x) .* cos(pi*y); tri = delaunay(x,y); % build a Delaunay triangulation trisurf(tri,x,y,z);
• Generate a signal whose regularity is C α (Sobolev).
alpha = 2; n = 100; y = randn(n,1); % gaussian noise fy = fft(y); fy = fftshift(fy); % filter the noise with |omega|ˆ-alpha h = (-n/2+1):(n/2); h = (abs(h)+1).ˆ(-alpha-0.5); fy = fy.*h’; fy = fftshift(fy); y = real( ifft(fy) ); y = (y-min(y))/(max(y)-min(y));
• Generate a signal whose regularity is nearly C α−1/2 .
alpha = 3; n = 300; x = rand(n,1); % uniform noise for i=1:alpha % integrate the noise alpha times x = cumsum(x - mean(x)); end
• Compute the PSNR between to signals x and y.
d = mean( mean( (x-y).ˆ2 ) ); m = max( max(x(:)),max(y(:)) ); PSNR = 10*log10( m/d );
8
• Evaluate a cubic spline at value t (can be a vector).
x = abs(t) ; I12 = (x>1)&(x<=2); I01 = (x<=1); y = I01.*( 2/3-x.ˆ2.*(1-x/2) ) + I12.*( 1/6*(2-x).ˆ3 );
• Perform spectral interpolation of a signal x (aka Fourier zero-padding). The original size is n and the final size is p
n f f x = = = = length(x); n0 = (n-1)/2; fft(x); % forward transform p/n*[f(1:n0+1); zeros(p-n,1); f(n0+2:n)]; real( ifft(f) ); % backward transform
• Compute the approximation error err= ||f −fM ||/||f || obtained when keeping the M best coefficients in an orthogonal basis.
% as an example we take the decomposition in the cosine basis M = 500; x = peaks(128); y = dct(x); % a sample function [tmp,I] = sort(abs(y(:))); y(I(1:end-M)) = 0; err = norm(y,’fro’)/norm(x,’fro’); % the relative error xx = idct(y); imagesc(xx); % the reconstructed function
• Perform a JPEG-like transform of an image x (replace dct by idct to compute the inverse transform).
bs = 8; % size of the blocks n = size(x,1); y = zeros(n,n); nb = n/bs; % n must be a multiple of bs for i=1:nb for j=1:nb xsel = ((i-1)*bs+1):(i*bs); ysel = ((j-1)*bs+1):(j*bs); y(xsel,ysel) = dct(x(xsel,ysel)); end end
• Extract interactively a part MM of an image M.
[n,p] = size(M); imagesc(M); axis image; axis off; sp = getrect; sp(1) = max(floor(sp(1)),1); % xmin sp(2) = max(floor(sp(2)),1); % ymin sp(3) = min(ceil(sp(1)+sp(3)),p); % xmax sp(4) = min(ceil(sp(2)+sp(4)),n); % ymax MM = M(sp(2):sp(4), sp(1):sp(3));
9
5
Graph Theory Tips
• Compute the shortest distance between all pair of nodes (D is the weighted adjacency matrix).
% non connected vectices must have Inf value N = length(D); for k=1:N D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1])); end D1 = D;
• Turn a triangulation into an adjacency matrix.
nvert = max(max(face)); nface = length(face); A = zeros(nvert); for i=1:nface A(face(i,1),face(i,2)) = 1; A(face(i,2),face(i,3)) = 1; A(face(i,3),face(i,1)) = 1; % make sure that all edges are symmetric A(face(i,2),face(i,1)) = 1; A(face(i,3),face(i,2)) = 1; A(face(i,1),face(i,3)) = 1; end
10
Matlab Tips and Tricks (ang)
Matlab Tips and Tricks
Gabriel Peyr´ e peyre@cmapx.polytechnique.fr August 10, 2004
First keep in mind that this is not a Matlab tutorial. This is just a list of tricks I have found useful while writing my toolboxes available on the Matlab Central repository http://www.mathworks.com/matlabcentral/ You can e-mail me if you have corrections about these pieces of code, or if you would like to add your own tips to those described in this document.
Contents
1 General Programming Tips 2 General Mathematical Tips 3 Advanced Mathematical Tips 4 Signal and Image Processing Tips 5 Graph Theory Tips 1 4 4 6 10
1
General Programming Tips
• Suppress entries in a vector.
x( 3:5 ) = [];
• Reverse a vector.
x = x(end:-1:1);
• Compute the running time of a function call.
1
tic; fft(rand(500)); disp( [’it takes ’ num2str(toc) ’s.’]
);
• Make a array full of a
% guess which one is the fastest ? tic; NaN*ones(2000,2000); toc; tic; repmat(NaN,2000,2000); toc;
• Turn an nD array into a vector.
x = x(:);
• Compute the maximum value of an nD array.
m = max(x(:));
• Access a matrix from a list of entries. Here, we have I = [I1; I2] and y(i) = M( I1(i), I2(i) )
J = sub2ind(size(M), I(1,:),I(2,:) y = M(J); );
• Create a function that take optional arguments in a struct.
function y = f(x,options) % parse the struct if nargin<2 options.null = 0; % force creation of options end if isfield(options, ’a’) options.a = 1; % default value end a = options.a; if isfield(options, ’b’) options.b = 1; % default value end b = options.b; % Here the body of the function ...
• Create a graphical waitbar.
n = 100; h = waitbar(0,’Waiting ...’); for i=1:n waitbar(i/n); % here perform some stuff end close(h);
• How to duplicate a character n times.
2
str = char( zeros(n,1)+’*’ );
• Output a string without carriage return.
fprintf(’Some Text’);
• Assign value v in a nD array at a position ind (lenth-n vector).
ind = num2cell(ind); x( ind{:} ) = v;
• Save the current figure as an image in e.g. EPS file format.
saveas(gcf, str, ’png’);
• Remove the ticks from a drawing.
set(gca, ’XTick’, []); set(gca, ’YTick’, []);
• Saving and loading an image.
saveas(gcf, ’my image’, ’png’); % save M = double( imread( ’my image.png’ ) ); % load
• Saving and loading a matrix M in a binary file.
[n,p] = size(M); % saving str = ’my file’; % name of the file fid = fopen(str,’wb’); if fid<0 error([’error writing to file ’, str]); end fwrite(fid,M,’double’); fclose(fid); % loading fid = fopen(str,’rb’); if fid<0 error([’error reading file ’,str]); end [M, cnt] = fread(fid,[n,p],’double’); fclose(fid); if cnt =n*p error([’Error reading file ’, str]); end
• Find the angle that makes a 2D vector x with the vector [1,0]
% just the angle theta = atan2(x(2),x(1)); % if you want to compute the full polar decomposition [theta,r] = cart2pol(x);
3
2
General Mathematical Tips
• Rescale the entries of a vector x so that it spans [0, 1]
m = min(x(:)); M = max(x(:)); x = (b-a) * (x-m)/(M-m) + a;
• Generate n points evenly sampled.
x = 0:1/(n-1):1; % faster than linspace
• Compute the L2 squared norm of a vector or matrix x.
m = sum(x(:).ˆ2);
• Subsample a vector x or an image M by a factor 2.
x = x(1:2:end); % useful for wavelet transform M = M(1:2:end,1:2:end);
• Compute centered finite differences.
D1 = [x(2:end),x(end)]; D2 = [x(1),x(1:end-1)]; y = (D1-D2)/2;
• Compute the prime number just before n
n = 150; P = primes(n); n = P(end);
• Compute J, the reverse of a permutation I, i.e. an array which contains the number 1:n in arbitrary order.
J(I) = 1:length(I);
• Shuffle an array x.
y = x( randperm(length(x)) );
3
Advanced Mathematical Tips
• Generate n points x sampled uniformly at random on a sphere.
4
% tensor product gaussian is isotropic x = randn(3,n); d = sqrt( x(1,:).ˆ2+x(2,:).ˆ2+x(2,:).ˆ2 ); x(1,:) = x(1,:)./d; x(2,:) = x(2,:)./d; x(3,:) = x(3,:)./d;
• Construct a polygon x whose ith sidelength is s(i). Here x(i) is the complex affix of the ith vertex.
theta = [0;cumsum(s)]; theta = theta/theta(end); theta = theta(1:(end-1)); x = exp(2i*pi*theta); L = abs(x(1)-x(2)); x = x*s(1)/L; % rescale the result
• Compute y, the inverse of an integer x modulo a prime p.
% use Bezout thm [u,y,d] = gcd(x,p); y = mod(y,p);
• Compute the curvilinear abscise s of a curve c. Here, c(:,i) is the ith point of the curve.
D = c(:,2:end)-c(:,1:(end-1)); s = zeros(size(c,2),1); s(2:end) = sqrt( D(1,:).ˆ2 + D(2,:).ˆ2 ); s = cumsum(s);
• Compute the 3D rotation matrix M around an axis v
% v S M taken from the OpenGL red book = v/norm(v,’fro’); = [0 -v(3) v(2); v(3) 0 -v(1); -v(2) v(1) 0]; = v*transp(v) + cos(alpha)*(eye(3) - v*transp(v)) + sin(alpha)*S;
• Compute a VanderMonde matrix M i.e. M(i,j)=x(i)ˆj for j=0:d.
n = length(x); % first method [J,I] = meshgrid(0:d,1:n); A = x(I).ˆJ; % second method, less elegant but faster A = ones(n); for j = 2:n A(:,j) = x.*A(:,j-1); end
• Threshold (i.e. set to 0) the entries below T.
5
% x % I
first solution = (abs(x)>=T) .* x; second one : nearly 2 times slower = find(abs(x)<T); x(I) = 0;
• Keep only the n biggest coefficients of a signal x (set the others to 0).
[,I] = sort(abs(x(:))); x( I(1:end-n) ) = 0;
• Draw a 3D sphere.
p = 20; % precision t = 0:1/(p-1):1; [th,ph] = meshgrid( t*pi,t*2*pi ); x = cos(th); y = sin(th).*cos(ph); z = sin(th).*sin(ph); surf(x,y,z, z.*0); % some pretty rendering options shading interp; lighting gouraud; camlight infinite; axis square; axis off;
• Project 3D points on a 2D plane (best fit plane). P(:,k) is the kth point.
for i=1:3 % substract mean P(i,:) = P(i,:) - mean(P(i,:)); end C = P*P’; % covariance matrix % project on the two most important eigenvectors [V,D] = eigs(C); Q = V(:,1:2)’*P;
4
Signal and Image Processing Tips
• Compute circular convolution of x and y.
% use the Fourier convolution thm z = real( ifft( fft(x).*fft(y) ) );
• Display the result of an FFT with the 0 frequency in the middle.
x = peaks(256); imagesc( real( fftshift( fft2(x) ) ) );
• Resize an image M (new size is (p1,q1)).
[p,q] = size(M); % the original image [X,Y] = meshgrid( (0:p-1)/(p-1), (0:q-1)/(q-1) ); % new sampling location [XI,YI] = meshgrid( (0:p1-1)/(p1-1) , (0:q1-1)/(q1-1) ); M1 = interp2( X,Y, M, XI,YI ,’cubic’); % the new image
6
• Build a 1D gaussian filter of variance s.
x = -1/2:1/(n-1):1/2; f = exp( -(x.ˆ2)/(2*sˆ2) ); f = f / sum(sum(f));
• Build a 2D gaussian filter of variance s.
x = -1/2:1/(n-1):1/2; [Y,X] = meshgrid(x,x); f = exp( -(X.ˆ2+Y.ˆ2)/(2*sˆ2) ); f = f / sum(f(:));
• Perform a 1D convolution of signal f and filter h with symmetric boundary conditions. The center of the filter is 0 for odd length filter, and 1/2 otherwise
n = length(x); p = length(h); if mod(p,2)==1 d1 = (p-1)/2; d2 = (p-1)/2; else d1 = p/2-1; d2 = p/2; end xx = [ x(d1:-1:1); x; x(end:-1:end-d2+1) ]; y = conv(xx,h); y = y( (2*d1+1):(2*d1+n) );
• Same but for 2D signals
n = length(x); p = length(h); if mod(p,2)==1 d1 = (p-1)/2; d2 = (p-1)/2; else d1 = p/2-1; d2 = p/2; end xx = [ x(d1:-1:1,:); x; x(end:-1:end-d2+1,:) ]; xx = [ xx(:,d1:-1:1), xx, xx(:,end:-1:end-d2+1) ]; y = conv2(xx,h); y = y( (2*d1+1):(2*d1+n), (2*d1+1):(2*d1+n) );
• Extract all 0th level curves from an image M an put these curves into a cell array c list.
c = contourc(M,[0,0]); k = 0; p = 1; while p < size(c, 2) % parse the result lc = c(2,p); % length of the curve cc = c(:,(p+1):(p+lc)); p = p+lc+1; k = k+1; c list{k} = cc; end
• Quick computation of the integral y of an image M along a 2D curve c (the
7
curve is assumed in [0, 1]2 )
cs = c*(n-1) + 1; % scale to [1,n] I = round(cs); J = sub2ind(size(M), I(1,:),I(2,:) ); y = sum(M(J));
• Draw the image of a disk and a square.
n = 100; x = -1:2/(n-1):1; [Y,X] = meshgrid(x,x); c = [0,0]; r = 0.4; % center and radius of the disk D = (X-c(1)).ˆ2 + (Y-c(2)).ˆ2 < rˆ2; imagesc(D); % a disk C = max(abs(X-c(1)),abs(Y-c(2)))<r; imagesc(C); % a square
• Draw a 2D function whose value z is known only at scattered 2D points (x,y).
n = 400; x = rand(n,1); y = rand(n,1); % this is an example of surface z = cos(pi*x) .* cos(pi*y); tri = delaunay(x,y); % build a Delaunay triangulation trisurf(tri,x,y,z);
• Generate a signal whose regularity is C α (Sobolev).
alpha = 2; n = 100; y = randn(n,1); % gaussian noise fy = fft(y); fy = fftshift(fy); % filter the noise with |omega|ˆ-alpha h = (-n/2+1):(n/2); h = (abs(h)+1).ˆ(-alpha-0.5); fy = fy.*h’; fy = fftshift(fy); y = real( ifft(fy) ); y = (y-min(y))/(max(y)-min(y));
• Generate a signal whose regularity is nearly C α−1/2 .
alpha = 3; n = 300; x = rand(n,1); % uniform noise for i=1:alpha % integrate the noise alpha times x = cumsum(x - mean(x)); end
• Compute the PSNR between to signals x and y.
d = mean( mean( (x-y).ˆ2 ) ); m = max( max(x(:)),max(y(:)) ); PSNR = 10*log10( m/d );
8
• Evaluate a cubic spline at value t (can be a vector).
x = abs(t) ; I12 = (x>1)&(x<=2); I01 = (x<=1); y = I01.*( 2/3-x.ˆ2.*(1-x/2) ) + I12.*( 1/6*(2-x).ˆ3 );
• Perform spectral interpolation of a signal x (aka Fourier zero-padding). The original size is n and the final size is p
n f f x = = = = length(x); n0 = (n-1)/2; fft(x); % forward transform p/n*[f(1:n0+1); zeros(p-n,1); f(n0+2:n)]; real( ifft(f) ); % backward transform
• Compute the approximation error err= ||f −fM ||/||f || obtained when keeping the M best coefficients in an orthogonal basis.
% as an example we take the decomposition in the cosine basis M = 500; x = peaks(128); y = dct(x); % a sample function [tmp,I] = sort(abs(y(:))); y(I(1:end-M)) = 0; err = norm(y,’fro’)/norm(x,’fro’); % the relative error xx = idct(y); imagesc(xx); % the reconstructed function
• Perform a JPEG-like transform of an image x (replace dct by idct to compute the inverse transform).
bs = 8; % size of the blocks n = size(x,1); y = zeros(n,n); nb = n/bs; % n must be a multiple of bs for i=1:nb for j=1:nb xsel = ((i-1)*bs+1):(i*bs); ysel = ((j-1)*bs+1):(j*bs); y(xsel,ysel) = dct(x(xsel,ysel)); end end
• Extract interactively a part MM of an image M.
[n,p] = size(M); imagesc(M); axis image; axis off; sp = getrect; sp(1) = max(floor(sp(1)),1); % xmin sp(2) = max(floor(sp(2)),1); % ymin sp(3) = min(ceil(sp(1)+sp(3)),p); % xmax sp(4) = min(ceil(sp(2)+sp(4)),n); % ymax MM = M(sp(2):sp(4), sp(1):sp(3));
9
5
Graph Theory Tips
• Compute the shortest distance between all pair of nodes (D is the weighted adjacency matrix).
% non connected vectices must have Inf value N = length(D); for k=1:N D = min(D,repmat(D(:,k),[1 N])+repmat(D(k,:),[N 1])); end D1 = D;
• Turn a triangulation into an adjacency matrix.
nvert = max(max(face)); nface = length(face); A = zeros(nvert); for i=1:nface A(face(i,1),face(i,2)) = 1; A(face(i,2),face(i,3)) = 1; A(face(i,3),face(i,1)) = 1; % make sure that all edges are symmetric A(face(i,2),face(i,1)) = 1; A(face(i,3),face(i,2)) = 1; A(face(i,1),face(i,3)) = 1; end
10
