Thread Subject: How to get the intersection point coordinate between polygon and segment

How to get the intersection point coordinate between polygon and segment;
for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
I want to return c=[0.5 0.5;0.4 0.4]
example2:
a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
I want to return c=[0.4 0;1 0]
example 3
a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
I want to return c=[0.4 0;0.5 0];
Thank you for your attention.Is there any mex implementation.and stable.I also write one but I think because of the eps I can't implement it well.
thank you for another time.

"zedong
" <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> How to get the intersection point coordinate between polygon and segment;
> for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
> a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> I want to return c=[0.5 0.5;0.4 0.4]
> example2:
> a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
> I want to return c=[0.4 0;1 0]
> example 3
> a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
> I want to return c=[0.4 0;0.5 0];
> Thank you for your attention.Is there any mex implementation.and stable.I also write one but I think because of the eps I can't implement it well.
> thank you for another time.

Did you look on the file exchange? Try it.
Here are some good key words to use:

curve intersections

John

"zedong
" <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> How to get the intersection point coordinate between polygon and segment;
> for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
> a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> I want to return c=[0.5 0.5;0.4 0.4]
> example2:
> a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
> I want to return c=[0.4 0;1 0]
> example 3
> a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
> I want to return c=[0.4 0;0.5 0];

Questions

(a) Can we assume the polygon is convex? If not, the problem is ill-defined.

(b) Do the rows of a give the polygon vertices in clockwise/counter-clockwise order, or is the order arbitrary? If the latter, then a sorting operation will have to be performed first.

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gi0gg5$878$1@fred.mathworks.com>...
> "zedong
> " <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> > How to get the intersection point coordinate between polygon and segment;
> > for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
> > a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> > I want to return c=[0.5 0.5;0.4 0.4]
> > example2:
> > a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
> > I want to return c=[0.4 0;1 0]
> > example 3
> > a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
> > I want to return c=[0.4 0;0.5 0];
>
> Questions
>
> (a) Can we assume the polygon is convex? If not, the problem is ill-defined.
>
> (b) Do the rows of a give the polygon vertices in clockwise/counter-clockwise order, or is the order arbitrary? If the latter, then a sorting operation will have to be performed first.
Thank you for your attention .we can assume that the polygon is convex.as a matter of fact it is a triangle.But I have many triangle.and many segment.so I must implement it use mex.

It is a convex.as a matter of fact it's just a triangle.and I can know that it's anti-clockwise.
But I have many segment and many triangle.I must use it in mex file.Thank you all.
"zedong
" <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> How to get the intersection point coordinate between polygon and segment;
> for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
> a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> I want to return c=[0.5 0.5;0.4 0.4]
> example2:
> a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
> I want to return c=[0.4 0;1 0]
> example 3
> a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
> I want to return c=[0.4 0;0.5 0];
> Thank you for your attention.Is there any mex implementation.and stable.I also write one but I think because of the eps I can't implement it well.
> thank you for another time.

> But I have many segment and many triangle.I must use it in mex file.Thank you all.

Why does it have to be a mex file? Certainly not just to avoid for loops.

Once you have an mfile that does it for one segment/triangle, you can do it for many using cellfun.

"zedong
" <zdongwu@gmail.com> wrote in message <gi0inj$t2m$1@fred.mathworks.com>...
> It is a convex.as a matter of fact it's just a triangle.and I can know that it's anti-clockwise.
> But I have many segment and many triangle.I must use it in mex file.Thank you all.

But a well written m-file can be as fast as a mex
function, especially if the mex function is not
carefully optimized.

John

"zedong
" <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> How to get the intersection point coordinate between polygon and segment;
> for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):
> a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> I want to return c=[0.5 0.5;0.4 0.4]
> example2:
> a=[0 0;1 0;0 1]; b=[0.4 0;2 0];
> I want to return c=[0.4 0;1 0]
> example 3
> a=[0 0;1 0;0 1]; b=[0.4 0;0.5 0];
> I want to return c=[0.4 0;0.5 0];
> Thank you for your attention.Is there any mex implementation.and stable.I also write one but I think because of the eps I can't implement it well.
> thank you for another time.

