intersectionpseudocodeline-intersection

How does the line-intersection routine in the Greiner-Hormann algorithm work?


As a part of the Greiner-Hormann algorithm for polygon-clipping (described here), there is this subroutine:

In image form:
The subroutine

And transcribed (attempted):

intersect(P1,P2,Q1,Q2,alphaP,alphaQ)
  WEC_P1 = <P1 - Q1 | (Q2 - Q1)⊥>
  WEC_P2 = <P2 - Q1 | (Q2 - Q1)⊥>
  if (WEC_P1*WEC_P2 <= 0)
    WEC_Q1 = <Q1 - P1 | (P2 - P1)⊥>
    WEC_Q2 = <Q2 - P1 | (P2 - P1)⊥>
      if (WEC_Q1*WEC_Q2 <= 0)
        alphaP = WEC_P1/(WEC_P1 - WEC_P2)
        alphaQ = WEC_Q1/(WEC_Q1 - WEC_Q2)
        return(true); exit
      end if
    end if
  return(false)
end intersect

Which I don't understand at all. The explanation in the document mentions these window edge coordinates, but I don't know what they are and couldn't find any info on them except this slideshow, that doesn't go very much in depth.

Could someone give me a more detailed explanation on what those "wec"s are, and how are they useful in finding the intersection point of two line segments?


Solution

  • My understanding of window edge coordinates, based on the presentation linked in the question, is that they are a signed distance from a point to a line. That is, a scalar whose absolute value is the distance from the point to the line, and is positive if the point and the origin are on the opposite sides of the line, negative if they are on the same side.

    enter image description here

    How you calculate that distance, and why it works, I can't explain it in formal terms, but this video by Casey M. can point you in the right direction. Use it to understand how the dot product relate to the lengths of the two input vectors and come back. You can use this image as a reference:

    enter image description here

    To check if two points are on the opposite sides of a line, just check whether their wecs have the same sign. Do that for both sets of points and line, and you know if the lines are intersecting.

    How do you find those alpha values? Well, for the line defined by p1 and p2, you are basically looking for a percentage that describes how much of the length of p1-p2 you should take to place the new point; and that's exactly what those instructions are doing.

    enter image description here

    Note that in the picture, wec_p1 > 0 and wec_p2 < 0 so (wec_p1 - wec_p2) > 0. Try to convince yourself that those last two assignments are equivalent to that a/b in the picture.

    Here's my version of the pseudocode, with type annotations:

    lines_intersect :: (p1: Vector2, p2: Vector2, q1: Vector2, q2: Vector2) -> bool, float, float {
      intersection_found: bool = false;
      p_alpha: float = 0;
      q_alpha: float = 0;
    
      q_normal: Vector2 = { q2.y - q1.y, q1.x - q2.x };
      p1_wec: float = dot_product(q_normal, p1 - q1);
      p2_wec: float = dot_product(q_normal, p2 - q1);
      
      if p1_wec * p2_wec <= 0 {
        p_normal: Vector2 = { p2.y - p1.y, p1.x - p2.x };
        q1_wec: float = dot_product(p_normal, q1 - p1);
        q2_wec: float = dot_product(p_normal, q2 - p1);
        
        if q1_wec * q2_wec <= 0 {
          intersection_found = true;
          p_alpha = p1_wec / (p1_wec - p2_wec);
          q_alpha = q1_wec / (q1_wec - q2_wec);
        }
      }
      
      return intersection_found, p_alpha, q_alpha;
    }