I'm porting some code from a book about ray tracing from C++ to Java, because I like java better and it seems like a decent way to learn how the code works. I ported the code for the Torus Primitive and all seemed well, until I tried to view it with an orthographic Z axis aligned camera and I got major errors. Image I can't find anything that should be wrong, perhaps because I don't entirely understand the Quartic formula, or perhaps I've missed something in the C++ implementation that needs to be changed to work with Java. I think the problem may be where solveCubic(..) is called, the original implementation ignored the returned number of roots, but I changed it to int shouldBeOne= so that I could look at it It is not one usually and we ignore more a root, and I don't know if that's a mistake or not. If this isn't a good place for this I'll have to head somewhere else.
public class Torus extends GeometricObject {
private double a;//outer circle center distance
private double b;//swept circle radius
private BBox bbox;
public Torus() {
this(2, 0.5);
}
public Torus(double a, double b) {
super();
this.a = a;
this.b = b;
bbox = new BBox(-a - b, a + b, -b, b, -a - b, a + b);
}
//...OTHER STUFF
public Normal computeNormal(Point3D p) {
Normal normal = new Normal();
double param_squared = a * a + b * b;
double x = p.x;
double y = p.y;
double z = p.z;
double sum_squared = x * x + y * y + z * z;
normal.x = 4.0 * x * (sum_squared - param_squared);
normal.y = 4.0 * y * (sum_squared - param_squared + 2.0 * a * a);
normal.z = 4.0 * z * (sum_squared - param_squared);
normal.normalize();
return (normal);
}
@Override
public boolean hit(Ray ray, ShadeRec s) {
if (!bbox.hit(ray)) {
return (false);
}
double x1 = ray.o.x;
double y1 = ray.o.y;
double z1 = ray.o.z;
double d1 = ray.d.x;
double d2 = ray.d.y;
double d3 = ray.d.z;
double[] coeffs = new double[5]; // coefficient array for the quartic equation
double[] roots = new double[4]; // solution array for the quartic equation
// define the coefficients of the quartic equation
double sum_d_sqrd = d1 * d1 + d2 * d2 + d3 * d3;
double e = x1 * x1 + y1 * y1 + z1 * z1 - a * a - b * b;
double f = x1 * d1 + y1 * d2 + z1 * d3;
double four_a_sqrd = 4.0 * a * a;
double E = e * e - four_a_sqrd * (b * b - y1 * y1); // constant term
double D = 4.0 * f * e + 2.0 * four_a_sqrd * y1 * d2;
double C = 2.0 * sum_d_sqrd * e + 4.0 * f * f + four_a_sqrd * d2 * d2;
double B = 4.0 * sum_d_sqrd * f;
double A = sum_d_sqrd * sum_d_sqrd; // coefficient of t^4
coeffs[0] = E;
coeffs[1] = D;
coeffs[2] = C;
coeffs[3] = B;
coeffs[4] = A;
boolean intersected = false;
// find roots of the quartic equation
int num_real_roots = Utility.solveQuartic(coeffs, roots);
double t = Utility.HUGE_VALUE;
if (num_real_roots == 0) {
return (false);
}
// find the smallest root greater than kEpsilon, if any
// the roots array is not sorted
for (int j = 0; j < num_real_roots; j++) {
if (roots[j] > Utility.EPSILON) {
intersected = true;
if (roots[j] < t) {
t = roots[j];
}
}
}
if (!intersected) {
return (false);
}
s.lastT = t;
s.localHitPosition.setTo(ray.o.add(ray.d.mul(t)));
s.normal.setTo(computeNormal(s.localHitPosition));
return (true);
}
//...OTHER STUFF
}
class Utility {
public static final double PI = Math.PI;
public static final double TWO_PI = 2 * PI;
public static final double PI_ON_180 = PI / 180.0;
public static final double INV_PI = 1.0 / PI;
public static final double INV_2_PI = 1.0 / TWO_PI;
public static final double EPSILON = 0.0001;
public static final double HUGE_VALUE = 1.0e10;
public static final boolean isZero(double x) {
return x > -EPSILON && x < EPSILON;
}
public static final int SolveQuadric(double c[], double s[]) {
double p, q, D;
/* normal form: x^2 + px + q = 0 */
p = c[ 1 ] / (2 * c[ 2 ]);
q = c[ 0 ] / c[ 2 ];
D = p * p - q;
if (isZero(D)) {
s[ 0 ] = - p;
return 1;
}
else if (D > 0) {
double sqrt_D = Math.sqrt(D);
s[ 0 ] = sqrt_D - p;
s[ 1 ] = - sqrt_D - p;
return 2;
}
else /* if (D < 0) */
return 0;
}
public static final int solveCubic(double c[], double s[]) {
int i, num;
double sub;
double A, B, C;
double sq_A, p, q;
double cb_p, D;
/* normal form: x^3 + Ax^2 + Bx + C = 0 */
A = c[ 2 ] / c[ 3 ];
B = c[ 1 ] / c[ 3 ];
C = c[ 0 ] / c[ 3 ];
/* substitute x = y - A/3 to eliminate quadric term:
x^3 +px + q = 0 */
sq_A = A * A;
p = 1.0/3 * (- 1.0/3 * sq_A + B);
q = 1.0/2 * (2.0/27 * A * sq_A - 1.0/3 * A * B + C);
/* use Cardano's formula */
cb_p = p * p * p;
D = q * q + cb_p;
if (isZero(D)) {
if (isZero(q)) { /* one triple solution */
s[ 0 ] = 0;
num = 1;
}
else { /* one single and one double solution */
double u = Math.cbrt(-q);
s[ 0 ] = 2 * u;
s[ 1 ] = - u;
num = 2;
}
}
else if (D < 0) { /* Casus irreducibilis: three real solutions */
double phi = 1.0/3 * Math.acos(-q / Math.sqrt(-cb_p));
double t = 2 * Math.sqrt(-p);
s[ 0 ] = t * Math.cos(phi);
s[ 1 ] = - t * Math.cos(phi + PI / 3);
s[ 2 ] = - t * Math.cos(phi - PI / 3);
num = 3;
}
else { /* one real solution */
double sqrt_D = Math.sqrt(D);
double u = Math.cbrt(sqrt_D - q);
double v = - Math.cbrt(sqrt_D + q);
s[ 0 ] = u + v;
num = 1;
}
/* resubstitute */
sub = 1.0/3 * A;
for (i = 0; i < num; ++i)
s[ i ] -= sub;
return num;
}
public static final int solveQuartic(double c[], double s[]) {
double[] coeffs=new double[4];
double z, u, v, sub;
double A, B, C, D;
double sq_A, p, q, r;
int i, num;
/* normal form: x^4 + Ax^3 + Bx^2 + Cx + D = 0 */
A = c[ 3 ] / c[ 4 ];
B = c[ 2 ] / c[ 4 ];
C = c[ 1 ] / c[ 4 ];
D = c[ 0 ] / c[ 4 ];
/* substitute x = y - A/4 to eliminate cubic term:
x^4 + px^2 + qx + r = 0 */
sq_A = A * A;
p = - 3.0/8 * sq_A + B;
q = 1.0/8 * sq_A * A - 1.0/2 * A * B + C;
r = - 3.0/256*sq_A*sq_A + 1.0/16*sq_A*B - 1.0/4*A*C + D;
if (isZero(r)) {
/* no absolute term: y(y^3 + py + q) = 0 */
coeffs[ 0 ] = q;
coeffs[ 1 ] = p;
coeffs[ 2 ] = 0;
coeffs[ 3 ] = 1;
num = solveCubic(coeffs, s);
s[ num++ ] = 0;
}
else {
/* solve the resolvent cubic ... */
coeffs[ 0 ] = 1.0/2 * r * p - 1.0/8 * q * q;
coeffs[ 1 ] = - r;
coeffs[ 2 ] = - 1.0/2 * p;
coeffs[ 3 ] = 1;
int shouldBeOne = solveCubic(coeffs, s);
/* ... and take the one real solution ... */
z = s[ 0 ];
/* ... to build two quadric equations */
u = z * z - r;
v = 2 * z - p;
if (isZero(u))
u = 0;
else if (u > 0)
u = Math.sqrt(u);
else
return 0;
if (isZero(v))
v = 0;
else if (v > 0)
v = Math.sqrt(v);
else
return 0;
coeffs[ 0 ] = z - u;
coeffs[ 1 ] = q < 0 ? -v : v;
coeffs[ 2 ] = 1;
num = SolveQuadric(coeffs, s);
coeffs[ 0 ]= z + u;
coeffs[ 1 ] = q < 0 ? v : -v;
coeffs[ 2 ] = 1;
int tempNum = 0;
double[] tempSols = new double[2];
tempNum = SolveQuadric(coeffs, tempSols);
for(int in =0; in<tempNum; in++){
s[in+num]=tempSols[in];
}
num+=tempNum;
}
/* resubstitute */
sub = 1.0/4 * A;
for (i = 0; i < num; ++i)
s[ i ] -= sub;
return num;
}
}