Polygon Challenge

General FreeBASIC programming questions.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Postby duke4e » Apr 26, 2007 21:13

counting_pine wrote:@duke4e:
So, now this thread has several suggestions on how to proceed, and a couple of solutions in code. Hopefully something there will help you.

If you don't mind me asking though, what do you need the algorithm for?

I'm prototyping for a QIX-kind of game. I wanted to take "modern" approach with freedom for map drawing and cutting. The only problem was that I've stucked on main part of idea - determinating the smaller part of map/polygon.

D.J.Peters wrote:Why i got a german message via email for this tread here?

Could be that someone of forum admins is experimenting with subscription to threads? Could be related to this thread http://www.freebasic.net/forum/viewtopic.php?t=8090

And last, but not the least - THANK YOU RICHARD FOR MAKING THIS POLYGON CHALLENGE SOLVABLE. Your code is very well commented and very useful to me. If I ever finish this game, you can count on your name on credits. Also, thanks everyone else who was involved in this thread!
counting_pine
Site Admin
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Postby counting_pine » Apr 27, 2007 2:07

I have a new solution for the problem. (Better late than never...)

This code has been much helped by (and originally written by) Richard, so thanks to him for showing us the simple genius of the 'trapezoid' method for calculating the area, and rotating the shape to get away from the "horizontal dividing line" problem.

Actually having to split the polygon in two, though, to be honest, seems like too much work to me.

I just keep track of the two areas, adding or subtracting to whichever is appropriate, depending on where the edge is, relative to the dividing line.
There is also the option to just keep track of the difference, like in my first version, because the math is a bit easier. The difference alone might or might not be enough for you, I don't know.

To match stakes with Richard's version, I've tried to comment well, and added drawing code.

Code: Select all

''comment this out to just get the difference of the areas (the math is easier)
#define getbothareas

screen 13

''input data (anticlockwise: positive area, clockwise: negative area)
const as integer n = 11
dim as double x(0 to n - 1), y(0 to n - 1)

x( 0) =  7: y( 0) = 11
x( 1) =  2: y( 1) =  9
x( 2) =  2: y( 2) =  8
x( 3) =  3: y( 3) =  5
x( 4) =  8: y( 4) =  5
x( 5) =  5: y( 5) =  2
x( 6) = 10: y( 6) =  2
x( 7) =  9: y( 7) =  9
x( 8) = 14: y( 8) =  4
x( 9) = 18: y( 9) =  6
x(10) = 16: y(10) = 10

''dividing line points
const as double lx1 =  0, ly1 = 9
const as double lx2 = 19, ly2 = 5

const as double ldx_ = lx2 - lx1, ldy_ = ly2 - ly1 ''calculate direction vector of line
const as double ldx = ldx_ / sqr(ldx_ ^ 2 + ldy_ ^ 2), ldy = ldy_ / sqr(ldx_ ^ 2 + ldy_ ^ 2) ''normalize it

''draw dividing line (pointing in light gray direction)
''(y coords reversed, so y=0 is at bottom of screen.)
line (10 * lx1, 199 - 10 * ly1) - (5 * (lx1 + lx2), 199 - 5 * (ly1 + ly2) ), 8
line (5 * (lx1 + lx2), 199 - 5 * (ly1 + ly2) ) - (10 * lx2, 199 - 10 * ly2), 7


#ifdef getbothareas
dim as double area(0 to 1) = {0, 0}
dim as double areaA, areaB
dim as integer A, B
#endif

''useful vars for following loop
dim as double x1, y1, x2, y2
dim as double x1_, y1_, x2_, y2_
dim as double iy, dy, mx
dim as double areadiff = 0


'' this is where the magic happens...

