Discrete Differential Geometry Lab-0: Combinatorial Surfaces

Posted on Edited on

为什么要开这个坑

学习Graphics\text{Graphics}，最少最少要学习33门课，渲染，模拟和几何各一门。几何方面的学习必不可少。此外，离散微分几何在物理模拟中也有重要应用（渲染我还不知道有没有），学这么一门课自然好处大大滴有。

这个系列用来记录课程的编程作业。

本次作业的要求是实现组合曲面中对于单形和复形的若干运算和操作，并且熟悉作业框架的各种接口。

1. buildVertexEdgeAdjacencyMatrix\text{1. buildVertexEdgeAdjacencyMatrix}

记住一点，在访问点/边/面的索引之前，都需要调用对应的require[Vertex|Edge|Face]Indices函数，否则会段错误。不知道这是什么奇葩设计（

框架使用Eigen\text{Eigen}数值线性代数库，使用稀疏矩阵的setFromTriplets方法可以很容易地从三元组列表构造稀疏矩阵。构建点边邻接矩阵自然是容易的，只需要枚举所有的半边，然后以对应的边的索引作为行号，对应的起点的索引作为列号即可。亲测，反过来是错的。

那么实现buildFaceEdgeAdjacencyMatrix函数也是同样的道理。

1
2
3
4
5
6
7
8
9
10
11
12
13
SparseMatrix<size_t> SimplicialComplexOperators::buildVertexEdgeAdjacencyMatrix() const {
    geometry->requireVertexIndices();
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    // Note: You can build an Eigen sparse matrix from triplets, then return it as a Geometry Central SparseMatrix.
    // See <https://eigen.tuxfamily.org/dox/group__TutorialSparse.html> for documentation.
    std::vector<Trip>tList;
    for(Halfedge he : mesh->halfedges())
        tList.emplace_back(geometry->edgeIndices[he.edge()], geometry->vertexIndices[he.tipVertex()], 1);
    SparseMatrix<size_t> SpM(mesh->nEdges(), mesh->nVertices());
    SpM.setFromTriplets(tList.begin(), tList.end());
    return SpM;
}

2. buildVertexVector\text{2. buildVertexVector}

这个就更简单了。。。但是注意一点，要构造的是顶点索引的onehot\text{onehot}向量（也就是每个分量非00即11，代表该顶点是否在选取的子集顶点中）。

1
2
3
4
5
6
7
8
Vector<size_t> SimplicialComplexOperators::buildVertexVector(const MeshSubset& subset) const {
    geometry->requireVertexIndices();
    Vector<size_t> Vec(mesh->nVertices());
    int cnt = 0;
    for(Vertex v : mesh->vertices())
        Vec[cnt++] = subset.vertices.count(geometry->vertexIndices[v]) > 0 ? 1 : 0;
    return Vec;
}

3. star\text{3. star}

实现St\text{St}算子，St(S)\text{St}(S)表示的曲面是：将所有与SS中顶点(0−0-单形)相邻的边(1−1-单形)和面(2−2-单形)都加入SS中，将所有与SS中的边相邻的面都加入SS中得到的曲面。

知道定义和接口以后就非常简单。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
MeshSubset SimplicialComplexOperators::star(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    MeshSubset st = subset;
    for(size_t idx : subset.edges)
    {
        Edge e = mesh->edge(idx);
        auto adjF = e.adjacentFaces();
        for(Face f : adjF)
        {
            size_t idxF = geometry->faceIndices[f];
            if(!st.faces.count(idxF)) st.faces.insert(idxF);
        }
    }
    for(size_t idx : st.vertices)
    {
        Vertex v = mesh->vertex(idx);
        auto adjE = v.adjacentEdges();
        for(Edge e : adjE)
        {
            size_t idxE = geometry->edgeIndices[e];
            if(!st.edges.count(idxE)) st.edges.insert(idxE);
        }
        auto adjF = v.adjacentFaces();
        for(Face f : adjF)
        {
            size_t idxF = geometry->faceIndices[f];
            if(!st.faces.count(idxF)) st.faces.insert(idxF);
        }
    }
    return st;
}

