-
Notifications
You must be signed in to change notification settings - Fork 0
/
patent2_InAndOut.m~
192 lines (172 loc) · 8.41 KB
/
patent2_InAndOut.m~
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
function patent2_InAndOut
% ��m�ļ���������ڲ�������������Ż���
% �����ƹ�������ͼ��ÿ��ͼ��ȡ10���㡣�����ÿ��ͼ��ĵ�Ӧ����Ȼ���������Ӧ�������һ��
% ��V*b=0���b����Ȼ����b�ֽ���ڲ�������A�������A�͵�Ӧ����������������
global COUNT;
n = 20;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% �Ա�ͼƬ-����ͼƬ����DIC��� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 0
height_ref = 750;
width_ref = 750;
pad = 75;
pathname = 'image/ori_0315/';
filename_ref = 'speckle_pattern_4000_pad_0111_20_15000_750.png';
filename_roi = 'ROI_750_750_75.png';
radius = 40; %DIC subset radius
spacing = 0;
% ���ROI��pad�б仯���ǵ���ncorr_auto_initseeds.m�����pad
mat_name = ['speckle_'];
for k=1:1:n
filename_cur = ['temp_',num2str(k),'.jpg'];
displacements = fun_dic(pathname, filename_ref, filename_roi ,filename_cur, radius, spacing, mat_name, k);
% save(['results/ori_ref_mat/speckle_',num2str(k),'.mat'],'displacements');
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ��DIC���������ȡ���Ƶ� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1
fid = fopen('speckle_ori.txt','w');
for k = 1:1:n
str=strcat('results/ori_ref_mat/speckle_',num2str(k),'.mat');
load(str);
for v_loc = 30:15:120
for u_loc = 30:15:120
% 19.2 = (1728-2*96)/(140+1)��ȥ���߿����Ƶ��ļ��
u_bias = displacements.plot_u(v_loc * 5, u_loc * 5);
v_bias = displacements.plot_v(v_loc * 5, u_loc * 5);
uu = u_loc * 5 + u_bias;
vv = v_loc * 5 + v_bias;
fprintf(fid, '%f %f\n',uu,vv);
end
end
end
fclose(fid);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ���´��ı��ж�����������Ͷ�Ӧ��ͼ������ %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1
M=load('board.txt'); %��ȡ��������
M=[M';ones(1,49)];
m_all=load('speckle_ori.txt'); %��ȡ�ؼ�������
m_one=ones(3,49,n);
for i=1:1:n
m_temp = m_all((i-1)*49+1:i*49,:);
m_one(:,:,i) = [m_temp';ones(1,49)];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ������ⵥӦ�����b %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = ones(3,3,n); %ÿ��ͼһ����Ӧ����
for i=1:1:n
H(:,:,i)=homography(m_one(:,:,i),M);
end
V=ones(2*n,6); %V*b=0
for i=1:n %��V
V(2*i-1,:)=[H(1,1,i)*H(1,2,i) H(1,1,i)*H(2,2,i)+H(2,1,i)*H(1,2,i) H(2,1,i)*H(2,2,i) ...
H(3,1,i)*H(1,2,i)+H(1,1,i)*H(3,2,i) H(3,1,i)*H(2,2,i)+H(2,1,i)*H(3,2,i) H(3,1,i)*H(3,2,i)];
p1=[H(1,1,i)^2 H(1,1,i)*H(2,1,i)+H(2,1,i)*H(1,1,i) H(2,1,i)^2 H(3,1,i)*H(1,1,i)+H(1,1,i)*H(3,1,i) H(3,1,i)*H(2,1,i)+H(2,1,i)*H(3,1,i) H(3,1,i)^2];
p2=[H(1,2,i)^2 H(1,2,i)*H(2,2,i)+H(2,2,i)*H(1,2,i) H(2,2,i)^2 H(3,2,i)*H(1,2,i)+H(1,2,i)*H(3,2,i) H(3,2,i)*H(2,2,i)+H(2,2,i)*H(3,2,i) H(3,2,i)^2];
V(2*i,:)=p1-p2;
end;
[u s v]=svd(V); %�������ֽ���b
b=v(:,6);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ���·ֽ��ڲ������� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cy=(b(2)*b(4)-b(1)*b(5))/(b(1)*b(3)-b(2)^2);
lamda=b(6)-(b(4)^2+(b(2)*b(4)-b(1)*b(5))/(b(1)*b(3)-b(2)^2)*(b(2)*b(4)-b(1)*b(5)))/b(1);
kx=sqrt(lamda/b(1));
ky=sqrt((lamda*b(1))/(b(1)*b(3)-b(2)^2));
gama=-(b(2)*kx^2*ky)/lamda;
cx=(gama*cy)/kx-(b(4)*kx^2)/lamda;
%A=[kx gama cx;0 ky cy;0 0 1]; %����ڲ�������
A=[kx 0 cx;0 ky cy;0 0 1] %����ڲ�������
mp=ones(3,49,n); %ͼ��������ϳ���ά��ʽ������ѭ��������
mp = m_one;
%para1=[kx gama cx ky cy ];
%A = [2500, 0, 512; 0, 2500, 512; 0, 0, 1]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ���������������� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inv(A)*H(:,1,1);
inv(A)*H(:,2,1);
q1=(norm(inv(A)*H(:,1,1))+norm(inv(A)*H(:,2,1)))/2;
q1=(norm(inv(A)*H(:,1,1),2)+norm(inv(A)*H(:,2,1),2))/2;
R=[];
Rm=[];
for i=1:n
q1=(norm(inv(A)*H(:,1,i))+norm(inv(A)*H(:,2,i)))/2; %����������������
r1=(1/q1)*inv(A)*H(:,1,i);
r2=(1/q1)*inv(A)*H(:,2,i);
r3=cross(r1,r2);
RL=[r1 r2 r3];
[u s v]=svd(RL);
RL=u*v'; %��ת����������
p=(1/q1)*inv(A)*H(:,3,i);
RT=[r1 r2 p];
R=[R;RT];
sq = sqrt(RL(3,2)^2+RL(3,3)^2);
Q1 = atan2(RL(3,2),RL(3,3));
Q2 = atan2(-RL(3,1),sq);
Q3 = atan2(RL(2,1),RL(1,1));
[Q1 Q2 Q3]';
p;
R_new=[Q1,Q2,Q3,p'];
% R_new=[rotationVectors',p'];
Rm=[Rm , R_new];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ��ʽ�����ͶӰ��� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
res = ones(3,49,n);
for i=1:n
RT=R([(i-1)*3+1 : (i-1)*3+3],:);
x=A*RT*M;
x=[x(1,:)./x(3,:) ; x(2,:)./x(3,:); x(3,:)./x(3,:)]; % ��֤���һ��Ϊ1
res(:,:,i) = mp(:,:,i)-x;
end
res_sum = zeros(1,98*n);
for i = 1:n
for j = 1:49
res_sum(1,2*((i-1)*49+j)-1) = res(1,j,i);
res_sum(1,2*((i-1)*49+j)) = res(2,j,i);
end
end
res_sum.*res_sum;
xxx=sum(res_sum.*res_sum);
initres = sqrt(sum(res_sum.*res_sum,2)/(n*49))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% �������map�Ż��������ڲ�����������ת�Ƕȣ�ƽ�������� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1
parpool(20);
COUNT = 0;
para=[Rm,A(1,1),A(1,3),A(2,2),A(2,3)];%�Ż��ڲ������������ƽ�ƾ��������ת�ǣ�
%options = optimset('Algorithm','levenberg-marquardt','InitDamping',1e2,'Display','iter');
options = optimoptions(@lsqnonlin,'Algorithm','levenberg-marquardt','Display','iter');
[x,nor,res,exitflag,output] = lsqnonlin( @fun_mapopt, para, [],[],options);
A=[x(n*6+1) 0 x(n*6+2); 0 x(n*6+3) x(n*6+4); 0,0,1];
disp('map�����Ż������ڲ���Ϊ');
disp(A);
res;
exitflag;
output;
%fun_saveresult(x);
end
if 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ������������ѱ궨�����Ż��������ڲ�����������ת�Ƕȣ�ƽ�������� %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
para=[Rm,A(1,1),A(1,3),A(2,2),A(2,3)];%�Ż��ڲ������������ƽ�ƾ��������ת�ǣ�
options = optimset('Algorithm','levenberg-marquardt','Display','iter');
[x,nor,res] = lsqnonlin( @fun2, para, [],[],options, mp, M, n);
A=[x(n*6+1) 0 x(n*6+2); 0 x(n*6+3) x(n*6+4); 0,0,1];
disp('�����ѱ궨���Ż������ڲ���Ϊ');
disp(A);
sqrt(nor/(49*n))
end