３点の座標から簡単に角度と回転方向を求める．
(2・3・N次元，外積を用いる方法)

2010/09/13(月)

３点の座標 (２次元または３次元) から三角形の辺の長さ・角度， ベクトルの内積・外積などを自動的に計算する Excel ファイルを， DLmarket 様にて委託ダウンロード販売開始しました．

1. 問題 (２次元の場合)
2. 解答
3. 解説
4. 注意
5. 余談 (Ｎ次元の外積について)
6. ２つのＮ次元ベクトルのなす角 (回転角) θと回転方向を求める．
   1. 角度だけを求める．(Ｎ次元，回転方向の情報なし)
   2. 角度＋回転方向を求める．(２次元)
   3. 角度＋回転方向を求める．(３次元)
   4. 角度＋回転方向を求める．(Ｎ次元)
7. ３点の座標から角度・距離・内積・外積などを計算する Excel ファイル
   • Triangle2D.xls (２次元用)
   • Triangle3D.xls (３次元用)
8. サンプルプログラム (Ｃ言語)
   1. 外積を計算する関数・マクロ
   2. ３点の座標から角度と回転方向を求める．(２次元)
9. 参考・関連図書
10. サイト内関連ページ
11. 外部へのリンク
12. 更新履歴

• ３次元の回転 (原点を通る任意方向回転軸，座標系に依存しないベクトル・回転行列表現) についてはこちらへ．
• 平面内のＮ個の座標列の回転方向を判別したい場合はこちらへ．

下図において，線分 CP が点 C を中心に角度θ (-180°＜ θ ＜ +180°) だけ回転した結果， 線分 CQ になったとする．３つの点の座標 C＝(Cx, Cy), P＝(Px, Py), Q＝(Qx, Qy) が与えられたとき，この回転が右回りか左回りかを判別せよ．
(別の言い方をすれば，「C → P → Q → C は右回りか左回りか？」)

Ｙ
↑
│          Ｑ
│        ／
│      ／θ
│    Ｃ────Ｐ
│
Ｏ────────→Ｘ

S �� (Px - Cx) * (Qy - Cy) - (Py - Cy) * (Qx - Cx)

とする．S＞0 なら左回り，S＜0 なら右回り，S＝0 ならば C，P，Q は一直線上にある．(注)

なお，この判別方法は，CP と CQ が同じ長さである必要はない．

θを求めたい場合はこちらへ．

この問題を見て，逆三角関数 tan^-1 (Ｃ言語では atan() や atan2()) を使って CP と CQ の角度をそれぞれ求め， 両者を比較しようと考えた方が多いのではないでしょうか． しかしこの問題では，角度そのものではなく角度差の符号を求めればよいので， 逆三角関数を使う方法よりも簡単で優れた，外積を使う方法を紹介します．

２つの２次元ベクトル A＝(Ax, Ay), B＝(Bx, By) の外積を次のように定義する．

A × B ≡ Ax * By − Ay * Bx

ここで OA から OB への回転角をθ (-180°＜ θ ＜ +180°) とする (下図)． ただし OA から OB に回転する方向が左回りならば θ＞０ とする． このとき，A×B ＝ |OA| * |OB| * sinθ である．したがって，

• A × B の符号はθの符号と同じ．正ならば左回り，負ならば右回り．
• |A × B| は，OA と OB を２辺とする平行四辺形の面積．(＝ △OAB の面積の２倍)

Ｙ
↑
│   Ｂ
│   / θ Ａ
│  /  ／
│ / ／
│/／
Ｏ───────→Ｘ

このように外積を使えば，２回の乗算だけですむので，(遅い) 逆三角関数サブルーチンを使う方法に比べて次のようにいいことずくめです．

• 高速．
• 計算誤差が少ない．
• 座標が整数ならば，外積も整数演算ですむ．その場合計算誤差は０．(ただしオーバーフローに注意)
• 逆三角関数では角度に不連続点があるので， ２つの角度を比較しようとすると場合分けが必要で，処理が複雑になる． また，角度のわずかな計算誤差が場合分けや判定結果に影響する場合もある． 外積を用いると解答に示したように場合分けは不要で，単一の数式を計算するだけでよい．

