pbd_constraints.h

• compareIntArrays
  2つの整数配列が等しいかどうかを比較
• sortRemoveDuplicates
  整数配列から重複する要素を削除し、要素をソートした新しい配列を返す
• calcConstraintTopology
  拘束トポロジーに関するデータIDを取得
• constraintAlias
  入力された拘束タイプの文字列に対応するエイリアスを返す
• constraintHash
  拘束タイプのエイリアスに基づいてハッシュ値を生成する
• isTriARAP
  拘束タイプがTriARAP(三角形の角度による最小二乗法)であるかどうかを判別
• isTetARAPVol
  TetARAPにおける体積維持に使用する拘束の種類を判別
• isTetARAP
  拘束タイプがTetARAP(四面体における剛体補間)であるかどうかを判別
• isNonLinearARAP
  拘束タイプが非線形ARAPであるかどうかを判別
• isTetFiber
  拘束タイプがTetFiber（特定方向に沿って圧縮することができる四面体拘束）であるかどうかを判別
• computeDistanceRestLength
  2ポイント間の距離を計算
• createDistanceConstraint
  入力ポイントから最も近いプリミティブ上の位置を計算し、restvector、restlengthなどを作成
• createAttachConstraint
  入力ポイントから最も近いプリミティブ上の位置を計算し、restvector、restlengthなどを作成
• projectToLine
  与えられた位置ベクトルを指定方向の直線上に射影
• fromTriangleSpaceXform
  三角形からワールド空間への変換行列を返す
• fromTriangleSpace
  3ポイントで構成される平面からワールド空間への回転ベクトルを返す
• computeDistanceLineRestLength
  与えられたポイントと、指定された直線間の距離に基づいてrestlengthを計算
• createAttachNormalConstraint
  対象プリミティブの法線情報を考慮したAttach拘束を作成
• pointPrimTargetPos
  UV座標と対象プリミティブのポイント位置からターゲット位置を計算
• pointPrimTargetPos
  pointPrimTargetPos関数のpattribをPとしてオーバーロードしたもの
• computePointPrimRestLength
  与えられたポイントとターゲットプリミティブ間の距離を計算
• remapPolyLineUVs
  ポリライン上のUV座標を、対象となる線分内のパラメトリック値に変換
• createStitchConstraint
  指定された条件に基づいて、2ポイント間にStitch(縫合)拘束を作成
• oppositepoint
  与えられたハーフエッジに対して、反対側にあるポリゴンのポイントを返す
• computeDihedralRestLength
  隣接する三角形の法線ベクトルから二面角を計算
• createDihedralConstraint
  与えられたジオメトリ内のポイントに接続された隣接三角形間におけるDihedral(二面角)拘束を作成
• createDihedralConstraintFromNewlyWeldedPrimitives
  新しくWeldされたプリミティブからDihedral(二面角)拘束を作成
• computeOrientRodlengths
  連続する2つのポイント間の各セグメントの回転と長さを計算
• computeBendTwistRestVector
  2ポイントのオリエンテーションを取得し、それらの間における最も近い回転を計算
• createOutEdgeConstraint
  ポリライン上のポイントに対して、Distance(距離)拘束やBend(曲げ)拘束を作成
• createStretchShearConstraint
  エッジに対して距離拘束と曲げ拘束を組み合わせたStretch-Shear拘束を作成
• createBendTwistConstraint
  エッジ上のポイントの向きが同じになるようなBend-Twist(曲げ・捻じれ)拘束を作成
• findBranchBendPoints
  ポリラインの分岐ポイントにおけるBend(曲げ)拘束を作成するために必要な2ポイントを見つける
• createBranchWeldConstraints
  枝分かれ構造を維持するためのBranchstitch拘束および曲げ・捻じれ(Bend-Twist)拘束を作成
• removeBranchWeldConstraints
  Branchstitch拘束および曲げ・捻じれ(Bend-Twist)拘束を削除
• findVolumePoints
  三角形上のポイントを体積と勾配を求めるのに適した順番にする
• computeVolume
  三角形の体積を計算
• computePressureRestLength
  Pressure(圧力)拘束に必要な体積を計算
• computePressureRestLength (volume)
  computePressureRestLength関数の計算モードをvolumeとしてオーバーロードしたもの
• createPressureConstraint
  体積を維持するPressure(圧力)拘束を作成
• createShapeMatchConstraint
  オブジェクトの形状を維持するShape Match拘束を作成
• computeTetVolumeRestLength
  四面体の体積を計算
• createTetVolumeConstraint
  四面体の体積を維持するTet Volume拘束を作成
• computeAngleRestLength
  3つの連続するポイントを受け取り、その間の角度を計算
• computeVertexCentroidDistance
  三角形の座標を受け取り、その座標から重心までの距離を計算
• computeTriangleBendRestLength
  三角形のインデックスを受け取り、その座標から重心までの距離を計算
• computeTriangleBendRestLengthFromAngle
  与えられた角度と仮定される三角形の辺の長さに基づいてrestlengthを計算
• createTriangleBendConstraints
  与えられたポイントとその隣接ポイントの間にBend(曲げ)拘束を作成
• createGlueConstraints
  指定された条件に基づいて、Glue(接着)拘束を作成
• createStrutConstraints
  指定された条件に基づいて、Strut(支柱)拘束を作成
• computeTetRestMatrix
  四面体のrestmatrix形状情報と体積を計算
• computeTetFiberRestLength
  四面体の体積と圧縮方向であるmaterialWを計算
• createTetFiberConstraint
  指定された条件に基づいて、Tet Fiber拘束を作成
• computeTriRestMatrix
  三角形のrestmatrix(形状情報)と面積を計算
• polardecomp2d
  2次元行列の極分解を計算
• computeTriAreaRestLength
  指定されたスケーリング係数に基いて三角形の面積を計算
• createTriStretchConstraint
  指定された条件に基づいて、Triangle Stretch拘束を作成
• createTetStretchConstraint
  指定された条件に基づいて、Tetrahedral Stretch拘束を作成
• orientedRestDifference
  ポイント間におけるオリエンテーションの差を計算
• logscaleStiffness
  ログスケール変換
• invlogscaleStiffness
  logscaleStiffnessの逆変換
• updateRestVectorOrient
  角度/軸の差分ベクトルをもとにrestvectorを更新し、その角度を返す
• squaredNorm2
  2x2行列におけるFrobeniusノルムの二乗を返す
• squaredNorm3
  3x3行列におけるFrobeniusノルムの二乗を返す
• updateRestMatrix2
  差分ベクトルをもとにrestvectorとrestlengthを更新しノルムの差を返す
• updateRestMatrix
  差分ベクトルをもとにrestmatrixとrestlengthを更新しノルムの差を返す
• computeRestVectorDifference
  四面体または三角形のcurrent-rest間におけるベクトル差分を計算
• computeRestMatrixDifference
  四面体のcurrent-rest間におけるmatrix差分を計算
• plasticDeformation
  塑性変形における塑性流動量を計算し、拘束の剛性値を更新する
• getSourcePoints
  拘束タイプがdistance/stretchshear/weld/ptprimであったときに、最初の配列要素を返す
• getTargetPoints
  拘束タイプがdistance/stretchshear/weld/ptprimであったときに、最初以外の配列要素を返す
• accumScaleValues
  アトリビュートを指定されたモード（最小・最大・平均・乗算）で累積させる
• falseColor
  入力値を指定した範囲に基づいてfalseを示すためのRGB値へ変換
• weldTarget
  接合の目標となるポイント番号の取得
• computeStress
  物体の変形や歪みによって生じる拘束に対するstress（応力）を計算
• maxConstraintStress
  拘束ジオメトリに含まれるstress（応力）の最大値を返す
• computeRestLengthDifference
  与えられた拘束タイプに基づいて、current-rest間における長さや角度の差分またはノルムを計算
