qfw_ccpmp.m

这是在MATLAB下编写的基于彩色条纹的相位测量轮廓术的源码,很实用.,在光学三维测量中

clear all
close all
clc

%%%%以下程序用来计算有颜色物体的三维信息
%%%读取白色参考面的信息
[filename,pathname] = uigetfile('*.bmp','打开白光参考面图象......');
fName = strcat(pathname,filename);
White = double(imread(fName));
%%滤波
% w = gausswin(9)*gausswin(3)';
% w1 = imfilter(White(:,:,1),w,'circular');
% w2 = imfilter(White(:,:,2),w,'circular');
% w3 = imfilter(White(:,:,3),w,'circular');
% White = cat(3,w1,w2,w3);
%%%%
% ma1 = max(max(w1));
% mi1 = min(min(w1));
% %
% ma2 = max(max(w2));
% mi2 = min(min(w2));
% %
% ma3 = max(max(w3));
% mi3 = min(min(w3));
% %%%%%
% ww1 = (w1 - mi1)./(ma1 - mi1);
% ww2 = (w2 - mi2)./(ma2 - mi2);
% ww3 = (w3 - mi3)./(ma3 - mi3);
% Wwhite = cat(3,ww1,ww2,ww3);
%%%%%%%%%%%%%%读取红色的补色的信息
for i = 1:3
    [filename,pathname] = uigetfile('*.bmp','打开保存的红色的补色图象......');
    fName = strcat(pathname,filename);
    CR = double(imread(fName));
%     cr1 = imfilter(CR(:,:,1),w,'circular');
%     cr2 = imfilter(CR(:,:,2),w,'circular');
%     cr3 = imfilter(CR(:,:,3),w,'circular');
%     CR = cat(3,cr1,cr2,cr3);
%     CR = (White - CR) + (White - CR).*(1-Wwhite);
    CRF{i} = CR;
end
%%%%%
CRFe1 = CRF{1};
CRFe2 = CRF{2};
CRFe3 = CRF{3};
%%%%%
CRFe1r = CRFe1(:,:,1);
CRFe1g = CRFe1(:,:,2);
CRFe1b = CRFe1(:,:,3);
%%%%%
CRFe2r = CRFe2(:,:,1);
CRFe2g = CRFe2(:,:,2);
CRFe2b = CRFe2(:,:,3);
%%%%%
CRFe3r = CRFe3(:,:,1);
CRFe3g = CRFe3(:,:,2);
CRFe3b = CRFe3(:,:,3);
%%%%%计算背景光强
CRrI1 = (CRFe1r + CRFe2r + CRFe3r)./3;
% CRgI1 = (CRFe1g + CRFe2g + CRFe3g)./3;
% CRbI1 = (CRFe1b + CRFe2b + CRFe3b)./3;
%%%%计算振幅调制度
CRrI2 = ((3*((CRFe1r - CRFe3r).^2) + (2*CRFe2r - CRFe1r - CRFe3r).^2).^(1/2))./3;
CRgI2 = ((3*((CRFe1g - CRFe3g).^2) + (2*CRFe2g - CRFe1g - CRFe3g).^2).^(1/2))./3;
CRbI2 = ((3*((CRFe1b - CRFe3b).^2) + (2*CRFe2b - CRFe1b - CRFe3b).^2).^(1/2))./3;
%%%%计算影响因子
Crg = sum(sum(CRgI2))/sum(sum(CRrI2));
Crb = sum(sum(CRbI2))/sum(sum(CRrI2));
% %%%%计算调制度
% CRrM = CRrI2./CRrI1;
% CRgM = CRgI2./CRgI1;
% CRbM = CRbI2./CRbI1;
% %%%%%计算影响因子
% Crg = sum(sum(CRgM))./sum(sum(CRrM));
% Crb = sum(sum(CRbM))./sum(sum(CRrM));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear CRFe1 CRFe2 CRFe3 CRFe1r CRFe1g CRFe1b CRFe2r CRFe2g CRFe2b CRFe3r CRFe3g CRFe3b
% clear CRgI1 CRbI1 CRgI2 CRbI2 w1 w2 w3 
%%%%%%%%%%%%%%读取绿色的补色的信息
for i = 1:3
    [filename,pathname] = uigetfile('*.bmp','打开保存的绿色的补色图象......');
    fName = strcat(pathname,filename);
    CR = double(imread(fName));
