cc 细分 网格细分

简介

pixel的大佬，就是cc细分论文的作者。

wiki的步骤是我见过的比较清晰的版本

Recursive evaluation
Catmull–Clark surfaces are defined recursively, using the following refinement scheme:[1]

Start with a mesh of an arbitrary polyhedron. All the vertices in this mesh shall be called original points.

For each face, add a face point
Set each face point to be the average of all original points for the respective face.
For each edge, add an edge point.
Set each edge point to be the average of the two neighbouring face points and its two original endpoints.
For each face point, add an edge for every edge of the face, connecting the face point to each edge point for the face.
For each original point P, take the average F of all n (recently created) face points for faces touching P, and take the average R of all n edge midpoints for (original) edges touching P, where each edge midpoint is the average of its two endpoint vertices (not to be confused with new "edge points" above). (Note that from the perspective of a vertex P, the number of edges neighboring P is also the number of adjacent faces, hence n). Move each original point to the point
({displaystyle {frac {F+2R+(n-3)P}{n}}}{displaystyle {frac {F+2R+(n-3)P}{n}}})
This is the barycenter of P, R and F with respective weights (n − 3), 2 and 1.
Connect each new face point to the new edge points of all original edges defining the original face.
Connect each new vertex point to the new edge points of all original edges incident on the original vertex.
Define new faces as enclosed by edges.
The new mesh will consist only of quadrilaterals, which in general will not be planar. The new mesh will generally look smoother than the old mesh.

Repeated subdivision results in smoother meshes. It can be shown that the limit surface obtained by this refinement process is at least ({displaystyle {mathcal {C}}^{1}}mathcal{C}^1) at extraordinary vertices and ({displaystyle {mathcal {C}}^{2}}{mathcal {C}}^{2}) everywhere else (when n indicates how many derivatives are continuous, we speak of ({displaystyle {mathcal {C}}^{n}}mathcal{C}^n) continuity). After one iteration, the number of extraordinary points on the surface remains constant.

The arbitrary-looking barycenter formula was chosen by Catmull and Clark based on the aesthetic appearance of the resulting surfaces rather than on a mathematical derivation, although Catmull and Clark do go to great lengths to rigorously show that the method converges to bicubic B-spline surfaces.[1]

参考

链接代码
http://rosettacode.org/wiki/Catmull–Clark_subdivision_surface#Python
https://blog.csdn.net/McQueen_LT/article/details/106102609

code

虽然 rosettacode 提供了 各种类型的细分代码，但是 个人觉得 都不是特别让人容易懂，使用openmesh重写了部分代码，结果很舒适。暂时有点缺憾，不能处理不封闭的网格。

#ifndef __CC_HPP__
#define __CC_HPP__
#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <map>
typedef OpenMesh::PolyMesh_ArrayKernelT<> MyMesh;

typedef struct FFACE{
    int faceId;
    int edgeId1;
    int vertexId;
    int edgeId2;
}FFACE;

class CC{
private:
    MyMesh mesh;
    std::map<int, OpenMesh::Vec3d> facePoints;
    std::map<int, OpenMesh::Vec3d> edgePoints;
    std::map<int, OpenMesh::Vec3d> vertexPoints;
    int times;
public:
    CC(int _times);
    ~CC();
    void genCube();
    void genFacePoint();
    void genEdgePoint();
    void genVertexPoint();
    void connectPoints();
    void genMesh(std::string name);
    int pipeLine();
};
#endif

/*
 * @Author: your name
 * @Date: 2020-09-02 13:49:22
 * @LastEditTime: 2020-09-03 14:35:38
 * @LastEditors: Please set LastEditors
 * @Description: In User Settings Edit
 * @FilePath: /cc/cc.cc
 */

#include "cc.hpp"
#include <iostream>
#include <unordered_map>

/**
 * @description: 构造函数 
 * @param {type} 
 * @return {type} 
 */
CC::CC(int _times):times(_times){
    genCube();
}

/**
 * @description: 析构函数 
 * @param {type} 
 * @return {type} 
 */
CC::~CC() {}

/**
 * @description: 一整套流程
 * @param {type} 
 * @return {int} 管线运行是否成功 
 */
