Mercurial > hg > human-echolocation

function outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy)
% EDB1edgeox - Calculates some plane- and edge-related geometrical parameters.
% EDB1edgeo calculates some plane- and edge-related geometrical parameters
% based on corners and planes in a .mat-file generated by EDB1readcad
% but only data that is independent of the source and receiver.
% The output is saved in a .mat-file.
%
% Input parameters:
%	cadgeofile (optional)	The .mat file with the corners and plane data.
%							If it is not specified, a file opening window is presented.
%	outputfile (optional)	The .mat file where the output data will be stored.
%							If it is not specified, an automatic file name is constructed
%							using the same filename stem as the cadgeofile, but ending
%                           with '_edgeo'.
%   open_or_closed_model    'o' or 'c' defining whether the model is open or closed.
%   int_or_ext_model        'i' or 'e' defining whether the model is interior or exterior.
%   specorder               The highest order of any reflection kind (specular and/or diffraction).
%	difforder (optional)	The highest order of diffraction. If it is 0 or 1 then the parameter
%							edgeseespartialedge is not calculated. Default: 1
%	nedgesubs (optional)	The number of segments that each edge will be
%							subdivided into for visibility/obstruction tests. Default: 2
%							NB! When nedgesubs = 2, only the two end points will be checked.
%	firstskipcorner (optional)	All edges including at least one corner with this number or
%                   		higher will be excluded from the calculations. Default: 1e6
%   planeseesplanestrategy (optional)   If this parameter is given the value 1, a plane-to-plane
%                           visibility check is done by checking the plane-midpoint to plane-midpoint
%                           for obstructions, which might be enough for some geometries.
%	SHOWTEXT (global)		An integer value; if it is 4 or higher, then messages will be printed on the
%							screen during calculations.
% Output parameters:
%	outputfile				The name of the outputfile, generated either automatically or specified as
%							an input parameter.
%
% Output parameters in the outputfile:
%    (Taken directly from the cadgeo-file:)
%		corners planecorners  planeeqs planenvecs ncornersperplanevec
%		minvals maxvals planehasindents indentingcorners
%    (Taken directly from the setup-file:)
%       int_or_ext_model
%	 (Calculated by EDB1edgeo:)
%	edgecorners             Matrix [nedges,2] with the corner numbers that
%	                        define each edge, in the reference order.
%   edgestartcoords         Matrix [nedges,3] with the coordinates of the
%                           start corner for each edge.
%   edgeendcoords           Matrix [nedges,3] with the coordinates of the
%                           end corner for each edge.
%   edgelengthvec           List [nedges,1] with the lengths of each edge.
%   planesatedge            Matrix [nedges,2] with the plane numbers of the
%                           two planes that connect to each edge. The first
%                           plane is considered as the reference plane for
%                           each edge. Rule: if the RH thumb is placed
%                           along the edge, from start to end, and the hand
%                           is rotated, the fingers should "come up
%                           through" the reference plane.
%   closwedangvec           List [nedges,1] with the angles in radians of
%                           each edge - for the closed part of each edge.
%   edgenvecs               Matrix [nedges,3] with the normal vectors of
%                           the reference plane for each edge.
%	offedges                List [nlength,1] where nlength = 0-nedges
%                           containing the edge numbers that should not
%                           used.
%   reflfactors             List [nplanes,1] with the reflection factors
%                           of each plane. The only supported values are
%                           +1 for rigid planes, 0 for perfectlt absorbing
%                           planes and -1 for pressure-release planes. The
%                           pressure-release cases might not be implemented
%                           fully.
%   planeisthin             List [nplanes,1] with the value 1 or 0 stating
%                           whether the plane is thin or not.
%   rearsideplane           List [nplanes,1] with the plane number that is at
%                           the rear side. If a plane is non-thin, the
%                           value in this list is zero.
%   planehasindents         List [nplanes,1] with the value 1 or 0 stating
%                           whether a plane has indents or not.
%   indentingedgepairs      Matrix [nindentingedgepairs,2] where each row
%                           gives two connected edges that form an indent
%                           in the shared plane. Consequently, those two
%                           edges can never see each other.
%   canplaneobstruct        List [nplanes,1] with the value 1 or 0 stating
%                           whether a plane has the potential to obstruct
%                           or not.
%   planeseesplane          A matrix, [nplanes,nplanes], with the value 1
%                           0,-1,-2 stating:
%                           1 A plane is front of another plane and could
%                             then potentially see it.
%                           0 A plane is completely behind another plane
%                             but they are not aligned
%                           -1 A plane is aligned with another plane
%                           -2 A plane is back-to-back with another plane.
%   edgesatplane            A matrix, [nplanes,nmaxedgesperplane] which for
%                           row N contains the edge numbers that connect to
%                           to plane N.
%   edgeseesplane           A matrix, [nplanes,nedges], which contains the
%                           values 0,1,-2 or -2, meaning:
%                           0   The edge is completely behind a plane
%                           -1  The edge can never see the plane and the
%                               edge is aligned with the plane
%                           -2  The edge can never see the plane and the
%                               edge belongs to the plane
%                           1   The edge has at least one point in front of the plane.
%
% Uses the functions EDB1coordtrans1
%                    EDB1infrontofplane, EDB1strpend, EDB1cross
%
% ----------------------------------------------------------------------------------------------
%   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) 20090709
%
% outputfile = EDB1edgeox(cadgeofile,outputfile,open_or_closed_model,int_or_ext_model,specorder,difforder,nedgesubs,firstskipcorner,planeseesplanestrategy);

global SHOWTEXT

geomacc = 1e-10;

if nargin < 9
    planeseesplanestrategy = 0;
	if nargin < 8
		firstskipcorner = 1e6;
		if nargin < 7
			nedgesubs = 2;
			if nargin < 6
				difforder = 1;
			else
				if isempty(difforder), difforder = 1; end
			end
		else
			if isempty(nedgesubs), nedgesubs = 2; end
		end
	else
		if isempty(firstskipcorner), firstskipcorner = 1e6; end
	end
end

%---------------------------------------------------------------
% If no cadgeofile was specified, present a file opening window

if isempty(cadgeofile)
	[cadgeofile,filepath] = uigetfile('*.mat','Please select the cadgeofile');
    [filepath,temp1,temp2] = fileparts(filepath);
	if ~isstr(cadgeofile)
		return
	end
    [temp1,filestem,fileext,temp2] = fileparts(cadgeofile);
    cadgeofile = [filepath,filestem,'.mat'];
else
    [filepath,filestem,fileext] = fileparts(cadgeofile);
    if ~isempty(filepath)
        cadgeofile = [[filepath,filesep],filestem,'.mat'];
    else
        cadgeofile = [filestem,'.mat'];
    end
end


%---------------------------------------------------------------
% If no output file was specified, construct an automatic file name

if isempty(outputfile)
	filestem = EDB1strpend(filestem,'_cadgeo');
	outputfile = [[filepath,filesep],filestem,'_eddata.mat'];
end

%---------------------------------------------------------------

eval(['load ',cadgeofile])

ncorners = length(cornernumbers);
[nplanes,maxncornersperplane] = size(planecorners);
minncornersperplane = min(ncornersperplanevec);
nedgesperplanevec = double(ncornersperplanevec);
ncornersperplanevec = double(ncornersperplanevec);

switchoffcorners = find( cornernumbers >= firstskipcorner );

