p5.jsでクォータニオンを使ってチューブを作る

はじめに

　チューブを作ろう。

やりたいこと

　ベジエ曲線の点列を作って座標系を滑らかに動かしチューブを作る。

ざっくりとした方針

　まず曲線上の点をevenlySpacingで均等にしておく。次に始点において次の点と差を取って接線ベクトルにして、それに直交するベクトルを取って正規直交基底のクォータニオンを作る。これを滑らかに変形していく。やり方は、次の点で同じように接線ベクトルを作り、それとの差が小さいならそのまま同じものを使う。ズレているならそれに応じて補正を掛ける。そうして最初まで戻る。
　閉じている場合、最初と最後が、一致してない可能性があるので、その場合は最初と最後でクォータニオンの商を取って、それに応じて全体を補正する。これで綺麗につながる。

コード全文

　getFrenetSerretSystem

const {Vecta, Quarternion} = fox3Dtools;
const {evenlySpacing} = foxApplications;

const getBezierPoints = (a0,a1,a2,b0,b1,b2,c0,c1,c2,d0,d1,d2, detail) => {
  const result = [];
  for(let k=0; k<=detail; k++){
    const t = k/detail;
    const ta = (1-t)*(1-t)*(1-t);
    const tb = 3*(1-t)*(1-t)*t;
    const tc = 3*(1-t)*t*t;
    const td = t*t*t;
    result.push(new Vecta(
      ta*a0 + tb*b0 + tc*c0 + td*d0,
      ta*a1 + tb*b1 + tc*c1 + td*d1,
      ta*a2 + tb*b2 + tc*c2 + td*d2
    ));
  }
  return result;
}

const getBezierPointsV = (va, vb, vc, vd, detail) => {
  return getBezierPoints(
    va.x, va.y, va.z, vb.x, vb.y, vb.z, vc.x, vc.y, vc.z, vd.x, vd.y, vd.z, detail
  );
} 

function setup() {
  createCanvas(400, 400, WEBGL);

  const points = [];
  points.push(...getBezierPoints(0,0,0,0,100,0,100,100,0,100,0,0, 200));
  points.push(...getBezierPoints(100,0,0,100,-100,0,100,-100,100,100,0,100, 200));
  points.push(...getBezierPoints(100,0,100,100,100,100,0,100,100,0,0,100, 200));
  points.push(...getBezierPoints(0,0,100,0,-100,100,-100,-100,100,-100,0,100, 200));
  points.push(...getBezierPoints(-100,0,100,-100,100,100,-100,100,0,-100,0,0, 200));
  points.push(...getBezierPoints(-100,0,0,-100,-100,0,0,-100,0,0,0,0, 200));

  // 300分割だとほぼ4だそうです。
  evenlySpacing(points, {closed:true, partition:300, showDetail:true});
  const fss = getFrenetSerretSystem(points, {closed:true, showDetail:true});

  const geom = new p5.Geometry();
  // 頂点ですかね
  const radius = 8;
  for(let k=0; k<300; k++){
    const p = points[k];
    const q = fss[k];
    const axes = q.getAxes();
    for(let i=0; i<12; i++){
      const t = TAU*i/12;
      const n = axes.x.mult(Math.cos(t), true).addScalar(axes.y, Math.sin(t));
      geom.vertexNormals.push(createVector(n.x, n.y, n.z));
      const v = p.addScalar(n, radius, true);
      geom.vertices.push(createVector(v.x, v.y, v.z));
      const g = Math.pow(0.5+0.5*Math.cos(t),0.5);
      geom.vertexColors.push(g,g,g,1);
    }
  }
  geom.faces = get2DFaces(12, 300);

  this._renderer.createBuffers("myTube", geom);

  // カリングチェック
  const gl = this._renderer.GL;
  gl.enable(gl.CULL_FACE);
  gl.cullFace(gl.FRONT); // p5はBACK描画である。自分のライブラリならBACKでないとまずい

  draw = () => {
    background(255);
    orbitControl(1,1,1,{freeRotation:true});
    lights();
    fill(255);
    noStroke();

    this._renderer.drawBuffers("myTube");
  }
}

