Ray direction is not computed by the code/algorithm, I noticed this while trying to apply the code.
Fortunately it does not seem to complex, though it could get a bit tricky with overflow potential.
It's likely that a user of this code will want to setup the ray based on a start coordinate and end coordinate.
To compute ray direction I think the correct code would be something like:
Ray.Direction = Ray.End - Ray.Start;
Ray.Direction = Ray.Direction * Ray.Direction // trying to compute delta x, delta y, delta z, no sure if this code is correct.
Also then I notice another problem/short coming, where the overlow potential begins, ray length must also be computed
Not sure how to express this with the ray class.
Ray.Length = std::sqrt( Ray.Direction )
Finally the direction must be normalized.
Ray.Direction = Ray.Direction / Ray.Length
Computing ray direction can become a bit tricky.
To prevent overflow might be save to write:
RayLength = DX * DX
RayLength += DY * DY
RayLength += DZ * DZ
RayLength = std:sqrt(RayLength)
So compute it in steps... to prevent any big multiplication overflows, maybe in this case it don't matter ;)
Also it's possible that the ray length or ray direction is already normalized, it depends on the user :) but some additional code or methods to do all of this would be appreciated and takes away any doubts.