Once upon a time we would have cut the polygon out of sheet material, cut along the fence line and then weighed or compared the two parts on a beam balance. Counting pixels is quicker but not as fast and accurate as vector maths solutions which take on the order of 50 usec for this example. Interestingly, balancing the polygon with the fence line along a knife edge could give a different result as we would be comparing moments not areas.
Counting_pine’s code measures the areas correctly. We get the same answers to the test examples. The area(a) and area(b) sub-totals are not both needed since they are only going to be subtracted later. Just accumulate one side while subtracting the other side from the running total. The sign of the final total tells which was greater. The rotation vector then does not need to be a unity vector because the scale is no longer relevant and the generated points will not be used. While I like the “mod n” to save copying the first to the end. I tend not to use it because I think in assembler and the mod function would slow down signal processing software too much. It constitutes a division within the inner loop. Unfortunately my RISC brain was unable to verify Counting_pine’s algorithm correctness and reliability by just reading his code. The use of #define inside a comment block to perform the interpolation function was the thing that threw me. I could not be sure where it would pop up or if the divisions contained would ever be zero. I felt uneasy about Counting_pine’s algorithm, too much like C. I think the level of abstraction was too deep for most readers. The use of Window(-) to scale the plot would have made it easier to play with. Maybe Counting_pine is too intelligent for some of us.
This requires that crossing points be available. That’s why I evaluated and kept them.
I wrote my code as a challenge and to ward off Alzheimer’s. To be certain the code was correct and worked I needed to prove each step. That is why it is written as a sequence of “hyphen lined” blocks, each understandable and testable by itself. The Chinese proverb “Teach a man to fish and you feed him for a lifetime” sums it up. It is always harder to teach how, than to do it once yourself. My aim was to teach by showing how it could be done, not to provide a free programming service. That is why I kept it simple, like me.
Since we know the fence line before we get the data, we could process the whole task as we read the data points in. We don’t even need to know how many points there are, just remember the first. We do not have to relocate the polygon, we just need to identify fence–edge intersections as we go.
So here is a new version that generates in-place sub-polygons without moving the input polygon. It uses my same tested algorithm but eliminates the blocks and uses just one loop. If you want the area comparison, the areas are there. It’s more flexible and useful than my original but the hidden mapping in this might throw some beginners. Again it is hopefully over commented.
Code: Select all
'===========================================================================
' Partition an irregular polygon with an infinite fence line. Generate
' two sub-polygons and compute areas. Copyright 2007, Richard. GNU free use.
'===========================================================================
' first specify the coordinates of two points on the fenceline
Data 0, 9 ' reference fence post position
Data 19, 5 ' secondary fence post position
' specify the number of data points following
Data 11
' the n datapoints here for testing are approximately the original example
' these points must be entered in anti-clockwise order around the polygon
Data 1, 7, 11 ' first point
Data 2, 2, 9
Data 3, 2, 8
Data 4, 3, 5
Data 5, 8, 5
Data 6, 5, 2
Data 7, 10, 2
Data 8, 9, 9
Data 9, 14, 4
Data 10, 18, 6
Data 11, 16, 10 ' last point
'---------------------------------------------------------------------------
Screen 12 ' Used by polygon plot routine during testing
Window (-1, 15)-(20, 0) ' values selected to display polygons
'
'---------------------------------------------------------------------------
' Where a polygon edge crosses the fence line, a new point is created. As
' data points are fed in, fence crossing points are detected and added to
' the sub-arrays as generated along with appropriate original data points.
'
'---------------------------------------------------------------------------
' Fence line has two posts, one a reference post, the other a secondary post
' If you stand at the reference post and look towards the secondary post
' then the polygon area on your left is the hi_side, to your right is the
' lo_side. If fence had been defined by y = mx + c then slope m could not
' be infinite. Whereas definition by two posts makes any fence possible.
'
'---------------------------------------------------------------------------
' How it works: If we reposition the polygon and fence together so that the
' fence lies along the +X-axis with the reference post at the origin then
' all points on or above the x-axis are part of the Hi-side sub-polygon and
' all points on or below the x-axis are part of the Lo-side sub-polygon.
' map() transforms from normal space to mapped space with fence on x-axis.
' unmap() transforms from mapped space back to normal space.
'
'---------------------------------------------------------------------------
' Understanding the math's: How vectors V(x,y) can effect points P(x,y).
' Addition of a vector translates the point by vector length and direction.
' A unity vector has a length of one. It is one radius of the unity circle.
' Multiplication by a unity vector rotates a point about the origin by
' the angle of the unity vector. Multiplication by a non-unity vector
' scales by vector length as it rotates. The conjugate of a vector is a
' vector with the negative of it's y-component.
' Multiplication by the conjugate of a vector rotates the point backwards.
' Relative to the origin, points are in effect just vectors.
'
'---------------------------------------------------------------------------
' Vectors are represented with complex numbers, (a + i b), where a is the
' real part and b the immaginary part. i squared is defined as minus one.
' For plane geometry the real component is the x value, the immaginary is y.
' Vector addition (a+ib) + (c+id) = ( (a+c) + i (b+d) )
' Vector multiplication, (a+ib) * (c+id) = ( (ac-bd) + i (ad+bc) )
'
'---------------------------------------------------------------------------
Declare Sub map ( As Double, As Double, Byref As Double, Byref As Double )
Declare Sub unmap ( As Double, As Double, Byref As Double, Byref As Double )
Declare Function Poly_area( As Integer,x()As Double,y()As Double) As Double
Declare Sub Poly_plot( As Integer, x() As Double, y() As Double, As Integer)
'---------------------------------------------------------------------------
Dim Shared As Double Tx, Ty ' used to translate by mapping subroutines
Dim Shared As Double Sx, Sy ' also used in plot routine
Read Tx, Ty ' input position of the reference fence post
Read Sx, Sy ' input position of the secondary fence post
'---------------------------------------------------------------------------
' generate a unity vector of the fence line between the two posts
Dim As Double dx, dy, radius
dx = Sx - Tx
dy = Sy - Ty
radius = Sqr(dx*dx+dy*dy) ' distance between fence posts
Dim Shared As Double Ux, Uy ' unity vector used for rotate
Ux = dx / radius ' convert to a unity length vector
Uy = dy / radius
' Tx, Ty become the translation vector to the origin
' Ux, Uy is the fence unity vector for rotation, (not the conjugate)
' all shared because they are used by the map, unmap and plot subroutines
'---------------------------------------------------------------------------
' identify number of polygon points n
Dim As Integer n, dummy ' dummy is an input data point sequence entry aid
Read n ' get number of points used to specify the polygon
Dim As Integer Nmax = 2 * (n + 1) ' estimate maximum storage needed
Dim As Double Px(n+1), Py(n+1) ' polygon points
' we need storage for the generated sub-polygons
Dim As Integer Hi_n = 0, Lo_n = 0 ' count number of points in following
Dim As Double Hi_x(Nmax), Hi_y(Nmax), Lo_x(Nmax), Lo_y(Nmax) ' sub-polygons
'---------------------------------------------------------------------------
' read in the n data points
For i As Integer = 1 To n
Read dummy, Px(i), Py(i) ' input array P(n)
Next i
'---------------------------------------------------------------------------
' start the feed of polygon points into the sausage machine
' The input array P(1..n) is not changed except that P(n+1) is overwriten.
Dim As Double tail_x, tail_y, head_x, head_y ' the two current mapped points
Dim As Double newx
Px(n+1) = Px(1) ' copy the first point to the end, close the polygon
Py(n+1) = Py(1)
map (Px(1), Py(1), head_x, head_y) ' generate initial mapped position
For i As Integer = 2 To n+1
tail_x = head_x ' prepare to advance to next point
tail_y = head_y ' avoid re-computing mapping of head as tail
map (Px(i), Py(i), head_x, head_y) ' generate next mapped position
' check to see if the tail point is part of either sub-ploygon
If tail_y >= 0 Then ' add point to Hi sub-polygon
Hi_n = Hi_n + 1
Hi_x(Hi_n) = Px(i-1)
Hi_y(Hi_n) = Py(i-1)
End If
If tail_y <= 0 Then ' add point to Lo sub-polygon
Lo_n = Lo_n + 1
Lo_x(Lo_n) = Px(i-1)
Lo_y(Lo_n) = Py(i-1)
End If
' now we test to see if the line between head and tail crosses x-axis
If (Head_y * Tail_y) < 0 Then ' if line jumped the x-axis then create
' calculate position of new point to insert in both sub-polygons
dx = Head_x - Tail_x ' dx, dy previously declared
dy = Head_y - Tail_y ' must be non-zero as it crosses the x-axis
newx = Tail_x - Tail_y * dx / dy ' linear interpolate with x = 0.
Hi_n = Hi_n + 1 ' add point to Hi sub-polygon
unmap( newx, 0 , Hi_x(Hi_n), Hi_y(Hi_n) )
Lo_n = Lo_n + 1 ' add point to Lo sub-polygon
unmap( newx, 0 , Lo_x(Lo_n), Lo_y(Lo_n) )
End If
Next i
'---------------------------------------------------------------------------
' output information
Print "Polygon =", Poly_area (n, Px(), Py())
Dim As Double Hi_area = Poly_area (Hi_n, Hi_x(), Hi_y())
Print "Hi_area =", Hi_area
Dim As Double Lo_area = Poly_area (Lo_n, Lo_x(), Lo_y())
Print "Lo_area =", Lo_area
Print "total area =", Lo_area + Hi_area
Poly_plot(n, Px(), Py(), 12) ' plot the polygon and fence line
Sleep 1000
Poly_plot(Hi_n, Hi_x(), Hi_y(), 14) ' plot the polygon and fence line
Sleep 1000
Poly_plot(n, Px(), Py(), 12) ' replot the polygon and fence line
Poly_plot(Lo_n, Lo_x(), Lo_y(), 14) ' plot the polygon and fence line
Sleep
'===========================================================================
' Functions and Subs follow
'===========================================================================
Sub map( in_x As Double, in_y As Double,_
Byref out_x As Double, Byref out_y As Double )
' transform from normal space to mapped space
Dim As Double x, y ' temporary
x = in_x - Tx ' translate to origin
y = in_y - Ty ' then do the rotate
out_x = Ux * x + Uy * y ' complex multiply by the conjugate Ux, -Uy
out_y = Ux * y - Uy * x ' of the unity vector
End Sub
'===========================================================================
Sub unmap( in_x As Double, in_y As Double,_
Byref out_x As Double, Byref out_y As Double )
' transform from mapped space to normal space
out_x = Tx + Ux * in_x - Uy * in_y ' complex multiply by the unity
out_y = Ty + Ux * in_y + Uy * in_x ' vector Ux,Uy then translate back
End Sub
'===========================================================================
Function Poly_area (n As Integer, x() As Double, y() As Double) As Double
'-----------------------------------------------------------------------
' Evaluate area of an irregular polygon of n points indexed from 1 to n
' the arrays x() and y() need to be sized at least one longer than n
' input data points must be in the correct order.
' anticlockwise -> positive area, clockwise -> negative area
' if edge lines cross then the included crossed area subtracts
'-----------------------------------------------------------------------
' Warning this changes the referenced array(n+1) data at its source
x(n + 1) = x(1) ' duplicate the first point onto the end so that
y(n + 1) = y(1) ' each polygon edge has a point at each end
'-----------------------------------------------------------------------
Dim As Double midx, dy, area = 0 ' each edge joins two points, evaluate
For i As Integer = 1 To n ' area between y axis and that edge
midx = (x(i + 1) + x(i)) ' twice midpoint of x values
dy = y(i + 1) - y(i) ' difference of y values
area = area + dy * midx ' accumulate areas of all point pairs
Next i
Poly_area = area / 2 ' correct for twice x midpoint values
End Function
'===========================================================================
Sub Poly_plot(n As Integer, x() As Double, y() As Double, colour As Integer)
Line(Tx, Ty)-( Sx, Sy), colour + 1
' Warning this changes the referenced array(n+1) data at its source
x(n + 1) = x(1)
y(n + 1) = y(1)
For i As Integer = 1 To n
Line (x(i),y(i))-(x(i+1),y(i+1)), colour
Circle (x(i),y(i)), 0.25 , colour
Next i
End Sub
'===========================================================================
' E N D o f P R O G R A M
'===========================================================================
