←Previous | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | Next→
          0, 0, 1, 0];
}

//o = nLetITranslateZ (o, f)
//  Array n
//  double f
function nLetITranslateZ (o, f) {
  o[ 0] = 1;
  o[ 1] = 0;
  o[ 2] = 0;
  o[ 3] = 0;
  o[ 4] = 0;
  o[ 5] = 1;
  o[ 6] = 0;
  o[ 7] = 0;
  o[ 8] = 0;
  o[ 9] = 0;
  o[10] = 1;
  o[11] = f;
  return o;
}

//o = nNewITranslateZ (f)
//  Array o
//  double f
function nNewITranslateZ (f) {
  return [1, 0, 0, 0,
          0, 1, 0, 0,
          0, 0, 1, f];
}

//o = nLetITranslateXYZ (o, d, e, f)
//  Array o
//  double d
//  double e
//  double f
function nLetITranslateXYZ (o, d, e, f) {
  o[ 0] = 1;
  o[ 1] = 0;
  o[ 2] = 0;
  o[ 3] = d;
  o[ 4] = 0;
  o[ 5] = 1;
  o[ 6] = 0;
  o[ 7] = e;
  o[ 8] = 0;
  o[ 9] = 0;
  o[10] = 1;
  o[11] = f;
  return o;
}

//o = nNewITranslateXYZ (d, e, f)
//  Array o
//  double d
//  double e
//  double f
function nNewITranslateXYZ (d, e, f) {
  return [1, 0, 0, d,
          0, 1, 0, e,
          0, 0, 1, f];
}

//n = nTranslateW (n, w)
//  Array n
//  Array w
function nTranslateW (n, w) {
  n[ 3] += w[0];
  n[ 7] += w[1];
  n[11] += w[2];
  return n;
}

//o = nLetTranslateW (o, n, w)
//  Array o
//  Array n
//  Array w
function nLetTranslateW (o, n, w) {
  o[ 0] = n[ 0];
  o[ 1] = n[ 1];
  o[ 2] = n[ 2];
  o[ 3] = n[ 3] + w[0];
  o[ 4] = n[ 4];
  o[ 5] = n[ 5];
  o[ 6] = n[ 6];
  o[ 7] = n[ 7] + w[1];
  o[ 8] = n[ 8];
  o[ 9] = n[ 9];
  o[10] = n[10];
  o[11] = n[11] + w[2];
  return o;
}

//o = nNewTranslateW (n, w)
//  Array o
//  Array n
//  Array w
function nNewTranslateW (n, w) {
  return [n[ 0], n[ 1], n[ 2], n[ 3] + w[0],
          n[ 4], n[ 5], n[ 6], n[ 7] + w[1],
          n[ 8], n[ 9], n[10], n[11] + w[2]];
}

//n = nTranslateX (n, d)
//  Array n
//  double d
function nTranslateX (n, d) {
  n[ 3] += d;
  return n;
}

//o = nLetTranslateX (o, n, d)
//  Array o
//  Array n
//  double d
function nLetTranslateX (o, n, d) {
  o[ 0] = n[ 0];
  o[ 1] = n[ 1];
  o[ 2] = n[ 2];
  o[ 3] = n[ 3] + d;
  o[ 4] = n[ 4];
  o[ 5] = n[ 5];
  o[ 6] = n[ 6];
  o[ 7] = n[ 7];
  o[ 8] = n[ 8];
  o[ 9] = n[ 9];
  o[10] = n[10];
  o[11] = n[11];
  return o;
}

//o = nNewTranslateX (n, d)
//  Array o
//  Array n
//  double d
function nNewTranslateX (n, d) {
  return [n[ 0], n[ 1], n[ 2], n[ 3] + d,
          n[ 4], n[ 5], n[ 6], n[ 7],
          n[ 8], n[ 9], n[10], n[11]];
}

//n = nTranslateY (n, e)
//  Array n
//  double e
function nTranslateY (n, e) {
  n[ 7] += e;
  return n;
}

