ある３次元ベクトルに直交するベクトルの算出方法

usao · #1

C++で，「ある３次元ベクトルに直行するベクトルを求める処理」を実装したいです．
そして，自分で考えて実装したのが，下記コードの　GetPerpendicularOf() 関数です．

この関数の実装に関して，
(1)誤った答えを出してしまうような入力が有り得るでしょうか？
(2)（何らかの意味で）より良い実装方法（アルゴリズム?）を御存じの方がおられましたら，ご教示願いたいです．

コード:

#include <iostream>
#include <limits>
#include <random>

//3次元ベクトル
struct Vec3
{
	Vec3( double x=0, double y=0, double z=0 ){	V[0]=x;	V[1]=y;	V[2]=z;	}
	double &operator[]( int i ){	return V[i];	}
	double operator[]( int i ) const {	return V[i];	}
	double V[3];	//要素
};

//3次元ベクトルVに直行する3次元ベクトルを返す．
//(無限に有り得るうちのどれか一つ．もちろん，ゼロベクトルではない．)
//※ただし，Vがゼロベクトルの時は，戻り値は何らかの単位ベクトルを返すこととする．
Vec3 GetPerpendicularOf( const Vec3 &V )
{
	//Vの3要素の絶対値を比較して{最大,真ん中,最小}の要素のindex {iMax,iMddle,iMin}を求む
	const double Abs[3] = {	fabs(V[0]), fabs(V[1]), fabs(V[2])	};
	int iMax = 0;
	if( Abs[iMax] < Abs[1] ){	iMax=1;	}
	if( Abs[iMax] < Abs[2] ){	iMax=2;	}

	int iMiddle = (iMax+1)%3;
	int iMin = (iMax+2)%3;
	if( Abs[iMiddle] < Abs[iMin] ){	std::swap( iMiddle, iMin );	}

	//Vに直行するベクトルRを作る
	Vec3 R;	//(0,0,0)
	if( Abs[iMiddle] <= std::numeric_limits<double>::epsilon() )
	{//Vの要素のうち，2つ以上が0と見なせるとき
		R[iMin] = 1;
	}
	else
	{//少なくともV[iMiddle]が0ではない値であるとき
		R[iMiddle] = V[iMin];
		R[iMin] = -V[iMiddle];
	}
	return R;
}

//
inline std::ostream &operator <<( std::ostream &OS, const Vec3 &V )
{	return ( OS << "[ " << V[0] << ", " << V[1] << ", " << V[2] << " ]" );	}

//
int main()
{
	//テストコード
	std::random_device seed_gen;
	std::mt19937 engine( seed_gen() );
	std::uniform_real_distribution<> dist( -10000, 10000 );

	for ( int i=0; i<32; ++i )
	{
		Vec3 V( dist(engine), dist(engine), dist(engine) );
		Vec3 R = GetPerpendicularOf( V );
		double IP = V[0]*R[0] + V[1]*R[1] + V[2]*R[2];	//innner-product

		std::cout << "V = " << V << std::endl;
		std::cout << "R = " << R << std::endl;
		std::cout << "IP = " << IP << std::endl;
		std::cout << std::endl;
	}

	//
	Vec3 R = GetPerpendicularOf( Vec3() );
	std::cout << "ResultOf(0,0,0) = " << R << std::endl;

	//
	std::cin.ignore();
	return 0;
}

[追加情報]
・何らかの課題ではありません．

以上，宜しくお願い申し上げます．

usao · #2

オフトピック

おおっと！
各所で　「直交」であるべき箇所が「直行」と誤変換されている！

かずま · #3

usao さんが書きました:(1)誤った答えを出してしまうような入力が有り得るでしょうか？

コード:

void test(Vec3 &V)
{
    Vec3 R = GetPerpendicularOf(V);
    double IP = V[0]*R[0] + V[1]*R[1] + V[2]*R[2];  //innner-product
    std::cout << "V = " << V << std::endl;
    std::cout << "R = " << R << std::endl;
    std::cout << "IP = " << IP << std::endl;
    std::cout << std::endl;
}

int main()
{
    test(Vec3(2e-16, 2.3e-16, 2.3e-16));
    test(Vec3(2e-16, 2.2e-16, 2.3e-16));
    test(Vec3(2.2e-16, 2e-16, 2.3e-16));
    test(Vec3(2e-16, 2.2e-16, 2.3e-16));
}

