[Algorithms] Weiler Atherton Implementation (2D polygon clipping)

Discussion:

p***@free.fr

2005-05-26 11:23:15 UTC

Hello,

I am implementing the Weiler Atherton algorithm as it really matches my needs
for some physic dev:
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes

You can find the original article here :
http://www.cs.drexel.edu/~david/Classes/CS430/HWs/p214-weiler.pdf

The general cases works really good but it is getting harder with the
pathological cases when some intersections points happens to be the same as
some polygon points:
- 2 non parallel edges intersect on one end of one edge
- 2 non parallel edges intersect on one end of one edge which is also an end of
the other edge
- 2 parallel edges intersect on two points which are end of edges

I tried different ways to treat them (insert intersection nodes or not) but
nothing really works for the moment.

The paper is quite "vague" about it :
"If care is taken in placement of intersections where the subject and clip
polygon contours are identical in the x-y plane, no degenerate polygons will be
produced by the clipping process"

Then the text makes two assertions :
"These two types of intersections will be found to alternate along any given
contour and the number of intersections will always be even"

I cant satisfy these two assertions with the different approaches i tried for
the pathological cases...

Anyone has experimented these issues ? Any help is welcome.

Thanks,
--
Pierre Loic Herve

Andrew Willmott

2005-05-30 17:52:16 UTC

Permalink

The short version: WA is not a practical, robust algorithm. This is not
a problem you can solve.

The long version: a friend of mine was working on 2D clipping algorithms
for his Master's thesis. He ran into exactly the same problems you
describe, and spent a lot of time trying to fix them. He eventually got
in touch with one of the paper authors, and established that their code
had the same set of problems. It was a nice algorithm in theory, but it
inevitably ran into robustness issues in practice. He went on to write
his own, which was based off a BSP tree approach, I think. (I might be
wrong about that -- this was all over ten years ago, and I can't find
his thesis online. It was "Object-precision methods for visible surface
determination", J. Williams, University of Auckland..)

Andrew

Post by p***@free.fr
Hello,
I am implementing the Weiler Atherton algorithm as it really matches my needs
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes
http://www.cs.drexel.edu/~david/Classes/CS430/HWs/p214-weiler.pdf
The general cases works really good but it is getting harder with the
pathological cases when some intersections points happens to be the same as
- 2 non parallel edges intersect on one end of one edge
- 2 non parallel edges intersect on one end of one edge which is also an end of
the other edge
- 2 parallel edges intersect on two points which are end of edges
I tried different ways to treat them (insert intersection nodes or not) but
nothing really works for the moment.
"If care is taken in placement of intersections where the subject and clip
polygon contours are identical in the x-y plane, no degenerate polygons will be
produced by the clipping process"
"These two types of intersections will be found to alternate along any given
contour and the number of intersections will always be even"
I cant satisfy these two assertions with the different approaches i tried for
the pathological cases...
Anyone has experimented these issues ? Any help is welcome.
Thanks,
--
Pierre Loic Herve
-------------------------------------------------------
SF.Net email is sponsored by: GoToMeeting - the easiest way to collaborate
online with coworkers and clients while avoiding the high cost of travel and
communications. There is no equipment to buy and you can meet as often as
you want. Try it free.http://ads.osdn.com/?ad_idt02&alloc_id135&op=click
_______________________________________________
GDAlgorithms-list mailing list
https://lists.sourceforge.net/lists/listinfo/gdalgorithms-list
http://sourceforge.net/mailarchive/forum.php?forum_ida88

Bretton Wade

2005-05-31 06:47:08 UTC

Permalink

(Despite my experience in this area, I continue to believe that object
space algorithms like this can be made to work. It's my masochistic
tendencies showing through.)

As a general rule, clipping algorithms on complex geometry are not
robust and accurate with floating point math.

In any large enough data set, you will find at least one unexpected
special case that you must solve. The last two assertions you claim to
have trouble with (types of intersections will alternate and the number
will always be even) are implications of the Jordan curve theorem (Any
continuous simple closed curve in the plane separates the plane into two
disjoint regions: the inside and the outside.) You have to make these
work in order to be sure the clipper is correct.

If this is for offline processing, use arbitrary precision arithmetic to
eliminate a major source of accuracy errors. It will be slow, but the
alternative is manually examining thousands of failure cases.

BSP based clipping is no better than the Weiler-Atherton algorithm in
the long run; it's just easier to code the operations. Unfortunately it
will also give you more polygons, meaning more opportunities for error.

I implemented the Sutherland and Hodgman algorithm a few years back for
offline processing of water polygons in Flight Simulator. Polys with
holes were clipped only to quads, and I had a @$%$ of a time getting
that robust across the entire world.

Here are a few suggestions:

1) Be careful about testing for equality. This comes up in the
collocation and collinear testing a lot. Epsilon testing is a generally
accepted method for comparison of floats [abs(a-b)<epsilon], but breaks
down if epsilon is meaningless in the range of numbers you are working
in. A better test is looking for least significant bits of deviation
(typecast the floats to ints, and compare the result to less than some
small enough integer tolerance).

