Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

spiral meshgrid in matlab

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

enter image description here

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)
like image 837
Tony Dong Avatar asked Apr 02 '13 04:04

Tony Dong


1 Answers

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:

enter image description here

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:

spiral2

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.

spiral3

  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)
like image 56
tom Avatar answered Sep 22 '22 02:09

tom