Surface.prejava

/*
Okay, maybe I need a interface
to encapsulate the constraint surface,
which might be a sphere or paraboloid.
*/

#include "macros.h"

import com.donhatchsw.util.VecMath;
import com.donhatchsw.util.MyMath;

public abstract class Surface
{
    // Return the offset of the point p from the surface.
    // This can mean different things to different surfaces;
    // for example, a Sphere might define it to be the offset normal to
    // the surface, but a Paraboloid might define it to be the offset
    // in the +z direction.
    // In any case, the answer should be suitable
    // for passing to rayIntersect, to generate other points
    // at the same offset.
    //
    public abstract double offset(double p[/*3*/]);

    /*
    // XXX TODO: not sure of requirements.
    //     what will this be used for?
    //     1. when creating point: picks intersection point, offset=0, in front of eye, closest to eye.
                                   record whether it's near or far surface.
           2. when starting to drag point:
                                   determine whether point is on near or far surface
    //     2. when dragging point: if there's only one intersection, pick it
                                      (does that guarantee it's on same side as original?
                                       maybe except for singular cases?)
                                      (for singular case: just don't do it.
                                       everything is measured from original mouse-down position,
                                       so we'll recover next time.)
                                   if two, pick the one on the same side as the original
                                   if none, reflect across the horizon, result should be two points,
                                        pick the one on the *opposite* side as the original
    Possible strategies of incremental movement:
        1. measure from original mouse-down position    (this seems best)
        2. measure from previous position
    Possible strategies when miss surface:
        1. fail
        2. force (knowing beforehand that it's a hit, except for roundoff error)
        2. clamp to horizon (will recover next time, assuming measuring from original mouse-down)
        3. reflect across horizon (as many times as necessary)
        4. simply don't do anything, will recover next time
           (assuming measuring from original mouse-down position)
           (same as "fail")
    Possible strategies when graze surface:
        1. return 1 point
        2. return 2 copies of same point?
    Possible strategies when eye is on surface:
        1. pretend eye is behind surface (problematic when reflecting)
        2. pretend eye is in front of surface (problematic when reflecting)
        blech.
    Possible strategies when infinite (or bigger than some large distance from origin):
        1. use homogeneous coords, that will prevent this from happening, I think
        2. clamp to large region
        3. simply don't do anything, will recover next time
           (assuming measuring from original mouse-down position)
    Possible strategies for coords:
        1. store x,y,z  (h computed when necessary)
        2. store x,y,z,h
        3. store x,y,z,w (h computed when necessary)
        4. store x,y,z,h,w (actual coords are x/w,y/w,z/w, h)
        5. store x,y,z,h,w (actual coords are x/w,y/w,z/w, h/w)
        5. store x,y,z,h,w (actual coords are x/w,y/w,z/w, h/w^2)
    */

    public final static int MISS_FAIL = 0;
    public final static int MISS_FORCE = 1; // use this when it's a known hit or graze, to prevent missing due to roundoff error.
    public final static int MISS_CLAMP = 2;
    public final static int MISS_REFLECT = 3;

    public String missStrategyToString(int missStrategy)
    {
        if (missStrategy == MISS_FAIL) return "MISS_FAIL";
        if (missStrategy == MISS_FORCE) return "MISS_FORCE";
        if (missStrategy == MISS_CLAMP) return "MISS_CLAMP";
        if (missStrategy == MISS_REFLECT) return "MISS_REFLECT";
        CHECK(false);
        return null;
    }

    // find the points on the line p + t*v that intersects the surface,
    // in order by increasing t.
    // this is the raw intersection function, implemented by the concrete subclass.
    // rayIntersect is implemented in terms of it.
    // fills in zero, or two t's (possibly infinite), into answer_t.
    // function return value is number of points filled in.
    // NOTE, the surface is assumed to be boundaryless, so this should NOT
    // ever return an odd number of t's-- if the line grazes the surface,
    // return the same point twice.
    public abstract int lineIntersect(double answer_t[/*2*/],
                                      double offset,
                                      double p[/*3*/],
                                      double v[/*3*/],
                                      int missStrategy); // FAIL or FORCE only!

    // assuming it's known that the ray misses the surface entirely,
    // return a vector pointing from p towards the horizon,
    // as close as possible to v.  need not be normalized.
    public abstract void closestHorizonDir(double answer[/*3*/],
                                           double offset,
                                           double p[/*3*/],
                                           double v[/*3*/]);


