# Stereo Digital Image Correlation in MATLAB
## 6 Conclusions
IndShift=-iftshift(-fix(SubSize/2):ceil(SubSize/2)-1);
7. u=IndShift(uid);
8. v=IndShift(vid);
9. else
10. u=NaN; v=NaN;function PD=ImgCorr(n,PD,FileNames,RefStrat,StopCritVal)
2. SubExtract=@(Mat,Xos,SubSize) Mat(Xos(2)-(SubSize-1)/2: Xos(2)+ (SubSize-1)/2,Xos(1)-(SubSize-1)/2:Xos(1)+(SubSize-1)/2);
3. for d=2:n, tic; % outer loop
4. G=im2 double(imread(FileNames{d}));
5. if all(PD(d).ImgFilt), G=imgaussfilt(G,PD(d).ImgFilt(1), 'FilterSize', PD(d).ImgFilt(2)); end
6. InterpCoef=griddedInterpolant({1:1:size(G,1),1:1:size(G,2)},G,'spline'); % Section 2.3.2
7. if any([RefStrat==1,d==2])
8. F=im2 double(imread(FileNames{d-1}));
9. if all(PD(d).ImgFilt), F=imagaussfilt(F,PD(d).ImgFilt(1), 'FilterSize',PD(d).ImgFilt(2)); end
10. [dFdx,dFdy]=imgradientxy(F,'prewitt');
11. PD(d).Xos(1,:)=PD(d-1). [PERSON](1,:)+fix(PD(d-1).P(1,:));
12. PD(d).Xos(2,:)=PD(d-1).Xos(2,:)+fix(PD(d-1).P(7,:));
13. [PD(d).P(1,:),PD(d).P(7,:)]=arrayfun(@(XosX,XosY,SubSize )PCM (F,G,SubSize,XosX,XosY,SubExtract),PD(d).Xos(1,:), PD(d).Xos(2,:),PD(d). SubSize); % Section 3.4.1
14. else
15. PD(d).P=PD(d-1).P;
16. end
17. P=NaN(size(PD(d).P)); C=NaN(size(PD(d).C)); Iter=NaN(size(PD(d).C)); StopVal=NaN(size(PD(d).C));
18. for q=1:size(PD(d).Xos,2) % inner loop (can be changed to parfor for parallel processing)
19. if (sum(ismam(PD(d).P(:,q)))==0)&&(sum(isman( PD(d).Xos(:,q)))==0)
20. [f,dfdx,dfdy,dX,dY]=SubShapeExtract(PD(d).SubSize(q), PD(d).SubShape(q,:),PD(d).Xos(:,q),F,dfdx,dFdy, SubExtract); % Section 3.2
21. [P(:,q),C(q),Iter(q),StopVal(q)]=SubCorr(InterpCoef,f,dfdx, dfdy,PD(d).SubSize(q),PD(d).SFOrder(q), PD(d).Xos(:,q),dX,dY,PD(d). P(:,q),StopCritVal); % Section 3.2
22. end
23. end
24. PD(d).P=P; PD(d).C=C; PD(d).IterIterIter; PD(d).StopVal=StopVal;
25. if rem(d-2,10)==0, fprintf('Image/Total| Time (s) | CC (min) | CC (mean) | Iter (max) \
'); end
26. fprintf(' %4.d/%4.d | %8.3f | %.6f | %.7f | %4.Of \
',d,n,toc,min (PD(d).C),nanmean(PD(d).C),max(PD(d).Iter));
27. end
1. function RD=CSTrans(n,RD,WorldCTs,ImCTs,RefStrat)
2. CamParams=estimateCameraParameters(ImgCTs,WorldCTs, 'NumRadialDistortionCoefficients',2); % Section 2.2
3. Qi=[CamParams.CameraParameters1. IntrinsicMatrix',[0; 0; 0]]*[[eye(3), [0;0;0]]; 0, 0, 0, 1];
4. Q2=[CamParams.CameraParameters2. IntrinsicMatrix',[0; 0; 0]]*[[CamParams.RotationOfCamera2', CamParams.TranslationOfCamera2']; 0, 0, 0, 1];
5. B=CamParams.FundamentalMatrix;
6. for d=1:n, tic
7. Xds1=RD.ProData1(d).Xos+[RD.ProcData1(d).P(1,:); RD.ProcData1(d). P(7,:)]; % Equation (49)
8. Xds2=RD.ProcData2(d).Xos+[RD.ProcData2(d).P(1,:); RD.ProcData2(d). P(7,:)]; % Equation (49)9. indValid=find((isnan(Xds1(1,:))+isnan(Xds1(2,:)) +isnan(Xds2(1,:))+isnan(Xds2(2,:)))==0);
10. if d==1|RefStrat=1
11. RD.DispTrans(d).Xow(:,indValid)=Triangulation(B,Q1,Q2, undistortPoints(RD.ProData1(d).Xos(:,indValid)',CamParams. CameraParameters1)',undistorPoints(RD.ProData2(d).Xos (:,indValid)',CamParams.CameraParameters2)');
12. else
13. RD.DispTrans(d).Xow=RD.DispTrans(d-1).Xow;
14. end
15. RD.DispTrans(d).Uw(:,indValid)=Triangulation(B,Q1,Q2, undistortPoints(Xds1(:,indValid)',CamParams.CameraParameters1)', undistortPoints(Xds2(:,indValid)',CamParams.CameraParameters2)')-RD. DispTrans(d).Xow(:,indValid);% Equation (51)
16. RD.DispTrans(d).CamParams=CamParams;
17. fprintf('CS transformation image: %d/%d\t{time:%.3f\
',d,n,toc};
18. end
1. function [ptsOut]=Triangulation(B,Q1,Q2,pts1,pts2)
2. for i=1:size(pts1,2)
3. T1 inv=[1, 0, pts1(1,i); 0, 1, pts1(2,i); 0, 0, 1]; % Equation (35)
4. T2 inv=[1, 0, pts2(1,i); 0, 1, pts2(2,i); 0, 0, 1]; % Equation (35)
5. B1=T2 inv'*B*T1 inv;% Equation (36)
6. [U,_,V]=svd(B1,0);
7. e1=V(:,3)./norm(V(1:,2,3));
8. e2=U(:,3)./norm(U(1:,2,3));
9. R1=[e1(1), e1(2), 0; -e1(2), e1(1), 0; 0, 0, 1]; % Equation (37)
10. R2=[e2(1), e2(2), 0; -e2(2), e2(1), 0; 0, 0, 1]; % Equation (37)
11.