• updateFromRest
  与えられた拘束タイプに基づいて、restlengthやrestvectorなどの値を更新
• computeStiffnessDropoff
  Stiffness Dropoffを考慮した場合の剛性値を計算
• restMetric
  塑性変形などの可視化に用いるためのrest状態におけるメトリックを計算
• constraintMetric
  指定タイプにおける拘束ジオメトリの変化量を計算
• maxConstraintMetric
  指定タイプにおける拘束ジオメトリの最大変化量を計算
• warpWeftScale
  マテリアルUV座標に基づいて、織り方向のスケーリング値を計算
• findwelds
  Weldに使用されているポイント配列を返す
• walkdist
  目標ポイントから指定プリミティブまでの最短距離やそのプリミティブにおけるUV座標を計算
• closestpttriangle
  指定したポイントpから三角形までの最も近いポイント座標を計算し、頂ポイントまたは辺の番号とバリセントリック座標を計算
• closestpttriangle
  closestpttriangleのラッパー関数であり、対応するポイント番号やハーフエッジを求める機能が追加されている
• ispolytri
  三角ポリゴンであるかどうか判別
• triwalkdist
  三角プリミティブを対象として、与えられたジオメトリ内で目標位置に最も近いポイントを取得
• extractRotation
  入力行列から回転成分を抽出
• accumVec
  カハンの加算アルゴリズムを利用してベクトルの和を計算
• accumMat3
  カハンの加算アルゴリズムを利用して3x3行列の和を計算
• computeCmAndRot
  与えられたポイント集合に対して、rest-currentにおける重心を計算し、その最適な回転を求める
• updateFromRestShapeMatch
  rest-current間における形状マッチングを行い、rest状態におけるポイントを更新
• plasticDeformationShapeMatch
  Shape Match拘束を対象として塑性変形における塑性流動量を計算し、拘束の剛性値を更新
• hasSmoothing
  与えられた拘束タイプが拘束解決時にスムージングが行えるか判別

compareIntArrays

2 つの整数配列が等しいかどうかを比較し、等しければ 1 を、そうでなければ 0 を返す

int compareIntArrays(const int a[], b[])
{
    // 配列aとbの大きさが異なっている場合
    if (len(a) != len(b))
        return 0;

    foreach(int idx; int i; a)
    {
        // 要素を比較
        if (i != b[idx])
            return 0;
    }

    return 1;
}

sortRemoveDuplicates

整数配列から重複する要素を削除し、要素をソートした新しい配列を返す

int [] sortRemoveDuplicates(const int a[])
{
    int b[];
    resize(b, len(a));  // 配列の大きさを合わせる
    int sorted[] = sort(a); // 並べ替え
    int last = 0, first = 1;
    int idx = 0;

    foreach(int c; sorted)
    {
        // 最初の要素でなく、
        // 取り出した要素が前回の取り出した要素と同じ場合
        if (!first && c == last)
            continue;

        b[idx] = c; // 配列を更新
        idx++; // インクリメント
        
        // 次のループのため更新
        first = 0;
        last = c;
    }

    resize(b, idx); // 配列の大きさを合わせる
    return b;
}

calcConstraintTopology

拘束トポロジーに関するデータIDを取得

int [] calcConstraintTopology(const int constraintgeo, pointgeo)
{
    // meta
    int topo[] = attribdataid(constraintgeo, "meta", "topology");
    append(topo, attribdataid(constraintgeo, "meta", "primitivelist"));

    // point
    append(topo, attribdataid(pointgeo, 'point', 'id'));
    append(topo, attribdataid(pointgeo, 'point', 'weld'));
    append(topo, attribdataid(pointgeo, 'point', 'branchweld'));

    return topo;
}

constraintAlias

入力された拘束タイプの文字列に対応するエイリアスを返す

string constraintAlias(const string type)
{
    // Handle alias conversion.
    if (type == "attach")
        return "pin";
    if (type == "stitch" || type == "branchstitch")
        return "distance";
    if (type == "attachnormal")
        return "distanceline";

    return type;
}

constraintHash

拘束タイプのエイリアスに基づいてハッシュ値を生成する

int constraintHash(const string type)
{
    return random_shash(constraintAlias(type));
}

isTriARAP

拘束タイプがTriARAP(三角形の角度による最小二乗法)であるかどうかを判別

1. triarap: 三角形に対するARAPアルゴリズム
2. triarapnl: 非線形変形を考慮した三角形に対するARAPアルゴリズム
3. triarapnorm: 法線情報を維持しながら変形を行う三角形に対するARAPアルゴリズム

int isTriARAP(const string type)
{
    return type == "triarap" || type == "triarapnl" || type == "triarapnorm";
}

isTetARAPVol

TetARAPにおける体積維持に使用する拘束の種類を判別

1. tetarapvol: 体積維持を考慮した四面体に対するARAPアルゴリズム
2. tetarapnlvol: 非線形変形と体積維持を考慮した四面体に対するARAPアルゴリズム
3. tetarapnormvol: 体積と法線を維持しながら変形を行う四面体に対するARAPアルゴリズム

isTetARAP

拘束タイプがTetARAP(四面体における剛体補間)であるかどうかを判別

1. tetarap: 四面体に対するARAPアルゴリズム
2. tetarapnl: 非線形変形を考慮した四面体に対するARAPアルゴリズム
3. tetarapnorm: 法線情報を維持しながら変形を行う四面体に対するARAPアルゴリズム

int isTetARAP(const string type)
{
    return type == "tetarap" || type == "tetarapnl" || type == "tetarapnorm" || isTetARAPVol(type);
}

isNonLinearARAP

拘束タイプが非線形ARAPであるかどうかを判別

int isNonLinearARAP(const string type)
{
    return type == "triarapnl" || type == "tetarapnl" || type == "tetarapnlvol";
}

isTetFiber

拘束タイプがTetFiber（特定方向に沿って圧縮することができる四面体拘束）であるかどうかを判別

1. tetfiber: 圧縮を考慮した四面体拘束
2. tetfibernorm: 圧縮と法線情報を考慮しながら変形を行う四面体拘束

int isTetFiber(const string type)
{
    return type == "tetfiber" || type == "tetfibernorm";
}

computeDistanceRestLength

2ポイント間の距離を計算

float
computeDistanceRestLength(const int geo; const int p0, p1)
{
    return distance(vector(point(geo, "P", p0)), vector(point(geo, "P", p1)));
}

createDistanceConstraint

隣接ポイントを結合し、polylineを作成

void createDistanceConstraint(const int geo; const int ptnum; const string srcgrp;
                              const int outgeo; const string outgrp)
{
    // 隣接ポイント
    int nbrs[] = neighbours(geo, ptnum);

    foreach(int n; nbrs)
    {
        // 隣接ポイントがptnumより小さく、
        // グループに属していない場合
        if (n <= ptnum || !inpointgroup(geo, srcgrp, n))
            continue;
        
        // polyline
        int prim = addprim(outgeo, "polyline", array(ptnum, n));
        // outgrp (prim) <- 1
        setprimgroup(outgeo, outgrp, prim, 1);
        // restlength (prim) <- computeDistanceRestLength
        setprimattrib(outgeo, "restlength", prim, computeDistanceRestLength(geo, ptnum, n));
    }
}

createAttachConstraint

入力ポイントから最も近いプリミティブ上の位置を計算し、Attach拘束を作成