2) Error check your input religiously. You want polygon chains that are
simple, have no collinear or collocated vertices, have area greater than
0, etc., etc. Do all this conditioning before you feed the clipper, not
while you feed it. Punt on any polygon that isn't well conditioned.

3) Skip the alternative winding order hole polygon approach if you can.
I inserted the hole polygon into the source polygons and updated the
test for simple-ness to allow collinear edges in the polygon. This
removes one major source of complexity from the algorithm.

4) Isolate your intersection routines so that they always order line
segments the same way. This ensures that you get the same intersection
point whether you are considering AB intersects CD, or CD intersects AB
(or DC intersects BA...)

5) To make the Jordan curve theorem hold, I found that I had to sort the
intersections along the clipping plane and fudge it if I got multiple
intersection points that appeared to be collocated. Basically the sort
routine had to enforce the alternating nature of the intersection
points.

If you keep going down this road... Good luck and please report back on
your results.

-----Original Message-----
From: gdalgorithms-list-***@lists.sourceforge.net
[mailto:gdalgorithms-list-***@lists.sourceforge.net] On Behalf Of
Andrew Willmott
Sent: Monday, May 30, 2005 12:51 PM
To: gdalgorithms-***@lists.sourceforge.net
Subject: Re: [Algorithms] Weiler Atherton Implementation (2D polygon
clipping)

The short version: WA is not a practical, robust algorithm. This is not
a problem you can solve.

The long version: a friend of mine was working on 2D clipping algorithms

