## Polygon Challenge

General FreeBASIC programming questions.
duke4e
Posts: 717
Joined: Dec 04, 2005 0:16
Location: Varazdin, Croatia, Europe
Contact:
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
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs
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 getbothareasscreen 13''input data (anticlockwise: positive area, clockwise: negative area)const as integer n = 11dim as double x(0 to n - 1), y(0 to n - 1)x( 0) =  7: y( 0) = 11x( 1) =  2: y( 1) =  9x( 2) =  2: y( 2) =  8x( 3) =  3: y( 3) =  5x( 4) =  8: y( 4) =  5x( 5) =  5: y( 5) =  2x( 6) = 10: y( 6) =  2x( 7) =  9: y( 7) =  9x( 8) = 14: y( 8) =  4x( 9) = 18: y( 9) =  6x(10) = 16: y(10) = 10''dividing line pointsconst as double lx1 =  0, ly1 = 9const as double lx2 = 19, ly2 = 5const as double ldx_ = lx2 - lx1, ldy_ = ly2 - ly1 ''calculate direction vector of lineconst 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) ), 8line (5 * (lx1 + lx2), 199 - 5 * (ly1 + ly2) ) - (10 * lx2, 199 - 10 * ly2), 7#ifdef getbothareasdim as double area(0 to 1) = {0, 0}dim as double areaA, areaBdim as integer A, B#endif''useful vars for following loopdim as double x1, y1, x2, y2dim as double x1_, y1_, x2_, y2_dim as double iy, dy, mxdim 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)   #endifprint "Difference = " & csng(areadiff)if areadiff > 0 then print "Right area heavier" else print "Left area heavier"sleep` 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:
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
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 fencelineData    0, 9        ' reference fence post position Data    19, 5       ' secondary fence post position' specify the number of data points followingData    11' the n datapoints here for testing are approximately the original example' these points must be entered in anti-clockwise order around the polygonData    1,  7, 11     ' first pointData    2,  2,  9Data    3,  2,  8Data    4,  3,  5Data    5,  8,  5Data    6,  5,  2Data    7, 10,  2Data    8,  9,  9Data    9, 14,  4Data   10, 18,  6Data   11, 16, 10     ' last point'---------------------------------------------------------------------------Screen 12           ' Used by polygon plot routine during testingWindow (-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 DoubleDeclare Sub Poly_plot( As Integer, x() As Double, y() As Double, As Integer)'---------------------------------------------------------------------------Dim Shared As Double Tx, Ty ' used to translate by mapping subroutinesDim Shared As Double Sx, Sy ' also used in plot routineRead Tx, Ty     ' input position of the reference fence postRead Sx, Sy     ' input position of the secondary fence post'---------------------------------------------------------------------------' generate a unity vector of the fence line between the two postsDim As Double dx, dy, radiusdx = Sx - Txdy = Sy - Tyradius = Sqr(dx*dx+dy*dy)   ' distance between fence postsDim Shared As Double Ux, Uy ' unity vector used for rotateUx = dx / radius        ' convert to a unity length vectorUy = 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 nDim As Integer n, dummy ' dummy is an input data point sequence entry aidRead n                  ' get number of points used to specify the polygonDim As Integer Nmax = 2 * (n + 1)       ' estimate maximum storage neededDim As Double Px(n+1), Py(n+1)          ' polygon points' we need storage for the generated sub-polygonsDim As Integer Hi_n = 0, Lo_n = 0  ' count number of points in followingDim As Double Hi_x(Nmax), Hi_y(Nmax), Lo_x(Nmax), Lo_y(Nmax) ' sub-polygons'---------------------------------------------------------------------------' read in the n data pointsFor 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 pointsDim As Double newxPx(n+1) = Px(1) ' copy the first point to the end, close the polygonPy(n+1) = Py(1)map (Px(1), Py(1), head_x, head_y)  ' generate initial mapped positionFor 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 IfNext i'---------------------------------------------------------------------------' output informationPrint "Polygon =", Poly_area (n, Px(), Py())Dim As Double Hi_area = Poly_area (Hi_n, Hi_x(), Hi_y())Print "Hi_area =", Hi_areaDim As Double Lo_area = Poly_area (Lo_n, Lo_x(), Lo_y())Print "Lo_area =", Lo_areaPrint "total area =", Lo_area + Hi_areaPoly_plot(n, Px(), Py(), 12)  ' plot the polygon and fence lineSleep 1000Poly_plot(Hi_n, Hi_x(), Hi_y(), 14)  ' plot the polygon and fence lineSleep 1000Poly_plot(n, Px(), Py(), 12)  ' replot the polygon and fence linePoly_plot(Lo_n, Lo_x(), Lo_y(), 14)  ' plot the polygon and fence lineSleep'===========================================================================' 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 vectorEnd 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 backEnd 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 valuesEnd 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 iEnd 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:
Damn, you should really write some vector tutorial for FB wiki ;)
counting_pine
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs
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 = 11dim as double x(0 to n - 1), y(0 to n - 1)x( 0) =  7: y( 0) = 11x( 1) =  2: y( 1) =  9x( 2) =  2: y( 2) =  8x( 3) =  3: y( 3) =  5x( 4) =  8: y( 4) =  5x( 5) =  5: y( 5) =  2x( 6) = 10: y( 6) =  2x( 7) =  9: y( 7) =  9x( 8) = 14: y( 8) =  4x( 9) = 18: y( 9) =  6x(10) = 16: y(10) = 10''dividing line pointsdim as double lx1 =  0, ly1 = 9dim as double lx2 = 19, ly2 = 5dim as double ldx_, ldy_, ldx, ldyldx_ = lx2 - lx1 ''calculate direction vector of lineldy_ = ly2 - ly1ldx = ldx_ / sqr(ldx_ ^ 2 + ldy_ ^ 2) ''make it a unit vectorldy = 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) / 2line (10 * lx1,   199 - 10 * ly1  )- _     (10 * lx1_5, 199 - 10 * ly1_5), 8line (10 * lx1_5, 199 - 10 * ly1_5)- _     (10 * lx2,   199 - 10 * ly2  ), 7''useful vars for following loopdim as double AreaDiff = 0          ''total difference between two areasdim as double Area0, Area1          ''left (0) and right (1) areasdim as double Area(0 to 1) = {0, 0} ''totals of left and right areasdim as double x1_, y1_, x2_, y2_    ''the line points from the arraydim as double x1,  y1,  x2,  y2     ''line points, after rotationdim as double dx, dy                ''displacement of x2 from x1; y2 from y1dim as double iy                    ''point at interception of x axisdim 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:
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? 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, y2End 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
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
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:
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. So, this case should return 0 or something similar...
Richard
Posts: 2958
Joined: Jan 15, 2007 20:44
Location: Australia
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
Posts: 6173
Joined: Jul 05, 2005 17:32
Location: Manchester, Lancs
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.