また，２つのベクトルの回転角 (上記のθ) を求める場合には， 後述するように内積と外積を併用すれば簡単です． この場合も，２つのベクトルそれぞれの角度を求め， それらの差を計算しようとすると場合分けが必要になり厄介です．

なお，外積は３次元で使われることが多く， 「外積」で検索してもヒットするのは多分ほとんどが３次元の話だと思います． ３次元ベクトル同士の外積は３次元ベクトルになりますが， ２次元ベクトル同士の外積はスカラーになります．
(３次元外積において，元の２つのベクトルがＸＹ平面内にある場合を考えれば， 外積のＸ，Ｙ成分は０になる．そのためＺ成分だけを考えればよく， それを２次元の外積と考えることもできる．)

上の解答および解説で，「外積の値が正ならば左回り」と書きましたが， それはＸ軸およびＹ軸の方向が上の図のようになっている場合です． もしＹ軸が下向き (あるいはＸ軸が左向き) になっていれば回転の向きが逆になります．(注の注)

そこで一般的な言い方をすると，次のようになります．

OX が90°回転して OY になる方向を正回転とすると， 外積が正ならば正回転，負ならば逆回転．

上の解説を読んで， 次のような疑問を持たれた方も多いでしょう．

あるいは，３次元外積だけを知っている人なら，次のように思われたかもしれません．

(a) のようになる理由は， ３次元の外積の定義を見るとわかりますが， 外積の１つの成分を計算するのに， 元の２つのベクトルの２つの座標軸の成分を使っていることにあります． 例えば３次元外積のＸ成分を計算するには， 元の２つのベクトルのＹおよびＺ成分を使います． ２つの座標軸の組合せごとに，外積の独立成分が１つ対応します． 外積の各独立成分は，Ｎ次元空間内の回転を， ２つの座標軸が張る平面に投影したものに対応していると考えられます．

したがってＮ次元ベクトル同士の外積には， _NC₂ 個の独立な成分があります． ３次元ではたまたま ₃C₂＝3 なのでベクトルとして表現できますし， ２次元ではたまたま ₂C₂＝1 なのでスカラーとして表現できます．

一般のＮ次元では，ベクトル A＝(A₁, A₂, …, A_N) と B＝(B₁, B₂, …, B_N) の外積は，次のような N×N 行列で表現できます．
((楔形の記号∧から) ウェッジ積，(考案者名から) グラスマン積ともいう．)

実はこれが本来の外積の定義であり， ベクトルでもスカラーでもありません． 上の説明では「行列で表現」と書きましたが， 正確に言うと外積は２階のテンソルと呼ばれる量です． Ｎ階のテンソルとは，大雑把に言えば，「Ｎ次元の行列」のようなものと考えてください． ０階のテンソルはスカラー，１階のテンソルはベクトルです．

さらに正確に言うと，外積は２階の 反対称 擬テンソルです． 「反対称」というのは，上に書いた行列の行と列を入れ替えると (つまり転置行列にすると)，符号が反転するということです． 「擬」とは，「右と左を入れ替える＝座標系の向きを反転させる (奇数本の座標軸を逆向きにする，あるいは２本の座標軸を交換する) と， 向きが反転する」ことを意味します．

上の「注意」で， 座標系の向きによって外積の符号に対応する回転の向きが変わるという意味のことを書きましたが， これは２次元外積が擬スカラーだからです． 同様に３次元外積は擬ベクトル (軸性ベクトルともいう) です．

軸性ベクトルに関する詳しい説明．

• ベクトルの外積　軸性ベクトルについて (OKWave)
• NO.1825　線形な（正直な？）エッセイ（？）其の１　2010.1.4.　DDT
  (Weekend Mathematics コロキウム室 NO.251)
• 有限要素法よもやま話 第86話 ベクトル余話 (FEMINGWAY 数理エッセイ)
• 有限要素法よもやま話 第87話 ベクトル余話・続き (FEMINGWAY 数理エッセイ)