    // find the points on the ray p + t*v (t>=0) that intersects the surface,
    // in order by increasing t.
    // if it misses, reflect v across the horizon and try again,
    // but return the points in order of decreasing t.
    // fills in zero, one, or two homogeneous points into answer.
    // function return value is number of points filled in.
    public int rayIntersect(double answer[/*2*/][/*4*/],
                            double offset,
                            double p[/*3*/],
                            double v[/*3*/],
                            int missStrategy)
    {
        System.out.println("    in Surface.rayIntersect");
        System.out.println("        offset = "+offset);
        System.out.println("        p = "+VecMath.toString(p));
        System.out.println("        v = "+VecMath.toString(v));
        System.out.println("        v normalized = "+VecMath.toString(VecMath.normalize(v)));
        System.out.println("        missStrategy = "+missStrategyToString(missStrategy));

        // when REFLECT, we keep track of change in v, requiring it decrease each time
        double tooBigDist2 = Double.POSITIVE_INFINITY;
        int nReflections = 0;

        while (true)
        {
            double ts[] = answer[1]; // use as scratch area initially
            int nAnswers = lineIntersect(ts, offset, p, v, MISS_FAIL);
            CHECK(nAnswers == 0
                || nAnswers == 2);
            int nValidAnswers = 0;
            FORI (i, nAnswers)
                if (ts[i] >= 0.
                 || (i==nAnswers-1 && missStrategy==MISS_FORCE)) // when FORCE, don't pass up last chance
                {
                    ts[nValidAnswers++] = ts[i];
                }
            FORI (i, nValidAnswers)
            {
                double t = ts[i];
                if (t == Double.POSITIVE_INFINITY)
                {
                    VecMath.copyvec(3, answer[i], v);
                    answer[i][3] = 0.;
                }
                else
                {
                    VecMath.vpsxv(3, answer[i], p, t, v);
                    answer[i][3] = 1.;
                }
                subNormalize(answer[i]);
            }
            if (nValidAnswers != 0)
            {
                if (nValidAnswers == 2)
                {
                    if (nReflections % 2 == 1)
                    {
                        double temp;
                        FORI (i, 4)
                            SWAP(answer[0][i], answer[1][i], temp);
                    }
                }
                System.out.println("    out Surface.rayIntersect (returning nValidAnswers="+nValidAnswers+")");
                return nValidAnswers;
            }
            else
            {
                //
                // nValidAnswers is 0.
                // what to do?
                //
                if (missStrategy == MISS_FAIL)
                {
                    System.out.println("    out Surface.rayIntersect (FAIL!)");
                    return 0;
                }
                else if (missStrategy == MISS_FORCE)
                {
                    CHECK(false); // can't happen
                }
                else // CLAMP or REFLECT
                {
                    // tail recurse, with either horizon ray
                    // or v reflected across horizon.

                    double vAdjusted[] = new double[3];
                    closestHorizonDir(vAdjusted, offset, p, v);

                    if (missStrategy == MISS_REFLECT)
                    {
                        // replace vAdjusted with v projected onto horizon dir
                        VecMath.vxs(vAdjusted, vAdjusted, VecMath.dot(vAdjusted,v) / VecMath.normsqrd(vAdjusted));
                        // replace vAdjusted with v reflected across horizon
                        VecMath.lerp(vAdjusted, vAdjusted, v, -1.);
                        double thisDist2 = VecMath.distsqrd(v, vAdjusted);
                        if (thisDist2 >= tooBigDist2) // didn't decrease; force the next one
                            missStrategy = MISS_FORCE;
                        tooBigDist2 = thisDist2;
                        System.out.println("        missed! reflecting");
                    }
                    else // CLAMP
                    {
                        missStrategy = MISS_FORCE;
                        System.out.println("        missed! forcing");
                    }

                    v = vAdjusted;

                    nReflections++;

                    // and tail recurse.
                }
            } // nValidAnswers==0
        } // while (true)
    } // rayIntersect


    // sample nPoints around the horizon seen by the given eye.
    public abstract double[][] horizon(double offset, double eye[/*3*/], int nPoints);




