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)".
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?
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?
Thank you again.
Kenan
"Amy Haskins" <amy.haskins.nospam@mathworks.com> wrote in message <kmbu9m$q70$1@newscl01ah.mathworks.com>...
> Hi Kenan,
>
> 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:
>
> [x3,y3] = polybool('subtraction',x0,y0,x1,y1); % subtraction
>
> 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:
>
> y2([1,4]) = y2([1,4])  .000001;
> [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % subtraction
>
> 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.
>
> This is just a natural limitation of doing computation geometry with floating point numbers, the results from polybool are not surprising.
>
> Amy
>
> "Kenan" wrote in message <kls6k5$osr$1@newscl01ah.mathworks.com>...
> >
> > I am using the function polybool from the Mapping Toolbox and my version of Matlab is R2013a.
> >
> > 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.
> >
> > px = .7;
> > py = .75;
> > lx = (5*2^.56)/20;
> > ly = (10*3^.59)/60;
> >
> > [x1,y1] = poly2cw([px; px+lx; px+lx; px; px],[py; py; py+ly; py+ly; py]); % rectangle
> > [x0,y0] = poly2cw([15;15;14;15]/19,[15;16;16;15]/19); % triangle
> > [x2,y2] = polybool('intersection',x1,y1,x0,y0); % intersection
> >
> > [x3,y3] = polybool('subtraction',x0,y0,x2,y2); % soustraction
> >
> > figure; plot(x0,y0,'b',x1,y1,'k:',x3,y3,'r+'); hold on
> > for i=1:length(x3)
> > text(x3(i)*1.003,y3(i)*1.003,num2str(i))
> > end
> >
> > 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).
> >
> >
