tp@0: function outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches)
tp@0: % EDB1ed2geox - Calculates 2nd- and higher-order edge-related geom. parameters.
tp@0: % EDB1ed2geox calculates some plane- and edge-related geometrical parameters
tp@0: % based on corners and planes in a .mat-file generated by EDB1readcad
tp@0: % but only data that is independent of the source and receiver.
tp@0: % The output is saved in a .mat-file.
tp@0: %
tp@0: % Input parameters:
tp@0: % eddatafile The .mat file where the output data will be stored.
tp@0: % If it is not specified, an automatic file name is constructed
tp@0: % using the same filename stem as the cadgeofile, but ending
tp@0: % with '_edgeo'.
tp@0: % sdatafile
tp@0: % rdatafile
tp@0: % doSbatches
tp@0: % doRbatches
tp@0: % specorder The highest order of any reflection kind (specular and/or diffraction).
tp@0: % difforder The highest order of diffraction. If it is 0 or 1 then the parameter
tp@0: % edgeseespartialedge is not calculated. Default: 1
tp@0: % nedgesubs
tp@0: % ndiff2batches
tp@0: % SHOWTEXT (global) An integer value; if it is 4 or higher, then messages will be printed on the
tp@0: % screen during calculations.
tp@0: % Output parameters:
tp@0: % outputfile The name of the outputfile, generated either automatically or specified as
tp@0: % an input parameter.
tp@0: %
tp@0: % Output parameters in the outputfile:
tp@0: % (Taken directly from the cadgeo-file:)
tp@0: % corners planecorners planeeqs planenvecs ncornersperplanevec
tp@0: % minvals maxvals
tp@0: % (Taken directly from the setup-file:)
tp@0: % int_or_ext_model
tp@0: % (Calculated by EDB1edgeo:)
tp@0: % edgecorners Matrix [nedges,2] with the corner numbers that
tp@0: % define each edge, in the reference order.
tp@0: % edgestartcoords Matrix [nedges,3] with the coordinates of the
tp@0: % start corner for each edge.
tp@0: % edgeendcoords Matrix [nedges,3] with the coordinates of the
tp@0: % end corner for each edge.
tp@0: % edgelengthvec List [nedges,1] with the lengths of each edge.
tp@0: % planesatedge Matrix [nedges,2] with the plane numbers of the
tp@0: % two planes that connect to each edge. The first
tp@0: % plane is considered as the reference plane for
tp@0: % each edge. Rule: if the RH thumb is placed
tp@0: % along the edge, from start to end, and the hand
tp@0: % is rotated, the fingers should "come up
tp@0: % through" the reference plane.
tp@0: % closwedangvec List [nedges,1] with the angles in radians of
tp@0: % each edge - for the closed part of each edge.
tp@0: % edgenvecs Matrix [nedges,3] with the normal vectors of
tp@0: % the reference plane for each edge.
tp@0: % offedges List [nlength,1] where nlength = 0-nedges
tp@0: % containing the edge numbers that should not
tp@0: % used.
tp@0: % reflfactors List [nplanes,1] with the reflection factors
tp@0: % of each plane. The only supported values are
tp@0: % +1 for rigid planes, 0 for perfectlt absorbing
tp@0: % planes and -1 for pressure-release planes. The
tp@0: % pressure-release cases might not be implemented
tp@0: % fully.
tp@0: % planeisthin List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether the plane is thin or not.
tp@0: % rearsideplane List [nplanes,1] with the plane number that is at
tp@0: % the rear side. If a plane is non-thin, the
tp@0: % value in this list is zero.
tp@0: % planehasindents List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether a plane has indents or not.
tp@0: % canplaneobstruct List [nplanes,1] with the value 1 or 0 stating
tp@0: % whether a plane has the potential to obstruct
tp@0: % or not.
tp@0: % planeseesplane A matrix, [nplanes,nplanes], with the value 1
tp@0: % 0,-1,-2 stating:
tp@0: % 1 A plane is front of another plane and could
tp@0: % then potentially see it.
tp@0: % 0 A plane is completely behind another plane
tp@0: % but they are not aligned
tp@0: % -1 A plane is aligned with another plane
tp@0: % -2 A plane is back-to-back with another plane.
tp@0: % edgesatplane A matrix, [nplanes,nmaxedgesperplane] which for
tp@0: % row N contains the edge numbers that connect to
tp@0: % to plane N.
tp@0: % edgeseesplane A matrix, [nplanes,nedges], which contains the
tp@0: % values 0,1,-2 or -2, meaning:
tp@0: % 0 The edge can never see the plane
tp@0: % -1 The edge can never see the plane and the
tp@0: % edge is aligned with the plane
tp@0: % -2 The edge can never see the plane and the
tp@0: % edge belongs to the plane
tp@0: % 1 The edge can potentially see the plane.
tp@0: %
tp@0: % (Calculated by EDB1edgeo, but only for diffraction order >= 2:)
tp@0: % edgeseespartialedge A sparse matrix, [nedges,nedges], with the
tp@0: % values 0, pos. integer or neg. integer:
tp@0: % 0 Two edges can not see each other (or one
tp@0: % of the edges is inactive)
tp@0: % pos.int. A value from 1 to 2^nedgesubs-1
tp@0: % indicating how much of the edges
tp@0: % see each other. A pos. value
tp@0: % indicates that the edge-to-edge
tp@0: % path runs along a plane.
tp@0: % NB! A complete visibility test is
tp@0: % not made; subsegment 1 on edge 1 is only
tp@0: % checked against subsegment 1 on
tp@0: % edge 2, not against all other
tp@0: % subsegments of edge 2.
tp@0: % neg.int. Same as for the pos.int. except
tp@0: % that the edge-to-edge path does not
tp@0: % run along a plane.
tp@0: %
tp@0: % reftoshortlistE
tp@0: % re1sho, re2sho
tp@0: % thetae1sho, thetae2sho
tp@0: % ze1sho, ze2sho
tp@0: % examplecombE
tp@0: % edgealignedwithedge A sparse matrix, [nedges,nedges], with the value 1
tp@0: % or 0 stating whether an (active) edge is aligned
tp@0: % with another (active) edge or not.
tp@0: % edgeperptoplane A sparse matrix, [nplanes,nedges], with the value 1
tp@0: % or 0 stating whether an (active) edge is perpendicular
tp@0: % to an active plane or not.
tp@0: % edgeplaneperptoplane1 A sparse matrix, [nplanes,nedges], with the value 1
tp@0: % or 0 stating whether an (active) edge has one
tp@0: % of its two defining planes perpendicular to
tp@0: % another plane with which it shares an edge.
tp@0: % edgeplaneperptoplane2 A sparse matrix, [nplanes,nedges], with the value 1
tp@0: % or 0 stating whether an (active) edge has one
tp@0: % of its two defining planes perpendicular to
tp@0: % another plane which has a flat edge (i.e., a 180 degree edge).
tp@0: %
tp@0: % Uses the functions EDB1coordtrans1, EDB1coordtrans2
tp@0: % EDB1infrontofplane, EDB1compress7, EDB1strpend, EDB1cross
tp@0: % EDB1checkobstr_edgetoedge
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0: %
tp@0: % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0: % it under the terms of the GNU General Public License as published by the Free Software
tp@0: % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0: %
tp@0: % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0: % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0: % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0: %
tp@0: % You should have received a copy of the GNU General Public License along with the
tp@0: % Edge Diffraction Toolbox. If not, see .
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20110822
tp@0: %
tp@0: % outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches);
tp@0:
tp@0: % 20110822 Bug found and fixed, giving the wrong amplitude for second- and higher-order
tp@0: % diffraction, but only for diffraction paths that do not travel along a surface.
tp@0:
tp@0: global SHOWTEXT
tp@0:
tp@0: geomacc = 1e-10;
tp@0:
tp@0: eval(['load ',eddatafile])
tp@0: edgestartcoordsnudge = edgedata.edgestartcoordsnudge;
tp@0: edgeendcoordsnudge = edgedata.edgeendcoordsnudge;
tp@0: edgenormvecs = edgedata.edgenormvecs;
tp@0:
tp@0: [sfilepath,sfilestem,temp1] = fileparts(sdatafile);
tp@0: sdatafile = [sfilepath,filesep,sfilestem];
tp@0: if doSbatches == 0
tp@0: eval(['load ',sdatafile,'.mat'])
tp@0: else
tp@0: eval(['load ',sdatafile,'_1.mat'])
tp@0: end
tp@0: [rfilepath,rfilestem,temp1] = fileparts(rdatafile);
tp@0: rdatafile = [rfilepath,filesep,rfilestem];
tp@0: if doRbatches == 0
tp@0: eval(['load ',rdatafile,'.mat'])
tp@0: else
tp@0: eval(['load ',rdatafile,'_1.mat'])
tp@0: end
tp@0:
tp@0: outputfile = eddatafile;
tp@0:
tp@0: nplanes = size(planecorners,1);
tp@0: nedges = size(edgecorners,1);
tp@0:
tp@0: zerosvec1 = zeros(nplanes,1);
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: % We check right away which edges that can not be seen by either the source
tp@0: % or the receiver. We don't need to check edgeseesedge for combinations
tp@0: % that involve any of these.
tp@0: %
tp@0: % Note that this is true only for difforder = 2! If difforder >= 3, then we
tp@0: % can not exclude any edges based on this! (It could be possible to exclude
tp@0: % edges, if we managed to "unwind" the
tp@0: % edge-sees-an-edge-which-sees-an-edge-which-can-be-seen-by-the-source.
tp@0:
tp@0: if difforder == 2
tp@0: edgesthatareseen = any( [vispartedgesfroms vispartedgesfromr].' );
tp@0: edgesnottoworryabout = [1:nedges];
tp@0: edgesnottoworryabout(edgesthatareseen) = [];
tp@0: clear vispartedgesfroms vispartedgesfromr visplanesfroms visplanesfromr
tp@0: else
tp@0: edgesnottoworryabout = [];
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % edgeperptoplane
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0: % Make matrices over:
tp@0: % 1. Which edges are perpendicular to planes ("edgeperptoplane")
tp@0: % 2. Which edges can give combinations with edge1-plane-edge1
tp@0: % where plane is perpendicular to one of the two planes that
tp@0: % defines the edge ("edgeplaneperptoplane1")
tp@0: % 3. Which edges can give combinations with edge1-plane-edge1
tp@0: % where plane is a part of a (unnecessarily) divided plane
tp@0: % ("edgeplaneperptoplane2")
tp@0: %
tp@0: % NB! Only active edges and active planes can get
tp@0: % such combinations.
tp@0: %
tp@0: % These are needed for diff-spec-diff combos, so only when
tp@0: % difforder >= 2 (and specorder >= 3).
tp@0:
tp@0: edgeperptoplane = [];
tp@0: edgeplaneperptoplane1 = [];
tp@0: edgeplaneperptoplane2 = [];
tp@0: if specorder >= 3
tp@0: edgeperptoplane = sparse(zeros(nplanes,nedges));
tp@0: edgeplaneperptoplane1 = sparse(zeros(nplanes,nedges));
tp@0: edgeplaneperptoplane2 = sparse(zeros(nplanes,nedges));
tp@0: if SHOWTEXT >= 4
tp@0: disp(' checking which edges are perpendicular to planes ...')
tp@0: end
tp@0:
tp@0: iv = find(edgeseesplane == 1);
tp@0: if ~isempty(iv)
tp@0: edgelist = floor((iv-1)/nplanes)+1;
tp@0: planelist = iv - (edgelist-1)*nplanes;
tp@0:
tp@0: % Check first which edges have vectors that are parallel to plane
tp@0: % normal vectors. These are "edgeperptoplane" cases. Clear them
tp@0: % afterwards from the edgelist and planelist.
tp@0:
tp@0: vec1 = edgenormvecs(edgelist,:);
tp@0: vec2 = planenvecs(planelist,:);
tp@0: ivperp = find( abs( sum( (vec1.*vec2).' ) )==1);
tp@0: if ~isempty(ivperp)
tp@0: edgeperptoplane(iv(ivperp)) = ones(size(ivperp));
tp@0: clear iv
tp@0: edgelist(ivperp) = [];
tp@0: planelist(ivperp) = [];
tp@0: vec2(ivperp,:) = [];
tp@0: else
tp@0: clear iv
tp@0: end
tp@0:
tp@0: % Now check which edges have plane1 or plane2 perpendicular to
tp@0: % other planes.
tp@0:
tp@0: if ~isempty(planelist)
tp@0:
tp@0: vec1 = planenvecs(planesatedge(edgelist,1),:);
tp@0: ivplaneperp1 = find( abs( sum( (vec1.*vec2).' ) )==0);
tp@0: vec1 = planenvecs(planesatedge(edgelist,2),:);
tp@0: ivplaneperp2 = find( abs( sum( (vec1.*vec2).' ) )==0);
tp@0: ivplaneperp = [ivplaneperp1 ivplaneperp2];
tp@0: if ~isempty(ivplaneperp)
tp@0: ivplaneperp = unique(ivplaneperp);
tp@0: edgelist = edgelist(ivplaneperp);
tp@0: planelist = planelist(ivplaneperp);
tp@0: edgeplane1 = planesatedge(edgelist,1);
tp@0: edgeplane2 = planesatedge(edgelist,2);
tp@0:
tp@0: % Check which of these other perpendicular planes that connnect
tp@0: % to one of the edge planes via a 90 deg. corner.
tp@0:
tp@0: ivtoclear = [];
tp@0: [ivfoundcombos,loc] = ismember([planelist edgeplane1],planesatedge,'rows');
tp@0: ivperp = find(ivfoundcombos);
tp@0: if ~isempty(ivperp)
tp@0: edgenumb = loc(ivperp);
tp@0: ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0: if ~isempty(ivfinalcomb)
tp@0: ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0: edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0: ivtoclear = ivfinalcomb;
tp@0: end
tp@0: end
tp@0: [ivfoundcombos,loc] = ismember([edgeplane1 planelist],planesatedge,'rows');
tp@0: ivperp = find(ivfoundcombos);
tp@0: if ~isempty(ivperp)
tp@0: edgenumb = loc(ivperp);
tp@0: ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0: if ~isempty(ivfinalcomb)
tp@0: ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0: edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0:
tp@0: ivtoclear = ivtoclear.';
tp@0: ivfinalcomb = ivfinalcomb.';
tp@0: ivtoclear = [ivtoclear ivfinalcomb];
tp@0: end
tp@0: end
tp@0:
tp@0: [ivfoundcombos,loc] = ismember([planelist edgeplane2],planesatedge,'rows');
tp@0: ivperp = find(ivfoundcombos);
tp@0: if ~isempty(ivperp)
tp@0: edgenumb = loc(ivperp);
tp@0: ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0: if ~isempty(ivfinalcomb)
tp@0: ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0: edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0:
tp@0: if size(ivtoclear,2) == size(ivfinalcomb,2)
tp@0: ivtoclear = [ivtoclear;ivfinalcomb];
tp@0: elseif size(ivtoclear,1) == size(ivfinalcomb,1),
tp@0: ivfinalcomb = ivfinalcomb.';
tp@0: ivtoclear = [ivtoclear ivfinalcomb];
tp@0: elseif isempty(ivtoclear)
tp@0: error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0:
tp@0: [ivfoundcombos,loc] = ismember([edgeplane2 planelist],planesatedge,'rows');
tp@0: ivperp = find(ivfoundcombos);
tp@0: if ~isempty(ivperp)
tp@0: edgenumb = loc(ivperp);
tp@0: ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
tp@0: if ~isempty(ivfinalcomb)
tp@0: ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
tp@0: edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0:
tp@0: if size(ivtoclear,2) == size(ivfinalcomb,2)
tp@0: ivtoclear = [ivtoclear;ivfinalcomb];
tp@0: elseif size(ivtoclear,1) == size(ivfinalcomb,1),
tp@0: ivfinalcomb = ivfinalcomb.';
tp@0: ivtoclear = [ivtoclear ivfinalcomb];
tp@0: elseif isempty(ivtoclear)
tp@0: error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: % Finally check which edges have plane1 or plane2 perpendicular to
tp@0: % two other planes that connnect via a 180 deg. corner.
tp@0:
tp@0: if ~isempty(ivtoclear)
tp@0: ivtoclear = unique(ivtoclear);
tp@0: edgelist(ivtoclear) = [];
tp@0: planelist(ivtoclear) = [];
tp@0: if ~isempty(edgelist)
tp@0:
tp@0: ivflatedges = find(closwedangvec==pi);
tp@0: if ~isempty(ivflatedges)
tp@0: planecombos = planesatedge(ivflatedges,:);
tp@0: for ii = 1:length(ivflatedges)
tp@0: iv1 = find(planelist == planecombos(ii,1));
tp@0: if ~isempty(iv1)
tp@0: ivreftomatrix = double(planecombos(ii,1)) + (edgelist(iv1)-1)*nplanes;
tp@0: edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0: end
tp@0: iv2 = find(planelist == planecombos(ii,2));
tp@0: if ~isempty(iv2)
tp@0: ivreftomatrix = planecombos(ii,2) + (edgelist(iv2)-1)*nplanes;
tp@0: edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0: end
tp@0:
tp@0: end
tp@0: end
tp@0:
tp@0: end
tp@0: end
tp@0:
tp@0: end
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % edgealignedwithedge
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0:
tp@0: edgealignedwithedge = [];
tp@0:
tp@0: if specorder >= 2
tp@0: edgealignedwithedge = sparse(eye(nedges));
tp@0: if SHOWTEXT >= 4
tp@0: disp(' checking which edges are aligned with other edges ...')
tp@0: end
tp@0:
tp@0: listofactiveedges = [1:nedges].';
tp@0: listofactiveedges(offedges) = [];
tp@0: vecstocheck = edgenormvecs(listofactiveedges,:);
tp@0: iv = find(vecstocheck(:,1)<0);
tp@0: vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0: iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)<0);
tp@0: vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0: iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)==0 & vecstocheck(:,3)<0);
tp@0: vecstocheck(iv,:) = -vecstocheck(iv,:);
tp@0: if ~isempty(listofactiveedges)
tp@0: [uniquevecs,II,JJ] = unique(vecstocheck,'rows');
tp@0: N = histc(JJ,[1:length(listofactiveedges)]);
tp@0: iv = find(N>1);
tp@0: for ii = 1:length(iv)
tp@0: iv2 = find(JJ==iv(ii));
tp@0: edgestocheck = listofactiveedges(iv2);
tp@0: ntocheck = length(edgestocheck);
tp@0: cornernumberstocheck = edgecorners(edgestocheck,:);
tp@0: expandedstartcornermatrix = [reshape(repmat(cornernumberstocheck(:,1).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,1),[ntocheck 1])];
tp@0: expandedendcornermatrix = [reshape(repmat(cornernumberstocheck(:,2).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,2),[ntocheck 1])];
tp@0: toexpandededgestocheck = repmat(edgestocheck, [ntocheck 1]);
tp@0: fromexpandededgestocheck = reshape(repmat(edgestocheck.',[ntocheck 1]),ntocheck^2,1);
tp@0:
tp@0: % Now we don't need to check edges against themselves, or edge1
tp@0: % vs edge2 if edge2 vs edge1 has already been tested.
tp@0:
tp@0: indmat = reshape([1:ntocheck^2],ntocheck,ntocheck);
tp@0: indmat = triu(indmat,1);
tp@0: indmat = indmat(find(indmat));
tp@0:
tp@0: expandedstartcornermatrix = expandedstartcornermatrix(indmat,:);
tp@0: expandedendcornermatrix = expandedendcornermatrix(indmat,:);
tp@0: fromexpandededgestocheck = fromexpandededgestocheck(indmat,:);
tp@0: toexpandededgestocheck = toexpandededgestocheck(indmat,:);
tp@0:
tp@0: % We make the easiest check first: if two edges with the same
tp@0: % direction vector share one corner, they must be aligned.
tp@0:
tp@0: A = [expandedstartcornermatrix expandedendcornermatrix];
tp@0: A = sort(A.').';
tp@0: A = diff(A.').';
tp@0: iv3 = find(sum(A.'==0).');
tp@0: validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
tp@0: if ~isempty(validcombs)
tp@0: indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
tp@0: edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0: indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
tp@0: edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0: fromexpandededgestocheck(iv3) = [];
tp@0: toexpandededgestocheck(iv3) = [];
tp@0: expandedstartcornermatrix(iv3,:) = [];
tp@0: expandedendcornermatrix(iv3,:) = [];
tp@0: end
tp@0:
tp@0: % For the other edge-pair combinations we must check if the two
tp@0: % edges are really aligned. We do this by checking if two
tp@0: % vectors are aligned: edge1-vector and the vector from
tp@0: % edge1start to edge2start.
tp@0:
tp@0: direcvec = corners(expandedstartcornermatrix(:,2),:) - corners(expandedstartcornermatrix(:,1),:);
tp@0: direcveclengths = sqrt( sum(direcvec.'.^2) ).';
tp@0: direcvec = direcvec./(direcveclengths(:,ones(1,3)) + eps*10);
tp@0: test = abs(sum( (direcvec.*edgenormvecs(fromexpandededgestocheck,:)).' ));
tp@0: iv3 = find(abs(test-1)< 1e-8);
tp@0: validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
tp@0: if ~isempty(validcombs)
tp@0: indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
tp@0: edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0: indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
tp@0: edgealignedwithedge(indexvals) = edgealignedwithedge(indexvals) + 1;
tp@0: fromexpandededgestocheck(iv3) = [];
tp@0: toexpandededgestocheck(iv3) = [];
tp@0: expandedstartcornermatrix(iv3,:) = [];
tp@0: expandedendcornermatrix(iv3,:) = [];
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % edgeseesedge
tp@0: %
tp@0: %---------------------------------------------------------------
tp@0: %
tp@0: % Make a matrix of which edges can see which edges
tp@0: %
tp@0: % Planes that are in front of the edge-planes might have visible edges.
tp@0: % If closwedang = 0, then all other edges are potentially visible.
tp@0: % If closwedang < pi, then all planes in front of plane1 or plane2 are OK
tp@0: % If closwedang > pi, then all planes in front of plane1 and plane2 are OK
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' checking which edges see which edges...')
tp@0: end
tp@0:
tp@0: edgeseesedge = int8(zeros(nedges,nedges));
tp@0: for ii = 1:nedges
tp@0: closwed = closwedangvec(ii);
tp@0: if closwed == 0
tp@0: edgeseesedge(:,ii) = ones(nedges,1);
tp@0: else
tp@0: plane1 = planesatedge(ii,1);
tp@0: plane2 = planesatedge(ii,2);
tp@0: okplanelist1 = find(planeseesplane(:,plane1)==1);
tp@0: okplanelist2 = find(planeseesplane(:,plane2)==1);
tp@0:
tp@0: if closwed < pi
tp@0: okplanes = [okplanelist1;okplanelist2];
tp@0: okplanes = sort(okplanes);
tp@0: else
tp@0: okplanes = zerosvec1;
tp@0: okplanes(okplanelist1) = okplanes(okplanelist1)+1;
tp@0: okplanes(okplanelist2) = okplanes(okplanelist2)+1;
tp@0: okplanes = find(okplanes==2);
tp@0: end
tp@0:
tp@0: okedges = edgesatplane(okplanes,:);
tp@0: [n1,n2] = size(okedges);
tp@0: okedges = reshape(okedges,n1*n2,1);
tp@0:
tp@0: if ~isempty(okedges)
tp@0: okedges = okedges(find(okedges~=0));
tp@0: edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1);
tp@0: end
tp@0:
tp@0: edgesatplane1 = edgesatplane(plane1,:);
tp@0: edgesatplane2 = edgesatplane(plane2,:);
tp@0: okedges = [edgesatplane1 edgesatplane2].';
tp@0: okedges = okedges(find(okedges~=0));
tp@0: edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1);
tp@0:
tp@0: end % else ..... if closwed == 0
tp@0:
tp@0: end % for ii = 1:nedges
tp@0:
tp@0: % Mask out all edge pairs that are aligned
tp@0: edgeseesedge = double(edgeseesedge) - double(edgealignedwithedge);
tp@0:
tp@0: % Make sure the edge-to-same-edge combos are masked out
tp@0: % edgeseesedge = sign(edgeseesedge.*(edgeseesedge.')).*(1-eye(nedges));
tp@0: edgeseesedge = int8( ( (edgeseesedge & (edgeseesedge.'))>0).*(1-eye(nedges)));
tp@0:
tp@0: % Mask out all edges that should be switched off
tp@0:
tp@0: mask = ones(nedges,nedges);
tp@0: mask(offedges,:) = mask(offedges,:)*0;
tp@0: mask(:,offedges.') = mask(:,offedges.')*0;
tp@0:
tp@0: edgeseesedge = edgeseesedge & mask;
tp@0: clear mask
tp@0:
tp@0:
tp@0: %-----------------------------------------------------------
tp@0: % Go through all edgeseesedge combinations. If two edges see
tp@0: % each other without belonging to the same plane, there should
tp@0: % be a gain factor of 2.
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' checking which edge-to-edge paths that run along planes')
tp@0: end
tp@0:
tp@0: edgetoedgestrength = 2*edgeseesedge;
tp@0: for ii = 1:nedges
tp@0: if sum(edgetoedgestrength(:,ii)) ~= 0
tp@0: plane1 = planesatedge(ii,1);
tp@0: plane2 = planesatedge(ii,2);
tp@0: edgesatplane1 = edgesatplane(plane1,:);
tp@0: edgesatplane2 = edgesatplane(plane2,:);
tp@0: sameplaneedges = [edgesatplane1 edgesatplane2].';
tp@0: sameplaneedges = sameplaneedges(find(sameplaneedges~=0));
tp@0: edgetoedgestrength(sameplaneedges,ii) = edgetoedgestrength(sameplaneedges,ii)*0 + 1;
tp@0: end % % if sum(edge
tp@0: end
tp@0:
tp@0: % an edge can not contribute to itself
tp@0:
tp@0: edgetoedgestrength = edgetoedgestrength.*(1-eye(nedges));
tp@0:
tp@0: % The edges belonging to the same planes, but that are inactive (90 degrees etc)
tp@0: % must be switched off here too.
tp@0: % 011021 This is redundant.
tp@0: edgetoedgestrength = edgetoedgestrength.*(edgeseesedge>0);
tp@0: edgetoedgestrength = int8(edgetoedgestrength);
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % CHECK OBSTRUCTION OF EDGE-TO-EDGE PATHS
tp@0: %
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % We construct two long lists: 'from-edges' and 'to-edges' for which
tp@0: % obstruction should be tested. These lists are made up of all the
tp@0: % combinations in edgeseesedge that have a value 1.
tp@0: % NB! We use symmetry - if edge 1 can see edge 2, then edge 2 can also see
tp@0: % edge 1.
tp@0: %
tp@0: % Both edges should be subdivided into nedgesubs segments. So, the long
tp@0: % lists must be expanded so that each from-edge/to-edge combination is
tp@0: % replaced by nedgesubs^2 as many.
tp@0: % For efficiency, we make a shortlist of the actual edge subdivisions and
tp@0: % we then expand long lists with pointers to this shortlist.
tp@0: %
tp@0: % The final matrix, edgeseespartialedge, will have an integer value which
tp@0: % tells whether two edges see each other:
tp@0: %
tp@0: % 0 Edges can not see each other, or one of the edges is inactive
tp@0: % (an inactive edge has an integer wedge index, or has one plane
tp@0: % which is TOTABS)
tp@0: % 2^nedgesubs-1 Two edges see each other completely.
tp@0: % Other integers Two edges see each other partly.
tp@0: % Neg. integer As above, but in addition, the edges do not belong
tp@0: % to the same plane.
tp@0: %
tp@0: % We check in the visedgesfroms and visedgesfromr lists to check which
tp@0: % combinations we don't need to check.
tp@0:
tp@0: if SHOWTEXT >= 3
tp@0: disp([' Checking for obstructing planes between edges and edges'])
tp@0: disp(' ')
tp@0: disp('NB!!! The value of nedgesubs is temporarily set to 1 for the edge-to-edge vis. test')
tp@0: end
tp@0:
tp@0: obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0: maxedgetoedgevisvalue = 2^(2*nedgesubs)-1;
tp@0:
tp@0: if obstructtestneeded
tp@0:
tp@0: nedgesubsorig = nedgesubs;
tp@0: nedgesubs = 1;
tp@0: edgeseesedge = triu(edgeseesedge);
tp@0:
tp@0: iv = full(find(edgeseesedge~=0));
tp@0:
tp@0: if ~isempty(iv)
tp@0: fromedges = floor(iv/nedges)+1;
tp@0: toedges = iv - (fromedges-1)*nedges;
tp@0: ncombs = length(fromedges);
tp@0:
tp@0: ivcancel = find(ismember(fromedges,edgesnottoworryabout));
tp@0: fromedges(ivcancel) = [];
tp@0: toedges(ivcancel) = [];
tp@0:
tp@0: ivcancel = find(ismember(toedges,edgesnottoworryabout));
tp@0: toedges(ivcancel) = [];
tp@0: fromedges(ivcancel) = [];
tp@0: clear ivcancel
tp@0:
tp@0: % Some of these fromedges-toedges combos might involve two
tp@0: % neighboring edges that form an indent. Those should be
tp@0: % removed before we start checking fro obstructions.
tp@0: %
tp@0: % We check the [fromedges toedges] matrix against the
tp@0: % indentingedgepairs matrix.
tp@0:
tp@0: [tf,loc]= ismember(indentingedgepairs,sort([fromedges toedges].').','rows');
tp@0: if ~isempty(loc)
tp@0: ivmatch = find(loc);
tp@0: fromedges(loc(ivmatch)) = [];
tp@0: toedges(loc(ivmatch)) = [];
tp@0: end
tp@0: ncombs = length(fromedges);
tp@0:
tp@0: [uniqueedges,iunique,junique] = unique([fromedges toedges]);
tp@0: clear fromedges toedges
tp@0:
tp@0: % The lists edgesubcoords are the shortlists.
tp@0:
tp@0: [edgesubcoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(uniqueedges,:),edgeendcoords(uniqueedges,:),edgelengthvec(uniqueedges,:),nedgesubs,1);
tp@0:
tp@0: % The two lists below contain pointers to the first edge segment in the
tp@0: % shortlist.
tp@0:
tp@0: fromcoords_refp1 = nedgesubs*(junique(1:ncombs)-1)+1;
tp@0: tocoords_refp1 = nedgesubs*(junique(ncombs+1:2*ncombs)-1)+1;
tp@0: clear junique
tp@0:
tp@0: % Now the original [fromedges toedges] could be recovered by:
tp@0: % uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)])
tp@0: % reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),6,2)
tp@0: % uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]);
tp@0: % npairs = length(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]));
tp@0: % reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),npairs/2,2)
tp@0:
tp@0: addmask = [0:nedgesubs-1];
tp@0: onesvec1 = ones(1,nedgesubs);
tp@0: ntot1 = ncombs*nedgesubs;
tp@0:
tp@0: % First, the from-list gets repetitions of the same values
tp@0: expandfrom_ref = fromcoords_refp1(:,onesvec1);
tp@0: clear fromcoords_refp1;
tp@0: expandfrom_ref = reshape(expandfrom_ref.',ntot1,1);
tp@0:
tp@0: % Second, the expanded from-list gets repetitions that are
tp@0: % incremented in steps of 1 (because the shortlists have the
tp@0: % edge segments placed after each other).
tp@0: expandfrom_ref = expandfrom_ref(:,onesvec1);
tp@0: expandfrom_ref = expandfrom_ref + addmask(ones(ntot1,1),:);
tp@0: expandfrom_ref = reshape(expandfrom_ref.',ncombs*nedgesubs^2,1);
tp@0:
tp@0: % Same thing for the to-list except that the first expansion
tp@0: % gives repetitions that are incremented in steps of 1.
tp@0: expandto_ref = tocoords_refp1(:,onesvec1);
tp@0: clear tocoords_refp1;
tp@0: expandto_ref = expandto_ref + addmask(ones(ncombs,1),:);
tp@0: expandto_ref = reshape(expandto_ref.',ntot1,1);
tp@0:
tp@0: % Second expansion for the to-list: repetitions
tp@0: expandto_ref = expandto_ref(:,onesvec1);
tp@0: expandto_ref = reshape(expandto_ref.',ncombs*nedgesubs^2,1);
tp@0:
tp@0: % Now we have expanded lists of all combinations, and we can
tp@0: % make lists with the coordinates, the edge numbers, and the weights of
tp@0: % the segments.
tp@0:
tp@0: % edgesubcoords(expandfrom_ref,:) will give all the starting points
tp@0: % "fromcoords"
tp@0: % edgesubcoords(expandto_ref,:) will give all the ending points
tp@0: % "tocoords"
tp@0: % uniqueedges(edgenumberlist(expandfrom_ref)) will give all
tp@0: % the starting edges, "fromedges"
tp@0: % uniqueedges(edgenumberlist(expandto_ref)) will give all
tp@0: % the ending edges, "toedges"
tp@0: % We can make a simple test of which planes that could be excluded
tp@0: % from the obstruction test by checking the edgeseesplane for the
tp@0: % relevant edges. It can be multiplied by the canplaneobstruct list
tp@0: % when calling the findobstructedpaths function.
tp@0: isplaneactiveforobstruction = (sum(edgeseesplane(:,uniqueedges).'==1)>0);
tp@0:
tp@0: global STARTPLANES ENDPLANES
tp@0: global BIGFROMCOORDS BIGTOCOORDS BIGSTARTPLANES BIGENDPLANES
tp@0: global REFTOFROMCOSHO REFTOTOCOSHO
tp@0:
tp@0: FROMCOORDSSHORTLIST = edgesubcoords;
tp@0: REFTOFROMCOSHO = uint32(expandfrom_ref);
tp@0: TOCOORDSSHORTLIST = edgesubcoords;
tp@0: REFTOTOCOSHO = uint32(expandto_ref);
tp@0: if nedges < 256
tp@0: bigfromedge = uint8(uniqueedges(edgenumberlist(expandfrom_ref)).');
tp@0: else
tp@0: bigfromedge = uint16(uniqueedges(edgenumberlist(expandfrom_ref)).');
tp@0: end
tp@0: clear expandfrom_ref
tp@0: bigfromedge = bigfromedge(:);
tp@0: if nedges < 256
tp@0: bigtoedge = uint8(uniqueedges(edgenumberlist(expandto_ref)).');
tp@0: else
tp@0: bigtoedge = uint16(uniqueedges(edgenumberlist(expandto_ref)).');
tp@0: end
tp@0: clear expandto_ref
tp@0: bigtoedge = bigtoedge(:);
tp@0: STARTPLANES = [planesatedge(bigfromedge,1) planesatedge(bigfromedge,2)];
tp@0: ENDPLANES = [planesatedge(bigtoedge,1) planesatedge(bigtoedge,2)];
tp@0: shouldplanebechecked = canplaneobstruct.*isplaneactiveforobstruction;
tp@0: nplanesperbatch = ceil(sum(shouldplanebechecked)/ndiff2batches);
tp@0: if nplanesperbatch > 0
tp@0: temp = cumsum(shouldplanebechecked);
tp@0: diff2batchlist = zeros(ndiff2batches,2);
tp@0: diff2batchlist(1,1) = 1;
tp@0: loopisfinished = 0;
tp@0: for ii = 1:ndiff2batches-1
tp@0: ivtemp = find(temp==nplanesperbatch*ii);
tp@0: if ~isempty(ivtemp)
tp@0: diff2batchlist(ii,2) = ivtemp(1);
tp@0: diff2batchlist(ii+1,1) = ivtemp(1)+1;
tp@0: else
tp@0: loopisfinished = 1;
tp@0: end
tp@0: end
tp@0: diff2batchlist(ii,2) = nplanes;
tp@0: ivtemp = find(diff2batchlist(:,1)>0);
tp@0: diff2batchlist = diff2batchlist(ivtemp,:);
tp@0: ndiff2batches = size(diff2batchlist,1);
tp@0:
tp@0: npathstocheck = size(REFTOFROMCOSHO,1);
tp@0: for ii = 1:ndiff2batches
tp@0: if npathstocheck > 0
tp@0: maskedchkpla = shouldplanebechecked;
tp@0: if ii > 1
tp@0: maskedchkpla(1:diff2batchlist(ii,1)-1) = maskedchkpla(1:diff2batchlist(ii,1)-1)*0;
tp@0: end
tp@0: if ii < ndiff2batches
tp@0: maskedchkpla(diff2batchlist(ii,2)+1:end) = maskedchkpla(diff2batchlist(ii,2)+1:end)*0;
tp@0: end
tp@0: nplanestocheck = sum(maskedchkpla);
tp@0: if SHOWTEXT >= 3,
tp@0: if round(ii/1)*1==ii
tp@0: disp([' Batch ',int2str(ii),' of ',int2str(ndiff2batches),': '])
tp@0: disp([' ',int2str(npathstocheck),' paths to check, ',int2str(nplanestocheck),' planes to check obstruction for'])
tp@0: end
tp@0: end
tp@0:
tp@0: % We create the BIGFROMCOORDS matrices and BIGTOCOORDS matrices
tp@0: % outside since they need to be exactly identical from batch to
tp@0: % batch as long as nplanestocheck is the same.
tp@0:
tp@0: npathstocheck = size(REFTOFROMCOSHO,1);
tp@0:
tp@0: ntot = nplanestocheck*npathstocheck;
tp@0:
tp@0: BIGFROMCOORDS = reshape(repmat(FROMCOORDSSHORTLIST(REFTOFROMCOSHO,:).',[nplanestocheck,1]),3,ntot).';
tp@0: BIGTOCOORDS = reshape(repmat(TOCOORDSSHORTLIST(REFTOTOCOSHO,:).',[nplanestocheck,1]),3,ntot).';
tp@0: BIGSTARTPLANES = reshape(repmat(STARTPLANES.',[nplanestocheck,1]),2,ntot).';
tp@0: BIGENDPLANES = reshape(repmat(ENDPLANES.',[nplanestocheck,1]),2,ntot).';
tp@0:
tp@0: [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_edgetoedge(maskedchkpla,planeseesplane,...
tp@0: planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0:
tp@0:
tp@0: if ~isempty(edgehits)
tp@0: % Look for the edge-to-edge combos that share a
tp@0: % plane. When such combos have hit an edge, it
tp@0: % must mean that they are planes with indents.
tp@0: wasedgehitnotindent = prod(double(diff(sort([STARTPLANES(edgehits,:) ENDPLANES(edgehits,:)].')))).';
tp@0: indentedgehits = find(wasedgehitnotindent==0);
tp@0: if ~isempty(indentedgehits)
tp@0: nonobstructedpaths = setdiff(nonobstructedpaths,edgehits(indentedgehits));
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: if ~isempty(nonobstructedpaths)
tp@0: REFTOFROMCOSHO = REFTOFROMCOSHO(nonobstructedpaths);
tp@0: REFTOTOCOSHO = REFTOTOCOSHO(nonobstructedpaths);
tp@0: STARTPLANES = STARTPLANES(nonobstructedpaths,:);
tp@0: ENDPLANES = ENDPLANES(nonobstructedpaths,:);
tp@0: bigfromedge = bigfromedge(nonobstructedpaths);
tp@0: bigtoedge = bigtoedge(nonobstructedpaths);
tp@0: npathstocheck = size(REFTOFROMCOSHO,1);
tp@0: else
tp@0: REFTOFROMCOSHO = [];
tp@0: REFTOTOCOSHO = [];
tp@0: STARTPLANES = [];
tp@0: ENDPLANES = [];
tp@0: bigfromedge = [];
tp@0: bigtoedge = [];
tp@0: npathstocheck = 0;
tp@0: end
tp@0:
tp@0: end % if npathstocheck > 0
tp@0:
tp@0: end % for ii = 1:ndiff2batches
tp@0: else
tp@0: nonobstructedpaths = 1;
tp@0:
tp@0: end %if nplanesperbatch > 0
tp@0:
tp@0: clear REFTOFROMCOSHO REFTOTOCOSHO STARTPLANES ENDPLANES
tp@0:
tp@0:
tp@0: if ~isempty(nonobstructedpaths)
tp@0:
tp@0: % Add the segment pieces together
tp@0: % NB! We don't try to implement the correct way to do it since
tp@0: % we would need nedgesubs^2 bits to represent all edge-segment
tp@0: % to edge-segment visibility possibilities.
tp@0: % Instead we just establish whether the two edges see each
tp@0: % other fully or partly.
tp@0:
tp@0: visiblesegmentscounter = ones(size(bigfromedge));
tp@0: test = [bigfromedge bigtoedge];
tp@0:
tp@0: ncombs = length(bigfromedge);
tp@0: dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0: ivremove = find(dtest==1);
tp@0:
tp@0: while ~isempty(ivremove)
tp@0: visiblesegmentscounter(ivremove+1) = visiblesegmentscounter(ivremove+1) + visiblesegmentscounter(ivremove);
tp@0: visiblesegmentscounter(ivremove) = [];
tp@0: bigfromedge(ivremove) = [];
tp@0: bigtoedge(ivremove,:) = [];
tp@0:
tp@0: test = [bigfromedge bigtoedge];
tp@0: ncombs = length(bigfromedge);
tp@0: dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0: ivremove = find(dtest==1);
tp@0:
tp@0: end
tp@0:
tp@0: indexvec = uint32((double(bigfromedge)-1)*nedges + double(bigtoedge));
tp@0: if maxedgetoedgevisvalue < 128
tp@0: edgeseespartialedge = int8(zeros(nedges,nedges));
tp@0: elseif maxedgetoedgevisvalue < 32768
tp@0: edgeseespartialedge = int16(zeros(nedges,nedges));
tp@0: else
tp@0: edgeseespartialedge = int32(zeros(nedges,nedges));
tp@0: end
tp@0:
tp@0: iv = find(visiblesegmentscounter==nedgesubs^2);
tp@0: edgeseespartialedge(indexvec(iv)) = maxedgetoedgevisvalue*ones(size(iv));
tp@0: iv = find(visiblesegmentscounter>0 & visiblesegmentscounter= 2 & ~isempty(edgeseespartialedge)
tp@0:
tp@0: if SHOWTEXT >= 4
tp@0: disp(' E2E')
tp@0: end
tp@0:
tp@0: zerosvec4 = zeros(nedges,nedges);
tp@0: Bigre1 = (zerosvec4);
tp@0: Bigthetae1 = (zerosvec4);
tp@0: Bigze1 = (zerosvec4);
tp@0: Bigre2 = (zerosvec4);
tp@0: Bigthetae2 = (zerosvec4);
tp@0: Bigze2 = (zerosvec4);
tp@0: clear zerosvec4
tp@0:
tp@0: for edge1 = 1:nedges
tp@0: if SHOWTEXT >= 4
tp@0: if round(edge1/1)*1 == edge1
tp@0: disp([' Edge no. ',int2str(edge1)])
tp@0: end
tp@0: end
tp@0: edge1coords = [edgestartcoords(edge1,:);edgeendcoords(edge1,:)];
tp@0: iv = find(edgeseespartialedge(edge1,:)~=0).';
tp@0:
tp@0: if ~isempty(iv),
tp@0: % First find the subset of edges which belong to the same plane as the
tp@0: % edge itself. These should be treated separately for higher accuracy
tp@0:
tp@0: % iv1 will be the edges on the reference plane. They should have theta = 0.
tp@0:
tp@0:
tp@0: refplane = planesatedge(edge1,1);
tp@0: ncornersperplanevec = double(ncornersperplanevec);
tp@0: iv1 = edgesatplane( refplane,1:ncornersperplanevec( refplane ));
tp@0: iv1 = iv1( find(iv1~= edge1)).';
tp@0:
tp@0: edgestart = edgestartcoordsnudge(iv1,:);
tp@0: edgeend = edgeendcoordsnudge(iv1,:);
tp@0:
tp@0: [rs,slask,zs,rr,slask,zr] = ...
tp@0: EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0: Bigre1(iv1,edge1) = rs;
tp@0: Bigthetae1(iv1,edge1) = 0*iv1;
tp@0: Bigze1(iv1,edge1) = zs;
tp@0: Bigre2(iv1,edge1) = rr;
tp@0: Bigthetae2(iv1,edge1) = 0*iv1;
tp@0: Bigze2(iv1,edge1) = zr;
tp@0:
tp@0: [samevalues,iivec,jjvec] = intersect(iv,iv1);
tp@0: % Bug found 20050421!! PS
tp@0: % % % Old wrong version:
tp@0: % % % % sameedges = [iivec jjvec];
tp@0: % % % if sum(sum(sameedges)) ~= 0
tp@0: % % % iv( sameedges(1,:).' ) = [];
tp@0: % % % end
tp@0: if ~isempty(samevalues)
tp@0: iv(iivec) = [];
tp@0: end
tp@0:
tp@0: % iv2 will be the edges on the non-reference plane.
tp@0: % They should have theta = 2*pi - closthetavec.
tp@0:
tp@0: if ~isempty(iv)
tp@0:
tp@0: if planesatedge(edge1,2) > 0
tp@0: secplane = planesatedge(edge1,2);
tp@0:
tp@0: iv2 = edgesatplane( secplane,1:ncornersperplanevec(secplane));
tp@0: iv2 = iv2( find(iv2~= edge1)).';
tp@0:
tp@0: if ~isempty(iv2)
tp@0: edgestart = edgestartcoordsnudge(iv2,:);
tp@0: edgeend = edgeendcoordsnudge(iv2,:);
tp@0:
tp@0: [rs,slask,zs,rr,slask,zr] = ...
tp@0: EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0: Bigre1(iv2,edge1) = rs;
tp@0: Bigthetae1(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
tp@0: Bigze1(iv2,edge1) = zs;
tp@0: Bigre2(iv2,edge1) = rr;
tp@0: Bigthetae2(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
tp@0: Bigze2(iv2,edge1) = zr;
tp@0:
tp@0: % sameedges = EDB1findsame(iv,iv2);
tp@0: % Bug found 20050421!!! PS
tp@0: % % % % Old wrong version:
tp@0: % % % % [samevalues,iivec,jjvec] = intersect(iv,iv1);
tp@0: % % % % sameedges = [iivec;jjvec];
tp@0: % % % % if sum(sum(sameedges)) ~= 0
tp@0: % % % % iv( sameedges(1,:).' ) = [];
tp@0: % % % % end
tp@0:
tp@0: [samevalues,iivec,jjvec] = intersect(iv,iv2);
tp@0: if ~isempty(samevalues)
tp@0: iv(iivec) = [];
tp@0: end
tp@0: end
tp@0: end
tp@0: end
tp@0: if ~isempty(iv)
tp@0:
tp@0: % Move the edge coordinates to be checked a short distance away
tp@0:
tp@0: edgestart = edgestartcoordsnudge(iv,:);
tp@0: edgeend = edgeendcoordsnudge(iv,:);
tp@0:
tp@0: [rs,thetas,zs,rr,thetar,zr] = ...
tp@0: EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
tp@0: Bigre1(iv,edge1) = rs;
tp@0: Bigthetae1(iv,edge1) = thetas;
tp@0: Bigze1(iv,edge1) = zs;
tp@0: Bigre2(iv,edge1) = rr;
tp@0: Bigthetae2(iv,edge1) = thetar;
tp@0: Bigze2(iv,edge1) = zr;
tp@0:
tp@0: end
tp@0: end
tp@0: end
tp@0:
tp@0: %-----------------------------------------------------------
tp@0: % Go through all edges that are in-plane with each other, that is, two
tp@0: % edges have at least one plane each that are in-plane with each other.
tp@0:
tp@0: % First we identify all planes that have at least one more co-planar
tp@0: % plane
tp@0:
tp@0: planehascoplanar = any(planeseesplane == -1);
tp@0:
tp@0: % Then we go through the list of edges and every edge with a connected
tp@0: % plane that has another co-planar plane must also have potential
tp@0: % in-plane edges.
tp@0:
tp@0: % For each edge, we check if there are other edges (that don't belong to the same plane!!)
tp@0: % with thetaangle = 0 or thetaangle = 2*pi. If there are other such
tp@0: % edges, those edge-to-edge combinations should be shut off.
tp@0: %
tp@0: % After all such edge-to-edge combinations have been shut off, we still
tp@0: % need another pass to cancel edge-to-edge paths that pass entirely
tp@0: % across other edges.
tp@0:
tp@0: ivedges = 1:nedges;
tp@0: ivedges(offedges) = [];
tp@0:
tp@0: for ii = ivedges
tp@0: if any( planehascoplanar(planesatedge(1,:)) )
tp@0:
tp@0: ivcancelcombs = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ( (Bigthetae1(:,ii) == 0) | (Bigthetae1(:,ii) == 2*pi) ) );
tp@0: if ~isempty(ivcancelcombs)
tp@0: edgeseespartialedge(ivcancelcombs,ii) = 0;
tp@0: edgeseespartialedge(ii,ivcancelcombs.') = 0;
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: % The only edge-to-edge combinations that we can easily discard are
tp@0: % those for which edges are aligned with each other (same zstart and
tp@0: % zend): only the closest of those should be kept!
tp@0:
tp@0: for ii = ivedges
tp@0: if any( planehascoplanar(planesatedge(1,:)) )
tp@0:
tp@0:
tp@0: ivotheredges = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ((Bigthetae1(:,ii) == pi)) );
tp@0: if ~isempty(ivotheredges)
tp@0:
tp@0: nudgeval = 1e-10*edgelengthvec(ii)*100;
tp@0: ivselectedges = find( abs(Bigze1(ivotheredges,ii)) 1
tp@0: disp(['Edge no. ',int2str(ii)])
tp@0:
tp@0: lengthok = abs([Bigze1(ivotheredges(ivselectedges),ii) Bigze2(ivotheredges(ivselectedges),ii)] - edgelengthvec(ii)) < nudgeval;
tp@0: ivstillok = find(any(lengthok.'));
tp@0:
tp@0: if ~isempty(ivstillok)
tp@0: ivotheredges = ivotheredges(ivselectedges(ivstillok));
tp@0: else
tp@0: ivotheredges = [];
tp@0: end
tp@0:
tp@0: meanradialdist = mean( [Bigre1(ivotheredges,ii) Bigre2(ivotheredges,ii)].' );
tp@0:
tp@0: ivshutoff = ivotheredges;
tp@0: noshutoff = find(meanradialdist == min(meanradialdist));
tp@0: ivshutoff(noshutoff) = [];
tp@0:
tp@0: if ~isempty(ivshutoff)
tp@0: edgeseespartialedge(ivshutoff,ii) = 0;
tp@0: edgeseespartialedge(ii,ivshutoff) = 0;
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: end
tp@0:
tp@0: %-----------------------------------------------------------------------
tp@0: % Find identical combinations
tp@0:
tp@0: if SHOWTEXT >= 3
tp@0: disp(' Looking for identical combinations')
tp@0: end
tp@0: [reftoshortlistE,re1sho,re2sho,thetae1sho,thetae2sho,...
tp@0: ze1sho,ze2sho,edgelengthsho,examplecombE] = EDB1compress7(Bigre1,Bigre2,...
tp@0: Bigthetae1,Bigthetae2,Bigze1,Bigze2, edgelengthvec(:,ones(1,nedges)).');
tp@0:
tp@0: else
tp@0: reftoshortlistE = [];
tp@0: re1sho = [];
tp@0: re2sho = [];
tp@0: thetae1sho = [];
tp@0: thetae2sho = [];
tp@0: ze1sho = [];
tp@0: ze2sho = [];
tp@0: examplecombE = [];
tp@0: end
tp@0:
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % SAVE THE VARIABLES
tp@0: %
tp@0: %----------------------------------------------------------------------------
tp@0: %
tp@0: % Also the variables from cadgeo??
tp@0: % Yes: planeeqs planenvecs ncornersperplanevec planecorners corners
tp@0: % minvals maxvals
tp@0:
tp@0: edgeseesplane = int8(edgeseesplane);
tp@0: planeseesplane = int8(planeseesplane);
tp@0:
tp@0: ncorners = size(corners,1);
tp@0: if ncorners < 256
tp@0: planecorners = uint8(planecorners);
tp@0: elseif ncorners < 65536
tp@0: planecorners = uint16(planecorners);
tp@0: end
tp@0: planeisthin = uint8(planeisthin);
tp@0:
tp@0: Varlist = [ ' edgecorners planesatedge closwedangvec planehasindents'];
tp@0: Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
tp@0: Varlist = [Varlist,' corners planecorners planeeqs planenvecs ncornersperplanevec'];
tp@0: Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
tp@0: Varlist = [Varlist,' minvals maxvals offedges edgestartcoordsnudge edgeendcoordsnudge'];
tp@0: Varlist = [Varlist,' reftoshortlistE re1sho re2sho thetae1sho thetae2sho ze1sho ze2sho examplecombE'];
tp@0: Varlist = [Varlist,' edgeseespartialedge int_or_ext_model'];
tp@0: Varlist = [Varlist,' edgealignedwithedge edgeperptoplane edgeplaneperptoplane1 edgeplaneperptoplane2'];
tp@0:
tp@0: eval(['save ',outputfile,Varlist])