for his Master's thesis. He ran into exactly the same problems you
describe, and spent a lot of time trying to fix them. He eventually got
in touch with one of the paper authors, and established that their code
had the same set of problems. It was a nice algorithm in theory, but it
inevitably ran into robustness issues in practice. He went on to write
his own, which was based off a BSP tree approach, I think. (I might be
wrong about that -- this was all over ten years ago, and I can't find
his thesis online. It was "Object-precision methods for visible surface
determination", J. Williams, University of Auckland..)

Andrew

Post by p***@free.fr
Hello,
I am implementing the Weiler Atherton algorithm as it really matches my

needs

Post by p***@free.fr
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes
http://www.cs.drexel.edu/~david/Classes/CS430/HWs/p214-weiler.pdf
The general cases works really good but it is getting harder with the
pathological cases when some intersections points happens to be the

same as

Post by p***@free.fr
- 2 non parallel edges intersect on one end of one edge
- 2 non parallel edges intersect on one end of one edge which is also

an end of

Post by p***@free.fr
the other edge
- 2 parallel edges intersect on two points which are end of edges
I tried different ways to treat them (insert intersection nodes or not)

but

Post by p***@free.fr
nothing really works for the moment.
"If care is taken in placement of intersections where the subject and

clip

Post by p***@free.fr
polygon contours are identical in the x-y plane, no degenerate polygons

will be

Post by p***@free.fr
produced by the clipping process"
"These two types of intersections will be found to alternate along any

given

Post by p***@free.fr
contour and the number of intersections will always be even"
I cant satisfy these two assertions with the different approaches i

tried for

Post by p***@free.fr
the pathological cases...
Anyone has experimented these issues ? Any help is welcome.
Thanks,
--
Pierre Loic Herve
-------------------------------------------------------
SF.Net email is sponsored by: GoToMeeting - the easiest way to

collaborate

Post by p***@free.fr
online with coworkers and clients while avoiding the high cost of

travel and

Post by p***@free.fr
communications. There is no equipment to buy and you can meet as often

Post by p***@free.fr
you want. Try it

free.http://ads.osdn.com/?ad_idt02&alloc_id135&op=click

Post by p***@free.fr
_______________________________________________
GDAlgorithms-list mailing list
https://lists.sourceforge.net/lists/listinfo/gdalgorithms-list
http://sourceforge.net/mailarchive/forum.php?forum_ida88

-------------------------------------------------------
This SF.Net email is sponsored by Yahoo.
Introducing Yahoo! Search Developer Network - Create apps using Yahoo!
Search APIs Find out how you can build Yahoo! directly into your own
Applications - visit
http://developer.yahoo.net/?fr=offad-ysdn-ostg-q22005
_______________________________________________
GDAlgorithms-list mailing list
GDAlgorithms-***@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/gdalgorithms-list
Archives:
http://sourceforge.net/mailarchive/forum.php?forum_id=6188

Willem de Boer

2005-05-31 07:35:59 UTC

Permalink

I recently had to implement a polygon scene clipper
algorithm. I used a beamtree to speed up poly clipping,
and a BSP tree to feed the beamtree a front-to-back
ordered list of polygons.

My findings were that you _do not_ want to convert to
triangles during your beamtree and BSP tree constructions,
as this will i) result in an exponential increase of
tree nodes and ii) will entice a lot of rounding-error
creepiness.

Bretton Wade wrote:
"As a general rule, clipping algorithms on complex
geometry are not robust and accurate with floating
point math. "

I can also confirm this.

Cheers,
Willem

Gino van den Bergen

2005-05-31 08:05:30 UTC

Permalink

Post by Bretton Wade
-----Original Message-----
Behalf Of Bretton Wade
Sent: Tuesday, May 31, 2005 10:45 AM
Subject: RE: [Algorithms] Weiler Atherton Implementation (2D
polygon clipping)

<snip>

Post by Bretton Wade
1) Be careful about testing for equality. This comes up in
the collocation and collinear testing a lot. Epsilon testing
is a generally accepted method for comparison of floats
[abs(a-b)<epsilon], but breaks down if epsilon is meaningless
in the range of numbers you are working in. A better test is
looking for least significant bits of deviation (typecast the
floats to ints, and compare the result to less than some
small enough integer tolerance).

I would like to add that you should try to relate your epsilon to the
actual input data. For instance, my (hyper)plane object maintains a
noise value that is used to determine whether a point lies on the plane.
Given a normal, an offset, and a noise value, a point p lies on the
plane if |dot(normal, p) - offset| <= noise. The noise value is not an
absolute epsilon. Rather, it is computed from the vertices that are
known to lie on the plane. Computed normals are usually not perfectly
orthogonal to the support plane of a polygon. Therefore, for these
vertices, the signed distance dot(normal, p) - offset will not be
exactly zero. The noise value is the maximum absolute signed distance of
a vertex that is known to lie on the plane. Furthermore, I add an
epsilon to this noise value that is related to the offset in order to
catch rounding errors. This espilon is taken to be offset *
std::numeric_limits<float>::epsilon() * C, where C is a small positive
constant (e.g. 10). In this way, planes that are further away from the
origin will have a larger noise value.

