山东大学 计算机图形学实验 二维网格剖分 Catmull-Clark算法

原理可以参考这里

想自己写的记得注意边界条件：度数<=3的顶点不做位置调整，只和一个面相邻的边点直接取边的中点

#include <iostream>
#include <fstream>
#include <vector>
#include <sstream>
#include <algorithm>
#include <cmath>
#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

struct V{
    int v_id;
    double x, y;
};

struct E{
    int e_id, V1, V2;
    //vector<int> nf;
};

class F{
public:
    int f_id;
    vector<int> es; // 属于该面的所有边的索引
    vector<int> vs; // 属于该面的所有点的索引

    F(){}

    F(const F& f){
        f_id = f.f_id;
        es.clear();
        vs.clear();
        for (auto &e: f.es)
            es.push_back(e);
        for (auto &v: f.vs)
            vs.push_back(v);
    }

    void operator = (const F& f){
        f_id = f.f_id;
        es.clear();
        vs.clear();
        for (auto &e: f.es)
            es.push_back(e);
        for (auto &v: f.vs)
            vs.push_back(v);
    }
};
 
void readFile(const string &filename, vector<V> &vs, vector<E> &es, vector<F> &fs){
    ifstream file(filename);
    string line;
    int i = 0, j = 0, k = 0;
    while (getline(file, line)){
        istringstream iss(line);
        string type;
        iss >> type;
        if (type == "v"){ // 顶点信息
            V v;
            i++;
            v.v_id = i;
            iss >> v.x >> v.y;
            vs.push_back(v);
        }
        else if (type == "e"){ // 边信息
            E e;
            j++;
            e.e_id = j;
            iss >> e.V1 >> e.V2;
            es.push_back(e);
        }
        else if (type == "f"){ // 面信息
            F f;
            k++;
            f.f_id = k;
            string vertexIndex;
            while (iss >> vertexIndex){
                istringstream vertexIss(vertexIndex);
                int vIndex;
                vertexIss >> vIndex;
                f.vs.push_back(vIndex);
            }
            for (auto &e : es){
                int count = 0;
                for (auto v : f.vs){
                    if (e.V1 == v || e.V2 == v){
                        count++;
                        if (count == 2){
                            f.es.push_back(e.e_id);
                            //e.nf.push_back(k);
                            break;
                        }
                    }

                }
            }
            fs.push_back(f);
        }
    }
}
void writeFile(const string &filename, const vector<V> &vs, const vector<E> &es, const vector<F> &fs){
    ofstream file(filename);
 
    // 写入顶点信息
    for (const auto &v : vs){
        file << "v " << v.x << " " << v.y << endl;
    }
 
    // 写入边信息
    for (const auto &e : es){
        file << "e " << e.V1 << " " << e.V2 << endl;
    }
 
    // 写入面信息
    for (const auto &f : fs){
        file << "f";
        for (int v : f.vs){
            file << " " << v; // 从1开始的索引转换为从0开始
        }
        file << endl;
    }
}

void draw(const vector<V> &vs, const vector<E> &es){
    cv::Mat image(600, 800, CV_8UC3, cv::Scalar(255, 255, 255));
    for (const auto &v : vs){
        cv::Point p(v.x * 200, v.y * 200);
        circle(image, p, 3, Scalar(0, 0, 0), -1);
    }
    for (const auto &e : es){
        cv::Point p1(vs[e.V1 - 1].x * 200, vs[e.V1 - 1].y * 200);
        cv::Point p2(vs[e.V2 - 1].x * 200, vs[e.V2 - 1].y * 200);
        cv::line(image, p1, p2, cv::Scalar(0, 0, 255), 1);
    }
    cv::imshow("画板", image);
    cv::waitKey(0);
}

