tp@0: function [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_pointtoedge(fromcoordsshortlist,reftoshortlist,tocoords,canplaneobstruct,planeseesplane,...
tp@0: planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0: % EDB1checkobstr_pointtoedge - Checks obstruction for a list of paths from points to edges.
tp@0: % Checks if the paths from N starting points (reftoshortlist -> fromcoordsshortlist) to
tp@0: % N ending points on edges (tocoords) pass through any of M planes (those marked by 1 in
tp@0: % canplaneobstruct), i.e., if they are obstructed.
tp@0: %
tp@0: % Input parameters:
tp@0: % fromcoordsshortlist Matrix with the unique coordinates of the N
tp@0: % starting points.
tp@0: % reftoshortlist List, [N,1], with the N indices to fromcoordsshortlist.
tp@0: % tocoords Matrix, [N,3], with the coordinates of the
tp@0: % N ending points (on edges) of the N paths to check.
tp@0: % canplaneobstruct,planeseesplane,planeeqs,planenvecs,minvals,maxvals,...
tp@0: % planecorners,corners,ncornersperplanevec,rearsideplane
tp@0: % Data that should have been passed on from the
tp@0: % eddatafile. See EDB1edgeo for more information.
tp@0: %
tp@0: % Output parameters:
tp@0: % nonobstructedpaths A list, [N_nobstructions,1] of the indices of paths that are not
tp@0: % obstructed. The indices refer to the input
tp@0: % lists fromcoords, tocoords, startplanes
tp@0: % endplanes
tp@0: % nobstructions The number of paths (of the N input paths) that
tp@0: % were not obstructed.
tp@0: % edgehits A list, [N_edgehits,1] of the indicies of paths that hit a
tp@0: % plane right at an edge.
tp@0: % cornerhits A list, [N_cornerhits,1] of the indicies of paths that hit a
tp@0: % plane right at a corner.
tp@0: %
tp@0: % Uses function EDB1poinpla
tp@0: %
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
tp@0: %
tp@0: % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
tp@0: % it under the terms of the GNU General Public License as published by the Free Software
tp@0: % Foundation, either version 3 of the License, or (at your option) any later version.
tp@0: %
tp@0: % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
tp@0: % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
tp@0: % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
tp@0: %
tp@0: % You should have received a copy of the GNU General Public License along with the
tp@0: % Edge Diffraction Toolbox. If not, see .
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20050922
tp@0: %
tp@0: % [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_pointtoedge(fromcoordsshortlist,reftoshortlist,tocoords,canplaneobstruct,planeseesplane,...
tp@0: % planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0:
tp@0: npaths = size(tocoords,1);
tp@0: nplanes = size(canplaneobstruct,2);
tp@0:
tp@0: % We only need to check planes for which canplaneobstruct = 1.
tp@0:
tp@0: if nplanes < 256
tp@0: maxlistofplanestocheck = uint8( find( sum(canplaneobstruct) ).' );
tp@0: elseif nplanes < 65536
tp@0: maxlistofplanestocheck = uint16( find( sum(canplaneobstruct) ).' );
tp@0: else
tp@0: maxlistofplanestocheck = uint32( find( sum(canplaneobstruct) ).' );
tp@0: end
tp@0:
tp@0: nplanestocheck = length(maxlistofplanestocheck);
tp@0:
tp@0: if nplanestocheck > 0
tp@0: ntot = nplanestocheck*npaths;
tp@0:
tp@0: nmaxcorners = double(max(ncornersperplanevec(maxlistofplanestocheck)));
tp@0: planecorners = planecorners(:,1:nmaxcorners);
tp@0:
tp@0: % For the obstruction test, we create expanded matrices/lists where the
tp@0: % plane-number is the fastest varying variable. We calculate
tp@0: % the direction vectors, from point to point, before the expanding.
tp@0:
tp@0: dirvec = tocoords - fromcoordsshortlist(reftoshortlist,:);
tp@0: clear tocoords
tp@0: dirveclengths = sqrt(sum(dirvec.'.^2)).' + eps*100;
tp@0: dirvec = dirvec./dirveclengths(:,ones(1,3));
tp@0:
tp@0: % Create the expanded matrices, by duplicating each path with the
tp@0: % number of planes that need to be checked obstruction against.
tp@0: % Immediately select only those individual combinations where
tp@0: % canplaneobstruct = 1. This selection is done both for bigplanelist
tp@0: % and for the reftoshortlist, and for the reftodirvec (which points to
tp@0: % the shorter dirvec matrix).
tp@0:
tp@0: iv1 = (find(canplaneobstruct(:,maxlistofplanestocheck)));
tp@0:
tp@0: bigplanelist = maxlistofplanestocheck(ceil((iv1)/npaths));
tp@0: if ntot < 65536
tp@0: iv1 = uint16(iv1);
tp@0: elseif ntot < 4.29e9
tp@0: iv1 = uint32(iv1);
tp@0: end
tp@0:
tp@0: reftoshortlist = reftoshortlist(:,ones(1,nplanestocheck));
tp@0: reftoshortlist = reftoshortlist(iv1);
tp@0:
tp@0: if npaths*nplanestocheck < 65536
tp@0: reftodirvec = uint16([1:npaths].');
tp@0: else
tp@0: reftodirvec = uint32([1:npaths].');
tp@0: end
tp@0: reftodirvec = reftodirvec(:,ones(1,nplanestocheck));
tp@0: reftodirvec = reftodirvec(iv1);
tp@0:
tp@0: %--------------------------------------------------------------------------
tp@0: % Tests
tp@0: % (Earlier, it was first tested if the IS and R were on the same side of a plane. This is checked at another place).
tp@0: % (Earlier, test 2 was: Are the vectors IS -> R parallel to the planes?
tp@0: % This shold not be needed now).
tp@0: %
tp@0: % Are the hitpoints inside the planes?
tp@0: % 1. Calculate the hit points' coordinates.
tp@0: % 2. Check that the hit point is not further away than the receiver
tp@0: % point, and that the hit point is not in the wrong direction!
tp@0: % 3. Check if the hit points are inside the plane that is checked.
tp@0:
tp@0: % Step 1
tp@0: udir_bigplanenvecs = ( planeeqs(bigplanelist,4) - sum( planenvecs(bigplanelist,:).'.*(fromcoordsshortlist(reftoshortlist,:).') ).' )./( sum( planenvecs(bigplanelist,:).'.*(dirvec(reftodirvec,:).') ).');
tp@0:
tp@0: % Step 2
tp@0: iv2 = find( udir_bigplanenvecs0 );
tp@0: clear dirveclengths
tp@0: iv1 = iv1(iv2);
tp@0: bigplanelist = bigplanelist(iv2);
tp@0: reftoshortlist = reftoshortlist(iv2);
tp@0: reftodirvec = reftodirvec(iv2);
tp@0: udir_bigplanenvecs = udir_bigplanenvecs(iv2,:);
tp@0: clear iv2
tp@0:
tp@0: if ~isempty(iv1)
tp@0: % Step 3
tp@0: xhit_list = fromcoordsshortlist(reftoshortlist,:) + udir_bigplanenvecs(:,ones(1,3)).*dirvec(reftodirvec,:);
tp@0:
tp@0: clear udir_bigplanenvecs dirvec reftoshortlist
tp@0:
tp@0:
tp@0: [hitvec,edgehit,cornerhit] = EDB1poinpla(xhit_list,bigplanelist,minvals,maxvals,planecorners,corners,ncornersperplanevec,planenvecs);
tp@0: clear xhit_list
tp@0:
tp@0: % We recycle canplaneobstruct since we need to compile how many hits
tp@0: % that occurred.
tp@0: canplaneobstruct = zeros(size(canplaneobstruct(:,maxlistofplanestocheck)));
tp@0: ivhits = find(hitvec);
tp@0: if ~isempty(ivhits)
tp@0: indexvec = iv1(ivhits);
tp@0: canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
tp@0: end
tp@0: edgehits = [];
tp@0: ivhits = find(edgehit);
tp@0: if ~isempty(ivhits)
tp@0: disp('WARNING: Edge hit. This is treated as an obstruction.')
tp@0: indexvec = iv1(ivhits);
tp@0: canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
tp@0: end
tp@0:
tp@0: cornerhits = [];
tp@0: ivhits = find(cornerhit);
tp@0: if ~isempty(ivhits)
tp@0: disp('WARNING: Corner hit. This is treated as an obstruction.')
tp@0: indexvec = iv1(ivhits);
tp@0: canplaneobstruct(iv1(ivhits)) = canplaneobstruct(iv1(ivhits))+1;
tp@0: end
tp@0:
tp@0: % Here was a bug: erroneous summing when there was a single plane that
tp@0: % could obstruct (PS 041127)
tp@0: if size(canplaneobstruct,2) > 1
tp@0: canplaneobstruct = 1-sign(sum(canplaneobstruct.'));
tp@0: else
tp@0: canplaneobstruct = 1-canplaneobstruct;
tp@0: end
tp@0: nonobstructedpaths = find(canplaneobstruct);
tp@0:
tp@0: nobstructions = npaths-length(nonobstructedpaths);
tp@0: else
tp@0: nonobstructedpaths = [1:npaths].';
tp@0: nobstructions = 0;
tp@0: edgehits = [];
tp@0: cornerhits = [];
tp@0: end
tp@0: else
tp@0: nonobstructedpaths = [1:npaths].';
tp@0: nobstructions = 0;
tp@0: edgehits = [];
tp@0: cornerhits = [];
tp@0: end