vrshoot
view libs/assimp/Subdivision.cpp @ 0:b2f14e535253
initial commit
| author | John Tsiombikas <nuclear@member.fsf.org> |
|---|---|
| date | Sat, 01 Feb 2014 19:58:19 +0200 |
| parents | |
| children |
line source
1 /*
2 Open Asset Import Library (assimp)
3 ----------------------------------------------------------------------
5 Copyright (c) 2006-2012, assimp team
6 All rights reserved.
8 Redistribution and use of this software in source and binary forms,
9 with or without modification, are permitted provided that the
10 following conditions are met:
12 * Redistributions of source code must retain the above
13 copyright notice, this list of conditions and the
14 following disclaimer.
16 * Redistributions in binary form must reproduce the above
17 copyright notice, this list of conditions and the
18 following disclaimer in the documentation and/or other
19 materials provided with the distribution.
21 * Neither the name of the assimp team, nor the names of its
22 contributors may be used to endorse or promote products
23 derived from this software without specific prior
24 written permission of the assimp team.
26 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
27 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
28 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
29 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
30 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
31 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
32 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
33 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
34 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
35 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
36 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
38 ----------------------------------------------------------------------
39 */
41 #include "AssimpPCH.h"
43 #include "Subdivision.h"
44 #include "SceneCombiner.h"
45 #include "SpatialSort.h"
46 #include "ProcessHelper.h"
47 #include "Vertex.h"
49 using namespace Assimp;
50 void mydummy() {}
52 // ------------------------------------------------------------------------------------------------
53 /** Subdivider stub class to implement the Catmull-Clarke subdivision algorithm. The
54 * implementation is basing on recursive refinement. Directly evaluating the result is also
55 * possible and much quicker, but it depends on lengthy matrix lookup tables. */
56 // ------------------------------------------------------------------------------------------------
57 class CatmullClarkSubdivider : public Subdivider
58 {
60 public:
62 void Subdivide (aiMesh* mesh, aiMesh*& out, unsigned int num, bool discard_input);
63 void Subdivide (aiMesh** smesh, size_t nmesh,
64 aiMesh** out, unsigned int num, bool discard_input);
66 // ---------------------------------------------------------------------------
67 /** Intermediate description of an edge between two corners of a polygon*/
68 // ---------------------------------------------------------------------------
69 struct Edge
70 {
71 Edge()
72 : ref(0)
73 {}
74 Vertex edge_point, midpoint;
75 unsigned int ref;
76 };
80 typedef std::vector<unsigned int> UIntVector;
81 typedef std::map<uint64_t,Edge> EdgeMap;
83 // ---------------------------------------------------------------------------
84 // Hashing function to derive an index into an #EdgeMap from two given
85 // 'unsigned int' vertex coordinates (!!distinct coordinates - same
86 // vertex position == same index!!).
87 // NOTE - this leads to rare hash collisions if a) sizeof(unsigned int)>4
88 // and (id[0]>2^32-1 or id[0]>2^32-1).
89 // MAKE_EDGE_HASH() uses temporaries, so INIT_EDGE_HASH() needs to be put
90 // at the head of every function which is about to use MAKE_EDGE_HASH().
91 // Reason is that the hash is that hash construction needs to hold the
92 // invariant id0<id1 to identify an edge - else two hashes would refer
93 // to the same edge.
94 // ---------------------------------------------------------------------------
95 #define MAKE_EDGE_HASH(id0,id1) (eh_tmp0__=id0,eh_tmp1__=id1,\
96 (eh_tmp0__<eh_tmp1__?std::swap(eh_tmp0__,eh_tmp1__):mydummy()),(uint64_t)eh_tmp0__^((uint64_t)eh_tmp1__<<32u))
99 #define INIT_EDGE_HASH_TEMPORARIES()\
100 unsigned int eh_tmp0__, eh_tmp1__;
102 private:
104 void InternSubdivide (const aiMesh* const * smesh,
105 size_t nmesh,aiMesh** out, unsigned int num);
106 };
109 // ------------------------------------------------------------------------------------------------
110 // Construct a subdivider of a specific type
111 Subdivider* Subdivider::Create (Algorithm algo)
112 {
113 switch (algo)
114 {
115 case CATMULL_CLARKE:
116 return new CatmullClarkSubdivider();
117 };
119 ai_assert(false);
120 return NULL; // shouldn't happen
121 }
123 // ------------------------------------------------------------------------------------------------
124 // Call the Catmull Clark subdivision algorithm for one mesh
125 void CatmullClarkSubdivider::Subdivide (
126 aiMesh* mesh,
127 aiMesh*& out,
128 unsigned int num,
129 bool discard_input
130 )
131 {
132 assert(mesh != out);
133 Subdivide(&mesh,1,&out,num,discard_input);
134 }
136 // ------------------------------------------------------------------------------------------------
137 // Call the Catmull Clark subdivision algorithm for multiple meshes
138 void CatmullClarkSubdivider::Subdivide (
139 aiMesh** smesh,
140 size_t nmesh,
141 aiMesh** out,
142 unsigned int num,
143 bool discard_input
144 )
145 {
146 ai_assert(NULL != smesh && NULL != out);
148 // course, both regions may not overlap
149 assert(smesh<out || smesh+nmesh>out+nmesh);
150 if (!num) {
152 // No subdivision at all. Need to copy all the meshes .. argh.
153 if (discard_input) {
154 for (size_t s = 0; s < nmesh; ++s) {
155 out[s] = smesh[s];
156 smesh[s] = NULL;
157 }
158 }
159 else {
160 for (size_t s = 0; s < nmesh; ++s) {
161 SceneCombiner::Copy(out+s,smesh[s]);
162 }
163 }
164 return;
165 }
167 std::vector<aiMesh*> inmeshes;
168 std::vector<aiMesh*> outmeshes;
169 std::vector<unsigned int> maptbl;
171 inmeshes.reserve(nmesh);
172 outmeshes.reserve(nmesh);
173 maptbl.reserve(nmesh);
175 // Remove pure line and point meshes from the working set to reduce the
176 // number of edge cases the subdivider is forced to deal with. Line and
177 // point meshes are simply passed through.
178 for (size_t s = 0; s < nmesh; ++s) {
179 aiMesh* i = smesh[s];
180 // FIX - mPrimitiveTypes might not yet be initialized
181 if (i->mPrimitiveTypes && (i->mPrimitiveTypes & (aiPrimitiveType_LINE|aiPrimitiveType_POINT))==i->mPrimitiveTypes) {
182 DefaultLogger::get()->debug("Catmull-Clark Subdivider: Skipping pure line/point mesh");
184 if (discard_input) {
185 out[s] = i;
186 smesh[s] = NULL;
187 }
188 else {
189 SceneCombiner::Copy(out+s,i);
190 }
191 continue;
192 }
194 outmeshes.push_back(NULL);inmeshes.push_back(i);
195 maptbl.push_back(s);
196 }
198 // Do the actual subdivision on the preallocated storage. InternSubdivide
199 // *always* assumes that enough storage is available, it does not bother
200 // checking any ranges.
201 ai_assert(inmeshes.size()==outmeshes.size()&&inmeshes.size()==maptbl.size());
202 if (inmeshes.empty()) {
203 DefaultLogger::get()->warn("Catmull-Clark Subdivider: Pure point/line scene, I can't do anything");
204 return;
205 }
206 InternSubdivide(&inmeshes.front(),inmeshes.size(),&outmeshes.front(),num);
207 for (unsigned int i = 0; i < maptbl.size(); ++i) {
208 ai_assert(outmeshes[i]);
209 out[maptbl[i]] = outmeshes[i];
210 }
212 if (discard_input) {
213 for (size_t s = 0; s < nmesh; ++s) {
214 delete smesh[s];
215 }
216 }
217 }
219 // ------------------------------------------------------------------------------------------------
220 // Note - this is an implementation of the standard (recursive) Cm-Cl algorithm without further
221 // optimizations (except we're using some nice LUTs). A description of the algorithm can be found
222 // here: http://en.wikipedia.org/wiki/Catmull-Clark_subdivision_surface
223 //
224 // The code is mostly O(n), however parts are O(nlogn) which is therefore the algorithm's
225 // expected total runtime complexity. The implementation is able to work in-place on the same
226 // mesh arrays. Calling #InternSubdivide() directly is not encouraged. The code can operate
227 // in-place unless 'smesh' and 'out' are equal (no strange overlaps or reorderings).
228 // Previous data is replaced/deleted then.
229 // ------------------------------------------------------------------------------------------------
230 void CatmullClarkSubdivider::InternSubdivide (
231 const aiMesh* const * smesh,
232 size_t nmesh,
233 aiMesh** out,
234 unsigned int num
235 )
236 {
237 ai_assert(NULL != smesh && NULL != out);
238 INIT_EDGE_HASH_TEMPORARIES();
240 // no subdivision requested or end of recursive refinement
241 if (!num) {
242 return;
243 }
245 UIntVector maptbl;
246 SpatialSort spatial;
248 // ---------------------------------------------------------------------
249 // 0. Offset table to index all meshes continuously, generate a spatially
250 // sorted representation of all vertices in all meshes.
251 // ---------------------------------------------------------------------
252 typedef std::pair<unsigned int,unsigned int> IntPair;
253 std::vector<IntPair> moffsets(nmesh);
254 unsigned int totfaces = 0, totvert = 0;
255 for (size_t t = 0; t < nmesh; ++t) {
256 const aiMesh* mesh = smesh[t];
258 spatial.Append(mesh->mVertices,mesh->mNumVertices,sizeof(aiVector3D),false);
259 moffsets[t] = IntPair(totfaces,totvert);
261 totfaces += mesh->mNumFaces;
262 totvert += mesh->mNumVertices;
263 }
265 spatial.Finalize();
266 const unsigned int num_unique = spatial.GenerateMappingTable(maptbl,ComputePositionEpsilon(smesh,nmesh));
269 #define FLATTEN_VERTEX_IDX(mesh_idx, vert_idx) (moffsets[mesh_idx].second+vert_idx)
270 #define FLATTEN_FACE_IDX(mesh_idx, face_idx) (moffsets[mesh_idx].first+face_idx)
272 // ---------------------------------------------------------------------
273 // 1. Compute the centroid point for all faces
274 // ---------------------------------------------------------------------
275 std::vector<Vertex> centroids(totfaces);
276 unsigned int nfacesout = 0;
277 for (size_t t = 0, n = 0; t < nmesh; ++t) {
278 const aiMesh* mesh = smesh[t];
279 for (unsigned int i = 0; i < mesh->mNumFaces;++i,++n)
280 {
281 const aiFace& face = mesh->mFaces[i];
282 Vertex& c = centroids[n];
284 for (unsigned int a = 0; a < face.mNumIndices;++a) {
285 c += Vertex(mesh,face.mIndices[a]);
286 }
288 c /= static_cast<float>(face.mNumIndices);
289 nfacesout += face.mNumIndices;
290 }
291 }
293 EdgeMap edges;
295 // ---------------------------------------------------------------------
296 // 2. Set each edge point to be the average of all neighbouring
297 // face points and original points. Every edge exists twice
298 // if there is a neighboring face.
299 // ---------------------------------------------------------------------
300 for (size_t t = 0; t < nmesh; ++t) {
301 const aiMesh* mesh = smesh[t];
303 for (unsigned int i = 0; i < mesh->mNumFaces;++i) {
304 const aiFace& face = mesh->mFaces[i];
306 for (unsigned int p =0; p< face.mNumIndices; ++p) {
307 const unsigned int id[] = {
308 face.mIndices[p],
309 face.mIndices[p==face.mNumIndices-1?0:p+1]
310 };
311 const unsigned int mp[] = {
312 maptbl[FLATTEN_VERTEX_IDX(t,id[0])],
313 maptbl[FLATTEN_VERTEX_IDX(t,id[1])]
314 };
316 Edge& e = edges[MAKE_EDGE_HASH(mp[0],mp[1])];
317 e.ref++;
318 if (e.ref<=2) {
319 if (e.ref==1) { // original points (end points) - add only once
320 e.edge_point = e.midpoint = Vertex(mesh,id[0])+Vertex(mesh,id[1]);
321 e.midpoint *= 0.5f;
322 }
323 e.edge_point += centroids[FLATTEN_FACE_IDX(t,i)];
324 }
325 }
326 }
327 }
329 // ---------------------------------------------------------------------
330 // 3. Normalize edge points
331 // ---------------------------------------------------------------------
332 {unsigned int bad_cnt = 0;
333 for (EdgeMap::iterator it = edges.begin(); it != edges.end(); ++it) {
334 if ((*it).second.ref < 2) {
335 ai_assert((*it).second.ref);
336 ++bad_cnt;
337 }
338 (*it).second.edge_point *= 1.f/((*it).second.ref+2.f);
339 }
341 if (bad_cnt) {
342 // Report the number of bad edges. bad edges are referenced by less than two
343 // faces in the mesh. They occur at outer model boundaries in non-closed
344 // shapes.
345 char tmp[512];
346 sprintf(tmp,"Catmull-Clark Subdivider: got %u bad edges touching only one face (totally %u edges). ",
347 bad_cnt,static_cast<unsigned int>(edges.size()));
349 DefaultLogger::get()->debug(tmp);
350 }}
352 // ---------------------------------------------------------------------
353 // 4. Compute a vertex-face adjacency table. We can't reuse the code
354 // from VertexTriangleAdjacency because we need the table for multiple
355 // meshes and out vertex indices need to be mapped to distinct values
356 // first.
357 // ---------------------------------------------------------------------
358 UIntVector faceadjac(nfacesout), cntadjfac(maptbl.size(),0), ofsadjvec(maptbl.size()+1,0); {
359 for (size_t t = 0; t < nmesh; ++t) {
360 const aiMesh* const minp = smesh[t];
361 for (unsigned int i = 0; i < minp->mNumFaces; ++i) {
363 const aiFace& f = minp->mFaces[i];
364 for (unsigned int n = 0; n < f.mNumIndices; ++n) {
365 ++cntadjfac[maptbl[FLATTEN_VERTEX_IDX(t,f.mIndices[n])]];
366 }
367 }
368 }
369 unsigned int cur = 0;
370 for (size_t i = 0; i < cntadjfac.size(); ++i) {
371 ofsadjvec[i+1] = cur;
372 cur += cntadjfac[i];
373 }
374 for (size_t t = 0; t < nmesh; ++t) {
375 const aiMesh* const minp = smesh[t];
376 for (unsigned int i = 0; i < minp->mNumFaces; ++i) {
378 const aiFace& f = minp->mFaces[i];
379 for (unsigned int n = 0; n < f.mNumIndices; ++n) {
380 faceadjac[ofsadjvec[1+maptbl[FLATTEN_VERTEX_IDX(t,f.mIndices[n])]]++] = FLATTEN_FACE_IDX(t,i);
381 }
382 }
383 }
385 // check the other way round for consistency
386 #ifdef _DEBUG
388 for (size_t t = 0; t < ofsadjvec.size()-1; ++t) {
389 for (unsigned int m = 0; m < cntadjfac[t]; ++m) {
390 const unsigned int fidx = faceadjac[ofsadjvec[t]+m];
391 ai_assert(fidx < totfaces);
392 for (size_t n = 1; n < nmesh; ++n) {
394 if (moffsets[n].first > fidx) {
395 const aiMesh* msh = smesh[--n];
396 const aiFace& f = msh->mFaces[fidx-moffsets[n].first];
398 bool haveit = false;
399 for (unsigned int i = 0; i < f.mNumIndices; ++i) {
400 if (maptbl[FLATTEN_VERTEX_IDX(n,f.mIndices[i])]==(unsigned int)t) {
401 haveit = true; break;
402 }
403 }
404 ai_assert(haveit);
405 break;
406 }
407 }
408 }
409 }
411 #endif
412 }
414 #define GET_ADJACENT_FACES_AND_CNT(vidx,fstartout,numout) \
415 fstartout = &faceadjac[ofsadjvec[vidx]], numout = cntadjfac[vidx]
417 typedef std::pair<bool,Vertex> TouchedOVertex;
418 std::vector<TouchedOVertex > new_points(num_unique,TouchedOVertex(false,Vertex()));
419 // ---------------------------------------------------------------------
420 // 5. Spawn a quad from each face point to the corresponding edge points
421 // the original points being the fourth quad points.
422 // ---------------------------------------------------------------------
423 for (size_t t = 0; t < nmesh; ++t) {
424 const aiMesh* const minp = smesh[t];
425 aiMesh* const mout = out[t] = new aiMesh();
427 for (unsigned int a = 0; a < minp->mNumFaces; ++a) {
428 mout->mNumFaces += minp->mFaces[a].mNumIndices;
429 }
431 // We need random access to the old face buffer, so reuse is not possible.
432 mout->mFaces = new aiFace[mout->mNumFaces];
434 mout->mNumVertices = mout->mNumFaces*4;
435 mout->mVertices = new aiVector3D[mout->mNumVertices];
437 // quads only, keep material index
438 mout->mPrimitiveTypes = aiPrimitiveType_POLYGON;
439 mout->mMaterialIndex = minp->mMaterialIndex;
441 if (minp->HasNormals()) {
442 mout->mNormals = new aiVector3D[mout->mNumVertices];
443 }
445 if (minp->HasTangentsAndBitangents()) {
446 mout->mTangents = new aiVector3D[mout->mNumVertices];
447 mout->mBitangents = new aiVector3D[mout->mNumVertices];
448 }
450 for(unsigned int i = 0; minp->HasTextureCoords(i); ++i) {
451 mout->mTextureCoords[i] = new aiVector3D[mout->mNumVertices];
452 mout->mNumUVComponents[i] = minp->mNumUVComponents[i];
453 }
455 for(unsigned int i = 0; minp->HasVertexColors(i); ++i) {
456 mout->mColors[i] = new aiColor4D[mout->mNumVertices];
457 }
459 mout->mNumVertices = mout->mNumFaces<<2u;
460 for (unsigned int i = 0, v = 0, n = 0; i < minp->mNumFaces;++i) {
462 const aiFace& face = minp->mFaces[i];
463 for (unsigned int a = 0; a < face.mNumIndices;++a) {
465 // Get a clean new face.
466 aiFace& faceOut = mout->mFaces[n++];
467 faceOut.mIndices = new unsigned int [faceOut.mNumIndices = 4];
469 // Spawn a new quadrilateral (ccw winding) for this original point between:
470 // a) face centroid
471 centroids[FLATTEN_FACE_IDX(t,i)].SortBack(mout,faceOut.mIndices[0]=v++);
473 // b) adjacent edge on the left, seen from the centroid
474 const Edge& e0 = edges[MAKE_EDGE_HASH(maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])],
475 maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a==face.mNumIndices-1?0:a+1])
476 ])]; // fixme: replace with mod face.mNumIndices?
478 // c) adjacent edge on the right, seen from the centroid
479 const Edge& e1 = edges[MAKE_EDGE_HASH(maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])],
480 maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[!a?face.mNumIndices-1:a-1])
481 ])]; // fixme: replace with mod face.mNumIndices?
483 e0.edge_point.SortBack(mout,faceOut.mIndices[3]=v++);
484 e1.edge_point.SortBack(mout,faceOut.mIndices[1]=v++);
486 // d= original point P with distinct index i
487 // F := 0
488 // R := 0
489 // n := 0
490 // for each face f containing i
491 // F := F+ centroid of f
492 // R := R+ midpoint of edge of f from i to i+1
493 // n := n+1
494 //
495 // (F+2R+(n-3)P)/n
496 const unsigned int org = maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])];
497 TouchedOVertex& ov = new_points[org];
499 if (!ov.first) {
500 ov.first = true;
502 const unsigned int* adj; unsigned int cnt;
503 GET_ADJACENT_FACES_AND_CNT(org,adj,cnt);
505 if (cnt < 3) {
506 ov.second = Vertex(minp,face.mIndices[a]);
507 }
508 else {
510 Vertex F,R;
511 for (unsigned int o = 0; o < cnt; ++o) {
512 ai_assert(adj[o] < totfaces);
513 F += centroids[adj[o]];
515 // adj[0] is a global face index - search the face in the mesh list
516 const aiMesh* mp = NULL;
517 size_t nidx;
519 if (adj[o] < moffsets[0].first) {
520 mp = smesh[nidx=0];
521 }
522 else {
523 for (nidx = 1; nidx<= nmesh; ++nidx) {
524 if (nidx == nmesh ||moffsets[nidx].first > adj[o]) {
525 mp = smesh[--nidx];
526 break;
527 }
528 }
529 }
531 ai_assert(adj[o]-moffsets[nidx].first < mp->mNumFaces);
532 const aiFace& f = mp->mFaces[adj[o]-moffsets[nidx].first];
533 # ifdef _DEBUG
534 bool haveit = false;
535 # endif
537 // find our original point in the face
538 for (unsigned int m = 0; m < f.mNumIndices; ++m) {
539 if (maptbl[FLATTEN_VERTEX_IDX(nidx,f.mIndices[m])] == org) {
541 // add *both* edges. this way, we can be sure that we add
542 // *all* adjacent edges to R. In a closed shape, every
543 // edge is added twice - so we simply leave out the
544 // factor 2.f in the amove formula and get the right
545 // result.
547 const Edge& c0 = edges[MAKE_EDGE_HASH(org,maptbl[FLATTEN_VERTEX_IDX(
548 nidx,f.mIndices[!m?f.mNumIndices-1:m-1])])];
549 // fixme: replace with mod face.mNumIndices?
551 const Edge& c1 = edges[MAKE_EDGE_HASH(org,maptbl[FLATTEN_VERTEX_IDX(
552 nidx,f.mIndices[m==f.mNumIndices-1?0:m+1])])];
553 // fixme: replace with mod face.mNumIndices?
554 R += c0.midpoint+c1.midpoint;
556 # ifdef _DEBUG
557 haveit = true;
558 # endif
559 break;
560 }
561 }
563 // this invariant *must* hold if the vertex-to-face adjacency table is valid
564 ai_assert(haveit);
565 }
567 const float div = static_cast<float>(cnt), divsq = 1.f/(div*div);
568 ov.second = Vertex(minp,face.mIndices[a])*((div-3.f) / div) + R*divsq + F*divsq;
569 }
570 }
571 ov.second.SortBack(mout,faceOut.mIndices[2]=v++);
572 }
573 }
574 }
576 // ---------------------------------------------------------------------
577 // 7. Apply the next subdivision step.
578 // ---------------------------------------------------------------------
579 if (num != 1) {
580 std::vector<aiMesh*> tmp(nmesh);
581 InternSubdivide (out,nmesh,&tmp.front(),num-1);
582 for (size_t i = 0; i < nmesh; ++i) {
583 delete out[i];
584 out[i] = tmp[i];
585 }
586 }
587 }