1 简介
内联全息图外推算法matlab代码
2 部分代码
close all
clear all
% addpath('C:/Program Files/MATLAB/R2010b/myfiles');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N0 = 500; % initial number of pixels
N = 1000; % final number of pixels
Iterations = 100; % number of iterations
wavelength = 532*10^(-9); % wavelength in meter
screen = 0.035; % screen size in meter
z = 0.075; % source-to-screen distance in meter
z0 = 0.252*10^(-3); % source-to-sample distance in meter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reading cropped normalized hologram, where
% normalized hologram = (hologram with object)/(hologram without object)
fid = fopen('00_norm_div_cropped.bin', 'r');
hologram = fread(fid, [N0, N0], 'real*4');
fclose(fid);
measured = sqrt(hologram);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% reading support mask in the object domain
fid = fopen('00_support_object.bin', 'r');
support = fread(fid, [N, N], 'real*4');
fclose(fid);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creating initial complex-valued field distribution in the detector plane
amplitude = ones(N,N);
phase = zeros(N,N);
spherical = spherical_wave(N, wavelength, screen, z);
% phase = angle(spherical);
field_detector = amplitude.*exp(i*phase);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% creating spherical wavefront term
plane0 = N*wavelength*z/screen;
spherical0 = spherical_wave(N, wavelength, plane0, z0);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% iterative reconstruction
N1 = (N - N0)/2;
N2 = (N + N0)/2;
for kk = 1:Iterations
fprintf('Iteration: %d\n', kk)
% replacing the updated amplitudes with the known measured amplitudes
for ii = N1+1:N2-1
for jj = N1+1:N2-1
amplitude(ii,jj) = measured(ii-N1,jj-N1);
end
end
field_detector = amplitude.*exp(i*phase);
% reconstruction of transmission function t
t = IFT2Dc((FT2Dc(field_detector)).*spherical0);
% filtering object function
object = t - 1;
object = object.*support;
object_amplitude = abs(object);
object_phase = angle(object);
% smoothing every 5 iterations
R = rem(kk,5);
if (R == 0)
object_amplitude = smooth2D(object_amplitude);
end;
object = object_amplitude.*exp(i*object_phase);
t = object + 1;
% calculating complex-valued wavefront in the detector plane
field_detector_updated = FT2Dc((IFT2Dc(t)).*conj(spherical0));
amplitude = abs(field_detector_updated);
phase = angle(field_detector_updated);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% showing the reconstructed object
figure, imshow(flipud(rot90(abs(object))), [],'colormap', 1-gray);
title('Reconstructed object (amplitude) / a.u.')
xlabel({'x / px'})
ylabel({'y / px'})
axis on
set(gca,'YDir','normal')
colorbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% showing the extrapolated hologram
for ii = N1+1:N2-1
for jj = N1+1:N2-1
amplitude(ii,jj) = measured(ii-N1,jj-N1);
end
end
hologram_extrapolated = amplitude.^2;
figure, imshow(flipud(rot90(hologram_extrapolated)), []);
title('Extrapolated hologram / a.u.')
xlabel({'x / px'})
ylabel({'y / px'})
axis on
set(gca,'YDir','normal')
colormap('gray')
colorbar;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% saving the extrapolated hologram as jpg image
h = hologram_extrapolated;
h = (h - min(min(h)))/(max(max(h)) - min(min(h)));
imwrite (rot90(h), 'hologram_extrapolated.jpg');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [out] = FT2Dc(in)
[Nx Ny] = size(in);
f1 = zeros(Nx,Ny);
for ii = 1:Nx
for jj = 1:Ny
f1(ii,jj) = exp(i*pi*(ii + jj));
end
end
FT = fft2(f1.*in);
out = f1.*FT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [out] = IFT2Dc(in)
[Nx Ny] = size(in);
f1 = zeros(Nx,Ny);
for ii = 1:Nx
for jj = 1:Ny
f1(ii, jj) = exp(-i*pi*(ii + jj));
end
end
FT = ifft2(f1.*in);
out = f1.*FT;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% smoothing 2d images
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The code is written by Tatiana Latychevskaia, 2013
function [out] = smooth2D(in)
[N N] = size(in);
f = zeros(N,N);
f(N/2+1,N/2+1) = 4;
f(N/2-1+1,N/2+1) = 1;
f(N/2+1+1,N/2+1) = 1;
f(N/2-1+1,N/2-1+1) = 1;
f(N/2+1,N/2-1+1) = 1;
f(N/2+1+1,N/2-1+1) = 1;
f(N/2-1+1,N/2+1+1) = 1;
f(N/2+1,N/2+1+1) = 1;
f(N/2+1+1,N/2+1+1) = 1;
out = (1/12)*abs(IFT2Dc(FT2Dc(in).*FT2Dc(f)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [p] = spherical_wave(N, wavelength, area, z)
p = zeros(N,N);
delta = area/N;
for ii = 1:N;
for jj = 1:N
x = delta*(ii - N/2 -1);
y = delta*(jj - N/2 -1);
p(ii,jj) = exp(i*pi*(x^2 + y^2)/(wavelength*z));
end
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 仿真结果
4 参考文献
[1]宁冉, 李重光, 楼宇丽,等. 基于多图像处理单元的Matlab计算全息图快速算法[J]. 激光与光电子学进展, 2019, 56(5):6.
博主简介:擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真,相关matlab代码问题可私信交流。
部分理论引用网络文献,若有侵权联系博主删除。