実行結果

コード:

V = [ 2e-016, 2.3e-016, 2.3e-016 ]
R = [ -2.3e-016, 0, 2e-016 ]
IP = 0

V = [ 2e-016, 2.2e-016, 2.3e-016 ]
R = [ 1, 0, 0 ]
IP = 2e-016

V = [ 2.2e-016, 2e-016, 2.3e-016 ]
R = [ 0, 1, 0 ]
IP = 2e-016

V = [ 2e-016, 2.2e-016, 2.3e-016 ]
R = [ 1, 0, 0 ]
IP = 2e-016

2e-16 は std::numeric_limits<double>::epsilon() より
小さいので、IP はゼロになったとも言えますが。

かずま · #4

すみません。4番目が 2番目と同じでした。
4番目は、次のつもりでした。

コード:

    test(Vec3(2.3e-16, 2.2e-16, 2e-16));

実行結果

コード:

V = [ 2.3e-016, 2.2e-016, 2e-016 ]
R = [ 0, 0, 1 ]
IP = 2e-016

usao · #5

ありがとうございます．

なるほど…　関数が

コード:

if( Abs[iMiddle] <= std::numeric_limits<double>::epsilon() )

として，勝手な方法で判断をしているせいで，関数使用側にその精度を押しつけてしまうということですね．
関数内は単純に

コード:

if( Abs[iMiddle] == 0 )  //本当に0のときだけは考える

とした方が良さそうですね．

かずま · #6

usao さんが書きました:(2)（何らかの意味で）より良い実装方法（アルゴリズム?）を御存じの方がおられましたら，ご教示願いたいです．

コード:

    int iMiddle = (iMax+1)%3;
    int iMin = (iMax+2)%3;

より

コード:

    int iMiddle = !iMax; // または iMiddle = (2 - iMax) >> 1;
    int iMin = 3 - iMax - iMiddle;

のほうが、割り算しないから速いと思うんですが、

コード:

int get();
void func(int, int);

void f1()
{
    int iMax = get();
    int iMiddle = (iMax + 1) % 3;
    int iMin = (iMax + 2) % 3;
    func(iMiddle, iMin);
}

void f2()
{
    int iMax = get();
    int iMiddle = !iMax;
    int iMin = 3 - iMiddle - iMax;
    func(iMiddle, iMin);
}

void f3()
{
    int iMax = get();
    int iMiddle = (2 - iMax) >> 1;
    int iMin = 3 - iMiddle - iMax;
    func(iMiddle, iMin);
}

gcc -O2 -c hoge.c; objdump -d hoge.o の結果の抜粋

コード:

0000000000000000 <f1>:
   0:   48 83 ec 28             sub    $0x28,%rsp
   4:   e8 00 00 00 00          callq  9 <f1+0x9>
   9:   44 8d 48 02             lea    0x2(%rax),%r9d
   d:   41 b8 56 55 55 55       mov    $0x55555556,%r8d
  13:   89 c1                   mov    %eax,%ecx
  15:   83 c1 01                add    $0x1,%ecx
  18:   44 89 c8                mov    %r9d,%eax
  1b:   41 f7 e8                imul   %r8d
  1e:   44 89 c8                mov    %r9d,%eax
  21:   c1 f8 1f                sar    $0x1f,%eax
  24:   29 c2                   sub    %eax,%edx
  26:   8d 04 52                lea    (%rdx,%rdx,2),%eax
  29:   41 29 c1                sub    %eax,%r9d
  2c:   89 c8                   mov    %ecx,%eax
  2e:   41 f7 e8                imul   %r8d
  31:   89 c8                   mov    %ecx,%eax
  33:   c1 f8 1f                sar    $0x1f,%eax
  36:   29 c2                   sub    %eax,%edx
  38:   8d 04 52                lea    (%rdx,%rdx,2),%eax
  3b:   44 89 ca                mov    %r9d,%edx
  3e:   29 c1                   sub    %eax,%ecx
  40:   48 83 c4 28             add    $0x28,%rsp
  44:   e9 00 00 00 00          jmpq   49 <f1+0x49>