I worked on this recently. If you want the equations for computing intersections between line segments and triangles, you can first google "ray triangle intersection" and check out some of the links. Alternatively, you can peek into one of my FEX submissions:

http://www.mathworks.com/matlabcentral/fileexchange/20637

This submission is specifically designed to work with RENDERED surfaces, but if you have triangle vertices stored in variables then you can look at how my code handles multiple ray-triangle intersections to figure out how to do it for yourself in an efficient way (I tried to use matrix operations over for loops wherever possible).

hth,
Ken

"zedong
" <zdongwu@gmail.com> wrote in message <gi0epq$lk3$1@fred.mathworks.com>...
> How to get the intersection point coordinate between polygon and segment;
> for example(a is a triangle a(i,:) is the ith coordinate of the triangle.b is an segment,b(i,:) is the side point coordinate):

I believe the function below does the job. It employs vert2con from the FEX

I also tested its speed on 10000 successive triangles using both a for loop and cellfun()

a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
A=cell(1,10000);B=A; A(:)={a}; B(:)={b};

>> tic, C=cellfun(@cutpoints,A,B,'UniformOutput',false); toc
Elapsed time is 11.256645 seconds.

>> tic, for ii=1:10000, c=cutpoints(a,b); end,toc
Elapsed time is 11.231848 seconds.

Not much difference, so I guess you'll have to be the judge of whether it's fast enough.


function c=cutpoints(a,b)
%find intersection points of line segment with polytope
%
% c=cutpoints(a,b)
%
%Rows of a define polytope vertices
%Rows of b are line segment end points
%Rows of c are intersection points


c=nan(2);

 xx=b(1,:).';
 yy=b(2,:).';
 dd=yy-xx;
  
 [AA,bb]=vert2con(a);
 
 CC=AA*dd; qq=bb-AA*xx;

 db=(CC==0);
 if any(qq(db)<0) %ray does not intersect polytope
    return
 end
 
 
 %%%Max lower bound
 lb=zeros(size(CC));
 
 ii=CC<0; lb(ii)=qq(ii)./CC(ii);
 
 tmin=max(max(lb),0);
 
  %%%Min upper bound
  
   ub=ones(size(CC));
   
   ii=CC>0; ub(ii)=qq(ii)./CC(ii);
   
   tmax=min(min(ub),1);
   
   c(1,:)=(xx+tmin*dd).';
   c(2,:)=(xx+tmax*dd).';

> I also tested its speed on 10000 successive triangles using both a for loop and cellfun()
>
> a=[0 0;1 0;0 1]; b=[0.4 0.4;1 1];
> A=cell(1,10000);B=A; A(:)={a}; B(:)={b};
>
> >> tic, C=cellfun(@cutpoints,A,B,'UniformOutput',false); toc
> Elapsed time is 11.256645 seconds.
>
> >> tic, for ii=1:10000, c=cutpoints(a,b); end,toc
> Elapsed time is 11.231848 seconds.
>
> Not much difference, so I guess you'll have to be the judge of whether it's fast enough.

%%%%%%%%%%%%%%%%


I found that the code becomes much faster, when I replace vert2con() with vert2con_special() below. The latter assumes a 2D convex polygon and with vertices sequentially ordered.

>> tic, C=cellfun(@cutpoints,A,B,'UniformOutput',false); toc
Elapsed time is 1.092368 seconds.
>> tic, for ii=1:10000, c=cutpoints(a,b); end,toc
Elapsed time is 1.107186 seconds.


For-loops and cellfun() are still about the same though...



function [A,b]=vert2con_special(a)

centroid=mean(a).';
R=[0 1; -1 0];

A=diff([a;a(1,:)])*R;

b=sum(A.*a,2);

ii=(A*centroid>=b);

b(ii)=-b(ii);
A(ii,:)=-A(ii,:);

I had some bugs in my previously posted code. Here it is again corrected (hopefully).


function c=cutpoints(a,b)
%find intersection points of line segment with polytope
%
% c=cutpoints(a,b)
%
%Rows of a define polytope vertices
%Rows of b are line segment end points
%Rows of c are intersection points