4. isComplex\text{4. isComplex}

判断给定的曲面是否是单纯复形。回忆单纯复形的定义，单纯复形本质上就是对子集操作封闭的集族，或者说子集系统。从几何角度来说，就是如果一个面在曲面中，那么该面的所有边以及顶点都在曲面中，并且如果一条边在曲面中，那么它的两个顶点都在曲面中。

这样我们就可以很容易地写出来。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
bool SimplicialComplexOperators::isComplex(const MeshSubset& subset) const {
    geometry->requireVertexIndices();
    geometry->requireEdgeIndices();
    for(size_t idx : subset.edges)
    {
        Edge e = mesh->edge(idx);
        if(!subset.vertices.count(geometry->vertexIndices[e.firstVertex()]) ||
            !subset.vertices.count(geometry->vertexIndices[e.secondVertex()]))
            return false;
    }
    for(size_t idx : subset.faces)
    {
        Face f = mesh->face(idx);
        auto Es = f.adjacentEdges();
        for(Edge e : Es)
            if (!subset.edges.count(geometry->edgeIndices[e]))
                return false;
    }
    return true;
}

5. isPureComplex\text{5. isPureComplex}

判断给定的曲面是否是纯(齐次)单纯复形，若是，还需判断纯单纯复形的维数。回忆概念的话，纯kk维单纯复形就是极大元集合的大小均相等的子集系统。

首先判断是不是单纯复形再进行后续检查。如果组合曲面中有面，那么就要求所有的边和顶点都是这些面的边和顶点，此时纯单纯复形的维数为22。如果组合曲面中没有面但是有边，就要求所有的顶点都是这些边的顶点，此时纯单纯复形的维数为11。如果既没有面也没有边，只有一群点，那纯单纯复形的维数就是00。

实际操作时，可以先用std::set保存所有边，然后遍历所有面并且删除边，检查最后得到的std::set是否为空。如果要检查顶点是否满足要求，有一个简单的办法是看顶点的度数。因为已经判断过是单纯复形，所以如果有不满足条件的顶点，一定是游离的顶点，度数为00。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
int SimplicialComplexOperators::isPureComplex(const MeshSubset& subset) const {
    if(!isComplex(subset)) return -1;
    if(!subset.faces.empty())
    {
        std::set<Edge> Es;
        for(size_t idx : subset.edges)
            Es.insert(mesh->edge(idx));
        for(size_t idx : subset.faces)
        {
            Face f = mesh->face(idx);
            for(Edge e : f.adjacentEdges())
                if(Es.count(e))Es.erase(e);
        }
        if(!Es.empty()) return -1;
        for(size_t idx : subset.vertices)
        {
            Vertex v = mesh->vertex(idx);
            if(!v.degree()) return -1;
        }
        return 2;
    }
    if(!subset.edges.empty())
    {
        for(size_t idx : subset.vertices)
        {
            Vertex v = mesh->vertex(idx);
            if(!v.degree()) return -1;
        }
        return 1;
    }
    return 0;
}

6. closure\text{6. closure}

实现闭包算子Cl\text{Cl}，Cl(S)\text{Cl}(S)表示的曲面是包括SS的最小单纯复形。

对于所有面，把面的所有边都加入曲面。然后检查新曲面的所有边，将边的所有顶点加入新曲面即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
MeshSubset SimplicialComplexOperators::closure(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    geometry->requireVertexIndices();
    MeshSubset st = subset;
    for(size_t idx : st.faces)
    {
        Face f = mesh->face(idx);
        auto adjE = f.adjacentEdges();
        for(Edge e : adjE)
        {
            size_t idxE = geometry->edgeIndices[e];
            if(!st.edges.count(idxE)) st.edges.insert(idxE);
        }
    }
    for(size_t idx : st.edges)
    {
        Edge e = mesh->edge(idx);
        if(!st.vertices.count(geometry->vertexIndices[e.firstVertex()]))
            st.vertices.insert(geometry->vertexIndices[e.firstVertex()]);
        if(!st.vertices.count(geometry->vertexIndices[e.secondVertex()]))
            st.vertices.insert(geometry->vertexIndices[e.secondVertex()]);
    }
    return st;
}

