@@ -145,67 +145,82 @@ static void buildTriangleAdjacencySparse(TriangleAdjacency2& adjacency, const un
145145 }
146146}
147147
148- static void computeBoundingSphere (float result[4 ], const float points[][ 3 ] , size_t count)
148+ static void computeBoundingSphere (float result[4 ], const float * points, size_t count, size_t points_stride, const float * radii, size_t radii_stride )
149149{
150150 assert (count > 0 );
151151
152+ size_t points_stride_float = points_stride / sizeof (float );
153+ size_t radii_stride_float = radii_stride / sizeof (float );
154+
152155 // find extremum points along all 3 axes; for each axis we get a pair of points with min/max coordinates
153156 size_t pmin[3 ] = {0 , 0 , 0 };
154157 size_t pmax[3 ] = {0 , 0 , 0 };
155158
156159 for (size_t i = 0 ; i < count; ++i)
157160 {
158- const float * p = points[i];
161+ const float * p = points + i * points_stride_float;
162+ float r = radii[i * radii_stride_float];
159163
160164 for (int axis = 0 ; axis < 3 ; ++axis)
161165 {
162- pmin[axis] = (p[axis] < points[pmin[axis]][axis]) ? i : pmin[axis];
163- pmax[axis] = (p[axis] > points[pmax[axis]][axis]) ? i : pmax[axis];
166+ float bmin = points[pmin[axis] * points_stride_float + axis] - radii[pmin[axis] * radii_stride_float];
167+ float bmax = points[pmax[axis] * points_stride_float + axis] + radii[pmax[axis] * radii_stride_float];
168+
169+ pmin[axis] = (p[axis] - r < bmin) ? i : pmin[axis];
170+ pmax[axis] = (p[axis] + r > bmax) ? i : pmax[axis];
164171 }
165172 }
166173
167174 // find the pair of points with largest distance
168- float paxisd2 = 0 ;
169175 int paxis = 0 ;
176+ float paxisdr = 0 ;
170177
171178 for (int axis = 0 ; axis < 3 ; ++axis)
172179 {
173- const float * p1 = points[pmin[axis]];
174- const float * p2 = points[pmax[axis]];
180+ const float * p1 = points + pmin[axis] * points_stride_float;
181+ const float * p2 = points + pmax[axis] * points_stride_float;
182+ float r1 = radii[pmin[axis] * radii_stride_float];
183+ float r2 = radii[pmax[axis] * radii_stride_float];
175184
176185 float d2 = (p2[0 ] - p1[0 ]) * (p2[0 ] - p1[0 ]) + (p2[1 ] - p1[1 ]) * (p2[1 ] - p1[1 ]) + (p2[2 ] - p1[2 ]) * (p2[2 ] - p1[2 ]);
186+ float dr = sqrtf (d2) + r1 + r2;
177187
178- if (d2 > paxisd2 )
188+ if (dr > paxisdr )
179189 {
180- paxisd2 = d2 ;
190+ paxisdr = dr ;
181191 paxis = axis;
182192 }
183193 }
184194
185195 // use the longest segment as the initial sphere diameter
186- const float * p1 = points[pmin[paxis]];
187- const float * p2 = points[pmax[paxis]];
196+ const float * p1 = points + pmin[paxis] * points_stride_float;
197+ const float * p2 = points + pmax[paxis] * points_stride_float;
198+ float r1 = radii[pmin[paxis] * radii_stride_float];
199+ float r2 = radii[pmax[paxis] * radii_stride_float];
200+
201+ float paxisd = sqrtf ((p2[0 ] - p1[0 ]) * (p2[0 ] - p1[0 ]) + (p2[1 ] - p1[1 ]) * (p2[1 ] - p1[1 ]) + (p2[2 ] - p1[2 ]) * (p2[2 ] - p1[2 ]));
202+ float paxisk = paxisd > 0 ? (paxisd + r2 - r1) / (2 * paxisd) : 0 .f ;
188203
189- float center[3 ] = {( p1[0 ] + p2[0 ]) / 2 , ( p1[1 ] + p2[1 ]) / 2 , ( p1[2 ] + p2[2 ]) / 2 };
190- float radius = sqrtf (paxisd2) / 2 ;
204+ float center[3 ] = {p1[0 ] + ( p2[0 ] - p1[ 0 ]) * paxisk, p1[1 ] + ( p2[1 ] - p1[ 1 ]) * paxisk, p1[2 ] + ( p2[2 ] - p1[ 2 ]) * paxisk };
205+ float radius = paxisdr / 2 ;
191206
192207 // iteratively adjust the sphere up until all points fit
193208 for (size_t i = 0 ; i < count; ++i)
194209 {
195- const float * p = points[i];
210+ const float * p = points + i * points_stride_float;
211+ float r = radii[i * radii_stride_float];
212+
196213 float d2 = (p[0 ] - center[0 ]) * (p[0 ] - center[0 ]) + (p[1 ] - center[1 ]) * (p[1 ] - center[1 ]) + (p[2 ] - center[2 ]) * (p[2 ] - center[2 ]);
214+ float d = sqrtf (d2);
197215
198- if (d2 > radius * radius)
216+ if (d + r > radius)
199217 {
200- float d = sqrtf (d2);
201- assert (d > 0 );
202-
203- float k = 0 .5f + (radius / d) / 2 ;
218+ float k = d > 0 ? (d + r - radius) / (2 * d) : 0 .f ;
204219
205- center[0 ] = center[ 0 ] * k + p[0 ] * ( 1 - k );
206- center[1 ] = center[ 1 ] * k + p[1 ] * ( 1 - k );
207- center[2 ] = center[ 2 ] * k + p[2 ] * ( 1 - k );
208- radius = (radius + d) / 2 ;
220+ center[0 ] += k * ( p[0 ] - center[ 0 ] );
221+ center[1 ] += k * ( p[1 ] - center[ 1 ] );
222+ center[2 ] += k * ( p[2 ] - center[ 2 ] );
223+ radius = (radius + d + r ) / 2 ;
209224 }
210225 }
211226
@@ -1004,15 +1019,17 @@ meshopt_Bounds meshopt_computeClusterBounds(const unsigned int* indices, size_t
10041019 if (triangles == 0 )
10051020 return bounds;
10061021
1022+ const float rzero = 0 .f ;
1023+
10071024 // compute cluster bounding sphere; we'll use the center to determine normal cone apex as well
10081025 float psphere[4 ] = {};
1009- computeBoundingSphere (psphere, corners[0 ], triangles * 3 );
1026+ computeBoundingSphere (psphere, corners[0 ][ 0 ] , triangles * 3 , sizeof ( float ) * 3 , &rzero, 0 );
10101027
10111028 float center[3 ] = {psphere[0 ], psphere[1 ], psphere[2 ]};
10121029
10131030 // treating triangle normals as points, find the bounding sphere - the sphere center determines the optimal cone axis
10141031 float nsphere[4 ] = {};
1015- computeBoundingSphere (nsphere, normals, triangles);
1032+ computeBoundingSphere (nsphere, normals[ 0 ] , triangles, sizeof ( float ) * 3 , &rzero, 0 );
10161033
10171034 float axis[3 ] = {nsphere[0 ], nsphere[1 ], nsphere[2 ]};
10181035 float axislength = sqrtf (axis[0 ] * axis[0 ] + axis[1 ] * axis[1 ] + axis[2 ] * axis[2 ]);
@@ -1122,6 +1139,33 @@ meshopt_Bounds meshopt_computeMeshletBounds(const unsigned int* meshlet_vertices
11221139 return meshopt_computeClusterBounds (indices, triangle_count * 3 , vertex_positions, vertex_count, vertex_positions_stride);
11231140}
11241141
1142+ meshopt_Bounds meshopt_computeSphereBounds (const float * positions, size_t count, size_t positions_stride, const float * radii, size_t radii_stride)
1143+ {
1144+ using namespace meshopt ;
1145+
1146+ assert (positions_stride >= 12 && positions_stride <= 256 );
1147+ assert (positions_stride % sizeof (float ) == 0 );
1148+ assert ((radii_stride >= 4 && radii_stride <= 256 ) || radii == NULL );
1149+ assert (radii_stride % sizeof (float ) == 0 );
1150+
1151+ meshopt_Bounds bounds = {};
1152+
1153+ if (count == 0 )
1154+ return bounds;
1155+
1156+ const float rzero = 0 .f ;
1157+
1158+ float psphere[4 ] = {};
1159+ computeBoundingSphere (psphere, positions, count, positions_stride, radii ? radii : &rzero, radii ? radii_stride : 0 );
1160+
1161+ bounds.center [0 ] = psphere[0 ];
1162+ bounds.center [1 ] = psphere[1 ];
1163+ bounds.center [2 ] = psphere[2 ];
1164+ bounds.radius = psphere[3 ];
1165+
1166+ return bounds;
1167+ }
1168+
11251169void meshopt_optimizeMeshlet (unsigned int * meshlet_vertices, unsigned char * meshlet_triangles, size_t triangle_count, size_t vertex_count)
11261170{
11271171 using namespace meshopt ;
0 commit comments