高校数学Ｂでは２〜３次元ベクトルしか習わないが，何次元だろうと (無限次元だろうと)，ベクトル A＝(Ax，Ay，Az，…) と B＝(Bx，By，Bz，…) の内積の定義は，

幾何学的定義：A・B ≡ |A| * |B| * cosθ
代数的定義  ：A・B ≡ Ax * Bx + Ay * By + Az * Bz + …

したがって A と B のなす角θは，

           A・B
cosθ ＝ ────
          |A||B|

                    Ax * Bx + Ay * By + Az * Bz + …
      ＝ ─────────────────────────
          √((Ax2 + Ay2 + Az2 + …)(Bx2 + By2 + Bz2 + …))

あとは逆三角関数 cos^-1 を使うだけでθ (0°≦ θ ≦ 180°) が求まる． ２次元の場合は Az および Bz 以後を０とすればよい．

内積を使うと，何次元のなす角でも求められるが，回転方向の情報は得られない．

謎の検索ワード集

２次元ベクトル A と B の内積および外積の定義を改めて書くと，

幾何学的定義：A ・ B ≡ |A| * |B| * cosθ
　　　　　　　A × B ≡ |A| * |B| * sinθ

代数的定義  ：A ・ B ≡ Ax * Bx + Ay * By
              A × B ≡ Ax * By - Ay * Bx

なので，両者を併用すると簡単に角度および回転方向 (-180°＜ θ ≦ +180°) が求まる． C/C++ の場合は atan2(y, x) 関数，C# の場合は Math.Atan2(y, x) を使うと，

θ ＝ atan2(A×B，A・B) (単位はラジアン，引数の順序にも注意)

Excel の ATAN2(x, y) 関数はＣの atan2(y, x) とは引数の順序が逆なので，

θ ＝ ATAN2(A・B，A×B)

■参考

• Manpage of ATAN2
  C/C++ 標準ライブラリの atan2() 関数のマニュアル．
• Math.Atan2 メソッド (Double, Double)
  C# の Math.Atan2(Double y, Double x)．
• atan2 (Wikipedia 英語版)
  atan2 関数の数学的説明，図解．

３次元ベクトル A と B の外積 (ベクトル積，クロス積ともいう) A×B の定義は，
大学１年レベルの力学，電磁気学，物理数学の教科書にデカデカと書いてあるはず (ただしほとんど右手系限定)．
最近では，3DCG やゲームプログラミングの本 (の一部？) でも解説されている．

「外積 求め方」などで検索して来る人が多いけど， 教科書読んだことがないんだったら何で「外積」という言葉を知ってるの？ (激謎)

幾何学的定義：|A × B| ≡ |A| * |B| * sinθ    (0°≦ θ ≦ 180°)
              回転に対する A × B の向きは，座標系が右手系ならば右ネジの進む方向，
              左手系ならば左ネジの進む方向．(A：親指方向，B：人差し指方向，A×B：中指方向)

代数的定義  ：A × B ≡ (Ay * Bz - Az * By, Az * Bx - Ax * Bz, Ax * By - Ay * Bx)

θを求める方法として，２次元の場合と同様に外積と内積を併用することも可能だが， ３次元では計算式が複雑になるだけでメリットが (たぶん) ないので， 6.1 の方法を使う方がいいと思う．

回転方向については，A×B 自体がその方向を表しているのでそのまま使ってもよいが， 単位ベクトルにする必要がある場合には正規化して A×B / |A×B| とすればよい．

注意：４次元以上の回転についてちゃんと考えたわけではないので， この節に書いてあることを読む前に，眉に唾をベットリとつけてほしい．(笑)

まず回転角 (なす角) は，何次元であろうが 6.1 の方法で求めることができる．(これは間違いない．)

「余談 (Ｎ次元の外積について)」 に書いたことと２〜３次元の回転から類推すれば， Ｎ次元でも次式が成立するはずである．

|A ∧ B| ＝ |A| * |B| * sinθ    (0°≦ θ ≦ 180°)

ただしＮ次元外積テンソル A∧B の「長さ」は次のように定義する．