//o = nLetTranslateY (o, n, e)
//  Array o
//  Array n
//  double e
function nLetTranslateY (o, n, e) {
  o[ 0] = n[ 0];
  o[ 1] = n[ 1];
  o[ 2] = n[ 2];
  o[ 3] = n[ 3];
  o[ 4] = n[ 4];
  o[ 5] = n[ 5];
  o[ 6] = n[ 6];
  o[ 7] = n[ 7] + e;
  o[ 8] = n[ 8];
  o[ 9] = n[ 9];
  o[10] = n[10];
  o[11] = n[11];
  return o;
}

//o = nNewTranslateY (n, e)
//  Array o
//  Array n
//  double e
function nNewTranslateY (n, e) {
  return [n[ 0], n[ 1], n[ 2], n[ 3],
          n[ 4], n[ 5], n[ 6], n[ 7] + e,
          n[ 8], n[ 9], n[10], n[11]];
}

//n = nTranslateZ (n, f)
//  Array n
//  double f
function nTranslateZ (n, f) {
  n[11] += f;
  return n;
}

//o = nLetTranslateZ (o, n, f)
//  Array o
//  Array n
//  double f
function nLetTranslateZ (o, n, f) {
  o[ 0] = n[ 0];
  o[ 1] = n[ 1];
  o[ 2] = n[ 2];
  o[ 3] = n[ 3];
  o[ 4] = n[ 4];
  o[ 5] = n[ 5];
  o[ 6] = n[ 6];
  o[ 7] = n[ 7];
  o[ 8] = n[ 8];
  o[ 9] = n[ 9];
  o[10] = n[10];
  o[11] = n[11] + f;
  return o;
}

//o = nNewTranslateZ (n, f)
//  Array o
//  Array n
//  double f
function nNewTranslateZ (n, f) {
  return [n[ 0], n[ 1], n[ 2], n[ 3],
          n[ 4], n[ 5], n[ 6], n[ 7],
          n[ 8], n[ 9], n[10], n[11] + f];
}

//n = nTranslateXYZ (n, d, e, f)
//  Array n
//  double d
//  double e
//  double f
function nTranslateXYZ (n, d, e, f) {
  n[ 3] += d;
  n[ 7] += e;
  n[11] += f;
  return n;
}

//o = nLetTranslateXYZ (o, n, d, e, f)
//  Array o
//  Array n
//  double d
//  double e
//  double f
function nLetTranslateXYZ (o, n, d, e, f) {
  o[ 0] = n[ 0];
  o[ 1] = n[ 1];
  o[ 2] = n[ 2];
  o[ 3] = n[ 3] + d;
  o[ 4] = n[ 4];
  o[ 5] = n[ 5];
  o[ 6] = n[ 6];
  o[ 7] = n[ 7] + e;
  o[ 8] = n[ 8];
  o[ 9] = n[ 9];
  o[10] = n[10];
  o[11] = n[11] + f;
  return o;
}

//o = nNewTranslateXYZ (n, d, e, f)
//  Array o
//  Array n
//  double d
//  double e
//  double f
function nNewTranslateXYZ (n, d, e, f) {
  return [n[ 0], n[ 1], n[ 2], n[ 3] + d,
          n[ 4], n[ 5], n[ 6], n[ 7] + e,
          n[ 8], n[ 9], n[10], n[11] + f];
}

//------------------------------------------------------------------------
//translationComponent 移動成分
//
//    [n11 n12 n13 n14]     [1 0 0 n14]
//    [n21 n22 n23 n24] ==> [0 1 0 n24]
//    [n31 n32 n33 n34]     [0 0 1 n34]
//    [ 0   0   0   1 ]     [0 0 0  1 ]
//

//o = nTranslationComponent (o)
//  Array o
function nTranslationComponent (o) {
  o[ 0] = 1;
  o[ 1] = 0;
  o[ 2] = 0;
  o[ 4] = 0;
  o[ 5] = 1;
  o[ 6] = 0;
  o[ 8] = 0;
  o[ 9] = 0;
  o[10] = 1;
  return o;
}