c=nan(2);

 xx=b(1,:).';
 yy=b(2,:).';
 dd=yy-xx;
  
 %[AA,bb]=vert2con(a);
 [AA,bb]=vert2con_special(a);
 
 
 CC=AA*dd; qq=bb-AA*xx;

 db=(CC==0);
 if any(qq(db)<0) %ray does not intersect polytope
    return
 end
 
 
 %%%Max lower bound
 
 ii=CC<0; lb=qq(ii)./CC(ii);
 
 tmin=max(lb);
 
 %%%Min upper bound
 
   
 ii=~ii; ub=qq(ii)./CC(ii);
   
 tmax=min(ub);
   
if tmin>1 | tmax<0, return, end %no intersection
    
tmin=max(tmin,0);
tmax=min(tmax,1);

   c(1,:)=(xx+tmin*dd).';
   c(2,:)=(xx+tmax*dd).';


%%%%%%%%%
function [A,b]=vert2con_special(a)

centroid=mean(a).';
R=[0 1; -1 0];

A=diff([a;a(1,:)])*R;

b=sum(A.*a,2);

ii=(A*centroid>=b);

b(ii)=-b(ii);
A(ii,:)=-A(ii,:);

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gi0rih$ljc$1@fred.mathworks.com>...
> I had some bugs in my previously posted code. Here it is again corrected (hopefully).
>
>
> function c=cutpoints(a,b)
> %find intersection points of line segment with polytope
> %
> % c=cutpoints(a,b)
> %
> %Rows of a define polytope vertices
> %Rows of b are line segment end points
> %Rows of c are intersection points
>
>
> c=nan(2);
>
> xx=b(1,:).';
> yy=b(2,:).';
> dd=yy-xx;
>
> %[AA,bb]=vert2con(a);
> [AA,bb]=vert2con_special(a);
>
>
> CC=AA*dd; qq=bb-AA*xx;
>
> db=(CC==0);
> if any(qq(db)<0) %ray does not intersect polytope
> return
> end
>
>
> %%%Max lower bound
>
> ii=CC<0; lb=qq(ii)./CC(ii);
>
> tmin=max(lb);
>
> %%%Min upper bound
>
>
> ii=~ii; ub=qq(ii)./CC(ii);
>
> tmax=min(ub);
>
> if tmin>1 | tmax<0, return, end %no intersection
>
> tmin=max(tmin,0);
> tmax=min(tmax,1);
>
> c(1,:)=(xx+tmin*dd).';
> c(2,:)=(xx+tmax*dd).';
>
>
> %%%%%%%%%
> function [A,b]=vert2con_special(a)
>
> centroid=mean(a).';
> R=[0 1; -1 0];
>
> A=diff([a;a(1,:)])*R;
>
> b=sum(A.*a,2);
>
> ii=(A*centroid>=b);
>
> b(ii)=-b(ii);
> A(ii,:)=-A(ii,:);


Thank you very much.Thank you for your kind help.You have done it very well.At many cases.It's right.There are still two problem:
(1) I have tested that if a=[0 0;1 0;0 1] b=[-1 0;1 0];
the right answer must be c=[0 0;1 0] of course.But your code return c=[-1 0;1 0].It maybe a small bug.
(2)I will tell you that why I say that we should implement it by mex file.If I write it in mex file.I have tested that 1million by 1 thousand it only need less than 1 second.

> Thank you very much.Thank you for your kind help.You have done it very well.At many cases.It's right.There are still two problem:
> (1) I have tested that if a=[0 0;1 0;0 1] b=[-1 0;1 0];
> the right answer must be c=[0 0;1 0] of course.But your code return c=[-1 0;1 0].It maybe a small bug.


I am not seeing this error. I get the correct answer c=[0 0;1 0]

I've posted my code again below, in case I accidentally gave you a buggy version before.

> (2)I will tell you that why I say that we should implement it by mex file.If I write it in mex file.I have tested that 1million by 1 thousand it only need less than 1 second.


Sorry, I didn't understand that. What does "1million by 1 thousand" mean?

You mean you can process 10^6*10^3=10^9 sets of a,b data in less than a second?

And if you already have a mex file that does so, why do you need help from us?



function c=cutpoints(a,b)
%find intersection points of line segment with polytope
%
% c=cutpoints(a,b)
%
%Rows of a define polytope vertices
%Rows of b are line segment end points
%Rows of c are intersection points


