Mercurial > hg > human-echolocation

annotate private/EDB1diffISESx.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,listoforders,...
tp@0 2 bigedgeweightlist] = EDB1diffISESx(eddatafile,S,R,...
tp@0 3 ivsinglediff,singlediffcol,startindicessinglediff,endindicessinglediff,...
tp@0 4 specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,...
tp@0 5 PointertoIRcombs,IRoriginsfrom)
tp@0 6 % EDB1diffISESx - Gives list of paths that includes a 1st-order diff. combination.
tp@0 7 % Gives the list of possible first-order diffraction paths, possibly with specular
tp@0 8 % reflections before and after.
tp@0 9 %
tp@0 10 % Input parameters:
tp@0 11 % eddatafile The name of the file that contains all edge related data.
tp@0 12 % This file will be loaded.
tp@0 13 % S,R,ivsinglediff,...
tp@0 14 % singlediffcol,startindicessinglediff,endindicessinglediff,specorder,visplanesfromR,...
tp@0 15 % vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,PointertoIRcombs,IRoriginsfrom
tp@0 16 % Data that should have been passed from the srdatafile
tp@0 17 % (S,R,visplanesfromR,vispartedgesfromS,vispartedgesfromR
tp@0 18 % )
tp@0 19 % from the ISEStreefile
tp@0 20 % (ivsinglediff
tp@0 21 % singlediffcol,startindicessinglediff,endindicessinglediff
tp@0 22 % ndecimaldivider,PointertoIRcombs,IRoriginsfrom)
tp@0 23 % and from the setupfile
tp@0 24 % (specorder,nedgesubs)
tp@0 25 % POTENTIALISES (global)
tp@0 26 % ISCOORDS (global)
tp@0 27 % ORIGINSFROM (global)
tp@0 28 % ISESVISIBILITY (global)
tp@0 29 % REFLORDER (global)
tp@0 30 %
tp@0 31 % Global parameters:
tp@0 32 % SHOWTEXT,JJ,JJnumbofchars See EDB1mainISES
tp@0 33 % POTENTIALISES,ISCOORDS
tp@0 34 %
tp@0 35 % Output parameters:
tp@0 36 % edgedifflist List [ncombs,1] of the edge number involved in each
tp@0 37 % spec-diff-spec combination.
tp@0 38 % startandendpoints Matrix [ncombs,2] of the relative start and end
tp@0 39 % points of each edge. The values, [0,1], indicate
tp@0 40 % which part of the edge that is visible.
tp@0 41 % prespeclist Matrix [ncombs,specorder-1] of the specular
tp@0 42 % reflections that precede every diffraction.
tp@0 43 % postspeclist Matrix [ncombs,specorder-1] of the specular
tp@0 44 % reflections that follow every diffraction.
tp@0 45 % validIScoords Matrix [ncombs,3] of the image source for each
tp@0 46 % multiple-spec that precedes the diffraction. If
tp@0 47 % there is no spec refl before the diffraction, the
tp@0 48 % value [0 0 0] is given.
tp@0 49 % validIRcoords Matrix [ncombs,3] of the image receiver for each
tp@0 50 % multiple-spec that follows the diffraction. If
tp@0 51 % there is no spec refl after the diffraction, the
tp@0 52 % value [0 0 0] is given.
tp@0 53 % listguide Matrix [nuniquecombs,3] which for each row gives
tp@0 54 % 1. The number of examples in edgefdifflist etc that
tp@0 55 % are the same type of spec-diff-spec comb.
tp@0 56 % 2. The first row number and 3. The last row number.
tp@0 57 % listoforders Matrix [nuniquecombs,2] which for each row gives
tp@0 58 % 1. The reflection order for the spec-diff-spec comb
tp@0 59 % in the same row in listguide.
tp@0 60 % 2. The order of the diffraction in the
tp@0 61 % spec-diff-spec comb.
tp@0 62 % bigedgeweightlist List [ncombs,1] of the visibility of each edge
tp@0 63 % expressed as a number 0 to 2^nedgesubs-1.
tp@0 64 %
tp@0 65 % Uses functions EDB1findis EDB1getedgepoints EDB1chkISvisible EDB1checkobstrpaths
tp@0 66 %
tp@0 67 % ----------------------------------------------------------------------------------------------
tp@0 68 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 69 %
tp@0 70 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 71 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 72 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 73 %
tp@0 74 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 75 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 76 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 77 %
tp@0 78 % You should have received a copy of the GNU General Public License along with the
tp@0 79 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 80 % ----------------------------------------------------------------------------------------------
tp@0 81 % Peter Svensson (svensson@iet.ntnu.no) 20061118
tp@0 82 %
tp@0 83 % [edgedifflist,startandendpoints,prespeclist,postspeclist,validIScoords,validIRcoords,listguide,listoforders,...
tp@0 84 % bigedgeweightlist] = EDB1diffISESx(eddatafile,S,R,...
tp@0 85 % ivsinglediff,singlediffcol,startindicessinglediff,endindicessinglediff,...
tp@0 86 % specorder,visplanesfromR,vispartedgesfromS,vispartedgesfromR,nedgesubs,ndecimaldivider,...
tp@0 87 % PointertoIRcombs,IRoriginsfrom)
tp@0 88
tp@0 89 global SHOWTEXT JJ JJnumbofchars
tp@0 90 global POTENTIALISES ISCOORDS ORIGINSFROM REFLORDER ISESVISIBILITY
tp@0 91
tp@0 92 eval(['load ',eddatafile])
tp@0 93 clear edgeseesplane cornerinfrontofplane
tp@0 94
tp@0 95 [nedges,slask] = size(planesatedge);
tp@0 96 [nplanes,slask] = size(planecorners);
tp@0 97
tp@0 98 edgedifflist = [];
tp@0 99 prespeclist = [];
tp@0 100 postspeclist = [];
tp@0 101 startandendpoints = [];
tp@0 102 bigedgeweightlist = [];
tp@0 103 validIScoords = [];
tp@0 104 validIRcoords = [];
tp@0 105
tp@0 106 [n1,n2] = size(POTENTIALISES);
tp@0 107 if n2 < specorder
tp@0 108 specorder = n2;
tp@0 109 end
tp@0 110
tp@0 111 maxvisibilityvalue = 2^nedgesubs-1;
tp@0 112 zerosvec1 = zeros(1,max([specorder-1 1]));
tp@0 113 zerosvec2 = zeros(1,3);
tp@0 114 listguide = zeros(specorder*2-1,3);
tp@0 115 listoforders = zeros(specorder*2-1,2);
tp@0 116
tp@0 117 obstructtestneeded = (sum(canplaneobstruct)~=0);
tp@0 118 onesvec = ones(1,nedgesubs);
tp@0 119 onesvec3 = ones(1,3);
tp@0 120
tp@0 121 % ###########################################
tp@0 122 % # #
tp@0 123 % # S - spec-spec- edge - R cases #
tp@0 124 % # #
tp@0 125 % # Prespec cases #
tp@0 126 % # #
tp@0 127 % ###########################################
tp@0 128 %
tp@0 129 % Possible edges for S-spec-spec-E-R are seen (at least partly) by the receiver.
tp@0 130 %
tp@0 131 % The visibility doesn't need to be checked since this was done in the
tp@0 132 % ISEStree.
tp@0 133
tp@0 134 % The vector masterivlist will always refer to the original data vector
tp@0 135 % i.e. POTENTIALISES, ORIGINSFROM, REFLORDER etc
tp@0 136 %
tp@0 137 % First we pick out those indices where there was a single diffraction, but
tp@0 138 % skip those with only diffraction (because we dealt with them already).
tp@0 139 % Also, select only those where the diffraction is the last in the sequence
tp@0 140 % of reflections.
tp@0 141
tp@0 142 for ii = 1:specorder, % NB!!! ii = 1 corresponds to zero spec. refl before the diff.
tp@0 143 % ii = 2 corresponds to one spec refl.
tp@0 144 % before the diff etc
tp@0 145
tp@0 146 if SHOWTEXT >= 3
tp@0 147 if ii > 1
tp@0 148 disp([' Checking for ',JJ(ii-1,1:JJnumbofchars(ii-1)),' spec refl before the edge diff'])
tp@0 149 else
tp@0 150 disp([' Checking for 0 spec refl before the edge diff'])
tp@0 151 end
tp@0 152 end
tp@0 153
tp@0 154 iv = uint32(startindicessinglediff(ii):endindicessinglediff(ii));
tp@0 155
tp@0 156 if ~isempty(iv)
tp@0 157
tp@0 158 ivkeep = find(singlediffcol(iv)==ii);
tp@0 159 iv = iv(ivkeep);
tp@0 160 masterivlist = ivsinglediff(iv);
tp@0 161 possibleedges = double(POTENTIALISES(masterivlist,ii)) - nplanes;
tp@0 162
tp@0 163 % Keep only combinations for which the receiver can see the edge
tp@0 164
tp@0 165 ivnotvisiblefromr = find(vispartedgesfromR(possibleedges)==0);
tp@0 166 if ~isempty(ivnotvisiblefromr)
tp@0 167 masterivlist(ivnotvisiblefromr) = [];
tp@0 168 possibleedges(ivnotvisiblefromr) = [];
tp@0 169 end
tp@0 170 % Pick out the pre-specs
tp@0 171
tp@0 172 nposs = length(masterivlist);
tp@0 173
tp@0 174 if nposs > 0
tp@0 175
tp@0 176 if ii > 1
tp@0 177 possiblecombs = POTENTIALISES(masterivlist,1:ii-1);
tp@0 178 reftoIScoords = ORIGINSFROM(masterivlist);
tp@0 179 edgeweightlist = bitand((ISESVISIBILITY(masterivlist)),(vispartedgesfromR(possibleedges)));
tp@0 180 prespeclist = [prespeclist;[possiblecombs zeros(nposs,specorder-ii)]];
tp@0 181 % NB! It is correct below that the indices for the IScoords should be
tp@0 182 % ORIGINSFROM(masterivlist), rather than masterivlist.
tp@0 183 % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0 184 % spec-spec-...-diff combinations and then
tp@0 185 % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0 186 % ends with a diff has no image source.
tp@0 187 validIScoords = [validIScoords;ISCOORDS(ORIGINSFROM(masterivlist),:)];
tp@0 188 else
tp@0 189 edgeweightlist = bitand((vispartedgesfromS(possibleedges)),(vispartedgesfromR(possibleedges)));
tp@0 190 prespeclist = [prespeclist;[zeros(nposs,max([specorder-1 1]))]];
tp@0 191 % For the case of no spec refl before the diffraction, we
tp@0 192 % let the IS get the S coordinates.
tp@0 193 validIScoords = [validIScoords;S(ones(nposs,1),:)];
tp@0 194 end
tp@0 195
tp@0 196 if SHOWTEXT >= 3
tp@0 197 disp([' ',int2str(nposs),' IS+edges valid'])
tp@0 198 end
tp@0 199
tp@0 200 edgedifflist = [edgedifflist;possibleedges];
tp@0 201 postspeclist = [postspeclist;zerosvec1(ones(nposs,1),:)];
tp@0 202 bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0 203 validIRcoords = [validIRcoords;R(ones(nposs,1),:)];
tp@0 204
tp@0 205 end
tp@0 206
tp@0 207 end
tp@0 208 end
tp@0 209
tp@0 210 iv = uint32(find(bigedgeweightlist==0));
tp@0 211 if ~isempty(iv)
tp@0 212 edgedifflist(iv) = [];
tp@0 213 bigedgeweightlist(iv) = [];
tp@0 214 prespeclist(iv,:) = [];
tp@0 215 postspeclist(iv,:) = [];
tp@0 216 validIScoords(iv,:) = [];
tp@0 217 validIRcoords(iv,:) = [];
tp@0 218 end
tp@0 219
tp@0 220 % #######################################################################
tp@0 221 % #######################################################################
tp@0 222 % #######################################################################
tp@0 223 % #######################################################################
tp@0 224 % #######################################################################
tp@0 225
tp@0 226 % #############################################
tp@0 227 % # #
tp@0 228 % # Build an IR tree #
tp@0 229 % # #
tp@0 230 % #############################################
tp@0 231 %
tp@0 232 % For all cases with a specular reflection after an edge diffraction
tp@0 233 % we build an IR tree, analogous to the IS tree
tp@0 234
tp@0 235 if specorder > 1
tp@0 236
tp@0 237 if SHOWTEXT >= 3
tp@0 238 disp([' For the edge+spec combs, an IR tree is built'])
tp@0 239 end
tp@0 240
tp@0 241 % Select all combinations with a single diffraction, anywhere except in
tp@0 242 % the last column, and store the useful data for these (which column
tp@0 243 % is the diffraction, what reflection order, how many specular
tp@0 244 % reflections?)
tp@0 245
tp@0 246 ivselect = uint32(find( singlediffcol~=REFLORDER(ivsinglediff) ));
tp@0 247 masterivlistorig = ivsinglediff(ivselect);
tp@0 248 singlediffcol_select = singlediffcol(ivselect);
tp@0 249 reflorder_select = REFLORDER(ivsinglediff(ivselect));
tp@0 250 clear ivsinglediff
tp@0 251 nspecular_select = uint32(double(reflorder_select) - double(singlediffcol_select));
tp@0 252 nspecular_select_not_empty = 1;
tp@0 253
tp@0 254 % Now we remove all those with a last reflection plane which can not be
tp@0 255 % seen by the receiver
tp@0 256
tp@0 257 if nplanes + nedges < 65536
tp@0 258 lastreflplane = uint16(zeros(length(ivselect),1));
tp@0 259 else
tp@0 260 lastreflplane = uint32(zeros(length(ivselect),1));
tp@0 261 end
tp@0 262 clear ivselect
tp@0 263 for ii = 2:specorder
tp@0 264 iv = find(reflorder_select==ii);
tp@0 265 lastreflplane(iv) = POTENTIALISES(masterivlistorig(iv),ii);
tp@0 266 end
tp@0 267 ivremove = uint32(find(visplanesfromR(lastreflplane)~=2 & visplanesfromR(lastreflplane)~=4 & visplanesfromR(lastreflplane)~=5));
tp@0 268 if ~isempty(ivremove)
tp@0 269 masterivlistorig(ivremove) = [];
tp@0 270 singlediffcol_select(ivremove) = [];
tp@0 271 reflorder_select(ivremove) = [];
tp@0 272 nspecular_select(ivremove) = [];
tp@0 273 lastreflplane(ivremove) = [];
tp@0 274 clear ivremove
tp@0 275 end
tp@0 276
tp@0 277 % Start by calculating all the first-order IR coordinates for
tp@0 278 % all the planes that are visible from the receiver
tp@0 279 IRcoordsstart = zeros(nplanes,3);
tp@0 280 iv = uint32(find(visplanesfromR==2 | visplanesfromR==4 | visplanesfromR==5));
tp@0 281 IRcoordsstart(iv,:) = EDB1findis(R,iv,planeeqs,1,onesvec3);
tp@0 282
tp@0 283 bigIRcoords = (zeros(size(ISCOORDS)));
tp@0 284 IRreflist = zeros(size(iv));
tp@0 285
tp@0 286 % IMPROVE: Not necessary to calculate the IR coordinates for all
tp@0 287 % combinations since many are the same!
tp@0 288
tp@0 289 % Now we make a temporary copy of POTENTIALISES which is displaced to
tp@0 290 % the right since this will make it easier to align the postspec
tp@0 291 % combos.
tp@0 292
tp@0 293 PotentialISESshift = POTENTIALISES;
tp@0 294 for ii = 2:specorder-1
tp@0 295 iv = uint32(find(reflorder_select==ii));
tp@0 296 PotentialISESshift(masterivlistorig(iv),:) = [zeros(length(iv),specorder-ii) PotentialISESshift(masterivlistorig(iv),1:ii)];
tp@0 297 end
tp@0 298
tp@0 299 % Go through ii, which is the number of specular reflections after
tp@0 300 % the diffraction.
tp@0 301
tp@0 302 if ~isempty(nspecular_select)
tp@0 303 for ii = 1:specorder-1
tp@0 304 if SHOWTEXT >= 3
tp@0 305 disp([' order ',int2str(ii)])
tp@0 306 end
tp@0 307
tp@0 308 % We save the index vectors from the different orders so they can
tp@0 309 % be reused in the next stage.
tp@0 310
tp@0 311 iv = uint32(find(nspecular_select == ii));
tp@0 312 listselect = masterivlistorig(iv);
tp@0 313 nlist = length(listselect);
tp@0 314 eval(['iv',JJ(ii,1:JJnumbofchars(ii)),' = iv;'])
tp@0 315
tp@0 316 if ii == 1
tp@0 317 bigIRcoords(listselect,:) = IRcoordsstart( PotentialISESshift(listselect,specorder),:);
tp@0 318
tp@0 319 elseif ii == 2
tp@0 320 bigIRcoords(listselect,:) = EDB1findis(IRcoordsstart(PotentialISESshift(listselect,specorder),:),...
tp@0 321 PotentialISESshift(listselect,specorder-1),planeeqs,nlist,onesvec3);
tp@0 322
tp@0 323 elseif ii >= 3
tp@0 324 ivcombo = double(PotentialISESshift(listselect,specorder));
tp@0 325 for jj = 1:ii-2
tp@0 326 ivcombo = ivcombo + double(PotentialISESshift(listselect,specorder-jj))*ndecimaldivider^(jj);
tp@0 327 end
tp@0 328 IRreflist = PointertoIRcombs(ivcombo);
tp@0 329 ivsimple = find(IRreflist~=0);
tp@0 330 ivnotsimple = find(IRreflist==0);
tp@0 331 newIRcoordssimple = EDB1findis(bigIRcoords(IRreflist(ivsimple),:),...
tp@0 332 PotentialISESshift(listselect(ivsimple),specorder-(ii-1)),planeeqs,size(ivsimple,1),onesvec3);
tp@0 333 newIRcoordsnotsimple = EDB1findis(IRcoordsstart(PotentialISESshift(listselect(ivnotsimple),specorder),:),...
tp@0 334 PotentialISESshift(listselect(ivnotsimple),specorder-1),planeeqs,size(ivnotsimple,1),onesvec3);
tp@0 335 for jj = 2:ii-1
tp@0 336 newIRcoordsnotsimple = EDB1findis(newIRcoordsnotsimple,...
tp@0 337 PotentialISESshift(listselect(ivnotsimple),specorder-jj),planeeqs,size(ivnotsimple,1),onesvec3);
tp@0 338 end
tp@0 339 newIRcoordstoadd = zeros(length(listselect),3);
tp@0 340 newIRcoordstoadd(ivsimple,:) = newIRcoordssimple;
tp@0 341 newIRcoordstoadd(ivnotsimple,:) = newIRcoordsnotsimple;
tp@0 342 bigIRcoords(listselect,:) = newIRcoordstoadd;
tp@0 343
tp@0 344 end
tp@0 345
tp@0 346 end
tp@0 347 clear nspecular_select
tp@0 348 else
tp@0 349 nspecular_select_not_empty = 0;
tp@0 350 end
tp@0 351 end
tp@0 352
tp@0 353 % #######################################################################
tp@0 354 % #######################################################################
tp@0 355 % #######################################################################
tp@0 356 % #######################################################################
tp@0 357 % #######################################################################
tp@0 358
tp@0 359 % ##################################################
tp@0 360 % # #
tp@0 361 % # S - (spec) - edge - spec - spec - R cases #
tp@0 362 % # #
tp@0 363 % # (Pre- and) Postspec cases #
tp@0 364 % # #
tp@0 365 % ##################################################
tp@0 366 %
tp@0 367 % Possible edges for S-E-spec-spec-R are seen (at least partly) by the receiver.
tp@0 368 %
tp@0 369 % The visibility must be checked, and also obstructions, but only
tp@0 370 % for the paths between the edge and the receiver
tp@0 371
tp@0 372 % The vector masterivlist will always refer to the original data vector
tp@0 373 % i.e. PotentialIR, OriginsfromIR, reflorderIR etc
tp@0 374 %
tp@0 375 % First we pick out those indices where there was a single diffraction,
tp@0 376 % which wasn't last in the sequence
tp@0 377
tp@0 378 if specorder > 1 & nspecular_select_not_empty == 1
tp@0 379
tp@0 380 % The ii-value refers to the number of specular reflections after the
tp@0 381 % diffraction
tp@0 382
tp@0 383 for ii = 1:specorder-1
tp@0 384
tp@0 385 if SHOWTEXT >= 3
tp@0 386 disp([' Checking for ',JJ(ii,1:JJnumbofchars(ii)),' spec refl after the edge diff'])
tp@0 387 end
tp@0 388
tp@0 389 % The selection of cases with ii specular reflection after one
tp@0 390 % diffraction have already been done, in the IR-building process.
tp@0 391
tp@0 392 eval(['masterivlist = masterivlistorig(iv',JJ(ii,1:JJnumbofchars(ii)),');'])
tp@0 393
tp@0 394 if nedges < 256
tp@0 395 possibleedges = uint8(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0 396 elseif nedges < 65536
tp@0 397 possibleedges = uint16(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0 398 else
tp@0 399 possibleedges = uint32(double(PotentialISESshift(masterivlist,specorder-ii)) - nplanes);
tp@0 400 end
tp@0 401
tp@0 402 % NB! possiblecombs will contain the specular sequences after the
tp@0 403 % diffraction but they appear in reversed order!
tp@0 404 possiblecombs = PotentialISESshift(masterivlist,specorder-ii+1:specorder);
tp@0 405
tp@0 406 % Specular combinations before the diffraction
tp@0 407 nmaxprespecs = specorder-1-ii;
tp@0 408 if nmaxprespecs == 0
tp@0 409 if nedges+nplanes < 256
tp@0 410 prespeccombs = uint8(zeros(length(masterivlist),1));
tp@0 411 elseif nedges+nplanes < 65536
tp@0 412 prespeccombs = uint16(zeros(length(masterivlist),1));
tp@0 413 else
tp@0 414 prespeccombs = uint32(zeros(length(masterivlist),1));
tp@0 415 end
tp@0 416 nmaxprespecs = 1;
tp@0 417 else
tp@0 418 prespeccombs = PotentialISESshift(masterivlist,1:nmaxprespecs);
tp@0 419 end
tp@0 420 % Visibility of the edge in each sequence, seen from the source
tp@0 421 edgeweightsfromsou = ISESVISIBILITY(masterivlist);
tp@0 422
tp@0 423 %--------------------------------------------------------------
tp@0 424 % Expand to take the edge subdivisions into account -
tp@0 425 % treat all the edge subdivisions as separate edges for the
tp@0 426 % visibility and obstruction test.
tp@0 427
tp@0 428 nposs = length(masterivlist);
tp@0 429 nposs = nposs*nedgesubs;
tp@0 430
tp@0 431 expandedposscombs = repmat(possiblecombs,[1,nedgesubs]);
tp@0 432 expandedposscombs = reshape(expandedposscombs.',ii,nposs).';
tp@0 433
tp@0 434 expandedprespecs = repmat(prespeccombs,[1,nedgesubs]);
tp@0 435 expandedprespecs = reshape(expandedprespecs.',nmaxprespecs,nposs).';
tp@0 436
tp@0 437 expandededgeweightsfromsou = repmat(edgeweightsfromsou,[1,nedgesubs]);
tp@0 438 expandededgeweightsfromsou = reshape(expandededgeweightsfromsou.',1,nposs).';
tp@0 439
tp@0 440 expandedpossedges = repmat(possibleedges,[1,nedgesubs]);
tp@0 441 expandedpossedges = reshape(expandedpossedges.',1,nposs).';
tp@0 442
tp@0 443 expandedmasterivlist = repmat(masterivlist,[1,nedgesubs]);
tp@0 444 expandedmasterivlist = reshape(expandedmasterivlist.',1,nposs).';
tp@0 445
tp@0 446 if SHOWTEXT >= 3
tp@0 447 disp([' ',int2str(nposs),' Edge+IR segments found initially '])
tp@0 448 end
tp@0 449
tp@0 450 % Go through, iteratively, and check if the path from edge to R passes
tp@0 451 % through all reflection planes along the way
tp@0 452
tp@0 453 for jj = ii:-1:1
tp@0 454
tp@0 455 if nposs > 0
tp@0 456
tp@0 457 if jj == ii
tp@0 458 fromcoords = full(bigIRcoords(expandedmasterivlist,:));
tp@0 459 [toedgecoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(possibleedges,:),edgeendcoords(possibleedges,:),edgelengthvec(possibleedges,:),nedgesubs);
tp@0 460 tocoords = toedgecoords;
tp@0 461 else
tp@0 462 eval(['tocoords = reflpoints',JJ(jj+1,1:JJnumbofchars(jj+1)),';'])
tp@0 463 startlistvec = [1:length(expandedmasterivlist)];
tp@0 464 ivref = IRoriginsfrom(expandedmasterivlist);
tp@0 465 for kk = jj:ii-2
tp@0 466 ivnoIRexist = find(ivref==0);
tp@0 467 if ~isempty(ivnoIRexist)
tp@0 468 ivIRexist = startlistvec;
tp@0 469 ivIRexist(ivnoIRexist) = [];
tp@0 470 ivref(ivIRexist) = IRoriginsfrom(ivref(ivIRexist));
tp@0 471 else
tp@0 472 ivref = IRoriginsfrom(ivref);
tp@0 473 end
tp@0 474 end
tp@0 475 ivnoIRexist = find(ivref==0);
tp@0 476 if isempty(ivnoIRexist)
tp@0 477 fromcoords = full(bigIRcoords(ivref,:));
tp@0 478 else
tp@0 479 ivIRexist = startlistvec;
tp@0 480 ivIRexist(ivnoIRexist) = [];
tp@0 481 fromcoords = zeros(length(ivref),3);
tp@0 482 fromcoords(ivIRexist,:) = full(bigIRcoords(ivref(ivIRexist),:));
tp@0 483
tp@0 484 ISEScutout = PotentialISESshift(expandedmasterivlist(ivnoIRexist),specorder-ii+1:specorder);
tp@0 485 newIRcoords = EDB1findis(R,ISEScutout(:,ii),planeeqs,1,onesvec3);
tp@0 486 for kk = 2:jj
tp@0 487 newIRcoords = EDB1findis(newIRcoords,ISEScutout(:,ii-kk+1),planeeqs,size(newIRcoords,1),onesvec3);
tp@0 488 end
tp@0 489 fromcoords(ivnoIRexist,:) = newIRcoords;
tp@0 490
tp@0 491 end
tp@0 492
tp@0 493 end
tp@0 494
tp@0 495 colno = ii-jj+1;
tp@0 496
tp@0 497 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(fromcoords,tocoords,planeeqs(expandedposscombs(:,colno),4),planenvecs(expandedposscombs(:,colno),:),minvals(expandedposscombs(:,colno),:),...
tp@0 498 maxvals(expandedposscombs(:,colno),:),planecorners(expandedposscombs(:,colno),:),corners,ncornersperplanevec(expandedposscombs(:,colno)));
tp@0 499 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 500 disp('WARNING! An edgehit or cornerhit occurred during the visibility test but this is not')
tp@0 501 disp(' handled correctly yet.')
tp@0 502 end
tp@0 503 eval(['reflpoints',JJ(jj,1:JJnumbofchars(jj)),' = reflpoints;'])
tp@0 504
tp@0 505 expandedmasterivlist = expandedmasterivlist(hitplanes);
tp@0 506 expandedposscombs = expandedposscombs(hitplanes,:);
tp@0 507 expandedpossedges = expandedpossedges(hitplanes);
tp@0 508 expandedprespecs = expandedprespecs(hitplanes,:);
tp@0 509 edgeweightlist = edgeweightlist(hitplanes);
tp@0 510 expandededgeweightsfromsou = expandededgeweightsfromsou(hitplanes);
tp@0 511 toedgecoords = toedgecoords(hitplanes,:);
tp@0 512
tp@0 513 if jj < ii
tp@0 514 for kk = jj+1:ii
tp@0 515 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(hitplanes,:);'])
tp@0 516 end
tp@0 517 end
tp@0 518
tp@0 519 nposs = length(expandedmasterivlist);
tp@0 520
tp@0 521 end
tp@0 522 if SHOWTEXT >= 3
tp@0 523 disp([' ',int2str(nposs),' Edge+IR segments survived the visibility test in refl plane ',JJ(jj,1:JJnumbofchars(jj))])
tp@0 524 if SHOWTEXT >= 5
tp@0 525 POTENTIALISES(expandedmasterivlist,:)
tp@0 526 end
tp@0 527 end
tp@0 528
tp@0 529 end
tp@0 530
tp@0 531 if obstructtestneeded & nposs > 0
tp@0 532
tp@0 533 % Check obstructions for all the paths: edge -> plane1 -> plane2 -> ...
tp@0 534 % -> edge (S to edge is not needed!)
tp@0 535 for jj = 1:ii+1
tp@0 536 if nposs > 0
tp@0 537 if jj==1
tp@0 538 fromcoords = R;
tp@0 539 startplanes = [];
tp@0 540 else
tp@0 541 startplanes = expandedposscombs(:,ii-jj+2);
tp@0 542 eval(['fromcoords = reflpoints',JJ(jj-1,1:JJnumbofchars(jj-1)),';'])
tp@0 543 end
tp@0 544 if jj == ii+1
tp@0 545 tocoords = toedgecoords;
tp@0 546 endplanes = [planesatedge(expandedpossedges,1) planesatedge(expandedpossedges,2)];
tp@0 547 else
tp@0 548 eval(['tocoords = reflpoints',JJ(jj,1:JJnumbofchars(jj)),';'])
tp@0 549 endplanes = expandedposscombs(:,ii-jj+1);
tp@0 550 end
tp@0 551
tp@0 552 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 553 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
tp@0 554 if ~isempty(edgehits) | ~isempty(cornerhits)
tp@0 555 disp('WARNING! An edgehit or cornerhit occurred during the obstruction test but this is not')
tp@0 556 disp(' handled correctly yet.')
tp@0 557 end
tp@0 558
tp@0 559 if nobstructions > 0
tp@0 560 expandedmasterivlist = expandedmasterivlist(nonobstructedpaths);
tp@0 561 expandedposscombs = expandedposscombs(nonobstructedpaths,:);
tp@0 562 expandedpossedges = expandedpossedges(nonobstructedpaths);
tp@0 563 expandedprespecs = expandedprespecs(nonobstructedpaths,:);
tp@0 564 edgeweightlist = edgeweightlist(nonobstructedpaths);
tp@0 565 expandededgeweightsfromsou = expandededgeweightsfromsou(nonobstructedpaths);
tp@0 566 toedgecoords = toedgecoords(nonobstructedpaths,:);
tp@0 567 nposs = length(expandedmasterivlist);
tp@0 568
tp@0 569 for kk = 1:ii
tp@0 570 eval(['reflpoints',JJ(kk,1:JJnumbofchars(kk)),' = reflpoints',JJ(kk,1:JJnumbofchars(kk)),'(nonobstructedpaths,:);'])
tp@0 571 end
tp@0 572
tp@0 573 end
tp@0 574
tp@0 575 end
tp@0 576 if SHOWTEXT >= 3
tp@0 577 disp([' ',int2str(nposs),' Edge+IR segments survived the obstruction test for path segment ',int2str(jj)])
tp@0 578 if SHOWTEXT >= 5
tp@0 579 POTENTIALISES(expandedmasterivlist,:)
tp@0 580 end
tp@0 581 end
tp@0 582 end
tp@0 583
tp@0 584 end
tp@0 585
tp@0 586 if nposs > 0
tp@0 587
tp@0 588 % Visibility of edge segments is given by bitand-ing the
tp@0 589 % edge visibility from the (image) source and from the receiver
tp@0 590
tp@0 591 edgeweightlist = bitand(edgeweightlist,expandededgeweightsfromsou);
tp@0 592 iv = find(edgeweightlist==0);
tp@0 593
tp@0 594 if ~isempty(iv)
tp@0 595 expandedpossedges(iv) = [];
tp@0 596 expandedposscombs(iv,:) = [];
tp@0 597 expandedprespecs(iv,:) = [];
tp@0 598 edgeweightlist(iv) = [];
tp@0 599 expandedmasterivlist(iv) = [];
tp@0 600 nposs = length(edgeweightlist);
tp@0 601
tp@0 602 end
tp@0 603
tp@0 604 % Add the newly found combinations to the outdata list
tp@0 605
tp@0 606 [nnew,slask] = size(expandedprespecs);
tp@0 607 if nnew == 1
tp@0 608 colstoclear = find( expandedprespecs==0 );
tp@0 609 else
tp@0 610 colstoclear = find(sum(expandedprespecs>0)==0);
tp@0 611 end
tp@0 612 if ~isempty(colstoclear)
tp@0 613 expandedprespecs(:,colstoclear) = [];
tp@0 614 end
tp@0 615
tp@0 616 [slask,ncolsexist] = size(prespeclist);
tp@0 617 [slask,ncolsnew] = size(expandedprespecs);
tp@0 618 ncolstoadd = ncolsexist - ncolsnew;
tp@0 619 if ncolstoadd > 0
tp@0 620 prespeclist = [prespeclist; [expandedprespecs zeros(nposs,ncolstoadd)]];
tp@0 621 else
tp@0 622 prespeclist = [prespeclist; expandedprespecs];
tp@0 623 end
tp@0 624 edgedifflist = [edgedifflist;expandedpossedges];
tp@0 625 postspeclist = [postspeclist;[expandedposscombs zeros(nposs,specorder-1-ii)]];
tp@0 626 bigedgeweightlist = [bigedgeweightlist;edgeweightlist];
tp@0 627 validIRcoords = [validIRcoords;bigIRcoords(expandedmasterivlist,:)];
tp@0 628 % NB! It is correct below that the indices for the ISCOORDS should be
tp@0 629 % ORIGINSFROM(masterivlist), rather than masterivlist.
tp@0 630 % The combinations in POTENTIALISES(masterivlist,:) all have
tp@0 631 % spec-spec-...-diff combinations and then
tp@0 632 % ISCOORDS(masterivlist,:) are zeros since a comb. that
tp@0 633 % ends with a diff has no image source.
tp@0 634 % If we have a spec-spec-...-diff-spec comb., then we must
tp@0 635 % use ORIGINSFROM(ORIGINSFROM(masterivlist)) etc.
tp@0 636 ivref = ORIGINSFROM(expandedmasterivlist);
tp@0 637 for kk = 1:ii
tp@0 638 ivref = ORIGINSFROM(ivref);
tp@0 639 end
tp@0 640 validIScoords = [validIScoords;ISCOORDS(ivref,:)];
tp@0 641
tp@0 642 end
tp@0 643 end
tp@0 644
tp@0 645 % clear PotentialISESshift IScoords bigIRcoords
tp@0 646 clear PotentialISESshift bigIRcoords
tp@0 647 end
tp@0 648
tp@0 649 % #######################################################################
tp@0 650 % #
tp@0 651 % # Pack the edge segments together because every little edge segment
tp@0 652 % # is present as a separate edge
tp@0 653 % # This can be done for all combinations at once.
tp@0 654 % #
tp@0 655 % #######################################################################
tp@0 656
tp@0 657 test = [prespeclist edgedifflist postspeclist];
tp@0 658
tp@0 659 if ~isempty(test)
tp@0 660
tp@0 661 ncombs = length(edgedifflist);
tp@0 662 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 663 ivremove = find(dtest==1);
tp@0 664
tp@0 665 while ~isempty(ivremove)
tp@0 666 bigedgeweightlist(ivremove+1) = double(bigedgeweightlist(ivremove+1)) + double(bigedgeweightlist(ivremove));
tp@0 667 bigedgeweightlist(ivremove) = [];
tp@0 668 edgedifflist(ivremove) = [];
tp@0 669 prespeclist(ivremove,:) = [];
tp@0 670 postspeclist(ivremove,:) = [];
tp@0 671 validIScoords(ivremove,:) = [];
tp@0 672 validIRcoords(ivremove,:) = [];
tp@0 673
tp@0 674 test = [prespeclist edgedifflist postspeclist];
tp@0 675 ncombs = length(edgedifflist);
tp@0 676 dtest = diff([0;prod( double(test(1:ncombs-1,:)==test(2:ncombs,:)).' ).']);
tp@0 677 ivremove = find(dtest==1);
tp@0 678 end
tp@0 679 end
tp@0 680
tp@0 681 % #######################################################################
tp@0 682 % #
tp@0 683 % # The weights of the visible edge segments should be
tp@0 684 % # translated into start and end points, together with the visibility
tp@0 685 % # weights from the receiver.
tp@0 686 % # This can be done for all combinations at once.
tp@0 687 % #
tp@0 688 % #######################################################################
tp@0 689
tp@0 690 % As a start we set the start and end values to 0 and 1, i.e. assuming full
tp@0 691 % visibility.
tp@0 692
tp@0 693 if specorder >= 1 & ~isempty(edgedifflist)
tp@0 694 startandendpoints = [startandendpoints;[zeros(length(edgedifflist),1) ones(length(edgedifflist),1)]];
tp@0 695
tp@0 696 totvisibility = bigedgeweightlist;
tp@0 697
tp@0 698 ivnovisibility = find(totvisibility==0);
tp@0 699 if ~isempty(ivnovisibility)
tp@0 700 disp('Clearing some cases.')
tp@0 701 edgedifflist(ivnovisibility) = [];
tp@0 702 prespeclist(ivnovisibility,:) = [];
tp@0 703 postspeclist(ivnovisibility,:) = [];
tp@0 704 bigedgeweightlist(ivnovisibility) = [];
tp@0 705 totvisibility(ivnovisibility) = [];
tp@0 706 validIScoords(ivnovisibility,:) = [];
tp@0 707 validIRcoords(ivnovisibility,:) = [];
tp@0 708 end
tp@0 709
tp@0 710 ivtoaccountfor = [1:length(totvisibility)].';
tp@0 711
tp@0 712 % If there are edges with full visibility we don't need to do anything to
tp@0 713 % the start and end values, but we remove the edges from the list to
tp@0 714 % process.
tp@0 715
tp@0 716 ivwholeedges = find( totvisibility == maxvisibilityvalue);
tp@0 717 if ~isempty(ivwholeedges)
tp@0 718 ivtoaccountfor(ivwholeedges) = [];
tp@0 719 end
tp@0 720
tp@0 721 if ~isempty(ivtoaccountfor)
tp@0 722 ncombs = length(ivtoaccountfor);
tp@0 723 bitpattern = zeros(ncombs,nedgesubs);
tp@0 724 for ii=1:nedgesubs
tp@0 725 bitpattern(:,ii) = bitget(totvisibility(ivtoaccountfor),ii);
tp@0 726 end
tp@0 727 dbit1 = diff([zeros(ncombs,1) bitpattern].').';
tp@0 728 dbit2 = [dbit1 -bitpattern(:,nedgesubs)];
tp@0 729
tp@0 730 nsegments = ceil((sum(abs(dbit1.')).')/2);
tp@0 731
tp@0 732 ivonesegments = find(nsegments==1);
tp@0 733 if ~isempty(ivonesegments)
tp@0 734
tp@0 735 nonesegments = length(ivonesegments);
tp@0 736 multvec = 2.^[0:nedgesubs];
tp@0 737 segstartpos = round(log(sum( ((dbit2(ivonesegments,:)== 1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 738 segendpos = round(log(sum( ((dbit2(ivonesegments,:)==-1).*multvec(ones(nonesegments,1),:)).').')/log(2))+1;
tp@0 739
tp@0 740 ivmodify = find(segstartpos==1);
tp@0 741 segstartpos(ivmodify) = ones(size(ivmodify))*1.5;
tp@0 742 ivmodify = find(segendpos>nedgesubs);
tp@0 743 segendpos(ivmodify) = ones(size(ivmodify))*(nedgesubs+0.5);
tp@0 744
tp@0 745 startandendpoints(ivtoaccountfor(ivonesegments),1) = (segstartpos-1.5)/(nedgesubs-1);
tp@0 746 startandendpoints(ivtoaccountfor(ivonesegments),2) = (segendpos-1.5)/(nedgesubs-1);
tp@0 747
tp@0 748 end
tp@0 749
tp@0 750 % If we have some two-or-more-subsegments cases, they will be
tp@0 751 % discovered by the if-condition below
tp@0 752
tp@0 753 if length(ivonesegments) < ncombs
tp@0 754 for nsegmentstocheck = 2:ceil(nedgesubs/2)
tp@0 755
tp@0 756 ivNsegments = find(nsegments==nsegmentstocheck);
tp@0 757 if ~isempty(ivNsegments)
tp@0 758 [n1,n2] = size(startandendpoints);
tp@0 759 if n2 < 2*nsegmentstocheck
tp@0 760 startandendpoints = [startandendpoints zeros(n1,2*nsegmentstocheck-n2)];
tp@0 761 end
tp@0 762 for jj = 1:length(ivNsegments)
tp@0 763 ivstartbits = find(dbit2(ivNsegments(jj),:) == 1);
tp@0 764 ivstartbits = (ivstartbits==1)*1.5 + (ivstartbits~=1).*ivstartbits;
tp@0 765 ivendbits = find(dbit2(ivNsegments(jj),:) == -1);
tp@0 766 ivendbits = (ivendbits>nedgesubs)*(nedgesubs+0.5) + (ivendbits<=nedgesubs).*ivendbits;
tp@0 767
tp@0 768 for kk = 1:nsegmentstocheck,
tp@0 769 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*2+1) = (ivstartbits(kk)-1.5)/(nedgesubs-1);
tp@0 770 startandendpoints(ivtoaccountfor(ivNsegments(jj)),(kk-1)*2+2) = (ivendbits(kk)-1.5)/(nedgesubs-1);
tp@0 771 end
tp@0 772 end
tp@0 773 end
tp@0 774
tp@0 775 end
tp@0 776 end
tp@0 777 end
tp@0 778 end
tp@0 779
tp@0 780 % #######################################################################
tp@0 781 % #
tp@0 782 % # Construct a list guide, which will tell which rows have only
tp@0 783 % # d, which rows have sd etc
tp@0 784 % #
tp@0 785 % #######################################################################
tp@0 786
tp@0 787 [n1,n2] = size(prespeclist);
tp@0 788 if n1 ~= 0
tp@0 789 if n2 == 0
tp@0 790 prespeclist = zeros(n1,1);
tp@0 791 end
tp@0 792
tp@0 793 if n2 > 1
tp@0 794 nprespecs = sum(prespeclist.' > 0).';
tp@0 795 else
tp@0 796 nprespecs = (prespeclist>0);
tp@0 797 end
tp@0 798 [n1,n2] = size(postspeclist);
tp@0 799 if n2 > 1
tp@0 800 npostspecs = sum(postspeclist.' > 0).';
tp@0 801 else
tp@0 802 npostspecs = (postspeclist>0);
tp@0 803 end
tp@0 804
tp@0 805 listofdiffcol = 1 + nprespecs;
tp@0 806 listofreflorder = listofdiffcol + npostspecs;
tp@0 807
tp@0 808 [B,ivec,jvec] = unique([listofreflorder listofdiffcol],'rows');
tp@0 809 nuniquecombs = length(ivec);
tp@0 810 ntotcombs = length(jvec);
tp@0 811
tp@0 812 listguide = zeros(nuniquecombs,3);
tp@0 813 listoforders = zeros(nuniquecombs,2);
tp@0 814 sortvec = zeros(ntotcombs,1);
tp@0 815 for ii = 1:length(ivec)
tp@0 816 ivfindcombs = find(jvec==ii);
tp@0 817 listguide(ii,1) = length(ivfindcombs);
tp@0 818 if ii > 1
tp@0 819 listguide(ii,2) = listguide(ii-1,3)+1;
tp@0 820 else
tp@0 821 listguide(ii,2) = 1;
tp@0 822 end
tp@0 823 listguide(ii,3) = listguide(ii,2)+listguide(ii,1)-1;
tp@0 824 listoforders(ii,:) = [B(ii,1) B(ii,2)];
tp@0 825
tp@0 826 sortvec(listguide(ii,2):listguide(ii,3)) = ivfindcombs;
tp@0 827
tp@0 828 end
tp@0 829
tp@0 830 prespeclist = prespeclist(sortvec,:);
tp@0 831 postspeclist = postspeclist(sortvec,:);
tp@0 832 validIScoords = validIScoords(sortvec,:);
tp@0 833 validIRcoords = validIRcoords(sortvec,:);
tp@0 834 edgedifflist = edgedifflist(sortvec,:);
tp@0 835 startandendpoints = startandendpoints(sortvec,:);
tp@0 836 bigedgeweightlist = bigedgeweightlist(sortvec,:);
tp@0 837
tp@0 838 % The prespeclist needs to be shifted back to the left because it is
tp@0 839 % "right-aligned"
tp@0 840
tp@0 841 [ncombs,ncols] = size(prespeclist);
tp@0 842 if ncols > 1
tp@0 843 if ncombs > 1
tp@0 844 iv = find( ( prespeclist(:,1)==0 ) & ( prespeclist(:,2)~=0 ) );
tp@0 845 if ~isempty(iv)
tp@0 846 prespeclist(iv,1:ncols-1) = prespeclist(iv,2:ncols);
tp@0 847 prespeclist(iv,ncols) = 0;
tp@0 848 end
tp@0 849
tp@0 850 for jj = 3:ncols
tp@0 851 iv = find( ( sum(prespeclist(:,1:jj-1).').' == 0 ) & ( prespeclist(:,jj)~=0 ) );
tp@0 852 if ~isempty(iv)
tp@0 853 prespeclist(iv,1:ncols-(jj-1)) = prespeclist(iv,jj:ncols);
tp@0 854 prespeclist(iv,ncols-(jj-2):ncols) = 0;
tp@0 855 end
tp@0 856 end
tp@0 857 else
tp@0 858 if prespeclist(1,1) == 0,
tp@0 859 ninitialzeros = find(diff(cumsum(prespeclist)));
tp@0 860 if ~isempty(ninitialzeros)
tp@0 861 ninitialzeros = ninitialzeros(1);
tp@0 862 prespeclist(1:ncols-ninitialzeros) = prespeclist(ninitialzeros+1:ncols);
tp@0 863 prespeclist(ncols-ninitialzeros+1:ncols) = 0;
tp@0 864 end
tp@0 865 end
tp@0 866 end
tp@0 867 end
tp@0 868
tp@0 869 else
tp@0 870 prespeclist = [];
tp@0 871 postspeclist = [];
tp@0 872 validIScoords = [];
tp@0 873 validIRcoords = [];
tp@0 874 edgedifflist = [];
tp@0 875 startandendpoints = [];
tp@0 876 bigedgeweightlist = [];
tp@0 877 listguide = [];
tp@0 878
tp@0 879 end