You should be careful comparing floats by casting them to integers. You
would like the float values 1.0 * 2^exp and 1.11111111111111111111111 *
2^(exp-1) (in base-2 representation) to be classified as equal. However,
a direct cast to integer (through *reinterpret_cast<int*>(&value)) will
show different values after masking away the least significant bits. It
is better to get rid of the N least significant bits by adding a bias
value that is 1.0 * 2^(exp + N) (2^N times the absolute maximum of the
two) to both values and compare their 32-bit float values. For N = 8,
adding the bias results in 1.00000001000000000000000 * 2^(exp + N) for
both values after rounding. (In order to make sure that you're testing
only 32 bits and not the 80-bit internal representation you could cast
them to integers using the pointer cast, or use the float-consistency
compiler option -Op.)

Gino

p***@free.fr

2005-05-31 13:53:40 UTC

Permalink

Thanks for your contributions.

I think i did progress a bit since the initial post as i managed to list two
main pathological cases and i *feel* like i have a solution for the first one
:-).
I still cant find anything for the second case and that is the bad news...

I did a page to explain where i am, you might have a look if you whish :
http://pierloic.free.fr/clip/clip.html

I dont have precision problems for the moment. I simply use fabs(a-b)<epsilon.
That's allright for now and i dont get equality test problems. It may become an
issue later...

I do some strong states check with asserts. Mainly on even count of
intersections and the alternates but that is after the graph is built.

So:

1) If anyone can think of something about those cases, that's *cool* !!!

2) Reading your posts, it sounds like i should give up with this...
Is there any other algorithm for the same specifications:
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes

Thanks
--
Pierre Loic

Willem de Boer

2005-06-01 04:54:25 UTC

Permalink

"Is there any other algorithm for the same specifications:"

A modified polygon beamtree with a BSP tree for front-to-back
ordering of polygons would do. I say modified because you want
to keep the polygons that are fully clipped by the beamtree,
whereas with the original beamtree algorithm you'd throw
those away.

---
Willem H. de Boer
Homepage: http://www.whdeboer.com

Post by Bretton Wade
-----Original Message-----
Sent: Tuesday, May 31, 2005 5:48 PM
Subject: RE: [Algorithms] Weiler Atherton Implementation (2D
polygon clipping)
Thanks for your contributions.
I think i did progress a bit since the initial post as i
managed to list two main pathological cases and i *feel* like
i have a solution for the first one :-).
I still cant find anything for the second case and that is
the bad news...
I did a page to explain where i am, you might have a look if
http://pierloic.free.fr/clip/clip.html
I dont have precision problems for the moment. I simply use
fabs(a-b)<epsilon.
That's allright for now and i dont get equality test
problems. It may become an issue later...
I do some strong states check with asserts. Mainly on even
count of intersections and the alternates but that is after
the graph is built.
1) If anyone can think of something about those cases, that's
*cool* !!!
2) Reading your posts, it sounds like i should give up with this...
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes
Thanks
--
Pierre Loic
-------------------------------------------------------
This SF.Net email is sponsored by Yahoo.
Introducing Yahoo! Search Developer Network - Create apps using Yahoo!
Search APIs Find out how you can build Yahoo! directly into
your own Applications - visit
http://developer.yahoo.net/?fr=fad-ysdn-ostg-q22005
_______________________________________________
GDAlgorithms-list mailing list
https://lists.sourceforge.net/lists/listinfo/gdalgorithms-list
http://sourceforge.net/mailarchive/forum.php?forum_ida88

Sam Martin

2005-06-01 07:19:39 UTC

Permalink

Hi,

I've also yet to see a float implementation of an algorithm that involves
line segment intersection that won't fail given sufficiently challenging
data. There are a few things I can suggest:

Try your best to ignore it if it's offline. The alternatives are all pretty
nasty. Try the brain dead approach of tweaking your epsilons and even adding
noise until you can get it to converge. If it works it'll save you a lot of
effort.

