Mercurial > hg > human-echolocation

annotate private/EDB1checkobstr_edgetoedge.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] = EDB1checkobstr_edgetoedge(canplaneobstruct,planeseesplane,...
tp@0 2 planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0 3 % EDB1checkobstr_edgetoedge - Checks obstruction for a list of edge-to-edge paths.
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 % Input parameters:
tp@0 9 % canplaneobstruct,planeseesplane,planeeqs,planenvecs,minvals,maxvals,...
tp@0 10 % planecorners,corners,ncornersperplanevec,rearsideplane
tp@0 11 % Data that should have been passed on from the
tp@0 12 % eddatafile. See EDB1edgeo for more information.
tp@0 13 % global:
tp@0 14 % FROMCOORDS Matrix, [N,3] or [1,3], with the coordinates of the
tp@0 15 % N (or 1) starting points of the N paths to check.
tp@0 16 % TOCOORDS Matrix, [N,3] or [1,3], with the coordinates of the
tp@0 17 % N (or 1) ending points of the N paths to check.
tp@0 18 % STARTPLANES Matrix, [N,1] or [N,2] or [0,0], with the plane
tp@0 19 % numbers that the starting points are positioned directly at.
tp@0 20 % The size is determined by:
tp@0 21 % [N,1] fromcoords are specular reflection points
tp@0 22 % [N,2] fromcoords are edge points
tp@0 23 % [0,0] fromcoords are image sources or receiver
tp@0 24 % positions.
tp@0 25 % ENDPLANES Matrix with the plane numbers that the ending points are
tp@0 26 % positioned directly at.
tp@0 27 %
tp@0 28 % Output parameters:
tp@0 29 % nonobstructedpaths A list, [N_nobstructions,1] of the indices of paths that are not
tp@0 30 % obstructed. The indices refer to the input
tp@0 31 % lists fromcoords, tocoords, startplanes
tp@0 32 % endplanes
tp@0 33 % nobstructions The number of paths (of the N input paths) that
tp@0 34 % were not obstructed.
tp@0 35 % edgehits A list, [N_edgehits,1] of the indicies of paths that hit a
tp@0 36 % plane right at an edge.
tp@0 37 % cornerhits A list, [N_cornerhits,1] of the indicies of paths that hit a
tp@0 38 % plane right at a corner.
tp@0 39 %
tp@0 40 % Uses no special functions.
tp@0 41 %
tp@0 42 % ----------------------------------------------------------------------------------------------
tp@0 43 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 44 %
tp@0 45 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 46 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 47 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 48 %
tp@0 49 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 50 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 51 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 52 %
tp@0 53 % You should have received a copy of the GNU General Public License along with the
tp@0 54 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 55 % ----------------------------------------------------------------------------------------------
tp@0 56 % Peter Svensson (svensson@iet.ntnu.no) 20050922
tp@0 57 %
tp@0 58 % [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_edgetoedge(canplaneobstruct,planeseesplane,...
tp@0 59 % planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0 60
tp@0 61 global FROMCOORDSSHORTLIST REFTOFROMCOSHO TOCOORDSSHORTLIST REFTOTOCOSHO STARTPLANES ENDPLANES
tp@0 62 global BIGTOCOORDS BIGFROMCOORDS BIGSTARTPLANES BIGENDPLANES
tp@0 63
tp@0 64 [nstartpoints,slask] = size(REFTOFROMCOSHO);
tp@0 65 [nendpoints,slask] = size(REFTOTOCOSHO);
tp@0 66 npaths = max(max([nstartpoints nendpoints]));
tp@0 67 nplanes = length(canplaneobstruct);
tp@0 68
tp@0 69 % The startplanes and endplanes can either be one-column or two-column
tp@0 70 % lists depending on if paths are between specular reflections or edges.
tp@0 71
tp@0 72 [slask,nplanecols1] = size(STARTPLANES);
tp@0 73 [slask,nplanecols2] = size(ENDPLANES);
tp@0 74
tp@0 75 % In addition, we only need to check planes for which canplaneobstruct = 1.
tp@0 76
tp@0 77 if nplanes < 256
tp@0 78 maxlistofplanestocheck = uint8(find( canplaneobstruct>0).');
tp@0 79 elseif nplanes < 65536
tp@0 80 maxlistofplanestocheck = uint16(find( canplaneobstruct>0).');
tp@0 81 else
tp@0 82 maxlistofplanestocheck = uint32(find( canplaneobstruct>0).');
tp@0 83 end
tp@0 84 nplanestocheck = length(maxlistofplanestocheck);
tp@0 85
tp@0 86 if nplanestocheck > 0
tp@0 87
tp@0 88 ntot = nplanestocheck*npaths;
tp@0 89
tp@0 90 % Make one long list of planes to check obstruction for. It will be aligned
tp@0 91 % with expanded lists of path startpoints and endpoints in addition to startplanes and
tp@0 92 % endplanes. We need to construct a ref. list because some of the
tp@0 93 % combinations will be thrown away and we must be able to refer back to the
tp@0 94 % rectangular, complete matrix of size [nplanestocheck,npaths].
tp@0 95
tp@0 96 bigplanelist = maxlistofplanestocheck(:,ones(1,npaths));
tp@0 97 bigplanelist = reshape(bigplanelist,ntot,1);
tp@0 98
tp@0 99 reftoorigmatrix = uint32([1:npaths*nplanestocheck].');
tp@0 100
tp@0 101 % Don't check the one or two planes that are involved in each path for obstruction
tp@0 102 % or, their respective rearside planes, if they happen to be thin.
tp@0 103
tp@0 104 if ~isempty(BIGSTARTPLANES)
tp@0 105 if nplanecols1 == 1
tp@0 106 iv = find(bigplanelist==BIGSTARTPLANES | bigplanelist==rearsideplane(BIGSTARTPLANES));
tp@0 107 bigplanelist(iv) = [];
tp@0 108 BIGSTARTPLANES(iv) = [];
tp@0 109 reftoorigmatrix(iv) = [];
tp@0 110 if ~isempty(BIGENDPLANES)
tp@0 111 BIGENDPLANES(iv,:) = [];
tp@0 112 end
tp@0 113 BIGTOCOORDS(iv,:) = [];
tp@0 114 BIGFROMCOORDS(iv,:) = [];
tp@0 115 else
tp@0 116
tp@0 117 iv = find(bigplanelist==BIGSTARTPLANES(:,1) | bigplanelist==BIGSTARTPLANES(:,2) | ...
tp@0 118 bigplanelist==rearsideplane(BIGSTARTPLANES(:,1)) | bigplanelist==rearsideplane(BIGSTARTPLANES(:,2)));
tp@0 119 bigplanelist(iv) = [];
tp@0 120 BIGSTARTPLANES(iv,:) = [];
tp@0 121 reftoorigmatrix(iv) = [];
tp@0 122 if ~isempty(BIGENDPLANES)
tp@0 123 BIGENDPLANES(iv,:) = [];
tp@0 124 end
tp@0 125 BIGTOCOORDS(iv,:) = [];
tp@0 126 BIGFROMCOORDS(iv,:) = [];
tp@0 127 end
tp@0 128 end
tp@0 129
tp@0 130 if ~isempty(BIGENDPLANES)
tp@0 131 if nplanecols2 == 1
tp@0 132 iv = find(bigplanelist==BIGENDPLANES | bigplanelist==rearsideplane(BIGENDPLANES));
tp@0 133 bigplanelist(iv) = [];
tp@0 134 if ~isempty(BIGSTARTPLANES)
tp@0 135 BIGSTARTPLANES(iv,:) = [];
tp@0 136 end
tp@0 137 BIGENDPLANES(iv,:) = [];
tp@0 138 reftoorigmatrix(iv) = [];
tp@0 139 BIGTOCOORDS(iv,:) = [];
tp@0 140 BIGFROMCOORDS(iv,:) = [];
tp@0 141 else
tp@0 142 iv = find(bigplanelist==BIGENDPLANES(:,1) | bigplanelist==BIGENDPLANES(:,2) | ...
tp@0 143 bigplanelist==rearsideplane(BIGENDPLANES(:,1)) | bigplanelist==rearsideplane(BIGENDPLANES(:,2)));
tp@0 144 bigplanelist(iv) = [];
tp@0 145 if ~isempty(BIGSTARTPLANES)
tp@0 146 BIGSTARTPLANES(iv,:) = [];
tp@0 147 end
tp@0 148 BIGENDPLANES(iv,:) = [];
tp@0 149 reftoorigmatrix(iv) = [];
tp@0 150 BIGTOCOORDS(iv,:) = [];
tp@0 151 BIGFROMCOORDS(iv,:) = [];
tp@0 152 end
tp@0 153 end
tp@0 154
tp@0 155 % [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(BIGFROMCOORDS,BIGTOCOORDS,planeeqs(bigplanelist,4),planenvecs(bigplanelist,:),minvals(bigplanelist,:),...
tp@0 156 % maxvals(bigplanelist,:),planecorners(bigplanelist,:),corners,ncornersperplanevec(bigplanelist));
tp@0 157 %
tp@0 158 % 20050922: The entire EDB1chkISvisible is pasted in here.
tp@0 159
tp@0 160 %##########################################################################
tp@0 161 %##########################################################################
tp@0 162 %##########################################################################
tp@0 163 %##########################################################################
tp@0 164
tp@0 165
tp@0 166 nsou = size(BIGFROMCOORDS,1);
tp@0 167 nplanes = length(bigplanelist);
tp@0 168
tp@0 169 [nrec,slask] = size(BIGTOCOORDS);
tp@0 170 if nrec == 1
tp@0 171 onesvec = uint8(ones(nplanes,1));
tp@0 172 BIGTOCOORDS = BIGTOCOORDS(onesvec,:);
tp@0 173 clear onesvec
tp@0 174 end
tp@0 175
tp@0 176 planenvecsexpanded = planenvecs(bigplanelist,:).';
tp@0 177
tp@0 178 %----------------------------------------------------------------
tp@0 179 % First we check if the IS and R are on the same side of the
tp@0 180 % respective planes. Then the path can not pass through the plane.
tp@0 181
tp@0 182 tempmatrix = corners(planecorners(bigplanelist,1),:);
tp@0 183 sseesplanes = sum( ((BIGFROMCOORDS - tempmatrix).').*planenvecsexpanded ).';
tp@0 184 rseesplanes = sum( ((BIGTOCOORDS - tempmatrix).').*planenvecsexpanded ).';
tp@0 185 iv1 = uint32(find( (sseesplanes>= (-eps*30) )~=(rseesplanes>= (-eps*30) ) ));
tp@0 186 clear sseesplanes rseesplanes
tp@0 187
tp@0 188
tp@0 189 %--------------------------------------------------------------------------
tp@0 190 % Next tests (if there were any planes that had IS and R at different sides)
tp@0 191 % 1. Are the vectors IS -> R parallel to the planes?
tp@0 192 % 2. Are the reflpoints inside the planes?
tp@0 193
tp@0 194 if ~isempty(iv1)
tp@0 195 npossible = length(iv1);
tp@0 196 % Below tempmatrix is the direction vector
tp@0 197 tempmatrix = BIGTOCOORDS - BIGFROMCOORDS;
tp@0 198 BIGTOCOORDS = [];
tp@0 199 dirveclengths = sqrt(sum(tempmatrix.'.^2)).' + eps*100;
tp@0 200 tempmatrix = tempmatrix./dirveclengths(:,ones(1,3));
tp@0 201
tp@0 202 % If test == 0, then the vector S -> R runs along the
tp@0 203 % plane.
tp@0 204
tp@0 205 iv1 = iv1(find(sum( planenvecsexpanded(:,iv1).*(tempmatrix(iv1,:).') ).'~=0));
tp@0 206
tp@0 207 if ~isempty(iv1)
tp@0 208
tp@0 209 % The last test is if the hitpoint is inside the plane
tp@0 210
tp@0 211 % Step 1: calculate the hit point using u = the line parameter (with
tp@0 212 % unit meters) from the source towards the receiver.
tp@0 213
tp@0 214 udir = ( planeeqs(bigplanelist(iv1),4) - sum( planenvecsexpanded(:,iv1).*(BIGFROMCOORDS(iv1,:).') ).' );
tp@0 215
tp@0 216 udir = udir./( sum( planenvecsexpanded(:,iv1).*(tempmatrix(iv1,:).') ).');
tp@0 217
tp@0 218 % Step 2: check that the hit point is at a shorter distance than
tp@0 219 % the receiver point, and that the hit point is not in the opposite
tp@0 220 % direction than the receiver.
tp@0 221
tp@0 222 iv2 = find( udir<(dirveclengths(iv1)*1.001) & udir>=0 );
tp@0 223 clear dirveclengths
tp@0 224 if length(iv2) < length(iv1)
tp@0 225 iv1 = iv1(iv2);
tp@0 226 udir = udir(iv2);
tp@0 227 clear iv2
tp@0 228 end
tp@0 229
tp@0 230 if ~isempty(iv1)
tp@0 231 % Step 3: calculate the actual xyz coordinates of the hit points
tp@0 232 % Now tempmatrix gets the values of the hit points
tp@0 233 tempmatrix = BIGFROMCOORDS(iv1,:) + udir(:,ones(1,3)).*tempmatrix(iv1,:);
tp@0 234
tp@0 235 clear udir
tp@0 236 BIGFROMCOORDS = [];
tp@0 237
tp@0 238 [hitvec,edgehitvec,cornerhitvec] = EDB1poinpla(tempmatrix,iv1,minvals(bigplanelist,:),maxvals(bigplanelist,:),planecorners(bigplanelist,:),corners,ncornersperplanevec(bigplanelist),planenvecsexpanded.');
tp@0 239
tp@0 240 hitplanes = [];
tp@0 241 reflpoints = [];
tp@0 242 edgehits = [];
tp@0 243 edgehitpoints = [];
tp@0 244 cornerhits = [];
tp@0 245 cornerhitpoints = [];
tp@0 246 if any(any(hitvec)) ~=0
tp@0 247 ivhit = find(hitvec==1);
tp@0 248 hitplanes = iv1( ivhit );
tp@0 249 reflpoints = tempmatrix(ivhit,:);
tp@0 250 end
tp@0 251 if ~isempty(edgehitvec)
tp@0 252 ivhit = find(edgehitvec==1);
tp@0 253 edgehits = iv1( ivhit );
tp@0 254 edgehitpoints = tempmatrix(ivhit,:);
tp@0 255 end
tp@0 256 if ~isempty(cornerhitvec)
tp@0 257 ivhit = find(cornerhitvec==1);
tp@0 258 cornerhits = iv1( ivhit );
tp@0 259 cornerhitpoints = tempmatrix(ivhit,:);
tp@0 260 end
tp@0 261 else
tp@0 262 hitplanes = [];
tp@0 263 reflpoints = [];
tp@0 264 edgehits = [];
tp@0 265 edgehitpoints = [];
tp@0 266 cornerhits = [];
tp@0 267 cornerhitpoints = [];
tp@0 268 end
tp@0 269 else
tp@0 270 hitplanes = [];
tp@0 271 reflpoints = [];
tp@0 272 edgehits = [];
tp@0 273 edgehitpoints = [];
tp@0 274 cornerhits = [];
tp@0 275 cornerhitpoints = [];
tp@0 276 end
tp@0 277 else
tp@0 278 hitplanes = [];
tp@0 279 reflpoints = [];
tp@0 280 edgehits = [];
tp@0 281 edgehitpoints = [];
tp@0 282 cornerhits = [];
tp@0 283 cornerhitpoints = [];
tp@0 284 end
tp@0 285
tp@0 286 %##########################################################################
tp@0 287 %##########################################################################
tp@0 288 %##########################################################################
tp@0 289 %##########################################################################
tp@0 290
tp@0 291 zerosvec1 = zeros(nplanestocheck,npaths);
tp@0 292
tp@0 293 obstructionoccurrence = zerosvec1;
tp@0 294 obstructionoccurrence(reftoorigmatrix(hitplanes)) = ones(size(hitplanes));
tp@0 295 if nplanestocheck > 1
tp@0 296 obstructionoccurrence = sum(obstructionoccurrence);
tp@0 297 end
tp@0 298 nonobstructedpaths = find(obstructionoccurrence==0);
tp@0 299
tp@0 300 if ~isempty(edgehits)
tp@0 301 edgehitoccurrence = zerosvec1;
tp@0 302 edgehitoccurrence(reftoorigmatrix(edgehits)) = ones(size(edgehits));
tp@0 303 if nplanestocheck > 1
tp@0 304 edgehitoccurrence = sum(edgehitoccurrence);
tp@0 305 end
tp@0 306 edgehits = find(edgehitoccurrence~=0);
tp@0 307 end
tp@0 308
tp@0 309 if ~isempty(cornerhits)
tp@0 310 cornerhitoccurrence = zerosvec1;
tp@0 311 cornerhitoccurrence(reftoorigmatrix(cornerhits)) = ones(size(cornerhits));
tp@0 312 if nplanestocheck > 1
tp@0 313 cornerhitoccurrence = sum(cornerhitoccurrence);
tp@0 314 end
tp@0 315 cornerhits = find(cornerhitoccurrence~=0);
tp@0 316 end
tp@0 317
tp@0 318 nobstructions = npaths-length(nonobstructedpaths);
tp@0 319 else
tp@0 320 npaths
tp@0 321 nonobstructedpaths = [1:npaths].'
tp@0 322 pause
tp@0 323 nobstructions = 0;
tp@0 324 edgehits = [];
tp@0 325 cornerhits = [];
tp@0 326 end