/
rhoRange.m
65 lines (44 loc) · 1.26 KB
/
rhoRange.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
function [rhoMin,rhoMax] = rhoRange(sR)
% compute range of the polar angle of a spherical region
if isempty(sR.N) || ...
(isscalar(sR.N) && ~isnull(sR.N.z) && sR.alpha == 0)
rhoMin = 0;
rhoMax = 2*pi;
return
end
% antipodal should not increase spherical region
sR.antipodal = false;
% discretisation
omega = linspace(0,2*pi,361);
% start with equator
v = vector3d('theta',pi/2,'rho',omega);
rho = v.rho(sR.checkInside(v));
% cylce through boundary
for i = 1:length(sR.N)
b = vector3d('theta',acos(sR.alpha(i)),'rho',omega);
rot = rotation.map(zvector,sR.N(i));
b = rot * b;
% remove points close to north and south
% as we can not determine rho very well
b(abs(b.theta-pi/2)>pi/2-1e-4)= [];
rho = [rho,b.rho(sR.checkInside(b))]; %#ok<AGROW>
end
rho = sort(mod(rho,2*pi));
% find all holes in the rho list
ind = find(abs(diff(rho))>20*degree);
% maybe there is also a hole in zero
if rho(end)-rho(1)<340*degree, ind = [ind, length(rho)]; end
if isempty(ind)
if isempty(rho)
rhoMin = 0;
rhoMax = 2*pi;
else
rhoMin = min(rho);
rhoMax = max(rho);
end
else
rhoMin = rho(mod(ind([end,1:end-1]),length(rho))+1);
rhoMax = rho(ind);
ind = rhoMin > rhoMax;
rhoMin(ind) = rhoMin(ind) - 2*pi;
end