Sutherland-Hodgman polygon clipping

From Rosetta Code
Task
Sutherland-Hodgman polygon clipping
You are encouraged to solve this task according to the task description, using any language you may know.

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:

and clip it by the rectangle defined by the points:

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.)

11l[edit]

Translation of: Python
F clip(subjectPolygon, clipPolygon)
   F inside(p, cp1, cp2)
      R (cp2.x - cp1.x) * (p.y - cp1.y) > (cp2.y - cp1.y) * (p.x - cp1.x)

   F computeIntersection(s, e, cp1, cp2)
      V dc = cp1 - cp2
      V dp = s - e
      V n1 = cp1.x * cp2.y - cp1.y * cp2.x
      V n2 = s.x * e.y - s.y * e.x
      V n3 = 1.0 / (dc.x * dp.y - dc.y * dp.x)
      R ((n1 * dp.x - n2 * dc.x) * n3, (n1 * dp.y - n2 * dc.y) * n3)

   V outputList = subjectPolygon
   V cp1 = clipPolygon.last

   L(clipVertex) clipPolygon
      V cp2 = clipVertex
      V inputList = outputList
      outputList.clear()
      V s = inputList.last

      L(subjectVertex) inputList
         V e = subjectVertex
         I inside(e, cp1, cp2)
            I !inside(s, cp1, cp2)
               outputList.append(computeIntersection(s, e, cp1, cp2))
            outputList.append(e)
         E I inside(s, cp1, cp2)
            outputList.append(computeIntersection(s, e, cp1, cp2))
         s = e
      cp1 = cp2
   R (outputList)