|A ∧ B| ≡ √ (Σ (Ai * Bj - Aj * Bi)2)
               i＜j

         ＝ √ ((1/2) ΣΣ (Ai * Bj - Aj * Bi)2)
                       i j

３次元と同様，A∧B 自体がその「回転方向」を表しているが， これを正規化すれば A∧B / |A∧B| となる．

３点の座標 (２次元または３次元) を入力するだけで，自動的に次の値を計算する Excel ファイルを，DLmarket 様にて委託ダウンロード販売しています．

• 三角形の辺の長さ，内角，面積，重心 (図心) 位置．(２次元，３次元)
• ２つのベクトルの内積，外積，長さ，なす角，回転方向．(２次元，３次元)
• 水平面 (XY 平面) に対する三角形 (を含む平面) の傾き (角度，方向)．(３次元のみ)
• 三角形の法線ベクトル (３点を含む面に垂直なベクトル)．(３次元のみ)

/*─────────────────────────────────────
機能  ：２つのベクトルの外積成分を計算する．
入力  ：(1) vector1，vector2：ベクトル (構造体のアドレスまたは配列)．
        (2) component1，component2：ベクトルの座標成分．ベクトルが構造体の
            場合はメンバ名，配列の場合は成分番号 (添字)．
戻り値：外積 vector1∧vector2 の (component1，component2) 成分．
2009/03/13(金) 作成
─────────────────────────────────────*/
// ベクトルが構造体の場合
#define CrossProductComponentS(vector1, vector2, component1, component2) \
  ((vector1)->component1 * (vector2)->component2 - \
   (vector1)->component2 * (vector2)->component1)

// ベクトルが配列の場合
#define CrossProductComponentA(vector1, vector2, component1, component2) \
  ((vector1)[component1] * (vector2)[component2] - \
   (vector1)[component2] * (vector2)[component1])

/*─────────────────────────────────────
機能  ：２つのｎ次元ベクトルのウェッジ積 (外積テンソル)を計算する．
        無駄な計算を避けるため，独立成分のみ求める．
入力  ：(1) vector1[0 〜 n-1]，vector2[0 〜 n-1]：ｎ次元ベクトル (配列)．
        (2) n (≧2)：ベクトルの次元．
出力  ：wedgeProduct[0 〜 nC2-1]：ウェッジ積 vector1∧vector2 の独立成分の配列．
注意  ：n＝3 の場合，wedgeProduct[] の要素の順序および符号は，通常の３次元
        外積の定義とは異なる．だから「３次元外積のプログラムを作成せよ」と
        いう学校の課題でこれを丸写ししてもダメよ．(笑)
2009/03/13(金) 作成
2013/01/14(月) 改名 (VectorN_CrossProduct → VectorN_WedgeProduct)
─────────────────────────────────────*/
typedef double VectorComponent_t; // ベクトルの成分の型

void VectorN_WedgeProduct(VectorComponent_t wedgeProduct[/* nC2 */],
                             const VectorComponent_t vector1[/* n */],
                             const VectorComponent_t vector2[/* n */],
                             unsigned n)
{
  unsigned i, j;

  for(i = 0;  i < n;  i++) {
    for(j = i;  ++j < n;  ) {
      *wedgeProduct++ = (VectorComponent_t)CrossProductComponentA(vector1, vector2, i, j);
    }
  }
}

#ifdef _MSC_VER
#define _USE_MATH_DEFINES // Visual C/C++ では，M_PI を使うのにこれが必要．
#endif /* _MSC_VER */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

// ラジアンを度に変換する．
#define RadianToDegree(radian)    ((180 / M_PI) * (radian))

// 度をラジアンに変換する．
#define DegreeToRadian(degree)    ((M_PI / 180) * (degree))

// ベクトルの成分の型
typedef double VectorComponent_t;

// ２次元ベクトル
typedef struct {
  VectorComponent_t x, y;
} Vector2_t;

// diff ← ベクトル v1 - v2
void Vector2_Diff(Vector2_t *diff, const Vector2_t *v1, const Vector2_t *v2)
{
  diff->x = v1->x - v2->x;
  diff->y = v1->y - v2->y;
}