for i as integer = 0 to n - 1
   
    ''this loop will find the area contained between each edge and the dividing line
    ''the points are rotated so the dividing line becomes the line x = 0
   
    x1_ = x(i): x2_ = x((i + 1) mod n)
    y1_ = y(i): y2_ = y((i + 1) mod n)
   
    line (10 * x1_, 199 - 10 * y1_) - _
         (10 * x2_, 199 - 10 * y2_), _
         56 + int(24 * i / n) ''draw edge (in red-violet order)
   
    ''rotate all the points
   
    x1 = ( ldy * x1_ - ldx * y1_) - (ldy * lx1 - ldx * ly1)
    x2 = ( ldy * x2_ - ldx * y2_) - (ldy * lx1 - ldx * ly1)
   
    y1 = ( ldx * x1_ + ldy * y1_)
    y2 = ( ldx * x2_ + ldy * y2_)
   
   
    ''if the edge intersects the dividing line, there will be two triangles made by it
   
    ''  |-/
    ''  |/
    '' /|
    ''/-|
   
    ''otherwise, there will be a tapezoid
   
    '' |----/
    '' |---/
    '' |--/
   
    ''dy will be the height of the trapezoid, or the height of the two triangles
    dy = y2 - y1
   
    #ifdef getbothareas
       
        if sgn(x1) * sgn(x2) >= 0 then
            ''point are same side of line: trapezoid
            if x1 + x2 > 0 then A = 1 else A = 0 ''which shape is the area added to?
           
            mx = (x1 + x2) / 2 ''area of trapezium:
            area(A) += mx * dy
        else
            ''points are opposite sides of line: two triangles
            if x1 > 0 then A = 1: B = 0 else A = 0: B = 1 ''which shape gets which area?
            iy = y1 + dy * x1 / (x1 - x2) ''find intersection point
           
            AreaA = x1 * (iy - y1) / 2 ''areas of two triangles:
            AreaB = x2 * (y2 - iy) / 2
            area(A) += AreaA
            area(B) += AreaB
        end if
       
    #else
       
        if sgn(x1) * sgn(x2) >= 0 then
            ''point are same side of line: trapezoid
            mx = abs(x1 + x2) / 2
            areadiff += mx * dy
        else
            ''points are opposite sides of line: two triangles
            mx = (x1 ^ 2 + x2 ^ 2) / (2 * abs(x2 - x1) )
            areadiff += mx * dy
        end if
       
    #endif
   
next i

''Report results:

#ifdef getbothareas
   
    print "Area to the right =  " & csng(area(1))
    print "Area to the left = " & csng(area(0))
    areadiff = area(1) - area(0)
   
#endif
print "Difference = " & csng(areadiff)
if areadiff > 0 then print "Right area heavier" else print "Left area heavier"

sleep

Image
Last edited by counting_pine on Apr 29, 2007 19:33, edited 1 time in total.
roook_ph
Posts: 402
Joined: Apr 01, 2006 20:50
Location: philippines
Contact:

Postby roook_ph » Apr 28, 2007 14:50

Dont know that much geometry or trigonometry but a simple basic code could solve that problem. Ill describe the solution in details . Paint the below red line inside the polygon with blue (or any color) paint the upper inside of the polygon with yellow. Now run pixel left to right top to bottom like so put it in array the bigger array has the bigger polygon
Richard
Posts: 2958
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Apr 29, 2007 0:47

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.

duke4e wrote:
I want to check which part of polygon is smaller.
It would be useful if program could also report all vertices which belong to that part of polygon.

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
'===========================================================================
Last edited by Richard on May 01, 2007 0:51, edited 1 time in total.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Postby duke4e » Apr 29, 2007 1:25

Damn, you should really write some vector tutorial for FB wiki ;)
counting_pine
Site Admin
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Postby counting_pine » Apr 29, 2007 16:02

I have tried to take some of your comments to heart. My aim was not to provide a 'free programming service', but to illustrate an idea, as best I could.

The mod, as much of my code, was written for brevity. I have removed it and replaced it with an if. It could be optimised more for speed, but I don't think it would help the simplicity.

I posted the first working solution to the problem, based on your first code, and only calculating the difference. I received no comments on it at all. I still don't understand why. I ifdeffed in some code to calculate both areas, in case that was the problem. I'm sorry if you found this confusing, but I hoped it would be simple enough for people to "preprocess" the code themselves and remove the unused code, if necessary.

In any case, I've removed the ifdefs and improved the code I added. I've expanded it a bit, added more comments, and made it clear why there are no divisions by zero. I've left the 'difference' version pretty much as it is.
Neither piece of code physically splits the polygon. I think this is a step that overcomplicates the procedure.

You have made some fantastic contributions to this thread, in ideas and code, and have produced multiple, very well written, working solutions, but if both algorithms are fully understood and well optimised, I think my solution will work better.

Code: Select all

screen 13

''input data (anticlockwise: positive area, clockwise: negative area)
const as integer n = 11
dim as double x(0 to n - 1), y(0 to n - 1)

x( 0) =  7: y( 0) = 11
x( 1) =  2: y( 1) =  9
x( 2) =  2: y( 2) =  8
x( 3) =  3: y( 3) =  5
x( 4) =  8: y( 4) =  5
x( 5) =  5: y( 5) =  2
x( 6) = 10: y( 6) =  2
x( 7) =  9: y( 7) =  9
x( 8) = 14: y( 8) =  4
x( 9) = 18: y( 9) =  6
x(10) = 16: y(10) = 10

''dividing line points
dim as double lx1 =  0, ly1 = 9
dim as double lx2 = 19, ly2 = 5
dim as double ldx_, ldy_, ldx, ldy

ldx_ = lx2 - lx1 ''calculate direction vector of line
ldy_ = ly2 - ly1

ldx = ldx_ / sqr(ldx_ ^ 2 + ldy_ ^ 2) ''make it a unit vector
ldy = ldy_ / sqr(ldx_ ^ 2 + ldy_ ^ 2)