void
createAttachConstraint(const int geo, attachgeo; const int ptnum, srcidx; const string attachgrp;
                      const int useclosestpt, useclosestprim, maxdistcheck; const float maxdist;
                      const int outgeo; const string outgrp)
{
    int targetprim = -1, targetpt = -1;
    vector targetuv, targetpos;
    float dist;

    // 座標取得
    vector pos = point(geo, "P", ptnum);

    // 近接ポイントの使用
    if (useclosestpt)
    {
        // 近接プリミティブの使用
        if (useclosestprim)
        {
            // 探索範囲の指定
            if (maxdistcheck)
                // 対象となるプリミティブの距離やパラメトリックUV情報の取得
                dist = xyzdist(attachgeo, attachgrp, pos, targetprim, targetuv, maxdist);
            else
                dist = xyzdist(attachgeo, attachgrp, pos, targetprim, targetuv);
        }
        else
        {
            // 探索範囲指定
            if (maxdistcheck)
                // 近接ポイント番号の取得
                targetpt = nearpoint(attachgeo, attachgrp, pos, maxdist);
            else
                targetpt = nearpoint(attachgeo, attachgrp, pos);
        }
    }
    else
    {
        // attachgrpに属するポイント番号の取得と
        // srcidxに対応するポイントをtargetptへ
        int tgtpts[] = expandpointgroup(attachgeo, attachgrp);
        if (srcidx < len(tgtpts))
            targetpt = tgtpts[srcidx];
    }
    
    // targetprim, targetptが更新されなかった場合
    if (targetprim < 0 && targetpt < 0)
        return;

    // sphereプリミティブの作成
    int prim = addprim(outgeo, "sphere", ptnum);
    matrix3 t = 0.01;
    // transform (prim) <- t
    setprimintrinsic(outgeo, "transform", prim, t);
    // outgrp (prim) <- 1
    setprimgroup(outgeo, outgrp, prim, 1);

    // targetprimが存在
    if (targetprim >= 0)
    {
        // target_prim(prim) <- targetprim
        setprimattrib(outgeo, "target_prim", prim, targetprim);

        // target_uv(prim) <- targetuv
        setprimattrib(outgeo, "target_uv", prim, targetuv);

        // 座標取得
        targetpos = primuv(attachgeo, "P", targetprim, targetuv);
    }
    else if (targetpt >= 0) // targetptのみが存在
    {
        // target_pt(prim) <- target_pt
        setprimattrib(outgeo, "target_pt", prim, targetpt);

        // 座標取得
        targetpos = point(attachgeo, "P", targetpt);

        // 距離取得
        dist = distance(pos, targetpos);
    }

    // restvector(prim) <- targetpos
    setprimattrib(outgeo, "restvector", prim, (vector4)targetpos);

    // restlength(prim) <- dist
    setprimattrib(outgeo, "restlength", prim, dist);
}

projectToLine

与えられた位置ベクトルを指定方向の直線上に射影

// dirは正規化されている必要がある
// Returns projection of pos onto the line, assuming dir is normalized.
vector projectToLine(const vector pos, orig, dir)
{
    return orig + dir * dot(pos - orig, dir);
}

fromTriangleSpaceXform

3ポイントで構成される平面からワールド空間への変換行列を返す
pbd_constraints.cl の toTriangleSpaceXformと一致させる必要がある

// Returns the matrix to transform from triangle to world space.
// This MUST match toTriangleSpaceXform in pbd_constraints.cl.
matrix3 fromTriangleSpaceXform(const vector p0, p1, p2)
{
    // 2ポイントからベクトルを作成
    vector e0 = p1 - p0; 
    vector e1 = p2 - p0;
    
    // 外積を用いてe0とe1と直行するベクトルを作成しz軸として使用
    // x軸、y軸は、z軸と同様の方法で作成
    vector e2 = cross(e1, e0); 
    vector z  = normalize(e2); // 正規化
    vector y = normalize(cross(e0, z));
    vector x = cross(y, z);
    return set(x, y, z);
}

fromTriangleSpace

3ポイントで構成される平面からワールド空間への回転ベクトルを返す

// 入力プリミティブ最初の3ポイントから構成される面を用いてワールド座標系への変換行列を作成
// 変換行列を用いてdirベクトルを回転

// Rotate the vector from triangle space using the first three points of the polygon
// specified in prim.  Invert from world space to triangle if directed.
vector
fromTriangleSpace(const string path; const int prim; const vector dir; const int invert)
{
    // ポリゴンを構成しているポイント3ポイントを取得し、ワールド座標系への変換行列を作成
    int pts[] = primpoints(path, prim);
    vector p0 = point(path, "P", pts[0]);
    vector p1 = point(path, "P", pts[1]);
    vector p2 = point(path, "P", pts[2]);

    // 変換行列
    matrix3 xform = fromTriangleSpaceXform(p0, p1, p2);

    // invertが有効の場合は、作成した変換行列を転置
    if (invert)
        return dir * transpose(xform);
    return dir * xform;
}

computeDistanceLineRestLength

与えられたポイントと指定された直線間の距離に基づいてrestlengthを計算

float
computeDistanceLineRestLength(const int geo; const int pt0; const vector orig, indir;
                              const string path; const int prim)
{
    vector pos = point(geo, "P", pt0);

    // Transform from triangle space if needed.
    // 三角面からなる空間の変換が必要な場合は、
    // pathとプリミティブ番号を指定し
    // fromTriangleSpaceを用いてdirベクトルを変換
    vector dir = indir;
    if (path != "" && prim >=0)
        dir = fromTriangleSpace(path, prim, dir, 0);

    // posと直線上に射影したポイントとの距離を返す
    return distance(pos, projectToLine(pos, orig, dir));
}

createAttachNormalConstraint

ターゲットジオメトリの法線を考慮したAttach拘束を作成

// Create an AttachNormal constraint that stores the line direction in restdir
// in triangle space.
// NOTE: Possibly more convenient would be to store restdir in world space, but also
// store a 3x3 matrix that would transform it back to triangle space.  Then on each
// timestep we would update restdir from the current target geometry with:
// v@restdir *= 3@restmatrix * curmatrix;
// where curmatrix is the current triangle-to-world transform, then:
// 3@restmatrix = invert(curmatrix);
// But this would require storing an extra matrix per constraint.
int
createAttachNormalConstraint(const int geo; const int ptnum;
                            const string targetpath;
                            const int targetprim; const vector targetuv;
                            const int outgeo; const string outgrp)
{
    vector pos = point(geo, "P", ptnum);
    // Orig is point on target attach geometry.
    vector orig = primuv(targetpath, "P", targetprim, targetuv);
    // Dir is ray from orig to point position.  This will be close to the normal
    // for convex geometry, but might vary quite a bit for non-convex. We use
    // this instead of the normal so the rest state is stable in the latter case.
    vector dir = normalize(pos - orig);

    // Transform to triangle space; the solver will transform back to world space
    // based on current target geometry point positions.
    // Only polygons.
    // ポリゴンのみ
    if (primintrinsic(targetpath, "typeid", targetprim) != 1) // typeid: 1 -> Poly
        return -1;
    // With at least three points to form a basis.
    // 少なくとも3ポイントから構成
    if (len(primpoints(targetpath, targetprim)) < 3)
        return -1;

    // dirベクトルの回転（invert）
    dir = fromTriangleSpace(targetpath, targetprim, dir, /*invert=*/ 1);

    // sphereプリミティブの作成
    int prim = addprim(outgeo, "sphere", ptnum);
    matrix3 t = 0.01;

    // transform (prim) <- t
    setprimintrinsic(outgeo, "transform", prim, t);
    // target_prim (prim) <- targetprim
    setprimattrib(outgeo, "target_prim", prim, targetprim);
    // target_uv (prim) <- targetuv
    setprimattrib(outgeo, "target_uv", prim, targetuv);
    // target_path (prim) <- targetpath
    setprimattrib(outgeo, "target_path", prim, targetpath);

    // Restlength is always zero to constraint directly to line to closest point.
    // 接ポイントは0
    // restlength (prim) <- 0.0 
    setprimattrib(outgeo, "restlength", prim, 0.0f);
    // restvector (prim) <- orig 接地ポイント
    setprimattrib(outgeo, "restvector", prim, (vector4) orig);
    // restdir (prim) <- dir 更新した方向ベクトル
    setprimattrib(outgeo, "restdir", prim, dir);
    // outgrp (prim) <- 1
    setprimgroup(outgeo, outgrp, prim, 1);

    return prim; // プリミティブ番号
}

