174 float c1 = m.
data[0].x*m.
data[1].y -
183 if (std::abs(c0) < std::numeric_limits<float>::epsilon())
187 const float s_inv3 = 1.0f/3.0f;
188 const float s_sqrt3 = sqrtf (3.0f);
191 float c2_over_3 = c2 * s_inv3;
192 float a_over_3 = (c1 - c2 * c2_over_3) * s_inv3;
196 float half_b = 0.5f * (c0 + c2_over_3 * (2.0f * c2_over_3 * c2_over_3 - c1));
198 float q = half_b * half_b + a_over_3 * a_over_3 * a_over_3;
203 float rho = sqrtf (-a_over_3);
204 float theta = std::atan2 (sqrtf (-q), half_b) * s_inv3;
205 float cos_theta = std::cos (theta);
206 float sin_theta = sin (theta);
207 roots.x = c2_over_3 + 2.f * rho * cos_theta;
208 roots.y = c2_over_3 - rho * (cos_theta + s_sqrt3 * sin_theta);
209 roots.z = c2_over_3 - rho * (cos_theta - s_sqrt3 * sin_theta);
212 if (roots.x >= roots.y)
213 swap (roots.x, roots.y);
214 if (roots.y >= roots.z)
216 swap (roots.y, roots.z);
217 if (roots.x >= roots.y)
218 swap (roots.x, roots.y);
229 evals = evecs.
data[0] = evecs.
data[1] = evecs.
data[2] = make_float3 (0.0f, 0.0f, 0.0f);
235 float3 scale_tmp = fmaxf (fmaxf (fabs (mat.
data[0]), fabs (mat.
data[1])), fabs (mat.
data[2]));
236 float scale = fmaxf (fmaxf (scale_tmp.x, scale_tmp.y), scale_tmp.z);
237 if (scale <= std::numeric_limits<float>::min())
241 scaledMat.
data[0] = mat.
data[0] / scale;
242 scaledMat.
data[1] = mat.
data[1] / scale;
243 scaledMat.
data[2] = mat.
data[2] / scale;
248 if ((evals.z - evals.x) <= std::numeric_limits<float>::epsilon())
251 evecs.
data[0] = make_float3 (1.0f, 0.0f, 0.0f);
252 evecs.
data[1] = make_float3 (0.0f, 1.0f, 0.0f);
253 evecs.
data[2] = make_float3 (0.0f, 0.0f, 1.0f);
255 else if ((evals.y - evals.x) <= std::numeric_limits<float>::epsilon())
263 tmp.
data[0].x -= evals.z;
264 tmp.
data[1].y -= evals.z;
265 tmp.
data[2].z -= evals.z;
267 float3 vec1 = cross (tmp.
data[0], tmp.
data[1]);
268 float3 vec2 = cross (tmp.
data[0], tmp.
data[2]);
269 float3 vec3 = cross (tmp.
data[1], tmp.
data[2]);
271 float len1 = dot (vec1, vec1);
272 float len2 = dot (vec2, vec2);
273 float len3 = dot (vec3, vec3);
275 if (len1 >= len2 && len1 >= len3)
276 evecs.
data[2] = vec1 / sqrtf (len1);
277 else if (len2 >= len1 && len2 >= len3)
278 evecs.
data[2] = vec2 / sqrtf (len2);
280 evecs.
data[2] = vec3 / sqrtf (len3);
285 else if ((evals.z - evals.y) <= std::numeric_limits<float>::epsilon())
292 tmp.
data[0].x -= evals.x;
293 tmp.
data[1].y -= evals.x;
294 tmp.
data[2].z -= evals.x;
296 float3 vec1 = cross (tmp.
data[0], tmp.
data[1]);
297 float3 vec2 = cross (tmp.
data[0], tmp.
data[2]);
298 float3 vec3 = cross (tmp.
data[1], tmp.
data[2]);
300 float len1 = dot (vec1, vec1);
301 float len2 = dot (vec2, vec2);
302 float len3 = dot (vec3, vec3);
304 if (len1 >= len2 && len1 >= len3)
305 evecs.
data[0] = vec1 / sqrtf (len1);
306 else if (len2 >= len1 && len2 >= len3)
307 evecs.
data[0] = vec2 / sqrtf (len2);
309 evecs.
data[0] = vec3 / sqrtf (len3);
320 tmp.
data[0].x -= evals.z;
321 tmp.
data[1].y -= evals.z;
322 tmp.
data[2].z -= evals.z;
324 float3 vec1 = cross (tmp.
data[0], tmp.
data[1]);
325 float3 vec2 = cross (tmp.
data[0], tmp.
data[2]);
326 float3 vec3 = cross (tmp.
data[1], tmp.
data[2]);
328 float len1 = dot (vec1, vec1);
329 float len2 = dot (vec2, vec2);
330 float len3 = dot (vec3, vec3);
333 unsigned int min_el = 2;
334 unsigned int max_el = 2;
335 if (len1 >= len2 && len1 >= len3)
338 evecs.
data[2] = vec1 / sqrtf (len1);
340 else if (len2 >= len1 && len2 >= len3)
343 evecs.
data[2] = vec2 / sqrtf (len2);
348 evecs.
data[2] = vec3 / sqrtf (len3);
354 tmp.
data[0].x -= evals.y;
355 tmp.
data[1].y -= evals.y;
356 tmp.
data[2].z -= evals.y;
358 vec1 = cross (tmp.
data[0], tmp.
data[1]);
359 vec2 = cross (tmp.
data[0], tmp.
data[2]);
360 vec3 = cross (tmp.
data[1], tmp.
data[2]);
362 len1 = dot (vec1, vec1);
363 len2 = dot (vec2, vec2);
364 len3 = dot (vec3, vec3);
365 if (len1 >= len2 && len1 >= len3)
368 evecs.
data[1] = vec1 / sqrtf (len1);
369 min_el = len1 <= mmax[min_el]? 1: min_el;
370 max_el = len1 > mmax[max_el]? 1: max_el;
372 else if (len2 >= len1 && len2 >= len3)
375 evecs.
data[1] = vec2 / sqrtf (len2);
376 min_el = len2 <= mmax[min_el]? 1: min_el;
377 max_el = len2 > mmax[max_el]? 1: max_el;
382 evecs.
data[1] = vec3 / sqrtf (len3);
383 min_el = len3 <= mmax[min_el]? 1: min_el;
384 max_el = len3 > mmax[max_el]? 1: max_el;
390 tmp.
data[0].x -= evals.x;
391 tmp.
data[1].y -= evals.x;
392 tmp.
data[2].z -= evals.x;
394 vec1 = cross (tmp.
data[0], tmp.
data[1]);
395 vec2 = cross (tmp.
data[0], tmp.
data[2]);
396 vec3 = cross (tmp.
data[1], tmp.
data[2]);
398 len1 = dot (vec1, vec1);
399 len2 = dot (vec2, vec2);
400 len3 = dot (vec3, vec3);
401 if (len1 >= len2 && len1 >= len3)
404 evecs.
data[0] = vec1 / sqrtf (len1);
405 min_el = len3 <= mmax[min_el]? 0: min_el;
406 max_el = len3 > mmax[max_el]? 0: max_el;
408 else if (len2 >= len1 && len2 >= len3)
411 evecs.
data[0] = vec2 / sqrtf (len2);
412 min_el = len3 <= mmax[min_el]? 0: min_el;
413 max_el = len3 > mmax[max_el]? 0: max_el;
418 evecs.
data[0] = vec3 / sqrtf (len3);
419 min_el = len3 <= mmax[min_el]? 0: min_el;
420 max_el = len3 > mmax[max_el]? 0: max_el;
423 unsigned mid_el = 3 - min_el - max_el;
424 evecs.
data[min_el] = normalize (cross (evecs.
data[(min_el+1)%3], evecs.
data[(min_el+2)%3] ));
425 evecs.
data[mid_el] = normalize (cross (evecs.
data[(mid_el+1)%3], evecs.
data[(mid_el+2)%3] ));
502 cov.
data[0] = make_float3 (0.0f, 0.0f, 0.0f);
503 cov.
data[1] = make_float3 (0.0f, 0.0f, 0.0f);
504 cov.
data[2] = make_float3 (0.0f, 0.0f, 0.0f);
506 cov = transform_reduce (begin, end,
517 cov.
data[0] /= (float) (end - begin);
518 cov.
data[1] /= (float) (end - begin);
519 cov.
data[2] /= (float) (end - begin);