Mercurial > hg > human-echolocation

annotate private/EDB1diff2ISES.m @ 18:2d5f50205527 jabuilder_int tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Escape the trailing backslash as well
author Chris Cannam
date Tue, 30 Sep 2014 16:23:00 +0100
parents 90220f7249fc
children
rev   line source
tp@0 1 function [edgedifflist,startandendpoints,prespeclist,midspeclist,postspeclist,validIScoords,validIRcoords,listguide,...
tp@0 2 listofallspecs] = EDB1diff2ISES(eddatafile,S,R,ivNdiffmatrix,...
tp@0 3 lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,...
tp@0 4 ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlistin,validEDIRcoords,edgeplaneperptoplane1,edgeplaneperptoplane2)
tp@0 5 % EDB1diff2ISES - Gives list of paths that includes a 2nd-order diff. combination.
tp@0 6 % EDB1diff2ISES gives the list of possible diffraction paths that includes one
tp@0 7 % second-order diffraction path, and possibly specular reflections before
tp@0 8 % and after.
tp@0 9 %
tp@0 10 % Input parameters:
tp@0 11 % eddatafile,S,R,ivNdiffmatrix,...
tp@0 12 % lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,...
tp@0 13 % ndecimaldivider,edgeplaneperptoplane1,edgeplaneperptoplane2
tp@0 14 % All these are taken from the ISEStreefile or the
tp@0 15 % setup file.
tp@0 16 % edgedifflistin List [nd1combs,1] of edges involved in all
tp@0 17 % first-order diffraction combs that have already
tp@0 18 % been found
tp@0 19 % postspeclistin Matrix [nd1combs,specorder-1] of all specular reflections
tp@0 20 % following the diffraction for all first-order diff combs
tp@0 21 % that have already been found
tp@0 22 % bigedgeweightlistin List [nd1combs,1] of visibility for edges involved in all
tp@0 23 % first-order diffraction combs that have already been found
tp@0 24 % validEDIRcoords Matrix [nd1combs,3] of image receiver coordinates for all
tp@0 25 % first-order diff combs that have already been found
tp@0 26 % GLOBAL parameters:
tp@0 27 % POTENTIALISES,ISCOORDS,ORIGINSFROM,ISESVISIBILITY,REFLORDER See EDB1findISEStree
tp@0 28 % SHOWTEXT JJ JJnumbofchars See EDB1mainISES
tp@0 29 %
tp@0 30 % Output parameters:
tp@0 31 % edgedifflist List [ncombs,2] of the edge numbers involved in each
tp@0 32 % spec-diff-diff-spec combination.
tp@0 33 % startandendpoints Matrix [ncombs,4] of the relative start and end
tp@0 34 % points of each edge. The values, [0,1], indicate
tp@0 35 % which part of the two edges that are visible.
tp@0 36 % prespeclist Matrix [ncombs,specorder-2] of the specular
tp@0 37 % reflections that precede every diffraction.
tp@0 38 % midspeclist Matrix [ncombs,specorder-2] of the specular
tp@0 39 % reflections inbetween the two diffractions.
tp@0 40 % postspeclist Matrix [ncombs,specorder-2] of the specular
tp@0 41 % reflections that follow every diffraction.
tp@0 42 % validIScoords Matrix [ncombs,3] of the image source for each
tp@0 43 % multiple-spec that precedes the diffraction. If
tp@0 44 % there is no spec refl before the diffraction, the
tp@0 45 % value [0 0 0] is given.
tp@0 46 % validIRcoords Matrix [ncombs,3] of the image receiver for each
tp@0 47 % multiple-spec that follows the diffraction. If
tp@0 48 % there is no spec refl after the diffraction, the
tp@0 49 % value [0 0 0] is given.
tp@0 50 % listguide Matrix [nuniquecombs,3] which for each row gives
tp@0 51 % 1. The number of examples in edgefdifflist etc that
tp@0 52 % are the same type of spec-diff-diff-spec comb.
tp@0 53 % 2. The first row number and 3. The last row number.
tp@0 54 % listofallspecs Matrix [nuniquecombs,3] which for each row gives
tp@0 55 % 1. The number of pre-specular reflections for the spec-diff-spec-diff-spec comb
tp@0 56 % in the same row in listguide.
tp@0 57 % 2. The number of mid-specular reflections for the spec-diff-spec-diff-spec comb
tp@0 58 % in the same row in listguide.
tp@0 59 % 3. The number of post-specular reflections for the spec-diff-spec-diff-spec comb
tp@0 60 % in the same row in listguide.
tp@0 61 %
tp@0 62 % Uses functions EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths
tp@0 63 %
tp@0 64 % ----------------------------------------------------------------------------------------------
tp@0 65 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 66 %
tp@0 67 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 68 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 69 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 70 %
tp@0 71 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 72 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 73 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 74 %
tp@0 75 % You should have received a copy of the GNU General Public License along with the
tp@0 76 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 77 % ----------------------------------------------------------------------------------------------
tp@0 78 % Peter Svensson (svensson@iet.ntnu.no) 20050202
tp@0 79 %
tp@0 80 % [edgedifflist,startandendpoints,prespeclist,midspeclist,postspeclist,validIScoords,validIRcoords,listguide,...
tp@0 81 % listofallspecs] = EDB1diff2ISES(eddatafile,S,R,...
tp@0 82 % ivNdiffmatrix,lengthNdiffmatrix,specorder,visplanesfromR,vispartedgesfromS,...
tp@0 83 % vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlist,validEDIRcoords,edgeplaneperptoplane1,edgeplaneperptoplane2)
tp@0 84
tp@0 85 global SHOWTEXT JJ JJnumbofchars
tp@0 86 global POTENTIALISES ISCOORDS ORIGINSFROM ISESVISIBILITY REFLORDER
tp@0 87
tp@0 88 eval(['load ',eddatafile])
tp@0 89
tp@0 90 [nedges,slask] = size(planesatedge);
tp@0 91 [nplanes,slask] = size(planecorners);
tp@0 92 multfac = 10^(ceil(log10(max([nedges nplanes]))));
tp@0 93
tp@0 94 edgedifflist = [];
tp@0 95 prespeclist = [];
tp@0 96 midspeclist = [];
tp@0 97 postspeclist = [];
tp@0 98 startandendpoints = [];
tp@0 99 bigedgeweightlist = [];
tp@0 100 validIScoords = [];
tp@0 101 validIRcoords = [];
tp@0 102
tp@0 103 [n1,n2] = size(POTENTIALISES);
tp@0 104 if n2 < specorder
tp@0 105 specorder = n2;
tp@0 106 end
tp@0 107
tp@0 108 maxvisibilityvalue = 2^nedgesubs-1;
tp@0 109 zerosvec1 = zeros(1,specorder-2);
tp@0 110 zerosvec2 = zeros(1,3);
tp@0 111 listguide = zeros(specorder*2-1,3);
tp@0 112 bitmultvec = 2.^[0:nedgesubs-1];
tp@0 113
tp@0 114 obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0 115 onesvec = ones(1,nedgesubs);
tp@0 116 onesvec3 = ones(1,3);
tp@0 117
tp@0 118 npostspecs = sum(double(postspeclistin.'>0)).';
tp@0 119 planesatedge = double(planesatedge);
tp@0 120
tp@0 121 if specorder >= 3
tp@0 122 coplanarsviaflatedge = sparse(zeros(nplanes,nplanes));
tp@0 123 listofflatedges = find(closwedangvec==pi);
tp@0 124 if ~isempty(listofflatedges)
tp@0 125 ivreftomatrix = planesatedge(listofflatedges,1) + (planesatedge(listofflatedges,2)-1)*nplanes;
tp@0 126 coplanarsviaflatedge(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 127 coplanarsviaflatedge = sign(coplanarsviaflatedge + coplanarsviaflatedge.');
tp@0 128 end
tp@0 129 end
tp@0 130
tp@0 131 if specorder >= 4
tp@0 132 cateyeplanecombs = sparse(zeros(nplanes,nplanes));
tp@0 133 listof90edges = find(closwedangvec==3*pi/2);
tp@0 134 if ~isempty(listof90edges)
tp@0 135 ivreftomatrix = planesatedge(listof90edges,1) + (planesatedge(listof90edges,2)-1)*nplanes;
tp@0 136 cateyeplanecombs(ivreftomatrix) = ones(size(ivreftomatrix));
tp@0 137 cateyeplanecombs = sign(cateyeplanecombs + cateyeplanecombs.');
tp@0 138 end
tp@0 139
tp@0 140 end
tp@0 141
tp@0 142 % ###########################################
tp@0 143 % # #
tp@0 144 % # S - edge - edge - R cases #
tp@0 145 % # #
tp@0 146 % ###########################################
tp@0 147 %
tp@0 148 % Possible edges for S-E-E-R are seen (at least partly) by the source and by the receiver.
tp@0 149 %
tp@0 150 % The visibility by the source is taken care of by the findISEStree
tp@0 151 % so we must check which ones are visible by the receiver.
tp@0 152 % Also, the active edge segment(s) must be selected but this is done
tp@0 153 % further down together with the S-spec-spec-edge-R cases
tp@0 154
tp@0 155 ivNdiff = ivNdiffmatrix(1:lengthNdiffmatrix(2),2);
tp@0 156 ivsinglediff = ivNdiffmatrix(1:lengthNdiffmatrix(1),1);
tp@0 157
tp@0 158 ivSEER = ivNdiff(find(REFLORDER(ivNdiff)==2));
tp@0 159
tp@0 160 possibleedgepairs = double(POTENTIALISES(ivSEER,1:2))-nplanes;
tp@0 161 ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0);
tp@0 162 if ~isempty(ivnotvisiblefromr)
tp@0 163 possibleedgepairs(ivnotvisiblefromr,:) = [];
tp@0 164 end
tp@0 165
tp@0 166 edgedifflist = [edgedifflist;possibleedgepairs];
tp@0 167 bigedgeweightlist = [bigedgeweightlist;[vispartedgesfromS(edgedifflist(:,1)) vispartedgesfromR(edgedifflist(:,2))]];
tp@0 168
tp@0 169 [nedgesadded,slask] = size(edgedifflist);
tp@0 170 zerosvec3 = zeros(nedgesadded,1);
tp@0 171 ndiffonly = nedgesadded;
tp@0 172
tp@0 173 prespeclist = [prespeclist; zerosvec3(:,ones(1,max(specorder-2,1)))];
tp@0 174 midspeclist = [midspeclist; zerosvec3(:,ones(1,max(specorder-2,1)))];
tp@0 175 postspeclist = [postspeclist; zerosvec3(:,ones(1,max(specorder-2,1)))];
tp@0 176 validIScoords = [validIScoords;S(ones(nedgesadded,1),:)];
tp@0 177 validIRcoords = [validIRcoords;R(ones(nedgesadded,1),:)];
tp@0 178
tp@0 179 if SHOWTEXT >= 3
tp@0 180 disp([' ',int2str(nedgesadded),' double diff valid'])
tp@0 181 end
tp@0 182
tp@0 183 % ###########################################
tp@0 184 % # #
tp@0 185 % # S - spec - edge - edge - R cases #
tp@0 186 % # #
tp@0 187 % # Prespec cases #
tp@0 188 % # #
tp@0 189 % ###########################################
tp@0 190 %
tp@0 191 % Possible edges for S-spec-E-E-R are seen (at least partly) by the receiver.
tp@0 192 %
tp@0 193 % The visibility doesn't need to be checked since the source-to-edge paths
tp@0 194 % were checked in the ISEStree, and the visibility from the receiver also
tp@0 195 % has been checked.
tp@0 196
tp@0 197 % The vector ivmultidiff will always refer to the original data vector
tp@0 198 % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc
tp@0 199 %
tp@0 200 % We should remove the combinations that involve an edge which the
tp@0 201 % receiver can not see, but since the edge number is in different columns
tp@0 202 % we do that in the for loop.
tp@0 203
tp@0 204 % The ii-loop will go through: spec-diff, spec-spec-diff
tp@0 205 % spec-spec-spec-diff etc
tp@0 206
tp@0 207 for ii = 1:specorder-2
tp@0 208
tp@0 209 if SHOWTEXT >= 3
tp@0 210 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before the double edge diff'])
tp@0 211 end
tp@0 212
tp@0 213 % Select the combinations where the reflection order == ii+2
tp@0 214 % which means ii specular reflections before the diffraction
tp@0 215
tp@0 216 iv = find(REFLORDER(ivNdiff) == ii+2 & POTENTIALISES(ivNdiff,ii+1)>nplanes & POTENTIALISES(ivNdiff,ii+2)>nplanes);
tp@0 217 masterivlist = ivNdiff(iv);
tp@0 218 possibleedgepairs = double(POTENTIALISES(masterivlist,ii+1:ii+2)) - nplanes;
tp@0 219
tp@0 220 % Keep only combinations for which the receiver can see the edge
tp@0 221
tp@0 222 ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0);
tp@0 223 if ~isempty(ivnotvisiblefromr)
tp@0 224 masterivlist(ivnotvisiblefromr) = [];
tp@0 225 possibleedgepairs(ivnotvisiblefromr,:) = [];
tp@0 226 end
tp@0 227
tp@0 228 possiblecombs = POTENTIALISES(masterivlist,1:ii);
tp@0 229
tp@0 230 edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgepairs(:,2))];
tp@0 231
tp@0 232 nposs = length(masterivlist);
tp@0 233
tp@0 234 if SHOWTEXT >= 3
tp@0 235 disp([' ',int2str(nposs),' IS+edge pairs valid'])
tp@0 236 end
tp@0 237
tp@0 238 edgedifflist = [edgedifflist;possibleedgepairs];
tp@0 239 prespeclist = [prespeclist;[possiblecombs zeros(nposs,specorder-2-ii)]];
tp@0 240 midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)];
tp@0 241 postspeclist = [postspeclist;zerosvec1(ones(nposs,1),:)];
tp@0 242 bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0 243
tp@0 244 % NB! It is correct below that the indices for the ISCOORDS should be
tp@0 245 % ORIGINSFROM(ORIGINSFROM(masterivlist)), rather than masterivlist.
tp@0 246 % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0 247 % spec-spec-...-diff-diff combinations and then
tp@0 248 % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0 249 % ends with a diff has no image source.
tp@0 250 % Also, two recursive references are needed since we need to get back
tp@0 251 % through the two last diffractions to reach the last specular
tp@0 252 % reflection.
tp@0 253
tp@0 254 validIScoords = [validIScoords;ISCOORDS(ORIGINSFROM(ORIGINSFROM(masterivlist)),:)];
tp@0 255 validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0 256
tp@0 257 end
tp@0 258
tp@0 259 % #######################################################################
tp@0 260 % #######################################################################
tp@0 261 % #######################################################################
tp@0 262 % #######################################################################
tp@0 263 % #######################################################################
tp@0 264
tp@0 265 % #############################################
tp@0 266 % # #
tp@0 267 % # S - edge - edge -spec - R cases #
tp@0 268 % # #
tp@0 269 % # Postspec cases #
tp@0 270 % # #
tp@0 271 % #############################################
tp@0 272 %
tp@0 273 % Possible edges for S-E-E-spec-spec-R are seen (at least partly) by the receiver.
tp@0 274 %
tp@0 275 % For some of the combos, we have the IR coords and the edge visibility
tp@0 276 % from the single-diffraction run. For potential combos that were not
tp@0 277 % included there, they are either invisible/obstructed or were not tested.
tp@0 278 % We can figure out which ones were tested because they can be found in the
tp@0 279 % POTENTIALISES under edge-spec-spec but not in the final valid list.
tp@0 280
tp@0 281 % The vector masterivlist will always refer to the original data vector
tp@0 282 % i.e. PotentialIR, OriginsfromIR, reflorderIR etc
tp@0 283 %
tp@0 284 % First we pick out those indices where there was a single diffraction, but
tp@0 285 % skip those with only diffraction (because we dealt with them already).
tp@0 286 % Also, select only those where the diffraction is the first in the sequence
tp@0 287 % of reflections.
tp@0 288
tp@0 289 for ii = 1:specorder-2
tp@0 290
tp@0 291 if SHOWTEXT >= 3
tp@0 292 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the double edge diff'])
tp@0 293 end
tp@0 294
tp@0 295 % Select the combinations where the reflection order == ii+2
tp@0 296 % (which means ii specular reflections in addition to the diffraction)
tp@0 297 % and where the first two columns in POTENTIALISES contain edges.
tp@0 298
tp@0 299 iv = find(REFLORDER(ivNdiff) == ii+2 & POTENTIALISES(ivNdiff,1)>nplanes & POTENTIALISES(ivNdiff,2)>nplanes);
tp@0 300 masterivlist = ivNdiff(iv);
tp@0 301 possibleedgepairs = double(POTENTIALISES(masterivlist,1:2)) - nplanes;
tp@0 302 possiblecombs = POTENTIALISES(masterivlist,3:2+ii);
tp@0 303 possibleweights = ISESVISIBILITY(masterivlist);
tp@0 304
tp@0 305 % Compare with those combinations that were found OK
tp@0 306 % in the first-order diffraction step (EDB1diffISES),
tp@0 307 % and that were input matrices to EDB1diff2ISES.
tp@0 308 % The index vector ivOK refers to these input matrices
tp@0 309 % (edgedifflistin,posspeclistin,validEDIRcoords,bigedgeweightlistin)
tp@0 310 % npostspecs is a list that was calculated inside EDB1diff2ISES but
tp@0 311 % refers to the input matrices.
tp@0 312
tp@0 313 ivOK = find(npostspecs==ii);
tp@0 314 if ~isempty(ivOK)
tp@0 315 patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:ii)];
tp@0 316
tp@0 317 % Find out which ones, of all the possible first-order diffraction combos
tp@0 318 % in POTENTIALISES, that were indeed tested and found
tp@0 319 % invisible/obstructed in EDB1diffISES.
tp@0 320
tp@0 321 ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == ii+1));
tp@0 322 patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+ii))];
tp@0 323 if ~isempty(patternOK) & ~isempty(patternALL)
tp@0 324 patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0 325 else
tp@0 326 if isempty(patternOK)
tp@0 327 patternNOTOK = patternALL;
tp@0 328 else % Then patternALL must be empty
tp@0 329 patternNOTOK = [];
tp@0 330 end
tp@0 331 end
tp@0 332
tp@0 333 % Now, the patterns in patternNOTOK can be removed from
tp@0 334 % masterivlist.
tp@0 335
tp@0 336 patterntocompare = [possibleedgepairs(:,2) possiblecombs(:,1:ii)];
tp@0 337
tp@0 338 ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0 339 masterivlist(ivtocancel) = [];
tp@0 340 possibleedgepairs(ivtocancel,:) = [];
tp@0 341 possiblecombs(ivtocancel,:) = [];
tp@0 342 possibleweights(ivtocancel,:) = [];
tp@0 343 patterntocompare(ivtocancel,:) = [];
tp@0 344
tp@0 345 [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0 346 ivmustbechecked = find(ivcompletelyOK==0);
tp@0 347 ivcompletelyOK = find(ivcompletelyOK);
tp@0 348
tp@0 349 if ~isempty(ivmustbechecked)
tp@0 350 masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0 351 ntocheckmore = length(masterlisttocheckmore);
tp@0 352
tp@0 353 %----------------------------------------------
tp@0 354 % Must carry out a visibility and obstruction check for the special
tp@0 355 % combinations here.
tp@0 356 % These combinations have a postspec-combination that hasn't been
tp@0 357 % encountered in the single diffraction cases, so no visibility
tp@0 358 % test has been made for these.
tp@0 359
tp@0 360 lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,2))-nplanes;
tp@0 361 newIRcoords = R;
tp@0 362 reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0 363 for jj = 1:ii
tp@0 364 reflplanes = POTENTIALISES(masterlisttocheckmore,3+ii-jj);
tp@0 365 reflplanesexpand(:,jj) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0 366 newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,1,onesvec3);
tp@0 367 newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0 368 eval(['newIRcoords',JJ(jj,1:JJnumbofchars(jj)),' = newIRcoordsexpand;'])
tp@0 369 end
tp@0 370 [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0 371 tocoords = toedgecoords;
tp@0 372 lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0 373 lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0 374 masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0 375 masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0 376
tp@0 377 ntocheckmore = length(masterlisttocheckmore);
tp@0 378 if SHOWTEXT >= 3
tp@0 379 disp([' ',int2str(ntocheckmore),' special edge+edge+IR combinations to check'])
tp@0 380 end
tp@0 381
tp@0 382 for jj = ii:-1:1
tp@0 383 if length(masterlisttocheckmore) > 0
tp@0 384 eval(['fromcoords = newIRcoords',JJ(jj,1:JJnumbofchars(jj)),';']);
tp@0 385 if jj < ii
tp@0 386 eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])
tp@0 387 end
tp@0 388
tp@0 389 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,jj),4),planenvecs(reflplanesexpand(:,jj),:),minvals(reflplanesexpand(:,jj),:),...
tp@0 390 maxvals(reflplanesexpand(:,jj),:),planecorners(reflplanesexpand(:,jj),:),corners,ncornersperplanevec(reflplanesexpand(:,jj)));
tp@0 391 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 392 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 393 disp(' handled correctly yet.')
tp@0 394 end
tp@0 395 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0 396
tp@0 397 masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0 398 edgeweightlist = edgeweightlist(hitplanes);
tp@0 399 lastedgenumbers = lastedgenumbers(hitplanes);
tp@0 400 reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0 401 toedgecoords = toedgecoords(hitplanes,:);
tp@0 402
tp@0 403 for kk = 1:ii
tp@0 404 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0 405 if kk > jj
tp@0 406 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0 407 end
tp@0 408 end
tp@0 409 ntocheckmore = length(masterlisttocheckmore);
tp@0 410
tp@0 411 end
tp@0 412
tp@0 413 if SHOWTEXT >= 3
tp@0 414 disp([' ',int2str(ntocheckmore),' of the special edge+edge+IR combinations survived the visibility test in refl plane ',int2str(jj)])
tp@0 415 end
tp@0 416 end
tp@0 417
tp@0 418 % Obstruction test of all the involved paths: R ->
tp@0 419 % reflplane1 -> reflplane2 -> ... -> last edge
tp@0 420
tp@0 421 if obstructtestneeded & ntocheckmore > 0
tp@0 422
tp@0 423 for jj = 1:ii+1
tp@0 424 if ntocheckmore > 0
tp@0 425 if jj == 1
tp@0 426 fromcoords = R;
tp@0 427 startplanes = [];
tp@0 428 else
tp@0 429 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0 430 startplanes = reflplanesexpand(:,jj-1);
tp@0 431 end
tp@0 432 if jj == ii+1,
tp@0 433 tocoords = toedgecoords;
tp@0 434 endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0 435 else
tp@0 436 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
tp@0 437 endplanes = reflplanesexpand(:,jj);
tp@0 438 end
tp@0 439
tp@0 440 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 441 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 442
tp@0 443 if nobstructions > 0
tp@0 444 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0 445 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 446 lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0 447 reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0 448 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 449 for kk = 1:ii
tp@0 450 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0 451 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0 452 end
tp@0 453
tp@0 454 end
tp@0 455 ntocheckmore = length(masterlisttocheckmore);
tp@0 456
tp@0 457 end
tp@0 458
tp@0 459 end
tp@0 460 if SHOWTEXT >= 3
tp@0 461 disp([' ',int2str(ntocheckmore),' of the special edge+edge+IR combinations survived the obstruction test'])
tp@0 462 end
tp@0 463
tp@0 464 end
tp@0 465
tp@0 466 % Add the found special combinations to the outdata list
tp@0 467
tp@0 468 edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,1:2))-nplanes];
tp@0 469 prespeclist = [prespeclist;zerosvec1(ones(ntocheckmore,1),:)];
tp@0 470 midspeclist = [midspeclist;zerosvec1(ones(ntocheckmore,1),:)];
tp@0 471 postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-2-ii)]];
tp@0 472 bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];
tp@0 473
tp@0 474 eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(ii,1:JJnumbofchars(ii)),'];']);
tp@0 475 validIScoords = [validIScoords;S(ones(ntocheckmore,1),:)];
tp@0 476
tp@0 477 end
tp@0 478
tp@0 479 masterivlist = masterivlist(ivcompletelyOK);
tp@0 480 possibleedgepairs = possibleedgepairs(ivcompletelyOK,:);
tp@0 481 possiblecombs = possiblecombs(ivcompletelyOK,:);
tp@0 482 possibleweights = possibleweights(ivcompletelyOK,:);
tp@0 483
tp@0 484 nposs = length(ivcompletelyOK);
tp@0 485
tp@0 486 if SHOWTEXT >= 3
tp@0 487 disp([' ',int2str(nposs),' Edge+edge+IR segments (non-special combinations) survived the obstruction test'])
tp@0 488 end
tp@0 489
tp@0 490 % Add the found "standard" combinations to the outdata list
tp@0 491
tp@0 492 edgedifflist = [edgedifflist;possibleedgepairs];
tp@0 493 prespeclist = [prespeclist;zerosvec1(ones(nposs,1),:)];
tp@0 494 midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)];
tp@0 495 postspeclist = [postspeclist;[possiblecombs zeros(nposs,specorder-2-ii)]];
tp@0 496 bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];
tp@0 497 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0 498 validIScoords = [validIScoords;S(ones(nposs,1),:)];
tp@0 499
tp@0 500 end
tp@0 501 end
tp@0 502
tp@0 503
tp@0 504 % #######################################################################
tp@0 505 % #######################################################################
tp@0 506 % #######################################################################
tp@0 507 % #######################################################################
tp@0 508 % #######################################################################
tp@0 509
tp@0 510 % ####################################################
tp@0 511 % # #
tp@0 512 % # S - spec - edge - edge - spec - R cases #
tp@0 513 % # #
tp@0 514 % # Pre and postspec cases #
tp@0 515 % # #
tp@0 516 % ####################################################
tp@0 517
tp@0 518
tp@0 519 % The ii- and jj-loops will go through all combinations of specular
tp@0 520 % reflection before and after the double diffraction.
tp@0 521 %
tp@0 522 % Unlike before, we will look in the list of already found combinations
tp@0 523 % (edgedifflist etc)
tp@0 524
tp@0 525 for ii = 1:specorder-3
tp@0 526
tp@0 527 for jj = 1:specorder-ii-2
tp@0 528
tp@0 529 if SHOWTEXT >= 3
tp@0 530 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 531 end
tp@0 532
tp@0 533 % Select the combinations where the reflection order == ii+jj+2
tp@0 534 % (which means ii+jj specular reflections in addition to the
tp@0 535 % diffraction), and where columns ii+1 and ii+2 of POTENTIALISES
tp@0 536 % contain edges.
tp@0 537
tp@0 538 iv = find(REFLORDER(ivNdiff) == ii+jj+2 & POTENTIALISES(ivNdiff,ii+1)>nplanes & POTENTIALISES(ivNdiff,ii+2)>nplanes);
tp@0 539 masterivlist = ivNdiff(iv);
tp@0 540 possibleedgepairs = double(POTENTIALISES(masterivlist,ii+1:ii+2)) - nplanes;
tp@0 541 possibleprespecs = POTENTIALISES(masterivlist,1:ii);
tp@0 542 possiblepostspecs = POTENTIALISES(masterivlist,ii+3:ii+2+jj);
tp@0 543 possibleweights = ISESVISIBILITY(masterivlist);
tp@0 544
tp@0 545 % Compare with those that have already been found OK
tp@0 546 ivOK = find(npostspecs==jj);
tp@0 547 if ~isempty(ivOK)
tp@0 548 patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:jj)];
tp@0 549 else
tp@0 550 patternOK = [];
tp@0 551 end
tp@0 552
tp@0 553 % Find out which ones have been checked and found invisible/obstructed
tp@0 554 ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == jj+1));
tp@0 555 patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+jj))];
tp@0 556 if ~isempty(patternOK) & ~isempty(patternALL)
tp@0 557 patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0 558 else
tp@0 559 if isempty(patternOK)
tp@0 560 patternNOTOK = patternALL;
tp@0 561 else % Then patternALL must be empty
tp@0 562 patternNOTOK = [];
tp@0 563 end
tp@0 564 end
tp@0 565
tp@0 566 patterntocompare = [possibleedgepairs(:,2) possiblepostspecs(:,1:jj)];
tp@0 567
tp@0 568 ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0 569 masterivlist(ivtocancel) = [];
tp@0 570 possibleedgepairs(ivtocancel,:) = [];
tp@0 571 possibleprespecs(ivtocancel,:) = [];
tp@0 572 possiblepostspecs(ivtocancel,:) = [];
tp@0 573 possibleweights(ivtocancel,:) = [];
tp@0 574 patterntocompare(ivtocancel,:) = [];
tp@0 575
tp@0 576 [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0 577 ivmustbechecked = find(ivcompletelyOK==0);
tp@0 578 ivcompletelyOK = find(ivcompletelyOK);
tp@0 579 if ~isempty(ivmustbechecked)
tp@0 580 masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0 581 ntocheckmore = length(masterlisttocheckmore);
tp@0 582
tp@0 583 %----------------------------------------------
tp@0 584 % Must carry out a visibility and obstruction check for the special
tp@0 585 % combinations here.
tp@0 586 %
tp@0 587 % NB! toedgecoords are the coordinates of the last edge in the sequence.
tp@0 588 % This name is because for post-specular reflections, the propagation
tp@0 589 % is viewed from the receiver towards the last edge!
tp@0 590
tp@0 591 lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,ii+2))-nplanes;
tp@0 592 newIRcoords = R;
tp@0 593 reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0 594 for kk = 1:jj
tp@0 595 reflplanes = POTENTIALISES(masterlisttocheckmore,3+ii+jj-kk);
tp@0 596 reflplanesexpand(:,kk) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0 597 newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,1,onesvec3);
tp@0 598 newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0 599 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoordsexpand;'])
tp@0 600 end
tp@0 601 [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0 602 tocoords = toedgecoords;
tp@0 603 lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0 604 lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0 605 masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0 606 masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0 607
tp@0 608 ntocheckmore = length(masterlisttocheckmore);
tp@0 609 if SHOWTEXT >= 3
tp@0 610 disp([' ',int2str(ntocheckmore),' special IS+edge+edge+IR combinations to check'])
tp@0 611 end
tp@0 612
tp@0 613 for kk = jj:-1:1
tp@0 614 if length(masterlisttocheckmore) > 0
tp@0 615 eval(['fromcoords = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),';']);
tp@0 616 if kk < jj
tp@0 617 eval(['tocoords = reflpoints',JJ(kk+1,1:JJnumbofchars(kk+1)),';'])
tp@0 618 end
tp@0 619
tp@0 620 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,kk),4),planenvecs(reflplanesexpand(:,kk),:),minvals(reflplanesexpand(:,kk),:),...
tp@0 621 maxvals(reflplanesexpand(:,kk),:),planecorners(reflplanesexpand(:,kk),:),corners,ncornersperplanevec(reflplanesexpand(:,kk)));
tp@0 622 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 623 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 624 disp(' handled correctly yet.')
tp@0 625 end
tp@0 626 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints;'])
tp@0 627
tp@0 628 masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0 629 edgeweightlist = edgeweightlist(hitplanes);
tp@0 630 lastedgenumbers = lastedgenumbers(hitplanes);
tp@0 631 reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0 632 toedgecoords = toedgecoords(hitplanes,:);
tp@0 633
tp@0 634 for ll = 1:jj
tp@0 635 eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0 636 if ll > kk
tp@0 637 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0 638 end
tp@0 639 end
tp@0 640 ntocheckmore = length(masterlisttocheckmore);
tp@0 641
tp@0 642 end
tp@0 643
tp@0 644 if SHOWTEXT >= 3
tp@0 645 disp([' ',int2str(ntocheckmore),' of the special IS+edge+edge+IR combinations survived the visibility test in refl plane ',int2str(jj)])
tp@0 646 end
tp@0 647
tp@0 648 end
tp@0 649
tp@0 650 % Obstruction test of all the involved paths: R ->
tp@0 651 % reflplane1 -> reflplane2 -> ... -> last edge
tp@0 652
tp@0 653 if obstructtestneeded & ntocheckmore > 0
tp@0 654
tp@0 655 for kk = 1:jj+1
tp@0 656 if ntocheckmore > 0
tp@0 657 if kk == 1
tp@0 658 fromcoords = R;
tp@0 659 startplanes = [];
tp@0 660 else
tp@0 661 eval(['fromcoords = reflpoints',JJ(kk-1,1:JJnumbofchars(kk-1)),';'])
tp@0 662 startplanes = reflplanesexpand(:,kk-1);
tp@0 663 end
tp@0 664 if kk == jj+1,
tp@0 665 tocoords = toedgecoords;
tp@0 666 endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0 667 else
tp@0 668 eval(['tocoords = reflpoints',JJ(kk,1:JJnumbofchars(kk)),';'])
tp@0 669 endplanes = reflplanesexpand(:,kk);
tp@0 670 end
tp@0 671
tp@0 672 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 673 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 674
tp@0 675 if nobstructions > 0
tp@0 676 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0 677 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 678 lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0 679 reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0 680 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 681 for ll = 1:jj
tp@0 682 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0 683 eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0 684 end
tp@0 685
tp@0 686 end
tp@0 687 ntocheckmore = length(masterlisttocheckmore);
tp@0 688
tp@0 689 end
tp@0 690
tp@0 691 end
tp@0 692 if SHOWTEXT >= 3
tp@0 693 disp([' ',int2str(ntocheckmore),' of the special IS+edge+edge+IR combinations survived the obstruction test'])
tp@0 694 end
tp@0 695
tp@0 696 end
tp@0 697
tp@0 698 % Add the found special combinations to the outdata list
tp@0 699
tp@0 700 edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,ii+1:ii+2))-nplanes];
tp@0 701 prespeclist = [prespeclist;[double(POTENTIALISES(masterlisttocheckmore,1:ii)) zeros(ntocheckmore,specorder-2-ii)]];
tp@0 702 midspeclist = [midspeclist;zerosvec1(ones(ntocheckmore,1),:)];
tp@0 703 postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-2-ii)]];
tp@0 704 bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];
tp@0 705 eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(jj,1:JJnumbofchars(jj)),'];']);
tp@0 706 % NB!! In the same way as earlier, we must a recursive reference
tp@0 707 % method to find the image source of the last specular reflection.
tp@0 708 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterlisttocheckmore)));
tp@0 709 for kk = 2:jj
tp@0 710 ivref = ORIGINSFROM(ivref);
tp@0 711 end
tp@0 712 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 713 end
tp@0 714
tp@0 715 masterivlist = masterivlist(ivcompletelyOK);
tp@0 716 possibleedgepairs = possibleedgepairs(ivcompletelyOK,:);
tp@0 717 possibleprespecs = possibleprespecs(ivcompletelyOK,:);
tp@0 718 possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0 719 possibleweights = possibleweights(ivcompletelyOK,:);
tp@0 720
tp@0 721 nposs = length(ivcompletelyOK);
tp@0 722
tp@0 723 if SHOWTEXT >= 3
tp@0 724 disp([' ',int2str(nposs),' IS+Edge+edge+IR segments (non-special) survived the obstruction test'])
tp@0 725 end
tp@0 726
tp@0 727 % Add the found "standard" combinations to the outdata list
tp@0 728
tp@0 729 edgedifflist = [edgedifflist;possibleedgepairs];
tp@0 730 prespeclist = [prespeclist; [possibleprespecs zeros(nposs,specorder-2-ii)]];
tp@0 731 midspeclist = [midspeclist;zerosvec1(ones(nposs,1),:)];
tp@0 732 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-2-jj)]];
tp@0 733 bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];
tp@0 734 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0 735 % NB!! In the same way as earlier, we must a recursive reference
tp@0 736 % method to find the image source of the last specular reflection.
tp@0 737 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist)));
tp@0 738 for kk = 2:jj
tp@0 739 ivref = ORIGINSFROM(ivref);
tp@0 740 end
tp@0 741 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 742
tp@0 743 end
tp@0 744
tp@0 745 end
tp@0 746
tp@0 747 % #######################################################################
tp@0 748 % #######################################################################
tp@0 749 % #######################################################################
tp@0 750 % #######################################################################
tp@0 751 % #######################################################################
tp@0 752
tp@0 753 % #######################################################################
tp@0 754 % #
tp@0 755 % # Midspec cases, with prespecs and postspecs
tp@0 756 % #
tp@0 757 % # S - spec - edge - spec - edge - R
tp@0 758 % #
tp@0 759 % #######################################################################
tp@0 760
tp@0 761 for iipre = 0:specorder-3
tp@0 762 for jjmid = 1:specorder-2-iipre
tp@0 763 for kkpost = 0:specorder-iipre-jjmid-2
tp@0 764
tp@0 765 if SHOWTEXT >= 3
tp@0 766 if iipre > 0
tp@0 767 JJpre = JJ(iipre,1:JJnumbofchars(iipre));
tp@0 768 else
tp@0 769 JJpre = '0';
tp@0 770 end
tp@0 771 if kkpost > 0
tp@0 772 JJpost = JJ(kkpost,1:JJnumbofchars(kkpost));
tp@0 773 else
tp@0 774 JJpost = '0';
tp@0 775 end
tp@0 776 disp([' Checking for ',JJpre,' spec refl before, ',JJ(jjmid,1:JJnumbofchars(jjmid)),' spec refl between the double edge diff, and'])
tp@0 777 disp([' ',JJpost,' spec refl after the double edge diff.'])
tp@0 778 end
tp@0 779 iv = find(REFLORDER(ivNdiff) == iipre+jjmid+kkpost+2 & POTENTIALISES(ivNdiff,iipre+1)>nplanes & POTENTIALISES(ivNdiff,iipre+jjmid+2)>nplanes);
tp@0 780 masterivlist = ivNdiff(iv);
tp@0 781 possibleedgepairs = double(POTENTIALISES(masterivlist,[iipre+1 iipre+jjmid+2])) - nplanes;
tp@0 782 nfast = 0;
tp@0 783
tp@0 784 if kkpost == 0
tp@0 785 possiblepostspecs = [];
tp@0 786
tp@0 787 % Keep only combinations for which the receiver can see the edge
tp@0 788 ivnotvisiblefromr = find(vispartedgesfromR(possibleedgepairs(:,2))==0);
tp@0 789 if ~isempty(ivnotvisiblefromr)
tp@0 790 masterivlist(ivnotvisiblefromr) = [];
tp@0 791 possibleedgepairs(ivnotvisiblefromr,:) = [];
tp@0 792 end
tp@0 793 else
tp@0 794 possiblepostspecs = POTENTIALISES(masterivlist,iipre+jjmid+3:iipre+jjmid+2+kkpost);
tp@0 795
tp@0 796 % Compare with those that have already been found OK
tp@0 797 ivOK = find(npostspecs==kkpost);
tp@0 798 if ~isempty(ivOK)
tp@0 799 patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:kkpost)];
tp@0 800 else
tp@0 801 patternOK = [];
tp@0 802 end
tp@0 803 % Find out which ones have been checked and found invisible/obstructed
tp@0 804 ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == kkpost+1));
tp@0 805 patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+kkpost))];
tp@0 806
tp@0 807 if ~isempty(patternOK) & ~isempty(patternALL)
tp@0 808 patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0 809 else
tp@0 810 if isempty(patternOK)
tp@0 811 patternNOTOK = patternALL;
tp@0 812 else % Then patternALL must be empty
tp@0 813 patternNOTOK = [];
tp@0 814 end
tp@0 815 end
tp@0 816
tp@0 817 patterntocompare = [possibleedgepairs(:,2) possiblepostspecs(:,1:kkpost)];
tp@0 818
tp@0 819 if ~isempty(patternNOTOK)
tp@0 820 ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0 821 masterivlist(ivtocancel) = [];
tp@0 822 possibleedgepairs(ivtocancel,:) = [];
tp@0 823 possiblepostspecs(ivtocancel,:) = [];
tp@0 824 patterntocompare(ivtocancel,:) = [];
tp@0 825 end
tp@0 826
tp@0 827 [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0 828 ivmustbechecked = find(ivcompletelyOK==0);
tp@0 829 ivcompletelyOK = find(ivcompletelyOK);
tp@0 830
tp@0 831 if ~isempty(ivmustbechecked) & SHOWTEXT > 0
tp@0 832 disp('WARNING!! For midspec and postspec case, all checks have not been implemented yet!!')
tp@0 833 end
tp@0 834
tp@0 835 masterivlist = masterivlist(ivcompletelyOK);
tp@0 836 possibleedgepairs = possibleedgepairs(ivcompletelyOK,:);
tp@0 837 possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0 838
tp@0 839 nposs = length(ivcompletelyOK);
tp@0 840
tp@0 841 end
tp@0 842
tp@0 843 if iipre > 0
tp@0 844 possibleprespecs = POTENTIALISES(masterivlist,1:iipre);
tp@0 845 else
tp@0 846 possibleprespecs = [];
tp@0 847 end
tp@0 848
tp@0 849 % NB!! possiblemidspecs is numbered in reverse order
tp@0 850 % since we view the propagation by starting to mirror the last edge
tp@0 851 % and move towards the first edge.
tp@0 852
tp@0 853 possiblemidspecs = POTENTIALISES(masterivlist,iipre+1+jjmid:-1:iipre+2);
tp@0 854
tp@0 855 if kkpost > 0
tp@0 856 edgeweightlist = [ISESVISIBILITY(masterivlist) bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))];
tp@0 857 else
tp@0 858 edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgepairs(:,2))];
tp@0 859 end
tp@0 860
tp@0 861 % Expand the various lists and matrices to represent the
tp@0 862 % sub-divided edges.
tp@0 863
tp@0 864 nposs = length(masterivlist);
tp@0 865 if nposs > 0
tp@0 866 if iipre == 1
tp@0 867 possibleprespecs = reshape(possibleprespecs(:,onesvec).',nposs*nedgesubs,1);
tp@0 868 elseif iipre > 1
tp@0 869 possibleprespecs = reshape(repmat(possibleprespecs.',nedgesubs,1),iipre,nposs*nedgesubs).';
tp@0 870 end
tp@0 871 if jjmid == 1
tp@0 872 possiblemidspecs = reshape(possiblemidspecs(:,onesvec).',nposs*nedgesubs,1);
tp@0 873 elseif jjmid > 1
tp@0 874 possiblemidspecs = reshape(repmat(possiblemidspecs.',nedgesubs,1),jjmid,nposs*nedgesubs).';
tp@0 875 end
tp@0 876 if kkpost == 1
tp@0 877 possiblepostspecs = reshape(possiblepostspecs(:,onesvec).',nposs*nedgesubs,1);
tp@0 878 elseif kkpost > 1
tp@0 879 possiblepostspecs = reshape(repmat(possiblepostspecs.',nedgesubs,1),kkpost,nposs*nedgesubs).';
tp@0 880 end
tp@0 881
tp@0 882 expandedmasterivlist = reshape(masterivlist(:,onesvec).',nposs*nedgesubs,1);
tp@0 883 if kkpost > 0
tp@0 884 expandedivcompletelyOK = reshape(ivcompletelyOK(:,onesvec).',nposs*nedgesubs,1);
tp@0 885 end
tp@0 886
tp@0 887 edgeweightlist = reshape(repmat(edgeweightlist.',[nedgesubs 1]),2,nposs*nedgesubs).';
tp@0 888
tp@0 889 for ll = 1:nedgesubs
tp@0 890 edgeweightlist(ll:nedgesubs:end,1) = double(bitget(edgeweightlist(ll:nedgesubs:end,1),ll))*bitmultvec(ll);
tp@0 891 edgeweightlist(ll:nedgesubs:end,2) = double(bitget(edgeweightlist(ll:nedgesubs:end,2),ll))*bitmultvec(ll);
tp@0 892 end
tp@0 893
tp@0 894 %----------------------------------------------
tp@0 895 % Must carry out a visibility and obstruction check for the
tp@0 896 % edge-spec-edge paths
tp@0 897 %
tp@0 898 % NB! toedgecoords are the coordinates of the first edge in the sequence
tp@0 899 % and fromedgecoords refers to the last edge, after the mid-specular
tp@0 900 % reflections.
tp@0 901 % This name is because for mid-specular reflections, the propagation
tp@0 902 % is viewed from the last edge towards the first edge!
tp@0 903
tp@0 904 [toedgecoords,firstedgeweightlist,slask] = EDB1getedgepoints(edgestartcoords(possibleedgepairs(:,2),:),edgeendcoords(possibleedgepairs(:,2),:),edgelengthvec(possibleedgepairs(:,1),:),nedgesubs);
tp@0 905 tocoords = toedgecoords;
tp@0 906 [fromedgecoords,lastedgeweightlist,slask] = EDB1getedgepoints(edgestartcoords(possibleedgepairs(:,1),:),edgeendcoords(possibleedgepairs(:,1),:),edgelengthvec(possibleedgepairs(:,2),:),nedgesubs);
tp@0 907 fromcoords = fromedgecoords;
tp@0 908
tp@0 909 possibleedgepairs = reshape(repmat(possibleedgepairs.',nedgesubs,1),2,nposs*nedgesubs).';
tp@0 910
tp@0 911 edgeimagecoords = fromedgecoords;
tp@0 912 for ll = 1:jjmid
tp@0 913 edgeimagecoords = EDB1findis(edgeimagecoords,possiblemidspecs(:,ll),planeeqs,size(fromedgecoords,1),onesvec3);
tp@0 914 eval(['bigedgeimagecoords',JJ(ll,1:JJnumbofchars(ll)),' = edgeimagecoords;'])
tp@0 915 end
tp@0 916
tp@0 917 % Some cases do not need to be checked, when jjmid = 2: the
tp@0 918 % cateye cases. For these, we will have doubles (both, e.g.
tp@0 919 % 3-7 and 7-3) and one should be tossed then. The non-doubles
tp@0 920 % can be kept in a "fast lane" so that visibility isn't
tp@0 921 % checked, but obstruction is.
tp@0 922
tp@0 923 if jjmid == 2
tp@0 924 specpattern = double(possiblemidspecs);
tp@0 925 ivreftomatrix = specpattern(:,1) + ( specpattern(:,2)-1)*nplanes;
tp@0 926 ivcateyes = find( cateyeplanecombs(ivreftomatrix) );
tp@0 927
tp@0 928 if ~isempty(ivcateyes),
tp@0 929 specpattern = specpattern(ivcateyes,:);
tp@0 930 fliporder = specpattern(:,2)<specpattern(:,1);
tp@0 931 ivfliporder = find(fliporder);
tp@0 932 specpattern(ivfliporder,:) = specpattern(ivfliporder,[2 1]);
tp@0 933 [uniquepatterns,ivec,jvec] = unique(specpattern,'rows');
tp@0 934
tp@0 935 countcases = histc(jvec,[1:max(jvec)]);
tp@0 936 ivtossone = find(fliporder & countcases(jvec)==nedgesubs*2);
tp@0 937
tp@0 938 ivfastlane = ivcateyes;
tp@0 939 ivfastlane(ivtossone) = [];
tp@0 940 if ~isempty(ivfastlane)
tp@0 941 expandedmasterivlistfast = expandedmasterivlist(ivfastlane);
tp@0 942 possibleedgepairsfast = possibleedgepairs(ivfastlane,:);
tp@0 943 fromedgecoordsfast = fromedgecoords(ivfastlane,:);
tp@0 944 toedgecoordsfast = toedgecoords(ivfastlane,:);
tp@0 945 if ~isempty(possibleprespecs)
tp@0 946 possibleprespecsfast = possibleprespecs(ivfastlane,:);
tp@0 947 end
tp@0 948 possiblemidspecsfast = possiblemidspecs(ivfastlane,:);
tp@0 949 if ~isempty(possiblepostspecs)
tp@0 950 possiblepostspecsfast = possiblepostspecs(ivfastlane,:);
tp@0 951 end
tp@0 952 edgeweightlistfast = edgeweightlist(ivfastlane,:);
tp@0 953 for mm = 1:jjmid
tp@0 954 eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'fast = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(ivfastlane,:);']);
tp@0 955 end
tp@0 956 reflpoints2fast = 0.5*(bigedgeimagecoords2fast-fromedgecoordsfast) + fromedgecoordsfast;
tp@0 957 reflpoints1fast = reflpoints2fast;
tp@0 958 nfast = length(edgeweightlistfast);
tp@0 959 end
tp@0 960
tp@0 961 if ~isempty(ivtossone)
tp@0 962 ivtossone = ivcateyes;
tp@0 963 expandedmasterivlist((ivtossone)) = [];
tp@0 964 edgeweightlist((ivtossone),:) = [];
tp@0 965 possibleedgepairs((ivtossone),:) = [];
tp@0 966 if ~isempty(possibleprespecs)
tp@0 967 possibleprespecs((ivtossone),:) = [];
tp@0 968 end
tp@0 969 possiblemidspecs((ivtossone),:) = [];
tp@0 970 if ~isempty(possiblepostspecs)
tp@0 971 possiblepostspecs((ivtossone),:) = [];
tp@0 972 end
tp@0 973 fromcoords((ivtossone),:) = [];
tp@0 974 tocoords((ivtossone),:) = [];
tp@0 975 for mm = 1:jjmid
tp@0 976 eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'((ivtossone),:) = [];']);
tp@0 977 end
tp@0 978 nposs = length(expandedmasterivlist);
tp@0 979 end
tp@0 980
tp@0 981 end
tp@0 982
tp@0 983 end
tp@0 984
tp@0 985 nposs = length(expandedmasterivlist);
tp@0 986 end
tp@0 987
tp@0 988 if SHOWTEXT >= 3
tp@0 989 if jjmid ~= 2
tp@0 990 disp([' ',int2str(nposs),' edge+spec+edge combos found initially:'])
tp@0 991 else
tp@0 992 disp([' ',int2str(nposs),' edge+spec+edge combos found initially, + ',int2str(nfast),' cateye combos'])
tp@0 993 end
tp@0 994 end
tp@0 995
tp@0 996 %--------------------------------------------------------------
tp@0 997 % Check the visibility through all the reflection planes
tp@0 998 %
tp@0 999
tp@0 1000 for ll = 1:jjmid
tp@0 1001 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = [];'])
tp@0 1002 end
tp@0 1003 for ll = jjmid:-1:1
tp@0 1004 if nposs > 0
tp@0 1005 eval(['fromcoords = bigedgeimagecoords',JJ(ll,1:JJnumbofchars(ll)),';'])
tp@0 1006 if ll < jjmid
tp@0 1007 eval(['tocoords = reflpoints',JJ(ll+1,1:JJnumbofchars(ll+1)),';'])
tp@0 1008 end
tp@0 1009
tp@0 1010 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(possiblemidspecs(:,ll),4),planenvecs(possiblemidspecs(:,ll),:),minvals(possiblemidspecs(:,ll),:),...
tp@0 1011 maxvals(possiblemidspecs(:,ll),:),planecorners(possiblemidspecs(:,ll),:),corners,ncornersperplanevec(possiblemidspecs(:,ll)));
tp@0 1012
tp@0 1013 % Make a special treatment for the cases with the
tp@0 1014 % specular reflection point right on an edge since some
tp@0 1015 % of these are special cases:
tp@0 1016 % "edgeplaneperptoplane1" indicates that edge-plane-edge
tp@0 1017 % travels along the edge's plane and is reflected at a
tp@0 1018 % 90 degree corner (which is an inactive edge).
tp@0 1019 % These are treated as good hits.
tp@0 1020 % "edgeplaneperptoplane2" indicates that edge-plane-edge
tp@0 1021 % has a specular reflection right at a flat edge
tp@0 1022 % between two coplanar planes.
tp@0 1023 % These come in pairs; one half-hit in the two
tp@0 1024 % coplanar planes. Only one in each pair should be
tp@0 1025 % kept.
tp@0 1026
tp@0 1027 if jjmid == 1 & ~isempty(edgehits)
tp@0 1028 edge1 = possibleedgepairs(edgehits,1);
tp@0 1029 edge2 = possibleedgepairs(edgehits,2);
tp@0 1030 midspec = possiblemidspecs(edgehits,1);
tp@0 1031 ivreftomatrix1 = double(midspec) + double(edge1-1)*nplanes;
tp@0 1032 ivreftomatrix2 = double(midspec) + double(edge2-1)*nplanes;
tp@0 1033
tp@0 1034 specialhit = edgeplaneperptoplane1(ivreftomatrix1).*edgeplaneperptoplane1(ivreftomatrix1);
tp@0 1035 ivspecial = find(specialhit);
tp@0 1036 if ~isempty(ivspecial)
tp@0 1037 hitplanes = [hitplanes;edgehits(ivspecial)];
tp@0 1038 reflpoints = [reflpoints;edgehitpoints(ivspecial,:)];
tp@0 1039 edgehits(ivspecial) = [];
tp@0 1040 edgehitpoints(ivspecial,:) = [];
tp@0 1041 ivreftomatrix1(ivspecial) = [];
tp@0 1042 ivreftomatrix2(ivspecial) = [];
tp@0 1043 end
tp@0 1044
tp@0 1045 specialhit = edgeplaneperptoplane2(ivreftomatrix1).*edgeplaneperptoplane2(ivreftomatrix1);
tp@0 1046 ivspecial = find(specialhit);
tp@0 1047 if ~isempty(ivspecial)
tp@0 1048 patternlist = double([possibleedgepairs(edgehits(ivspecial),1) possiblemidspecs(edgehits(ivspecial),1) possibleedgepairs(edgehits(ivspecial),1)]);
tp@0 1049 [uniquepatterns,ivec,jvec] = unique(patternlist,'rows');
tp@0 1050 keeppattern = zeros(size(ivec));
tp@0 1051 for mm = 1:length(ivec)
tp@0 1052 ivreftomatrix = uniquepatterns(mm,2) + (uniquepatterns(mm+1:end,2)-1)*nplanes;
tp@0 1053 coplanarindicator = coplanarsviaflatedge(ivreftomatrix);
tp@0 1054 ivcoplanars = find(uniquepatterns(mm+1:end,1)==uniquepatterns(mm,1) & uniquepatterns(mm+1:end,3)==uniquepatterns(mm,3) & coplanarindicator);
tp@0 1055 if ~isempty(ivcoplanars)
tp@0 1056 keeppattern(mm) = 1;
tp@0 1057 end
tp@0 1058 end
tp@0 1059 ivkeepers = find(keeppattern(jvec));
tp@0 1060 hitplanes = [hitplanes;edgehits(ivspecial(ivkeepers))];
tp@0 1061 reflpoints = [reflpoints;edgehitpoints(ivspecial(ivkeepers),:)];
tp@0 1062
tp@0 1063 end
tp@0 1064 end
tp@0 1065
tp@0 1066 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints;'])
tp@0 1067
tp@0 1068 expandedmasterivlist = expandedmasterivlist(hitplanes);
tp@0 1069 edgeweightlist = edgeweightlist(hitplanes,:);
tp@0 1070 possibleedgepairs = possibleedgepairs(hitplanes,:);
tp@0 1071 if ~isempty(possibleprespecs)
tp@0 1072 possibleprespecs = possibleprespecs(hitplanes,:);
tp@0 1073 end
tp@0 1074 possiblemidspecs = possiblemidspecs(hitplanes,:);
tp@0 1075 if ~isempty(possiblepostspecs)
tp@0 1076 possiblepostspecs = possiblepostspecs(hitplanes,:);
tp@0 1077 end
tp@0 1078 fromedgecoords = fromedgecoords(hitplanes,:);
tp@0 1079 toedgecoords = toedgecoords(hitplanes,:);
tp@0 1080 if kkpost > 0
tp@0 1081 expandedivcompletelyOK = expandedivcompletelyOK(hitplanes);
tp@0 1082 end
tp@0 1083
tp@0 1084 for mm = 1:jjmid
tp@0 1085 eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(hitplanes,:);']);
tp@0 1086 if mm > ll
tp@0 1087 eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)),' = reflpoints',JJ(mm,1:JJnumbofchars(mm)),'(hitplanes,:);']);
tp@0 1088 end
tp@0 1089 end
tp@0 1090 nposs = length(expandedmasterivlist);
tp@0 1091
tp@0 1092 if SHOWTEXT >= 3
tp@0 1093 if jjmid ~= 2
tp@0 1094 disp([' ',int2str(nposs),' edge+spec+edge combos survived the visibility test in reflection plane ',int2str(ll)])
tp@0 1095 else
tp@0 1096 disp([' ',int2str(nposs),' edge+spec+edge combos survived the visibility test in reflection plane ',int2str(ll),' + ',int2str(nfast),' cateye combos'])
tp@0 1097 end
tp@0 1098 end
tp@0 1099 end
tp@0 1100 end
tp@0 1101
tp@0 1102 %--------------------------------------------------------------
tp@0 1103 % Check for obstructions for all the paths, starting from edge number 2
tp@0 1104 % towards edge number 1.
tp@0 1105 %
tp@0 1106 % Reinsert the "fast lane" cases, i.e., the cateye reflections.
tp@0 1107
tp@0 1108 if jjmid == 2
tp@0 1109 if nfast > 0,
tp@0 1110 expandedmasterivlist = [expandedmasterivlist;expandedmasterivlistfast];
tp@0 1111 edgeweightlist = [edgeweightlist;edgeweightlistfast];
tp@0 1112 possibleedgepairs = [possibleedgepairs;possibleedgepairsfast];
tp@0 1113 if ~isempty(possibleprespecs)
tp@0 1114 possibleprespecs = [possibleprespecs;possibleprespecsfast];
tp@0 1115 end
tp@0 1116 possiblemidspecs = [possiblemidspecs;possiblemidspecsfast];
tp@0 1117 if ~isempty(possiblepostspecs)
tp@0 1118 possiblepostspecs = [possiblepostspecs;possiblepostspecsfast];
tp@0 1119 end
tp@0 1120 fromedgecoords = [fromedgecoords;fromedgecoordsfast];
tp@0 1121 toedgecoords = [toedgecoords;toedgecoordsfast];
tp@0 1122 for mm = 1:jjmid
tp@0 1123 eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = [bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),';bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'fast];']);
tp@0 1124 eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)), ' = [reflpoints',JJ(mm,1:JJnumbofchars(mm)), ';reflpoints',JJ(mm,1:JJnumbofchars(mm)),'fast];']);
tp@0 1125 end
tp@0 1126 nposs = length(expandedmasterivlist);
tp@0 1127 end
tp@0 1128 end
tp@0 1129
tp@0 1130 if nposs > 0 & obstructtestneeded
tp@0 1131
tp@0 1132 for ll = 1:jjmid+1
tp@0 1133 if nposs > 0
tp@0 1134 if ll == 1
tp@0 1135 fromcoords = fromedgecoords;
tp@0 1136 startplanes = [planesatedge(possibleedgepairs(:,2),1) planesatedge(possibleedgepairs(:,2),2)];
tp@0 1137 else
tp@0 1138 eval(['fromcoords = reflpoints',JJ(ll-1,1:JJnumbofchars(ll-1)),';'])
tp@0 1139 startplanes = possiblemidspecs(:,ll-1);
tp@0 1140 end
tp@0 1141 if ll == jjmid+1,
tp@0 1142 tocoords = toedgecoords;
tp@0 1143 endplanes = [planesatedge(possibleedgepairs(:,1),1) planesatedge(possibleedgepairs(:,1),2)];
tp@0 1144 else
tp@0 1145 eval(['tocoords = reflpoints',JJ(ll,1:JJnumbofchars(ll)),';'])
tp@0 1146 endplanes = possiblemidspecs(:,ll);
tp@0 1147 end
tp@0 1148 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 1149 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 1150
tp@0 1151 if ~isempty(edgehits)
tp@0 1152 nonobstructedpaths = setdiff(nonobstructedpaths,edgehits);
tp@0 1153 nobstructions = nposs - length(nonobstructedpaths);
tp@0 1154 end
tp@0 1155
tp@0 1156 if nobstructions > 0
tp@0 1157 expandedmasterivlist = expandedmasterivlist(nonobstructedpaths);
tp@0 1158
tp@0 1159 edgeweightlist = edgeweightlist(nonobstructedpaths,:);
tp@0 1160 possibleedgepairs = possibleedgepairs(nonobstructedpaths,:);
tp@0 1161 if ~isempty(possibleprespecs)
tp@0 1162 possibleprespecs = possibleprespecs(nonobstructedpaths,:);
tp@0 1163 end
tp@0 1164 possiblemidspecs = possiblemidspecs(nonobstructedpaths,:);
tp@0 1165 if ~isempty(possiblepostspecs)
tp@0 1166 possiblepostspecs = possiblepostspecs(nonobstructedpaths,:);
tp@0 1167 end
tp@0 1168 fromedgecoords = fromedgecoords(nonobstructedpaths,:);
tp@0 1169 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 1170 if kkpost > 0
tp@0 1171 expandedivcompletelyOK = expandedivcompletelyOK(nonobstructedpaths,:);
tp@0 1172 end
tp@0 1173
tp@0 1174 for mm = 1:jjmid
tp@0 1175 eval(['bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),' = bigedgeimagecoords',JJ(mm,1:JJnumbofchars(mm)),'(nonobstructedpaths,:);']);
tp@0 1176 if mm > ll
tp@0 1177 eval(['reflpoints',JJ(mm,1:JJnumbofchars(mm)),' = reflpoints',JJ(mm,1:JJnumbofchars(mm)),'(nonobstructedpaths,:);']);
tp@0 1178 end
tp@0 1179 end
tp@0 1180
tp@0 1181 end
tp@0 1182 nposs = length(expandedmasterivlist);
tp@0 1183
tp@0 1184 end
tp@0 1185 end
tp@0 1186
tp@0 1187 if SHOWTEXT >= 3
tp@0 1188 if jjmid == 2
tp@0 1189 disp([' ',int2str(nposs),' edge+spec+edge combos survived the obstruction test (including cateye cases)'])
tp@0 1190 else
tp@0 1191 disp([' ',int2str(nposs),' edge+spec+edge combos survived the obstruction test'])
tp@0 1192 end
tp@0 1193 end
tp@0 1194
tp@0 1195 end
tp@0 1196
tp@0 1197 if nposs > 0
tp@0 1198 edgedifflist = [edgedifflist;possibleedgepairs];
tp@0 1199 if iipre == 0
tp@0 1200 possibleprespecs = zeros(nposs,1);
tp@0 1201 end
tp@0 1202 if specorder <= 4
tp@0 1203 if specorder == 3
tp@0 1204 prespeclist = [prespeclist;[possibleprespecs]];
tp@0 1205 else
tp@0 1206 prespeclist = [prespeclist;[possibleprespecs zeros(nposs,1)]];
tp@0 1207 end
tp@0 1208 else
tp@0 1209 [n1,n2] = size(prespeclist);
tp@0 1210 [n3,n4] = size(possibleprespecs);
tp@0 1211 if n1 > 0
tp@0 1212 % Error found 20050202 PS
tp@0 1213 % Old wrong version prespeclist = [prespeclist;[possibleprespecs zeros(nposs,n4-n2)]];
tp@0 1214 prespeclist = [prespeclist;[possibleprespecs zeros(nposs,n2-n4)]];
tp@0 1215 else
tp@0 1216 % Error found 20050202 PS
tp@0 1217 % Old wrong version prespeclist = [possibleprespecs zeros(nposs,n4-n2)];
tp@0 1218 prespeclist = [possibleprespecs zeros(nposs,n2-n4)];
tp@0 1219 end
tp@0 1220 end
tp@0 1221 if jjmid == specorder-2
tp@0 1222 midspeclist = [midspeclist;possiblemidspecs];
tp@0 1223 else
tp@0 1224 midspeclist = [midspeclist;[possiblemidspecs zeros(nposs,specorder-2-jjmid) ]];
tp@0 1225 end
tp@0 1226 if kkpost == 0
tp@0 1227 possiblepostspecs = zeros(nposs,1);
tp@0 1228 end
tp@0 1229 if specorder <= 4
tp@0 1230 if specorder == 3
tp@0 1231 postspeclist = [postspeclist;[possiblepostspecs]];
tp@0 1232 else
tp@0 1233 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,1)]];
tp@0 1234 end
tp@0 1235 else
tp@0 1236 if kkpost == 0
tp@0 1237 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-3)]];
tp@0 1238 else
tp@0 1239 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-2-kkpost)]];
tp@0 1240 end
tp@0 1241 end
tp@0 1242 bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0 1243
tp@0 1244 % NB! It is correct below that the indices for the ISCOORDS should be
tp@0 1245 % ORIGINSFROM(ORIGINSFROM(masterivlist)), rather than masterivlist.
tp@0 1246 % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0 1247 % spec-spec-...-diff-diff combinations and then
tp@0 1248 % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0 1249 % ends with a diff has no image source.
tp@0 1250 % Also, two recursive references are needed since we need to get back
tp@0 1251 % through the two last diffractions to reach the last specular
tp@0 1252 % reflection.
tp@0 1253
tp@0 1254 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(expandedmasterivlist)));
tp@0 1255 for kk = 2:jjmid+kkpost
tp@0 1256 ivref = ORIGINSFROM(ivref);
tp@0 1257 end
tp@0 1258 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 1259 if kkpost == 0
tp@0 1260 validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0 1261 else
tp@0 1262 if jjmid ~= 2
tp@0 1263 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(expandedivcompletelyOK)),:)];
tp@0 1264 else
tp@0 1265 error(['ERROR: Calculation of IR coords for midspec = 2 and postspec not implemented yet!']);
tp@0 1266 end
tp@0 1267 end
tp@0 1268 end
tp@0 1269 end
tp@0 1270 end
tp@0 1271 end
tp@0 1272
tp@0 1273 % #######################################################################
tp@0 1274 % #
tp@0 1275 % # Pack the edge segments together because every little edge segment
tp@0 1276 % # is present as a separate edge
tp@0 1277 % # This can be done for all combinations at once.
tp@0 1278 % #
tp@0 1279 % #######################################################################
tp@0 1280
tp@0 1281 test = [prespeclist edgedifflist(:,1) midspeclist edgedifflist(:,2) postspeclist];
tp@0 1282
tp@0 1283 if ~isempty(test)
tp@0 1284
tp@0 1285 [ncombs,slask] = size(edgedifflist);
tp@0 1286 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 1287 ivremove = find(dtest==1);
tp@0 1288
tp@0 1289 while ~isempty(ivremove)
tp@0 1290 bigedgeweightlist(ivremove+1,:) = double(bigedgeweightlist(ivremove+1,:)) + double(bigedgeweightlist(ivremove,:));
tp@0 1291 bigedgeweightlist(ivremove,:) = [];
tp@0 1292 edgedifflist(ivremove,:) = [];
tp@0 1293 prespeclist(ivremove,:) = [];
tp@0 1294 midspeclist(ivremove,:) = [];
tp@0 1295 postspeclist(ivremove,:) = [];
tp@0 1296 validIScoords(ivremove,:) = [];
tp@0 1297 validIRcoords(ivremove,:) = [];
tp@0 1298
tp@0 1299 test = [prespeclist edgedifflist(:,1) midspeclist edgedifflist(:,2) postspeclist];
tp@0 1300 [ncombs,slask] = size(edgedifflist);
tp@0 1301 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 1302 ivremove = find(dtest==1);
tp@0 1303
tp@0 1304 end
tp@0 1305
tp@0 1306 end
tp@0 1307
tp@0 1308 % #######################################################################
tp@0 1309 % #
tp@0 1310 % # The weights of the visible edge segments should be
tp@0 1311 % # translated into start and end points, together with the visibility
tp@0 1312 % # weights from the receiver.
tp@0 1313 % # This can be done for all combinations at once!
tp@0 1314 % #
tp@0 1315 % #######################################################################
tp@0 1316
tp@0 1317 % As a start we set the start and end values to 0 and 1, i.e. assuming full
tp@0 1318 % visibility.
tp@0 1319
tp@0 1320 startandendpoints = [startandendpoints;...
tp@0 1321 [zeros(length(edgedifflist),1),ones(length(edgedifflist),1),...
tp@0 1322 zeros(length(edgedifflist),1),ones(length(edgedifflist),1)]];
tp@0 1323
tp@0 1324 % Treat edge 1
tp@0 1325
tp@0 1326 ivtoaccountfor = [1:size(edgedifflist,1)].';
tp@0 1327
tp@0 1328 ivwholeedge1 = find( bigedgeweightlist(:,1) == maxvisibilityvalue);
tp@0 1329 if ~isempty(ivwholeedge1)
tp@0 1330 ivtoaccountfor(ivwholeedge1) = [];
tp@0 1331 end
tp@0 1332
tp@0 1333 if ~isempty(ivtoaccountfor)
tp@0 1334 ncombs = length(ivtoaccountfor);
tp@0 1335 bitpattern = zeros(ncombs,nedgesubs);
tp@0 1336 for ii=1:nedgesubs
tp@0 1337 bitpattern(:,ii) = bitget(bigedgeweightlist(ivtoaccountfor,1),ii);
tp@0 1338 end
tp@0 1339 dbit1 = diff([zeros(ncombs,1) bitpattern].').';
tp@0 1340 dbit2 = [dbit1 -bitpattern(:,nedgesubs)];
tp@0 1341
tp@0 1342 nsegments = ceil((sum(abs(dbit1.')).')/2);
tp@0 1343 ivonesegments = find(nsegments==1);
tp@0 1344
tp@0 1345 if ~isempty(ivonesegments)
tp@0 1346
tp@0 1347 nonesegments = length(ivonesegments);
tp@0 1348 multvec = 2.^[0:nedgesubs];
tp@0 1349 segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 1350 segendpos = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 1351
tp@0 1352 ivmodify = find(segstartpos==1);
tp@0 1353 segstartpos(ivmodify) = ones(size(ivmodify))*1.5;
tp@0 1354 ivmodify = find(segendpos>nedgesubs);
tp@0 1355 segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5);
tp@0 1356
tp@0 1357 startandendpoints(ivtoaccountfor(ivonesegments),1) = (segstartpos-1.5)/(nedgesubs-1);
tp@0 1358 startandendpoints(ivtoaccountfor(ivonesegments),2) = (segendpos-1.5)/(nedgesubs-1);
tp@0 1359
tp@0 1360 end
tp@0 1361
tp@0 1362 % If we have some two-or-more-subsegments cases, they will be
tp@0 1363 % discovered by the if-condition below
tp@0 1364
tp@0 1365 if length(ivonesegments) < ncombs
tp@0 1366 for nsegmentstocheck = 2:ceil(nedgesubs/2)
tp@0 1367 disp(['Checking for ',int2str(nsegmentstocheck),' sub-segments'])
tp@0 1368 ivNsegments = find(nsegments==nsegmentstocheck);
tp@0 1369 if ~isempty(ivNsegments)
tp@0 1370 [n1,n2] = size(startandendpoints);
tp@0 1371 if n2 < 4*nsegmentstocheck
tp@0 1372 startandendpoints = [startandendpoints zeros(n1,4*nsegmentstocheck-n2)];
tp@0 1373 end
tp@0 1374 for jj = 1:length(ivNsegments)
tp@0 1375 ivstartbits = find(dbit2(ivNsegments(jj),:) == 1);
tp@0 1376 ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits;
tp@0 1377 ivendbits = find(dbit2(ivNsegments(jj),:) == -1);
tp@0 1378 ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits;
tp@0 1379
tp@0 1380 for kk = 1:nsegmentstocheck,
tp@0 1381 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+1) = (ivstartbits(kk)-1.5)/(nedgesubs-1);
tp@0 1382 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+2) = (ivendbits(kk)-1.5)/(nedgesubs-1);
tp@0 1383 end
tp@0 1384 end
tp@0 1385 end
tp@0 1386
tp@0 1387 end
tp@0 1388 end
tp@0 1389 end
tp@0 1390
tp@0 1391 % Treat edge 2
tp@0 1392
tp@0 1393 ivtoaccountfor = [1:size(edgedifflist,1)].';
tp@0 1394
tp@0 1395 ivwholeedge2 = find( bigedgeweightlist(:,2) == maxvisibilityvalue);
tp@0 1396 if ~isempty(ivwholeedge2)
tp@0 1397 ivtoaccountfor(ivwholeedge2) = [];
tp@0 1398 end
tp@0 1399
tp@0 1400 if ~isempty(ivtoaccountfor)
tp@0 1401 ncombs = length(ivtoaccountfor);
tp@0 1402 bitpattern = zeros(ncombs,nedgesubs);
tp@0 1403 for ii=1:nedgesubs
tp@0 1404 bitpattern(:,ii) = bitget(bigedgeweightlist(ivtoaccountfor,2),ii);
tp@0 1405 end
tp@0 1406 dbit1 = diff([zeros(ncombs,1) bitpattern].').';
tp@0 1407 dbit2 = [dbit1 -bitpattern(:,nedgesubs)];
tp@0 1408
tp@0 1409 nsegments = ceil((sum(abs(dbit1.')).')/2);
tp@0 1410 ivonesegments = find(nsegments==1);
tp@0 1411
tp@0 1412 if ~isempty(ivonesegments)
tp@0 1413
tp@0 1414 nonesegments = length(ivonesegments);
tp@0 1415 multvec = 2.^[0:nedgesubs];
tp@0 1416 segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 1417 segendpos = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 1418
tp@0 1419 ivmodify = find(segstartpos==1);
tp@0 1420 segstartpos(ivmodify) = ones(size(ivmodify))*1.5;
tp@0 1421 ivmodify = find(segendpos>nedgesubs);
tp@0 1422 segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5);
tp@0 1423
tp@0 1424 startandendpoints(ivtoaccountfor(ivonesegments),3) = (segstartpos-1.5)/(nedgesubs-1);
tp@0 1425 startandendpoints(ivtoaccountfor(ivonesegments),4) = (segendpos-1.5)/(nedgesubs-1);
tp@0 1426
tp@0 1427 end
tp@0 1428
tp@0 1429 % If we have some two-or-more-subsegments cases, they will be
tp@0 1430 % discovered by the if-condition below
tp@0 1431
tp@0 1432 if length(ivonesegments) < ncombs
tp@0 1433 for nsegmentstocheck = 2:ceil(nedgesubs/2)
tp@0 1434
tp@0 1435 ivNsegments = find(nsegments==nsegmentstocheck);
tp@0 1436 if ~isempty(ivNsegments)
tp@0 1437 [n1,n2] = size(startandendpoints);
tp@0 1438 if n2 < 4*nsegmentstocheck
tp@0 1439 startandendpoints = [startandendpoints zeros(n1,4*nsegmentstocheck-n2)];
tp@0 1440 end
tp@0 1441 for jj = 1:length(ivNsegments)
tp@0 1442 ivstartbits = find(dbit2(ivNsegments(jj),:) == 1);
tp@0 1443 ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits;
tp@0 1444 ivendbits = find(dbit2(ivNsegments(jj),:) == -1);
tp@0 1445 ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits;
tp@0 1446
tp@0 1447 for kk = 1:nsegmentstocheck,
tp@0 1448 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+3) = (ivstartbits(kk)-1.5)/(nedgesubs-1);
tp@0 1449 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*4+4) = (ivendbits(kk)-1.5)/(nedgesubs-1);
tp@0 1450 end
tp@0 1451 end
tp@0 1452 end
tp@0 1453
tp@0 1454 end
tp@0 1455 end
tp@0 1456 end
tp@0 1457
tp@0 1458 % #######################################################################
tp@0 1459 % #
tp@0 1460 % # Construct a list guide, which will tell which rows have only
tp@0 1461 % # dd, which rows have sdd etc
tp@0 1462 % # Syntax: dd,sdd,ssdd,sssdd,...,dds,ddss,ddsss,...
tp@0 1463 % # (should also continue with sdds, sddss,...
tp@0 1464 % #
tp@0 1465 % #######################################################################
tp@0 1466
tp@0 1467 [n1,n2] = size(prespeclist);
tp@0 1468 if n2 > 1
tp@0 1469 nprespecs = sum(prespeclist.' > 0).';
tp@0 1470 else
tp@0 1471 nprespecs = (prespeclist>0);
tp@0 1472 end
tp@0 1473 [n1,n2] = size(midspeclist);
tp@0 1474 if n2 > 1
tp@0 1475 nmidspecs = sum(midspeclist.' > 0).';
tp@0 1476 else
tp@0 1477 nmidspecs = (midspeclist>0);
tp@0 1478 end
tp@0 1479 [n1,n2] = size(postspeclist);
tp@0 1480 if n2 > 1
tp@0 1481 npostspecs = sum(postspeclist.' > 0).';
tp@0 1482 else
tp@0 1483 npostspecs = (postspeclist>0);
tp@0 1484 end
tp@0 1485
tp@0 1486 [B,ivec,jvec] = unique([nprespecs nmidspecs npostspecs],'rows');
tp@0 1487 nuniquecombs = length(ivec);
tp@0 1488 ntotcombs = length(jvec);
tp@0 1489
tp@0 1490 listguide = zeros(nuniquecombs,3);
tp@0 1491 listofallspecs = zeros(nuniquecombs,3);
tp@0 1492 sortvec = zeros(ntotcombs,1);
tp@0 1493 for ii = 1:length(ivec)
tp@0 1494 ivfindcombs = find(jvec==ii);
tp@0 1495 listguide(ii,1) = length(ivfindcombs);
tp@0 1496 if ii > 1
tp@0 1497 listguide(ii,2) = listguide(ii-1,3)+1;
tp@0 1498 else
tp@0 1499 listguide(ii,2) = 1;
tp@0 1500 end
tp@0 1501 listguide(ii,3) = listguide(ii,2)+listguide(ii,1)-1;
tp@0 1502 listofallspecs(ii,:) = [B(ii,:)];
tp@0 1503
tp@0 1504 sortvec(listguide(ii,2):listguide(ii,3)) = ivfindcombs;
tp@0 1505
tp@0 1506 end
tp@0 1507
tp@0 1508 prespeclist = prespeclist(sortvec,:);
tp@0 1509 midspeclist = midspeclist(sortvec,:);
tp@0 1510 postspeclist = postspeclist(sortvec,:);
tp@0 1511 validIScoords = validIScoords(sortvec,:);
tp@0 1512 validIRcoords = validIRcoords(sortvec,:);
tp@0 1513 edgedifflist = edgedifflist(sortvec,:);
tp@0 1514 startandendpoints = startandendpoints(sortvec,:);