Mercurial > hg > human-echolocation
view private/EDB1chkISvisiblex.m @ 18:2d5f50205527 jabuilder_int tip
Escape the trailing backslash as well
author | Chris Cannam |
---|---|
date | Tue, 30 Sep 2014 16:23:00 +0100 |
parents | 90220f7249fc |
children |
line wrap: on
line source
function [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisiblex(bigplanelist,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec) % EDB1chkISvisiblex - Checks if paths from a set of IS to a set of R pass through resp. refl. planes. Extended version % EDB1chkISvisiblex checks if the paths between a number of IS and a single R, or a number of R, % pass through their respective reflecting planes. % % Input parameters: % bigplanelist List, [nIS,1], of the plane numbers in a long list of % combinations to check. % planeeqs_lastvalue List [nplanes,1] of the fourth value of the plane % equations. % planenvecs,minvals, maxvals, planecorners, corners, ncornersperplanevec % Data that should have been taken from the corresponding % variables in the eddatafile,see EDB1edgeo for more information. % NB!! The matrices planeeqs, minvals, maxvals % have been rearranged so that they have nIS rows, and each % row contain the data for one specific plane: the reflecting % plane for the IS. % BIGFROMCOORDSSHORTLIST,REFTOFROMSHORTLIST,BIGTOCOORDSSHORTLIST,REFTOTOSHORTLIST % (GLOBAL) % % Output parameters: % hitplanes List, [nhits,1], of the planes that have been hit. The values % in this list refer to the index numbers of the input list ISlist % hitpoints List, [nhits,3] of the hitpoint coordinates, for the planes that have been hit % edgehits List, [nedgehits,1] of the planes that were hit right at an edge. % These combinations were marked as hit in the list hitplanes, % but the extra information in edgehits can possibly be used. % edgehitpoints List, [nedgehits,3] of the hitpoint coordinates for the % planes that were hit at an edge. % cornerhits List, [ncornerhits,1] of the planes that were hit right at a corner. % These combinations were marked as hit in the list hitplanes, % but the extra information in edgehits can possibly be used. % cornerhitpoints List, [ncornerhits,3] of the hitpoint coordinates for the % planes that were hit at a corner. % % Uses the function EDB1poinplax. % % ---------------------------------------------------------------------------------------------- % This file is part of the Edge Diffraction Toolbox by Peter Svensson. % % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify % it under the terms of the GNU General Public License as published by the Free Software % Foundation, either version 3 of the License, or (at your option) any later version. % % The Edge Diffraction Toolbox is distributed in the hope that it will be useful, % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. % % You should have received a copy of the GNU General Public License along with the % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>. % ---------------------------------------------------------------------------------------------- % Peter Svensson (svensson@iet.ntnu.no) 20050116 % % [hitplanes,hitpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisiblex(bigplanelist,planeeqs_lastvalue,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec) global BIGFROMCOORDSSHORTLIST REFTOFROMSHORTLIST BIGTOCOORDSSHORTLIST REFTOTOSHORTLIST if isempty(bigplanelist) hitplanes = []; hitpoints = []; edgehits = []; edgehitpoints = []; cornerhits = []; cornerhitpoints = []; return end nsou = size(REFTOFROMSHORTLIST,1); npaths = size(bigplanelist,1); if size(planeeqs_lastvalue,2) > 1 error(['Fix the call to EDB1chkISvisible!!!']) end [nrec,slask] = size(REFTOTOSHORTLIST); planenvecs = planenvecs.'; %---------------------------------------------------------------- % First we check if the IS and R are on the same side of the % respective planes. Then the path can not pass through the plane. if nrec == 1 iv1 = uint32(find( (sum( ((BIGFROMCOORDSSHORTLIST(REFTOFROMSHORTLIST,:) - corners(planecorners(bigplanelist,1),:)).').*planenvecs(:,bigplanelist) ).' >= (-eps*30) )~=( sum( ((BIGTOCOORDSSHORTLIST(ones(npaths,1),:) - corners(planecorners(bigplanelist,1),:)).').*planenvecs(:,bigplanelist) ).'>= (-eps*30) ) )); else iv1 = uint32(find( (sum( ((BIGFROMCOORDSSHORTLIST(REFTOFROMSHORTLIST,:) - corners(planecorners(bigplanelist,1),:)).').*planenvecs(:,bigplanelist) ).' >= (-eps*30) )~=( sum( ((BIGTOCOORDSSHORTLIST(REFTOTOSHORTLIST,:) - corners(planecorners(bigplanelist,1),:)).').*planenvecs(:,bigplanelist) ).'>= (-eps*30) ) )); end %-------------------------------------------------------------------------- % Next tests (if there were any planes that had IS and R at different sides) % 1. Are the vectors IS -> R parallel to the planes? % 2. Are the hitpoints inside the planes? if ~isempty(iv1) npossible = length(iv1); % Below tempmatrix is the direction vector if nrec == 1 tempmatrix = BIGTOCOORDSSHORTLIST(ones(npaths,1),:) - BIGFROMCOORDSSHORTLIST(REFTOFROMSHORTLIST,:); else tempmatrix = BIGTOCOORDSSHORTLIST(REFTOTOSHORTLIST,:) - BIGFROMCOORDSSHORTLIST(REFTOFROMSHORTLIST,:); end dirveclengths = sqrt(sum(tempmatrix.'.^2)).' + eps*100; tempmatrix = tempmatrix./dirveclengths(:,ones(1,3)); % If test == 0, then the vector S -> R runs along the % plane. iv1 = iv1(find(sum( planenvecs(:,bigplanelist(iv1)).*(tempmatrix(iv1,:).') ).'~=0)); if ~isempty(iv1) % The last test is if the hitpoint is inside the plane % Step 1: calculate the hit point using u = the line parameter (with % unit meters) from the source towards the receiver. ISlist = BIGFROMCOORDSSHORTLIST(REFTOFROMSHORTLIST(iv1),:); udir = ( planeeqs_lastvalue(bigplanelist(iv1)) - sum( planenvecs(:,bigplanelist(iv1)).*(ISlist.') ).' ); clear planeeqs_lastvalue udir = udir./( sum( planenvecs(:,bigplanelist(iv1)).*(tempmatrix(iv1,:).') ).'); % Step 2: check that the hit point is at a shorter distance than % the receiver point, and that the hit point is not in the opposite % direction than the receiver. iv2 = find(udir < -eps*1e4); if ~isempty(iv2) disp(['Lowest value (should always be > 0)']) min(udir(iv2)) error('WARNING!!! Check EDB1chkISvisiblex; some hit point is in the wrong direction from the receiver. Check the plane warping tolerances etc.') end iv2 = find( (udir - dirveclengths(iv1)) > eps*1e4); if ~isempty(iv2) disp(['Lowest value (should always be > eps*1e4):']) min(udir(iv2) - dirveclengths(iv1(iv2))) disp(['Problematic hit point is related to plane ',int2str(bigplanelist(iv1(iv2(1)))),' and IS w coords: ']) ISlist(iv2(1),:) error('WARNING!!! Check EDB1chkISvisiblex; some hit point was at a shorter distance than the IS-receiver distance. Check the plane warping tolerances etc.') end % Step 3: calculate the actual xyz coordinates of the hit points % Now tempmatrix gets the values of the hit points tempmatrix = ISlist + udir(:,ones(1,3)).*tempmatrix(iv1,:); clear ISlist udir [hitvec,edgehitvec,cornerhitvec] = EDB1poinplax(bigplanelist,tempmatrix,iv1,minvals,maxvals,planecorners,corners,ncornersperplanevec,planenvecs.'); hitplanes = []; hitpoints = []; edgehits = []; edgehitpoints = []; cornerhits = []; cornerhitpoints = []; if any(any(hitvec)) ~=0 ivhit = find(hitvec==1); hitplanes = iv1( ivhit ); hitpoints = tempmatrix(ivhit,:); end if ~isempty(edgehitvec) ivhit = find(edgehitvec==1); edgehits = iv1( ivhit ); edgehitpoints = tempmatrix(ivhit,:); end if ~isempty(cornerhitvec) ivhit = find(cornerhitvec==1); cornerhits = iv1( ivhit ); cornerhitpoints = tempmatrix(ivhit,:); end else hitplanes = []; hitpoints = []; edgehits = []; edgehitpoints = []; cornerhits = []; cornerhitpoints = []; end else hitplanes = []; hitpoints = []; edgehits = []; edgehitpoints = []; cornerhits = []; cornerhitpoints = []; end