QEM三维模型抽稀算法C++Wraper实例

去年研究了QEM算法，下面简要记录下研究过程

一，QEM三维模型简化算法，介绍见

1.QEM(Quadric Error Mactrics，二次误差测度)模型简化算法，具有兼顾执行效率和模型质量的优点（https://www.cnblogs.com/kekec/archive/2011/12/07/2279320.html）

2.C++源码：http://mgarland.org/software/qslim.html

二，Wraper代码实现

新建CLR C++项目“QSlimWraper”，包含Mesh.h，QSlim.h，QSlim.cpp代码

新建Wraper.cpp，实现qslim_init()，qslim_run()函数，见下面

 1 #include "Mesh.h"
 2 #include "QSlim.h"  
 3  
 4 
 5 using namespace System::Runtime::InteropServices;
 6 static void qslim_init()
 7 {
 8     int i;
 9     cerr << "Reading input ..." << endl;
10     cerr << "Cleaning up initial input ..." << endl;
11     int initialVertCount = M0.vertCount();
12     int initialEdgeCount = M0.edgeCount();
13     int initialFaceCount = M0.faceCount();
14     for (i = 0; i < M0.faceCount(); i++)
15         if (!M0.face(i)->plane().isValid())
16             M0.killFace(M0.face(i));
17     M0.removeDegeneracy(M0.allFaces());
18     for (i = 0; i < M0.vertCount(); i++)
19     {
20         if (M0.vertex(i)->edgeUses().length() == 0)
21             M0.vertex(i)->kill();
22     }
23     cerr << "Input model summary:" << endl;
24     cerr << "    Vertices    : " << initialVertCount << endl;
25     cerr << "    Edges       : " << initialEdgeCount << endl;
26     int man = 0, non = 0, bndry = 0, bogus = 0;
27     for (i = 0; i < M0.edgeCount(); i++)
28         switch (M0.edge(i)->faceUses().length())
29         {
30         case 0:
31             bogus++;
32             break;
33         case 1:
34             bndry++;
35             break;
36         case 2:
37             man++;
38             break;
39         default:
40             non++;
41             break;
42         }
43     if (bogus)
44         cerr << "        Bogus       : " << bogus << endl;
45     cerr << "        Boundary    : " << bndry << endl;
46     cerr << "        Manifold    : " << man << endl;
47     cerr << "        Nonmanifold : " << non << endl;
48 
49     cerr << "    Faces       : " << initialFaceCount << endl;
50 }
51 
52 static void qslim_run()
53 {
54     decimate_init(M0, pair_selection_tolerance); 
55     while (M0.validFaceCount > face_target&& decimate_min_error() < error_tolerance)
56     {
57         decimate_contract(M0); 
58     }
59 }

View Code