然后可以根据定义很快地完成link\text{link}。Lk(S)=Cl(St(S))\St(Cl(S))\text{Lk}(S)=\text{Cl}(\text{St}(S)) \backslash \text{St}(\text{Cl}(S))

1
2
3
4
5
6
7
8
9
10
11
MeshSubset SimplicialComplexOperators::link(const MeshSubset& subset) const {
    MeshSubset st1 = closure(star(subset));
    MeshSubset st2 = star(closure(subset));
    for(size_t idx : st2.vertices)
        st1.vertices.erase(idx);
    for(size_t idx : st2.faces)
        st1.faces.erase(idx);
    for(size_t idx : st2.edges)
        st1.edges.erase(idx);
    return st1;
}

7. boundary\text{7. boundary}

首先取闭包（闭包不改变边界，而且性质更好），然后取出所有度数小于等于11的点和度数小于等于11的边，再取闭包即可。

这么做是合理的，因为边界中一定不会出现面，所以只需要考虑哪些边和顶点应该在边界上。根据定义，度数小于等于11的顶点（这些是端点或者游离的顶点）和边都是满足要求的，除此以外，边界上度数的顶点难以根据原曲面直接判断，但是边界上的边很好判断，因此先把所有边都取出来，然后取闭包就知道边界上有哪些点了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
MeshSubset SimplicialComplexOperators::boundary(const MeshSubset& subset) const {
    // TODO
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    MeshSubset cl = closure(subset);
    MeshSubset st;
    for(size_t idx : cl.vertices)
    {
        Vertex v = mesh->vertex(idx);
        size_t deg = 0;
        for(Edge e : v.adjacentEdges())
        {
            if(cl.edges.count(geometry->edgeIndices[e]))
            {
                deg++;
                if(deg > 1) break;
            }
        }
        if(deg == 1) st.addVertex(idx);
    }
    for(size_t idx : cl.edges)
    {
        Edge e = mesh->edge(idx);
        size_t deg = 0;
        for(Face f : e.adjacentFaces())
        {
            if(cl.faces.count(geometry->faceIndices[f]))
            {
                deg++;
                if(deg > 1) break;
            }
        }
        if(deg == 1) st.addEdge(idx);
    }
    return closure(st);
}

8. 全部代码\text{8. 全部代码}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
using Trip = Eigen::Triplet<size_t>;
/*
 * Construct the unsigned vertex-edge adjacency matrix A0.
 *
 * Input:
 * Returns: The sparse vertex-edge adjacency matrix which gets stored in the global variable A0.
 */
SparseMatrix<size_t> SimplicialComplexOperators::buildVertexEdgeAdjacencyMatrix() const {
    geometry->requireVertexIndices();
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    // Note: You can build an Eigen sparse matrix from triplets, then return it as a Geometry Central SparseMatrix.
    // See <https://eigen.tuxfamily.org/dox/group__TutorialSparse.html> for documentation.
    std::vector<Trip>tList;
    for(Halfedge he : mesh->halfedges())
        tList.emplace_back(geometry->edgeIndices[he.edge()], geometry->vertexIndices[he.tipVertex()], 1);
    SparseMatrix<size_t> SpM(mesh->nEdges(), mesh->nVertices());
    SpM.setFromTriplets(tList.begin(), tList.end());
    return SpM;
}

/*
 * Construct the unsigned face-edge adjacency matrix A1.
 *
 * Input:
 * Returns: The sparse face-edge adjacency matrix which gets stored in the global variable A1.
 */
