MATLAB: Gerchberg–Saxton algorithm

digital image processingimage processingiteration

Hello! I have the following code:
for k=1:1:20;
G_pr=absP.*exp(i.*theta); 
g_pr=ifft2(ifftshift(G_pr)); 
absPhase=abs(angle(g_pr)); 
maxPh=max(max(absPhase)); 
minPh=min(min(absPhase)); 
g_pr(absPhase>=(minPh+0.2*(maxPh-minPh)))=0; 
g_pr=real(g_pr); 
gg=255*g_pr/(max(max(g_pr))); 
figure(1),imshow(uint8(gg)); title(num2str(k)) 
G=fftshift(fft2(g_pr)); 
G=G./abs(G); 
theta=angle(G); 
end
The first theta is the phase of my model image (angle(model)) However, this code diverges instead of converge, Does someone knows why?
Thank you

Best Answer

I'm not sure where in your code is the error. The following code works perfectly for me (adopted from wikipedia):
A = fftshift(ifft2(fftshift(Target)));
for i=1:25
  B = abs(Source) .* exp(1i*angle(A));
  C = fftshift(fft2(fftshift(B)));
  D = abs(Target) .* exp(1i*angle(C));
  A = fftshift(ifft2(fftshift(D)));
    imagesc(abs(C)) %Present current pattern
    title(sprintf('%d',i));
    pause(0.5)
end
Before running the code, make sure 'Source' contains your input beam, for example:
Source = exp(-1/2*(xx0.^2+yy0.^2)/sigma^2);
And 'Target' contains your requested pattern.
The phase mask can be presented at the end of the for loop:
imagesc(angle(A))
Related Question