%%
clear all
clc
%%
lamd=632.8e-9; %设定入射光波长
R=50; %设定牛顿环曲率 R越大 条纹越宽
rm=25.5e-3; %设定干涉条纹区域 图像大小可以根据这个来定
x=0:0.0001:rm;
y=rm:-0.0001:0;
[X,Y]=meshgrid(x,y);
r2=X.^2+Y.^2;
phi=2*pi*(r2/R+lamd/2)/lamd; %相位差
I=4*cos(phi./2).^2; %第一象限干涉光强
N=255; %设定灰度等级
Ir2=(I/4.0)*N; %最大光强为最大灰度
% Ir1=fliplr(Ir2); %矩阵对称操作 flipur左右转
% Ir3=flipud(Ir1); %矩阵对称操作 flipud上下转
% Ir4=flipud(Ir2);
% Ir=[Ir1 Ir2;Ir3 Ir4]; %构造图像矩阵
I0=Ir2;
figure,image(I0,'XData',[-0.015,0.015],'YData',[0.015,-0.015]); %画干涉条纹
figure,image(I0);
colormap(gray(N));
% axis square
title('牛顿环'); %为图形添加标题
%% 二值化图牛顿
Img00=I0;
figure;imshow(mat2gray(Img00));
%% 计算x-方向载频大小
[height, width] = size(Img00);
xprofile = Img00(floor(height/2) , :);
fx0 = TemperalFrequency(xprofile, width);
%% 计算y-方向载频大小
yprofile = Img00(:, floor(width/2));
fy0 = -TemperalFrequency(yprofile, height);
%% 产生虚拟光栅
x=linspace(-1,1,width);
y=x'; %单引号求反
u=ones(size(y))*x; %定义x
v=y*ones(size(x)); %定义y
I1 =0.5*cos(fx0*pi*u + fy0*pi*v);
I2 =0.5*cos(fx0*pi*u + fy0*pi*v + pi/2);
I3 =0.5*cos(fx0*pi*u + fy0*pi*v + pi); %二维的移相,可去背景光强
I4 =0.5*cos(fx0*pi*u + fy0*pi*v + 1.5*pi);
I1=I1(1:256,1:256);
I2=I2(1:256,1:256); %规定横纵坐标的大小,取值 (大小自己定吗?)
I3=I3(1:256,1:256);
I4=I4(1:256,1:256);
figure, imshow(I1);
% subplot(2,2,1);imshow(I1);
% subplot(2,2,2);imshow(I2);
% subplot(2,2,3);imshow(I3);
% subplot(2,2,4);imshow(I4);
%% 虚拟光栅与干涉条纹形成莫尔条纹
Imgc1=double(I1).*double(Img00);
Imgc2=double(I2).*double(Img00);
Imgc3=double(I3).*double(Img00);
Imgc4=double(I4).*double(Img00);
figure,imshow(mat2gray(Imgc1));
figure,imshow(mat2gray(Imgc2));
figure,imshow(mat2gray(Imgc3)); %显示出四幅相移莫尔条纹
figure,imshow(mat2gray(Imgc4));
% subplot(2,2,1);imshow(mat2gray(Imgc1));
% subplot(2,2,2);imshow(mat2gray(Imgc2));
% subplot(2,2,3);imshow(mat2gray(Imgc3));
% subplot(2,2,4);imshow(mat2gray(Imgc4));
%% 提取低频信息
% % %设置滤波半径
D1=0.8; %%%%调试到最佳结果
%傅里叶变换并中心化
Fuv1=fftshift(fft2(Imgc1));
a1=abs(Fuv1);
figure,imshow(log(a1),[]);
x=1:256;y=1:256;
figure,mesh(x,y,a1(x,y));colormap(jet);colorbar;
%滤波阶数
n=3;
[M,N]=size(Fuv1);
%确定傅里叶变换的原点
xo=floor(M/2);
yo=floor(N/2);
%分别求4个不同滤波半径的ILPF(理想低通滤波器)
for i=1:M
for j=1:N
D=(sqrt((i-112)^2+(j-146)^2)); % 滤波区域,调试
h1(i,j)=1/(1+(D/D1)^(2*n));
%D=sqrt((i-168)^2+(j-168)^2);
%h1(i,j)=1/(1+(D/D1)^(2*n));
%D01=(sqrt((i-195)^2+(j-163)^2));
%D02=(sqrt((i-158)^2+(j-190)^2));
%h01(i,j)=1/(1+(D01/D1)^(2*n));
%h02(i,j)=1/(1+(D02/D1)^(2*n));
%h1=h01+h02;
end
end
%滤波矩阵点乘
Fuv2=fftshift(fft2(Imgc2));
Fuv3=fftshift(fft2(Imgc3));
Fuv4=fftshift(fft2(Imgc4));
Guv1=h1.*Fuv1;c=real(Guv1);figure,imshow(c);x=1:256;y=1:256;
figure,mesh(x,y,c(x,y));colormap(jet);colorbar;
Guv2=h1.*Fuv2;
Guv3=h1.*Fuv3;
Guv4=h1.*Fuv4;
%傅里叶逆变换
g1=ifft2(ifftshift(Guv1));
b1=mat2gray(real(g1));
g2=ifft2(ifftshift(Guv2));
b2=mat2gray(real(g2));
g3=ifft2(ifftshift(Guv3));
b3=mat2gray(real(g3));
g4=ifft2(ifftshift(Guv4));
b4=mat2gray(real(g4));
figure,imshow(b1);figure,imshow(b2);figure,imshow(b3);figure,imshow(b4);
%% 截取图像
d1=mat2gray(real(g1));
d2=mat2gray(real(g2));
d3=mat2gray(real(g3));
d4=mat2gray(real(g4));
figure,imshow(d1);figure,imshow(d2);figure,imshow(d3);figure,imshow(d4);
subplot(2,2,1);imshow(d1);
subplot(2,2,2);imshow(d2);
subplot(2,2,3);imshow(d3);
subplot(2,2,4);imshow(d4);
%% 求解相位
H=(double(g4)-double(g2))./(double(g1)-double(g3));
xiangwei=atan(H);
xiangwei=real(xiangwei);
figure,imshow(mat2gray(xiangwei),[]);
%% 寻找误差点,并绕过误差点解包裹
[m,n]=size(xiangwei);
im_unwrapped=unwrap5(xiangwei,m,n);
z=real(im_unwrapped);
figure;imshow(mat2gray(z),[]);
%% 求前36项zernike多项式系数wws36
[wws36 idx a1]=zernikecoeff36(im_unwrapped);
%% 用前36项zernike多项式进行波面拟合
[z_fit]=zernike36_fit(wws36,idx,a1);
%% 求波相差的PV值及RMS值
[PV RMS]=PV_RMS(z_fit,idx);
评论0