【物理应用】内联全息图外推算法matlab代码-CFANZ编程社区

1 简介

内联全息图外推算法matlab代码

2 部分代码

close allclear all% addpath('C:/Program Files/MATLAB/R2010b/myfiles');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%N0 = 500;                  % initial number of pixelsN = 1000;                  % final number of pixelsIterations = 100;          % number of iterationswavelength = 532*10^(-9);  % wavelength in meterscreen = 0.035;            % screen size in meterz = 0.075;                 % source-to-screen distance in meterz0 = 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 planeamplitude = 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 termplane0 = N*wavelength*z/screen;spherical0 = spherical_wave(N, wavelength, plane0, z0);%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% iterative reconstruction    N1 = (N - N0)/2;    N2 = (N + N0)/2;for kk = 1:Iterationsfprintf('Iteration: %d\n', kk)    % replacing the updated amplitudes with the known measured amplitudes  for ii = N1+1:N2-1for jj = N1+1:N2-1  amplitude(ii,jj) = measured(ii-N1,jj-N1);endend  field_detector = amplitude.*exp(i*phase);% reconstruction of transmission function tt = IFT2Dc((FT2Dc(field_detector)).*spherical0);% filtering object functionobject = t - 1; object = object.*support;object_amplitude = abs(object);object_phase = angle(object);% smoothing every 5 iterationsR = 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 planefield_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 hologramfor ii = N1+1:N2-1for jj = N1+1:N2-1  amplitude(ii,jj) = measured(ii-N1,jj-N1);endend  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));    endendFT = 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));    endendFT = ifft2(f1.*in);out = f1.*FT;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% smoothing 2d images%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% The code is written by Tatiana Latychevskaia, 2013function [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));    endend;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%