Mercurial > hg > human-echolocation

annotate private/EDB1diffNISES.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,postspeclist,validIScoords,validIRcoords,listguide,...
tp@0 2 listoforders] = EDB1diffNISES(eddatafile,S,R,...
tp@0 3 ivNdiffmatrix,lengthNdiffmatrix,Ndifforder,specorder,visplanesfromR,vispartedgesfromS,...
tp@0 4 vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlistin,validEDIRcoords)
tp@0 5 % EDB1diffNISES - Gives a list of paths that includes an N-order diffraction path.
tp@0 6 % Gives the list of visible N-diffraction paths, possibly with specular reflections before and after.
tp@0 7 %
tp@0 8 % Input parameters:
tp@0 9 % eddatafile,S,R,ivsinglediff,...
tp@0 10 % singlediffcol,startindicessinglediff,endindicessinglediff,specorder,visplanesfromR,...
tp@0 11 % vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,PointertoIRcombs,IRoriginsfrom;
tp@0 12 %
tp@0 13 % GLOBAL parameters:
tp@0 14 % SHOWTEXT JJ JJnumbofchars See EDB1mainISES
tp@0 15 % POTENTIALISES,ISCOORDS,ORIGINSFROM,ISESVISIBILITY,REFLORDER See EDB1findISEStree
tp@0 16 %
tp@0 17 % Output parameters:
tp@0 18 % edgedifflist List [ncombs,N] of the edge numbers involved in each
tp@0 19 % spec-diff-...-diff-spec combination.
tp@0 20 % startandendpoints Matrix [ncombs,4] of the relative start and end
tp@0 21 % points of each edge. The values, [0,1], indicate
tp@0 22 % which part of the two edges that are visible.
tp@0 23 % prespeclist Matrix [ncombs,specorder-N] of the specular
tp@0 24 % reflections that precede every diffraction.
tp@0 25 % postspeclist Matrix [ncombs,specorder-N] of the specular
tp@0 26 % reflections that follow every diffraction.
tp@0 27 % validIScoords Matrix [ncombs,3] of the image source for each
tp@0 28 % multiple-spec that precedes the diffraction. If
tp@0 29 % there is no spec refl before the diffraction, the
tp@0 30 % value [0 0 0] is given.
tp@0 31 % validIRcoords Matrix [ncombs,3] of the image receiver for each
tp@0 32 % multiple-spec that follows the diffraction. If
tp@0 33 % there is no spec refl after the diffraction, the
tp@0 34 % value [0 0 0] is given.
tp@0 35 % listguide Matrix [nuniquecombs,3] which for each row gives
tp@0 36 % 1. The number of examples in edgefdifflist etc that
tp@0 37 % are the same type of spec-diff-...-diff-spec comb.
tp@0 38 % 2. The first row number and 3. The last row number.
tp@0 39 % listoforders Matrix [nuniquecombs,2] which for each row gives
tp@0 40 % 1. The reflection order for the spec-diff-spec comb
tp@0 41 % in the same row in listguide.
tp@0 42 % 2. The order of the diffraction in the
tp@0 43 % spec-diff-...-diff-spec comb.
tp@0 44 %
tp@0 45 % Uses functions EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths
tp@0 46 %
tp@0 47 % ----------------------------------------------------------------------------------------------
tp@0 48 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 49 %
tp@0 50 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 51 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 52 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 53 %
tp@0 54 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 55 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 56 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 57 %
tp@0 58 % You should have received a copy of the GNU General Public License along with the
tp@0 59 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 60 % ----------------------------------------------------------------------------------------------
tp@0 61 % Peter Svensson (svensson@iet.ntnu.no) 20031018
tp@0 62 %
tp@0 63 % [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,...
tp@0 64 % listguide] = EDB1diff3ISES(eddatafile,S,R,...
tp@0 65 % ivNdiffmatrix,lengthNdiffmatrix,Ndifforder,specorder,visplanesfromR,vispartedgesfromS,...
tp@0 66 % vispartedgesfromR,nedgesubs,ndecimaldivider,edgedifflistin,postspeclistin,bigedgeweightlist,validEDIRcoords);
tp@0 67
tp@0 68 global SHOWTEXT JJ JJnumbofchars
tp@0 69 global POTENTIALISES ISCOORDS ORIGINSFROM ISESVISIBILITY REFLORDER
tp@0 70
tp@0 71 eval(['load ',eddatafile])
tp@0 72
tp@0 73 [nedges,slask] = size(planesatedge);
tp@0 74 [nplanes,slask] = size(planecorners);
tp@0 75 multfac = 10^(ceil(log10(max([nedges nplanes]))));
tp@0 76
tp@0 77 edgedifflist = [];
tp@0 78 prespeclist = [];
tp@0 79 postspeclist = [];
tp@0 80 startandendpoints = [];
tp@0 81 bigedgeweightlist = [];
tp@0 82 validIScoords = [];
tp@0 83 validIRcoords = [];
tp@0 84
tp@0 85 [n1,n2] = size(POTENTIALISES);
tp@0 86 if n2 < specorder
tp@0 87 specorder = n2;
tp@0 88 end
tp@0 89
tp@0 90 maxvisibilityvalue = 2^nedgesubs-1;
tp@0 91 zerosvec1 = zeros(1,specorder-Ndifforder);
tp@0 92 zerosvec2 = zeros(1,3);
tp@0 93 listguide = zeros(specorder*2-1,3);
tp@0 94
tp@0 95 obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0 96 onesvec = ones(1,nedgesubs);
tp@0 97 onesvec3 = ones(1,3);
tp@0 98
tp@0 99 npostspecs = sum(double(postspeclistin.'>0)).';
tp@0 100
tp@0 101 % ###########################################
tp@0 102 % # #
tp@0 103 % # S - edge - edge - edge -R cases #
tp@0 104 % # #
tp@0 105 % ###########################################
tp@0 106 %
tp@0 107 % Possible edges for S-E-E-E-R are seen (at least partly) by the source and by the receiver.
tp@0 108 %
tp@0 109 % The visibility by the source is taken care of by the findISEStree
tp@0 110 % so we must check which ones are visible by the receiver.
tp@0 111 % Also, the active edge segment(s) must be selected but this is done
tp@0 112 % further down together with the S-spec-spec-edge-R cases
tp@0 113
tp@0 114 ivNdiff = ivNdiffmatrix(1:lengthNdiffmatrix(Ndifforder),Ndifforder);
tp@0 115 ivsinglediff = ivNdiffmatrix(1:lengthNdiffmatrix(1),1);
tp@0 116
tp@0 117 ivSEEER = ivNdiff(find(REFLORDER(ivNdiff)==Ndifforder));
tp@0 118
tp@0 119 possibleedgecombs = double(POTENTIALISES(ivSEEER,1:Ndifforder))-nplanes;
tp@0 120 ivnotvisiblefromr = find(vispartedgesfromR(possibleedgecombs(:,Ndifforder))==0);
tp@0 121 if ~isempty(ivnotvisiblefromr)
tp@0 122 possibleedgecombs(ivnotvisiblefromr,:) = [];
tp@0 123 end
tp@0 124
tp@0 125 edgedifflist = [edgedifflist;possibleedgecombs];
tp@0 126 bigedgeweightlist = [bigedgeweightlist;[vispartedgesfromS(edgedifflist(:,1)) vispartedgesfromR(edgedifflist(:,Ndifforder))]];
tp@0 127
tp@0 128 [nedgesadded,slask] = size(edgedifflist);
tp@0 129 zerosvec3 = zeros(nedgesadded,1);
tp@0 130 ndiffonly = nedgesadded;
tp@0 131
tp@0 132 prespeclist = [prespeclist;zerosvec3(:,ones(1,max(specorder-Ndifforder,1)))];
tp@0 133 postspeclist = [postspeclist;zerosvec3(:,ones(1,max(specorder-Ndifforder,1)))];
tp@0 134 validIScoords = [validIScoords;S(ones(nedgesadded,1),:)];
tp@0 135 validIRcoords = [validIRcoords;R(ones(nedgesadded,1),:)];
tp@0 136
tp@0 137 if SHOWTEXT >= 3
tp@0 138 disp([' ',int2str(nedgesadded),' Nth diff valid'])
tp@0 139 end
tp@0 140
tp@0 141 % ###########################################
tp@0 142 % ###########################################
tp@0 143 % ###########################################
tp@0 144 % ###########################################
tp@0 145 % # #
tp@0 146 % # Prespec cases #
tp@0 147 % # #
tp@0 148 % ###########################################
tp@0 149 %
tp@0 150 % Possible edges for S-spec-spec-E-R are seen (at least partly) by the receiver.
tp@0 151 %
tp@0 152 % The visibility doesn't need to be checked since the source-to-edge paths
tp@0 153 % were checked in the ISEStree, and the visibility from the receiver also
tp@0 154 % has been checked.
tp@0 155
tp@0 156 % The vector ivmultidiff will always refer to the original data vector
tp@0 157 % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc
tp@0 158 %
tp@0 159 % We should remove the combinations which involve an edge that the
tp@0 160 % receiver can not see, but since the edge number is in different columns
tp@0 161 % we do that in the for loop.
tp@0 162
tp@0 163 % The ii-loop will go through: spec-diff, spec-spec-diff
tp@0 164 % spec-spec-spec-diff etc
tp@0 165
tp@0 166 for ii = 1:specorder-Ndifforder
tp@0 167
tp@0 168 if SHOWTEXT >= 3
tp@0 169 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before the Nth edge diff'])
tp@0 170 end
tp@0 171
tp@0 172 % Select the combinations where the reflection order == ii+Ndifforder
tp@0 173 % (which means ii specular reflections in addition to the diffraction)
tp@0 174 % and where the last Ndifforder columns contain edges.
tp@0 175
tp@0 176 iv = find((REFLORDER(ivNdiff)==ii+Ndifforder) & prod( double( POTENTIALISES(ivNdiff,ii+1:ii+Ndifforder)>nplanes).' ).' );
tp@0 177 masterivlist = ivNdiff(iv);
tp@0 178 possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0 179
tp@0 180 % Keep only combinations for which the receiver can see the edge
tp@0 181
tp@0 182 ivnotvisiblefromr = find(vispartedgesfromR(possibleedgecombs(:,Ndifforder))==0);
tp@0 183 if ~isempty(ivnotvisiblefromr)
tp@0 184 masterivlist(ivnotvisiblefromr) = [];
tp@0 185 possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0 186 end
tp@0 187
tp@0 188 possibleprespecs = POTENTIALISES(masterivlist,1:ii);
tp@0 189
tp@0 190 reftoIScoords = ORIGINSFROM(masterivlist);
tp@0 191 edgeweightlist = [ISESVISIBILITY(masterivlist) vispartedgesfromR(possibleedgecombs(:,Ndifforder))];
tp@0 192
tp@0 193 nposs = length(masterivlist);
tp@0 194
tp@0 195 if SHOWTEXT >= 3
tp@0 196 disp([' ',int2str(nposs),' IS+edge multiples valid'])
tp@0 197 end
tp@0 198
tp@0 199 edgedifflist = [edgedifflist;possibleedgecombs];
tp@0 200 prespeclist = [prespeclist;[possibleprespecs zeros(nposs,specorder-Ndifforder-ii)]];
tp@0 201 postspeclist = [postspeclist;zeros(nposs,specorder-Ndifforder)];
tp@0 202 bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0 203 % The ref. to ISCOORDS must be nested the correct number of time to
tp@0 204 % come back to the image source coords of the last spec. refl.
tp@0 205 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist)));
tp@0 206 for jj = 4:Ndifforder
tp@0 207 ivref = ORIGINSFROM(ivref);
tp@0 208 end
tp@0 209 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 210 validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0 211
tp@0 212 end
tp@0 213
tp@0 214 % #######################################################################
tp@0 215 % #######################################################################
tp@0 216 % #######################################################################
tp@0 217 % #######################################################################
tp@0 218 % #######################################################################
tp@0 219
tp@0 220 % #############################################
tp@0 221 % # #
tp@0 222 % # S - edge - edge -spec - R cases #
tp@0 223 % # #
tp@0 224 % # Postspec cases #
tp@0 225 % # #
tp@0 226 % #############################################
tp@0 227 %
tp@0 228 % Possible edges for S-E-E-spec-spec-R are seen (at least partly) by the receiver.
tp@0 229 %
tp@0 230 % For some of the combos, we have the IR coords and the edge visibility
tp@0 231 % from the single-diffraction run. For potential combos that were not
tp@0 232 % included there, they are either invisible/obstructed or were not tested.
tp@0 233 % We can figure out which ones were tested because they can be found in the
tp@0 234 % POTENTIALISES under edge-spec-spec but not in the final valid list.
tp@0 235
tp@0 236 % The vector masterivlist will always refer to the original data vector
tp@0 237 % i.e. PotentialIR, ORIGINSFROMIR, reflorderIR etc
tp@0 238 %
tp@0 239 % First we pick out those indices where there was a single diffraction, but
tp@0 240 % skip those with only diffraction (because we dealt with them already).
tp@0 241 % Also, select only those where the diffraction is the first in the sequence
tp@0 242 % of reflections.
tp@0 243
tp@0 244 for ii = 1:specorder-Ndifforder
tp@0 245
tp@0 246 if SHOWTEXT >= 3
tp@0 247 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the multiple edge diff'])
tp@0 248 end
tp@0 249
tp@0 250 % Select the combinations where the reflection order == ii+Ndifforder
tp@0 251 % (which means ii specular reflections in addition to the diffraction)
tp@0 252 % and where the first Ndifforder columns contain edges.
tp@0 253
tp@0 254 iv = find((REFLORDER(ivNdiff)==ii+Ndifforder) & prod( double( POTENTIALISES(ivNdiff,1:Ndifforder)>nplanes).' ).' );
tp@0 255 masterivlist = ivNdiff(iv);
tp@0 256 possibleedgecombs = double(POTENTIALISES(masterivlist,1:Ndifforder)) - nplanes;
tp@0 257 possiblepostspecs = POTENTIALISES(masterivlist,Ndifforder+1:specorder);
tp@0 258 possibleweights = ISESVISIBILITY(masterivlist);
tp@0 259
tp@0 260 % Compare with those that have already been found OK
tp@0 261 ivOK = find(npostspecs==ii);
tp@0 262 if ~isempty(ivOK)
tp@0 263 patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:ii)];
tp@0 264 else
tp@0 265 patternOK = [];
tp@0 266 end
tp@0 267 % Find out which ones have been checked and found invisible/obstructed
tp@0 268 ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == ii+1));
tp@0 269 patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+ii))];
tp@0 270 if ~isempty(patternOK)
tp@0 271 patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0 272 else
tp@0 273 patternNOTOK = patternALL;
tp@0 274 end
tp@0 275 patterntocompare = [possibleedgecombs(:,Ndifforder) possiblepostspecs(:,1:ii)];
tp@0 276
tp@0 277 ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0 278 masterivlist(ivtocancel) = [];
tp@0 279 possibleedgecombs(ivtocancel,:) = [];
tp@0 280 possiblepostspecs(ivtocancel,:) = [];
tp@0 281 possibleweights(ivtocancel,:) = [];
tp@0 282 patterntocompare(ivtocancel,:) = [];
tp@0 283
tp@0 284 [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0 285 ivmustbechecked = find(ivcompletelyOK==0);
tp@0 286 ivcompletelyOK = find(ivcompletelyOK);
tp@0 287 if ~isempty(ivmustbechecked)
tp@0 288 masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0 289 ntocheckmore = length(masterlisttocheckmore);
tp@0 290
tp@0 291 %----------------------------------------------
tp@0 292 % Must carry out a visibility and obstruction check for the special
tp@0 293 % combinations here.
tp@0 294 % These combinations have a postspec-combination that hasn't been
tp@0 295 % encountered in the single diffraction cases, so no visibility
tp@0 296 % test has been made for these.
tp@0 297
tp@0 298 lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,Ndifforder))-nplanes;
tp@0 299 newIRcoords = R;
tp@0 300 reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0 301 for jj = 1:ii
tp@0 302 reflplanes = POTENTIALISES(masterlisttocheckmore,Ndifforder+1+ii-jj);
tp@0 303 reflplanesexpand(:,jj) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0 304 newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,size(newIRcoords,1),onesvec3);
tp@0 305 newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0 306 eval(['newIRcoords',JJ(jj,1:JJnumbofchars(jj)),' = newIRcoordsexpand;'])
tp@0 307 end
tp@0 308 [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0 309 tocoords = toedgecoords;
tp@0 310 lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0 311 lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0 312 masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0 313 masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0 314
tp@0 315 ntocheckmore = length(masterlisttocheckmore);
tp@0 316 if SHOWTEXT >= 3
tp@0 317 disp([' ',int2str(ntocheckmore),' special N-edge - IR combinations to check'])
tp@0 318 end
tp@0 319
tp@0 320 for jj = ii:-1:1
tp@0 321 if length(masterlisttocheckmore) > 0
tp@0 322 eval(['fromcoords = newIRcoords',JJ(jj,1:JJnumbofchars(jj)),';']);
tp@0 323 if jj < ii
tp@0 324 eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])
tp@0 325 end
tp@0 326
tp@0 327 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,jj),4),planenvecs(reflplanesexpand(:,jj),:),minvals(reflplanesexpand(:,jj),:),...
tp@0 328 maxvals(reflplanesexpand(:,jj),:),planecorners(reflplanesexpand(:,jj),:),corners,ncornersperplanevec(reflplanesexpand(:,jj)));
tp@0 329 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 330 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 331 disp(' implemented.')
tp@0 332 end
tp@0 333 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0 334
tp@0 335 masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0 336 edgeweightlist = edgeweightlist(hitplanes);
tp@0 337 lastedgenumbers = lastedgenumbers(hitplanes);
tp@0 338 reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0 339 toedgecoords = toedgecoords(hitplanes,:);
tp@0 340
tp@0 341 for kk = 1:ii
tp@0 342 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0 343 if kk > jj
tp@0 344 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);']);
tp@0 345 end
tp@0 346 end
tp@0 347 ntocheckmore = length(masterlisttocheckmore);
tp@0 348
tp@0 349 end
tp@0 350
tp@0 351 if SHOWTEXT >= 3
tp@0 352 disp([' ',int2str(ntocheckmore),' of the special N-edge - IR combinations survived the visibility test in refl plane ',int2str(jj)])
tp@0 353 end
tp@0 354 end
tp@0 355
tp@0 356 if obstructtestneeded & ntocheckmore > 0
tp@0 357
tp@0 358 for jj = 1:ii+1
tp@0 359 if ntocheckmore > 0
tp@0 360 if jj == 1
tp@0 361 fromcoords = R;
tp@0 362 startplanes = [];
tp@0 363 else
tp@0 364 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0 365 startplanes = reflplanesexpand(:,jj-1);
tp@0 366 end
tp@0 367 if jj == ii+1,
tp@0 368 tocoords = toedgecoords;
tp@0 369 endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0 370 else
tp@0 371 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
tp@0 372 endplanes = reflplanesexpand(:,jj);
tp@0 373 end
tp@0 374
tp@0 375 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 376 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 377 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 378 disp('WARNING! An edgehit or cornerhit occurred during the obstruction test but this is not')
tp@0 379 disp(' implemented.')
tp@0 380 end
tp@0 381
tp@0 382 if nobstructions > 0
tp@0 383 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0 384 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 385 lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0 386 reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0 387 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 388 for kk = 1:ii
tp@0 389 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0 390 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);']);
tp@0 391 end
tp@0 392
tp@0 393 end
tp@0 394 ntocheckmore = length(masterlisttocheckmore);
tp@0 395
tp@0 396 end
tp@0 397
tp@0 398 end
tp@0 399 if SHOWTEXT >= 3
tp@0 400 disp([' ',int2str(ntocheckmore),' of the special N-edge - IR combinations survived the obstruction test'])
tp@0 401 end
tp@0 402
tp@0 403 end
tp@0 404
tp@0 405 % Add the found special combinations to the outdata list
tp@0 406
tp@0 407 edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,1:Ndifforder))-nplanes];
tp@0 408 prespeclist = [prespeclist;zerosvec1(ones(ntocheckmore,1),:)];
tp@0 409 postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0 410 bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];
tp@0 411 eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(ii,1:JJnumbofchars(ii)),'];']);
tp@0 412 validIScoords = [validIScoords;S(ones(ntocheckmore,1),:)];
tp@0 413
tp@0 414 end
tp@0 415
tp@0 416 masterivlist = masterivlist(ivcompletelyOK);
tp@0 417 possibleedgecombs = possibleedgecombs(ivcompletelyOK,:);
tp@0 418 possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0 419 possibleweights = possibleweights(ivcompletelyOK,:);
tp@0 420
tp@0 421 nposs = length(ivcompletelyOK);
tp@0 422
tp@0 423 if SHOWTEXT >= 3
tp@0 424 disp([' ',int2str(nposs),' Edge,Nth-order + IR segments (non-special) survived the obstruction test'])
tp@0 425 end
tp@0 426
tp@0 427 % Add the found "standard" combinations to the outdata list
tp@0 428
tp@0 429 edgedifflist = [edgedifflist;possibleedgecombs];
tp@0 430 prespeclist = [prespeclist;zeros(nposs,specorder-Ndifforder)];
tp@0 431 postspeclist = [postspeclist;possiblepostspecs ];
tp@0 432 bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];
tp@0 433 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0 434 validIScoords = [validIScoords;S(ones(nposs,1),:)];
tp@0 435
tp@0 436 end
tp@0 437
tp@0 438
tp@0 439 % #######################################################################
tp@0 440 % #######################################################################
tp@0 441 % #######################################################################
tp@0 442 % #######################################################################
tp@0 443 % #######################################################################
tp@0 444
tp@0 445 % #############################################
tp@0 446 % # #
tp@0 447 % # S - spec - edge - spec - R cases #
tp@0 448 % # #s
tp@0 449 % # Pre and postspec cases #
tp@0 450 % # #
tp@0 451 % #############################################
tp@0 452
tp@0 453 % The ii- and jj-loops will go through all combinations of spec before and after
tp@0 454 % the Ndifforder diffractions.
tp@0 455 %
tp@0 456 % Unlike before, we will look in the list of already found combinations
tp@0 457 % (edgedifflist etc)
tp@0 458
tp@0 459
tp@0 460 for ii = 1:specorder-Ndifforder-1
tp@0 461
tp@0 462 for jj = 1:specorder-ii-Ndifforder
tp@0 463
tp@0 464 if SHOWTEXT >= 3
tp@0 465 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl before and ',JJ(jj,1:JJnumbofchars(jj)),' spec refl after the N-order edge diff'])
tp@0 466 end
tp@0 467
tp@0 468 % Select the combinations where the reflection order ==
tp@0 469 % ii+jj+Ndifforder
tp@0 470 % (which means ii+jj specular reflections in addition to the
tp@0 471 % diffraction), and where columns ii+1:ii+Ndifforder of POTENTIALISES
tp@0 472 % contain edges.
tp@0 473
tp@0 474 iv = find((REFLORDER(ivNdiff) == ii+jj+Ndifforder) & prod( double( POTENTIALISES(ivNdiff,ii+1:ii+Ndifforder)>nplanes).' ).' );
tp@0 475 masterivlist = ivNdiff(iv);
tp@0 476 possibleedgecombs = double(POTENTIALISES(masterivlist,ii+1:ii+Ndifforder)) - nplanes;
tp@0 477 possibleprespecs = POTENTIALISES(masterivlist,1:ii);
tp@0 478 possiblepostspecs = POTENTIALISES(masterivlist,ii+Ndifforder+1:ii+Ndifforder+jj);
tp@0 479 possibleweights = ISESVISIBILITY(masterivlist);
tp@0 480
tp@0 481 % Compare with those that have already been found OK
tp@0 482 ivOK = find(npostspecs==jj);
tp@0 483 if ~isempty(ivOK)
tp@0 484 patternOK = [edgedifflistin(ivOK) postspeclistin(ivOK,1:jj)];
tp@0 485 else
tp@0 486 patternOK = [];
tp@0 487 end
tp@0 488
tp@0 489 % Find out which ones have been checked and found invisible/obstructed
tp@0 490 ivallcombs = ivsinglediff(find( POTENTIALISES(ivsinglediff,1)>nplanes & REFLORDER(ivsinglediff) == jj+1));
tp@0 491 patternALL = [double(POTENTIALISES(ivallcombs,1))-nplanes double(POTENTIALISES(ivallcombs,2:1+jj))];
tp@0 492 if ( (~isempty(patternOK)) & (~isempty(patternALL)) )
tp@0 493 patternNOTOK = setdiff(patternALL,patternOK,'rows');
tp@0 494 else
tp@0 495 if ( (~isempty(patternOK)) )
tp@0 496 patternNOTOK = patternOK;
tp@0 497 else
tp@0 498 patternNOTOK = patternALL;
tp@0 499 end
tp@0 500 end
tp@0 501 patterntocompare = [possibleedgecombs(:,Ndifforder) possiblepostspecs(:,1:jj)];
tp@0 502
tp@0 503 ivtocancel = find(ismember(patterntocompare,patternNOTOK,'rows'));
tp@0 504 masterivlist(ivtocancel) = [];
tp@0 505 possibleedgecombs(ivtocancel,:) = [];
tp@0 506 possibleprespecs(ivtocancel,:) = [];
tp@0 507 possiblepostspecs(ivtocancel,:) = [];
tp@0 508 possibleweights(ivtocancel,:) = [];
tp@0 509 patterntocompare(ivtocancel,:) = [];
tp@0 510
tp@0 511 [ivcompletelyOK,ivreftoindata] = ismember(patterntocompare,patternOK,'rows');
tp@0 512 ivmustbechecked = find(ivcompletelyOK==0);
tp@0 513 ivcompletelyOK = find(ivcompletelyOK);
tp@0 514 if ~isempty(ivmustbechecked)
tp@0 515 masterlisttocheckmore = masterivlist(ivmustbechecked);
tp@0 516 ntocheckmore = length(masterlisttocheckmore);
tp@0 517
tp@0 518 %----------------------------------------------
tp@0 519 % Must carry out a visibility and obstruction check for the special
tp@0 520 % combinations here.
tp@0 521
tp@0 522 lastedgenumbers = double(POTENTIALISES(masterlisttocheckmore,ii+Ndifforder))-nplanes;
tp@0 523 newIRcoords = R;
tp@0 524 reflplanesexpand = zeros(ntocheckmore*nedgesubs,ii);
tp@0 525 for kk = 1:jj
tp@0 526 reflplanes = POTENTIALISES(masterlisttocheckmore,Ndifforder+1+ii+jj-kk);
tp@0 527 reflplanesexpand(:,kk) = reshape(reflplanes(:,onesvec).',ntocheckmore*nedgesubs,1);
tp@0 528 newIRcoords = EDB1findis(newIRcoords,reflplanes,planeeqs,size(newIRcoords,1),onesvec3);
tp@0 529 newIRcoordsexpand = reshape(repmat(newIRcoords.',nedgesubs,1),3,ntocheckmore*nedgesubs).';
tp@0 530 eval(['newIRcoords',JJ(kk,1:JJnumbofchars(kk)),' = newIRcoordsexpand;'])
tp@0 531 end
tp@0 532 [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(lastedgenumbers,:),edgeendcoords(lastedgenumbers,:),edgelengthvec(lastedgenumbers,:),nedgesubs);
tp@0 533 tocoords = toedgecoords;
tp@0 534 lastedgenumbers = lastedgenumbers(:,onesvec);
tp@0 535 lastedgenumbers = reshape(lastedgenumbers.',ntocheckmore*nedgesubs,1);
tp@0 536 masterlisttocheckmore = masterlisttocheckmore(:,onesvec);
tp@0 537 masterlisttocheckmore = reshape(masterlisttocheckmore.',ntocheckmore*nedgesubs,1);
tp@0 538
tp@0 539 ntocheckmore = length(masterlisttocheckmore);
tp@0 540 if SHOWTEXT >= 3
tp@0 541 disp([' ',int2str(ntocheckmore),' special IS - N-edge - IR combinations to check'])
tp@0 542 end
tp@0 543
tp@0 544 for kk = jj:-1:1
tp@0 545 if length(masterlisttocheckmore) > 0
tp@0 546 eval(['fromcoords = newIRcoords',JJ(kk,1:JJnumbofchars(kk)),';']);
tp@0 547 if kk < jj
tp@0 548 eval(['tocoords = reflpoints',JJ(kk+1,1:JJnumbofchars(kk+1)),';'])
tp@0 549 end
tp@0 550
tp@0 551 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(reflplanesexpand(:,kk),4),planenvecs(reflplanesexpand(:,kk),:),minvals(reflplanesexpand(:,kk),:),...
tp@0 552 maxvals(reflplanesexpand(:,kk),:),planecorners(reflplanesexpand(:,kk),:),corners,ncornersperplanevec(reflplanesexpand(:,kk)));
tp@0 553 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 554 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 555 disp(' implemented.')
tp@0 556 end
tp@0 557 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints;'])
tp@0 558
tp@0 559 masterlisttocheckmore = masterlisttocheckmore(hitplanes);
tp@0 560 edgeweightlist = edgeweightlist(hitplanes);
tp@0 561 lastedgenumbers = lastedgenumbers(hitplanes);
tp@0 562 reflplanesexpand = reflplanesexpand(hitplanes,:);
tp@0 563 toedgecoords = toedgecoords(hitplanes,:);
tp@0 564
tp@0 565 for ll = 1:jj
tp@0 566 eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0 567 if ll > kk
tp@0 568 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(hitplanes,:);']);
tp@0 569 end
tp@0 570 end
tp@0 571 ntocheckmore = length(masterlisttocheckmore);
tp@0 572
tp@0 573 end
tp@0 574
tp@0 575 if SHOWTEXT >= 3
tp@0 576 disp([' ',int2str(ntocheckmore),' of the special IS - N-edge - IR combinations survived the visibility test in refl plane ',int2str(kk)])
tp@0 577 end
tp@0 578
tp@0 579 end
tp@0 580
tp@0 581 % Obstruction test of all the involved paths: R ->
tp@0 582 % reflplane1 -> reflplane2 -> ... -> last edge
tp@0 583
tp@0 584 if obstructtestneeded & ntocheckmore > 0
tp@0 585
tp@0 586 for kk = 1:jj+1
tp@0 587 if ntocheckmore > 0
tp@0 588 if kk == 1
tp@0 589 fromcoords = R;
tp@0 590 startplanes = [];
tp@0 591 else
tp@0 592 eval(['fromcoords = reflpoints',JJ(kk-1,1:JJnumbofchars(kk-1)),';'])
tp@0 593 startplanes = reflplanesexpand(:,kk-1);
tp@0 594 end
tp@0 595 if kk == jj+1,
tp@0 596 tocoords = toedgecoords;
tp@0 597 endplanes = [planesatedge(lastedgenumbers,1) planesatedge(lastedgenumbers,2)];
tp@0 598 else
tp@0 599 eval(['tocoords = reflpoints',JJ(kk,1:JJnumbofchars(kk)),';'])
tp@0 600 endplanes = reflplanesexpand(:,kk);
tp@0 601 end
tp@0 602
tp@0 603 [nonobstructedpaths,nobstructions] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 604 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 605
tp@0 606 if nobstructions > 0
tp@0 607 masterlisttocheckmore = masterlisttocheckmore(nonobstructedpaths);
tp@0 608 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 609 lastedgenumbers = lastedgenumbers(nonobstructedpaths);
tp@0 610 reflplanesexpand = reflplanesexpand(nonobstructedpaths,:);
tp@0 611 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 612 for ll = 1:jj
tp@0 613 eval(['reflpoints',JJ(ll,1:JJnumbofchars(ll)),' = reflpoints',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0 614 eval(['newIRcoords',JJ(ll,1:JJnumbofchars(ll)),' = newIRcoords',JJ(ll,1:JJnumbofchars(ll)),'(nonobstructedpaths,:);']);
tp@0 615 end
tp@0 616
tp@0 617 end
tp@0 618 ntocheckmore = length(masterlisttocheckmore);
tp@0 619
tp@0 620 end
tp@0 621
tp@0 622 end
tp@0 623 if SHOWTEXT >= 3
tp@0 624 disp([' ',int2str(ntocheckmore),' of the special IS - N-edge - IR combinations survived the obstruction test'])
tp@0 625 end
tp@0 626
tp@0 627 end
tp@0 628
tp@0 629 % Add the found special combinations to the outdata list
tp@0 630
tp@0 631 edgedifflist = [edgedifflist;double(POTENTIALISES(masterlisttocheckmore,ii+1:ii+Ndifforder))-nplanes];
tp@0 632 prespeclist = [prespeclist;[double(POTENTIALISES(masterlisttocheckmore,1:ii)) zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0 633 postspeclist = [postspeclist;[reflplanesexpand zeros(ntocheckmore,specorder-Ndifforder-ii)]];
tp@0 634 bigedgeweightlist = [bigedgeweightlist;[ ISESVISIBILITY(masterlisttocheckmore) edgeweightlist]];
tp@0 635 eval(['validIRcoords = [validIRcoords;newIRcoords',JJ(jj,1:JJnumbofchars(jj)),'];']);
tp@0 636 % NB!! In the same way as earlier, we must a recursive reference
tp@0 637 % method to find the image source of the last specular reflection.
tp@0 638 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterlisttocheckmore)));
tp@0 639 for kk = 2:jj
tp@0 640 ivref = ORIGINSFROM(ivref);
tp@0 641 end
tp@0 642 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 643
tp@0 644 end
tp@0 645
tp@0 646 masterivlist = masterivlist(ivcompletelyOK);
tp@0 647 possibleedgecombs = possibleedgecombs(ivcompletelyOK,:);
tp@0 648 possibleprespecs = possibleprespecs(ivcompletelyOK,:);
tp@0 649 possiblepostspecs = possiblepostspecs(ivcompletelyOK,:);
tp@0 650 possibleweights = possibleweights(ivcompletelyOK,:);
tp@0 651
tp@0 652 nposs = length(ivcompletelyOK);
tp@0 653
tp@0 654 if SHOWTEXT >= 3
tp@0 655 disp([' ',int2str(nposs),' IS + Edge (Nth-order) + IR segments (non-special combinations) survived the obstruction test'])
tp@0 656 end
tp@0 657
tp@0 658 % Add the found "standard" combinations to the outdata list
tp@0 659
tp@0 660 edgedifflist = [edgedifflist;possibleedgecombs];
tp@0 661 prespeclist = [prespeclist; [possibleprespecs zeros(nposs,specorder-Ndifforder-ii)]];
tp@0 662 postspeclist = [postspeclist;[possiblepostspecs zeros(nposs,specorder-Ndifforder-jj)]];
tp@0 663 bigedgeweightlist = [bigedgeweightlist;[possibleweights bigedgeweightlistin(ivOK(ivreftoindata(ivcompletelyOK)))]];
tp@0 664 validIRcoords = [validIRcoords;validEDIRcoords(ivOK(ivreftoindata(ivcompletelyOK)),:)];
tp@0 665 % NB!! In the same way as earlier, we must a recursive reference
tp@0 666 % method to find the image source of the last specular reflection.
tp@0 667 ivref = ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(ORIGINSFROM(masterivlist)))));
tp@0 668 for kk = 6:jj+Ndifforder
tp@0 669 ivref = ORIGINSFROM(ivref);
tp@0 670 end
tp@0 671 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 672
tp@0 673 end
tp@0 674 end
tp@0 675
tp@0 676
tp@0 677 % #######################################################################
tp@0 678 % #
tp@0 679 % # The weights of the visible edge segments should be
tp@0 680 % # translated into start and end points, together with the visibility
tp@0 681 % # weights from the receiver.
tp@0 682 % # This can be done for all combinations at once!
tp@0 683 % #
tp@0 684 % #######################################################################
tp@0 685
tp@0 686 % As a start we set the start and end values to 0 and 1, i.e. assuming full
tp@0 687 % visibility.
tp@0 688
tp@0 689 addvec = [];
tp@0 690 for ii = 1:Ndifforder
tp@0 691 addvec = [addvec,zeros(length(edgedifflist),1),ones(length(edgedifflist),1)];
tp@0 692 end
tp@0 693 startandendpoints = [startandendpoints addvec];
tp@0 694
tp@0 695 % #######################################################################
tp@0 696 % #
tp@0 697 % # Construct a list guide, which will tell which rows have only
tp@0 698 % # dd, which rows have sdd etc
tp@0 699 % # Syntax: dd,sdd,ssdd,sssdd,...,dds,ddss,ddsss,...
tp@0 700 % # (should also continue with sdds, sddss,...
tp@0 701 % #
tp@0 702 % #######################################################################
tp@0 703
tp@0 704 [n1,n2] = size(prespeclist);
tp@0 705 if n2 > 1
tp@0 706 nprespecs = sum(prespeclist.' > 0).';
tp@0 707 else
tp@0 708 nprespecs = (prespeclist>0);
tp@0 709 end
tp@0 710 [n1,n2] = size(postspeclist);
tp@0 711 if n2 > 1
tp@0 712 npostspecs = sum(postspeclist.' > 0).';
tp@0 713 else
tp@0 714 npostspecs = (postspeclist>0);
tp@0 715 end
tp@0 716
tp@0 717 listofdiffcol = 1 + nprespecs;
tp@0 718 listofreflorder = listofdiffcol + npostspecs + Ndifforder - 1;
tp@0 719
tp@0 720 [B,ivec,jvec] = unique([listofreflorder listofdiffcol],'rows');
tp@0 721 nuniquecombs = length(ivec);
tp@0 722 ntotcombs = length(jvec);
tp@0 723
tp@0 724 listguide = zeros(nuniquecombs,3);
tp@0 725 listoforders = zeros(nuniquecombs,2);
tp@0 726 sortvec = zeros(ntotcombs,1);
tp@0 727 for ii = 1:length(ivec)
tp@0 728 ivfindcombs = find(jvec==ii);
tp@0 729 listguide(ii,1) = length(ivfindcombs);
tp@0 730 if ii > 1
tp@0 731 listguide(ii,2) = listguide(ii-1,3)+1;
tp@0 732 else
tp@0 733 listguide(ii,2) = 1;
tp@0 734 end
tp@0 735 listguide(ii,3) = listguide(ii,2)+listguide(ii,1)-1;
tp@0 736 listoforders(ii,:) = [B(ii,1) B(ii,2)];
tp@0 737
tp@0 738 sortvec(listguide(ii,2):listguide(ii,3)) = ivfindcombs;
tp@0 739
tp@0 740 end
tp@0 741
tp@0 742 prespeclist = prespeclist(sortvec,:);
tp@0 743 postspeclist = postspeclist(sortvec,:);
tp@0 744 validIScoords = validIScoords(sortvec,:);
tp@0 745 validIRcoords = validIRcoords(sortvec,:);
tp@0 746 edgedifflist = edgedifflist(sortvec,:);
tp@0 747 startandendpoints = startandendpoints(sortvec,:);