%camera center: C = [8;-8;8;1] %rotation: %around x-axis, Pi/6: Rx = [[1;0;0], [0; cos(pi/6); sin(pi/6)],[0; -sin(pi/6); cos(pi/6)]]; %around y-axis, Pi/4: Ry = [[cos(pi/4); 0; -sin(pi/4)],[0;1;0],[sin(pi/4);0;cos(pi/4)]]; %overall rotation: R = Rx * Ry %intrinsic params (f = 3, mx=my=1, px=py=0): K = [[3;0;0],[0;3;0],[0;0;1]] %now, make projection matrix: P = K * [R,(-R*C(1:3))] X = [[0,0,0,0,-4,-4,-4,-4];[0,0,4,4,0,4,4,0];[0,4,0,4,4,4,0,0];[1,1,1,1,1,1,1,1]] %The corresponding 2D points on the image plane: x=P*X trials = 10000; %Exercise e,f,g eMeasure = zeros(1,8); eEstimation = zeros(1,8); eResidual = zeros(1,8); Pcov = zeros(12, trials); for k=1:trials %copy x to xn (x with noise) xn = x; % Exercise a % add normally distributed noise to x,y but not w xn(1:2,1:8) = xn(1:2,1:8)+0.07*randn(2,8); % Exercise b [Pnoise, r] = DLTforCamCalib(xn',X'); Pcov(1,k) = P(1,1); Pcov(2,k) = P(1,2); Pcov(3,k) = P(1,3); Pcov(4,k) = P(1,4); Pcov(5,k) = P(2,1); Pcov(6,k) = P(2,2); Pcov(7,k) = P(2,3); Pcov(8,k) = P(2,4); Pcov(9,k) = P(3,1); Pcov(10,k) = P(3,2); Pcov(11,k) = P(3,3); Pcov(12,k) = P(3,4); %Exercise c xe=Pnoise*X; hold on; for i=1:8 xn(:,i)=xn(:,i)/xn(3,i); x(:,i)=x(:,i)/x(3,i); xe(:,i)=xe(:,i)/xe(3,i); end %Exercise e,f,g eMeasure = eMeasure + ((xn(1,:)-x(1,:)) .* (xn(1,:)-x(1,:)) + (xn(2,:)-x(2,:)) .* (xn(2,:)-x(2,:))); eEstimation = eEstimation + ((x(1,:)-xe(1,:)) .* (x(1,:)-xe(1,:)) + (x(2,:)-xe(2,:)) .* (x(2,:)-xe(2,:))); eResidual = eResidual + ((xn(1,:)-xe(1,:)) .* (xn(1,:)-xe(1,:)) + (xn(2,:)-xe(2,:)) .* (xn(2,:)-xe(2,:))); %Exercise d plot3(xn(1,1:8), xn(2,1:8),f*ones(1,8),'.r'); %plot3(x(1,1:8), x(2,1:8),f*ones(1,8),'ob'); %plot3(xe(1,1:8), xe(2,1:8),f*ones(1,8),'.g'); end %Exercise e,f,g eMeasure = eMeasure ./ trials; eEstimation = eEstimation ./ trials; eResidual = eResidual ./ trials; %determine mean and covariance: meanP = zeros(12, 1); for i=1:12 for j=1:trials meanP(i,1) = meanP(i,1) + Pcov(i,j); end meanP(i,1) = meanP(i,1) / trials; end %(same as mean(covP');) covP = zeros(12,12); for i=1:12 for j=1:12 for k=1:trials covP(i,j) = covP(i,j) + (Pcov(i,k) - meanP(i,1)) * (Pcov(j,k) - meanP(j,1)); end covP(i,j) = covP(i,j) / trials; end end %(same as cov(covP');) %plotLine(xe(1:4,1),xe(1:4,2)) %plotLine(xe(1:4,1),xe(1:4,3)) %plotLine(xe(1:4,2),xe(1:4,4)) %plotLine(xe(1:4,4),xe(1:4,3)) %plotLine(xe(1:4,4),xe(1:4,6)) %plotLine(xe(1:4,6),xe(1:4,7)) %plotLine(xe(1:4,7),xe(1:4,3)) %plotLine(xe(1:4,6),xe(1:4,5)) %plotLine(xe(1:4,5),xe(1:4,8)) %plotLine(xe(1:4,8),xe(1:4,7)) %plotLine(xe(1:4,8),xe(1:4,1)) %plotLine(xe(1:4,5),xe(1:4,2)) %Exercise e for i=1:8 eMeasure = eMeasure + ((xn(1,i)-xe(1,i))^2 + (xn(2,i)-xe(2,i))^2); end eMeasure = eMeasure / 8;