Testing three points for collinearity in a numerically robust way

You can test that three 2D points, a, b, and c, fall on a line by noting the slope of line segment (a,b) will have to be the same as that of (b,c), or by noting that the area of the triangle they define will be zero if and only if the three points are collinear. I chose the former method because the math seemed more concise: this answer contains the formula.

The trouble with the above is that although the mathematics is perfectly correct when we translate the test into code we find it behaves poorly given the fundamental imprecision of floating point types. Consider the following C++:

using point = std::tuple<float, float>;

bool are_collinear(point a, point b, point c, float eps) {
    auto [a_x, a_y] = a;
    auto [b_x, b_y] = b;
    auto [c_x, c_y] = c;

    auto test = (b_x - a_x) * (c_y - a_y) - (c_x - a_x) * (b_y - a_y);
    return std::abs(test) < eps;
}

int main()
{
    point a = { 28.8171,77.9103 };
    point b = { 55.7515,75.5051 };
    point c = { 122.831,69.8003 };

    std::cout << "are_collinear(a, b, c, 0.01) => "
        << (are_collinear(a, b, c, 0.001) ? "yes\n" : "no\n"); // no

    std::cout << "are_collinear(a, b, c, 0.1) => "
        << (are_collinear(a, b, c, 0.1) ? "yes\n" : "no\n"); // no

    std::cout << "are_collinear(a, b, c, 10) => "
        << (are_collinear(a, b, c, 10) ? "yes\n" : "no\n"); // yes
}

The three points in that code look like this:

What is happening is that test in are_collinear(...) works out to approximately 7.685 which is way outside of a typical value we would think of as an acceptable error. The problem here is that the terms in the formula are products of relative lengths in absolute coordinates e.g. (b_x - a_x) * (c_y - a_y). In order to make this function behave better we need to normalize the coordinates in some way.

Below I scale the coordinates such that the longest side of triangle (a,b,c) is 1 unit in length:

using point = std::tuple<float, float>;

float distance(point a, point b) {
    auto [a_x, a_y] = a;
    auto [b_x, b_y] = b;
    auto x_diff = a_x - b_x;
    auto y_diff = a_y - b_y;
    return std::sqrt(x_diff * x_diff + y_diff * y_diff);
}

float max_distance(point a, point b, point c) {
    auto d1 = distance(a, b);
    auto d2 = distance(b, c);
    auto d3 = distance(a, c);
    return std::max(d1, std::max(d2, d3));
}

bool are_collinear(point a, point b, point c, float eps) {
    auto scale = max_distance(a, b, c);
    if (scale == 0)
        return true; // they are the same point.

    auto [a_x, a_y] = a;
    auto [b_x, b_y] = b;
    auto [c_x, c_y] = c;

    a_x /= scale;
    a_y /= scale;
    b_x /= scale;
    b_y /= scale;
    c_x /= scale;
    c_y /= scale;

    auto test = (b_x - a_x) * (c_y - a_y) - (c_x - a_x) * (b_y - a_y);
    return std::abs(test) < eps;
}

int main()
{
    point a = { 28.8171,77.9103 };
    point b = { 55.7515,75.5051 };
    point c = { 122.831,69.8003 };

    std::cout << "are_collinear(a, b, c, 0.01) => "
        << (are_collinear(a, b, c, 0.001) ? "yes\n" : "no\n"); // yes

    std::cout << "are_collinear(a, b, c, 0.1) => "
        << (are_collinear(a, b, c, 0.1) ? "yes\n" : "no\n"); // yes

    std::cout << "are_collinear(a, b, c, 10) => "
        << (are_collinear(a, b, c, 10) ? "yes\n" : "no\n"); // yes
}

The above fixes the problem but I have two problems with it (1) it is arbitrary in the sense that I made up the scaling factor by intuition and (2) I chose the slope method originally because the code was a one-liner; the new version is substantially longer so choosing this test for concision no longer makes sense.

Is there a better way to perform this test than my second version and are there any problems with that code that I am missing?