// wは帯に垂直な方向の分割数、hは帯の分割の長さ。
// セグメントの個数。だからh=8の場合分点は0,1,2,...,8と9つできる。
// ループの場合は8が0に一致するから8つだが、切れるなら9つである。wも同様。
function get2DFaces(w=8, h=200){
  const result = [];
  // WEBGLで考えるがp5でも通用する。私のライブラリは逆向きなのでそれで考えさせて
  // いただく。いずれ供用したいので。ただどっちでも結果は同じ。
  // 進行方向に対していわゆる右ねじの向きに点が並び上昇していく。
  // p5だと時計回り（左ねじ）。
  for(let y=0; y<h; y++){
    for(let x=0; x<w; x++){
      const ld = x + w*y;
      const rd = (x+1)%w + w*y;
      const lu = x + w*((y+1)%h);
      const ru = (x+1)%w + w*((y+1)%h);
      result.push([ld, rd, ru], [ld, ru, lu]);
    }
  }
  return result;
}

// 戻り値：pointsと同じ長さのクォータニオンの列
// これも洗練させて（console.logは除くかdetailに含めてshowDetailで表示するように
// するなどして）供用できるようにしましょ。
function getFrenetSerretSystem(points, options = {}){
  const {closed = false, showDetail = false} = options;
  const z0 = points[1].sub(points[0], true).normalize();
  const x0 = Vecta.getOrtho(z0);
  const y0 = z0.cross(x0, true);

  const result = [Quarternion.getFromAxes(x0, y0, z0)];
  const currentTangent = z0.copy();

  // closedの場合に最初の点と重複させる（意図的に）
  for(let i=1; i <= points.length; i++){
    // non-closedの場合は直前で切る
    if(!closed && i === points.length-1) break;
    // closedの場合は最終的に始点を参照する
    const nextTangent = points[(i+1)%points.length].sub(points[i%points.length], true).normalize();
    const axis = currentTangent.cross(nextTangent, true);
    if(axis.magSq() < Number.EPSILON){
      result.push(result[result.length-1].copy());
      currentTangent.set(nextTangent);
      continue;
    }else{
      axis.normalize();
      // axisのまわりに回転させる角度の計算
      const differenceAngle = currentTangent.angleBetween(nextTangent, axis);
      // globalRotateを使う。q0は直前でしたね。バカ...
      const q0 = result[result.length-1];
      const q1 = q0.globalRotate(axis, differenceAngle, true);

      result.push(q1);
      currentTangent.set(nextTangent);
    }
  }

  if(!closed){
    // 同じものを重複して置く
    result.push(result[result.length-1].copy());
  }

  if(closed){
    const diff = result[0].conj(true).multQ(result[result.length-1]);
    // x,y成分が0（ほぼ0）ですね。OKです。角度...
    if(showDetail){ console.log(`difference of begin~end: ${diff.show()}`); }

    // ああそうか
    // まず角度の算出にはtheta/2なので2倍すると
    // そんでもってマイナスで正解だと
    // OKです
    // TAUに近い場合におかしくなるのを防ぐために、wが負の時は-1を掛ける
    if(diff.w < 0){ diff.mult(-1); }
    const diffAngle = Math.atan2(diff.z, diff.w)*2;
    if(showDetail){ console.log(`angle difference: ${diffAngle}`); }
    // そういうわけでこれを...600で割って...
    // i番を-i*diffAngle/600だけzの周りに回す。ローカルで。
    for(let i=0; i<=points.length; i++){
      result[i].localRotate(0,0,1,-i*diffAngle/points.length);
    }
    // 末尾は要らないので弾く
    result.pop();
  }
  return result;
}

実行結果

おわりに

　ここまでお読みいただいてありがとうございました。