新建QSlimWraper namespace，实现CLR，见下面

  1 namespace QSlimWraper
  2 {
  3     public ref struct Point3dWrap
  4     {
  5         float X;
  6         float Y;
  7         float Z;
  8         int R;//red
  9         int G;//green
 10         int B;//blue
 11         Point3dWrap()
 12         {
 13             X = 0;
 14             Y = 0;
 15             Z = 0;
 16         }
 17         ~Point3dWrap()
 18         {
 19         }
 20         Point3dWrap(float x, float y, float z)
 21         {
 22             this->X = x;
 23             this->Y = y;
 24             this->Z = z;
 25         }
 26         Point3dWrap(float x, float y, float z, int r, int g, int b)
 27         {
 28             this->X = x;
 29             this->Y = y;
 30             this->Z = z;
 31             this->R = r;
 32             this->G = g;
 33             this->B = b;
 34         }
 35     };
 36 
 37     public ref struct TriangleWrap
 38     {
 39         int pIndex0;
 40         int pIndex1;
 41         int pIndex2;
 42         int R;//red
 43         int G;//green
 44         int B;//blue
 45 
 46         TriangleWrap()
 47         {
 48             this->pIndex0 = -1;
 49             this->pIndex1 = -1;
 50             this->pIndex2 = -1;
 51         }
 52         ~TriangleWrap()
 53         {
 54         }
 55         TriangleWrap(int p0index, int p1index, int p2index)
 56         {
 57             this->pIndex0 = p0index;
 58             this->pIndex1 = p1index;
 59             this->pIndex2 = p2index;
 60         }
 61         TriangleWrap(int p0index, int p1index, int p2index, int r, int g, int b)
 62         {
 63             this->pIndex0 = p0index;
 64             this->pIndex1 = p1index;
 65             this->pIndex2 = p2index;
 66             this->R = r;
 67             this->G = g;
 68             this->B = b;
 69         }
 70     };
 71 
 72     public ref class QSlimWrap
 73     {
 74     public:
 75         QSlimWrap();
 76         ~QSlimWrap();
 77         int ReadFile(char* file);
 78         int OutputFile(char* newFile);
 79 
 80         int InportData(cli::array<Point3dWrap^>^ vertices, cli::array<TriangleWrap^>^ faces, int fTarget);
 81         int OutputData([Out] cli::array<Point3dWrap^>^% vertices, [Out] cli::array<TriangleWrap^>^% faces);
 82 
 83         int TestParam1(cli::array<Point3dWrap^>^ datas);
 84         bool Uninitialize();
 85          
 86 
 87     private:
 88         Mesh *mesh;
 89         int faceTarget;
 90         void InitM0();
 91         void ReplaceM();
 92         int Run();
 93     };
 94 
 95     QSlimWrap::QSlimWrap()
 96     {
 97         mesh = new Mesh;
 98     }
 99 
100     QSlimWrap::~QSlimWrap()
101     {
102         Uninitialize();
103     }
104 
105     bool QSlimWrap::Uninitialize()
106     { 
107         DisposeQS();
108         if (mesh)
109         {
110             delete mesh;
111             mesh = NULL;
112         }
113         else
114             return false;
115         return true;
116     }
117 
118     int QSlimWrap::TestParam1(cli::array<Point3dWrap^>^ datas)
119     {
120         return 12;
121     }
122 
123     int QSlimWrap::ReadFile(char* file)
124     {
125         Mesh m;
126         PlyManager::ReadFile(m, file);//result_decix
127         mesh->Vertices = m.Vertices;
128         mesh->Faces = m.Faces;
129         faceTarget = 2 * m.Faces.size() / 3;
130         Run();
131         return 1;
132     }
133 
134     int QSlimWrap::OutputFile(char* newFile)
135     {
136         Mesh m;
137         m.Vertices = mesh->Vertices;
138         m.Faces = mesh->Faces;
139         PlyManager::Output(m, newFile);
140         cerr << "Output File:" << newFile << endl;
141         return 0;
142     }
143 
144     void QSlimWrap::InitM0()
145     { 
146         for (size_t i = 0;i < mesh->Vertices.size();i++)
147         {
148             Point3d p = mesh->Vertices[i];
149             Vec3 v(p.X, p.Y, p.Z);
150             //M0.in_Vertex(v);
151             M0.in_VertexEx(v, p.R, p.G, p.B);
152         }
153         for (size_t i = 0;i < mesh->Faces.size();i++)
154         {
155             Triangle t = mesh->Faces[i];
156             //M0.in_Face(t.P0Index,t.P1Index,t.P2Index);
157             M0.in_FaceEx(t.P0Index, t.P1Index, t.P2Index, t.R, t.G, t.B);
158         } 
159     }
160 
161     void QSlimWrap::ReplaceM()
162     { 
163         mesh->Vertices.swap(std::vector<Point3d>());
164         mesh->Faces.swap(std::vector<Triangle>());
165         mesh->Vertices.reserve(M0.vertCount());
166         mesh->Faces.reserve(M0.faceCount());
167         int* map = new int[M0.vertCount()];
168         for (int i = 0;i < M0.vertCount();i++)
169             map[i] = -1;
170         int number = 0, number1 = 0;
171         for (int i = 0;i < M0.vertCount();i++)
172         { 
173             if (M0.vertex(i)->isValid())
174             { 
175                 real* data = M0.vertex(i)->raw();
176                 //Point3d p((float)data[0], (float)data[1], (float)data[2]);
177                 Point3d p((float)data[0], (float)data[1], (float)data[2], M0.vertex(i)->R, M0.vertex(i)->G, M0.vertex(i)->B);
178                 map[i] = mesh->AddVertex(p);
179                 number++;
180             }
181             else 
182                 number1++;
183         }
184          cerr << "ReplaceM M0 Vertices valid  number: " << number << "  not valid  number : " << number1 << "  total Count:" << M0.vertCount()<< endl;
185          number=0, number1=0;
186         for (int i = 0;i < M0.faceCount();i++)
187         {
188             if (M0.face(i)->isValid())
189             {
190                 Vertex* v0 = M0.face(i)->vertex(0);
191                 Vertex* v1 = M0.face(i)->vertex(1);
192                 Vertex* v2 = M0.face(i)->vertex(2);
193                 //Triangle t(map[v0->uniqID], map[v1->uniqID], map[v2->uniqID]); 
194                 Triangle t(map[v0->uniqID], map[v1->uniqID], map[v2->uniqID], M0.face(i)->R, M0.face(i)->G, M0.face(i)->B);
195                 mesh->AddFace(t);
196                 number++;
197             }
198             else
199                 number1++;
200         } 
201         cerr << "ReplaceM M0 faces valid  number: " << number << "  not valid  number : " << number1 << "  total Count:" << M0.faceCount() << endl;
202         delete[] map;
203     }
204 
205     int QSlimWrap::Run()
206     {
207         InitM0();
208         //cerr << "InportData runing    fTarget =" << face_target << endl;
209         //qslim_init();
210         face_target = faceTarget;//2 * m.Faces.size() / 3;
211         error_tolerance = HUGE;
212         will_use_plane_constraint = true;
213         will_use_vertex_constraint = false;
214         will_preserve_boundaries = true;
215         will_preserve_mesh_quality = true;
216         will_constrain_boundaries = true;
217         boundary_constraint_weight = 1.0;
218         will_weight_by_area = false;
219         placement_policy = 1;
220         pair_selection_tolerance = 0.0;
221         qslim_run();
222         ReplaceM();
223         return 1;
224     }
225 
226     int QSlimWrap::InportData(cli::array<Point3dWrap^>^ vertices, cli::array<TriangleWrap^>^ faces, int fTarget)
227     { 
228         int vCount = vertices->LongLength;
229         int fCount = faces->LongLength;
230         faceTarget = fTarget;
231         float v1 = 0, v2 = 0, v3 = 0;
232         int r = 0, g = 0, b = 0;
233         long i1 = 0, i2 = 0, i3 = 0;
234 
235         for (long i = 0;i < vCount;i++)
236         {
237             Point3dWrap^ ver = (Point3dWrap^)vertices[i];
238             Point3d p3d(ver->X, ver->Y, ver->Z, ver->R, ver->G, ver->B);
239             mesh->AddVertex(p3d); 
240         }
241         delete[] vertices;
242 
243         for (long j = 0;j < fCount;j++)
244         {
245             TriangleWrap^ tg = faces[j];
246             Triangle t(tg->pIndex0, tg->pIndex1, tg->pIndex2, tg->R, tg->G, tg->B);
247             mesh->AddFace(t); 
248         }
249         delete[] faces;
250         cerr << "InportData run before Vertices Count:" << mesh->Vertices.size() << "  face_target:" << faceTarget << endl;
251         Run();
252         cerr << "InportData run completed Vertices Count:" << mesh->Vertices.size() << "  faces Count: "<< mesh->Faces.size() <<endl; 
253         return 1;
254     }
255      
256 
257     int QSlimWrap::OutputData([Out] cli::array<Point3dWrap^>^% vertices, [Out]  cli::array<TriangleWrap^>^% faces)
258     { 
259         int vCount = mesh->Vertices.size(); 
260         vertices = gcnew cli::array<Point3dWrap^>(vCount);
261         for (int i = 0;i < vCount;i++)
262         {
263             Point3d v = mesh->Vertices[i];
264             Point3dWrap^ ver = gcnew Point3dWrap();
265             ver->X = v.X, ver->Y = v.Y, ver->Z = v.Z, ver->R = v.R, ver->G = v.G, ver->B = v.B;
266             vertices[i] = ver;
267         }
268         mesh->Vertices.clear();
269 
270         int fCount = mesh->Faces.size();
271         faces = gcnew cli::array<TriangleWrap^>(fCount);
272         for (int i = 0;i < fCount;i++)
273         {
274             Triangle t = mesh->Faces[i]; 
275             TriangleWrap^ tw = gcnew TriangleWrap(t.P0Index, t.P1Index, t.P2Index, t.R, t.G, t.B);
276             faces[i] = tw;
277         }
278         
279         mesh->Faces.clear();
280         cerr << "OutputData vertices Count:" << vCount << "  faces Count:" << fCount << endl;
281         return 1;
282     } 
283 }