//o = nLetTranslationComponent (o, n)
//  Array o
//  Array n
function nLetTranslationComponent (o, n) {
  o[ 0] = 1;
  o[ 1] = 0;
  o[ 2] = 0;
  o[ 3] = n[ 3];
  o[ 4] = 0;
  o[ 5] = 1;
  o[ 6] = 0;
  o[ 7] = n[ 7];
  o[ 8] = 0;
  o[ 9] = 0;
  o[10] = 1;
  o[11] = n[11];
  return o;
}

//o = nNewTranslationComponent (n)
//  Array o
//  Array n
function nNewTranslationComponent (n) {
  return [1, 0, 0, n[ 3],
          0, 1, 0, n[ 7],
          0, 0, 1, n[11]];
}



//========================================================================
//多角形
//
//  vertexArray  [v0,v1,…]  頂点の配列
//  faceArray    [f0,f1,…]  面の配列
//
//  頂点
//    v  [vx,vy,vz]  座標の配列
//
//  面
//    f  [vn0,vn1,…]  頂点の番号の配列
//

//linecap 線の端の形
//
//    butt
//      __
//         │
//         │
//         ●
//         │
//         │
//       ̄ ̄
//    round
//      ___
//            \
//             │
//         ●  │
//             │
//            /
//       ̄ ̄ ̄
//    square
//      ____
//             │
//             │
//         ●  │
//             │
//             │
//       ̄ ̄ ̄ ̄
//
let pgStrokeLinecap = "butt";

//linejoin 線の繋ぎ方
//
//    bevel
//        ───────┐
//                      │\
//                      ● \
//                        \
//        ──○          /
//          /          /
//                    /
//                  /
//                /
//
//    miter
//        ───────┬────○
//                      │      /
//                      ●    /
//                        \/
//        ──○          /
//          /          /
//                    /
//                  /
//                /
//
//    round
//        ───────┐-_
//                      │  \
//                      ●   |
//                        \/
//        ──○          /
//          /          /
//                    /
//                  /
//                /
//
let pgStrokeLinejoin = "bevel";

//  linewidth 線の幅

//  miterlimit 線の繋ぎ目の太さの限度
//    θ<=2*asin(1/miterlimit)でmiterからbevelに切り替わる
//      perl -e "use Math::Trig;for$miterlimit(1,sqrt(2),2,2.6131259297527,3.86370330515627){print $miterlimit,' ',rad2deg(2*asin(1/$miterlimit)),qq'\n';}"
//      1 180
//      1.4142135623731 90
//      2 60
//      2.6131259297527 45.000000000001
//      3.86370330515627 30
let pgStrokeMiterlimit = 4;

