Mercurial > hg > human-echolocation

annotate private/EDB1ed2geox.m @ 18:2d5f50205527 jabuilder_int tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Escape the trailing backslash as well
author Chris Cannam
date Tue, 30 Sep 2014 16:23:00 +0100
parents 90220f7249fc
children
rev   line source
tp@0 1 function outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches)
tp@0 2 % EDB1ed2geox - Calculates 2nd- and higher-order edge-related geom. parameters.
tp@0 3 % EDB1ed2geox 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 % eddatafile The .mat file where the output data will be stored.
tp@0 10 % If it is not specified, an automatic file name is constructed
tp@0 11 % using the same filename stem as the cadgeofile, but ending
tp@0 12 % with '_edgeo'.
tp@0 13 % sdatafile
tp@0 14 % rdatafile
tp@0 15 % doSbatches
tp@0 16 % doRbatches
tp@0 17 % specorder The highest order of any reflection kind (specular and/or diffraction).
tp@0 18 % difforder 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
tp@0 21 % ndiff2batches
tp@0 22 % SHOWTEXT (global) An integer value; if it is 4 or higher, then messages will be printed on the
tp@0 23 % screen during calculations.
tp@0 24 % Output parameters:
tp@0 25 % outputfile The name of the outputfile, generated either automatically or specified as
tp@0 26 % an input parameter.
tp@0 27 %
tp@0 28 % Output parameters in the outputfile:
tp@0 29 % (Taken directly from the cadgeo-file:)
tp@0 30 % corners planecorners planeeqs planenvecs ncornersperplanevec
tp@0 31 % minvals maxvals
tp@0 32 % (Taken directly from the setup-file:)
tp@0 33 % int_or_ext_model
tp@0 34 % (Calculated by EDB1edgeo:)
tp@0 35 % edgecorners Matrix [nedges,2] with the corner numbers that
tp@0 36 % define each edge, in the reference order.
tp@0 37 % edgestartcoords Matrix [nedges,3] with the coordinates of the
tp@0 38 % start corner for each edge.
tp@0 39 % edgeendcoords Matrix [nedges,3] with the coordinates of the
tp@0 40 % end corner for each edge.
tp@0 41 % edgelengthvec List [nedges,1] with the lengths of each edge.
tp@0 42 % planesatedge Matrix [nedges,2] with the plane numbers of the
tp@0 43 % two planes that connect to each edge. The first
tp@0 44 % plane is considered as the reference plane for
tp@0 45 % each edge. Rule: if the RH thumb is placed
tp@0 46 % along the edge, from start to end, and the hand
tp@0 47 % is rotated, the fingers should "come up
tp@0 48 % through" the reference plane.
tp@0 49 % closwedangvec List [nedges,1] with the angles in radians of
tp@0 50 % each edge - for the closed part of each edge.
tp@0 51 % edgenvecs Matrix [nedges,3] with the normal vectors of
tp@0 52 % the reference plane for each edge.
tp@0 53 % offedges List [nlength,1] where nlength = 0-nedges
tp@0 54 % containing the edge numbers that should not
tp@0 55 % used.
tp@0 56 % reflfactors List [nplanes,1] with the reflection factors
tp@0 57 % of each plane. The only supported values are
tp@0 58 % +1 for rigid planes, 0 for perfectlt absorbing
tp@0 59 % planes and -1 for pressure-release planes. The
tp@0 60 % pressure-release cases might not be implemented
tp@0 61 % fully.
tp@0 62 % planeisthin List [nplanes,1] with the value 1 or 0 stating
tp@0 63 % whether the plane is thin or not.
tp@0 64 % rearsideplane List [nplanes,1] with the plane number that is at
tp@0 65 % the rear side. If a plane is non-thin, the
tp@0 66 % value in this list is zero.
tp@0 67 % planehasindents List [nplanes,1] with the value 1 or 0 stating
tp@0 68 % whether a plane has indents or not.
tp@0 69 % canplaneobstruct List [nplanes,1] with the value 1 or 0 stating
tp@0 70 % whether a plane has the potential to obstruct
tp@0 71 % or not.
tp@0 72 % planeseesplane A matrix, [nplanes,nplanes], with the value 1
tp@0 73 % 0,-1,-2 stating:
tp@0 74 % 1 A plane is front of another plane and could
tp@0 75 % then potentially see it.
tp@0 76 % 0 A plane is completely behind another plane
tp@0 77 % but they are not aligned
tp@0 78 % -1 A plane is aligned with another plane
tp@0 79 % -2 A plane is back-to-back with another plane.
tp@0 80 % edgesatplane A matrix, [nplanes,nmaxedgesperplane] which for
tp@0 81 % row N contains the edge numbers that connect to
tp@0 82 % to plane N.
tp@0 83 % edgeseesplane A matrix, [nplanes,nedges], which contains the
tp@0 84 % values 0,1,-2 or -2, meaning:
tp@0 85 % 0 The edge can never see the plane
tp@0 86 % -1 The edge can never see the plane and the
tp@0 87 % edge is aligned with the plane
tp@0 88 % -2 The edge can never see the plane and the
tp@0 89 % edge belongs to the plane
tp@0 90 % 1 The edge can potentially see the plane.
tp@0 91 %
tp@0 92 % (Calculated by EDB1edgeo, but only for diffraction order >= 2:)
tp@0 93 % edgeseespartialedge A sparse matrix, [nedges,nedges], with the
tp@0 94 % values 0, pos. integer or neg. integer:
tp@0 95 % 0 Two edges can not see each other (or one
tp@0 96 % of the edges is inactive)
tp@0 97 % pos.int. A value from 1 to 2^nedgesubs-1
tp@0 98 % indicating how much of the edges
tp@0 99 % see each other. A pos. value
tp@0 100 % indicates that the edge-to-edge
tp@0 101 % path runs along a plane.
tp@0 102 % NB! A complete visibility test is
tp@0 103 % not made; subsegment 1 on edge 1 is only
tp@0 104 % checked against subsegment 1 on
tp@0 105 % edge 2, not against all other
tp@0 106 % subsegments of edge 2.
tp@0 107 % neg.int. Same as for the pos.int. except
tp@0 108 % that the edge-to-edge path does not
tp@0 109 % run along a plane.
tp@0 110 %
tp@0 111 % reftoshortlistE
tp@0 112 % re1sho, re2sho
tp@0 113 % thetae1sho, thetae2sho
tp@0 114 % ze1sho, ze2sho
tp@0 115 % examplecombE
tp@0 116 % edgealignedwithedge A sparse matrix, [nedges,nedges], with the value 1
tp@0 117 % or 0 stating whether an (active) edge is aligned
tp@0 118 % with another (active) edge or not.
tp@0 119 % edgeperptoplane A sparse matrix, [nplanes,nedges], with the value 1
tp@0 120 % or 0 stating whether an (active) edge is perpendicular
tp@0 121 % to an active plane or not.
tp@0 122 % edgeplaneperptoplane1 A sparse matrix, [nplanes,nedges], with the value 1
tp@0 123 % or 0 stating whether an (active) edge has one
tp@0 124 % of its two defining planes perpendicular to
tp@0 125 % another plane with which it shares an edge.
tp@0 126 % edgeplaneperptoplane2 A sparse matrix, [nplanes,nedges], with the value 1
tp@0 127 % or 0 stating whether an (active) edge has one
tp@0 128 % of its two defining planes perpendicular to
tp@0 129 % another plane which has a flat edge (i.e., a 180 degree edge).
tp@0 130 %
tp@0 131 % Uses the functions EDB1coordtrans1, EDB1coordtrans2
tp@0 132 % EDB1infrontofplane, EDB1compress7, EDB1strpend, EDB1cross
tp@0 133 % EDB1checkobstr_edgetoedge
tp@0 134 %
tp@0 135 % ----------------------------------------------------------------------------------------------
tp@0 136 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 137 %
tp@0 138 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 139 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 140 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 141 %
tp@0 142 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 143 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 144 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 145 %
tp@0 146 % You should have received a copy of the GNU General Public License along with the
tp@0 147 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 148 % ----------------------------------------------------------------------------------------------
tp@0 149 % Peter Svensson (svensson@iet.ntnu.no) 20110822
tp@0 150 %
tp@0 151 % outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches);
tp@0 152
tp@0 153 % 20110822 Bug found and fixed, giving the wrong amplitude for second- and higher-order
tp@0 154 % diffraction, but only for diffraction paths that do not travel along a surface.
tp@0 155
tp@0 156 global SHOWTEXT
tp@0 157
tp@0 158 geomacc = 1e-10;
tp@0 159
tp@0 160 eval(['load ',eddatafile])
tp@0 161 edgestartcoordsnudge = edgedata.edgestartcoordsnudge;
tp@0 162 edgeendcoordsnudge = edgedata.edgeendcoordsnudge;
tp@0 163 edgenormvecs = edgedata.edgenormvecs;
tp@0 164
tp@0 165 [sfilepath,sfilestem,temp1] = fileparts(sdatafile);
tp@0 166 sdatafile = [sfilepath,filesep,sfilestem];
tp@0 167 if doSbatches == 0
tp@0 168 eval(['load ',sdatafile,'.mat'])
tp@0 169 else
tp@0 170 eval(['load ',sdatafile,'_1.mat'])
tp@0 171 end
tp@0 172 [rfilepath,rfilestem,temp1] = fileparts(rdatafile);
tp@0 173 rdatafile = [rfilepath,filesep,rfilestem];
tp@0 174 if doRbatches == 0
tp@0 175 eval(['load ',rdatafile,'.mat'])
tp@0 176 else
tp@0 177 eval(['load ',rdatafile,'_1.mat'])
tp@0 178 end
tp@0 179
tp@0 180 outputfile = eddatafile;
tp@0 181
tp@0 182 nplanes = size(planecorners,1);
tp@0 183 nedges = size(edgecorners,1);
tp@0 184
tp@0 185 zerosvec1 = zeros(nplanes,1);
tp@0 186
tp@0 187 %---------------------------------------------------------------
tp@0 188 % We check right away which edges that can not be seen by either the source
tp@0 189 % or the receiver. We don't need to check edgeseesedge for combinations
tp@0 190 % that involve any of these.
tp@0 191 %
tp@0 192 % Note that this is true only for difforder = 2! If difforder >= 3, then we
tp@0 193 % can not exclude any edges based on this! (It could be possible to exclude
tp@0 194 % edges, if we managed to "unwind" the
tp@0 195 % edge-sees-an-edge-which-sees-an-edge-which-can-be-seen-by-the-source.
tp@0 196
tp@0 197 if difforder == 2
tp@0 198 edgesthatareseen = any( [vispartedgesfroms vispartedgesfromr].' );
tp@0 199 edgesnottoworryabout = [1:nedges];
tp@0 200 edgesnottoworryabout(edgesthatareseen) = [];
tp@0 201 clear vispartedgesfroms vispartedgesfromr visplanesfroms visplanesfromr
tp@0 202 else
tp@0 203 edgesnottoworryabout = [];
tp@0 204 end
tp@0 205
tp@0 206 %---------------------------------------------------------------
tp@0 207 %
tp@0 208 % edgeperptoplane
tp@0 209 %
tp@0 210 %---------------------------------------------------------------
tp@0 211 % Make matrices over:
tp@0 212 % 1. Which edges are perpendicular to planes ("edgeperptoplane")
tp@0 213 % 2. Which edges can give combinations with edge1-plane-edge1
tp@0 214 % where plane is perpendicular to one of the two planes that
tp@0 215 % defines the edge ("edgeplaneperptoplane1")
tp@0 216 % 3. Which edges can give combinations with edge1-plane-edge1
tp@0 217 % where plane is a part of a (unnecessarily) divided plane
tp@0 218 % ("edgeplaneperptoplane2")
tp@0 219 %
tp@0 220 % NB! Only active edges and active planes can get
tp@0 221 % such combinations.
tp@0 222 %
tp@0 223 % These are needed for diff-spec-diff combos, so only when
tp@0 224 % difforder >= 2 (and specorder >= 3).
tp@0 225
tp@0 226 edgeperptoplane = [];
tp@0 227 edgeplaneperptoplane1 = [];
tp@0 228 edgeplaneperptoplane2 = [];
tp@0 229 if specorder >= 3
tp@0 230 edgeperptoplane = sparse(zeros(nplanes,nedges));
tp@0 231 edgeplaneperptoplane1 = sparse(zeros(nplanes,nedges));
tp@0 232 edgeplaneperptoplane2 = sparse(zeros(nplanes,nedges));
tp@0 233 if SHOWTEXT >= 4
tp@0 234 disp(' checking which edges are perpendicular to planes ...')
tp@0 235 end
tp@0 236
tp@0 237 iv = find(edgeseesplane == 1);
tp@0 238 if ~isempty(iv)
tp@0 239 edgelist = floor((iv-1)/nplanes)+1;
tp@0 240 planelist = iv - (edgelist-1)*nplanes;
tp@0 241
tp@0 242 % Check first which edges have vectors that are parallel to plane
tp@0 243 % normal vectors. These are "edgeperptoplane" cases. Clear them
tp@0 244 % afterwards from the edgelist and planelist.
tp@0 245
tp@0 246 vec1 = edgenormvecs(edgelist,:);
tp@0 247 vec2 = planenvecs(planelist,:);
tp@0 248 ivperp = find( abs( sum( (vec1.*vec2).' ) )==1);
tp@0 249 if ~isempty(ivperp)
tp@0 250 edgeperptoplane(iv(ivperp)) = ones(size(ivperp));
tp@0 251 clear iv
tp@0 252 edgelist(ivperp) = [];
tp@0 253 planelist(ivperp) = [];
tp@0 254 vec2(ivperp,:) = [];
tp@0 255 else
tp@0 256 clear iv
tp@0 257 end
tp@0 258
tp@0 259 % Now check which edges have plane1 or plane2 perpendicular to
tp@0 260 % other planes.
tp@0 261
tp@0 262 if ~isempty(planelist)
tp@0 263
tp@0 264 vec1 = planenvecs(planesatedge(edgelist,1),:);
tp@0 265 ivplaneperp1 = find( abs( sum( (vec1.*vec2).' ) )==0);
tp@0 266 vec1 = planenvecs(planesatedge(edgelist,2),:);
tp@0 267 ivplaneperp2 = find( abs( sum( (vec1.*vec2).' ) )==0);
tp@0 268 ivplaneperp = [ivplaneperp1 ivplaneperp2];
tp@0 269 if ~isempty(ivplaneperp)
tp@0 270 ivplaneperp = unique(ivplaneperp);
tp@0 271 edgelist = edgelist(ivplaneperp);
tp@0 272 planelist = planelist(ivplaneperp);
tp@0 273 edgeplane1 = planesatedge(edgelist,1);
tp@0 274 edgeplane2 = planesatedge(edgelist,2);
tp@0 275
tp@0 276 % Check which of these other perpendicular planes that connnect
tp@0 277 % to one of the edge planes via a 90 deg. corner.
tp@0 278
tp@0 279 ivtoclear = [];
tp@0 280 [ivfoundcombos,loc] = ismember([planelist edgeplane1],planesatedge,'rows');
tp@0 281 ivperp = find(ivfoundcombos);
tp@0 282 if ~isempty(ivperp)
tp@0 283 edgenumb = loc(ivperp);
tp@0 284 ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0 285 if ~isempty(ivfinalcomb)
tp@0 286 ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0 287 edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 288 ivtoclear = ivfinalcomb;
tp@0 289 end
tp@0 290 end
tp@0 291 [ivfoundcombos,loc] = ismember([edgeplane1 planelist],planesatedge,'rows');
tp@0 292 ivperp = find(ivfoundcombos);
tp@0 293 if ~isempty(ivperp)
tp@0 294 edgenumb = loc(ivperp);
tp@0 295 ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0 296 if ~isempty(ivfinalcomb)
tp@0 297 ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0 298 edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 299
tp@0 300 ivtoclear = ivtoclear.';
tp@0 301 ivfinalcomb = ivfinalcomb.';
tp@0 302 ivtoclear = [ivtoclear ivfinalcomb];
tp@0 303 end
tp@0 304 end
tp@0 305
tp@0 306 [ivfoundcombos,loc] = ismember([planelist edgeplane2],planesatedge,'rows');
tp@0 307 ivperp = find(ivfoundcombos);
tp@0 308 if ~isempty(ivperp)
tp@0 309 edgenumb = loc(ivperp);
tp@0 310 ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0 311 if ~isempty(ivfinalcomb)
tp@0 312 ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0 313 edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 314
tp@0 315 if size(ivtoclear,2) == size(ivfinalcomb,2)
tp@0 316 ivtoclear = [ivtoclear;ivfinalcomb];
tp@0 317 elseif size(ivtoclear,1) == size(ivfinalcomb,1),
tp@0 318 ivfinalcomb = ivfinalcomb.';
tp@0 319 ivtoclear = [ivtoclear ivfinalcomb];
tp@0 320 elseif isempty(ivtoclear)
tp@0 321 error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')
tp@0 322 end
tp@0 323 end
tp@0 324 end
tp@0 325
tp@0 326
tp@0 327 [ivfoundcombos,loc] = ismember([edgeplane2 planelist],planesatedge,'rows');
tp@0 328 ivperp = find(ivfoundcombos);
tp@0 329 if ~isempty(ivperp)
tp@0 330 edgenumb = loc(ivperp);
tp@0 331 ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0 332 if ~isempty(ivfinalcomb)
tp@0 333 ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0 334 edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 335
tp@0 336 if size(ivtoclear,2) == size(ivfinalcomb,2)
tp@0 337 ivtoclear = [ivtoclear;ivfinalcomb];
tp@0 338 elseif size(ivtoclear,1) == size(ivfinalcomb,1),
tp@0 339 ivfinalcomb = ivfinalcomb.';
tp@0 340 ivtoclear = [ivtoclear ivfinalcomb];
tp@0 341 elseif isempty(ivtoclear)
tp@0 342 error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')
tp@0 343 end
tp@0 344 end
tp@0 345 end
tp@0 346
tp@0 347 % Finally check which edges have plane1 or plane2 perpendicular to
tp@0 348 % two other planes that connnect via a 180 deg. corner.
tp@0 349
tp@0 350 if ~isempty(ivtoclear)
tp@0 351 ivtoclear = unique(ivtoclear);
tp@0 352 edgelist(ivtoclear) = [];
tp@0 353 planelist(ivtoclear) = [];
tp@0 354 if ~isempty(edgelist)
tp@0 355
tp@0 356 ivflatedges = find(closwedangvec==pi);
tp@0 357 if ~isempty(ivflatedges)
tp@0 358 planecombos = planesatedge(ivflatedges,:);
tp@0 359 for ii = 1:length(ivflatedges)
tp@0 360 iv1 = find(planelist == planecombos(ii,1));
tp@0 361 if ~isempty(iv1)
tp@0 362 ivreftomatrix = double(planecombos(ii,1)) + (edgelist(iv1)-1)*nplanes;
tp@0 363 edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 364 end
tp@0 365 iv2 = find(planelist == planecombos(ii,2));
tp@0 366 if ~isempty(iv2)
tp@0 367 ivreftomatrix = planecombos(ii,2) + (edgelist(iv2)-1)*nplanes;
tp@0 368 edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 369 end
tp@0 370
tp@0 371 end
tp@0 372 end
tp@0 373
tp@0 374 end
tp@0 375 end
tp@0 376
tp@0 377 end
tp@0 378 end
tp@0 379
tp@0 380 end
tp@0 381
tp@0 382 end
tp@0 383
tp@0 384 %---------------------------------------------------------------
tp@0 385 %
tp@0 386 % edgealignedwithedge
tp@0 387 %
tp@0 388 %---------------------------------------------------------------
tp@0 389
tp@0 390 edgealignedwithedge = [];
tp@0 391
tp@0 392 if specorder >= 2
tp@0 393 edgealignedwithedge = sparse(eye(nedges));
tp@0 394 if SHOWTEXT >= 4
tp@0 395 disp(' checking which edges are aligned with other edges ...')
tp@0 396 end
tp@0 397
tp@0 398 listofactiveedges = [1:nedges].';
tp@0 399 listofactiveedges(offedges) = [];
tp@0 400 vecstocheck = edgenormvecs(listofactiveedges,:);
tp@0 401 iv = find(vecstocheck(:,1)<0);
tp@0 402 vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0 403 iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)<0);
tp@0 404 vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0 405 iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)==0 & vecstocheck(:,3)<0);
tp@0 406 vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0 407 if ~isempty(listofactiveedges)
tp@0 408 [uniquevecs,II,JJ] = unique(vecstocheck,'rows');
tp@0 409 N = histc(JJ,[1:length(listofactiveedges)]);
tp@0 410 iv = find(N>1);
tp@0 411 for ii = 1:length(iv)
tp@0 412 iv2 = find(JJ==iv(ii));
tp@0 413 edgestocheck = listofactiveedges(iv2);
tp@0 414 ntocheck = length(edgestocheck);
tp@0 415 cornernumberstocheck = edgecorners(edgestocheck,:);
tp@0 416 expandedstartcornermatrix = [reshape(repmat(cornernumberstocheck(:,1).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,1),[ntocheck 1])];
tp@0 417 expandedendcornermatrix = [reshape(repmat(cornernumberstocheck(:,2).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,2),[ntocheck 1])];
tp@0 418 toexpandededgestocheck = repmat(edgestocheck, [ntocheck 1]);
tp@0 419 fromexpandededgestocheck = reshape(repmat(edgestocheck.',[ntocheck 1]),ntocheck^2,1);
tp@0 420
tp@0 421 % Now we don't need to check edges against themselves, or edge1
tp@0 422 % vs edge2 if edge2 vs edge1 has already been tested.
tp@0 423
tp@0 424 indmat = reshape([1:ntocheck^2],ntocheck,ntocheck);
tp@0 425 indmat = triu(indmat,1);
tp@0 426 indmat = indmat(find(indmat));
tp@0 427
tp@0 428 expandedstartcornermatrix = expandedstartcornermatrix(indmat,:);
tp@0 429 expandedendcornermatrix = expandedendcornermatrix(indmat,:);
tp@0 430 fromexpandededgestocheck = fromexpandededgestocheck(indmat,:);
tp@0 431 toexpandededgestocheck = toexpandededgestocheck(indmat,:);
tp@0 432
tp@0 433 % We make the easiest check first: if two edges with the same
tp@0 434 % direction vector share one corner, they must be aligned.
tp@0 435
tp@0 436 A = [expandedstartcornermatrix expandedendcornermatrix];
tp@0 437 A = sort(A.').';
tp@0 438 A = diff(A.').';
tp@0 439 iv3 = find(sum(A.'==0).');
tp@0 440 validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
tp@0 441 if ~isempty(validcombs)
tp@0 442 indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
tp@0 443 edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0 444 indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
tp@0 445 edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0 446 fromexpandededgestocheck(iv3) = [];
tp@0 447 toexpandededgestocheck(iv3) = [];
tp@0 448 expandedstartcornermatrix(iv3,:) = [];
tp@0 449 expandedendcornermatrix(iv3,:) = [];
tp@0 450 end
tp@0 451
tp@0 452 % For the other edge-pair combinations we must check if the two
tp@0 453 % edges are really aligned. We do this by checking if two
tp@0 454 % vectors are aligned: edge1-vector and the vector from
tp@0 455 % edge1start to edge2start.
tp@0 456
tp@0 457 direcvec = corners(expandedstartcornermatrix(:,2),:) - corners(expandedstartcornermatrix(:,1),:);
tp@0 458 direcveclengths = sqrt( sum(direcvec.'.^2) ).';
tp@0 459 direcvec = direcvec./(direcveclengths(:,ones(1,3)) + eps*10);
tp@0 460 test = abs(sum( (direcvec.*edgenormvecs(fromexpandededgestocheck,:)).' ));
tp@0 461 iv3 = find(abs(test-1)< 1e-8);
tp@0 462 validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
tp@0 463 if ~isempty(validcombs)
tp@0 464 indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
tp@0 465 edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0 466 indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
tp@0 467 edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0 468 fromexpandededgestocheck(iv3) = [];
tp@0 469 toexpandededgestocheck(iv3) = [];
tp@0 470 expandedstartcornermatrix(iv3,:) = [];
tp@0 471 expandedendcornermatrix(iv3,:) = [];
tp@0 472 end
tp@0 473 end
tp@0 474 end
tp@0 475
tp@0 476 end
tp@0 477
tp@0 478 %---------------------------------------------------------------
tp@0 479 %
tp@0 480 % edgeseesedge
tp@0 481 %
tp@0 482 %---------------------------------------------------------------
tp@0 483 %
tp@0 484 % Make a matrix of which edges can see which edges
tp@0 485 %
tp@0 486 % Planes that are in front of the edge-planes might have visible edges.
tp@0 487 % If closwedang = 0, then all other edges are potentially visible.
tp@0 488 % If closwedang < pi, then all planes in front of plane1 or plane2 are OK
tp@0 489 % If closwedang > pi, then all planes in front of plane1 and plane2 are OK
tp@0 490
tp@0 491 if SHOWTEXT >= 4
tp@0 492 disp(' checking which edges see which edges...')
tp@0 493 end
tp@0 494
tp@0 495 edgeseesedge = int8(zeros(nedges,nedges));
tp@0 496 for ii = 1:nedges
tp@0 497 closwed = closwedangvec(ii);
tp@0 498 if closwed == 0
tp@0 499 edgeseesedge(:,ii) = ones(nedges,1);
tp@0 500 else
tp@0 501 plane1 = planesatedge(ii,1);
tp@0 502 plane2 = planesatedge(ii,2);
tp@0 503 okplanelist1 = find(planeseesplane(:,plane1)==1);
tp@0 504 okplanelist2 = find(planeseesplane(:,plane2)==1);
tp@0 505
tp@0 506 if closwed < pi
tp@0 507 okplanes = [okplanelist1;okplanelist2];
tp@0 508 okplanes = sort(okplanes);
tp@0 509 else
tp@0 510 okplanes = zerosvec1;
tp@0 511 okplanes(okplanelist1) = okplanes(okplanelist1)+1;
tp@0 512 okplanes(okplanelist2) = okplanes(okplanelist2)+1;
tp@0 513 okplanes = find(okplanes==2);
tp@0 514 end
tp@0 515
tp@0 516 okedges = edgesatplane(okplanes,:);
tp@0 517 [n1,n2] = size(okedges);
tp@0 518 okedges = reshape(okedges,n1*n2,1);
tp@0 519
tp@0 520 if ~isempty(okedges)
tp@0 521 okedges = okedges(find(okedges~=0));
tp@0 522 edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1);
tp@0 523 end
tp@0 524
tp@0 525 edgesatplane1 = edgesatplane(plane1,:);
tp@0 526 edgesatplane2 = edgesatplane(plane2,:);
tp@0 527 okedges = [edgesatplane1 edgesatplane2].';
tp@0 528 okedges = okedges(find(okedges~=0));
tp@0 529 edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1);
tp@0 530
tp@0 531 end % else ..... if closwed == 0
tp@0 532
tp@0 533 end % for ii = 1:nedges
tp@0 534
tp@0 535 % Mask out all edge pairs that are aligned
tp@0 536 edgeseesedge = double(edgeseesedge) - double(edgealignedwithedge);
tp@0 537
tp@0 538 % Make sure the edge-to-same-edge combos are masked out
tp@0 539 % edgeseesedge = sign(edgeseesedge.*(edgeseesedge.')).*(1-eye(nedges));
tp@0 540 edgeseesedge = int8( ( (edgeseesedge & (edgeseesedge.'))>0).*(1-eye(nedges)));
tp@0 541
tp@0 542 % Mask out all edges that should be switched off
tp@0 543
tp@0 544 mask = ones(nedges,nedges);
tp@0 545 mask(offedges,:) = mask(offedges,:)*0;
tp@0 546 mask(:,offedges.') = mask(:,offedges.')*0;
tp@0 547
tp@0 548 edgeseesedge = edgeseesedge & mask;
tp@0 549 clear mask
tp@0 550
tp@0 551
tp@0 552 %-----------------------------------------------------------
tp@0 553 % Go through all edgeseesedge combinations. If two edges see
tp@0 554 % each other without belonging to the same plane, there should
tp@0 555 % be a gain factor of 2.
tp@0 556
tp@0 557 if SHOWTEXT >= 4
tp@0 558 disp(' checking which edge-to-edge paths that run along planes')
tp@0 559 end
tp@0 560
tp@0 561 edgetoedgestrength = 2*edgeseesedge;
tp@0 562 for ii = 1:nedges
tp@0 563 if sum(edgetoedgestrength(:,ii)) ~= 0
tp@0 564 plane1 = planesatedge(ii,1);
tp@0 565 plane2 = planesatedge(ii,2);
tp@0 566 edgesatplane1 = edgesatplane(plane1,:);
tp@0 567 edgesatplane2 = edgesatplane(plane2,:);
tp@0 568 sameplaneedges = [edgesatplane1 edgesatplane2].';
tp@0 569 sameplaneedges = sameplaneedges(find(sameplaneedges~=0));
tp@0 570 edgetoedgestrength(sameplaneedges,ii) = edgetoedgestrength(sameplaneedges,ii)*0 + 1;
tp@0 571 end % % if sum(edge
tp@0 572 end
tp@0 573
tp@0 574 % an edge can not contribute to itself
tp@0 575
tp@0 576 edgetoedgestrength = edgetoedgestrength.*(1-eye(nedges));
tp@0 577
tp@0 578 % The edges belonging to the same planes, but that are inactive (90 degrees etc)
tp@0 579 % must be switched off here too.
tp@0 580 % 011021 This is redundant.
tp@0 581 edgetoedgestrength = edgetoedgestrength.*(edgeseesedge>0);
tp@0 582 edgetoedgestrength = int8(edgetoedgestrength);
tp@0 583
tp@0 584 %----------------------------------------------------------------------------
tp@0 585 %
tp@0 586 % CHECK OBSTRUCTION OF EDGE-TO-EDGE PATHS
tp@0 587 %
tp@0 588 %----------------------------------------------------------------------------
tp@0 589 %
tp@0 590 % We construct two long lists: 'from-edges' and 'to-edges' for which
tp@0 591 % obstruction should be tested. These lists are made up of all the
tp@0 592 % combinations in edgeseesedge that have a value 1.
tp@0 593 % NB! We use symmetry - if edge 1 can see edge 2, then edge 2 can also see
tp@0 594 % edge 1.
tp@0 595 %
tp@0 596 % Both edges should be subdivided into nedgesubs segments. So, the long
tp@0 597 % lists must be expanded so that each from-edge/to-edge combination is
tp@0 598 % replaced by nedgesubs^2 as many.
tp@0 599 % For efficiency, we make a shortlist of the actual edge subdivisions and
tp@0 600 % we then expand long lists with pointers to this shortlist.
tp@0 601 %
tp@0 602 % The final matrix, edgeseespartialedge, will have an integer value which
tp@0 603 % tells whether two edges see each other:
tp@0 604 %
tp@0 605 % 0 Edges can not see each other, or one of the edges is inactive
tp@0 606 % (an inactive edge has an integer wedge index, or has one plane
tp@0 607 % which is TOTABS)
tp@0 608 % 2^nedgesubs-1 Two edges see each other completely.
tp@0 609 % Other integers Two edges see each other partly.
tp@0 610 % Neg. integer As above, but in addition, the edges do not belong
tp@0 611 % to the same plane.
tp@0 612 %
tp@0 613 % We check in the visedgesfroms and visedgesfromr lists to check which
tp@0 614 % combinations we don't need to check.
tp@0 615
tp@0 616 if SHOWTEXT >= 3
tp@0 617 disp([' Checking for obstructing planes between edges and edges'])
tp@0 618 disp(' ')
tp@0 619 disp('NB!!! The value of nedgesubs is temporarily set to 1 for the edge-to-edge vis. test')
tp@0 620 end
tp@0 621
tp@0 622 obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0 623 maxedgetoedgevisvalue = 2^(2*nedgesubs)-1;
tp@0 624
tp@0 625 if obstructtestneeded
tp@0 626
tp@0 627 nedgesubsorig = nedgesubs;
tp@0 628 nedgesubs = 1;
tp@0 629 edgeseesedge = triu(edgeseesedge);
tp@0 630
tp@0 631 iv = full(find(edgeseesedge~=0));
tp@0 632
tp@0 633 if ~isempty(iv)
tp@0 634 fromedges = floor(iv/nedges)+1;
tp@0 635 toedges = iv - (fromedges-1)*nedges;
tp@0 636 ncombs = length(fromedges);
tp@0 637
tp@0 638 ivcancel = find(ismember(fromedges,edgesnottoworryabout));
tp@0 639 fromedges(ivcancel) = [];
tp@0 640 toedges(ivcancel) = [];
tp@0 641
tp@0 642 ivcancel = find(ismember(toedges,edgesnottoworryabout));
tp@0 643 toedges(ivcancel) = [];
tp@0 644 fromedges(ivcancel) = [];
tp@0 645 clear ivcancel
tp@0 646
tp@0 647 % Some of these fromedges-toedges combos might involve two
tp@0 648 % neighboring edges that form an indent. Those should be
tp@0 649 % removed before we start checking fro obstructions.
tp@0 650 %
tp@0 651 % We check the [fromedges toedges] matrix against the
tp@0 652 % indentingedgepairs matrix.
tp@0 653
tp@0 654 [tf,loc]= ismember(indentingedgepairs,sort([fromedges toedges].').','rows');
tp@0 655 if ~isempty(loc)
tp@0 656 ivmatch = find(loc);
tp@0 657 fromedges(loc(ivmatch)) = [];
tp@0 658 toedges(loc(ivmatch)) = [];
tp@0 659 end
tp@0 660 ncombs = length(fromedges);
tp@0 661
tp@0 662 [uniqueedges,iunique,junique] = unique([fromedges toedges]);
tp@0 663 clear fromedges toedges
tp@0 664
tp@0 665 % The lists edgesubcoords are the shortlists.
tp@0 666
tp@0 667 [edgesubcoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(uniqueedges,:),edgeendcoords(uniqueedges,:),edgelengthvec(uniqueedges,:),nedgesubs,1);
tp@0 668
tp@0 669 % The two lists below contain pointers to the first edge segment in the
tp@0 670 % shortlist.
tp@0 671
tp@0 672 fromcoords_refp1 = nedgesubs*(junique(1:ncombs)-1)+1;
tp@0 673 tocoords_refp1 = nedgesubs*(junique(ncombs+1:2*ncombs)-1)+1;
tp@0 674 clear junique
tp@0 675
tp@0 676 % Now the original [fromedges toedges] could be recovered by:
tp@0 677 % uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)])
tp@0 678 % reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),6,2)
tp@0 679 % uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]);
tp@0 680 % npairs = length(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]));
tp@0 681 % reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),npairs/2,2)
tp@0 682
tp@0 683 addmask = [0:nedgesubs-1];
tp@0 684 onesvec1 = ones(1,nedgesubs);
tp@0 685 ntot1 = ncombs*nedgesubs;
tp@0 686
tp@0 687 % First, the from-list gets repetitions of the same values
tp@0 688 expandfrom_ref = fromcoords_refp1(:,onesvec1);
tp@0 689 clear fromcoords_refp1;
tp@0 690 expandfrom_ref = reshape(expandfrom_ref.',ntot1,1);
tp@0 691
tp@0 692 % Second, the expanded from-list gets repetitions that are
tp@0 693 % incremented in steps of 1 (because the shortlists have the
tp@0 694 % edge segments placed after each other).
tp@0 695 expandfrom_ref = expandfrom_ref(:,onesvec1);
tp@0 696 expandfrom_ref = expandfrom_ref + addmask(ones(ntot1,1),:);
tp@0 697 expandfrom_ref = reshape(expandfrom_ref.',ncombs*nedgesubs^2,1);
tp@0 698
tp@0 699 % Same thing for the to-list except that the first expansion
tp@0 700 % gives repetitions that are incremented in steps of 1.
tp@0 701 expandto_ref = tocoords_refp1(:,onesvec1);
tp@0 702 clear tocoords_refp1;
tp@0 703 expandto_ref = expandto_ref + addmask(ones(ncombs,1),:);
tp@0 704 expandto_ref = reshape(expandto_ref.',ntot1,1);
tp@0 705
tp@0 706 % Second expansion for the to-list: repetitions
tp@0 707 expandto_ref = expandto_ref(:,onesvec1);
tp@0 708 expandto_ref = reshape(expandto_ref.',ncombs*nedgesubs^2,1);
tp@0 709
tp@0 710 % Now we have expanded lists of all combinations, and we can
tp@0 711 % make lists with the coordinates, the edge numbers, and the weights of
tp@0 712 % the segments.
tp@0 713
tp@0 714 % edgesubcoords(expandfrom_ref,:) will give all the starting points
tp@0 715 % "fromcoords"
tp@0 716 % edgesubcoords(expandto_ref,:) will give all the ending points
tp@0 717 % "tocoords"
tp@0 718 % uniqueedges(edgenumberlist(expandfrom_ref)) will give all
tp@0 719 % the starting edges, "fromedges"
tp@0 720 % uniqueedges(edgenumberlist(expandto_ref)) will give all
tp@0 721 % the ending edges, "toedges"
tp@0 722 % We can make a simple test of which planes that could be excluded
tp@0 723 % from the obstruction test by checking the edgeseesplane for the
tp@0 724 % relevant edges. It can be multiplied by the canplaneobstruct list
tp@0 725 % when calling the findobstructedpaths function.
tp@0 726 isplaneactiveforobstruction = (sum(edgeseesplane(:,uniqueedges).'==1)>0);
tp@0 727
tp@0 728 global STARTPLANES ENDPLANES
tp@0 729 global BIGFROMCOORDS BIGTOCOORDS BIGSTARTPLANES BIGENDPLANES
tp@0 730 global REFTOFROMCOSHO REFTOTOCOSHO
tp@0 731
tp@0 732 FROMCOORDSSHORTLIST = edgesubcoords;
tp@0 733 REFTOFROMCOSHO = uint32(expandfrom_ref);
tp@0 734 TOCOORDSSHORTLIST = edgesubcoords;
tp@0 735 REFTOTOCOSHO = uint32(expandto_ref);
tp@0 736 if nedges < 256
tp@0 737 bigfromedge = uint8(uniqueedges(edgenumberlist(expandfrom_ref)).');
tp@0 738 else
tp@0 739 bigfromedge = uint16(uniqueedges(edgenumberlist(expandfrom_ref)).');
tp@0 740 end
tp@0 741 clear expandfrom_ref
tp@0 742 bigfromedge = bigfromedge(:);
tp@0 743 if nedges < 256
tp@0 744 bigtoedge = uint8(uniqueedges(edgenumberlist(expandto_ref)).');
tp@0 745 else
tp@0 746 bigtoedge = uint16(uniqueedges(edgenumberlist(expandto_ref)).');
tp@0 747 end
tp@0 748 clear expandto_ref
tp@0 749 bigtoedge = bigtoedge(:);
tp@0 750 STARTPLANES = [planesatedge(bigfromedge,1) planesatedge(bigfromedge,2)];
tp@0 751 ENDPLANES = [planesatedge(bigtoedge,1) planesatedge(bigtoedge,2)];
tp@0 752 shouldplanebechecked = canplaneobstruct.*isplaneactiveforobstruction;
tp@0 753 nplanesperbatch = ceil(sum(shouldplanebechecked)/ndiff2batches);
tp@0 754 if nplanesperbatch > 0
tp@0 755 temp = cumsum(shouldplanebechecked);
tp@0 756 diff2batchlist = zeros(ndiff2batches,2);
tp@0 757 diff2batchlist(1,1) = 1;
tp@0 758 loopisfinished = 0;
tp@0 759 for ii = 1:ndiff2batches-1
tp@0 760 ivtemp = find(temp==nplanesperbatch*ii);
tp@0 761 if ~isempty(ivtemp)
tp@0 762 diff2batchlist(ii,2) = ivtemp(1);
tp@0 763 diff2batchlist(ii+1,1) = ivtemp(1)+1;
tp@0 764 else
tp@0 765 loopisfinished = 1;
tp@0 766 end
tp@0 767 end
tp@0 768 diff2batchlist(ii,2) = nplanes;
tp@0 769 ivtemp = find(diff2batchlist(:,1)>0);
tp@0 770 diff2batchlist = diff2batchlist(ivtemp,:);
tp@0 771 ndiff2batches = size(diff2batchlist,1);
tp@0 772
tp@0 773 npathstocheck = size(REFTOFROMCOSHO,1);
tp@0 774 for ii = 1:ndiff2batches
tp@0 775 if npathstocheck > 0
tp@0 776 maskedchkpla = shouldplanebechecked;
tp@0 777 if ii > 1
tp@0 778 maskedchkpla(1:diff2batchlist(ii,1)-1) = maskedchkpla(1:diff2batchlist(ii,1)-1)*0;
tp@0 779 end
tp@0 780 if ii < ndiff2batches
tp@0 781 maskedchkpla(diff2batchlist(ii,2)+1:end) = maskedchkpla(diff2batchlist(ii,2)+1:end)*0;
tp@0 782 end
tp@0 783 nplanestocheck = sum(maskedchkpla);
tp@0 784 if SHOWTEXT >= 3,
tp@0 785 if round(ii/1)*1==ii
tp@0 786 disp([' Batch ',int2str(ii),' of ',int2str(ndiff2batches),': '])
tp@0 787 disp([' ',int2str(npathstocheck),' paths to check, ',int2str(nplanestocheck),' planes to check obstruction for'])
tp@0 788 end
tp@0 789 end
tp@0 790
tp@0 791 % We create the BIGFROMCOORDS matrices and BIGTOCOORDS matrices
tp@0 792 % outside since they need to be exactly identical from batch to
tp@0 793 % batch as long as nplanestocheck is the same.
tp@0 794
tp@0 795 npathstocheck = size(REFTOFROMCOSHO,1);
tp@0 796
tp@0 797 ntot = nplanestocheck*npathstocheck;
tp@0 798
tp@0 799 BIGFROMCOORDS = reshape(repmat(FROMCOORDSSHORTLIST(REFTOFROMCOSHO,:).',[nplanestocheck,1]),3,ntot).';
tp@0 800 BIGTOCOORDS = reshape(repmat(TOCOORDSSHORTLIST(REFTOTOCOSHO,:).',[nplanestocheck,1]),3,ntot).';
tp@0 801 BIGSTARTPLANES = reshape(repmat(STARTPLANES.',[nplanestocheck,1]),2,ntot).';
tp@0 802 BIGENDPLANES = reshape(repmat(ENDPLANES.',[nplanestocheck,1]),2,ntot).';
tp@0 803
tp@0 804 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_edgetoedge(maskedchkpla,planeseesplane,...
tp@0 805 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 806
tp@0 807
tp@0 808 if ~isempty(edgehits)
tp@0 809 % Look for the edge-to-edge combos that share a
tp@0 810 % plane. When such combos have hit an edge, it
tp@0 811 % must mean that they are planes with indents.
tp@0 812 wasedgehitnotindent = prod(double(diff(sort([STARTPLANES(edgehits,:) ENDPLANES(edgehits,:)].')))).';
tp@0 813 indentedgehits = find(wasedgehitnotindent==0);
tp@0 814 if ~isempty(indentedgehits)
tp@0 815 nonobstructedpaths = setdiff(nonobstructedpaths,edgehits(indentedgehits));
tp@0 816 end
tp@0 817
tp@0 818 end
tp@0 819
tp@0 820 if ~isempty(nonobstructedpaths)
tp@0 821 REFTOFROMCOSHO = REFTOFROMCOSHO(nonobstructedpaths);
tp@0 822 REFTOTOCOSHO = REFTOTOCOSHO(nonobstructedpaths);
tp@0 823 STARTPLANES = STARTPLANES(nonobstructedpaths,:);
tp@0 824 ENDPLANES = ENDPLANES(nonobstructedpaths,:);
tp@0 825 bigfromedge = bigfromedge(nonobstructedpaths);
tp@0 826 bigtoedge = bigtoedge(nonobstructedpaths);
tp@0 827 npathstocheck = size(REFTOFROMCOSHO,1);
tp@0 828 else
tp@0 829 REFTOFROMCOSHO = [];
tp@0 830 REFTOTOCOSHO = [];
tp@0 831 STARTPLANES = [];
tp@0 832 ENDPLANES = [];
tp@0 833 bigfromedge = [];
tp@0 834 bigtoedge = [];
tp@0 835 npathstocheck = 0;
tp@0 836 end
tp@0 837
tp@0 838 end % if npathstocheck > 0
tp@0 839
tp@0 840 end % for ii = 1:ndiff2batches
tp@0 841 else
tp@0 842 nonobstructedpaths = 1;
tp@0 843
tp@0 844 end %if nplanesperbatch > 0
tp@0 845
tp@0 846 clear REFTOFROMCOSHO REFTOTOCOSHO STARTPLANES ENDPLANES
tp@0 847
tp@0 848
tp@0 849 if ~isempty(nonobstructedpaths)
tp@0 850
tp@0 851 % Add the segment pieces together
tp@0 852 % NB! We don't try to implement the correct way to do it since
tp@0 853 % we would need nedgesubs^2 bits to represent all edge-segment
tp@0 854 % to edge-segment visibility possibilities.
tp@0 855 % Instead we just establish whether the two edges see each
tp@0 856 % other fully or partly.
tp@0 857
tp@0 858 visiblesegmentscounter = ones(size(bigfromedge));
tp@0 859 test = [bigfromedge bigtoedge];
tp@0 860
tp@0 861 ncombs = length(bigfromedge);
tp@0 862 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 863 ivremove = find(dtest==1);
tp@0 864
tp@0 865 while ~isempty(ivremove)
tp@0 866 visiblesegmentscounter(ivremove+1) = visiblesegmentscounter(ivremove+1) + visiblesegmentscounter(ivremove);
tp@0 867 visiblesegmentscounter(ivremove) = [];
tp@0 868 bigfromedge(ivremove) = [];
tp@0 869 bigtoedge(ivremove,:) = [];
tp@0 870
tp@0 871 test = [bigfromedge bigtoedge];
tp@0 872 ncombs = length(bigfromedge);
tp@0 873 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 874 ivremove = find(dtest==1);
tp@0 875
tp@0 876 end
tp@0 877
tp@0 878 indexvec = uint32((double(bigfromedge)-1)*nedges + double(bigtoedge));
tp@0 879 if maxedgetoedgevisvalue < 128
tp@0 880 edgeseespartialedge = int8(zeros(nedges,nedges));
tp@0 881 elseif maxedgetoedgevisvalue < 32768
tp@0 882 edgeseespartialedge = int16(zeros(nedges,nedges));
tp@0 883 else
tp@0 884 edgeseespartialedge = int32(zeros(nedges,nedges));
tp@0 885 end
tp@0 886
tp@0 887 iv = find(visiblesegmentscounter==nedgesubs^2);
tp@0 888 edgeseespartialedge(indexvec(iv)) = maxedgetoedgevisvalue*ones(size(iv));
tp@0 889 iv = find(visiblesegmentscounter>0 & visiblesegmentscounter<nedgesubs^2);
tp@0 890 edgeseespartialedge(indexvec(iv)) = maxedgetoedgevisvalue/2*ones(size(iv));
tp@0 891
tp@0 892 if maxedgetoedgevisvalue < 128
tp@0 893 edgeseespartialedge = int8(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');
tp@0 894 elseif maxedgetoedgevisvalue < 32768
tp@0 895 edgeseespartialedge = int16(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');
tp@0 896 else
tp@0 897 edgeseespartialedge = int32(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');
tp@0 898 end
tp@0 899 else
tp@0 900 edgeseespartialedge = sparse(zeros(nedges,nedges));
tp@0 901 end
tp@0 902 else
tp@0 903 edgeseespartialedge = [];
tp@0 904 end
tp@0 905 else
tp@0 906 if maxedgetoedgevisvalue < 128
tp@0 907 edgeseespartialedge = int8((edgeseesedge*maxedgetoedgevisvalue));
tp@0 908 elseif maxedgetoedgevisvalue < 32768
tp@0 909 edgeseespartialedge = int16((edgeseesedge*maxedgetoedgevisvalue));
tp@0 910 else
tp@0 911 edgeseespartialedge = int32((edgeseesedge*maxedgetoedgevisvalue));
tp@0 912 end
tp@0 913 disp(' No obstruction tests needed!')
tp@0 914
tp@0 915 end %(obstructtestneeded)
tp@0 916
tp@0 917 % For the special case of thin planes that are all in the same plane, only
tp@0 918 % edge pairs that belong to the same plane should be included. That is
tp@0 919 % when edgeseespartialedge is negative, then we should set
tp@0 920 % edgeseespartialedge to zero.
tp@0 921
tp@0 922 % Bug detected 20110822 PS
tp@0 923 % Old: allthinplanes = sign(1-sum(~planeisthin));
tp@0 924 allthinplanes = sign(1-sign(sum(~planeisthin)));
tp@0 925
tp@0 926 % Bug detected 20110822 PS
tp@0 927 % Old: allplanesinsameplane = 1 - sum(sum(planeseesplane==1))
tp@0 928 allplanesinsameplane = 1 - sign(sum(sum(planeseesplane==1)));
tp@0 929
tp@0 930 multfactor = 1 - allthinplanes*allplanesinsameplane;
tp@0 931
tp@0 932 iv = find(edgetoedgestrength==2);
tp@0 933 if ~isempty(iv)
tp@0 934 if maxedgetoedgevisvalue < 128
tp@0 935 edgeseespartialedge(iv) = int8(-double(edgeseespartialedge(iv))*multfactor);
tp@0 936 elseif maxedgetoedgevisvalue < 32768
tp@0 937 edgeseespartialedge(iv) = int16(-double(edgeseespartialedge(iv))*multfactor);
tp@0 938 else
tp@0 939 edgeseespartialedge(iv) = int32(-double(edgeseespartialedge(iv))*multfactor);
tp@0 940 end
tp@0 941 end
tp@0 942
tp@0 943 clear edgeseesedge edgetoedgestrength
tp@0 944
tp@0 945 %----------------------------------------------------------------------------
tp@0 946 %
tp@0 947 % CYLINDRICAL EDGE-TO-EDGE PARAMETERS
tp@0 948 %
tp@0 949 %----------------------------------------------------------------------------
tp@0 950 %
tp@0 951 % For each edge, calculate the cylindrical coordinates of the starting
tp@0 952 % and ending points of all other edges relative to the first edge.
tp@0 953
tp@0 954 if difforder >= 2 & ~isempty(edgeseespartialedge)
tp@0 955
tp@0 956 if SHOWTEXT >= 4
tp@0 957 disp(' E2E')
tp@0 958 end
tp@0 959
tp@0 960 zerosvec4 = zeros(nedges,nedges);
tp@0 961 Bigre1 = (zerosvec4);
tp@0 962 Bigthetae1 = (zerosvec4);
tp@0 963 Bigze1 = (zerosvec4);
tp@0 964 Bigre2 = (zerosvec4);
tp@0 965 Bigthetae2 = (zerosvec4);
tp@0 966 Bigze2 = (zerosvec4);
tp@0 967 clear zerosvec4
tp@0 968
tp@0 969 for edge1 = 1:nedges
tp@0 970 if SHOWTEXT >= 4
tp@0 971 if round(edge1/1)*1 == edge1
tp@0 972 disp([' Edge no. ',int2str(edge1)])
tp@0 973 end
tp@0 974 end
tp@0 975 edge1coords = [edgestartcoords(edge1,:);edgeendcoords(edge1,:)];
tp@0 976 iv = find(edgeseespartialedge(edge1,:)~=0).';
tp@0 977
tp@0 978 if ~isempty(iv),
tp@0 979 % First find the subset of edges which belong to the same plane as the
tp@0 980 % edge itself. These should be treated separately for higher accuracy
tp@0 981
tp@0 982 % iv1 will be the edges on the reference plane. They should have theta = 0.
tp@0 983
tp@0 984
tp@0 985 refplane = planesatedge(edge1,1);
tp@0 986 ncornersperplanevec = double(ncornersperplanevec);
tp@0 987 iv1 = edgesatplane( refplane,1:ncornersperplanevec( refplane ));
tp@0 988 iv1 = iv1( find(iv1~= edge1)).';
tp@0 989
tp@0 990 edgestart = edgestartcoordsnudge(iv1,:);
tp@0 991 edgeend = edgeendcoordsnudge(iv1,:);
tp@0 992
tp@0 993 [rs,slask,zs,rr,slask,zr] = ...
tp@0 994 EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0 995 Bigre1(iv1,edge1) = rs;
tp@0 996 Bigthetae1(iv1,edge1) = 0*iv1;
tp@0 997 Bigze1(iv1,edge1) = zs;
tp@0 998 Bigre2(iv1,edge1) = rr;
tp@0 999 Bigthetae2(iv1,edge1) = 0*iv1;
tp@0 1000 Bigze2(iv1,edge1) = zr;
tp@0 1001
tp@0 1002 [samevalues,iivec,jjvec] = intersect(iv,iv1);
tp@0 1003 % Bug found 20050421!! PS
tp@0 1004 % % % Old wrong version:
tp@0 1005 % % % % sameedges = [iivec jjvec];
tp@0 1006 % % % if sum(sum(sameedges)) ~= 0
tp@0 1007 % % % iv( sameedges(1,:).' ) = [];
tp@0 1008 % % % end
tp@0 1009 if ~isempty(samevalues)
tp@0 1010 iv(iivec) = [];
tp@0 1011 end
tp@0 1012
tp@0 1013 % iv2 will be the edges on the non-reference plane.
tp@0 1014 % They should have theta = 2*pi - closthetavec.
tp@0 1015
tp@0 1016 if ~isempty(iv)
tp@0 1017
tp@0 1018 if planesatedge(edge1,2) > 0
tp@0 1019 secplane = planesatedge(edge1,2);
tp@0 1020
tp@0 1021 iv2 = edgesatplane( secplane,1:ncornersperplanevec(secplane));
tp@0 1022 iv2 = iv2( find(iv2~= edge1)).';
tp@0 1023
tp@0 1024 if ~isempty(iv2)
tp@0 1025 edgestart = edgestartcoordsnudge(iv2,:);
tp@0 1026 edgeend = edgeendcoordsnudge(iv2,:);
tp@0 1027
tp@0 1028 [rs,slask,zs,rr,slask,zr] = ...
tp@0 1029 EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0 1030 Bigre1(iv2,edge1) = rs;
tp@0 1031 Bigthetae1(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
tp@0 1032 Bigze1(iv2,edge1) = zs;
tp@0 1033 Bigre2(iv2,edge1) = rr;
tp@0 1034 Bigthetae2(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
tp@0 1035 Bigze2(iv2,edge1) = zr;
tp@0 1036
tp@0 1037 % sameedges = EDB1findsame(iv,iv2);
tp@0 1038 % Bug found 20050421!!! PS
tp@0 1039 % % % % Old wrong version:
tp@0 1040 % % % % [samevalues,iivec,jjvec] = intersect(iv,iv1);
tp@0 1041 % % % % sameedges = [iivec;jjvec];
tp@0 1042 % % % % if sum(sum(sameedges)) ~= 0
tp@0 1043 % % % % iv( sameedges(1,:).' ) = [];
tp@0 1044 % % % % end
tp@0 1045
tp@0 1046 [samevalues,iivec,jjvec] = intersect(iv,iv2);
tp@0 1047 if ~isempty(samevalues)
tp@0 1048 iv(iivec) = [];
tp@0 1049 end
tp@0 1050 end
tp@0 1051 end
tp@0 1052 end
tp@0 1053 if ~isempty(iv)
tp@0 1054
tp@0 1055 % Move the edge coordinates to be checked a short distance away
tp@0 1056
tp@0 1057 edgestart = edgestartcoordsnudge(iv,:);
tp@0 1058 edgeend = edgeendcoordsnudge(iv,:);
tp@0 1059
tp@0 1060 [rs,thetas,zs,rr,thetar,zr] = ...
tp@0 1061 EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0 1062 Bigre1(iv,edge1) = rs;
tp@0 1063 Bigthetae1(iv,edge1) = thetas;
tp@0 1064 Bigze1(iv,edge1) = zs;
tp@0 1065 Bigre2(iv,edge1) = rr;
tp@0 1066 Bigthetae2(iv,edge1) = thetar;
tp@0 1067 Bigze2(iv,edge1) = zr;
tp@0 1068
tp@0 1069 end
tp@0 1070 end
tp@0 1071 end
tp@0 1072
tp@0 1073 %-----------------------------------------------------------
tp@0 1074 % Go through all edges that are in-plane with each other, that is, two
tp@0 1075 % edges have at least one plane each that are in-plane with each other.
tp@0 1076
tp@0 1077 % First we identify all planes that have at least one more co-planar
tp@0 1078 % plane
tp@0 1079
tp@0 1080 planehascoplanar = any(planeseesplane == -1);
tp@0 1081
tp@0 1082 % Then we go through the list of edges and every edge with a connected
tp@0 1083 % plane that has another co-planar plane must also have potential
tp@0 1084 % in-plane edges.
tp@0 1085
tp@0 1086 % For each edge, we check if there are other edges (that don't belong to the same plane!!)
tp@0 1087 % with thetaangle = 0 or thetaangle = 2*pi. If there are other such
tp@0 1088 % edges, those edge-to-edge combinations should be shut off.
tp@0 1089 %
tp@0 1090 % After all such edge-to-edge combinations have been shut off, we still
tp@0 1091 % need another pass to cancel edge-to-edge paths that pass entirely
tp@0 1092 % across other edges.
tp@0 1093
tp@0 1094 ivedges = 1:nedges;
tp@0 1095 ivedges(offedges) = [];
tp@0 1096
tp@0 1097 for ii = ivedges
tp@0 1098 if any( planehascoplanar(planesatedge(1,:)) )
tp@0 1099
tp@0 1100 ivcancelcombs = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ( (Bigthetae1(:,ii) == 0) | (Bigthetae1(:,ii) == 2*pi) ) );
tp@0 1101 if ~isempty(ivcancelcombs)
tp@0 1102 edgeseespartialedge(ivcancelcombs,ii) = 0;
tp@0 1103 edgeseespartialedge(ii,ivcancelcombs.') = 0;
tp@0 1104 end
tp@0 1105
tp@0 1106 end
tp@0 1107
tp@0 1108 end
tp@0 1109
tp@0 1110 % The only edge-to-edge combinations that we can easily discard are
tp@0 1111 % those for which edges are aligned with each other (same zstart and
tp@0 1112 % zend): only the closest of those should be kept!
tp@0 1113
tp@0 1114 for ii = ivedges
tp@0 1115 if any( planehascoplanar(planesatedge(1,:)) )
tp@0 1116
tp@0 1117
tp@0 1118 ivotheredges = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ((Bigthetae1(:,ii) == pi)) );
tp@0 1119 if ~isempty(ivotheredges)
tp@0 1120
tp@0 1121 nudgeval = 1e-10*edgelengthvec(ii)*100;
tp@0 1122 ivselectedges = find( abs(Bigze1(ivotheredges,ii))<nudgeval | abs(Bigze2(ivotheredges,ii))<nudgeval);
tp@0 1123 if length(ivselectedges) > 1
tp@0 1124 disp(['Edge no. ',int2str(ii)])
tp@0 1125
tp@0 1126 lengthok = abs([Bigze1(ivotheredges(ivselectedges),ii) Bigze2(ivotheredges(ivselectedges),ii)] - edgelengthvec(ii)) < nudgeval;
tp@0 1127 ivstillok = find(any(lengthok.'));
tp@0 1128
tp@0 1129 if ~isempty(ivstillok)
tp@0 1130 ivotheredges = ivotheredges(ivselectedges(ivstillok));
tp@0 1131 else
tp@0 1132 ivotheredges = [];
tp@0 1133 end
tp@0 1134
tp@0 1135 meanradialdist = mean( [Bigre1(ivotheredges,ii) Bigre2(ivotheredges,ii)].' );
tp@0 1136
tp@0 1137 ivshutoff = ivotheredges;
tp@0 1138 noshutoff = find(meanradialdist == min(meanradialdist));
tp@0 1139 ivshutoff(noshutoff) = [];
tp@0 1140
tp@0 1141 if ~isempty(ivshutoff)
tp@0 1142 edgeseespartialedge(ivshutoff,ii) = 0;
tp@0 1143 edgeseespartialedge(ii,ivshutoff) = 0;
tp@0 1144 end
tp@0 1145
tp@0 1146 end
tp@0 1147
tp@0 1148 end
tp@0 1149
tp@0 1150 end
tp@0 1151
tp@0 1152 end
tp@0 1153
tp@0 1154 %-----------------------------------------------------------------------
tp@0 1155 % Find identical combinations
tp@0 1156
tp@0 1157 if SHOWTEXT >= 3
tp@0 1158 disp(' Looking for identical combinations')
tp@0 1159 end
tp@0 1160 [reftoshortlistE,re1sho,re2sho,thetae1sho,thetae2sho,...
tp@0 1161 ze1sho,ze2sho,edgelengthsho,examplecombE] = EDB1compress7(Bigre1,Bigre2,...
tp@0 1162 Bigthetae1,Bigthetae2,Bigze1,Bigze2, edgelengthvec(:,ones(1,nedges)).');
tp@0 1163
tp@0 1164 else
tp@0 1165 reftoshortlistE = [];
tp@0 1166 re1sho = [];
tp@0 1167 re2sho = [];
tp@0 1168 thetae1sho = [];
tp@0 1169 thetae2sho = [];
tp@0 1170 ze1sho = [];
tp@0 1171 ze2sho = [];
tp@0 1172 examplecombE = [];
tp@0 1173 end
tp@0 1174
tp@0 1175 %----------------------------------------------------------------------------
tp@0 1176 %
tp@0 1177 % SAVE THE VARIABLES
tp@0 1178 %
tp@0 1179 %----------------------------------------------------------------------------
tp@0 1180 %
tp@0 1181 % Also the variables from cadgeo??
tp@0 1182 % Yes: planeeqs planenvecs ncornersperplanevec planecorners corners
tp@0 1183 % minvals maxvals
tp@0 1184
tp@0 1185 edgeseesplane = int8(edgeseesplane);
tp@0 1186 planeseesplane = int8(planeseesplane);
tp@0 1187
tp@0 1188 ncorners = size(corners,1);
tp@0 1189 if ncorners < 256
tp@0 1190 planecorners = uint8(planecorners);
tp@0 1191 elseif ncorners < 65536
tp@0 1192 planecorners = uint16(planecorners);
tp@0 1193 end
tp@0 1194 planeisthin = uint8(planeisthin);
tp@0 1195
tp@0 1196 Varlist = [ ' edgecorners planesatedge closwedangvec planehasindents'];
tp@0 1197 Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
tp@0 1198 Varlist = [Varlist,' corners planecorners planeeqs planenvecs ncornersperplanevec'];
tp@0 1199 Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
tp@0 1200 Varlist = [Varlist,' minvals maxvals offedges edgestartcoordsnudge edgeendcoordsnudge'];
tp@0 1201 Varlist = [Varlist,' reftoshortlistE re1sho re2sho thetae1sho thetae2sho ze1sho ze2sho examplecombE'];
tp@0 1202 Varlist = [Varlist,' edgeseespartialedge int_or_ext_model'];
tp@0 1203 Varlist = [Varlist,' edgealignedwithedge edgeperptoplane edgeplaneperptoplane1 edgeplaneperptoplane2'];
tp@0 1204
tp@0 1205 eval(['save ',outputfile,Varlist])