3D File and Example Code

Example code of single wavelength reconstruction (to estimate the distance):

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
clc;
close all;
clear;
set(0,'defaultfigurecolor','w');%background color white
I=imread('G.tiff'); %Hologram
I0=imread('BG.tiff'); %Background
I = double(I); 
I0 = double(I0); 
[N,M]=size(I);
figure(1);imshow(I,[]);title('Original Hologram');       
I1=I;
Rect=round(getrect);                                                    
Xmin=max(Rect(1),1);                                                   
Ymin=max(Rect(2),1);       
Xmax=(Rect(1)+Rect(3));
Ymax=(Rect(2)+Rect(4));
hold on;
plot([Xmin Xmax Xmax Xmin Xmin],[Ymin Ymin Ymax Ymax Ymin],'r');       
hold off;
I1=I(Ymin:Ymax,Xmin:Xmax);   
[N,M]=size(I1);   
I0=I0(Ymin:Ymax,Xmin:Xmax);   
lamda=525*10^-9;
for kk = 0:1
 %d=(0.6+0.4*kk/10)*10^-3;   
 d=0.88*10^-3;%Estimated Distance
k=1.5*pi/lamda;
dx=2.2*10^-6;                                  
dy=2.2*10^-6;
dfx=1/(M*dx);                                     
dfy=1/(N*dy);
fx=ones(N,1)*(-M/2:M/2-1)*dfx;
fy=(ones(M,1)*(-N/2:N/2-1)*dfy)';
    hol_fft111=fftshift(fft2(I1));     
    hol_fft=fftshift(fft2(I1));                          
    obj_fft=hol_fft.*exp(1i*k*(-d)*sqrt(1-lamda^2*(fx.^2+fy.^2))); 
    obj=ifft2(ifftshift( obj_fft));
    amp_obj=abs(obj);
    obj_phase=angle(obj);
 figure(kk+3);imshow( amp_obj,[]);title('Intensity');  
 %figure(3);imshow( amp_obj,[]);title('Intensity');  
 %figure(4);imshow( obj_phase,[]);title('Phase');             
 
end
phase1=obj_phase;
AAA=I1./I0;
amp_hol=I1./I0;
 hol_phase=0;
 for kk = 1:10     
    hol=amp_hol.*exp(1i.*hol_phase);                  
    hol_fft=fftshift(fft2(hol));                        
    obj_fft=hol_fft.*exp(1i*k*(-d)*sqrt(1-lamda^2*(fx.^2+fy.^2))); 
    obj=ifft2(ifftshift( obj_fft));
    amp_obj=abs(obj);
    obj_phase=angle(obj);
    %obj_phase(amp_obj>0.58)=-obj_phase(amp_obj>0.58);     
    amp_obj(amp_obj>0.88)=1;
    %amp_obj=kkk;
    obj_1=amp_obj.*exp(1i.*obj_phase);
    obj_1_fft=fftshift(fft2(obj_1)); 
    hol_fft=obj_1_fft.*exp(1i*k*(d)*sqrt(1-lamda^2*(fx.^2+fy.^2))); 
    hol=ifft2(ifftshift(  hol_fft));
    hol_phase=angle(hol);
 end
figure(5);imshow(abs(obj),[]);title('Intensity'); 
figure(6);imshow(angle(obj),[]);title('Phase');  
I=abs(obj);
P=angle(obj);
mesh(P);