白光干涉三维形貌恢复算法

这可能是白光干涉领域目前全网唯一的开源代码，废话不多说，直接上。质心法、相移法等较为简单的就不附上来了。

傅里叶变换法

这里展示的是通过重心法的方式来取包络峰值，也可以采取拟合的方法提取包络峰值

%傅里叶变换法
clear all;clc;
st=0;
en=116; %干涉图片总张数
tic;
disp('Step1 Starting...')
steplength=72;%扫描步长,单位nm
for ii=1:(en-st)
    n=st+ii;
    file=['F:\WLI\pics\test\1\',num2str(n) '.bmp'];%读取干涉图路径
    I=im2double(im2gray(imread(file)));
    [e,r]=size(I); 
    for jj=1:e
        for k=1:r
            I2(jj,k,ii)=I(jj,k); %方便接下来的循环处理
        end
    end
end

I2=double(I2);
disp('Step2 Starting...');
h=zeros(e,r);
t=1:en;

parfor jj=1:e
    for k=1:r
        y=zeros(1,en);
        max1=double(0),min1=double(0),sum=double(0);   
        for ii=1:en
            y(ii)=I2(jj,k,ii);
        end
        y1=fftshift(fft(y)); %傅里叶变化+频移
        y1(1:length(y1)/2+1)=0;%移除零频率和负频
        res=abs(ifft(y1)); %傅里叶逆变换取包络
        temp1=double(0),temp2=double(0);  %重心法提取包络峰值，也可以用曲线拟合包络
        for ii=1:length(res)
             mm=res(ii)^2;
             temp1=temp1+mm*ii;   
             temp2=temp2+mm; 
        end
        h(jj,k)=steplength*temp1/temp2; %扫描步长为72nm
    end
end

toc;


figure(1);surf(h);xlabel('Pixels');shading interp;colormap('default');toc;


                

空间频域分析法(FDA)

 de Groot提出的算法

%%FDA算法
clear all;close all;clc;
st=0;
en=116; %干涉图片张数
steplength=0.07;%扫描波长 单位微米
l=0.56 %中心波长 单位微米

tic;
disp('Step1 Starting...')
    for ii=1:(en-st)
        n=st+ii;
       file=['F:\WLI\pics\test\1\',num2str(n) '.bmp'];
        I=im2double(im2gray(imread(file)));
        [e,r]=size(I); 
        for jj=1:e
            for k=1:r
                I2(jj,k,ii)=I(jj,k); %方便接下来的循环处理
            end
        end   
        
    end

    I2=double(I2);
    disp('Step2 Starting...');
    h=zeros(e,r);

    t=1:en;
    parfor jj=1:e
        for k=1:r
            y=zeros(1,en);
            N=en;
            for ii=1:en
                y(ii)=I2(jj,k,ii);
            end
            y1=fftshift(fft(y)); %傅里叶变换+频移
            y2=y1(N/2+2:N);     %取正频分量

            fs=4*(steplength/l)*pi/steplength; 

            N=length(y2);
            t=fs*(1:N)/N;
                fhi=unwrap(angle(y2)); %求出相位并解包裹
                x=round(2*pi*N/(l*fs)); %找到中心波长对应包络的峰值位置
                xx=[t(x-3) t(x-2) t(x-1) t(x) t(x+1) t(x+2) t(x+3)]
                yy=[fhi(x-3) fhi(x-2) fhi(x-1) fhi(x) fhi(x+1) fhi(x+2) fhi(x+3)];
                f=polyfit(xx,yy,1); %线性拟合
                hc=(f(1)/2);
                h(jj,k)=1000*(hc);

        end
    end
    disp('Step3 Starting...');


 figure(1);surf(h);xlabel('X(pixel)');ylabel('Y(pixel)');zlabel('高度(nm)','FontSize',12);
 shading interp;colormap('default');toc;