Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Speedup: Fast polar intersection

Subject: Speedup: Fast polar intersection

From: Johannes Korsawe

Date: 11 Apr, 2011 13:37:06

Message: 1 of 9

Dear community,

i have the following task.

Let xy be a closed curve in 2d around [0,0] with size(xy)=[n1 2].
Let a be a vector of angles in between 0 and 2*pi and r = [cos(a),sin(a)] a vector of polar directions with size(r) = [n2 2];

I am looking for a fast (!) algorithm to compute the norm of the outmost intersections (there could be multiple) of xy with the rays r (consider r beeing prolongated infinitely from [0,0]). Using different words, i want to compute approximately the polar coordinates of the curve xy according to the angles.

I have already written a vectorized function:

nw=size(r,1);tol = 1e-6;
b=diff(xy,1)';a=xy(1:end-1,:)';
sraw=repmat(b(2,:).*a(1,:)-b(1,:).*a(2,:),nw,1);
b2=[b(2,:);-b(1,:)];a2=[-a(2,:);a(1,:)];
det=r*b2;
% due to numerical stability set ...
det(det==0)=tol;
s=sraw./det;
traw=r*a2./det;
t=traw./det;
ind=-tol<=t & t<=1+tol & s>-tol;
s(~ind)=0;
result=max(s,[],2);

Is there any idea to make this faster???

Best regards,
Johannes

Subject: Speedup: Fast polar intersection

From: Roger Stafford

Date: 11 Apr, 2011 19:23:05

Message: 2 of 9

"Johannes Korsawe" wrote in message <inv061$f7p$1@fred.mathworks.com>...
> Dear community,
>
> i have the following task.
>
> Let xy be a closed curve in 2d around [0,0] with size(xy)=[n1 2].
> Let a be a vector of angles in between 0 and 2*pi and r = [cos(a),sin(a)] a vector of polar directions with size(r) = [n2 2];
>
> I am looking for a fast (!) algorithm to compute the norm of the outmost intersections (there could be multiple) of xy with the rays r (consider r beeing prolongated infinitely from [0,0]). Using different words, i want to compute approximately the polar coordinates of the curve xy according to the angles.
>
> I have already written a vectorized function:
>
> nw=size(r,1);tol = 1e-6;
> b=diff(xy,1)';a=xy(1:end-1,:)';
> sraw=repmat(b(2,:).*a(1,:)-b(1,:).*a(2,:),nw,1);
> b2=[b(2,:);-b(1,:)];a2=[-a(2,:);a(1,:)];
> det=r*b2;
> % due to numerical stability set ...
> det(det==0)=tol;
> s=sraw./det;
> traw=r*a2./det;
> t=traw./det;
> ind=-tol<=t & t<=1+tol & s>-tol;
> s(~ind)=0;
> result=max(s,[],2);
>
> Is there any idea to make this faster???
>
> Best regards,
> Johannes
- - - - - - - - - -
  I notice that your computation combines every possible r ray with every possible segment in your closed curve xy. That would make the algorithm of order O(n1*n2), which would be very large if n1 and n2 were both large.

  I don't have time - I have to do blasted income taxes - to work out the details, but here is the outline of an idea to make it of order O(n1+n2*log(n2)).

  Steps:

1. Do 'atan2' on each point in xy and then apply 'unwrap' so that successive angles never exceed pi in absolute value. If your curve is actually closed, the difference between the first and last angle should be plus or minus 2*pi (or else multiples of that if the curve crosses over itself,) though in between, the angles may range over a larger interval.

2. Do 'atan2' to find each angle in r and then sort these angles in ascending order. That's your O(n2*log(n2)) part. At this point they will range from -pi to +pi.

3. Add and/or subtract multiples of 2*pi of these r angles sufficient to cover the full range of angles that appear in the "unwrapped" xy angles. This is to enable the comparisons that are to be made in step 4.

4. Run through each successive xy segment angle interval, picking up angles from the step 3 r-list that lie in between and compute the polar coordinates of their intersections with the segment. This will probably require some kind of while-loop which is able to go up or down angle-wise as you traverse the r-angles in the interval.

5. Finally run through the r-angles, (taken modulo 2*pi,) and in cases where more than one segment is crossed, find the maximum norm case. This should leave you with one set of polar coordinates for each original r-angle.

  (If you have trouble carrying out this ambitious scheme I might be able to give further help after my taxes are done, or, even better, someone else might have come up with a superior solution.)

Roger Stafford

Subject: Speedup: Fast polar intersection

From: Johannes Korsawe

Date: 13 Apr, 2011 07:23:05

Message: 3 of 9

Roger,

thank you for your remarks.

Following the thoughts of your posts, i came to the idea to sort the angles of the xy curve and insert the rays' angles into this sort via hist or something.

I will try and see. Any more hints are well appreciated. (<-- is this correct English? or is it "much appreciated"?)

Thanks,
Johannes

Subject: Speedup: Fast polar intersection

From: Johannes Korsawe

Date: 13 Apr, 2011 08:34:04

Message: 4 of 9

And here it is:
Much faster than the version before. But it fails, if there are mutiple cuts of the rays r with the curve xy. Luckily, this case happens seldom in my applications, so i can take care of that with the slower variant then and still have a significant speedup.

But how to take care of this inside this idea of histc???

nk = size(xy,1);nw = size(r,1);
scale=zeros(nw,1);
wr=unwrap(atan2(r(:,2),r(:,1)));
wk=unwrap(atan2(kurve(:,2),kurve(:,1)));
if isempty(find(sign(diff(wk))<0)),
   wk=[wk-2*pi;wk;wk+2*pi];
   [n,bin]=histc(wr,wk);
   bin=rem(bin,nk);kurve2=[kurve;kurve(2,:)];
   for i=1:nw,
       k = kurve2(bin(i),:);
       dk=kurve2(bin(i)+1,:)-k;
       scale(i)=(dk(2)*k(1)-dk(1)*k(2))/(r(i,1)*dk(2)-r(i,2)*dk(1));
   end
else, % multiple cuts
          % take care of that seperately
end
 

Best regards,
Johannes

Subject: Speedup: Fast polar intersection

From: Roger Stafford

Date: 15 Apr, 2011 06:24:04

Message: 5 of 9

"Johannes Korsawe" wrote in message <inv061$f7p$1@fred.mathworks.com>...
> Dear community,
>
> i have the following task.
>
> Let xy be a closed curve in 2d around [0,0] with size(xy)=[n1 2].
> Let a be a vector of angles in between 0 and 2*pi and r = [cos(a),sin(a)] a vector of polar directions with size(r) = [n2 2];
>
> I am looking for a fast (!) algorithm to compute the norm of the outmost intersections (there could be multiple) of xy with the rays r (consider r beeing prolongated infinitely from [0,0]). Using different words, i want to compute approximately the polar coordinates of the curve xy according to the angles.
>
> I have already written a vectorized function:
>
> nw=size(r,1);tol = 1e-6;
> b=diff(xy,1)';a=xy(1:end-1,:)';
> sraw=repmat(b(2,:).*a(1,:)-b(1,:).*a(2,:),nw,1);
> b2=[b(2,:);-b(1,:)];a2=[-a(2,:);a(1,:)];
> det=r*b2;
> % due to numerical stability set ...
> det(det==0)=tol;
> s=sraw./det;
> traw=r*a2./det;
> t=traw./det;
> ind=-tol<=t & t<=1+tol & s>-tol;
> s(~ind)=0;
> result=max(s,[],2);
>
> Is there any idea to make this faster???
>
> Best regards,
> Johannes
- - - - - - - -
  I finally got back to your problem, Johannes. I realize now that there is no need to do an 'unwrap'. Each segment of your xy polygon can be dealt with independently of the others. However, when two successive vertices have atan2 angles that differ by more than pi, it is necessary to make an appropriate adjustment that allows the equivalent of the lower of the two angles being increased by 2*pi. This is accomplished indirectly by the manipulations you see in the code below with the start and end indices, c1 and c2, so that s(c1:c2) can address all the angles in the sorted 'as' angle vector that intersect with the line segment between the vertices.

  The assumption here is that 'a' is a column vector of angles that lie between -pi and +pi, and that in the two-column xy array the first and last points coincide. The vector r is the result of the computation and r.*cos(a),r.*sin(a) give you the coordinates of the most distant segment intersection along each 'a' ray.

  As you are probably aware, when the two vertices of a line segment lie very nearly on one of the 'a' rays, there is a reduction in accuracy. If they lie exactly on the ray, a NaN can be produced from the indeterminacy of the quotient in the calculation of r when both the numerator and denominator become zero. That is inherent in the problem.

  I hope this helps with your speed problem.

m = size(xy,1);
n = size(a,1);
[as,p] = sort(a);
ca = cos(as); sa = sin(as);
x = xy(:,1); y = xy(:,2);
dx = diff(x); dy = diff(y);
axy = atan2(y,x);
[~,b] = histc(axy,[-inf;as;inf]);
s = repmat((1:n).',2,1);
r = zeros(n,1);
for k = 1:m-1 % Do this once for each line segment
 b1 = b(k); b2 = b(k+1);
 if abs(axy(k)-axy(k+1)) < pi
  c1 = min(b1,b2);
  c2 = b1+b2-c1-1;
 else % In this case we have to play tricks with c1 and c2
  c1 = max(b1,b2);
  c2 = b1+b2-c1+n-1; % Adding n is the equivalent of adding 2*pi
 end
 if c1<=c2 % Don't bother with empty sets
  t = s(c1:c2);
  u = (x(k)*dy(k)-y(k)*dx(k))./(ca(t)*dy(k)-sa(t)*dx(k));
  r(t) = max(r(t),u);
 end
end
r(p) = r;

Roger Stafford

Subject: Speedup: Fast polar intersection

From: Rune Allnor

Date: 15 Apr, 2011 09:38:47

Message: 6 of 9

On Apr 11, 3:37 pm, "Johannes Korsawe"
<johannes.korsawe.nos...@volkswagen.de> wrote:
> Dear community,
>
> i have the following task.
>
> Let xy be a closed curve in 2d around [0,0] with size(xy)=[n1 2].
> Let a be a vector of angles in between 0 and 2*pi and r = [cos(a),sin(a)] a vector of polar directions with size(r) = [n2 2];
>
> I am looking for a fast (!) algorithm to compute the norm of the outmost intersections (there could be multiple) of xy with the rays r (consider r beeing prolongated infinitely from [0,0]). Using different words, i want to compute approximately the polar coordinates of the curve xy according to the angles.

This is a classical problem in optical ray tracing.
Very fast algorithms exist, but they require the data
to be stored in particular data structures for the
most efficient algorithms to work.

So you have two choises: Either stick with the slow
matlab, or re-formulate your problem to one of optical
raytracing and use some of the efficient raytracer
engines.

Rune

Subject: Speedup: Fast polar intersection

From: Roger Stafford

Date: 16 Apr, 2011 03:05:20

Message: 7 of 9

"Roger Stafford" wrote in message <io8oa4$gc5$1@fred.mathworks.com>...
> I finally got back to your problem, Johannes. I realize now that there is no ....
- - - - - - - -
  I vectorized the logic that was in the previous for-loop. Perhaps it will help the speed. At least it looks a little neater.

n = size(a,1);
r = zeros(n,1);
[as,p] = sort(a);
ca = cos(as); sa = sin(as);
s = repmat((1:n).',2,1);
m = size(xy,1);
x = xy(:,1); y = xy(:,2);
dx = diff(x); dy = diff(y);
num = x(1:m-1).*dy-y(1:m-1).*dx;
axy = atan2(y,x);
[~,b] = histc(axy,[-inf;as;inf]);
t1 = abs(diff(axy)) >= pi;
b1 = b(1:m-1); b2 = b(2:m);
t2 = ((b1<b2)==t1);
c1 = t2.*(b2-b1)+b1;
c2 = t2.*(b1-b2)+b2+t1.*n-1;
m1 = 1:m-1;
for k = m1(c1<=c2)
 t = s(c1(k):c2(k));
 u = num(k)./(ca(t)*dy(k)-sa(t)*dx(k));
 r(t) = max(r(t),u);
end
r(p) = r;

Roger Stafford

Subject: Speedup: Fast polar intersection

From: Johannes Korsawe

Date: 18 Apr, 2011 09:22:04

Message: 8 of 9

Roger,

thank you for your help and thoughts.

In comparing your two versions, the vectorized version took considerably longer than the first version. Sometimes the jit-acceleration inside MATLAb is the winner.

But even the faster version did not beat my own implementation after your ideas, IF there are not multiple cuts of a ray through the polyline. This could also depend from the test scenario. Maybe if i have a longer list of rays than list of polylines' nodes, it is the other way round. I will check this and give report...

So by now for my applications, my version of choice is using my own implementation and using your first version as a fallback if i detect several crossings of a ray.

Thank you again and best regards,
Johannes

Subject: Speedup: Fast polar intersection

From: Johannes Korsawe

Date: 18 Apr, 2011 14:02:20

Message: 9 of 9

Roger,

i made some more tests.

The remarks in my last posts only apply to small sizes of points in curve xy and number of rays in r (~300).

If you use more rays (~3000), your first version beats my version clearly. The same applies if you use a larger numbers of points (~1000) in the discretized curve xy.

Just your second version is always slower than your first.

So thats the outcome if noone else has additional ideas. Thank you so far!

Best,
Johannes

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us