%     cr1 = imfilter(CR(:,:,1),w,'circular');
%     cr2 = imfilter(CR(:,:,2),w,'circular');
%     cr3 = imfilter(CR(:,:,3),w,'circular');
%     CR = cat(3,cr1,cr2,cr3);
%     CR = (White - CR) + (White - CR).*(1 - Wwhite);
    CRF{i} = CR;
end
%%%%%
CGFe1 = CRF{1};
CGFe2 = CRF{2};
CGFe3 = CRF{3};
%%%%%
CGFe1r = CGFe1(:,:,1);
CGFe1g = CGFe1(:,:,2);
CGFe1b = CGFe1(:,:,3);
%%%%%
CGFe2r = CGFe2(:,:,1);
CGFe2g = CGFe2(:,:,2);
CGFe2b = CGFe2(:,:,3);
%%%%%
CGFe3r = CGFe3(:,:,1);
CGFe3g = CGFe3(:,:,2);
CGFe3b = CGFe3(:,:,3);
%%%%%计算背景光强
% CGrI1 = (CGFe1r + CGFe2r + CGFe3r)./3;
CGgI1 = (CGFe1g + CGFe2g + CGFe3g)./3;
% CGbI1 = (CGFe1b + CGFe2b + CGFe3b)./3;
%%%%计算振幅
CGrI2 = ((3*((CGFe1r - CGFe3r).^2) + (2*CGFe2r - CGFe1r - CGFe3r).^2).^(1/2))./3;
CGgI2 = ((3*((CGFe1g - CGFe3g).^2) + (2*CGFe2g - CGFe1g - CGFe3g).^2).^(1/2))./3;
CGbI2 = ((3*((CGFe1b - CGFe3b).^2) + (2*CGFe2b - CGFe1b - CGFe3b).^2).^(1/2))./3;
%%%计算影响因子
Cgr = sum(sum(CGrI2))/sum(sum(CGgI2));
Cgb = sum(sum(CGbI2))/sum(sum(CGgI2));
% %%%%计算调制度
% CGrM = CGrI2./CGrI1;
% CGgM = CGgI2./CGgI1;
% CGbM = CGbI2./CGbI1;
% %%%%%计算影响因子
% Cgr = sum(sum(CGrM))./sum(sum(CGgM));
% Cgb = sum(sum(CGbM))./sum(sum(CGgM));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear CGFe1 CGFe2 CGFe3 CGFe1r CGFe1g CGFe1b CGFe2r CGFe2g CGFe2b CGFe3r CGFe3g CGFe3b
% clear CGrI1 CGbI1 CGrI2 CGbI2    
%%%%%%%%%%%%%%读取蓝色的补色的信息
for i = 1:3
    [filename,pathname] = uigetfile('*.bmp','打开保存的蓝色的补色图象......');
    fName = strcat(pathname,filename);
    CR = double(imread(fName));
%     cr1 = imfilter(CR(:,:,1),w,'circular');
%     cr2 = imfilter(CR(:,:,2),w,'circular');
%     cr3 = imfilter(CR(:,:,3),w,'circular');
%     CR = cat(3,cr1,cr2,cr3);
%     CR = (White - CR) + (White - CR).*(1 - Wwhite);
    CRF{i} = CR;