0000000000000050 <f2>:
  50:   48 83 ec 28             sub    $0x28,%rsp
  54:   e8 00 00 00 00          callq  59 <f2+0x9>
  59:   31 c9                   xor    %ecx,%ecx
  5b:   85 c0                   test   %eax,%eax
  5d:   ba 03 00 00 00          mov    $0x3,%edx
  62:   0f 94 c1                sete   %cl
  65:   29 ca                   sub    %ecx,%edx
  67:   29 c2                   sub    %eax,%edx
  69:   48 83 c4 28             add    $0x28,%rsp
  6d:   e9 00 00 00 00          jmpq   72 <f2+0x22>

0000000000000080 <f3>:
  80:   48 83 ec 28             sub    $0x28,%rsp
  84:   e8 00 00 00 00          callq  89 <f3+0x9>
  89:   b9 02 00 00 00          mov    $0x2,%ecx
  8e:   ba 03 00 00 00          mov    $0x3,%edx
  93:   29 c1                   sub    %eax,%ecx
  95:   d1 f9                   sar    %ecx
  97:   29 ca                   sub    %ecx,%edx
  99:   29 c2                   sub    %eax,%edx
  9b:   48 83 c4 28             add    $0x28,%rsp
  9f:   e9 00 00 00 00          jmpq   a4 <f3+0x24>

VC++ で、cl -O2 -c hoge.c; dumpbin -disasm hoge.obj の結果の抜粋

コード:

_f1:
  00000000: 56                 push        esi
  00000001: E8 00 00 00 00     call        _get
  00000006: 8B C8              mov         ecx,eax
  00000008: BE 03 00 00 00     mov         esi,3
  0000000D: 8D 41 02           lea         eax,[ecx+2]
  00000010: 99                 cdq
  00000011: F7 FE              idiv        eax,esi
  00000013: 8D 41 01           lea         eax,[ecx+1]
  00000016: 52                 push        edx
  00000017: 99                 cdq
  00000018: F7 FE              idiv        eax,esi
  0000001A: 52                 push        edx
  0000001B: E8 00 00 00 00     call        _func
  00000020: 83 C4 08           add         esp,8
  00000023: 5E                 pop         esi
  00000024: C3                 ret

_f2:
  00000000: E8 00 00 00 00     call        _get
  00000005: 33 D2              xor         edx,edx
  00000007: B9 03 00 00 00     mov         ecx,3
  0000000C: 85 C0              test        eax,eax
  0000000E: 0F 94 C2           sete        dl
  00000011: 2B CA              sub         ecx,edx
  00000013: 2B C8              sub         ecx,eax
  00000015: 51                 push        ecx
  00000016: 52                 push        edx
  00000017: E8 00 00 00 00     call        _func
  0000001C: 83 C4 08           add         esp,8
  0000001F: C3                 ret

_f3:
  00000000: E8 00 00 00 00     call        _get
  00000005: BA 02 00 00 00     mov         edx,2
  0000000A: B9 03 00 00 00     mov         ecx,3
  0000000F: 2B D0              sub         edx,eax
  00000011: D1 FA              sar         edx,1
  00000013: 2B CA              sub         ecx,edx
  00000015: 2B C8              sub         ecx,eax
  00000017: 51                 push        ecx
  00000018: 52                 push        edx
  00000019: E8 00 00 00 00     call        _func
  0000001E: 83 C4 08           add         esp,8
  00000021: C3                 ret

usao · #7

引き続きのご指摘，ありがとうございます．
割り算不要はその通りだと思いましたので，
自分なりにもっと直接的に（というか，素直に）indexを求めるように変えてみました．

コード:

Vec3 GetPerpendicularOf( const Vec3 &V )
{
	//Vの3要素のうち，絶対値が最大ではない2つの要素のindex {iNotMax1,iNotMax2}を求む
	const double Abs[3] = {	fabs(V[0]), fabs(V[1]), fabs(V[2])	};
	int iNotMax1 = 1;
	int iNotMax2 = 2;
	{
		int iMax = 0;
		if( Abs[0] < Abs[1] ){	iNotMax1 = 0;	iMax = 1;	}
		if( Abs[iMax] < Abs[2] ){	iNotMax2 = iMax;	}
	}

	//Vに直交するベクトルRを作る
	Vec3 R;	//(0,0,0)
	if( Abs[iNotMax1]==0  &&  Abs[iNotMax2]==0 )  //「==0」で，２つの要素が0であることを判定するようにした
	{
		R[iNotMax1] = 1;
	}
	else
	{
		R[iNotMax1] =  V[iNotMax2];
		R[iNotMax2] = -V[iNotMax1];
	}
	return R;
}

