|
tp@0
|
1 function outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy)
|
|
tp@0
|
2 % EDB1edgeox - Calculates some plane- and edge-related geometrical parameters.
|
|
tp@0
|
3 % EDB1edgeo calculates some plane- and edge-related geometrical parameters
|
|
tp@0
|
4 % based on corners and planes in a .mat-file generated by EDB1readcad
|
|
tp@0
|
5 % but only data that is independent of the source and receiver.
|
|
tp@0
|
6 % The output is saved in a .mat-file.
|
|
tp@0
|
7 %
|
|
tp@0
|
8 % Input parameters:
|
|
tp@0
|
9 % cadgeofile (optional) The .mat file with the corners and plane data.
|
|
tp@0
|
10 % If it is not specified, a file opening window is presented.
|
|
tp@0
|
11 % outputfile (optional) The .mat file where the output data will be stored.
|
|
tp@0
|
12 % If it is not specified, an automatic file name is constructed
|
|
tp@0
|
13 % using the same filename stem as the cadgeofile, but ending
|
|
tp@0
|
14 % with '_edgeo'.
|
|
tp@0
|
15 % open_or_closed_model 'o' or 'c' defining whether the model is open or closed.
|
|
tp@0
|
16 % int_or_ext_model 'i' or 'e' defining whether the model is interior or exterior.
|
|
tp@0
|
17 % specorder The highest order of any reflection kind (specular and/or diffraction).
|
|
tp@0
|
18 % difforder (optional) The highest order of diffraction. If it is 0 or 1 then the parameter
|
|
tp@0
|
19 % edgeseespartialedge is not calculated. Default: 1
|
|
tp@0
|
20 % nedgesubs (optional) The number of segments that each edge will be
|
|
tp@0
|
21 % subdivided into for visibility/obstruction tests. Default: 2
|
|
tp@0
|
22 % NB! When nedgesubs = 2, only the two end points will be checked.
|
|
tp@0
|
23 % firstskipcorner (optional) All edges including at least one corner with this number or
|
|
tp@0
|
24 % higher will be excluded from the calculations. Default: 1e6
|
|
tp@0
|
25 % planeseesplanestrategy (optional) If this parameter is given the value 1, a plane-to-plane
|
|
tp@0
|
26 % visibility check is done by checking the plane-midpoint to plane-midpoint
|
|
tp@0
|
27 % for obstructions, which might be enough for some geometries.
|
|
tp@0
|
28 % SHOWTEXT (global) An integer value; if it is 4 or higher, then messages will be printed on the
|
|
tp@0
|
29 % screen during calculations.
|
|
tp@0
|
30 % Output parameters:
|
|
tp@0
|
31 % outputfile The name of the outputfile, generated either automatically or specified as
|
|
tp@0
|
32 % an input parameter.
|
|
tp@0
|
33 %
|
|
tp@0
|
34 % Output parameters in the outputfile:
|
|
tp@0
|
35 % (Taken directly from the cadgeo-file:)
|
|
tp@0
|
36 % corners planecorners planeeqs planenvecs ncornersperplanevec
|
|
tp@0
|
37 % minvals maxvals planehasindents indentingcorners
|
|
tp@0
|
38 % (Taken directly from the setup-file:)
|
|
tp@0
|
39 % int_or_ext_model
|
|
tp@0
|
40 % (Calculated by EDB1edgeo:)
|
|
tp@0
|
41 % edgecorners Matrix [nedges,2] with the corner numbers that
|
|
tp@0
|
42 % define each edge, in the reference order.
|
|
tp@0
|
43 % edgestartcoords Matrix [nedges,3] with the coordinates of the
|
|
tp@0
|
44 % start corner for each edge.
|
|
tp@0
|
45 % edgeendcoords Matrix [nedges,3] with the coordinates of the
|
|
tp@0
|
46 % end corner for each edge.
|
|
tp@0
|
47 % edgelengthvec List [nedges,1] with the lengths of each edge.
|
|
tp@0
|
48 % planesatedge Matrix [nedges,2] with the plane numbers of the
|
|
tp@0
|
49 % two planes that connect to each edge. The first
|
|
tp@0
|
50 % plane is considered as the reference plane for
|
|
tp@0
|
51 % each edge. Rule: if the RH thumb is placed
|
|
tp@0
|
52 % along the edge, from start to end, and the hand
|
|
tp@0
|
53 % is rotated, the fingers should "come up
|
|
tp@0
|
54 % through" the reference plane.
|
|
tp@0
|
55 % closwedangvec List [nedges,1] with the angles in radians of
|
|
tp@0
|
56 % each edge - for the closed part of each edge.
|
|
tp@0
|
57 % edgenvecs Matrix [nedges,3] with the normal vectors of
|
|
tp@0
|
58 % the reference plane for each edge.
|
|
tp@0
|
59 % offedges List [nlength,1] where nlength = 0-nedges
|
|
tp@0
|
60 % containing the edge numbers that should not
|
|
tp@0
|
61 % used.
|
|
tp@0
|
62 % reflfactors List [nplanes,1] with the reflection factors
|
|
tp@0
|
63 % of each plane. The only supported values are
|
|
tp@0
|
64 % +1 for rigid planes, 0 for perfectlt absorbing
|
|
tp@0
|
65 % planes and -1 for pressure-release planes. The
|
|
tp@0
|
66 % pressure-release cases might not be implemented
|
|
tp@0
|
67 % fully.
|
|
tp@0
|
68 % planeisthin List [nplanes,1] with the value 1 or 0 stating
|
|
tp@0
|
69 % whether the plane is thin or not.
|
|
tp@0
|
70 % rearsideplane List [nplanes,1] with the plane number that is at
|
|
tp@0
|
71 % the rear side. If a plane is non-thin, the
|
|
tp@0
|
72 % value in this list is zero.
|
|
tp@0
|
73 % planehasindents List [nplanes,1] with the value 1 or 0 stating
|
|
tp@0
|
74 % whether a plane has indents or not.
|
|
tp@0
|
75 % indentingedgepairs Matrix [nindentingedgepairs,2] where each row
|
|
tp@0
|
76 % gives two connected edges that form an indent
|
|
tp@0
|
77 % in the shared plane. Consequently, those two
|
|
tp@0
|
78 % edges can never see each other.
|
|
tp@0
|
79 % canplaneobstruct List [nplanes,1] with the value 1 or 0 stating
|
|
tp@0
|
80 % whether a plane has the potential to obstruct
|
|
tp@0
|
81 % or not.
|
|
tp@0
|
82 % planeseesplane A matrix, [nplanes,nplanes], with the value 1
|
|
tp@0
|
83 % 0,-1,-2 stating:
|
|
tp@0
|
84 % 1 A plane is front of another plane and could
|
|
tp@0
|
85 % then potentially see it.
|
|
tp@0
|
86 % 0 A plane is completely behind another plane
|
|
tp@0
|
87 % but they are not aligned
|
|
tp@0
|
88 % -1 A plane is aligned with another plane
|
|
tp@0
|
89 % -2 A plane is back-to-back with another plane.
|
|
tp@0
|
90 % edgesatplane A matrix, [nplanes,nmaxedgesperplane] which for
|
|
tp@0
|
91 % row N contains the edge numbers that connect to
|
|
tp@0
|
92 % to plane N.
|
|
tp@0
|
93 % edgeseesplane A matrix, [nplanes,nedges], which contains the
|
|
tp@0
|
94 % values 0,1,-2 or -2, meaning:
|
|
tp@0
|
95 % 0 The edge is completely behind a plane
|
|
tp@0
|
96 % -1 The edge can never see the plane and the
|
|
tp@0
|
97 % edge is aligned with the plane
|
|
tp@0
|
98 % -2 The edge can never see the plane and the
|
|
tp@0
|
99 % edge belongs to the plane
|
|
tp@0
|
100 % 1 The edge has at least one point in front of the plane.
|
|
tp@0
|
101 %
|
|
tp@0
|
102 % Uses the functions EDB1coordtrans1
|
|
tp@0
|
103 % EDB1infrontofplane, EDB1strpend, EDB1cross
|
|
tp@0
|
104 %
|
|
tp@0
|
105 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
106 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
107 %
|
|
tp@0
|
108 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
109 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
110 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
111 %
|
|
tp@0
|
112 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
113 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
114 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
115 %
|
|
tp@0
|
116 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
117 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
118 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
119 % Peter Svensson (svensson@iet.ntnu.no) 20090709
|
|
tp@0
|
120 %
|
|
tp@0
|
121 % outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy);
|
|
tp@0
|
122
|
|
tp@0
|
123 global SHOWTEXT
|
|
tp@0
|
124
|
|
tp@0
|
125 geomacc = 1e-10;
|
|
tp@0
|
126
|
|
tp@0
|
127 if nargin < 9
|
|
tp@0
|
128 planeseesplanestrategy = 0;
|
|
tp@0
|
129 if nargin < 8
|
|
tp@0
|
130 firstskipcorner = 1e6;
|
|
tp@0
|
131 if nargin < 7
|
|
tp@0
|
132 nedgesubs = 2;
|
|
tp@0
|
133 if nargin < 6
|
|
tp@0
|
134 difforder = 1;
|
|
tp@0
|
135 else
|
|
tp@0
|
136 if isempty(difforder), difforder = 1; end
|
|
tp@0
|
137 end
|
|
tp@0
|
138 else
|
|
tp@0
|
139 if isempty(nedgesubs), nedgesubs = 2; end
|
|
tp@0
|
140 end
|
|
tp@0
|
141 else
|
|
tp@0
|
142 if isempty(firstskipcorner), firstskipcorner = 1e6; end
|
|
tp@0
|
143 end
|
|
tp@0
|
144 end
|
|
tp@0
|
145
|
|
tp@0
|
146 %---------------------------------------------------------------
|
|
tp@0
|
147 % If no cadgeofile was specified, present a file opening window
|
|
tp@0
|
148
|
|
tp@0
|
149 if isempty(cadgeofile)
|
|
tp@0
|
150 [cadgeofile,filepath] = uigetfile('*.mat','Please select the cadgeofile');
|
|
tp@0
|
151 [filepath,temp1,temp2] = fileparts(filepath);
|
|
tp@0
|
152 if ~isstr(cadgeofile)
|
|
tp@0
|
153 return
|
|
tp@0
|
154 end
|
|
tp@0
|
155 [temp1,filestem,fileext,temp2] = fileparts(cadgeofile);
|
|
tp@0
|
156 cadgeofile = [filepath,filestem,'.mat'];
|
|
tp@0
|
157 else
|
|
tp@0
|
158 [filepath,filestem,fileext] = fileparts(cadgeofile);
|
|
tp@0
|
159 if ~isempty(filepath)
|
|
tp@0
|
160 cadgeofile = [[filepath,filesep],filestem,'.mat'];
|
|
tp@0
|
161 else
|
|
tp@0
|
162 cadgeofile = [filestem,'.mat'];
|
|
tp@0
|
163 end
|
|
tp@0
|
164 end
|
|
tp@0
|
165
|
|
tp@0
|
166
|
|
tp@0
|
167 %---------------------------------------------------------------
|
|
tp@0
|
168 % If no output file was specified, construct an automatic file name
|
|
tp@0
|
169
|
|
tp@0
|
170 if isempty(outputfile)
|
|
tp@0
|
171 filestem = EDB1strpend(filestem,'_cadgeo');
|
|
tp@0
|
172 outputfile = [[filepath,filesep],filestem,'_eddata.mat'];
|
|
tp@0
|
173 end
|
|
tp@0
|
174
|
|
tp@0
|
175 %---------------------------------------------------------------
|
|
tp@0
|
176
|
|
tp@0
|
177 eval(['load ',cadgeofile])
|
|
tp@0
|
178
|
|
tp@0
|
179 ncorners = length(cornernumbers);
|
|
tp@0
|
180 [nplanes,maxncornersperplane] = size(planecorners);
|
|
tp@0
|
181 minncornersperplane = min(ncornersperplanevec);
|
|
tp@0
|
182 nedgesperplanevec = double(ncornersperplanevec);
|
|
tp@0
|
183 ncornersperplanevec = double(ncornersperplanevec);
|
|
tp@0
|
184
|
|
tp@0
|
185 switchoffcorners = find( cornernumbers >= firstskipcorner );
|
|
tp@0
|
186
|
|
tp@0
|
187 onesvec2 = ones(nplanes,1);
|
|
tp@0
|
188 zerosvec1 = zeros(nplanes,1);
|
|
tp@0
|
189 zerosvec2 = int8(zeros(nplanes,nplanes));
|
|
tp@0
|
190
|
|
tp@0
|
191 %----------------------------------------------------------------------------
|
|
tp@0
|
192 %
|
|
tp@0
|
193 % GEOMETRICAL ACOUSTICS PARAMETERS
|
|
tp@0
|
194 %
|
|
tp@0
|
195 %----------------------------------------------------------------------------
|
|
tp@0
|
196 %
|
|
tp@0
|
197 % Go through all planeabstypes. If it is:
|
|
tp@0
|
198 % 'SOFT' or 'soft' then reflfactors should be -1.
|
|
tp@0
|
199 % 'TOTA' or 'tota' then reflfactors should be 0 (short for TOTABS)
|
|
tp@0
|
200 % Otherwise it is assumed to be rigid, so reflfactors = 1.
|
|
tp@0
|
201
|
|
tp@0
|
202 if SHOWTEXT >= 4
|
|
tp@0
|
203 disp(' Check abstypes')
|
|
tp@0
|
204 end
|
|
tp@0
|
205
|
|
tp@0
|
206 [slask,nchars] = size(planeabstypes);
|
|
tp@0
|
207 reflfactors = ones(nplanes,1);
|
|
tp@0
|
208 if nchars >= 4
|
|
tp@0
|
209 planeabstypes = lower(full(planeabstypes(:,1:min([4,nchars]))));
|
|
tp@0
|
210
|
|
tp@0
|
211 comptxt = 'soft';
|
|
tp@0
|
212 ivpotential = find(planeabstypes(:,1)==comptxt(1));
|
|
tp@0
|
213 if ~isempty(ivpotential)
|
|
tp@0
|
214 comptxt = 'oft';
|
|
tp@0
|
215 compmat = comptxt(ones(length(ivpotential),1),:);
|
|
tp@0
|
216 ivsoft = ivpotential(find(prod( (planeabstypes(ivpotential,2:4).*compmat).' ).'));
|
|
tp@0
|
217 reflfactors(ivsoft) = -1*ones(size(ivsoft));
|
|
tp@0
|
218 end
|
|
tp@0
|
219
|
|
tp@0
|
220 comptxt = 'tota';
|
|
tp@0
|
221 ivpotential = find(planeabstypes(:,1)==comptxt(1));
|
|
tp@0
|
222 if ~isempty(ivpotential)
|
|
tp@0
|
223 comptxt = 'ota';
|
|
tp@0
|
224 compmat = comptxt(ones(length(ivpotential),1),:);
|
|
tp@0
|
225 ivtotabs = ivpotential(find(prod( double(planeabstypes(ivpotential,2:4)==compmat).' ).'));
|
|
tp@0
|
226 reflfactors(ivtotabs) = zeros(size(ivtotabs));
|
|
tp@0
|
227 end
|
|
tp@0
|
228 end
|
|
tp@0
|
229
|
|
tp@0
|
230 %----------------------------------------------------------------------------
|
|
tp@0
|
231 %
|
|
tp@0
|
232 % EDGE PARAMETERS
|
|
tp@0
|
233 %
|
|
tp@0
|
234 %----------------------------------------------------------------------------
|
|
tp@0
|
235 %
|
|
tp@0
|
236 % Check that the planecorners matrix contains no zeros. In such case
|
|
tp@0
|
237 % add the circular repetition of coordinates.
|
|
tp@0
|
238
|
|
tp@0
|
239 if SHOWTEXT >= 4
|
|
tp@0
|
240 disp(' Circular')
|
|
tp@0
|
241 end
|
|
tp@0
|
242 if sum(sum( planecorners == 0 )) ~= 0
|
|
tp@0
|
243 for ii = 1:nplanes
|
|
tp@0
|
244 iv = find( planecorners(ii,:) ~= 0);
|
|
tp@0
|
245 ncornersatplane = length(iv);
|
|
tp@0
|
246 if ncornersatplane ~= maxncornersperplane
|
|
tp@0
|
247 pattern = planecorners(ii,iv);
|
|
tp@0
|
248 nrepeatings = ceil(maxncornersperplane/ncornersatplane);
|
|
tp@0
|
249 for jj = 1:nrepeatings-1
|
|
tp@0
|
250 pattern = [pattern planecorners(ii,iv)];
|
|
tp@0
|
251 end
|
|
tp@0
|
252 planecorners(ii,:) = pattern(1:maxncornersperplane);
|
|
tp@0
|
253 end
|
|
tp@0
|
254 end
|
|
tp@0
|
255 end
|
|
tp@0
|
256
|
|
tp@0
|
257 %--------------------------------------------------------------------------------
|
|
tp@0
|
258 %
|
|
tp@0
|
259 % Find all edges
|
|
tp@0
|
260 %
|
|
tp@0
|
261 % Go through all planes. All consecutive pairs of corners form an edge.
|
|
tp@0
|
262 % The list planesatedge gives the one or two planes that connect
|
|
tp@0
|
263 % to each edge. If the second plane number, for a certain edge
|
|
tp@0
|
264 % is zero, then that specific edge is a plate edge.
|
|
tp@0
|
265
|
|
tp@0
|
266 edgesatplane = zeros( size(planecorners) );
|
|
tp@0
|
267
|
|
tp@0
|
268 if SHOWTEXT >= 4
|
|
tp@0
|
269 disp(' Defining edges...')
|
|
tp@0
|
270 end
|
|
tp@0
|
271
|
|
tp@0
|
272 nedgesguess = sum(ncornersperplanevec);
|
|
tp@0
|
273 edgecorners = zeros(nedgesguess,2);
|
|
tp@0
|
274 tempplanesatedge = zeros(nedgesguess,1);
|
|
tp@0
|
275 nedges = 0;
|
|
tp@0
|
276 thirdpoint = [];
|
|
tp@0
|
277 edgecounter = 0;
|
|
tp@0
|
278 multfac = round(10^(floor(log10(nedgesguess))+2));
|
|
tp@0
|
279
|
|
tp@0
|
280 % First go through all planes, and construct edges from the list of corners
|
|
tp@0
|
281 % that define every plane. This way, many edges will occur twice, unless
|
|
tp@0
|
282 % a thin plane has no rear side (then the edge occurs a single time)
|
|
tp@0
|
283 % or if two thin planes are connected (then the edge occurs four times)
|
|
tp@0
|
284
|
|
tp@0
|
285 for ii = 1:nplanes
|
|
tp@0
|
286 for jj = 1:nedgesperplanevec(ii)
|
|
tp@0
|
287 edgecounter = edgecounter + 1;
|
|
tp@0
|
288 corner1 = planecorners(ii,jj);
|
|
tp@0
|
289 if jj + 1 <= ncornersperplanevec(ii)
|
|
tp@0
|
290 corner2 = planecorners(ii,jj+1);
|
|
tp@0
|
291 else
|
|
tp@0
|
292 corner2 = planecorners(ii,1);
|
|
tp@0
|
293 end
|
|
tp@0
|
294 edgecorners(edgecounter,:) = [corner1 corner2];
|
|
tp@0
|
295 tempplanesatedge(edgecounter) = ii;
|
|
tp@0
|
296 end
|
|
tp@0
|
297
|
|
tp@0
|
298 end
|
|
tp@0
|
299
|
|
tp@0
|
300 if edgecounter < nedgesguess
|
|
tp@0
|
301 edgecorners(edgecounter+1:nedgesguess,:) = [];
|
|
tp@0
|
302 tempplanesatedge(edgecounter+1:nedgesguess,:) = [];
|
|
tp@0
|
303 end
|
|
tp@0
|
304
|
|
tp@0
|
305 % To find all duplicates, the edge list is sorted in numerical order
|
|
tp@0
|
306 % If an edge is defined twice or four times, then this is an OK edge.
|
|
tp@0
|
307 % If an edge is defined a single time, it is connected to a plane which is not connected, and then the
|
|
tp@0
|
308 % edge should be switched off.
|
|
tp@0
|
309 % If an edge is defined three times, there must be something like a reflector array where one of
|
|
tp@0
|
310 % the rear planes is missing.
|
|
tp@0
|
311
|
|
tp@0
|
312 [test,flipvec] = sort(edgecorners.');
|
|
tp@0
|
313 test = test.';
|
|
tp@0
|
314 flipvec = flipvec(1,:).';
|
|
tp@0
|
315
|
|
tp@0
|
316 test = test(:,1)*multfac + test(:,2);
|
|
tp@0
|
317 [test,sortvec] = sort(test);
|
|
tp@0
|
318 tempplanesatedge = tempplanesatedge(sortvec);
|
|
tp@0
|
319 flipvec = flipvec(sortvec);
|
|
tp@0
|
320 planesatedge = [tempplanesatedge.*(flipvec==1) tempplanesatedge.*(flipvec==2)];
|
|
tp@0
|
321
|
|
tp@0
|
322 % Check if there are loners; remove and give a warning.
|
|
tp@0
|
323
|
|
tp@0
|
324 dtest = diff([0;test;test(length(test))+100]);
|
|
tp@0
|
325 ntest = length(dtest);
|
|
tp@0
|
326 loners = find(dtest(1:ntest-1)~=0 & dtest(2:ntest)~=0);
|
|
tp@0
|
327
|
|
tp@0
|
328 if ~isempty(loners)
|
|
tp@0
|
329 disp('WARNING! Some edges had only one plane connected and they will be excluded.')
|
|
tp@0
|
330 disp(' Check if they should be corrected. ')
|
|
tp@0
|
331 disp(' Remember that reflectors must have a plane at the rear side too')
|
|
tp@0
|
332 disp(' Set SHOWTEXT to 4 or higher to get a list of these edges')
|
|
tp@0
|
333 if SHOWTEXT >= 4
|
|
tp@0
|
334 disp(' The edges are (EDB1 corner numbers for the two edge point are given):')
|
|
tp@0
|
335 nloners = length(loners);
|
|
tp@0
|
336 for jjdisp = 1:nloners
|
|
tp@0
|
337 indlon = loners(jjdisp);
|
|
tp@0
|
338 disp([' ',int2str( floor(test(indlon)/multfac) ),' ',int2str( test(indlon)-floor(test(indlon)/multfac)*multfac )])
|
|
tp@0
|
339 disp([' CATT: ',int2str( cornernumbers(floor(test(indlon)/multfac) )),' ',int2str( cornernumbers(test(indlon)-floor(test(indlon)/multfac)*multfac ))])
|
|
tp@0
|
340 end
|
|
tp@0
|
341 end
|
|
tp@0
|
342
|
|
tp@0
|
343 test(loners) = [];
|
|
tp@0
|
344 planesatedge(loners,:) = [];
|
|
tp@0
|
345
|
|
tp@0
|
346 end
|
|
tp@0
|
347
|
|
tp@0
|
348 ntest = length(test);
|
|
tp@0
|
349 if ntest >= 2
|
|
tp@0
|
350 if test(ntest) ~= test(ntest-1)
|
|
tp@0
|
351 test(ntest) = [];
|
|
tp@0
|
352 planesatedge(ntest,:) = [];
|
|
tp@0
|
353 end
|
|
tp@0
|
354 end
|
|
tp@0
|
355
|
|
tp@0
|
356 % Check if there are triplets
|
|
tp@0
|
357
|
|
tp@0
|
358 triptest = test(1:2:length(test))-test(2:2:length(test));
|
|
tp@0
|
359 triplets = find(triptest~=0);
|
|
tp@0
|
360
|
|
tp@0
|
361 if ~isempty(triplets)
|
|
tp@0
|
362 disp('ERROR: Triplets: some edges are defined three times which means that two thin planes have a correct')
|
|
tp@0
|
363 disp(' definition but some error on the rear side. You must correct these. ')
|
|
tp@0
|
364 disp(' Only the first one can be detected; it is:')
|
|
tp@0
|
365 disp(' (EDB1 corner numbers for the two edge points are given)')
|
|
tp@0
|
366 [floor(test(triplets(1)*2-1)/multfac) test(triplets(1)*2-1)-floor(test(triplets(1)*2-1)/multfac)*multfac]
|
|
tp@0
|
367 save ~/Documents/Temp/test.mat
|
|
tp@0
|
368 error
|
|
tp@0
|
369 end
|
|
tp@0
|
370
|
|
tp@0
|
371 edgecounter = length(test);
|
|
tp@0
|
372
|
|
tp@0
|
373 iv1 = find(planesatedge(:,1)~=0);
|
|
tp@0
|
374 iv2 = 1:edgecounter;
|
|
tp@0
|
375 iv2(iv1) = [];
|
|
tp@0
|
376
|
|
tp@0
|
377 edgecorners = [floor(test(iv1)/multfac) test(iv1)-floor(test(iv1)/multfac)*multfac];
|
|
tp@0
|
378 planesatedge = planesatedge(iv1,:) + planesatedge(iv2,:);
|
|
tp@0
|
379
|
|
tp@0
|
380 [nedges,slask] = size(edgecorners);
|
|
tp@0
|
381
|
|
tp@0
|
382 zerosvec3 = zeros(nedges,1);
|
|
tp@0
|
383 zerosvec5 = zeros(nedges,3);
|
|
tp@0
|
384
|
|
tp@0
|
385 if SHOWTEXT >= 4
|
|
tp@0
|
386 disp([' ',int2str(nedges),' edges found'])
|
|
tp@0
|
387 end
|
|
tp@0
|
388 thirdpoint = zerosvec5;
|
|
tp@0
|
389 for ii = 1:nedges
|
|
tp@0
|
390 refplane = planesatedge(ii,1);
|
|
tp@0
|
391 secondplane = planesatedge(ii,2);
|
|
tp@0
|
392 if secondplane ~= 0
|
|
tp@0
|
393 edgeco1 = corners(edgecorners(ii,1),:);
|
|
tp@0
|
394 edgeco2 = corners(edgecorners(ii,2),:);
|
|
tp@0
|
395 edgemidpoint = edgeco1 + (edgeco2-edgeco1)/2;
|
|
tp@0
|
396 intoplanevec = EDB1cross( (edgeco2-edgeco1).',(planenvecs(secondplane,:)).').';
|
|
tp@0
|
397 inplanepoint = edgemidpoint + intoplanevec*0.1;
|
|
tp@0
|
398 if sum(abs(inplanepoint)) == 0
|
|
tp@0
|
399 inplanepoint = edgemidpoint + intoplanevec*0.2;
|
|
tp@0
|
400 end
|
|
tp@0
|
401 thirdpoint(ii,:) = inplanepoint;
|
|
tp@0
|
402 end
|
|
tp@0
|
403 end
|
|
tp@0
|
404
|
|
tp@0
|
405
|
|
tp@0
|
406 %---------------------------------------------------------------
|
|
tp@0
|
407 % Calculate the closwedang values for all edges
|
|
tp@0
|
408
|
|
tp@0
|
409 if SHOWTEXT >= 4
|
|
tp@0
|
410 disp(' wedge angles...')
|
|
tp@0
|
411 end
|
|
tp@0
|
412 closwedangvec = zerosvec3;
|
|
tp@0
|
413
|
|
tp@0
|
414 ivec = find( planesatedge(:,1).*planesatedge(:,2) ~= 0);
|
|
tp@0
|
415 xwedgevec = [corners(edgecorners(ivec,1),:) corners(edgecorners(ivec,2),:)];
|
|
tp@0
|
416 nvec1vec = planenvecs( planesatedge(ivec,1),: );
|
|
tp@0
|
417 xsouvec = thirdpoint(ivec,:);
|
|
tp@0
|
418
|
|
tp@0
|
419 for jj = 1:length(ivec)
|
|
tp@0
|
420
|
|
tp@0
|
421 ii = ivec(jj);
|
|
tp@0
|
422 [slask,thetas,slask] = EDB1coordtrans1(xsouvec(jj,:),[xwedgevec(jj,1:3);xwedgevec(jj,4:6)],nvec1vec(jj,:));
|
|
tp@0
|
423 if thetas == 0
|
|
tp@0
|
424 closwedangvec(ii) = 0;
|
|
tp@0
|
425 else
|
|
tp@0
|
426 closwedangvec(ii) = 2*pi - thetas;
|
|
tp@0
|
427 end
|
|
tp@0
|
428
|
|
tp@0
|
429 end
|
|
tp@0
|
430
|
|
tp@0
|
431 %-------------------------------------------------------------------
|
|
tp@0
|
432 % Now we check for duplicates of the edge definitions
|
|
tp@0
|
433 % Some edge definitions might need to be swapped: if there are duplicate edges,
|
|
tp@0
|
434 % and one of them has closwedangvec = 0
|
|
tp@0
|
435 % then there is a mistake in the edge definition, so a swap will be made.
|
|
tp@0
|
436
|
|
tp@0
|
437 edgecosinglelist = edgecorners(:,1)*multfac + edgecorners(:,2);
|
|
tp@0
|
438 nsteps = diff([0;edgecosinglelist]);
|
|
tp@0
|
439 ivduplicates = find(nsteps == 0);
|
|
tp@0
|
440 nduplicates = length(ivduplicates);
|
|
tp@0
|
441 if nduplicates > 0
|
|
tp@0
|
442 if SHOWTEXT >= 4
|
|
tp@0
|
443 disp([' ',int2str(nduplicates),' edges formed by connected thin plates'])
|
|
tp@0
|
444 end
|
|
tp@0
|
445 for ii = 1:nduplicates
|
|
tp@0
|
446 comb1 = ivduplicates(ii)-1;
|
|
tp@0
|
447 comb2 = comb1+1;
|
|
tp@0
|
448 if closwedangvec(comb1) == 0 | closwedangvec(comb2) == 0, % edge definitions should be swapped
|
|
tp@0
|
449 plane1 = planesatedge(comb1,2);
|
|
tp@0
|
450 plane2 = planesatedge(comb2,2);
|
|
tp@0
|
451 planesatedge(comb1,2) = plane2;
|
|
tp@0
|
452 planesatedge(comb2,2) = plane1;
|
|
tp@0
|
453 temp1 = thirdpoint(comb1,:);
|
|
tp@0
|
454 temp2 = thirdpoint(comb2,:);
|
|
tp@0
|
455 thirdpoint(comb1,:) = temp2;
|
|
tp@0
|
456 thirdpoint(comb2,:) = temp1;
|
|
tp@0
|
457 if SHOWTEXT >= 4
|
|
tp@0
|
458 disp([' Swapping an edge definition'])
|
|
tp@0
|
459 disp([' Planes ',int2str(planesatedge(comb1,1)),' and ',int2str(planesatedge(comb1,2))])
|
|
tp@0
|
460 disp([' CATT: ',int2str(planenumbers(planesatedge(comb1,1))),' and ',int2str(planenumbers(planesatedge(comb1,2)))])
|
|
tp@0
|
461 disp([' Planes ',int2str(planesatedge(comb2,1)),' and ',int2str(planesatedge(comb2,2))])
|
|
tp@0
|
462 disp([' CATT: ',int2str(planenumbers(planesatedge(comb2,1))),' and ',int2str(planenumbers(planesatedge(comb2,2)))])
|
|
tp@0
|
463 end
|
|
tp@0
|
464
|
|
tp@0
|
465 xwedge = [corners(edgecorners(comb1,1),:);corners(edgecorners(comb1,2),:)];
|
|
tp@0
|
466 nvec1 = planenvecs( planesatedge(comb1,1),: );
|
|
tp@0
|
467 xsou = thirdpoint(comb1,:);
|
|
tp@0
|
468 [slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
|
|
tp@0
|
469 if thetas == 0
|
|
tp@0
|
470 closwedangvec(comb1) = 0;
|
|
tp@0
|
471 else
|
|
tp@0
|
472 closwedangvec(comb1) = 2*pi - thetas;
|
|
tp@0
|
473 end
|
|
tp@0
|
474
|
|
tp@0
|
475 xwedge = [corners(edgecorners(comb2,1),:);corners(edgecorners(comb2,2),:)];
|
|
tp@0
|
476 nvec1 = planenvecs( planesatedge(comb2,1),: );
|
|
tp@0
|
477 xsou = thirdpoint(comb2,:);
|
|
tp@0
|
478 [slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
|
|
tp@0
|
479 if thetas == 0
|
|
tp@0
|
480 closwedangvec(comb2) = 0;
|
|
tp@0
|
481 else
|
|
tp@0
|
482 closwedangvec(comb2) = 2*pi - thetas;
|
|
tp@0
|
483 end
|
|
tp@0
|
484
|
|
tp@0
|
485 end % if closwedangvec(comb1) == 0 |
|
|
tp@0
|
486
|
|
tp@0
|
487 end % for ii = 1:nduplicates
|
|
tp@0
|
488 end
|
|
tp@0
|
489
|
|
tp@0
|
490 if nplanes < 256
|
|
tp@0
|
491 planesatedge = uint8(planesatedge);
|
|
tp@0
|
492 elseif nplanes < 65535
|
|
tp@0
|
493 planesatedge = uint16(planesatedge);
|
|
tp@0
|
494 else
|
|
tp@0
|
495 planesatedge = uint32(planesatedge);
|
|
tp@0
|
496 end
|
|
tp@0
|
497
|
|
tp@0
|
498 %-------------------------------------------------------------------
|
|
tp@0
|
499 % Find which edges each plane is connected to
|
|
tp@0
|
500
|
|
tp@0
|
501 if SHOWTEXT >= 4
|
|
tp@0
|
502 disp(' find edgesatplane')
|
|
tp@0
|
503 end
|
|
tp@0
|
504
|
|
tp@0
|
505 edgesatplane = zeros(nplanes,double(max(ncornersperplanevec)));
|
|
tp@0
|
506
|
|
tp@0
|
507 for ii= 1:nplanes
|
|
tp@0
|
508 iv = find(planesatedge==ii);
|
|
tp@0
|
509 iv = iv - floor((iv-1)/nedges)*nedges;
|
|
tp@0
|
510 if ~isempty(iv)
|
|
tp@0
|
511 edgesatplane(ii,1:length(iv)) = iv.';
|
|
tp@0
|
512 end
|
|
tp@0
|
513 end
|
|
tp@0
|
514
|
|
tp@0
|
515 % Check how many edges are defined for each plane. There must be at
|
|
tp@0
|
516 % least three edges per plane.
|
|
tp@0
|
517
|
|
tp@0
|
518 tempnedgesperplanevec = sum(edgesatplane.'~=0).';
|
|
tp@0
|
519
|
|
tp@0
|
520 iv = find(tempnedgesperplanevec<3);
|
|
tp@0
|
521
|
|
tp@0
|
522 if ~isempty(iv)
|
|
tp@0
|
523
|
|
tp@0
|
524 disp(' ')
|
|
tp@0
|
525 if SHOWTEXT >= 4
|
|
tp@0
|
526 for ii = 1:length(iv)
|
|
tp@0
|
527 disp([' Plane ',int2str(iv(ii)),' has only ',int2str(tempnedgesperplanevec(iv(ii))),' edges connected to it'])
|
|
tp@0
|
528 end
|
|
tp@0
|
529 error(['ERROR: The planes above have less than three edges.'])
|
|
tp@0
|
530 else
|
|
tp@0
|
531 error(['ERROR: Some planes have less than three edges. Set SHOWTEXT >= 4 to see a list of those planes.'])
|
|
tp@0
|
532 end
|
|
tp@0
|
533
|
|
tp@0
|
534 end
|
|
tp@0
|
535
|
|
tp@0
|
536
|
|
tp@0
|
537
|
|
tp@0
|
538 %-------------------------------------------------------------------
|
|
tp@0
|
539 % Make a special list of thin planes
|
|
tp@0
|
540
|
|
tp@0
|
541 planeisthin = zerosvec1;
|
|
tp@0
|
542 rearsideplane = zerosvec1;
|
|
tp@0
|
543
|
|
tp@0
|
544 iv = find(closwedangvec==0);
|
|
tp@0
|
545 if ~isempty(iv)
|
|
tp@0
|
546 planenumberlist = planesatedge(iv,1);
|
|
tp@0
|
547 planeisthin(planenumberlist) = planeisthin(planenumberlist) + 1;
|
|
tp@0
|
548 rearsideplane(planenumberlist) = planesatedge(iv,2);
|
|
tp@0
|
549 end
|
|
tp@0
|
550
|
|
tp@0
|
551 iv = find(closwedangvec == 0 & planesatedge(:,2)~=0);
|
|
tp@0
|
552 if ~isempty(iv)
|
|
tp@0
|
553 planenumberlist = planesatedge(iv,2);
|
|
tp@0
|
554 planeisthin(planenumberlist) = planeisthin(planenumberlist) + 1;
|
|
tp@0
|
555 rearsideplane(planenumberlist) = planesatedge(iv,1);
|
|
tp@0
|
556 end
|
|
tp@0
|
557
|
|
tp@0
|
558 planeisthin = sign(planeisthin);
|
|
tp@0
|
559 listofthinplanes = find(planeisthin);
|
|
tp@0
|
560
|
|
tp@0
|
561 %---------------------------------------------------------------
|
|
tp@0
|
562 % Closwedangvec should be calculated into nyvec
|
|
tp@0
|
563
|
|
tp@0
|
564 nyvec = pi./(2*pi - closwedangvec);
|
|
tp@0
|
565 integerny = ( abs(nyvec - round(nyvec)) < 1e-10 );
|
|
tp@0
|
566
|
|
tp@0
|
567 %-----------------------------------------------------------
|
|
tp@0
|
568 % Construct some other useful edge variables
|
|
tp@0
|
569
|
|
tp@0
|
570 if difforder >= 1
|
|
tp@0
|
571 edgestartcoords = zerosvec5;
|
|
tp@0
|
572 edgestartcoordsnudge = zerosvec5;
|
|
tp@0
|
573 edgeendcoords = zerosvec5;
|
|
tp@0
|
574 edgeendcoordsnudge = zerosvec5;
|
|
tp@0
|
575 edgemidcoords = zerosvec5;
|
|
tp@0
|
576 edgenvecs = zerosvec5;
|
|
tp@0
|
577
|
|
tp@0
|
578 edgestartcoords = [corners(edgecorners(:,1),:)];
|
|
tp@0
|
579 edgeendcoords = [corners(edgecorners(:,2),:)];
|
|
tp@0
|
580 edgemidcoords = (edgestartcoords+edgeendcoords)/2;
|
|
tp@0
|
581 edgenvecs = planenvecs(planesatedge(:,1),:);
|
|
tp@0
|
582
|
|
tp@0
|
583 edgenormvecs = edgeendcoords - edgestartcoords;
|
|
tp@0
|
584
|
|
tp@0
|
585 edgestartcoordsnudge = edgestartcoords + 1e-10*edgenormvecs;
|
|
tp@0
|
586 edgeendcoordsnudge = edgeendcoords - 1e-10*edgenormvecs;
|
|
tp@0
|
587
|
|
tp@0
|
588 edgelengthvec = sqrt(sum( ((edgenormvecs).^2).' )).';
|
|
tp@0
|
589 edgenormvecs = edgenormvecs./edgelengthvec(:,ones(1,3));
|
|
tp@0
|
590
|
|
tp@0
|
591 else
|
|
tp@0
|
592 edgestartcoords = [];
|
|
tp@0
|
593 edgeendcoords = [];
|
|
tp@0
|
594 edgenvecs = [];
|
|
tp@0
|
595 edgelengthvec = [];
|
|
tp@0
|
596 end
|
|
tp@0
|
597
|
|
tp@0
|
598
|
|
tp@0
|
599 %---------------------------------------------------------------
|
|
tp@0
|
600 %
|
|
tp@0
|
601 % BACK TO SOME PURE PLANE RELATED PARAMETERS
|
|
tp@0
|
602 %
|
|
tp@0
|
603 %---------------------------------------------------------------
|
|
tp@0
|
604 %
|
|
tp@0
|
605 % First, make a list of which planes that are absorbing/non-reflective
|
|
tp@0
|
606
|
|
tp@0
|
607 planeisabsorbing = zerosvec1;
|
|
tp@0
|
608
|
|
tp@0
|
609 listofactiveplanes = find(reflfactors~=0);
|
|
tp@0
|
610 listofabsorbingplanes = find(reflfactors == 0);
|
|
tp@0
|
611 planeisabsorbing(listofabsorbingplanes) = ones(size(listofabsorbingplanes));
|
|
tp@0
|
612
|
|
tp@0
|
613 %---------------------------------------------------------------
|
|
tp@0
|
614 % Help variables: planealignedwplane and planeconnectstoplane
|
|
tp@0
|
615 %
|
|
tp@0
|
616 % Preparatory for the check of which planes see which plane: check
|
|
tp@0
|
617 % if any planes are aligned with each other. If two planes are aligned
|
|
tp@0
|
618 % then planealignedwplane = 1 for that combination.
|
|
tp@0
|
619 % However, if two planes are back-to back (thin planes), then
|
|
tp@0
|
620 % planealignedwplane = 2 for that combination.
|
|
tp@0
|
621 %
|
|
tp@0
|
622 % Also make a matrix of which planes connect to which planes and via
|
|
tp@0
|
623 % which edges.
|
|
tp@0
|
624
|
|
tp@0
|
625 planealignedwplane = zerosvec2;
|
|
tp@0
|
626
|
|
tp@0
|
627 % First, check which planes are aligned with each other
|
|
tp@0
|
628
|
|
tp@0
|
629 for ii = 1:nplanes-1
|
|
tp@0
|
630 if planeisabsorbing(ii) ~= 1
|
|
tp@0
|
631 oneplaneeq = planeeqs(ii,:);
|
|
tp@0
|
632 diffvec1 = oneplaneeq(ones(nplanes-ii,1),:) - planeeqs(ii+1:nplanes,:);
|
|
tp@0
|
633 diffvec2 = oneplaneeq(ones(nplanes-ii,1),:) + planeeqs(ii+1:nplanes,:);
|
|
tp@0
|
634 diffvec = min( [sum( diffvec1.'.^2 ).' sum( diffvec2.'.^2 ).'].' ).';
|
|
tp@0
|
635 alignedplanes = find(diffvec < geomacc) + ii;
|
|
tp@0
|
636 if ~isempty(alignedplanes)
|
|
tp@0
|
637 planealignedwplane(alignedplanes,ii) = int8(double(planealignedwplane(alignedplanes,ii)) + double((planeisabsorbing(alignedplanes)~=1)));
|
|
tp@0
|
638 end
|
|
tp@0
|
639 end
|
|
tp@0
|
640 end
|
|
tp@0
|
641
|
|
tp@0
|
642 planealignedwplane = (sparse(sign(double(planealignedwplane) + double(planealignedwplane).')));
|
|
tp@0
|
643
|
|
tp@0
|
644 % Second, check which planes are back to back
|
|
tp@0
|
645
|
|
tp@0
|
646 if ~isempty(listofthinplanes)
|
|
tp@0
|
647 for ii = 1:length(listofthinplanes)
|
|
tp@0
|
648 plane = listofthinplanes(ii);
|
|
tp@0
|
649 rearplane = rearsideplane(plane);
|
|
tp@0
|
650 if planeisabsorbing(plane) ~= 1 & planeisabsorbing(rearplane) ~= 1
|
|
tp@0
|
651 planealignedwplane(plane,rearplane) = 2;
|
|
tp@0
|
652 planealignedwplane(rearplane,plane) = 2;
|
|
tp@0
|
653 end
|
|
tp@0
|
654 end
|
|
tp@0
|
655 end
|
|
tp@0
|
656
|
|
tp@0
|
657 planeconnectstoplane = zerosvec2;
|
|
tp@0
|
658 clear zerosvec2
|
|
tp@0
|
659
|
|
tp@0
|
660 for ii = 1:nplanes
|
|
tp@0
|
661 if planeisabsorbing(ii) ~= 1
|
|
tp@0
|
662 edgelist = edgesatplane(ii,:);
|
|
tp@0
|
663 edgelist = edgesatplane(ii,1:double(ncornersperplanevec(ii)));
|
|
tp@0
|
664 ivec = planesatedge(edgelist,:);
|
|
tp@0
|
665 ivec = reshape(ivec.',length(edgelist)*2,1);
|
|
tp@0
|
666 ivec = ivec(find(ivec~=ii));
|
|
tp@0
|
667 okplanes = find(planeisabsorbing(ivec)~=1);
|
|
tp@0
|
668 ivec = ivec(okplanes);
|
|
tp@0
|
669 if ~isempty(ivec)
|
|
tp@0
|
670 planeconnectstoplane(ii,ivec) = edgelist(okplanes);
|
|
tp@0
|
671 end
|
|
tp@0
|
672 end
|
|
tp@0
|
673 end
|
|
tp@0
|
674
|
|
tp@0
|
675 %---------------------------------------------------------------
|
|
tp@0
|
676 %
|
|
tp@0
|
677 % planeseesplane
|
|
tp@0
|
678 %
|
|
tp@0
|
679 %---------------------------------------------------------------
|
|
tp@0
|
680 % Check which planes see which planes:
|
|
tp@0
|
681 % 1. If one of the planes has reflfactors = 0, then the two can't see each other.
|
|
tp@0
|
682 % 2. Reflective planes that are aligned with each other can not reach each other
|
|
tp@0
|
683 % 3. Reflective planes that are not aligned but have the same normal vectors
|
|
tp@0
|
684 % can not reach each other
|
|
tp@0
|
685 % 4. Planes that have all its corners behind another plane can not be seen by that plane.
|
|
tp@0
|
686 %
|
|
tp@0
|
687 % planeseesplane = 0 means that either, one of the planes is non-reflective or, that
|
|
tp@0
|
688 % two planes never can see each other, but they are not aligned with each other.
|
|
tp@0
|
689 % planeseesplane = 1 means that two reflective planes might see each other, but there could
|
|
tp@0
|
690 % be obstructions
|
|
tp@0
|
691 % planeseesplane = -1 means that two reflective planes are aligned (but they are not back-to-back)
|
|
tp@0
|
692 % and thus they can never see each other.
|
|
tp@0
|
693 % planeseesplane = -2 means that two reflective planes are back-to-back with each other and thus
|
|
tp@0
|
694 % they can never see each other.
|
|
tp@0
|
695
|
|
tp@0
|
696 if SHOWTEXT >= 4
|
|
tp@0
|
697 disp(' Check which planes see which planes')
|
|
tp@0
|
698 end
|
|
tp@0
|
699
|
|
tp@0
|
700 % First an auxilliary matrix which can be used for edgeseesplane later:
|
|
tp@0
|
701 % cornerinfrontofplane, size [nplanes,ncorners]
|
|
tp@0
|
702 % Values are the same as given by EDB1infrontofplane:
|
|
tp@0
|
703 % 1 means that a point is in front of the plane
|
|
tp@0
|
704 % 0 means that a point is aligned with a plane
|
|
tp@0
|
705 % -1 means that a point is behind a plane
|
|
tp@0
|
706 % All corners that belong to a plane will have the value 0.
|
|
tp@0
|
707
|
|
tp@0
|
708 % Corner number is given by the col no.
|
|
tp@0
|
709 if ncorners < 256
|
|
tp@0
|
710 cornernumb = uint8([1:ncorners]);
|
|
tp@0
|
711 elseif ncorners < 65536
|
|
tp@0
|
712 cornernumb = uint16([1:ncorners]);
|
|
tp@0
|
713 else
|
|
tp@0
|
714 cornernumb = uint32([1:ncorners]);
|
|
tp@0
|
715 end
|
|
tp@0
|
716
|
|
tp@0
|
717 % Plane numbers is given by the row no.
|
|
tp@0
|
718 if nplanes < 256
|
|
tp@0
|
719 planenumb = uint8([1:nplanes].');
|
|
tp@0
|
720 elseif nplanes < 65536
|
|
tp@0
|
721 planenumb = uint16([1:nplanes].');
|
|
tp@0
|
722 else
|
|
tp@0
|
723 planenumb = uint32([1:nplanes].');
|
|
tp@0
|
724 end
|
|
tp@0
|
725
|
|
tp@0
|
726 cornerinfrontofplane = int8(EDB1infrontofplane(corners,planenvecs,planecorners,[],cornernumb,planenumb));
|
|
tp@0
|
727 clear planenumb cornernumb
|
|
tp@0
|
728 cornerinfrontofplane = reshape(cornerinfrontofplane,nplanes,ncorners);
|
|
tp@0
|
729
|
|
tp@0
|
730 planeseesplane = int8(ones(nplanes,nplanes));
|
|
tp@0
|
731
|
|
tp@0
|
732 % 1. If one of the planes has reflfactors = 0, then the two can't see each other.
|
|
tp@0
|
733
|
|
tp@0
|
734 if ~isempty(listofabsorbingplanes)
|
|
tp@0
|
735 zerosvec = zeros(length(listofabsorbingplanes),nplanes);
|
|
tp@0
|
736 planeseesplane(listofabsorbingplanes,:) = zerosvec;
|
|
tp@0
|
737 planeseesplane(:,listofabsorbingplanes) = zerosvec.';
|
|
tp@0
|
738 end
|
|
tp@0
|
739
|
|
tp@0
|
740 % 2. Reflective planes that are aligned with each other can not reach each other
|
|
tp@0
|
741
|
|
tp@0
|
742 ivec = find( planealignedwplane == 1);
|
|
tp@0
|
743 if ~isempty(ivec)
|
|
tp@0
|
744 planeseesplane(ivec) = - 1*ones(size(ivec));
|
|
tp@0
|
745 end
|
|
tp@0
|
746
|
|
tp@0
|
747 ivec = find( planealignedwplane == 2);
|
|
tp@0
|
748 if ~isempty(ivec)
|
|
tp@0
|
749 planeseesplane(ivec) = - 2*ones(size(ivec));
|
|
tp@0
|
750 end
|
|
tp@0
|
751
|
|
tp@0
|
752 % 3. Reflective planes that have the same normal vectors can not reach each other
|
|
tp@0
|
753
|
|
tp@0
|
754 numvec = [1:nplanes];
|
|
tp@0
|
755 for ii = 1:nplanes
|
|
tp@0
|
756 if planeisabsorbing(ii) ~= 1
|
|
tp@0
|
757 nvec1 = planenvecs(ii,:);
|
|
tp@0
|
758 ivec = find(planealignedwplane(ii,:)==0 & planeisabsorbing.'==0 & numvec>ii);
|
|
tp@0
|
759 if ~isempty(ivec)
|
|
tp@0
|
760 similarvec = abs( nvec1(ones(length(ivec),1),:) - planenvecs(ivec,:)) < geomacc;
|
|
tp@0
|
761 similarvec = prod( double(similarvec.') ).';
|
|
tp@0
|
762 ivec2 = ivec(find(similarvec==1));
|
|
tp@0
|
763 if ~isempty(ivec2)
|
|
tp@0
|
764 zerosvec = zeros(size(ivec2));
|
|
tp@0
|
765 planeseesplane(ii,ivec2) = zerosvec;
|
|
tp@0
|
766 planeseesplane(ivec2,ii) = zerosvec.';
|
|
tp@0
|
767 end
|
|
tp@0
|
768 end
|
|
tp@0
|
769 end
|
|
tp@0
|
770 end
|
|
tp@0
|
771
|
|
tp@0
|
772 % 3.5 A plane does not see itself
|
|
tp@0
|
773
|
|
tp@0
|
774 iv = [0:nplanes-1]*nplanes + [1:nplanes];
|
|
tp@0
|
775 planeseesplane(iv) = zeros(size(iv));
|
|
tp@0
|
776
|
|
tp@0
|
777 % 4. Planes that have all its corners behind another plane can not be seen by that plane.
|
|
tp@0
|
778 % First, we construct a list referring to the complete matrix of size [nplanes,nplanes]
|
|
tp@0
|
779 % The index vector is called iv. Then we find which combinations that
|
|
tp@0
|
780 % could be cleared. After having cleared a number of combinations (in
|
|
tp@0
|
781 % iv, fromplanenumb and toplanenumb) we pick out out only the non-zero
|
|
tp@0
|
782 % indices in iv.
|
|
tp@0
|
783
|
|
tp@0
|
784 if nplanes*nplanes < 128
|
|
tp@0
|
785 iv = int8([1:nplanes*nplanes].');
|
|
tp@0
|
786 elseif nplanes*nplanes < 32768
|
|
tp@0
|
787 iv = int16([1:nplanes*nplanes].');
|
|
tp@0
|
788 else
|
|
tp@0
|
789 iv = int32([1:nplanes*nplanes].');
|
|
tp@0
|
790 end
|
|
tp@0
|
791 iv = reshape(iv,nplanes,nplanes);
|
|
tp@0
|
792
|
|
tp@0
|
793 % If there are any absorbing planes, we zero all plane-to-plane
|
|
tp@0
|
794 % combinations that involve such a plane
|
|
tp@0
|
795
|
|
tp@0
|
796 if ~isempty(listofabsorbingplanes)
|
|
tp@0
|
797 nabsorbingplanes = length(listofabsorbingplanes);
|
|
tp@0
|
798 iv(listofabsorbingplanes,:) = zeros(nabsorbingplanes,nplanes);
|
|
tp@0
|
799 iv(:,listofabsorbingplanes) = zeros(nplanes,nabsorbingplanes);
|
|
tp@0
|
800 end
|
|
tp@0
|
801
|
|
tp@0
|
802 % Check connecting planes first: if two planes are connected and their
|
|
tp@0
|
803 % shared edge has a closwedangvec > pi, then the two planes can see each
|
|
tp@0
|
804 % other. So, we can zero the plane-pair combinations where closwedangvec
|
|
tp@0
|
805 % <= pi but larger than zero.
|
|
tp@0
|
806
|
|
tp@0
|
807 edgecombstozero = find(closwedangvec<=pi & closwedangvec > 0);
|
|
tp@0
|
808 if ~isempty(edgecombstozero)
|
|
tp@0
|
809 ivzerocombs = (double(planesatedge(edgecombstozero,1))-1)*nplanes + double(planesatedge(edgecombstozero,2));
|
|
tp@0
|
810 iv(ivzerocombs) = zeros(size(ivzerocombs));
|
|
tp@0
|
811 ivzerocombs = (double(planesatedge(edgecombstozero,2))-1)*nplanes + double(planesatedge(edgecombstozero,1));
|
|
tp@0
|
812 iv(ivzerocombs) = zeros(size(ivzerocombs));
|
|
tp@0
|
813 end
|
|
tp@0
|
814
|
|
tp@0
|
815 % Now we should keep only the index numbers in iv that are still non-zero
|
|
tp@0
|
816 % We also take the chance to zero planeseesplane for the combinations that
|
|
tp@0
|
817 % could be cleared.
|
|
tp@0
|
818 ivclear = uint32(find(iv==0));
|
|
tp@0
|
819 planeseesplane(ivclear) = zeros(size(ivclear));
|
|
tp@0
|
820 iv(ivclear) = [];
|
|
tp@0
|
821
|
|
tp@0
|
822 % Of the remaining plane-to-plane combinations, we pick out the ones
|
|
tp@0
|
823 % for which we need to check if the "to-plane" is in front of the "from-plane"
|
|
tp@0
|
824 iv = iv(find(planeseesplane(iv)==1 & planeconnectstoplane(iv)==0));
|
|
tp@0
|
825
|
|
tp@0
|
826 % We create full matrices, fromplanenumb and toplanenumb, and then keep
|
|
tp@0
|
827 % only the combinations that remain in iv
|
|
tp@0
|
828 % to keep fromplanenumb and toplanenumb aligned with iv.
|
|
tp@0
|
829
|
|
tp@0
|
830 % To-plane numbers is given by the row no.
|
|
tp@0
|
831 if nplanes < 256
|
|
tp@0
|
832 toplanenumb = uint8([1:nplanes].');
|
|
tp@0
|
833 elseif nplanes < 65536
|
|
tp@0
|
834 toplanenumb = uint16([1:nplanes].');
|
|
tp@0
|
835 else
|
|
tp@0
|
836 toplanenumb = uint32([1:nplanes].');
|
|
tp@0
|
837 end
|
|
tp@0
|
838 toplanenumb = reshape(toplanenumb(:,ones(1,nplanes)),nplanes*nplanes,1);
|
|
tp@0
|
839 toplanenumb = toplanenumb(iv);
|
|
tp@0
|
840
|
|
tp@0
|
841 % From-plane numbers is given by the col no.
|
|
tp@0
|
842 if nplanes < 256
|
|
tp@0
|
843 fromplanenumb = uint8([1:nplanes]);
|
|
tp@0
|
844 elseif nplanes < 65536
|
|
tp@0
|
845 fromplanenumb = uint16([1:nplanes]);
|
|
tp@0
|
846 else
|
|
tp@0
|
847 fromplanenumb = uint32([1:nplanes]);
|
|
tp@0
|
848 end
|
|
tp@0
|
849 fromplanenumb = reshape(fromplanenumb(ones(nplanes,1),:),nplanes*nplanes,1);
|
|
tp@0
|
850 fromplanenumb = fromplanenumb(iv);
|
|
tp@0
|
851
|
|
tp@0
|
852 % Now, if a plane has *all* its corners behind another plane, those two
|
|
tp@0
|
853 % planes can not see each other.
|
|
tp@0
|
854 % We check the corners of the "to-plane" (columns) and check if they are in front of
|
|
tp@0
|
855 % the "from-plane" (rows).
|
|
tp@0
|
856 % The ivreftocoplamatrix is the index number to the
|
|
tp@0
|
857 % cornerinfrontofplanematrix, which has the size [nplanes,ncorners].
|
|
tp@0
|
858 % NB! In order to save some space, we don't construct the ivreftocoplamatrix specifically.
|
|
tp@0
|
859
|
|
tp@0
|
860 % First we check the first 3 corners since all planes have at least 3
|
|
tp@0
|
861 % corners.
|
|
tp@0
|
862
|
|
tp@0
|
863 % Plane corner 1
|
|
tp@0
|
864 cornersbehind = cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,1 ))-1 )*nplanes)<1;
|
|
tp@0
|
865
|
|
tp@0
|
866 % Plane corner 2
|
|
tp@0
|
867 cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,2 ))-1 )*nplanes) < 1);
|
|
tp@0
|
868
|
|
tp@0
|
869 % Plane corner 3
|
|
tp@0
|
870 cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,3 ))-1 )*nplanes) < 1);
|
|
tp@0
|
871
|
|
tp@0
|
872 if minncornersperplane == 4 & maxncornersperplane == 4
|
|
tp@0
|
873 % Plane corner 4
|
|
tp@0
|
874 cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,4 ))-1 )*nplanes) < 1);
|
|
tp@0
|
875 elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
|
|
tp@0
|
876 ivall = uint32([1:length(fromplanenumb)].');
|
|
tp@0
|
877 if minncornersperplane == 3
|
|
tp@0
|
878 iv3cornerplanes = find(ncornersperplanevec(toplanenumb)==3);
|
|
tp@0
|
879 ivall(iv3cornerplanes) = [];
|
|
tp@0
|
880 end
|
|
tp@0
|
881 if maxncornersperplane == 4
|
|
tp@0
|
882 % Plane corner 4
|
|
tp@0
|
883 cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
|
|
tp@0
|
884 else
|
|
tp@0
|
885 ivmorethan4 = find(ncornersperplanevec(toplanenumb(ivall))>4);
|
|
tp@0
|
886 temp = ivall(ivmorethan4);
|
|
tp@0
|
887 ivall(ivmorethan4) = [];
|
|
tp@0
|
888 ivmorethan4 = temp;
|
|
tp@0
|
889 cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
|
|
tp@0
|
890 for ii = 5:maxncornersperplane
|
|
tp@0
|
891 ivselection = find(ncornersperplanevec(toplanenumb(ivmorethan4))>=ii);
|
|
tp@0
|
892 cornersbehind(ivmorethan4(ivselection)) = cornersbehind(ivmorethan4(ivselection)) &(cornerinfrontofplane(double(fromplanenumb(ivmorethan4(ivselection))) + ( double(planecorners( toplanenumb(ivmorethan4(ivselection)),ii ))-1 )*nplanes) < 1);
|
|
tp@0
|
893 end
|
|
tp@0
|
894 clear ivselection
|
|
tp@0
|
895 end
|
|
tp@0
|
896 end
|
|
tp@0
|
897 clear toplanenumb fromplanenumb
|
|
tp@0
|
898
|
|
tp@0
|
899 ivclear = iv(find(cornersbehind==1));
|
|
tp@0
|
900 clear iv
|
|
tp@0
|
901
|
|
tp@0
|
902 planeseesplane(ivclear) = zeros(size(ivclear));
|
|
tp@0
|
903 clear ivclear
|
|
tp@0
|
904
|
|
tp@0
|
905 temp1 = (planeseesplane~=0);
|
|
tp@0
|
906 temp2 = (planeseesplane.'~=0);
|
|
tp@0
|
907 planeseesplane = int8(double(planeseesplane).*(temp1.*temp2));
|
|
tp@0
|
908
|
|
tp@0
|
909 % New addition:
|
|
tp@0
|
910 % If the user has asked for it, check plane-mid-point to plane-mid-point
|
|
tp@0
|
911 % for obstruction. If there are obstructions, set planeseesplane to 0 for
|
|
tp@0
|
912 % that combination.
|
|
tp@0
|
913
|
|
tp@0
|
914 if planeseesplanestrategy == 1
|
|
tp@0
|
915 ivorig = uint32(find(planeseesplane==1));
|
|
tp@0
|
916 iv = ivorig;
|
|
tp@0
|
917 fromplane = ceil(double(iv)/nplanes);
|
|
tp@0
|
918 toplane = uint16(double(iv) - (fromplane-1)*nplanes);
|
|
tp@0
|
919 fromplane = uint16(fromplane);
|
|
tp@0
|
920 listofplanesthatcanobscure = find( sum(planeseesplane==1) );
|
|
tp@0
|
921
|
|
tp@0
|
922 iv4planes = find(ncornersperplanevec == 4);
|
|
tp@0
|
923 planemidpoints = zeros(nplanes,3);
|
|
tp@0
|
924
|
|
tp@0
|
925 if any(ncornersperplanevec~=4)
|
|
tp@0
|
926 disp('WARNING: planeseesplanestrategy 1 not implemented for models with non-4-corner planes')
|
|
tp@0
|
927 planeseesplanestrategy = 0;
|
|
tp@0
|
928 else
|
|
tp@0
|
929 xvalues = corners(planecorners.',1);
|
|
tp@0
|
930 xvalues = reshape(xvalues,4,length(xvalues)/4);
|
|
tp@0
|
931 yvalues = corners(planecorners.',2);
|
|
tp@0
|
932 yvalues = reshape(yvalues,4,length(yvalues)/4);
|
|
tp@0
|
933 zvalues = corners(planecorners.',3);
|
|
tp@0
|
934 zvalues = reshape(zvalues,4,length(zvalues)/4);
|
|
tp@0
|
935 planemidpoints = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
|
|
tp@0
|
936 end
|
|
tp@0
|
937 end
|
|
tp@0
|
938
|
|
tp@0
|
939 if planeseesplanestrategy == 1
|
|
tp@0
|
940 for ii = listofplanesthatcanobscure
|
|
tp@0
|
941 planeobsc = ii;
|
|
tp@0
|
942 ivcheckvis1 = uint32( planeobsc + (double(fromplane)-1)*nplanes);
|
|
tp@0
|
943 ivcheckvis2 = uint32( planeobsc + (double(toplane)-1)*nplanes);
|
|
tp@0
|
944 iv3 = find(toplane~=planeobsc & fromplane~=planeobsc & ((planeseesplane(ivcheckvis1) == 1) | (planeseesplane(ivcheckvis2) == 1)) & (planeseesplane(ivcheckvis1) >= 0) & (planeseesplane(ivcheckvis2) >= 0) );
|
|
tp@0
|
945 tempcanplaneobstruct = zerosvec1.';
|
|
tp@0
|
946 tempcanplaneobstruct(planeobsc) = 1;
|
|
tp@0
|
947 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(planemidpoints(fromplane(iv3),:),planemidpoints(toplane(iv3),:),[],[],tempcanplaneobstruct,planeseesplane,...
|
|
tp@0
|
948 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
|
|
tp@0
|
949 if length(iv3) > length(nonobstructedpaths)
|
|
tp@0
|
950 iv3(nonobstructedpaths) = [];
|
|
tp@0
|
951 iv(iv3) = [];
|
|
tp@0
|
952 toplane(iv3) = [];
|
|
tp@0
|
953 fromplane(iv3) = [];
|
|
tp@0
|
954 end
|
|
tp@0
|
955 end
|
|
tp@0
|
956 planeseesplane(ivorig) = 0;
|
|
tp@0
|
957 planeseesplane(iv) = 1;
|
|
tp@0
|
958 clear iv ivorig toplane fromplane ivcheckvis ivcheckvis2
|
|
tp@0
|
959 end
|
|
tp@0
|
960
|
|
tp@0
|
961 %---------------------------------------------------------------
|
|
tp@0
|
962 %
|
|
tp@0
|
963 % EDGE RELATED PARAMETERS AGAIN
|
|
tp@0
|
964 %
|
|
tp@0
|
965 %---------------------------------------------------------------
|
|
tp@0
|
966 %
|
|
tp@0
|
967 % Which edges should be switched off?
|
|
tp@0
|
968 % (1) Go through the list edgecorners. If any edge contains one of the
|
|
tp@0
|
969 % corners, whose numbers are in switchoffcorners, then that edge should be
|
|
tp@0
|
970 % stored in the list offedges.
|
|
tp@0
|
971 % (2) If an edge has an integer ny-value, then it should be switched off.
|
|
tp@0
|
972 % (3) If one of the two connecting planes has reflfactors = 0, then it should
|
|
tp@0
|
973 % be switched off.
|
|
tp@0
|
974
|
|
tp@0
|
975 % disp(' switching off possibly unwanted edges...')
|
|
tp@0
|
976
|
|
tp@0
|
977 offedges = zerosvec3;
|
|
tp@0
|
978
|
|
tp@0
|
979 % (1) Go through the list edgecorners. If any edge contains one of the
|
|
tp@0
|
980 % corners, whose numbers are in switchoffcorners, then that edge should be
|
|
tp@0
|
981 % stored in the list offedges.
|
|
tp@0
|
982
|
|
tp@0
|
983 if ~isempty(switchoffcorners)
|
|
tp@0
|
984 for ii = 1:nedges
|
|
tp@0
|
985 corner1 = edgecorners(ii,1);
|
|
tp@0
|
986 corner2 = edgecorners(ii,2);
|
|
tp@0
|
987 remove = sum( corner1 == switchoffcorners ) + sum( corner2 == switchoffcorners );
|
|
tp@0
|
988 offedges(ii) = sign(remove);
|
|
tp@0
|
989 end
|
|
tp@0
|
990 end
|
|
tp@0
|
991
|
|
tp@0
|
992 % (2) If an edge has an integer ny-value, then it should be switched off.
|
|
tp@0
|
993
|
|
tp@0
|
994 offedges = offedges + integerny;
|
|
tp@0
|
995
|
|
tp@0
|
996 % (3) If one of the two connecting planes has reflfactors = 0, then it should
|
|
tp@0
|
997 % be switched off.
|
|
tp@0
|
998
|
|
tp@0
|
999 reflectingedges = reshape(reflfactors(planesatedge),nedges,2);
|
|
tp@0
|
1000 reflectingedges = reflectingedges(:,1).*reflectingedges(:,2);
|
|
tp@0
|
1001 offedges = offedges + (1-reflectingedges);
|
|
tp@0
|
1002
|
|
tp@0
|
1003 offedges = find(offedges ~= 0);
|
|
tp@0
|
1004
|
|
tp@0
|
1005 %--------------------------------------------------------------------------------
|
|
tp@0
|
1006 %
|
|
tp@0
|
1007 % edgeseesplane
|
|
tp@0
|
1008 %
|
|
tp@0
|
1009 %--------------------------------------------------------------------------------
|
|
tp@0
|
1010 %
|
|
tp@0
|
1011 % Now check which planes each edge can see.
|
|
tp@0
|
1012 % 1. Set all combinations with reflfactors = 0 to -3
|
|
tp@0
|
1013 % 2. Switch off all combinations with offedges
|
|
tp@0
|
1014 % 3. For all edges that are aligned with a plane, set edgeseesplane to:
|
|
tp@0
|
1015 % -1 if the edge doesn't belong to the plane
|
|
tp@0
|
1016 % -2 if the edge does belong to the plane
|
|
tp@0
|
1017 % 4. If at least one edge end point is in front of a plane, then the
|
|
tp@0
|
1018 % plane is potentially visible (edgeseesplane = 1)
|
|
tp@0
|
1019 % 5. If at least one corner point is in front of:
|
|
tp@0
|
1020 % *one* of the two edge planes, when closwedang < pi
|
|
tp@0
|
1021 % *both* edge planes, when closwedangvec > pi
|
|
tp@0
|
1022
|
|
tp@0
|
1023 if SHOWTEXT >= 4
|
|
tp@0
|
1024 disp(' Check which edges see which planes')
|
|
tp@0
|
1025 end
|
|
tp@0
|
1026
|
|
tp@0
|
1027 edgeseesplane = int8(ones(nplanes,nedges));
|
|
tp@0
|
1028
|
|
tp@0
|
1029 % % % 20120414 PS: Why has this first stage been shut off? CHECK
|
|
tp@0
|
1030 % 1. Set all edges belonging to planes with reflfactors = 0 to
|
|
tp@0
|
1031 % edgeseesplane = -3.
|
|
tp@0
|
1032
|
|
tp@0
|
1033 % % if ~isempty(listofabsorbingplanes)
|
|
tp@0
|
1034 % % edgeseesplane(listofabsorbingplanes,:) = -3*ones(length(listofabsorbingplanes),nedges);
|
|
tp@0
|
1035 % % end
|
|
tp@0
|
1036
|
|
tp@0
|
1037 % 2. Switch off all combinations with offedges
|
|
tp@0
|
1038
|
|
tp@0
|
1039 if ~isempty(offedges)
|
|
tp@0
|
1040 edgeseesplane(:,offedges) = zeros(nplanes,length(offedges));
|
|
tp@0
|
1041 end
|
|
tp@0
|
1042
|
|
tp@0
|
1043 % 3. For all edges that belong to a plane, set edgeseesplane to -2
|
|
tp@0
|
1044 % Also, for all other edges that are aligned with a plane, set
|
|
tp@0
|
1045 % edgeseesplane to -1
|
|
tp@0
|
1046
|
|
tp@0
|
1047 edgelist = [1:nedges].';
|
|
tp@0
|
1048 edgelist(offedges) = [];
|
|
tp@0
|
1049
|
|
tp@0
|
1050 plane1list = planesatedge(edgelist,1);
|
|
tp@0
|
1051 plane2list = planesatedge(edgelist,2);
|
|
tp@0
|
1052 indexvec = uint32( double(plane1list) + (edgelist-1)*nplanes);
|
|
tp@0
|
1053 edgeseesplane(indexvec) = -2*(reflfactors(plane1list)~=0);
|
|
tp@0
|
1054 indexvec = uint32( double(plane2list) + (edgelist-1)*nplanes);
|
|
tp@0
|
1055 edgeseesplane(indexvec) = -2*(reflfactors(plane2list)~=0);
|
|
tp@0
|
1056
|
|
tp@0
|
1057 for ii = 1:length(edgelist)
|
|
tp@0
|
1058
|
|
tp@0
|
1059 aligners1 = find(planealignedwplane(:,plane1list(ii))==1);
|
|
tp@0
|
1060 if ~isempty(aligners1)
|
|
tp@0
|
1061 edgeseesplane(aligners1,edgelist(ii)) = -1*ones(size(aligners1));
|
|
tp@0
|
1062 end
|
|
tp@0
|
1063
|
|
tp@0
|
1064 aligners2 = find(planealignedwplane(:,plane2list(ii))==1);
|
|
tp@0
|
1065 if ~isempty(aligners2)
|
|
tp@0
|
1066 edgeseesplane(aligners2,edgelist(ii)) = -1*ones(size(aligners2));
|
|
tp@0
|
1067 end
|
|
tp@0
|
1068
|
|
tp@0
|
1069 end
|
|
tp@0
|
1070
|
|
tp@0
|
1071 % 4. If at least one edge end point is in front of a plane, then the
|
|
tp@0
|
1072 % plane is potentially visible.
|
|
tp@0
|
1073 % Do this for all plane-edge combinations for which edgeseesplane = 1
|
|
tp@0
|
1074 % up til now.
|
|
tp@0
|
1075 % 5. If at least one corner point is in front of:
|
|
tp@0
|
1076 % *one* of the two edge planes, when closwedang < pi
|
|
tp@0
|
1077 % *both* edge planes, when closwedangvec > pi
|
|
tp@0
|
1078
|
|
tp@0
|
1079 ntot = nplanes*nedges;
|
|
tp@0
|
1080 ivclear = find(edgeseesplane~=1);
|
|
tp@0
|
1081 % Edge number is given by the col no.
|
|
tp@0
|
1082 if nedges < 256
|
|
tp@0
|
1083 edgenumb = uint8([1:nedges]);
|
|
tp@0
|
1084 elseif nedges < 65536
|
|
tp@0
|
1085 edgenumb = uint16([1:nedges]);
|
|
tp@0
|
1086 else
|
|
tp@0
|
1087 edgenumb = uint32([1:nedges]);
|
|
tp@0
|
1088 end
|
|
tp@0
|
1089 edgenumb = reshape(edgenumb(ones(nplanes,1),:),nplanes*nedges,1);
|
|
tp@0
|
1090 edgenumb(ivclear) = [];
|
|
tp@0
|
1091 % Plane numbers is given by the row no.
|
|
tp@0
|
1092 if nplanes < 256
|
|
tp@0
|
1093 planenumb = uint8([1:nplanes].');
|
|
tp@0
|
1094 elseif nplanes < 65536
|
|
tp@0
|
1095 planenumb = uint16([1:nplanes].');
|
|
tp@0
|
1096 else
|
|
tp@0
|
1097 planenumb = uint32([1:nplanes].');
|
|
tp@0
|
1098 end
|
|
tp@0
|
1099 planenumb = reshape(planenumb(:,ones(1,nedges)),nplanes*nedges,1);
|
|
tp@0
|
1100 planenumb(ivclear) = [];
|
|
tp@0
|
1101
|
|
tp@0
|
1102 if nplanes*nedges < 128
|
|
tp@0
|
1103 iv = int8([1:nplanes*nedges].');
|
|
tp@0
|
1104 elseif nplanes*nedges < 32768
|
|
tp@0
|
1105 iv = int16([1:nplanes*nedges].');
|
|
tp@0
|
1106 else
|
|
tp@0
|
1107 iv = int32([1:nplanes*nedges].');
|
|
tp@0
|
1108 end
|
|
tp@0
|
1109 iv(ivclear) = [];
|
|
tp@0
|
1110
|
|
tp@0
|
1111 if ~isempty(edgenumb)
|
|
tp@0
|
1112 % For the "edgecorners in front of plane"-test
|
|
tp@0
|
1113 % we can use the cornerinfrontofplane matrix, but we need to construct
|
|
tp@0
|
1114 % index vectors, based on the planenumb and edgenumb values.
|
|
tp@0
|
1115
|
|
tp@0
|
1116 edgeinfrontofplane = cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,1) )-1 )*nplanes );
|
|
tp@0
|
1117 edgeinfrontofplane = (edgeinfrontofplane==1) | (cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,2) )-1 )*nplanes )==1);
|
|
tp@0
|
1118 edgeseesplane(iv) = int8(double(edgeseesplane(iv)).*double(edgeinfrontofplane));
|
|
tp@0
|
1119 clear edgeinfrontofplane
|
|
tp@0
|
1120 end
|
|
tp@0
|
1121
|
|
tp@0
|
1122 if ~isempty(edgenumb)
|
|
tp@0
|
1123 % For the "planecorners in front of edge-plane"-test
|
|
tp@0
|
1124 % we can use the cornerinfrontofplane matrix, but we need to construct
|
|
tp@0
|
1125 % index vectors, based on the planenumb and edgenumb values.
|
|
tp@0
|
1126 %
|
|
tp@0
|
1127 % First we split up the iv into two halves: the one with edges < pi and the
|
|
tp@0
|
1128 % one with edges > pi;
|
|
tp@0
|
1129
|
|
tp@0
|
1130 ivsmaller = uint32(find(closwedangvec(edgenumb)<pi));
|
|
tp@0
|
1131 ivlarger = uint32([1:length(edgenumb)].');
|
|
tp@0
|
1132 ivlarger(ivsmaller) = [];
|
|
tp@0
|
1133
|
|
tp@0
|
1134 if ~isempty(ivsmaller)
|
|
tp@0
|
1135 % First the edges smaller than pi (closwedang < pi): at least one plane corner needs
|
|
tp@0
|
1136 % to be in front of one or the other plane.
|
|
tp@0
|
1137 % Plane corner 1
|
|
tp@0
|
1138 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),1) ) + ( double(planecorners( planenumb(ivsmaller),1 ))-1 )*nplanes;
|
|
tp@0
|
1139 planecornerinfrontofedge = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1140 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),2) ) + ( double(planecorners( planenumb(ivsmaller),1 ))-1 )*nplanes;
|
|
tp@0
|
1141 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1142 % Plane corner 2
|
|
tp@0
|
1143 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),1) ) + ( double(planecorners( planenumb(ivsmaller),2 ))-1 )*nplanes;
|
|
tp@0
|
1144 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1145 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),2) ) + ( double(planecorners( planenumb(ivsmaller),2 ))-1 )*nplanes;
|
|
tp@0
|
1146 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1147 % Plane corner 3
|
|
tp@0
|
1148 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),1) ) + ( double(planecorners( planenumb(ivsmaller),3 ))-1 )*nplanes;
|
|
tp@0
|
1149 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1150 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),2) ) + ( double(planecorners( planenumb(ivsmaller),3 ))-1 )*nplanes;
|
|
tp@0
|
1151 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1152
|
|
tp@0
|
1153 if minncornersperplane == 4 & maxncornersperplane == 4
|
|
tp@0
|
1154 % Plane corner 4
|
|
tp@0
|
1155 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),1) ) + ( double(planecorners( planenumb(ivsmaller),4 ))-1 )*nplanes;
|
|
tp@0
|
1156 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1157 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller),2) ) + ( double(planecorners( planenumb(ivsmaller),4 ))-1 )*nplanes;
|
|
tp@0
|
1158 planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1159 elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
|
|
tp@0
|
1160 ivall = uint32(1:length(ivsmaller));
|
|
tp@0
|
1161 if minncornersperplane == 3
|
|
tp@0
|
1162 iv3cornerplanes = find(ncornersperplanevec(planenumb(ivsmaller))==3);
|
|
tp@0
|
1163 ivall(iv3cornerplanes) = [];
|
|
tp@0
|
1164 end
|
|
tp@0
|
1165 if maxncornersperplane == 4
|
|
tp@0
|
1166 % Plane corner 4
|
|
tp@0
|
1167 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),1) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1168 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1169 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),2) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1170 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1171 else
|
|
tp@0
|
1172 ivmorethan4 = find(ncornersperplanevec(planenumb(ivsmaller(ivall)))>4);
|
|
tp@0
|
1173 temp = ivall(ivmorethan4);
|
|
tp@0
|
1174 ivall(ivmorethan4) = [];
|
|
tp@0
|
1175 ivmorethan4 = temp;
|
|
tp@0
|
1176 % Plane corner 4
|
|
tp@0
|
1177 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),1) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1178 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1179 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivall)),2) ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1180 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1181 for ii = 5:maxncornersperplane
|
|
tp@0
|
1182 % Plane corner 5,6,7,etc
|
|
tp@0
|
1183 ivselection = find(ncornersperplanevec(planenumb(ivsmaller(ivmorethan4)))>=ii);
|
|
tp@0
|
1184 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivmorethan4)),1) ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
|
|
tp@0
|
1185 planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1186 ivreftocoplamatrix = double( planesatedge(edgenumb(ivsmaller(ivmorethan4)),2) ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
|
|
tp@0
|
1187 planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1188 end
|
|
tp@0
|
1189 clear ivselection
|
|
tp@0
|
1190 end
|
|
tp@0
|
1191 end
|
|
tp@0
|
1192 clear ivreftocoplamatrix
|
|
tp@0
|
1193
|
|
tp@0
|
1194 ivclear = ivsmaller(find(planecornerinfrontofedge==0));
|
|
tp@0
|
1195 edgeseesplane(iv(ivclear)) = 0;
|
|
tp@0
|
1196 end
|
|
tp@0
|
1197
|
|
tp@0
|
1198 if ~isempty(ivlarger)
|
|
tp@0
|
1199 % Second the edges larger than pi (closwedang > pi): at least one plane corner needs
|
|
tp@0
|
1200 % to be in front of *both* edge planes.
|
|
tp@0
|
1201 % Plane corner 1
|
|
tp@0
|
1202 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
|
|
tp@0
|
1203 planecornerinfrontofedge = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1204 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
|
|
tp@0
|
1205 planecornerinfrontofedge = planecornerinfrontofedge & (cornerinfrontofplane(ivreftocoplamatrix)==1);
|
|
tp@0
|
1206 % Plane corner 2
|
|
tp@0
|
1207 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
|
|
tp@0
|
1208 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1209 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
|
|
tp@0
|
1210 planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1211 % Plane corner 3
|
|
tp@0
|
1212 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
|
|
tp@0
|
1213 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1214 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
|
|
tp@0
|
1215 planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1216
|
|
tp@0
|
1217 if minncornersperplane == 4 & maxncornersperplane == 4
|
|
tp@0
|
1218 % Plane corner 4
|
|
tp@0
|
1219 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),1) ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
|
|
tp@0
|
1220 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1221 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger),2) ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
|
|
tp@0
|
1222 planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1223 elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
|
|
tp@0
|
1224 ivall = uint32(1:length(ivlarger));
|
|
tp@0
|
1225 if minncornersperplane == 3
|
|
tp@0
|
1226 iv3cornerplanes = find(ncornersperplanevec(planenumb(ivlarger))==3);
|
|
tp@0
|
1227 ivall(iv3cornerplanes) = [];
|
|
tp@0
|
1228 end
|
|
tp@0
|
1229 if maxncornersperplane == 4
|
|
tp@0
|
1230 % Plane corner 4
|
|
tp@0
|
1231 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),1) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1232 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1233 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),2) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1234 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1235 else
|
|
tp@0
|
1236 ivmorethan4 = find(ncornersperplanevec(planenumb(ivlarger(ivall)))>4);
|
|
tp@0
|
1237 temp = ivall(ivmorethan4);
|
|
tp@0
|
1238 ivall(ivmorethan4) = [];
|
|
tp@0
|
1239 ivmorethan4 = temp;
|
|
tp@0
|
1240 % Plane corner 4
|
|
tp@0
|
1241 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),1) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1242 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1243 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivall)),2) ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
|
|
tp@0
|
1244 planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1245 for ii = 5:maxncornersperplane
|
|
tp@0
|
1246 % Plane corner 5,6,7,etc
|
|
tp@0
|
1247 ivselection = find(ncornersperplanevec(planenumb(ivlarger(ivmorethan4)))>=ii);
|
|
tp@0
|
1248 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivmorethan4)),1) ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
|
|
tp@0
|
1249 temp = cornerinfrontofplane(ivreftocoplamatrix)==1;
|
|
tp@0
|
1250 ivreftocoplamatrix = double( planesatedge(edgenumb(ivlarger(ivmorethan4)),2) ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
|
|
tp@0
|
1251 planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
|
|
tp@0
|
1252 end
|
|
tp@0
|
1253 clear ivselection
|
|
tp@0
|
1254 end
|
|
tp@0
|
1255 end
|
|
tp@0
|
1256
|
|
tp@0
|
1257 ivclear = ivlarger(find(planecornerinfrontofedge==0));
|
|
tp@0
|
1258 edgeseesplane(iv(ivclear)) = 0;
|
|
tp@0
|
1259 end
|
|
tp@0
|
1260 end
|
|
tp@0
|
1261
|
|
tp@0
|
1262 clear ivlarger ivsmaller ivreftocoplamatrix edgenumb planenumb iv
|
|
tp@0
|
1263
|
|
tp@0
|
1264 %----------------------------------------------------------------------------
|
|
tp@0
|
1265 %
|
|
tp@0
|
1266 % indentingedgepairs
|
|
tp@0
|
1267 %
|
|
tp@0
|
1268 %--------------------------------------------------------------------------------
|
|
tp@0
|
1269 %
|
|
tp@0
|
1270 % If there are any planes with indents, make a list of all consecutive edge pairs that
|
|
tp@0
|
1271 % form an indent, because such edge pairs could never see each other.
|
|
tp@0
|
1272
|
|
tp@0
|
1273 indentingedgepairs = [];
|
|
tp@0
|
1274
|
|
tp@0
|
1275 if any(planehasindents)
|
|
tp@0
|
1276 iv = find(indentingcorners);
|
|
tp@0
|
1277 indentingedgepairs = zeros(length(iv),2);
|
|
tp@0
|
1278 indentingplanenumbers = mod(iv,nplanes);
|
|
tp@0
|
1279 ivzeros = find(indentingplanenumbers==0);
|
|
tp@0
|
1280 if ~isempty(ivzeros)
|
|
tp@0
|
1281 indentingplanenumbers(ivzeros) = nplanes;
|
|
tp@0
|
1282 end
|
|
tp@0
|
1283 conumbers = (iv - indentingplanenumbers)/nplanes+1;
|
|
tp@0
|
1284
|
|
tp@0
|
1285 % We expand planecorners temporarily, cyclically
|
|
tp@0
|
1286 planecorners = [planecorners planecorners(:,1:2)];
|
|
tp@0
|
1287
|
|
tp@0
|
1288 for ii = 1:length(iv)
|
|
tp@0
|
1289 edge1cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii):conumbers(ii)+1));
|
|
tp@0
|
1290 edge2cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii)+1:conumbers(ii)+2));
|
|
tp@0
|
1291
|
|
tp@0
|
1292 [slask,edge1] = ismember(edge1cornerpair,edgecorners,'rows');
|
|
tp@0
|
1293 [slask,edge2] = ismember(edge2cornerpair,edgecorners,'rows');
|
|
tp@0
|
1294
|
|
tp@0
|
1295 edgepair = sort([edge1 edge2]);
|
|
tp@0
|
1296 indentingedgepairs(ii,:) = edgepair;
|
|
tp@0
|
1297
|
|
tp@0
|
1298 end
|
|
tp@0
|
1299
|
|
tp@0
|
1300 % Remove the temporary expansion
|
|
tp@0
|
1301
|
|
tp@0
|
1302 planecorners = planecorners(:,1:size(planecorners,2)-2);
|
|
tp@0
|
1303 end
|
|
tp@0
|
1304
|
|
tp@0
|
1305 %----------------------------------------------------------------------------
|
|
tp@0
|
1306 %
|
|
tp@0
|
1307 % canplaneobstruct
|
|
tp@0
|
1308 %
|
|
tp@0
|
1309 %--------------------------------------------------------------------------------
|
|
tp@0
|
1310 %
|
|
tp@0
|
1311 % Check which planes that are potentially obstructing:
|
|
tp@0
|
1312 %
|
|
tp@0
|
1313 % If **all** closwedangvec are < pi and we have an exterior problem =>
|
|
tp@0
|
1314 % the scatterer is without indents.
|
|
tp@0
|
1315 %
|
|
tp@0
|
1316 % If **all** closwedangvec are > pi and we have an interior problem =>
|
|
tp@0
|
1317 % the room is without any indents so there can be no obstructions.
|
|
tp@0
|
1318 %
|
|
tp@0
|
1319 % For exterior problems, all active planes can potentially obstruct.
|
|
tp@0
|
1320 %
|
|
tp@0
|
1321 % If these simple checks give a results that conflict with the input
|
|
tp@0
|
1322 % parameters open_or_closed_model or int_or_ext_model, give an
|
|
tp@0
|
1323 % error message.
|
|
tp@0
|
1324
|
|
tp@0
|
1325 canplaneobstruct = [];
|
|
tp@0
|
1326 if sum(closwedangvec<pi) == 0
|
|
tp@0
|
1327 if SHOWTEXT >= 2
|
|
tp@0
|
1328 disp(' ')
|
|
tp@0
|
1329 disp(' The model is an interior problem, and the room has no indents')
|
|
tp@0
|
1330 disp(' so no obstruction tests will be necessary.')
|
|
tp@0
|
1331 end
|
|
tp@0
|
1332 if int_or_ext_model == 'e'
|
|
tp@0
|
1333 disp( 'ERROR: In the setup file this model was defined as an exterior problem, but it is clearly an interior problem.')
|
|
tp@0
|
1334 error(' Please check the orientation of the planes.')
|
|
tp@0
|
1335 end
|
|
tp@0
|
1336 canplaneobstruct = zeros(1,nplanes);
|
|
tp@0
|
1337 elseif sum(closwedangvec>0) == 0
|
|
tp@0
|
1338 if SHOWTEXT >= 2
|
|
tp@0
|
1339 disp(' ')
|
|
tp@0
|
1340 disp(' The model is an exterior problem (scattering)')
|
|
tp@0
|
1341 disp(' with only thin planes.')
|
|
tp@0
|
1342 end
|
|
tp@0
|
1343 if int_or_ext_model == 'i'
|
|
tp@0
|
1344 disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
|
|
tp@0
|
1345 error(' Please check the orientation of the planes.')
|
|
tp@0
|
1346 end
|
|
tp@0
|
1347 canplaneobstruct = ones(1,nplanes);
|
|
tp@0
|
1348 elseif sum(closwedangvec>pi) == 0
|
|
tp@0
|
1349 if SHOWTEXT >= 2
|
|
tp@0
|
1350 disp(' ')
|
|
tp@0
|
1351 disp(' The model is an exterior problem (scattering)')
|
|
tp@0
|
1352 disp(' and the scatterer has no indents.')
|
|
tp@0
|
1353 end
|
|
tp@0
|
1354 if int_or_ext_model == 'i'
|
|
tp@0
|
1355 disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
|
|
tp@0
|
1356 error(' Please check the orientation of the planes.')
|
|
tp@0
|
1357 end
|
|
tp@0
|
1358 canplaneobstruct = ones(1,nplanes);
|
|
tp@0
|
1359 end
|
|
tp@0
|
1360
|
|
tp@0
|
1361 % For an interior problem we know for sure that if a plane
|
|
tp@0
|
1362 % has all other points in front of itself, it can not obstruct.
|
|
tp@0
|
1363 % NB! We need to check only points that do not belong to planes
|
|
tp@0
|
1364 % that are aligned with the plane, or planes that are non-active.
|
|
tp@0
|
1365
|
|
tp@0
|
1366 if int_or_ext_model == 'e'
|
|
tp@0
|
1367 canplaneobstruct = (reflfactors.'~=0);
|
|
tp@0
|
1368 elseif isempty(canplaneobstruct)
|
|
tp@0
|
1369
|
|
tp@0
|
1370 canplaneobstruct = (reflfactors.'~=0);
|
|
tp@0
|
1371 maskvec1 = canplaneobstruct.';
|
|
tp@0
|
1372 listofactiveplanes = find(canplaneobstruct);
|
|
tp@0
|
1373 zerosvec = zeros(nplanes,1);
|
|
tp@0
|
1374
|
|
tp@0
|
1375 for ii = 1:length(listofactiveplanes)
|
|
tp@0
|
1376 iplane = listofactiveplanes(ii);
|
|
tp@0
|
1377 maskvec2 = zerosvec;
|
|
tp@0
|
1378 maskvec2(iplane) = 1;
|
|
tp@0
|
1379 otherplanestocheck = find((double(planeseesplane(:,iplane))+maskvec2)~=1 & planeseesplane(:,iplane)~=-1 & maskvec1~=0);
|
|
tp@0
|
1380 if ~isempty(otherplanestocheck),
|
|
tp@0
|
1381 cornerstocheck = unique(planecorners(otherplanestocheck,:));
|
|
tp@0
|
1382 cornerstocheck = setdiff(unique(cornerstocheck),planecorners(iplane,1:ncornersperplanevec(iplane)));
|
|
tp@0
|
1383 pointinfront = EDB1infrontofplane(corners(cornerstocheck,:),planenvecs(iplane,:),corners(planecorners(iplane,1),:),corners(planecorners(iplane,2),:));
|
|
tp@0
|
1384 if isempty(find(pointinfront==-1))
|
|
tp@0
|
1385 canplaneobstruct(iplane) = 0;
|
|
tp@0
|
1386 end
|
|
tp@0
|
1387 else
|
|
tp@0
|
1388 canplaneobstruct(iplane) = 0;
|
|
tp@0
|
1389 end
|
|
tp@0
|
1390
|
|
tp@0
|
1391 end
|
|
tp@0
|
1392 end
|
|
tp@0
|
1393
|
|
tp@0
|
1394 %---------------------------------------------------------------
|
|
tp@0
|
1395 % For each thin plane, make sure that all edges have the same normal
|
|
tp@0
|
1396 % vector. If not, switch the direction of that edge, and all the other
|
|
tp@0
|
1397 % relevant parameters.
|
|
tp@0
|
1398
|
|
tp@0
|
1399 if difforder >= 1
|
|
tp@0
|
1400 ivthin = find(planeisthin);
|
|
tp@0
|
1401
|
|
tp@0
|
1402 for ii = 1:length(ivthin)
|
|
tp@0
|
1403 plane = ivthin(ii);
|
|
tp@0
|
1404 if rearsideplane(plane) > plane
|
|
tp@0
|
1405 edgelist = unique(edgesatplane(plane,1:nedgesperplanevec(plane)));
|
|
tp@0
|
1406 iv = find(closwedangvec(edgelist,:)==0);
|
|
tp@0
|
1407 if ~isempty(iv)
|
|
tp@0
|
1408 edgelist = edgelist(iv);
|
|
tp@0
|
1409 nedgestemp = length(edgelist);
|
|
tp@0
|
1410 edgenveclist = edgenvecs(edgelist,:);
|
|
tp@0
|
1411 refnvec = edgenveclist(1,:);
|
|
tp@0
|
1412 nvecdiff = edgenveclist(2:nedgestemp,:) - refnvec(ones(nedgestemp-1,1),:);
|
|
tp@0
|
1413 nvecdiff = sum(abs(nvecdiff.')).';
|
|
tp@0
|
1414 ivswitch = find(nvecdiff)+1;
|
|
tp@0
|
1415 if ~isempty(ivswitch)
|
|
tp@0
|
1416 edgenumber = edgelist(ivswitch);
|
|
tp@0
|
1417 edgecorners(edgenumber,:) = [edgecorners(edgenumber,2) edgecorners(edgenumber,1)];
|
|
tp@0
|
1418 planesatedge(edgenumber,:) = [planesatedge(edgenumber,2) planesatedge(edgenumber,1)];
|
|
tp@0
|
1419 edgenvecs(edgenumber,:) = -edgenvecs(edgenumber,:);
|
|
tp@0
|
1420 edgestartcoords(edgenumber,:) = corners(edgecorners(edgenumber,1),:);
|
|
tp@0
|
1421 edgeendcoords(edgenumber,:) = corners(edgecorners(edgenumber,2),:);
|
|
tp@0
|
1422 tempvec = edgeendcoordsnudge(edgenumber,:);
|
|
tp@0
|
1423 edgeendcoordsnudge(edgenumber,:) = edgestartcoordsnudge(edgenumber,:);
|
|
tp@0
|
1424 edgestartcoordsnudge(edgenumber,:) = tempvec;
|
|
tp@0
|
1425 end
|
|
tp@0
|
1426 end
|
|
tp@0
|
1427 end
|
|
tp@0
|
1428
|
|
tp@0
|
1429 end
|
|
tp@0
|
1430 end
|
|
tp@0
|
1431
|
|
tp@0
|
1432 %----------------------------------------------------------------------------
|
|
tp@0
|
1433 %
|
|
tp@0
|
1434 % SAVE THE VARIABLES
|
|
tp@0
|
1435 %
|
|
tp@0
|
1436 %----------------------------------------------------------------------------
|
|
tp@0
|
1437 %
|
|
tp@0
|
1438 % Also the variables from cadgeo??
|
|
tp@0
|
1439 % Yes: planeeqs planenvecs ncornersperplanevec planecorners corners
|
|
tp@0
|
1440 % minvals maxvals
|
|
tp@0
|
1441
|
|
tp@0
|
1442 edgeseesplane = int8(edgeseesplane);
|
|
tp@0
|
1443 planeseesplane = int8(planeseesplane);
|
|
tp@0
|
1444
|
|
tp@0
|
1445 if ncorners < 256
|
|
tp@0
|
1446 planecorners = uint8(planecorners);
|
|
tp@0
|
1447 elseif ncorners < 65536
|
|
tp@0
|
1448 planecorners = uint16(planecorners);
|
|
tp@0
|
1449 end
|
|
tp@0
|
1450 planeisthin = uint8(planeisthin);
|
|
tp@0
|
1451
|
|
tp@0
|
1452 if max(ncornersperplanevec) <= 255
|
|
tp@0
|
1453 ncornersperplanevec = uint8(ncornersperplanevec);
|
|
tp@0
|
1454 else
|
|
tp@0
|
1455 ncornersperplanevec = uint16(ncornersperplanevec);
|
|
tp@0
|
1456 end
|
|
tp@0
|
1457
|
|
tp@0
|
1458 Varlist = [ ' edgecorners planesatedge closwedangvec planehasindents indentingedgepairs'];
|
|
tp@0
|
1459 Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
|
|
tp@0
|
1460 Varlist = [Varlist,' corners planecorners planeeqs planenvecs ncornersperplanevec cornerinfrontofplane'];
|
|
tp@0
|
1461 Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
|
|
tp@0
|
1462 Varlist = [Varlist,' minvals maxvals offedges '];
|
|
tp@0
|
1463 Varlist = [Varlist,' int_or_ext_model'];
|
|
tp@0
|
1464 if difforder >= 1
|
|
tp@0
|
1465 edgedata.edgestartcoordsnudge = edgestartcoordsnudge;
|
|
tp@0
|
1466 edgedata.edgeendcoordsnudge = edgeendcoordsnudge;
|
|
tp@0
|
1467 edgedata.edgenormvecs = edgenormvecs;
|
|
tp@0
|
1468 Varlist = [Varlist, ' edgedata '];
|
|
tp@0
|
1469 end
|
|
tp@0
|
1470 eval(['save ',outputfile,Varlist])
|