//pg = pgPolyline (pl, linewidth)
//  折れ線から多角形を作る
//
//  pl  折れ線
//    vertexArray
//      頂点の配列
//      頂点毎に、座標を列挙する
//    edgeArray
//      辺の配列
//      頂点毎に、辺で繋がっている頂点の番号を列挙する
//      1つの頂点が3つ以上の辺に所属しているときは、辺で繋がっている頂点の番号を左回りに書く
//
//  例
//
//    pl = {
//    pl$vertexArray: [[-1, 0],  //0
//                     [0, 0],  //1
//                     [1, 0],  //2
//                     [0, 1]],  //3
//    pl$edgeArray: [[1],  //0
//                   [0, 2, 3],  //1
//                   [1],  //2
//                   [1]]  //3
//      }
//
//  右隣の頂点を辿って一周する
//
//               3
//              ●
//              │
//            ↓│↑
//          ←  │  ←
//     0●───●───●2
//          →   1  →
//
//
//    -180と180は線の端
//        ──────★
//                    │
//      ◎──────●
//                    │
//        ──────★
//
//    -135→-135/2-90=-157.5
//        ───────┬────☆
//                  /\│      /
//      ○────/──●    /
//              /    /│\/
//        ──★──/─┘/
//          /    /    /
//              /    /
//            ○    /
//                /
//
//    -90→-90/2-90=-135
//        ──┬─┬─☆
//            │  │  │
//      ○──┼─●─┤
//            │  │  │
//        ──★─┼─┤
//            │  │  │
//            │  │  │
//                ○
//
//    -45→-45/2-90=-112.5
//        ────┐☆
//                │/\
//      ○────●    \
//              /│\    \
//        ───★┘  \    \
//                \    \
//                  \    ○
//                    \
//
//    0→0/2-90=-90
//        ────☆────
//                │
//      ○────●────○
//                │
//        ────★────
//
//    45→45/2-90=-67.5
//                  /    ○
//                /    /
//        ───☆┐  /
//              \│/    /
//      ○────●    /
//                │\/
//        ────┘★
//
//    90→90/2-90=-45
//                ○
//            │  │  │
//            │  │  │
//        ──☆─┼─┤
//            │  │  │
//      ○──┼─●─┤
//            │  │  │
//        ──┴─┴─★
//
//    135→135/2-90=-22.5
//                \
//            ○    \
//              \    \
//          \    \    \
//        ──☆──\─┐\
//              \    \│/\
//      ○────\──●    \
//                  \/│      \
//        ───────┴────★
function pgPolyline (pl, linewidth) {
  const step = rad (15);
  let va = pl.pl$vertexArray;
  let ea = pl.pl$edgeArray;
  let wh = linewidth * 0.5;  //太さの半分
  let va2 = [];
  let fa2 = [];
  let wa = [];  //頂点と辺の処理済みフラグ
  for (let i = 0; i < ea.length; i++) {
    wa[i] = [];
    for (let j = 0; j < ea[i].length; j++) {
      wa[i][j] = false;
    }
  }
  let vn1 = 0;  //頂点0から
  let en12 = 0;  //枝0へ進む
  let vstart = 0;  //面の頂点の開始番号
loop:
  for (;;) {
    wa[vn1][en12] = true;
    //頂点vn1と枝en12が選択されている
    let vn2 = ea[vn1][en12];  //vn2=頂点vn1の枝en12の先にある頂点
    let en21 = 0;  //en21=頂点vn2から見た頂点vn1の枝
    while (ea[vn2][en21] != vn1) {
      en21++;
    }
    let en23 = en21 + 1 < ea[vn2].length ? en21 + 1 : 0;  //en23=頂点vn2から見た頂点vn1の枝の次の枝=頂点vn2から見た頂点vn3の枝
    let vn3 = ea[vn2][en23];  //vn3=頂点vn2の枝en23の先にある頂点
    let v1 = va[vn1];
    let v2 = va[vn2];
    let v3 = va[vn3];
    let h12 = vHat (vNewSubtractV (v2, v1));  //1→2の単位ベクトル
    let h23 = vHat (vNewSubtractV (v3, v2));  //2→3の単位ベクトル
    let a12 = vAzimuth (h12);  //1→2の方位角
    let a23 = vAzimuth (h23);  //2→3の方位角
    let a13 = a23 - a12 - floor ((a23 - a12 + PI) / TWOPI) * TWOPI;  //1→2→3で曲がる角度。-PI~PI
    if (a13 <= -PI || PI <= a13) {  //線の端
      if (pgStrokeLinecap == "square") {
        va2.push (vTranslateV (vRotate (vScaleS ([1, -1], wh), a12), v2));
        va2.push (vTranslateV (vRotate (vScaleS ([1, 1], wh), a12), v2));
      } else {
        let t1 = -PI / 2;
        let t2 = PI / 2;
        let n = pgStrokeLinecap == "round" ? max2 (1, roundn ((t2 - t1) / step)) : 1;
        for (let i = 0; i <= n; i++) {
          let t = (n - i) / n * t1 + i / n * t2;
          va2.push (vTranslateV (vRotate (vScaleS ([cos (t), sin (t)], wh), a12), v2));
        }
      }
    } else if (a13 == 0) {  //直進は角を省略
    } else if (a13 < 0 ||  //右折の内側は1点で曲がる
               (pgStrokeLinejoin == "miter" &&  //左折の外側で1点で曲がる
                2 * asin (1 / pgStrokeMiterlimit) < PI - a13)) {  //極端に尖らない
      va2.push (vTranslateV (vRotate (vScaleS ([-1 / tan (a13 * 0.5 - PI / 2), -1], wh), a12), v2));
    } else {  //左折の外側で2点以上で曲がる
      let t1 = -PI / 2;
      let t2 = -PI / 2 + a13;
      let n = pgStrokeLinejoin == "round" ? max2 (1, roundn ((t2 - t1) / step)) : 1;
      for (let i = 0; i <= n; i++) {
        let t = (n - i) / n * t1 + i / n * t2;
        va2.push (vTranslateV (vRotate (vScaleS ([cos (t), sin (t)], wh), a12), v2));
      }
    }
    vn1 = vn2;
    en12 = en23;
    if (wa[vn1][en12]) {  //処理済みの頂点と辺が出てきた
      //面を追加する
      let f2 = [];
      for (let i = vstart; i < va2.length; i++) {
        f2.push (i);
      }
      fa2.push (f2);
      vstart = va2.length;
      //未処理の頂点と辺を探す
      for (let i = 0; i < ea.length; i++) {
        for (let j = 0; j < ea[i].length; j++) {
          if (!wa[i][j]) {  //未処理の頂点と辺が見つかった
            vn1 = i;
            en12 = j;
            continue loop;
          }
        }
      }
      //未処理の頂点と辺が残っていない
      break;
    }
  }  //loop: for
  return { pg$vertexArray: va2, pg$faceArray: fa2 };
}