pointPrimTargetPos

UV座標と対象プリミティブのポイント位置からターゲット位置を計算

// Compute the target position directly from the uv coords and the primitive's
// point positions.  (We can't use primuv because we don't know the primitive
// number, and don't want to try to track it.)
vector pointPrimTargetPos(const int geo; const int pts[]; const float u, v; const string pattrib)
{
    vector p0 = point(geo, pattrib, pts[0]);
    vector p1 = point(geo, pattrib, pts[1]);
    // Only line segments, triangles and quads.
    int npts = len(pts);

    if (npts == 2)  // ライン
    {
        float w0 = (1 - u);
        float w1 = u;
        return w0 * p0 + w1 * p1;
    }
    else if (npts == 3)  // 三角ポリゴン
    {
        vector p2 = point(geo, pattrib, pts[2]);
        float w0 = (1 - u - v);
        float w1 = u;
        float w2 = v;
        return w0 * p0 + w1 * p1 + w2 * p2;
    }
    else if (npts == 4) // 四角ポリゴン
    {
        vector p2 = point(geo, pattrib, pts[2]);
        vector p3 = point(geo, pattrib, pts[3]);
        float u1 = 1 - u;
        float v1 = 1 - v;
        float w0 = (u1 * v1);
        float w1 = (u1 * v);
        float w2 = (u * v);
        float w3 = (u * v1);
        return w0 * p0 + w1 * p1 + w2 * p2 + w3 * p3;
    }
    return 0;
}

pointPrimTargetPos

pointPrimTargetPos関数のpattribを"P" としてオーバーロードしたもの

// オーバーロード
vector pointPrimTargetPos(const int geo; const int pts[]; const float u, v)
{
    return pointPrimTargetPos(geo, pts, u, v, "P");
}

computePointPrimRestLength

与えられたポイントとターゲットプリミティブ間の距離を計算

// pointPrimTargetPosで求めた接地ポイントとの距離を計算

float computePointPrimRestLength(const int geo; const int pt; const int tgtpts[]; const vector4 restvector)
{
    int npts = len(tgtpts);
    // Only line segments, triangles, and quads.
    // ライン、三角ポリゴン、四角ポリゴンのみ
    if (npts < 2 || npts > 4)
        return 0;

    // 入力ジオメトリの座標
    vector pos = point(geo, "P", pt);
    // 接地ポイントの座標
    vector tgtpos = pointPrimTargetPos(geo, tgtpts, restvector.x, restvector.y);
    // 距離計算
    return distance(pos, tgtpos);
}

remapPolyLineUVs

ポリライン上のUV座標を、対象となる線分内のパラメトリック値に変換

// Map from uv's over an entire polyline to a single parametric value along the
// line segment in which the uv point falls.  Targetpts should contain all the line
// points as input, but will have the segment points as output.
void
remapPolyLineUVs(const string geo; const int targetprim; int targetpts[]; vector targetuv)
{
    vector2 uv = set(targetuv.x, 0);

    // Convert to segment based mapping
    // ポリライン全体のUV座標を線分ベースのマッピングに変換
    // 正規化された値からセグメントへ (0-1 range -> actual point number)
    uv = primuvconvert(geo, uv, targetprim, PRIMUV_UNIT_TO_REAL);
    int pt0 = floor(uv.x);

    // If on last point we need to back up 1.
    if (pt0 >= len(targetpts) - 1)
        pt0 -= 1;

    uv.x -= pt0;

    targetpts = array(targetpts[pt0], targetpts[pt0 + 1]);
    targetuv = set(uv.x, 0, 0);
}

createStitchConstraint

指定された条件に基づいて、2ポイント間にStitch(縫合)拘束を作成

void
createStitchConstraint(const int geo; const int ptnum, srcidx; const string targetgrp;
                      const int useclosestpt, useclosestprim, maxdistcheck; const float maxdist;
                      const int outgeo; const string outgrp)
{
    int targetprim = -1, targetpt = -1;
    vector targetuv, targetpos;
    float dist;
    int tgtpts[];

    // 座標取得
    vector pos = point(geo, "P", ptnum);

    // 近接ポイントの使用
    if (useclosestpt)
    {
        // 近接プリミティブの使用
        if (useclosestprim)
        {
            // 探索範囲の指定
            if (maxdistcheck)
                // 対象プリミティブの距離やパラメトリックUV情報の取得
                dist = xyzdist(geo, targetgrp, pos, targetprim, targetuv, maxdist);
            else
                dist = xyzdist(geo, targetgrp, pos, targetprim, targetuv);
        
            // targetprimの値が更新されない場合
            if (targetprim < 0)
                return;

            // Poly
            if (primintrinsic(geo, "typeid", targetprim) != 1)
                return;

            // 対象プリミティブ上のポイント取得
            tgtpts = primpoints(geo, targetprim);

            // polyline, remap to line segment.
            if (!primintrinsic(geo, "closed", targetprim))

                // UV座標を線分ベース（point number）へ変換
                // PRIMUV_UNIT_TO_REAL
                remapPolyLineUVs(sprintf("opinput:%d", geo), targetprim, tgtpts, targetuv);

            // Line segment, triangle, or quad polygons only.
            int npts = len(tgtpts);
            if (npts < 2 || npts > 4)
                return;
        }
        else
        {
            // 探索範囲の指定
            if (maxdistcheck)
                // 近接ポイント番号の取得
                targetpt = nearpoint(geo, targetgrp, pos, maxdist);
            else
                targetpt = nearpoint(geo, targetgrp, pos);
        }
    }
    else
    {
        // targetgrpに属するポイント番号の取得
        tgtpts = expandpointgroup(geo, targetgrp);

        // srcidxに対応するポイントをtargetptへ
        if (srcidx < len(tgtpts))
            targetpt = tgtpts[srcidx];
    }

    // targetprim, targetptが更新されなかった場合
    if (targetprim < 0 && targetpt < 0)
        return;

    int prim = -1;

    // targetprimが存在
    if (targetprim >= 0)
    {
        int pts[];
        append(pts, ptnum);
        append(pts, tgtpts);

        // UV座標のセット
        vector4 restvec = set(targetuv.x, targetuv.y, 0, 0);

        // polylineの作成
        prim = addprim(outgeo, "polyline", pts);

        // すでにxyzdistで距離は求めているが再度計算
        // Just to ensure these match exactly, even though we already know dist.
        dist = computePointPrimRestLength(geo, ptnum, tgtpts, restvec);

        // restlength (prim) <- dist
        setprimattrib(outgeo, "restlength", prim, dist);

        // restvector (prim) <- restvec 
        setprimattrib(outgeo, "restvector", prim, restvec);

        // type: ptprim
        setprimattrib(outgeo, "type", prim, "ptprim");
    }
    else if (targetpt >= 0) // targetptが存在
    {
        // polylineの作成
        prim = addprim(outgeo, "polyline", array(ptnum, targetpt));

        // restlength <- computeDistanceRestLength
        setprimattrib(outgeo, "restlength", prim, computeDistanceRestLength(geo, ptnum, targetpt));

        // type: stitch
        setprimattrib(outgeo, "type", prim, "stitch");
    }

    // outgrp (prim) <- 1
    setprimgroup(outgeo, outgrp, prim, 1);
}

oppositepoint

与えられたハーフエッジに対して、反対側にあるポリゴンのポイントを返す