SparseMatrix<size_t> SimplicialComplexOperators::buildFaceEdgeAdjacencyMatrix() const {
    geometry->requireVertexIndices();
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    std::vector<Trip>tList;
    for(Halfedge he : mesh->halfedges())
        tList.emplace_back(geometry->faceIndices[he.face()], geometry->edgeIndices[he.edge()], 1);
    SparseMatrix<size_t> SpM(mesh->nFaces(), mesh->nEdges());
    SpM.setFromTriplets(tList.begin(), tList.end());
    return SpM;
}

/*
 * Construct a vector encoding the vertices in the selected subset of simplices.
 *
 * Input: Selected subset of simplices.
 * Returns: Vector of length |V|, where |V| = # of vertices in the mesh.
 */
Vector<size_t> SimplicialComplexOperators::buildVertexVector(const MeshSubset& subset) const {
    geometry->requireVertexIndices();
    Vector<size_t> Vec(mesh->nVertices());
    int cnt = 0;
    for(Vertex v : mesh->vertices())
        Vec[cnt++] = subset.vertices.count(geometry->vertexIndices[v]) > 0 ? 1 : 0;
    return Vec;
}

/*
 * Construct a vector encoding the edges in the selected subset of simplices.
 *
 * Input: Selected subset of simplices.
 * Returns: Vector of length |E|, where |E| = # of edges in mesh.
 */
Vector<size_t> SimplicialComplexOperators::buildEdgeVector(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    Vector<size_t> Vec(mesh->nEdges());
    int cnt = 0;
    for(Edge e : mesh->edges())
        Vec[cnt++] = subset.edges.count(geometry->edgeIndices[e]) > 0 ? 1 : 0;
    return Vec;
}

/*
 * Construct a vector encoding the faces in the selected subset of simplices.
 *
 * Input: Selected subset of simplices.
 * Returns: Vector of length |F|, where |F| = # of faces in mesh.
 */
Vector<size_t> SimplicialComplexOperators::buildFaceVector(const MeshSubset& subset) const {
    geometry->requireFaceIndices();
    Vector<size_t> Vec(mesh->nFaces());
    int cnt = 0;
    for(Face f : mesh->faces())
        Vec[cnt++] = subset.faces.count(geometry->faceIndices[f]) > 0 ? 1 : 0;
    return Vec;
}

/*
 * Compute the simplicial star St(S) of the selected subset of simplices.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: The star of the given subset.
 */
MeshSubset SimplicialComplexOperators::star(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    MeshSubset st = subset;
    for(size_t idx : subset.edges)
    {
        Edge e = mesh->edge(idx);
        auto adjF = e.adjacentFaces();
        for(Face f : adjF)
        {
            size_t idxF = geometry->faceIndices[f];
            if(!st.faces.count(idxF)) st.faces.insert(idxF);
        }
    }
    for(size_t idx : st.vertices)
    {
        Vertex v = mesh->vertex(idx);
        auto adjE = v.adjacentEdges();
        for(Edge e : adjE)
        {
            size_t idxE = geometry->edgeIndices[e];
            if(!st.edges.count(idxE)) st.edges.insert(idxE);
        }
        auto adjF = v.adjacentFaces();
        for(Face f : adjF)
        {
            size_t idxF = geometry->faceIndices[f];
            if(!st.faces.count(idxF)) st.faces.insert(idxF);
        }
    }
    return st;
}


/*
 * Compute the closure Cl(S) of the selected subset of simplices.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: The closure of the given subset.
 */
MeshSubset SimplicialComplexOperators::closure(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    geometry->requireVertexIndices();
    MeshSubset st = subset;
    for(size_t idx : st.faces)
    {
        Face f = mesh->face(idx);
        auto adjE = f.adjacentEdges();
        for(Edge e : adjE)
        {
            size_t idxE = geometry->edgeIndices[e];
            if(!st.edges.count(idxE)) st.edges.insert(idxE);
        }
    }
    for(size_t idx : st.edges)
    {
        Edge e = mesh->edge(idx);
        if(!st.vertices.count(geometry->vertexIndices[e.firstVertex()]))
            st.vertices.insert(geometry->vertexIndices[e.firstVertex()]);
        if(!st.vertices.count(geometry->vertexIndices[e.secondVertex()]))
            st.vertices.insert(geometry->vertexIndices[e.secondVertex()]);
    }
    return st;
}