// ベクトル v1 と v2 の内積 (v1・v2) を返す．
VectorComponent_t Vector2_DotProduct(const Vector2_t *v1, const Vector2_t *v2)
{
  return v1->x * v2->x + v1->y * v2->y;
}

// ベクトル v1 と v2 の外積 (v1×v2) を返す．
VectorComponent_t Vector2_CrossProduct(const Vector2_t *v1, const Vector2_t *v2)
{
  return v1->x * v2->y - v1->y * v2->x;
  // return CrossProductComponentS(v1, v2, x, y); でもよい．
}

// ベクトル C→P と C→Q のなす角θおよび回転方向を求める．
int main(void)
{
  Vector2_t c, p, q;   // 入力データ
  Vector2_t cp;        // ベクトル C→P
  Vector2_t cq;        // ベクトル C→Q
  VectorComponent_t s; // 外積：(C→P) × (C→Q)
  VectorComponent_t t; // 内積：(C→P) ・ (C→Q)
  double theta;        // θ (ラジアン)

  // c，p，q を所望の値に設定する．
  c.x = 1;
  c.y = 2;
  p.x = 2;
  p.y = 3;
  q.x = 3;
  q.y = 5;

  // 回転方向および角度θを計算する．
  Vector2_Diff(&cp, &p, &c);          // cp ← p - c
  Vector2_Diff(&cq, &q, &c);          // cq ← q - c
  s = Vector2_CrossProduct(&cp, &cq); // s ← cp × cq
  t = Vector2_DotProduct(&cp, &cq);   // t ← cp ・ cq
  theta = atan2(s, t);

  // 結果を出力する．
  printf("Ｃ＝(%G, %G)\n", c.x, c.y);
  printf("Ｐ＝(%G, %G)\n", p.x, p.y);
  printf("Ｑ＝(%G, %G)\n", q.x, q.y);
  printf("θ＝%Gラジアン (%G度)\n", theta, RadianToDegree(theta));
  if(s == 0.0) {
    printf("３点 Ｃ，Ｐ，Ｑ は一直線上にある。\n");
  } else {
    printf("Ｃ→Ｐ→Ｑは%s回り。\n", (s > 0.0) ? "左" : "右");
  }

  return EXIT_SUCCESS;
}

実行結果：

Ｃ＝(1, 2)
Ｐ＝(2, 3)
Ｑ＝(3, 5)
θ＝0.197396ラジアン (11.3099度)
Ｃ→Ｐ→Ｑは左回り。

• 幾何学・CG のアルゴリズム集
• ３次元の回転 (原点を通る任意方向回転軸，座標系に依存しないベクトル・回転行列表現)
• 多角形の面積，重心(図心)， 断面Ｎ次モーメントの公式と，向き (頂点列の回転方向) の判別方法
  より一般の，Ｎ個の座標列の回転方向を判別したい場合はこちらへどうぞ．
• DirectX Graphics に関するメモ (Direct3D，DXGI ほか)
• OpenCV に関するメモ

• ゲームつくろー！ 衝突判定編 基礎の基礎編 その１ 内積と外積の使い方 (○×（まるぺけ）つくろーどっとコム)
• ２次元ベクトルの外積の効用 (線形代数学の教科内容の改善に向けて)
  ２次元外積のわかりやすい解説．
• 外積について (高校生のための微分幾何)
• ベクトルの外積に関する考察 ―その概念の広がりと教材性に関して―
  (今岡 光範，広島大学大学院教育学研究科紀要，2008)
• n次元ベクトルの外積の定義 (OKWave)
• Levi-Civita symbol and cross product vector/tensor
  レビ＝チビタの記号 (エディントンのイプシロン) と外積ベクトル／テンソルの関係
• 右と左
  • 時計回り・反時計回り (Wikipedia)
  • 右回り、左回り？
    このページと似た問題．平面上の３点を順番に巡るとき， それが右回りか左回りかを判別せよという問題．
  • なぜ鏡には左右が反対に映るに上下は反対に映らないのか