int oppositepoint(const int geo; const int hedge)
{
    // hedgeの後に続くハーフエッジを取得し、ゴールポイント番号を取得
    return hedge_dstpoint(geo, hedge_next(geo, hedge));
}

computeDihedralRestLength

隣接する三角形の法線ベクトルから二面角を計算

// 指定した4つのポイントが形成する二面角の初期角度を計算
int computeDihedralRestLength(const int geo;
                              const int pt0, pt1, pt2, pt3;
                              float restlength)
{
    vector p0 = point(geo, "P", pt0);
    vector p1 = point(geo, "P", pt1);
    vector p2 = point(geo, "P", pt2);
    vector p3 = point(geo, "P", pt3);

    vector e = p3 - p2;
    float elen = length(e);

    // 長さが短すぎるものは0
    if (elen < 1e-6)
        return 0;

    float invElen = 1 / elen;

    // Find initial rest angle.
    // 平面の法線ベクトルを計算
    vector n1 = cross(p3 - p0, p2 - p0);
    vector n2 = cross(p2 - p1, p3 - p1);

    // 2つの法線ベクトルから二面角を計算
    float d = dot(normalize(n1), normalize(n2));
    d = clamp(d, -1, 1);
    float phi = acos(d);

    // We want to xpress phi a -PI..PI
    // 二面角は-PiからPiまで
    if (dot(cross(n1, n2), e) < 0)
        phi = -phi;

    // restlength <- 角度に変換した二面角
    restlength = degrees(phi);
    return 1;
}

createDihedralConstraint

与えられたジオメトリ内のポイントに接続された隣接三角形間におけるDihedral(二面角)拘束を作成

void
createDihedralConstraint(const int geo; const int ptnum; const string srcgrp;
                        const int outgeo; const string outgrp)
{

    // 指定したポイントに属するプリミティブの配列を取得
    int prims[] = pointprims(geo, ptnum);

    // 取得したプリミティブごとに処理
    foreach(int prim; prims)
    {
        // プリミティブに属するポイントの配列を取得
        int primpts[] = primpoints(geo, prim);

        // Process triangles where this point is minimum point number,
        // and all points are in src group.
        // そのプリミティブが三角形であり、
        // 指定ポイントが最小のポイント番号であり、
        // 他の2つのポイントがsrcgrpというグループに属している
        if (len(primpts) != 3)
            continue;
        if (ptnum != min(primpts))
            continue;
        if (!inpointgroup(geo, srcgrp, primpts[1]) ||
            !inpointgroup(geo, srcgrp, primpts[2]))
            continue;

        // そのプリミティブに属する最初のハーフエッジ番号を取得
        int starthedge = primhedge(geo, prim);

        // ハーフエッジ番号が負ならreturn
        if (starthedge < 0)
            return;

        // Ignore open curves as the hedge function won't give us
        // the setup we expect with them.
        // 閉じているか
        if (primintrinsic(geo, 'closed', prim) == 0)
            return;

        int h = starthedge;
        int p = prim;
        int hasgrp = strlen(outgrp) > 0;

        while (1)
        {
            // Giving half edge h, add a bend polygon.
            // 隣接するハーフエッジを取得
            int oh = hedge_nextequiv(geo, h);

           // Skip boundary edges, invalid hedges, and non-manifold
           // edges
            if (h != oh && oh >= 0 && h == hedge_nextequiv(geo, oh))
            {
                int op = hedge_prim(geo, oh);
                // Always build in ascending direction.
                if (op >= 0 && p < op)
                {
                    // ハーフエッジと隣接するハーフエッジが形成する4つのポイント番号を取得
                    int pt0 = oppositepoint(geo, h);
                    int pt1 = oppositepoint(geo, oh);
                    int pt2 = hedge_srcpoint(geo, h);
                    int pt3 = hedge_dstpoint(geo, h);

                   // ポイントがすべてsrcgrpに属している
                    if (inpointgroup(geo, srcgrp, pt0) &&
                        inpointgroup(geo, srcgrp, pt1) &&
                        inpointgroup(geo, srcgrp, pt2) &&
                        inpointgroup(geo, srcgrp, pt3))
                    {
                        // 二面角の計算
                        float restlength;
                        if (computeDihedralRestLength(geo, pt0, pt1, pt2, pt3,
                                                      restlength))
                        {
                            // polyline
                            int newprim = addprim(outgeo, 'polyline', array(pt0, pt1, pt2, pt3));
                            // restlength (prim) <- restlength
                            setprimattrib(outgeo, "restlength", newprim, restlength);
                            
                            // グループが存在
                            if (hasgrp)
                                setprimgroup(outgeo, outgrp, newprim, 1);
                        }
                    }
                }
            }

            // 次のハーフエッジ取得
            int nh = hedge_next(geo, h);

            // Stop the loop when we complete the polygon
            // Closed polygons won't complete.
            // nhが最初のハーフエッジと同じか、無効ならループから抜ける
            if (nh == starthedge || nh < 0)
                break;

            // 次のループへ
            h = nh;
        }
    }
}

createDihedralConstraintFromNewlyWeldedPrimitives

新しくWeldされたプリミティブからDihedral(二面角)拘束を作成

// この関数の目的は、ソフトボディシミュレーションに使用できる、
// 新しく溶接されたプリミティブの間の二面体拘束を作成することです
// この拘束は、複雑なオブジェクトが外力を受けて変形する際に、その形状を保持するのに役立ちます

void
createDihedralConstraintFromNewlyWeldedPrimitives(const int geo; const int oldgeo; const int prim;
                                                  const int outgeo; const string outgrp)
{
    // ポイント番号の取得
    int primpts[] = primpoints(geo, prim);

    // プリミティブが3つのポイントを持っていた場合は何もしない
    if (len(primpts) != 3)
      return;

    // ハーフエッジの取得
    int starthedge = primhedge(geo, prim);

    // ハーフエッジがなければ何もしない
    if (starthedge < 0)
      return;

    // Ignore open curves as the hedge function won't give us
    // the setup we expect with them.
    // 閉じているか
    if (primintrinsic(geo, 'closed', prim) == 0)
      return;

    int h = starthedge;
    int p = prim;
    int hasgrp = strlen(outgrp) > 0; // グループを持っていた場合

    // hedge_next を用いて、プリミティブ上すべてのハーフエッジをループし、
    // 各ハーフエッジ (h) に対して、別のプリミティブ (op) 上と同等の半辺 (oh) を hedge_nextequiv 関数で探索
    while (1)
    {
      // Giving half edge h, add a bend polygon.
      // starthedgeに対応する次のハーフエッジを取得
      int oh = hedge_nextequiv(geo, h);

      // Skip boundary edges, invalid hedges, and non-manifold
      // edges
      // 境界のエッジ、無効なハーフエッジ、非多様体エッジはスキップ
      if (h != oh && oh >= 0 && h == hedge_nextequiv(geo, oh))
      {
        // ハーフエッジ (oh)に含まれるプリミティブ (op) を取得
        int op = hedge_prim(geo, oh);

        // Always build in ascending direction.
        // 常に方向は昇順
        // プリミティブインデックス (p) が他のプリミティブインデックス (op) よりも小さい場合のみ、
        // ベンドポリゴン (4ポイントを持つポリライン) を構築
        if (op >= 0 && p < op)
        **{
          // 二面角を計算するためのポイントを取得**
          int pt0 **=** oppositepoint(geo, h);
          int pt1 = oppositepoint(geo, oh);
          int pt2 = hedge_srcpoint(geo, h);
          int pt3 = hedge_dstpoint(geo, h);

          // See if we became welded.  Because hedges are
          // linear vertex numbers, these will match in the
          // two and we can find the corresponding half edge
          // in the other geometry directly.
          // So, if our h & oh don't have matching points
          // in the old geo, we know we are a newly welded edge.
          // ハーフエッジが新しくウェルドされたかどうかは、
          // 両方のジオメトリ（geoとoldgeo）のソースポイントとデスティネーションポイントを比較
          // もし一致しない場合は、何らかの操作（ブール演算など）によりマージされたことを意味
          if (hedge_srcpoint(oldgeo, h) != hedge_dstpoint(oldgeo, oh) ||
              hedge_dstpoint(oldgeo, h) != hedge_srcpoint(oldgeo, oh))
          {
              float restlength;

              // computeDihedralRestLength というカスタム関数を使用して、ベンドポリゴンの二面角を計算
              // この関数は、変形時に元の形状を維持するために、ベンドポリゴンをどれだけ伸縮させるべきかを計算
              if (computeDihedralRestLength(geo, pt0, pt1, pt2, pt3, restlength))
              {
                  // 出力ジオメトリ（outgeo）にタイプ 'polyline' の新しいプリミティブを追加し、
                  // その属性 'restlength' を計算されたものと同じに設定
                  int newprim = addprim(outgeo, 'polyline', array(pt0, pt1, pt2, pt3));
                  setprimattrib(outgeo, "restlength", newprim, restlength);

                  // 出力グループ名 (outgrp) が与えられている場合は、
                  // この新しいプリミティブもそのグループに追加
                  if (hasgrp)
                      setprimgroup(outgeo, outgrp, newprim, 1);
              }
            }
          }
        }
  
      // 次のハーフエッジを取得
      int nh = hedge_next(geo, h);
    
      // Stop the loop when we complete the polygon
      // Closed polygons won't complete.
      // プリミティブ上のハーフエッジを一周したらループから抜ける
      if (nh == starthedge || nh < 0)
        break;

      // 次のstarthedgeとする
      h = nh;
    }
}

