mathematics - The correct way to transform a ray with a matrix?

Playing with XNA Triangle Picking Sample I found out that it does not work well if you scale the world matrix of the objects you want to pick. When I dug into the implementation I found this comment in the RayIntersectsModel method:

        // The input ray is in world space, but our model data is stored in object
        // space. We would normally have to transform all the model data by the
        // modelTransform matrix, moving it into world space before we test it
        // against the ray. That transform can be slow if there are a lot of
        // triangles in the model, however, so instead we do the opposite.

        // Transforming our ray by the inverse modelTransform moves it into object
        // space, where we can test it directly against our model data. Since there
        // is only one ray but typically many triangles, doing things this way
        // around can be much faster.

After the comment they actually transformed the ray:

ray.Position = Vector3.Transform(ray.Position, inverseTransform);
ray.Direction = Vector3.TransformNormal(ray.Direction, inverseTransform);

With this implementation, picking suffered from "short-sightedness" if you scaled the models: it could only pick those objects, that were close enough to it. Even the ray-boundingSphere intersection test, which implementation is hardcoded into XNA, failed in the same way.

I fixed this by "doing the wrong thing" - I actually started transforming every vertex by the model's world matrix and to fix the boundingSphere test I added this code:

Quaternion rot;
Vector3 scale, trans;
modelTransform.Decompose(out scale, out rot, out trans);

float maxScale = Math.Max(Math.Max(scale.X, scale.Y), scale.Z);

boundingSphere.Center = Vector3.Transform(boundingSphere.Center, modelTransform);
boundingSphere.Radius *= maxScale;

This obviously is not optimal and I wanted to know if there is a way to actually transform the ray back to the model's space with the model's inverted matrix, while making it work for scaled matrices?

SOLUTION: Thanks to Nathan's answer I found a way to fix the ray scaling - just renormalize the ray direction:

ray.Position = Vector3.Transform(ray.Position, inverseTransform);
ray.Direction = Vector3.TransformNormal(ray.Direction, inverseTransform);
//ADD THE FOLLOWING LINE:
ray.Direction.Normalize();

SOLUTION UPDATE: As I tested the app, I found that Nathan was indeed completely right and another change was necessary. Here is the full code for the correct RayIntersectsModel() method:

static float? RayIntersectsModel(Ray ray, Model model, Matrix modelTransform,
                                out bool insideBoundingSphere,

                                out Vector3 vertex1, out Vector3 vertex2,
                                out Vector3 vertex3)
{
    vertex1 = vertex2 = vertex3 = Vector3.Zero;
...
    Matrix inverseTransform = Matrix.Invert(modelTransform);
    // STORE WORLDSPACE RAY.
    Ray oldRay = ray;

    ray.Position = Vector3.Transform(ray.Position, inverseTransform);

    ray.Direction = Vector3.TransformNormal(ray.Direction, inverseTransform);
    ray.Direction.Normalize();

    // Look up our custom collision data from the Tag property of the model.
    Dictionary tagData = (Dictionary)model.Tag;

    if (tagData == null)
    {
        throw new InvalidOperationException(
            "Model.Tag is not set correctly. Make sure your model " +

            "was built using the custom TrianglePickingProcessor.");
    }

    // Start off with a fast bounding sphere test.
    BoundingSphere boundingSphere = (BoundingSphere)tagData["BoundingSphere"];

    if (boundingSphere.Intersects(ray) == null)
    {
        // If the ray does not intersect the bounding sphere, we cannot
        // possibly have picked this model, so there is no need to even

        // bother looking at the individual triangle data.
        insideBoundingSphere = false;

        return null;
    }
    else
    {
        // The bounding sphere test passed, so we need to do a full
        // triangle picking test.
        insideBoundingSphere = true;


        // Keep track of the closest triangle we found so far,
        // so we can always return the closest one.
        float? closestIntersection = null;

        // Loop over the vertex data, 3 at a time (3 vertices = 1 triangle).
        Vector3[] vertices = (Vector3[])tagData["Vertices"];

        for (int i = 0; i < vertices.Length; i += 3)
        {

            // Perform a ray to triangle intersection test.
            float? intersection;

            RayIntersectsTriangle(ref ray,
                                    ref vertices[i],
                                    ref vertices[i + 1],
                                    ref vertices[i + 2],
                                    out intersection);

            // Does the ray intersect this triangle?

            if (intersection != null)
            {
                // RECOMPUTE DISTANCE IN WORLD SPACE:
                Vector3 vertexA = Vector3.Transform(vertices[i], modelTransform);
                Vector3 vertexB = Vector3.Transform(vertices[i+1], modelTransform);
                Vector3 vertexC = Vector3.Transform(vertices[i+2], modelTransform);

                RayIntersectsTriangle(ref oldRay,
                                    ref vertexA,
                                    ref vertexB,

                                    ref vertexC,
                                    out intersection);

                // If so, is it closer than any other previous triangle?
                if ((closestIntersection == null) ||
                    (intersection < closestIntersection))
                {
                    // Store the distance to this triangle.
                    closestIntersection = intersection;


                    // Store the three vertex positions in world space.
                    vertex1 = vertexA;
                    vertex2 = vertexB;
                    vertex3 = vertexC;
                }
            }
        }

        return closestIntersection;
    }

}

Answer

Transforming the ray position and direction by the inverse model transformation is correct. However, many ray-intersection routines assume that the ray direction is a unit vector. If the model transformation involves scaling, the ray direction won't be a unit vector afterward, and should likely be renormalized.

However, the distance along the ray returned by the intersection routines will then be measured in model space, and won't represent the distance in world space. If it's a uniform scale, you can simply multiply the returned distance by the scale factor to convert it back to world-space distance. For non-uniform scaling it's trickier; probably the best way is to transform the intersection point back to world space and then re-measure the distance from the ray origin there.