【物理应用】内联全息图外推算法matlab代码

close all clear all % addpath( ’C:/Program Files/MATLAB/R2010b/myfiles’ ); %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% N 0 = 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 z 0 = 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, [N 0 , N 0 ], ’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 plane 0 = N*wavelength*z/screen; spherical 0 = spherical_wave(N, wavelength, plane 0 , z 0 ); %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % iterative reconstruction N1 = (N - N 0 )/ 2 ; N2 = (N + N 0 )/ 2 ; for kk = 1 :Iterations fprintf( ’Iteration: %d ’ , 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)).*spherical 0 ); % 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(spherical 0 )); amplitude = abs(field_detector_updated); phase = angle(field_detector_updated); end %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% % showing the reconstructed object figure, imshow(flipud(rot9 0 (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(rot9 0 (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 (rot9 0 (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 2 d 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; % %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %%% %