void catmullClark(vector<V> &vs, vector<E> &es, vector<F> &fs){
    vector<E> new_es;
    vector<F> new_fs;
    int v_cnt = vs.size();
    int n = v_cnt;
    int e_cnt = es.size();
    int f_cnt = fs.size();
    int m = 0;
    int k = 0;
    // 面点
    for (auto &f : fs){
        double x = 0, y = 0;
        for (auto& v_index : f.vs){
            x += vs[v_index - 1].x; 
            y += vs[v_index - 1].y;
        }
        x /= 1.0 * f.vs.size();
        y /= 1.0 * f.vs.size();

        //circle(image, Point(x * 200, y * 200), 6, Scalar(0, 0, 0), -1);

        vs.push_back(V{v_cnt + f.f_id, x, y});
    }
 
    // 边点
    for (auto &e : es){
        double i = 2;
        double x = vs[e.V1 - 1].x + vs[e.V2 - 1].x; 
        double y = vs[e.V1 - 1].y + vs[e.V2 - 1].y;
        double x0 = x, y0 = y;
        for (auto &f : fs){
            for (auto& e_index : f.es){//查找包含e的面
                if (e_index == e.e_id){ 
// v_cnt + f.f_id – 1为对应面点的下标
                    x += vs[v_cnt + f.f_id - 1].x; 
                    y += vs[v_cnt + f.f_id - 1].y;
                    ++m;
                    new_es.push_back(E{m, v_cnt + f_cnt + e.e_id, v_cnt + f.f_id});
                    i += 1.0;
                }
            }
        }
        if (i < 3.5) //面数<=1，取中点
            x = x0, y = y0, i = 2.0;
        x = x / i;
        y = y / i;
        V edgeCenter = V{v_cnt + f_cnt + e.e_id, x, y};
        //边点和两端点连边
        ++m;
        new_es.push_back(E{m, v_cnt + f_cnt + e.e_id, e.V1});
        ++m;
        new_es.push_back(E{m, v_cnt + f_cnt + e.e_id, e.V2});

        vs.push_back(edgeCenter);//边点加入点集
    }
    for (auto& f : fs){//连接剩余的边，得到新面
        for (auto& v_index : f.vs){
            k++;
            F new_f1;
            new_f1.f_id = k;
            //面点-边点-原顶点-边点构成一个新面
            //面和原顶点加入f的点集
            new_f1.vs.push_back(v_index);
            new_f1.vs.push_back(v_cnt + f.f_id);
            for (auto& e_index : f.es) {
                //查找得到对应的边点
                if (es[e_index - 1].V1 == v_index || es[e_index - 1].V2 == v_index){
                    new_f1.vs.push_back(v_cnt + f_cnt + e_index);
                }
            }
            for (auto &e:new_es){
                //在新边集中查找，得到f的边集
                int c = 0;
                for (auto &v:new_f1.vs){
                    if (e.V1 == v || e.V2 == v){
                        ++c;
                    }
                    if (c == 2){
                        new_f1.es.push_back(e.e_id);
                        break;
                    }
                }
            }
            // cout<<endl;
            new_fs.push_back(new_f1);
        }
    }
    vector<V> modiv;
    for (auto &v : vs){
        n--;
        if (n < 0)
            break;
        //更新原顶点
        double x = 0, y = 0;
        double count = 0, c2 = 0;
        for (auto& f : fs){
            for (auto& v_index : f.vs){//计算面点部分
                if (v_index == v.v_id){
                    x += vs[v_cnt + f.f_id - 1].x;
                    y += vs[v_cnt + f.f_id - 1].y;
                    count += 1.0;
                }
            }
        }
        for (auto& e : es) {//计算边点部分
            if (e.V1 == v.v_id || e.V2 == v.v_id){
                x += 2.0 * vs[v_cnt + f_cnt + e.e_id - 1].x; 
                y += 2.0 * vs[v_cnt + f_cnt + e.e_id - 1].y;
            }
        }
        if (count < 4) continue;//度数小于4，跳过
        x /= 4.0 * count;
        y /= 4.0 * count;
        x += v.x / 4.0, y += v.y / 4.0;

        modiv.push_back(V{v.v_id, x, y});
    }
    for (auto &v:modiv){
        vs[v.v_id - 1].x = v.x;
        vs[v.v_id - 1].y = v.y;
    }
    es.clear();
    for (auto& e:new_es)
        es.push_back(e);
    fs.clear();
    for (auto& f:new_fs)
        fs.push_back(f);
//cv::imshow("画板", image);
 //cv::waitKey(0);
}
 
int main(){
    vector<V> vs;
    vector<E> es;
    vector<F> fs;
    int op;
    cin >> op;
    readFile("dat", vs, es, fs);
    for (int i = 0; i < op; i++){
        catmullClark(vs, es, fs);
       
    }
    writeFile("output", vs, es, fs);
    draw(vs, es);
    vs.clear();
    fs.clear();
    es.clear();
    return 0;
}