onesvec2 = ones(nplanes,1);
zerosvec1 = zeros(nplanes,1);
zerosvec2 = int8(zeros(nplanes,nplanes));

%----------------------------------------------------------------------------
%
%		GEOMETRICAL ACOUSTICS PARAMETERS
%
%----------------------------------------------------------------------------
%
% Go through all planeabstypes. If it is:
%	'SOFT' or 'soft'  then reflfactors should be -1.
%	'TOTA' or 'tota' then reflfactors should be 0 (short for TOTABS)
%   Otherwise it is assumed to be rigid, so reflfactors = 1.

if SHOWTEXT >= 4
	disp('         Check abstypes')
end

[slask,nchars] = size(planeabstypes);
reflfactors = ones(nplanes,1);
if nchars >= 4
    planeabstypes = lower(full(planeabstypes(:,1:min([4,nchars]))));

    comptxt = 'soft';
    ivpotential = find(planeabstypes(:,1)==comptxt(1));
    if ~isempty(ivpotential)
        comptxt = 'oft';
        compmat = comptxt(ones(length(ivpotential),1),:);
        ivsoft   = ivpotential(find(prod( (planeabstypes(ivpotential,2:4).*compmat).' ).'));
        reflfactors(ivsoft)   = -1*ones(size(ivsoft));
    end

    comptxt = 'tota';
    ivpotential = find(planeabstypes(:,1)==comptxt(1));
    if ~isempty(ivpotential)
        comptxt = 'ota';
        compmat = comptxt(ones(length(ivpotential),1),:);
        ivtotabs =  ivpotential(find(prod( double(planeabstypes(ivpotential,2:4)==compmat).' ).'));
        reflfactors(ivtotabs) = zeros(size(ivtotabs));
    end
end

%----------------------------------------------------------------------------
%
%		EDGE PARAMETERS
%
%----------------------------------------------------------------------------
%
% Check that the planecorners matrix contains no zeros. In such case
% add the circular repetition of coordinates.

if SHOWTEXT >= 4
	disp('         Circular')
end
if sum(sum( planecorners == 0 )) ~= 0
	for ii = 1:nplanes
		iv = find( planecorners(ii,:) ~= 0);
		ncornersatplane = length(iv);
		if ncornersatplane ~= maxncornersperplane
			pattern = planecorners(ii,iv);
			nrepeatings = ceil(maxncornersperplane/ncornersatplane);
			for jj = 1:nrepeatings-1
				pattern = [pattern planecorners(ii,iv)];
			end
			planecorners(ii,:) = pattern(1:maxncornersperplane);
		end
	end
end

%--------------------------------------------------------------------------------
%
% Find all edges
%
% Go through all planes. All consecutive pairs of corners form an edge.
% The list planesatedge gives the one or two planes that connect
% to each edge. If the second plane number, for a certain edge
% is zero, then that specific edge is a plate edge.

edgesatplane = zeros( size(planecorners) );

if SHOWTEXT >= 4
	disp('         Defining edges...')
end

nedgesguess = sum(ncornersperplanevec);
edgecorners = zeros(nedgesguess,2);
tempplanesatedge = zeros(nedgesguess,1);
nedges = 0;
thirdpoint = [];
edgecounter = 0;
multfac = round(10^(floor(log10(nedgesguess))+2));

% First go through all planes, and construct edges from the list of corners
% that define every plane. This way, many edges will occur twice, unless
% a thin plane has no rear side (then the edge occurs a single time)
% or if two thin planes are connected (then the edge occurs four times)

for ii = 1:nplanes
   for jj = 1:nedgesperplanevec(ii)
      edgecounter = edgecounter + 1;
      corner1 = planecorners(ii,jj);
      if jj + 1 <= ncornersperplanevec(ii)
         corner2 = planecorners(ii,jj+1);
      else
         corner2 = planecorners(ii,1);
      end
	  edgecorners(edgecounter,:) = [corner1 corner2];
	  tempplanesatedge(edgecounter) = ii;
   end

end

if edgecounter < nedgesguess
	edgecorners(edgecounter+1:nedgesguess,:) = [];
	tempplanesatedge(edgecounter+1:nedgesguess,:) = [];
end

% To find all duplicates, the edge list is sorted in numerical order
% If an edge is defined twice or four times, then this is an OK edge.
% If an edge is defined a single time, it is connected to a plane which is not connected, and then the
% edge should be switched off.
% If an edge is defined three times, there must be something like a reflector array where one of
% the rear planes is missing.

[test,flipvec] = sort(edgecorners.');
test = test.';
flipvec = flipvec(1,:).';

test = test(:,1)*multfac + test(:,2);
[test,sortvec] = sort(test);
tempplanesatedge = tempplanesatedge(sortvec);
flipvec = flipvec(sortvec);
planesatedge = [tempplanesatedge.*(flipvec==1) tempplanesatedge.*(flipvec==2)];

% Check if there are loners; remove and give a warning.

dtest = diff([0;test;test(length(test))+100]);
ntest = length(dtest);
loners = find(dtest(1:ntest-1)~=0 & dtest(2:ntest)~=0);

if ~isempty(loners)
	disp('WARNING! Some edges had only one plane connected and they will be excluded.')
	disp('         Check if they should be corrected. ')
	disp('         Remember that reflectors must have a plane at the rear side too')
	disp('         Set SHOWTEXT to 4 or higher to get a list of these edges')
	if SHOWTEXT >= 4
		disp('         The edges are (EDB1 corner numbers for the two edge point are given):')
		nloners = length(loners);
		for jjdisp = 1:nloners
			indlon = loners(jjdisp);
			disp(['      ',int2str( floor(test(indlon)/multfac) ),' ',int2str( test(indlon)-floor(test(indlon)/multfac)*multfac )])
			disp(['      CATT: ',int2str( cornernumbers(floor(test(indlon)/multfac) )),' ',int2str( cornernumbers(test(indlon)-floor(test(indlon)/multfac)*multfac ))])
		end
	end

	test(loners) = [];
	planesatedge(loners,:) = [];

end

ntest = length(test);
if ntest >= 2
    if test(ntest) ~= test(ntest-1)
	    test(ntest) = [];
	    planesatedge(ntest,:) = [];
    end
end

% Check if there are triplets

triptest = test(1:2:length(test))-test(2:2:length(test));
triplets = find(triptest~=0);

if ~isempty(triplets)
	disp('ERROR: Triplets: some edges are defined three times which means that two thin planes have a correct')
	disp('       definition but some error on the rear side. You must correct these. ')
	disp('       Only the first one can be detected; it is:')
	disp('       (EDB1 corner numbers for the two edge points are given)')
	[floor(test(triplets(1)*2-1)/multfac) test(triplets(1)*2-1)-floor(test(triplets(1)*2-1)/multfac)*multfac]
    save ~/Documents/Temp/test.mat
	error
end

edgecounter = length(test);

iv1 = find(planesatedge(:,1)~=0);
iv2 = 1:edgecounter;
iv2(iv1) = [];

edgecorners = [floor(test(iv1)/multfac) test(iv1)-floor(test(iv1)/multfac)*multfac];
planesatedge = planesatedge(iv1,:) + planesatedge(iv2,:);

[nedges,slask] = size(edgecorners);

zerosvec3 = zeros(nedges,1);
zerosvec5 = zeros(nedges,3);

if SHOWTEXT >= 4
	disp(['         ',int2str(nedges),' edges found'])
end
thirdpoint = zerosvec5;
for ii = 1:nedges
	refplane = planesatedge(ii,1);
	secondplane = planesatedge(ii,2);
	if secondplane ~= 0
		edgeco1 = corners(edgecorners(ii,1),:);
		edgeco2 = corners(edgecorners(ii,2),:);
		edgemidpoint = edgeco1 + (edgeco2-edgeco1)/2;
		intoplanevec = EDB1cross( (edgeco2-edgeco1).',(planenvecs(secondplane,:)).').';
		inplanepoint = edgemidpoint + intoplanevec*0.1;
		if sum(abs(inplanepoint)) == 0
			inplanepoint = edgemidpoint + intoplanevec*0.2;
		end
		thirdpoint(ii,:) = inplanepoint;
	end
end


%---------------------------------------------------------------
% Calculate the closwedang values for all edges

if SHOWTEXT >= 4
	disp('         wedge angles...')
end
closwedangvec = zerosvec3;

ivec = find( planesatedge(:,1).*planesatedge(:,2) ~= 0);
xwedgevec = [corners(edgecorners(ivec,1),:) corners(edgecorners(ivec,2),:)];
nvec1vec  = planenvecs( planesatedge(ivec,1),: );
xsouvec   = thirdpoint(ivec,:);

for jj = 1:length(ivec)

   ii = ivec(jj);
   [slask,thetas,slask] = EDB1coordtrans1(xsouvec(jj,:),[xwedgevec(jj,1:3);xwedgevec(jj,4:6)],nvec1vec(jj,:));
   if thetas == 0
      closwedangvec(ii) = 0;
   else
      closwedangvec(ii) = 2*pi - thetas;
   end

end

%-------------------------------------------------------------------
% Now we check for duplicates of the edge definitions
% Some edge definitions might need to be swapped: if there are duplicate edges,
% and one of them has closwedangvec = 0
% then there is a mistake in the edge definition, so a swap will be made.

edgecosinglelist = edgecorners(:,1)*multfac + edgecorners(:,2);
nsteps = diff([0;edgecosinglelist]);
ivduplicates = find(nsteps == 0);
nduplicates = length(ivduplicates);
if nduplicates > 0
	if SHOWTEXT >= 4
		disp(['         ',int2str(nduplicates),' edges formed by connected thin plates'])
	end
	for ii = 1:nduplicates
		comb1 = ivduplicates(ii)-1;
		comb2 = comb1+1;
		if closwedangvec(comb1) == 0 | closwedangvec(comb2) == 0,		% edge definitions should be swapped
			plane1 = planesatedge(comb1,2);
			plane2 = planesatedge(comb2,2);
			planesatedge(comb1,2) = plane2;
			planesatedge(comb2,2) = plane1;
			temp1 = thirdpoint(comb1,:);
			temp2 = thirdpoint(comb2,:);
			thirdpoint(comb1,:) = temp2;
			thirdpoint(comb2,:) = temp1;
			if SHOWTEXT >= 4
				disp(['         Swapping an edge definition'])
				disp(['            Planes ',int2str(planesatedge(comb1,1)),' and ',int2str(planesatedge(comb1,2))])
				disp(['            CATT:  ',int2str(planenumbers(planesatedge(comb1,1))),' and ',int2str(planenumbers(planesatedge(comb1,2)))])
				disp(['            Planes ',int2str(planesatedge(comb2,1)),' and ',int2str(planesatedge(comb2,2))])
				disp(['            CATT:  ',int2str(planenumbers(planesatedge(comb2,1))),' and ',int2str(planenumbers(planesatedge(comb2,2)))])
			end

   			xwedge = [corners(edgecorners(comb1,1),:);corners(edgecorners(comb1,2),:)];
   			nvec1   = planenvecs( planesatedge(comb1,1),: );
   			xsou   = thirdpoint(comb1,:);
   			[slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
   			if thetas == 0
    	  		closwedangvec(comb1) = 0;
   			else
    	  		closwedangvec(comb1) = 2*pi - thetas;
   			end

   			xwedge = [corners(edgecorners(comb2,1),:);corners(edgecorners(comb2,2),:)];
   			nvec1   = planenvecs( planesatedge(comb2,1),: );
   			xsou   = thirdpoint(comb2,:);
   			[slask,thetas,slask] = EDB1coordtrans1(xsou,xwedge,nvec1);
   			if thetas == 0
    	  		closwedangvec(comb2) = 0;
   			else
    	  		closwedangvec(comb2) = 2*pi - thetas;
   			end

		end   % if closwedangvec(comb1) == 0 |

	end		% for ii = 1:nduplicates
end

if nplanes < 256
    planesatedge = uint8(planesatedge);
elseif nplanes < 65535
    planesatedge = uint16(planesatedge);
else
    planesatedge = uint32(planesatedge);
end

%-------------------------------------------------------------------
% Find which edges each plane is connected to

if SHOWTEXT >= 4
	disp('         find edgesatplane')
end

edgesatplane = zeros(nplanes,double(max(ncornersperplanevec)));

for ii= 1:nplanes
	iv = find(planesatedge==ii);
	iv = iv - floor((iv-1)/nedges)*nedges;
	if ~isempty(iv)
		edgesatplane(ii,1:length(iv)) = iv.';
	end
end

% Check how many edges are defined for each plane. There must be at
% least three edges per plane.

tempnedgesperplanevec = sum(edgesatplane.'~=0).';

iv = find(tempnedgesperplanevec<3);

if ~isempty(iv)

    disp(' ')
    if SHOWTEXT >= 4
        for ii = 1:length(iv)
            disp(['         Plane ',int2str(iv(ii)),' has only ',int2str(tempnedgesperplanevec(iv(ii))),' edges connected to it'])
        end
        error(['ERROR: The planes above have less than three edges.'])
    else
        error(['ERROR: Some planes have less than three edges. Set SHOWTEXT >= 4 to see a list of those planes.'])
    end

end


%-------------------------------------------------------------------
% Make a special list of thin planes

planeisthin = zerosvec1;
rearsideplane = zerosvec1;

iv = find(closwedangvec==0);
if ~isempty(iv)
	planenumberlist = planesatedge(iv,1);
	planeisthin(planenumberlist) = planeisthin(planenumberlist) + 1;
	rearsideplane(planenumberlist) = planesatedge(iv,2);
end

iv = find(closwedangvec == 0 & planesatedge(:,2)~=0);
if ~isempty(iv)
	planenumberlist = planesatedge(iv,2);
	planeisthin(planenumberlist)   = planeisthin(planenumberlist) + 1;
	rearsideplane(planenumberlist) = planesatedge(iv,1);
end

planeisthin = sign(planeisthin);
listofthinplanes = find(planeisthin);

%---------------------------------------------------------------
% Closwedangvec should be calculated into nyvec

nyvec = pi./(2*pi - closwedangvec);
integerny = ( abs(nyvec - round(nyvec)) < 1e-10 );

%-----------------------------------------------------------
% Construct some other useful edge variables

if difforder >= 1
	edgestartcoords      = zerosvec5;
	edgestartcoordsnudge = zerosvec5;
	edgeendcoords        = zerosvec5;
	edgeendcoordsnudge   = zerosvec5;
	edgemidcoords        = zerosvec5;
	edgenvecs            = zerosvec5;

    edgestartcoords = [corners(edgecorners(:,1),:)];
	edgeendcoords   = [corners(edgecorners(:,2),:)];
	edgemidcoords   = (edgestartcoords+edgeendcoords)/2;
	edgenvecs       = planenvecs(planesatedge(:,1),:);

    edgenormvecs = edgeendcoords - edgestartcoords;

	edgestartcoordsnudge = edgestartcoords + 1e-10*edgenormvecs;
	edgeendcoordsnudge   = edgeendcoords   - 1e-10*edgenormvecs;

    edgelengthvec = sqrt(sum( ((edgenormvecs).^2).' )).';
	edgenormvecs = edgenormvecs./edgelengthvec(:,ones(1,3));

else
    edgestartcoords = [];
    edgeendcoords = [];
    edgenvecs = [];
    edgelengthvec = [];
end


%---------------------------------------------------------------
%
%		BACK TO SOME PURE PLANE RELATED PARAMETERS
%
%---------------------------------------------------------------
%
% First, make a list of which planes that are absorbing/non-reflective

planeisabsorbing = zerosvec1;

listofactiveplanes = find(reflfactors~=0);
listofabsorbingplanes = find(reflfactors == 0);
planeisabsorbing(listofabsorbingplanes) = ones(size(listofabsorbingplanes));

%---------------------------------------------------------------
% Help variables: planealignedwplane and planeconnectstoplane
%
% Preparatory for the check of which planes see which plane: check
% if any planes are aligned with each other. If two planes are aligned
% then planealignedwplane = 1 for that combination.
% However, if two planes are back-to back (thin planes), then
% planealignedwplane = 2 for that combination.
%
% Also make a matrix of which planes connect to which planes and via
% which edges.

planealignedwplane = zerosvec2;

% First, check which planes are aligned with each other

for ii = 1:nplanes-1
	if planeisabsorbing(ii) ~= 1
		oneplaneeq = planeeqs(ii,:);
		diffvec1 = oneplaneeq(ones(nplanes-ii,1),:) - planeeqs(ii+1:nplanes,:);
		diffvec2 = oneplaneeq(ones(nplanes-ii,1),:) + planeeqs(ii+1:nplanes,:);
		diffvec = min( [sum( diffvec1.'.^2 ).'   sum( diffvec2.'.^2 ).'].' ).';
		alignedplanes = find(diffvec < geomacc) + ii;
		if ~isempty(alignedplanes)
			planealignedwplane(alignedplanes,ii) = int8(double(planealignedwplane(alignedplanes,ii)) + double((planeisabsorbing(alignedplanes)~=1)));
		end
	end
end

planealignedwplane = (sparse(sign(double(planealignedwplane) + double(planealignedwplane).')));

% Second, check which planes are back to back

if ~isempty(listofthinplanes)
	for ii = 1:length(listofthinplanes)
		plane = listofthinplanes(ii);
		rearplane = rearsideplane(plane);
		if planeisabsorbing(plane) ~= 1 & planeisabsorbing(rearplane) ~= 1
			planealignedwplane(plane,rearplane) = 2;
			planealignedwplane(rearplane,plane) = 2;
		end
	end
end

planeconnectstoplane = zerosvec2;
clear zerosvec2

for ii = 1:nplanes
	if planeisabsorbing(ii) ~= 1
		edgelist = edgesatplane(ii,:);
 		edgelist = edgesatplane(ii,1:double(ncornersperplanevec(ii)));
		ivec = planesatedge(edgelist,:);
		ivec = reshape(ivec.',length(edgelist)*2,1);
		ivec = ivec(find(ivec~=ii));
		okplanes = find(planeisabsorbing(ivec)~=1);
		ivec = ivec(okplanes);
		if ~isempty(ivec)
			planeconnectstoplane(ii,ivec) = edgelist(okplanes);
		end
	end
end

%---------------------------------------------------------------
%
%			planeseesplane
%
%---------------------------------------------------------------
% Check which planes see which planes:
% 1. If one of the planes has reflfactors = 0, then the two can't see each other.
% 2. Reflective planes that are aligned with each other can not reach each other
% 3. Reflective planes that are not aligned but have the same normal vectors
%    can not reach each other
% 4. Planes that have all its corners behind another plane can not be seen by that plane.
%
% planeseesplane = 0 means that either, one of the planes is non-reflective or, that
%					 two planes never can see each other, but they are not aligned with each other.
% planeseesplane = 1 means that two reflective planes might see each other, but there could
%                    be obstructions
% planeseesplane = -1 means that two reflective planes are aligned (but they are not back-to-back)
%					  and thus they can never see each other.
% planeseesplane = -2 means that two reflective planes are back-to-back with each other and thus
%					  they can never see each other.

if SHOWTEXT >= 4
	disp('         Check which planes see which planes')
end

% First an auxilliary matrix which can be used for edgeseesplane later:
% cornerinfrontofplane, size [nplanes,ncorners]
% Values are the same as given by EDB1infrontofplane:
%   1 means that a point is in front of the plane
%   0 means that a point is aligned with a plane
%   -1 means that a point is behind a plane
% All corners that belong to a plane will have the value 0.

% Corner number is given by the col no.
if ncorners < 256
    cornernumb = uint8([1:ncorners]);
elseif ncorners < 65536
    cornernumb = uint16([1:ncorners]);
else
    cornernumb = uint32([1:ncorners]);
end

% Plane numbers is given by the row no.
if nplanes < 256
    planenumb = uint8([1:nplanes].');
elseif nplanes < 65536
    planenumb = uint16([1:nplanes].');
else
    planenumb = uint32([1:nplanes].');
end

cornerinfrontofplane = int8(EDB1infrontofplane(corners,planenvecs,planecorners,[],cornernumb,planenumb));
clear planenumb cornernumb
cornerinfrontofplane = reshape(cornerinfrontofplane,nplanes,ncorners);

planeseesplane = int8(ones(nplanes,nplanes));

% 1. If one of the planes has reflfactors = 0, then the two can't see each other.

if ~isempty(listofabsorbingplanes)
	zerosvec = zeros(length(listofabsorbingplanes),nplanes);
	planeseesplane(listofabsorbingplanes,:) = zerosvec;
	planeseesplane(:,listofabsorbingplanes) = zerosvec.';
end

% 2. Reflective planes that are aligned with each other can not reach each other

ivec = find( planealignedwplane == 1);
if ~isempty(ivec)
	planeseesplane(ivec) = - 1*ones(size(ivec));
end

ivec = find( planealignedwplane == 2);
if ~isempty(ivec)
	planeseesplane(ivec) = - 2*ones(size(ivec));
end

% 3. Reflective planes that have the same normal vectors can not reach each other

numvec = [1:nplanes];
for ii = 1:nplanes
	if planeisabsorbing(ii) ~= 1
		nvec1 = planenvecs(ii,:);
		ivec = find(planealignedwplane(ii,:)==0 & planeisabsorbing.'==0 & numvec>ii);
		if ~isempty(ivec)
			similarvec = abs( nvec1(ones(length(ivec),1),:) - planenvecs(ivec,:)) < geomacc;
			similarvec = prod(  double(similarvec.') ).';
			ivec2 = ivec(find(similarvec==1));
			if ~isempty(ivec2)
				zerosvec = zeros(size(ivec2));
				planeseesplane(ii,ivec2) = zerosvec;
				planeseesplane(ivec2,ii) = zerosvec.';
			end
		end
	end
end

% 3.5 A plane does not see itself

iv = [0:nplanes-1]*nplanes + [1:nplanes];
planeseesplane(iv) = zeros(size(iv));

% 4. Planes that have all its corners behind another plane can not be seen by that plane.
%    First, we construct a list referring to the complete matrix of size [nplanes,nplanes]
%    The index vector is called iv. Then we find which combinations that
%    could be cleared. After having cleared a number of combinations (in
%    iv, fromplanenumb and toplanenumb) we pick out out only the non-zero
%    indices in iv.

if nplanes*nplanes < 128
    iv = int8([1:nplanes*nplanes].');
elseif nplanes*nplanes < 32768
    iv = int16([1:nplanes*nplanes].');
else
    iv = int32([1:nplanes*nplanes].');
end
iv = reshape(iv,nplanes,nplanes);

% If there are any absorbing planes, we zero all plane-to-plane
% combinations that involve such a plane

if ~isempty(listofabsorbingplanes)
    nabsorbingplanes = length(listofabsorbingplanes);
    iv(listofabsorbingplanes,:) = zeros(nabsorbingplanes,nplanes);
    iv(:,listofabsorbingplanes) = zeros(nplanes,nabsorbingplanes);
end

% Check connecting planes first: if two planes are connected and their
% shared edge has a closwedangvec > pi, then the two planes can see each
% other. So, we can zero the plane-pair combinations where closwedangvec
% <= pi but larger than zero.

edgecombstozero = find(closwedangvec<=pi & closwedangvec > 0);
if ~isempty(edgecombstozero)
    ivzerocombs = (double(planesatedge(edgecombstozero,1))-1)*nplanes + double(planesatedge(edgecombstozero,2));
    iv(ivzerocombs) = zeros(size(ivzerocombs));
    ivzerocombs = (double(planesatedge(edgecombstozero,2))-1)*nplanes + double(planesatedge(edgecombstozero,1));
    iv(ivzerocombs) = zeros(size(ivzerocombs));
end

% Now we should keep only the index numbers in iv that are still non-zero
% We also take the chance to zero planeseesplane for the combinations that
% could be cleared.
ivclear = uint32(find(iv==0));
planeseesplane(ivclear) = zeros(size(ivclear));
iv(ivclear) = [];

% Of the remaining plane-to-plane combinations, we pick out the ones
% for which we need to check if the "to-plane" is in front of the "from-plane"
iv = iv(find(planeseesplane(iv)==1 & planeconnectstoplane(iv)==0));

% We create full matrices, fromplanenumb and toplanenumb, and then keep
% only the combinations that remain in iv
% to keep fromplanenumb and toplanenumb aligned with iv.

% To-plane numbers is given by the row no.
if nplanes < 256
    toplanenumb = uint8([1:nplanes].');
elseif nplanes < 65536
    toplanenumb = uint16([1:nplanes].');
else
    toplanenumb = uint32([1:nplanes].');
end
toplanenumb = reshape(toplanenumb(:,ones(1,nplanes)),nplanes*nplanes,1);
toplanenumb = toplanenumb(iv);

% From-plane numbers is given by the col no.
if nplanes < 256
    fromplanenumb = uint8([1:nplanes]);
elseif nplanes < 65536
    fromplanenumb = uint16([1:nplanes]);
else
    fromplanenumb = uint32([1:nplanes]);
end
fromplanenumb = reshape(fromplanenumb(ones(nplanes,1),:),nplanes*nplanes,1);
fromplanenumb = fromplanenumb(iv);

% Now, if a plane has *all* its corners behind another plane, those two
% planes can not see each other.
% We check the corners of the "to-plane" (columns) and check if they are in front of
% the "from-plane" (rows).
% The ivreftocoplamatrix is the index number to the
% cornerinfrontofplanematrix, which has the size [nplanes,ncorners].
% NB! In order to save some space, we don't construct the ivreftocoplamatrix specifically.

% First we check the first 3 corners since all planes have at least 3
% corners.

% Plane corner 1
cornersbehind = cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,1 ))-1 )*nplanes)<1;

% Plane corner 2
cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,2 ))-1 )*nplanes) < 1);

% Plane corner 3
cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,3 ))-1 )*nplanes) < 1);

if minncornersperplane == 4 & maxncornersperplane == 4
%   Plane corner 4
    cornersbehind = cornersbehind &(cornerinfrontofplane(double(fromplanenumb) + ( double(planecorners( toplanenumb,4 ))-1 )*nplanes) < 1);
elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
    ivall = uint32([1:length(fromplanenumb)].');
    if minncornersperplane == 3
        iv3cornerplanes = find(ncornersperplanevec(toplanenumb)==3);
        ivall(iv3cornerplanes) = [];
    end
    if maxncornersperplane == 4
% Plane corner 4
        cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
    else
        ivmorethan4 = find(ncornersperplanevec(toplanenumb(ivall))>4);
        temp = ivall(ivmorethan4);
        ivall(ivmorethan4) = [];
        ivmorethan4 = temp;
        cornersbehind(ivall) = cornersbehind(ivall) &(cornerinfrontofplane(double(fromplanenumb(ivall)) + ( double(planecorners( toplanenumb(ivall),4 ))-1 )*nplanes) < 1);
        for ii = 5:maxncornersperplane
            ivselection = find(ncornersperplanevec(toplanenumb(ivmorethan4))>=ii);
            cornersbehind(ivmorethan4(ivselection)) = cornersbehind(ivmorethan4(ivselection)) &(cornerinfrontofplane(double(fromplanenumb(ivmorethan4(ivselection))) + ( double(planecorners( toplanenumb(ivmorethan4(ivselection)),ii ))-1 )*nplanes) < 1);
        end
        clear ivselection
    end
end
clear toplanenumb fromplanenumb

ivclear = iv(find(cornersbehind==1));
clear iv

planeseesplane(ivclear) = zeros(size(ivclear));
clear ivclear

temp1 = (planeseesplane~=0);
temp2 = (planeseesplane.'~=0);
planeseesplane = int8(double(planeseesplane).*(temp1.*temp2));

% New addition:
% If the user has asked for it, check plane-mid-point to plane-mid-point
% for obstruction. If there are obstructions, set planeseesplane to 0 for
% that combination.

if planeseesplanestrategy == 1
	ivorig = uint32(find(planeseesplane==1));
    iv = ivorig;
	fromplane = ceil(double(iv)/nplanes);
	toplane = uint16(double(iv) - (fromplane-1)*nplanes);
    fromplane = uint16(fromplane);
	listofplanesthatcanobscure = find( sum(planeseesplane==1) );

    iv4planes = find(ncornersperplanevec == 4);
    planemidpoints = zeros(nplanes,3);

    if any(ncornersperplanevec~=4)
        disp('WARNING: planeseesplanestrategy 1 not implemented for models with non-4-corner planes')
        planeseesplanestrategy = 0;
    else
        xvalues = corners(planecorners.',1);
        xvalues = reshape(xvalues,4,length(xvalues)/4);
        yvalues = corners(planecorners.',2);
        yvalues = reshape(yvalues,4,length(yvalues)/4);
        zvalues = corners(planecorners.',3);
        zvalues = reshape(zvalues,4,length(zvalues)/4);
        planemidpoints = [mean(xvalues).' mean(yvalues).' mean(zvalues).'];
    end
end

if planeseesplanestrategy == 1
	for ii = listofplanesthatcanobscure
        planeobsc = ii;
        ivcheckvis1 = uint32( planeobsc + (double(fromplane)-1)*nplanes);
        ivcheckvis2 = uint32( planeobsc + (double(toplane)-1)*nplanes);
        iv3 = find(toplane~=planeobsc & fromplane~=planeobsc & ((planeseesplane(ivcheckvis1) == 1) | (planeseesplane(ivcheckvis2) == 1)) & (planeseesplane(ivcheckvis1) >= 0) & (planeseesplane(ivcheckvis2) >= 0) );
        tempcanplaneobstruct = zerosvec1.';
        tempcanplaneobstruct(planeobsc) = 1;
        [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstrpaths(planemidpoints(fromplane(iv3),:),planemidpoints(toplane(iv3),:),[],[],tempcanplaneobstruct,planeseesplane,...
            planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
        if length(iv3) > length(nonobstructedpaths)
            iv3(nonobstructedpaths) = [];
            iv(iv3) = [];
	        toplane(iv3) = [];
	        fromplane(iv3) = [];
        end
	end
    planeseesplane(ivorig) = 0;
    planeseesplane(iv) = 1;
    clear iv ivorig toplane fromplane ivcheckvis ivcheckvis2
end

%---------------------------------------------------------------
%
%		EDGE RELATED PARAMETERS AGAIN
%
%---------------------------------------------------------------
%
% Which edges should be switched off?
% (1) Go through the list edgecorners. If any edge contains one of the
%     corners, whose numbers are in switchoffcorners, then that edge should be
%     stored in the list offedges.
% (2) If an edge has an integer ny-value, then it should be switched off.
% (3) If one of the two connecting planes has reflfactors = 0, then it should
%     be switched off.

% disp('   switching off possibly unwanted edges...')

offedges = zerosvec3;

% (1) Go through the list edgecorners. If any edge contains one of the
%     corners, whose numbers are in switchoffcorners, then that edge should be
%     stored in the list offedges.

if ~isempty(switchoffcorners)
	for ii = 1:nedges
	   corner1 = edgecorners(ii,1);
	   corner2 = edgecorners(ii,2);
	   remove = sum( corner1 == switchoffcorners ) + sum( corner2 == switchoffcorners );
	   offedges(ii) = sign(remove);
	end
end

% (2) If an edge has an integer ny-value, then it should be switched off.

offedges = offedges + integerny;

% (3) If one of the two connecting planes has reflfactors = 0, then it should
%     be switched off.

reflectingedges = reshape(reflfactors(planesatedge),nedges,2);
reflectingedges = reflectingedges(:,1).*reflectingedges(:,2);
offedges = offedges + (1-reflectingedges);

offedges = find(offedges ~= 0);

%--------------------------------------------------------------------------------
%
%			edgeseesplane
%
%--------------------------------------------------------------------------------
%
% Now check which planes each edge can see.
%    1. Set all combinations with reflfactors = 0 to -3
%    2. Switch off all combinations with offedges
% 	 3. For all edges that are aligned with a plane, set edgeseesplane to:
%			-1 if the edge doesn't belong to the plane
%			-2 if the edge does belong to the plane
%    4. If at least one edge end point is in front of a plane, then the
%       plane is potentially visible (edgeseesplane = 1)
%    5. If at least one corner point is in front of:
%           *one* of the two edge planes, when closwedang < pi
%           *both* edge planes, when closwedangvec > pi

if SHOWTEXT >= 4
	disp('         Check which edges see which planes')
end

edgeseesplane = int8(ones(nplanes,nedges));

% % % 20120414 PS: Why has this first stage been shut off? CHECK
%    1. Set all edges belonging to planes with reflfactors = 0 to
%    edgeseesplane = -3.

% % if ~isempty(listofabsorbingplanes)
% % 	edgeseesplane(listofabsorbingplanes,:) = -3*ones(length(listofabsorbingplanes),nedges);
% % end

%    2. Switch off all combinations with offedges

if ~isempty(offedges)
	edgeseesplane(:,offedges) = zeros(nplanes,length(offedges));
end

% 	3.  For all edges that belong to a plane, set edgeseesplane to -2
%    	Also, for all other edges that are aligned with a plane, set
%    	edgeseesplane to -1

edgelist = [1:nedges].';
edgelist(offedges) = [];

plane1list = planesatedge(edgelist,1);
plane2list = planesatedge(edgelist,2);
indexvec = uint32( double(plane1list) + (edgelist-1)*nplanes);
edgeseesplane(indexvec) = -2*(reflfactors(plane1list)~=0);
indexvec = uint32( double(plane2list) + (edgelist-1)*nplanes);
edgeseesplane(indexvec) = -2*(reflfactors(plane2list)~=0);

for ii = 1:length(edgelist)

	aligners1 = find(planealignedwplane(:,plane1list(ii))==1);
	if ~isempty(aligners1)
		edgeseesplane(aligners1,edgelist(ii)) = -1*ones(size(aligners1));
	end

	aligners2 = find(planealignedwplane(:,plane2list(ii))==1);
	if ~isempty(aligners2)
		edgeseesplane(aligners2,edgelist(ii)) = -1*ones(size(aligners2));
	end

end

%    4. If at least one edge end point is in front of a plane, then the
%       plane is potentially visible.
%       Do this for all plane-edge combinations for which edgeseesplane = 1
%       up til now.
%    5. If at least one corner point is in front of:
%           *one* of the two edge planes, when closwedang < pi
%           *both* edge planes, when closwedangvec > pi

ntot = nplanes*nedges;
ivclear = find(edgeseesplane~=1);
% Edge number is given by the col no.
if nedges < 256
    edgenumb = uint8([1:nedges]);
elseif nedges < 65536
    edgenumb = uint16([1:nedges]);
else
    edgenumb = uint32([1:nedges]);
end
edgenumb = reshape(edgenumb(ones(nplanes,1),:),nplanes*nedges,1);
edgenumb(ivclear) = [];
% Plane numbers is given by the row no.
if nplanes < 256
    planenumb = uint8([1:nplanes].');
elseif nplanes < 65536
    planenumb = uint16([1:nplanes].');
else
    planenumb = uint32([1:nplanes].');
end
planenumb = reshape(planenumb(:,ones(1,nedges)),nplanes*nedges,1);
planenumb(ivclear) = [];

if nplanes*nedges < 128
    iv = int8([1:nplanes*nedges].');
elseif nplanes*nedges < 32768
    iv = int16([1:nplanes*nedges].');
else
    iv = int32([1:nplanes*nedges].');
end
iv(ivclear) = [];

if ~isempty(edgenumb)
	% For the "edgecorners in front of plane"-test
	% we can use the cornerinfrontofplane matrix, but we need to construct
	% index vectors, based on the planenumb and edgenumb values.

	edgeinfrontofplane =                      cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,1) )-1 )*nplanes );
	edgeinfrontofplane = (edgeinfrontofplane==1) | (cornerinfrontofplane( double(planenumb) + ( double( edgecorners(edgenumb,2) )-1 )*nplanes )==1);
	edgeseesplane(iv) = int8(double(edgeseesplane(iv)).*double(edgeinfrontofplane));
	clear edgeinfrontofplane
end

if ~isempty(edgenumb)
	% For the "planecorners in front of edge-plane"-test
	% we can use the cornerinfrontofplane matrix, but we need to construct
	% index vectors, based on the planenumb and edgenumb values.
	%
	% First we split up the iv into two halves: the one with edges < pi and the
	% one with edges > pi;

	ivsmaller = uint32(find(closwedangvec(edgenumb)<pi));
	ivlarger = uint32([1:length(edgenumb)].');
	ivlarger(ivsmaller) = [];

    if ~isempty(ivsmaller)
		% First the edges smaller than pi (closwedang < pi): at least one plane corner needs
		% to be in front of one or the other plane.
		% Plane corner 1
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),1)  ) + ( double(planecorners( planenumb(ivsmaller),1 ))-1 )*nplanes;
		planecornerinfrontofedge =                            cornerinfrontofplane(ivreftocoplamatrix)==1;
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),2)  ) + ( double(planecorners( planenumb(ivsmaller),1 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
		% Plane corner 2
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),1)  ) + ( double(planecorners( planenumb(ivsmaller),2 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),2)  ) + ( double(planecorners( planenumb(ivsmaller),2 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
		% Plane corner 3
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),1)  ) + ( double(planecorners( planenumb(ivsmaller),3 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),2)  ) + ( double(planecorners( planenumb(ivsmaller),3 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);

		if minncornersperplane == 4 & maxncornersperplane == 4
            % Plane corner 4
			ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),1)  ) + ( double(planecorners( planenumb(ivsmaller),4 ))-1 )*nplanes;
			planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
			ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller),2)  ) + ( double(planecorners( planenumb(ivsmaller),4 ))-1 )*nplanes;
			planecornerinfrontofedge = planecornerinfrontofedge | (cornerinfrontofplane(ivreftocoplamatrix)==1);
		elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
            ivall = uint32(1:length(ivsmaller));
            if minncornersperplane == 3
                iv3cornerplanes = find(ncornersperplanevec(planenumb(ivsmaller))==3);
                ivall(iv3cornerplanes) = [];
            end
            if maxncornersperplane == 4
                % Plane corner 4
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivall)),1)  ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivall)),2)  ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
            else
                ivmorethan4 = find(ncornersperplanevec(planenumb(ivsmaller(ivall)))>4);
                temp = ivall(ivmorethan4);
                ivall(ivmorethan4) = [];
                ivmorethan4 = temp;
                % Plane corner 4
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivall)),1)  ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivall)),2)  ) + ( double(planecorners( planenumb(ivsmaller(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
                for ii = 5:maxncornersperplane
                    % Plane corner 5,6,7,etc
                    ivselection = find(ncornersperplanevec(planenumb(ivsmaller(ivmorethan4)))>=ii);
					ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivmorethan4)),1)  ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
					planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
					ivreftocoplamatrix = double(  planesatedge(edgenumb(ivsmaller(ivmorethan4)),2)  ) + ( double(planecorners( planenumb(ivsmaller(ivmorethan4)),4 ))-1 )*nplanes;
					planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (cornerinfrontofplane(ivreftocoplamatrix)==1);
                end
                clear ivselection
            end
		end
		clear ivreftocoplamatrix

		ivclear = ivsmaller(find(planecornerinfrontofedge==0));
		edgeseesplane(iv(ivclear)) = 0;
    end

    if ~isempty(ivlarger)
		% Second the edges larger than pi (closwedang > pi): at least one plane corner needs
		% to be in front of *both* edge planes.
		% Plane corner 1
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),1)  ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
		planecornerinfrontofedge =                            cornerinfrontofplane(ivreftocoplamatrix)==1;
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),2)  ) + ( double(planecorners( planenumb(ivlarger),1 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge & (cornerinfrontofplane(ivreftocoplamatrix)==1);
		% Plane corner 2
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),1)  ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
		temp =                                                cornerinfrontofplane(ivreftocoplamatrix)==1;
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),2)  ) + ( double(planecorners( planenumb(ivlarger),2 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
		% Plane corner 3
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),1)  ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
		temp =                                                cornerinfrontofplane(ivreftocoplamatrix)==1;
		ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),2)  ) + ( double(planecorners( planenumb(ivlarger),3 ))-1 )*nplanes;
		planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));

		if minncornersperplane == 4 & maxncornersperplane == 4
            % Plane corner 4
			ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),1)  ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
			temp =                                                cornerinfrontofplane(ivreftocoplamatrix)==1;
			ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger),2)  ) + ( double(planecorners( planenumb(ivlarger),4 ))-1 )*nplanes;
			planecornerinfrontofedge = planecornerinfrontofedge | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
		elseif not(minncornersperplane == 3 & maxncornersperplane == 3)
            ivall = uint32(1:length(ivlarger));
            if minncornersperplane == 3
                iv3cornerplanes = find(ncornersperplanevec(planenumb(ivlarger))==3);
                ivall(iv3cornerplanes) = [];
            end
            if maxncornersperplane == 4
                % Plane corner 4
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivall)),1)  ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
				temp =                                                              cornerinfrontofplane(ivreftocoplamatrix)==1;
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivall)),2)  ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
            else
                ivmorethan4 = find(ncornersperplanevec(planenumb(ivlarger(ivall)))>4);
                temp = ivall(ivmorethan4);
                ivall(ivmorethan4) = [];
                ivmorethan4 = temp;
              % Plane corner 4
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivall)),1)  ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
				temp =                                                              cornerinfrontofplane(ivreftocoplamatrix)==1;
				ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivall)),2)  ) + ( double(planecorners( planenumb(ivlarger(ivall)),4 ))-1 )*nplanes;
				planecornerinfrontofedge(ivall) = planecornerinfrontofedge(ivall) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
                for ii = 5:maxncornersperplane
                    % Plane corner 5,6,7,etc
                    ivselection = find(ncornersperplanevec(planenumb(ivlarger(ivmorethan4)))>=ii);
					ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivmorethan4)),1)  ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
					temp =                                                          cornerinfrontofplane(ivreftocoplamatrix)==1;
					ivreftocoplamatrix = double(  planesatedge(edgenumb(ivlarger(ivmorethan4)),2)  ) + ( double(planecorners( planenumb(ivlarger(ivmorethan4)),4 ))-1 )*nplanes;
					planecornerinfrontofedge(ivmorethan4) = planecornerinfrontofedge(ivmorethan4) | (temp & (cornerinfrontofplane(ivreftocoplamatrix)==1));
                end
                clear ivselection
            end
		end

		ivclear = ivlarger(find(planecornerinfrontofedge==0));
		edgeseesplane(iv(ivclear)) = 0;
    end