computeOrientRodlengths

連続する2つのポイント間の各セグメントの回転と長さを計算

// inertia: パーティクルの回転ヘアー拘束の抵抗力

// Compute point orients and store the rodlengths in the inertia attribute.
// vellum constraint: compute_orient_rodlengths
// computeOrientRodlengths(0, @primnum, "__constraintsrc", geoself())

void
computeOrientRodlengths(const int geo; const int primnum; const string srcgrp;
                        const int outgeo)
{ 
    // Ignore anything but open polylines.
    // 開いているポリラインのみ処理対象とする
    if (primintrinsic(geo, "typename", primnum) != "Poly" ||
        primintrinsic(geo, "closed", primnum) == 1)
        return;

    // Only test group if it's not all points.
    // グループを持っているか
    int hasgrp = npointsgroup(geo, srcgrp) < npoints(geo);

    // Check if the points have an incoming orient attribute. If os, we assume
    // that it provides a stable basis in which to calculate our rod-aligned orients.
    // ジオメトリがorientアトリビュートを持っているかどうか
    int hasporient = haspointattrib(geo, "orient");
    vector from = {0, 0, 1};

    // Iterate over pts in vertex order.
    // 頂ポイント順にイテレーション
    int pts[] = primpoints(geo, primnum); // 各プリミティブ上のポイントを取得
    int npts = len(pts); // 頂ポイントの総数
    vector4 orients[];  
    float rodlens[];
    resize(orients, npts - 1);  // 配列をリサイズ
    resize(rodlens, npts - 1);　// 配列をリサイズ

    int lastpt = pts[npts - 1];  // 最後のポイントインデックス
    int loop = 0;

  　// 連続する2ポイント間の各セグメントについて回転と長さを計算
    for(int i=0; i < npts - 1; i++)
    {
        // 
        vector d = point(geo, "P", pts[i + 1]) - point(geo, "P", pts[i]);
        vector to = normalize(d);

        // ポイントがorientを持っていた場合は向きを変換
        if (hasporient)
        {
            // Transform vectors back to rest orientation.
            vector4 porient = point(geo, "orient", pts[i]);
            to = qrotate(qinvert(porient), to);
        }
  
        // {0, 0, 1} から toベクトルへのクォータニオンを作成
        vector4 dq = dihedral(from, to);

        // 始ポイント以外は、前のorient分に傾ける
        if (i == 0)
            orients[i] = dq;
        else
            orients[i] = qmultiply(dq, orients[i-1]);
  
        // rodlens <- 2ポイント間の距離
        rodlens[i] = length(d);

        from = to;  // 次のループで使用するベクトル方向

        // Check if this is a loop.
        // ポリラインがループしているかどうか
        loop |= (pts[i] == lastpt); 
    }

    for(int i=0; i < npts - 1; i++)
    {
        // orientを持っていた場合は、ワールド空間に戻す
        if (hasporient)
        {
            // Rotate new orientation back to current world orientation.
            vector4 porient = point(geo, "orient", pts[i]);
            orients[i] = qmultiply(porient, orients[i]);
        }

        // Make sure to do the above rotation in the orients array
        // before possibly skipping this point based on group membership,
        // otherwise if we then grab the rotation for the last point from
        // the previous vertex, it could be incorrect if the previous
        // vertex wasn't in srcgrp. (Exiting too early caused bug 97655.)
        // 最後のポイントの回転を前の頂ポイントから取得する際に，
        // 前の頂ポイントが srcgrp に含まれていない場合に不正確な回転が発生する
        if (hasgrp && !inpointgroup(geo, srcgrp, pts[i]))
            continue;

        // Set new orientation, possibly overwriting input orientation.
        // orient <- orients[i] of pts[i]
        setpointattrib(geoself(), "orient", pts[i], orients[i]);

        float inertia = point(geo, "inertia", pts[i]);

        // Don't overwrite pinned inertia
        // pin固定されている場合（inertiaが0）、上書きしない
        if (inertia == 0.0)
            continue;

        // 距離をinertiaとして設定
        // inertia <- rodlens[i] of pts[i]
        setpointattrib(geoself(), "inertia", pts[i], rodlens[i]);
    }

    // Set inertia and orient for last point if not already seen as part of a loop.
    // ポイント総数が1以上で、最終ポイントがsrcgrpに属しており、かつポリラインがループしていない
    if (npts > 1 && inpointgroup(geo, srcgrp, lastpt) && !loop)
    {
        float inertia = point(geo, "inertia", lastpt);
        // Don't overwrite pinned inertia
        if (inertia != 0.0)
            setpointattrib(geoself(), "inertia", lastpt, rodlens[npts - 2]);
        // Copy previous orient so final bend/twist rest relative orient is
        // identity, which reduces flipping when final vertex orientation is
        // otherwise unconstrained.
        setpointattrib(geoself(), "orient", lastpt, orients[npts - 2]);
    }
}

computeBendTwistRestVector

2ポイントのオリエンテーションを取得し、それらの間における最も近い回転を計算

computeBendTwistRestVector(const int geo; const int pt0, pt1)
{
    // Discrete Darbeaux vector is just closest rotation from
    // first orientation to second.

    // pt0の共役
    vector4 q0conj = point(geo, "orient", pt0) * {-1, -1, -1, 1};
    vector4 q1 = point(geo, "orient", pt1);

    // q0からq1への回転
    vector4 restDarbeaux = qmultiply(q0conj, q1);

    // quaternion: (x, y, z, w)
    // wの範囲が-1から1であることを利用し単位クォータニオンを加算・減算
    vector4 omegaplus = restDarbeaux + {0, 0, 0, 1}; // 正の方向
    vector4 omegaminus = restDarbeaux - {0, 0, 0, 1}; // 負の方向

    // もしomegaminusのdot積がomegaplusのdot積より大きい場合、
    // 最も近い回転ではないので-1を掛けて反転
    if (dot(omegaminus, omegaminus) > dot(omegaplus, omegaplus))
        restDarbeaux *= -1;

    return restDarbeaux;
}

