Fast Voxel Traversal Algorithm with negative direction

Question

I'm trying to implement fast voxel traversal algorithm and calculate T and M according to this answer (T is tDelta, M is tMax). All is good if the two components of the direction vector V are positive. But if at least one of them is negative, it's work wrong.

Green point is start, red is end. All seems correct.

And now from bigger to less position.

Traversal method:

private IEnumerable<Vector2> GetCrossedCells(Vector2 pPoint1, Vector2 pPoint2)
{
    Vector2 V = pPoint2 - pPoint1; // direction & distance vector
    if (V != Vector2.Zero)
    {
        Vector2 U = Vector2.Normalize(V); // direction unit vector
        Vector2 S = new Vector2(Math.Sign(U.X), Math.Sign(U.Y)); // sign vector
        Vector2 P = pPoint1; // position
        Vector2 G = new Vector2((int) Math.Floor(P.X / CELL_SIZE), (int) Math.Floor(P.Y / CELL_SIZE)); // grid coord
        Vector2 T = new Vector2(Math.Abs(CELL_SIZE / U.X), Math.Abs(CELL_SIZE / U.Y));
        Vector2 M = new Vector2(
            Single.IsInfinity(T.X) ? Single.PositiveInfinity : T.X * (1.0f - (P.X / CELL_SIZE) % 1),
            Single.IsInfinity(T.Y) ? Single.PositiveInfinity : T.Y * (1.0f - (P.Y / CELL_SIZE) % 1));

        Vector2 D = Vector2.Zero;
        bool isCanMoveByX = S.X != 0;
        bool isCanMoveByY = S.Y != 0;

        while (isCanMoveByX || isCanMoveByY)
        {
            yield return G;

            D = new Vector2(
                S.X > 0 ? (float) (Math.Floor(P.X / CELL_SIZE) + 1) * CELL_SIZE - P.X :
                S.X < 0 ? (float) (Math.Ceiling(P.X / CELL_SIZE) - 1) * CELL_SIZE - P.X :
                0,
                S.Y > 0 ? (float) (Math.Floor(P.Y / CELL_SIZE) + 1) * CELL_SIZE - P.Y :
                S.Y < 0 ? (float) (Math.Ceiling(P.Y / CELL_SIZE) - 1) * CELL_SIZE - P.Y :
                0);

            if (Math.Abs(V.X) <= Math.Abs(D.X))
            {
                D.X = V.X;
                isCanMoveByX = false;
            }

            if (Math.Abs(V.Y) <= Math.Abs(D.Y))
            {
                D.Y = V.Y;
                isCanMoveByY = false;
            }

            if (M.X <= M.Y)
            {
                M.X += T.X;
                G.X += S.X;
                if (isCanMoveByY)
                {
                    D.Y = U.Y / U.X * D.X; // U.X / U.Y = D.X / D.Y => U.X * D.Y = U.Y * D.X
                }
            }
            else
            {
                M.Y += T.Y;
                G.Y += S.Y;
                if (isCanMoveByX)
                {
                    D.X = U.X / U.Y * D.Y;
                }
            }

            V -= D;
            P += D;
        }
    }
}

In debug I can see that for example M.Y > M.X when should be the opposite if S.X < 0 or S.Y < 0.

Tell me please what my code work wrong for negative directions?

Answer 1

So, I solve it. I make code cleaner and problem is gone.

private IEnumerable<Vector2> GetCrossedCells(Vector2 pPoint1, Vector2 pPoint2)
{
    if (pPoint1 != pPoint2)
    {
        Vector2 V = (pPoint2 - pPoint1) / CELL_SIZE; // direction & distance vector
        Vector2 U = Vector2.Normalize(V); // direction unit vector
        Vector2 S = new Vector2(Math.Sign(U.X), Math.Sign(U.Y)); // sign vector
        Vector2 P = pPoint1 / CELL_SIZE; // position in grid coord system
        Vector2 G = new Vector2((int) Math.Floor(P.X), (int) Math.Floor(P.Y)); // grid coord
        Vector2 T = new Vector2(Math.Abs(CELL_SIZE / U.X), Math.Abs(CELL_SIZE / U.Y));
        Vector2 D = new Vector2(
            S.X > 0 ? 1 - P.X % 1 : S.X < 0 ? P.X % 1 : 0,
            S.Y > 0 ? 1 - P.Y % 1 : S.Y < 0 ? P.Y % 1 : 0);
        Vector2 M = new Vector2(
            Single.IsInfinity(T.X) || S.X == 0 ? Single.PositiveInfinity : T.X * D.X,
            Single.IsInfinity(T.Y) || S.Y == 0 ? Single.PositiveInfinity : T.Y * D.Y);

        bool isCanMoveByX = S.X != 0;
        bool isCanMoveByY = S.Y != 0;

        while (isCanMoveByX || isCanMoveByY)
        {
            yield return G;

            D = new Vector2(
                S.X > 0 ? (float) Math.Floor(P.X) + 1 - P.X :
                S.X < 0 ? (float) Math.Ceiling(P.X) - 1 - P.X :
                0,
                S.Y > 0 ? (float) Math.Floor(P.Y) + 1 - P.Y :
                S.Y < 0 ? (float) Math.Ceiling(P.Y) - 1 - P.Y :
                0);

            if (Math.Abs(V.X) <= Math.Abs(D.X))
            {
                D.X = V.X;
                isCanMoveByX = false;
            }

            if (Math.Abs(V.Y) <= Math.Abs(D.Y))
            {
                D.Y = V.Y;
                isCanMoveByY = false;
            }

            if (M.X <= M.Y)
            {
                M.X += T.X;
                G.X += S.X;
                if (isCanMoveByY)
                {
                    D.Y = U.Y / U.X * D.X; // U.X / U.Y = D.X / D.Y => U.X * D.Y = U.Y * D.X
                }
            }
            else
            {
                M.Y += T.Y;
                G.Y += S.Y;
                if (isCanMoveByX)
                {
                    D.X = U.X / U.Y * D.Y;
                }
            }

            V -= D;
            P += D;
        }
    }
}

Update

I'm began from removing redundant divisions on GRID_CELL and then notice mistake in M calculation. There are using Frac() function in answer to the question, a link to which I provided. I'm calculate it like (1 - P % 1), but that is a case for S > 0, and there are should be (P % 1) if S < 0, and Inf for S = 0.

Update 2

Also there should be

Vector2 D = new Vector2(
    S.X > 0 ? (float) Math.Floor(P.X) + 1 - P.X :
    S.X < 0 ? (float) Math.Ceiling(P.X) - 1 - P.X :
    0,
    S.Y > 0 ? (float) Math.Floor(P.Y) + 1 - P.Y :
    S.Y < 0 ? (float) Math.Ceiling(P.Y) - 1 - P.Y :
    0);

Instead of

Vector2 D = new Vector2(
    S.X > 0 ? 1 - P.X % 1 : S.X < 0 ? P.X % 1 : 0,
    S.Y > 0 ? 1 - P.Y % 1 : S.Y < 0 ? P.Y % 1 : 0);

Because M will be infinity in case S < 0 and P haven't fractional part.