tp@0: function [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0:     planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0: % EDB1checkobstrpaths - Checks obstruction for a list of paths, or check if specular reflections are valid.
tp@0: % Checks if the paths from N starting points (fromcoords) to
tp@0: % N ending points (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: % Instead of testing for N to N points, it is possible to specify
tp@0: % 1 starting point and N ending points or the other way around.
tp@0: % Also, it is possible to specify, for each individual starting or ending point
tp@0: % one or two planes that the respective point is directly at
tp@0: % because then, those planes will be excluded from the list of planes to check for.
tp@0: %
tp@0: % Input parameters:
tp@0: %   fromcoords          Matrix, [N,3] or [1,3], with the coordinates of the
tp@0: %                       N (or 1) starting points of the N paths to check.
tp@0: %   tocoords            Matrix, [N,3] or [1,3], with the coordinates of the
tp@0: %                       N (or 1) ending points of the N paths to check.
tp@0: %   startplanes         Matrix, [N,1] or [N,2] or [0,0], with the plane
tp@0: %                       numbers that the starting points are positioned directly at.
tp@0: %                       The size is determined by:
tp@0: %                       [N,1] fromcoords are specular reflection points
tp@0: %                       [N,2] fromcoords are edge points
tp@0: %                       [0,0] fromcoords are image sources or receiver
tp@0: %                       positions.
tp@0: %   endplanes           Matrix with the plane numbers that the ending points are
tp@0: %                       positioned directly at.
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 EDB1chkISvisiblex
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 <http://www.gnu.org/licenses/>.                 
tp@0: % ----------------------------------------------------------------------------------------------
tp@0: % Peter Svensson (svensson@iet.ntnu.no) 20060621
tp@0: %
tp@0: % [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(fromcoords,tocoords,startplanes,endplanes,canplaneobstruct,planeseesplane,...
tp@0: %    planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane)
tp@0: 
tp@0: global BIGFROMCOORDSSHORTLIST REFTOFROMSHORTLIST BIGTOCOORDSSHORTLIST REFTOTOSHORTLIST
tp@0: 
tp@0: [nstartpoints,slask] = size(fromcoords);
tp@0: [nendpoints,slask] = size(tocoords);
tp@0: npaths = max(max([nstartpoints nendpoints]));
tp@0: nplanes = length(canplaneobstruct);
tp@0: 
tp@0: % The startplanes and endplanes can either be one-column or two-column
tp@0: % lists depending on if paths are between specular reflections or edges.
tp@0: 
tp@0: [nplanes1,nplanecols1] = size(startplanes);
tp@0: [nplanes2,nplanecols2] = size(endplanes);
tp@0: 
tp@0: % If startplanes and/or endplanes have been specified, then we should only
tp@0: % need to check planes that are seen by the start- or end-plane. 
tp@0: 
tp@0: if ~isempty(startplanes) & ~isempty(endplanes)
tp@0:     if nplanecols1 == 1 & nplanecols2 == 1
tp@0:         planesareseen = ( sum( planeseesplane(endplanes,:) > 0 ) + sum( planeseesplane(startplanes,:) > 0 ) ) > 0;
tp@0:     elseif nplanecols1 == 2 & nplanecols2 == 2
tp@0:         planesareseen = ( sum( planeseesplane(endplanes(:,1),:) > 0 ) + sum( planeseesplane(endplanes(:,2),:) > 0 ) + ...
tp@0:             sum( planeseesplane(startplanes(:,1),:) > 0 ) + sum( planeseesplane(startplanes(:,2),:) > 0 ) ) > 0;
tp@0:     elseif nplanecols1 == 2 & nplanecols2 == 1
tp@0:         planesareseen = ( sum( planeseesplane(endplanes,:) > 0 ) + ...
tp@0:             sum( planeseesplane(startplanes(:,1),:) > 0 ) + sum( planeseesplane(startplanes(:,2),:) > 0 ) ) > 0;
tp@0:     elseif nplanecols1 == 1 & nplanecols2 == 2
tp@0:         planesareseen = ( sum( planeseesplane(endplanes(:,1),:) > 0 ) + sum( planeseesplane(endplanes(:,2),:) > 0 ) + ...
tp@0:            sum( planeseesplane(startplanes,:) > 0 ) ) > 0;
tp@0:     end
tp@0: else
tp@0:     if ~isempty(startplanes)
tp@0:         if nplanes1 > 1
tp@0:             planesareseen = sum( planeseesplane(startplanes,:) > 0 ) > 0;       
tp@0:         else
tp@0:             planesareseen = planeseesplane(startplanes,:) > 0;
tp@0:         end
tp@0:     elseif ~isempty(endplanes)
tp@0:         if nplanes2 > 1
tp@0:             planesareseen = sum( planeseesplane(endplanes,:) > 0 ) > 0;      
tp@0:         else
tp@0:             planesareseen = planeseesplane(endplanes,:) > 0;                  
tp@0:         end
tp@0:     else
tp@0:         planesareseen = ones(1,nplanes);
tp@0:     end
tp@0: end
tp@0: 
tp@0: % In addition, we only need to check planes for which canplaneobstruct = 1.
tp@0: 
tp@0: if nplanes <= 255
tp@0:     maxlistofplanestocheck = uint8(find((planesareseen & canplaneobstruct)>0).');
tp@0: elseif nplanes <= 65535
tp@0:     maxlistofplanestocheck = uint16(find((planesareseen & canplaneobstruct)>0).');    
tp@0: else
tp@0:     maxlistofplanestocheck = uint32(find((planesareseen & canplaneobstruct)>0).');
tp@0: end
tp@0: 
tp@0: nplanestocheck = length(maxlistofplanestocheck);
tp@0: 
tp@0: if nplanestocheck > 0
tp@0: 	onesvec1 = ones(1,nplanestocheck);
tp@0: 	onesvec2 = ones(npaths*nplanestocheck,1);
tp@0: 	ntot = nplanestocheck*npaths;
tp@0: 
tp@0:     % Make one long list of planes to check obstruction for. It will be aligned
tp@0: 	% with expanded lists of path startpoints and endpoints in addition to startplanes and
tp@0: 	% endplanes. We need to construct a ref. list because some of the
tp@0: 	% combinations will be thrown away and we must be able to refer back to the
tp@0: 	% rectangular, complete matrix of size [nplanestocheck,npaths].
tp@0: 
tp@0:     bigplanelist = maxlistofplanestocheck(:,ones(1,npaths));
tp@0: 	bigplanelist = reshape(bigplanelist,ntot,1);
tp@0: 	
tp@0: 	reftoorigmatrix = uint32([1:npaths*nplanestocheck].');
tp@0:     
tp@0:     if nstartpoints == 1 & nendpoints ~= 1,       
tp@0:         speedyalgorithm = 1;
tp@0:     else
tp@0:         speedyalgorithm = 0;        
tp@0:     end
tp@0: 
tp@0:     % Make one long matrix of path endpoints and path endplanes
tp@0: 	    
tp@0: 	if nendpoints > 1
tp@0:         if speedyalgorithm == 0
tp@0:             BIGTOCOORDSSHORTLIST = tocoords;
tp@0:             if npaths < 65536
tp@0:                 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0:             else
tp@0:                 REFTOTOSHORTLIST = uint32([1:npaths]);                
tp@0:             end
tp@0:           REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0:         end
tp@0:         
tp@0:         if isempty(endplanes)
tp@0:             bigendplanes = [];
tp@0:         else
tp@0:             if nplanecols2 == 1
tp@0:                 bigendplanes = endplanes(:,onesvec1);
tp@0:                 bigendplanes = reshape(bigendplanes.',ntot,1);
tp@0:             else               
tp@0:                 bigendplanes = reshape(repmat(endplanes.',[nplanestocheck,1]),2,ntot).';
tp@0:             end
tp@0:         end
tp@0: 	else
tp@0:         if speedyalgorithm == 0
tp@0:             BIGTOCOORDSSHORTLIST = tocoords;
tp@0:             REFTOTOSHORTLIST = int8(onesvec2);
tp@0:         end
tp@0:         if isempty(endplanes)
tp@0:             bigendplanes = [];
tp@0:         else
tp@0:             if length(endplanes) > 2
tp@0:                 error('ERROR: Some strange error here!')    
tp@0:             end
tp@0:             if nplanecols2 == 1
tp@0:                 bigendplanes = endplanes(onesvec2,:);
tp@0:             else
tp@0:                 bigendplanes = [endplanes(1) endplanes(2)];
tp@0:                 bigendplanes = bigendplanes(onesvec2,:);
tp@0:             end
tp@0:         end
tp@0: 	end
tp@0: 	
tp@0: 	% Make one long matrix of path startplanes and startpoints
tp@0: 	
tp@0: 	if nstartpoints > 1
tp@0:         if speedyalgorithm == 0
tp@0:             BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0:             if npaths < 65536
tp@0:                 REFTOFROMSHORTLIST = uint16([1:npaths]);
tp@0:             else
tp@0:                 REFTOFROMSHORTLIST = uint32([1:npaths]);                
tp@0:             end
tp@0:           REFTOFROMSHORTLIST = reshape(REFTOFROMSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0:         end
tp@0:         if isempty(startplanes)
tp@0:             bigstartplanes = [];
tp@0:         else
tp@0:             if nplanecols1 == 1
tp@0:                 bigstartplanes = startplanes(:,onesvec1);
tp@0:                 bigstartplanes = reshape(bigstartplanes.',ntot,1);
tp@0:             else
tp@0:                 bigstartplanes = reshape(repmat(startplanes.',[nplanestocheck,1]),2,ntot).';
tp@0:             end  
tp@0:         end    
tp@0: 	else
tp@0:         if speedyalgorithm == 0
tp@0:             BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0:             REFTOFROMSHORTLIST = int8(onesvec2);
tp@0:         end
tp@0:         if isempty(startplanes)
tp@0:             bigstartplanes = [];
tp@0:         else
tp@0:             if nplanecols1 == 1
tp@0:                 bigstartplanes = startplanes(onesvec2,:);
tp@0:             else
tp@0:                 bigstartplanes = [startplanes(1) startplanes(2)];
tp@0:                 bigstartplanes = bigstartplanes(onesvec2,:);
tp@0:             end
tp@0:         end
tp@0: 	end
tp@0: 
tp@0:     if speedyalgorithm == 0
tp@0:         clear onesvec2
tp@0:     end
tp@0:     
tp@0:     % Now we can clear all combinations where neither the startplane or the
tp@0:     % endplane sees the obstructing plane!
tp@0:     %
tp@0:     % After a number of combinations have been cleared, the fromcoords and
tp@0:     % tocoords can be expanded and then pruned.
tp@0:     %
tp@0:     % This is the speedy algorithm...
tp@0: 
tp@0:     if nstartpoints == 1 & nendpoints ~= 1,       
tp@0:         if size(bigendplanes,2) == 1
tp@0:             iv = uint32(find(planeseesplane(uint32(double(bigendplanes) + (double(bigplanelist)-1)*nplanes))==0));
tp@0:             bigplanelist(iv) = [];
tp@0:             bigendplanes(iv,:) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             
tp@0:             BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0:             REFTOFROMSHORTLIST = int8(onesvec2(1:length(reftoorigmatrix)));
tp@0:             clear onesvec2
tp@0:             BIGTOCOORDSSHORTLIST = tocoords;
tp@0:             if npaths < 65536
tp@0:                 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0:             else
tp@0:                 REFTOTOSHORTLIST = uint32([1:npaths]);                
tp@0:             end
tp@0:             REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             clear iv
tp@0:         else
tp@0:             ref1to_planeseesplane = uint32(double(bigendplanes,1) + (double(bigplanelist)-1)*nplanes);
tp@0:             ref2to_planeseesplane = uint32(double(bigendplanes,2) + (double(bigplanelist)-1)*nplanes);
tp@0:             iv = uint32(find(planeseesplane(ref1to_planeseesplane)==0 & planeseesplane(ref2to_planeseesplane)==0));
tp@0:             clear refto_planeseesplane
tp@0:             bigplanelist(iv) = [];
tp@0:             bigendplanes(iv,:) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             
tp@0:             BIGFROMCOORDSSHORTLIST = fromcoords;
tp@0:             REFTOFROMSHORTLIST = int8(onesvec2(1:length(reftoorigmatrix)));
tp@0:             clear onesvec2
tp@0:             
tp@0:             BIGTOCOORDSSHORTLIST = tocoords;
tp@0:             if npaths < 65536
tp@0:                 REFTOTOSHORTLIST = uint16([1:npaths]);
tp@0:             else
tp@0:                 REFTOTOSHORTLIST = uint32([1:npaths]);                
tp@0:             end
tp@0:             REFTOTOSHORTLIST = reshape(REFTOTOSHORTLIST(ones(nplanestocheck,1),:),ntot,1);
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             clear iv
tp@0:         end
tp@0:     end
tp@0:     if nstartpoints ~= 1 & nendpoints ~= 1,       
tp@0:         
tp@0:     end
tp@0: 
tp@0:     
tp@0: 	% Don't check the one or two planes that are involved in each path for obstruction
tp@0: 	% or, their respective rearside planes, if they happen to be thin.
tp@0: 	
tp@0: 	if ~isempty(bigstartplanes)
tp@0:         if nplanecols1 == 1
tp@0:             iv = find(bigplanelist==bigstartplanes | bigplanelist==rearsideplane(bigstartplanes));
tp@0:             clear bigstartplanes
tp@0:             bigplanelist(iv) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             if ~isempty(bigendplanes)
tp@0:                 bigendplanes(iv,:) = [];
tp@0:             end
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             REFTOFROMSHORTLIST(iv) = [];
tp@0:         else
tp@0:             iv = find(bigplanelist==bigstartplanes(:,1) | bigplanelist==bigstartplanes(:,2) | ...
tp@0:                       bigplanelist==rearsideplane(bigstartplanes(:,1)) | bigplanelist==rearsideplane(bigstartplanes(:,2)));
tp@0:             clear bigstartplanes
tp@0:             bigplanelist(iv) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             if ~isempty(bigendplanes)
tp@0:                 bigendplanes(iv,:) = [];
tp@0:             end
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             REFTOFROMSHORTLIST(iv) = [];
tp@0:         end
tp@0: 	end
tp@0: 	
tp@0: 	if ~isempty(bigendplanes)
tp@0:         if nplanecols2 == 1
tp@0:             
tp@0:             iv = find(bigplanelist==bigendplanes | bigplanelist==rearsideplane(bigendplanes));
tp@0:             
tp@0:             clear bigendplanes
tp@0:             bigplanelist(iv) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             REFTOFROMSHORTLIST(iv) = [];
tp@0:         else
tp@0:             iv = find(bigplanelist==bigendplanes(:,1)           | bigplanelist==bigendplanes(:,2) | ...
tp@0:                       bigplanelist==rearsideplane(bigendplanes(:,1)) | bigplanelist==rearsideplane(bigendplanes(:,2)));
tp@0:             clear bigendplanes
tp@0:             bigplanelist(iv) = [];
tp@0:             reftoorigmatrix(iv) = [];
tp@0:             REFTOTOSHORTLIST(iv) = [];
tp@0:             REFTOFROMSHORTLIST(iv) = [];
tp@0:         end
tp@0: 	end
tp@0: 
tp@0: [hitplanes,reflpoints,edgehits,edgehitpoints,cornerhits,cornerhitpoints] = EDB1chkISvisiblex(bigplanelist,planeeqs(:,4),planenvecs,minvals,...
tp@0:         maxvals,planecorners,corners,ncornersperplanevec);
tp@0:    
tp@0:     zerosvec1 = zeros(nplanestocheck,npaths);
tp@0: 
tp@0:     obstructionoccurrence = zerosvec1;
tp@0: 	obstructionoccurrence(reftoorigmatrix(hitplanes)) = ones(size(hitplanes));	
tp@0: 	if nplanestocheck > 1
tp@0:         obstructionoccurrence = sum(obstructionoccurrence);
tp@0: 	end
tp@0: 	nonobstructedpaths = find(obstructionoccurrence==0);
tp@0:     
tp@0:     if ~isempty(edgehits)
tp@0:         edgehitoccurrence = zerosvec1;
tp@0:         edgehitoccurrence(reftoorigmatrix(edgehits)) = ones(size(edgehits));
tp@0:         if nplanestocheck > 1
tp@0:             edgehitoccurrence = sum(edgehitoccurrence);
tp@0:         end
tp@0:         edgehits = find(edgehitoccurrence~=0);
tp@0:     end
tp@0: 
tp@0:     if ~isempty(cornerhits)
tp@0:         cornerhitoccurrence = zerosvec1;
tp@0:         cornerhitoccurrence(reftoorigmatrix(cornerhits)) = ones(size(cornerhits));
tp@0:         if nplanestocheck > 1
tp@0:             cornerhitoccurrence = sum(cornerhitoccurrence);
tp@0:         end
tp@0:         cornerhits = find(cornerhitoccurrence~=0);
tp@0:     end
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