createOutEdgeConstraint

ポリライン上のポイントに対して、Distance(距離)拘束やBend(曲げ)拘束を作成

// Create a single constraint of the specified type along the (first and only) output edge from
// the provided point, assuming we're dealing with points on polylines.
// ポリライン上のポイントを処理対象として想定し、
// Constraintを作成（Type: stretchshear or bendtwist）
void
createOutEdgeConstraint(const int geo; const int ptnum; const string srcgrp; const string type;
                        const int outgeo; const string outgrp)
{
    // 隣接しているポイント番号の取得
    int nbrs[] = neighbours(geo, ptnum);

    foreach(int n; nbrs)
    {
        // Only create if the in-group neighbor is the destination of an edge.
        // srcgrpに所属していない、またはハーフエッジが無効
        if (!inpointgroup(geo, srcgrp, n) || pointhedge(geo, ptnum, n) < 0)
            continue;

        // polylineの作成
        int prim = addprim(outgeo, "polyline", array(ptnum, n));

        // stretchshear
        if (type == "stretchshear")
            // 長さを保持
            // restlength (prim) <- computeDistanceRestLength
            setprimattrib(outgeo, "restlength", prim, computeDistanceRestLength(geo, ptnum, n));
        else // bendtwist
            // 角度を保持
            // restvector (prim) <- computeBendTwistRestVector
            setprimattrib(outgeo, "restvector", prim, computeBendTwistRestVector(geo, ptnum, n));
        
        // outgrp <- 1
        setprimgroup(outgeo, outgrp, prim, 1);

        // Only one constraint per point.
        return;
    }
}

createStretchShearConstraint

エッジに対して距離拘束と曲げ拘束を組み合わせたStretch-Shear拘束を作成

void
createStretchShearConstraint(const int geo; const int ptnum; const string srcgrp;
                            const int outgeo; const string outgrp)
{
    createOutEdgeConstraint(geo, ptnum, srcgrp, "stretchshear", outgeo, outgrp);
}

createBendTwistConstraint

エッジ上のポイントの向きが同じになるようなBend-Twist(曲げ・捻じれ)拘束を作成

// Constraintを作成（Type: bendtwist）

void
createBendTwistConstraint(const int geo; const int ptnum; const string srcgrp;
                          const int outgeo; const string outgrp)
{
    createOutEdgeConstraint(geo, ptnum, srcgrp, "bendtwist", outgeo, outgrp);
}

findBranchBendPoints

ポリラインの分岐ポイントにおけるBend(曲げ)拘束を作成するために必要な2ポイントを見つける

// bend or twist constraintを作成するためのポイントを探す

int
findBranchBendPoints(const int geo; const int ptnum; const int weld;
                    int pt0, pt1)
{
    // weldの隣接ポイントを取得
    int nbrs[] = neighbours(geo, weld);
    
    foreach(int n; nbrs)
    {
        // Check if there's an edge pointing from n to the weld.
        // 隣接ポイントを含むハーフエッジがある場合
        if (pointhedge(geo, n, weld) >= 0)
        {
            // We need to create a constraint from n to this point also.
            pt0 = n;
            pt1 = ptnum;
            return 1;
        }
    }
    // If only one neighbor and there's no edge with weld as the destination,
    // then punt and make the constraint between the weld and ptnum
    // 隣接ポイントが1つだけで、目的のエッジが見つからなかった場合は、
    // pt0にweldを代入し、pt1にptnumを代入する
    if (len(nbrs) == 1)
    {
        pt0 = weld;
        pt1 = ptnum;
        return 1;
    }
    return 0;
}

createBranchWeldConstraints

枝分かれ構造を維持するためのBranchstitch拘束および曲げ・捻じれ(Bend-Twist)拘束を作成

void
createBranchWeldConstraints(const int geo; const int ptnum;
                            const string srcgrp; const float maxbranchangle;
                            const int outgeo; const string outbendgrp)
{
    // 結合ポイント
    int weld = point(geo, "branchweld", ptnum);

    // weldが無効、またはsrcgrpに属していない場合はreturn
    if (weld < 0 || !inpointgroup(geo, srcgrp, weld))
        return;

    // Stitch points together with stiff distance constraint.
    // TODO - do we need to parameterize this somehow?
    // プリミティブの作成（type: branchstitch）
    int cprim = addprim(outgeo, "polyline", array(weld, ptnum));
    setprimattrib(outgeo, "restlength", cprim, 0); // restlength <- 0
    setprimattrib(outgeo, "stiffness", cprim, 1e20);  // stiffness <- 1e20
    setprimattrib(outgeo, "dampingratio", cprim, 0.001); // dampingratio <- 0.001
    setprimattrib(outgeo, "type", cprim, "branchstitch"); // type <- branchstitch

    // Find the appropriate points for creating a bend/twist
    // constraint on the geometry.
    // 捻じれの回転角を計算するために
    // 曲げ・捻じれ拘束に適しているpt0, pt1を探す
    int pt0, pt1;
    if (!findBranchBendPoints(geo, ptnum, weld, pt0, pt1))
        return;

    // Compute Darbeaux vector for rotation.
    // 回転角の計算
    vector4 restvec = computeBendTwistRestVector(geo, pt0, pt1);

    // Skip creation if greater than max branching angle.
    // ラジアンへ変換
    float maxang = radians(maxbranchangle);

    // Compare against axis/angle represenation of quaternion.
    // 回転量がmaxang以上であれば、return
    if (length(qconvert(restvec)) > maxang)
        return;

    // polylineの作成
    cprim = addprim(outgeo, "polyline", array(pt0, pt1));

    // restvector (vector4) <- restvec
    setprimattrib(outgeo, "restvector", cprim, restvec);

    // outbendgrp (prim) <- 1
    setprimgroup(outgeo, outbendgrp, cprim, 1);
}

removeBranchWeldConstraints

Branchstitch拘束および曲げ・捻じれ(Bend-Twist)拘束を削除

// Called with the primitive representing the branchweld stitch/distance constraint.
void
removeBranchWeldConstraints(const int congeo, ptgeo; const int primnum; const int outgeo)
{
    // プリミティブ上のポイント取得
    int pts[] = primpoints(congeo, primnum);

    // If pt is not still branchwelded to weld, we need to delete constraints.
    int weld = pts[0];
    int pt = pts[1];
    int bw = point(ptgeo, "branchweld", pt);

    // branchweldに対応するポイント番号を取得
    if (bw >= 0)
        bw = idtopoint(ptgeo, bw);

    // Still welded, nothing to do.
    // まだweldされていた場合は、return
    if (weld == bw)
        return;

    // Remove this constraint stitching them together.
    // weldされていなかったものに関しては、プリミティブを削除する
    removeprim(outgeo, primnum, 0);

    // Look for bend/twist constraints joining the bend points,
    // assume it was created for branchwelding and remove.
    
    int bendpt0, bendpt1;
    if (!findBranchBendPoints(ptgeo, pt, weld, bendpt0, bendpt1))
        return;

    // ptに対応するプリミティブの取得
    int prims[] = pointprims(congeo, pt);

    // typeがbendtwistであるもの（曲げ/捻じれ拘束）に関しても削除する
    foreach(int prim; prims)
    {
        string type = prim(congeo, "type", prim);

        // typeがbendtwistのみ
        if (type != "bendtwist")
            continue;

        // // プリミティブ上のポイント取得
        int primpts[] = primpoints(congeo, prim);

        if (bendpt0 == primpts[0] && bendpt1 == primpts[1])
            removeprim(outgeo, prim, 0);
    }
}

findVolumePoints

三角形上のポイントを体積と勾配を求めるのに適した順番にする

