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