Mercurial > hg > human-echolocation
diff private/EDB1diff2ISES.m @ 0:90220f7249fc
Initial commit
| author | Timos Papadopoulos <tp@isvr.soton.ac.uk> |
|---|---|
| date | Sat, 16 Nov 2013 18:25:24 +0000 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/private/EDB1diff2ISES.m Sat Nov 16 18:25:24 2013 +0000 @@ -0,0 +1,1514 @@ +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,:);