vrshoot
view libs/assimp/SpatialSort.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 ---------------------------------------------------------------------------
3 Open Asset Import Library (assimp)
4 ---------------------------------------------------------------------------
6 Copyright (c) 2006-2012, assimp team
8 All rights reserved.
10 Redistribution and use of this software in source and binary forms,
11 with or without modification, are permitted provided that the following
12 conditions are met:
14 * Redistributions of source code must retain the above
15 copyright notice, this list of conditions and the
16 following disclaimer.
18 * Redistributions in binary form must reproduce the above
19 copyright notice, this list of conditions and the
20 following disclaimer in the documentation and/or other
21 materials provided with the distribution.
23 * Neither the name of the assimp team, nor the names of its
24 contributors may be used to endorse or promote products
25 derived from this software without specific prior
26 written permission of the assimp team.
28 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
29 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
30 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
31 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
32 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
33 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
34 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
35 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
36 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
37 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
38 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
39 ---------------------------------------------------------------------------
40 */
42 /** @file Implementation of the helper class to quickly find vertices close to a given position */
44 #include "AssimpPCH.h"
45 #include "SpatialSort.h"
47 using namespace Assimp;
49 // CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8.
50 #ifndef CHAR_BIT
51 # define CHAR_BIT 8
52 #endif
54 // ------------------------------------------------------------------------------------------------
55 // Constructs a spatially sorted representation from the given position array.
56 SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions,
57 unsigned int pElementOffset)
59 // define the reference plane. We choose some arbitrary vector away from all basic axises
60 // in the hope that no model spreads all its vertices along this plane.
61 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
62 {
63 mPlaneNormal.Normalize();
64 Fill(pPositions,pNumPositions,pElementOffset);
65 }
67 // ------------------------------------------------------------------------------------------------
68 SpatialSort :: SpatialSort()
69 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
70 {
71 mPlaneNormal.Normalize();
72 }
74 // ------------------------------------------------------------------------------------------------
75 // Destructor
76 SpatialSort::~SpatialSort()
77 {
78 // nothing to do here, everything destructs automatically
79 }
81 // ------------------------------------------------------------------------------------------------
82 void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions,
83 unsigned int pElementOffset,
84 bool pFinalize /*= true */)
85 {
86 mPositions.clear();
87 Append(pPositions,pNumPositions,pElementOffset,pFinalize);
88 }
90 // ------------------------------------------------------------------------------------------------
91 void SpatialSort :: Finalize()
92 {
93 std::sort( mPositions.begin(), mPositions.end());
94 }
96 // ------------------------------------------------------------------------------------------------
97 void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions,
98 unsigned int pElementOffset,
99 bool pFinalize /*= true */)
100 {
101 // store references to all given positions along with their distance to the reference plane
102 const size_t initial = mPositions.size();
103 mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2));
104 for( unsigned int a = 0; a < pNumPositions; a++)
105 {
106 const char* tempPointer = reinterpret_cast<const char*> (pPositions);
107 const aiVector3D* vec = reinterpret_cast<const aiVector3D*> (tempPointer + a * pElementOffset);
109 // store position by index and distance
110 float distance = *vec * mPlaneNormal;
111 mPositions.push_back( Entry( a+initial, *vec, distance));
112 }
114 if (pFinalize) {
115 // now sort the array ascending by distance.
116 Finalize();
117 }
118 }
120 // ------------------------------------------------------------------------------------------------
121 // Returns an iterator for all positions close to the given position.
122 void SpatialSort::FindPositions( const aiVector3D& pPosition,
123 float pRadius, std::vector<unsigned int>& poResults) const
124 {
125 const float dist = pPosition * mPlaneNormal;
126 const float minDist = dist - pRadius, maxDist = dist + pRadius;
128 // clear the array in this strange fashion because a simple clear() would also deallocate
129 // the array which we want to avoid
130 poResults.erase( poResults.begin(), poResults.end());
132 // quick check for positions outside the range
133 if( mPositions.size() == 0)
134 return;
135 if( maxDist < mPositions.front().mDistance)
136 return;
137 if( minDist > mPositions.back().mDistance)
138 return;
140 // do a binary search for the minimal distance to start the iteration there
141 unsigned int index = (unsigned int)mPositions.size() / 2;
142 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
143 while( binaryStepSize > 1)
144 {
145 if( mPositions[index].mDistance < minDist)
146 index += binaryStepSize;
147 else
148 index -= binaryStepSize;
150 binaryStepSize /= 2;
151 }
153 // depending on the direction of the last step we need to single step a bit back or forth
154 // to find the actual beginning element of the range
155 while( index > 0 && mPositions[index].mDistance > minDist)
156 index--;
157 while( index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist)
158 index++;
160 // Mow start iterating from there until the first position lays outside of the distance range.
161 // Add all positions inside the distance range within the given radius to the result aray
162 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
163 const float pSquared = pRadius*pRadius;
164 while( it->mDistance < maxDist)
165 {
166 if( (it->mPosition - pPosition).SquareLength() < pSquared)
167 poResults.push_back( it->mIndex);
168 ++it;
169 if( it == mPositions.end())
170 break;
171 }
173 // that's it
174 }
176 namespace {
178 // Binary, signed-integer representation of a single-precision floating-point value.
179 // IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are
180 // ordered the same way when their bits are reinterpreted as sign-magnitude integers."
181 // This allows us to convert all floating-point numbers to signed integers of arbitrary size
182 // and then use them to work with ULPs (Units in the Last Place, for high-precision
183 // computations) or to compare them (integer comparisons are faster than floating-point
184 // comparisons on many platforms).
185 typedef signed int BinFloat;
187 // --------------------------------------------------------------------------------------------
188 // Converts the bit pattern of a floating-point number to its signed integer representation.
189 BinFloat ToBinary( const float & pValue) {
191 // If this assertion fails, signed int is not big enough to store a float on your platform.
192 // Please correct the declaration of BinFloat a few lines above - but do it in a portable,
193 // #ifdef'd manner!
194 BOOST_STATIC_ASSERT( sizeof(BinFloat) >= sizeof(float));
196 #if defined( _MSC_VER)
197 // If this assertion fails, Visual C++ has finally moved to ILP64. This means that this
198 // code has just become legacy code! Find out the current value of _MSC_VER and modify
199 // the #if above so it evaluates false on the current and all upcoming VC versions (or
200 // on the current platform, if LP64 or LLP64 are still used on other platforms).
201 BOOST_STATIC_ASSERT( sizeof(BinFloat) == sizeof(float));
203 // This works best on Visual C++, but other compilers have their problems with it.
204 const BinFloat binValue = reinterpret_cast<BinFloat const &>(pValue);
205 #else
206 // On many compilers, reinterpreting a float address as an integer causes aliasing
207 // problems. This is an ugly but more or less safe way of doing it.
208 union {
209 float asFloat;
210 BinFloat asBin;
211 } conversion;
212 conversion.asBin = 0; // zero empty space in case sizeof(BinFloat) > sizeof(float)
213 conversion.asFloat = pValue;
214 const BinFloat binValue = conversion.asBin;
215 #endif
217 // floating-point numbers are of sign-magnitude format, so find out what signed number
218 // representation we must convert negative values to.
219 // See http://en.wikipedia.org/wiki/Signed_number_representations.
221 // Two's complement?
222 if( (-42 == (~42 + 1)) && (binValue & 0x80000000))
223 return BinFloat(1 << (CHAR_BIT * sizeof(BinFloat) - 1)) - binValue;
224 // One's complement?
225 else if( (-42 == ~42) && (binValue & 0x80000000))
226 return BinFloat(-0) - binValue;
227 // Sign-magnitude?
228 else if( (-42 == (42 | (-0))) && (binValue & 0x80000000)) // -0 = 1000... binary
229 return binValue;
230 else
231 return binValue;
232 }
234 } // namespace
236 // ------------------------------------------------------------------------------------------------
237 // Fills an array with indices of all positions indentical to the given position. In opposite to
238 // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units.
239 void SpatialSort::FindIdenticalPositions( const aiVector3D& pPosition,
240 std::vector<unsigned int>& poResults) const
241 {
242 // Epsilons have a huge disadvantage: they are of constant precision, while floating-point
243 // values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but
244 // if you apply it to 0.001, it is enormous.
246 // The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs
247 // tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of
248 // logarithmic precision - around 1, they are 1÷(2^24) and around 10000, they are 0.00125.
250 // For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The
251 // incoming vertex positions might have already been transformed, probably using rather
252 // inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify
253 // identical vertex positions.
254 static const int toleranceInULPs = 4;
255 // An interesting point is that the inaccuracy grows linear with the number of operations:
256 // multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs
257 // plus 0.5 ULPs for the multiplication.
258 // To compute the distance to the plane, a dot product is needed - that is a multiplication and
259 // an addition on each number.
260 static const int distanceToleranceInULPs = toleranceInULPs + 1;
261 // The squared distance between two 3D vectors is computed the same way, but with an additional
262 // subtraction.
263 static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1;
265 // Convert the plane distance to its signed integer representation so the ULPs tolerance can be
266 // applied. For some reason, VC won't optimize two calls of the bit pattern conversion.
267 const BinFloat minDistBinary = ToBinary( pPosition * mPlaneNormal) - distanceToleranceInULPs;
268 const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs;
270 // clear the array in this strange fashion because a simple clear() would also deallocate
271 // the array which we want to avoid
272 poResults.erase( poResults.begin(), poResults.end());
274 // do a binary search for the minimal distance to start the iteration there
275 unsigned int index = (unsigned int)mPositions.size() / 2;
276 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
277 while( binaryStepSize > 1)
278 {
279 // Ugly, but conditional jumps are faster with integers than with floats
280 if( minDistBinary > ToBinary(mPositions[index].mDistance))
281 index += binaryStepSize;
282 else
283 index -= binaryStepSize;
285 binaryStepSize /= 2;
286 }
288 // depending on the direction of the last step we need to single step a bit back or forth
289 // to find the actual beginning element of the range
290 while( index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance) )
291 index--;
292 while( index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance))
293 index++;
295 // Now start iterating from there until the first position lays outside of the distance range.
296 // Add all positions inside the distance range within the tolerance to the result aray
297 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
298 while( ToBinary(it->mDistance) < maxDistBinary)
299 {
300 if( distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength()))
301 poResults.push_back(it->mIndex);
302 ++it;
303 if( it == mPositions.end())
304 break;
305 }
307 // that's it
308 }
310 // ------------------------------------------------------------------------------------------------
311 unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int>& fill,float pRadius) const
312 {
313 fill.resize(mPositions.size(),UINT_MAX);
314 float dist, maxDist;
316 unsigned int t=0;
317 const float pSquared = pRadius*pRadius;
318 for (size_t i = 0; i < mPositions.size();) {
319 dist = mPositions[i].mPosition * mPlaneNormal;
320 maxDist = dist + pRadius;
322 fill[mPositions[i].mIndex] = t;
323 const aiVector3D& oldpos = mPositions[i].mPosition;
324 for (++i; i < fill.size() && mPositions[i].mDistance < maxDist
325 && (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i)
326 {
327 fill[mPositions[i].mIndex] = t;
328 }
329 ++t;
330 }
332 #ifdef _DEBUG
334 // debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1
335 for (size_t i = 0; i < fill.size(); ++i) {
336 ai_assert(fill[i]<mPositions.size());
337 }
339 #endif
340 return t;
341 }