Algoritmos geométricos
Algoritmo de Sutherland-Hodgman
Sutherland-Hodgman
Empezando por el conjunto inicial de vertices del polígono, primero recorta el polígono contra una frontera para producir una nueva secuencia de vertices, con esta nueva secuencia se recorta contra otra frontera y así sucesivamente con las restantes.
Los polígonos cóncavos se pueden desplegar con líneas ajenas cuando el polígono recortado debe tener dos o más secciones separadas. Lo cual requiere medidas adicionales en estos casos como por ejemplo dividir el polígono cóncavo en varios convexos y procesarlos por separado.
Recorte de polígonos: algoritmo de Sutherland-Hodgman
- Sirve para cualquier ventana poligonal convexa, no sólo ventanas rectangulares
- Fácilmente generalizable a 3D
- Cada borde de la ventana (convexa) define dos semiplanos (semiespacios en 3D):
uno interior (in) y otro exterior (out) - Input: lista de puntos (vértices consecutivos del polígono original)
- Output: lista de puntos (vértices del polígono recortado)
- Algoritmo:
- Repetir para todos los bordes de la ventana
- Repetir para todos los lados (parejas de puntos,
, de la lista) del polígono:
(I: intersección de lado del polígono con el borde de la ventana)- Si S in y P in , introducir P en la lista
- Si S in y P out, introducir I en la lista
- Si S out y P out, no introducir nada en la lista
- Si S out y P in , introducir I y P en la lista
- Cerrar el polígono (repetir el primer punto al final de la lista) si es necesario
- Repetir para todos los lados (parejas de puntos,
- Repetir para todos los bordes de la ventana
- ¡OJO!
- Los bordes de la ventana se pueden recorrer en cualquier orden
- Los vértices del polígono se han de recorrer e introducir en la lista en orden
(Los mismos vértices pueden pertenecer a polígonos distintos)
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:
- [(50,150),(200,50),(350,150),(350,300),(250,300),(200,250),(150,350),(100,250),(100,200)]
and clip it by the rectangle defined by the points:
- [(100,100),(300,100),(300,300),(100,300)]
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.)
Contents[hide] |
[edit]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;
- 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)
}
[edit]BBC BASIC
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
[edit]C
Most of the code is actually storage util routines, such is C. Prints out nodes, and writes test.eps file in current dir.
#include#include #include 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 250175 300 125 300 100 250 100 200 100 116.667 125 100 275 100 300 116.667 300 300 250 300test.eps written
[edit]C#
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 ////// This represents a line segment ///
private class Edge
{
public Edge(Point from, Point to)
{
this.From = from;
this.To = to;
}
public readonly Point From;
public readonly Point To;
}
#endregion
///
/// This clips the subject polygon against the clip polygon (gets the intersection of the two polygons)
///
///
/// Based on the psuedocode from:
/// http://en.wikipedia.org/wiki/Sutherland%E2%80%93Hodgman
///
/// Can be concave or convex
/// Must be convex
///
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
///
/// This iterates through the edges of the polygon, always clockwise
///
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
}
}
///
/// Returns the intersection of the two lines (line segments are passed in, but they are treated like infinite lines)
///
///
/// Got this here:
/// http://stackoverflow.com/questions/14480124/how-do-i-detect-triangle-and-rectangle-intersection
///
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");
}
///
/// Tells if the test point lies on the left side of the edge line
///
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:
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">
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); } ////// This is just a wrapper to the color converter (why can't they have a method off the color class with all /// the others?) ///
private static Color ColorFromHex(string hexValue)
{
if (hexValue.StartsWith("#"))
{
return (Color)ColorConverter.ConvertFromString(hexValue);
}
else
{
return (Color)ColorConverter.ConvertFromString("#" + hexValue);
}
}
#endregion
}
}
[edit]D
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.
[edit]Fortran
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
[edit]Go
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)
[edit]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
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 <- span=""> createPen Solid 3 blue penP <- span=""> 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)]
[edit]J
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
[edit]Java
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); } } }
[edit]JavaScript
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 (j in clipPolygon) { var cp2 = clipPolygon[j]; var inputList = outputList; outputList = []; s = inputList[inputList.length - 1]; //last on the input list for (i in inputList) { var 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[edit]Lua
No extra credit.
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)
[edit]MATLAB / Octave
%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
[edit]OCaml
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 ()
[edit]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
- 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)
No hay comentarios:
Publicar un comentario