end

clear ivlarger ivsmaller ivreftocoplamatrix edgenumb planenumb iv

%----------------------------------------------------------------------------
%
%           indentingedgepairs
%
%--------------------------------------------------------------------------------
%
% If there are any planes with indents, make a list of all consecutive edge pairs that
% form an indent, because such edge pairs could never see each other.

indentingedgepairs = [];

if any(planehasindents)
    iv = find(indentingcorners);
    indentingedgepairs = zeros(length(iv),2);
    indentingplanenumbers = mod(iv,nplanes);
    ivzeros = find(indentingplanenumbers==0);
    if ~isempty(ivzeros)
       indentingplanenumbers(ivzeros) = nplanes;
    end
    conumbers = (iv - indentingplanenumbers)/nplanes+1;

    % We expand planecorners temporarily, cyclically
    planecorners = [planecorners planecorners(:,1:2)];

    for ii = 1:length(iv)
        edge1cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii):conumbers(ii)+1));
        edge2cornerpair = sort(planecorners(indentingplanenumbers(ii),conumbers(ii)+1:conumbers(ii)+2));

        [slask,edge1] = ismember(edge1cornerpair,edgecorners,'rows');
        [slask,edge2] = ismember(edge2cornerpair,edgecorners,'rows');

        edgepair = sort([edge1 edge2]);
        indentingedgepairs(ii,:) = edgepair;

    end

    % Remove the temporary expansion

    planecorners = planecorners(:,1:size(planecorners,2)-2);