    // adjust homogeneous coords so that max abs value is in interval (.5,1].
    // this protects against overflow, without introducing any roundoff error.
    private void subNormalize(double v[])
    {
        double maxAbs = ABS(v[0]);
        for (int i = 1; i < v.length; ++i) // skip [0]
        {
            double thisAbs = ABS(v[i]);
            maxAbs = MAX(maxAbs, thisAbs);
        }
        double scaleFactor = 1.;
        if (maxAbs > 1.)
            do
            {
                scaleFactor *= .5;
                maxAbs *= .5;
            } while (maxAbs > 1.);
        else if (maxAbs <= .5)
            do
            {
                scaleFactor *= 2.;
                maxAbs *= 2.;
            } while (maxAbs <= .5);
        for (int i = 0; i < v.length; ++i)
            v[i] *= scaleFactor;
    } // subNormalize

    private static int solveQuadratic(double answers[], double A, double B, double C, int missStrategy)
    {
        // numerical recipes in c, p. 184
        double discr = B*B - 4*A*C;
        if (discr < 0.)
        {
            if (missStrategy == MISS_FORCE)
            {
                System.out.println("HEY! fudged discr "+discr+" to zero");
                discr = 0.;
            }
            else // MISS_FAIL
            {
                return 0;
            }
        }
        double q = -.5 * (B + (B < 0 ? -1 : 1) * Math.sqrt(discr));
        double tThis = q/A;
        double tThat = C/q;
        answers[0] = MIN(tThis, tThat);
        answers[1] = MAX(tThis, tThat);
        return 2;
    } // solveQuadratic


    // z = zScale * (x^2 + y^2)
    public static class Paraboloid extends Surface
    {
        public Paraboloid(double zScale)
        {
            this.zScale = zScale;
        }

        // required by Surface interface
        @Override public double offset(double p[/*3*/])
        {
            return p[2] - zScale * (p[0]*p[0] + p[1]*p[1]);
        }

        // required by Surface interface
        @Override public void closestHorizonDir(double answer[], double offset, double p[/*3*/], double v[/*3*/])
        {
            System.out.println("        in Surface.Paraboloid.closestHorizonDir");
            if (this.zScale == 0.)
            {
                // project v onto xy plane, that's it.
                // could just zero out z component,
                // but then we risk getting zero,
                // so use robust normalization instead.

                double basis[][] = new double[2][3];
                basis[0][2] = 1.; // z axis
                VecMath.copyvec(basis[1], v);
                VecMath.extendAndOrthogonalize(2, 2, basis, basis);
                VecMath.copyvec(answer, basis[1]);
                System.out.println("        out Surface.Paraboloid.closestHorizonDir (zScale=0)");
                return;
            }

            // XXX TODO: the following method is not robust if very close to surface
            double px = p[0];
            double py = p[1];
            double pz = p[2];
            double vx = v[0];
            double vy = v[1];
            double vz = v[2];
            // The horizon curve, projected onto the xy plane,
            // is, magically, a circle of radius sqrt(-(z-zScale*(x^2+y^2)-offset)/zScale)
            // centered at the eye.
            double horizonRadius2 = -(pz-zScale*(px*px+py*py)-offset)/zScale;
            CHECK(GEQ(horizonRadius2, 0., 1e-3)); // if we were inside sphere, the above would have succeeded
            double horizonRadius = Math.sqrt(MAX(horizonRadius2,0.));
            double vLenXY = MyMath.hypot(vx,vy); // length of v projected to xy plane
            // XXX TODO: not robust if v is vertical!
            double horizonPoint[] = new double[] {
                px + vx*(horizonRadius/vLenXY),
                py + vy*(horizonRadius/vLenXY),
                0. // will be filled in in a moment
            };
            horizonPoint[2] = this.zScale*(SQR(horizonPoint[0])+SQR(horizonPoint[1])) + offset;

            if (false)
            {
                // TODO: the following is wrong wrong wrong. we don't get the right answer by just orthogonalizing v against the normal.

                // At this point we could just return horizonPoint-p,
                // but that's not robust if p is very close to surface.
                // so instead, we return v orthogonalized against the normal
                // at horizonPoint.
                //
                // Surface normal (not necessarily unit length)
                // of the paraboloid at x,y,(x^2+y^2)*zScale
                // is -2*zScale*x,-2*zScale*y,1.
                //
                double basis[][] = new double[2][3];
                basis[0][0] = -2*this.zScale*horizonPoint[0];
                basis[0][1] = -2*this.zScale*horizonPoint[1];
                basis[0][1] = 1.;
                VecMath.copyvec(basis[1], v);
                VecMath.extendAndOrthogonalize(2, 2, basis, basis);
                VecMath.copyvec(answer, basis[1]);
            }
            else
            {
                // do it simple non-robust way for now
                VecMath.vmv(answer, horizonPoint, p);
            }
            System.out.println("        out Surface.Paraboloid.closestHorizonDir (normal)");
        } // closestHorizonDir

