Mercurial > hg > human-echolocation

annotate private/EDB1chkISvisible.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 [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(ISlist,R,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec)
tp@0 2 % EDB1chkISvisible - Checks if paths from a set of IS to a set of R pass through their refl. planes.
tp@0 3 % EDB1chkISvisible checks if the paths between a number of IS and a single R, or a number of R,
tp@0 4 % pass through their respective reflecting planes.
tp@0 5 %
tp@0 6 % Input parameters:
tp@0 7 % ISlist Matrix, [nIS,3], of the nIS image source coordinates
tp@0 8 % R Receiver coordinates, [1,3] or [nIS,3]
tp@0 9 % planeeqs_lastvalue List [nIS,1] of the fourth value of the plane
tp@0 10 % equations.
tp@0 11 % planenvecs,minvals, maxvals, planecorners, corners, ncornersperplanevec
tp@0 12 % Data that should have been taken from the corresponding
tp@0 13 % variables in the eddatafile,see EDB1edgeo for more information.
tp@0 14 % NB!! The matrices planeeqs, minvals, maxvals
tp@0 15 % have been rearranged so that they have nIS rows, and each
tp@0 16 % row contain the data for one specific plane: the reflecting
tp@0 17 % plane for the IS.
tp@0 18 %
tp@0 19 % Output parameters:
tp@0 20 % hitplanes List, [nhits,1], of the planes that have been hit. The values
tp@0 21 % in this list refer to the index numbers of the input list ISlist
tp@0 22 % hitpoints List, [nhits,3] of the hitpoint coordinates, for the planes that have been hit
tp@0 23 % edgehits List, [nedgehits,1] of the planes that were hit right at an edge.
tp@0 24 % These combinations were marked as hit in the list hitplanes,
tp@0 25 % but the extra information in edgehits can possibly be used.
tp@0 26 % edgehitpoints List, [nedgehits,3] of the hitpoint coordinates for the
tp@0 27 % planes that were hit at an edge.
tp@0 28 % cornerhits List, [ncornerhits,1] of the planes that were hit right at a corner.
tp@0 29 % These combinations were marked as hit in the list hitplanes,
tp@0 30 % but the extra information in edgehits can possibly be used.
tp@0 31 % cornerhitpoints List, [ncornerhits,3] of the hitpoint coordinates for the
tp@0 32 % planes that were hit at a corner.
tp@0 33 %
tp@0 34 % Uses the function EDB1poinpla.
tp@0 35 %
tp@0 36 % ----------------------------------------------------------------------------------------------
tp@0 37 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0 38 %
tp@0 39 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0 40 % it under the terms of the GNU General Public License as published by the Free Software
tp@0 41 % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0 42 %
tp@0 43 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0 44 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0 45 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0 46 %
tp@0 47 % You should have received a copy of the GNU General Public License along with the
tp@0 48 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
tp@0 49 % ----------------------------------------------------------------------------------------------
tp@0 50 % Peter Svensson (svensson@iet.ntnu.no) 20050116
tp@0 51 %
tp@0 52 % [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisible(ISlist,R,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec)
tp@0 53
tp@0 54 nsou = size(ISlist,1);
tp@0 55 nplanes = size(planeeqs_lastvalue,1);
tp@0 56
tp@0 57 if size(planeeqs_lastvalue,2) > 1
tp@0 58 error(['Fix the call to EDB1chkISvisible!!!'])
tp@0 59 end
tp@0 60
tp@0 61 [nrec,slask] = size(R);
tp@0 62 if nrec == 1
tp@0 63 onesvec = uint8(ones(nplanes,1));
tp@0 64 R = R(onesvec,:);
tp@0 65 clear onesvec
tp@0 66 end
tp@0 67
tp@0 68 planenvecs = planenvecs.';
tp@0 69
tp@0 70 %----------------------------------------------------------------
tp@0 71 % First we check if the IS and R are on the same side of the
tp@0 72 % respective planes. Then the path can not pass through the plane.
tp@0 73
tp@0 74 tempmatrix = corners(planecorners(:,1),:);
tp@0 75 sseesplanes = sum( ((ISlist - tempmatrix).').*planenvecs ).';
tp@0 76 rseesplanes = sum( ((R - tempmatrix).').*planenvecs ).';
tp@0 77 iv1 = uint32(find( (sseesplanes>= (-eps*30) )~=(rseesplanes>= (-eps*30) ) ));
tp@0 78 clear sseesplanes rseesplanes
tp@0 79
tp@0 80
tp@0 81 %--------------------------------------------------------------------------
tp@0 82 % Next tests (if there were any planes that had IS and R at different sides)
tp@0 83 % 1. Are the vectors IS -> R parallel to the planes?
tp@0 84 % 2. Are the hitpoints inside the planes?
tp@0 85
tp@0 86 if ~isempty(iv1)
tp@0 87 npossible = length(iv1);
tp@0 88 % Below tempmatrix is the direction vector
tp@0 89 tempmatrix = R - ISlist;
tp@0 90 clear R
tp@0 91 dirveclengths = sqrt(sum(tempmatrix.'.^2)).' + eps*100;
tp@0 92 tempmatrix = tempmatrix./dirveclengths(:,ones(1,3));
tp@0 93
tp@0 94 % If test == 0, then the vector S -> R runs along the
tp@0 95 % plane.
tp@0 96
tp@0 97 iv1 = iv1(find(sum( planenvecs(:,iv1).*(tempmatrix(iv1,:).') ).'~=0));
tp@0 98
tp@0 99 if ~isempty(iv1)
tp@0 100
tp@0 101 % The last test is if the hitpoint is inside the plane
tp@0 102
tp@0 103 % Step 1: calculate the hit point using u = the line parameter (with
tp@0 104 % unit meters) from the source towards the receiver.
tp@0 105 udir = ( planeeqs_lastvalue(iv1) - sum( planenvecs(:,iv1).*(ISlist(iv1,:).') ).' );
tp@0 106 clear planeeqs_lastvalue
tp@0 107 udir = udir./( sum( planenvecs(:,iv1).*(tempmatrix(iv1,:).') ).');
tp@0 108
tp@0 109 % Step 2: check that the hit point is at a shorter distance than
tp@0 110 % the receiver point, and that the hit point is not in the opposite
tp@0 111 % direction than the receiver.
tp@0 112
tp@0 113 iv2 = find( udir<(dirveclengths(iv1)*1.001) & udir>=0 );
tp@0 114 clear dirveclengths
tp@0 115 if length(iv2) < length(iv1)
tp@0 116 iv1 = iv1(iv2);
tp@0 117 udir = udir(iv2);
tp@0 118 clear iv2
tp@0 119 end
tp@0 120
tp@0 121 if ~isempty(iv1)
tp@0 122 % Step 3: calculate the actual xyz coordinates of the hit points
tp@0 123 % Now tempmatrix gets the values of the hit points
tp@0 124 tempmatrix = ISlist(iv1,:) + udir(:,ones(1,3)).*tempmatrix(iv1,:);
tp@0 125
tp@0 126 clear ISlist udir
tp@0 127
tp@0 128 [hitvec,edgehitvec,cornerhitvec] = EDB1poinpla(tempmatrix,iv1,minvals,maxvals,planecorners,corners,ncornersperplanevec,planenvecs.');
tp@0 129
tp@0 130 hitplanes = [];
tp@0 131 hitpoints = [];
tp@0 132 edgehits = [];
tp@0 133 edgehitpoints = [];
tp@0 134 cornerhits = [];
tp@0 135 cornerhitpoints = [];
tp@0 136 if any(any(hitvec)) ~=0
tp@0 137 ivhit = find(hitvec==1);
tp@0 138 hitplanes = iv1( ivhit );
tp@0 139 hitpoints = tempmatrix(ivhit,:);
tp@0 140 end
tp@0 141 if ~isempty(edgehitvec)
tp@0 142 ivhit = find(edgehitvec==1);
tp@0 143 edgehits = iv1( ivhit );
tp@0 144 edgehitpoints = tempmatrix(ivhit,:);
tp@0 145 end
tp@0 146 if ~isempty(cornerhitvec)
tp@0 147 ivhit = find(cornerhitvec==1);
tp@0 148 cornerhits = iv1( ivhit );
tp@0 149 cornerhitpoints = tempmatrix(ivhit,:);
tp@0 150 end
tp@0 151 else
tp@0 152 hitplanes = [];
tp@0 153 hitpoints = [];
tp@0 154 edgehits = [];
tp@0 155 edgehitpoints = [];
tp@0 156 cornerhits = [];
tp@0 157 cornerhitpoints = [];
tp@0 158 end
tp@0 159 else
tp@0 160 hitplanes = [];
tp@0 161 hitpoints = [];
tp@0 162 edgehits = [];
tp@0 163 edgehitpoints = [];
tp@0 164 cornerhits = [];
tp@0 165 cornerhitpoints = [];
tp@0 166 end
tp@0 167 else
tp@0 168 hitplanes = [];
tp@0 169 hitpoints = [];
tp@0 170 edgehits = [];
tp@0 171 edgehitpoints = [];
tp@0 172 cornerhits = [];
tp@0 173 cornerhitpoints = [];
tp@0 174 end