Traversing through a line in a loop

Question

I want to find the crosses of a line (given by the start and end points) with the grid lines (x=0,1,2,3... or y=0,1,2,3,...). Thus, whether x or y of the cross points is an integer.

I came up with the following code for finding the cross points of the line given by 5.5,3.5 and 10.5,20.5 points with the grids (lines parallel to x and y axes).

#include <stdio.h>
#include <math.h>

float ax[2] = {5.0, 10.5}; // x coordinates of start and end points
float ay[2] = {3.5, 20.5}; // y coordinates of start and end points
int l, k = 0;

int main()
{
    if (ax[k] < ax[k + 1])
    {
        for (l = ceil(ax[k]); l < ax[k + 1]; l++)
        {
            printf("%d,%f\n", l, (ay[k + 1] - ay[k]) / (ax[k + 1] - ax[k]) * l + ay[k] - (ay[k + 1] - ay[k]) / (ax[k + 1] - ax[k]) * ax[k]);
        }
    }

    if (ax[k] > ax[k + 1])
    {
        for (l = ceil(ax[k + 1]); l < ax[k]; l++)
        {
            printf("%d,%f\n", l, (ay[k + 1] - ay[k]) / (ax[k + 1] - ax[k]) * l + ay[k] - (ay[k + 1] - ay[k]) / (ax[k + 1] - ax[k]) * ax[k]);
        }
    }

    if (ay[k] < ay[k + 1])
    {
        for (l = ceil(ay[k]); l < ay[k + 1]; l++)
        {
            printf("%f,%d\n", (ax[k + 1] - ax[k]) / (ay[k + 1] - ay[k]) * l + ax[k] - (ax[k + 1] - ax[k]) / (ay[k + 1] - ay[k]) * ay[k], l);
        }
    }

    if (ay[k] > ay[k + 1])
    {
        for (l = ceil(ay[k]); l < ay[k + 1]; l++)
        {
            printf("%f,%d\n", (ax[k + 1] - ax[k]) / (ay[k + 1] - ay[k]) * l + ax[k] - (ax[k + 1] - ax[k]) / (ay[k + 1] - ay[k]) * ay[k], l);
        }
    }

    return 0;
}

but I need to get the crosses in order (scanning from the start to the end of the line). I can sort the crosses found, but I think it would be easier to do this in a single loop.

The expected output is

5,3.499999
5.161765,4
5.485294,5
5.808824,6
6,6.590909
6.132353,7
6.455883,8
6.779412,9
7,9.681819
7.102942,10
7.426471,11
7.750000,12
8,12.772727
8.073530,13
8.397058,14
8.720589,15
9,15.863635
9.044118,16
9.367647,17
9.691177,18
10,18.954544
10.014706,19
10.338236,20

As you can see, the cross points have been found in order of distance from the start point (5.0,3.5).

Answer 1

Here's the full solution to your problem. Let's say we're currently at coordinate(x,y), then we calculate the offsets of next integers from the x, y and call them dx, dy. Now, the next nearest point will be of the two possible points:

1. (x+dx, Some-corresponding-Y) or
2. (SomeX, y+dx)

Now our challenge is to find the which one of the closest one. The slope (Dy/Dx) will actually help us find the next nearest one on the line.

If dx offset gives the next nearest point on the line, the corresponding offset of y coordinate, dy, can be calculated by dy = minv * dx. Conversely, if dy is the offset which is nearest, then we can calculate dx with dx = m * dy. So, we actually write an if statement to find which of the two possible alternatives is the nearest one. That's it. Once we find the next nearest point, we repeat the same process, until we reach the end point. But, beware of the limitations of the floating point calculations below. When Dx or Dy is nearly equal to 0 then then no.of valid points will become infinite. So, you need to define an upper bound on the maximum points you want to calculate. So, I didn't handle those cases (Dx ~ 0 or Dy ~ 0) in the code below.

#include <math.h>
#include <stdio.h>
#include <stdbool.h>
#include <stdlib.h>
#include <float.h>

float ax[2] = {5.0, 10.5}; // x coordinates of start and end points
float ay[2] = {3.5, 20.5}; // y coordinates of start and end points

double add(double x, double dx)
{
    return x + dx;
}

double sub(double x, double dx)
{
    return x - dx;
}

typedef double (*nextint) (double);
typedef double (*next) (double, double);

int main()
{
    double x, x1, x2, dx, Dx;
    double y, y1, y2, dy, Dy;
    double m = 0, minv = 0;

    bool dirX, dirY;

    nextint nextintx, nextinty;
    next nextx, nexty;

    x1 = ax[0];
    x2 = ax[1];
    y1 = ay[0];
    y2 = ay[1];

    dirX = x2 > x1;
    dirY = y2 > y1;

    // Tries to find the next integer of x/y
    // depending on their direction
    nextintx = dirX ? ceil : floor;
    nextinty = dirY ? ceil : floor;

    Dx = fabs(x2 - x1);
    Dy = fabs(y2 - y1);
    
    if (Dx == 0 && Dy == 0) {
    //Only a single point. Not a line segment
            if (ceil(x1) == x1 || ceil(y1) == y1)
                    printf("%f, %f\n", x1, y1);
            return 0;
    }

    if (Dx && Dy) {
            m = Dy/Dx; // slope
            minv = Dx/Dy; // Inverse of slope
    }

    // Tries to find the next point on x/y depending
    // on their direction (+ve or -ve)
    nextx = dirX ? add : sub;
    nexty = dirY ? add : sub;

    x = x1;
    y = y1;

#define morex(x)  (dirX ? (x <= x2) : (x >= x2))
#define morey(y)  (dirY ? (y <= y2) : (y >= y2))
#define more(x,y) (morex(x)) && (morey(y))

    while (more(x,y)) {
        // dx, dy values track the offsets from
        // current (x, y) coordinates respectively where
        // one of the two (x+dx, m*(x+dx)) or (minv*(y+dy), y+dy)
        // can be the desired coordinates(with integeral x/y)
        dx = fabs(nextintx(x) - x);
        dy = fabs(nextinty(y) - y);

        // We found the required (x,y) pair
        // and update their
        if (dx == 0 || dy == 0) {
            // Print the coordinates
            printf("%f, %f\n", x, y);

            // Possible offset for x/y to get integers
            dx = (dx == 0) ? 1 : dx;
            dy = (dy == 0) ? 1 : dy;
        } 

        // This is the main logic of this program.
        // The next possible pair can occur at either (x+dx, someY)
        // or (someX, y+dy). The below condition essentially checks
        // which is the closest one to the current (x,y) cooridnate.
        if (dx * Dy <= dy * Dx) {
            x  = nextx(x, dx);
            dy = dx * m;
            y = nexty(y, dy);
        } else {
            y = nexty(y, dy);
            dx = dy * minv;
            x = nextx(x, dx);
        }
    }
}