I'm trying to produce some computer generated holograms by using MATLAB. I used equally spaced mesh grid to initialize the spatial grid, and I got the following image
This pattern is sort of what I need except the center region. The fringe should be sharp but blurred. I think it might be the problem of the mesh grid. I tried generate a grid in polar coordinates and the map it into Cartesian coordinates by using MATLAB's pol2cart
function. Unfortunately, it doesn't work as well. One may suggest that using fine grids. It doesn't work too. I think if I can generate a spiral mesh grid, perhaps the problem is solvable. In addition, the number of the spiral arms could, in general, be arbitrary, could anyone give me a hint on this?
I've attached the code (My final projects are not exactly the same, but it has a similar problem).
clc; clear all; close all;
%% initialization
tic
lambda = 1.55e-6;
k0 = 2*pi/lambda;
c0 = 3e8;
eta0 = 377;
scale = 0.25e-6;
NELEMENTS = 1600;
GoldenRatio = (1+sqrt(5))/2;
g = 2*pi*(1-1/GoldenRatio);
pntsrc = zeros(NELEMENTS, 3);
phisrc = zeros(NELEMENTS, 1);
for idxe = 1:NELEMENTS
pntsrc(idxe, :) = scale*sqrt(idxe)*[cos(idxe*g), sin(idxe*g), 0];
phisrc(idxe) = angle(-sin(idxe*g)+1i*cos(idxe*g));
end
phisrc = 3*phisrc/2; % 3 arms (topological charge ell=3)
%% post processing
sigma = 1;
polfilter = [0, 0, 1i*sigma; 0, 0, -1; -1i*sigma, 1, 0]; % cp filter
xboundl = -100e-6; xboundu = 100e-6;
yboundl = -100e-6; yboundu = 100e-6;
xf = linspace(xboundl, xboundu, 100);
yf = linspace(yboundl, yboundu, 100);
zf = -400e-6;
[pntobsx, pntobsy] = meshgrid(xf, yf);
% how to generate a right mesh grid such that we can generate a decent result?
pntobs = [pntobsx(:), pntobsy(:), zf*ones(size(pntobsx(:)))];
% arbitrary mesh may result in "wrong" results
NPNTOBS = size(pntobs, 1);
nxp = length(xf);
nyp = length(yf);
%% observation
Eobs = zeros(NPNTOBS, 3);
matlabpool open local 12
parfor nobs = 1:NPNTOBS
rp = pntobs(nobs, :);
Erad = [0; 0; 0];
for idx = 1:NELEMENTS
rs = pntsrc(idx, :);
p = exp(sigma*1i*2*phisrc(idx))*[1 -sigma*1i 0]/2; % simplified here
u = rp - rs;
r = sqrt(u(1)^2+u(2)^2+u(3)^2); %norm(u);
u = u/r; % unit vector
ut = [u(2)*p(3)-u(3)*p(2),...
u(3)*p(1)-u(1)*p(3), ...
u(1)*p(2)-u(2)*p(1)]; % cross product: u cross p
Erad = Erad + ... % u cross p cross u, do not use the built-in func
c0*k0^2/4/pi*exp(1i*k0*r)/r*eta0*...
[ut(2)*u(3)-ut(3)*u(2);...
ut(3)*u(1)-ut(1)*u(3); ...
ut(1)*u(2)-ut(2)*u(1)];
end
Eobs(nobs, :) = Erad; % filter neglected here
end
matlabpool close
Eobs = Eobs/max(max(sum(abs(Eobs), 2))); % normailized
%% source, gaussian beam
E0 = 1;
w0 = 80e-6;
theta = 0; % may be titled
RotateX = [1, 0, 0; ...
0, cosd(theta), -sind(theta); ...
0, sind(theta), cosd(theta)];
Esrc = zeros(NPNTOBS, 3);
for nobs = 1:NPNTOBS
rp = RotateX*[pntobs(nobs, 1:2).'; 0];
z = rp(3);
r = sqrt(sum(abs(rp(1:2)).^2));
zR = pi*w0^2/lambda;
wz = w0*sqrt(1+z^2/zR^2);
Rz = z^2+zR^2;
zetaz = atan(z/zR);
gaussian = E0*w0/wz*exp(-r^2/wz^2-1i*k0*z-1i*k0*0*r^2/Rz/2+1i*zetaz);% ...
Esrc(nobs, :) = (polfilter*gaussian*[1; -1i; 0]).'/sqrt(2)/2;
end
Esrc = [Esrc(:, 2), Esrc(:, 3), Esrc(:, 1)];
Esrc = Esrc/max(max(sum(abs(Esrc), 2))); % normailized
toc
%% visualization
fringe = Eobs + Esrc; % I'll have a different formula in my code
normEsrc = reshape(sum(abs(Esrc).^2, 2), [nyp nxp]);
normEobs = reshape(sum(abs(Eobs).^2, 2), [nyp nxp]);
normFringe = reshape(sum(abs(fringe).^2, 2), [nyp nxp]);
close all;
xf0 = linspace(xboundl, xboundu, 500);
yf0 = linspace(yboundl, yboundu, 500);
[xfi, yfi] = meshgrid(xf0, yf0);
data = interp2(xf, yf, normFringe, xfi, yfi);
figure; surf(xfi, yfi, data,'edgecolor','none');
% tri = delaunay(xfi, yfi); trisurf(tri, xfi, yfi, data, 'edgecolor','none');
xlim([xboundl, xboundu])
ylim([yboundl, yboundu])
% colorbar
view(0,90)
colormap(hot)
axis equal
axis off
title('fringe thereo. ', ...
'fontsize', 18)
I didn't read your code because it is too long to do such a simple thing. I wrote mine and here is the result:
the code is
%spiral.m
function val = spiral(x,y)
r = sqrt( x*x + y*y);
a = atan2(y,x)*2+r;
x = r*cos(a);
y = r*sin(a);
val = exp(-x*x*y*y);
val = 1/(1+exp(-1000*(val)));
endfunction
%show.m
n=300;
l = 7;
A = zeros(n);
for i=1:n
for j=1:n
A(i,j) = spiral( 2*(i/n-0.5)*l,2*(j/n-0.5)*l);
end
end
imshow(A) %don't know if imshow is in matlab. I used octave.
the key for the sharpnes is line
val = 1/(1+exp(-1000*(val)));
It is logistic function. The number 1000 defines how sharp your image will be. So lower it for more blurry image or higher it for sharper.
I hope this answers your question ;)
Edit: It is real fun to play with. Here is another spiral:
function val = spiral(x,y)
s= 0.5;
r = sqrt( x*x + y*y);
a = atan2(y,x)*2+r*r*r;
x = r*cos(a);
y = r*sin(a);
val = 0;
if (abs(x)<s )
val = s-abs(x);
endif
if(abs(y)<s)
val =max(s-abs(y),val);
endif
%val = 1/(1+exp(-1*(val)));
endfunction
Edit2: Fun, fun, fun! Here the arms do not get thinner.
function val = spiral(x,y)
s= 0.1;
r = sqrt( x*x + y*y);
a = atan2(y,x)*2+r*r; % h
x = r*cos(a);
y = r*sin(a);
val = 0;
s = s*exp(r);
if (abs(x)<s )
val = s-abs(x);
endif
if(abs(y)<s)
val =max(s-abs(y),val);
endif
val = val/s;
val = 1/(1+exp(-10*(val)));
endfunction
Damn your question I really need to study for my exam, arghhh!
Edit3:
I vectorised the code and it runs much faster.
%spiral.m
function val = spiral(x,y)
s= 2;
r = sqrt( x.*x + y.*y);
a = atan2(y,x)*8+exp(r);
x = r.*cos(a);
y = r.*sin(a);
val = 0;
s = s.*exp(-0.1*r);
val = r;
val = (abs(x)<s ).*(s-abs(x));
val = val./s;
% val = 1./(1.+exp(-1*(val)));
endfunction
%show.m
n=1000;
l = 3;
A = zeros(n);
[X,Y] = meshgrid(-l:2*l/n:l);
A = spiral(X,Y);
imshow(A)
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With