逐行解释matlab代码clc;
close all;
clear all;
load ('E:\镜面实验数据\20240816\Calib_Results.mat','fc','cc','KK');
data=load('E:\镜面实验数据\20240816\pmxs.txt');
R = load('E:\镜面实验数据\20240816\R.txt');
Point = load('E:\镜面实验数据\20240816\P.txt');
phase_shift_num=4;
l_imgs = cell(1,phase_shift_num);
for i=1
for j=1:phase_shift_num
l_imgs{i,j} = double(imread(['E:\镜面实验数据\20240816\wuti2\' num2str(phase_shift_num*(i-1)+j),'.bmp']));
end
end
%%求极线参数
a=data(1,1);b=data(1,2);c=data(1,3);d=data(1,4);
D=abs(d)/(sqrt(a*a+b*b+c*c));
n=[a,b,c]';
I=[1,0,0;0,1,0;0,0,1];
D1=I-2*n*n';
D2=2*D*n;
N0 = [-1,0,0;0,1,0;0,0,1];
I0=[-1;1;1];
R=R;T=Point*R;
fx=fc(1,1);fy=fc(2,1);u0=cc(1,1);v0=cc(2,1); %内参
M=KK*[R,T'];
M0=M(1:3,1:3);m0=M(:,4);
% m=M0*D2+m0-M0*D1*inv(M0)*m0;
% mx=[0,-m(3,1),m(2,1);m(3,1),0,-m(1,1);-m(2,1),m(1,1),0];
% F=mx*M0*D1*inv(M0);
M01 = inv(M0);
m=m0-M0*N0*M01*m0;
mx=[0,-m(3,1),m(2,1);m(3,1),0,-m(1,1);-m(2,1),m(1,1),0];
F = mx*M0*N0*M01;
c1 = -inv(M0)*m0 ; c2 = [-1,0,0;0,1,0;0,0,1]*c1;
K = R*[0;0;1];
V = c2 - c1;
% 将向量 V 单位化
r1 = V / norm(V);
% r2 = R*[0;1;0];
r2 = cross(K,r1);
r3 = cross(r1,r2);
R1 = [r1';r2';r3']';
R2 = -R1*c2;
M2 = KK*[R1 R2];
left_M2 = M2(:,1:3);
Ma = left_M2*inv(M0);
Ma = double(Ma);
% 假设 Ma 是每张图像对应的 3x3 仿射变换矩阵,这里只是一个示例
Mz = cell(1, phase_shift_num); % 创建一个单元数组来存储每个相移图像的变换矩阵
for j = 1:phase_shift_num
Mz{j} = Ma;% 计算或加载 Ma{j}
end
% 初始化校正后的图像数组
l_imgs_corrected = cell(1, phase_shift_num);
% 对每张图像进行校正
for j = 1:phase_shift_num
% 读取原始图像
img = l_imgs{1, j};
% 计算校正后图像的边界框
s1 = size(img);
corners = [0, 0, s1(2), s1(2); 0, s1(1), 0, s1(1)];
corners_x = p2t(Mz{j}, corners);
minx = floor(min(corners_x(1,:)));
maxx = ceil(max(corners_x(1,:)));
miny = floor(min(corners_x(2,:)));
maxy = ceil(max(corners_x(2,:)));
bb1 = [minx; miny; maxx; maxy];
% 应用校正变换
[img_corrected, alpha] = imwarp(img, Mz{j}, 'Bilinear', bb1);
% 存储校正后的图像
l_imgs_corrected{j} = double(img_corrected);
end
ROIL=[];
imgL=[];
for i=1:phase_shift_num
imgL(:,:,i)= l_imgs_corrected{1,i};
end
ROILg=max(imgL,[],3);
figure,imshow(ROILg,[])
%闭运算
freq = 64;
ROI1=imclose(ROILg,strel('disk',10));
ROI1 = ROI1<80;
ROI1 = ~ROI1;
[ROIR,ROIL]=findLargeBlocks_D(ROI1, 8);
figure,imshow(ROI1,[])
figure,imshow(ROIL,[])
figure,imshow(ROIR,[])
l_uph1=calculatePhase(l_imgs_corrected,phase_shift_num,freq,10);
l_uph1 = double(l_uph1);
figure,imshow(l_uph1,[]),title('相机截断相位');
l_uph1ROIL=l_uph1.*ROIL;
figure,imshow(l_uph1ROIL,[]),title('相机截断相位');
[LPx,~]=gradient(l_uph1ROIL);
figure,imshow(LPx,[]),title('左相机截断相位x梯度');
BWL = (LPx >= -4 & LPx <= -2.9) | (LPx >= 2.9 & LPx <= 4);
BWLXX=imclose(BWL,strel('rectangle',[3 100]));
% figure,imshow(BWLXX,[])
BWL=bwareaopen(BWL,10,8);
% BWL=imclose(BWL|~ROIL,strel('rectangle',[15 1]));
BWL=imdilate(BWL,strel('rectangle',[5 1]));%图形膨胀
BWL=imclose(BWL&ROIL,strel('disk',8));
BWL=imclose(BWL|~ROIL,strel('rectangle',[15 3]));
BWL=BWL&ROIL;
figure,imshow(BWL+ROIL*2,[]),title('左相机截断相位x梯度');
figure,imshow(BWL,[]),title('左相机截断相位x梯度');
LL=bwlabel(~BWL&ROIL);
LLA=regionprops(LL);
LLy=[];
for i=1:size(LLA,1)
LLy(i)=LLA(i).Centroid(2);
end
figure,imshow(LL,[]),title('左相机');
l_uph1ROIR=l_uph1.*ROIR;
figure,imshow(l_uph1ROIR,[]),title('相机截断相位');
[RPx,~]=gradient(l_uph1ROIR);
figure,imshow(RPx,[]),title('左相机截断相位x梯度');
BWR = (RPx >= -4 & RPx <= -2.9) | (RPx >= 2.9 & RPx <= 4);
BWRXX=imclose(BWR,strel('rectangle',[3 100]));
BWR=bwareaopen(BWR,10,8);
BWR=imdilate(BWR,strel('rectangle',[5 1]));%图形膨胀
BWR=imclose(BWR&ROIR,strel('disk',8));
BWR=imclose(BWR|~ROIR,strel('rectangle',[15 3]));
BWR=BWR&ROIR;
figure,imshow(BWR+ROIR*2,[]),title('左相机截断相位x梯度');
figure,imshow(BWR,[]),title('左相机截断相位x梯度');
LR=bwlabel(~BWR&ROIR);
LRA=regionprops(LR);
LRy=[];
for i=1:size(LRA,1)
LRy(i)=LRA(i).Centroid(2);
end
figure,imshow(LR,[]),title('右相机');
for i=1:size(LLy,2)
LineL=LL(round(LLy(i)),:);
LineR=LR(round(LLy(i)),:);
PmL{i,1}=unique(LineL);%数组中的唯一值
PmL{i,2}=unique(LineR);
end
for i=1:size(LRy,2)
LineL=LL(round(LRy(i)),:);
LineR=LR(round(LRy(i)),:);
PmR{i,1}=unique(LineL);
PmR{i,2}=unique(LineR);
end