Intersection points of line bisector with rectangle

I've been trying to wrap my head around this the whole day...

Basically, I have the coordinates of two points that will always be inside a rectangle. I also know the position of the corners of the rectangle. Those two entry points are given at runtime.

I need an algorithm to find the 2 points where the bisector line made by the line segment between the given points intersects th开发者_JS百科at rectangle.

Some details:

In the above image, A and B are given by their coordinates: A(x1, y1) and B(x2, y2). Basically, I'll need to find position of C and D. Red X is the center of the AB segment. This point (let's call it center) will have to be on the CD line.

What I've did:

• found the center:

  center.x = (A.x+B.x)/2;
  center.y = (A.y+B.y)/2;
  

• found CD slope:

  AB_slope =  A.y - B.y / A.x - B.x;
  CD_slope = -1/AB_slope;
  

Knowing the center and CD slope gave me the equation of CD and such, I've attempted to find a solution by trying the position of the points on all the 4 borders of the rectangle. However, for some reason it doesn't work: every time I have a solution let's say for C, D is plotted outside or vice-versa.

Here are the equations I'm using:

• knowing x:

  y = (CD_slope * (x - center.x)) + center.y;
  if y > 0 && y < 512: #=> solution found!
  

• knowing y:

  x = (y - center.y + CD_slope*center.x)/CD_slope;
  if x > 0 && x < 512: #=> solution found!
  

From this, I could also end up with another segment (let's say I've found C and I know the center), but geometry failed on me to find the extension of this segment till it intersects the other side of the rectangle.

Updated to include coding snippet

(see comments in main function)

typedef struct { double x; double y; } Point;

Point calculate_center(Point p1, Point p2) {
    Point point;
    point.x = (p1.x+p2.x)/2;
    point.y = (p1.y+p2.y)/2;
    return point;
}

double calculate_pslope(Point p1, Point p2) {
    double dy = p1.y - p2.y;
    double dx = p1.x - p2.x;
    double slope = dy/dx; // this is p1 <-> p2 slope

    return -1/slope;
}

int calculate_y_knowing_x(double pslope, Point center, double x, Point *point) {
    double min= 0.00;
    double max= 512.00;
    double y = (pslope * (x - center.x)) + center.y;

    if(y >= min && y <= max) {
        point->x = corner;
        point->y = y;
        printf("++> found Y for X, point is P(%f, %f)\n", point->x, point->y);
        return 1;
    }
    return 0;
}

int calculate_x_knowing_y(double pslope, Point center, double y, Point *point) {
    double min= 0.00;
    double max= 512.00;
    double x = (y - center.y + pslope*center.x)/pslope;

    if(x >= min && x <= max) {
        point->x = x;
        point->y = y;
        printf("++> found X for Y, point is: P(%f, %f)\n", point->x, point->y);
        return 1;
    }
    return 0;
}

int main(int argc, char **argv) {
    Point A, B;

    // parse argv and define A and B
    // this code is omitted here, let's assume:
    // A.x = 175.00;
    // A.y = 420.00;
    // B.x = 316.00;
    // B.y = 62.00;

    Point C;
    Point D;

    Point center;
    double pslope;

    center = calculate_center(A, B);
    pslope = calculate_pslope(A, B);

    // Here's where the fun happens:
    // I'll need to find the right succession of calls to calculate_*_knowing_* 
    // for 4 cases: x=0, X=512 #=> call calculate_y_knowing_x
    // y=0, y=512 #=> call calculate_x_knowing_y
    // and do this 2 times for both C and D points.
    // Also, if point C is found, point D should not be on the same side (thus C != D)

    // for the given A and B points the succession is:
    calculate_y_knowing_x(pslope, center, 0.00, C);
    calculate_y_knowing_x(pslope, center, 512.00, D);
    // will yield: C(0.00, 144.308659), D(512.00, 345.962291)

    // But if A(350.00, 314.00) and B(106.00, 109.00)
    // the succesion should be:
    // calculate_y_knowing_x(pslope, center, 0.00, C);
    // calculate_x_knowing_y(pslope, center, 512.00, D);
    // to yield C(0.00, 482.875610) and D(405.694672, 0.00)


    return 0;
}

This is C code.

Notes:

• The image was drawn by hand.
• The coordinate system is rotated 90° CCW but should not have an impact on the solution
• I'm looking for an algorithm in C, but I can read other programming languages
• This is a 2D problem

You have the equation for CD (in the form (y - y0) = m(x - x0)) which you can transform into the form y = mx + c. You can also transform it into the form x = (1/m)y - (c/m).

You then simply need to find solutions for when x=0, x=512, y=0, y=512.

The following code should do the trick:

typedef struct { float x; float y; } Point;
typedef struct { Point point[2]; } Line;
typedef struct { Point origin; float width; float height; } Rect;
typedef struct { Point origin; Point direction; } Vector;

Point SolveVectorForX(Vector vector, float x)
{
    Point solution;
    solution.x = x;
    solution.y = vector.origin.y +
        (x - vector.origin.x)*vector.direction.y/vector.direction.x;
    return solution;
}

Point SolveVectorForY(Vector vector, float y)
{
    Point solution;
    solution.x = vector.origin.x +
        (y - vector.origin.y)*vector.direction.x/vector.direction.y;
    solution.y = y;
    return solution;
}

Line FindLineBisectorIntersectionWithRect(Rect rect, Line AB)
{
    Point A = AB.point[0];
    Point B = AB.point[1];
    int pointCount = 0;
    int testEdge = 0;
    Line result;
    Vector CD;

    // CD.origin = midpoint of line AB
    CD.origin.x = (A.x + B.x)/2.0;
    CD.origin.y = (A.y + B.y)/2.0;

    // CD.direction = negative inverse of AB.direction (perpendicular to AB)
    CD.direction.x = (B.y - A.y);
    CD.direction.y = (A.x - B.x);

    // for each edge of the rectangle, check:
    // 1. that an intersection with CD is possible (avoid division by zero)
    // 2. that the intersection point falls within the endpoints of the edge
    // 3. if both check out, use that point as one of the solution points
    while ((++testEdge <= 4) && (pointCount < 2))
    {
        Point point;

        switch (testEdge)
        {
            case 1: // check minimum x edge of rect
                if (CD.direction.x == 0) { continue; }
                point = SolveVectorForX(CD, rect.origin.x);
                if (point.y < rect.origin.y) { continue; }
                if (point.y > (rect.origin.y + rect.height)) { continue; }
                break;

            case 2: // check maximum x edge of rect
                if (CD.direction.x == 0) { continue; }
                point = SolveVectorForX(CD, rect.origin.x + rect.width);
                if (point.y < rect.origin.y) { continue; }
                if (point.y > (rect.origin.y + rect.height)) { continue; }
                break;

            case 3: // check minimum y edge of rect
                if (CD.direction.y == 0) { continue; }
                point = SolveVectorForY(CD, rect.origin.y);
                if (point.x < rect.origin.x) { continue; }
                if (point.x > (rect.origin.x + rect.width)) { continue; }
                break;

            case 4: // check maximum y edge of rect
                if (CD.direction.y == 0) { continue; }
                point = SolveVectorForY(CD, rect.origin.y + rect.height);
                if (point.x < rect.origin.x) { continue; }
                if (point.x > (rect.origin.x + rect.width)) { continue; }
                break;
        };

        // if we made it here, this point is one of the solution points
        result.point[pointCount++] = point;
    }

    // pointCount should always be 2
    assert(pointCount == 2);

    return result;
}

We start from the center point C and the direction of AB, D:

C.x = (A.x+B.x) / 2
C.y = (A.y+B.y) / 2
D.x = (A.x-B.x) / 2
D.y = (A.y-B.y) / 2

then if P is a point on the line, CP must be perpendicular to D. The equation of the line is:

DotProduct(P-C, D) = 0

or

CD = C.x*D.x + C.y*D.y
P.x * D.x + P.y * D.y - CD = 0

for each of the four edges of the square, we have an equation:

P.x=0 -> P.y = CD / D.y
P.y=0 -> P.x = CD / D.x
P.x=512 -> P.y = (CD - 512*D.x) / D.y
P.y=512 -> P.x = (CD - 512*D.y) / D.x

Except for degenerate cases where 2 points coincide, only 2 of these 4 points will have both P.x and P.y between 0 and 512. You'll also have to check for the special cases D.x = 0 or D.y = 0.