View Code

三，Wraper验证

新建C# 项目，引用QSlimWraper.dll,验证代码如下：

 1 static void Main(string[] args)
 2         {
 3 
 4             QSlimWraper wrap = new QSlimWraper();
 5             unsafe
 6             {
 7                 //Point3d[] test = new Point3d[10];
 8                 //for (int i = 0; i < 10; i++)
 9                 //{
10                 //    test[i].X = i + 1;
11                 //}
12 
13                 string str = "E:\resultXXX2.ply";
14                 IntPtr intPtrStr = (IntPtr)Marshal.StringToHGlobalAnsi(str);
15                 sbyte* sbyteStr = (sbyte*)intPtrStr;
16                 int res = wrap.ReadFile(sbyteStr);
17 
18                 str = "E:\resultXXX2XXX.ply";
19                 intPtrStr = (IntPtr)Marshal.StringToHGlobalAnsi(str);
20                 sbyteStr = (sbyte*)intPtrStr;
21                 int res1 = wrap.OutputFile(sbyteStr);
22 
23                 //Point3dWrap[] pts = new Point3dWrap[100];
24                 //for (int i = 0; i < 100; i++)
25                 //{
26                 //    Point3dWrap aa = new Point3dWrap();
27                 //    aa.X = i + 1;
28                 //    pts[i] = aa;
29                 //}
30                 //TriangleWrap[] ts = new TriangleWrap[300];
31                 //for (int i = 0; i < 100; i++)
32                 //{
33                 //    TriangleWrap aa = new TriangleWrap();
34                 //    aa.pIndex0 = i;
35                 //    aa.pIndex1 = i + 100;
36                 //    aa.pIndex2 = i + 200;
37                 //    ts[i] = aa;
38                 //}
39                 ////int iii = wrap.TestParam1(test);
40                 //int face_target = 2 * ts.Length / 3;
41                 //wrap.InportData(pts, ts, face_target);
42                 //Point3dWrap[] point3Ds;
43                 //TriangleWrap[] triangles;
44                 //wrap.OutputData(out point3Ds, out triangles);
45                 wrap.Dispose();
46             }
47         }

View Code

四，示例效果，19800 surfaces to 1000 surfaces

如需要原码可联系