以下は、Multiple View Geometry で説明されている三角形分割法を実装する Matlab 関数です。立方体の 8 つのコーナー ポイント (OriginalPoints) を取り、それらを 2 つのスクリーンに投影します。スクリーン ポイントは、CameraPoints と ProjectorPoints にあります。目標は、直接線形変換を使用して、これら 2 つのビューから元のポイントを再構築することです。残念ながら、結果は意味不明です。
SO に関する別のトピックがあり、そこから A の行の正規化が行われますが、結果は改善されませんでした。
3 つの一般的な質問:
- SVD (原点 = 重心および平均距離 = sqrt(2)) の前にデータを正規化する場合、(x,y) または (x, y, w) の平均距離を計算しますか? 不均一または均一？
- すべての操作の後、可能な場合は同次座標 w を 1.0 に戻しますが、これは正しいですか?
- カメラ行列は (射影変換までではなく) 完全にわかっているため、結果のポイントは、ユークリッド変換と同じように元のポイントと変わらないはずですよね?
% Demos the triangulation method using the direct linear transform
% algorithm
function testTriangulation()
close all;
camWidth = 800;
camHeight = 600;
projectorWidth = 5080;
projectorHeight = 760;
% camera matrices
CameraMatrix = [
-1.81066 0.00000 0.00000 0.00000;
0.00000 2.41421 0.00000 0.00000;
0.00000 0.00000 1.00000 1.50000
];
ProjectorMatrix = [
-0.36118 0.00000 0.00000 0.00000
0.00000 2.00875 1.33916 0.00000
0.00000 -0.55470 0.83205 1.80278
];
% generating original points and displaying them
OriginalPoints = [
-0.2 -0.2 -0.2 -0.2 0.2 0.2 0.2 0.2;
-0.2 -0.2 0.2 0.2 -0.2 -0.2 0.2 0.2;
-0.2 0.2 -0.2 0.2 -0.2 0.2 -0.2 0.2;
1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
];
scatter3(OriginalPoints(1, :), OriginalPoints(2, :), OriginalPoints(3,:));
title('Original Points');
% projecting and normalizing camera points
CameraPoints = CameraMatrix * OriginalPoints;
for i = 1:3
CameraPoints(i, :) = CameraPoints(i, :) ./ CameraPoints(3,:);
end
% figure();
% scatter(CameraPoints(1,:), CameraPoints(2,:));
% title('Projected Camera Points');
% transforming camera points to screen coordinates
camXOffset = camWidth / 2;
camYOffset = camHeight / 2;
CameraPoints(1,:) = CameraPoints(1,:) * camXOffset + camXOffset;
CameraPoints(2,:) = CameraPoints(2,:) * camYOffset + camYOffset;
figure();
scatter(CameraPoints(1,:), CameraPoints(2,:));
title('Projected Camera Points in Screen Coordinates');
% projecting and normalizing projector points
ProjectorPoints = ProjectorMatrix * OriginalPoints;
for i = 1:3
ProjectorPoints(i, :) = ProjectorPoints(i, :) ./ ProjectorPoints(3,:);
end
% figure();
% scatter(ProjectorPoints(1,:), ProjectorPoints(2,:));
% title('Projected Projector Points');
% transforming projector points to screen coordinates
projectorXOffset = projectorWidth / 2;
projectorYOffset = projectorHeight / 2;
ProjectorPoints(1,:) = ProjectorPoints(1,:) * projectorXOffset + projectorXOffset;
ProjectorPoints(2,:) = ProjectorPoints(2,:) * projectorYOffset + projectorYOffset;
figure();
scatter(ProjectorPoints(1,:), ProjectorPoints(2,:));
title('Projected Projector Points in Screen Coordinates');
% normalizing camera points (i.e. origin is
% the centroid and average distance is sqrt(2)).
camCenter = mean(CameraPoints, 2);
CameraPoints(1,:) = CameraPoints(1,:) - camCenter(1);
CameraPoints(2,:) = CameraPoints(2,:) - camCenter(2);
avgDistance = 0.0;
for i = 1:length(CameraPoints(1,:))
avgDistance = avgDistance + norm(CameraPoints(:, i));
end
avgDistance = avgDistance / length(CameraPoints(1, :));
scaleFactor = sqrt(2) / avgDistance;
CameraPoints(1:2, :) = CameraPoints(1:2, :) * scaleFactor;
% normalizing projector points (i.e. origin is
% the centroid and average distance is sqrt(2)).
projectorCenter = mean(ProjectorPoints, 2);
ProjectorPoints(1,:) = ProjectorPoints(1,:) - projectorCenter(1);
ProjectorPoints(2,:) = ProjectorPoints(2,:) - projectorCenter(2);
avgDistance = 0.0;
for i = 1:length(ProjectorPoints(1,:))
avgDistance = avgDistance + norm(ProjectorPoints(:, i));
end
avgDistance = avgDistance / length(ProjectorPoints(1, :));
scaleFactor = sqrt(2) / avgDistance;
ProjectorPoints(1:2, :) = ProjectorPoints(1:2, :) * scaleFactor;
% triangulating points
TriangulatedPoints = zeros(size(OriginalPoints));
for i = 1:length(OriginalPoints(1, :))
A = [
CameraPoints(1, i) * CameraMatrix(3, :) - CameraMatrix(1, :);
CameraPoints(2, i) * CameraMatrix(3, :) - CameraMatrix(2, :);
ProjectorPoints(1, i) * ProjectorMatrix(3,:) - ProjectorMatrix(1,:);
ProjectorPoints(2, i) * ProjectorMatrix(3,:) - ProjectorMatrix(2,:)
];
% normalizing rows of A, as suggested in a topic on StackOverflow
A(1,:) = A(1,:)/norm(A(1,:));
A(2,:) = A(2,:)/norm(A(2,:));
A(3,:) = A(3,:)/norm(A(3,:));
A(4,:) = A(4,:)/norm(A(4,:));
[U S V] = svd(A);
TriangulatedPoints(:, i) = V(:, end);
end
for i = 1:4
TriangulatedPoints(i, :) = TriangulatedPoints(i, :) ./ TriangulatedPoints(4, :);
end
figure()
scatter3(TriangulatedPoints(1, :), TriangulatedPoints(2, :), TriangulatedPoints(3, :));
title('Triangulated Points');