1 /*
2 ---------------------------------------------------------------------------
3 Open Asset Import Library (assimp)
4 ---------------------------------------------------------------------------
5
6 Copyright (c) 2006-2017, assimp team
7
8
9 All rights reserved.
10
11 Redistribution and use of this software in source and binary forms,
12 with or without modification, are permitted provided that the following
13 conditions are met:
14
15 * Redistributions of source code must retain the above
16 copyright notice, this list of conditions and the
17 following disclaimer.
18
19 * Redistributions in binary form must reproduce the above
20 copyright notice, this list of conditions and the
21 following disclaimer in the documentation and/or other
22 materials provided with the distribution.
23
24 * Neither the name of the assimp team, nor the names of its
25 contributors may be used to endorse or promote products
26 derived from this software without specific prior
27 written permission of the assimp team.
28
29 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
30 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
31 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
32 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
33 OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
34 SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
35 LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
36 DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
37 THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
38 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
39 OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
40 ---------------------------------------------------------------------------
41 */
42
43 /** @file Implementation of the helper class to quickly find vertices close to a given position */
44
45 #include "SpatialSort.h"
46 #include <assimp/ai_assert.h>
47
48 using namespace Assimp;
49
50 // CHAR_BIT seems to be defined under MVSC, but not under GCC. Pray that the correct value is 8.
51 #ifndef CHAR_BIT
52 # define CHAR_BIT 8
53 #endif
54
55 // ------------------------------------------------------------------------------------------------
56 // Constructs a spatially sorted representation from the given position array.
SpatialSort(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset)57 SpatialSort::SpatialSort( const aiVector3D* pPositions, unsigned int pNumPositions,
58 unsigned int pElementOffset)
59
60 // define the reference plane. We choose some arbitrary vector away from all basic axises
61 // in the hope that no model spreads all its vertices along this plane.
62 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
63 {
64 mPlaneNormal.Normalize();
65 Fill(pPositions,pNumPositions,pElementOffset);
66 }
67
68 // ------------------------------------------------------------------------------------------------
SpatialSort()69 SpatialSort :: SpatialSort()
70 : mPlaneNormal(0.8523f, 0.34321f, 0.5736f)
71 {
72 mPlaneNormal.Normalize();
73 }
74
75 // ------------------------------------------------------------------------------------------------
76 // Destructor
~SpatialSort()77 SpatialSort::~SpatialSort()
78 {
79 // nothing to do here, everything destructs automatically
80 }
81
82 // ------------------------------------------------------------------------------------------------
Fill(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset,bool pFinalize)83 void SpatialSort::Fill( const aiVector3D* pPositions, unsigned int pNumPositions,
84 unsigned int pElementOffset,
85 bool pFinalize /*= true */)
86 {
87 mPositions.clear();
88 Append(pPositions,pNumPositions,pElementOffset,pFinalize);
89 }
90
91 // ------------------------------------------------------------------------------------------------
Finalize()92 void SpatialSort :: Finalize()
93 {
94 std::sort( mPositions.begin(), mPositions.end());
95 }
96
97 // ------------------------------------------------------------------------------------------------
Append(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset,bool pFinalize)98 void SpatialSort::Append( const aiVector3D* pPositions, unsigned int pNumPositions,
99 unsigned int pElementOffset,
100 bool pFinalize /*= true */)
101 {
102 // store references to all given positions along with their distance to the reference plane
103 const size_t initial = mPositions.size();
104 mPositions.reserve(initial + (pFinalize?pNumPositions:pNumPositions*2));
105 for( unsigned int a = 0; a < pNumPositions; a++)
106 {
107 const char* tempPointer = reinterpret_cast<const char*> (pPositions);
108 const aiVector3D* vec = reinterpret_cast<const aiVector3D*> (tempPointer + a * pElementOffset);
109
110 // store position by index and distance
111 ai_real distance = *vec * mPlaneNormal;
112 mPositions.push_back( Entry( static_cast<unsigned int>(a+initial), *vec, distance));
113 }
114
115 if (pFinalize) {
116 // now sort the array ascending by distance.
117 Finalize();
118 }
119 }
120
121 // ------------------------------------------------------------------------------------------------
122 // Returns an iterator for all positions close to the given position.
FindPositions(const aiVector3D & pPosition,ai_real pRadius,std::vector<unsigned int> & poResults) const123 void SpatialSort::FindPositions( const aiVector3D& pPosition,
124 ai_real pRadius, std::vector<unsigned int>& poResults) const
125 {
126 const ai_real dist = pPosition * mPlaneNormal;
127 const ai_real minDist = dist - pRadius, maxDist = dist + pRadius;
128
129 // clear the array
130 poResults.clear();
131
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;
139
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;
149
150 binaryStepSize /= 2;
151 }
152
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++;
159
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 ai_real 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 }
172
173 // that's it
174 }
175
176 namespace {
177
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 ai_int BinFloat;
186
187 // --------------------------------------------------------------------------------------------
188 // Converts the bit pattern of a floating-point number to its signed integer representation.
ToBinary(const ai_real & pValue)189 BinFloat ToBinary( const ai_real & pValue) {
190
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 static_assert( sizeof(BinFloat) >= sizeof(ai_real), "sizeof(BinFloat) >= sizeof(ai_real)");
195
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 static_assert( sizeof(BinFloat) == sizeof(ai_real), "sizeof(BinFloat) == sizeof(ai_real)");
202
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 ai_real 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
216
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.
220
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 }
233
234 } // namespace
235
236 // ------------------------------------------------------------------------------------------------
237 // Fills an array with indices of all positions identical to the given position. In opposite to
238 // FindPositions(), not an epsilon is used but a (very low) tolerance of four floating-point units.
FindIdenticalPositions(const aiVector3D & pPosition,std::vector<unsigned int> & poResults) const239 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.
245
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.
249
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;
264
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;
269
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.resize( 0 );
273
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;
284
285 binaryStepSize /= 2;
286 }
287
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++;
294
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 }
306
307 // that's it
308 }
309
310 // ------------------------------------------------------------------------------------------------
GenerateMappingTable(std::vector<unsigned int> & fill,ai_real pRadius) const311 unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int>& fill, ai_real pRadius) const
312 {
313 fill.resize(mPositions.size(),UINT_MAX);
314 ai_real dist, maxDist;
315
316 unsigned int t=0;
317 const ai_real pSquared = pRadius*pRadius;
318 for (size_t i = 0; i < mPositions.size();) {
319 dist = mPositions[i].mPosition * mPlaneNormal;
320 maxDist = dist + pRadius;
321
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 }
331
332 #ifdef ASSIMP_BUILD_DEBUG
333
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 }
338
339 #endif
340 return t;
341 }
342