【JavaScript】 Catmull-Clark subdivision surface の実装例

はじめに

Catmull-Clark subdivision surfaceのErkamanの実装を見つけたので、解説します。
該当ソースはindex.jsの_catmullClark関数です。1行ずつではありませんが、全て解説します。

0. 引数

positionsは点配列でArray<Array<number>>やArray<float32Array>です。
cellsは面配列で、Array<Array<number>>で面の頂点のインデックスの配列です。

index.js: 23～27行目

/*
Implement Catmull-Clark subvision, as it is described on Wikipedia
 */
function _catmullClark(positions, cells) {

1. 面、辺、点配列の構築
面、辺、点配列の構築し、互いに辿りやすいようにします。
originalPointsに点配列を、facesに面配列を、edgesに辺配列を構築していきます。
originalPointsを構築する際に、点の座標.point、点を含む面のインデックス配列.faces、点を含む辺のインデックス配列.edgesも構築します。
facesを構築する際に、面の中心点.facePoint、面の頂点.pointsも構築します。面の中心点は「面の新しい点」と呼ぶことにします。
edgesを構築する際に、辺の端点のインデックス配列.points、辺を含む面のインデックス配列.facesを構築します。

index.js: 28～166行目

    // original points, indexed by their indices.
    // For every point, we store adjacent faces and adjacent edges.
    originalPoints = [];

    // original faces, in their original order.
    // For every face, we store the edges, the points, and the face point.
    faces = [];

    // original edges. Indexed by the sorted indices of their vertices
    // So the edge whose edge vertices has index `6` and `2`, will be 
    // indexed by the array [2,6]
    edges = [];


    /*
    First we collect all the information that we need to run the algorithm.
    Each point must know its adjacent edges and faces.
    Each face must know its edges and points.
    Each edge must know its adjacent faces and points.

    We collect all this information in this loop.
     */
    for (var iCell = 0; iCell < cells.length; ++iCell) {

        var cellPositions = cells[iCell];
        var facePoints = [];

        // initialize:
        faces[iCell] = {};


        // go through all the points of the face.
        for (var j = 0; j < cellPositions.length; ++j) {

            var positionIndex = cellPositions[j];

            var pointObject;

            /*
           On the fly, Create an object for every point.
             */
            if (typeof originalPoints[positionIndex] === 'undefined') {
                // create the object on the fly.
                var v = positions[positionIndex];

                var vec = vec3.fromValues(v[0], v[1], v[2]);
                pointObject = {
                    point: vec,
                    faces: [],
                    edges: new Set(),

                };

                originalPoints[positionIndex] = pointObject;
            } else {
                pointObject = originalPoints[positionIndex];
            }

            // every point should have a reference to its face.
            pointObject.faces.push(faces[iCell]);

            facePoints.push(pointObject);
        }

        // every face should now its points.
        faces[iCell].points = facePoints;

        var avg = vec3.fromValues(0, 0, 0);

        // now compute the facepoint(see Wikipedia).
        for (var i = 0; i < faces[iCell].points.length; ++i) {
            var v = faces[iCell].points[i].point;
            vec3.add(avg, v, avg);
        }
        vec3.scale(avg, avg, 1.0 / faces[iCell].points.length);
        faces[iCell].facePoint = avg;

        var faceEdges = [];

        // go through all the edges of the face.
        for (var iEdge = 0; iEdge < cellPositions.length; ++iEdge) {

            var edge;

            if (cellPositions.length == 3) { // for triangles
                if (iEdge == 0) {
                    edge = [cellPositions[0], cellPositions[1]];
                } else if (iEdge == 1) {
                    edge = [cellPositions[1], cellPositions[2]];
                } else if (iEdge == 2) {
                    edge = [cellPositions[2], cellPositions[0]];
                }
            } else { // for quads.
                if (iEdge == 0) {
                    edge = [cellPositions[0], cellPositions[1]];
                } else if (iEdge == 1) {
                    edge = [cellPositions[1], cellPositions[2]];
                } else if (iEdge == 2) {
                    edge = [cellPositions[2], cellPositions[3]];
                } else if (iEdge == 3) {
                    edge = [cellPositions[3], cellPositions[0]];
                }
            }

            // every edge is represented by the sorted indices of its vertices.
            // (the sorting ensures that [1,2] and [2,1] are considered the same edge, which they are )
            edge = _sort(edge);

            var edgeObject;
            // on the fly, create an edge object.
            if (typeof edges[edge] === 'undefined') {

                edgeObject = {
                    points: [originalPoints[edge[0]], originalPoints[edge[1]]],
                    faces: []

                };

                edges[edge] = edgeObject;
            } else {
                edgeObject = edges[edge];
            }

            // every edge should know its adjacent faces.
            edgeObject.faces.push(faces[iCell]);

            // every point should know its adjacent edges.
            edgeObject.points[0].edges.add(edgeObject);
            edgeObject.points[1].edges.add(edgeObject);


            faceEdges.push(edgeObject);
        }

        // every face should know its edges.
        faces[iCell].edges = faceEdges;
    }

2. 辺の新しい点と辺の中点を求める

辺の含む面の新しい点がA,F、辺の端点をM,Eとすると、辺の新しい点は です。
あとで使用するので、ここで辺の中点を求めているようです。

index.js: 167～210行目

    // Compute the edge points and the midpoints of every edge.
    for (key in edges) {

        var edge = edges[key];

        var avg = vec3.fromValues(0, 0, 0);
        var count = 0;

        // add face points of edge.
        for (var i = 0; i < edge.faces.length; ++i) {
            var facePoint = edge.faces[i].facePoint;
            vec3.add(avg, facePoint, avg);
            ++count;
        }

        // sum together the two endpoints.
        for (var i = 0; i < edge.points.length; ++i) {
            var endPoint = edge.points[i].point;
            vec3.add(avg, endPoint, avg);
            ++count;
        }

        // finally, compute edge point.
        vec3.scale(avg, avg, 1.0 / count);
        edge.edgePoint = avg;

        /*
         Next we compute the midpoint.
         */
        count = 0;
        var avg2 = vec3.fromValues(0, 0, 0);

        for (var i = 0; i < edge.points.length; ++i) {
            var endPoint = edge.points[i].point;
            vec3.add(avg2, endPoint, avg2);
            ++count;
        }
        vec3.scale(avg2, avg2, 1.0 / count);


        edge.midPoint = avg2;
    }

3. 新しい点の座標を求めます。

元の点をP、Pを含む面の新しい点の平均をF、Pに接続する辺の中点の平均をRとすると、新しい点の座標は以下です。
nはPを含む面の数とします。(nはPに接続する辺の数でもあります。)
この時、新しい点はです。

プログラムでは上記のアルゴリズムとは、ほんの少し新しい点の求め方を変えていますが勿論同じです。

index.js: 211～238行目

/*
     Each original point is moved to the position
     (F + 2R + (n-3)P) / n. See the wikipedia article for more details.
     */
    for (var i = 0; i < positions.length; ++i) {

        var point = originalPoints[i];
        var n = point.faces.length;
        var newPoint = vec3.fromValues(0, 0, 0);

        for (var j = 0; j < point.faces.length; ++j) {
            var facePoint = point.faces[j].facePoint;
            vec3.add(newPoint, newPoint, facePoint);
        }

        for (var edge of point.edges) {
            _mad(newPoint, newPoint, edge.midPoint, 2);
        }
        vec3.scale(newPoint, newPoint, 1.0 / n);

        _mad(newPoint, newPoint, point.point, n - 3);

        vec3.scale(newPoint, newPoint, 1.0 / n);

        point.newPoint = newPoint

    }

4. 新しい面を作成する

戻り値である新しい頂点配列、新しい面配列はnewPositions、newCellsです。
面配列をたどって、その面の頂点分だけ新しい面配列の要素を作成します。
新しい面は、該当点の新しい点と該当点に接続する辺の新しい点2つと、その面の新しい点の4点からなる面です。
新しい面の点の並び順に関しては、ソースをご覧ください。

index.js: 239～288行目

    newPositions = [];
    newCells = [];

    var index = 0;

    /*
     We create indices on the fly by using this method.
     The index of every vertex is stored in the vertex, in a property named `index`.
     */
    function getIndex(p) {
        if (!("index" in p)) {
            p.index = index++;
            newPositions.push([p[0], p[1], p[2]]);

        }
        return p.index;

    }

    /*
     We go through all faces.
     Triangle face we subdivide into 3 new quads.
     Quad faces we subdivide into 4 new quads.

     */
    for (var iFace = 0; iFace < faces.length; ++iFace) {

        var face = faces[iFace];

        for (var iPoint = 0; iPoint < face.points.length; ++iPoint) {
            var point = face.points[iPoint];

            var a = point.newPoint;
            var b = face.edges[(iPoint + 0) % face.edges.length].edgePoint;
            var c = face.facePoint;
            var d = face.edges[(iPoint + face.edges.length - 1) % face.edges.length].edgePoint;

            var ia = getIndex(a);
            var ib = getIndex(b);
            var ic = getIndex(c);
            var id = getIndex(d);

            newCells.push([id, ia, ib, ic]);
        }
    }


    return {positions: newPositions, cells: newCells};

}

最後に

このプログラム書いている人上手ですね。勉強になります。

参考リンク

Catmull-Clark subdivision surface のアルゴリズム解説