Polygon clipping by a line

[krtxmrtz]

(I wonder if this would be more appropriate for the game development codebank or the math forum).

This is a subroutine for clipping a general (convex or concave) polygon by a (infinite) straight line defined by 2 points.

The polygon is a closed contour with N vertices, such that points 1 and N are the same. A polygon UDT is used (see code below) which includes the coordinates of the vertices, the number of branches and pointers to the first point of each branch. Initially the polygon is supposed to have a single branch, but as a result of clipping it may be split into a number of individual closed contours.

The UDT definitions follow below. The subroutine comes next in a separate post as it's too long to fit in one post. Two screen shots show a demo polygon just before and after the clipping operation. A demo project is also included.

'Data type for the polygon vertices
Private Type punkt
X As Single
Y As Single
End Type
'Data type for the entire polygon
Private Type Polygon
vertex() As punkt
n As Integer 'Total number of vertices (first and last are the same)
up() As Integer 'Pointers to the vertices where the polygons start:
'when the polygon is plotted, we must go to them with "pen up"
nups As Integer 'Number of sub-polygons (branches): initially 1
'but after the polygon has been clipped, new polygons may be created
End Type

[krtxmrtz]

Private Sub Clip(plgn As Polygon, p() As punkt, lowerhalf As Boolean)
'
'*********************************************************
'This subroutine takes a general polygon (may be concave)
'and clips it by a straight line
'The boolean variable "lowerhalf" selects which half is
'to be clipped
'
'The polygon must not intersect itself (the code does not
'allow for this possibility, though it shouldn't be too
'difficult to add some checking code to the calling program)
'
'The algorithm has been developed from that outlined in
'the book: "Algorithms for Graphics and Image Processing"
'by Theo Pavlidis (Computer Science Press, 1982)
'ISBN: 0-914894-65-X
'*********************************************************
'
Dim A As Single, B As Single, C As Single
Dim U1 As Single, U2 As Single
Dim i As Integer, j As Integer, k As Integer
Dim k0 As Integer, jcumm As Integer
'Initial value of n (n will vary as the polygon is being clipped)
Dim nini As Integer
Dim kold As Integer
Dim u As Single, v As Single
Dim xa As Single, ya As Single, wa As Single, cro As Single
Dim intr As punkt
Dim auxpol As Polygon
'Array of pointers to the next vertex of the polygon
'Initially, nxt(i) = i+1 for all i (except nxt(plgn.n) = 1) but
'as the clipping operation proceeds it will have to be updated
Dim nxt() As Integer
Dim tempo() As Integer
Dim dumbo As Integer
Dim thesign As Integer
Dim dummy() As Single, dumdum As Single
Dim zeichen As Integer

'Keep a copy of the initial number of vertices
nini = plgn.n

'Initialize the array nxt
ReDim nxt(1 To plgn.n)
For i = 1 To plgn.n
Select Case i
Case plgn.n
nxt(i) = 1
Case Else
nxt(i) = i + 1
End Select
Next

'Straight line equation:
'Ax + By + C = 0
'Calculate the parameters A, B and C from the
'coordinates of both ends of the line segment
If p(1).X = p(2).X Then
'Vertical line
A = 1
B = 0
C = -p(1).X
Else
'Any other slope
A = (p(2).Y - p(1).Y) / (p(2).X - p(1).X)
B = -1
C = p(1).Y - A * p(1).X
End If

'Assign the value of the parameter zeichen
'according to which of the 2 halves of the
'the polygon is to be clipped
If lowerhalf Then
zeichen = 1
Else
zeichen = -1
End If
'For the first vertex of the polygon, calculate the quantity U1
'which represents its position relative to the line:
'If U1>o the point lies below the line (or at its left side if
'the line is vertical)
'If U1 is negative the point is located above it (or to the right)
U1 = A * plgn.vertex(1).X + B * plgn.vertex(1).Y + C
If U1 * zeichen <= 0 Then nxt(1) = 0
'Loop over the other sides of the polygon
For i = 2 To plgn.n
'Calculate the position of vertex i
U2 = A * plgn.vertex(i).X + B * plgn.vertex(i).Y + C
'Check the position of the previous vertex and of the current vertex
If U1 * zeichen <= 0 And U2 * zeichen <= 0 Then
'Set the current vertex pointer to 0 if both vertices
'are located on the side to be clipped
nxt(i) = 0
Else
'Either the vertices are on the opposite sides or
'both lie on the side opposite to that to be clipped
If U1 * zeichen < 0 Or U2 * zeichen < 0 Then
'If on opposite sides, calculate the coordinates of the
'interception between the line and the polygon side from
'the previous vertex to the current one
u = plgn.vertex(i - 1).X - plgn.vertex(i).X
v = plgn.vertex(i - 1).Y - plgn.vertex(i).Y
cro = plgn.vertex(i - 1).X * plgn.vertex(i).Y - plgn.vertex(i - 1).Y * plgn.vertex(i).X
xa = -(u * C + B * cro)
ya = -(v * C - A * cro)
wa = v * B + u * A
If wa <> 0 Then
'Coordinates of the interception
intr.X = xa / wa
intr.Y = ya / wa
Else
MsgBox "Error"
Exit Sub
End If
'The interception point is to be counted as a new
'point, so increment the total number of points
plgn.n = plgn.n + 1
ReDim Preserve plgn.vertex(1 To plgn.n)
ReDim Preserve nxt(1 To plgn.n)
plgn.vertex(plgn.n).X = intr.X
plgn.vertex(plgn.n).Y = intr.Y
If U1 * zeichen < 0 Then
'If the first point is on the side to be clipped
'(the second point is on the opposite side) then
'set its pointer to point to the second point
nxt(plgn.n) = i
Else
'Otherwise, set the previous point's pointer to point
'to the interception point and the second point's
'pointer to zero (as it lies on the side to be clipped)
nxt(i - 1) = plgn.n
nxt(i) = 0
nxt(plgn.n) = 0
End If
End If
'Set the current second vertex as the first
'for the comparison in the next loop
U1 = U2
End If
Next