        @Override public int lineIntersect(double answer_t[/*2*/],
                                           double offset,
                                           double p[/*3*/],
                                           double v[/*3*/],
                                           int missStrategy) // FAIL or FORCE only!
        {
            double px = p[0];
            double py = p[1];
            double pz = p[2];
            double vx = v[0];
            double vy = v[1];
            double vz = v[2];

            if (vx == 0. && vy == 0.)
            {
                // solving quadratic would blow up...
                // just return the point on the offset surface at the same x,y az p,
                // and also the infinite point at the far end of the paraboloid.
                // XXX TODO: actually, does it blow up properly? i.e. returning the correct infinities. if so, no need to special case it.
                if (vz < 0.)
                {
                    answer_t[0] = (offset-pz) / -vz;
                    answer_t[1] = Double.POSITIVE_INFINITY;
                }
                else
                {
                    answer_t[0] = Double.NEGATIVE_INFINITY;
                    answer_t[1] = (offset-pz) / -vz;
                }
                return 2;
            }

            // want t such that p+t*v is on the offset paraboloid, i.e.
            //     pz+t*vz = zScale * ((px+t*vx)^2 + (py+t*vy)^2) + offset
            //             = zScale * (px^2 + 2*t*px*vx + t^2*vx^2 + py^2 + 2*t*py*vy + t^2*vy^2) + offset
            //             = zScale * ((px^2+py^2) + 2*t*(px*vx + py*vy) + t^2*(vx^2 + vy^2)) + offset
            //     zScale*(vx^2+vy^2)*t^2 + (2*zScale*(px*vx+py*vy)-vz)*t + (zScale*(px^2+py^2)-pz+offset) == 0
            // This is a quadratic equation:
            double A = zScale * (vx*vx + vy*vy);
            double B = 2 * zScale * (px*vx + py*vy) - vz;
            double C = zScale * (px*px + py*py) - pz + offset; // note, could be computed outside loop since it doesn't depend on v... but that's awkward now that we're inside a required virtual function
            return solveQuadratic(answer_t, A, B, C, missStrategy);

        } // lineIntersect

        // required by Surface interface
        @Override public double[][] horizon(double offset, double eye[/*3*/], int nPoints)
        {
            double x = eye[0];
            double y = eye[1];
            double z = eye[2];

            double points[][] = new double[nPoints][3];

            if (zScale == 0.)
            {
                // horizon is infinite; draw something in the *direction* of the horizon.

                double radius = 10.; // somewhat arbitrary

                FORI (iPoint, nPoints)
                {
                    double point[] = points[iPoint];
                    double angle = 2*Math.PI*iPoint/nPoints;
                    point[0] = x + radius * Math.cos(angle);
                    point[1] = y + radius * Math.sin(angle);
                    point[2] = z;
                }
            }
            else
            {
                // The horizon curve, projected onto the xy plane,
                // is, magically, a circle of radius sqrt(-(z-zScale*(x^2+y^2)-offset)/zScale)
                // centered at the eye.
                double radius2 = -(z-zScale*(x*x+y*y)-offset)/zScale;
                double radius = Math.sqrt(MAX(radius2,0.));

                FORI (iPoint, nPoints)
                {
                    double point[] = points[iPoint];
                    double angle = 2*Math.PI*iPoint/nPoints;
                    point[0] = x + radius * Math.cos(angle);
                    point[1] = y + radius * Math.sin(angle);
                    point[2] = zScale*(SQR(point[0])+SQR(point[1])) + offset;
                }
            }
            return points;
        } // horizon

        public double zScale;
    } // class Paraboloid

    // ||p - center|| == radius
    public static class Sphere extends Surface
    {
        public Sphere(double center[/*3*/], double radius)
        {
            this.center = center;
            this.radius = radius;
        } // ctor

        // required by Surface interface
        @Override public double offset(double p[/*3*/])
        {
            return VecMath.dist(this.center, p) - this.radius;
        }