c=nan(2);

 xx=b(1,:).';
 yy=b(2,:).';
 dd=yy-xx;
  
 %[AA,bb]=vert2con(a);
 [AA,bb]=vert2con_special(a);
 
 
 CC=AA*dd; qq=bb-AA*xx;

 db=(CC==0);
 if any(qq(db)<0) %ray does not intersect polytope
    return
 end
 
 
 %%%Max lower bound
 
 ii=CC<0; lb=qq(ii)./CC(ii);
 
 tmin=max(lb);
 
 %%%Min upper bound
 
   
 ii=CC>0; ub=qq(ii)./CC(ii);
   
 tmax=min(ub);
   
if tmin>1 | tmax<0 | tmax<tmin, return, end %no intersection
    
tmin=max(tmin,0);
tmax=min(tmax,1);

   c(1,:)=(xx+tmin*dd).';
   c(2,:)=(xx+tmax*dd).';

"Matt" <mjacobson.removethis@xorantech.com> wrote in message <gi38m1$h17$1@fred.mathworks.com>...
> > Thank you very much.Thank you for your kind help.You have done it very well.At many cases.It's right.There are still two problem:
> > (1) I have tested that if a=[0 0;1 0;0 1] b=[-1 0;1 0];
> > the right answer must be c=[0 0;1 0] of course.But your code return c=[-1 0;1 0].It maybe a small bug.
>
>
> I am not seeing this error. I get the correct answer c=[0 0;1 0]
>
> I've posted my code again below, in case I accidentally gave you a buggy version before.
>
> > (2)I will tell you that why I say that we should implement it by mex file.If I write it in mex file.I have tested that 1million by 1 thousand it only need less than 1 second.
>
>
> Sorry, I didn't understand that. What does "1million by 1 thousand" mean?
>
> You mean you can process 10^6*10^3=10^9 sets of a,b data in less than a second?
>
> And if you already have a mex file that does so, why do you need help from us?
>
>
>
> function c=cutpoints(a,b)
> %find intersection points of line segment with polytope
> %
> % c=cutpoints(a,b)
> %
> %Rows of a define polytope vertices
> %Rows of b are line segment end points
> %Rows of c are intersection points
>
>
> c=nan(2);
>
> xx=b(1,:).';
> yy=b(2,:).';
> dd=yy-xx;
>
> %[AA,bb]=vert2con(a);
> [AA,bb]=vert2con_special(a);
>
>
> CC=AA*dd; qq=bb-AA*xx;
>
> db=(CC==0);
> if any(qq(db)<0) %ray does not intersect polytope
> return
> end
>
>
> %%%Max lower bound
>
> ii=CC<0; lb=qq(ii)./CC(ii);
>
> tmin=max(lb);
>
> %%%Min upper bound
>
>
> ii=CC>0; ub=qq(ii)./CC(ii);
>
> tmax=min(ub);
>
> if tmin>1 | tmax<0 | tmax<tmin, return, end %no intersection
>
> tmin=max(tmin,0);
> tmax=min(tmax,1);
>
> c(1,:)=(xx+tmin*dd).';
> c(2,:)=(xx+tmax*dd).';
 Thank you very much.This time it's really right.Could you explain your compute method.The computation idea to me. I can change it into mex function.Thank you very much.

> Thank you very much.This time it's really right.Could you explain your compute method.The computation idea to me. I can change it into mex function.Thank you very much.

Sure. Any convex polygon/polyhedron, and in particular a triangle, can be expressed as a set of linear inequalities. In matrix-vector form, point z is in the polygon if it satisfies

A*z<=b

In my code, the subfunction vert2con_special takes the given triangle vertices and finds the corresponding [A,b].

Additionally, the parametric equation for a line segment between x and y is

z(t)=x+t*(y-x), 0<=t<=1.

Substituting z(t) into the polygon inequalities gives
  
(A*(y-x))*t<=(b-A*x)

You must find the minimum and maximum scalar t satisfying all of these inequalities. Once you have tmin and tmax, the intersection points are given by z(tmin) and z(tmax)

NOTE: The requirement 0<=t<=1 ensures that z(t) will lie on the line segment between x and y. You must therefore threshold tmin and tmax to this interval appropriately.