Ideal point spread function

Contents

ideal1

% calculate point spread function
N = 256;
k = -N/2: N/2-1;
dx = 1/16;
x = k*dx;
D = N*dx;
circ = fftprep(imread('circ16.tif'),N);
subplot(1,2,1);
rgb = cat(3,circ,circ,circ);
imagesc(x,x,rgb);
axis square
%imwrite(circ,'circ.jpg');
[p1, vref] = psf(circ);
xf = k/D;
view = logim(p1,3);
subplot(1,2,2);
rgb = cat(3,view,view,view);
imagesc(xf,xf,rgb);
axis square
%imwrite(view,'circ_psf.jpg');

% compare FFT to analytic curve
close all;
d = 32;
a = D/d;
s = 5;
idx = find(abs(xf)<s);
u = linspace(-s,s,1+50*s);
v = somb(2*u).^2;
subplot(2,1,1);
plot(u,v,'k',xf(idx),p1(129,idx),'o');
xlabel('radial distance');
ylabel('value');
subplot(2,1,2);
semilogy(u,v,'k',xf(idx),p1(129,idx),'o');
xlabel('radial distance');
ylabel('value');

ideal 2

% calculate point spread function
N = 256;
circ = fftprep(imread('circ16.tif'),N);
[p1, vref] = psf(circ);
view = logim(p1,8);
subplot(1,2,1);
imshow(view);
%imwrite(view,'circ_psf5.jpg');

% compare FFT to analytic curve (unit diameter)
k = -N/2: N/2-1;
dx = 1/32;
D = N*dx;
x = ones(N,1)*k/D;
y = k'*ones(1,N)/D;
r = sqrt(x.*x+y.*y);
v = somb(r).^2;
z = logim(v,8);
subplot(1,2,2);
imshow(z);
%imwrite(z,'circ_exact.jpg');

shift

N=256;
k = -N/2:N/2-1;
dx = 1/16;
x = k*dx;

x = ones(N,1)*x;
phase = exp(-i*2*pi*x);
p2 = psf(circ.*phase);
rgb = cat(3,logim(p2,3),zeros(N,N),logim(p1,3));
close all;
imshow(rgb);
%imwrite(rgb,'shift.jpg');

coma

N = 256;
k = -N/2:N/2-1;
dx = 1/16;
x = k*dx;
x = ones(N,1)*x;
y = x';
rsq = x.*x + y.*y;

W = 5; % wave of coma

phase = exp(-i*2*pi*W*rsq.*x);
coma = psf(circ.*phase);
subplot(1,2,1);
imshow(logim(coma,3));
%imwrite(logim(coma,3),'coma.jpg');

phase = exp(-i*2*pi*W*(rsq-1).*x);
coma2 = psf(circ.*phase);
subplot(1,2,2);
imshow(logim(coma2,3));
%imwrite(logim(coma2,3),'coma2.jpg');

coma video

close all;
[coma, vref] = psf(circ);
%mov = avifile('coma1.avi','fps',5);
mov = avifile('coma2.avi','fps',5,'compression','Cinepak');
out = logim(coma,3);
rgb = cat(3,out,out,out);
f = im2frame(rgb);
mov = addframe(mov,f);
waves = 3;

for w=1:60
   A = waves*w/60;
   phase = exp(-i*2*pi*A*(rsq-1).*x);
   coma = psf(circ.*phase,vref);
   out = logim(coma,3);
   rgb = cat(3,out,out,out);
   f = im2frame(rgb);
   mov = addframe(mov,f);
end

mov = close(mov);

Published with MATLAB® 7.6