Calculate intersection points of line segments

[BeKinect]

I'm trying to calculate the intersection points of two line segments.

So I got n points stored in a vector called "boundingbox" and m points stored in a second vector called "matrix".
First I create a line segment by two points of boundingbox. The second line segment is created by two points of matrix. Than the intersection points are calculated, if they exist. The last step of my code is to check, if the found intersection point is on both line segments. So heres my code:

  int IsPointInBoundingBox(double x1, double y1, double x2, double y2, double px, double py)
{
    if(min(x1,x2) <= px && px <= max(x1,x2) && min(y1,y2) <= py  && py <= max(y1,y2))
    {
        return 1;
    }
    else
        return 0;
}

for(int a=fillingStart; a<boundingbox.size()-1; a += 2)
{
    double Ax=boundingbox[a][0];
    double Ay=boundingbox[a][1];
    double Bx=boundingbox[a+1][0];
    double By=boundingbox[a+1][1];
			
    double B1=Bx-Ax;
    double A1=By-Ay;
    double C1=A1*Ax+B1*Ay;
				
    for(size_t j=0; j<(matrix.size()-1); j++)
    {
	double Cx=matrix[j][0];
	double Cy=matrix[j][1];
	double Dx=matrix[j+1][0];
	double Dy=matrix[j+1][1];
				
	double B2=Dx-Cx;
	double A2=Dy-Cy;
	double C2=A2*Cx + B2*Cy;

	double det=A1*B2-A2*B1;
						
	if(det!=0) { // det == 0 -> Lines are parallel
	  double point_x=(B2*C1-B1*C2)/det;
	  double point_y=(A1*C2-A2*C1)/det;
	  if(IsPointInBoundingBox(Ax, Ay, Bx, By, point_x, point_y) == 1 && 
	     IsPointInBoundingBox(Cx, Cy, Dx, Dy, point_x, point_y) == 1)
	  {
	    intersects.push_back(std::array<double, 3>());
	    int q = intersects.size()-1;
	    intersects[q][0]=point_x;
	    intersects[q][1]=point_y;
	    intersects[q][2]=matrix[j][2];
	  }
	}
     }
}

The code to find the intersections is working. But the function IsPointInBoundingBox, does not work correctly. When I execute my program, there are no points in the csv file I print to. The function never become true. I already tried to calculate the intersection points manually and found points matching the criteria. I just don't know what I'm doing wrong.

[Half Life G Man]

bool IsPointInBoundingBox(double _x1, double _y1, double _x2, double _y2, double px, double py)
{
    double x1 = min(_x1, _x2);
    double x2 = max(_x1, _x2);
    double y1 = min(_y1, _y2);
    double y2 = max(_y1, _y2);

    if((px >= x1 && px <= x2) && (py >= y1 && py <= y2))
    {
        return true;
    }
    else return false;
}

[brown ogum]

I tested the function "IsPointInBoundingBox" and it works properly.
I don't seem to understand where your problem is

#include <iostream>

using namespace std;

void IsPointInBoundingBox(double x1, double y1, double x2, double y2, double px, double py)
{
    if(min(x1,x2) <= px && px <= max(x1,x2) && min(y1,y2) <= py  && py <= max(y1,y2))
    	cout<<1;
    else
    	cout<<0;
}

int main()
{
	double x1=0, y1=0;
	double x2=2, y2=2;
	double px=2.1, py=1;
	
	IsPointInBoundingBox(x1,y1,x2,y2,px,py);
	return 0;
}

[BeKinect]

@Half Life G Man: I tested your function, but still no intersection points are found. I now tested a second version of my IsPointInBoundingBox function, wich kind of worked.

  int IsPointInBoundingBox(double x1, double y1, double x2, double y2, double px, double py)
{
  if( (px+0.1) >= min(x1,x2) && (px-0.1) <= max(x1,x2) && (py+0.1) >= min(y1,y2) && (py-0.1) <= max(y1,y2) )
  {
      return 1;
   }
   else
     return 0;
}

This version of the function finds points depending on the bounderies. But I would have to do this for every set of data points manually. What I would need is a automated version.

[BeKinect]

@brown ogum: I think the function is to strict. In my answer to Half Life G Man I tested a version with boundaries and found intersection points, but in this way I also found points in an area where I do not need them. So i would have to adjust the boundaries every time.

[Duthomhas]

I'm not really sure what your matrix is or how it is used, but I'm not going to care. (Sorry.)

Here's some code to compute the intersection of two lines given by AB and CD:

#include <cmath>

struct point
{
  double x, y;
  point( double x = 0, double y = 0 ): x(x), y(y) { }
};

struct intersection: public point
{
  bool parallel;
  bool collinear;
  enum extends_type { none, AB, BA, CD, DC } extends;
  