かずま · #8

今、飲み屋でこんなのを思いつきました。

コード:

if (V[0]==0 && V[1]==0 && V[2]==0) return Vec3(1);
if (V[0]==0) return Vec3(0, V[2], -V[1]);
if (V[1]==0) return Vec3(V[2], 0, -V[0]);
return Vec3(V[1], -V[0]);

PC が無いのでチェックできません。

かずま · #9

確認しました。

コード:

#include <iostream>

struct Vec3 {
    Vec3(double x=0, double y=0, double z=0) { V[0]=x; V[1]=y; V[2]=z; }
    double &operator[](int i) { return V[i]; }
    double operator[](int i) const { return V[i]; }
    double V[3];    //要素
};
 
Vec3 GetPerpendicularOf(const Vec3 &V)
{
    return V[0] ? (V[1] ? Vec3(V[1],-V[0]) : Vec3(V[2],0,-V[0]))
                : (V[1] || V[2] ? Vec3(0,V[2],-V[1]) : Vec3(1));
}

inline std::ostream &operator <<(std::ostream &OS, const Vec3 &V)
{
    return OS << "[ " << V[0] << ", " << V[1] << ", " << V[2] << " ]";
}
 
void test(const Vec3 &V)
{
    Vec3 R = GetPerpendicularOf(V);
    double IP = V[0]*R[0] + V[1]*R[1] + V[2]*R[2];  //innner-product
    std::cout << "V = " << V << std::endl;
    std::cout << "R = " << R << std::endl;
    std::cout << "IP = " << IP << std::endl;
    std::cout << std::endl;
}

int main()
{
    test(Vec3(0, 0, 0));
    test(Vec3(0, 0, 1));
    test(Vec3(0, 2, 0));
    test(Vec3(3, 0, 0));
    test(Vec3(4, 4, 0));
    test(Vec3(5, 0, 5));
    test(Vec3(0, 6, 6));
    test(Vec3(7, 7, 7));
}

実行結果

コード:

V = [ 0, 0, 0 ]
R = [ 1, 0, 0 ]
IP = 0

V = [ 0, 0, 1 ]
R = [ 0, 1, -0 ]
IP = 0

V = [ 0, 2, 0 ]
R = [ 0, 0, -2 ]
IP = 0

V = [ 3, 0, 0 ]
R = [ 0, 0, -3 ]
IP = 0

V = [ 4, 4, 0 ]
R = [ 4, -4, 0 ]
IP = 0

V = [ 5, 0, 5 ]
R = [ 5, 0, -5 ]
IP = 0

V = [ 0, 6, 6 ]
R = [ 0, 6, -6 ]
IP = 0

V = [ 7, 7, 7 ]
R = [ 7, -7, 0 ]
IP = 0

-0 がちょっとあれですが。

usao · #10

無い頭でいろいろ考えてた経緯の中で，要素の大小関係をごちゃごちゃやっていたのですが…そんなことをする必要はなくて，

【非ゼロな要素を１つ見つけてそれを別の要素とswapして一方に-1をかける，残りの1要素は0にする】
という操作さえ行えば完了
ということで良いのか……な？

てことは，

コード:

	return V[0] ? Vec3(V[1],-V[0])
		: (V[1] || V[2] ? Vec3(0,V[2],-V[1]) : Vec3(1));

でいいのかもしれませんね．

usao · #11

テストと考察とをしてみた感じ，
No.10のコードで目的の物が算出できていると思われるので，ここで締めさせていただきます([解決]を付す)．

かずまさん
数々のご指摘とアイデア頂き本当にありがとうございました．

ある３次元ベクトルに直交するベクトルの算出方法

ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

Re: ある３次元ベクトルに直交するベクトルの算出方法

		8月 2025
日	月	火	水	木	金	土
					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