螺旋meshgrid在matlab

Question

我想通過使用MATLAB來產生一些計算機生成的全息圖。我用等距網格初始化空間網格，我得到了下面的圖像

這種模式是那種我所需要的，除了中心區域。邊緣應該清晰但模糊。我認爲這可能是網狀網格的問題。我嘗試使用MATLAB的pol2cart函數在極座標中生成網格，並將其映射到笛卡爾座標中。不幸的是，它不工作。有人可能會建議使用細網格。它也行不通。我想如果我能生成一個螺旋網狀網格，問題也許是可以解決的。另外，螺旋形武器的數量一般來說可以是任意的，任何人都可以給我一個暗示嗎？

我附上了代碼（我的最終項目不完全相同，但它有類似的問題）。

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) 

Answer 1

我沒看過你的代碼，因爲做這麼簡單的事情太長了。我寫我的，這裏是結果：

代碼

%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. 

關鍵的sharpnes是線

val = 1/(1+exp(-1000*(val))); 

這是logistic function。數字1000定義了圖像的銳利程度。所以降低它的模糊圖像或更高的清晰度。

我希望這回答了你的問題;）

編輯：這是一起玩真正的樂趣。這裏是另一個螺旋：

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：好玩，有趣，有趣！這裏的手臂不會變薄。

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 

該死你的問題我真的需要學習爲我的考試，arghhh！

EDIT3：

我矢量化的代碼，它運行得更快。

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

Answer 2

對不起，不能發佈的數字。但這可能有幫助。我寫它的幅度空間調製器的實驗...

R=70;   % radius of curvature of fresnel lens (in pixel units) 
A=0;    % oblique incidence by linear grating (1=oblique 0=collinear) 
B=1;    % expanding by fresnel lens (1=yes 0=no) 
L=7;   % topological charge 
Lambda=30;  % linear grating fringe spacing (in pixels) 
aspect=1/2;  % fraction of fringe period that is white/clear 
xsize=1024;  % resolution (xres x yres number data pts calculated) 
ysize=768;  % 

% define the X and Y ranges (defined to skip zero) 
xvec = linspace(-xsize/2, xsize/2, xsize);  % list of x values 
yvec = linspace(-ysize/2, ysize/2, ysize);  % list of y values 

% define the meshes - matrices linear in one dimension 
[xmesh, ymesh] = meshgrid(xvec, yvec); 

% calculate the individual phase components 
vortexPh = atan2(ymesh,xmesh);  % the vortex phase 
linPh = -2*pi*ymesh;   % a phase of linear grating 
radialPh = (xmesh.^2+ymesh.^2);  % a phase of defocus 

% combine the phases with appropriate scales (phases are additive) 
% the 'pi' at the end causes inversion of the pattern 
Ph = L*vortexPh + A*linPh/Lambda + B*radialPh/R^2; 

% transmittance function (the real part of exp(I*Ph)) 
T = cos(Ph); 
% the binary version 
binT = T > cos(pi*aspect); 

% plot the pattern 
% imagesc(binT) 
imagesc(T) 
colormap(gray)