はじめに

　こないだの記事
　p5.jsで円内の領域の上の均等分布を実装する
を応用して、三角形内の均等分布を実行する記事です。

コード全文

　triangle distribution

/*
  アフィン変換。
  https://qiita.com/inaba_darkfox/items/c9df0726f4ed597fa82c
  p5.jsで円内の領域の上の均等分布を実装する
*/

function setup() {
  createCanvas(400, 400);
  const m0 = [1,2,4,8,5,4,1,1,0];
  const m1 = [3,2,1,8,6,3,2,5,7];
  const m2 = multM(m0,m1);
  console.log(m2); // 27, 34, 35, 72, 66, 51, 11, 8, 4 正解。
  const m3 = invert3x3(m0);
  const m4 = multM(m0,m3);
  const m5 = multM(m3,m0);
  console.log(m4); // 1, 0, 0, 0, 1, 0, 0, 0, 1
  console.log(m5); // 1, 0, 0, 0, 1, 0, 0, 0, 1
  // 数学ってすごいよね。

  // 本題。
  const p0 = createVector(50,80);
  const p1 = createVector(350,34);
  const p2 = createVector(176,360);

  // 行列0,1,1,0,0,1,1,1,1を作る。
  // ベクトル的には(0,0),(1,0),(1,1)のなす三角形。
  // 逆行列を作りp0.x, p1.x, p2.x, p0.y, p1.y, p2.y, 1, 1, 1に左から掛ける。
  // 上記の三角形にp0,p1,p2が移るわけね
  // 得られた行列をベクトルx,y,1に左から掛けるだけ
  const m = invert3x3([0,1,1,0,0,1,1,1,1]);
  const a = multM([p0.x, p1.x, p2.x, p0.y, p1.y, p2.y, 1, 1, 1], m);
  console.log(m);
  console.log(a);

  const vectors = [];
  for(let i=0; i<8000; i++){
    const x = sqrt(random());
    const y = random(x);
    const v = createVector(x,y,1);
    const w = multV(a, v);
    vectors.push(w);
  }

  background(0);
  noStroke();
  colorMode(HSB,1);
  for(const v of vectors){
    fill(0.55+random(0.2), random(), 1);
    circle(v.x, v.y, 2);
  }
}

function multV(m, v){
  // ベクトルを行列に作用させる
  return createVector(
    m[0]*v.x + m[1]*v.y + m[2]*v.z,
    m[3]*v.x + m[4]*v.y + m[5]*v.z,
    m[6]*v.x + m[7]*v.y + m[8]*v.z
  );
}

function multM(m, n){
  // m,nは行ベース。通常の掛け算。
  const result = new Array(9).fill(0);
  for(let k=0; k<3; k++){
    for(let i=0; i<3; i++){
      result[3*k + i] = m[3*k]*n[i] + m[3*k+1]*n[i+3] + m[3*k+2]*n[i+6];
    }
  }
  return result;
}

function invert3x3(n){
  // nは長さ9の配列。結果も長さ9の配列で行ベース。
  const result = new Array(9).fill(0);
  const det = n[0]*n[4]*n[8] + n[1]*n[5]*n[6] + n[2]*n[3]*n[7] - n[2]*n[4]*n[6] - n[1]*n[3]*n[8] - n[0]*n[5]*n[7];
  const indices = [4,8,5,7, 2,7,1,8, 1,5,2,4,
                   5,6,3,8, 0,8,2,6, 2,3,0,5,
                   3,7,4,6, 1,6,0,7, 0,4,1,3];
  for(let i=0; i<9; i++){
    const offset = i*4;
    const a0 = indices[offset];
    const a1 = indices[offset+1];
    const a2 = indices[offset+2];
    const a3 = indices[offset+3];
    result[i] = (n[a0] * n[a1] - n[a2] * n[a3]) / det;
  }
  return result;
}

実行結果

説明

　上記の記事でランダム分布を三角形に対して実行しました。これによると、三つの頂点を

(0,0),~~~(1,0),~~~(1,1)

としたとき、

  for(let i=0; i<8000; i++){
    const x = sqrt(random());
    const y = random(x);
    circle(x,y,1);
  }

とすることで、その中の点をランダムで均等に取得できます。そこで、3つの頂点を任意に選び、

(a,b),~~~(c,d),~~~(e,f)

とします。これらを上の3つの頂点に移してみることを考えます。これは実は2次元の範囲では実現できず、ひとつ次元を挙げる必要があります。すなわち、

\begin{pmatrix} a & c & e\\ b & d & f  \\ 1 & 1 & 1 \end{pmatrix} = A\begin{pmatrix}0&1&1 \\ 0&0&1 \\ 1&1&1\end{pmatrix}

のように、第三成分として1を追加し、3次の正方行列で移すことを考えます。すると、これらは正則行列なので、逆行列を取ることで$A$を決定できます。このとき、この$A$に、先ほどの計算で得られる小さい三角形の中の点$(x,y)$に対して、

A\begin{pmatrix}x \\ y \\ 1\end{pmatrix}

を考えると、

\begin{align}
A\begin{pmatrix}x \\ y \\ 1\end{pmatrix} &= (1-x)\cdot A\begin{pmatrix}0 \\ 0 \\ 1\end{pmatrix} + (x-y)\cdot A\begin{pmatrix}1 \\ 0 \\ 1\end{pmatrix} + y\cdot A\begin{pmatrix}1 \\ 1 \\ 1\end{pmatrix} \\
&= (1-x) \begin{pmatrix}a \\ b \\ 1\end{pmatrix} + (x-y) \begin{pmatrix}c \\ d \\ 1\end{pmatrix} + y \begin{pmatrix}e \\ f \\ 1\end{pmatrix}
\end{align}

となりますね。重心座標の概念を知っている人ならわかると思いますが、これで三角形内部の点になっています。そういうことですね。

行列の計算

　これを実現するために、いくつか行列のメソッドを即席で用意しています。3x3行列の逆行列を取るもの、それらの間の積を取るもの、ベクトルに掛けるもの、以上の3つです。それらを効果的に駆使して、実行しました。どうやら合っているようです。良かった～。

おわりに

　次元を上げれば四面体の内部の点も均等に取ることができると思います。ちょっと大変かもしれないですが...。
　ここまでお読みいただいてありがとうございました。

　行列計算の詳細に関しては割愛しました。知らない人はググって調べてください。

追記

　sqrtを使わない方法があるようです。
　三角形内の一様乱数
　３次元空間、複数三角形内に均一に、点をばらまく
どういう理屈なんだろう...ざっくりいうと、これでいいということです。

  for(let i=0; i<8000; i++){
    const r0 = random();
    const r1 = random();
    const v = createVector(max(r0,r1),min(r0,r1),1);
    const w = multV(a, v);
    vectors.push(w);
  }

結果は同じです。ひぇーすごい。暇があったら考えてみたいですね（おそらくですがどっちが大きくなるかが等確率なので、単純に正方形分布を半分こしているのだと思います。うまいやりかたですね）。

　あとどうでもいいですが、行列$a$は成分から割と簡単に出せるので、きちんと整形した方がそれなりに速くなるかもですね。この記事では過程を重視したので、ある程度の冗長さはやむを得ないと思って書きました。