/*
 * Compute the link Lk(S) of the selected subset of simplices.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: The link of the given subset.
 */
MeshSubset SimplicialComplexOperators::link(const MeshSubset& subset) const {
    MeshSubset st1 = closure(star(subset));
    MeshSubset st2 = star(closure(subset));
    for(size_t idx : st2.vertices)
        st1.vertices.erase(idx);
    for(size_t idx : st2.faces)
        st1.faces.erase(idx);
    for(size_t idx : st2.edges)
        st1.edges.erase(idx);
    return st1;
}

/*
 * Return true if the selected subset is a simplicial complex, false otherwise.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: True if given subset is a simplicial complex, false otherwise.
 */
bool SimplicialComplexOperators::isComplex(const MeshSubset& subset) const {
    geometry->requireVertexIndices();
    geometry->requireEdgeIndices();
    for(size_t idx : subset.edges)
    {
        Edge e = mesh->edge(idx);
        if(!subset.vertices.count(geometry->vertexIndices[e.firstVertex()]) ||
            !subset.vertices.count(geometry->vertexIndices[e.secondVertex()]))
            return false;
    }
    for(size_t idx : subset.faces)
    {
        Face f = mesh->face(idx);
        auto Es = f.adjacentEdges();
        for(Edge e : Es)
            if (!subset.edges.count(geometry->edgeIndices[e]))
                return false;
    }
    return true;
}

/*
 * Check if the given subset S is a pure simplicial complex. If so, return the degree of the complex. Otherwise, return
 * -1.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: int representing the degree of the given complex (-1 if not pure)
 */
int SimplicialComplexOperators::isPureComplex(const MeshSubset& subset) const {
    if(!isComplex(subset)) return -1; 
    if(!subset.faces.empty())
    {
        std::set<Edge> Es;
        for(size_t idx : subset.edges)
            Es.insert(mesh->edge(idx));
        for(size_t idx : subset.faces)
        {
            Face f = mesh->face(idx);
            for(Edge e : f.adjacentEdges())
                if(Es.count(e))Es.erase(e);
        }
        if(!Es.empty()) return -1;
        for(size_t idx : subset.vertices)
        {
            Vertex v = mesh->vertex(idx);
            if(!v.degree()) return -1;
        }
        return 2;
    }
    if(!subset.edges.empty())
    {
        for(size_t idx : subset.vertices)
        {
            Vertex v = mesh->vertex(idx);
            if(!v.degree()) return -1;
        }
        return 1;
    }
    return 0;
}

/*
 * Compute the set of simplices contained in the boundary bd(S) of the selected subset S of simplices.
 *
 * Input: A MeshSubset object containing the indices of the currently active vertices, edges, and faces, respectively.
 * Returns: The boundary of the given subset.
 */
MeshSubset SimplicialComplexOperators::boundary(const MeshSubset& subset) const {
    geometry->requireEdgeIndices();
    geometry->requireFaceIndices();
    MeshSubset cl = closure(subset);
    MeshSubset st;
    for(size_t idx : cl.vertices)
    {
        Vertex v = mesh->vertex(idx);
        size_t deg = 0;
        for(Edge e : v.adjacentEdges())
        {
            if(cl.edges.count(geometry->edgeIndices[e]))
            {
                deg++;
                if(deg > 1) break;
            }
        }
        if(deg == 1) st.addVertex(idx);
    }
    for(size_t idx : cl.edges)
    {
        Edge e = mesh->edge(idx);
        size_t deg = 0;
        for(Face f : e.adjacentFaces())
        {
            if(cl.faces.count(geometry->faceIndices[f]))
            {
                deg++;
                if(deg > 1) break;
            }
        }
        if(deg == 1) st.addEdge(idx);
    }
    return closure(st);
}