//s = pgToString (pg)
//  文字列にする
//  配列の配列を文字列にするとき内側の配列を囲む括弧がないと外側の配列の要素の区切りがどこだかわからない
function pgToString (pg) {
  let va = pg.pg$vertexArray;
  let fa = pg.pg$faceArray;
  let s = "";
  s += "{va:[";
  for (let i = 0; i < va.length; i++) {
    if (0 < i) {
      s += ",";
    }
    s += "[" + va[i].join (",") + "]";
  }
  s += "],fa:[";
  for (let i = 0; i < fa.length; i++) {
    if (0 < i) {
      s += ",";
    }
    s += "[" + fa[i].join (",") + "]";
  }
  s += "]}";
  return s;
}

const RE_NUM = "([-+]?\\d+(?:\\.\\d*)?(?:[Ee][-+]?\\d+)?)";

//pg = pgTurtle (str)
//  タートルで多角形を作る
//  引数の角度の単位は度
//  初期状態は u j0,0 a0
//  a a     aを向く
//  cl a,r  半径rで左にa曲がる。0<a,0<r
//  cr a,r  半径rで右にa曲がる。0<a,0<r
//  d       ペンを下ろす
//  f d     d進む。0<d
//  j x,y   x,yまで移動する
//  l a     左にa曲がる。0<a
//  m x,y   x,yだけ移動する
//  r a     右にa曲がる。0<a
//  u       多角形を閉じてペンを上げる
function pgTurtle (str) {
  const step = 5;
  const eps = 1e-8;
  let tx = 0;  //現在の位置
  let ty = 0;
  let ta = 0;  //現在の方位角
  let td = false;  //false=ペンが上がっている,true=ペンが下りている
  let va = [];
  let fa = [];
  let vn = -1;
  let fn = -1;
  let rex = new RegExp ("\\s*(?:" +
                        "(a\\s*" + RE_NUM + ")|" +  //1,2  a a
                        "(cl\\s*" + RE_NUM + "\\s*,\\s*" + RE_NUM + ")|" +  //3,4,5  cl a,r
                        "(cr\\s*" + RE_NUM + "\\s*,\\s*" + RE_NUM + ")|" +  //6,7,8  cr a,r
                        "(d\\s*)|" +  //9  d
                        "(f\\s*" + RE_NUM + ")|" +  //10,11  f d
                        "(j\\s*" + RE_NUM + "\\s*,\\s*" + RE_NUM + ")|" +  //12,13,14  j x,y
                        "(l\\s*" + RE_NUM + ")|" +  //15,16  l a
                        "(m\\s*" + RE_NUM + "\\s*,\\s*" + RE_NUM + ")|" +  //17,18,19  m x,y
                        "(r\\s*" + RE_NUM + ")|" +  //20,21  r a
                        "(u\\s*)" +  //22  u
                        ")", "g");
  let mat;
  while (mat = rex.exec (str)) {
    if (mat[1]) {  //a a  aを向く
      ta = parseFloat (mat[2]);
      ta -= floor (ta / 360) * 360;
    } else if (mat[3]) {  //cl a,r  半径rで左にa曲がる。0<a,0<r
      let a = parseFloat (mat[4]);
      let r = parseFloat (mat[5]);
      let ox = tx - r * sin (rad (ta));  //円の中心
      let oy = ty + r * cos (rad (ta));
      let b1 = ta - 90;  //開始角
      ta += a;
      let b2 = ta - 90;  //終了角
      ta -= floor (ta / 360) * 360;
      if (td) {
        let n = max2 (1, roundn ((b2 - b1) / step));  //分割数
        for (let i = 0; i < n; i++) {
          let b = (n - i) / n * b1 + i / n * b2;
          tx = ox + r * cos (rad (b));
          ty = oy + r * sin (rad (b));
          va[++vn] = [tx, ty];
          fa[fn].push (vn);
        }
      }
      tx = ox + r * cos (rad (b2));
      ty = oy + r * sin (rad (b2));
    } else if (mat[6]) {  //cr a,r  半径rで右にa曲がる。0<a,0<r
      let a = parseFloat (mat[7]);
      let r = parseFloat (mat[8]);
      let ox = tx + r * sin (rad (ta));  //円の中心
      let oy = ty - r * cos (rad (ta));
      let b1 = ta + 90;  //開始角
      ta -= a;
      let b2 = ta + 90;  //終了角
      ta -= floor (ta / 360) * 360;
      if (td) {
        let n = max2 (1, roundn ((b1 - b2) / step));  //分割数
        for (let i = 0; i < n; i++) {
          let b = (n - i) / n * b1 + i / n * b2;
          tx = ox + r * cos (rad (b));
          ty = oy + r * sin (rad (b));
          va[++vn] = [tx, ty];
          fa[fn].push (vn);
        }
      }
      tx = ox + r * cos (rad (b2));
      ty = oy + r * sin (rad (b2));
    } else if (mat[9]) {  //d  ペンを下ろす
      if (!td) {
        td = true;
        fa[++fn] = [];
      }
    } else if (mat[10]) {  //f d  d進む。0<d
      if (td) {
        va[++vn] = [tx, ty];
        fa[fn].push (vn);
      }
      let d = parseFloat (mat[11]);
      tx += d * cos (rad (ta));
      ty += d * sin (rad (ta));
    } else if (mat[12]) {  //j x,y  x,yまで移動する
      if (td) {
        va[++vn] = [tx, ty];
        fa[fn].push (vn);
      }
      tx = parseFloat (mat[13]);
      ty = parseFloat (mat[14]);
    } else if (mat[15]) {  //l a  左にa曲がる。0<a
      ta += parseFloat (mat[16]);
      ta -= floor (ta / 360) * 360;
    } else if (mat[17]) {  //m x,y  x,yだけ移動する
      if (td) {
        va[++vn] = [tx, ty];
        fa[fn].push (vn);
      }
      tx += parseFloat (mat[18]);
      ty += parseFloat (mat[19]);
    } else if (mat[20]) {  //r a  右にa曲がる。0<a
      ta -= parseFloat (mat[21]);
      ta -= floor (ta / 360) * 360;
    } else if (mat[22]) {  //u  多角形を閉じてペンを上げる
      if (td) {
        let v = va[fa[fn][0]];
        let vx = v[0];
        let vy = v[1];
        if (eps < abs (tx - vx) || eps < abs (ty - vy)) {
          va[++vn] = [tx, ty];
          fa[fn].push (vn);
        }
        tx = vx;
        ty = vy;
        ta = 0;
        td = false;
      }
    }
  }
  if (td) {
    let v = va[fa[fn][0]];
    let vx = v[0];
    let vy = v[1];
    if (eps < abs (tx - vx) || eps < abs (ty - vy)) {
      va[++vn] = [tx, ty];
      fa[fn].push (vn);
    }
    tx = vx;
    ty = vy;
    ta = 0;
    td = false;
  }
  return { pg$vertexArray: va, pg$faceArray: fa };
}



