1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
|
// This file is part of meshoptimizer library; see meshoptimizer.h for version/license details
#include "meshoptimizer.h"
#include <assert.h>
#include <float.h>
#include <string.h>
// This work is based on:
// Fabian Giesen. Decoding Morton codes. 2009
namespace meshopt
{
// "Insert" two 0 bits after each of the 10 low bits of x
inline unsigned int part1By2(unsigned int x)
{
x &= 0x000003ff; // x = ---- ---- ---- ---- ---- --98 7654 3210
x = (x ^ (x << 16)) & 0xff0000ff; // x = ---- --98 ---- ---- ---- ---- 7654 3210
x = (x ^ (x << 8)) & 0x0300f00f; // x = ---- --98 ---- ---- 7654 ---- ---- 3210
x = (x ^ (x << 4)) & 0x030c30c3; // x = ---- --98 ---- 76-- --54 ---- 32-- --10
x = (x ^ (x << 2)) & 0x09249249; // x = ---- 9--8 --7- -6-- 5--4 --3- -2-- 1--0
return x;
}
static void computeOrder(unsigned int* result, const float* vertex_positions_data, size_t vertex_count, size_t vertex_positions_stride)
{
size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
float minv[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
float maxv[3] = {-FLT_MAX, -FLT_MAX, -FLT_MAX};
for (size_t i = 0; i < vertex_count; ++i)
{
const float* v = vertex_positions_data + i * vertex_stride_float;
for (int j = 0; j < 3; ++j)
{
float vj = v[j];
minv[j] = minv[j] > vj ? vj : minv[j];
maxv[j] = maxv[j] < vj ? vj : maxv[j];
}
}
float extent = 0.f;
extent = (maxv[0] - minv[0]) < extent ? extent : (maxv[0] - minv[0]);
extent = (maxv[1] - minv[1]) < extent ? extent : (maxv[1] - minv[1]);
extent = (maxv[2] - minv[2]) < extent ? extent : (maxv[2] - minv[2]);
float scale = extent == 0 ? 0.f : 1.f / extent;
// generate Morton order based on the position inside a unit cube
for (size_t i = 0; i < vertex_count; ++i)
{
const float* v = vertex_positions_data + i * vertex_stride_float;
int x = int((v[0] - minv[0]) * scale * 1023.f + 0.5f);
int y = int((v[1] - minv[1]) * scale * 1023.f + 0.5f);
int z = int((v[2] - minv[2]) * scale * 1023.f + 0.5f);
result[i] = part1By2(x) | (part1By2(y) << 1) | (part1By2(z) << 2);
}
}
static void computeHistogram(unsigned int (&hist)[1024][3], const unsigned int* data, size_t count)
{
memset(hist, 0, sizeof(hist));
// compute 3 10-bit histograms in parallel
for (size_t i = 0; i < count; ++i)
{
unsigned int id = data[i];
hist[(id >> 0) & 1023][0]++;
hist[(id >> 10) & 1023][1]++;
hist[(id >> 20) & 1023][2]++;
}
unsigned int sumx = 0, sumy = 0, sumz = 0;
// replace histogram data with prefix histogram sums in-place
for (int i = 0; i < 1024; ++i)
{
unsigned int hx = hist[i][0], hy = hist[i][1], hz = hist[i][2];
hist[i][0] = sumx;
hist[i][1] = sumy;
hist[i][2] = sumz;
sumx += hx;
sumy += hy;
sumz += hz;
}
assert(sumx == count && sumy == count && sumz == count);
}
static void radixPass(unsigned int* destination, const unsigned int* source, const unsigned int* keys, size_t count, unsigned int (&hist)[1024][3], int pass)
{
int bitoff = pass * 10;
for (size_t i = 0; i < count; ++i)
{
unsigned int id = (keys[source[i]] >> bitoff) & 1023;
destination[hist[id][pass]++] = source[i];
}
}
} // namespace meshopt
void meshopt_spatialSortRemap(unsigned int* destination, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
{
using namespace meshopt;
assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
assert(vertex_positions_stride % sizeof(float) == 0);
meshopt_Allocator allocator;
unsigned int* keys = allocator.allocate<unsigned int>(vertex_count);
computeOrder(keys, vertex_positions, vertex_count, vertex_positions_stride);
unsigned int hist[1024][3];
computeHistogram(hist, keys, vertex_count);
unsigned int* scratch = allocator.allocate<unsigned int>(vertex_count);
for (size_t i = 0; i < vertex_count; ++i)
destination[i] = unsigned(i);
// 3-pass radix sort computes the resulting order into scratch
radixPass(scratch, destination, keys, vertex_count, hist, 0);
radixPass(destination, scratch, keys, vertex_count, hist, 1);
radixPass(scratch, destination, keys, vertex_count, hist, 2);
// since our remap table is mapping old=>new, we need to reverse it
for (size_t i = 0; i < vertex_count; ++i)
destination[scratch[i]] = unsigned(i);
}
void meshopt_spatialSortTriangles(unsigned int* destination, const unsigned int* indices, size_t index_count, const float* vertex_positions, size_t vertex_count, size_t vertex_positions_stride)
{
using namespace meshopt;
assert(index_count % 3 == 0);
assert(vertex_positions_stride >= 12 && vertex_positions_stride <= 256);
assert(vertex_positions_stride % sizeof(float) == 0);
(void)vertex_count;
size_t face_count = index_count / 3;
size_t vertex_stride_float = vertex_positions_stride / sizeof(float);
meshopt_Allocator allocator;
float* centroids = allocator.allocate<float>(face_count * 3);
for (size_t i = 0; i < face_count; ++i)
{
unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2];
assert(a < vertex_count && b < vertex_count && c < vertex_count);
const float* va = vertex_positions + a * vertex_stride_float;
const float* vb = vertex_positions + b * vertex_stride_float;
const float* vc = vertex_positions + c * vertex_stride_float;
centroids[i * 3 + 0] = (va[0] + vb[0] + vc[0]) / 3.f;
centroids[i * 3 + 1] = (va[1] + vb[1] + vc[1]) / 3.f;
centroids[i * 3 + 2] = (va[2] + vb[2] + vc[2]) / 3.f;
}
unsigned int* remap = allocator.allocate<unsigned int>(face_count);
meshopt_spatialSortRemap(remap, centroids, face_count, sizeof(float) * 3);
// support in-order remap
if (destination == indices)
{
unsigned int* indices_copy = allocator.allocate<unsigned int>(index_count);
memcpy(indices_copy, indices, index_count * sizeof(unsigned int));
indices = indices_copy;
}
for (size_t i = 0; i < face_count; ++i)
{
unsigned int a = indices[i * 3 + 0], b = indices[i * 3 + 1], c = indices[i * 3 + 2];
unsigned int r = remap[i];
destination[r * 3 + 0] = a;
destination[r * 3 + 1] = b;
destination[r * 3 + 2] = c;
}
}
|