''draw dividing line (pointing in light gray direction)
''all coordinates ares caled by 10.
''y coords are reversed, so y=0 is at bottom of screen.
dim as double lx1_5 = (lx1 + lx2) / 2, ly1_5 = (ly1 + ly2) / 2

line (10 * lx1,   199 - 10 * ly1  )- _
     (10 * lx1_5, 199 - 10 * ly1_5), 8

line (10 * lx1_5, 199 - 10 * ly1_5)- _
     (10 * lx2,   199 - 10 * ly2  ), 7

''useful vars for following loop

dim as double AreaDiff = 0          ''total difference between two areas

dim as double Area0, Area1          ''left (0) and right (1) areas
dim as double Area(0 to 1) = {0, 0} ''totals of left and right areas

dim as double x1_, y1_, x2_, y2_    ''the line points from the array
dim as double x1,  y1,  x2,  y2     ''line points, after rotation
dim as double dx, dy                ''displacement of x2 from x1; y2 from y1
dim as double iy                    ''point at interception of x axis
dim as double mx                    ''average width of trapezoid


''here is the main loop that gets the edges, draws them and finds the areas:

for i as integer = 0 to n - 1
   
    ''this loop will find the area contained between each edge and the dividing line
    ''the points are rotated so the dividing line becomes the line x = 0
   
    'get the points for edge n
    x1_ = x(i)
    y1_ = y(i)
    if i + 1 = n then
        x2_ = x(0)
        y2_ = y(0)
    else
        x2_ = x(i + 1)
        y2_ = y(i + 1)
    end if
   
    line (10 * x1_, 199 - 10 * y1_)- _ ''draw edge (in red-violet order)
         (10 * x2_, 199 - 10 * y2_), _
         56 + int(24 * i / n)
   
    ''rotate all the points
   
    x1 = ( ldy * x1_ - ldx * y1_) - (ldy * lx1 - ldx * ly1)
    x2 = ( ldy * x2_ - ldx * y2_) - (ldy * lx1 - ldx * ly1)
   
    y1 = ( ldx * x1_ + ldy * y1_)
    y2 = ( ldx * x2_ + ldy * y2_)
   
   
    ''if the edge intersects the dividing line, there will be two triangles made by it
   
    ''   |-/
    ''   |/
    ''  /|
    '' /-|
    ''/--|
   
    ''otherwise, there will be a tapezoid
   
    '' |----/
    '' |---/
    '' |--/
   
    ''dy will be the height of the trapezoid, or the height of the two triangles
    ''when dy is positive, the areas on the right of the line will be positive
    ''and areas on the left will be negative.
    ''things will be the other way round when dy is negative.
    dy = y2 - y1
   
   
    ''find the total of both areas - one will be subtracted from the other at the end
    ''if you want to skip this bit, comment out or remove the whole of the following IF block:
   
    if x1 >= 0 and x2 >= 0 then ''both points are to the right of the line
       
        mx = (x1 + x2) / 2  ''average width of trapezoid (positive)
        Area1 = mx * dy     ''area of trapezoid
        Area(1) += Area1    ''add area of trapezoid to right hand area
       
    elseif x1 <= 0 and x2 <= 0 then ''both points are to the left of the line
       
        mx = (x1 + x2) / 2  ''average width of trapezoid (negative)
        Area0 = mx * dy     ''area of trapezoid
        Area(0) += Area0    ''add area of trapezoid to left hand area
       
    else
        ''points are opposite sides of the line: two triangles
       
        dx = x2 - x1 ''displacement of x2 from x1
        ''(never zero because points are opposite sides of the line)
       
        iy = y1 - x1 * dy / dx ''point where edge intercepts x axis
       
        if x1 > 0 then ''point 1 on right (so point 2 on left)
           
            Area1 = x1 * (iy - y1) / 2 ''right hand area is first triangle
            Area0 = x2 * (y2 - iy) / 2 ''left hand area is second triangle
           
            Area(0) += Area0
            Area(1) += Area1
           
        else ''point 1 on left (so point 2 on right)
           
            Area0 = x1 * (iy - y1) / 2 ''left hand area is first triangle
            Area1 = x2 * (y2 - iy) / 2 ''right hand area is second triangle
           
            Area(0) += Area0
            Area(1) += Area1
           
        end if
       
    end if
   
   
    ''find just the difference in areas and keep track of the total
    ''the math is basically the same, but it is simplified a lot.  If you
    ''understand the previous section, you should be able to derive this on your own
   
    if (x1 >= 0 and x2 >= 0) or (x1 <= 0 and x2 <= 0) then ''points are same side of line: trapezoid
        mx = abs(x1 + x2) / 2
        AreaDiff += mx * dy
    else ''points are opposite sides of line: two triangles
        dx = abs(x2 - x1)
        mx = (x1^2 + x2^2) / (2 * dx)
        AreaDiff += mx * dy
    end if
   