int CC::pipeLine(){
    for(int i=0; i<times; i++){
        genFacePoint();
        genEdgePoint();
        genVertexPoint();
        connectPoints();
        char a[20];
        sprintf(a, "output%d.off", i);
        genMesh(a);
        // 清空变量
        std::cout << "[DEBUG] 迭代了第 " << i << " 次。" << std::endl;
    }
    
    return 0;
}


/**
 * @description: 生成一个立方体（四边形网格） 
 * @param {type} 
 * @return {type} 
 */
void CC::genCube()
{
    MyMesh::VertexHandle vhandle[9];
    vhandle[0] = mesh.add_vertex(MyMesh::Point(-1, -1, 1));
    vhandle[1] = mesh.add_vertex(MyMesh::Point(1, -1, 1));
    vhandle[2] = mesh.add_vertex(MyMesh::Point(1, 1, 1));
    vhandle[3] = mesh.add_vertex(MyMesh::Point(-1, 1, 1));
    vhandle[4] = mesh.add_vertex(MyMesh::Point(-1, -1, -1));
    vhandle[5] = mesh.add_vertex(MyMesh::Point(1, -1, -1));
    vhandle[6] = mesh.add_vertex(MyMesh::Point(1, 1, -1));
    vhandle[7] = mesh.add_vertex(MyMesh::Point(-1, 1, -1));

    std::vector<MyMesh::VertexHandle> face_vhandles;
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[3]);
    mesh.add_face(face_vhandles);

    face_vhandles.clear();
    face_vhandles.push_back(vhandle[7]);
    face_vhandles.push_back(vhandle[6]);
    face_vhandles.push_back(vhandle[5]);
    face_vhandles.push_back(vhandle[4]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[4]);
    face_vhandles.push_back(vhandle[5]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[1]);
    face_vhandles.push_back(vhandle[5]);
    face_vhandles.push_back(vhandle[6]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[3]);
    face_vhandles.push_back(vhandle[2]);
    face_vhandles.push_back(vhandle[6]);
    face_vhandles.push_back(vhandle[7]);
    mesh.add_face(face_vhandles);
    face_vhandles.clear();
    face_vhandles.push_back(vhandle[0]);
    face_vhandles.push_back(vhandle[3]);
    face_vhandles.push_back(vhandle[7]);
    face_vhandles.push_back(vhandle[4]);
    mesh.add_face(face_vhandles);
}

/**
 * @description: 生成所有面点
 * @param {type} 
 * @return {type} 
 */
void CC::genFacePoint()
{
    facePoints.clear();
    for (const auto &fh : mesh.faces())
    {
        OpenMesh::Vec3d facePoint(0, 0, 0);
        int facePointsNumber = 0;
        for (const auto &fvh : mesh.fv_range(fh))
        {
            OpenMesh::DefaultTraits::Point point = mesh.point(fvh);
            facePoint += point;
            facePointsNumber++;
        }
        facePoint /= facePointsNumber;
        facePoints[fh.idx()] = facePoint;
    }
}

/**
 * @description: 生成所有的边点暂时不考虑hole 就是全是日子结构的
 * @param {type} 
 * @return {type} 
 */
void CC::genEdgePoint()
{
    edgePoints.clear();
    for (auto e_it = mesh.edges_begin(); e_it != mesh.edges_end(); ++e_it)
    {
        // 得到边所代表的半边
        OpenMesh::HalfedgeHandle heh1 = mesh.halfedge_handle(*e_it, 0); // 默认一个方向的半边
        OpenMesh::Vec3d edgePoint(0, 0, 0);
        int edgePointsNumber = 0;
        OpenMesh::DefaultTraits::Point pointV = mesh.point(mesh.from_vertex_handle(heh1)); // 这条（半）边的起点
        OpenMesh::DefaultTraits::Point pointW = mesh.point(mesh.to_vertex_handle(heh1));   // 这条（半）边的终点
        OpenMesh::FaceHandle fh1 = mesh.face_handle(heh1);
        OpenMesh::Vec3d fh1p = facePoints[fh1.idx()];                 // 得到面点
        OpenMesh::HalfedgeHandle heh = mesh.opposite_halfedge_handle(heh1); // 逆半边
        OpenMesh::FaceHandle fh2 = mesh.face_handle(heh);
        OpenMesh::Vec3d fh2p = facePoints[fh2.idx()]; // 得到面点
        edgePoints[heh1.idx()] = (fh1p + fh2p + pointV + pointW) / 4.0;
    }
}