• 軸性ベクトル (擬ベクトル)
  • ベクトル解析 (物理のかぎしっぽ)
    • 軸性ベクトルと極性ベクトル
    • 擬テンソル
  • ベクトルの外積　軸性ベクトルについて (OKWave)
  • NO.1825　線形な（正直な？）エッセイ（？）其の１　2010.1.4.　DDT
    (Weekend Mathematics コロキウム室 NO.251)
  • 何でＢは軸性ベクトルでＥは極性ベクトルなの？ (OKWave)
• Wikipedia
  • ドット積
  • クロス積 (Cross product)
  • 外積
  • 外積代数
  • 直積 (Outer product)
  • テンソル
  • 擬ベクトル
  • 擬スカラー
  • 向き (Orientation (mathematics))
  • Chirality (mathematics)
  • Cartesian coordinate system
    ⇒ Orientation and handedness
  • エディントンのイプシロン (Levi-Civita symbol)
    英語版には ε_ijk を３次元行列 (３次行列ではない) で表現した図がある．
• ３次元の回転
  • 3D空間における回転の表現形式
  • ７０秒で分る、使える、四元数・４元数・クォータニオン・ Quaternionで回転
• ベクトル (++C++; // 未確認飛行 C)
  ベクトルに使う文字，外積．
• 四次元への扉 (岡田好一ホームページ)
  四次元の超立体について．
• 座標変換とスピノール
  Ｎ次元の回転について．
• 高次元回転行列の補間とその応用
• スピンの回転について教えて下さい

このページの主な更新は Blog でお知らせします．

1. 2006/10/11(水) OKWave に寄せられた質問に回答したことをまとめ， 「２次元幾何学・CG のアルゴリズム集」のページで公開．
2. 2006/10/12(木) 解説を少し修正．
3. 2006/10/27(金) 解説を少し修正，リンクを追加．
4. 2006/10/28(土) 「２次元幾何学・CG のアルゴリズム集」のページから分離独立．
5. 2006/11/05(日) 右回り、左回り？のリンクを追加．
6. 2006/12/13(水) 時々，「なす角」や「Ｃ言語」等で検索して来る人がいるので， 下記を追加．
   • 角度θを求めたい場合
   • サンプルプログラム (Ｃ言語)
7. 2006/12/14(木) サイト内関連ページを追加．
8. 2006/12/24(日) 解説を少し修正・追記．
9. 2007/04/15(日) 外積の独立成分の意味を (１行だけ) 追記，記号を変更．
10. 2007/10/07(日) 相変わらず「なす角」や「３次元 回転」などで検索して来る人が多いので， ２次元ベクトルの回転角 (なす角＋回転方向) の求め方を更新，３次元，Ｎ次元 (なす角，回転方向) の場合を追加．
11. 2007/10/08(月) タイトルを「３点の座標から簡単に角度と回転方向を求める． (2・3・N次元，外積を用いる方法)」に変更．
12. 2008/08/17(日) 関数名変更 (Vector2{Inner,Outer}Product() → Vector2{Dot,Cross}Product())．
13. 2009/02/12(木) 「外積 Ｃ言語」とか「外積 アルゴリズム」とか， あまりにも自分の頭で考えずに検索して来る人が多いので， サンプルプログラムを穴埋めクイズ化．(笑)
14. 2009/03/13(金) 外積を計算する関数・マクロを追加．
15. 2010/06/12(土) 参考・関連図書を追加．
16. 2010/09/13(月) ３点の座標から角度・距離・内積・外積などを計算する Excel ファイルを委託販売開始．
17. 2011/01/22(土) 図 (sin cos tan) を追加．
18. 2011/11/03(木) 参考・関連図書を追加．
19. 2012/08/25(土) サンプルプログラムを改定．
20. 2012/12/29(土) 内積と外積に関するリンクと用語を少し整理．
21. 2013/01/14(月) 穴埋めクイズ化解除．
22. 2014/01/11(土) 図 (atan2(A×B, A・B)) を追加．
23. 2016/04/03(日)
    • DLmarket の販売手数料率変更による価格変更．
    • 角度の記号をθに統一．(φはフォントによって字形が変わるので廃止．)