end

%----------------------------------------------------------------------------
%
%           canplaneobstruct
%
%--------------------------------------------------------------------------------
%
% Check which planes that are potentially obstructing:
%
% If **all** closwedangvec are < pi and we have an exterior problem =>
%   the scatterer is without indents.
%
% If **all** closwedangvec are > pi and we have an interior problem =>
%   the room is without any indents so there can be no obstructions.
%
% For exterior problems, all active planes can potentially obstruct.
%
% If these simple checks give a results that conflict with the input
% parameters open_or_closed_model or int_or_ext_model, give an
% error message.

canplaneobstruct = [];
if sum(closwedangvec<pi) == 0
    if SHOWTEXT >= 2
        disp(' ')
        disp('   The model is an interior problem, and the room has no indents')
        disp('   so no obstruction tests will be necessary.')
    end
    if int_or_ext_model == 'e'
        disp( 'ERROR: In the setup file this model was defined as an exterior problem, but it is clearly an interior problem.')
        error('       Please check the orientation of the planes.')
    end
    canplaneobstruct = zeros(1,nplanes);
elseif sum(closwedangvec>0) == 0
    if SHOWTEXT >= 2
        disp(' ')
        disp('   The model is an exterior problem (scattering)')
        disp('   with only thin planes.')
    end
    if int_or_ext_model == 'i'
        disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
        error('       Please check the orientation of the planes.')
    end
    canplaneobstruct = ones(1,nplanes);
