Mercurial > hg > human-echolocation
view private/EDB1diff2ISES.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 [edgedifflist,startandendpoints,prespeclist,midspeclist,postspeclist,validIScoords,validIRcoords,listguide,... listofallspecs] = EDB1diff2ISES(eddatafile,S,R,ivNdiffmatrix,... lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,... ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlistin,validEDIRcoords,edgeplaneperptoplane1,edgeplaneperptoplane2) % EDB1diff2ISES - Gives list of paths that includes a 2nd-order diff. combination. % EDB1diff2ISES gives the list of possible diffraction paths that includes one % second-order diffraction path, and possibly specular reflections before % and after. % % Input parameters: % eddatafile,S,R,ivNdiffmatrix,... % lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,... % ndecimaldivider,edgeplaneperptoplane1,edgeplaneperptoplane2 % All these are taken from the ISEStreefile or the % setup file. % edgedifflistin List [nd1combs,1] of edges involved in all % first-order diffraction combs that have already % been found % postspeclistin Matrix [nd1combs,specorder-1] of all specular reflections % following the diffraction for all first-order diff combs % that have already been found % bigedgeweightlistin List [nd1combs,1] of visibility for edges involved in all % first-order diffraction combs that have already been found % validEDIRcoords Matrix [nd1combs,3] of image receiver coordinates for all % first-order diff combs that have already been found % GLOBAL parameters: % POTENTIALISES,ISCOORDS,ORIGINSFROM,ISESVISIBILITY,REFLORDER See EDB1findISEStree % SHOWTEXT JJ JJnumbofchars See EDB1mainISES % % Output parameters: % edgedifflist List [ncombs,2] of the edge numbers involved in each % spec-diff-diff-spec combination. % startandendpoints Matrix [ncombs,4] of the relative start and end % points of each edge. The values, [0,1], indicate % which part of the two edges that are visible. % prespeclist Matrix [ncombs,specorder-2] of the specular % reflections that precede every diffraction. % midspeclist Matrix [ncombs,specorder-2] of the specular % reflections inbetween the two diffractions. % postspeclist Matrix [ncombs,specorder-2] of the specular % reflections that follow every diffraction. % validIScoords Matrix [ncombs,3] of the image source for each % multiple-spec that precedes the diffraction. If % there is no spec refl before the diffraction, the % value [0 0 0] is given. % validIRcoords Matrix [ncombs,3] of the image receiver for each % multiple-spec that follows the diffraction. If % there is no spec refl after the diffraction, the % value [0 0 0] is given. % listguide Matrix [nuniquecombs,3] which for each row gives % 1. The number of examples in edgefdifflist etc that % are the same type of spec-diff-diff-spec comb. % 2. The first row number and 3. The last row number. % listofallspecs Matrix [nuniquecombs,3] which for each row gives % 1. The number of pre-specular reflections for the spec-diff-spec-diff-spec comb % in the same row in listguide. % 2. The number of mid-specular reflections for the spec-diff-spec-diff-spec comb % in the same row in listguide. % 3. The number of post-specular reflections for the spec-diff-spec-diff-spec comb % in the same row in listguide. % % Uses functions EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths % % ---------------------------------------------------------------------------------------------- % 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) 20050202 % % [edgedifflist,startandendpoints,prespeclist,midspeclist,postspeclist,validIScoords,validIRcoords,listguide,... % listofallspecs] = EDB1diff2ISES(eddatafile,S,R,... % ivNdiffmatrix,lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,... % vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlist,validEDIRcoords,edgeplaneperptoplane1,edgeplaneperptoplane2) global SHOWTEXT JJ JJnumbofchars global POTENTIALISES ISCOORDS ORIGINSFROM ISESVISIBILITY REFLORDER eval(['load ',eddatafile]) [nedges,slask] = size(planesatedge); [nplanes,slask] = size(planecorners); multfac = 10^(ceil(log10(max([nedges nplanes])))); edgedifflist = []; prespeclist = []; midspeclist = []; postspeclist = []; startandendpoints = []; bigedgeweightlist = []; validIScoords = []; validIRcoords = []; [n1,n2] = size(POTENTIALISES); if n2 < specorder specorder = n2; end maxvisibilityvalue = 2^nedgesubs-1; zerosvec1 = zeros(1,specorder-2); zerosvec2 = zeros(1,3); listguide = zeros(specorder*2-1,3); bitmultvec = 2.^[0:nedgesubs-1]; obstructtestneeded = (sum(canplaneobstruct)~=0); onesvec = ones(1,nedgesubs); onesvec3 = ones(1,3); npostspecs = sum(double(postspeclistin.'>0)).'; planesatedge = double(planesatedge); if specorder >= 3 coplanarsviaflatedge = sparse(zeros(nplanes,nplanes)); listofflatedges = find(closwedangvec==pi); if ~isempty(listofflatedges) ivreftomatrix = planesatedge(listofflatedges,1) + (planesatedge(listofflatedges,2)-1)*nplanes; coplanarsviaflatedge(ivreftomatrix) = ones(size(ivreftomatrix)); coplanarsviaflatedge = sign(coplanarsviaflatedge + coplanarsviaflatedge.'); end end if specorder >= 4 cateyeplanecombs = sparse(zeros(nplanes,nplanes)); listof90edges = find(closwedangvec==3*pi/2); if ~isempty(listof90edges) ivreftomatrix = planesatedge(listof90edges,1) + (planesatedge(listof90edges,2)-1)*nplanes; cateyeplanecombs(ivreftomatrix) = ones(size(ivreftomatrix)); cateyeplanecombs = sign(cateyeplanecombs + cateyeplanecombs.'); end end % ########################################### % # # % # S - edge - edge - R cases # % # # % ########################################### % % Possible edges for S-E-E-R are seen (at least partly) by the source and by the receiver. % % The visibility by the source is taken care of by the findISEStree % so we must check which ones are visible by the receiver. % Also, the active edge segment(s) must be selected but this is done % further down together with the S-spec-spec-edge-R cases ivNdiff = ivNdiffmatrix(1:lengthNdiffmatrix(2),2); ivsinglediff = ivNdiffmatrix(1:lengthNdiffmatrix(1),1); ivSEER = ivNdiff(find(REFLORDER(ivNdiff)==2)); possibleedgepairs = double(POTENTIALISES(ivSEER,1:2))-nplanes; ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0); if ~isempty(ivnotvisiblefromr) possibleedgepairs(ivnotvisiblefromr,:) = []; end edgedifflist = [edgedifflist;possibleedgepairs]; bigedgeweightlist = [bigedgeweightlist;[vispartedgesfromS(edgedifflist(:,1)) vispartedgesfromR(edgedifflist(:,2))]]; [nedgesadded,slask] = size(edgedifflist); zerosvec3 = zeros(nedgesadded,1); ndiffonly = nedgesadded; prespeclist = [prespeclist; zerosvec3(:,ones(1,max(specorder-2,1)))]; midspeclist = [midspeclist; zerosvec3(:,ones(1,max(specorder-2,1)))]; postspeclist = [postspeclist; zerosvec3(:,ones(1,max(specorder-2,1)))]; validIScoords = [validIScoords;S(ones(nedgesadded,1),:)]; validIRcoords = [validIRcoords;R(ones(nedgesadded,1),:)]; if SHOWTEXT >= 3 disp([' ',int2str(nedgesadded),' double diff valid']) end % ########################################### % # # % # S - spec - edge - edge - R cases # % # # % # Prespec cases # % # # % ########################################### % % Possible edges for S-spec-E-E-R are seen (at least partly) by the receiver. % % The visibility doesn't need to be checked since the source-to-edge paths % were checked in the ISEStree, and the visibility from the receiver also % has been checked. % The vector ivmultidiff will always refer to the original data vector % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc % % We should remove the combinations that involve an edge which the % receiver can not see, but since the edge number is in different columns % we do that in the for loop. % The ii-loop will go through: spec-diff, spec-spec-diff % spec-spec-spec-diff etc for ii = 1:specorder-2 if SHOWTEXT >= 3 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before the double edge diff']) end % Select the combinations where the reflection order == ii+2 % which means ii specular reflections before the diffraction iv = find(REFLORDER(ivNdiff) == ii+2 & POTENTIALISES(ivNdiff,ii+1)>nplanes & POTENTIALISES(ivNdiff,ii+2)>nplanes); masterivlist = ivNdiff(iv); possibleedgepairs = double(POTENTIALISES(masterivlist,ii+1:ii+2)) - nplanes; % Keep only combinations for which the receiver can see the edge ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0); if ~isempty(ivnotvisiblefromr) masterivlist(ivnotvisiblefromr) = []; possibleedgepairs(ivnotvisiblefromr,:) = []; end possiblecombs = POTENTIALISES(masterivlist,1:ii); edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgepairs(:,2))]; nposs = length(masterivlist); if SHOWTEXT >= 3 disp([' ',int2str(nposs),' IS+edge pairs valid']) end edgedifflist = [edgedifflist;possibleedgepairs]; prespeclist = [prespeclist;[possiblecombs zeros(nposs,specorder-2-ii)]]; midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)]; postspeclist = [postspeclist;zerosvec1(ones(nposs,1),:)]; bigedgeweightlist = [bigedgeweightlist;edgeweightlist]; % NB! It is correct below that the indices for the ISCOORDS should be % ORIGINSFROM(ORIGINSFROM(masterivlist)), rather than masterivlist. % The combinations in POTENTIALISES(masterivlist,:) all have % spec-spec-...-diff-diff combinations and then % ISCOORDS(masterivlist,:) are zeros since a comb. that % ends with a diff has no image source. % Also, two recursive references are needed since we need to get back % through the two last diffractions to reach the last specular % reflection. validIScoords = [validIScoords;ISCOORDS(ORIGINSFROM(ORIGINSFROM(masterivlist)),:)]; validIRcoords = [validIRcoords;R(ones(nposs,1),:)]; end % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % ############################################# % # # % # S - edge - edge -spec - R cases # % # # % # Postspec cases # % # # % ############################################# % % Possible edges for S-E-E-spec-spec-R are seen (at least partly) by the receiver. % % For some of the combos, we have the IR coords and the edge visibility % from the single-diffraction run. For potential combos that were not % included there, they are either invisible/obstructed or were not tested. % We can figure out which ones were tested because they can be found in the % POTENTIALISES under edge-spec-spec but not in the final valid list. % The vector masterivlist will always refer to the original data vector % i.e. PotentialIR, OriginsfromIR, reflorderIR etc % % First we pick out those indices where there was a single diffraction, but % skip those with only diffraction (because we dealt with them already). % Also, select only those where the diffraction is the first in the sequence % of reflections. for ii = 1:specorder-2 if SHOWTEXT >= 3 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the double edge diff']) end % Select the combinations where the reflection order == ii+2 % (which means ii specular reflections in addition to the diffraction) % and where the first two columns in POTENTIALISES contain edges. iv = find(REFLORDER(ivNdiff) == ii+2 & POTENTIALISES(ivNdiff,1)>nplanes & POTENTIALISES(ivNdiff,2)>nplanes); masterivlist = ivNdiff(iv); possibleedgepairs = double(POTENTIALISES(masterivlist,1:2)) - nplanes; possiblecombs = POTENTIALISES(masterivlist,3:2+ii); possibleweights = ISESVISIBILITY(masterivlist); % Compare with those combinations that were found OK % in the first-order diffraction step (EDB1diffISES), % and that were input matrices to EDB1diff2ISES. % The index vector ivOK refers to these input matrices % (edgedifflistin,posspeclistin,validEDIRcoords,bigedgeweightlistin) % npostspecs is a list that was calculated inside EDB1diff2ISES but % refers to the input matrices. ivOK = find(npostspecs==ii); if ~isempty(ivOK) patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:ii)]; % Find out which ones, of all the possible first-order diffraction combos % in POTENTIALISES, that were indeed tested and found % invisible/obstructed in EDB1diffISES. ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == ii+1)); patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+ii))]; if ~isempty(patternOK) & ~isempty(patternALL) patternNOTOK = setdiff(patternALL,patternOK,'rows'); else if isempty(patternOK) patternNOTOK = patternALL; else % Then patternALL must be empty patternNOTOK = []; end end % Now, the patterns in patternNOTOK can be removed from % masterivlist. patterntocompare = [possibleedgepairs(:,2) possiblecombs(:,1:ii)]; ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows')); masterivlist(ivtocancel) = []; possibleedgepairs(ivtocancel,:) = []; possiblecombs(ivtocancel,:) = []; possibleweights(ivtocancel,:) = []; patterntocompare(ivtocancel,:) = []; [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows'); ivmustbechecked = find(ivcompletelyOK==0); ivcompletelyOK = find(ivcompletelyOK); if ~isempty(ivmustbechecked) masterlisttocheckmore = masterivlist(ivmustbechecked); ntocheckmore = length(masterlisttocheckmore); %---------------------------------------------- % Must carry out a visibility and obstruction check for the special % combinations here. % These combinations have a postspec-combination that hasn't been % encountered in the single diffraction cases, so no visibility % test has been made for these. lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,2))-nplanes; newIRcoords = R; reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii); for jj = 1:ii reflplanes = POTENTIALISES(masterlisttocheckmore,3+ii-jj); reflplanesexpand(:,jj) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1); newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,1,onesvec3); newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).'; eval(['newIRcoords',JJ(jj,1:JJnumbofchars(jj)),' = newIRcoordsexpand;']) end [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs); tocoords = toedgecoords; lastedgenumbers = lastedgenumbers(:,onesvec); lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1); masterlisttocheckmore = masterlisttocheckmore(:,onesvec); masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1); ntocheckmore = length(masterlisttocheckmore); if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' special edge+edge+IR combinations to check']) end for jj = ii:-1:1 if length(masterlisttocheckmore) > 0 eval(['fromcoords = newIRcoords',JJ(jj,1:JJnumbofchars(jj)),';']); if jj < ii eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';']) end [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,jj),4),planenvecs(reflplanesexpand(:,jj),:),minvals(reflplanesexpand(:,jj),:),... maxvals(reflplanesexpand(:,jj),:),planecorners(reflplanesexpand(:,jj),:),corners,ncornersperplanevec(reflplanesexpand(:,jj))); if ~isempty(edgehits) | ~isempty(cornerhits) disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not') disp(' handled correctly yet.') end eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;']) masterlisttocheckmore = masterlisttocheckmore(hitplanes); edgeweightlist = edgeweightlist(hitplanes); lastedgenumbers = lastedgenumbers(hitplanes); reflplanesexpand = reflplanesexpand(hitplanes,:); toedgecoords = toedgecoords(hitplanes,:); for kk = 1:ii eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']); if kk > jj eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']); end end ntocheckmore = length(masterlisttocheckmore); end if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' of the special edge+edge+IR combinations survived the visibility test in refl plane ',int2str(jj)]) end end % Obstruction test of all the involved paths: R -> % reflplane1 -> reflplane2 -> ... -> last edge if obstructtestneeded & ntocheckmore > 0 for jj = 1:ii+1 if ntocheckmore > 0 if jj == 1 fromcoords = R; startplanes = []; else eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';']) startplanes = reflplanesexpand(:,jj-1); end if jj == ii+1, tocoords = toedgecoords; endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)]; else eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';']) endplanes = reflplanesexpand(:,jj); end [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,... planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane); if nobstructions > 0 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths); edgeweightlist = edgeweightlist(nonobstructedpaths); lastedgenumbers = lastedgenumbers(nonobstructedpaths); reflplanesexpand = reflplanesexpand(nonobstructedpaths,:); toedgecoords = toedgecoords(nonobstructedpaths,:); for kk = 1:ii eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']); eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']); end end ntocheckmore = length(masterlisttocheckmore); end end if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' of the special edge+edge+IR combinations survived the obstruction test']) end end % Add the found special combinations to the outdata list edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,1:2))-nplanes]; prespeclist = [prespeclist;zerosvec1(ones(ntocheckmore,1),:)]; midspeclist = [midspeclist;zerosvec1(ones(ntocheckmore,1),:)]; postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-2-ii)]]; bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]]; eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(ii,1:JJnumbofchars(ii)),'];']); validIScoords = [validIScoords;S(ones(ntocheckmore,1),:)]; end masterivlist = masterivlist(ivcompletelyOK); possibleedgepairs = possibleedgepairs(ivcompletelyOK,:); possiblecombs = possiblecombs(ivcompletelyOK,:); possibleweights = possibleweights(ivcompletelyOK,:); nposs = length(ivcompletelyOK); if SHOWTEXT >= 3 disp([' ',int2str(nposs),' Edge+edge+IR segments (non-special combinations) survived the obstruction test']) end % Add the found "standard" combinations to the outdata list edgedifflist = [edgedifflist;possibleedgepairs]; prespeclist = [prespeclist;zerosvec1(ones(nposs,1),:)]; midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)]; postspeclist = [postspeclist;[possiblecombs zeros(nposs,specorder-2-ii)]]; bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]]; validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)]; validIScoords = [validIScoords;S(ones(nposs,1),:)]; end end % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % #################################################### % # # % # S - spec - edge - edge - spec - R cases # % # # % # Pre and postspec cases # % # # % #################################################### % The ii- and jj-loops will go through all combinations of specular % reflection before and after the double diffraction. % % Unlike before, we will look in the list of already found combinations % (edgedifflist etc) for ii = 1:specorder-3 for jj = 1:specorder-ii-2 if SHOWTEXT >= 3 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before and ',JJ(jj,1:JJnumbofchars(jj)),' spec refl after the double edge diff']) end % Select the combinations where the reflection order == ii+jj+2 % (which means ii+jj specular reflections in addition to the % diffraction), and where columns ii+1 and ii+2 of POTENTIALISES % contain edges. iv = find(REFLORDER(ivNdiff) == ii+jj+2 & POTENTIALISES(ivNdiff,ii+1)>nplanes & POTENTIALISES(ivNdiff,ii+2)>nplanes); masterivlist = ivNdiff(iv); possibleedgepairs = double(POTENTIALISES(masterivlist,ii+1:ii+2)) - nplanes; possibleprespecs = POTENTIALISES(masterivlist,1:ii); possiblepostspecs = POTENTIALISES(masterivlist,ii+3:ii+2+jj); possibleweights = ISESVISIBILITY(masterivlist); % Compare with those that have already been found OK ivOK = find(npostspecs==jj); if ~isempty(ivOK) patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:jj)]; else patternOK = []; end % Find out which ones have been checked and found invisible/obstructed ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == jj+1)); patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+jj))]; if ~isempty(patternOK) & ~isempty(patternALL) patternNOTOK = setdiff(patternALL,patternOK,'rows'); else if isempty(patternOK) patternNOTOK = patternALL; else % Then patternALL must be empty patternNOTOK = []; end end patterntocompare = [possibleedgepairs(:,2) possiblepostspecs(:,1:jj)]; ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows')); masterivlist(ivtocancel) = []; possibleedgepairs(ivtocancel,:) = []; possibleprespecs(ivtocancel,:) = []; possiblepostspecs(ivtocancel,:) = []; possibleweights(ivtocancel,:) = []; patterntocompare(ivtocancel,:) = []; [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows'); ivmustbechecked = find(ivcompletelyOK==0); ivcompletelyOK = find(ivcompletelyOK); if ~isempty(ivmustbechecked) masterlisttocheckmore = masterivlist(ivmustbechecked); ntocheckmore = length(masterlisttocheckmore); %---------------------------------------------- % Must carry out a visibility and obstruction check for the special % combinations here. % % NB! toedgecoords are the coordinates of the last edge in the sequence. % This name is because for post-specular reflections, the propagation % is viewed from the receiver towards the last edge! lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,ii+2))-nplanes; newIRcoords = R; reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii); for kk = 1:jj reflplanes = POTENTIALISES(masterlisttocheckmore,3+ii+jj-kk); reflplanesexpand(:,kk) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1); newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,1,onesvec3); newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).'; eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoordsexpand;']) end [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs); tocoords = toedgecoords; lastedgenumbers = lastedgenumbers(:,onesvec); lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1); masterlisttocheckmore = masterlisttocheckmore(:,onesvec); masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1); ntocheckmore = length(masterlisttocheckmore); if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' special IS+edge+edge+IR combinations to check']) end for kk = jj:-1:1 if length(masterlisttocheckmore) > 0 eval(['fromcoords = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),';']); if kk < jj eval(['tocoords = reflpoints',JJ(kk+1,1:JJnumbofchars(kk+1)),';']) end [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,kk),4),planenvecs(reflplanesexpand(:,kk),:),minvals(reflplanesexpand(:,kk),:),... maxvals(reflplanesexpand(:,kk),:),planecorners(reflplanesexpand(:,kk),:),corners,ncornersperplanevec(reflplanesexpand(:,kk))); if ~isempty(edgehits) | ~isempty(cornerhits) disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not') disp(' handled correctly yet.') end eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints;']) masterlisttocheckmore = masterlisttocheckmore(hitplanes); edgeweightlist = edgeweightlist(hitplanes); lastedgenumbers = lastedgenumbers(hitplanes); reflplanesexpand = reflplanesexpand(hitplanes,:); toedgecoords = toedgecoords(hitplanes,:); for ll = 1:jj eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']); if ll > kk eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']); end end ntocheckmore = length(masterlisttocheckmore); end if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' of the special IS+edge+edge+IR combinations survived the visibility test in refl plane ',int2str(jj)]) end end % Obstruction test of all the involved paths: R -> % reflplane1 -> reflplane2 -> ... -> last edge if obstructtestneeded & ntocheckmore > 0 for kk = 1:jj+1 if ntocheckmore > 0 if kk == 1 fromcoords = R; startplanes = []; else eval(['fromcoords = reflpoints',JJ(kk-1,1:JJnumbofchars(kk-1)),';']) startplanes = reflplanesexpand(:,kk-1); end if kk == jj+1, tocoords = toedgecoords; endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)]; else eval(['tocoords = reflpoints',JJ(kk,1:JJnumbofchars(kk)),';']) endplanes = reflplanesexpand(:,kk); end [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,... planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane); if nobstructions > 0 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths); edgeweightlist = edgeweightlist(nonobstructedpaths); lastedgenumbers = lastedgenumbers(nonobstructedpaths); reflplanesexpand = reflplanesexpand(nonobstructedpaths,:); toedgecoords = toedgecoords(nonobstructedpaths,:); for ll = 1:jj eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']); eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']); end end ntocheckmore = length(masterlisttocheckmore); end end if SHOWTEXT >= 3 disp([' ',int2str(ntocheckmore),' of the special IS+edge+edge+IR combinations survived the obstruction test']) end end % Add the found special combinations to the outdata list edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,ii+1:ii+2))-nplanes]; prespeclist = [prespeclist;[double(POTENTIALISES(masterlisttocheckmore,1:ii)) zeros(ntocheckmore,specorder-2-ii)]]; midspeclist = [midspeclist;zerosvec1(ones(ntocheckmore,1),:)]; postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-2-ii)]]; bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]]; eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(jj,1:JJnumbofchars(jj)),'];']); % NB!! In the same way as earlier, we must a recursive reference % method to find the image source of the last specular reflection. ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterlisttocheckmore))); for kk = 2:jj ivref = ORIGINSFROM(ivref); end validIScoords = [validIScoords;ISCOORDS(ivref,:)]; end masterivlist = masterivlist(ivcompletelyOK); possibleedgepairs = possibleedgepairs(ivcompletelyOK,:); possibleprespecs = possibleprespecs(ivcompletelyOK,:); possiblepostspecs = possiblepostspecs(ivcompletelyOK,:); possibleweights = possibleweights(ivcompletelyOK,:); nposs = length(ivcompletelyOK); if SHOWTEXT >= 3 disp([' ',int2str(nposs),' IS+Edge+edge+IR segments (non-special) survived the obstruction test']) end % Add the found "standard" combinations to the outdata list edgedifflist = [edgedifflist;possibleedgepairs]; prespeclist = [prespeclist; [possibleprespecs zeros(nposs,specorder-2-ii)]]; midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)]; postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-2-jj)]]; bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]]; validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)]; % NB!! In the same way as earlier, we must a recursive reference % method to find the image source of the last specular reflection. ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist))); for kk = 2:jj ivref = ORIGINSFROM(ivref); end validIScoords = [validIScoords;ISCOORDS(ivref,:)]; end end % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % ####################################################################### % # % # Midspec cases, with prespecs and postspecs % # % # S - spec - edge - spec - edge - R % # % ####################################################################### for iipre = 0:specorder-3 for jjmid = 1:specorder-2-iipre for kkpost = 0:specorder-iipre-jjmid-2 if SHOWTEXT >= 3 if iipre > 0 JJpre = JJ(iipre,1:JJnumbofchars(iipre)); else JJpre = '0'; end if kkpost > 0 JJpost = JJ(kkpost,1:JJnumbofchars(kkpost)); else JJpost = '0'; end disp([' Checking for ',JJpre,' spec refl before, ',JJ(jjmid,1:JJnumbofchars(jjmid)),' spec refl between the double edge diff, and']) disp([' ',JJpost,' spec refl after the double edge diff.']) end iv = find(REFLORDER(ivNdiff) == iipre+jjmid+kkpost+2 & POTENTIALISES(ivNdiff,iipre+1)>nplanes & POTENTIALISES(ivNdiff,iipre+jjmid+2)>nplanes); masterivlist = ivNdiff(iv); possibleedgepairs = double(POTENTIALISES(masterivlist,[iipre+1 iipre+jjmid+2])) - nplanes; nfast = 0; if kkpost == 0 possiblepostspecs = []; % Keep only combinations for which the receiver can see the edge ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0); if ~isempty(ivnotvisiblefromr) masterivlist(ivnotvisiblefromr) = []; possibleedgepairs(ivnotvisiblefromr,:) = []; end else possiblepostspecs = POTENTIALISES(masterivlist,iipre+jjmid+3:iipre+jjmid+2+kkpost); % Compare with those that have already been found OK ivOK = find(npostspecs==kkpost); if ~isempty(ivOK) patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:kkpost)]; else patternOK = []; end % Find out which ones have been checked and found invisible/obstructed ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == kkpost+1)); patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+kkpost))]; if ~isempty(patternOK) & ~isempty(patternALL) patternNOTOK = setdiff(patternALL,patternOK,'rows'); else if isempty(patternOK) patternNOTOK = patternALL; else % Then patternALL must be empty patternNOTOK = []; end end patterntocompare = [possibleedgepairs(:,2) possiblepostspecs(:,1:kkpost)]; if ~isempty(patternNOTOK) ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows')); masterivlist(ivtocancel) = []; possibleedgepairs(ivtocancel,:) = []; possiblepostspecs(ivtocancel,:) = []; patterntocompare(ivtocancel,:) = []; end [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows'); ivmustbechecked = find(ivcompletelyOK==0); ivcompletelyOK = find(ivcompletelyOK); if ~isempty(ivmustbechecked) & SHOWTEXT > 0 disp('WARNING!! For midspec and postspec case, all checks have not been implemented yet!!') end masterivlist = masterivlist(ivcompletelyOK); possibleedgepairs = possibleedgepairs(ivcompletelyOK,:); possiblepostspecs = possiblepostspecs(ivcompletelyOK,:); nposs = length(ivcompletelyOK); end if iipre > 0 possibleprespecs = POTENTIALISES(masterivlist,1:iipre); else possibleprespecs = []; end % NB!! possiblemidspecs is numbered in reverse order % since we view the propagation by starting to mirror the last edge % and move towards the first edge. possiblemidspecs = POTENTIALISES(masterivlist,iipre+1+jjmid:-1:iipre+2); if kkpost > 0 edgeweightlist = [ISESVISIBILITY(masterivlist) bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]; else edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgepairs(:,2))]; end % Expand the various lists and matrices to represent the % sub-divided edges. nposs = length(masterivlist); if nposs > 0 if iipre == 1 possibleprespecs = reshape(possibleprespecs(:,onesvec).',nposs*nedgesubs,1); elseif iipre > 1 possibleprespecs = reshape(repmat(possibleprespecs.',nedgesubs,1),iipre,nposs*nedgesubs).'; end if jjmid == 1 possiblemidspecs = reshape(possiblemidspecs(:,onesvec).',nposs*nedgesubs,1); elseif jjmid > 1 possiblemidspecs = reshape(repmat(possiblemidspecs.',nedgesubs,1),jjmid,nposs*nedgesubs).'; end if kkpost == 1 possiblepostspecs = reshape(possiblepostspecs(:,onesvec).',nposs*nedgesubs,1); elseif kkpost > 1 possiblepostspecs = reshape(repmat(possiblepostspecs.',nedgesubs,1),kkpost,nposs*nedgesubs).'; end expandedmasterivlist = reshape(masterivlist(:,onesvec).',nposs*nedgesubs,1); if kkpost > 0 expandedivcompletelyOK = reshape(ivcompletelyOK(:,onesvec).',nposs*nedgesubs,1); end edgeweightlist = reshape(repmat(edgeweightlist.',[nedgesubs 1]),2,nposs*nedgesubs).'; for ll = 1:nedgesubs edgeweightlist(ll:nedgesubs:end,1) = double(bitget(edgeweightlist(ll:nedgesubs:end,1),ll))*bitmultvec(ll); edgeweightlist(ll:nedgesubs:end,2) = double(bitget(edgeweightlist(ll:nedgesubs:end,2),ll))*bitmultvec(ll); end %---------------------------------------------- % Must carry out a visibility and obstruction check for the % edge-spec-edge paths % % NB! toedgecoords are the coordinates of the first edge in the sequence % and fromedgecoords refers to the last edge, after the mid-specular % reflections. % This name is because for mid-specular reflections, the propagation % is viewed from the last edge towards the first edge! [toedgecoords,firstedgeweightlist,slask] = EDB1getedgepoints(edgestartcoords(possibleedgepairs(:,2),:),edgeendcoords(possibleedgepairs(:,2),:),edgelengthvec(possibleedgepairs(:,1),:),nedgesubs); tocoords = toedgecoords; [fromedgecoords,lastedgeweightlist,slask] = EDB1getedgepoints(edgestartcoords(possibleedgepairs(:,1),:),edgeendcoords(possibleedgepairs(:,1),:),edgelengthvec(possibleedgepairs(:,2),:),nedgesubs); fromcoords = fromedgecoords; possibleedgepairs = reshape(repmat(possibleedgepairs.',nedgesubs,1),2,nposs*nedgesubs).'; edgeimagecoords = fromedgecoords; for ll = 1:jjmid edgeimagecoords = EDB1findis(edgeimagecoords,possiblemidspecs(:,ll),planeeqs,size(fromedgecoords,1),onesvec3); eval(['bigedgeimagecoords',JJ(ll,1:JJnumbofchars(ll)),' = edgeimagecoords;']) end % Some cases do not need to be checked, when jjmid = 2: the % cateye cases. For these, we will have doubles (both, e.g. % 3-7 and 7-3) and one should be tossed then. The non-doubles % can be kept in a "fast lane" so that visibility isn't % checked, but obstruction is. if jjmid == 2 specpattern = double(possiblemidspecs); ivreftomatrix = specpattern(:,1) + ( specpattern(:,2)-1)*nplanes; ivcateyes = find( cateyeplanecombs(ivreftomatrix) ); if ~isempty(ivcateyes), specpattern = specpattern(ivcateyes,:); fliporder = specpattern(:,2)<specpattern(:,1); ivfliporder = find(fliporder); specpattern(ivfliporder,:) = specpattern(ivfliporder,[2 1]); [uniquepatterns,ivec,jvec] = unique(specpattern,'rows'); countcases = histc(jvec,[1:max(jvec)]); ivtossone = find(fliporder & countcases(jvec)==nedgesubs*2); ivfastlane = ivcateyes; ivfastlane(ivtossone) = []; if ~isempty(ivfastlane) expandedmasterivlistfast = expandedmasterivlist(ivfastlane); possibleedgepairsfast = possibleedgepairs(ivfastlane,:); fromedgecoordsfast = fromedgecoords(ivfastlane,:); toedgecoordsfast = toedgecoords(ivfastlane,:); if ~isempty(possibleprespecs) possibleprespecsfast = possibleprespecs(ivfastlane,:); end possiblemidspecsfast = possiblemidspecs(ivfastlane,:); if ~isempty(possiblepostspecs) possiblepostspecsfast = possiblepostspecs(ivfastlane,:); end edgeweightlistfast = edgeweightlist(ivfastlane,:); for mm = 1:jjmid eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'fast = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(ivfastlane,:);']); end reflpoints2fast = 0.5*(bigedgeimagecoords2fast-fromedgecoordsfast) + fromedgecoordsfast; reflpoints1fast = reflpoints2fast; nfast = length(edgeweightlistfast); end if ~isempty(ivtossone) ivtossone = ivcateyes; expandedmasterivlist((ivtossone)) = []; edgeweightlist((ivtossone),:) = []; possibleedgepairs((ivtossone),:) = []; if ~isempty(possibleprespecs) possibleprespecs((ivtossone),:) = []; end possiblemidspecs((ivtossone),:) = []; if ~isempty(possiblepostspecs) possiblepostspecs((ivtossone),:) = []; end fromcoords((ivtossone),:) = []; tocoords((ivtossone),:) = []; for mm = 1:jjmid eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'((ivtossone),:) = [];']); end nposs = length(expandedmasterivlist); end end end nposs = length(expandedmasterivlist); end if SHOWTEXT >= 3 if jjmid ~= 2 disp([' ',int2str(nposs),' edge+spec+edge combos found initially:']) else disp([' ',int2str(nposs),' edge+spec+edge combos found initially, + ',int2str(nfast),' cateye combos']) end end %-------------------------------------------------------------- % Check the visibility through all the reflection planes % for ll = 1:jjmid eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = [];']) end for ll = jjmid:-1:1 if nposs > 0 eval(['fromcoords = bigedgeimagecoords',JJ(ll,1:JJnumbofchars(ll)),';']) if ll < jjmid eval(['tocoords = reflpoints',JJ(ll+1,1:JJnumbofchars(ll+1)),';']) end [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(possiblemidspecs(:,ll),4),planenvecs(possiblemidspecs(:,ll),:),minvals(possiblemidspecs(:,ll),:),... maxvals(possiblemidspecs(:,ll),:),planecorners(possiblemidspecs(:,ll),:),corners,ncornersperplanevec(possiblemidspecs(:,ll))); % Make a special treatment for the cases with the % specular reflection point right on an edge since some % of these are special cases: % "edgeplaneperptoplane1" indicates that edge-plane-edge % travels along the edge's plane and is reflected at a % 90 degree corner (which is an inactive edge). % These are treated as good hits. % "edgeplaneperptoplane2" indicates that edge-plane-edge % has a specular reflection right at a flat edge % between two coplanar planes. % These come in pairs; one half-hit in the two % coplanar planes. Only one in each pair should be % kept. if jjmid == 1 & ~isempty(edgehits) edge1 = possibleedgepairs(edgehits,1); edge2 = possibleedgepairs(edgehits,2); midspec = possiblemidspecs(edgehits,1); ivreftomatrix1 = double(midspec) + double(edge1-1)*nplanes; ivreftomatrix2 = double(midspec) + double(edge2-1)*nplanes; specialhit = edgeplaneperptoplane1(ivreftomatrix1).*edgeplaneperptoplane1(ivreftomatrix1); ivspecial = find(specialhit); if ~isempty(ivspecial) hitplanes = [hitplanes;edgehits(ivspecial)]; reflpoints = [reflpoints;edgehitpoints(ivspecial,:)]; edgehits(ivspecial) = []; edgehitpoints(ivspecial,:) = []; ivreftomatrix1(ivspecial) = []; ivreftomatrix2(ivspecial) = []; end specialhit = edgeplaneperptoplane2(ivreftomatrix1).*edgeplaneperptoplane2(ivreftomatrix1); ivspecial = find(specialhit); if ~isempty(ivspecial) patternlist = double([possibleedgepairs(edgehits(ivspecial),1) possiblemidspecs(edgehits(ivspecial),1) possibleedgepairs(edgehits(ivspecial),1)]); [uniquepatterns,ivec,jvec] = unique(patternlist,'rows'); keeppattern = zeros(size(ivec)); for mm = 1:length(ivec) ivreftomatrix = uniquepatterns(mm,2) + (uniquepatterns(mm+1:end,2)-1)*nplanes; coplanarindicator = coplanarsviaflatedge(ivreftomatrix); ivcoplanars = find(uniquepatterns(mm+1:end,1)==uniquepatterns(mm,1) & uniquepatterns(mm+1:end,3)==uniquepatterns(mm,3) & coplanarindicator); if ~isempty(ivcoplanars) keeppattern(mm) = 1; end end ivkeepers = find(keeppattern(jvec)); hitplanes = [hitplanes;edgehits(ivspecial(ivkeepers))]; reflpoints = [reflpoints;edgehitpoints(ivspecial(ivkeepers),:)]; end end eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints;']) expandedmasterivlist = expandedmasterivlist(hitplanes); edgeweightlist = edgeweightlist(hitplanes,:); possibleedgepairs = possibleedgepairs(hitplanes,:); if ~isempty(possibleprespecs) possibleprespecs = possibleprespecs(hitplanes,:); end possiblemidspecs = possiblemidspecs(hitplanes,:); if ~isempty(possiblepostspecs) possiblepostspecs = possiblepostspecs(hitplanes,:); end fromedgecoords = fromedgecoords(hitplanes,:); toedgecoords = toedgecoords(hitplanes,:); if kkpost > 0 expandedivcompletelyOK = expandedivcompletelyOK(hitplanes); end for mm = 1:jjmid eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(hitplanes,:);']); if mm > ll eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)),' = reflpoints',JJ(mm,1:JJnumbofchars(mm)),'(hitplanes,:);']); end end nposs = length(expandedmasterivlist); if SHOWTEXT >= 3 if jjmid ~= 2 disp([' ',int2str(nposs),' edge+spec+edge combos survived the visibility test in reflection plane ',int2str(ll)]) else disp([' ',int2str(nposs),' edge+spec+edge combos survived the visibility test in reflection plane ',int2str(ll),' + ',int2str(nfast),' cateye combos']) end end end end %-------------------------------------------------------------- % Check for obstructions for all the paths, starting from edge number 2 % towards edge number 1. % % Reinsert the "fast lane" cases, i.e., the cateye reflections. if jjmid == 2 if nfast > 0, expandedmasterivlist = [expandedmasterivlist;expandedmasterivlistfast]; edgeweightlist = [edgeweightlist;edgeweightlistfast]; possibleedgepairs = [possibleedgepairs;possibleedgepairsfast]; if ~isempty(possibleprespecs) possibleprespecs = [possibleprespecs;possibleprespecsfast]; end possiblemidspecs = [possiblemidspecs;possiblemidspecsfast]; if ~isempty(possiblepostspecs) possiblepostspecs = [possiblepostspecs;possiblepostspecsfast]; end fromedgecoords = [fromedgecoords;fromedgecoordsfast]; toedgecoords = [toedgecoords;toedgecoordsfast]; for mm = 1:jjmid eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = [bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),';bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'fast];']); eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)), ' = [reflpoints',JJ(mm,1:JJnumbofchars(mm)), ';reflpoints',JJ(mm,1:JJnumbofchars(mm)),'fast];']); end nposs = length(expandedmasterivlist); end end if nposs > 0 & obstructtestneeded for ll = 1:jjmid+1 if nposs > 0 if ll == 1 fromcoords = fromedgecoords; startplanes = [planesatedge(possibleedgepairs(:,2),1) planesatedge(possibleedgepairs(:,2),2)]; else eval(['fromcoords = reflpoints',JJ(ll-1,1:JJnumbofchars(ll-1)),';']) startplanes = possiblemidspecs(:,ll-1); end if ll == jjmid+1, tocoords = toedgecoords; endplanes = [planesatedge(possibleedgepairs(:,1),1) planesatedge(possibleedgepairs(:,1),2)]; else eval(['tocoords = reflpoints',JJ(ll,1:JJnumbofchars(ll)),';']) endplanes = possiblemidspecs(:,ll); end [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,... planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane); if ~isempty(edgehits) nonobstructedpaths = setdiff(nonobstructedpaths,edgehits); nobstructions = nposs - length(nonobstructedpaths); end if nobstructions > 0 expandedmasterivlist = expandedmasterivlist(nonobstructedpaths); edgeweightlist = edgeweightlist(nonobstructedpaths,:); possibleedgepairs = possibleedgepairs(nonobstructedpaths,:); if ~isempty(possibleprespecs) possibleprespecs = possibleprespecs(nonobstructedpaths,:); end possiblemidspecs = possiblemidspecs(nonobstructedpaths,:); if ~isempty(possiblepostspecs) possiblepostspecs = possiblepostspecs(nonobstructedpaths,:); end fromedgecoords = fromedgecoords(nonobstructedpaths,:); toedgecoords = toedgecoords(nonobstructedpaths,:); if kkpost > 0 expandedivcompletelyOK = expandedivcompletelyOK(nonobstructedpaths,:); end for mm = 1:jjmid eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(nonobstructedpaths,:);']); if mm > ll eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)),' = reflpoints',JJ(mm,1:JJnumbofchars(mm)),'(nonobstructedpaths,:);']); end end end nposs = length(expandedmasterivlist); end end if SHOWTEXT >= 3 if jjmid == 2 disp([' ',int2str(nposs),' edge+spec+edge combos survived the obstruction test (including cateye cases)']) else disp([' ',int2str(nposs),' edge+spec+edge combos survived the obstruction test']) end end end if nposs > 0 edgedifflist = [edgedifflist;possibleedgepairs]; if iipre == 0 possibleprespecs = zeros(nposs,1); end if specorder <= 4 if specorder == 3 prespeclist = [prespeclist;[possibleprespecs]]; else prespeclist = [prespeclist;[possibleprespecs zeros(nposs,1)]]; end else [n1,n2] = size(prespeclist); [n3,n4] = size(possibleprespecs); if n1 > 0 % Error found 20050202 PS % Old wrong version prespeclist = [prespeclist;[possibleprespecs zeros(nposs,n4-n2)]]; prespeclist = [prespeclist;[possibleprespecs zeros(nposs,n2-n4)]]; else % Error found 20050202 PS % Old wrong version prespeclist = [possibleprespecs zeros(nposs,n4-n2)]; prespeclist = [possibleprespecs zeros(nposs,n2-n4)]; end end if jjmid == specorder-2 midspeclist = [midspeclist;possiblemidspecs]; else midspeclist = [midspeclist;[possiblemidspecs zeros(nposs,specorder-2-jjmid) ]]; end if kkpost == 0 possiblepostspecs = zeros(nposs,1); end if specorder <= 4 if specorder == 3 postspeclist = [postspeclist;[possiblepostspecs]]; else postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,1)]]; end else if kkpost == 0 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-3)]]; else postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-2-kkpost)]]; end end bigedgeweightlist = [bigedgeweightlist;edgeweightlist]; % NB! It is correct below that the indices for the ISCOORDS should be % ORIGINSFROM(ORIGINSFROM(masterivlist)), rather than masterivlist. % The combinations in POTENTIALISES(masterivlist,:) all have % spec-spec-...-diff-diff combinations and then % ISCOORDS(masterivlist,:) are zeros since a comb. that % ends with a diff has no image source. % Also, two recursive references are needed since we need to get back % through the two last diffractions to reach the last specular % reflection. ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(expandedmasterivlist))); for kk = 2:jjmid+kkpost ivref = ORIGINSFROM(ivref); end validIScoords = [validIScoords;ISCOORDS(ivref,:)]; if kkpost == 0 validIRcoords = [validIRcoords;R(ones(nposs,1),:)]; else if jjmid ~= 2 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(expandedivcompletelyOK)),:)]; else error(['ERROR: Calculation of IR coords for midspec = 2 and postspec not implemented yet!']); end end end end end end % ####################################################################### % # % # Pack the edge segments together because every little edge segment % # is present as a separate edge % # This can be done for all combinations at once. % # % ####################################################################### test = [prespeclist edgedifflist(:,1) midspeclist edgedifflist(:,2) postspeclist]; if ~isempty(test) [ncombs,slask] = size(edgedifflist); dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']); ivremove = find(dtest==1); while ~isempty(ivremove) bigedgeweightlist(ivremove+1,:) = double(bigedgeweightlist(ivremove+1,:)) + double(bigedgeweightlist(ivremove,:)); bigedgeweightlist(ivremove,:) = []; edgedifflist(ivremove,:) = []; prespeclist(ivremove,:) = []; midspeclist(ivremove,:) = []; postspeclist(ivremove,:) = []; validIScoords(ivremove,:) = []; validIRcoords(ivremove,:) = []; test = [prespeclist edgedifflist(:,1) midspeclist edgedifflist(:,2) postspeclist]; [ncombs,slask] = size(edgedifflist); dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']); ivremove = find(dtest==1); end end % ####################################################################### % # % # The weights of the visible edge segments should be % # translated into start and end points, together with the visibility % # weights from the receiver. % # This can be done for all combinations at once! % # % ####################################################################### % As a start we set the start and end values to 0 and 1, i.e. assuming full % visibility. startandendpoints = [startandendpoints;... [zeros(length(edgedifflist),1),ones(length(edgedifflist),1),... zeros(length(edgedifflist),1),ones(length(edgedifflist),1)]]; % Treat edge 1 ivtoaccountfor = [1:size(edgedifflist,1)].'; ivwholeedge1 = find( bigedgeweightlist(:,1) == maxvisibilityvalue); if ~isempty(ivwholeedge1) ivtoaccountfor(ivwholeedge1) = []; end if ~isempty(ivtoaccountfor) ncombs = length(ivtoaccountfor); bitpattern = zeros(ncombs,nedgesubs); for ii=1:nedgesubs bitpattern(:,ii) = bitget(bigedgeweightlist(ivtoaccountfor,1),ii); end dbit1 = diff([zeros(ncombs,1) bitpattern].').'; dbit2 = [dbit1 -bitpattern(:,nedgesubs)]; nsegments = ceil((sum(abs(dbit1.')).')/2); ivonesegments = find(nsegments==1); if ~isempty(ivonesegments) nonesegments = length(ivonesegments); multvec = 2.^[0:nedgesubs]; segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1; segendpos = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1; ivmodify = find(segstartpos==1); segstartpos(ivmodify) = ones(size(ivmodify))*1.5; ivmodify = find(segendpos>nedgesubs); segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5); startandendpoints(ivtoaccountfor(ivonesegments),1) = (segstartpos-1.5)/(nedgesubs-1); startandendpoints(ivtoaccountfor(ivonesegments),2) = (segendpos-1.5)/(nedgesubs-1); end % If we have some two-or-more-subsegments cases, they will be % discovered by the if-condition below if length(ivonesegments) < ncombs for nsegmentstocheck = 2:ceil(nedgesubs/2) disp(['Checking for ',int2str(nsegmentstocheck),' sub-segments']) ivNsegments = find(nsegments==nsegmentstocheck); if ~isempty(ivNsegments) [n1,n2] = size(startandendpoints); if n2 < 4*nsegmentstocheck startandendpoints = [startandendpoints zeros(n1,4*nsegmentstocheck-n2)]; end for jj = 1:length(ivNsegments) ivstartbits = find(dbit2(ivNsegments(jj),:) == 1); ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits; ivendbits = find(dbit2(ivNsegments(jj),:) == -1); ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits; for kk = 1:nsegmentstocheck, startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+1) = (ivstartbits(kk)-1.5)/(nedgesubs-1); startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+2) = (ivendbits(kk)-1.5)/(nedgesubs-1); end end end end end end % Treat edge 2 ivtoaccountfor = [1:size(edgedifflist,1)].'; ivwholeedge2 = find( bigedgeweightlist(:,2) == maxvisibilityvalue); if ~isempty(ivwholeedge2) ivtoaccountfor(ivwholeedge2) = []; end if ~isempty(ivtoaccountfor) ncombs = length(ivtoaccountfor); bitpattern = zeros(ncombs,nedgesubs); for ii=1:nedgesubs bitpattern(:,ii) = bitget(bigedgeweightlist(ivtoaccountfor,2),ii); end dbit1 = diff([zeros(ncombs,1) bitpattern].').'; dbit2 = [dbit1 -bitpattern(:,nedgesubs)]; nsegments = ceil((sum(abs(dbit1.')).')/2); ivonesegments = find(nsegments==1); if ~isempty(ivonesegments) nonesegments = length(ivonesegments); multvec = 2.^[0:nedgesubs]; segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1; segendpos = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1; ivmodify = find(segstartpos==1); segstartpos(ivmodify) = ones(size(ivmodify))*1.5; ivmodify = find(segendpos>nedgesubs); segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5); startandendpoints(ivtoaccountfor(ivonesegments),3) = (segstartpos-1.5)/(nedgesubs-1); startandendpoints(ivtoaccountfor(ivonesegments),4) = (segendpos-1.5)/(nedgesubs-1); end % If we have some two-or-more-subsegments cases, they will be % discovered by the if-condition below if length(ivonesegments) < ncombs for nsegmentstocheck = 2:ceil(nedgesubs/2) ivNsegments = find(nsegments==nsegmentstocheck); if ~isempty(ivNsegments) [n1,n2] = size(startandendpoints); if n2 < 4*nsegmentstocheck startandendpoints = [startandendpoints zeros(n1,4*nsegmentstocheck-n2)]; end for jj = 1:length(ivNsegments) ivstartbits = find(dbit2(ivNsegments(jj),:) == 1); ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits; ivendbits = find(dbit2(ivNsegments(jj),:) == -1); ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits; for kk = 1:nsegmentstocheck, startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+3) = (ivstartbits(kk)-1.5)/(nedgesubs-1); startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+4) = (ivendbits(kk)-1.5)/(nedgesubs-1); end end end end end end % ####################################################################### % # % # Construct a list guide, which will tell which rows have only % # dd, which rows have sdd etc % # Syntax: dd,sdd,ssdd,sssdd,...,dds,ddss,ddsss,... % # (should also continue with sdds, sddss,... % # % ####################################################################### [n1,n2] = size(prespeclist); if n2 > 1 nprespecs = sum(prespeclist.' > 0).'; else nprespecs = (prespeclist>0); end [n1,n2] = size(midspeclist); if n2 > 1 nmidspecs = sum(midspeclist.' > 0).'; else nmidspecs = (midspeclist>0); end [n1,n2] = size(postspeclist); if n2 > 1 npostspecs = sum(postspeclist.' > 0).'; else npostspecs = (postspeclist>0); end [B,ivec,jvec] = unique([nprespecs nmidspecs npostspecs],'rows'); nuniquecombs = length(ivec); ntotcombs = length(jvec); listguide = zeros(nuniquecombs,3); listofallspecs = zeros(nuniquecombs,3); sortvec = zeros(ntotcombs,1); for ii = 1:length(ivec) ivfindcombs = find(jvec==ii); listguide(ii,1) = length(ivfindcombs); if ii > 1 listguide(ii,2) = listguide(ii-1,3)+1; else listguide(ii,2) = 1; end listguide(ii,3) = listguide(ii,2)+listguide(ii,1)-1; listofallspecs(ii,:) = [B(ii,:)]; sortvec(listguide(ii,2):listguide(ii,3)) = ivfindcombs; end prespeclist = prespeclist(sortvec,:); midspeclist = midspeclist(sortvec,:); postspeclist = postspeclist(sortvec,:); validIScoords = validIScoords(sortvec,:); validIRcoords = validIRcoords(sortvec,:); edgedifflist = edgedifflist(sortvec,:); startandendpoints = startandendpoints(sortvec,:);