end
%%%%%
CBFe1 = CRF{1};
CBFe2 = CRF{2};
CBFe3 = CRF{3};
%%%%%
CBFe1r = CBFe1(:,:,1);
CBFe1g = CBFe1(:,:,2);
CBFe1b = CBFe1(:,:,3);
%%%%%
CBFe2r = CBFe2(:,:,1);
CBFe2g = CBFe2(:,:,2);
CBFe2b = CBFe2(:,:,3);
%%%%%
CBFe3r = CBFe3(:,:,1);
CBFe3g = CBFe3(:,:,2);
CBFe3b = CBFe3(:,:,3);
%%%%%计算背景光强
% CBrI1 = (CBFe1r + CBFe2r + CBFe3r)./3;
% CBgI1 = (CBFe1g + CBFe2g + CBFe3g)./3;
CBbI1 = (CBFe1b + CBFe2b + CBFe3b)./3;
%%%%计算振幅
CBrI2 = ((3*((CBFe1r - CBFe3r).^2) + (2*CBFe2r - CBFe1r - CBFe3r).^2).^(1/2))./3;
CBgI2 = ((3*((CBFe1g - CBFe3g).^2) + (2*CBFe2g - CBFe1g - CBFe3g).^2).^(1/2))./3;
CBbI2 = ((3*((CBFe1b - CBFe3b).^2) + (2*CBFe2b - CBFe1b - CBFe3b).^2).^(1/2))./3;
%%%%计算影响因子
Cbr = sum(sum(CBrI2))/sum(sum(CBbI2));
Cbg = sum(sum(CBgI2))/sum(sum(CBbI2));
%%%%%计算相位修正的强度比值
K1g = sum(sum(CGgI1))/sum(sum(CRrI1));
K1b = sum(sum(CBbI1))/sum(sum(CRrI1));
%%%%%%
K2g = sum(sum(CGgI2))/sum(sum(CRrI2));
K2b = sum(sum(CBbI2))/sum(sum(CRrI2));
% %%%%计算调制度
% CBrM = CBrI2./CBrI1;
% CBgM = CBgI2./CBgI1;
% CBbM = CBbI2./CBbI1;
% %%%%%计算影响因子
% Cbr = sum(sum(CBrM))./sum(sum(CBbM));
% Cbg = sum(sum(CBgM))./sum(sum(CBbM));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% clear CBFe1 CBFe2 CBFe3 CBFe1r CBFe1g CBFe1b CBFe2r CBFe2g CBFe2b CBFe3r CBFe3g CBFe3b
% clear CBrI1 CBgI1 CBrI2 CBgI2 cr1 cr2 cr3 CRF CR 
%%%%%%以下提取物体进行处理
[filename,pathname] = uigetfile('*.bmp','打开投射物体的图象......');
fName = strcat(pathname,filename);
CObjectImage = double(imread(fName));
% O1 = imfilter(CObjectImage(:,:,1),w,'circular');
% O2 = imfilter(CObjectImage(:,:,2),w,'circular');
% O3 = imfilter(CObjectImage(:,:,3),w,'circular');
% CObjectImage = cat(3,O1,O2,O3);
%%%%%%%去除背景色
[filename,pathname] = uigetfile('*.bmp','打开投射白光的物体图象......');
fName = strcat(pathname,filename);
White0 = double(imread(fName));
% O1 = imfilter(White0(:,:,1),w,'circular');
% O2 = imfilter(White0(:,:,2),w,'circular');
% O3 = imfilter(White0(:,:,3),w,'circular');
% White0 = cat(3,O1,O2,O3);
%%%%
% ma1 = max(max(O1));
% mi1 = min(min(O1));
% %
% ma2 = max(max(O2));
% mi2 = min(min(O2));
% %
% ma3 = max(max(O3));
% mi3 = min(min(O3));
% %%%%%
% OO1 = (O1 - mi1)./(ma1 - mi1);
% OO2 = (O2 - mi2)./(ma2 - mi2);
% OO3 = (O3 - mi3)./(ma3 - mi3);
% Wohite0 = cat(3,OO1,OO2,OO3);
%%%%
ObjectImage = (White0 - CObjectImage);
% ObjectImage = (White0 - CObjectImage) + (White0 - CObjectImage).*(1 - Wohite0);
%%%%%通道信息分离
RObject = ObjectImage(:,:,1);
GObject = ObjectImage(:,:,2);
BObject = ObjectImage(:,:,3);
%%%%%
CRObject = (RObject - GObject.*Cgr - BObject.*Cbr)./(1 - Cgr - Cbr);
CGObject = (GObject - RObject.*Crg - BObject.*Cbg)./(1 - Crg - Cbg);
CBObject = (BObject - GObject.*Cgb - RObject.*Crb)./(1 - Crb - Cgb);
%%%%%
w = gausswin(9)*gausswin(3)';
CRObject = imfilter(CRObject,w,'circular');
CGObject = imfilter(CGObject,w,'circular');
CBObject = imfilter(CBObject,w,'circular');

