Function to find the closest points between two line segements

Question

The following function is supposed to return the points p1 on line segment a, and p2 line segment b, such that |p1 - p2| is minimized.

point_pair get_closest_points(line_segment const& a, line_segment const& b)
{
    auto const intersect = intersection(extension(a), extension(b));
    if(line_segment::valid(intersect.a))
    {
        if(line_segment::valid(intersect.b))
        {
            auto const loc_a = point_at(a, intersect.a);
            auto const loc_b = point_at(b, intersect.b);
            return point_pair{loc_a, loc_b};
        }
        else
        {
            auto const loc_b = point_at(extension(b), line_segment::clamp(intersect.b));
            auto const proj = project(extension(a), loc_b);
            auto const loc_a = point_at(extension(a), proj);
            return point_pair{loc_a, loc_b};
        }
    }
    else
    {
        if(line_segment::valid(intersect.b))
        {
            auto const loc_a = point_at(extension(a), line_segment::clamp(intersect.a));
            auto const proj = project(extension(b), loc_a);
            auto const loc_b = point_at(extension(b), proj);
            return point_pair{loc_a, loc_b};
        }
        else
        {
            auto const loc_a = point_at(extension(a), line_segment::clamp(intersect.a));
            auto const proj_a_on_b = project(extension(b), loc_a);
            if(line_segment::valid(proj_a_on_b))
            {
                return point_pair{loc_a, point_at(extension(b), proj_a_on_b)};
            }
            else
            {
                auto const loc_b = point_at(extension(b), line_segment::clamp(intersect.b));
                auto const proj_b_on_a = project(extension(a), loc_b);
                if(line_segment::valid(proj_b_on_a))
                {
                    return point_pair{point_at(extension(a), proj_b_on_a), loc_b};
                }
                else
                {
                    return point_pair{loc_a, loc_b};
                }
            }
        }
    }
}

Helper functions used:

• line extension(line_segment a) converts a line segment and extends it to a line

• line_intersection intersection(line a, line b), computes the parameters of a and b that corresponds to the closest points on the lines.

• bool line_segment::valid(line_parameter t), checks that t is a valid parameter for a line_segment (should be between 0 and 1)

• point point_at(line a, line_parameter t) returns the point on a that corresponds to t

• line_parameter line_segment::clamp(line_parameter t) clamps t to the valid range

• line_parameter project(line l, point p) projects p on l

From my testing, the function appears to work. Questions:

1. Is there some case that is not handled properly? As a precondition, a and b are assumed be non-parallel and to have non-zero length.

2. Can I get rid of some branches?

Answer 1

You should've written a pseudo code explaining what the algo does and why it works. That would've helped reviewers to understand your code easier and also perhaps you would've been able to find out that the algorithm doesn't provide the right answer in some cases. Sometimes it will give a pair that isn't minimal or even one of the points might not belong to its' segment. The following cases are not handled correctly:

• In the branch where intersect.a is valid and intersect.b is not valid, your loc_a taken as a projection of a point from b to extension(a) - it doesn't necessarily fall into a. An example of this scenario is segment_a = {(-1,0),(1,0)} segment_b = {(100,100), (101,101)}

• In case both intersect.a and intersect.b are not valid, and projections of loc_a and loc_b onto the other extending lines do not fall into their respective segments - then it is not necessary that {loc_a,loc_b} is the optimal pair. It might be the case that you need to switch loc_a or loc_b for the other endvertex of their respective segment. An example of this scenario is segment_a = {(1,0),(2,0)} segment_b = {(100,100), (101,101)}: {(2,0), (100, 100)} is optimal but the algorith returns {(1,0),(100,100)}.

Let me write an example of algorithm that is both much clearer and has no if/else branching depth.

// finds the closest point to b on the segment a
point get_closet_point(line_segment const& a, point const& b);

// computes distance between the two points
double distance(point_pair const& a);

point_pair get_closest_points(line_segment const& a, line_segment const& b)
{
    auto const intersect = intersection(extension(a), extension(b));
    
    if(line_segment::valid(intersect.a) && line_segment::valid(intersect.b))
    {
            auto const loc_a = point_at(a, intersect.a);
            auto const loc_b = point_at(b, intersect.b);
            return point_pair{loc_a, loc_b};
    }
    
    std::optional<point_pair> candidate_a, candidate_b;

    if(!line_segment::valid(intersect.a))
    {
        auto const loc_a = point_at(extension(a), line_segment::clamp(intersect.a));
        auto const loc_b = get_closet_point(b, loc_a );
        candidate_a = {loc_a, loc_b};
    }
   
    if(!line_segment::valid(intersect.b))
    {
        auto const loc_b = point_at(extension(b), line_segment::clamp(intersect.b));
        auto const loc_a = get_closet_point(a, loc_b );
        candidate_b = {loc_a, loc_b};
    }

    if(candidate_a.has_value() && candidate_b.has_value())
    {
          if(distance(*candidate_a) < distance(*candidate_b))
          {
              return *candidate_a;
          }
          else
          {
              return *candidate_b;
          }
    }
    
    if(candidate_a)
    {
         return *candidate_a;
    }
    else
    {
         return *candidate_b;
    }
}

The idea of the algorithm is as follows: if the closest points between the extending lines belong to the segments then it is clearly the answer.

Otherwise, one needs to consider only the endvertices of the segments. For each of these 4 points, search for closest point in the other segment. We get 4 candidates of point pairs, now one just needs to check the 4 candidates and find the closest pair.

It happens that it is unnecessary to check all 4 candidates. It is sufficient to consider only the segment' endvertices that are the closest to the extending lines' closest points. In case the extending lines' closest point is inside the segment then there is no need to include its' envertices into candidates.

It is not exactly easy to prove that it works but it is doable.

---

For general writing. Hope the example is helpful for you to see how code can be simplified and flattened. A repetitive piece of code can be made into a function if can be given a clear name. Or into a local lambda/closure if it has no uses outside the function or cannot be given understandable name.

The function line_intersection intersection(line a, line b). The name is very misleading. It doesn't do what intersection is supposed to do.

Also, it accepts the extending lines but the output is usable for the segments. I find this behavior odd. Hypothetically, getting the extending lines can remove the information about its' original line segments. It's like you rely on behavior that isn't inherently promised.

Also, some people to do things this way, but I generally prefer member function over global free function as it is much easier to find and use them for associated classes. If you want to keep them free, it's better to put them into a namespace that clarifies to whom they belong and it will be much easier to search for them.