//========================================================================
//多面体
//
//  ph$vertexArray  [v0,v1,…]  頂点の配列。頂点は座標の配列
//  ph$faceArray    [f0,f1,…]  面の配列。面は頂点の番号の配列
//  ph$openArray    [o0,o1,…]  面の開放フラグの配列。0=閉じる,1=閉じない
//  ph$normalArray  [n0,n1,…]  面の法線ベクトルの配列
//  ph$reciprocalMatrix         単位立方体変換行列。外接直方体を(0,0,0)-(1,1,1)の単位立方体に変換する行列。外接直方体の辺は軸に並行とは限らない
//
//  座標系
//    ここでは、右が+x、奥が+y、上が+zの右手座標系を用いる
//            +z
//             | +y
//             |/
//      -x-----+-----+x
//            /|
//          -y |
//            -z
//


//------------------------------------------------------------------------
//コンストラクタ

//ph = phNew (va, fa, oa)
//  頂点の配列と面の配列と面の開放フラグの配列から多面体を作る
//  頂点の配列と面の配列は必須。面の開放フラグは省略可
function phNew (va, fa, oa) {
  return phMakeFields ({
    ph$vertexArray: va,
    ph$faceArray: fa,
    ph$openArray: oa || null
    });
}

//s = pgArea ([[x[0],y[0]],…,[x[n-1],y[n-1]])
//  凸多角形の面積
//  頂点は左回りであること
//  右回りのときは面積に負符号を付けた値が返る(左回りか右回りかの判別に使えるように絶対値をとらずそのまま返す)
//
//  凸多角形(x[k],y[k]) (k=0..n-1)の面積
//    x[n]=x[0],y[n]=y[0]とする
//    平行移動しても面積は変わらないので第一象限で考える
//    各頂点からx軸に垂線を下ろして辺と垂線とx軸で台形を作る
//    頂点の順序が左回りの場合
//      上(+y側)の辺はkが増えるとxが減少するので、上の辺を含む台形の面積は、
//        (x[k]-x[k+1])*(y[k]+y[k+1])/2 (k=0..n-1)
//      下(-y側)の辺はkが増えるとxが増加するので、下の辺を含む台形の面積は、
//        (x[k+1]-x[k])*(y[k]+y[k+1])/2 (k=0..n-1)
//      下の辺を含む台形の面積の式は上の辺を含む台形の面積の式の符号を逆にしたものなので、
//      上の辺を含む台形の面積の式をすべての辺について合計すると、
//      上の辺を含む台形の面積の合計から下の辺を含む台形の面積の合計を引いた値、すなわち、多角形の面積が求まる
//        S=1/2*sum[k=0..n-1]{(x[k]-x[k+1])*(y[k]+y[k+1])}
//    頂点の順序が右回りの場合
//      左回りと符号が逆になる
//    与えられた頂点の順序が左回りか右回りかわからないとき、左回りの場合のSが正ならば左回り、負ならば右回りである
//
function pgArea (va) {
  let s = 0;
  let v = va[0];
  let x1 = v[0];
  let y1 = v[1];
  for (let i = va.length - 1; 0 <= i; i--) {
    v = va[i];
    let x0 = v[0];
    let y0 = v[1];
    s += (x0 - x1) * (y0 + y1);
    x1 = x0;
    y1 = y0;
  }
  return s * 0.5;
}

