Sutherland-Hodgman polygon clipping: Difference between revisions

From Rosetta Code
← Older editNewer edit →
Content deleted Content added
→‎{{header|MATLAB}}: Added PureBasic
Line 359: Line 359:
</pre>
</pre>
[[File:Sutherland-Hodgman_haskell.png]]
[[File:Sutherland-Hodgman_haskell.png]]

=={{header|J}}==
'''Solution:'''
<lang j> 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
)</lang>

'''Example use:'''
<lang j> 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</lang>


=={{header|MATLAB}}==
=={{header|MATLAB}}==
Line 456: Line 499:
>> patch(clippedSubject(:,1),clippedSubject(:,2),0);</lang>
>> patch(clippedSubject(:,1),clippedSubject(:,2),0);</lang>
[[File:Sutherland-Hodgman_MATLAB.png]]
[[File:Sutherland-Hodgman_MATLAB.png]]

=={{header|PureBasic}}==
=={{header|PureBasic}}==
{{trans|Go}}
{{trans|Go}}

Revision as of 23:51, 22 June 2011

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.

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

Ada

<lang Ada>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;</lang>

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

D

No extra credit.

Translation of: Go

<lang d>import std.stdio, std.typecons;

alias Tuple!(double, "x", double, "y") Point;

Point[] clip(const Point[] pts, const Point[] clips) { 

   bool isInside(Point p, Point cp1, Point cp2) {
       return (cp2.x - cp1.x) * (p.y - cp1.y) >
              (cp2.y - cp1.y) * (p.x - cp1.x);
   }

   Point intersection(Point cp1, Point cp2, Point s, Point e) {
       auto dc = Point(cp1.x - cp2.x, cp1.y - cp2.y);
       auto dp = Point(s.x - e.x, s.y - e.y);
       auto n1 = cp1.x * cp2.y - cp1.y * cp2.x;
       auto n2 = s.x * e.y - s.y * e.x;
       auto n3 = 1 / (dc.x * dp.y - dc.y * dp.x);
       return Point((n1 * dp.x - n2 * dc.x) * n3,
                    (n1 * dp.y - n2 * dc.y) * n3);
   }

   auto result = pts.dup; // save input
   Point cp1 = clips[$ - 1]; // not const
   foreach (cp2; clips) {
       auto inputList = result;
       result = [];
       auto s = inputList[$ - 1];
       foreach (e; inputList) {
           if (isInside(e, cp1, cp2)) {
               if (!isInside(s, cp1, cp2))
                   result ~= intersection(cp1, cp2, s, e);
               result ~= e;
           } else if (isInside(s, cp1, cp2)) {
               result ~= intersection(cp1, cp2, s, e);
           }
           s = e;
       }
       cp1 = cp2;
   }
   return result;

}


void main() { 

   alias Point P;
   auto subjectPolygon = [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)];

   auto clipPolygon = [P(100, 100), P(300, 100),
                       P(300, 300), P(100, 300)];

   foreach (p; clip(subjectPolygon, clipPolygon))
       writeln(p);

}</lang> Output: 

Tuple!(double,"x",double,"y")(100, 116.667)
Tuple!(double,"x",double,"y")(125, 100)
Tuple!(double,"x",double,"y")(275, 100)
Tuple!(double,"x",double,"y")(300, 116.667)
Tuple!(double,"x",double,"y")(300, 300)
Tuple!(double,"x",double,"y")(250, 300)
Tuple!(double,"x",double,"y")(200, 250)
Tuple!(double,"x",double,"y")(175, 300)
Tuple!(double,"x",double,"y")(125, 300)
Tuple!(double,"x",double,"y")(100, 250)

Go

No extra credit today. <lang go>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)

}</lang> Output: 

[{100 116.66667} {125 100} {275 100} {300 116.66667} {300 300} {250 300} {200 250} {175 300} {125 300} {100 250}]

Haskell

<lang haskell> 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

</lang>

Print the resulting list of points and display the polygons in a window.
<lang haskell> 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

</lang> 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)]

J

Solution: <lang j> 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 )</lang>

Example use: <lang j> 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</lang>

MATLAB

<lang MATLAB>%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 = [];
       previousVertex = inputList(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</lang>

Output: <lang MATLAB>>> 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)],'blue') >> hold on >> plot([clipPolygon(:,1);clipPolygon(1,1)],[clipPolygon(:,2);clipPolygon(1,2)],'r') >> patch(clippedSubject(:,1),clippedSubject(:,2),0);</lang>

PureBasic

Translation of: Go

<lang PureBasic>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</lang> Sample 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)

Tcl

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

  1. Does the path (x0,y0)->(x1,y1)->(x2,y2) turn clockwise
  2. 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 

   }

}

  1. Calculate the point of intersection of two lines
  2. 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

}

  1. 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 {}

}

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

}

  1. 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 {}

}

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

}</lang>

The specifics of the task:

Library: Tk

<lang tcl>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}</lang>

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

Retrieved from "https://rosettacode.org/wiki/Sutherland-Hodgman_polygon_clipping?oldid=143748"