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