I haven't had much luck with arbitrary arithmetic purely because I haven't
seen a scheme that's efficient enough to use all the time. I may not be the
best authority on this, but I've found it to be something of an all or
nothing approach. If you intend to take the result of an arbitrary
arithmetic-clipped shape and cast it back to regular floats you can still
suffer from topological changes, face collapses, intersecting edges, etc. If
anyone knows any way round this I'd be very interested to hear about it.

I have used the GPC library for clipping and csg operations before
(http://www.cs.man.ac.uk/aig/staff/alan/software/). However, it's another
float implementation and therefore also suffers from numerical precision
problems. The upside is it doesn't fail very often and it's fairly easy to
hack so it doesn't crash when it fails and still returns something half
sensible. That may sound like a tremendous cop-out but it's at least a
usable solution!

However, if you really need something that is rock solid under any
configuration and efficient enough to use at runtime I can suggest you look
up snap rounding: i.e. http://citeseer.csail.mit.edu/guibas95rounding.html.
It's purely a solution for calculating line segment intersections, but it
does so completely robustly and I feel there's some scope in adapting this
work into more complex algorithms. For my purposes I managed to marry snap
rounding with constrained delaunay triangulations (similar in style to:
http://ligwww.epfl.ch/Publications/pdf/Kallmann_and_al_Geometric_Modeling_03.pdf).
The result is entirely robust and not a mile away from a general clipping or
csg implementation, but it wasn't a light undertaking by any means.

Hope that helps!

Regards,
Sam

Message: 5
Date: Tue, 31 May 2005 17:47:56 +0200
Subject: RE: [Algorithms] Weiler Atherton Implementation (2D polygon
clipping)
Thanks for your contributions.
I think i did progress a bit since the initial post as i managed to list =
two
main pathological cases and i *feel* like i have a solution for the first=
one
:-).
I still cant find anything for the second case and that is the bad news..=
.
http://pierloic.free.fr/clip/clip.html
I dont have precision problems for the moment. I simply use fabs(a-b)<eps=
ilon.
That's allright for now and i dont get equality test problems. It may bec=
ome an
issue later...
I do some strong states check with asserts. Mainly on even count of
intersections and the alternates but that is after the graph is built.
1) If anyone can think of something about those cases, that's *cool* !!!
2) Reading your posts, it sounds like i should give up with this...
- concave with holes polygon clipped with concave with holes polygon
- the algorithm outputs 2 lists: inside polygons and outside polygons
- ability to use results as inputs for several passes
Thanks
--
Pierre Loic
--__--__--
_______________________________________________
GDAlgorithms-list mailing list
https://lists.sourceforge.net/lists/listinfo/gdalgorithms-list
End of GDAlgorithms-list Digest

Gino van den Bergen

2005-06-01 07:23:10 UTC

Permalink

I'm using the BSP tree approach, so I'm clipping convex polygons against
(hyper)planes, which is easier than your case. However, in my case it
helped to treat edges (actually half-edges) as half-open. That is, the
edge connecting points p and q is the set { (1 - t) * p + t * q : 0 < t
<= 1 }. A single iteration of the clipper does something like this

if (bottomSign * topSign <= 0 && bottomSign != 0)
{
if (topSign == 0)
{
// Degenerate case: top vertex lies on the plane

if (bottomSign * topSignOfNextVertex < 0)
{
// Add top vertex
}
}
else
{
// Edge intersects at interior. Add intersection
point.
}
}

In this way, there is only one degenerate case, and I always end up with
zero or two intersection points.

NB: This only works for convex polygons. Using finite-precision numbers,
a polygon having more than three vertices has a very slim chance of
being absolutely convex. (Each internal edge has a slight crease.) In
order to deal with this issue, my planes are fattened by a noise factor,
so a sign of zero actually means distance(point, plane) <= plane.noise.