OpenGL 三角网格简化：QEM 算法实现
#include <GL/glut.h>
#include
using namespace std;
struct Vertex { float x, y, z; };
struct Triangle { int v1, v2, v3; };
vector
void loadMesh(const string& filename) { ifstream file(filename); if (!file.is_open()) { cout << "Failed to open file: " << filename << endl; return; }
string line;
while (getline(file, line)) {
if (line.empty())
continue;
stringstream ss(line);
char type;
ss >> type;
if (type == 'v') {
Vertex vertex;
ss >> vertex.x >> vertex.y >> vertex.z;
vertices.push_back(vertex);
}
else if (type == 'f') {
Triangle triangle;
ss >> triangle.v1 >> triangle.v2 >> triangle.v3;
triangles.push_back(triangle);
}
}
file.close();
}
void computeQuadricError(Vertex& v) { v.x = v.y = v.z = 0.0f; for (const auto& triangle : triangles) { const Vertex& v1 = vertices[triangle.v1 - 1]; const Vertex& v2 = vertices[triangle.v2 - 1]; const Vertex& v3 = vertices[triangle.v3 - 1];
float a = v2.y - v1.y;
float b = v1.x - v2.x;
float c = -(a * v1.x + b * v1.y);
float d = sqrt(a * a + b * b);
float distance = (a * v.x + b * v.y + c) / d;
if (distance > 0) {
v.x += (distance * a / d);
v.y += (distance * b / d);
v.z += (distance * c / d);
}
}
}
void simplifyMesh(float reductionRate) { int targetNumVertices = (1 - reductionRate) * vertices.size();
while (vertices.size() > targetNumVertices) {
vector<float> errors(vertices.size(), 0.0f);
vector<int> collapseMap(vertices.size(), -1);
for (int i = 0; i < triangles.size(); i++) {
const auto& triangle = triangles[i];
const Vertex& v1 = vertices[triangle.v1 - 1];
const Vertex& v2 = vertices[triangle.v2 - 1];
const Vertex& v3 = vertices[triangle.v3 - 1];
float q[4][4] = {
{ v1.x * v1.x, v1.x * v1.y, v1.x * v1.z, v1.x },
{ v1.x * v1.y, v1.y * v1.y, v1.y * v1.z, v1.y },
{ v1.x * v1.z, v1.y * v1.z, v1.z * v1.z, v1.z },
{ v1.x, v1.y, v1.z, 1 }
};
float determinant = q[0][0] * (q[1][1] * q[2][2] - q[1][2] * q[2][1]) -
q[0][1] * (q[1][0] * q[2][2] - q[2][0] * q[1][2]) +
q[0][2] * (q[1][0] * q[2][1] - q[2][0] * q[1][1]);
float qInverse[4][4] = {
{ (q[1][1] * q[2][2] - q[1][2] * q[2][1]) / determinant, (q[0][2] * q[2][1] - q[0][1] * q[2][2]) / determinant, (q[0][1] * q[1][2] - q[0][2] * q[1][1]) / determinant, (q[0][1] * q[1][2] * q[2][3] - q[0][2] * q[1][1] * q[2][3] - q[0][1] * q[1][3] * q[2][2] + q[0][3] * q[1][1] * q[2][2] + q[0][2] * q[1][3] * q[2][1] - q[0][3] * q[1][2] * q[2][1]) / determinant },
{ (q[1][2] * q[2][0] - q[1][0] * q[2][2]) / determinant, (q[0][0] * q[2][2] - q[0][2] * q[2][0]) / determinant, (q[0][2] * q[1][0] - q[0][0] * q[1][2]) / determinant, (q[0][0] * q[1][2] * q[2][3] - q[0][2] * q[1][0] * q[2][3] + q[0][0] * q[1][3] * q[2][2] - q[0][3] * q[1][0] * q[2][2] - q[0][2] * q[1][3] * q[2][0] + q[0][3] * q[1][2] * q[2][0]) / determinant },
{ (q[1][0] * q[2][1] - q[1][1] * q[2][0]) / determinant, (q[0][1] * q[2][0] - q[0][0] * q[2][1]) / determinant, (q[0][0] * q[1][1] - q[0][1] * q[1][0]) / determinant, (q[0][1] * q[1][0] * q[2][3] - q[0][0] * q[1][1] * q[2][3] - q[0][0] * q[1][3] * q[2][1] + q[0][3] * q[1][0] * q[2][1] + q[0][1] * q[1][3] * q[2][0] - q[0][3] * q[1][1] * q[2][0]) / determinant },
{ 0, 0, 0, 1 }
};
for (int j = 0; j < 3; j++) {
int vIndex = triangle.v1;
if (j == 1)
vIndex = triangle.v2;
else if (j == 2)
vIndex = triangle.v3;
Vertex& v = vertices[vIndex - 1];
float qe = v.x * (qInverse[0][0] * v.x + qInverse[1][0] * v.y + qInverse[2][0] * v.z + qInverse[3][0]) +
v.y * (qInverse[0][1] * v.x + qInverse[1][1] * v.y + qInverse[2][1] * v.z + qInverse[3][1]) +
v.z * (qInverse[0][2] * v.x + qInverse[1][2] * v.y + qInverse[2][2] * v.z + qInverse[3][2]) +
1 * (qInverse[0][3] * v.x + qInverse[1][3] * v.y + qInverse[2][3] * v.z + qInverse[3][3]);
errors[vIndex - 1] += qe;
}
}
float minError = numeric_limits<float>::max();
int collapseVertex = -1;
for (int i = 0; i < vertices.size(); i++) {
if (collapseMap[i] != -1)
continue;
if (errors[i] < minError) {
minError = errors[i];
collapseVertex = i;
}
}
if (collapseVertex == -1)
break;
Vertex& v1 = vertices[collapseVertex];
Vertex v2;
float minCost = numeric_limits<float>::max();
for (int i = 0; i < vertices.size(); i++) {
if (i == collapseVertex || collapseMap[i] != -1)
continue;
Vertex& v = vertices[i];
Vertex temp = v1;
float q[4][4] = {
{ v.x * v.x, v.x * v.y, v.x * v.z, v.x },
{ v.x * v.y, v.y * v.y, v.y * v.z, v.y },
{ v.x * v.z, v.y * v.z, v.z * v.z, v.z },
{ v.x, v.y, v.z, 1 }
};
float determinant = q[0][0] * (q[1][1] * q[2][2] - q[1][2] * q[2][1]) -
q[0][1] * (q[1][0] * q[2][2] - q[2][0] * q[1][2]) +
q[0][2] * (q[1][0] * q[2][1] - q[2][0] * q[1][1]);
float qInverse[4][4] = {
{ (q[1][1] * q[2][2] - q[1][2] * q[2][1]) / determinant, (q[0][2] * q[2][1] - q[0][1] * q[2][2]) / determinant, (q[0][1] * q[1][2] - q[0][2] * q[1][1]) / determinant, (q[0][1] * q[1][2] * q[2][3] - q[0][2] * q[1][1] * q[2][3] - q[0][1] * q[1][3] * q[2][2] + q[0][3] * q[1][1] * q[2][2] + q[0][2] * q[1][3] * q[2][1] - q[0][3] * q[1][2] * q[2][1]) / determinant },
{ (q[1][2] * q[2][0] - q[1][0] * q[2][2]) / determinant, (q[0][0] * q[2][2] - q[0][2] * q[2][0]) / determinant, (q[0][2] * q[1][0] - q[0][0] * q[1][2]) / determinant, (q[0][0] * q[1][2] * q[2][3] - q[0][2] * q[1][0] * q[2][3] + q[0][0] * q[1][3] * q[2][2] - q[0][3] * q[1][0] * q[2][2] - q[0][2] * q[1][3] * q[2][0] + q[0][3] * q[1][2] * q[2][0]) / determinant },
{ (q[1][0] * q[2][1] - q[1][1] * q[2][0]) / determinant, (q[0][1] * q[2][0] - q[0][0] * q[2][1]) / determinant, (q[0][0] * q[1][1] - q[0][1] * q[1][0]) / determinant, (q[0][1] * q[1][0] * q[2][3] - q[0][0] * q[1][1] * q[2][3] - q[0][0] * q[1][3] * q[2][1] + q[0][3] * q[1][0] * q[2][1] + q[0][1] * q[1][3] * q[2][0] - q[0][3] * q[1][1] * q[2][0]) / determinant },
{ 0, 0, 0, 1 }
};
for (const auto& triangle : triangles) {
if (triangle.v1 == collapseVertex + 1) {
v2 = vertices[triangle.v2 - 1];
break;
}
else if (triangle.v2 == collapseVertex + 1) {
v2 = vertices[triangle.v3 - 1];
break;
}
else if (triangle.v3 == collapseVertex + 1) {
v2 = vertices[triangle.v1 - 1];
break;
}
}
float e = v2.x * (qInverse[0][0] * v2.x + qInverse[1][0] * v2.y + qInverse[2][0] * v2.z + qInverse[3][0]) +
v2.y * (qInverse[0][1] * v2.x + qInverse[1][1] * v2.y + qInverse[2][1] * v2.z + qInverse[3][1]) +
v2.z * (qInverse[0][2] * v2.x + qInverse[1][2] * v2.y + qInverse[2][2] * v2.z + qInverse[3][2]) +
1 * (qInverse[0][3] * v2.x + qInverse[1][3] * v2.y + qInverse[2][3] * v2.z + qInverse[3][3]);
float d = v1.x * (qInverse[0][0] * v1.x + qInverse[1][0] * v1.y + qInverse[2][0] * v1.z + qInverse[3][0]) +
v1.y * (qInverse[0][1] * v1.x + qInverse[1][1] * v1.y + qInverse[2][1] * v1.z + qInverse[3][1]) +
v1.z * (qInverse[0][2] * v1.x + qInverse[1][2] * v1.y + qInverse[2][2] * v1.z + qInverse[3][2]) +
1 * (qInverse[0][3] * v1.x + qInverse[1][3] * v1.y + qInverse[2][3] * v1.z + qInverse[3][3]);
float cost = e - (2 * d) + minError;
if (cost < minCost) {
minCost = cost;
temp = v;
}
}
collapseMap[collapseVertex] = temp.x;
for (auto& triangle : triangles) {
if (triangle.v1 == collapseVertex + 1)
triangle.v1 = temp.x + 1;
if (triangle.v2 == collapseVertex + 1)
triangle.v2 = temp.x + 1;
if (triangle.v3 == collapseVertex + 1)
triangle.v3 = temp.x + 1;
}
}
vector<Vertex> newVertices;
for (int i = 0; i < vertices.size(); i++) {
if (vertices[i].z != -1) {
vertices[i].z = -1;
newVertices.push_back(vertices[i]);
}
}
vertices = newVertices;
triangles.clear();
}
void display() { glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
gluLookAt(0.0, 0.0, 5.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
glColor3f(1.0, 1.0, 1.0);
glBegin(GL_TRIANGLES);
for (const auto& triangle : triangles) {
const Vertex& v1 = vertices[triangle.v1 - 1];
const Vertex& v2 = vertices[triangle.v2 - 1];
const Vertex& v3 = vertices[triangle.v3 - 1];
glVertex3f(v1.x, v1.y, v1.z);
glVertex3f(v2.x, v2.y, v2.z);
glVertex3f(v3.x, v3.y, v3.z);
}
glEnd();
glutSwapBuffers();
}
void reshape(int width, int height) { glViewport(0, 0, width, height);
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
gluPerspective(60.0, (float)width / height, 1.0, 100.0);
}
void keyboard(unsigned char key, int x, int y) { if (key == 27) { exit(0); } }
int main(int argc, char** argv) { if (argc < 2) { cout << "Usage: " << argv[0] << " <mesh_file.obj>" << endl; return 0; }
string filename = argv[1];
loadMesh(filename);
float reductionRate;
cout << "Enter reduction rate (0~0.99): ";
cin >> reductionRate;
simplifyMesh(reductionRate);
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
glutInitWindowSize(800, 600);
glutCreateWindow("Triangle Mesh Simplification");
glClearColor(0.0, 0.0, 0.0, 1.0);
glEnable(GL_DEPTH_TEST);
glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutKeyboardFunc(keyboard);
glutMainLoop();
return 0;
}
原文地址: https://www.cveoy.top/t/topic/o8Fd 著作权归作者所有。请勿转载和采集!