|
tp@0
|
1 function [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_pointtoedge(fromcoordsshortlist,reftoshortlist,tocoords,canplaneobstruct,planeseesplane,...
|
|
tp@0
|
2 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
|
|
tp@0
|
3 % EDB1checkobstr_pointtoedge - Checks obstruction for a list of paths from points to edges.
|
|
tp@0
|
4 % Checks if the paths from N starting points (reftoshortlist -> fromcoordsshortlist) to
|
|
tp@0
|
5 % N ending points on edges (tocoords) pass through any of M planes (those marked by 1 in
|
|
tp@0
|
6 % canplaneobstruct), i.e., if they are obstructed.
|
|
tp@0
|
7 %
|
|
tp@0
|
8 % Input parameters:
|
|
tp@0
|
9 % fromcoordsshortlist Matrix with the unique coordinates of the N
|
|
tp@0
|
10 % starting points.
|
|
tp@0
|
11 % reftoshortlist List, [N,1], with the N indices to fromcoordsshortlist.
|
|
tp@0
|
12 % tocoords Matrix, [N,3], with the coordinates of the
|
|
tp@0
|
13 % N ending points (on edges) of the N paths to check.
|
|
tp@0
|
14 % canplaneobstruct,planeseesplane,planeeqs,planenvecs,minvals,maxvals,...
|
|
tp@0
|
15 % planecorners,corners,ncornersperplanevec,rearsideplane
|
|
tp@0
|
16 % Data that should have been passed on from the
|
|
tp@0
|
17 % eddatafile. See EDB1edgeo for more information.
|
|
tp@0
|
18 %
|
|
tp@0
|
19 % Output parameters:
|
|
tp@0
|
20 % nonobstructedpaths A list, [N_nobstructions,1] of the indices of paths that are not
|
|
tp@0
|
21 % obstructed. The indices refer to the input
|
|
tp@0
|
22 % lists fromcoords, tocoords, startplanes
|
|
tp@0
|
23 % endplanes
|
|
tp@0
|
24 % nobstructions The number of paths (of the N input paths) that
|
|
tp@0
|
25 % were not obstructed.
|
|
tp@0
|
26 % edgehits A list, [N_edgehits,1] of the indicies of paths that hit a
|
|
tp@0
|
27 % plane right at an edge.
|
|
tp@0
|
28 % cornerhits A list, [N_cornerhits,1] of the indicies of paths that hit a
|
|
tp@0
|
29 % plane right at a corner.
|
|
tp@0
|
30 %
|
|
tp@0
|
31 % Uses function EDB1poinpla
|
|
tp@0
|
32 %
|
|
tp@0
|
33 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
34 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
35 %
|
|
tp@0
|
36 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
37 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
38 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
39 %
|
|
tp@0
|
40 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
41 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
42 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
43 %
|
|
tp@0
|
44 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
45 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
46 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
47 % Peter Svensson (svensson@iet.ntnu.no) 20050922
|
|
tp@0
|
48 %
|
|
tp@0
|
49 % [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_pointtoedge(fromcoordsshortlist,reftoshortlist,tocoords,canplaneobstruct,planeseesplane,...
|
|
tp@0
|
50 % planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
|
|
tp@0
|
51
|
|
tp@0
|
52 npaths = size(tocoords,1);
|
|
tp@0
|
53 nplanes = size(canplaneobstruct,2);
|
|
tp@0
|
54
|
|
tp@0
|
55 % We only need to check planes for which canplaneobstruct = 1.
|
|
tp@0
|
56
|
|
tp@0
|
57 if nplanes < 256
|
|
tp@0
|
58 maxlistofplanestocheck = uint8( find( sum(canplaneobstruct) ).' );
|
|
tp@0
|
59 elseif nplanes < 65536
|
|
tp@0
|
60 maxlistofplanestocheck = uint16( find( sum(canplaneobstruct) ).' );
|
|
tp@0
|
61 else
|
|
tp@0
|
62 maxlistofplanestocheck = uint32( find( sum(canplaneobstruct) ).' );
|
|
tp@0
|
63 end
|
|
tp@0
|
64
|
|
tp@0
|
65 nplanestocheck = length(maxlistofplanestocheck);
|
|
tp@0
|
66
|
|
tp@0
|
67 if nplanestocheck > 0
|
|
tp@0
|
68 ntot = nplanestocheck*npaths;
|
|
tp@0
|
69
|
|
tp@0
|
70 nmaxcorners = double(max(ncornersperplanevec(maxlistofplanestocheck)));
|
|
tp@0
|
71 planecorners = planecorners(:,1:nmaxcorners);
|
|
tp@0
|
72
|
|
tp@0
|
73 % For the obstruction test, we create expanded matrices/lists where the
|
|
tp@0
|
74 % plane-number is the fastest varying variable. We calculate
|
|
tp@0
|
75 % the direction vectors, from point to point, before the expanding.
|
|
tp@0
|
76
|
|
tp@0
|
77 dirvec = tocoords - fromcoordsshortlist(reftoshortlist,:);
|
|
tp@0
|
78 clear tocoords
|
|
tp@0
|
79 dirveclengths = sqrt(sum(dirvec.'.^2)).' + eps*100;
|
|
tp@0
|
80 dirvec = dirvec./dirveclengths(:,ones(1,3));
|
|
tp@0
|
81
|
|
tp@0
|
82 % Create the expanded matrices, by duplicating each path with the
|
|
tp@0
|
83 % number of planes that need to be checked obstruction against.
|
|
tp@0
|
84 % Immediately select only those individual combinations where
|
|
tp@0
|
85 % canplaneobstruct = 1. This selection is done both for bigplanelist
|
|
tp@0
|
86 % and for the reftoshortlist, and for the reftodirvec (which points to
|
|
tp@0
|
87 % the shorter dirvec matrix).
|
|
tp@0
|
88
|
|
tp@0
|
89 iv1 = (find(canplaneobstruct(:,maxlistofplanestocheck)));
|
|
tp@0
|
90
|
|
tp@0
|
91 bigplanelist = maxlistofplanestocheck(ceil((iv1)/npaths));
|
|
tp@0
|
92 if ntot < 65536
|
|
tp@0
|
93 iv1 = uint16(iv1);
|
|
tp@0
|
94 elseif ntot < 4.29e9
|
|
tp@0
|
95 iv1 = uint32(iv1);
|
|
tp@0
|
96 end
|
|
tp@0
|
97
|
|
tp@0
|
98 reftoshortlist = reftoshortlist(:,ones(1,nplanestocheck));
|
|
tp@0
|
99 reftoshortlist = reftoshortlist(iv1);
|
|
tp@0
|
100
|
|
tp@0
|
101 if npaths*nplanestocheck < 65536
|
|
tp@0
|
102 reftodirvec = uint16([1:npaths].');
|
|
tp@0
|
103 else
|
|
tp@0
|
104 reftodirvec = uint32([1:npaths].');
|
|
tp@0
|
105 end
|
|
tp@0
|
106 reftodirvec = reftodirvec(:,ones(1,nplanestocheck));
|
|
tp@0
|
107 reftodirvec = reftodirvec(iv1);
|
|
tp@0
|
108
|
|
tp@0
|
109 %--------------------------------------------------------------------------
|
|
tp@0
|
110 % Tests
|
|
tp@0
|
111 % (Earlier, it was first tested if the IS and R were on the same side of a plane. This is checked at another place).
|
|
tp@0
|
112 % (Earlier, test 2 was: Are the vectors IS -> R parallel to the planes?
|
|
tp@0
|
113 % This shold not be needed now).
|
|
tp@0
|
114 %
|
|
tp@0
|
115 % Are the hitpoints inside the planes?
|
|
tp@0
|
116 % 1. Calculate the hit points' coordinates.
|
|
tp@0
|
117 % 2. Check that the hit point is not further away than the receiver
|
|
tp@0
|
118 % point, and that the hit point is not in the wrong direction!
|
|
tp@0
|
119 % 3. Check if the hit points are inside the plane that is checked.
|
|
tp@0
|
120
|
|
tp@0
|
121 % Step 1
|
|
tp@0
|
122 udir_bigplanenvecs = ( planeeqs(bigplanelist,4) - sum( planenvecs(bigplanelist,:).'.*(fromcoordsshortlist(reftoshortlist,:).') ).' )./( sum( planenvecs(bigplanelist,:).'.*(dirvec(reftodirvec,:).') ).');
|
|
tp@0
|
123
|
|
tp@0
|
124 % Step 2
|
|
tp@0
|
125 iv2 = find( udir_bigplanenvecs<dirveclengths(reftodirvec) & udir_bigplanenvecs>0 );
|
|
tp@0
|
126 clear dirveclengths
|
|
tp@0
|
127 iv1 = iv1(iv2);
|
|
tp@0
|
128 bigplanelist = bigplanelist(iv2);
|
|
tp@0
|
129 reftoshortlist = reftoshortlist(iv2);
|
|
tp@0
|
130 reftodirvec = reftodirvec(iv2);
|
|
tp@0
|
131 udir_bigplanenvecs = udir_bigplanenvecs(iv2,:);
|
|
tp@0
|
132 clear iv2
|
|
tp@0
|
133
|
|
tp@0
|
134 if ~isempty(iv1)
|
|
tp@0
|
135 % Step 3
|
|
tp@0
|
136 xhit_list = fromcoordsshortlist(reftoshortlist,:) + udir_bigplanenvecs(:,ones(1,3)).*dirvec(reftodirvec,:);
|
|
tp@0
|
137
|
|
tp@0
|
138 clear udir_bigplanenvecs dirvec reftoshortlist
|
|
tp@0
|
139
|
|
tp@0
|
140
|
|
tp@0
|
141 [hitvec,edgehit,cornerhit] = EDB1poinpla(xhit_list,bigplanelist,minvals,maxvals,planecorners,corners,ncornersperplanevec,planenvecs);
|
|
tp@0
|
142 clear xhit_list
|
|
tp@0
|
143
|
|
tp@0
|
144 % We recycle canplaneobstruct since we need to compile how many hits
|
|
tp@0
|
145 % that occurred.
|
|
tp@0
|
146 canplaneobstruct = zeros(size(canplaneobstruct(:,maxlistofplanestocheck)));
|
|
tp@0
|
147 ivhits = find(hitvec);
|
|
tp@0
|
148 if ~isempty(ivhits)
|
|
tp@0
|
149 indexvec = iv1(ivhits);
|
|
tp@0
|
150 canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
|
|
tp@0
|
151 end
|
|
tp@0
|
152 edgehits = [];
|
|
tp@0
|
153 ivhits = find(edgehit);
|
|
tp@0
|
154 if ~isempty(ivhits)
|
|
tp@0
|
155 disp('WARNING: Edge hit. This is treated as an obstruction.')
|
|
tp@0
|
156 indexvec = iv1(ivhits);
|
|
tp@0
|
157 canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
|
|
tp@0
|
158 end
|
|
tp@0
|
159
|
|
tp@0
|
160 cornerhits = [];
|
|
tp@0
|
161 ivhits = find(cornerhit);
|
|
tp@0
|
162 if ~isempty(ivhits)
|
|
tp@0
|
163 disp('WARNING: Corner hit. This is treated as an obstruction.')
|
|
tp@0
|
164 indexvec = iv1(ivhits);
|
|
tp@0
|
165 canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
|
|
tp@0
|
166 end
|
|
tp@0
|
167
|
|
tp@0
|
168 % Here was a bug: erroneous summing when there was a single plane that
|
|
tp@0
|
169 % could obstruct (PS 041127)
|
|
tp@0
|
170 if size(canplaneobstruct,2) > 1
|
|
tp@0
|
171 canplaneobstruct = 1-sign(sum(canplaneobstruct.'));
|
|
tp@0
|
172 else
|
|
tp@0
|
173 canplaneobstruct = 1-canplaneobstruct;
|
|
tp@0
|
174 end
|
|
tp@0
|
175 nonobstructedpaths = find(canplaneobstruct);
|
|
tp@0
|
176
|
|
tp@0
|
177 nobstructions = npaths-length(nonobstructedpaths);
|
|
tp@0
|
178 else
|
|
tp@0
|
179 nonobstructedpaths = [1:npaths].';
|
|
tp@0
|
180 nobstructions = 0;
|
|
tp@0
|
181 edgehits = [];
|
|
tp@0
|
182 cornerhits = [];
|
|
tp@0
|
183 end
|
|
tp@0
|
184 else
|
|
tp@0
|
185 nonobstructedpaths = [1:npaths].';
|
|
tp@0
|
186 nobstructions = 0;
|
|
tp@0
|
187 edgehits = [];
|
|
tp@0
|
188 cornerhits = [];
|
|
tp@0
|
189 end
|