/**
 * @description: 生成新的顶点
 * For each original point P, take the average F of all n (recently created) face points for faces touching P,
 * and take the average R of all n edge midpoints for (original) edges touching P, 
 * where each edge midpoint is the average of its two endpoint vertices (not to be confused with new "edge points" above). 
 * (Note that from the perspective of a vertex P, the number of edges neighboring P is also the number of adjacent faces, hence n). 
 * Move each original point to the point
 * @param {type} 
 * @return {type} 
 */
void CC::genVertexPoint()
{
    vertexPoints.clear();
    // 原始点接触的面的所有的面点的均值
    for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); v_it++)
    {
        OpenMesh::Vec3d originPoint = (OpenMesh::Vec3d)mesh.point(*v_it);
        OpenMesh::Vec3d facePoint(0, 0, 0);
        int faceNumber = 0;
        for (auto vf_it = mesh.vf_iter(*v_it); vf_it.is_valid(); ++vf_it)
        { //这个顶点所带有的面迭代器
            facePoint += facePoints[(*vf_it).idx()];
            faceNumber++;
        }
        facePoint /= faceNumber;

        // 原始点接触的边的中间点的值的均值 * 2
        OpenMesh::Vec3d edgePoint(0, 0, 0);
        int edgeNumber = 0;
        for(auto vv_it = mesh.vv_begin(*v_it); vv_it != mesh.vv_end(*v_it); vv_it++){ // 为啥还有孤立点？
            OpenMesh::Vec3d point =(OpenMesh::Vec3d)mesh.point(*vv_it);
            edgePoint += (point + originPoint) / 2;
            edgeNumber++;
        }
        if(faceNumber != edgeNumber){
            std::cout << "[debug] is not equal f " << faceNumber << " e " << edgeNumber << std::endl; 
        }
        std::cout << "[DEBUG] n " << edgeNumber << std::endl; 
        edgePoint /= edgeNumber;
        
        OpenMesh::Vec3d newPoint = (edgeNumber - 3.0) / (edgeNumber) * originPoint + facePoint * 1.0 / edgeNumber + 2.0 * edgePoint / edgeNumber;
        std::cout << "[DEBUG] oldPoint " << originPoint << " newPoint " << newPoint << std::endl;
        vertexPoints[(*v_it).idx()] = newPoint;
    }
}

/**
 * @description: 连接面点和边点   （要严格的封闭的四边形）
 * @param {type} 
 * @return {type} 
 */
