85375563DSCM_v1.0.rar_cropddt_亚像素 拟合_图像匹配_整像素搜索_曲面拟合法

function varargout = subpixelEstimation_new(oxy,txy,cType,speType,np,f,targetImg) %give the u,v displacement result in subpixel accuracy %oxy ---- original coordinate %txy ---- target coordinate found by pixel level processing %cType ---- definition for calculation of c in this stage %speType ---- subpixel estimation method %np ---- number of points for fitting or interpolation of c distribution %f ---- source pattern %targetImg ---- target image % this function return the [u v] in subpixel accuracy %此处的坐标系是与整像素搜索相反的 switch (speType) case 1 % 二维显式6个点,第6个点左上 % [ex ey cc p]=calculateSp1(oxy,txy,cType,np,f,targetImg); [ex ey cc]=calculateSp1(oxy,txy,cType,np,f,targetImg); case 2 % 二维显式9个点 [ex ey cc]=calculateSp2(oxy,txy,cType,np,f,targetImg); case 3 % 二维高斯5个点 [ex ey cc]=calculateSp3(oxy,txy,cType,np,f,targetImg); case 4 % 二维高斯9个点,MA 简化+polyfit,错误 [ex ey cc]=calculateSp4(oxy,txy,cType,np,f,targetImg); case 5 % 一维曲线5个点 [ex ey cc]=calculateSp5(oxy,txy,cType,np,f,targetImg); case 6 % 一维曲线9个点 [ex ey cc]=calculateSp6(oxy,txy,cType,np,f,targetImg); case 7 % 二维隐式9个点 % [ex ey cc p]=calculateSp7(oxy,txy,cType,np,f,targetImg); [ex ey cc]=calculateSp7(oxy,txy,cType,np,f,targetImg); case 8 % 二维高斯9个点,原始 NaN [ex ey cc]=calculateSp8(oxy,txy,cType,np,f,targetImg); case 9 % 一维曲线9个点 polyfit [ex ey cc]=calculateSp9(oxy,txy,cType,np,f,targetImg); case 10 % 二维高斯9个点 简化， 错误 [ex ey cc]=calculateSp10(oxy,txy,cType,np,f,targetImg); case 11 % 二维隐式9个点 polyfit [ex ey cc]=calculateSp11(oxy,txy,cType,np,f,targetImg); case 12 % 二维隐式9个点 复杂 [ex ey cc]=calculateSp12(oxy,txy,cType,np,f,targetImg); case 13 % 二维显式6个点,第6个点左下 [ex ey cc]=calculateSp13(oxy,txy,cType,np,f,targetImg); case 14 % 二维显式6个点,第6个点右上 [ex ey cc]=calculateSp14(oxy,txy,cType,np,f,targetImg); case 15 % 二维显式6个点,第6个点右下 [ex ey cc]=calculateSp15(oxy,txy,cType,np,f,targetImg); otherwise end %%%%%%%%%%%%%% End of Method %%%%%%%%%%%%%%%%%%% %% u = txy(2) - oxy(2) + ey; v = txy(1) - oxy(1) + ex; varargout{1} = u; varargout{2} = v; varargout{3} = cc; % varargout{4} = p; %% 函数--相关系数拟合法 function varargout = calculateSp1(oxy,txy,cType,np,f,targetImg) ws = size(f,1); %np = 3; hfnp = floor(np/2);%np是用于进行相关系数拟合的点数 [pointsX pointsY]=meshgrid(txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp);%pointsX1是图像行向坐标 pointsX1(1:2)=pointsX(1:2); pointsX1(3:5)=pointsX(4:6); pointsX1(6)=pointsX(8); pointsY1(1:2,1)=pointsY(1:2); pointsY1(3:5,1)=pointsY(4:6); pointsY1(6,1)=pointsY(8); % [cValues p]= calculateC(cType,f,pointsX1,pointsY1,targetImg); [cValues]= calculateC(cType,f,pointsX1,pointsY1,targetImg); cValues=cValues(:); % c0 = cValues(1); % c1 = cValues(2); % c2 = cValues(4); % c3 = cValues(5); % c4 = cValues(6); % c5 = cValues(8); c0 = cValues(1); c1 = cValues(2); c2 = cValues(3); c3 = cValues(4); c4 = cValues(5); c5 = cValues(6); a0 = c1/2 - c3 + c5/2; a1 = c2/2 - c3 + c4/2; a2 = c0 - c1 - c2 + c3; a3 = c5/2 - c1/2; a4 = c4/2 - c2/2; a5 = c3; %% x轴是竖向 % a0 = c2/2 - c3 + c4/2; % a1 = c1/2 - c3 + c5/2; % a2 = c0 - c1 - c2 + c3; % a3 = c4/2 - c2/2; % a4 = c5/2 - c1/2; % a5 = c3; %% x轴是横向 aa = 4*a0*a1-a2^2; ex=(a2*a4-2*a1*a3)/aa; ey=(a2*a3-2*a0*a4)/aa; % ex=e(2); % ey=e(1); cc = a0*ex*ex + a1*ey*ey + a2*ex*ey + a3*ex + a4*ey + a5; varargout{1}=ex; varargout{2}=ey; varargout{3}=cc; % varargout{4}=p; %% function varargout = calculateSp2(oxy,txy,cType,np,f,targetImg) ws = size(f,1); %np = 3; hfnp = floor(np/2);%np是用于进行相关系数拟合的点数 [pointsX1 pointsY1]=meshgrid( txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp ); cValues = calculateC(cType,f,pointsX1,pointsY1,targetImg); [pointsX pointsY]=meshgrid(-hfnp:hfnp, -hfnp:hfnp); xa1=pointsX(:); ya1= pointsY(:); A=[xa1.^2 ya1.^2 xa1.*ya1 xa1 ya1 ones(size(xa1))]; C=cValues(:); % para=inv(A'*A)*A'*C; A1=Tr(A); A2=Ni(A1*A); para=A2*A1*C; a=[2*para(1) para(3);para(3) 2*para(2)]; b=Ni(a); c=[-para(4);-para(5)]; e=b*c; ex=e(1); ey=e(2); cc = para(1)*ex*ex + para(2)*ey*ey + para(3)*ex*ey + para(4)*ex + para(5)*ey + para(6); varargout{1}=ex; varargout{2}=ey; varargout{3}=cc; %% function varargout = calculateSp3(oxy,txy,cType,np,f,targetImg) ws = size(f,1); %np = 3; hfnp = floor(np/2);%np是用于进行相关系数拟合的点数 [pointsX1 pointsY1]=meshgrid( txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp ); cValues = calculateC(cType,f,pointsX1,pointsY1,targetImg); % [pointsX pointsY]=meshgrid(-hfnp:hfnp, -hfnp:hfnp); % xa(1)=pointsX(2); % xa(2:4)=pointsX(4:6); % xa(5)=pointsX(8); % xa1=Tr(xa); % pointsY= pointsY(:); % ya(1)= pointsY(2); % ya(2:4)= pointsY(4:6); % ya(5)= pointsY(8); % ya1=Tr(ya); % A=[xa1.^2 ya1.^2 xa1 ya1 ones(size(xa1))]; CValues=log(abs(cValues(:)));%???有问题 c0 = cValues(2); c1 = cValues(4); c2 = cValues(5); c3 = cValues(6); c4 = cValues(8); a0 = c0/2 - c2 + c4/2; a1 = c1/2 - c2 + c3/2; a2 = c4/2 - c0/2; a3 = c3/2 - c1/2; a4 = c2; sigma2 = 1/(-2*a0); miu1 = a2*sigma2; miu2 = a3*sigma2; ex=miu1; ey=miu2; cc = exp(a0*ex*ex+a1*ey*ey+a2*ex+a3*ey+a4); varargout{1}=ex; varargout{2}=ey; varargout{3}=cc;%% function varargout = calculateSp4(oxy,txy,cType,np,f,targetImg) ws = size(f,1); %np = 3; hfnp = floor(np/2); [pointsX pointsY] = meshgrid( txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp ); cValues = calculateC(cType,f,pointsX,pointsY,targetImg); [pointsX pointsY] = meshgrid( txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp ); cValues = calculateC(cType,f,pointsX,pointsY,targetImg); cValues = log(cValues); c0 = cValues(1); c1 = cValues(2); c2 = cValues(3); c3 = cValues(4); c4 = cValues(5); c5 = cValues(6); c6 = cValues(7); c7 = cValues(8); c8 = cValues(9); a0 = 1/4*c0-1/2*c1+1/4*c2-1/2*c3+c4-1/2*c5+1/4*c6-1/2*c7+1/4*c8; a1 = -1/4*c0+1/2*c1-1/4*c2+1/4*c6-1/2*c7+1/4*c8; a2 = -1/4*c0+1/4*c2+1/2*c3-1/2*c5-1/4*c6+1/4*c8; a3 = 1/2*c3-c4+1/2*c5; a4 = 1/2*c1-c4+1/2*c7; a5 = 1/4*c0-1/4*c2-1/4*c6+1/4*c8; a6 = -1/2*c3+1/2*c5; a7 = -1/2*c1+1/2*c7; a8 = c4; %% c = a0*X^2*Y^2 + a1*X^2*Y + a2*X*Y^2 + a3*X^2 + a4*Y^2 + ... %% +a5*X*Y + a6*X + a7*Y + a8 %% Cx = 0 %% Cy = 0 %%% Initial value from the 1D method %%% x = -hfnp:hfnp; nx1 = 1 + hfnp*np; ny1 = 1 + hfnp; nx = [nx1:nx1+np-1]; ny = [ny1:np:ny1+(np-1)*np]; yx = cValues(nx);%%%the C on x direction yy = cValues(ny);%%%the C on y direction p = polyfit(x,yx,2); ex = -0.5*p(2)/p(1); p = polyfit(x,yy,2); ey = -0.5*p(2)/p(1); c(1) = 2*a0*ex*ex + 2*a2*ex + 2*a4; c(2) = a1*ex*ex + a5*ex + a7; ey = -c(2)/c(1); %cross search on y c(1) = 2*a0*ey*ey + 2*a1*ey + 2*a3; c(2) = a2*ey*ey + a5*ey + a6; ex = -c(2)/c(1); cc = 0; varargout{1}=ex; varargout{2}=ey; varargout{3}=cc; %% function varargout = calculateSp5(oxy,txy,cType,np,f,targetImg) ws = size(f,1); %np = 3; hfnp = floor(np/2);%np是用于进行相关系数拟合的点数 [pointsX pointsY] = meshgrid(txy(1)-hfnp:txy(1)+hfnp, txy(2)-hfnp:txy(2)+hfnp ); cValues = calculateC(cType,f,pointsX,pointsY,targetImg); C=cValues(:); c0 = cValues(2); c1 = cValues(4); c2 = cValues(5); c3 = cVa