Mercurial > hg > human-echolocation

annotate private/EDB1checkobstrpaths.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 [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 2 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0 3 % EDB1checkobstrpaths - Checks obstruction for a list of paths, or check if specular reflections are valid.
tp@0 4 % Checks if the paths from N starting points (fromcoords) to
tp@0 5 % N ending points (tocoords) pass through any of M planes (those marked by 1 in
tp@0 6 % canplaneobstruct), i.e., if they are obstructed.
tp@0 7 %
tp@0 8 % Instead of testing for N to N points, it is possible to specify
tp@0 9 % 1 starting point and N ending points or the other way around.
tp@0 10 % Also, it is possible to specify, for each individual starting or ending point
tp@0 11 % one or two planes that the respective point is directly at
tp@0 12 % because then, those planes will be excluded from the list of planes to check for.
tp@0 13 %
tp@0 14 % Input parameters:
tp@0 15 % fromcoords Matrix, [N,3] or [1,3], with the coordinates of the
tp@0 16 % N (or 1) starting points of the N paths to check.
tp@0 17 % tocoords Matrix, [N,3] or [1,3], with the coordinates of the
tp@0 18 % N (or 1) ending points of the N paths to check.
tp@0 19 % startplanes Matrix, [N,1] or [N,2] or [0,0], with the plane
tp@0 20 % numbers that the starting points are positioned directly at.
tp@0 21 % The size is determined by:
tp@0 22 % [N,1] fromcoords are specular reflection points
tp@0 23 % [N,2] fromcoords are edge points
tp@0 24 % [0,0] fromcoords are image sources or receiver
tp@0 25 % positions.
tp@0 26 % endplanes Matrix with the plane numbers that the ending points are
tp@0 27 % positioned directly at.
tp@0 28 % canplaneobstruct,planeseesplane,planeeqs,planenvecs,minvals,maxvals,..
tp@0 29 % planecorners,corners,ncornersperplanevec,rearsideplane
tp@0 30 % Data that should have been passed on from the
tp@0 31 % eddatafile. See EDB1edgeo for more information.
tp@0 32 %
tp@0 33 % Output parameters:
tp@0 34 % nonobstructedpaths A list, [N_nobstructions,1] of the indices of paths that are not
tp@0 35 % obstructed. The indices refer to the input
tp@0 36 % lists fromcoords, tocoords, startplanes
tp@0 37 % endplanes
tp@0 38 % nobstructions The number of paths (of the N input paths) that
tp@0 39 % were not obstructed.
tp@0 40 % edgehits A list, [N_edgehits,1] of the indicies of paths that hit a
tp@0 41 % plane right at an edge.
tp@0 42 % cornerhits A list, [N_cornerhits,1] of the indicies of paths that hit a
tp@0 43 % plane right at a corner.
tp@0 44 %
tp@0 45 % Uses function EDB1chkISvisiblex
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) 20060621
tp@0 62 %
tp@0 63 % [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0 64 % planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0 65
tp@0 66 global BIGFROMCOORDSSHORTLIST REFTOFROMSHORTLIST BIGTOCOORDSSHORTLIST REFTOTOSHORTLIST
tp@0 67
tp@0 68 [nstartpoints,slask] = size(fromcoords);
tp@0 69 [nendpoints,slask] = size(tocoords);
tp@0 70 npaths = max(max([nstartpoints nendpoints]));
tp@0 71 nplanes = length(canplaneobstruct);
tp@0 72
tp@0 73 % The startplanes and endplanes can either be one-column or two-column
tp@0 74 % lists depending on if paths are between specular reflections or edges.
tp@0 75
tp@0 76 [nplanes1,nplanecols1] = size(startplanes);
tp@0 77 [nplanes2,nplanecols2] = size(endplanes);
tp@0 78
tp@0 79 % If startplanes and/or endplanes have been specified, then we should only
tp@0 80 % need to check planes that are seen by the start- or end-plane.
tp@0 81
tp@0 82 if ~isempty(startplanes) & ~isempty(endplanes)
tp@0 83 if nplanecols1 == 1 & nplanecols2 == 1
tp@0 84 planesareseen = ( sum( planeseesplane(endplanes,:) > 0 ) + sum( planeseesplane(startplanes,:) > 0 ) ) > 0;
tp@0 85 elseif nplanecols1 == 2 & nplanecols2 == 2
tp@0 86 planesareseen = ( sum( planeseesplane(endplanes(:,1),:) > 0 ) + sum( planeseesplane(endplanes(:,2),:) > 0 ) + ...
tp@0 87 sum( planeseesplane(startplanes(:,1),:) > 0 ) + sum( planeseesplane(startplanes(:,2),:) > 0 ) ) > 0;
tp@0 88 elseif nplanecols1 == 2 & nplanecols2 == 1
tp@0 89 planesareseen = ( sum( planeseesplane(endplanes,:) > 0 ) + ...
tp@0 90 sum( planeseesplane(startplanes(:,1),:) > 0 ) + sum( planeseesplane(startplanes(:,2),:) > 0 ) ) > 0;
tp@0 91 elseif nplanecols1 == 1 & nplanecols2 == 2
tp@0 92 planesareseen = ( sum( planeseesplane(endplanes(:,1),:) > 0 ) + sum( planeseesplane(endplanes(:,2),:) > 0 ) + ...
tp@0 93 sum( planeseesplane(startplanes,:) > 0 ) ) > 0;
tp@0 94 end
tp@0 95 else
tp@0 96 if ~isempty(startplanes)
tp@0 97 if nplanes1 > 1
tp@0 98 planesareseen = sum( planeseesplane(startplanes,:) > 0 ) > 0;
tp@0 99 else
tp@0 100 planesareseen = planeseesplane(startplanes,:) > 0;
tp@0 101 end
tp@0 102 elseif ~isempty(endplanes)
tp@0 103 if nplanes2 > 1
tp@0 104 planesareseen = sum( planeseesplane(endplanes,:) > 0 ) > 0;
tp@0 105 else
tp@0 106 planesareseen = planeseesplane(endplanes,:) > 0;
tp@0 107 end
tp@0 108 else
tp@0 109 planesareseen = ones(1,nplanes);
tp@0 110 end
tp@0 111 end
tp@0 112
tp@0 113 % In addition, we only need to check planes for which canplaneobstruct = 1.
tp@0 114
tp@0 115 if nplanes <= 255
tp@0 116 maxlistofplanestocheck = uint8(find((planesareseen & canplaneobstruct)>0).');
tp@0 117 elseif nplanes <= 65535
tp@0 118 maxlistofplanestocheck = uint16(find((planesareseen & canplaneobstruct)>0).');
tp@0 119 else
tp@0 120 maxlistofplanestocheck = uint32(find((planesareseen & canplaneobstruct)>0).');
tp@0 121 end
tp@0 122
tp@0 123 nplanestocheck = length(maxlistofplanestocheck);
tp@0 124
tp@0 125 if nplanestocheck > 0
tp@0 126 onesvec1 = ones(1,nplanestocheck);
tp@0 127 onesvec2 = ones(npaths*nplanestocheck,1);
tp@0 128 ntot = nplanestocheck*npaths;
tp@0 129
tp@0 130 % Make one long list of planes to check obstruction for. It will be aligned
tp@0 131 % with expanded lists of path startpoints and endpoints in addition to startplanes and
tp@0 132 % endplanes. We need to construct a ref. list because some of the
tp@0 133 % combinations will be thrown away and we must be able to refer back to the
tp@0 134 % rectangular, complete matrix of size [nplanestocheck,npaths].
tp@0 135
tp@0 136 bigplanelist = maxlistofplanestocheck(:,ones(1,npaths));
tp@0 137 bigplanelist = reshape(bigplanelist,ntot,1);
tp@0 138
tp@0 139 reftoorigmatrix = uint32([1:npaths*nplanestocheck].');
tp@0 140
tp@0 141 if nstartpoints == 1 & nendpoints ~= 1,
tp@0 142 speedyalgorithm = 1;
tp@0 143 else
tp@0 144 speedyalgorithm = 0;
tp@0 145 end
tp@0 146
tp@0 147 % Make one long matrix of path endpoints and path endplanes
tp@0 148
tp@0 149 if nendpoints > 1
tp@0 150 if speedyalgorithm == 0
tp@0 151 BIGTOCOORDSSHORTLIST = tocoords;
tp@0 152 if npaths < 65536
tp@0 153 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0 154 else
tp@0 155 REFTOTOSHORTLIST = uint32([1:npaths]);
tp@0 156 end
tp@0 157 REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0 158 end
tp@0 159
tp@0 160 if isempty(endplanes)
tp@0 161 bigendplanes = [];
tp@0 162 else
tp@0 163 if nplanecols2 == 1
tp@0 164 bigendplanes = endplanes(:,onesvec1);
tp@0 165 bigendplanes = reshape(bigendplanes.',ntot,1);
tp@0 166 else
tp@0 167 bigendplanes = reshape(repmat(endplanes.',[nplanestocheck,1]),2,ntot).';
tp@0 168 end
tp@0 169 end
tp@0 170 else
tp@0 171 if speedyalgorithm == 0
tp@0 172 BIGTOCOORDSSHORTLIST = tocoords;
tp@0 173 REFTOTOSHORTLIST = int8(onesvec2);
tp@0 174 end
tp@0 175 if isempty(endplanes)
tp@0 176 bigendplanes = [];
tp@0 177 else
tp@0 178 if length(endplanes) > 2
tp@0 179 error('ERROR: Some strange error here!')
tp@0 180 end
tp@0 181 if nplanecols2 == 1
tp@0 182 bigendplanes = endplanes(onesvec2,:);
tp@0 183 else
tp@0 184 bigendplanes = [endplanes(1) endplanes(2)];
tp@0 185 bigendplanes = bigendplanes(onesvec2,:);
tp@0 186 end
tp@0 187 end
tp@0 188 end
tp@0 189
tp@0 190 % Make one long matrix of path startplanes and startpoints
tp@0 191
tp@0 192 if nstartpoints > 1
tp@0 193 if speedyalgorithm == 0
tp@0 194 BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0 195 if npaths < 65536
tp@0 196 REFTOFROMSHORTLIST = uint16([1:npaths]);
tp@0 197 else
tp@0 198 REFTOFROMSHORTLIST = uint32([1:npaths]);
tp@0 199 end
tp@0 200 REFTOFROMSHORTLIST = reshape(REFTOFROMSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0 201 end
tp@0 202 if isempty(startplanes)
tp@0 203 bigstartplanes = [];
tp@0 204 else
tp@0 205 if nplanecols1 == 1
tp@0 206 bigstartplanes = startplanes(:,onesvec1);
tp@0 207 bigstartplanes = reshape(bigstartplanes.',ntot,1);
tp@0 208 else
tp@0 209 bigstartplanes = reshape(repmat(startplanes.',[nplanestocheck,1]),2,ntot).';
tp@0 210 end
tp@0 211 end
tp@0 212 else
tp@0 213 if speedyalgorithm == 0
tp@0 214 BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0 215 REFTOFROMSHORTLIST = int8(onesvec2);
tp@0 216 end
tp@0 217 if isempty(startplanes)
tp@0 218 bigstartplanes = [];
tp@0 219 else
tp@0 220 if nplanecols1 == 1
tp@0 221 bigstartplanes = startplanes(onesvec2,:);
tp@0 222 else
tp@0 223 bigstartplanes = [startplanes(1) startplanes(2)];
tp@0 224 bigstartplanes = bigstartplanes(onesvec2,:);
tp@0 225 end
tp@0 226 end
tp@0 227 end
tp@0 228
tp@0 229 if speedyalgorithm == 0
tp@0 230 clear onesvec2
tp@0 231 end
tp@0 232
tp@0 233 % Now we can clear all combinations where neither the startplane or the
tp@0 234 % endplane sees the obstructing plane!
tp@0 235 %
tp@0 236 % After a number of combinations have been cleared, the fromcoords and
tp@0 237 % tocoords can be expanded and then pruned.
tp@0 238 %
tp@0 239 % This is the speedy algorithm...
tp@0 240
tp@0 241 if nstartpoints == 1 & nendpoints ~= 1,
tp@0 242 if size(bigendplanes,2) == 1
tp@0 243 iv = uint32(find(planeseesplane(uint32(double(bigendplanes) + (double(bigplanelist)-1)*nplanes))==0));
tp@0 244 bigplanelist(iv) = [];
tp@0 245 bigendplanes(iv,:) = [];
tp@0 246 reftoorigmatrix(iv) = [];
tp@0 247
tp@0 248 BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0 249 REFTOFROMSHORTLIST = int8(onesvec2(1:length(reftoorigmatrix)));
tp@0 250 clear onesvec2
tp@0 251 BIGTOCOORDSSHORTLIST = tocoords;
tp@0 252 if npaths < 65536
tp@0 253 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0 254 else
tp@0 255 REFTOTOSHORTLIST = uint32([1:npaths]);
tp@0 256 end
tp@0 257 REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0 258 REFTOTOSHORTLIST(iv) = [];
tp@0 259 clear iv
tp@0 260 else
tp@0 261 ref1to_planeseesplane = uint32(double(bigendplanes,1) + (double(bigplanelist)-1)*nplanes);
tp@0 262 ref2to_planeseesplane = uint32(double(bigendplanes,2) + (double(bigplanelist)-1)*nplanes);
tp@0 263 iv = uint32(find(planeseesplane(ref1to_planeseesplane)==0 & planeseesplane(ref2to_planeseesplane)==0));
tp@0 264 clear refto_planeseesplane
tp@0 265 bigplanelist(iv) = [];
tp@0 266 bigendplanes(iv,:) = [];
tp@0 267 reftoorigmatrix(iv) = [];
tp@0 268
tp@0 269 BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0 270 REFTOFROMSHORTLIST = int8(onesvec2(1:length(reftoorigmatrix)));
tp@0 271 clear onesvec2
tp@0 272
tp@0 273 BIGTOCOORDSSHORTLIST = tocoords;
tp@0 274 if npaths < 65536
tp@0 275 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0 276 else
tp@0 277 REFTOTOSHORTLIST = uint32([1:npaths]);
tp@0 278 end
tp@0 279 REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0 280 REFTOTOSHORTLIST(iv) = [];
tp@0 281 clear iv
tp@0 282 end
tp@0 283 end
tp@0 284 if nstartpoints ~= 1 & nendpoints ~= 1,
tp@0 285
tp@0 286 end
tp@0 287
tp@0 288
tp@0 289 % Don't check the one or two planes that are involved in each path for obstruction
tp@0 290 % or, their respective rearside planes, if they happen to be thin.
tp@0 291
tp@0 292 if ~isempty(bigstartplanes)
tp@0 293 if nplanecols1 == 1
tp@0 294 iv = find(bigplanelist==bigstartplanes | bigplanelist==rearsideplane(bigstartplanes));
tp@0 295 clear bigstartplanes
tp@0 296 bigplanelist(iv) = [];
tp@0 297 reftoorigmatrix(iv) = [];
tp@0 298 if ~isempty(bigendplanes)
tp@0 299 bigendplanes(iv,:) = [];
tp@0 300 end
tp@0 301 REFTOTOSHORTLIST(iv) = [];
tp@0 302 REFTOFROMSHORTLIST(iv) = [];
tp@0 303 else
tp@0 304 iv = find(bigplanelist==bigstartplanes(:,1) | bigplanelist==bigstartplanes(:,2) | ...
tp@0 305 bigplanelist==rearsideplane(bigstartplanes(:,1)) | bigplanelist==rearsideplane(bigstartplanes(:,2)));
tp@0 306 clear bigstartplanes
tp@0 307 bigplanelist(iv) = [];
tp@0 308 reftoorigmatrix(iv) = [];
tp@0 309 if ~isempty(bigendplanes)
tp@0 310 bigendplanes(iv,:) = [];
tp@0 311 end
tp@0 312 REFTOTOSHORTLIST(iv) = [];
tp@0 313 REFTOFROMSHORTLIST(iv) = [];
tp@0 314 end
tp@0 315 end
tp@0 316
tp@0 317 if ~isempty(bigendplanes)
tp@0 318 if nplanecols2 == 1
tp@0 319
tp@0 320 iv = find(bigplanelist==bigendplanes | bigplanelist==rearsideplane(bigendplanes));
tp@0 321
tp@0 322 clear bigendplanes
tp@0 323 bigplanelist(iv) = [];
tp@0 324 reftoorigmatrix(iv) = [];
tp@0 325 REFTOTOSHORTLIST(iv) = [];
tp@0 326 REFTOFROMSHORTLIST(iv) = [];
tp@0 327 else
tp@0 328 iv = find(bigplanelist==bigendplanes(:,1) | bigplanelist==bigendplanes(:,2) | ...
tp@0 329 bigplanelist==rearsideplane(bigendplanes(:,1)) | bigplanelist==rearsideplane(bigendplanes(:,2)));
tp@0 330 clear bigendplanes
tp@0 331 bigplanelist(iv) = [];
tp@0 332 reftoorigmatrix(iv) = [];
tp@0 333 REFTOTOSHORTLIST(iv) = [];
tp@0 334 REFTOFROMSHORTLIST(iv) = [];
tp@0 335 end
tp@0 336 end
tp@0 337
tp@0 338 [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisiblex(bigplanelist,planeeqs(:,4),planenvecs,minvals,...
tp@0 339 maxvals,planecorners,corners,ncornersperplanevec);
tp@0 340
tp@0 341 zerosvec1 = zeros(nplanestocheck,npaths);
tp@0 342
tp@0 343 obstructionoccurrence = zerosvec1;
tp@0 344 obstructionoccurrence(reftoorigmatrix(hitplanes)) = ones(size(hitplanes));
tp@0 345 if nplanestocheck > 1
tp@0 346 obstructionoccurrence = sum(obstructionoccurrence);
tp@0 347 end
tp@0 348 nonobstructedpaths = find(obstructionoccurrence==0);
tp@0 349
tp@0 350 if ~isempty(edgehits)
tp@0 351 edgehitoccurrence = zerosvec1;
tp@0 352 edgehitoccurrence(reftoorigmatrix(edgehits)) = ones(size(edgehits));
tp@0 353 if nplanestocheck > 1
tp@0 354 edgehitoccurrence = sum(edgehitoccurrence);
tp@0 355 end
tp@0 356 edgehits = find(edgehitoccurrence~=0);
tp@0 357 end
tp@0 358
tp@0 359 if ~isempty(cornerhits)
tp@0 360 cornerhitoccurrence = zerosvec1;
tp@0 361 cornerhitoccurrence(reftoorigmatrix(cornerhits)) = ones(size(cornerhits));
tp@0 362 if nplanestocheck > 1
tp@0 363 cornerhitoccurrence = sum(cornerhitoccurrence);
tp@0 364 end
tp@0 365 cornerhits = find(cornerhitoccurrence~=0);
tp@0 366 end
tp@0 367
tp@0 368 nobstructions = npaths-length(nonobstructedpaths);
tp@0 369 else
tp@0 370 nonobstructedpaths = [1:npaths].';
tp@0 371 nobstructions = 0;
tp@0 372 edgehits = [];
tp@0 373 cornerhits = [];
tp@0 374 end