sciwizeh 62 Posting Pro in Training

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;
}

}