%%%%%
a1 = CGObject*(2*K2b+K1b)-CRObject*(2*K2b*K1g+K1b*K2g)-CBObject*(K1g-K2g);
a2 = sqrt(3)*(CBObject*(K2g+K1g)-CRObject*K2g*K1b-CGObject*K1b);
%%%%%
% a1 = K1b*K2g - K1g*K2b + (K1b*K2g + K1g*K2b + K2b*K2g).*((K1g*K2b - K1b*K2g).*CRObject - K2b.*CGObject + K2g.*CBObject);
% a2 = (sqrt(3)*K2b*K2g).*((K1g*K2b + K1b*K2g).*CRObject - K2b.*CGObject - K2g.*BObject);
ObjWrapPhai = atan2(a1,a2);
%%%%%
% maR = max(max(CRObject));
% miR = min(min(CRObject));
% CRObject = (CRObject - miR)*255./(maR - miR);
% %
% maG = max(max(CGObject));
% miG = min(min(CGObject));
% CGObject = (CGObject - miG)*255./(maG - miG);
% %
% maB = max(max(CGObject));
% miB = min(min(CGObject));
% CGObject = (CGObject - miB)*255./(maB - miB);
% %%%%%
% a1 = sqrt(3).*(CRObject - CBObject);
% a2 = 2*CGObject -CBObject - CRObject;
% ObjWrapPhai = atan2(a1,(a2+eps));
[Row,Column] = size(ObjWrapPhai);
mask = ones(Row,Column);
% [filename,pathname] = uigetfile('*.bmp','打开模板图象......');
% fName = strcat(pathname,filename);
% mask = double(imread(fName));
% maa = max(max(mask));
% mii = min(min(mask));
% mask = (mask-mii)./(maa-mii);
ObjUnwrapPhai = zlf_dunwrap(ObjWrapPhai,mask,[460,500]);
figure;mesh(ObjUnwrapPhai)
zz = ObjUnwrapPhai;
[xs,ys] = meshgrid(1:Column,1:Row);
aa = polyval2d(polyfit2d(xs,ys,zz,1,1),xs,ys);
figure;mesh(zz - aa)
clear BObject CBObject CBbI1 CBbI2 CGObject CGgI1 CGgI2 CObjectImage CRObject CRrI1 CRrI2 GObject
clear CBFe1 CBFe1b CBFe1g CBFe1b CBFe2 CBFe2b CBFe2g CBFe2r
clear CBFe3 CBFe3b CBFe3g CBFe3r CBgI2 CBrI2 
clear CGFe1 CGFe1b CGFe1g CGFe1r CGFe2 CGFe2b CGFe2g CGFe2r
clear CGFe3 CGFe3b CGFe3g CGFe3r CGbI2 CGrI2 
clear CR CRF CRFe1 CRFe1b CRFe1g CRFe1r CRFe2 CRFe2b CRFe2g CRFe2r CRFe3 CRFe3b CRFe3g CRFe3r
clear CRbI2 CRgI2 Column a1 a2 a3 fName filename i maa mask mii pathfile w xs ys