TO BE CONTINUED IN MY NEXT POST

[krtxmrtz]

REST OF THE CODE
'An array "tempo" is defined with a number of elements
'equal to the number of added points, which will be used
'below to rearrange the array "nxt"
'Here "tempo" is initialized
ReDim tempo(nini + 1 To plgn.n)
For i = nini + 1 To plgn.n
tempo(i) = i
Next

'Determine the contour orientation through the calculation
'of the signed area
'If it's negative, the contoun is clockwise oriented
'If it's positive, it's counter-clockwise oriented
thesign = Int(Sgn(SignedPolygonArea(po)))

'Array dummy is a copy of the Y coordinates (or of the
'X coordinates, if the line is horizontal) of the new
'added vertices
'It will be used for the rearrangement of array tempo
'Because the actual polygon vertices can't be rearranged
'we need to reaarange a copy of them in order for the
'sorting algorithm to work properly
ReDim dummy(nini + 1 To plgn.n)
For i = nini + 1 To plgn.n
Select Case A
Case 0
'Horitzontal line: compare the x coordinates
dummy(i) = plgn.vertex(i).X
Case Else
'Any other slope: compare the y coordinates
dummy(i) = plgn.vertex(i).Y
End Select
Next

'Sort array tempo by increasing (or decreasing) values of
'the y (or x if the line is horizontal) coordinates of the
'corresponding added polygon vertices
For i = nini + 1 To plgn.n - 1
For j = i + 1 To plgn.n
'Use "thesign" to take into account the contour orientation
If thesign * dummy(i) > thesign * dummy(j) Then
dumbo = tempo(i)
tempo(i) = tempo(j)
tempo(j) = dumbo
dumdum = dummy(i)
dummy(i) = dummy(j)
dummy(j) = dumdum
End If
Next
Next

'Elements of array nxt are taken in consecutive pairs
'One (and only one) of each pair will be zero; this one
'is then set to point to the other element
'to the new added points by increasing values of y
'(or of x is the line is horizontal)
For i = 1 To (plgn.n - nini) \ 2
j = nini + 2 * i
If nxt(tempo(j - 1)) = 0 Then
nxt(tempo(j - 1)) = tempo(j)
Else
nxt(tempo(j)) = tempo(j - 1)
End If
Next

'Use a new polygon structure to make up the clipped polygon
'Initialize the number of branches to zero
plgn.nups = 0
'Initialize counter for the vertices in a specific branch
j = 0
'Initialize counter for the total number of vertices
jcumm = 0
'Allow for plgn.n points, the total number
'of polygon points including the added ones
ReDim auxpol.vertex(1 To plgn.n)
For i = 1 To plgn.n
'Vertices with nxt(i)=0 are located
'on the side to be clipped and are
'neglected
If nxt(i) <> 0 Then
'Vertex is on the side to be kept
If j = 0 Then
'First point of the current branch
'Update the branch counter
plgn.nups = plgn.nups + 1
'Update the array of pointers to
'the first point of each branch
ReDim Preserve plgn.up(1 To plgn.nups)
'Assign value of pointer
plgn.up(plgn.nups) = jcumm + 1
'Initialize k to the first
'point of the branch
k = i
'Set k0 to the index of the
'first point of the branch
k0 = i
End If
Do
j = j + 1
auxpol.vertex(jcumm + j) = plgn.vertex(k)
'Save copy of k
kold = k
'Update index k
k = nxt(k)
'Make pointer for the current index = zero so
'it will be skipped past in the following loops
'(kold must be used as k has just been modified)
nxt(kold) = 0
'Make up the branch by collecting all points
'until the first one is reached
Loop While k <> k0
'Add one more point to the current branch and make it
'the same as the first one (i.e. close the contour)
j = j + 1
auxpol.vertex(jcumm + j) = plgn.vertex(k0)
'Increment total point counter
jcumm = jcumm + j
'Initialize the number of points for the next branch
j = 0
End If
Next

'Set total number of points
plgn.n = jcumm
'And re-build the polygon after clipping
ReDim plgn.vertex(1 To plgn.n)
For i = 1 To plgn.n
plgn.vertex(i) = auxpol.vertex(i)
Next
'Free some memory
Erase auxpol.vertex

End Sub

[krtxmrtz]

Finally, here's this function for calculationg the area of the polygon, called from the above clipping routine:

Private Function SignedPolygonArea(pol As Polygon) As Single
'
'*********************************************
'This function has been taken from
'http://vb-helper.com/howto_polygon_area.html
'and has been only slightly modified to allow
'for the "polygon" data type
'
'In the current project it is used ONLY to
'determine the orientation (CW or CCW) of
'the only branch of the original polygon
'*********************************************
'
' Return the polygon's area in "square" pixels.
' Add the areas of the trapezoids defined by the
' polygon's edges dropped to the X-axis. When the
' program considers a bottom edge of a polygon, the
' calculation gives a negative area so the space
' between the polygon and the axis is subtracted,
' leaving the polygon's area. This method gives odd
' results for non-simple polygons.
'
' The value will be negative if the polyogn is
' oriented clockwise.
Dim vk As Integer
Dim area As Single
' Get the areas.
For vk = 1 To pol.n - 1
area = area + _
(pol.vertex(vk + 1).X - pol.vertex(vk).X) * _
(pol.vertex(vk + 1).Y + pol.vertex(vk).Y) / 2
Next
' Return the result.
SignedPolygonArea = area
End Function