Bug in lineToBorderPoints?

I am using the lineToBorderPoints function to find where a straight line, detected by hough, hough peaks and houghlines, crosses the border of my image. There appears to be a bug in the output of the function in that the order in which crossing coordinates are given changes depending on which image edges are crossed.

Here is a section of my code:

p=polyfit([lines.x1(1),lines.x2(1)],[lines.y1(1),lines.y2(1)],1)
coefs=[-p(1),1,-p(2)]

horizon_line=lineToBorderPoints(coefs1,size(output));

The output is 4 coordinates. An example image is below:

horizon_line =

    0.5000  295.9742  640.5000  183.8093

Given the MATLAB image coordinate system (x ascending left to right; y ascending top to bottom) we can see here that the output order is horizon_line=[x1, y1, x2, y2].

There are however specific instances where the output order changes, the images and outputs for whch follow:

1. Clockwise line rotation, line intersecting with top, but not bottom, image edge

horizon_line =

  640.5000  431.0721   43.0812    0.5000

2. Anti-clockwise line rotation, line intersecting with bottom edge (example shown), or bottom and top edge (example not shown)

  horizon_line =

  640.5000   73.9493   43.6489  480.5000

This can be easily corrected by switching the order under the circumstances that x1 is greater than x3 but unfortunately this also obscures in my case the fact that the detected horizon line sometimes goes eyond the vertical so there can be real instances where x1 should be greater than x3.

Has anyone else encountered this problem? Have they found solutions? Is this a bug in the MATLAB function? If so, what is the protocol is these circumstances?