V subjectp = [(50.0, 150.0), (200.0, 50.0), (350.0, 150.0), (350.0, 300.0), (250.0, 300.0), (200.0, 250.0), (150.0, 350.0), (100.0, 250.0), (100.0, 200.0)]
V clipp = [(100.0, 100.0), (300.0, 100.0), (300.0, 300.0), (100.0, 300.0)]
print_elements(clip(subjectp, clipp), sep' "\n")
Output:
(100, 116.667)
(125, 100)
(275, 100)
(300, 116.667)
(300, 300)
(250, 300)
(200, 250)
(175, 300)
(125, 300)
(100, 250)

Ada[edit]

with Ada.Containers.Doubly_Linked_Lists;
with Ada.Text_IO;

procedure Main is
   package FIO is new Ada.Text_IO.Float_IO (Float);

   type Point is record
      X, Y : Float;
   end record;

   function "-" (Left, Right : Point) return Point is
   begin
      return (Left.X - Right.X, Left.Y - Right.Y);
   end "-";

   type Edge is array (1 .. 2) of Point;

   package Point_Lists is new Ada.Containers.Doubly_Linked_Lists
     (Element_Type => Point);
   use type Point_Lists.List;
   subtype Polygon is Point_Lists.List;

   function Inside (P : Point; E : Edge) return Boolean is
   begin
      return (E (2).X - E (1).X) * (P.Y - E (1).Y) >
             (E (2).Y - E (1).Y) * (P.X - E (1).X);
   end Inside;

   function Intersecton (P1, P2 : Point; E : Edge) return Point is
      DE : Point := E (1) - E (2);
      DP : Point := P1 - P2;
      N1 : Float := E (1).X * E (2).Y - E (1).Y * E (2).X;
      N2 : Float := P1.X * P2.Y - P1.Y * P2.X;
      N3 : Float := 1.0 / (DE.X * DP.Y - DE.Y * DP.X);
   begin
      return ((N1 * DP.X - N2 * DE.X) * N3, (N1 * DP.Y - N2 * DE.Y) * N3);
   end Intersecton;

   function Clip (P, C : Polygon) return Polygon is
      use Point_Lists;
      A, B, S, E : Cursor;
      Inputlist  : List;
      Outputlist : List := P;
      AB         : Edge;
   begin
      A := C.First;
      B := C.Last;
      while A /= No_Element loop
         AB        := (Element (B), Element (A));
         Inputlist := Outputlist;
         Outputlist.Clear;
         S := Inputlist.Last;
         E := Inputlist.First;
         while E /= No_Element loop
            if Inside (Element (E), AB) then
               if not Inside (Element (S), AB) then
                  Outputlist.Append
                    (Intersecton (Element (S), Element (E), AB));
               end if;
               Outputlist.Append (Element (E));
            elsif Inside (Element (S), AB) then
               Outputlist.Append
                 (Intersecton (Element (S), Element (E), AB));
            end if;
            S := E;
            E := Next (E);
         end loop;
         B := A;
         A := Next (A);
      end loop;
      return Outputlist;
   end Clip;

   procedure Print (P : Polygon) is
      use Point_Lists;
      C : Cursor := P.First;
   begin
      Ada.Text_IO.Put_Line ("{");
      while C /= No_Element loop
         Ada.Text_IO.Put (" (");
         FIO.Put (Element (C).X, Exp => 0);
         Ada.Text_IO.Put (',');
         FIO.Put (Element (C).Y, Exp => 0);
         Ada.Text_IO.Put (')');
         C := Next (C);
         if C /= No_Element then
            Ada.Text_IO.Put (',');
         end if;
         Ada.Text_IO.New_Line;
      end loop;
      Ada.Text_IO.Put_Line ("}");
   end Print;

   Source  : Polygon;
   Clipper : Polygon;
   Result  : Polygon;
begin
   Source.Append ((50.0, 150.0));
   Source.Append ((200.0, 50.0));
   Source.Append ((350.0, 150.0));
   Source.Append ((350.0, 300.0));
   Source.Append ((250.0, 300.0));
   Source.Append ((200.0, 250.0));
   Source.Append ((150.0, 350.0));
   Source.Append ((100.0, 250.0));
   Source.Append ((100.0, 200.0));
   Clipper.Append ((100.0, 100.0));
   Clipper.Append ((300.0, 100.0));
   Clipper.Append ((300.0, 300.0));
   Clipper.Append ((100.0, 300.0));
   Result := Clip (Source, Clipper);
   Print (Result);
end Main;
Output:
{
 (100.00000,116.66667),
 (125.00000,100.00000),
 (275.00000,100.00000),
 (300.00000,116.66667),
 (300.00000,300.00000),
 (250.00000,300.00000),
 (200.00000,250.00000),
 (175.00000,300.00000),
 (125.00000,300.00000),
 (100.00000,250.00000)
}

ATS[edit]

(*------------------------------------------------------------------*)
(* Sutherland-Hodgman polygon clipping. *)

#include "share/atspre_staload.hats"

#define NIL list_nil ()
#define ::  list_cons

(*------------------------------------------------------------------*)

typedef coordinate = double

fn {tk : tkind}
coord_make_g0int (x : g0int tk)
    :<> coordinate =
  g0i2f x

fn {tk : tkind}
coord_make_g0float (x : g0float tk)
    :<> coordinate =
  g0f2f x

overload coord with coord_make_g0int
overload coord with coord_make_g0float

(*------------------------------------------------------------------*)

datatype point = point of (coordinate, coordinate)
datatype closedpoly = closedpoly of (arrszref point)

fn
fprint_coordinate (outf : FILEref,
                   x    : coordinate)
    : void =
  let
    val _ = $extfcall (int, "fprintf", outf, "%g", x)
  in
  end

fn
fprint_point (outf : FILEref,
              pt   : point) =
  let

    val+ point (x, y) = pt
  in
    fprint! (outf, "(");
    fprint_coordinate (outf, x);
    fprint! (outf, ",");
    fprint_coordinate (outf, y);
    fprint! (outf, ")")
  end

overload fprint with fprint_point

fn
fprint_closedpoly
          (outf : FILEref,
           poly : closedpoly)
    : void =
  let
    val+ closedpoly points = poly
    val n = size points
    var i : size_t
  in
    for (i := i2sz 0; i <> n; i := succ i)
      fprint! (outf, points[i], "---");
    fprint! (outf, "cycle")
  end

fn
print_closedpoly (poly : closedpoly) =
  fprint_closedpoly (stdout_ref, poly)

overload fprint with fprint_closedpoly

fn
closedpoly_make_list
          (points : List point)
    : closedpoly =
  closedpoly (arrszref_make_list<point> points)

(*------------------------------------------------------------------*)

fn
evaluate_line (x1 : coordinate,
               y1 : coordinate,
               x2 : coordinate,
               y2 : coordinate,
               x  : coordinate)
    :<> coordinate =
  let
    val dy = y2 - y1
    and dx = x2 - x1
    val slope = dy / dx
    and intercept = ((dx * y1) - (dy * x1)) / dx
  in
    (slope * x) + intercept
  end

fn
intersection_of_lines
          (x1 : coordinate,
           y1 : coordinate,
           x2 : coordinate,
           y2 : coordinate,
           x3 : coordinate,
           y3 : coordinate,
           x4 : coordinate,
           y4 : coordinate)
    :<> point =
  if x1 = x2 then
    point (x1, evaluate_line (x3, y3, x4, y4, x1))
  else if x3 = x4 then
    point (x3, evaluate_line (x1, y1, x2, y2, x3))
  else
    let
      val denominator =
        ((x1 - x2) * (y3 - y4)) - ((y1 - y2) * (x3 - x4))
      and x1y2_y1x2 = (x1 * y2) - (y1 * x2)
      and x3y4_y3x4 = (x3 * y4) - (y3 * x4)

      val xnumerator = (x1y2_y1x2 * (x3 - x4)) - ((x1 - x2) * x3y4_y3x4)
      and ynumerator = (x1y2_y1x2 * (y3 - y4)) - ((y1 - y2) * x3y4_y3x4)
    in
      point (xnumerator / denominator,
             ynumerator / denominator)
    end

fn
intersection_of_edges
          (e1 : @(point, point),
           e2 : @(point, point))
    :<> point =
  let
    val+ @(point (x1, y1), point (x2, y2)) = e1
    and  @(point (x3, y3), point (x4, y4)) = e2
  in
    intersection_of_lines (x1, y1, x2, y2, x3, y3, x4, y4)
  end

fn
point_is_left_of_edge
          (pt   : point,
           edge : @(point, point))
    :<> bool =
  let
    val+ point (x, y) = pt
    and  @(point (x1, y1), point (x2, y2)) = edge
  in
    (* Outer product of the vectors (x1,y1)-->(x,y) and
       (x1,y1)-->(x2,y2). *)
    ((x - x1) * (y2 - y1)) - ((x2 - x1) * (y - y1)) < coord 0
  end

fn
clip_subject_edge
          (subject_edge : @(point, point),
           clip_edge    : @(point, point),
           accum        : List0 point)
    : List0 point =
  let
    macdef left_of = point_is_left_of_edge
    macdef intersection =
      intersection_of_edges (subject_edge, clip_edge)

    val @(s1, s2) = subject_edge
    val s2_is_inside = s2 \left_of clip_edge
    and s1_is_inside = s1 \left_of clip_edge
  in
    case+ (s2_is_inside, s1_is_inside) of
    | (true, true) => s2 :: accum
    | (true, false) => s2 :: intersection :: accum
    | (false, true) => intersection :: accum
    | (false, false) => accum
  end

fun
for_each_subject_edge
          (i              : size_t,
           subject_points : arrszref point,
           clip_edge      : @(point, point),
           accum          : List0 point)
    : arrszref point =
  let
    val n = size subject_points
  in
    if i = n then
      arrszref_make_rlist accum
    else
      let
        val s2 = subject_points[i]
        and s1 =
          begin
            if i = 0 then
              subject_points[pred n]
            else
              subject_points[pred i]
          end
        val accum = clip_subject_edge (@(s1, s2), clip_edge, accum)
      in
        for_each_subject_edge (succ i, subject_points, clip_edge,
                               accum)
      end
  end

fun
for_each_clip_edge
          (i              : size_t,
           subject_points : arrszref point,
           clip_points    : arrszref point)
    : arrszref point =
  let
    val n = size clip_points
  in
    if i = n then
      subject_points
    else
      let
        val c2 = clip_points[i]
        and c1 =
          begin
            if i = 0 then
              clip_points[pred n]
            else
              clip_points[pred i]
          end

        val subject_points =
          for_each_subject_edge
            (i2sz 0, subject_points, @(c1, c2), NIL)
      in
        for_each_clip_edge (succ i, subject_points, clip_points)
      end
  end

fn
clip_closedpoly_closedpoly
          (subject_poly : closedpoly,
           clip_poly    : closedpoly)
    : closedpoly =
  let
    val+ closedpoly subject_points = subject_poly
    and  closedpoly clip_points = clip_poly
    val result_points =
      for_each_clip_edge (i2sz 0, subject_points, clip_points)
  in
    closedpoly result_points
  end

overload clip with clip_closedpoly_closedpoly

(*------------------------------------------------------------------*)
(* A function to create an EPS file. *)

(* The EPS code is based on that which is generated by the C
   implementation of this task. *)

fn
write_eps (outf         : FILEref,
           subject_poly : closedpoly,
           clip_poly    : closedpoly,
           result_poly  : closedpoly)
    : void =
  let
    fn
    moveto (pt : point)
        : void =
      let
        val+ point (x, y) = pt
      in
        fprintln! (outf, x, " ", y, " moveto");
      end

    fn
    lineto (pt : point)
        : void =
      let
        val+ point (x, y) = pt
      in
        fprintln! (outf, x, " ", y, " lineto");
      end

    fn
    setrgbcolor (rgb : string)
        : void =
      fprintln! (outf, rgb, " setrgbcolor")

    fn closepath () : void = fprintln! (outf, "closepath")
    fn fill () : void = fprintln! (outf, "fill")
    fn stroke () : void = fprintln! (outf, "stroke")
    fn gsave () : void = fprintln! (outf, "gsave")
    fn grestore () : void = fprintln! (outf, "grestore")

    fn
    showpoly (poly       : closedpoly,
              line_color : string,
              fill_color : string)
        : void =
      let
        val+ closedpoly p = poly
        val n = size p

        var i : size_t
      in
        moveto p[0];
        for (i := i2sz 1; i <> n; i := succ i)
          lineto p[i];
        closepath ();
        setrgbcolor line_color;
        gsave ();
        setrgbcolor fill_color;
        fill ();
        grestore ();
        stroke ()
      end
  in
    fprintln! (outf, "%!PS-Adobe-3.0 EPSF-3.0");
    fprintln! (outf, "%%BoundingBox: 40 40 360 360");
    fprintln! (outf, "0 setlinewidth ");
    showpoly (clip_poly, ".5 0 0", "1 .7 .7");
    showpoly (subject_poly, "0 .2 .5", ".4 .7 1");
    fprintln! (outf, "2 setlinewidth");
    fprintln! (outf, "[10 8] 0 setdash");
    showpoly (result_poly, ".5 0 .5", ".7 .3 .8");
    fprintln! (outf, "%%EOF")
  end

fn
write_eps_to_file
          (outfile      : string,
           subject_poly : closedpoly,
           clip_poly    : closedpoly,
           result_poly  : closedpoly)
    : void =
  let
    val outf = fileref_open_exn (outfile, file_mode_w)
  in
    write_eps (outf, subject_poly, clip_poly, result_poly);
    fileref_close outf
  end

(*------------------------------------------------------------------*)

implement
main0 () =
  let
    val outf = stdout_ref
    
    val subject_poly =
      closedpoly_make_list
        $list (point (coord 50, coord 150),
               point (coord 200, coord 50),
               point (coord 350, coord 150),
               point (coord 350, coord 300),
               point (coord 250, coord 300),
               point (coord 200, coord 250),
               point (coord 150, coord 350),
               point (coord 100, coord 250),
               point (coord 100, coord 200))
    val clip_poly =
      closedpoly_make_list
        $list (point (coord 100, coord 100),
               point (coord 300, coord 100),
               point (coord 300, coord 300),
               point (coord 100, coord 300))

    val result_poly = clip (subject_poly, clip_poly)
  in
    fprintln! (outf, result_poly);
    write_eps_to_file ("sutherland-hodgman.eps",
                       subject_poly, clip_poly, result_poly);
    fprintln! (outf, "Wrote sutherland-hodgman.eps")
  end

(*------------------------------------------------------------------*)
Output:
$ patscc -O3 -DATS_MEMALLOC_GCBDW sutherland-hodgman.dats -lgc && ./a.out
(100,116.667)---(125,100)---(275,100)---(300,116.667)---(300,300)---(250,300)---(200,250)---(175,300)---(125,300)---(100,250)---cycle
Wrote sutherland-hodgman.eps

The polygons of the task problem.

Here is a simple fixed-point version of the same program:

(*------------------------------------------------------------------*)
(* Sutherland-Hodgman polygon clipping (fixed-point version). *)

#include "share/atspre_staload.hats"

#define NIL list_nil ()
#define ::  list_cons

implement g0int2float<intknd,ldblknd> x = $UNSAFE.cast x
implement g0int2float<llintknd,ldblknd> x = $UNSAFE.cast x

(*------------------------------------------------------------------*)

abst@ype coordinate = llint

extern castfn llint2coordinate : llint -<> coordinate
extern castfn coordinate2llint : coordinate -<> llint

overload lli2coord with llint2coordinate
overload coord2lli with coordinate2llint

#define SCALE 262144            (* 18 fraction bits. *)

fn {tk : tkind}
coordinate_make_g0int (x : g0int tk)
    :<> coordinate =
  let
    val x : llint = g0i2i x
  in
    llint2coordinate (x * g0i2i SCALE)
  end

fn
add_coordinate (x : coordinate, y : coordinate)
    :<> coordinate =
  lli2coord (coord2lli x + coord2lli y)

fn
sub_coordinate (x : coordinate, y : coordinate)
    :<> coordinate =
  lli2coord (coord2lli x - coord2lli y)

fn
mul_coordinate (x : coordinate, y : coordinate)
    :<> coordinate =
  lli2coord ((coord2lli x * coord2lli y) / g0i2i SCALE)

fn
div_coordinate (x : coordinate, y : coordinate)
    :<> coordinate =
  lli2coord ((coord2lli x * g0i2i SCALE) / coord2lli y)

fn
eq_coordinate (x : coordinate, y : coordinate)
    :<> bool =
  coord2lli x = coord2lli y

fn
lt_coordinate (x : coordinate, y : coordinate)
    :<> bool =
  coord2lli x < coord2lli y

overload coord with coordinate_make_g0int
overload + with add_coordinate
overload - with sub_coordinate
overload * with mul_coordinate
overload / with div_coordinate
overload = with eq_coordinate
overload < with lt_coordinate

fn
fprint_coordinate (outf : FILEref,
                   x    : coordinate)
    : void =
  let
    val x : ldouble = g0i2f (coord2lli x)
    val x = x / g0i2f SCALE
    val _ = $extfcall (int, "fprintf", outf, "%Lg", x)
  in
  end

(*------------------------------------------------------------------*)

datatype point = point of (coordinate, coordinate)
datatype closedpoly = closedpoly of (arrszref point)

fn
fprint_point (outf : FILEref,
              pt   : point) =
  let

    val+ point (x, y) = pt
  in
    fprint! (outf, "(");
    fprint_coordinate (outf, x);
    fprint! (outf, ",");
    fprint_coordinate (outf, y);
    fprint! (outf, ")")
  end

overload fprint with fprint_point

fn
fprint_closedpoly
          (outf : FILEref,
           poly : closedpoly)
    : void =
  let
    val+ closedpoly points = poly
    val n = size points
    var i : size_t
  in
    for (i := i2sz 0; i <> n; i := succ i)
      fprint! (outf, points[i], "---");
    fprint! (outf, "cycle")
  end

fn
print_closedpoly (poly : closedpoly) =
  fprint_closedpoly (stdout_ref, poly)

overload fprint with fprint_closedpoly

fn
closedpoly_make_list
          (points : List point)
    : closedpoly =
  closedpoly (arrszref_make_list<point> points)

(*------------------------------------------------------------------*)

fn
evaluate_line (x1 : coordinate,
               y1 : coordinate,
               x2 : coordinate,
               y2 : coordinate,
               x  : coordinate)
    :<> coordinate =
  let
    val dy = y2 - y1
    and dx = x2 - x1
    val slope = dy / dx
    and intercept = ((dx * y1) - (dy * x1)) / dx
  in
    (slope * x) + intercept
  end

fn
intersection_of_lines
          (x1 : coordinate,
           y1 : coordinate,
           x2 : coordinate,
           y2 : coordinate,
           x3 : coordinate,
           y3 : coordinate,
           x4 : coordinate,
           y4 : coordinate)
    :<> point =
  if x1 = x2 then
    point (x1, evaluate_line (x3, y3, x4, y4, x1))
  else if x3 = x4 then
    point (x3, evaluate_line (x1, y1, x2, y2, x3))
  else
    let
      val denominator =
        ((x1 - x2) * (y3 - y4)) - ((y1 - y2) * (x3 - x4))
      and x1y2_y1x2 = (x1 * y2) - (y1 * x2)
      and x3y4_y3x4 = (x3 * y4) - (y3 * x4)

      val xnumerator = (x1y2_y1x2 * (x3 - x4)) - ((x1 - x2) * x3y4_y3x4)
      and ynumerator = (x1y2_y1x2 * (y3 - y4)) - ((y1 - y2) * x3y4_y3x4)
    in
      point (xnumerator / denominator,
             ynumerator / denominator)
    end

fn
intersection_of_edges
          (e1 : @(point, point),
           e2 : @(point, point))
    :<> point =
  let
    val+ @(point (x1, y1), point (x2, y2)) = e1
    and  @(point (x3, y3), point (x4, y4)) = e2
  in
    intersection_of_lines (x1, y1, x2, y2, x3, y3, x4, y4)
  end

fn
point_is_left_of_edge
          (pt   : point,
           edge : @(point, point))
    :<> bool =
  let
    val+ point (x, y) = pt
    and  @(point (x1, y1), point (x2, y2)) = edge
  in
    (* Outer product of the vectors (x1,y1)-->(x,y) and
       (x1,y1)-->(x2,y2). *)
    ((x - x1) * (y2 - y1)) - ((x2 - x1) * (y - y1)) < coord 0
  end

fn
clip_subject_edge
          (subject_edge : @(point, point),
           clip_edge    : @(point, point),
           accum        : List0 point)
    : List0 point =
  let
    macdef left_of = point_is_left_of_edge
    macdef intersection =
      intersection_of_edges (subject_edge, clip_edge)

    val @(s1, s2) = subject_edge
    val s2_is_inside = s2 \left_of clip_edge
    and s1_is_inside = s1 \left_of clip_edge
  in
    case+ (s2_is_inside, s1_is_inside) of
    | (true, true) => s2 :: accum
    | (true, false) => s2 :: intersection :: accum
    | (false, true) => intersection :: accum
    | (false, false) => accum
  end

fun
for_each_subject_edge
          (i              : size_t,
           subject_points : arrszref point,
           clip_edge      : @(point, point),
           accum          : List0 point)
    : arrszref point =
  let
    val n = size subject_points
  in
    if i = n then
      arrszref_make_rlist accum
    else
      let
        val s2 = subject_points[i]
        and s1 =
          begin
            if i = 0 then
              subject_points[pred n]
            else
              subject_points[pred i]
          end
        val accum = clip_subject_edge (@(s1, s2), clip_edge, accum)
      in
        for_each_subject_edge (succ i, subject_points, clip_edge,
                               accum)
      end
  end

fun
for_each_clip_edge
          (i              : size_t,
           subject_points : arrszref point,
           clip_points    : arrszref point)
    : arrszref point =
  let
    val n = size clip_points
  in
    if i = n then
      subject_points
    else
      let
        val c2 = clip_points[i]
        and c1 =
          begin
            if i = 0 then
              clip_points[pred n]
            else
              clip_points[pred i]
          end

        val subject_points =
          for_each_subject_edge
            (i2sz 0, subject_points, @(c1, c2), NIL)
      in
        for_each_clip_edge (succ i, subject_points, clip_points)
      end
  end

fn
clip_closedpoly_closedpoly
          (subject_poly : closedpoly,
           clip_poly    : closedpoly)
    : closedpoly =
  let
    val+ closedpoly subject_points = subject_poly
    and  closedpoly clip_points = clip_poly
    val result_points =
      for_each_clip_edge (i2sz 0, subject_points, clip_points)
  in
    closedpoly result_points
  end

overload clip with clip_closedpoly_closedpoly

(*------------------------------------------------------------------*)
(* A function to create an EPS file. *)

(* The EPS code is based on that which is generated by the C
   implementation of this task. *)

fn
write_eps (outf         : FILEref,
           subject_poly : closedpoly,
           clip_poly    : closedpoly,
           result_poly  : closedpoly)
    : void =
  let
    fn
    moveto (pt : point)
        : void =
      let
        val+ point (x, y) = pt
      in
        fprint_coordinate (outf, x);
        fprint! (outf, " ");
        fprint_coordinate (outf, y);
        fprintln! (outf, " moveto")
      end

    fn
    lineto (pt : point)
        : void =
      let
        val+ point (x, y) = pt
      in
        fprint_coordinate (outf, x);
        fprint! (outf, " ");
        fprint_coordinate (outf, y);
        fprintln! (outf, " lineto")
      end

    fn
    setrgbcolor (rgb : string)
        : void =
      fprintln! (outf, rgb, " setrgbcolor")

    fn closepath () : void = fprintln! (outf, "closepath")
    fn fill () : void = fprintln! (outf, "fill")
    fn stroke () : void = fprintln! (outf, "stroke")
    fn gsave () : void = fprintln! (outf, "gsave")
    fn grestore () : void = fprintln! (outf, "grestore")

    fn
    showpoly (poly       : closedpoly,
              line_color : string,
              fill_color : string)
        : void =
      let
        val+ closedpoly p = poly
        val n = size p

        var i : size_t
      in
        moveto p[0];
        for (i := i2sz 1; i <> n; i := succ i)
          lineto p[i];
        closepath ();
        setrgbcolor line_color;
        gsave ();
        setrgbcolor fill_color;
        fill ();
        grestore ();
        stroke ()
      end
  in
    fprintln! (outf, "%!PS-Adobe-3.0 EPSF-3.0");
    fprintln! (outf, "%%BoundingBox: 40 40 360 360");
    fprintln! (outf, "0 setlinewidth ");
    showpoly (clip_poly, ".5 0 0", "1 .7 .7");
    showpoly (subject_poly, "0 .2 .5", ".4 .7 1");
    fprintln! (outf, "2 setlinewidth");
    fprintln! (outf, "[10 8] 0 setdash");
    showpoly (result_poly, ".5 0 .5", ".7 .3 .8");
    fprintln! (outf, "%%EOF")
  end

fn
write_eps_to_file
          (outfile      : string,
           subject_poly : closedpoly,
           clip_poly    : closedpoly,
           result_poly  : closedpoly)
    : void =
  let
    val outf = fileref_open_exn (outfile, file_mode_w)
  in
    write_eps (outf, subject_poly, clip_poly, result_poly);
    fileref_close outf
  end

(*------------------------------------------------------------------*)

implement
main0 () =
  let
    val outf = stdout_ref
    
    val subject_poly =
      closedpoly_make_list
        $list (point (coord 50, coord 150),
               point (coord 200, coord 50),
               point (coord 350, coord 150),
               point (coord 350, coord 300),
               point (coord 250, coord 300),
               point (coord 200, coord 250),
               point (coord 150, coord 350),
               point (coord 100, coord 250),
               point (coord 100, coord 200))
    val clip_poly =
      closedpoly_make_list
        $list (point (coord 100, coord 100),
               point (coord 300, coord 100),
               point (coord 300, coord 300),
               point (coord 100, coord 300))

    val result_poly = clip (subject_poly, clip_poly)
  in
    fprintln! (outf, result_poly);
    write_eps_to_file ("sutherland-hodgman.eps",
                       subject_poly, clip_poly, result_poly);
    fprintln! (outf, "Wrote sutherland-hodgman.eps")
  end

(*------------------------------------------------------------------*)
Output:
(100,116.667)---(125,100)---(275,100)---(300,116.666)---(300,300)---(250,300)---(200,250)---(175,300)---(125,300)---(100,250)---cycle
Wrote sutherland-hodgman.eps

BBC BASIC[edit]

      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

Suthhodg bbc.gif

C[edit]

Most of the code is actually storage util routines, such is C. Prints out nodes, and writes test.eps file in current dir.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

typedef struct { double x, y; } vec_t, *vec;

inline double dot(vec a, vec b)
{
	return a->x * b->x + a->y * b->y;
}

inline double cross(vec a, vec b)
{
	return a->x * b->y - a->y * b->x;
}

inline vec vsub(vec a, vec b, vec res)
{
	res->x = a->x - b->x;
	res->y = a->y - b->y;
	return res;
}

/* tells if vec c lies on the left side of directed edge a->b
 * 1 if left, -1 if right, 0 if colinear
 */
int left_of(vec a, vec b, vec c)
{
	vec_t tmp1, tmp2;
	double x;
	vsub(b, a, &tmp1);
	vsub(c, b, &tmp2);
	x = cross(&tmp1, &tmp2);
	return x < 0 ? -1 : x > 0;
}

int line_sect(vec x0, vec x1, vec y0, vec y1, vec res)
{
	vec_t dx, dy, d;
	vsub(x1, x0, &dx);
	vsub(y1, y0, &dy);
	vsub(x0, y0, &d);
	/* x0 + a dx = y0 + b dy ->
	   x0 X dx = y0 X dx + b dy X dx ->
	   b = (x0 - y0) X dx / (dy X dx) */
	double dyx = cross(&dy, &dx);
	if (!dyx) return 0;
	dyx = cross(&d, &dx) / dyx;
	if (dyx <= 0 || dyx >= 1) return 0;

	res->x = y0->x + dyx * dy.x;
	res->y = y0->y + dyx * dy.y;
	return 1;
}

/* === polygon stuff === */
typedef struct { int len, alloc; vec v; } poly_t, *poly;

poly poly_new()
{
	return (poly)calloc(1, sizeof(poly_t));
}

void poly_free(poly p)
{
	free(p->v);
	free(p);
}

void poly_append(poly p, vec v)
{
	if (p->len >= p->alloc) {
		p->alloc *= 2;
		if (!p->alloc) p->alloc = 4;
		p->v = (vec)realloc(p->v, sizeof(vec_t) * p->alloc);
	}
	p->v[p->len++] = *v;
}

/* this works only if all of the following are true:
 *   1. poly has no colinear edges;
 *   2. poly has no duplicate vertices;
 *   3. poly has at least three vertices;
 *   4. poly is convex (implying 3).
*/
int poly_winding(poly p)
{
	return left_of(p->v, p->v + 1, p->v + 2);
}

void poly_edge_clip(poly sub, vec x0, vec x1, int left, poly res)
{
	int i, side0, side1;
	vec_t tmp;
	vec v0 = sub->v + sub->len - 1, v1;
	res->len = 0;

	side0 = left_of(x0, x1, v0);
	if (side0 != -left) poly_append(res, v0);

	for (i = 0; i < sub->len; i++) {
		v1 = sub->v + i;
		side1 = left_of(x0, x1, v1);
		if (side0 + side1 == 0 && side0)
			/* last point and current straddle the edge */
			if (line_sect(x0, x1, v0, v1, &tmp))
				poly_append(res, &tmp);
		if (i == sub->len - 1) break;
		if (side1 != -left) poly_append(res, v1);
		v0 = v1;
		side0 = side1;
	}
}

poly poly_clip(poly sub, poly clip)
{
	int i;
	poly p1 = poly_new(), p2 = poly_new(), tmp;

	int dir = poly_winding(clip);
	poly_edge_clip(sub, clip->v + clip->len - 1, clip->v, dir, p2);
	for (i = 0; i < clip->len - 1; i++) {
		tmp = p2; p2 = p1; p1 = tmp;
		if(p1->len == 0) {
			p2->len = 0;
			break;
		}
		poly_edge_clip(p1, clip->v + i, clip->v + i + 1, dir, p2);
	}

	poly_free(p1);
	return p2;
}

int main()
{
	int i;
	vec_t c[] = {{100,100}, {300,100}, {300,300}, {100,300}};
	//vec_t c[] = {{100,300}, {300,300}, {300,100}, {100,100}};
	vec_t s[] = {	{50,150}, {200,50}, {350,150},
			{350,300},{250,300},{200,250},
			{150,350},{100,250},{100,200}};
#define clen (sizeof(c)/sizeof(vec_t))
#define slen (sizeof(s)/sizeof(vec_t))
	poly_t clipper = {clen, 0, c};
	poly_t subject = {slen, 0, s};

	poly res = poly_clip(&subject, &clipper);

	for (i = 0; i < res->len; i++)
		printf("%g %g\n", res->v[i].x, res->v[i].y);

	/* long and arduous EPS printout */
	FILE * eps = fopen("test.eps", "w");
	fprintf(eps, "%%!PS-Adobe-3.0\n%%%%BoundingBox: 40 40 360 360\n"
		"/l {lineto} def /m{moveto} def /s{setrgbcolor} def"
		"/c {closepath} def /gs {fill grestore stroke} def\n");
	fprintf(eps, "0 setlinewidth %g %g m ", c[0].x, c[0].y);
	for (i = 1; i < clen; i++)
		fprintf(eps, "%g %g l ", c[i].x, c[i].y);
	fprintf(eps, "c .5 0 0 s gsave 1 .7 .7 s gs\n");

	fprintf(eps, "%g %g m ", s[0].x, s[0].y);
	for (i = 1; i < slen; i++)
		fprintf(eps, "%g %g l ", s[i].x, s[i].y);
	fprintf(eps, "c 0 .2 .5 s gsave .4 .7 1 s gs\n");

	fprintf(eps, "2 setlinewidth [10 8] 0 setdash %g %g m ",
		res->v[0].x, res->v[0].y);
	for (i = 1; i < res->len; i++)
		fprintf(eps, "%g %g l ", res->v[i].x, res->v[i].y);
	fprintf(eps, "c .5 0 .5 s gsave .7 .3 .8 s gs\n");

	fprintf(eps, "%%%%EOF");
	fclose(eps);
	printf("test.eps written\n");

	return 0;
}
Output:
200 250

175 300 125 300 100 250 100 200 100 116.667 125 100 275 100 300 116.667 300 300 250 300

test.eps written
Poly-clip-C.png

C#[edit]

This was written in .net 4.0 using wpf

Worker class:

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows;

namespace Sutherland
{
    public static class SutherlandHodgman
    {
        #region Class: Edge

        /// <summary>
        /// This represents a line segment
        /// </summary>
        private class Edge
        {
            public Edge(Point from, Point to)
            {
                this.From = from;
                this.To = to;
            }

            public readonly Point From;
            public readonly Point To;
        }

        #endregion

        /// <summary>
        /// This clips the subject polygon against the clip polygon (gets the intersection of the two polygons)
        /// </summary>
        /// <remarks>
        /// Based on the psuedocode from:
        /// http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman
        /// </remarks>
        /// <param name="subjectPoly">Can be concave or convex</param>
        /// <param name="clipPoly">Must be convex</param>
        /// <returns>The intersection of the two polygons (or null)</returns>
        public static Point[] GetIntersectedPolygon(Point[] subjectPoly, Point[] clipPoly)
        {
            if (subjectPoly.Length < 3 || clipPoly.Length < 3)
            {
                throw new ArgumentException(string.Format("The polygons passed in must have at least 3 points: subject={0}, clip={1}", subjectPoly.Length.ToString(), clipPoly.Length.ToString()));
            }

            List<Point> outputList = subjectPoly.ToList();

            //	Make sure it's clockwise
            if (!IsClockwise(subjectPoly))
            {
                outputList.Reverse();
            }

            //	Walk around the clip polygon clockwise
            foreach (Edge clipEdge in IterateEdgesClockwise(clipPoly))
            {
                List<Point> inputList = outputList.ToList();		//	clone it
                outputList.Clear();

                if (inputList.Count == 0)
                {
                    //	Sometimes when the polygons don't intersect, this list goes to zero.  Jump out to avoid an index out of range exception
                    break;
                }

                Point S = inputList[inputList.Count - 1];

                foreach (Point E in inputList)
                {
                    if (IsInside(clipEdge, E))
                    {
                        if (!IsInside(clipEdge, S))
                        {
                            Point? point = GetIntersect(S, E, clipEdge.From, clipEdge.To);
                            if (point == null)
                            {
                                throw new ApplicationException("Line segments don't intersect");		//	may be colinear, or may be a bug
                            }
                            else
                            {
                                outputList.Add(point.Value);
                            }
                        }

                        outputList.Add(E);
                    }
                    else if (IsInside(clipEdge, S))
                    {
                        Point? point = GetIntersect(S, E, clipEdge.From, clipEdge.To);
                        if (point == null)
                        {
                            throw new ApplicationException("Line segments don't intersect");		//	may be colinear, or may be a bug
                        }
                        else
                        {
                            outputList.Add(point.Value);
                        }
                    }

                    S = E;
                }
            }

            //	Exit Function
            return outputList.ToArray();
        }

        #region Private Methods

        /// <summary>
        /// This iterates through the edges of the polygon, always clockwise
        /// </summary>
        private static IEnumerable<Edge> IterateEdgesClockwise(Point[] polygon)
        {
            if (IsClockwise(polygon))
            {
                #region Already clockwise

                for (int cntr = 0; cntr < polygon.Length - 1; cntr++)
                {
                    yield return new Edge(polygon[cntr], polygon[cntr + 1]);
                }

                yield return new Edge(polygon[polygon.Length - 1], polygon[0]);

                #endregion
            }
            else
            {
                #region Reverse

                for (int cntr = polygon.Length - 1; cntr > 0; cntr--)
                {
                    yield return new Edge(polygon[cntr], polygon[cntr - 1]);
                }

                yield return new Edge(polygon[0], polygon[polygon.Length - 1]);

                #endregion
            }
        }

        /// <summary>
        /// Returns the intersection of the two lines (line segments are passed in, but they are treated like infinite lines)
        /// </summary>
        /// <remarks>
        /// Got this here:
        /// http://stackoverflow.com/questions/14480124/how-do-i-detect-triangle-and-rectangle-intersection
        /// </remarks>
        private static Point? GetIntersect(Point line1From, Point line1To, Point line2From, Point line2To)
        {
            Vector direction1 = line1To - line1From;
            Vector direction2 = line2To - line2From;
            double dotPerp = (direction1.X * direction2.Y) - (direction1.Y * direction2.X);

            // If it's 0, it means the lines are parallel so have infinite intersection points
            if (IsNearZero(dotPerp))
            {
                return null;
            }

            Vector c = line2From - line1From;
            double t = (c.X * direction2.Y - c.Y * direction2.X) / dotPerp;
            //if (t < 0 || t > 1)
            //{
            //    return null;		//	lies outside the line segment
            //}

            //double u = (c.X * direction1.Y - c.Y * direction1.X) / dotPerp;
            //if (u < 0 || u > 1)
            //{
            //    return null;		//	lies outside the line segment
            //}

            //	Return the intersection point
            return line1From + (t * direction1);
        }

        private static bool IsInside(Edge edge, Point test)
        {
            bool? isLeft = IsLeftOf(edge, test);
            if (isLeft == null)
            {
                //	Colinear points should be considered inside
                return true;
            }

            return !isLeft.Value;
        }
        private static bool IsClockwise(Point[] polygon)
        {
            for (int cntr = 2; cntr < polygon.Length; cntr++)
            {
                bool? isLeft = IsLeftOf(new Edge(polygon[0], polygon[1]), polygon[cntr]);
                if (isLeft != null)		//	some of the points may be colinear.  That's ok as long as the overall is a polygon
                {
                    return !isLeft.Value;
                }
            }

            throw new ArgumentException("All the points in the polygon are colinear");
        }

        /// <summary>
        /// Tells if the test point lies on the left side of the edge line
        /// </summary>
        private static bool? IsLeftOf(Edge edge, Point test)
        {
            Vector tmp1 = edge.To - edge.From;
            Vector tmp2 = test - edge.To;

            double x = (tmp1.X * tmp2.Y) - (tmp1.Y * tmp2.X);		//	dot product of perpendicular?

            if (x < 0)
            {
                return false;
            }
            else if (x > 0)
            {
                return true;
            }
            else
            {
                //	Colinear points;
                return null;
            }
        }

        private static bool IsNearZero(double testValue)
        {
            return Math.Abs(testValue) <= .000000001d;
        }

        #endregion
    }
}

Window code:

<Window x:Class="Sutherland.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        Title="Sutherland Hodgman" Background="#B0B0B0" ResizeMode="CanResizeWithGrip" Width="525" Height="450">
    <Grid Margin="4">
        <Grid.RowDefinitions>
            <RowDefinition Height="1*"/>
            <RowDefinition Height="auto"/>
        </Grid.RowDefinitions>

        <Border Grid.Row="0" CornerRadius="4" BorderBrush="#707070" Background="#FFFFFF" BorderThickness="2">
            <Canvas Name="canvas"/>
        </Border>

        <UniformGrid Grid.Row="1" Rows="1" Margin="0,4,0,0">
            <Button Name="btnTriRect" Content="Triangle - Rectangle" Margin="4,0" Click="btnTriRect_Click"/>
            <Button Name="btnConvex" Content="Concave - Convex" Click="btnConvex_Click"/>
        </UniformGrid>
    </Grid>
</Window>
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows;
using System.Windows.Controls;
using System.Windows.Data;
using System.Windows.Documents;
using System.Windows.Input;
using System.Windows.Media;
using System.Windows.Media.Imaging;
using System.Windows.Navigation;
using System.Windows.Shapes;

namespace Sutherland
{
    public partial class MainWindow : Window
    {
        #region Declaration Section

        private Random _rand = new Random();

        private Brush _subjectBack = new SolidColorBrush(ColorFromHex("30427FCF"));
        private Brush _subjectBorder = new SolidColorBrush(ColorFromHex("427FCF"));
        private Brush _clipBack = new SolidColorBrush(ColorFromHex("30D65151"));
        private Brush _clipBorder = new SolidColorBrush(ColorFromHex("D65151"));
        private Brush _intersectBack = new SolidColorBrush(ColorFromHex("609F18CC"));
        private Brush _intersectBorder = new SolidColorBrush(ColorFromHex("9F18CC"));

        #endregion

        #region Constructor

        public MainWindow()
        {
            InitializeComponent();
        }

        #endregion

        #region Event Listeners

        private void btnTriRect_Click(object sender, RoutedEventArgs e)
        {
            try
            {
                double width = canvas.ActualWidth;
                double height = canvas.ActualHeight;

                Point[] poly1 = new Point[] {
				    new Point(_rand.NextDouble() * width, _rand.NextDouble() * height),
				    new Point(_rand.NextDouble() * width, _rand.NextDouble() * height),
				    new Point(_rand.NextDouble() * width, _rand.NextDouble() * height) };

                Point rectPoint = new Point(_rand.NextDouble() * (width * .75d), _rand.NextDouble() * (height * .75d));		//	don't let it start all the way at the bottom right
                Rect rect = new Rect(
                    rectPoint,
                    new Size(_rand.NextDouble() * (width - rectPoint.X), _rand.NextDouble() * (height - rectPoint.Y)));

                Point[] poly2 = new Point[] { rect.TopLeft, rect.TopRight, rect.BottomRight, rect.BottomLeft };

                Point[] intersect = SutherlandHodgman.GetIntersectedPolygon(poly1, poly2);

                canvas.Children.Clear();
                ShowPolygon(poly1, _subjectBack, _subjectBorder, 1d);
                ShowPolygon(poly2, _clipBack, _clipBorder, 1d);
                ShowPolygon(intersect, _intersectBack, _intersectBorder, 3d);
            }
            catch (Exception ex)
            {
                MessageBox.Show(ex.ToString(), this.Title, MessageBoxButton.OK, MessageBoxImage.Error);
            }
        }
        private void btnConvex_Click(object sender, RoutedEventArgs e)
        {
            try
            {
                Point[] poly1 = new Point[] { new Point(50, 150), new Point(200, 50), new Point(350, 150), new Point(350, 300), new Point(250, 300), new Point(200, 250), new Point(150, 350), new Point(100, 250), new Point(100, 200) };
                Point[] poly2 = new Point[] { new Point(100, 100), new Point(300, 100), new Point(300, 300), new Point(100, 300) };

                Point[] intersect = SutherlandHodgman.GetIntersectedPolygon(poly1, poly2);

                canvas.Children.Clear();
                ShowPolygon(poly1, _subjectBack, _subjectBorder, 1d);
                ShowPolygon(poly2, _clipBack, _clipBorder, 1d);
                ShowPolygon(intersect, _intersectBack, _intersectBorder, 3d);
            }
            catch (Exception ex)
            {
                MessageBox.Show(ex.ToString(), this.Title, MessageBoxButton.OK, MessageBoxImage.Error);
            }
        }

        #endregion

        #region Private Methods

        private void ShowPolygon(Point[] points, Brush background, Brush border, double thickness)
        {
            if (points == null || points.Length == 0)
            {
                return;
            }

            Polygon polygon = new Polygon();
            polygon.Fill = background;
            polygon.Stroke = border;
            polygon.StrokeThickness = thickness;

            foreach (Point point in points)
            {
                polygon.Points.Add(point);
            }

            canvas.Children.Add(polygon);
        }

        /// <summary>
        /// This is just a wrapper to the color converter (why can't they have a method off the color class with all
        /// the others?)
        /// </summary>
        private static Color ColorFromHex(string hexValue)
        {
            if (hexValue.StartsWith("#"))
            {
                return (Color)ColorConverter.ConvertFromString(hexValue);
            }
            else
            {
                return (Color)ColorConverter.ConvertFromString("#" + hexValue);
            }
        }

        #endregion
    }
}

PolyIntersect.png

C++[edit]

#include <iostream>
#include <span>
#include <vector>

struct vec2 {
    float x = 0.0f, y = 0.0f;

    constexpr vec2 operator+(vec2 other) const {
        return vec2{x + other.x, y + other.y};
    }

    constexpr vec2 operator-(vec2 other) const {
        return vec2{x - other.x, y - other.y};
    }
};

constexpr vec2 operator*(vec2 a, float b) { return vec2{a.x * b, a.y * b}; }

constexpr float dot(vec2 a, vec2 b) { return a.x * b.x + a.y * b.y; }

constexpr float cross(vec2 a, vec2 b) { return a.x * b.y - b.x * a.y; }

// check if a point is on the LEFT side of an edge
constexpr bool is_inside(vec2 point, vec2 a, vec2 b) {
    return (cross(a - b, point) + cross(b, a)) < 0.0f;
}

// calculate intersection point
constexpr vec2 intersection(vec2 a1, vec2 a2, vec2 b1, vec2 b2) {
    return ((b1 - b2) * cross(a1, a2) - (a1 - a2) * cross(b1, b2)) *
           (1.0f / cross(a1 - a2, b1 - b2));
}

// Sutherland-Hodgman clipping
std::vector<vec2> suther_land_hodgman(
    std::span<vec2 const> subject_polygon, std::span<vec2 const> clip_polygon) {
    if (clip_polygon.empty() || subject_polygon.empty()) {
        return {};
    }

    std::vector<vec2> ring{subject_polygon.begin(), subject_polygon.end()};

    vec2 p1 = clip_polygon[clip_polygon.size() - 1];

    std::vector<vec2> input;

    for (vec2 p2 : clip_polygon) {
        input.clear();
        input.insert(input.end(), ring.begin(), ring.end());
        vec2 s = input[input.size() - 1];

        ring.clear();

        for (vec2 e : input) {
            if (is_inside(e, p1, p2)) {
                if (!is_inside(s, p1, p2)) {
                    ring.push_back(intersection(p1, p2, s, e));
                }

                ring.push_back(e);
            } else if (is_inside(s, p1, p2)) {
                ring.push_back(intersection(p1, p2, s, e));
            }

            s = e;
        }

        p1 = p2;
    }

    return ring;
}

int main(int argc, char **argv) {
    // subject polygon
    vec2 subject_polygon[] = {{50, 150},  {200, 50},  {350, 150},
                              {350, 300}, {250, 300}, {200, 250},
                              {150, 350}, {100, 250}, {100, 200}};

    // clipping polygon
    vec2 clip_polygon[] = {{100, 100}, {300, 100}, {300, 300}, {100, 300}};

    // apply clipping
    std::vector<vec2> clipped_polygon =
        suther_land_hodgman(subject_polygon, clip_polygon);

    // print clipped polygon points
    std::cout << "Clipped polygon points:" << std::endl;
    for (vec2 p : clipped_polygon) {
        std::cout << "(" << p.x << ", " << p.y << ")" << std::endl;
    }

    return EXIT_SUCCESS;
}
Output:
Clipped polygon points:
(100, 116.667)
(125, 100)
(275, 100)
(300, 116.667)
(300, 300)
(250, 300)
(200, 250)
(175, 300)
(125, 300)
(100, 250)

Common Lisp[edit]

Translation of: Scheme
;;; Sutherland-Hodgman polygon clipping.

(defun evaluate-line (x1 y1 x2 y2 x)
  ;; Given the straight line between (x1,y1) and (x2,y2), evaluate it
  ;; at x.
  (let ((dy (- y2 y1))
        (dx (- x2 x1)))
    (let ((slope (/ dy dx))
          (intercept (/ (- (* dx y1) (* dy x1)) dx)))
      (+ (* slope x) intercept))))

(defun intersection-of-lines (x1 y1 x2 y2 x3 y3 x4 y4)
  ;; Given the line between (x1,y1) and (x2,y2), and the line between
  ;; (x3,y3) and (x4,y4), find their intersection.
  (cond ((= x1 x2) (list x1 (evaluate-line x3 y3 x4 y4 x1)))
        ((= x3 x4) (list x3 (evaluate-line x1 y1 x2 y2 x3)))
        (t (let ((denominator (- (* (- x1 x2) (- y3 y4))
                                 (* (- y1 y2) (- x3 x4))))
                 (x1*y2-y1*x2 (- (* x1 y2) (* y1 x2)))
                 (x3*y4-y3*x4 (- (* x3 y4) (* y3 x4))))
             (let ((xnumerator (- (* x1*y2-y1*x2 (- x3 x4))
                                  (* (- x1 x2) x3*y4-y3*x4)))
                   (ynumerator (- (* x1*y2-y1*x2 (- y3 y4))
                                  (* (- y1 y2) x3*y4-y3*x4))))
               (list (/ xnumerator denominator)
                     (/ ynumerator denominator)))))))

(defun intersection-of-edges (e1 e2)
  ;;
  ;; A point is a list of two coordinates, and an edge is a list of
  ;; two points.
  ;;
  (let ((point1 (car e1))
        (point2 (cadr e1))
        (point3 (car e2))
        (point4 (cadr e2)))
    (let ((x1 (car point1))
          (y1 (cadr point1))
          (x2 (car point2))
          (y2 (cadr point2))
          (x3 (car point3))
          (y3 (cadr point3))
          (x4 (car point4))
          (y4 (cadr point4)))
      (intersection-of-lines x1 y1 x2 y2 x3 y3 x4 y4))))

(defun point-is-left-of-edge-p (pt edge)
  (let ((x (car pt))
        (y (cadr pt))
        (x1 (caar edge))
        (y1 (cadar edge))
        (x2 (caadr edge))
        (y2 (cadadr edge)))
    ;; Outer product of the vectors (x1,y1)-->(x,y) and
    ;; (x1,y1)-->(x2,y2)
    (< (- (* (- x x1) (- y2 y1))
          (* (- x2 x1) (- y y1)))
       0)))

(defun clip-subject-edge (subject-edge clip-edge accum)
  (flet ((intersect ()
           (intersection-of-edges subject-edge clip-edge)))
    (let ((s1 (car subject-edge))
          (s2 (cadr subject-edge)))
      (let ((s2-is-inside (point-is-left-of-edge-p s2 clip-edge))
            (s1-is-inside (point-is-left-of-edge-p s1 clip-edge)))
        (if s2-is-inside
            (if s1-is-inside
                (cons s2 accum)
                (cons s2 (cons (intersect) accum)))
            (if s1-is-inside
                (cons (intersect) accum)
                accum))))))

(defun for-each-subject-edge (i subject-points clip-edge accum)
  (let ((n (length subject-points))
        (accum '()))
    (loop for i from 0 to (1- n)
          do (let ((s2 (aref subject-points i))
                   (s1 (aref subject-points
                             (- (if (zerop i) n i) 1))))
               (setf accum (clip-subject-edge (list s1 s2)
                                              clip-edge accum))))
    (coerce (reverse accum) 'vector)))

(defun for-each-clip-edge (i subject-points clip-points)
  (let ((n (length clip-points)))
    (loop for i from 0 to (1- n)
          do (let ((c2 (aref clip-points i))
                   (c1 (aref clip-points (- (if (zerop i) n i) 1))))
               (setf subject-points
                     (for-each-subject-edge 0 subject-points
                                            (list c1 c2) '()))))
    subject-points))

(defun clip (subject-points clip-points)
  (for-each-clip-edge 0 subject-points clip-points))

(defun write-eps (outf subject-points clip-points result-points)
  (flet ((x (pt) (coerce (car pt) 'float))
         (y (pt) (coerce (cadr pt) 'float))
         (code (s)
           (princ s outf)
           (terpri outf)))
    (flet ((moveto (pt)
             (princ (x pt) outf)
             (princ " " outf)
             (princ (y pt) outf)
             (princ " moveto" outf)
             (terpri outf))
           (lineto (pt)
             (princ (x pt) outf)
             (princ " " outf)
             (princ (y pt) outf)
             (princ " lineto" outf)
             (terpri outf))
           (setrgbcolor (rgb)
             (princ rgb outf)
             (princ " setrgbcolor" outf)
             (terpri outf))
           (closepath () (code "closepath"))
           (fill-it () (code "fill"))
           (stroke () (code "stroke"))
           (gsave () (code "gsave"))
           (grestore () (code "grestore")))
      (flet ((showpoly (poly line-color fill-color)
               (let ((n (length poly)))
                 (moveto (aref poly 0))
                 (loop for i from 1 to (1- n)
                       do (lineto (aref poly i)))
                 (closepath)
                 (setrgbcolor line-color)
                 (gsave)
                 (setrgbcolor fill-color)
                 (fill-it)
                 (grestore)
                 (stroke))))

        (code "%!PS-Adobe-3.0 EPSF-3.0")
        (code "%%BoundingBox: 40 40 360 360")
        (code "0 setlinewidth")
        (showpoly clip-points ".5 0 0" "1 .7 .7")
        (showpoly subject-points "0 .2 .5" ".4 .7 1")
        (code "2 setlinewidth")
        (code "[10 8] 0 setdash")
        (showpoly result-points ".5 0 .5" ".7 .3 .8")
        (code "%%EOF")))))

(defun write-eps-to-file (outfile subject-points clip-points
                          result-points)
  (with-open-file (outf outfile :direction :output
                                :if-exists :supersede
                                :if-does-not-exist :create)
    (write-eps outf subject-points clip-points result-points)))

(defvar subject-points
  #((50 150)
    (200 50)
    (350 150)
    (350 300)
    (250 300)
    (200 250)
    (150 350)
    (100 250)
    (100 200)))

(defvar clip-points
  #((100 100)
    (300 100)
    (300 300)
    (100 300)))

(defvar result-points (clip subject-points clip-points))

(princ result-points)
(terpri)
(write-eps-to-file "sutherland-hodgman.eps"
                   subject-points clip-points result-points)
(princ "Wrote sutherland-hodgman.eps")
(terpri)
Output:
$ clisp sutherland-hodgman.lisp
#((100 350/3) (125 100) (275 100) (300 350/3) (300 300) (250 300) (200 250) (175 300) (125 300) (100 250))
Wrote sutherland-hodgman.eps

The polygons as generated by Common Lisp.

D[edit]

import std.stdio, std.array, std.range, std.typecons, std.algorithm;

struct Vec2 { // To be replaced with Phobos code.
    double x, y;

    Vec2 opBinary(string op="-")(in Vec2 other)
    const pure nothrow @safe @nogc {
        return Vec2(this.x - other.x, this.y - other.y);
    }

    typeof(x) cross(in Vec2 other) const pure nothrow @safe @nogc {
        return this.x * other.y - this.y * other.x;
    }
}

immutable(Vec2)[] clip(in Vec2[] subjectPolygon, in Vec2[] clipPolygon)
pure /*nothrow*/ @safe in {
    assert(subjectPolygon.length > 1);
    assert(clipPolygon.length > 1);
    // Probably clipPolygon needs to be convex and probably
    // its vertices need to be listed in a direction.
} out(result) {
    assert(result.length > 1);
} body {
    alias Edge = Tuple!(Vec2,"p", Vec2,"q");

    static enum isInside = (in Vec2 p, in Edge cle)
    pure nothrow @safe @nogc =>
        (cle.q.x - cle.p.x) * (p.y - cle.p.y) >
        (cle.q.y - cle.p.y) * (p.x - cle.p.x);

    static Vec2 intersection(in Edge se, in Edge cle)
    pure nothrow @safe @nogc {
        immutable dc = cle.p - cle.q;
        immutable dp = se.p - se.q;
        immutable n1 = cle.p.cross(cle.q);
        immutable n2 = se.p.cross(se.q);
        immutable n3 = 1.0 / dc.cross(dp);
        return Vec2((n1 * dp.x - n2 * dc.x) * n3,
                    (n1 * dp.y - n2 * dc.y) * n3);
    }

    // How much slower is this compared to lower-level code?
    static enum edges = (in Vec2[] poly) pure nothrow @safe @nogc =>
        // poly[$ - 1 .. $].chain(poly).zip!Edge(poly);
        poly[$ - 1 .. $].chain(poly).zip(poly).map!Edge;

    immutable(Vec2)[] result = subjectPolygon.idup; // Not nothrow.

    foreach (immutable clipEdge; edges(clipPolygon)) {
        immutable inputList = result;
        result.destroy;
        foreach (immutable inEdge; edges(inputList)) {
            if (isInside(inEdge.q, clipEdge)) {
                if (!isInside(inEdge.p, clipEdge))
                    result ~= intersection(inEdge, clipEdge);
                result ~= inEdge.q;
            } else if (isInside(inEdge.p, clipEdge))
                result ~= intersection(inEdge, clipEdge);
        }
    }

    return result;
}

// Code adapted from the C version.
void saveEPSImage(in string fileName, in Vec2[] subjPoly,
                  in Vec2[] clipPoly, in Vec2[] clipped)
in {
    assert(!fileName.empty);
    assert(subjPoly.length > 1);
    assert(clipPoly.length > 1);
    assert(clipped.length > 1);
} body {
    auto eps = File(fileName, "w");

    // The image bounding box is hard-coded, not computed.
    eps.writeln(
"%%!PS-Adobe-3.0
%%%%BoundingBox: 40 40 360 360
/l {lineto} def
/m {moveto} def
/s {setrgbcolor} def
/c {closepath} def
/gs {fill grestore stroke} def
");

    eps.writef("0 setlinewidth %g %g m ", clipPoly[0].tupleof);
    foreach (immutable cl; clipPoly[1 .. $])
        eps.writef("%g %g l ", cl.tupleof);
    eps.writefln("c 0.5 0 0 s gsave 1 0.7 0.7 s gs");

    eps.writef("%g %g m ", subjPoly[0].tupleof);
    foreach (immutable s; subjPoly[1 .. $])
        eps.writef("%g %g l ", s.tupleof);
    eps.writefln("c 0 0.2 0.5 s gsave 0.4 0.7 1 s gs");

    eps.writef("2 setlinewidth [10 8] 0 setdash %g %g m ",
               clipped[0].tupleof);
    foreach (immutable c; clipped[1 .. $])
        eps.writef("%g %g l ", c.tupleof);
    eps.writefln("c 0.5 0 0.5 s gsave 0.7 0.3 0.8 s gs");

    eps.writefln("%%%%EOF");
    eps.close;
    writeln(fileName, " written.");
}

void main() {
    alias V = Vec2;
    immutable subjectPolygon = [V(50, 150), V(200, 50), V(350, 150),
                                V(350, 300), V(250, 300), V(200, 250),
                                V(150, 350), V(100, 250), V(100, 200)];
    immutable clippingPolygon = [V(100, 100), V(300, 100),
                                 V(300, 300), V(100, 300)];
    immutable clipped = subjectPolygon.clip(clippingPolygon);
    writefln("%(%s\n%)", clipped);
    saveEPSImage("sutherland_hodgman_clipping_out.eps",
                 subjectPolygon, clippingPolygon, clipped);
}
Output:
immutable(Vec2)(100, 116.667)
immutable(Vec2)(125, 100)
immutable(Vec2)(275, 100)
immutable(Vec2)(300, 116.667)
immutable(Vec2)(300, 300)
immutable(Vec2)(250, 300)
immutable(Vec2)(200, 250)
immutable(Vec2)(175, 300)
immutable(Vec2)(125, 300)
immutable(Vec2)(100, 250)
sutherland_hodgman_clipping_out.eps written.

It also outputs an EPS file, the same as the C entry.

Elixir[edit]

Translation of: Ruby
defmodule SutherlandHodgman do
  defp inside(cp1, cp2, p), do: (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x)
  
  defp intersection(cp1, cp2, s, e) do
    {dcx, dcy} = {cp1.x-cp2.x, cp1.y-cp2.y}
    {dpx, dpy} = {s.x-e.x, s.y-e.y}
    n1 = cp1.x*cp2.y - cp1.y*cp2.x
    n2 = s.x*e.y - s.y*e.x
    n3 = 1.0 / (dcx*dpy - dcy*dpx)
    %{x: (n1*dpx - n2*dcx) * n3, y: (n1*dpy - n2*dcy) * n3}
  end
  
  def polygon_clipping(subjectPolygon, clipPolygon) do
    Enum.chunk([List.last(clipPolygon) | clipPolygon], 2, 1)
    |> Enum.reduce(subjectPolygon, fn [cp1,cp2],acc ->
         Enum.chunk([List.last(acc) | acc], 2, 1)
         |> Enum.reduce([], fn [s,e],outputList ->
              case {inside(cp1, cp2, e), inside(cp1, cp2, s)} do
                {true,  true} -> [e | outputList]
                {true, false} -> [e, intersection(cp1,cp2,s,e) | outputList]
                {false, true} -> [intersection(cp1,cp2,s,e) | outputList]
                _             -> outputList
              end
            end)
         |> Enum.reverse
       end)
  end
end

subjectPolygon = [[50, 150], [200, 50], [350, 150], [350, 300], [250, 300],
                  [200, 250], [150, 350], [100, 250], [100, 200]]
                 |> Enum.map(fn [x,y] -> %{x: x, y: y} end)

clipPolygon = [[100, 100], [300, 100], [300, 300], [100, 300]]
              |> Enum.map(fn [x,y] -> %{x: x, y: y} end)

SutherlandHodgman.polygon_clipping(subjectPolygon, clipPolygon)
|> Enum.each(&IO.inspect/1)
Output:
%{x: 100.0, y: 116.66666666666667}
%{x: 125.00000000000001, y: 100.0}
%{x: 275.0, y: 100.0}
%{x: 300.0, y: 116.66666666666667}
%{x: 300.0, y: 299.99999999999994}
%{x: 250.0, y: 300.0}
%{x: 200, y: 250}
%{x: 175.0, y: 300.0}
%{x: 125.0, y: 300.0}
%{x: 100.0, y: 250.0}

Fortran[edit]

Infos: The polygons are fortran type with an allocatable array "vertex" that contains the vertices and an integer n that is the size of the polygon. For any polygon, the first vertex and the last vertex have to be the same. As you will see, in the main function, we allocate the vertex array of the result polygon with its maximal size.

module SutherlandHodgmanUtil
  ! functions and type needed for Sutherland-Hodgman algorithm

  ! -------------------------------------------------------- !
  type polygon
    !type for polygons
    ! when you define a polygon, the first and the last vertices have to be the same
    integer :: n
    double precision, dimension(:,:), allocatable :: vertex
  end type polygon
  
  contains 
  
  ! -------------------------------------------------------- !
  subroutine sutherlandHodgman( ref, clip, outputPolygon )
    ! Sutherland Hodgman algorithm for 2d polygons
  
    ! -- parameters of the subroutine --
    type(polygon) :: ref, clip, outputPolygon
  
    ! -- variables used is the subroutine
    type(polygon) :: workPolygon               ! polygon clipped step by step 
    double precision, dimension(2) :: y1,y2    ! vertices of edge to clip workPolygon
    integer :: i  
  
    ! allocate workPolygon with the maximal possible size
    !   the sum of the size of polygon ref and clip
    allocate(workPolygon%vertex( ref%n+clip%n , 2 ))
    
    !  initialise the work polygon with clip
    workPolygon%n = clip%n
    workPolygon%vertex(1:workPolygon%n,:) = clip%vertex(1:workPolygon%n,:)

    do i=1,ref%n-1 ! for each edge i of the polygon ref
      y1(:) = ref%vertex(i,:)   !  vertex 1 of edge i
      y2(:) = ref%vertex(i+1,:) !  vertex 2 of edge i
  
      ! clip the work polygon by edge i
      call edgeClipping( workPolygon, y1, y2, outputPolygon)
      ! workPolygon <= outputPolygon
      workPolygon%n = outputPolygon%n
      workPolygon%vertex(1:workPolygon%n,:) = outputPolygon%vertex(1:workPolygon%n,:)

    end do 
    deallocate(workPolygon%vertex)
  end subroutine sutherlandHodgman
  
  ! -------------------------------------------------------- !
  subroutine edgeClipping( poly, y1, y2, outputPoly )
    ! make the clipping  of the polygon by the line (x1x2)
    
    type(polygon) :: poly, outputPoly
    double precision, dimension(2) :: y1, y2, x1, x2, intersecPoint
    integer ::  i, c
    
    c = 0 ! counter for the output polygon
    
    do i=1,poly%n-1 ! for each edge i of poly
      x1(:) = poly%vertex(i,:)   ! vertex 1 of edge i
      x2(:) = poly%vertex(i+1,:) ! vertex 2 of edge i
      
      if ( inside(x1, y1, y2) ) then ! if vertex 1 in inside clipping region
        if ( inside(x2, y1, y2) ) then ! if vertex 2 in inside clipping region
          ! add the vertex 2 to the output polygon
          c = c+1
          outputPoly%vertex(c,:) = x2(:)

        else ! vertex i+1 is outside
          intersecPoint = intersection(x1, x2, y1,y2)
          c = c+1
          outputPoly%vertex(c,:) = intersecPoint(:)
        end if
      else ! vertex i is outside
        if ( inside(x2, y1, y2) ) then
          intersecPoint = intersection(x1, x2, y1,y2)
          c = c+1
          outputPoly%vertex(c,:) = intersecPoint(:)
          
          c = c+1
          outputPoly%vertex(c,:) = x2(:)
        end if
      end if
    end do
    
    if (c .gt. 0) then
      ! if the last vertice is not equal to the first one
      if ( (outputPoly%vertex(1,1) .ne. outputPoly%vertex(c,1)) .or. & 
           (outputPoly%vertex(1,2) .ne. outputPoly%vertex(c,2)))  then
        c=c+1
        outputPoly%vertex(c,:) = outputPoly%vertex(1,:)
      end if
    end if
    ! set the size of the outputPolygon
    outputPoly%n = c
  end subroutine edgeClipping
  
  ! -------------------------------------------------------- !
  function intersection( x1, x2, y1, y2)
    ! computes the intersection between segment [x1x2] 
    ! and line the line (y1y2) 

    ! -- parameters of the function --
    double precision, dimension(2) :: x1, x2, &  ! points of the segment
                                      y1, y2     ! points of the line
    
    double precision, dimension(2) :: intersection, vx, vy, x1y1 
    double precision :: a
  
    vx(:) = x2(:) - x1(:) 
    vy(:) = y2(:) - y1(:)

    ! if the vectors are colinear
    if ( crossProduct(vx,vy) .eq. 0.d0) then
      x1y1(:) = y1(:) - x1(:)
      ! if the the segment [x1x2] is included in the line (y1y2)
      if ( crossProduct(x1y1,vx) .eq. 0.d0) then
        ! the intersection is the last point of the segment
        intersection(:) = x2(:)
      end if
    else ! the vectors are not colinear
      ! we want to find the inersection between [x1x2]
      ! and (y1,y2).
      ! mathematically, we want to find a in [0;1] such
      ! that :
      !     x1 + a vx = y1 + b vy        
      ! <=> a vx = x1y1 + b vy
      ! <=> a vx^vy = x1y1^vy      , ^ is cross product
      ! <=> a = x1y1^vy / vx^vy
     
      x1y1(:) = y1(:) - x1(:) 
      ! we compute a
      a = crossProduct(x1y1,vy)/crossProduct(vx,vy)
      ! if a is not in [0;1]
      if ( (a .gt. 1.d0) .or. (a .lt. 0)) then
        ! no intersection
      else
        intersection(:) = x1(:) + a*vx(:)
      end if
    end if

  end function intersection
  
  
  ! -------------------------------------------------------- !
  function inside( p, y1, y2)
    ! function that tells is the point p is at left of the line (y1y2)
    
    double precision, dimension(2) :: p, y1, y2, v1, v2
    logical :: inside
    v1(:) = y2(:) -  y1(:)
    v2(:) = p(:)  -  y1(:)  
    if ( crossProduct(v1,v2) .ge. 0.d0) then
      inside = .true.
    else 
      inside = .false.
    end if
   
   contains 
  end function inside

  ! -------------------------------------------------------- !
  function dotProduct( v1, v2)
    ! compute the dot product of vectors v1 and v2
    double precision, dimension(2) :: v1
    double precision, dimension(2) :: v2
    double precision :: dotProduct
    dotProduct = v1(1)*v2(1) + v1(2)*v2(2)
  end function dotProduct

  ! -------------------------------------------------------- !
  function crossProduct( v1, v2)
    ! compute the crossproduct of vectors v1 and v2
    double precision, dimension(2) :: v1
    double precision, dimension(2) :: v2
    double precision :: crossProduct
    crossProduct = v1(1)*v2(2) - v1(2)*v2(1)
  end function crossProduct

end module SutherlandHodgmanUtil

program main
  
  ! load the module for S-H algorithm
  use SutherlandHodgmanUtil, only : polygon, &
                                    sutherlandHodgman, &
                                    edgeClipping

  type(polygon) :: p1, p2, res
  integer :: c, n 
  double precision, dimension(2) :: y1, y2
  
  ! when you define a polygon, the first and the last vertices have to be the same

  ! first polygon
  p1%n = 10
  allocate(p1%vertex(p1%n,2))
  p1%vertex(1,1)=50.d0
  p1%vertex(1,2)=150.d0
  
  p1%vertex(2,1)=200.d0
  p1%vertex(2,2)=50.d0
  
  p1%vertex(3,1)= 350.d0
  p1%vertex(3,2)= 150.d0
  
  p1%vertex(4,1)= 350.d0
  p1%vertex(4,2)= 300.d0
  
  p1%vertex(5,1)= 250.d0
  p1%vertex(5,2)= 300.d0
  
  p1%vertex(6,1)= 200.d0
  p1%vertex(6,2)= 250.d0
  
  p1%vertex(7,1)= 150.d0
  p1%vertex(7,2)= 350.d0
  
  p1%vertex(8,1)= 100.d0
  p1%vertex(8,2)= 250.d0
  
  p1%vertex(9,1)= 100.d0
  p1%vertex(9,2)= 200.d0
  
  p1%vertex(10,1)=  50.d0
  p1%vertex(10,2)= 150.d0
 
  y1 = (/ 100.d0, 300.d0 /)
  y2 = (/ 300.d0, 300.d0 /)
  
  ! second polygon
  p2%n = 5
  allocate(p2%vertex(p2%n,2))

  p2%vertex(1,1)= 100.d0
  p2%vertex(1,2)= 100.d0
  
  p2%vertex(2,1)= 300.d0
  p2%vertex(2,2)= 100.d0
  
  p2%vertex(3,1)= 300.d0
  p2%vertex(3,2)= 300.d0
  
  p2%vertex(4,1)= 100.d0
  p2%vertex(4,2)= 300.d0
  
  p2%vertex(5,1)= 100.d0
  p2%vertex(5,2)= 100.d0
 
  allocate(res%vertex(p1%n+p2%n,2))
  call sutherlandHodgman( p2, p1, res)
  write(*,*) "Suterland-Hodgman"
  do c=1, res%n
    write(*,*) res%vertex(c,1), res%vertex(c,2)
  end do
  deallocate(res%vertex)

end program main

Output:

  Suterland-Hodgman
  300.00000000000000        300.00000000000000     
  250.00000000000000        300.00000000000000     
  200.00000000000000        250.00000000000000     
  175.00000000000000        300.00000000000000     
  125.00000000000000        300.00000000000000     
  100.00000000000000        250.00000000000000     
  100.00000000000000        200.00000000000000     
  100.00000000000000        200.00000000000000     
  100.00000000000000        116.66666666666667     
  125.00000000000000        100.00000000000000     
  275.00000000000000        100.00000000000000     
  300.00000000000000        116.66666666666666     
  300.00000000000000        300.00000000000000

FreeBASIC[edit]

FreeBASIC has inbuilt gfx graphics (a main feature), but I have no access to graphics uploads. So no extra credits.

Type Point
      As Double x,y
End Type

Type Line
      As Point s,f'start/finish
End Type

Function isleft(L As Line,p As Point) As Long
      Return  -Sgn((L.s.x-L.f.x)*(p.y-L.f.y)-(p.x-L.f.x)*(L.s.y-L.f.y))
End Function

Function segmentintersections(L1 As Line,L2 As Line) As Long
      If isleft(L2,L1.s) = isleft(L2,L1.f) Then Return 0 
      If isleft(L1,L2.s) = isleft(L1,L2.f) Then Return 0
      Return 1
End Function

Function allintersections(l1 As Line,l2 As Line,_out As Point) As Long
      Const tolerance=.01
      Var p1=l1.s, p2=l1.f, p3=l2.s, p4=l2.f
      Var x12=p1.x-p2.x, x34=p3.x-p4.x, y12=p1.y-p2.y, y34=p3.y-p4.y
      Var c=x12*y34-y12*x34
      If Abs(c)<tolerance Then Return 0
      Var a=p1.x*p2.y-p1.y*p2.x, b=p3.x*p4.y-p3.y*p4.x
      _out.x = (a*x34-b*x12)/c
      _out.y = (a*y34-b*y12)/c
      Return 1
End Function

Dim As Point p1(...)={(50,150),(200,50),(350,150),(350,300),(250,300),(200,250), _
                     (150,350),(100,250),(100,200)}
                     
Dim As Point p2(...)={(100,100),(300,100),(300,300),(100,300)}
'get the line segments around the polygons
Dim As Line L1(...)={(p1(0),p1(1)),(p1(1),p1(2)),(p1(2),p1(3)),(p1(3),p1(4)),(p1(4),p1(5)),_
                    (p1(5),p1(6)),(p1(6),p1(7)),(p1(7),p1(8)),(p1(8),p1(0))}

Dim As Line L2(...)={(p2(0),p2(1)),(p2(1),p2(2)),(p2(2),p2(3)),(p2(3),p2(0))}

'would normally draw these lines now, but not here.
Dim As Point x
For n1 As Long=Lbound(L1) To Ubound(L1)
      For n2 As Long=Lbound(L2) To Ubound(L2)
            If allintersections(L1(n1),L2(n2),x) And segmentintersections(L1(n1),L2(n2)) Then 
                  Print x.x,x.y
            End If
      Next
Next

Sleep
Output:
 125           100
 100           116.6666666666667
 275           100
 300           116.6666666666667
 300           300
 250           300
 175           300
 125           300
 100           250
 100           200


Go[edit]

No extra credit today.

package main

import "fmt"

type point struct {
    x, y float32
}

var subjectPolygon = []point{{50, 150}, {200, 50}, {350, 150}, {350, 300},
    {250, 300}, {200, 250}, {150, 350}, {100, 250}, {100, 200}}

var clipPolygon = []point{{100, 100}, {300, 100}, {300, 300}, {100, 300}}

func main() {
    var cp1, cp2, s, e point
    inside := func(p point) bool {
        return (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x)
    }
    intersection := func() (p point) {
        dcx, dcy := cp1.x-cp2.x, cp1.y-cp2.y
        dpx, dpy := s.x-e.x, s.y-e.y
        n1 := cp1.x*cp2.y - cp1.y*cp2.x
        n2 := s.x*e.y - s.y*e.x
        n3 := 1 / (dcx*dpy - dcy*dpx)
        p.x = (n1*dpx - n2*dcx) * n3
        p.y = (n1*dpy - n2*dcy) * n3
        return
    }
    outputList := subjectPolygon
    cp1 = clipPolygon[len(clipPolygon)-1]
    for _, cp2 = range clipPolygon { // WP clipEdge is cp1,cp2 here
        inputList := outputList
        outputList = nil
        s = inputList[len(inputList)-1]
        for _, e = range inputList {
            if inside(e) {
                if !inside(s) {
                    outputList = append(outputList, intersection())
                }
                outputList = append(outputList, e)
            } else if inside(s) {
                outputList = append(outputList, intersection())
            }
            s = e
        }
        cp1 = cp2
    }
    fmt.Println(outputList)
}
Output:
[{100 116.66667} {125 100} {275 100} {300 116.66667} {300 300} {250 300} {200 250} {175 300} {125 300} {100 250}]

(You can try it online)

Haskell[edit]

module SuthHodgClip (clipTo) where

import Data.List

type   Pt a = (a, a)
type   Ln a = (Pt a, Pt a)
type Poly a = [Pt a]

-- Return a polygon from a list of points.
polyFrom ps = last ps : ps

-- Return a list of lines from a list of points.
linesFrom pps@(_:ps) = zip pps ps

-- Return true if the point (x,y) is on or to the left of the oriented line
-- defined by (px,py) and (qx,qy).
(.|) :: (Num a, Ord a) => Pt a -> Ln a -> Bool
(x,y) .| ((px,py),(qx,qy)) = (qx-px)*(y-py) >= (qy-py)*(x-px)

-- Return the intersection of two lines.
(><) :: Fractional a => Ln a -> Ln a -> Pt a
((x1,y1),(x2,y2)) >< ((x3,y3),(x4,y4)) =
    let (r,s) = (x1*y2-y1*x2, x3*y4-y3*x4)
        (t,u,v,w) = (x1-x2, y3-y4, y1-y2, x3-x4)
        d = t*u-v*w 
    in ((r*w-t*s)/d, (r*u-v*s)/d)

-- Intersect the line segment (p0,p1) with the clipping line's left halfspace,
-- returning the point closest to p1.  In the special case where p0 lies outside
-- the halfspace and p1 lies inside we return both the intersection point and
-- p1.  This ensures we will have the necessary segment along the clipping line.
(-|) :: (Fractional a, Ord a) => Ln a -> Ln a -> [Pt a]
ln@(p0, p1) -| clipLn =
    case (p0 .| clipLn, p1 .| clipLn) of
      (False, False) -> []
      (False, True)  -> [isect, p1]
      (True,  False) -> [isect]
      (True,  True)  -> [p1]
    where isect = ln >< clipLn

-- Intersect the polygon with the clipping line's left halfspace.
(<|) :: (Fractional a, Ord a) => Poly a -> Ln a -> Poly a
poly <| clipLn = polyFrom $ concatMap (-| clipLn) (linesFrom poly)

-- Intersect a target polygon with a clipping polygon.  The latter is assumed to
-- be convex.
clipTo :: (Fractional a, Ord a) => [Pt a] -> [Pt a] -> [Pt a]
targPts `clipTo` clipPts = 
    let targPoly = polyFrom targPts
        clipLines = linesFrom (polyFrom clipPts)
    in foldl' (<|) targPoly clipLines

Print the resulting list of points and display the polygons in a window.

import Graphics.HGL
import SuthHodgClip

targPts = [( 50,150), (200, 50), (350,150), (350,300), (250,300), 
           (200,250), (150,350), (100,250), (100,200)] :: [(Float,Float)]
clipPts = [(100,100), (300,100), (300,300), (100,300)] :: [(Float,Float)]

toInts = map (\(a,b) -> (round a, round b))
complete xs = last xs : xs

drawSolid w c = drawInWindow w . withRGB c . polygon
drawLines w p = drawInWindow w . withPen p . polyline . toInts . complete

blue  = RGB 0x99 0x99 0xff
green = RGB 0x99 0xff 0x99
pink  = RGB 0xff 0x99 0x99
white = RGB 0xff 0xff 0xff

main = do
  let resPts = targPts `clipTo` clipPts
      sz = 400
      win = [(0,0), (sz,0), (sz,sz), (0,sz)]
  runWindow "Sutherland-Hodgman Polygon Clipping" (sz,sz) $ \w -> do
         print $ toInts resPts
         penB <- createPen Solid 3 blue
         penP <- createPen Solid 5 pink
         drawSolid w white win
         drawLines w penB targPts
         drawLines w penP clipPts
         drawSolid w green $ toInts resPts
         getKey w
Output:
[(100,200),(100,200),(100,117),(125,100),(275,100),(300,117),(300,300),(250,300),(200,250),(175,300),(125,300),(100,250),(100,200)]

Sutherland-Hodgman haskell.png

J[edit]

Solution:

NB. assumes counterclockwise orientation.
NB. determine whether point y is inside edge x.
isinside=:0< [:-/ .* {.@[ -~"1 {:@[,:]
 
NB. (p0,:p1) intersection (p2,:p3)
intersection=:|:@[ (+/ .* (,-.)) [:{. ,.&(-~/) %.~ -&{:
 
SutherlandHodgman=:4 :0 NB. clip S-H subject
  clip=.2 ]\ (,{.) x
  subject=.y
  for_edge. clip do.
    S=.{:input=.subject
    subject=.0 2$0
    for_E. input do.
      if. edge isinside E do.
        if. -.edge isinside S do.
          subject=.subject,edge intersection S,:E end.
        subject=.subject,E
      elseif. edge isinside S do.
        subject=.subject,edge intersection S,:E end.
      S=.E
    end.
  end.
  subject
)
Example use:
   subject=: 50 150,200 50,350 150,350 300,250 300,200 250,150 350,100 250,:100 200
   clip=: 100 100,300 100,300 300,:100 300
   clip SutherlandHodgman subject
100 116.667
125     100
275     100
300 116.667
300     300
250     300
200     250
175     300
125     300
100     250

Java[edit]

Works with: Java version 7
import java.awt.*;
import java.awt.geom.Line2D;
import java.util.*;
import java.util.List;
import javax.swing.*;

public class SutherlandHodgman extends JFrame {

    SutherlandHodgmanPanel panel;

    public static void main(String[] args) {
        JFrame f = new SutherlandHodgman();
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setVisible(true);
    }

    public SutherlandHodgman() {
        Container content = getContentPane();
        content.setLayout(new BorderLayout());
        panel = new SutherlandHodgmanPanel();
        content.add(panel, BorderLayout.CENTER);
        setTitle("SutherlandHodgman");
        pack();
        setLocationRelativeTo(null);
    }
}

class SutherlandHodgmanPanel extends JPanel {
    List<double[]> subject, clipper, result;

    public SutherlandHodgmanPanel() {
        setPreferredSize(new Dimension(600, 500));

        // these subject and clip points are assumed to be valid
        double[][] subjPoints = {{50, 150}, {200, 50}, {350, 150}, {350, 300},
        {250, 300}, {200, 250}, {150, 350}, {100, 250}, {100, 200}};

        double[][] clipPoints = {{100, 100}, {300, 100}, {300, 300}, {100, 300}};

        subject = new ArrayList<>(Arrays.asList(subjPoints));
        result  = new ArrayList<>(subject);
        clipper = new ArrayList<>(Arrays.asList(clipPoints));

        clipPolygon();
    }

    private void clipPolygon() {
        int len = clipper.size();
        for (int i = 0; i < len; i++) {

            int len2 = result.size();
            List<double[]> input = result;
            result = new ArrayList<>(len2);

            double[] A = clipper.get((i + len - 1) % len);
            double[] B = clipper.get(i);

            for (int j = 0; j < len2; j++) {

                double[] P = input.get((j + len2 - 1) % len2);
                double[] Q = input.get(j);

                if (isInside(A, B, Q)) {
                    if (!isInside(A, B, P))
                        result.add(intersection(A, B, P, Q));
                    result.add(Q);
                } else if (isInside(A, B, P))
                    result.add(intersection(A, B, P, Q));
            }
        }
    }

    private boolean isInside(double[] a, double[] b, double[] c) {
        return (a[0] - c[0]) * (b[1] - c[1]) > (a[1] - c[1]) * (b[0] - c[0]);
    }

    private double[] intersection(double[] a, double[] b, double[] p, double[] q) {
        double A1 = b[1] - a[1];
        double B1 = a[0] - b[0];
        double C1 = A1 * a[0] + B1 * a[1];

        double A2 = q[1] - p[1];
        double B2 = p[0] - q[0];
        double C2 = A2 * p[0] + B2 * p[1];

        double det = A1 * B2 - A2 * B1;
        double x = (B2 * C1 - B1 * C2) / det;
        double y = (A1 * C2 - A2 * C1) / det;

        return new double[]{x, y};
    }

    @Override
    public void paintComponent(Graphics g) {
        super.paintComponent(g);
        Graphics2D g2 = (Graphics2D) g;
        g2.translate(80, 60);
        g2.setStroke(new BasicStroke(3));
        g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING,
                RenderingHints.VALUE_ANTIALIAS_ON);

        drawPolygon(g2, subject, Color.blue);
        drawPolygon(g2, clipper, Color.red);
        drawPolygon(g2, result, Color.green);
    }

    private void drawPolygon(Graphics2D g2, List<double[]> points, Color color) {
        g2.setColor(color);
        int len = points.size();
        Line2D line = new Line2D.Double();
        for (int i = 0; i < len; i++) {
            double[] p1 = points.get(i);
            double[] p2 = points.get((i + 1) % len);
            line.setLine(p1[0], p1[1], p2[0], p2[1]);
            g2.draw(line);
        }
    }
}

JavaScript[edit]

Solution:

<html>
    <head>
	<script>
        function clip (subjectPolygon, clipPolygon) {
            
            var cp1, cp2, s, e;
            var inside = function (p) {
                return (cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0]);
            };
            var intersection = function () {
                var dc = [ cp1[0] - cp2[0], cp1[1] - cp2[1] ],
                    dp = [ s[0] - e[0], s[1] - e[1] ],
                    n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0],
                    n2 = s[0] * e[1] - s[1] * e[0], 
                    n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0]);
                return [(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3];
            };
            var outputList = subjectPolygon;
            cp1 = clipPolygon[clipPolygon.length-1];
            for (var j in clipPolygon) {
                cp2 = clipPolygon[j];
                var inputList = outputList;
                outputList = [];
                s = inputList[inputList.length - 1]; //last on the input list
                for (var i in inputList) {
                    e = inputList[i];
                    if (inside(e)) {
                        if (!inside(s)) {
                            outputList.push(intersection());
                        }
                        outputList.push(e);
                    }
                    else if (inside(s)) {
                        outputList.push(intersection());
                    }
                    s = e;
                }
                cp1 = cp2;
            }
            return outputList
        }

        function drawPolygon(context, polygon, strokeStyle, fillStyle) {
            context.strokeStyle = strokeStyle;
            context.fillStyle = fillStyle;
            context.beginPath();
            context.moveTo(polygon[0][0],polygon[0][1]); //first vertex
            for (var i = 1; i < polygon.length ; i++)
                context.lineTo(polygon[i][0],polygon[i][1]);
            context.lineTo(polygon[0][0],polygon[0][1]); //back to start
            context.fill();
            context.stroke();
            context.closePath();
        }

        window.onload = function () {
	        var context = document.getElementById('canvas').getContext('2d');
	        var subjectPolygon = [[50, 150], [200, 50], [350, 150], [350, 300], [250, 300], [200, 250], [150, 350], [100, 250], [100, 200]],
	            clipPolygon = [[100, 100], [300, 100], [300, 300], [100, 300]];
	        var clippedPolygon = clip(subjectPolygon, clipPolygon);
	        drawPolygon(context, clipPolygon, '#888','#88f');
	        drawPolygon(context, subjectPolygon, '#888','#8f8');
	        drawPolygon(context, clippedPolygon, '#000','#0ff');
    	}
        </script>
    <body>
    	<canvas id='canvas' width='400' height='400'></canvas>
    </body>
</html>

You can see it running here

Julia[edit]

using Luxor

isinside(p, a, b) = (b.x - a.x) * (p.y - a.y) > (b.y - a.y) * (p.x - a.x)

function intersection(a, b, s, f)
    dc = [a.x - b.x, a.y - b.y]
    dp = [s.x - f.x, s.y - f.y]
    n1 = a.x * b.y - a.y * b.x
    n2 = s.x * f.y - s.y * f.x
    n3 = 1.0 / (dc[1] * dp[2] - dc[2] * dp[1])
    Point((n1 * dp[1] - n2 * dc[1]) * n3, (n1 * dp[2] - n2 * dc[2]) * n3)
end

function clipSH(spoly, cpoly)
    outarr = spoly
    q = cpoly[end]
    for p in cpoly
        inarr = outarr
        outarr = Point[]
        s = inarr[end]
        for vtx in inarr
            if isinside(vtx, q, p)
                if !isinside(s, q, p)
                    push!(outarr, intersection(q, p, s, vtx))
                end
                push!(outarr, vtx)
            elseif isinside(s, q, p)
                push!(outarr, intersection(q, p, s, vtx))
            end
            s = vtx
        end
        q = p
    end
    outarr
end

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

clipp = [Point(100, 100), Point(300, 100), Point(300, 300), Point(100, 300)]

Drawing(400, 400, "intersecting-polygons.png")
background("white")
sethue("red")
poly(subjectp, :stroke, close=true)
sethue("blue")
poly(clipp, :stroke, close=true)
clipped = clipSH(subjectp, clipp)
sethue("gold")
poly(clipped, :fill, close=true)
finish()
preview()
println(clipped)
Output:
Point[Point(100.0, 116.667), Point(125.0, 100.0), Point(275.0, 100.0), Point(300.0, 116.667),
    Point(300.0, 300.0), Point(250.0, 300.0), Point(200.0, 250.0), Point(175.0, 300.0), 
    Point(125.0, 300.0), Point(100.0, 250.0)]

Kotlin[edit]

Translation of: Java
// version 1.1.2

import java.awt.*
import java.awt.geom.Line2D
import javax.swing.*
 
class SutherlandHodgman : JPanel() {
    private val subject = listOf(
        doubleArrayOf( 50.0, 150.0), doubleArrayOf(200.0,  50.0), doubleArrayOf(350.0, 150.0), 
        doubleArrayOf(350.0, 300.0), doubleArrayOf(250.0, 300.0), doubleArrayOf(200.0, 250.0), 
        doubleArrayOf(150.0, 350.0), doubleArrayOf(100.0, 250.0), doubleArrayOf(100.0, 200.0)
    )

    private val clipper = listOf(
        doubleArrayOf(100.0, 100.0), doubleArrayOf(300.0, 100.0), 
        doubleArrayOf(300.0, 300.0), doubleArrayOf(100.0, 300.0)
    )

    private var result = subject.toMutableList()

    init {
        preferredSize = Dimension(600, 500)
        clipPolygon()
    } 

    private fun clipPolygon() {
        val len = clipper.size
        for (i in 0 until len) { 
            val len2 = result.size
            val input = result
            result = mutableListOf<DoubleArray>() 
            val a = clipper[(i + len - 1) % len]
            val b = clipper[i]
 
            for (j in 0 until len2) {
                val p = input[(j + len2 - 1) % len2]
                val q = input[j]
 
                if (isInside(a, b, q)) {
                    if (!isInside(a, b, p)) result.add(intersection(a, b, p, q))
                    result.add(q)
                } 
                else if (isInside(a, b, p)) result.add(intersection(a, b, p, q))
            }
        }
    } 

    private fun isInside(a: DoubleArray, b: DoubleArray, c: DoubleArray) =
        (a[0] - c[0]) * (b[1] - c[1]) > (a[1] - c[1]) * (b[0] - c[0])

    private fun intersection(a: DoubleArray, b: DoubleArray, 
                             p: DoubleArray, q: DoubleArray): DoubleArray {
        val a1 = b[1] - a[1]
        val b1 = a[0] - b[0]
        val c1 = a1 * a[0] + b1 * a[1]
 
        val a2 = q[1] - p[1]
        val b2 = p[0] - q[0]
        val c2 = a2 * p[0] + b2 * p[1]
 
        val d = a1 * b2 - a2 * b1
        val x = (b2 * c1 - b1 * c2) / d
        val y = (a1 * c2 - a2 * c1) / d
 
        return doubleArrayOf(x, y)
    }

    override fun paintComponent(g: Graphics) {
        super.paintComponent(g)
        val g2 = g as Graphics2D         
        g2.translate(80, 60)
        g2.stroke = BasicStroke(3.0f)
        g2.setRenderingHint(RenderingHints.KEY_ANTIALIASING,
                            RenderingHints.VALUE_ANTIALIAS_ON) 
        drawPolygon(g2, subject, Color.blue)
        drawPolygon(g2, clipper, Color.red)
        drawPolygon(g2, result, Color.green)
    }

    private fun drawPolygon(g2: Graphics2D, points: List<DoubleArray>, color: Color) {
        g2.color = color
        val len = points.size
        val line = Line2D.Double()
        for (i in 0 until len) {
            val p1 = points[i]
            val p2 = points[(i + 1) % len]
            line.setLine(p1[0], p1[1], p2[0], p2[1])
            g2.draw(line)
        }
    }
}

fun main(args: Array<String>) {
    SwingUtilities.invokeLater {
        val f = JFrame()
        with(f) {
            defaultCloseOperation = JFrame.EXIT_ON_CLOSE
            add(SutherlandHodgman(), BorderLayout.CENTER)
            title = "Sutherland-Hodgman"
            pack()
            setLocationRelativeTo(null)
            isVisible = true
        }
    }
}

Lua[edit]

No extra credit.

Translation of: Go
subjectPolygon = {
  {50, 150}, {200, 50}, {350, 150}, {350, 300},
  {250, 300}, {200, 250}, {150, 350}, {100, 250}, {100, 200}
}
 
clipPolygon = {{100, 100}, {300, 100}, {300, 300}, {100, 300}}

function inside(p, cp1, cp2)
  return (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x)
end

function intersection(cp1, cp2, s, e)
  local dcx, dcy = cp1.x-cp2.x, cp1.y-cp2.y
  local dpx, dpy = s.x-e.x, s.y-e.y
  local n1 = cp1.x*cp2.y - cp1.y*cp2.x
  local n2 = s.x*e.y - s.y*e.x
  local n3 = 1 / (dcx*dpy - dcy*dpx)
  local x = (n1*dpx - n2*dcx) * n3
  local y = (n1*dpy - n2*dcy) * n3
  return {x=x, y=y}
end

function clip(subjectPolygon, clipPolygon)
  local outputList = subjectPolygon
  local cp1 = clipPolygon[#clipPolygon]
  for _, cp2 in ipairs(clipPolygon) do  -- WP clipEdge is cp1,cp2 here
    local inputList = outputList
    outputList = {}
    local s = inputList[#inputList]
    for _, e in ipairs(inputList) do
      if inside(e, cp1, cp2) then
        if not inside(s, cp1, cp2) then
          outputList[#outputList+1] = intersection(cp1, cp2, s, e)
        end
        outputList[#outputList+1] = e
      elseif inside(s, cp1, cp2) then
        outputList[#outputList+1] = intersection(cp1, cp2, s, e)
      end
      s = e
    end
    cp1 = cp2
  end
  return outputList
end

function main()
  local function mkpoints(t)
    for i, p in ipairs(t) do
      p.x, p.y = p[1], p[2]
    end
  end
  mkpoints(subjectPolygon)
  mkpoints(clipPolygon)

  local outputList = clip(subjectPolygon, clipPolygon)

  for _, p in ipairs(outputList) do
    print(('{%f, %f},'):format(p.x, p.y))
  end
end

main()
Output:
{100.000000, 116.666667},
{125.000000, 100.000000},
{275.000000, 100.000000},
{300.000000, 116.666667},
{300.000000, 300.000000},
{250.000000, 300.000000},
{200.000000, 250.000000},
{175.000000, 300.000000},
{125.000000, 300.000000},
{100.000000, 250.000000},

(You can also see it live)

Mathematica/Wolfram Language[edit]

Geometry is built in to the Wolfram Language.

p1 = Polygon[{{50, 150}, {200, 50}, {350, 150}, {350, 300}, {250, 300}, {200, 250}, {150, 350}, {100, 250}, {100, 200}}];
p2 = Polygon[{{100, 100}, {300, 100}, {300, 300}, {100, 300}}];
RegionIntersection[p1, p2]
Graphics[{Red, p1, Blue, p2, Green, RegionIntersection[p1, p2]}]
Output:
Polygon[{{125, 100}, {100, 350/3}, {100, 200}, {100, 250}, {125, 300}, {175, 300}, {200, 250}, {250, 300}, {300, 300}, {300, 350/3}, {275, 100}}]

MATLAB / Octave[edit]

%The inputs are a table of x-y pairs for the verticies of the subject
%polygon and boundary polygon. (x values in column 1 and y values in column
%2) The output is a table of x-y pairs for the clipped version of the 
%subject polygon.

function clippedPolygon = sutherlandHodgman(subjectPolygon,clipPolygon)

%% Helper Functions

    %computerIntersection() assumes the two lines intersect
    function intersection = computeIntersection(line1,line2)

        %this is an implementation of
        %http://en.wikipedia.org/wiki/Line-line_intersection
        
        intersection = zeros(1,2);

        detL1 = det(line1);
        detL2 = det(line2);

        detL1x = det([line1(:,1),[1;1]]);
        detL1y = det([line1(:,2),[1;1]]);

        detL2x = det([line2(:,1),[1;1]]);
        detL2y = det([line2(:,2),[1;1]]);

        denominator = det([detL1x detL1y;detL2x detL2y]);

        intersection(1) = det([detL1 detL1x;detL2 detL2x]) / denominator;
        intersection(2) = det([detL1 detL1y;detL2 detL2y]) / denominator;

    end %computeIntersection

    %inside() assumes the boundary is oriented counter-clockwise
    function in = inside(point,boundary)
       
        pointPositionVector = [diff([point;boundary(1,:)]) 0];
        boundaryVector = [diff(boundary) 0];
        crossVector = cross(pointPositionVector,boundaryVector);
        
        if ( crossVector(3) <= 0 )
            in = true;
        else
            in = false;
        end
        
    end %inside

%% Sutherland-Hodgman Algorithm

    clippedPolygon = subjectPolygon;
    numVerticies = size(clipPolygon,1);
    clipVertexPrevious = clipPolygon(end,:);
    
    for clipVertex = (1:numVerticies)
       
        clipBoundary = [clipPolygon(clipVertex,:) ; clipVertexPrevious];
        
        inputList = clippedPolygon;
        
        clippedPolygon = [];
        if ~isempty(inputList),
            previousVertex = inputList(end,:);
        end
        
        for subjectVertex = (1:size(inputList,1))

            if ( inside(inputList(subjectVertex,:),clipBoundary) )
                
                if( not(inside(previousVertex,clipBoundary)) )  
                    subjectLineSegment = [previousVertex;inputList(subjectVertex,:)];
                    clippedPolygon(end+1,1:2) = computeIntersection(clipBoundary,subjectLineSegment);
                end
                
                clippedPolygon(end+1,1:2) = inputList(subjectVertex,:);
                
            elseif( inside(previousVertex,clipBoundary) )
                    subjectLineSegment = [previousVertex;inputList(subjectVertex,:)];
                    clippedPolygon(end+1,1:2) = computeIntersection(clipBoundary,subjectLineSegment);                            
            end
            
            previousVertex = inputList(subjectVertex,:);
            clipVertexPrevious = clipPolygon(clipVertex,:);
            
        end %for subject verticies                
    end %for boundary verticies
end %sutherlandHodgman
Output:
>> subject = [[50;200;350;350;250;200;150;100;100],[150;50;150;300;300;250;350;250;200]];
>> clipPolygon = [[100;300;300;100],[100;100;300;300]];
>> clippedSubject = sutherlandHodgman(subject,clipPolygon);
>> plot([subject(:,1);subject(1,1)],[subject(:,2);subject(1,2)],[0,0,1])
>> hold on
>> plot([clipPolygon(:,1);clipPolygon(1,1)],[clipPolygon(:,2);clipPolygon(1,2)],'r')
>> patch(clippedSubject(:,1),clippedSubject(:,2),0);
>> axis square

Sutherland-Hodgman MATLAB.png

Mercury[edit]

Translation of: ATS
Works with: Mercury version 22.01.1
:- module sutherland_hodgman_task.

:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.

:- implementation.
:- import_module exception.
:- import_module float.
:- import_module list.
:- import_module pair.
:- import_module string.

:- type plane_point == pair(float).
:- func xcoord(plane_point) = float.
:- func ycoord(plane_point) = float.
:- func plane_point(float, float) = plane_point.
:- pred write_plane_point(plane_point::in, io::di, io::uo) is det.
:- pred write_plane_point_list(list(plane_point)::in, string::in,
                               io::di, io::uo) is det.
xcoord(Pt) = fst(Pt).
ycoord(Pt) = snd(Pt).
plane_point(X, Y) = pair(X, Y).
write_plane_point(Pt, !IO) :-
  write_string("(", !IO),
  write_float(xcoord(Pt), !IO),
  write_string(", ", !IO),
  write_float(ycoord(Pt), !IO),
  write_string(")", !IO).
write_plane_point_list(Pts, Separator, !IO) :-
  write_list(Pts, Separator, write_plane_point, !IO).

:- type plane_edge == pair(plane_point).
:- func point0(plane_edge) = plane_point.
:- func point1(plane_edge) = plane_point.
:- func plane_edge(plane_point, plane_point) = plane_edge.
point0(Edge) = fst(Edge).
point1(Edge) = snd(Edge).
plane_edge(Pt0, Pt1) = pair(Pt0, Pt1).

:- func evaluate_line(float, float, float, float, float) = float.
evaluate_line(X1, Y1, X2, Y2, X) = Y :-
  %% Given the line (X1,Y1)--(X2,Y2), evaluate it at X.
  Dy = Y2 - Y1,
  Dx = X2 - X1,
  Slope = Dy / Dx,
  Intercept = ((Dx * Y1) - (Dy * X1)) / Dx,
  Y = (Slope * X) + Intercept.

:- func intersection_of_lines(float, float, float, float,
                              float, float, float, float)
   = plane_point.
intersection_of_lines(X1, Y1, X2, Y2, X3, Y3, X4, Y4) = Pt :-
  %% Given the lines (X1,Y1)--(X2,Y2) and (X3,Y3)--(X3,Y4), find their
  %% point of intersection.
  (if (X1 = X2)
   then (Pt = plane_point(X1, evaluate_line(X3, Y3, X4, Y4, X1)))
   else if (X3 = X4)
   then (Pt = plane_point(X3, evaluate_line(X1, Y1, X2, Y2, X3)))
   else (Pt = plane_point(X, Y),
         X = Xnumerator / Denominator,
         Y = Ynumerator / Denominator,
         Denominator =
         ((X1 - X2) * (Y3 - Y4)) - ((Y1 - Y2) * (X3 - X4)),
         Xnumerator =
         (X1Y2_Y1X2 * (X3 - X4)) - ((X1 - X2) * X3Y4_Y3X4),
         Ynumerator =
         (X1Y2_Y1X2 * (Y3 - Y4)) - ((Y1 - Y2) * X3Y4_Y3X4),
         X1Y2_Y1X2 = (X1 * Y2) - (Y1 * X2),
         X3Y4_Y3X4 = (X3 * Y4) - (Y3 * X4))).

:- func intersection_of_edges(plane_edge, plane_edge) = plane_point.
intersection_of_edges(E1, E2) = Pt :-
  %% Given two edges, find their point of intersection (on the
  %% assumption that there is such an intersection).
  Pt = intersection_of_lines(X1, Y1, X2, Y2, X3, Y3, X4, Y4),
  Pt1 = point0(E1), Pt2 = point1(E1),
  Pt3 = point0(E2), Pt4 = point1(E2),
  X1 = xcoord(Pt1), Y1 = ycoord(Pt1),
  X2 = xcoord(Pt2), Y2 = ycoord(Pt2),
  X3 = xcoord(Pt3), Y3 = ycoord(Pt3),
  X4 = xcoord(Pt4), Y4 = ycoord(Pt4).

:- pred point_is_left_of_edge(plane_point::in, plane_edge::in)
   is semidet.
point_is_left_of_edge(Pt, Edge) :-
  %% Is Pt left of Edge?
  (OP < 0.0),
  %% OP = outer product of the vectors (x1,y1)-->(x,y) and
  %% (x1,y1)-->(x2,y2). *)
  OP = ((X - X1) * (Y2 - Y1)) - ((X2 - X1) * (Y - Y1)),
  Pt1 = point0(Edge), Pt2 = point1(Edge),
  X1 = xcoord(Pt1), Y1 = ycoord(Pt1),
  X2 = xcoord(Pt2), Y2 = ycoord(Pt2),
  X = xcoord(Pt), Y = ycoord(Pt).

:- func clip_subject_edge(plane_edge, plane_edge,
                          list(plane_point)) = list(plane_point).
clip_subject_edge(Subject_edge, Clip_edge, Accum0) = Accum :-
  S1 = point0(Subject_edge), S2 = point1(Subject_edge),
  (if (point_is_left_of_edge(S2, Clip_edge))
   then (if (point_is_left_of_edge(S1, Clip_edge))
         then (Accum = [S2 | Accum0])
         else (Accum = [S2, Intersection | Accum0],
               Intersection =
               intersection_of_edges(Subject_edge, Clip_edge)))
   else (if (point_is_left_of_edge(S1, Clip_edge))
         then (Accum = [Intersection | Accum0],
               Intersection =
               intersection_of_edges(Subject_edge, Clip_edge))
         else (Accum = Accum0))).

:- func plane_points_to_plane_edges(list(plane_point))
   = list(plane_edge).
plane_points_to_plane_edges(Pts) = Edges :-
  plane_points_to_plane_edges_(Pt_first, Pts, [], Edges),
  Pt_first = det_head(Pts).

:- pred plane_points_to_plane_edges_(plane_point::in,
                                     list(plane_point)::in,
                                     list(plane_edge)::in,
                                     list(plane_edge)::out) is det.
%% Convert a list of points to a list of edges.
plane_points_to_plane_edges_(Pt_first, [Pt0, Pt1 | Rest],
                             Edges0, Edges) :-
  plane_points_to_plane_edges_(Pt_first, [Pt1 | Rest],
                               [plane_edge(Pt0, Pt1) | Edges0],
                               Edges).
plane_points_to_plane_edges_(Pt_first, [Pt_last], Edges0, Edges) :-
  Edges = [plane_edge(Pt_last, Pt_first) | reverse(Edges0)].
plane_points_to_plane_edges_(_, [], _, _) :-
  throw("list(plane_point) was expected to have length >= 2").

:- pred for_each_subject_edge(list(plane_edge)::in, plane_edge::in,
                              list(plane_point)::in,
                              list(plane_point)::out) is det.
for_each_subject_edge([], _, Accum0, Accum) :-
  (Accum = reverse(Accum0)).
for_each_subject_edge([Subject_edge | Rest], Clip_edge,
                      Accum0, Accum) :-
  Accum1 = clip_subject_edge(Subject_edge, Clip_edge, Accum0),
  for_each_subject_edge(Rest, Clip_edge, Accum1, Accum).

:- func clip_subject_with_clip_edge(list(plane_point), plane_edge)
   = list(plane_point).
clip_subject_with_clip_edge(Subject_pts, Clip_edge) = Pts :-
  for_each_subject_edge(Subject_edges, Clip_edge, [], Pts),
  Subject_edges = plane_points_to_plane_edges(Subject_pts).

:- pred for_each_clip_edge(list(plane_point)::in,
                           list(plane_point)::out,
                           list(plane_edge)::in) is det.
for_each_clip_edge(Subject_pts0, Subject_pts, []) :-
  (Subject_pts = Subject_pts0).
for_each_clip_edge(Subject_pts0, Subject_pts,
                   [Clip_edge | Rest]) :-
  Subject_pts1 = clip_subject_with_clip_edge(Subject_pts0, Clip_edge),
  for_each_clip_edge(Subject_pts1, Subject_pts, Rest).

:- func clip(list(plane_point), list(plane_point))
   = list(plane_point).
clip(Subject_pts, Clip_pts) = Result_pts :-
  for_each_clip_edge(Subject_pts, Result_pts, Clip_edges),
  Clip_edges = plane_points_to_plane_edges(Clip_pts).

:- pred moveto(text_output_stream::in, plane_point::in,
               io::di, io::uo) is det.
moveto(Stream, Pt, !IO) :-
  write_float(Stream, xcoord(Pt), !IO),
  write_string(Stream, " ", !IO),
  write_float(Stream, ycoord(Pt), !IO),
  write_string(Stream, " moveto\n", !IO).

:- pred lineto(plane_point::in, io::di, io::uo) is det.
lineto(Pt, !IO) :-
  write_float(xcoord(Pt), !IO),
  write_string(" ", !IO),
  write_float(ycoord(Pt), !IO),
  write_string(" lineto\n", !IO).

:- pred setrgbcolor(text_output_stream::in,
                    string::in, io::di, io::uo) is det.
setrgbcolor(Stream, Color, !IO) :-
  write_string(Stream, Color, !IO),
  write_string(Stream, " setrgbcolor\n", !IO).

:- pred write_polygon(text_output_stream::in,
                      list(plane_point)::in,
                      string::in, string::in,
                      io::di, io::uo) is det.
write_polygon(Stream, Pts, Line_color, Fill_color, !IO) :-
  if ([First_pt | Rest] = Pts)
  then (moveto(Stream, First_pt, !IO),
        write_list(Stream, Rest, "", lineto, !IO),
        write_string(Stream, "closepath\n", !IO),
        setrgbcolor(Stream, Line_color, !IO),
        write_string(Stream, "gsave\n", !IO),
        setrgbcolor(Stream, Fill_color, !IO),
        write_string(Stream, "fill\n", !IO),
        write_string(Stream, "grestore\n", !IO),
        write_string(Stream, "stroke\n", !IO))
  else true.

:- pred write_eps(text_output_stream::in,
                  list(plane_point)::in,
                  list(plane_point)::in,
                  list(plane_point)::in,
                  io::di, io::uo) is det.
write_eps(Stream, Subject_pts, Clip_pts, Result_pts, !IO) :-
  write_string(Stream, "%!PS-Adobe-3.0 EPSF-3.0\n", !IO),
  write_string(Stream, "%%BoundingBox: 40 40 360 360\n", !IO),
  write_string(Stream, "0 setlinewidth\n", !IO),
  write_polygon(Stream, Clip_pts, ".5 0 0", "1 .7 .7", !IO),
  write_polygon(Stream, Subject_pts, "0 .2 .5", ".4 .7 1", !IO),
  write_string(Stream, "2 setlinewidth\n", !IO),
  write_string(Stream, "[10 8] 0 setdash\n", !IO),
  write_polygon(Stream, Result_pts, ".5 0 .5", ".7 .3 .8", !IO),
  write_string(Stream, "%%EOF\n", !IO).

:- pred write_eps_to_file(string::in,
                          list(plane_point)::in,
                          list(plane_point)::in,
                          list(plane_point)::in,
                          io::di, io::uo) is det.
write_eps_to_file(Filename, Subject_pts, Clip_pts, Result_pts, !IO) :-
  open_output(Filename, Open_result, !IO),
  (if (Open_result = ok(Outp))
   then write_eps(Outp, Subject_pts, Clip_pts, Result_pts, !IO)
   else throw("Failed to open " ++ Filename ++ " for output.")).

main(!IO) :-
  Subject_pts = [plane_point(50.0, 150.0),
                 plane_point(200.0, 50.0),
                 plane_point(350.0, 150.0),
                 plane_point(350.0, 300.0),
                 plane_point(250.0, 300.0),
                 plane_point(200.0, 250.0),
                 plane_point(150.0, 350.0),
                 plane_point(100.0, 250.0),
                 plane_point(100.0, 200.0)],
  Clip_pts = [plane_point(100.0, 100.0),
              plane_point(300.0, 100.0),
              plane_point(300.0, 300.0),
              plane_point(100.0, 300.0)],
  Result_pts = clip(Subject_pts, Clip_pts),
  write_plane_point_list(Result_pts, "\n", !IO), nl(!IO),
  EPSF = "sutherland-hodgman.eps",
  write_eps_to_file(EPSF, Subject_pts, Clip_pts, Result_pts, !IO),
  write_string("Wrote " ++ EPSF, !IO), nl(!IO).

%%% local variables:
%%% mode: mercury
%%% prolog-indent-width: 2
%%% end:
Output:
$ mmc sutherland_hodgman_task.m && ./sutherland_hodgman_task
(100.0, 116.66666666666669)
(124.99999999999999, 100.0)
(275.0, 100.0)
(300.0, 116.66666666666667)
(300.0, 300.0)
(250.0, 300.0)
(200.0, 250.0)
(175.0, 300.0)
(125.0, 300.0)
(100.0, 250.0)
Wrote sutherland-hodgman.eps

The polygons as generated by a Mercury program.

Modula-2[edit]

Translation of: ATS
Works with: GNU Modula-2 version 13.0.0 20220926 (experimental)
(* Sutherland-Hodgman polygon clipping, for ISO Modula-2. *)

MODULE Sutherland_Hodgman_Task;

IMPORT STextIO, SRealIO;
IMPORT TextIO, RealIO;
IMPORT IOChan, StreamFile;

TYPE PlanePoint =
     RECORD
       x : REAL;
       y : REAL;
     END;

     PlaneEdge =
     RECORD
       pt0 : PlanePoint;        (* The start point. *)
       pt1 : PlanePoint;        (* The end point. *)
     END;

PROCEDURE evaluate_line (x1, y1, x2, y2, x : REAL) : REAL;
  VAR dy, dx, slope, intercept : REAL;
BEGIN
  dy := y2 - y1;
  dx := x2 - x1;
  slope := dy / dx;
  intercept := ((dx * y1) - (dy * x1)) / dx;
  RETURN (slope * x) + intercept
END evaluate_line;

PROCEDURE intersection_of_lines
            (x1, y1, x2, y2, x3, y3, x4, y4 : REAL) : PlanePoint;
  VAR intersection : PlanePoint;
      denominator, xnumerator, ynumerator : REAL;
      x1y2_y1x2, x3y4_y3x4 : REAL;
BEGIN
  IF x1 = x2 THEN
    intersection.x := x1;
    intersection.y := evaluate_line (x3, y3, x4, y4, x1);
  ELSIF x3 = x4 THEN
    intersection.x := x3;
    intersection.y := evaluate_line (x1, y1, x2, y2, x3);
  ELSE
    denominator := ((x1 - x2) * (y3 - y4)) - ((y1 - y2) * (x3 - x4));
    x1y2_y1x2 := (x1 * y2) - (y1 * x2);
    x3y4_y3x4 := (x3 * y4) - (y3 * x4);
    xnumerator := (x1y2_y1x2 * (x3 - x4)) - ((x1 - x2) * x3y4_y3x4);
    ynumerator := (x1y2_y1x2 * (y3 - y4)) - ((y1 - y2) * x3y4_y3x4);
    intersection.x := xnumerator / denominator;
    intersection.y := ynumerator / denominator;
  END;
  RETURN intersection;
END intersection_of_lines;

PROCEDURE intersection_of_edges
            (e1, e2 : PlaneEdge) : PlanePoint;
BEGIN
  RETURN intersection_of_lines (e1.pt0.x, e1.pt0.y,
                                e1.pt1.x, e1.pt1.y,
                                e2.pt0.x, e2.pt0.y,
                                e2.pt1.x, e2.pt1.y);
END intersection_of_edges;

PROCEDURE point_is_left_of_edge
            (pt   : PlanePoint;
             edge : PlaneEdge) : BOOLEAN;
  VAR x, y, x1, y1, x2, y2, op : REAL;
BEGIN
  x := pt.x;
  y := pt.y;
  x1 := edge.pt0.x;
  y1 := edge.pt0.y;
  x2 := edge.pt1.x;
  y2 := edge.pt1.y;

  (* Outer product of the vectors (x1,y1)-->(x,y) and
     (x1,y1)-->(x2,y2). *)
  op := ((x - x1) * (y2 - y1)) - ((x2 - x1) * (y - y1));

  RETURN (op < 0.0);
END point_is_left_of_edge;

PROCEDURE clip_subject_edge
            (subject_edge : PlaneEdge;
             clip_edge    : PlaneEdge;
             VAR n        : CARDINAL;
             VAR points   : ARRAY OF PlanePoint);
  VAR s1, s2 : PlanePoint;
      s2_is_inside, s1_is_inside : BOOLEAN;
BEGIN
  s1 := subject_edge.pt0;
  s2 := subject_edge.pt1;
  s2_is_inside := point_is_left_of_edge (s2, clip_edge);
  s1_is_inside := point_is_left_of_edge (s1, clip_edge);
  IF s2_is_inside THEN
    IF s1_is_inside THEN
      points[n] := s2;
      n := n + 1;
    ELSE
      points[n] := intersection_of_edges (subject_edge, clip_edge);
      n := n + 1;
      points[n] := s2;
      n := n + 1;
    END;
  ELSIF s1_is_inside THEN
    points[n] := intersection_of_edges (subject_edge, clip_edge);
    n := n + 1;
  END;
END clip_subject_edge;

PROCEDURE for_each_subject_edge
            (nsubject       : CARDINAL;
             subject_points : ARRAY OF PlanePoint;
             clip_edge      : PlaneEdge;
             VAR n          : CARDINAL;
             VAR points     : ARRAY OF PlanePoint);
  VAR subject_edge : PlaneEdge;
      i, j         : CARDINAL;
BEGIN
  n := 0;
  FOR i := 0 TO nsubject - 1 DO
    IF i = 0 THEN
      j := nsubject - 1;
    ELSE
      j := i - 1;
    END;
    subject_edge.pt1 := subject_points[i];
    subject_edge.pt0 := subject_points[j];
    clip_subject_edge (subject_edge, clip_edge, n, points);
  END;
END for_each_subject_edge;

PROCEDURE clip (VAR nsubject       : CARDINAL;
                VAR subject_points : ARRAY OF PlanePoint;
                nclip              : CARDINAL;
                clip_points        : ARRAY OF PlanePoint;
                VAR workspace      : ARRAY OF PlanePoint);
  VAR clip_edge   : PlaneEdge;
      i, j, nwork : CARDINAL;
BEGIN
  FOR i := 0 TO nclip - 1 DO
    IF i = 0 THEN
      j := nclip - 1;
    ELSE
      j := i - 1;
    END;
    clip_edge.pt1 := clip_points[i];
    clip_edge.pt0 := clip_points[j];
    for_each_subject_edge (nsubject, subject_points, clip_edge,
                           nwork, workspace);
    FOR j := 0 TO nwork - 1 DO
      subject_points[j] := workspace[j];
    END;
    nsubject := nwork;
  END;
END clip;

PROCEDURE set_point
            (VAR points : ARRAY OF PlanePoint;
             i          : CARDINAL;
             x, y       : REAL);
BEGIN
  points[i].x := x;
  points[i].y := y;
END set_point;

PROCEDURE write_polygon
            (cid        : IOChan.ChanId;
             npoly      : CARDINAL;
             polygon    : ARRAY OF PlanePoint;
             line_color : ARRAY OF CHAR;
             fill_color : ARRAY OF CHAR);
  VAR i : CARDINAL;
BEGIN
  RealIO.WriteReal (cid, polygon[0].x, 10);
  TextIO.WriteString (cid, ' ');
  RealIO.WriteReal (cid, polygon[0].y, 10);
  TextIO.WriteString (cid, ' moveto');
  TextIO.WriteLn (cid);
  FOR i := 1 TO npoly - 1 DO
    RealIO.WriteReal (cid, polygon[i].x, 10);
    TextIO.WriteString (cid, ' ');
    RealIO.WriteReal (cid, polygon[i].y, 10);
    TextIO.WriteString (cid, ' lineto');
    TextIO.WriteLn (cid);
  END;
  TextIO.WriteString (cid, 'closepath');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, line_color);
  TextIO.WriteString (cid, ' setrgbcolor');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, 'gsave');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, fill_color);
  TextIO.WriteString (cid, ' setrgbcolor');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, 'fill');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, 'grestore');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, 'stroke');
  TextIO.WriteLn (cid);
END write_polygon;

PROCEDURE write_eps
            (cid             : IOChan.ChanId;
             nsubject        : CARDINAL;
             subject_polygon : ARRAY OF PlanePoint;
             nclip           : CARDINAL;
             clip_polygon    : ARRAY OF PlanePoint;
             nresult         : CARDINAL;
             result_polygon  : ARRAY OF PlanePoint);
BEGIN
  TextIO.WriteString (cid, '%!PS-Adobe-3.0 EPSF-3.0');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, '%%BoundingBox: 40 40 360 360');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, '0 setlinewidth');
  TextIO.WriteLn (cid);
  write_polygon (cid, nclip, clip_polygon,
                 '.5 0 0', '1 .7 .7');
  write_polygon (cid, nsubject, subject_polygon,
                 '0 .2 .5', '.4 .7 1');
  TextIO.WriteString (cid, '2 setlinewidth');
  TextIO.WriteLn (cid);
  TextIO.WriteString (cid, '[10 8] 0 setdash');
  TextIO.WriteLn (cid);
  write_polygon (cid, nresult, result_polygon,
                 '.5 0 .5', '.7 .3 .8');
  TextIO.WriteString (cid, '%%EOF');
  TextIO.WriteLn (cid);
END write_eps;

PROCEDURE write_eps_to_file
            (filename        : ARRAY OF CHAR;
             nsubject        : CARDINAL;
             subject_polygon : ARRAY OF PlanePoint;
             nclip           : CARDINAL;
             clip_polygon    : ARRAY OF PlanePoint;
             nresult         : CARDINAL;
             result_polygon  : ARRAY OF PlanePoint);
VAR cid          : IOChan.ChanId;
    open_results : StreamFile.OpenResults;
BEGIN
  StreamFile.Open (cid, filename,
                   StreamFile.write,
                   open_results);
  write_eps (cid,
             nsubject, subject_polygon,
             nclip, clip_polygon,
             nresult, result_polygon);
  StreamFile.Close (cid);
END write_eps_to_file;

CONST NMax = 100;

VAR subject_polygon : ARRAY [0 .. NMax - 1] OF PlanePoint;
    clip_polygon    : ARRAY [0 .. NMax - 1] OF PlanePoint;
    workspace       : ARRAY [0 .. NMax - 1] OF PlanePoint;
    result_polygon  : ARRAY [0 .. NMax - 1] OF PlanePoint;
    nsubject, nclip, nresult, i : CARDINAL;

BEGIN
  nsubject := 9;
  set_point (subject_polygon, 0, 50.0, 150.0);
  set_point (subject_polygon, 1, 200.0, 50.0);
  set_point (subject_polygon, 2, 350.0, 150.0);
  set_point (subject_polygon, 3, 350.0, 300.0);
  set_point (subject_polygon, 4, 250.0, 300.0);
  set_point (subject_polygon, 5, 200.0, 250.0);
  set_point (subject_polygon, 6, 150.0, 350.0);
  set_point (subject_polygon, 7, 100.0, 250.0);
  set_point (subject_polygon, 8, 100.0, 200.0);

  nclip := 4;
  set_point (clip_polygon, 0, 100.0, 100.0);
  set_point (clip_polygon, 1, 300.0, 100.0);
  set_point (clip_polygon, 2, 300.0, 300.0);
  set_point (clip_polygon, 3, 100.0, 300.0);

  FOR i := 0 TO nsubject - 1 DO
    result_polygon[i] := subject_polygon[i];
  END;
  nresult := nsubject;

  clip (nresult, result_polygon, nclip, clip_polygon,
        workspace);

  FOR i := 0 TO nsubject - 1 DO
    STextIO.WriteString ('(');
    SRealIO.WriteReal (result_polygon[i].x, 8);
    STextIO.WriteString (', ');
    SRealIO.WriteReal (result_polygon[i].y, 8);
    STextIO.WriteString (')');
    STextIO.WriteLn;
  END;

  write_eps_to_file ('sutherland-hodgman.eps',
                     nsubject, subject_polygon,
                     nclip, clip_polygon,
                     nresult, result_polygon);
  STextIO.WriteString ('Wrote sutherland-hodgman.eps');
  STextIO.WriteLn;
END Sutherland_Hodgman_Task.
Output:
gm2 sutherland_hodgman_task.mod && ./a.out
(100.0000, 116.6667)
(125.0000, 100.0000)
(275.0000, 100.0000)
(300.0000, 116.6667)
(300.0000, 300.0000)
(250.0000, 300.0000)
(200.0000, 250.0000)
(175.0000, 300.0000)
(125.0000, 300.0000)
Wrote sutherland-hodgman.eps

Sutherland-Hodgman task polygons from Modula-2.

Nim[edit]

Translation of: D
import sequtils, strformat

type
  Vec2 = tuple[x, y: float]
  Edge = tuple[p, q: Vec2]
  Polygon = seq[Vec2]


func `-`(a, b: Vec2): Vec2 = (a.x - b.x, a.y - b.y)

func cross(a, b: Vec2): float = a.x * b.y - a.y * b.x

func isInside(p: Vec2; edge: Edge): bool =
  (edge.q.x - edge.p.x) * (p.y - edge.p.y) > (edge.q.y - edge.p.y) * (p.x - edge.p.x)


func intersection(sEdge, cEdge: Edge): Vec2 =
  let
    dc = cEdge.p - cEdge.q
    dp = sEdge.p - sEdge.q
    n1 = cEdge.p.cross(cEdge.q)
    n2 = sEdge.p.cross(sEdge.q)
    n3 = 1 / dc.cross(dp)
  result = ((n1 * dp.x - n2 * dc.x) * n3, (n1 * dp.y - n2 * dc.y) * n3)


func edges(poly: Polygon): seq[Edge] =
  (poly[^1] & poly).zip(poly)


func clip(subjectPolygon, clipPolygon: Polygon): Polygon =
  assert subjectPolygon.len > 1
  assert clipPolygon.len > 1

  result = subjectPolygon
  for clipEdge in clipPolygon.edges:
    let inputList = move(result)
    result.reset()
    for inEdge in inputList.edges:
      if inEdge.q.isInside(clipEdge):
        if not inEdge.p.isInside(clipEdge):
          result.add intersection(inEdge, clipEdge)
        result.add inEdge.q
      elif inEdge.p.isInside(clipEdge):
        result.add intersection(inEdge, clipEdge)


proc saveEpsImage(filename: string; subject, clip, clipped: Polygon) =
  let eps = open(filename, fmWrite)
  eps.write "%%!PS-Adobe-3.0\n%%%%BoundingBox: 40 40 360 360\n/l {lineto} def\n/m {moveto} def\n",
            "/s {setrgbcolor} def\n/c {closepath} def\n/gs {fill grestore stroke} def\n"

  eps.write &"0 setlinewidth {clip[0].x} {clip[0].y} m "
  for i in 1..clip.high:
    eps.write &"{clip[i].x} {clip[i].y} l "
  eps.writeLine "c 0.5 0 0 s gsave 1 0.7 0.7 s gs"

  eps.write &"{subject[0].x} {subject[0].y} m "
  for i in 1..subject.high:
    eps.write &"{subject[i].x} {subject[i].y} l "
  eps.writeLine "c 0 0.2 0.5 s gsave 0.4 0.7 1 s gs"

  eps.write &"2 setlinewidth [10 8] 0 setdash {clipped[0].x} {clipped[0].y} m "
  for i in 1..clipped.high:
    eps.write &"{clipped[i].x} {clipped[i].y} l "
  eps.writeLine &"c 0.5 0 0.5 s gsave 0.7 0.3 0.8 s gs"

  eps.writeLine "%%%%EOF"
  eps.close()
  echo &"File “{filename}” written."


when isMainModule:

  let
    subjectPolygon = @[(50.0, 150.0), (200.0, 50.0), (350.0, 150.0),
                      (350.0, 300.0), (250.0, 300.0), (200.0, 250.0),
                      (150.0, 350.0), (100.0, 250.0), (100.0, 200.0)]
    clippingPolygon = @[(100.0, 100.0), (300.0, 100.0), (300.0, 300.0), (100.0, 300.0)]

    clipped = subjectPolygon.clip(clippingPolygon)

  for point in clipped:
    echo &"({point.x:.3f}, {point.y:.3f})"
  saveEpsImage("sutherland_hodgman_clipping_out.eps", subjectPolygon, clippingPolygon, clipped)
Output:
(100.000, 116.667)
(125.000, 100.000)
(275.000, 100.000)
(300.000, 116.667)
(300.000, 300.000)
(250.000, 300.000)
(200.000, 250.000)
(175.000, 300.000)
(125.000, 300.000)
(100.000, 250.000)
File “sutherland_hodgman_clipping_out.eps” written.

OCaml[edit]

let is_inside (x,y) ((ax,ay), (bx,by)) =
  (bx -. ax) *. (y -. ay) > (by -. ay) *. (x -. ax)

let intersection (sx,sy) (ex,ey) ((ax,ay), (bx,by)) =
  let dc_x, dc_y = (ax -. bx, ay -. by) in
  let dp_x, dp_y = (sx -. ex, sy -. ey) in
  let n1 = ax *. by -. ay *. bx in
  let n2 = sx *. ey -. sy *. ex in
  let n3 = 1.0 /. (dc_x *. dp_y -. dc_y *. dp_x) in
  ((n1 *. dp_x -. n2 *. dc_x) *. n3,
   (n1 *. dp_y -. n2 *. dc_y) *. n3)

let last lst = List.hd (List.rev lst)

let polygon_iter_edges poly f init =
  if poly = [] then init else
    let p0 = List.hd poly in
    let rec aux acc = function
      | p1 :: p2 :: tl -> aux (f (p1, p2) acc) (p2 :: tl)
      | p :: [] -> f (p, p0) acc
      | [] -> acc
    in
    aux init poly

let poly_clip subject_polygon clip_polygon =
  polygon_iter_edges clip_polygon (fun clip_edge input_list ->
    fst (
      List.fold_left (fun (out, s) e ->

        match (is_inside e clip_edge), (is_inside s clip_edge) with
        | true, false -> (e :: (intersection s e clip_edge) :: out), e
        | true, true -> (e :: out), e
        | false, true -> ((intersection s e clip_edge) :: out), e
        | false, false -> (out, e)

      ) ([], last input_list) input_list)

  ) subject_polygon

let () =
  let subject_polygon =
    [ ( 50.0, 150.0); (200.0,  50.0); (350.0, 150.0);
      (350.0, 300.0); (250.0, 300.0); (200.0, 250.0);
      (150.0, 350.0); (100.0, 250.0); (100.0, 200.0); ] in

  let clip_polygon =
    [ (100.0, 100.0); (300.0, 100.0); (300.0, 300.0); (100.0, 300.0) ] in

  List.iter (fun (x,y) ->
      Printf.printf " (%g, %g)\n" x y;
    ) (poly_clip subject_polygon clip_polygon)
Output:
 (100, 116.667)
 (125, 100)
 (275, 100)
 (300, 116.667)
 (300, 300)
 (250, 300)
 (200, 250)
 (175, 300)
 (125, 300)
 (100, 250)

We can display the result in a window using the Graphics module:

let subject_polygon =
  [ ( 50.0, 150.0); (200.0,  50.0); (350.0, 150.0);
    (350.0, 300.0); (250.0, 300.0); (200.0, 250.0);
    (150.0, 350.0); (100.0, 250.0); (100.0, 200.0); ]

let clip_polygon =
  [ (100.0, 100.0); (300.0, 100.0); (300.0, 300.0); (100.0, 300.0) ]

let () =
  Graphics.open_graph " 400x400";
  let to_grid poly =
    let round x = int_of_float (floor (x +. 0.5)) in
    Array.map
      (fun (x, y) -> (round x, round y))
      (Array.of_list poly)
  in
  let draw_poly fill stroke poly =
    let p = to_grid poly in
    Graphics.set_color fill;
    Graphics.fill_poly p;
    Graphics.set_color stroke;
    Graphics.draw_poly p;
  in
  draw_poly Graphics.red Graphics.blue subject_polygon;
  draw_poly Graphics.cyan Graphics.blue clip_polygon;
  draw_poly Graphics.magenta Graphics.blue (poly_clip subject_polygon clip_polygon);
  let _ = Graphics.wait_next_event [Graphics.Button_down; Graphics.Key_pressed] in
  Graphics.close_graph ()

SuthHodgClip OCaml.png

Perl[edit]

Translation of: Raku
use strict;
use warnings;

sub intersection {
    my($L11, $L12, $L21, $L22) = @_;
    my ($d1x, $d1y) = ($$L11[0] - $$L12[0], $$L11[1] - $$L12[1]);
    my ($d2x, $d2y) = ($$L21[0] - $$L22[0], $$L21[1] - $$L22[1]);
    my $n1 = $$L11[0] * $$L12[1] - $$L11[1] * $$L12[0];
    my $n2 = $$L21[0] * $$L22[1] - $$L21[1] * $$L22[0];
    my $n3 = 1 / ($d1x * $d2y - $d2x * $d1y);
    [($n1 * $d2x - $n2 * $d1x) * $n3, ($n1 * $d2y - $n2 * $d1y) * $n3]
}

sub is_inside {
    my($p1, $p2, $p3) = @_;
    ($$p2[0] - $$p1[0]) * ($$p3[1] - $$p1[1]) > ($$p2[1] - $$p1[1]) * ($$p3[0] - $$p1[0])
}

sub sutherland_hodgman {
    my($polygon, $clip) = @_;
    my @output = @$polygon;
    my $clip_point1 = $$clip[-1];
    for my $clip_point2 (@$clip) {
        my @input = @output;
        @output = ();
        my $start = $input[-1];
        for my $end (@input) {
            if (is_inside($clip_point1, $clip_point2, $end)) {
                push @output, intersection($clip_point1, $clip_point2, $start, $end)
                  unless is_inside($clip_point1, $clip_point2, $start);
                push @output, $end;
            } elsif (is_inside($clip_point1, $clip_point2, $start)) {
                push @output, intersection($clip_point1, $clip_point2, $start, $end);
            }
            $start = $end;
        }
        $clip_point1 = $clip_point2;
    }
    @output
}

my @polygon = ([50,  150], [200, 50],  [350, 150], [350, 300], [250, 300],
              [200, 250], [150, 350], [100, 250], [100, 200]);

my @clip    = ([100, 100], [300, 100], [300, 300], [100, 300]);

my @clipped = sutherland_hodgman(\@polygon, \@clip);

print "Clipped polygon:\n";
print '(' . join(' ', @$_) . ') ' for @clipped;
Output:
Clipped polygon:
(100 116.666666666667) (125 100) (275 100) (300 116.666666666667) (300 300) (250 300) (200 250) (175 300) (125 300) (100 250)

Phix[edit]

Library: Phix/pGUI
Library: Phix/online

You can run this online here.

--
-- demo\rosetta\Sutherland_Hodgman_polygon_clipping.exw
-- ====================================================
--
with javascript_semantics
enum X,Y

function inside(sequence cp1, cp2, p)
    return (cp2[X]-cp1[X])*(p[Y]-cp1[Y])>(cp2[Y]-cp1[Y])*(p[X]-cp1[X])
end function

function intersect(sequence cp1, cp2, s, e)
    atom {dcx,dcy} = {cp1[X]-cp2[X],cp1[Y]-cp2[Y]},
         {dpx,dpy} = {s[X]-e[X],s[Y]-e[Y]},
         n1 = cp1[X]*cp2[Y]-cp1[Y]*cp2[X],
         n2 = s[X]*e[Y]-s[Y]*e[X],
         n3 = 1/(dcx*dpy-dcy*dpx)
    return {(n1*dpx-n2*dcx)*n3,(n1*dpy-n2*dcy)*n3}
end function

function sutherland_hodgman(sequence subjectPolygon, clipPolygon)
    sequence cp1 = clipPolygon[$],
      outputList = subjectPolygon
    for i=1 to length(clipPolygon) do
        sequence cp2 = clipPolygon[i],
           inputList = outputList,
                   s = inputList[$]
        outputList = {}
        for j=1 to length(inputList) do
            sequence e = inputList[j]
            if inside(cp1,cp2,e) then
                if not inside(cp1,cp2,s) then
                    outputList = append(outputList,intersect(cp1,cp2,s,e))
                end if
                outputList = append(outputList,e)
            elsif inside(cp1,cp2,s) then
                outputList = append(outputList,intersect(cp1,cp2,s,e))
            end if
            s = e
        end for
        cp1 = cp2
    end for
    return outputList
end function

constant subjectPolygon = {{50, 150}, {200, 50}, {350, 150}, {350, 300},
                           {250, 300}, {200, 250}, {150, 350}, {100, 250},
                           {100, 200}},
         clipPolygon = {{100, 100}, {300, 100}, {300, 300}, {100, 300}}

sequence clippedPolygon = sutherland_hodgman(subjectPolygon,clipPolygon)

include pGUI.e

Ihandle dlg, canvas
cdCanvas cddbuffer, cdcanvas

procedure draw_poly(sequence poly)
    cdCanvasBegin(cddbuffer,CD_FILL)
    for i=1 to length(poly) do
        atom {x,y} = poly[i]
        cdCanvasVertex(cddbuffer,x,y)
    end for
    cdCanvasEnd(cddbuffer)
end procedure

function redraw_cb(Ihandle /*ih*/)
    cdCanvasActivate(cddbuffer)
    cdCanvasClear(cddbuffer)
    cdCanvasSetForeground(cddbuffer, CD_CYAN)
    draw_poly(subjectPolygon)
    cdCanvasSetForeground(cddbuffer, CD_MAGENTA)
    draw_poly(clipPolygon)
    cdCanvasSetForeground(cddbuffer, CD_ORANGE)
    draw_poly(clippedPolygon)
    cdCanvasFlush(cddbuffer)
    return IUP_DEFAULT
end function

function map_cb(Ihandle ih)
    cdcanvas = cdCreateCanvas(CD_IUP, ih)
    cddbuffer = cdCreateCanvas(CD_DBUFFER, cdcanvas)
    cdCanvasSetBackground(cddbuffer, CD_WHITE)
    cdCanvasSetForeground(cddbuffer, CD_GRAY)
    return IUP_DEFAULT
end function

procedure main()
    IupOpen()
    canvas = IupCanvas("RASTERSIZE=400x400")
    IupSetCallbacks(canvas, {"MAP_CB", Icallback("map_cb"),
                             "ACTION", Icallback("redraw_cb")})
    dlg = IupDialog(canvas, "RESIZE=NO")
    IupSetAttribute(dlg, "TITLE", "Sutherland-Hodgman polygon clipping")
    IupShow(dlg)
    if platform()!=JS then
        IupMainLoop()
        IupClose()
    end if
end procedure

main()

PHP[edit]

<?php
function clip ($subjectPolygon, $clipPolygon) {

    function inside ($p, $cp1, $cp2) {
        return ($cp2[0]-$cp1[0])*($p[1]-$cp1[1]) > ($cp2[1]-$cp1[1])*($p[0]-$cp1[0]);
    }
    
    function intersection ($cp1, $cp2, $e, $s) {
        $dc = [ $cp1[0] - $cp2[0], $cp1[1] - $cp2[1] ];
        $dp = [ $s[0] - $e[0], $s[1] - $e[1] ];
        $n1 = $cp1[0] * $cp2[1] - $cp1[1] * $cp2[0];
        $n2 = $s[0] * $e[1] - $s[1] * $e[0];
        $n3 = 1.0 / ($dc[0] * $dp[1] - $dc[1] * $dp[0]);

        return [($n1*$dp[0] - $n2*$dc[0]) * $n3, ($n1*$dp[1] - $n2*$dc[1]) * $n3];
    }
    
    $outputList = $subjectPolygon;
    $cp1 = end($clipPolygon);
    foreach ($clipPolygon as $cp2) {
        $inputList = $outputList;
        $outputList = [];
        $s = end($inputList);
        foreach ($inputList as $e) {
            if (inside($e, $cp1, $cp2)) {
                if (!inside($s, $cp1, $cp2)) {
                    $outputList[] = intersection($cp1, $cp2, $e, $s);
                }
                $outputList[] = $e;
            }
            else if (inside($s, $cp1, $cp2)) {
                $outputList[] = intersection($cp1, $cp2, $e, $s);
            }
            $s = $e;
        }
        $cp1 = $cp2;
    }
    return $outputList;
}

$subjectPolygon = [[50, 150], [200, 50], [350, 150], [350, 300], [250, 300], [200, 250], [150, 350], [100, 250], [100, 200]];
$clipPolygon = [[100, 100], [300, 100], [300, 300], [100, 300]];
$clippedPolygon = clip($subjectPolygon, $clipPolygon);

echo json_encode($clippedPolygon);
echo "\n";
?>

PureBasic[edit]

Translation of: Go
Structure point_f
  x.f
  y.f
EndStructure

Procedure isInside(*p.point_f, *cp1.point_f, *cp2.point_f)  
  If (*cp2\x - *cp1\x) * (*p\y - *cp1\y) > (*cp2\y - *cp1\y) * (*p\x - *cp1\x)
    ProcedureReturn 1
  EndIf 
EndProcedure

Procedure intersection(*cp1.point_f, *cp2.point_f, *s.point_f, *e.point_f, *newPoint.point_f)
  Protected.point_f dc, dp
  Protected.f n1, n2, n3
  dc\x = *cp1\x - *cp2\x: dc\y = *cp1\y - *cp2\y
  dp\x = *s\x - *e\x: dp\y = *s\y - *e\y
  n1 = *cp1\x * *cp2\y - *cp1\y * *cp2\x
  n2 = *s\x * *e\y - *s\y * *e\x
  n3 = 1 / (dc\x * dp\y - dc\y * dp\x)
  *newPoint\x = (n1 * dp\x - n2 * dc\x) * n3: *newPoint\y = (n1 * dp\y - n2 * dc\y) * n3
EndProcedure

Procedure clip(List vPolygon.point_f(), List vClippedBy.point_f(), List vClippedPolygon.point_f())
  Protected.point_f cp1, cp2, s, e, newPoint
  CopyList(vPolygon(), vClippedPolygon())
  If LastElement(vClippedBy())
    cp1 = vClippedBy()
    
    NewList vPreClipped.point_f()
    ForEach vClippedBy()
      cp2 = vClippedBy()
      CopyList(vClippedPolygon(), vPreClipped())
      ClearList(vClippedPolygon())
      If LastElement(vPreClipped())
        s = vPreClipped()
        ForEach vPreClipped()
          e = vPreClipped()
          If isInside(e, cp1, cp2)
            If Not isInside(s, cp1, cp2)
              intersection(cp1, cp2, s, e, newPoint)
              AddElement(vClippedPolygon()): vClippedPolygon() = newPoint
            EndIf 
            AddElement(vClippedPolygon()): vClippedPolygon() = e
          ElseIf isInside(s, cp1, cp2)
            intersection(cp1, cp2, s, e, newPoint)
            AddElement(vClippedPolygon()): vClippedPolygon() = newPoint
          EndIf 
          s = e
        Next
      EndIf 
      cp1 = cp2
    Next 
  EndIf 
EndProcedure

DataSection
  Data.f 50,150, 200,50, 350,150, 350,300, 250,300, 200,250, 150,350, 100,250, 100,200 ;subjectPolygon's vertices (x,y)
  Data.f 100,100, 300,100, 300,300, 100,300 ;clipPolygon's vertices (x,y)
EndDataSection

NewList subjectPolygon.point_f()
For i = 1 To 9
  AddElement(subjectPolygon())
  Read.f subjectPolygon()\x
  Read.f subjectPolygon()\y
Next 

NewList clipPolygon.point_f()
For i = 1 To 4
  AddElement(clipPolygon())
  Read.f clipPolygon()\x
  Read.f clipPolygon()\y
Next 

NewList newPolygon.point_f()
clip(subjectPolygon(), clipPolygon(), newPolygon())
If OpenConsole()
  ForEach newPolygon()
    PrintN("(" + StrF(newPolygon()\x, 2) + ", " + StrF(newPolygon()\y, 2) + ")")
  Next
   
  Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()
  CloseConsole()
EndIf
Output:
(100.00, 116.67)
(125.00, 100.00)
(275.00, 100.00)
(300.00, 116.67)
(300.00, 300.00)
(250.00, 300.00)
(200.00, 250.00)
(175.00, 300.00)
(125.00, 300.00)
(100.00, 250.00)

Python[edit]

def clip(subjectPolygon, clipPolygon):
   def inside(p):
      return(cp2[0]-cp1[0])*(p[1]-cp1[1]) > (cp2[1]-cp1[1])*(p[0]-cp1[0])
      
   def computeIntersection():
      dc = [ cp1[0] - cp2[0], cp1[1] - cp2[1] ]
      dp = [ s[0] - e[0], s[1] - e[1] ]
      n1 = cp1[0] * cp2[1] - cp1[1] * cp2[0]
      n2 = s[0] * e[1] - s[1] * e[0] 
      n3 = 1.0 / (dc[0] * dp[1] - dc[1] * dp[0])
      return [(n1*dp[0] - n2*dc[0]) * n3, (n1*dp[1] - n2*dc[1]) * n3]

   outputList = subjectPolygon
   cp1 = clipPolygon[-1]
   
   for clipVertex in clipPolygon:
      cp2 = clipVertex
      inputList = outputList
      outputList = []
      s = inputList[-1]

      for subjectVertex in inputList:
         e = subjectVertex
         if inside(e):
            if not inside(s):
               outputList.append(computeIntersection())
            outputList.append(e)
         elif inside(s):
            outputList.append(computeIntersection())
         s = e
      cp1 = cp2
   return(outputList)

Racket[edit]

Shameless rewrite of haskell version.

#lang racket

(module sutherland-hodgman racket
  (provide clip-to)
  (provide make-edges)
  (provide (struct-out point))

  (struct point (x y) #:transparent)
  (struct edge (p1 p2) #:transparent)
  (struct polygon (points edges) #:transparent)

  (define (make-edges points)
    (let ([points-shifted
	   (match points
	     [(list a b ...) (append b (list a))])])
      (map edge points points-shifted)))

  (define (is-point-left? pt ln)
    (match-let ([(point x y) pt]
                [(edge (point px py) (point qx qy)) ln])
               (>= (* (- qx px) (- y py))
                   (* (- qy py) (- x px)))))

  ;; Return the intersection of two lines
  (define (isect-lines l1 l2)
    (match-let ([(edge (point x1 y1) (point x2 y2)) l1]
                [(edge (point x3 y3) (point x4 y4)) l2])
               (let* ([r (- (* x1 y2) (* y1 x2))] [s (- (* x3 y4) (* y3 x4))]
                      [t (- x1 x2)] [u (- y3 y4)] [v (- y1 y2)] [w (- x3 x4)]
                      [d (- (* t u) (* v w))])
                 (point (/ (- (* r w) (* t s)) d)
                        (/ (- (* r u) (* v s)) d)))))
  
  ;; Intersect the line segment (p0,p1) with the clipping line's left halfspace, 
  ;; returning the point closest to p1.  In the special case where p0 lies outside  
  ;; the halfspace and p1 lies inside we return both the intersection point and p1.  
  ;; This ensures we will have the necessary segment along the clipping line.

  (define (intersect segment clip-line)
    (define (isect) (isect-lines segment clip-line))

    (match-let ([(edge p0 p1) segment])
               (match/values (values (is-point-left? p0 clip-line) (is-point-left? p1 clip-line))
                             [(#f #f) '()]
                             [(#f #t) (list (isect) p1)]
                             [(#t #f) (list (isect))]
                             [(#t #t) (list p1)])))

  ;; Intersect the polygon with the clipping line's left halfspace
  (define (isect-polygon poly-edges clip-line)
    (for/fold ([p '()]) ([e  poly-edges])
      (append p (intersect e clip-line))))

  ;; Intersect a subject polygon with a clipping polygon.  The latter is assumed to be convex.
  (define (clip-to sp-pts cp-edges)
    (for/fold ([out-poly sp-pts]) ([clip-line cp-edges])
      (isect-polygon (make-edges out-poly) clip-line))))

Testing code (Couldn't find a way to attach image with polygons)

(require racket/gui)  
(require 'sutherland-hodgman)

(define (make-points pt-list)
    (for/list ([p pt-list])
      (make-object point% (point-x p) (point-y p))))

(define subject-poly-points 
    (list (point 50 150)  (point 200 50)  (point 350 150) 
          (point 350 300) (point 250 300) (point 200 250) 
          (point 150 350) (point 100 250) (point 100 200)))

(define clip-poly-points
    (list (point 100 100)
          (point 300 100)
          (point 300 300)
          (point 100 300)))

(define clip-poly-edges 
    (make-edges clip-poly-points))

(define (run)
  (let* ([frame (new frame% [label "Sutherland-Hodgman racket demo"]
		     [width 320]
		     [height 320])]
	 [canvas (new canvas% [parent frame])]
	 [dc (send canvas get-dc)]
         [clipped-poly (clip-to subject-poly-points clip-poly-edges)])
    
    (send frame show #t)
    (sleep/yield 1)

    (send dc set-pen (make-pen 
                        #:color (send the-color-database find-color "Blue")
                        #:width 3))
    (send dc draw-polygon (make-points subject-poly-points))
    (send dc set-pen (make-pen 
                        #:color (send the-color-database find-color "Red")
                        #:width 4
                        #:style 'long-dash))
    (send dc draw-polygon (make-points clip-poly-points))
    (send dc set-pen (make-pen 
                        #:color (send the-color-database find-color "Green")))
    (send dc set-brush (make-brush 
                        #:color (send the-color-database find-color "Green")
                        #:style 'solid))
    (send dc draw-polygon (make-points clipped-poly))
    clipped-poly))

(run)

Output:

(list
 (point 300 300)
 (point 250 300)
 (point 200 250)
 (point 175 300)
 (point 125 300)
 (point 100 250)
 (point 100 200)
 (point 100 200)
 (point 100 350/3)
 (point 125 100)
 (point 275 100)
 (point 300 350/3))

Raku[edit]

(formerly Perl 6)

Works with: Rakudo version 2019.11
Translation of: Sidef
sub intersection ($L11, $L12, $L21, $L22) {
    my ($Δ1x, $Δ1y) = $L11 »-« $L12;
    my ($Δ2x, $Δ2y) = $L21 »-« $L22;
    my $n1 = $L11[0] * $L12[1] - $L11[1] * $L12[0];
    my $n2 = $L21[0] * $L22[1] - $L21[1] * $L22[0];
    my $n3 = 1 / ($Δ1x * $Δ2y - $Δ2x * $Δ1y);
    (($n1 * $Δ2x - $n2 * $Δ1x) * $n3, ($n1 * $Δ2y - $n2 * $Δ1y) * $n3)
}


sub is-inside ($p1, $p2, $p3) {
    ($p2[0] - $p1[0]) * ($p3[1] - $p1[1]) > ($p2[1] - $p1[1]) * ($p3[0] - $p1[0])
}

sub sutherland-hodgman (@polygon, @clip) {
    my @output = @polygon;
    my $clip-point1 = @clip.tail;
    for @clip -> $clip-point2 {
        my @input = @output;
        @output = ();
        my $start = @input.tail;
        for @input -> $end {
            if is-inside($clip-point1, $clip-point2, $end) {
                @output.push: intersection($clip-point1, $clip-point2, $start, $end)
                  unless is-inside($clip-point1, $clip-point2, $start);
                @output.push: $end;
            } elsif is-inside($clip-point1, $clip-point2, $start) {
                @output.push: intersection($clip-point1, $clip-point2, $start, $end);
            }
            $start = $end;
        }
        $clip-point1 = $clip-point2;
    }
    @output
}

my @polygon = (50,  150), (200, 50),  (350, 150), (350, 300), (250, 300),
              (200, 250), (150, 350), (100, 250), (100, 200);

my @clip    = (100, 100), (300, 100), (300, 300), (100, 300);

my @clipped = sutherland-hodgman(@polygon, @clip);

say "Clipped polygon: ", @clipped;

# Output an SVG as well as it is easier to visualize
use SVG;
my $outfile = 'Sutherland-Hodgman-polygon-clipping-perl6.svg'.IO.open(:w);
$outfile.say: SVG.serialize(
    svg => [
        :400width, :400height,
        :rect[ :400width, :400height, :fill<white> ],
        :text[ :10x, :20y, "Polygon (blue)" ],
        :text[ :10x, :35y, "Clip port (green)" ],
        :text[ :10x, :50y, "Clipped polygon (red)" ],
        :polyline[ :points(@polygon.join: ','), :style<stroke:blue>,  :fill<blue>,  :opacity<.3> ],
        :polyline[ :points(   @clip.join: ','), :style<stroke:green>, :fill<green>, :opacity<.3> ],
        :polyline[ :points(@clipped.join: ','), :style<stroke:red>,   :fill<red>,   :opacity<.5> ],
    ],
);
Output:
Clipped polygon: [(100 116.666667) (125 100) (275 100) (300 116.666667) (300 300) (250 300) (200 250) (175 300) (125 300) (100 250)]

Also see output image: offsite SVG image

RATFOR[edit]

Translation of: ATS
Works with: ratfor77 version public domain 1.0
# Sutherland-Hodgman polygon clipping.
#
# On a POSIX platform, the program can be compiled with f2c and run
# somewhat as follows:
#
#    ratfor77 sutherland-hodgman.r > sutherland-hodgman.f
#    f2c -C sutherland-hodgman.f
#    cc sutherland-hodgman.c -lf2c
#    ./a.out
#
# With gfortran, a little differently:
#
#    ratfor77 sutherland-hodgman.r > sutherland-hodgman.f
#    gfortran -fcheck=all -std=legacy sutherland-hodgman.f
#    ./a.out

function evalln (x1, y1, x2, y2, x)
  #
  # Given the line (x1,y1)--(x2,y2), evaluate it at x.
  #
  implicit none
  real evalln
  real x1, y1, x2, y2, x
  real dy, dx, slope, xcept
  dy = y2 - y1
  dx = x2 - x1
  slope = dy / dx
  xcept = ((dx * y1) - (dy * x1)) / dx
  evalln = (slope * x) + xcept
end

subroutine xsctln (x1, y1, x2, y2, x3, y3, x4, y4, x, y)
  #
  # Given lines (x1,y1)--(x2,y2) and (x3,y3)--(x4,y4), find their
  # intersection (x,y).
  #
  implicit none
  real x1, y1, x2, y2, x3, y3, x4, y4, x, y
  real evalln
  real denom, xnumer, ynumer, xyyx12, xyyx34
  if (x1 == x2)
    {
      x = x1
      y = evalln (x3, y3, x4, y4, x1)
    }
  else if (x3 == x4)
    {
      x = x3
      y = evalln (x1, y1, x2, y2, x3)
    }
  else
    {
      denom = ((x1 - x2) * (y3 - y4)) - ((y1 - y2) * (x3 - x4))
      xyyx12 = (x1 * y2) - (y1 * x2)
      xyyx34 = (x3 * y4) - (y3 * x4)
      xnumer = (xyyx12 * (x3 - x4)) - ((x1 - x2) * xyyx34)
      ynumer = (xyyx12 * (y3 - y4)) - ((y1 - y2) * xyyx34)
      x = xnumer / denom
      y = ynumer / denom
    }
end

function ptleft (x, y, x1, y1, x2, y2)
  #
  # Is the point (x,y) left of the edge (x1,y1)--(x2,y2)?
  #
  implicit none
  logical ptleft
  real x, y, x1, y1, x2, y2
  ptleft = (((x - x1) * (y2 - y1)) - ((x2 - x1) * (y - y1)) < 0)
end

subroutine clpsbe (xs1, ys1, xs2, ys2, xc1, yc1, xc2, yc2, n, pts)
  #
  # Clip subject edge (xs1,ys1)--(xs2,ys2) with clip edge
  # (xc1,yc1)--(xc2,yc2).
  #
  implicit none
  real xs1, ys1, xs2, ys2, xc1, yc1, xc2, yc2
  integer n
  real pts(2,*), x, y
  logical ptleft, s2left, s1left
  s2left = ptleft (xs2, ys2, xc1, yc1, xc2, yc2)
  s1left = ptleft (xs1, ys1, xc1, yc1, xc2, yc2)
  if (s2left)
    {
      if (s1left)
        {
          n = n + 1
          pts(1,n) = xs2
          pts(2,n) = ys2
        }
      else
        {
          call xsctln (xs1, ys1, xs2, ys2, xc1, yc1, xc2, yc2, x, y)
          n = n + 1
          pts(1,n) = x
          pts(2,n) = y
          n = n + 1
          pts(1,n) = xs2
          pts(2,n) = ys2
        }
    }
  else if (s1left)
    {
      call xsctln (xs1, ys1, xs2, ys2, xc1, yc1, xc2, yc2, x, y)
      n = n + 1
      pts(1,n) = x
      pts(2,n) = y
    }
end

subroutine sublp (nsub, ptssub, xc1, yc1, xc2, yc2, n, pts)
  #
  # Loop over the subject points in ptssub, clipping the edges
  # therein. Produce a result in pts.
  #
  implicit none
  integer nsub, n
  real ptssub(2,*), pts(2,*)
  real xc1, yc1, xc2, yc2
  real xs1, ys1, xs2, ys2
  integer i, j
  for (i = 1; i <= nsub; i = i + 1)
    {
      xs2 = ptssub(1,i)
      ys2 = ptssub(2,i)
      j = i - 1
      if (j == 0) j = nsub
      xs1 = ptssub(1,j)
      ys1 = ptssub(2,j)
      call clpsbe (xs1, ys1, xs2, ys2, xc1, yc1, xc2, yc2, n, pts)
    }
end

subroutine clip (nsub, ptssub, nclp, ptsclp, ptswrk)
  #
  # Loop over the clip points in ptsclp, clipping the subject stored
  # in ptssub. Use ptswrk as workspace.
  #
  implicit none
  integer nsub, nclp
  real ptssub(2,*), ptsclp(2,*), ptswrk(2,*)
  integer i, j, nwrk
  real xc1, yc1, xc2, yc2
  for (i = 1; i <= nclp; i = i + 1)
    {
      xc2 = ptsclp(1,i)
      yc2 = ptsclp(2,i)
      j = i - 1
      if (j == 0) j = nclp
      xc1 = ptsclp(1,j)
      yc1 = ptsclp(2,j)
      nwrk = 0
      call sublp (nsub, ptssub, xc1, yc1, xc2, yc2, nwrk, ptswrk)

      # Copy the new subject over the old subject.
      for (j = 1; j <= nwrk; j = j + 1)
        {
          ptssub(1,j) = ptswrk(1,j)
          ptssub(2,j) = ptswrk(2,j)
        }
      nsub = nwrk
    }
end

subroutine wrtpts (eps, n, pts, linclr, filclr)
  #
  # Write a polygon as PostScript code.
  #
  implicit none
  character*10 linclr, filclr
  integer eps, n, i
  real pts(2,*)

10 format (F12.6, ' ', F12.6, ' moveto')
20 format (F12.6, ' ', F12.6, ' lineto')
30 format ('closepath')
40 format ('gsave')
50 format ('grestore')
60 format ('fill')
70 format ('stroke')
80 format (A10, ' setrgbcolor')

  write (eps, 10) pts(1,1), pts(2,1)
  for (i = 2; i <= n; i = i + 1)
    write (eps, 20) pts(1,i), pts(2,i)
  write (eps, 30)
  write (eps, 80) linclr
  write (eps, 40)
  write (eps, 80) filclr
  write (eps, 60)
  write (eps, 50)
  write (eps, 70)
end

subroutine wrteps (eps, nsub, ptssub, nclp, ptsclp, nres, ptsres)
  #
  # Write an Encapsulated PostScript file.
  #
  implicit none
  integer eps
  integer nsub, nclp, nres
  real ptssub(2,*), ptsclp(2,*), ptsres(2,*)

10 format ('%!PS-Adobe-3.0 EPSF-3.0')
20 format ('%%BoundingBox: 40 40 360 360')
30 format ('0 setlinewidth ')
40 format ('2 setlinewidth')
50 format ('[10 8] 0 setdash')
60 format ('%%EOF')

  write (eps, 10)
  write (eps, 20)
  write (eps, 30)
  call wrtpts (eps, nclp, ptsclp, '.5 0 0    ', '1 .7 .7   ')
  call wrtpts (eps, nsub, ptssub, '0 .2 .5   ', '.4 .7 1   ')
  write (eps, 40)
  write (eps, 50)
  call wrtpts (eps, nres, ptsres, '.5 0 .5   ', '.7 .3 .8  ')
  write (eps, 60)
end

define(NMAX,100)              # Maximum number of points in a polygon.
define(EPSF,9)                # Unit number for the EPS file.

program shclip
  implicit none
  integer nsub, nclp, nres
  real ptssub(2,NMAX), ptsclp(2,NMAX), ptsres(2,NMAX), ptswrk(2,NMAX)
  integer i
  integer eps

  nsub = 9
  ptssub(1,1) = 50; ptssub(2,1) = 150
  ptssub(1,2) = 200; ptssub(2,2) = 50
  ptssub(1,3) = 350; ptssub(2,3) = 150
  ptssub(1,4) = 350; ptssub(2,4) = 300
  ptssub(1,5) = 250; ptssub(2,5) = 300
  ptssub(1,6) = 200; ptssub(2,6) = 250
  ptssub(1,7) = 150; ptssub(2,7) = 350
  ptssub(1,8) = 100; ptssub(2,8) = 250
  ptssub(1,9) = 100; ptssub(2,9) = 200

  nclp = 4
  ptsclp(1,1) = 100; ptsclp(2,1) = 100
  ptsclp(1,2) = 300; ptsclp(2,2) = 100
  ptsclp(1,3) = 300; ptsclp(2,3) = 300
  ptsclp(1,4) = 100; ptsclp(2,4) = 300

  # Copy the subject points to the "result" array.
  for (i = 1; i <= nsub; i = i + 1)
    {
      ptsres(1,i) = ptssub(1,i)
      ptsres(2,i) = ptssub(2,i)
    }
  nres = nsub

  call clip (nres, ptsres, nclp, ptsclp, ptswrk)
  for (i = 1; i <= nres; i = i + 1)
    write (*, 1000) ptsres(1,i), ptsres(2,i)
1000 format ('(', F8.4, ', ', F8.4, ')')

  eps = EPSF
  open (unit=eps, file='sutherland-hodgman.eps')
  call wrteps (eps, nsub, ptssub, nclp, ptsclp, nres, ptsres)
  write (*, 1010)
1010 format ('Wrote sutherland-hodgman.eps')
end
Output:
(100.0000, 116.6667)
(125.0000, 100.0000)
(275.0000, 100.0000)
(300.0000, 116.6667)
(300.0000, 300.0000)
(250.0000, 300.0000)
(200.0000, 250.0000)
(175.0000, 300.0000)
(125.0000, 300.0000)
(100.0000, 250.0000)
Wrote sutherland-hodgman.eps

The polygons as generated by Ratfor.

Ruby[edit]

Translation of: Go
Point = Struct.new(:x,:y) do
  def to_s; "(#{x}, #{y})" end
end

def sutherland_hodgman(subjectPolygon, clipPolygon)
  # These inner functions reduce the argument passing to
  # "inside" and "intersection".
  cp1, cp2, s, e = nil
  inside = proc do |p|
    (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x)
  end
  intersection = proc do
    dcx, dcy = cp1.x-cp2.x, cp1.y-cp2.y
    dpx, dpy = s.x-e.x, s.y-e.y
    n1 = cp1.x*cp2.y - cp1.y*cp2.x
    n2 = s.x*e.y - s.y*e.x
    n3 = 1.0 / (dcx*dpy - dcy*dpx)
    Point[(n1*dpx - n2*dcx) * n3, (n1*dpy - n2*dcy) * n3]
  end
  
  outputList = subjectPolygon
  cp1 = clipPolygon.last
  for cp2 in clipPolygon
    inputList = outputList
    outputList = []
    s = inputList.last
    for e in inputList
      if inside[e]
        outputList << intersection[] unless inside[s]
        outputList << e
      elsif inside[s]
        outputList << intersection[]
      end
      s = e
    end
    cp1 = cp2
  end
  outputList
end

subjectPolygon = [[50, 150], [200, 50], [350, 150], [350, 300],
                  [250, 300], [200, 250], [150, 350], [100, 250],
                  [100, 200]].collect{|pnt| Point[*pnt]}
 
clipPolygon = [[100, 100], [300, 100], [300, 300], [100, 300]].collect{|pnt| Point[*pnt]}
 
puts sutherland_hodgman(subjectPolygon, clipPolygon)
Output:
(100.0, 116.66666666666667)
(125.00000000000001, 100.0)
(275.0, 100.0)
(300.0, 116.66666666666667)
(300.0, 299.99999999999994)
(250.0, 300.0)
(200, 250)
(175.0, 300.0)
(125.0, 300.0)
(100.0, 250.0)

Rust[edit]

Translation of: Ruby
#[derive(Debug, Clone)]
struct Point {
    x: f64,
    y: f64,
}

#[derive(Debug, Clone)]
struct Polygon(Vec<Point>);

fn is_inside(p: &Point, cp1: &Point, cp2: &Point) -> bool {
    (cp2.x - cp1.x) * (p.y - cp1.y) > (cp2.y - cp1.y) * (p.x - cp1.x)
}

fn compute_intersection(cp1: &Point, cp2: &Point, s: &Point, e: &Point) -> Point {
    let dc = Point {
        x: cp1.x - cp2.x,
        y: cp1.y - cp2.y,
    };
    let dp = Point {
        x: s.x - e.x,
        y: s.y - e.y,
    };
    let n1 = cp1.x * cp2.y - cp1.y * cp2.x;
    let n2 = s.x * e.y - s.y * e.x;
    let n3 = 1.0 / (dc.x * dp.y - dc.y * dp.x);
    Point {
        x: (n1 * dp.x - n2 * dc.x) * n3,
        y: (n1 * dp.y - n2 * dc.y) * n3,
    }
}

fn sutherland_hodgman_clip(subject_polygon: &Polygon, clip_polygon: &Polygon) -> Polygon {
    let mut result_ring = subject_polygon.0.clone();
    let mut cp1 = clip_polygon.0.last().unwrap();
    for cp2 in &clip_polygon.0 {
        let input = result_ring;
        let mut s = input.last().unwrap();
        result_ring = vec![];
        for e in &input {
            if is_inside(e, cp1, cp2) {
                if !is_inside(s, cp1, cp2) {
                    result_ring.push(compute_intersection(cp1, cp2, s, e));
                }
                result_ring.push(e.clone());
            } else if is_inside(s, cp1, cp2) {
                result_ring.push(compute_intersection(cp1, cp2, s, e));
            }
            s = e;
        }
        cp1 = cp2;
    }
    Polygon(result_ring)
}

fn main() {
    let _p = |x: f64, y: f64| Point { x, y };
    let subject_polygon = Polygon(vec![
        _p(50.0, 150.0), _p(200.0, 50.0), _p(350.0, 150.0), _p(350.0, 300.0), _p(250.0, 300.0),
        _p(200.0, 250.0), _p(150.0, 350.0), _p(100.0, 250.0), _p(100.0, 200.0),
    ]);
    let clip_polygon = Polygon(vec![
        _p(100.0, 100.0),_p(300.0, 100.0),_p(300.0, 300.0),_p(100.0, 300.0),
    ]);
    let result = sutherland_hodgman_clip(&subject_polygon, &clip_polygon);
    println!("{:?}", result);
}
Output:
Polygon([
    Point { x: 100, y: 116.66666666666667 }, Point { x: 125.00000000000001, y: 100 }, Point { x: 275, y: 100 },
    Point { x: 300, y: 116.66666666666667 }, Point { x: 300, y: 299.99999999999994 }, Point { x: 250, y: 300 },
    Point { x: 200, y: 250 }, Point { x: 175, y: 300 }, Point { x: 125, y: 300 }, Point { x: 100, y: 250 }])

Scala[edit]

From Java snippet.

import javax.swing.{ JFrame, JPanel }

object SutherlandHodgman extends JFrame with App {
    import java.awt.BorderLayout

    setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE)
    setVisible(true)
    val content = getContentPane()
    content.setLayout(new BorderLayout())
    content.add(SutherlandHodgmanPanel, BorderLayout.CENTER)
    setTitle("SutherlandHodgman")
    pack()
    setLocationRelativeTo(null)
}

object SutherlandHodgmanPanel extends JPanel {
    import java.awt.{ Color, Graphics, Graphics2D }

    setPreferredSize(new java.awt.Dimension(600, 500))

    // subject and clip points are assumed to be valid
    val subject = Seq((50D, 150D), (200D, 50D), (350D, 150D), (350D, 300D), (250D, 300D), (200D, 250D), (150D, 350D), (100D, 250D), (100D, 200D))
    val clipper = Seq((100D, 100D), (300D, 100D), (300D, 300D), (100D, 300D))
    var result = subject

    val len = clipper.size
    for (i <- 0 until len) {
        val len2 = result.size
        val input = result
        result = Seq()

        val A = clipper((i + len - 1) % len)
        val B = clipper(i)

        for (j <- 0 until len2) {
            val P = input((j + len2 - 1) % len2)
            val Q = input(j)

            if (inside(A, B, Q)) {
                if (!inside(A, B, P))
                    result = result :+ intersection(A, B, P, Q)
                result = result :+ Q
            }
            else if (inside(A, B, P))
                result = result :+ intersection(A, B, P, Q)
        }
    }

    override def paintComponent(g: Graphics) {
        import java.awt.RenderingHints._

        super.paintComponent(g)
        val g2 = g.asInstanceOf[Graphics2D]
        g2.translate(80, 60)
        g2.setStroke(new java.awt.BasicStroke(3))
        g2.setRenderingHint(KEY_ANTIALIASING, VALUE_ANTIALIAS_ON)
        g2.draw_polygon(subject, Color.blue)
        g2.draw_polygon(clipper, Color.red)
        g2.draw_polygon(result, Color.green)
    }

    private def inside(a: (Double, Double), b: (Double, Double), c: (Double, Double)) =
        (a._1 - c._1) * (b._2 - c._2) > (a._2 - c._2) * (b._1 - c._1)

    private def intersection(a: (Double, Double), b: (Double, Double), p: (Double, Double), q: (Double, Double)) = {
        val A1 = b._2 - a._2
        val B1 = a._1 - b._1
        val C1 = A1 * a._1 + B1 * a._2
        val A2 = q._2 - p._2
        val B2 = p._1 - q._1
        val C2 = A2 * p._1 + B2 * p._2

        val det = A1 * B2 - A2 * B1
        ((B2 * C1 - B1 * C2) / det, (A1 * C2 - A2 * C1) / det)
    }

    private implicit final class Polygon_drawing(g: Graphics2D) {
        def draw_polygon(points: Seq[(Double, Double)], color: Color) {
            g.setColor(color)
            val len = points.length
            val line = new java.awt.geom.Line2D.Double()
            for (i <- 0 until len) {
                val p1 = points(i)
                val p2 = points((i + 1) % len)
                line.setLine(p1._1, p1._2, p2._1, p2._2)
                g.draw(line)
            }
        }
    }
}

Scheme[edit]

Translation of: ATS
Works with: Guile version at least 2.0
Works with: Gauche Scheme version 0.9.12
Works with: CHICKEN Scheme version 5.3.0
Works with: Gambit Scheme version 4.9.4
;;; Sutherland-Hodgman polygon clipping.

(define (evaluate-line x1 y1 x2 y2 x)
  ;; Given the straight line between (x1,y1) and (x2,y2), evaluate it
  ;; at x.
  (let ((dy (- y2 y1))
        (dx (- x2 x1)))
    (let ((slope (/ dy dx))
          (intercept (/ (- (* dx y1) (* dy x1)) dx)))
      (+ (* slope x) intercept))))

(define (intersection-of-lines x1 y1 x2 y2 x3 y3 x4 y4)
  ;; Given the line between (x1,y1) and (x2,y2), and the line between
  ;; (x3,y3) and (x4,y4), find their intersection.
  (cond ((= x1 x2) (list x1 (evaluate-line x3 y3 x4 y4 x1)))
        ((= x3 x4) (list x3 (evaluate-line x1 y1 x2 y2 x3)))
        (else (let ((denominator (- (* (- x1 x2) (- y3 y4))
                                    (* (- y1 y2) (- x3 x4))))
                    (x1*y2-y1*x2 (- (* x1 y2) (* y1 x2)))
                    (x3*y4-y3*x4 (- (* x3 y4) (* y3 x4))))
                (let ((xnumerator (- (* x1*y2-y1*x2 (- x3 x4))
                                     (* (- x1 x2) x3*y4-y3*x4)))
                      (ynumerator (- (* x1*y2-y1*x2 (- y3 y4))
                                     (* (- y1 y2) x3*y4-y3*x4))))
                  (list (/ xnumerator denominator)
                        (/ ynumerator denominator)))))))

(define (intersection-of-edges e1 e2)
  ;;
  ;; A point is a list of two coordinates, and an edge is a list of
  ;; two points.
  ;;
  ;; I am not using any SRFI-9 records, or the like, that define
  ;; actual new types, although I would do so if writing a more
  ;; serious implementation. Also, I am not using any pattern matcher.
  ;; A pattern matcher would make this code less tedious with
  ;; "cadaddaddr" notations.
  (let ((point1 (car e1))
        (point2 (cadr e1))
        (point3 (car e2))
        (point4 (cadr e2)))
    (let ((x1 (car point1))
          (y1 (cadr point1))
          (x2 (car point2))
          (y2 (cadr point2))
          (x3 (car point3))
          (y3 (cadr point3))
          (x4 (car point4))
          (y4 (cadr point4)))
      (intersection-of-lines x1 y1 x2 y2 x3 y3 x4 y4))))

(define (point-is-left-of-edge? pt edge)
  (let ((x (car pt))
        (y (cadr pt))
        (x1 (caar edge))
        (y1 (cadar edge))
        (x2 (caadr edge))
        (y2 (cadadr edge)))
    ;; Outer product of the vectors (x1,y1)-->(x,y) and
    ;; (x1,y1)-->(x2,y2)
    (negative? (- (* (- x x1) (- y2 y1))
                  (* (- x2 x1) (- y y1))))))

(define (clip-subject-edge subject-edge clip-edge accum)
  (define left-of? point-is-left-of-edge?)
  (define (intersection)
    (intersection-of-edges subject-edge clip-edge))
  (let ((s1 (car subject-edge))
        (s2 (cadr subject-edge)))
    (let ((s2-is-inside? (left-of? s2 clip-edge))
          (s1-is-inside? (left-of? s1 clip-edge)))
      (if s2-is-inside?
          (if s1-is-inside?
              (cons s2 accum)
              (cons s2 (cons (intersection) accum)))
          (if s1-is-inside?
              (cons (intersection) accum)
              accum)))))

(define (for-each-subject-edge i subject-points clip-edge accum)
  (define n (vector-length subject-points))
  (if (= i n)
      (list->vector (reverse accum))
      (let ((s2 (vector-ref subject-points i))
            (s1 (vector-ref subject-points
                            (- (if (zero? i) n i) 1))))
        (let ((accum (clip-subject-edge (list s1 s2)
                                        clip-edge accum)))
          (for-each-subject-edge (+ i 1) subject-points
                                 clip-edge accum)))))

(define (for-each-clip-edge i subject-points clip-points)
  (define n (vector-length clip-points))
  (if (= i n)
      subject-points
      (let ((c2 (vector-ref clip-points i))
            (c1 (vector-ref clip-points (- (if (zero? i) n i) 1))))
        (let ((subject-points
               (for-each-subject-edge 0 subject-points
                                      (list c1 c2) '())))
          (for-each-clip-edge (+ i 1) subject-points clip-points)))))

(define (clip subject-points clip-points)
  (for-each-clip-edge 0 subject-points clip-points))

(define (write-eps subject-points clip-points result-points)

  ;; I use only some of the most basic output procedures. Schemes tend
  ;; to include more advanced means to write output, often resembling
  ;; those of Common Lisp.

  (define (x pt) (exact->inexact (car pt)))
  (define (y pt) (exact->inexact (cadr pt)))

  (define (moveto pt)
    (display (x pt))
    (display " ")
    (display (y pt))
    (display " moveto")
    (newline))

  (define (lineto pt)
    (display (x pt))
    (display " ")
    (display (y pt))
    (display " lineto")
    (newline))

  (define (setrgbcolor rgb)
    (display rgb)
    (display " setrgbcolor")
    (newline))

  (define (simple-word word)
    (lambda ()
      (display word)
      (newline)))

  (define closepath (simple-word "closepath"))
  (define fill (simple-word "fill"))
  (define stroke (simple-word "stroke"))
  (define gsave (simple-word "gsave"))
  (define grestore (simple-word "grestore"))

  (define (showpoly poly line-color fill-color)
    (define n (vector-length poly))
    (moveto (vector-ref poly 0))
    (do ((i 1 (+ i 1)))
        ((= i n))
      (lineto (vector-ref poly i)))
    (closepath)
    (setrgbcolor line-color)
    (gsave)
    (setrgbcolor fill-color)
    (fill)
    (grestore)
    (stroke))

  (define (code s)
    (display s)
    (newline))

  (code "%!PS-Adobe-3.0 EPSF-3.0")
  (code "%%BoundingBox: 40 40 360 360")
  (code "0 setlinewidth")
  (showpoly clip-points ".5 0 0" "1 .7 .7")
  (showpoly subject-points "0 .2 .5" ".4 .7 1")
  (code "2 setlinewidth")
  (code "[10 8] 0 setdash")
  (showpoly result-points ".5 0 .5" ".7 .3 .8")
  (code "%%EOF"))

(define (write-eps-to-file outfile subject-points clip-points
                           result-points)
  (with-output-to-file outfile
    (lambda ()
      (write-eps subject-points clip-points result-points))))

(define subject-points
  #((50 150)
    (200 50)
    (350 150)
    (350 300)
    (250 300)
    (200 250)
    (150 350)
    (100 250)
    (100 200)))

(define clip-points
  #((100 100)
    (300 100)
    (300 300)
    (100 300)))

(define result-points (clip subject-points clip-points))

(display result-points)
(newline)
(write-eps-to-file "sutherland-hodgman.eps"
                   subject-points clip-points result-points)
(display "Wrote sutherland-hodgman.eps")
(newline)
Output:
#((100 350/3) (125 100) (275 100) (300 350/3) (300 300) (250 300) (200 250) (175 300) (125 300) (100 250))
Wrote sutherland-hodgman.eps

The polygons of the task

Sidef[edit]

Translation of: Ruby
class Point(x, y) {
    method to_s {
        "(#{'%.2f' % x}, #{'%.2f' % y})"
    }
}

func sutherland_hodgman(subjectPolygon, clipPolygon) {
  var inside = { |cp1, cp2, p|
    ((cp2.x-cp1.x)*(p.y-cp1.y)) > ((cp2.y-cp1.y)*(p.x-cp1.x))
  }

  var intersection = { |cp1, cp2, s, e|
    var (dcx, dcy) = (cp1.x-cp2.x, cp1.y-cp2.y)
    var (dpx, dpy) = (s.x-e.x, s.y-e.y)
    var n1 = (cp1.x*cp2.y - cp1.y*cp2.x)
    var n2 = (s.x*e.y - s.y*e.x)
    var n3 = (1 / (dcx*dpy - dcy*dpx))
    Point((n1*dpx - n2*dcx) * n3, (n1*dpy - n2*dcy) * n3)
  }

  var outputList = subjectPolygon
  var cp1 = clipPolygon.last
  for cp2 in clipPolygon {
    var inputList = outputList
    outputList = []
    var s = inputList.last
    for e in inputList {
      if (inside(cp1, cp2, e)) {
        outputList << intersection(cp1, cp2, s, e) if !inside(cp1, cp2, s)
        outputList << e
      }
      elsif(inside(cp1, cp2, s)) {
        outputList << intersection(cp1, cp2, s, e)
      }
      s = e
    }
    cp1 = cp2
  }
  outputList
}

var subjectPolygon = [
    [50,  150], [200,  50], [350, 150], [350, 300],
    [250, 300], [200, 250], [150, 350], [100, 250],
    [100, 200]
].map{|pnt| Point(pnt...) }

var clipPolygon = [
    [100, 100], [300, 100],
    [300, 300], [100, 300]
].map{|pnt| Point(pnt...) }

sutherland_hodgman(subjectPolygon, clipPolygon).each { .say }
Output:
(100.00, 116.67)
(125.00, 100.00)
(275.00, 100.00)
(300.00, 116.67)
(300.00, 300.00)
(250.00, 300.00)
(200.00, 250.00)
(175.00, 300.00)
(125.00, 300.00)
(100.00, 250.00)

Standard ML[edit]

Translation of: ATS
Works with: MLton version 20210117
(* Sutherland-Hodgman polygon clipping. *)

fun evaluate_line (x1 : real, y1 : real,
                   x2 : real, y2 : real,
                   x : real) =
    let
      val dy = y2 - y1
      and dx = x2 - x1
      val slope = dy / dx
      and intercept = ((dx * y1) - (dy * x1)) / dx
    in
      (slope * x) + intercept
    end

fun intersection_of_lines (x1 : real, y1 : real,
                           x2 : real, y2 : real,
                           x3 : real, y3 : real,
                           x4 : real, y4 : real) =
    if Real.== (x1, x2) then
      (x1, evaluate_line (x3, y3, x4, y4, x1))
    else if Real.== (x3, x4) then
      (x3, evaluate_line (x1, y1, x2, y2, x3))
    else
      let
        val denominator =
            ((x1 - x2) * (y3 - y4)) - ((y1 - y2) * (x3 - x4))
        and x1y2_y1x2 = (x1 * y2) - (y1 * x2)
        and x3y4_y3x4 = (x3 * y4) - (y3 * x4)

        val xnumerator =
            (x1y2_y1x2 * (x3 - x4)) - ((x1 - x2) * x3y4_y3x4)
        and ynumerator =
            (x1y2_y1x2 * (y3 - y4)) - ((y1 - y2) * x3y4_y3x4)
      in
        (xnumerator / denominator,
         ynumerator / denominator)
      end

fun intersection_of_edges (((x1, y1), (x2, y2)),
                           ((x3, y3), (x4, y4))) =
    intersection_of_lines (x1, y1, x2, y2, x3, y3, x4, y4)

fun point_is_left_of_edge ((x, y), ((x1, y1), (x2, y2))) =
    (* Outer product of the vectors (x1,y1)-->(x,y) and
       (x1,y1)-->(x2,y2). *)
    ((x - x1) * (y2 - y1)) - ((x2 - x1) * (y - y1)) < 0.0

fun clip_subject_edge (subject_edge, clip_edge, accum) =
    let
      fun intersection () =
          intersection_of_edges (subject_edge, clip_edge)

      val (s1, s2) = subject_edge
      val s2_is_inside = point_is_left_of_edge (s2, clip_edge)
      and s1_is_inside = point_is_left_of_edge (s1, clip_edge)
    in
      case (s2_is_inside, s1_is_inside) of
          (true, true) => s2 :: accum
        | (true, false) => s2 :: intersection () :: accum
        | (false, true) => intersection () :: accum
        | (false, false) => accum
    end

fun for_each_subject_edge (i, subject_points, clip_edge, accum) =
    let
      val n = Array.length subject_points
    in
      if i = n then
        Array.fromList (rev accum)
      else
        let
          val s2 = Array.sub (subject_points, i)
          and s1 = (if i = 0 then
                      Array.sub (subject_points, n - 1)
                    else
                      Array.sub (subject_points, i - 1))
          val accum = clip_subject_edge ((s1, s2), clip_edge, accum)
        in
          for_each_subject_edge (i + 1, subject_points, clip_edge,
                                 accum)
        end
    end

fun for_each_clip_edge (i, subject_points, clip_points) =
    let
      val n = Array.length clip_points
    in
      if i = n then
        subject_points
      else
        let
          val c2 = Array.sub (clip_points, i)
          and c1 = (if i = 0 then
                      Array.sub (clip_points, n - 1)
                    else
                      Array.sub (clip_points, i - 1))
          val subject_points =
              for_each_subject_edge (0, subject_points, (c1, c2), [])
        in
          for_each_clip_edge (i + 1, subject_points, clip_points)
        end
    end

fun clip (subject_points, clip_points) =
    for_each_clip_edge (0, subject_points, clip_points)

fun write_eps (outf, subject_points, clip_points, result_points) =
    (* The EPS code that will be generated is based on that which is
       generated by the C implementation of this task. *)
    let
      fun moveto (x, y) =
          (TextIO.output (outf, Real.toString x);
           TextIO.output (outf, " ");
           TextIO.output (outf, Real.toString y);
           TextIO.output (outf, " moveto\n"))
      fun lineto (x, y) =
          (TextIO.output (outf, Real.toString x);
           TextIO.output (outf, " ");
           TextIO.output (outf, Real.toString y);
           TextIO.output (outf, " lineto\n"))
      fun setrgbcolor rgb =
          (TextIO.output (outf, rgb);
           TextIO.output (outf, " setrgbcolor\n"))
      fun closepath () = TextIO.output (outf, "closepath\n")
      fun fill () = TextIO.output (outf, "fill\n")
      fun stroke () = TextIO.output (outf, "stroke\n")
      fun gsave () = TextIO.output (outf, "gsave\n")
      fun grestore () = TextIO.output (outf, "grestore\n")
      fun showpoly (poly, line_color, fill_color) =
          let
            val n = Array.length poly
          in
            moveto (Array.sub (poly, 0));
            Array.app lineto poly;
            closepath ();
            setrgbcolor line_color;
            gsave ();
            setrgbcolor fill_color;
            fill ();
            grestore ();
            stroke ()
          end
    in
      TextIO.output (outf, "%!PS-Adobe-3.0 EPSF-3.0\n");
      TextIO.output (outf, "%%BoundingBox: 40 40 360 360\n");
      TextIO.output (outf, "0 setlinewidth\n");
      showpoly (clip_points, ".5 0 0", "1 .7 .7");
      showpoly (subject_points, "0 .2 .5", ".4 .7 1");
      TextIO.output (outf, "2 setlinewidth\n");
      TextIO.output (outf, "[10 8] 0 setdash\n");
      showpoly (result_points, ".5 0 .5", ".7 .3 .8");
      TextIO.output (outf, "%%EOF\n")
    end

fun write_eps_to_file (outfile, subject_points, clip_points,
                       result_points) =
    let
      val outf = TextIO.openOut outfile
    in
      write_eps (outf, subject_points, clip_points, result_points);
      TextIO.closeOut outf
    end

val subject_points =
    Array.fromList
      [(50.0, 150.0),
       (200.0, 50.0),
       (350.0, 150.0),
       (350.0, 300.0),
       (250.0, 300.0),
       (200.0, 250.0),
       (150.0, 350.0),
       (100.0, 250.0),
       (100.0, 200.0)]

val clip_points =
    Array.fromList
      [(100.0, 100.0),
       (300.0, 100.0),
       (300.0, 300.0),
       (100.0, 300.0)]

val result_points = clip (subject_points, clip_points)

fun print_point (x, y) =
    (TextIO.print " (";
     TextIO.print (Real.toString x);
     TextIO.print " ";
     TextIO.print (Real.toString y);
     TextIO.print ")")
  ;

Array.app print_point result_points;
TextIO.print "\n";
write_eps_to_file ("sutherland-hodgman.eps",
                   subject_points, clip_points, result_points);
TextIO.print "Wrote sutherland-hodgman.eps\n";

(*
local variables:
mode: SML
sml-indent-level: 2
end:
*)
Output:
$ mlton sutherland-hodgman.sml && ./sutherland-hodgman
 (100 116.666666667) (125 100) (275 100) (300 116.666666667) (300 300) (250 300) (200 250) (175 300) (125 300) (100 250)
Wrote sutherland-hodgman.eps

The polygons of the task.

Swift[edit]

Translation of: Rust
struct Point {
  var x: Double
  var y: Double
}

struct Polygon {
  var points: [Point]

  init(points: [Point]) {
    self.points = points
  }

  init(points: [(Double, Double)]) {
    self.init(points: points.map({ Point(x: $0.0, y: $0.1) }))
  }
}

func isInside(_ p1: Point, _ p2: Point, _ p3: Point) -> Bool {
  (p3.x - p2.x) * (p1.y - p2.y) > (p3.y - p2.y) * (p1.x - p2.x)
}

func computeIntersection(_ p1: Point, _ p2: Point, _ s: Point, _ e: Point) -> Point {
  let dc = Point(x: p1.x - p2.x, y: p1.y - p2.y)
  let dp = Point(x: s.x - e.x, y: s.y - e.y)
  let n1 = p1.x * p2.y - p1.y * p2.x
  let n2 = s.x * e.y - s.y * e.x
  let n3 = 1.0 / (dc.x * dp.y - dc.y * dp.x)

  return Point(x: (n1 * dp.x - n2 * dc.x) * n3, y: (n1 * dp.y - n2 * dc.y) * n3)
}

func sutherlandHodgmanClip(subjPoly: Polygon, clipPoly: Polygon) -> Polygon {
  var ring = subjPoly.points
  var p1 = clipPoly.points.last!

  for p2 in clipPoly.points {
    let input = ring
    var s = input.last!

    ring = []

    for e in input {
      if isInside(e, p1, p2) {
        if !isInside(s, p1, p2) {
          ring.append(computeIntersection(p1, p2, s, e))
        }

        ring.append(e)
      } else if isInside(s, p1, p2) {
        ring.append(computeIntersection(p1, p2, s, e))
      }

      s = e
    }

    p1 = p2
  }

  return Polygon(points: ring)
}

let subj = Polygon(points: [
  (50.0, 150.0),
  (200.0, 50.0),
  (350.0, 150.0),
  (350.0, 300.0),
  (250.0, 300.0),
  (200.0, 250.0),
  (150.0, 350.0),
  (100.0, 250.0),
  (100.0, 200.0)
])

let clip = Polygon(points: [
  (100.0, 100.0),
  (300.0, 100.0),
  (300.0, 300.0),
  (100.0, 300.0)
])

print(sutherlandHodgmanClip(subjPoly: subj, clipPoly: clip))
Output:
Polygon(points: [Point(x: 100.0, y: 116.66666666666667), Point(x: 125.00000000000001, y: 100.0), Point(x: 275.0, y: 100.0), Point(x: 300.0, y: 116.66666666666667), Point(x: 300.0, y: 299.99999999999994), Point(x: 250.0, y: 300.0), Point(x: 200.0, y: 250.0), Point(x: 175.0, y: 300.0), Point(x: 125.0, y: 300.0), Point(x: 100.0, y: 250.0)])

Tcl[edit]

# Find intersection of an arbitrary polygon with a convex one.
package require Tcl 8.6

#	Does the path (x0,y0)->(x1,y1)->(x2,y2) turn clockwise
#	or counterclockwise?
proc cw {x0 y0 x1 y1 x2 y2} {
    set dx1 [expr {$x1 - $x0}]; set dy1 [expr {$y1 - $y0}]
    set dx2 [expr {$x2 - $x0}]; set dy2 [expr {$y2 - $y0}]
    # (0,0,$dx1*$dy2 - $dx2*$dy1) is the crossproduct of
    # ($x1-$x0,$y1-$y0,0) and ($x2-$x0,$y2-$y0,0). 
    # Its z-component is positive if the turn
    # is clockwise, negative if the turn is counterclockwise.
    set pr1 [expr {$dx1 * $dy2}]
    set pr2 [expr {$dx2 * $dy1}]
    if {$pr1 > $pr2} {
	# Clockwise
	return 1
    } elseif {$pr1 < $pr2} {
	# Counter-clockwise
	return -1
    } elseif {$dx1*$dx2 < 0 || $dy1*$dy2 < 0} {
	# point 0 is the middle point
	return 0
    } elseif {($dx1*$dx1 + $dy1*$dy1) < ($dx2*$dx2 + $dy2+$dy2)} {
	# point 1 is the middle point
	return 0
    } else {
	# point 2 lies on the segment joining points 0 and 1
	return 1
    }
}

#	Calculate the point of intersection of two lines
#	containing the line segments (x1,y1)-(x2,y2) and (x3,y3)-(x4,y4)
proc intersect {x1 y1 x2 y2 x3 y3 x4 y4} {
    set d [expr {($y4 - $y3) * ($x2 - $x1) - ($x4 - $x3) * ($y2 - $y1)}]
    set na [expr {($x4 - $x3) * ($y1 - $y3) - ($y4 - $y3) * ($x1 - $x3)}]
    if {$d == 0} {
	return {}
    }
    set r [list \
	    [expr {$x1 + $na * ($x2 - $x1) / $d}] \
	    [expr {$y1 + $na * ($y2 - $y1) / $d}]]
    return $r
}

#	Coroutine that yields the elements of a list in pairs
proc pairs {list} {
    yield [info coroutine]
    foreach {x y} $list {
	yield [list $x $y]
    }
    return {}
}

#	Coroutine to clip one segment of a polygon against a line.
proc clipsegment {inside0 cx0 cy0 cx1 cy1 sx0 sy0 sx1 sy1} {
    set inside1 [expr {[cw $cx0 $cy0 $cx1 $cy1 $sx1 $sy1] > 0}]
    if {$inside1} {
	if {!$inside0} {
	    set int [intersect $cx0 $cy0 $cx1 $cy1 \
		    $sx0 $sy0 $sx1 $sy1]
	    if {[llength $int] >= 0} {
		yield $int
	    }
	}
	yield [list $sx1 $sy1]
    } else {
	if {$inside0} {
	    set int [intersect $cx0 $cy0 $cx1 $cy1 \
		    $sx0 $sy0 $sx1 $sy1]
	    if {[llength $int] >= 0} {
		yield $int
	    }
	}
    }
    return $inside1
}

#	Coroutine to perform one step of Sutherland-Hodgman polygon clipping
proc clipstep {source cx0 cy0 cx1 cy1} {
    yield [info coroutine]
    set pt0 [{*}$source]
    if {[llength $pt0] == 0} {
	return
    }
    lassign $pt0 sx0 sy0
    set inside0 [expr {[cw $cx0 $cy0 $cx1 $cy1 $sx0 $sy0] > 0}]
    set finished 0
    while {!$finished} {
	set thispt [{*}$source]
	if {[llength $thispt] == 0} {
	    set thispt $pt0
	    set finished 1
	}
	lassign $thispt sx1 sy1
	set inside0 [clipsegment $inside0 \
		$cx0 $cy0 $cx1 $cy1 $sx0 $sy0 $sx1 $sy1]
	set sx0 $sx1
	set sy0 $sy1
    }
    return {}
}

#	Perform Sutherland-Hodgman polygon clipping
proc clippoly {cpoly spoly} {
    variable clipindx
    set source [coroutine clipper[incr clipindx] pairs $spoly]
    set cx0 [lindex $cpoly end-1]
    set cy0 [lindex $cpoly end]
    foreach {cx1 cy1} $cpoly {
	set source [coroutine clipper[incr clipindx] \
		clipstep $source $cx0 $cy0 $cx1 $cy1]
	set cx0 $cx1; set cy0 $cy1
    }
    set result {}
    while {[llength [set pt [{*}$source]]] > 0} {
	lappend result {*}$pt
    }
    return $result
}

The specifics of the task:

Library: Tk
package require Tk

grid [canvas .c -width 400 -height 400 -background \#ffffff]
proc demonstrate {cpoly spoly} {
    set rpoly [clippoly $cpoly $spoly]
    puts $rpoly
    .c create polygon $cpoly -outline \#ff9999 -fill {} -width 5 
    .c create polygon $spoly -outline \#9999ff -fill {} -width 3
    .c create polygon $rpoly -fill \#99ff99 -outline black -width 1
}

demonstrate {100 100 300 100 300 300 100 300} \
    {50 150 200 50 350 150 350 300 250 300 200 250 150 350 100 250 100 200}
Output:
300 116 300 300 250 300 200 250 175 300 125 300 100 250 100 200 100 200 100 116 124 100 275 100

Sutherland-Hodgman.gif

TypeScript[edit]

interface XYCoords  {
    x : number;
    y : number;
}

const inside = ( cp1 : XYCoords, cp2 : XYCoords, p : XYCoords) : boolean => {
    return (cp2.x-cp1.x)*(p.y-cp1.y) > (cp2.y-cp1.y)*(p.x-cp1.x);
};

const intersection = ( cp1 : XYCoords ,cp2 : XYCoords ,s : XYCoords, e : XYCoords ) : XYCoords => {
    const dc = {
	x:  cp1.x - cp2.x,
	y : cp1.y - cp2.y
	 },
    dp = { x: s.x - e.x,
	   y : s.y - e.y
	 },
    n1 = cp1.x * cp2.y - cp1.y * cp2.x,
    n2 = s.x * e.y - s.y * e.x, 
    n3 = 1.0 / (dc.x * dp.y - dc.y * dp.x);
    return { x : (n1*dp.x - n2*dc.x) * n3,
	     y : (n1*dp.y - n2*dc.y) * n3
	   };
};

export const sutherland_hodgman = ( subjectPolygon : Array<XYCoords>,
				   clipPolygon : Array<XYCoords> ) : Array<XYCoords> => {
    
    let cp1 : XYCoords = clipPolygon[clipPolygon.length-1];
    let cp2 : XYCoords;
    let s : XYCoords;
    let e : XYCoords;
    
    let outputList : Array<XYCoords> = subjectPolygon;
    
    for( var j in clipPolygon ) {
        cp2 = clipPolygon[j];
        var inputList = outputList;
        outputList = [];
        s = inputList[inputList.length - 1]; // last on the input list
        for (var i in inputList) {
	    e = inputList[i];
	    if (inside(cp1,cp2,e)) {
                if (!inside(cp1,cp2,s)) {
		    outputList.push(intersection(cp1,cp2,s,e));
                }
                outputList.push(e);
	    }
	    else if (inside(cp1,cp2,s)) {
                outputList.push(intersection(cp1,cp2,s,e));
	    }
	    s = e;
        }
        cp1 = cp2;
    }
    return outputList
}


Wren[edit]

Library: DOME
Library: Wren-polygon
import "graphics" for Canvas, Color
import "dome" for Window
import "./polygon" for Polygon

class SutherlandHodgman {
    construct new(width, height, subject, clipper) {
        Window.title = "Sutherland-Hodgman"
        Window.resize(width, height)
        Canvas.resize(width, height)
        _subject = subject
        _result = subject.toList
        _clipper = clipper
    }

    init() {
        clipPolygon()
        System.print("Clipped polygon points:")
        for (p in _result) {
            p[1] = (1000*p[1]).round/1000
            System.print(p)
        }
        // display all 3 polygons
        Polygon.quick(_subject).drawfill(Color.blue)
        Polygon.quick(_clipper).drawfill(Color.red)
        Polygon.quick(_result ).drawfill(Color.green)
    }

    clipPolygon() {
        var len = _clipper.count
        for (i in 0...len) {
            var len2 = _result.count
            var input = _result
            _result = []
            var a = _clipper[(i + len - 1) % len]
            var b = _clipper[i]

            for (j in 0...len2) {
                var p = input[(j + len2 - 1) % len2]
                var q = input[j]

                if (isInside(a, b, q)) {
                    if (!isInside(a, b, p)) _result.add(intersection(a, b, p, q))
                    _result.add(q)
                } else if (isInside(a, b, p)) _result.add(intersection(a, b, p, q))
            }
        }
    }

    isInside(a, b, c) {  (a[0] - c[0]) * (b[1] - c[1]) > (a[1] - c[1]) * (b[0] - c[0]) }

    intersection(a, b, p, q) {
        var a1 = b[1] - a[1]
        var b1 = a[0] - b[0]
        var c1 = a1 * a[0] + b1 * a[1]

        var a2 = q[1] - p[1]
        var b2 = p[0] - q[0]
        var c2 = a2 * p[0] + b2 * p[1]

        var d = a1 * b2 - a2 * b1
        var x = (b2 * c1 - b1 * c2) / d
        var y = (a1 * c2 - a2 * c1) / d

        return [x, y]
    }
    update() {}

    draw(alpha) {}
}

var subject = [
    [ 50, 150], [200,  50], [350, 150], 
    [350, 300], [250, 300], [200, 250],
    [150, 350], [100, 250], [100, 200]
]
var clipper = [
    [100, 100], [300, 100],
    [300, 300], [100, 300]
]
var Game = SutherlandHodgman.new(600, 500, subject, clipper)
Output:
Clipped polygon points:
[100, 116.667]
[125, 100]
[275, 100]
[300, 116.667]
[300, 300]
[250, 300]
[200, 250]
[175, 300]
[125, 300]
[100, 250]

Yabasic[edit]

Translation of: BBC BASIC
open window 400, 400
backcolor 0,0,0
clear window

DPOL = 8
DREC = 3
CX = 1 : CY = 2

dim poligono(DPOL, 2)
dim rectang(DREC, 2)
dim clipped(DPOL + DREC, 2)

for n = 0 to DPOL : read poligono(n, CX), poligono(n, CY) : next n
DATA 50,150, 200,50, 350,150, 350,300, 250,300, 200,250, 150,350, 100,250, 100,200
for n = 0 to DREC : read rectang(n, CX), rectang(n, CY) : next n
DATA 100,100, 300,100, 300,300, 100,300


color 255,0,0
dibuja(poligono(), DPOL)
color 0,0,255
dibuja(rectang(), DREC)

nvert = FNsutherland_hodgman(poligono(), rectang(), clipped(), DPOL + DREC)
color 250,250,0
dibuja(clipped(), nvert - 1)


sub dibuja(figura(), i)
	local n
	
	print
	new curve
	for n = 0 to i
		line to figura(n, CX), figura(n, CY)
		print figura(n, CX), ", ", figura(n, CY)
	next n
	close curve
end sub


sub FNsutherland_hodgman(subj(), clip(), out(), n)
	local i, j, o, tclip, p1(2), p2(2), s(2), e(2), p(2), inp(n, 2)
	
	FOR o = 0 TO arraysize(subj(), 1) : out(o, CX) = subj(o, CX) : out(o, CY) = subj(o, CY) : NEXT o
	
	tclip = arraysize(clip(),1)
	p1(CX) = clip(tclip, CX) : p1(CY) = clip(tclip, CY)
	
	FOR i = 0 TO tclip
	    p2(CX) = clip(i, CX) : p2(CY) = clip(i, CY)
	    FOR n = 0 TO o - 1 : inp(n, CX) = out(n, CX) : inp(n, CY) = out(n, CY) : NEXT n : o = 0
	  	IF n >= 2 THEN
	            s(CX) = inp(n - 1, CX) : s(CY) = inp(n - 1, CY)

	    	    FOR j = 0 TO n - 1
	      		e(CX) = inp(j, CX) : e(CY) = inp(j, CY)
	      		IF FNside(e(), p1(), p2()) THEN
	        		IF NOT FNside(s(), p1(), p2()) THEN
	          			PROCintersection(p1(), p2(), s(), e(), p())
	          			out(o, CX) = round(p(CX)) : out(o, CY) = round(p(CY))
	          			o = o + 1
	        		ENDIF
	        		out(o, CX) = round(e(CX)) : out(o, CY) = round(e(CY))
	        		o = o + 1
	      		ELSE
	        		IF FNside(s(), p1(), p2()) THEN
	          			PROCintersection(p1(), p2(), s(), e(), p())
	          			out(o, CX) = round(p(CX)) : out(o, CY) = round(p(CY))
	          			o = o + 1
	        		ENDIF
	      		ENDIF
	      		s(CX) = e(CX) : s(CY) = e(CY)
	    	    NEXT j
	  	ENDIF
	  	p1(CX) = p2(CX) : p1(CY) = p2(CY)
	NEXT i
	return o
end sub


sub FNside(p(), p1(), p2())
	return  (p2(CX) - p1(CX)) * (p(CY) - p1(CY)) > (p2(CY) - p1(CY)) * (p(CX) - p1(CX))
end sub


sub PROCintersection(p1(), p2(), p3(), p4(), p())
	LOCAL a(2), b(2), k, l, m
	
	a(CX) = p1(CX) - p2(CX) : a(CY) = p1(CY) - p2(CY)
	b(CX) = p3(CX) - p4(CX) : b(CY) = p3(CY) - p4(CY)
	k = p1(CX) * p2(CY) - p1(CY) * p2(CX)
	l = p3(CX) * p4(CY) - p3(CY) * p4(CX)
	m = 1 / (a(CX) * b(CY) - a(CY) * b(CX))
	p(CX) =  m * (k * b(CX) - l * a(CX))
	p(CY) =  m * (k * b(CY) - l * a(CY))
	
end sub


sub round(n)
	return int(n + .5)
end sub

zkl[edit]

Translation of: C
Translation of: Wikipedia

Uses the PPM class from http://rosettacode.org/wiki/Bitmap/Bresenham%27s_line_algorithm#zkl

class P{	// point
   fcn init(_x,_y){ var [const] x=_x.toFloat(), y=_y.toFloat() }
   fcn __opSub(p) { self(x - p.x, y - p.y) }
   fcn cross(p)   { x*p.y - y*p.x          }
   fcn toString   { "(%7.2f,%7.2f)".fmt(x,y) }
   var [const,proxy] ps=fcn{ T(x.toInt(),y.toInt()) };    // property
}
fcn shClipping(clip,polygon){
   inputList,outputList,clipEdge:=List(), polygon.copy(), List(Void,clip[-1]);
   foreach p in (clip){
      clipEdge.del(0).append(p);
      inputList.clear().extend(outputList);
      outputList.clear();
      S:=inputList[-1];
      foreach E in (inputList){
         if(leftOf(clipEdge,E)){
	    if(not leftOf(clipEdge,S))
	       outputList.append(intersection(S,E,clipEdge));
	    outputList.append(E);
	 }
	 else if(leftOf(clipEdge,S))
	         outputList.append(intersection(S,E,clipEdge));
	 S=E;
      }
   }
   outputList
}
fcn leftOf(line,p){ //-->True (p is left of line), direction of line matters
   p1,p2:=line;		// line is (p1,p2)
   (p2-p1).cross(p-p2)>0;
}
fcn intersection(p1,p2, line){	//-->Point of intersection or False
   p3,p4:=line;
   dx,dy,d:=p2-p1, p3-p4, p1-p3;
   // x0 + a dx = y0 + b dy ->
   // x0 X dx = y0 X dx + b dy X dx ->
   // b = (x0 - y0) X dx / (dy X dx)
   dyx:=dy.cross(dx);
   if(not dyx) return(False);  // parallel lines, could just throw on next line
   dyx=d.cross(dx)/dyx;
   P(p3.x + dyx*dy.x, p3.y + dyx*dy.y);
}
fcn drawPolygon(ppm,listOfPoints,rgb){
   foreach n in (listOfPoints.len()-1){
      ppm.line(listOfPoints[n].ps.xplode(),listOfPoints[n+1].ps.xplode(),rgb);
   }
   ppm.line(listOfPoints[0].ps.xplode(),listOfPoints[-1].ps.xplode(),rgb);
}
ppm:=PPM(400,400);
clip:=T( P(100,100), P(300,100), P(300,300), P(100,300) );
polygon:=T( P( 50,150),P(200, 50),P(350,150),
	    P(350,300),P(250,300),P(200,250),
	    P(150,350),P(100,250),P(100,200) );
drawPolygon(ppm,polygon,0x0000ff);	// blue: polygon
ppm.flood(200,200,0x000030);
drawPolygon(ppm,clip,0xff0000);		// red:  clip region

clipped:=shClipping(clip,polygon);
drawPolygon(ppm,clipped,0x00ff00);	// green: clipped polygon
ppm.flood(200,200,0x003000);		// which is the clipped region anyway
clipped.apply('wrap(p){ ppm.cross(p.ps.xplode(),0x00ff00) }); // mark vertices

ppm.writeJPGFile("sutherland_hodgman.zkl.jpg");

println("Clipped polygon has ",clipped.len()," points:");
clipped.pump(Console.println);
Output:

Until local image uploading is re-enabled, see this image.

Clipped polygon has 10 points:
( 100.00, 116.67)
( 125.00, 100.00)
( 275.00, 100.00)
( 300.00, 116.67)
( 300.00, 300.00)
( 250.00, 300.00)
( 200.00, 250.00)
( 175.00, 300.00)
( 125.00, 300.00)
( 100.00, 250.00)