        // required by Surface interface
        @Override public void closestHorizonDir(double answer[], double offset, double p[/*3*/], double v[/*3*/])
        {
            System.out.println("        in Surface.Sphere.closestHorizonDir");
            double r = this.radius + offset;
            double r2 = r*r;
            double P2 = VecMath.distsqrd(this.center, p);
            double centerToHorizonChordMidpoint2 = r2*r2/P2;
            double centerToHorizonChordMidpoint = Math.sqrt(centerToHorizonChordMidpoint2);
            //PRINT(centerToHorizonChordMidpoint2);
            //PRINT(Math.sqrt(centerToHorizonChordMidpoint2));
            double horizonChordHalfLength2 = r2 - centerToHorizonChordMidpoint2;

            // orthogonalize v against p-center to get horizon chord dir
            double basis[][] = new double[2][3];
            VecMath.vmv(basis[0], p, this.center);
            VecMath.copyvec(basis[1], v);
            VecMath.extendAndOrthogonalize(2, 2, basis, basis);
            // basis[1] is now horizon chord dir

            if (horizonChordHalfLength2 <= 0.)
            {
                // if equal to zero, p is right on the surface.
                // if less than zero, that means we're inside the sphere,
                // which means we shouldn't have gotten this far but we did
                // due to roundoff error, so we still consider p
                // to be right on the surface.
                CHECK_GE(horizonChordHalfLength2, -SQR(1e-3));
                VecMath.copyvec(answer, basis[1]); // horizon chord dir
                System.out.println("        out Surface.Sphere.closestHorizonDir (FUDGED!)");
                return;
            }
            double horizonChordHalfLength = Math.sqrt(horizonChordHalfLength2);

            // answer is v orthogonalized against centerToHorizonPoint
            VecMath.sxvpsxv(basis[0],
                            centerToHorizonChordMidpoint, basis[0],
                            horizonChordHalfLength, basis[1]);
            // basis[0] is now centerToHorizonPoint
            VecMath.copyvec(basis[1], v);
            VecMath.extendAndOrthogonalize(2, 2, basis, basis);
            VecMath.copyvec(answer, basis[1]);
            System.out.println("        out Surface.Sphere.closestHorizonDir (normal)");
        } // closestHorizonDir

        @Override public int lineIntersect(double answer_t[/*2*/],
                                           double offset,
                                           double p[/*3*/],
                                           double v[/*3*/],
                                           int missStrategy) // FAIL or FORCE only!
        {
            double r = this.radius + offset;
            double r2 = r*r;
            double P2 = VecMath.distsqrd(this.center, p);
            double pToHorizon2 = P2 - r2; // distance squared between p (the eye) and horizon

            // let P = p-c.
            // want t such that ||P + t*v|| == r
            // i.e. (P dot P) + t^2 (v dot v) + 2*t*(P dot v) == r^2
            // This is a quadratic equation with:
            //      A = v dot v
            //      B = 2 (P dot v)
            //      C = (P dot P) - r^2 (happens to be pToHorizon2)
            double A = VecMath.dot(v, v);
            double B = 2 * (VecMath.dot(p,v) - VecMath.dot(this.center,v)); // 2*((p-c) dot v)
            double C = pToHorizon2; // note, doesn't depend on v

            return solveQuadratic(answer_t, A, B, C, missStrategy);

        } // lineIntersect

        // required by Surface interface
        @Override public double[][] horizon(double offset, double eye[/*3*/], int nPoints)
        {
            //System.out.println("    in Sphere.horizon");

            double dist2 = VecMath.distsqrd(this.center, eye);
            double r = this.radius+offset;
            double r2 = r*r;

            double horizonRadius = r * Math.sqrt(MAX(1. - r2/dist2, 0.));
            double horizonCenter[] = VecMath.lerp(this.center,
                                                  eye,
                                                  r2/dist2);
            double basis[][] = new double[3][3];
            VecMath.vmv(basis[0], this.center, eye);
            VecMath.extendAndOrthogonalize(1, 3, basis, basis);

            double points[][] = new double[nPoints][3];
            FORI (iPoint, nPoints)
            {
                double angle = 2*Math.PI*iPoint/nPoints;
                VecMath.vpsxvpsxv(points[iPoint],
                                  horizonCenter,
                                  horizonRadius*Math.cos(angle), basis[1],
                                  horizonRadius*Math.sin(angle), basis[2]);
            }
            //System.out.println("    out Sphere.horizon");
            return points;
        } // horizon

        public double center[];
        public double radius;
    } // class Sphere

    public static class Box extends Surface
    {
        public Box(double xRadius, double yRadius, double zRadius)
        {
            this.r = new double[] {xRadius, yRadius, zRadius};
        }

