Not correctly solving certain quartics #5
Description
Under some situations, I have found that the solver is outputting all zeros, below are some coefficients that caused it for me. I have tried putting some of these into wolfram alpha, which gives non-zero solutions.
I was using the solver to render a torus in a ray-tracing application for an assignment.
I am currently working around this issue by detecting the all-zeros condition, and retrying with slightly different numbers. As this is relatively rare, it hasn't been causing a noticable performance impact.
1
-198.995
15006.6
-508127
6.52025e+06
1
-200.246
15131.5
-511320
6.52025e+06
1
-202.193
15327.2
-516281
6.52025e+06
1
-202.205
15328.4
-516313
6.52025e+06
1
-202.215
15329.4
-516336
6.52025e+06
1
-202.268
15334.8
-516476
6.52025e+06
1
-202.289
15336.8
-516525
6.52025e+06
1
-202.35
15343
-516680
6.52025e+06
1
-202.389
15347
-516781
6.52025e+06
1
-202.388
15346.8
-516775
6.52025e+06
1
-202.427
15350.8
-516876
6.52025e+06
1
-202.427
15350.8
-516876
6.52025e+06
code producing the unwanted behaviour:
//find closest intersection with a ray originating at posn, and traveling in unit direction d
float Torus::intersect(glm::vec3 posn, glm::vec3 d)
{
glm::vec3 p = posn - this->center;
//checking for intersection with bounding sphere for optimization reasons
float len = glm::length(posn);
glm::vec3 finalPos = posn + d*len;
if (glm::length(finalPos) > (r1 + d1))
{
return -1.0;
}
float a = glm::dot(d, d);
float b = 2 * glm::dot(p, d);
float c = glm::dot(p, p) + powf(d1, 2) - powf(r1, 2);
float alpha = powf(d.x, 2) + powf(d.z, 2);
float beta = 2 * (p.x*d.x + p.z*d.z);
float gamma = powf(p.x, 2) + powf(p.z, 2);
Eigen::VectorXd polynomial(5);
polynomial[0] = powf(a, 2);
polynomial[1] = 2 * a*b;
polynomial[2] = 2 * a*c + powf(b, 2) - 4 * powf(d1, 2) *alpha;
polynomial[3] = 2 * b*c - 4 * powf(d1, 2) *beta;
polynomial[4] = powf(c, 2) - 4 * powf(d1, 2) *gamma;
Eigen::VectorXd *realRoots = new Eigen::VectorXd();
Eigen::VectorXd *compRoots = new Eigen::VectorXd();
rpoly_plus_plus::FindPolynomialRootsJenkinsTraub(polynomial, realRoots, compRoots);
float closest = -1.0;
bool nonzero = 0;
for (int i = 0; (*realRoots).size() > i; i++)
{
double d = (*realRoots).data()[i];
double c = (*compRoots).data()[i];
nonzero = nonzero || (d != 0 || c != 0);
if (fabs(c) < 1e-7 && d > 1e-7 && (d < closest || closest == -1.0))
{
closest = d;
}
}
if (!nonzero)
{
std::cout << polynomial << "\n" << "\n";
//return(this->intersect(glm::vec3(posn.x+0.0001, posn.y, posn.z), d));
}
//std::cout<<closest<<"\n";
return closest;
}
program output(set to display surface normals for debugging purposes), the pixel sized holes represent the locations where the solver returned zeros:
