Vectorize my code in azimuth calculation.

2014 9 23
2 个回答
回答已采纳
更新时间：2014 9 25
7 次查看（30 天）

0 个投票

Hi I have lots of points(in Latitude & Longitude) and I want to calculate the azimuth of each pairs of them. I use code as attached file 'AzimuthCalc.m'. when I run my function, it's so time consuming and it's elapsed time is 109.371341 seconds. I'm sure the problem caused by my non-vectorize code.
So, would you please check my function and help me to improve it?
Thanks a lot Mani

2 个评论
显示 无 隐藏 无

Adam
Adam 2014-9-23
Run:
profile on
your code
profile off
profile viewer
and find out where your code is slow. No point in guessing when you can use facts!
Mani Ahmadian
Mani Ahmadian 2014-9-23
Dear Adam,
Thanks so much. I guess the problem is here(nested for):
for ii=1:length(lat)
for jj=ii+1:length(lat)
AzimuthPoints(ii,jj)=azimuth('rh',lat(ii,1),long(ii,1),lat(jj,1),long(jj,1));
end
end
If I can change my simple code to vectorize one, problem will solve. So, I'm trying to find a way to do that, so plz help me.
Thanks
Mnai

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Guillaume
Guillaume 2014-9-23

0 个投票

I don't have the mapping toolbox, so I'm just relying on the documentation to write my answer.
For a start, you can replace the inner loop with:
AzimuthPoints(ii,ii+1:end)=azimuth('rh',lat(ii),long(ii),lat(ii+1:end),long(ii+1:end));
I believe the following may work to replace both loops (and the AzimuthPoints declaration)
[lat1, lat2] = meshgrid(lat);
[long1, long2] = meshgrid(long);
AzimuthPoints = azimuth('rh', lat1, long1, lat2, long2);
The downside of this method is that it's going to calculate the whole AzimuthPoints matrix, not just the upper triangle, so it may actually be slower.
You could do:
lat1 = triu(lat1, 1);
%same with lat2, long1 and long2
to replace the lower triangle and main diagonal with zero, on the assumption that calculating
azimuth('rh', 0, 0, 0, 0)
is faster.

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Mani Ahmadian
Mani Ahmadian 2014-9-24
Dear Guillaume;
Your answer shocked me, It's great and is very useful for me. I run it and elapsed time is about 0.222 seconds. This time is acceptable for my project and I used your offer in my code.
Another stage is to find all azimuths are between 180+-30 (150<=azimuths and azimuths<=210). I use my code as bellow to find all azimuths are true, it's so fast as I need:
TrueAzimuth=zeros(size(AzimuthPoints));
TrueAzimuth(AzimuthPoints()<=UserAz+AzTol & AzimuthPoints()>=UserAz-AzTol)=1;
As final part I should pick all lat/long values that it's corresponding value in TrueAzimuth is 1. For example if TrueAzimuth(1,3)=1, I pick lat/long values of point#1 and point#3 for future calculations. To select points, I use my simple code as bellow, But it's so slow.
for ii=1:length(lat)
for jj=1:length(long)
if TrueAzimuth(ii,jj)==1
TrueLat=vertcat(TrueLat, lat(ii), lat(jj));
TrueLong=vertcat(TrueLong, long(ii), long(jj));
end
end
end
Please help me to vectorize it, too.
Thanks a lot.
Mani
Guillaume
Guillaume 2014-9-24
First, TrueAzimuth is just
TrueAzimuth = AzimuthPoints<=UserAz+AzTol & AzimuthPoints>=UserAz-AzTol; %No need to predeclare.
Then, I believe the following will give you exactly the same result as your loop:
TrueLat = [lat1(TrueAzimuth) lat2(TruAzimuth)]';
TrueLat = TrueLat(:);
%same with TrueLong
But the following seems to make more sense to me, if you're going to make further processing:
TrueLat = [lat1(TrueAzimuth) lat2(TrueAzimuth)]; %keep as 2 column matrix
Mani Ahmadian
Mani Ahmadian 2014-9-24
I use your code, but because we're useing TtueAzimuth(consist of 0,1) to index lat1,long1,... I receive an error message: Subscript indices must either be real positive integers or logicals.
How can I solve it?
Thanks so much
Mani
Guillaume
Guillaume 2014-9-24
TrueAzimuth should be 0s and 1s but of type logical. If you used my code to generate it, it shoul already be logical, but if not:
TrueAzimuth = logical(TrueAzimuth);
Mani Ahmadian
Mani Ahmadian 2014-9-25
Thanks, It works correctly.

请先登录，再进行评论。

更多回答（1 个）

Roger Stafford
Roger Stafford 2014-9-24

0 个投票

You could try this. It calls on 'azimuth' in the same sequence as in your for-loop code, but it is vectorized. Whether it is faster depends on whether the 'tril' and 'repmat' calls take less time than the overhead in the for-loops and individual 'azimuth' calls in your code.
n = length(lat);
p = tril(repmat(1:n,n,1),-1);
p = p(p>0);
q = tril(repmat((1:n)',1,n),-1);
q = q(q>0);
AzimuthPoints = zeros(n);
AzimuthPoints(p+n*(q-1)) = azimuth('rh',lat(p),long(p),lat(q),long(q));

类别

在 帮助中心 和 File Exchange 中查找有关 Get Started with Mapping Toolbox 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by