// Returns array of neighboring points on triangles that can
// be iterated through to calculate volume and gradients
// using cross products.
int []
findVolumePoints(const int geo; const int ptnum)
{
    // プリミティブ取得
    int prims[] = pointprims(geo, ptnum);
    int volpts[];

    foreach(int prim; prims)
    {
        // ポイント番号の取得
        int pts[] = primpoints(geo, prim);

        // Only triangles.
        // 三角形のみ
        if (len(pts) != 3)
            continue;

        // ptnumにマッチするポイント番号を取得
        int ppos = find(pts, ptnum);
        vector n = 0;

        // ppos以外のポイントをvolptsに追加
        if (ppos == 0)
            append(volpts, array(pts[2], pts[1]));
        if (ppos == 1)
            append(volpts, array(pts[0], pts[2]));
        if (ppos == 2)
            append(volpts, array(pts[1], pts[0]));
    }
    return volpts;
}

computeVolume

三角形の体積を計算

// Compute the volume for all triangles specified in the
// volumepts array for each point.  The volume is only
// counted if the current points is the lowest index of the tri.
float
computeVolume(const int geo; const int ptnum; const int volumepts[])
{
    int ntris = len(volumepts) / 2; // 三角形の数

    float volume = 0;
    vector pos = point(geo, "P", ptnum);

    for(int i = 0; i < ntris; i++)
    {
        int pt1 = volumepts[i * 2];
        int pt2 = volumepts[i * 2 + 1];

        // ポイントが最もインデックスの小さいポイントである場合にのみ体積を加算
        if (ptnum < pt1 && ptnum < pt2)
            volume += dot(pos,
                          cross(point(geo, "P", pt1), point(geo, "P", pt2)));
    }
    return volume / 6; // 体積
}

computePressureRestLength

Pressure(圧力)拘束に必要な体積を計算

float
computePressureRestLength(const int geo; const int inpts[]; const string mode)
{
    float restvol = 0;
    // Remove any duplicates that can arise from welding.
    // 重複しているポイントを削除しソート
    int pts[] = sortRemoveDuplicates(inpts);
    int volumepts[];

    foreach(int pt; pts)
    {
        if (mode == "volume")
        {
            // ポイントがvolume attributeを持っている場合はそのまま加算
            restvol += point(geo, "volume", pt);
        }
        else
        {
            // volumepts attributeを持っている場合は体積計算に利用
            // 持っていない場合はfindVolumePoints関数を用いてvolumeptsを見つける
            if (mode == "volumepts")
                volumepts = point(geo, "volumepts", pt);
            else // Use topology.
                volumepts = findVolumePoints(geo, pt);

            // 体積計算
            restvol += computeVolume(geo, pt, volumepts);
        }
    }
    return restvol;
}

computePressureRestLength (volume)

computePressureRestLength関数の計算モードをvolumeとしてオーバーロードしたもの

float
computePressureRestLength(const int geo; const int inpts[])
{
    // modeをvolumeとし体積計算
    return computePressureRestLength(geo, inpts, "volume");
}

createPressureConstraint

体積を維持するPressure(圧力)拘束を作成

// Create a pressure constraint over the provided group of points.
void
createPressureConstraint(const int geo; const string srcgrp;
                        const int outgeo; const string outgrp)
{
    // srcgrpに属するポイント番号の取得
    int pts[] = expandpointgroup(geo, srcgrp);

    // Don't create empty constraint.
    // srcgrpに属するポイントがない場合はreturn
    if (len(pts) == 0)
        return;
    
    // polylineの作成
    int prim = addprim(outgeo, "polyline", pts);
    
    // volume attributeを用いて体積計算
    float restvol = computePressureRestLength(geo, pts, "volume");
    
    // restlength (prim) <- restvole
    setprimattrib(outgeo, "restlength", prim, restvol);
    
    // outgrp <- 1
    setprimgroup(outgeo, outgrp, prim, 1);
}

createShapeMatchConstraint

オブジェクトの形状を維持するShape Match拘束を作成

// srcgrpに属するポイントからShapematch Constraintを作成
// point attrs: rest <- P
// grp: outgrp <- 1

// Create a shapematch constraint over the provided group of points.
void
createShapeMatchConstraint(const int geo; const string srcgrp;
                          const int outgeo; const string outgrp)
{
    // srcgrpに属するポイント番号の取得
    int pts[] = expandpointgroup(geo, srcgrp);

    // Don't create empty constraint.
    // srcgrpに属するポイントがない場合はreturn
    if (len(pts) == 0)
        return;

    // Set rest for each point from current geo position.
    // 現在のポイントポジションからrest attributeを設定
    foreach(int pt; pts)
    {
        vector pos = point(geo, "P", pt);
        
        // rest (pt) <- pos
        setpointattrib(outgeo, "rest", pt, pos);
    }

    // polylineの作成
    int prim = addprim(outgeo, "polyline", pts);

    // outgrp <- 1
    setprimgroup(outgeo, outgrp, prim, 1);
}

computeTetVolumeRestLength

四面体の体積を計算

float
computeTetVolumeRestLength(const int geo; const int pt0, pt1, pt2, pt3)
{
    vector p0 = point(geo, "P", pt0);
    vector p1 = point(geo, "P", pt1);
    vector p2 = point(geo, "P", pt2);
    vector p3 = point(geo, "P", pt3);

    return dot(cross(p0 - p1, p0 - p2), p0 - p3) / 6;
}

createTetVolumeConstraint

四面体の体積を維持するTet Volume拘束を作成

void
createTetVolumeConstraint(const int geo; const int ptnum; const string srcgrp;
                          const int outgeo; const string outgrp)
{
    // プリミティブの取得
    int prims[] = pointprims(geo, ptnum);

    foreach(int prim; prims)
    {
        // Process tets where this point is minimum point number,
        // and all points are in src group.
        // Typeid: 21 -> Tetrahedron のみ処理対象
        if(primintrinsic(geo, "typeid", prim) != 21)
            continue;

        // プリミティブからポイント取得
        int pts[] = primpoints(geo, prim);

        // 4つのポイントを持ち、
        // 指定されたポイント番号が最小であり、
        // 残りのポイントがすべてsrcgrpに属している場合のみ処理対象
        if (len(pts) != 4 ||
            ptnum != min(pts) ||
            !inpointgroup(geo, srcgrp, pts[1]) ||
            !inpointgroup(geo, srcgrp, pts[2]) ||
            !inpointgroup(geo, srcgrp, pts[3]))
            continue;

        // Tetrahedron Primitiveの作成
        int cprim = addprim(outgeo, "tet", pts);

        // 体積計算
        float volume = computeTetVolumeRestLength(geo, pts[0], pts[1], pts[2], pts[3]);
        
        // restlength (prim) <- volume
        setprimattrib(outgeo, "restlength", cprim, volume);

        // outgrp <- 1
        setprimgroup(outgeo, outgrp, cprim, 1);
    }
}

computeAngleRestLength

3つの連続するポイントを受け取り、その間の角度を計算

float
computeAngleRestLength(const int geo;
                      const int pt0, pt1, pt2)
{
    vector p0 = point(geo, "P", pt0);
    vector p1 = point(geo, "P", pt1);
    vector p2 = point(geo, "P", pt2);

    vector n1 = normalize(p1 - p0);
    vector n2 = normalize(p2 - p1);

    // n1とn2の角度を計算
    // theta = acos(dot)
    return degrees(acos(clamp(dot(n1, n2), -1, 1)));
}

computeVertexCentroidDistance

三角形の座標を受け取り、その座標から重心までの距離を計算

// Compute distance from vertex to centroid for triangle
float
computeVertexCentroidDistance(const vector p0, p1, p2)
{
    vector c = (p0 + p1 + p2) / 3;

    // Distance from vertex to centroid.
    return length(p1 - c);
}

computeTriangleBendRestLength

三角形のインデックスを受け取り、その座標から重心までの距離を計算