1 /*
2 ---------------------------------------------------------------------------
3 Open Asset Import Library (assimp)
4 ---------------------------------------------------------------------------
5
6 Copyright (c) 2006-2021, assimp team
7
8 All rights reserved.
9
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:
13
14 * Redistributions of source code must retain the above
15 copyright notice, this list of conditions and the
16 following disclaimer.
17
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.
22
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.
27
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 */
41
42 /** @file Implementation of the helper class to quickly find vertices close to a given position */
43
44 #include <assimp/SpatialSort.h>
45 #include <assimp/ai_assert.h>
46
47 using namespace Assimp;
48
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
53
54 const aiVector3D PlaneInit(0.8523f, 0.34321f, 0.5736f);
55
56 // ------------------------------------------------------------------------------------------------
57 // Constructs a spatially sorted representation from the given position array.
58 // define the reference plane. We choose some arbitrary vector away from all basic axes
59 // in the hope that no model spreads all its vertices along this plane.
SpatialSort(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset)60 SpatialSort::SpatialSort(const aiVector3D *pPositions, unsigned int pNumPositions, unsigned int pElementOffset) :
61 mPlaneNormal(PlaneInit),
62 mFinalized(false) {
63 mPlaneNormal.Normalize();
64 Fill(pPositions, pNumPositions, pElementOffset);
65 }
66
67 // ------------------------------------------------------------------------------------------------
SpatialSort()68 SpatialSort::SpatialSort() :
69 mPlaneNormal(PlaneInit),
70 mFinalized(false) {
71 mPlaneNormal.Normalize();
72 }
73
74 // ------------------------------------------------------------------------------------------------
75 // Destructor
~SpatialSort()76 SpatialSort::~SpatialSort() {
77 // empty
78 }
79
80 // ------------------------------------------------------------------------------------------------
Fill(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset,bool pFinalize)81 void SpatialSort::Fill(const aiVector3D *pPositions, unsigned int pNumPositions,
82 unsigned int pElementOffset,
83 bool pFinalize /*= true */) {
84 mPositions.clear();
85 mFinalized = false;
86 Append(pPositions, pNumPositions, pElementOffset, pFinalize);
87 mFinalized = pFinalize;
88 }
89
90 // ------------------------------------------------------------------------------------------------
CalculateDistance(const aiVector3D & pPosition) const91 ai_real SpatialSort::CalculateDistance(const aiVector3D &pPosition) const {
92 return (pPosition - mCentroid) * mPlaneNormal;
93 }
94
95 // ------------------------------------------------------------------------------------------------
Finalize()96 void SpatialSort::Finalize() {
97 const ai_real scale = 1.0f / mPositions.size();
98 for (unsigned int i = 0; i < mPositions.size(); i++) {
99 mCentroid += scale * mPositions[i].mPosition;
100 }
101 for (unsigned int i = 0; i < mPositions.size(); i++) {
102 mPositions[i].mDistance = CalculateDistance(mPositions[i].mPosition);
103 }
104 std::sort(mPositions.begin(), mPositions.end());
105 mFinalized = true;
106 }
107
108 // ------------------------------------------------------------------------------------------------
Append(const aiVector3D * pPositions,unsigned int pNumPositions,unsigned int pElementOffset,bool pFinalize)109 void SpatialSort::Append(const aiVector3D *pPositions, unsigned int pNumPositions,
110 unsigned int pElementOffset,
111 bool pFinalize /*= true */) {
112 ai_assert(!mFinalized && "You cannot add positions to the SpatialSort object after it has been finalized.");
113 // store references to all given positions along with their distance to the reference plane
114 const size_t initial = mPositions.size();
115 mPositions.reserve(initial + pNumPositions);
116 for (unsigned int a = 0; a < pNumPositions; a++) {
117 const char *tempPointer = reinterpret_cast<const char *>(pPositions);
118 const aiVector3D *vec = reinterpret_cast<const aiVector3D *>(tempPointer + a * pElementOffset);
119 mPositions.push_back(Entry(static_cast<unsigned int>(a + initial), *vec));
120 }
121
122 if (pFinalize) {
123 // now sort the array ascending by distance.
124 Finalize();
125 }
126 }
127
128 // ------------------------------------------------------------------------------------------------
129 // Returns an iterator for all positions close to the given position.
FindPositions(const aiVector3D & pPosition,ai_real pRadius,std::vector<unsigned int> & poResults) const130 void SpatialSort::FindPositions(const aiVector3D &pPosition,
131 ai_real pRadius, std::vector<unsigned int> &poResults) const {
132 ai_assert(mFinalized && "The SpatialSort object must be finalized before FindPositions can be called.");
133 const ai_real dist = CalculateDistance(pPosition);
134 const ai_real minDist = dist - pRadius, maxDist = dist + pRadius;
135
136 // clear the array
137 poResults.clear();
138
139 // quick check for positions outside the range
140 if (mPositions.size() == 0)
141 return;
142 if (maxDist < mPositions.front().mDistance)
143 return;
144 if (minDist > mPositions.back().mDistance)
145 return;
146
147 // do a binary search for the minimal distance to start the iteration there
148 unsigned int index = (unsigned int)mPositions.size() / 2;
149 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
150 while (binaryStepSize > 1) {
151 if (mPositions[index].mDistance < minDist)
152 index += binaryStepSize;
153 else
154 index -= binaryStepSize;
155
156 binaryStepSize /= 2;
157 }
158
159 // depending on the direction of the last step we need to single step a bit back or forth
160 // to find the actual beginning element of the range
161 while (index > 0 && mPositions[index].mDistance > minDist)
162 index--;
163 while (index < (mPositions.size() - 1) && mPositions[index].mDistance < minDist)
164 index++;
165
166 // Mow start iterating from there until the first position lays outside of the distance range.
167 // Add all positions inside the distance range within the given radius to the result array
168 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
169 const ai_real pSquared = pRadius * pRadius;
170 while (it->mDistance < maxDist) {
171 if ((it->mPosition - pPosition).SquareLength() < pSquared)
172 poResults.push_back(it->mIndex);
173 ++it;
174 if (it == mPositions.end())
175 break;
176 }
177
178 // that's it
179 }
180
181 namespace {
182
183 // Binary, signed-integer representation of a single-precision floating-point value.
184 // IEEE 754 says: "If two floating-point numbers in the same format are ordered then they are
185 // ordered the same way when their bits are reinterpreted as sign-magnitude integers."
186 // This allows us to convert all floating-point numbers to signed integers of arbitrary size
187 // and then use them to work with ULPs (Units in the Last Place, for high-precision
188 // computations) or to compare them (integer comparisons are faster than floating-point
189 // comparisons on many platforms).
190 typedef ai_int BinFloat;
191
192 // --------------------------------------------------------------------------------------------
193 // Converts the bit pattern of a floating-point number to its signed integer representation.
ToBinary(const ai_real & pValue)194 BinFloat ToBinary(const ai_real &pValue) {
195
196 // If this assertion fails, signed int is not big enough to store a float on your platform.
197 // Please correct the declaration of BinFloat a few lines above - but do it in a portable,
198 // #ifdef'd manner!
199 static_assert(sizeof(BinFloat) >= sizeof(ai_real), "sizeof(BinFloat) >= sizeof(ai_real)");
200
201 #if defined(_MSC_VER)
202 // If this assertion fails, Visual C++ has finally moved to ILP64. This means that this
203 // code has just become legacy code! Find out the current value of _MSC_VER and modify
204 // the #if above so it evaluates false on the current and all upcoming VC versions (or
205 // on the current platform, if LP64 or LLP64 are still used on other platforms).
206 static_assert(sizeof(BinFloat) == sizeof(ai_real), "sizeof(BinFloat) == sizeof(ai_real)");
207
208 // This works best on Visual C++, but other compilers have their problems with it.
209 const BinFloat binValue = reinterpret_cast<BinFloat const &>(pValue);
210 //::memcpy(&binValue, &pValue, sizeof(pValue));
211 //return binValue;
212 #else
213 // On many compilers, reinterpreting a float address as an integer causes aliasing
214 // problems. This is an ugly but more or less safe way of doing it.
215 union {
216 ai_real asFloat;
217 BinFloat asBin;
218 } conversion;
219 conversion.asBin = 0; // zero empty space in case sizeof(BinFloat) > sizeof(float)
220 conversion.asFloat = pValue;
221 const BinFloat binValue = conversion.asBin;
222 #endif
223
224 // floating-point numbers are of sign-magnitude format, so find out what signed number
225 // representation we must convert negative values to.
226 // See http://en.wikipedia.org/wiki/Signed_number_representations.
227 const BinFloat mask = BinFloat(1) << (CHAR_BIT * sizeof(BinFloat) - 1);
228
229 // Two's complement?
230 const bool DefaultValue = ((-42 == (~42 + 1)) && (binValue & mask));
231 const bool OneComplement = ((-42 == ~42) && (binValue & mask));
232
233 if (DefaultValue)
234 return mask - binValue;
235 // One's complement?
236 else if (OneComplement)
237 return BinFloat(-0) - binValue;
238 // Sign-magnitude? -0 = 1000... binary
239 return binValue;
240 }
241
242 } // namespace
243
244 // ------------------------------------------------------------------------------------------------
245 // Fills an array with indices of all positions identical to the given position. In opposite to
246 // 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) const247 void SpatialSort::FindIdenticalPositions(const aiVector3D &pPosition, std::vector<unsigned int> &poResults) const {
248 ai_assert(mFinalized && "The SpatialSort object must be finalized before FindIdenticalPositions can be called.");
249 // Epsilons have a huge disadvantage: they are of constant precision, while floating-point
250 // values are of log2 precision. If you apply e=0.01 to 100, the epsilon is rather small, but
251 // if you apply it to 0.001, it is enormous.
252
253 // The best way to overcome this is the unit in the last place (ULP). A precision of 2 ULPs
254 // tells us that a float does not differ more than 2 bits from the "real" value. ULPs are of
255 // logarithmic precision - around 1, they are 1*(2^24) and around 10000, they are 0.00125.
256
257 // For standard C math, we can assume a precision of 0.5 ULPs according to IEEE 754. The
258 // incoming vertex positions might have already been transformed, probably using rather
259 // inaccurate SSE instructions, so we assume a tolerance of 4 ULPs to safely identify
260 // identical vertex positions.
261 static const int toleranceInULPs = 4;
262 // An interesting point is that the inaccuracy grows linear with the number of operations:
263 // multiplying to numbers, each inaccurate to four ULPs, results in an inaccuracy of four ULPs
264 // plus 0.5 ULPs for the multiplication.
265 // To compute the distance to the plane, a dot product is needed - that is a multiplication and
266 // an addition on each number.
267 static const int distanceToleranceInULPs = toleranceInULPs + 1;
268 // The squared distance between two 3D vectors is computed the same way, but with an additional
269 // subtraction.
270 static const int distance3DToleranceInULPs = distanceToleranceInULPs + 1;
271
272 // Convert the plane distance to its signed integer representation so the ULPs tolerance can be
273 // applied. For some reason, VC won't optimize two calls of the bit pattern conversion.
274 const BinFloat minDistBinary = ToBinary(CalculateDistance(pPosition)) - distanceToleranceInULPs;
275 const BinFloat maxDistBinary = minDistBinary + 2 * distanceToleranceInULPs;
276
277 // clear the array in this strange fashion because a simple clear() would also deallocate
278 // the array which we want to avoid
279 poResults.resize(0);
280
281 // do a binary search for the minimal distance to start the iteration there
282 unsigned int index = (unsigned int)mPositions.size() / 2;
283 unsigned int binaryStepSize = (unsigned int)mPositions.size() / 4;
284 while (binaryStepSize > 1) {
285 // Ugly, but conditional jumps are faster with integers than with floats
286 if (minDistBinary > ToBinary(mPositions[index].mDistance))
287 index += binaryStepSize;
288 else
289 index -= binaryStepSize;
290
291 binaryStepSize /= 2;
292 }
293
294 // depending on the direction of the last step we need to single step a bit back or forth
295 // to find the actual beginning element of the range
296 while (index > 0 && minDistBinary < ToBinary(mPositions[index].mDistance))
297 index--;
298 while (index < (mPositions.size() - 1) && minDistBinary > ToBinary(mPositions[index].mDistance))
299 index++;
300
301 // Now start iterating from there until the first position lays outside of the distance range.
302 // Add all positions inside the distance range within the tolerance to the result array
303 std::vector<Entry>::const_iterator it = mPositions.begin() + index;
304 while (ToBinary(it->mDistance) < maxDistBinary) {
305 if (distance3DToleranceInULPs >= ToBinary((it->mPosition - pPosition).SquareLength()))
306 poResults.push_back(it->mIndex);
307 ++it;
308 if (it == mPositions.end())
309 break;
310 }
311
312 // that's it
313 }
314
315 // ------------------------------------------------------------------------------------------------
GenerateMappingTable(std::vector<unsigned int> & fill,ai_real pRadius) const316 unsigned int SpatialSort::GenerateMappingTable(std::vector<unsigned int> &fill, ai_real pRadius) const {
317 ai_assert(mFinalized && "The SpatialSort object must be finalized before GenerateMappingTable can be called.");
318 fill.resize(mPositions.size(), UINT_MAX);
319 ai_real dist, maxDist;
320
321 unsigned int t = 0;
322 const ai_real pSquared = pRadius * pRadius;
323 for (size_t i = 0; i < mPositions.size();) {
324 dist = (mPositions[i].mPosition - mCentroid) * mPlaneNormal;
325 maxDist = dist + pRadius;
326
327 fill[mPositions[i].mIndex] = t;
328 const aiVector3D &oldpos = mPositions[i].mPosition;
329 for (++i; i < fill.size() && mPositions[i].mDistance < maxDist && (mPositions[i].mPosition - oldpos).SquareLength() < pSquared; ++i) {
330 fill[mPositions[i].mIndex] = t;
331 }
332 ++t;
333 }
334
335 #ifdef ASSIMP_BUILD_DEBUG
336
337 // debug invariant: mPositions[i].mIndex values must range from 0 to mPositions.size()-1
338 for (size_t i = 0; i < fill.size(); ++i) {
339 ai_assert(fill[i] < mPositions.size());
340 }
341
342 #endif
343 return t;
344 }
345