http://www.mathworks.com/matlabcentral/newsreader/view_thread/328722
MATLAB Central Newsreader  Polybool gives subtraction with wrong vertices
Feed for thread: Polybool gives subtraction with wrong vertices
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 01 May 2013 22:57:09 +0000
Polybool gives subtraction with wrong vertices
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328722#903570
Kenan
Hello all,<br>
<br>
I am using the function polybool from the Mapping Toolbox and my version of Matlab is R2013a.<br>
<br>
I have problems when I perform a simple subtraction of two polygonal regions (a rectangle described by [x1 y1] and a triangle described by [x0 y0]). See attached code and figures.<br>
<br>
px = .7;<br>
py = .75;<br>
lx = (5*2^.56)/20;<br>
ly = (10*3^.59)/60;<br>
<br>
[x1,y1] = poly2cw([px; px+lx; px+lx; px; px],[py; py; py+ly; py+ly; py]); % rectangle<br>
[x0,y0] = poly2cw([15;15;14;15]/19,[15;16;16;15]/19); % triangle<br>
[x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection<br>
<br>
[x3,y3] = polybool('subtraction',x0,y0,x2,y2); % soustraction<br>
<br>
figure; plot(x0,y0,'b',x1,y1,'k:',x3,y3,'r+'); hold on<br>
for i=1:length(x3)<br>
text(x3(i)*1.003,y3(i)*1.003,num2str(i))<br>
end<br>
<br>
As you can see, there is an extra vertex that should not be there (vertex 23). So instead of being 12,23,34,45,56; the subtraction should be 13,34,45,56. (6 being equal to 1).<br>
<br>
<br>
This wouldn't be such a bad issue if it was all, but it causes the rest of the algorithm to behave abnormally (after performing a constrained Delaunay triangulation, the function isInterior gives an incorrect result)<br>
<br>
First I check if the subtraction [x3 y3] is a connected polygon with ~isShapeMultipart(x3,y3): it is.<br>
<br>
If I perform a non constrained Delaunay triangulation, the result contains triangles outside of [x3 y3] (in the general case) and there is no direct way of finding which ones. Hence it is easier to perform a constrained Delaunay triangulation. If I understand the word constrain correctly, Matlab performs the triangulation making sure the vertices in the constrain are vertices of the triangulation. However, it does not mean that these vertices will be the only vertices of the triangulation. I need someone to confirm this statement, or explain to me what is a constrain exactly. So I construct the constraint matrix for the Delaunay triangulation and perform the triangulation:<br>
<br>
C3 = zeros(length(x3)1,2);<br>
for j=1:2<br>
for i2=1:length(x3)1<br>
C3(i2,j) = i2+j1;<br>
end<br>
end<br>
C3(end,end)=1;<br>
DT31 = delaunayTriangulation(x3,y3,C3);<br>
<br>
Finally, the statement inside31 = isInterior(DT31) gives 0 0 0 meaning that all 3 triangles in DT31 are outside of the constraint C3.<br>
<br>
The code for constructing the constrain matrix C3 could be adapted to the weird result of polybool but I would lose robustness in my code. Plus, I would rather correct the weird result the earliest possible (ie just after [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % soustraction) rather than making corrections afterwards.<br>
<br>
Thank you for your precious help and/or remarks.<br>
<br>
See here for illustrations: <a href="http://stackoverflow.com/questions/16326810/polyboolgivessubtractionwithwrongvertices">http://stackoverflow.com/questions/16326810/polyboolgivessubtractionwithwrongvertices</a>

Tue, 07 May 2013 18:57:09 +0000
Re: Polybool gives subtraction with wrong vertices
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328722#903900
Kenan
NOTE: I mistook 'vertex' with 'edge'.<br>
<br>
I received an answer from the technical support: "the algorithms are kind of sensitive to these rounding errors. I do not know whether this can be improved, but I will inform development, maybe they can make polybool more robust against these problems."

Tue, 07 May 2013 22:13:10 +0000
Re: Polybool gives subtraction with wrong vertices
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328722#903911
Amy Haskins
Hi Kenan,<br>
<br>
In regards to your question, it looks like you're running into some rounding errors here. If you subtract the rectangle from the large triangle directly, you will get the result you would expect:<br>
<br>
[x3,y3] = polybool('subtraction',x0,y0,x1,y1); % subtraction<br>
<br>
When you do this as two steps, there is a very small amount of rounding error in the calculation which causes the extra line in your results. If we shift the bottom vertex down even slightly, the correct answer is given:<br>
<br>
y2([1,4]) = y2([1,4])  .000001;<br>
[x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction<br>
<br>
In a world of infinite percision, the twostep approach would work just fine. However, since we are dealing with floating point, there is about a 50/50 chance that each edges of the smaller triangle will be ever so slightly inside the larger one. <br>
<br>
This is just a natural limitation of doing computation geometry with floating point numbers, the results from polybool are not surprising.<br>
<br>
Amy<br>
<br>
"Kenan" wrote in message <kls6k5$osr$1@newscl01ah.mathworks.com>...<br>
> <br>
> I am using the function polybool from the Mapping Toolbox and my version of Matlab is R2013a.<br>
> <br>
> I have problems when I perform a simple subtraction of two polygonal regions (a rectangle described by [x1 y1] and a triangle described by [x0 y0]). See attached code and figures.<br>
> <br>
> px = .7;<br>
> py = .75;<br>
> lx = (5*2^.56)/20;<br>
> ly = (10*3^.59)/60;<br>
> <br>
> [x1,y1] = poly2cw([px; px+lx; px+lx; px; px],[py; py; py+ly; py+ly; py]); % rectangle<br>
> [x0,y0] = poly2cw([15;15;14;15]/19,[15;16;16;15]/19); % triangle<br>
> [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection<br>
> <br>
> [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % soustraction<br>
> <br>
> figure; plot(x0,y0,'b',x1,y1,'k:',x3,y3,'r+'); hold on<br>
> for i=1:length(x3)<br>
> text(x3(i)*1.003,y3(i)*1.003,num2str(i))<br>
> end<br>
> <br>
> As you can see, there is an extra vertex that should not be there (vertex 23). So instead of being 12,23,34,45,56; the subtraction should be 13,34,45,56. (6 being equal to 1).<br>
> <br>
>

Wed, 08 May 2013 18:54:08 +0000
Re: Polybool gives subtraction with wrong vertices
http://www.mathworks.com/matlabcentral/newsreader/view_thread/328722#904001
Kenan
Thank you Amy. Yes, I had noticed I was obtaining the correct result when using "triangle MINUS rectangle" instead of "triangle MINUS (rectangle INTER triangle)".<br>
However, I was curious as to why there was a difference between the two methods (that are on paper the same). I understand it comes from rounding errors. However, correct me if I'm wrong, but there is no way of knowing *for sure* that the direct method will never fail either, right?<br>
<br>
Talking about computational geometry and rounding errors, do you think I will be luckier if I use functions from the PDE Toolbox (initmesh instead of delaunayTriangulation for instance), or the rounding errors would be the same?<br>
<br>
Thank you again.<br>
<br>
Kenan<br>
<br>
"Amy Haskins" <amy.haskins.nospam@mathworks.com> wrote in message <kmbu9m$q70$1@newscl01ah.mathworks.com>...<br>
> Hi Kenan,<br>
> <br>
> In regards to your question, it looks like you're running into some rounding errors here. If you subtract the rectangle from the large triangle directly, you will get the result you would expect:<br>
> <br>
> [x3,y3] = polybool('subtraction',x0,y0,x1,y1); % subtraction<br>
> <br>
> When you do this as two steps, there is a very small amount of rounding error in the calculation which causes the extra line in your results. If we shift the bottom vertex down even slightly, the correct answer is given:<br>
> <br>
> y2([1,4]) = y2([1,4])  .000001;<br>
> [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction<br>
> <br>
> In a world of infinite percision, the twostep approach would work just fine. However, since we are dealing with floating point, there is about a 50/50 chance that each edges of the smaller triangle will be ever so slightly inside the larger one. <br>
> <br>
> This is just a natural limitation of doing computation geometry with floating point numbers, the results from polybool are not surprising.<br>
> <br>
> Amy<br>
> <br>
> "Kenan" wrote in message <kls6k5$osr$1@newscl01ah.mathworks.com>...<br>
> > <br>
> > I am using the function polybool from the Mapping Toolbox and my version of Matlab is R2013a.<br>
> > <br>
> > I have problems when I perform a simple subtraction of two polygonal regions (a rectangle described by [x1 y1] and a triangle described by [x0 y0]). See attached code and figures.<br>
> > <br>
> > px = .7;<br>
> > py = .75;<br>
> > lx = (5*2^.56)/20;<br>
> > ly = (10*3^.59)/60;<br>
> > <br>
> > [x1,y1] = poly2cw([px; px+lx; px+lx; px; px],[py; py; py+ly; py+ly; py]); % rectangle<br>
> > [x0,y0] = poly2cw([15;15;14;15]/19,[15;16;16;15]/19); % triangle<br>
> > [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection<br>
> > <br>
> > [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % soustraction<br>
> > <br>
> > figure; plot(x0,y0,'b',x1,y1,'k:',x3,y3,'r+'); hold on<br>
> > for i=1:length(x3)<br>
> > text(x3(i)*1.003,y3(i)*1.003,num2str(i))<br>
> > end<br>
> > <br>
> > As you can see, there is an extra vertex that should not be there (vertex 23). So instead of being 12,23,34,45,56; the subtraction should be 13,34,45,56. (6 being equal to 1).<br>
> > <br>
> >