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