        // required by Surface interface
        @Override public double offset(double p[/*3*/])
        {
            CHECK(false); // XXX IMPLEMENT ME
            return 0.;
        }

        // required by Surface interface
        @Override public void closestHorizonDir(double answer[], double offset, double p[/*3*/], double v[/*3*/])
        {
            CHECK(false); // XXX IMPLEMENT ME
        }

        @Override public int lineIntersect(double answer_t[/*2*/],
                                           double offset,
                                           double p[/*3*/],
                                           double v[/*3*/],
                                           int missStrategy) // FAIL or FORCE only!
        {
            CHECK(false); // XXX IMPLEMENT ME
            return 0;
        }

        // required by Surface interface
        @Override public double[][] horizon(double offset, double eye[/*3*/], int nPoints)
        {
            CHECK(false); // XXX IMPLEMENT ME!
            return null;
        } // horizon

        public double r[/*3*/];
    } // class Box


    //
    // Given a Surface,
    // this class returns another Surface
    // representing the original surface
    // transformed in the given way.
    // The transform must be a rigid xform
    // with optional uniform scale.
    //
    public static class XformedSurface extends Surface
    {
        public XformedSurface(Surface originalSurface,
                              double s,     // scale
                              double r[][], // rotation/reflection
                              double t[])   // translation
        {
            this.originalSurface = originalSurface;
            this.s = s;
            this.r = r!=null ? r : VecMath.identitymat(3);
            this.t = t!=null ? t : new double[3]; // zeros

            CHECK_EQ(r.length, 3);
            CHECK_EQ(r[0].length, 3);
            CHECK(VecMath.equals(VecMath.mxm(r, VecMath.transpose(r)), VecMath.identitymat(3), 1e-6));
            CHECK_EQ(t.length, 3);

            this.M = new double[][] {
                {s*r[0][0], s*r[0][1], s*r[0][2]},
                {s*r[1][0], s*r[1][1], s*r[1][2]},
                {s*r[2][0], s*r[2][1], s*r[2][2]},
                {t[0],      t[1],      t[2],    },
            };
            this.invM = VecMath.invertmat(this.M);
            CHECK_EQ(this.invM.length, 4);
            CHECK_EQ(this.invM[0].length, 3);
        } // ctor

        // required by Surface interface
        @Override public double offset(double p[/*3*/])
        {
            // xform vert from world space into original surface's space
            double originalP[] = VecMath.vxm(p, invM);

            // compute offset in original
            double originalOffset = this.originalSurface.offset(originalP);

            // offset is the original offset, scaled back into world space.
            double offset = this.s * originalOffset;

            return offset;
        } // computeNormalAndOffset

        // required by Surface interface
        @Override public double[][] horizon(double offset, double eye[/*3*/], int nPoints)
        {
            double originalEye[] = VecMath.vxm(eye, this.invM);
            double originalOffset = offset / s;
            double horizon[][] = this.originalSurface.horizon(originalOffset, originalEye, nPoints);
            FORI (iPoint, nPoints)
                horizon[iPoint] = VecMath.vxm(horizon[iPoint], this.M);
            return horizon;
        } // horizon

        // required by Surface interface
        @Override public void closestHorizonDir(double answer[], double offset, double p[/*3*/], double v[/*3*/])
        {
            CHECK(false); // XXX IMPLEMENT ME
        }

        // required by Surface interface
        @Override public int lineIntersect(double answer_t[/*2*/],
                                           double offset,
                                           double p[/*3*/],
                                           double v[/*3*/],
                                           int missStrategy) // FAIL or FORCE only!
        {
            double invS = 1./s;
            double originalOffset = offset * invS;
            double originalP[] = VecMath.vxm(p, invM);
            double originalV[] = VecMath.mxv(r, v); // apply inverse of r
            VecMath.vxs(originalV, originalV, invS); // and inverse of s
            return this.originalSurface.lineIntersect(answer_t,
                                                      originalOffset,
                                                      originalP,
                                                      originalV,
                                                      missStrategy);
        }

        // private since public access would let them get out of sync
        private Surface originalSurface;
        private double s;
        private double r[/*3*/][/*3*/];
        private double t[/*3*/];
        private double M[/*4*/][/*4*/];    // s r t
        private double invM[/*4*/][/*4*/]; // inverse of M
    } // class XformedSurface

} // class Surface