1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/assimp/Subdivision.cpp	Sat Feb 01 19:58:19 2014 +0200
     1.3 @@ -0,0 +1,587 @@
     1.4 +/*
     1.5 +Open Asset Import Library (assimp)
     1.6 +----------------------------------------------------------------------
     1.7 +
     1.8 +Copyright (c) 2006-2012, assimp team
     1.9 +All rights reserved.
    1.10 +
    1.11 +Redistribution and use of this software in source and binary forms, 
    1.12 +with or without modification, are permitted provided that the 
    1.13 +following conditions are met:
    1.14 +
    1.15 +* Redistributions of source code must retain the above
    1.16 +  copyright notice, this list of conditions and the
    1.17 +  following disclaimer.
    1.18 +
    1.19 +* Redistributions in binary form must reproduce the above
    1.20 +  copyright notice, this list of conditions and the
    1.21 +  following disclaimer in the documentation and/or other
    1.22 +  materials provided with the distribution.
    1.23 +
    1.24 +* Neither the name of the assimp team, nor the names of its
    1.25 +  contributors may be used to endorse or promote products
    1.26 +  derived from this software without specific prior
    1.27 +  written permission of the assimp team.
    1.28 +
    1.29 +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
    1.30 +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
    1.31 +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
    1.32 +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
    1.33 +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
    1.34 +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
    1.35 +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
    1.36 +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
    1.37 +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
    1.38 +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
    1.39 +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    1.40 +
    1.41 +----------------------------------------------------------------------
    1.42 +*/
    1.43 +
    1.44 +#include "AssimpPCH.h"
    1.45 +
    1.46 +#include "Subdivision.h"
    1.47 +#include "SceneCombiner.h"
    1.48 +#include "SpatialSort.h"
    1.49 +#include "ProcessHelper.h"
    1.50 +#include "Vertex.h"
    1.51 +
    1.52 +using namespace Assimp;
    1.53 +void mydummy() {}
    1.54 +
    1.55 +// ------------------------------------------------------------------------------------------------
    1.56 +/** Subdivider stub class to implement the Catmull-Clarke subdivision algorithm. The 
    1.57 + *  implementation is basing on recursive refinement. Directly evaluating the result is also
    1.58 + *  possible and much quicker, but it depends on lengthy matrix lookup tables. */
    1.59 +// ------------------------------------------------------------------------------------------------
    1.60 +class CatmullClarkSubdivider : public Subdivider
    1.61 +{
    1.62 +
    1.63 +public:
    1.64 +
    1.65 +	void Subdivide (aiMesh* mesh, aiMesh*& out, unsigned int num, bool discard_input);
    1.66 +	void Subdivide (aiMesh** smesh, size_t nmesh,
    1.67 +		aiMesh** out, unsigned int num, bool discard_input);
    1.68 +
    1.69 +	// ---------------------------------------------------------------------------
    1.70 +	/** Intermediate description of an edge between two corners of a polygon*/
    1.71 +	// ---------------------------------------------------------------------------
    1.72 +	struct Edge
    1.73 +	{
    1.74 +		Edge()
    1.75 +			: ref(0)
    1.76 +		{}
    1.77 +		Vertex edge_point, midpoint;
    1.78 +		unsigned int ref;
    1.79 +	};
    1.80 +
    1.81 +
    1.82 +
    1.83 +	typedef std::vector<unsigned int> UIntVector;
    1.84 +	typedef std::map<uint64_t,Edge> EdgeMap;
    1.85 +
    1.86 +	// ---------------------------------------------------------------------------
    1.87 +	// Hashing function to derive an index into an #EdgeMap from two given
    1.88 +	// 'unsigned int' vertex coordinates (!!distinct coordinates - same 
    1.89 +	// vertex position == same index!!). 
    1.90 +	// NOTE - this leads to rare hash collisions if a) sizeof(unsigned int)>4
    1.91 +	// and (id[0]>2^32-1 or id[0]>2^32-1).
    1.92 +	// MAKE_EDGE_HASH() uses temporaries, so INIT_EDGE_HASH() needs to be put
    1.93 +	// at the head of every function which is about to use MAKE_EDGE_HASH().
    1.94 +	// Reason is that the hash is that hash construction needs to hold the
    1.95 +	// invariant id0<id1 to identify an edge - else two hashes would refer
    1.96 +	// to the same edge.
    1.97 +	// ---------------------------------------------------------------------------
    1.98 +#define MAKE_EDGE_HASH(id0,id1) (eh_tmp0__=id0,eh_tmp1__=id1,\
    1.99 +	(eh_tmp0__<eh_tmp1__?std::swap(eh_tmp0__,eh_tmp1__):mydummy()),(uint64_t)eh_tmp0__^((uint64_t)eh_tmp1__<<32u))
   1.100 +
   1.101 +
   1.102 +#define INIT_EDGE_HASH_TEMPORARIES()\
   1.103 +	unsigned int eh_tmp0__, eh_tmp1__;
   1.104 +
   1.105 +private:
   1.106 +
   1.107 +	void InternSubdivide (const aiMesh* const * smesh, 
   1.108 +		size_t nmesh,aiMesh** out, unsigned int num);
   1.109 +};
   1.110 +
   1.111 +
   1.112 +// ------------------------------------------------------------------------------------------------
   1.113 +// Construct a subdivider of a specific type
   1.114 +Subdivider* Subdivider::Create (Algorithm algo)
   1.115 +{
   1.116 +	switch (algo)
   1.117 +	{
   1.118 +	case CATMULL_CLARKE:
   1.119 +		return new CatmullClarkSubdivider();
   1.120 +	};
   1.121 +
   1.122 +	ai_assert(false);
   1.123 +	return NULL; // shouldn't happen
   1.124 +}
   1.125 +
   1.126 +// ------------------------------------------------------------------------------------------------
   1.127 +// Call the Catmull Clark subdivision algorithm for one mesh
   1.128 +void  CatmullClarkSubdivider::Subdivide (
   1.129 +	aiMesh* mesh, 
   1.130 +	aiMesh*& out,
   1.131 +	unsigned int num,
   1.132 +	bool discard_input
   1.133 +	)
   1.134 +{
   1.135 +	assert(mesh != out);
   1.136 +	Subdivide(&mesh,1,&out,num,discard_input);
   1.137 +}
   1.138 +
   1.139 +// ------------------------------------------------------------------------------------------------
   1.140 +// Call the Catmull Clark subdivision algorithm for multiple meshes
   1.141 +void CatmullClarkSubdivider::Subdivide (
   1.142 +	aiMesh** smesh, 
   1.143 +	size_t nmesh,
   1.144 +	aiMesh** out, 
   1.145 +	unsigned int num,
   1.146 +	bool discard_input
   1.147 +	)
   1.148 +{
   1.149 +	ai_assert(NULL != smesh && NULL != out);
   1.150 +
   1.151 +	// course, both regions may not overlap
   1.152 +	assert(smesh<out || smesh+nmesh>out+nmesh);
   1.153 +	if (!num) {
   1.154 +		
   1.155 +		// No subdivision at all. Need to copy all the meshes .. argh.
   1.156 +		if (discard_input) {
   1.157 +			for (size_t s = 0; s < nmesh; ++s) {
   1.158 +				out[s] = smesh[s];
   1.159 +				smesh[s] = NULL;
   1.160 +			}
   1.161 +		}
   1.162 +		else {
   1.163 +			for (size_t s = 0; s < nmesh; ++s) {
   1.164 +				SceneCombiner::Copy(out+s,smesh[s]);
   1.165 +			}
   1.166 +		}
   1.167 +		return;
   1.168 +	}
   1.169 +
   1.170 +	std::vector<aiMesh*> inmeshes;
   1.171 +	std::vector<aiMesh*> outmeshes;
   1.172 +	std::vector<unsigned int> maptbl;
   1.173 +
   1.174 +	inmeshes.reserve(nmesh);
   1.175 +	outmeshes.reserve(nmesh);
   1.176 +	maptbl.reserve(nmesh);
   1.177 +
   1.178 +	// Remove pure line and point meshes from the working set to reduce the 
   1.179 +	// number of edge cases the subdivider is forced to deal with. Line and 
   1.180 +	// point meshes are simply passed through.
   1.181 +	for (size_t s = 0; s < nmesh; ++s) {
   1.182 +		aiMesh* i = smesh[s];
   1.183 +		// FIX - mPrimitiveTypes might not yet be initialized
   1.184 +		if (i->mPrimitiveTypes && (i->mPrimitiveTypes & (aiPrimitiveType_LINE|aiPrimitiveType_POINT))==i->mPrimitiveTypes) {
   1.185 +			DefaultLogger::get()->debug("Catmull-Clark Subdivider: Skipping pure line/point mesh");
   1.186 +
   1.187 +			if (discard_input) {
   1.188 +				out[s] = i;
   1.189 +				smesh[s] = NULL;
   1.190 +			}
   1.191 +			else {
   1.192 +				SceneCombiner::Copy(out+s,i);
   1.193 +			}
   1.194 +			continue;
   1.195 +		}
   1.196 +
   1.197 +		outmeshes.push_back(NULL);inmeshes.push_back(i);
   1.198 +		maptbl.push_back(s);
   1.199 +	}
   1.200 +
   1.201 +	// Do the actual subdivision on the preallocated storage. InternSubdivide 
   1.202 +	// *always* assumes that enough storage is available, it does not bother
   1.203 +	// checking any ranges.
   1.204 +	ai_assert(inmeshes.size()==outmeshes.size()&&inmeshes.size()==maptbl.size());
   1.205 +	if (inmeshes.empty()) {
   1.206 +		DefaultLogger::get()->warn("Catmull-Clark Subdivider: Pure point/line scene, I can't do anything");
   1.207 +		return;
   1.208 +	}
   1.209 +	InternSubdivide(&inmeshes.front(),inmeshes.size(),&outmeshes.front(),num);
   1.210 +	for (unsigned int i = 0; i < maptbl.size(); ++i) {
   1.211 +		ai_assert(outmeshes[i]);
   1.212 +		out[maptbl[i]] = outmeshes[i];
   1.213 +	}
   1.214 +
   1.215 +	if (discard_input) {
   1.216 +		for (size_t s = 0; s < nmesh; ++s) {
   1.217 +			delete smesh[s];
   1.218 +		}
   1.219 +	}
   1.220 +}
   1.221 +
   1.222 +// ------------------------------------------------------------------------------------------------
   1.223 +// Note - this is an implementation of the standard (recursive) Cm-Cl algorithm without further 
   1.224 +// optimizations (except we're using some nice LUTs). A description of the algorithm can be found
   1.225 +// here: http://en.wikipedia.org/wiki/Catmull-Clark_subdivision_surface
   1.226 +//
   1.227 +// The code is mostly O(n), however parts are O(nlogn) which is therefore the algorithm's
   1.228 +// expected total runtime complexity. The implementation is able to work in-place on the same
   1.229 +// mesh arrays. Calling #InternSubdivide() directly is not encouraged. The code can operate
   1.230 +// in-place unless 'smesh' and 'out' are equal (no strange overlaps or reorderings). 
   1.231 +// Previous data is replaced/deleted then.
   1.232 +// ------------------------------------------------------------------------------------------------
   1.233 +void CatmullClarkSubdivider::InternSubdivide (
   1.234 +	const aiMesh* const * smesh, 
   1.235 +	size_t nmesh,
   1.236 +	aiMesh** out, 
   1.237 +	unsigned int num
   1.238 +	)
   1.239 +{
   1.240 +	ai_assert(NULL != smesh && NULL != out);
   1.241 +	INIT_EDGE_HASH_TEMPORARIES();
   1.242 +
   1.243 +	// no subdivision requested or end of recursive refinement
   1.244 +	if (!num) {
   1.245 +		return;
   1.246 +	}
   1.247 +
   1.248 +	UIntVector maptbl;
   1.249 +	SpatialSort spatial;
   1.250 +
   1.251 +	// ---------------------------------------------------------------------
   1.252 +	// 0. Offset table to index all meshes continuously, generate a spatially
   1.253 +	// sorted representation of all vertices in all meshes.
   1.254 +	// ---------------------------------------------------------------------
   1.255 +	typedef std::pair<unsigned int,unsigned int> IntPair;
   1.256 +	std::vector<IntPair> moffsets(nmesh);
   1.257 +	unsigned int totfaces = 0, totvert = 0;
   1.258 +	for (size_t t = 0; t < nmesh; ++t) {
   1.259 +		const aiMesh* mesh = smesh[t];
   1.260 +
   1.261 +		spatial.Append(mesh->mVertices,mesh->mNumVertices,sizeof(aiVector3D),false);
   1.262 +		moffsets[t] = IntPair(totfaces,totvert);
   1.263 +
   1.264 +		totfaces += mesh->mNumFaces;
   1.265 +		totvert  += mesh->mNumVertices;
   1.266 +	}
   1.267 +
   1.268 +	spatial.Finalize();
   1.269 +	const unsigned int num_unique = spatial.GenerateMappingTable(maptbl,ComputePositionEpsilon(smesh,nmesh));
   1.270 +
   1.271 +
   1.272 +#define FLATTEN_VERTEX_IDX(mesh_idx, vert_idx) (moffsets[mesh_idx].second+vert_idx)
   1.273 +#define   FLATTEN_FACE_IDX(mesh_idx, face_idx) (moffsets[mesh_idx].first+face_idx)
   1.274 +
   1.275 +	// ---------------------------------------------------------------------
   1.276 +	// 1. Compute the centroid point for all faces
   1.277 +	// ---------------------------------------------------------------------
   1.278 +	std::vector<Vertex> centroids(totfaces);
   1.279 +	unsigned int nfacesout = 0;
   1.280 +	for (size_t t = 0, n = 0; t < nmesh; ++t) {
   1.281 +		const aiMesh* mesh = smesh[t];
   1.282 +		for (unsigned int i = 0; i < mesh->mNumFaces;++i,++n)
   1.283 +		{
   1.284 +			const aiFace& face = mesh->mFaces[i];
   1.285 +			Vertex& c = centroids[n];
   1.286 +
   1.287 +			for (unsigned int a = 0; a < face.mNumIndices;++a) {
   1.288 +				c += Vertex(mesh,face.mIndices[a]);
   1.289 +			}
   1.290 +
   1.291 +			c /= static_cast<float>(face.mNumIndices);
   1.292 +			nfacesout += face.mNumIndices;
   1.293 +		}
   1.294 +	}
   1.295 +	
   1.296 +	EdgeMap edges;
   1.297 +
   1.298 +	// ---------------------------------------------------------------------
   1.299 +	// 2. Set each edge point to be the average of all neighbouring 
   1.300 +	// face points and original points. Every edge exists twice
   1.301 +	// if there is a neighboring face.
   1.302 +	// ---------------------------------------------------------------------
   1.303 +	for (size_t t = 0; t < nmesh; ++t) {
   1.304 +		const aiMesh* mesh = smesh[t];
   1.305 +
   1.306 +		for (unsigned int i = 0; i < mesh->mNumFaces;++i)	{
   1.307 +			const aiFace& face = mesh->mFaces[i];
   1.308 +
   1.309 +			for (unsigned int p =0; p< face.mNumIndices; ++p) {
   1.310 +				const unsigned int id[] = { 
   1.311 +					face.mIndices[p], 
   1.312 +					face.mIndices[p==face.mNumIndices-1?0:p+1]
   1.313 +				};
   1.314 +				const unsigned int mp[] = { 
   1.315 +					maptbl[FLATTEN_VERTEX_IDX(t,id[0])], 
   1.316 +					maptbl[FLATTEN_VERTEX_IDX(t,id[1])]
   1.317 +				};
   1.318 +
   1.319 +				Edge& e = edges[MAKE_EDGE_HASH(mp[0],mp[1])];
   1.320 +				e.ref++;
   1.321 +				if (e.ref<=2) {
   1.322 +					if (e.ref==1) { // original points (end points) - add only once
   1.323 +						e.edge_point = e.midpoint = Vertex(mesh,id[0])+Vertex(mesh,id[1]);
   1.324 +						e.midpoint *= 0.5f;
   1.325 +					}
   1.326 +					e.edge_point += centroids[FLATTEN_FACE_IDX(t,i)];
   1.327 +				}
   1.328 +			}
   1.329 +		}
   1.330 +	}
   1.331 +
   1.332 +	// ---------------------------------------------------------------------
   1.333 +	// 3. Normalize edge points
   1.334 +	// ---------------------------------------------------------------------
   1.335 +	{unsigned int bad_cnt = 0;
   1.336 +	for (EdgeMap::iterator it = edges.begin(); it != edges.end(); ++it) {
   1.337 +		if ((*it).second.ref < 2) {
   1.338 +			ai_assert((*it).second.ref);
   1.339 +			++bad_cnt;
   1.340 +		}
   1.341 +		(*it).second.edge_point *= 1.f/((*it).second.ref+2.f);
   1.342 +	}
   1.343 +
   1.344 +	if (bad_cnt) {
   1.345 +		// Report the number of bad edges. bad edges are referenced by less than two
   1.346 +		// faces in the mesh. They occur at outer model boundaries in non-closed
   1.347 +		// shapes.
   1.348 +		char tmp[512];
   1.349 +		sprintf(tmp,"Catmull-Clark Subdivider: got %u bad edges touching only one face (totally %u edges). ",
   1.350 +			bad_cnt,static_cast<unsigned int>(edges.size()));
   1.351 +
   1.352 +		DefaultLogger::get()->debug(tmp);
   1.353 +	}}
   1.354 +
   1.355 +	// ---------------------------------------------------------------------
   1.356 +	// 4. Compute a vertex-face adjacency table. We can't reuse the code
   1.357 +	// from VertexTriangleAdjacency because we need the table for multiple
   1.358 +	// meshes and out vertex indices need to be mapped to distinct values 
   1.359 +	// first.
   1.360 +	// ---------------------------------------------------------------------
   1.361 +	UIntVector faceadjac(nfacesout), cntadjfac(maptbl.size(),0), ofsadjvec(maptbl.size()+1,0); {
   1.362 +	for (size_t t = 0; t < nmesh; ++t) {
   1.363 +		const aiMesh* const minp = smesh[t];
   1.364 +		for (unsigned int i = 0; i < minp->mNumFaces; ++i) {
   1.365 +			
   1.366 +			const aiFace& f = minp->mFaces[i];
   1.367 +			for (unsigned int n = 0; n < f.mNumIndices; ++n) {
   1.368 +				++cntadjfac[maptbl[FLATTEN_VERTEX_IDX(t,f.mIndices[n])]];
   1.369 +			}
   1.370 +		}
   1.371 +	}
   1.372 +	unsigned int cur = 0;
   1.373 +	for (size_t i = 0; i < cntadjfac.size(); ++i) {
   1.374 +		ofsadjvec[i+1] = cur;
   1.375 +		cur += cntadjfac[i];
   1.376 +	}
   1.377 +	for (size_t t = 0; t < nmesh; ++t) {
   1.378 +		const aiMesh* const minp = smesh[t];
   1.379 +		for (unsigned int i = 0; i < minp->mNumFaces; ++i) {
   1.380 +			
   1.381 +			const aiFace& f = minp->mFaces[i];
   1.382 +			for (unsigned int n = 0; n < f.mNumIndices; ++n) {
   1.383 +				faceadjac[ofsadjvec[1+maptbl[FLATTEN_VERTEX_IDX(t,f.mIndices[n])]]++] = FLATTEN_FACE_IDX(t,i);
   1.384 +			}
   1.385 +		}
   1.386 +	} 
   1.387 +	
   1.388 +	// check the other way round for consistency
   1.389 +#ifdef _DEBUG
   1.390 +
   1.391 +	for (size_t t = 0; t < ofsadjvec.size()-1; ++t) {
   1.392 +		for (unsigned int m = 0; m <  cntadjfac[t]; ++m) {
   1.393 +			const unsigned int fidx = faceadjac[ofsadjvec[t]+m];
   1.394 +			ai_assert(fidx < totfaces);
   1.395 +			for (size_t n = 1; n < nmesh; ++n) {
   1.396 +			
   1.397 +				if (moffsets[n].first > fidx) {
   1.398 +					const aiMesh* msh = smesh[--n];
   1.399 +					const aiFace& f = msh->mFaces[fidx-moffsets[n].first];
   1.400 +					
   1.401 +					bool haveit = false;
   1.402 +					for (unsigned int i = 0; i < f.mNumIndices; ++i) {
   1.403 +						if (maptbl[FLATTEN_VERTEX_IDX(n,f.mIndices[i])]==(unsigned int)t) {
   1.404 +							haveit = true; break;
   1.405 +						}
   1.406 +					}
   1.407 +					ai_assert(haveit);
   1.408 +					break;
   1.409 +				}
   1.410 +			}
   1.411 +		}
   1.412 +	}
   1.413 +
   1.414 +#endif
   1.415 +	}
   1.416 +
   1.417 +#define GET_ADJACENT_FACES_AND_CNT(vidx,fstartout,numout) \
   1.418 +	fstartout = &faceadjac[ofsadjvec[vidx]], numout = cntadjfac[vidx]
   1.419 +
   1.420 +	typedef std::pair<bool,Vertex> TouchedOVertex;
   1.421 +	std::vector<TouchedOVertex > new_points(num_unique,TouchedOVertex(false,Vertex()));
   1.422 +	// ---------------------------------------------------------------------
   1.423 +	// 5. Spawn a quad from each face point to the corresponding edge points
   1.424 +	// the original points being the fourth quad points.
   1.425 +	// ---------------------------------------------------------------------
   1.426 +	for (size_t t = 0; t < nmesh; ++t) {
   1.427 +		const aiMesh* const minp = smesh[t];
   1.428 +		aiMesh* const mout = out[t] = new aiMesh();
   1.429 +
   1.430 +		for (unsigned int a  = 0; a < minp->mNumFaces; ++a) {
   1.431 +			mout->mNumFaces += minp->mFaces[a].mNumIndices;
   1.432 +		}
   1.433 +
   1.434 +		// We need random access to the old face buffer, so reuse is not possible.
   1.435 +		mout->mFaces = new aiFace[mout->mNumFaces];
   1.436 +
   1.437 +		mout->mNumVertices = mout->mNumFaces*4;
   1.438 +		mout->mVertices = new aiVector3D[mout->mNumVertices];
   1.439 +
   1.440 +		// quads only, keep material index
   1.441 +		mout->mPrimitiveTypes = aiPrimitiveType_POLYGON;
   1.442 +		mout->mMaterialIndex = minp->mMaterialIndex;
   1.443 +
   1.444 +		if (minp->HasNormals()) {
   1.445 +			mout->mNormals = new aiVector3D[mout->mNumVertices];
   1.446 +		}
   1.447 +
   1.448 +		if (minp->HasTangentsAndBitangents()) {
   1.449 +			mout->mTangents = new aiVector3D[mout->mNumVertices];
   1.450 +			mout->mBitangents = new aiVector3D[mout->mNumVertices];
   1.451 +		}
   1.452 +
   1.453 +		for(unsigned int i = 0; minp->HasTextureCoords(i); ++i) {
   1.454 +			mout->mTextureCoords[i] = new aiVector3D[mout->mNumVertices];
   1.455 +			mout->mNumUVComponents[i] = minp->mNumUVComponents[i];
   1.456 +		}
   1.457 +
   1.458 +		for(unsigned int i = 0; minp->HasVertexColors(i); ++i) {
   1.459 +			mout->mColors[i] = new aiColor4D[mout->mNumVertices];
   1.460 +		}
   1.461 +
   1.462 +		mout->mNumVertices = mout->mNumFaces<<2u;
   1.463 +		for (unsigned int i = 0, v = 0, n = 0; i < minp->mNumFaces;++i)	{
   1.464 +
   1.465 +			const aiFace& face = minp->mFaces[i];
   1.466 +			for (unsigned int a = 0; a < face.mNumIndices;++a)	{
   1.467 +
   1.468 +				// Get a clean new face.
   1.469 +				aiFace& faceOut = mout->mFaces[n++];
   1.470 +				faceOut.mIndices = new unsigned int [faceOut.mNumIndices = 4];
   1.471 +
   1.472 +				// Spawn a new quadrilateral (ccw winding) for this original point between:
   1.473 +				// a) face centroid
   1.474 +				centroids[FLATTEN_FACE_IDX(t,i)].SortBack(mout,faceOut.mIndices[0]=v++); 
   1.475 +
   1.476 +				// b) adjacent edge on the left, seen from the centroid
   1.477 +				const Edge& e0 = edges[MAKE_EDGE_HASH(maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])],
   1.478 +					maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a==face.mNumIndices-1?0:a+1])
   1.479 +					])];  // fixme: replace with mod face.mNumIndices? 
   1.480 +
   1.481 +				// c) adjacent edge on the right, seen from the centroid
   1.482 +				const Edge& e1 = edges[MAKE_EDGE_HASH(maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])],
   1.483 +					maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[!a?face.mNumIndices-1:a-1])
   1.484 +					])];  // fixme: replace with mod face.mNumIndices? 
   1.485 +
   1.486 +				e0.edge_point.SortBack(mout,faceOut.mIndices[3]=v++);
   1.487 +				e1.edge_point.SortBack(mout,faceOut.mIndices[1]=v++);
   1.488 +
   1.489 +				// d= original point P with distinct index i
   1.490 +				// F := 0
   1.491 +				// R := 0
   1.492 +				// n := 0
   1.493 +				// for each face f containing i
   1.494 +				//    F := F+ centroid of f
   1.495 +				//    R := R+ midpoint of edge of f from i to i+1
   1.496 +				//    n := n+1
   1.497 +				//
   1.498 +				// (F+2R+(n-3)P)/n         
   1.499 +				const unsigned int org = maptbl[FLATTEN_VERTEX_IDX(t,face.mIndices[a])];
   1.500 +				TouchedOVertex& ov = new_points[org];
   1.501 +
   1.502 +				if (!ov.first) {
   1.503 +					ov.first = true;
   1.504 +
   1.505 +					const unsigned int* adj; unsigned int cnt;
   1.506 +					GET_ADJACENT_FACES_AND_CNT(org,adj,cnt);
   1.507 +
   1.508 +					if (cnt < 3) {
   1.509 +						ov.second = Vertex(minp,face.mIndices[a]);
   1.510 +					}
   1.511 +					else {
   1.512 +
   1.513 +						Vertex F,R;
   1.514 +						for (unsigned int o = 0; o < cnt; ++o) {
   1.515 +							ai_assert(adj[o] < totfaces);
   1.516 +							F += centroids[adj[o]];
   1.517 +
   1.518 +							// adj[0] is a global face index - search the face in the mesh list
   1.519 +							const aiMesh* mp = NULL;
   1.520 +							size_t nidx;
   1.521 +
   1.522 +							if (adj[o] < moffsets[0].first) {
   1.523 +								mp = smesh[nidx=0];
   1.524 +							}
   1.525 +							else {
   1.526 +								for (nidx = 1; nidx<= nmesh; ++nidx) {
   1.527 +									if (nidx == nmesh ||moffsets[nidx].first > adj[o]) {
   1.528 +										mp = smesh[--nidx];
   1.529 +										break;
   1.530 +									}
   1.531 +								}
   1.532 +							}
   1.533 +
   1.534 +							ai_assert(adj[o]-moffsets[nidx].first < mp->mNumFaces);
   1.535 +							const aiFace& f = mp->mFaces[adj[o]-moffsets[nidx].first];
   1.536 +#				ifdef _DEBUG
   1.537 +							bool haveit = false;
   1.538 +#				endif
   1.539 +
   1.540 +							// find our original point in the face
   1.541 +							for (unsigned int m = 0; m < f.mNumIndices; ++m) {
   1.542 +								if (maptbl[FLATTEN_VERTEX_IDX(nidx,f.mIndices[m])] == org) {
   1.543 +
   1.544 +									// add *both* edges. this way, we can be sure that we add
   1.545 +									// *all* adjacent edges to R. In a closed shape, every
   1.546 +									// edge is added twice - so we simply leave out the
   1.547 +									// factor 2.f in the amove formula and get the right
   1.548 +									// result.
   1.549 +
   1.550 +									const Edge& c0 = edges[MAKE_EDGE_HASH(org,maptbl[FLATTEN_VERTEX_IDX(
   1.551 +										nidx,f.mIndices[!m?f.mNumIndices-1:m-1])])];
   1.552 +									// fixme: replace with mod face.mNumIndices? 
   1.553 +
   1.554 +									const Edge& c1 = edges[MAKE_EDGE_HASH(org,maptbl[FLATTEN_VERTEX_IDX(
   1.555 +										nidx,f.mIndices[m==f.mNumIndices-1?0:m+1])])];
   1.556 +									// fixme: replace with mod face.mNumIndices? 
   1.557 +									R += c0.midpoint+c1.midpoint;
   1.558 +
   1.559 +#						ifdef _DEBUG
   1.560 +									haveit = true;
   1.561 +#						endif
   1.562 +									break;
   1.563 +								}
   1.564 +							}
   1.565 +
   1.566 +							// this invariant *must* hold if the vertex-to-face adjacency table is valid
   1.567 +							ai_assert(haveit);
   1.568 +						}
   1.569 +
   1.570 +						const float div = static_cast<float>(cnt), divsq = 1.f/(div*div);
   1.571 +						ov.second = Vertex(minp,face.mIndices[a])*((div-3.f) / div) + R*divsq + F*divsq;  
   1.572 +					}
   1.573 +				}
   1.574 +				ov.second.SortBack(mout,faceOut.mIndices[2]=v++);
   1.575 +			}
   1.576 +		}
   1.577 +	}
   1.578 +
   1.579 +	// ---------------------------------------------------------------------
   1.580 +	// 7. Apply the next subdivision step. 
   1.581 +	// ---------------------------------------------------------------------
   1.582 +	if (num != 1) {
   1.583 +		std::vector<aiMesh*> tmp(nmesh);
   1.584 +		InternSubdivide (out,nmesh,&tmp.front(),num-1);
   1.585 +		for (size_t i = 0; i < nmesh; ++i) {
   1.586 +			delete out[i];
   1.587 +			out[i] = tmp[i];
   1.588 +		}
   1.589 +	}
   1.590 +}