1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/libs/assimp/SpatialSort.cpp	Sat Feb 01 19:58:19 2014 +0200
     1.3 @@ -0,0 +1,342 @@
     1.4 +/*
     1.5 +---------------------------------------------------------------------------
     1.6 +Open Asset Import Library (assimp)
     1.7 +---------------------------------------------------------------------------
     1.8 +
     1.9 +Copyright (c) 2006-2012, assimp team
    1.10 +
    1.11 +All rights reserved.
    1.12 +
    1.13 +Redistribution and use of this software in source and binary forms, 
    1.14 +with or without modification, are permitted provided that the following 
    1.15 +conditions are met:
    1.16 +
    1.17 +* Redistributions of source code must retain the above
    1.18 +  copyright notice, this list of conditions and the
    1.19 +  following disclaimer.
    1.20 +
    1.21 +* Redistributions in binary form must reproduce the above
    1.22 +  copyright notice, this list of conditions and the
    1.23 +  following disclaimer in the documentation and/or other
    1.24 +  materials provided with the distribution.
    1.25 +
    1.26 +* Neither the name of the assimp team, nor the names of its
    1.27 +  contributors may be used to endorse or promote products
    1.28 +  derived from this software without specific prior
    1.29 +  written permission of the assimp team.
    1.30 +
    1.31 +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
    1.32 +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
    1.33 +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
    1.34 +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 
    1.35 +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
    1.36 +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 
    1.37 +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
    1.38 +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
    1.39 +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
    1.40 +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 
    1.41 +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
    1.42 +---------------------------------------------------------------------------
    1.43 +*/
    1.44 +
    1.45 +/** @file Implementation of the helper class to quickly find vertices close to a given position */
    1.46 +
    1.47 +#include "AssimpPCH.h"
    1.48 +#include "SpatialSort.h"
    1.49 +
    1.50 +using namespace Assimp;
    1.51 +
    1.52 +// CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8.
    1.53 +#ifndef CHAR_BIT
    1.54 +#	define CHAR_BIT 8
    1.55 +#endif
    1.56 +
    1.57 +// ------------------------------------------------------------------------------------------------
    1.58 +// Constructs a spatially sorted representation from the given position array.
    1.59 +SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions, 
    1.60 +	unsigned int pElementOffset)
    1.61 +
    1.62 +	// define the reference plane. We choose some arbitrary vector away from all basic axises 
    1.63 +	// in the hope that no model spreads all its vertices along this plane.
    1.64 +	: mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
    1.65 +{
    1.66 +	mPlaneNormal.Normalize();
    1.67 +	Fill(pPositions,pNumPositions,pElementOffset);
    1.68 +}
    1.69 +
    1.70 +// ------------------------------------------------------------------------------------------------
    1.71 +SpatialSort :: SpatialSort()
    1.72 +: mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
    1.73 +{
    1.74 +	mPlaneNormal.Normalize();
    1.75 +}
    1.76 +
    1.77 +// ------------------------------------------------------------------------------------------------
    1.78 +// Destructor
    1.79 +SpatialSort::~SpatialSort()
    1.80 +{
    1.81 +	// nothing to do here, everything destructs automatically
    1.82 +}
    1.83 +
    1.84 +// ------------------------------------------------------------------------------------------------
    1.85 +void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions, 
    1.86 +	unsigned int pElementOffset,
    1.87 +	bool pFinalize /*= true */)
    1.88 +{
    1.89 +	mPositions.clear();
    1.90 +	Append(pPositions,pNumPositions,pElementOffset,pFinalize);
    1.91 +}
    1.92 +
    1.93 +// ------------------------------------------------------------------------------------------------
    1.94 +void SpatialSort :: Finalize()
    1.95 +{
    1.96 +	std::sort( mPositions.begin(), mPositions.end());
    1.97 +}
    1.98 +
    1.99 +// ------------------------------------------------------------------------------------------------
   1.100 +void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions, 
   1.101 +	unsigned int pElementOffset,
   1.102 +	bool pFinalize /*= true */)
   1.103 +{
   1.104 +	// store references to all given positions along with their distance to the reference plane
   1.105 +	const size_t initial = mPositions.size();
   1.106 +	mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2));
   1.107 +	for( unsigned int a = 0; a < pNumPositions; a++)
   1.108 +	{
   1.109 +		const char* tempPointer = reinterpret_cast<const char*> (pPositions);
   1.110 +		const aiVector3D* vec   = reinterpret_cast<const aiVector3D*> (tempPointer + a * pElementOffset);
   1.111 +
   1.112 +		// store position by index and distance
   1.113 +		float distance = *vec * mPlaneNormal;
   1.114 +		mPositions.push_back( Entry( a+initial, *vec, distance));
   1.115 +	}
   1.116 +
   1.117 +	if (pFinalize) {
   1.118 +		// now sort the array ascending by distance.
   1.119 +		Finalize();
   1.120 +	}
   1.121 +}
   1.122 +
   1.123 +// ------------------------------------------------------------------------------------------------
   1.124 +// Returns an iterator for all positions close to the given position.
   1.125 +void SpatialSort::FindPositions( const aiVector3D& pPosition, 
   1.126 +	float pRadius, std::vector<unsigned int>& poResults) const
   1.127 +{
   1.128 +	const float dist = pPosition * mPlaneNormal;
   1.129 +	const float minDist = dist - pRadius, maxDist = dist + pRadius;
   1.130 +
   1.131 +	// clear the array in this strange fashion because a simple clear() would also deallocate
   1.132 +    // the array which we want to avoid
   1.133 +	poResults.erase( poResults.begin(), poResults.end());
   1.134 +
   1.135 +	// quick check for positions outside the range
   1.136 +	if( mPositions.size() == 0)
   1.137 +		return;
   1.138 +	if( maxDist < mPositions.front().mDistance)
   1.139 +		return;
   1.140 +	if( minDist > mPositions.back().mDistance)
   1.141 +		return;
   1.142 +
   1.143 +	// do a binary search for the minimal distance to start the iteration there
   1.144 +	unsigned int index = (unsigned int)mPositions.size() / 2;
   1.145 +	unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
   1.146 +	while( binaryStepSize > 1)
   1.147 +	{
   1.148 +		if( mPositions[index].mDistance < minDist)
   1.149 +			index += binaryStepSize;
   1.150 +		else
   1.151 +			index -= binaryStepSize;
   1.152 +
   1.153 +		binaryStepSize /= 2;
   1.154 +	}
   1.155 +
   1.156 +	// depending on the direction of the last step we need to single step a bit back or forth
   1.157 +	// to find the actual beginning element of the range
   1.158 +	while( index > 0 && mPositions[index].mDistance > minDist)
   1.159 +		index--;
   1.160 +	while( index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist)
   1.161 +		index++;
   1.162 +	
   1.163 +	// Mow start iterating from there until the first position lays outside of the distance range.
   1.164 +	// Add all positions inside the distance range within the given radius to the result aray
   1.165 +	std::vector<Entry>::const_iterator it = mPositions.begin() + index;
   1.166 +	const float pSquared = pRadius*pRadius;
   1.167 +	while( it->mDistance < maxDist)
   1.168 +	{
   1.169 +		if( (it->mPosition - pPosition).SquareLength() < pSquared)
   1.170 +			poResults.push_back( it->mIndex);
   1.171 +		++it;
   1.172 +		if( it == mPositions.end())
   1.173 +			break;
   1.174 +	}
   1.175 +
   1.176 +	// that's it
   1.177 +}
   1.178 +
   1.179 +namespace {
   1.180 +
   1.181 +	// Binary, signed-integer representation of a single-precision floating-point value.
   1.182 +	// IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are
   1.183 +	//	ordered the same way when their bits are reinterpreted as sign-magnitude integers."
   1.184 +	// This allows us to convert all floating-point numbers to signed integers of arbitrary size
   1.185 +	//	and then use them to work with ULPs (Units in the Last Place, for high-precision
   1.186 +	//	computations) or to compare them (integer comparisons are faster than floating-point
   1.187 +	//	comparisons on many platforms).
   1.188 +	typedef signed int BinFloat;
   1.189 +
   1.190 +	// --------------------------------------------------------------------------------------------
   1.191 +	// Converts the bit pattern of a floating-point number to its signed integer representation.
   1.192 +	BinFloat ToBinary( const float & pValue) {
   1.193 +
   1.194 +		// If this assertion fails, signed int is not big enough to store a float on your platform.
   1.195 +		//	Please correct the declaration of BinFloat a few lines above - but do it in a portable,
   1.196 +		//	#ifdef'd manner!
   1.197 +		BOOST_STATIC_ASSERT( sizeof(BinFloat) >= sizeof(float));
   1.198 +
   1.199 +		#if defined( _MSC_VER)
   1.200 +			// If this assertion fails, Visual C++ has finally moved to ILP64. This means that this
   1.201 +			//	code has just become legacy code! Find out the current value of _MSC_VER and modify
   1.202 +			//	the #if above so it evaluates false on the current and all upcoming VC versions (or
   1.203 +			//	on the current platform, if LP64 or LLP64 are still used on other platforms).
   1.204 +			BOOST_STATIC_ASSERT( sizeof(BinFloat) == sizeof(float));
   1.205 +
   1.206 +			// This works best on Visual C++, but other compilers have their problems with it.
   1.207 +			const BinFloat binValue = reinterpret_cast<BinFloat const &>(pValue);
   1.208 +		#else
   1.209 +			// On many compilers, reinterpreting a float address as an integer causes aliasing
   1.210 +			// problems. This is an ugly but more or less safe way of doing it.
   1.211 +			union {
   1.212 +				float		asFloat;
   1.213 +				BinFloat	asBin;
   1.214 +			} conversion;
   1.215 +			conversion.asBin	= 0; // zero empty space in case sizeof(BinFloat) > sizeof(float)
   1.216 +			conversion.asFloat	= pValue;
   1.217 +			const BinFloat binValue = conversion.asBin;
   1.218 +		#endif
   1.219 +
   1.220 +		// floating-point numbers are of sign-magnitude format, so find out what signed number
   1.221 +		//	representation we must convert negative values to.
   1.222 +		// See http://en.wikipedia.org/wiki/Signed_number_representations.
   1.223 +
   1.224 +		// Two's complement?
   1.225 +		if( (-42 == (~42 + 1)) && (binValue & 0x80000000))
   1.226 +			return BinFloat(1 << (CHAR_BIT * sizeof(BinFloat) - 1)) - binValue;
   1.227 +		// One's complement?
   1.228 +		else if( (-42 == ~42) && (binValue & 0x80000000))
   1.229 +			return BinFloat(-0) - binValue;
   1.230 +		// Sign-magnitude?
   1.231 +		else if( (-42 == (42 | (-0))) && (binValue & 0x80000000)) // -0 = 1000... binary
   1.232 +			return binValue;
   1.233 +		else
   1.234 +			return binValue;
   1.235 +	}
   1.236 +
   1.237 +} // namespace
   1.238 +
   1.239 +// ------------------------------------------------------------------------------------------------
   1.240 +// Fills an array with indices of all positions indentical to the given position. In opposite to
   1.241 +// FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units.
   1.242 +void SpatialSort::FindIdenticalPositions( const aiVector3D& pPosition, 
   1.243 +	std::vector<unsigned int>& poResults) const
   1.244 +{
   1.245 +	// Epsilons have a huge disadvantage: they are of constant precision, while floating-point
   1.246 +	//	values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but
   1.247 +	//	if you apply it to 0.001, it is enormous.
   1.248 +
   1.249 +	// The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs
   1.250 +	//	tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of
   1.251 +	//	logarithmic precision - around 1, they are 1÷(2^24) and around 10000, they are 0.00125.
   1.252 +
   1.253 +	// For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The
   1.254 +	//	incoming vertex positions might have already been transformed, probably using rather
   1.255 +	//	inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify
   1.256 +	//	identical vertex positions.
   1.257 +	static const int toleranceInULPs = 4;
   1.258 +	// An interesting point is that the inaccuracy grows linear with the number of operations:
   1.259 +	//	multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs
   1.260 +	//	plus 0.5 ULPs for the multiplication.
   1.261 +	// To compute the distance to the plane, a dot product is needed - that is a multiplication and
   1.262 +	//	an addition on each number.
   1.263 +	static const int distanceToleranceInULPs = toleranceInULPs + 1;
   1.264 +	// The squared distance between two 3D vectors is computed the same way, but with an additional
   1.265 +	//	subtraction.
   1.266 +	static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1;
   1.267 +
   1.268 +	// Convert the plane distance to its signed integer representation so the ULPs tolerance can be
   1.269 +	//	applied. For some reason, VC won't optimize two calls of the bit pattern conversion.
   1.270 +	const BinFloat minDistBinary = ToBinary( pPosition * mPlaneNormal) - distanceToleranceInULPs;
   1.271 +	const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs;
   1.272 +
   1.273 +	// clear the array in this strange fashion because a simple clear() would also deallocate
   1.274 +    // the array which we want to avoid
   1.275 +	poResults.erase( poResults.begin(), poResults.end());
   1.276 +
   1.277 +	// do a binary search for the minimal distance to start the iteration there
   1.278 +	unsigned int index = (unsigned int)mPositions.size() / 2;
   1.279 +	unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
   1.280 +	while( binaryStepSize > 1)
   1.281 +	{
   1.282 +		// Ugly, but conditional jumps are faster with integers than with floats
   1.283 +		if( minDistBinary > ToBinary(mPositions[index].mDistance))
   1.284 +			index += binaryStepSize;
   1.285 +		else
   1.286 +			index -= binaryStepSize;
   1.287 +
   1.288 +		binaryStepSize /= 2;
   1.289 +	}
   1.290 +
   1.291 +	// depending on the direction of the last step we need to single step a bit back or forth
   1.292 +	// to find the actual beginning element of the range
   1.293 +	while( index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance) )
   1.294 +		index--;
   1.295 +	while( index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance))
   1.296 +		index++;
   1.297 +
   1.298 +	// Now start iterating from there until the first position lays outside of the distance range.
   1.299 +	// Add all positions inside the distance range within the tolerance to the result aray
   1.300 +	std::vector<Entry>::const_iterator it = mPositions.begin() + index;
   1.301 +	while( ToBinary(it->mDistance) < maxDistBinary)
   1.302 +	{
   1.303 +		if( distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength()))
   1.304 +			poResults.push_back(it->mIndex);
   1.305 +		++it;
   1.306 +		if( it == mPositions.end())
   1.307 +			break;
   1.308 +	}
   1.309 +
   1.310 +	// that's it
   1.311 +}
   1.312 +
   1.313 +// ------------------------------------------------------------------------------------------------
   1.314 +unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int>& fill,float pRadius) const
   1.315 +{
   1.316 +	fill.resize(mPositions.size(),UINT_MAX);
   1.317 +	float dist, maxDist;
   1.318 +
   1.319 +	unsigned int t=0;
   1.320 +	const float pSquared = pRadius*pRadius;
   1.321 +	for (size_t i = 0; i < mPositions.size();) {
   1.322 +		dist = mPositions[i].mPosition * mPlaneNormal;
   1.323 +		maxDist = dist + pRadius;
   1.324 +
   1.325 +		fill[mPositions[i].mIndex] = t;
   1.326 +		const aiVector3D& oldpos = mPositions[i].mPosition;
   1.327 +		for (++i; i < fill.size() && mPositions[i].mDistance < maxDist 
   1.328 +			&& (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i) 
   1.329 +		{
   1.330 +			fill[mPositions[i].mIndex] = t;
   1.331 +		}
   1.332 +		++t;
   1.333 +	}
   1.334 +
   1.335 +#ifdef _DEBUG
   1.336 +
   1.337 +	// debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1
   1.338 +	for (size_t i = 0; i < fill.size(); ++i) {
   1.339 +		ai_assert(fill[i]<mPositions.size());
   1.340 +	}
   1.341 +
   1.342 +#endif
   1.343 +	return t;
   1.344 +}
   1.345 +