heigth2d.m

来自「Matlab code for encoding an unwrapping p」· M 代码 · 共 49 行

%Construct test surface
x=[-4:0.08:4];
y=[-4:0.08:4];
pixelx=length(x)    ;
pixely=length(y)    ;
 hinitial=zeros(length(x),length(y));

    for i=1:length(x)
        for j=1:length(y)
                  hinitial(i,j)=30*(exp(-0.5*(x(i)+1).^2-0.5*(y(j)+1).^2) + 1.2*exp(-(x(i)-1).^2-(y(j)-1).^2));   %  fonction construite par 2 distributions gaussiennes
        end
    end
    
%       Plot original function hinitial

% figure(1)
% surf(x,y,hinitial)
%title('Original function h(x,y)')

%       Construct intensity function phi

lambda=10;
phi=(cos(2*pi.*hinitial/lambda)).^2;
%       Plot intensity function phi

figure(1)
imagesc(phi);colormap(gray);%title(['Intensity function phi(x) for lambda=',num2str(lambda)])

G1  =   mod(hinitial,2*pi)       ;
% %----ajout de bruit additif------------
moyenne   =0.6      ;
variance  =1        ;%pour une variance du bruit d閜asant 3 la recontruction n'est plus robuste
n = moyenne + (sqrt(variance) * randn(pixelx,pixely));%bruit de distribution normale
                                                %moyenne=0 variance=1 
G2  =   G1  +   n   ;% ajout d'un bruit additif

phder = zstep(G2)   ;
%---- plot --------------
figure(2);
subplot(141);imagesc(hinitial);title('Gaussienne')
colorbar
subplot(142);imagesc(G1);title('interferogram')
colorbar
subplot(143);imagesc(G2);title('ajout de bruit additif')
subplot(144);imagesc(phder);title('phase d閞oul閑')
colorbar
figure(3);
subplot(121);mesh(x,y,hinitial);title('Relief synth閠ique')
subplot(122);mesh(x,y,phder);title('Relief reconstruit')