The   Sutherland-Hodgman clipping algorithm   finds the polygon that is the intersection between an arbitrary polygon (the “subject polygon”) and a convex polygon (the “clip polygon”).

It is used in computer graphics (especially 2D graphics) to reduce the complexity of a scene being displayed by eliminating parts of a polygon that do not need to be displayed.



Task

Take the closed polygon defined by the points:

 



[
(
50
,
150
)
,
(
200
,
50
)
,
(
350
,
150
)
,
(
350
,
300
)
,
(
250
,
300
)
,
(
200
,
250
)
,
(
150
,
350
)
,
(
100
,
250
)
,
(
100
,
200
)
]


{\displaystyle [(50,150),(200,50),(350,150),(350,300),(250,300),(200,250),(150,350),(100,250),(100,200)]}


and clip it by the rectangle defined by the points:

 



[
(
100
,
100
)
,
(
300
,
100
)
,
(
300
,
300
)
,
(
100
,
300
)
]


{\displaystyle [(100,100),(300,100),(300,300),(100,300)]}


Print the sequence of points that define the resulting clipped polygon.



Extra credit

Display all three polygons on a graphical surface, using a different color for each polygon and filling the resulting polygon.

(When displaying you may use either a north-west or a south-west origin, whichever is more convenient for your display mechanism.)




      VDU 23,22,200;200;8,16,16,128
      VDU 23,23,2;0;0;0;
 
      DIM SubjPoly{(8) x, y}
      DIM ClipPoly{(3) x, y}
      FOR v% = 0 TO 8 : READ SubjPoly{(v%)}.x, SubjPoly{(v%)}.y : NEXT
      DATA 50,150,200,50,350,150,350,300,250,300,200,250,150,350,100,250,100,200
      FOR v% = 0 TO 3 : READ ClipPoly{(v%)}.x, ClipPoly{(v%)}.y : NEXT
      DATA 100,100, 300,100, 300,300, 100,300
 
      GCOL 4 : PROCplotpoly(SubjPoly{()}, 9)
      GCOL 1 : PROCplotpoly(ClipPoly{()}, 4)
      nvert% = FNsutherland_hodgman(SubjPoly{()}, ClipPoly{()}, Clipped{()})
      GCOL 2 : PROCplotpoly(Clipped{()}, nvert%)
      END
 
      DEF FNsutherland_hodgman(subj{()}, clip{()}, RETURN out{()})
      LOCAL i%, j%, n%, o%, p1{}, p2{}, s{}, e{}, p{}, inp{()}
      DIM p1{x,y}, p2{x,y}, s{x,y}, e{x,y}, p{x,y}
      n% = DIM(subj{()},1) + DIM(clip{()},1)
      DIM inp{(n%) x, y}, out{(n%) x,y}
      FOR o% = 0 TO DIM(subj{()},1) : out{(o%)} = subj{(o%)} : NEXT
      p1{} = clip{(DIM(clip{()},1))}
      FOR i% = 0 TO DIM(clip{()},1)
        p2{} = clip{(i%)}
        FOR n% = 0 TO o% - 1 : inp{(n%)} = out{(n%)} : NEXT : o% = 0
        IF n% >= 2 THEN
          s{} = inp{(n% - 1)}
          FOR j% = 0 TO n% - 1
            e{} = inp{(j%)}
            IF FNside(e{}, p1{}, p2{}) THEN
              IF NOT FNside(s{}, p1{}, p2{}) THEN
                PROCintersection(p1{}, p2{}, s{}, e{}, p{})
                out{(o%)} = p{}
                o% += 1
              ENDIF
              out{(o%)} = e{}
              o% += 1
            ELSE
              IF FNside(s{}, p1{}, p2{}) THEN
                PROCintersection(p1{}, p2{}, s{}, e{}, p{})
                out{(o%)} = p{}
                o% += 1
              ENDIF
            ENDIF
            s{} = e{}
          NEXT
        ENDIF
        p1{} = p2{}
      NEXT i%
      = o%
 
      REM Which side of the line p1-p2 is the point p?
      DEF FNside(p{}, p1{}, p2{})
      =  (p2.x - p1.x) * (p.y - p1.y) > (p2.y - p1.y) * (p.x - p1.x)
 
      REM Find the intersection of two lines p1-p2 and p3-p4
      DEF PROCintersection(p1{}, p2{}, p3{}, p4{}, p{})
      LOCAL a{}, b{}, k, l, m : DIM a{x,y}, b{x,y}
      a.x = p1.x - p2.x : a.y = p1.y - p2.y
      b.x = p3.x - p4.x : b.y = p3.y - p4.y
      k = p1.x * p2.y - p1.y * p2.x
      l = p3.x * p4.y - p3.y * p4.x
      m = 1 / (a.x * b.y - a.y * b.x)
      p.x =  m * (k * b.x - l * a.x)
      p.y =  m * (k * b.y - l * a.y)
      ENDPROC
 
      REM plot a polygon
      DEF PROCplotpoly(poly{()}, n%)
      LOCAL i%
      MOVE poly{(0)}.x, poly{(0)}.y
      FOR i% = 1 TO n%-1
        DRAW poly{(i%)}.x, poly{(i%)}.y
      NEXT
      DRAW poly{(0)}.x, poly{(0)}.y
      ENDPROC