  intersection( point A, point B, point C, point D, double epsilon = 0.000000001 ): 
    parallel(false),
    collinear(false),
    extends(none)
  {
    // Uses the determinant of the two lines. For more information, refer to one of the following:
    // https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection#Given_two_points_on_each_line
    // http://www.faqs.org/faqs/graphics/algorithms-faq/ (Subject 1.03)
    
    double d = denominator( A, B, C, D );
    
    if (std::abs( d ) < epsilon)
    {
      parallel  = true;
      collinear = abs(numerator( A, C, C, D )) < epsilon;
      return;
    }
    
    double r = numerator( A, C, C, D ) / d;
    double s = numerator( A, C, A, B ) / d;

    x = A.x + r * (B.x - A.x);
    y = A.y + r * (B.y - A.y);
    
    extends = (extends_type)( (r>1) + 2*(r<0) + 3*(s>1) + 4*(s<0) );
  }
  
  private:
    inline double numerator  ( point A, point C, point E, point F ) { return (A.y - C.y) * (F.x - E.x) - (A.x - C.x) * (F.y - E.y); }
    inline double denominator( point A, point B, point C, point D ) { return (B.x - A.x) * (D.y - C.y) - (B.y - A.y) * (D.x - C.x); }
};

You can use it conveniently enough:

point A{ ... };
point B{ ... };
...
auto p = intersection( A, B, C, D );

if (!p.parallel and !p.extends)
  cout << "intersection at (" << p.x << "," << p.y << ")\n";

Hope this helps.

[BeKinect]

@Duoas

I'm not really sure what your matrix is or how it is used, but I'm not going to care. (Sorry.)

Matrix is just a simple vector, which contains the coordinates of a 3d contour.
std::vector<std::array<double, 3>> matrix;
It's not the best solution, but I got it from a former project and I'm not allowed to change it ("Never change a running system"), so I have to deal with it.

To your code. It looks good. The only thing i do not understand is Line 39
extends = (extends_type)( (r>1) + 2*(r<0) + 3*(s>1) + 4*(s<0) );

As I get it, this line tells us if the point is not between the segments, where the point is found. Since I'm just searching for points inside the segments I could ignore this line. Is this correct?

[Duthomhas]

Matrix is just a simple vector, which contains the coordinates of a 3d contour.

Ah. The code I posted above only works over 2D line segments...

Are you looking for a 3D intersection algorithm?

The only thing i do not understand is Line 39

This is an integer-to-enum C-style cast.

As I get it, this line tells us if the point is not between the segments, where the point is found.

Yes, it tells you off of which end of which line segment the point is found.

For example, if 'extends' is 'AB', then the point is on the ray A→B (past B).

Since I'm just searching for points inside the segments I could ignore this line. Is this correct?

Yes.

Be aware that collinear lines may intersect, but my code does not calculate the intersection, as it may involve more than one point. You'll need to handle that case specially if it matters.

Hope this helps.

[BeKinect]

Ah. The code I posted above only works over 2D line segments...
Are you looking for a 3D intersection algorithm?

I store 3d informations, but since I'm working with layers, I just have to deal with 2D intersection. So first I search for intersection in one layer, than in the next and so on..

[48T7M4Gy]

Another method using determinants.
See: http://mathworld.wolfram.com/Line-LineIntersection.html

The article shows how 3D intersections are calculated - somewhat more complex.

#include <iostream>

double det( const double, const double, const double, const double );

int main()
{
    double x = 0;
    double y = 0;
    
    double x1, y1, x2, y2, x3, y3, x4, y4;
    
    x1 = 15; y1 = 10;
    x2 = 15; y2 = 20;
    x3 = 7;  y3 = 28;
    x4 = 45; y4 = 9;
    
    double denominator = det(x1 - x2, y1 - y2, x3 - x4, y3 - y4);
    
    if(denominator != 0)
    {
        x = det( det(x1, y1, x2, y2), x1 - x2, det(x3, y3, x4, y4), x3 - x4 ) / denominator;
        y = det( det(x1, y1, x2, y2), y1 - y2, det(x3, y3, x4, y4), y3 - y4 ) / denominator;

        std::cout << "Intersection X: " << x << "   Y: " << y << '\n';
    }
    else
        std::cout << "Lines don't intersect.\n";
    
   return 0;
}

double det( const double a, const double b, const double c, const double d)
{
    return a * d - b * c;
}

Intersection X: 15   Y: 24
Program ended with exit code: 0

[Duthomhas]

It is the same method (using determinants).

[48T7M4Gy]

Probably unwittingly, but the comment about similarity is profoundly philosophical.

[48T7M4Gy]

However, there is something endearing about the efficiency in number of lines and clarity of the latter submission of the two. ;)

[48T7M4Gy]

BTW collinearity is a concept suited to points rather than lines. A group of collinear lines, if that could be true, is a single line which (only esotericaly) intersects itself at every point in its length. ;)

[48T7M4Gy]

Of course there is another case which has been overlooked so far. And that is where the intersection is only valid where it occurs between 2 line segments and not their unlimited extension, in which case only one of the two lines need to be considered for the requisite test.

[48T7M4Gy]

How about that fern is only one-way.

[48T7M4Gy]

Whether the requirements apply to extended lines or line segments is obviously not clear. Separately there are additional cases where only one line can be extended or if the directions of extension are limited.