void CC::connectPoints(){
    if (!mesh.has_vertex_status())  mesh.request_vertex_status();
	if (!mesh.has_face_status())    mesh.request_face_status();
	if (!mesh.has_edge_status())    mesh.request_edge_status();
    // 先将旧顶点移动到新顶点
    std::unordered_map<int, OpenMesh::VertexHandle> vertexHandle;
    for (auto v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); v_it++){
        mesh.set_point(*v_it, (OpenMesh::DefaultTraits::Point)vertexPoints[(*v_it).idx()]);
        vertexHandle[(*v_it).idx()] = *v_it;
    }
    std::vector<MyMesh::VertexHandle> facePointsHandle;
    std::vector<MyMesh::FaceHandle> faceHandle;
    // 再加入所有的要加入的面点
    for (const auto &fh : mesh.faces()) {
        facePointsHandle.push_back(mesh.add_vertex(MyMesh::Point(facePoints[fh.idx()])));
        faceHandle.push_back(fh);
    }
    std::unordered_map<int, MyMesh::VertexHandle> edgesHandle;
    // 再加入所有的要加入的边点
    for (const auto &eh : mesh.edges()){
        edgesHandle[mesh.halfedge_handle(eh, 0).idx()] = mesh.add_vertex(MyMesh::Point(edgePoints[mesh.halfedge_handle(eh, 0).idx()]));
    }

    std::vector<std::vector<MyMesh::VertexHandle>> v;
    int i = 0;
    for (const auto &fh : mesh.faces()){
        std::vector<MyMesh::VertexHandle> handlerVector(9);
        handlerVector[0] = facePointsHandle[i];
        i++;
        int pointsNumber = 1;
        for (const auto &feh : mesh.fe_range(fh)){
            auto fehh = mesh.halfedge_handle(feh, 0);
            handlerVector[pointsNumber] = edgesHandle[fehh.idx()];
            pointsNumber++;
        }
        for (const auto &fvh : mesh.fv_range(fh)){
            handlerVector[pointsNumber] = vertexHandle[fvh.idx()];
            pointsNumber++;
        }
        v.push_back(handlerVector);
    }
    // 开始删除面
    // for(int i=0; i<faceHandle.size(); i++){
    //     mesh.delete_face(faceHandle[i], false);
    // }
    for(int i=0; i<v.size(); i++){
        mesh.delete_face(faceHandle[i], true);
        std::vector<MyMesh::VertexHandle> handlerVector;
        for(int j=0; j<v[i].size(); j++){
            handlerVector.push_back(v[i][j]);
        }
        std::vector<MyMesh::VertexHandle> face_vhandles;
        face_vhandles.push_back(handlerVector[0]);
        face_vhandles.push_back(handlerVector[1]);
        face_vhandles.push_back(handlerVector[5]);
        face_vhandles.push_back(handlerVector[2]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[0]);
        face_vhandles.push_back(handlerVector[2]);
        face_vhandles.push_back(handlerVector[6]);
        face_vhandles.push_back(handlerVector[3]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[0]);
        face_vhandles.push_back(handlerVector[3]);
        face_vhandles.push_back(handlerVector[7]);
        face_vhandles.push_back(handlerVector[4]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();

        face_vhandles.push_back(handlerVector[0]);
        face_vhandles.push_back(handlerVector[4]);
        face_vhandles.push_back(handlerVector[8]);
        face_vhandles.push_back(handlerVector[1]);
        mesh.add_face(face_vhandles);
        face_vhandles.clear();
    }

	mesh.garbage_collection();
	if (mesh.has_vertex_status())  mesh.release_vertex_status();
	if (mesh.has_face_status())    mesh.release_face_status();
	if (mesh.has_edge_status())    mesh.release_edge_status();
}



/**
 * @description: 输出新网格 
 * @param {string} name 文件名称 默认 output.obj 
 * @return {type} 
 */
void CC::genMesh(std::string name){
    if(name == ""){
        name = "output.off";
    }
    try{
        if( !OpenMesh::IO::write_mesh(mesh, name) ){
            std::cerr << "Cannot write mesh to file 'output.off'" << std::endl;
            return;
        }
    }
    catch(std::exception &e){
        std::cerr << e.what() << std::endl;
        return;
    }
}

image

附赠Makefile

LOCAL_LIBRARY +=  -L/home/ling/lee/lib/openmesh/static
LOCAL_LDFLAGS += -lm -lpthread -ldl -lOpenMeshCore -lOpenMeshTools
LOCAL_CFLAGS += -I/home/ling/lee/include -I./include -D OM_STATIC_BUILD

CC := g++ -g
TARGETS1 = genCube
SRCS1 = mian1.cc cc.cc
OBJS1 = $(patsubst %.cc, %.o, $(SRCS1))

CFLAGS += $(LOCAL_CFLAGS)
LDFLAGS += $(LOCAL_LIBRARY) $(LOCAL_LDFLAGS)

$(info $(OBJS))
$(info $(TARGETS))

all: $(TARGETS1)


$(TARGETS1):$(OBJS1)
	$(CC)  -o $@ $^  $(LDFLAGS) $(CFLAGS)

$(OBJS1): %.o:%.cc
	$(CC) -c $< -o $@ $(CFLAGS)

clean :
	@rm -rf $(TARGETS1) $(OBJS1)

#.SUFFIXES:
.PHONY : all clean

附赠 main.cc

/*
 * @Author: your name
 * @Date: 2020-09-02 09:25:40
 * @LastEditTime: 2020-09-03 11:09:19
 * @LastEditors: Please set LastEditors
 * @Description: In User Settings Edit
 * @FilePath: /cc/mian1.cc
 */
#include <iostream>
#include <OpenMesh/Core/IO/MeshIO.hh>
#include <OpenMesh/Core/Mesh/PolyMesh_ArrayKernelT.hh>
#include "cc.hpp"


int main(){
    CC c(5);
    c.pipeLine();
    return 0;
}