next i

''Report results:

print "Area to the right = " & csng(Area(1))
print "Area to the left =  " & csng(Area(0))
print "Difference (method 1) = " & csng(Area(1) - Area(0))

print "Difference (method 2) = " & csng(AreaDiff)
if AreaDiff > 0 then print "Right area heavier" else print "Left area heavier"

sleep
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Postby duke4e » Apr 29, 2007 22:32

Driven by sucess of this thread, and some very useful code by Richard and counting_pine, I must ask the (probably) final question. How to deal with polygons which have holes?
Image

My first guess it that "normal vectors" whould be out, since

Code: Select all

x(n + 1) = x(1)
y(n + 1) = y(1)
wouldn't make much sense. And something similar to code I've posted page before whould be needed - to deal with 2 vector pairs which reprecent a line.

Code: Select all

Type Wall_Type
    As Single x1, y1, x2, y2
End Type


Now, the main problem whould be how whould program recognize a hole in polygon?

Since the start of this thread, I've learned a massive ammount of new and very useful stuff and tricks, but I'm still "too new" in this vector/polygonal mess to figure this out. So, any advices where to begin?
Richard
Posts: 2958
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Apr 29, 2007 22:45

Counting_pine wrote:
I received no comments on it at all. I still don't understand why.

Although you did post the first working solution to the problem, you received no comments when you originally posted your algorithm because no one trusted or understood how it worked. It was too advanced for the readership. You overestimated the math’s skills of your readers. They could not cut and paste or modify it for their own applications. At the time I was still working on my algorithm and so I deliberately avoided studying or commenting on yours. I waited until I had finished mine and had time to test and understand yours. I really needed a proven reference that generated the sub-polygons to test against.

I am not going to compete with you for a “best” algorithm because we have different aims, we are in different events. You have focussed on comparison of the areas, I have focussed on the generation of two sub-polygons. We are both the leader, but each in our own field. Duke4e mentioned both requirements in the original post.

Your latest post is much more readable and therefore understandable. The removal of #define has made a big difference.
Richard
Posts: 2958
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Apr 29, 2007 22:56

Holes in polygons are separate polygons themselves. Holes subtract from the outer polygon area.
Run the algorithm on the outer polygon, then run it again on the hole, using the same fence line. Subtract the hole areas from the outer areas. There you have it.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:

Postby duke4e » Apr 30, 2007 0:17

Damn, you're right! Such a simple solution! Thanks!

edit: I've found a flaw in both codes - both programs assume that fence has infinite lenght. While they should actually have some kind of recognition of it's lenght.
Image

So, this case should return 0 or something similar...
Richard
Posts: 2958
Joined: Jan 15, 2007 20:44
Location: Australia

Postby Richard » Apr 30, 2007 10:36

There is no flaw in either code, you just moved the goal posts. You had not fully specified the problem. If the fence line has ends then your fence posts can mark them. The question that then rises is how to determine if the fence post is inside or outside the polygon. You have the maths needed to solve that riddle. This is an age old exercise from computer graphics.

1. You can understand how to solve it like this. First translate and rotate all the polygon points to put the fence post in question at the origin with the fence along the x-axis. You know the length of the fence, (it could still be infinite), find its end x-value on the positive x-axis. This could be your other fence post. Then go through all the polygon edges and count the number that cross the x-axis between the origin and the far end of the fence. If this count is odd then the point lies within the polygon. Even says that it crosses the polygon. All the math's for this are present in my first example, (Posted: Apr 25, 2007 13:03).

2. Once you can tell if the end fence post point is within the polygon you have the math's needed to tell if a polygon has a hole in it. Take any point from the hole and test it like the end fence post in the previous example. There are exceptions to this rule which can get very messy to test, consider:

3. Now if two polygons partly overlap. That is a problem to think about!
counting_pine
Site Admin
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs

Postby counting_pine » Apr 30, 2007 16:37

The latest flaw is inherent in my algorithm - it copes fine with holes, the same way Richard's does (as I mentioned somewhere in my first post), but I think a finite dividing line unavoidably requires splitting up the polygon.
Mine doesn't split polygons, but instead works out where each edge is in relation to the line. If some of the edges move away from the line boundaries, it has no way of working out which side of the line the polygon is on.

My algorithm (though it may not have seemed it) was a simplification of the one Richard subsequently used, and is designed to work effectively for the initial problem. But, as is evident, there are some more general problems that Richard's code can be modified to work with, but mine can't.

Thanks, Richard, for your constructive criticism, and I hope someone finds my algorithm useful. If you or anyone else has any questions about it, I will try to answer them as clearly as I can.

Return to “General”

Who is online

Users browsing this forum: MSN [Bot] and 1 guest