任意の平面に投影した点（垂線の足）の座標の求め方

三次元座標上に存在する点を「法線ベクトルで規定される平面」に投影する方法を教えてください．

H = P - (n・AP)ｎ（ｎは単位ベクトル）を用いて計算していますが，上手くいきません．

m1～m9のマーカーの三次元座標データがあります．

load marker

マーカーの並びは以下の通りです．

% m1------m2------m3
% |  \    |    /  |
% |    \  |  /    |
% m4------m5------m6
% |    /  |  \    |
% |  /    |    \  |
% m7------m8------m9

m5を注目マーカーとします．

【行いたい操作】

ｍ5と周りの８個のマーカーによって作られる8個の平面（例．m5,m1,m2によって作られる平面）の「平均の面」にm1～m9を投影する．

問題点が分からないため，間違っている箇所を教えていただきたいです．

m5を中心とし，周りの各マーカーに向かう8個のベクトル　m5_1～m5_9　を作成する．

m5_1 = m1 - m5;         % m5からm1に向かうベクトル
m5_2 = m2 - m5;         % m5からm2に向かうベクトル
m5_3 = m3 - m5;         % m5からm3に向かうベクトル
m5_4 = m4 - m5;         % m5からm4に向かうベクトル
m5_6 = m6 - m5;         % m5からm6に向かうベクトル
m5_7 = m7 - m5;         % m5からm7に向かうベクトル
m5_8 = m8 - m5;         % m5からm8に向かうベクトル
m5_9 = m9 - m5;         % m5からm9に向かうベクトル

隣り合うベクトルの法線ベクトルにより，８個の面を規定する．

c1 = cross(m5_2,m5_1);      % m5_1とm5_2の法線ベクトルによって規定される面         m1-----------m2----------m3
c2 = cross(m5_3,m5_2);      % m5_2とm5_3の法線ベクトルによって規定される面          |  \   c1   |    c2   / |
c3 = cross(m5_6,m5_3);      % m5_3とm5_6の法線ベクトルによって規定される面          |    \      |      /    |
c4 = cross(m5_9,m5_6);      % m5_6とm5_9の法線ベクトルによって規定される面          |  c8   \   |   /   c3  |
c5 = cross(m5_8,m5_9);      % m5_9とm5_8の法線ベクトルによって規定される面         m4-----------m5----------m6
c6 = cross(m5_7,m5_8);      % m5_8とm5_7の法線ベクトルによって規定される面          |  c7   /   |   \   c4  |  
c7 = cross(m5_4,m5_7);      % m5_7とm5_4の法線ベクトルによって規定される面          |    /      |      \    | 
c8 = cross(m5_1,m5_4);      % m5_4とm5_1の法線ベクトルによって規定される面          |  /   c6   |    c5   \ |
                            %                                                   m7----------m8----------m9

c1～c8の単位ベクトルを算出する．

unit_c1 = c1./sqrt(sum(c1.^2));
unit_c2 = c2./sqrt(sum(c2.^2));
unit_c3 = c3./sqrt(sum(c3.^2));
unit_c4 = c4./sqrt(sum(c4.^2));
unit_c5 = c5./sqrt(sum(c5.^2));
unit_c6 = c6./sqrt(sum(c6.^2));
unit_c7 = c7./sqrt(sum(c7.^2));
unit_c8 = c8./sqrt(sum(c8.^2));

8個の面の平均面を規定する法線ベクトルを算出する．

ave_c = (unit_c1 + unit_c2 + unit_c3 + unit_c4 + unit_c5 + unit_c6 + unit_c7 + unit_c8)./8;     % 各面を規定する法線ベクトルの平均
unit_ave_c = ave_c./sqrt(sum(ave_c.^2));                                                        % 単位ベクトル化

unit_ave_cで規定される面に m1～m9を投影する．

p1 = m1 - (dot(m5_1,unit_ave_c).*unit_ave_c);       % m1の投影点
p2 = m2 - (dot(m5_2,unit_ave_c).*unit_ave_c);       % m2の投影点
p3 = m3 - (dot(m5_3,unit_ave_c).*unit_ave_c);       % m3の投影点
p4 = m4 - (dot(m5_4,unit_ave_c).*unit_ave_c);       % m4の投影点
p6 = m6 - (dot(m5_6,unit_ave_c).*unit_ave_c);       % m6の投影点
p7 = m7 - (dot(m5_7,unit_ave_c).*unit_ave_c);       % m7の投影点
p8 = m8 - (dot(m5_8,unit_ave_c).*unit_ave_c);       % m8の投影点
p9 = m9 - (dot(m5_9,unit_ave_c).*unit_ave_c);       % m9の投影点

p1～p9がunit_ave_cで規定される面に投影されている場合，「m5」から「p1～p9」に向かう８個のベクトルはすべてunit_ave_cと垂直になるはずですが，どうも上手くいきません．

m5_p1 = p1 - m5;        % m5からp1に向かうベクトル
m5_p2 = p2 - m5;        % m5からp2に向かうベクトル
m5_p3 = p3 - m5;        % m5からp3に向かうベクトル
m5_p4 = p4 - m5;        % m5からp4に向かうベクトル
m5_p6 = p6 - m5;        % m5からp6に向かうベクトル
m5_p7 = p7 - m5;        % m5からp7に向かうベクトル
m5_p8 = p8 - m5;        % m5からp8に向かうベクトル
m5_p9 = p9 - m5;        % m5からp9に向かうベクトル
s1 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p1のドット積

s1 = -1.7347e-17

s2 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p2のドット積

s2 = -1.7347e-17

s3 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p3のドット積

s3 = -1.7347e-17

s4 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p4のドット積

s4 = -1.7347e-17

s6 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p6のドット積

s6 = -1.7347e-17

s7 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p7のドット積

s7 = -1.7347e-17

s8 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p8のドット積

s8 = -1.7347e-17

s9 = dot(unit_ave_c,m5_p1)        % unit_ave_cとm5_p9のドット積

s9 = -1.7347e-17

[Keita Abe]

ドット積の値が 10^-17のオーダーなので、丸め誤差や桁落ちで厳密に0にならなかったのではないでしょうか？実用上はeps=10^-10などと決めて、ドット積の絶対値がepsを下回ることで垂直判定をしてはいかがでしょうか。

[知真 梶山]

ご回答ありがとうございます．

おっしゃる通り，丸め誤差で厳密に０にならないようです．

教えていただいた方法を使わせていただきます．