//ph = phMakeFields (ph)
//  足りないフィールドを作る
//  頂点の配列と面の配列は必須
//  既存のフィールドを再構築させたいときはフィールドをnullにするかdeleteしてから呼び出すこと
function phMakeFields (ph) {
  const eps = 1e-4;
  //開放フラグを作る
  if (!ph.ph$openArray) {
    let fa = ph.ph$faceArray;
    let oa = [];
    for (let i = 0; i < fa.length; i++) {
      oa[i] = 0;
    }
    ph.ph$openArray = oa;
  }
  //法線ベクトルを作る
  if (!ph.ph$normalArray) {
    let va = ph.ph$vertexArray;
    let fa = ph.ph$faceArray;
    let na = [];
    for (let i = 0; i < fa.length; i++) {
      let f = fa[i];
      let va1 = [];
      for (let j = 0; j < f.length; j++) {
        va1[j] = va[f[j]].concat ();
      }
      //最初の2つの頂点と、同一直線上にない3つ目の頂点を選ぶ
      let a = va1[0];
      let b = va1[1];
      let c;
      for (let j = 2; j < va1.length; j++) {
        c = va1[j];
        if (eps < wAreaOfTriangle (a, b, c)) {  //三角形の面積が0でなければ同一直線上にない
          break;
        }
      }
      //三角形の法線ベクトルを求める
      //  この時点では凸多角形でないとき法線ベクトルが多面体の内側を向いてしまっている可能性がある
      let n = wNewNormal (a, b, c);
      let nx = n[0];
      let ny = n[1];
      let nz = n[2];
      //法線ベクトルが+zを向くように面を回転させてxy平面と平行にする
      //                    z
      //                    |
      //                    |  nx n
      //                    +-----+
      //                    |    /
      //                  nz|   /
      //        --          |  / hypot(nx,nz)
      //          ~~--      |t/
      //              ~~--  |/
      //                  ~~+-----+----------- x
      //                      ~~-- t
      //                          ~~--
      //                              ~~-
      //  頂点と法線ベクトルを+x→+zにtxz=atan2(nx,nz)回転させて法線ベクトルのx成分を0にする
      //  txz=atan2(nx,nz)
      //  cos(txz)=nz/hypot(nx,nz)
      //  sin(txz)=nx/hypot(nx,nz)
      let n2x;
      let n2y;
      let n2z;
      let va2 = [];
      let r = hypot2 (nx, nz);
      if (r < eps) {
        n2x = nx;
        n2y = ny;
        n2z = nz;
        for (let j = 0; j < va1.length; j++) {
          va2[j] = va1[j].concat ();
        }
      } else {
        let cxz = nz / r;
        let sxz = nx / r;
        for (let j = 0; j < va1.length; j++) {
          let v = va1[j];
          va2[j] = [cxz * v[0] - sxz * v[2],
                    v[1],
                    sxz * v[0] + cxz * v[2]];
        }
        n2x = cxz * nx - sxz * nz;
        n2y = ny;
        n2z = sxz * nx + cxz * nz;
      }
      //  頂点と法線ベクトルを+y→+zにtyz=atan2(ny,nz)回転させて法線ベクトルのy成分を0にする
      //  tyz=atan2(ny,nz)
      //  cos(tyz)=nz/hypot(ny,nz)
      //  sin(tyz)=ny/hypot(ny,nz)
      let n3x;
      let n3y;
      let n3z;
      let va3 = [];
      let r2 = hypot2 (n2y, n2z);
      if (r2 < eps) {
        n3x = n2x;
        n3y = n2y;
        n3z = n2z;
        for (let j = 0; j < va2.length; j++) {
          va3[j] = va2[j].concat ();
        }
      } else {
        let cyz = n2z / r2;
        let syz = n2y / r2;
        for (let j = 0; j < va2.length; j++) {
          let v2 = va2[j];
          va3[j] = [v2[0],
                    cyz * v2[1] - syz * v2[2],
                    syz * v2[1] + cyz * v2[2]];
        }
        n3x = n2x;
        n3y = cyz * n2y - syz * n2z;
        n3z = syz * n2y + cyz * n2z;
      }
      //xy平面上で多角形の面積を計算して負になったときは法線ベクトルを逆向きにする
      if (pgArea (va3) < 0) {  //多角形の面積が負
        wReverse (n);  //法線ベクトルを逆向きにする
      }
      na[i] = n;
    }
    ph.ph$normalArray = na;
  }
  //単位立方体変換行列を作る
  if (!ph.ph$reciprocalMatrix) {
    let va = ph.ph$vertexArray;
    let xmin = Infinity;
    let ymin = Infinity;
    let zmin = Infinity;
    let xmax = -Infinity;
    let ymax = -Infinity;
    let zmax = -Infinity;
    for (let i = 0; i < va.length; i++) {
      let v = va[i];
      let x = v[0];
      let y = v[1];
      let z = v[2];
      if (x < xmin) {
        xmin = x;
      } else if (xmax < x) {
        xmax = x;
      }
←Previous | 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | Next→