elseif sum(closwedangvec>pi) == 0
    if SHOWTEXT >= 2
        disp(' ')
        disp('   The model is an exterior problem (scattering)')
        disp('   and the scatterer has no indents.')
    end
    if int_or_ext_model == 'i'
        disp( 'ERROR: In the setup file this model was defined as an interior problem, but it is clearly an exterior problem.')
        error('       Please check the orientation of the planes.')
    end
    canplaneobstruct = ones(1,nplanes);
end

% For an interior problem we know for sure that if a plane
% has all other points in front of itself, it can not obstruct.
% NB! We need to check only points that do not belong to planes
% that are aligned with the plane, or planes that are non-active.

if int_or_ext_model == 'e'
    canplaneobstruct = (reflfactors.'~=0);
elseif isempty(canplaneobstruct)

    canplaneobstruct = (reflfactors.'~=0);
    maskvec1 = canplaneobstruct.';
    listofactiveplanes = find(canplaneobstruct);
    zerosvec = zeros(nplanes,1);

    for ii = 1:length(listofactiveplanes)
        iplane = listofactiveplanes(ii);
        maskvec2 = zerosvec;
        maskvec2(iplane) = 1;
        otherplanestocheck = find((double(planeseesplane(:,iplane))+maskvec2)~=1 & planeseesplane(:,iplane)~=-1 & maskvec1~=0);
        if ~isempty(otherplanestocheck),
            cornerstocheck = unique(planecorners(otherplanestocheck,:));
            cornerstocheck = setdiff(unique(cornerstocheck),planecorners(iplane,1:ncornersperplanevec(iplane)));
            pointinfront = EDB1infrontofplane(corners(cornerstocheck,:),planenvecs(iplane,:),corners(planecorners(iplane,1),:),corners(planecorners(iplane,2),:));
            if isempty(find(pointinfront==-1))
                canplaneobstruct(iplane) = 0;
            end
        else
            canplaneobstruct(iplane) = 0;
        end

    end
end

%---------------------------------------------------------------
% For each thin plane, make sure that all edges have the same normal
% vector. If not, switch the direction of that edge, and all the other
% relevant parameters.

if difforder >= 1
	ivthin = find(planeisthin);

	for ii = 1:length(ivthin)
        plane = ivthin(ii);
        if rearsideplane(plane) > plane
            edgelist = unique(edgesatplane(plane,1:nedgesperplanevec(plane)));
            iv = find(closwedangvec(edgelist,:)==0);
            if ~isempty(iv)
                edgelist = edgelist(iv);
                nedgestemp = length(edgelist);
                edgenveclist = edgenvecs(edgelist,:);
                refnvec = edgenveclist(1,:);
                nvecdiff = edgenveclist(2:nedgestemp,:) - refnvec(ones(nedgestemp-1,1),:);
                nvecdiff = sum(abs(nvecdiff.')).';
                ivswitch = find(nvecdiff)+1;
                if ~isempty(ivswitch)
                    edgenumber = edgelist(ivswitch);
                    edgecorners(edgenumber,:)  = [edgecorners(edgenumber,2) edgecorners(edgenumber,1)];
                    planesatedge(edgenumber,:) = [planesatedge(edgenumber,2) planesatedge(edgenumber,1)];
                    edgenvecs(edgenumber,:) = -edgenvecs(edgenumber,:);
                    edgestartcoords(edgenumber,:) = corners(edgecorners(edgenumber,1),:);
                    edgeendcoords(edgenumber,:)   = corners(edgecorners(edgenumber,2),:);
                    tempvec = edgeendcoordsnudge(edgenumber,:);
                    edgeendcoordsnudge(edgenumber,:) = edgestartcoordsnudge(edgenumber,:);
                    edgestartcoordsnudge(edgenumber,:) = tempvec;
                end
            end
        end

	end
end

%----------------------------------------------------------------------------
%
%		SAVE THE VARIABLES
%
%----------------------------------------------------------------------------
%
% Also the variables from cadgeo??
% Yes: planeeqs planenvecs ncornersperplanevec planecorners corners
%      minvals maxvals

edgeseesplane = int8(edgeseesplane);
planeseesplane = int8(planeseesplane);

if ncorners < 256
    planecorners = uint8(planecorners);
elseif ncorners < 65536
    planecorners = uint16(planecorners);
end
planeisthin = uint8(planeisthin);

if max(ncornersperplanevec) <= 255
    ncornersperplanevec = uint8(ncornersperplanevec);
else
    ncornersperplanevec = uint16(ncornersperplanevec);
end

Varlist =        [ ' edgecorners planesatedge closwedangvec planehasindents indentingedgepairs'];
Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
Varlist = [Varlist,' corners planecorners  planeeqs planenvecs ncornersperplanevec cornerinfrontofplane'];
Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
Varlist = [Varlist,' minvals maxvals offedges '];
Varlist = [Varlist,'  int_or_ext_model'];
if difforder >= 1
    edgedata.edgestartcoordsnudge = edgestartcoordsnudge;
    edgedata.edgeendcoordsnudge = edgeendcoordsnudge;
    edgedata.edgenormvecs = edgenormvecs;
    Varlist = [Varlist, '  edgedata '];
end
eval(['save ',outputfile,Varlist])
author	Chris Cannam
date	Tue, 30 Sep 2014 16:23:00 +0100
parents	90220f7249fc
children