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