function outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches)
% EDB1ed2geox - Calculates 2nd- and higher-order edge-related geom. parameters.
% EDB1ed2geox 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:
%	eddatafile	            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'.
%   sdatafile
%   rdatafile
%   doSbatches
%   doRbatches
%   specorder               The highest order of any reflection kind (specular and/or diffraction).
%	difforder 	            The highest order of diffraction. If it is 0 or 1 then the parameter
%							edgeseespartialedge is not calculated. Default: 1
%   nedgesubs
%   ndiff2batches
%	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
%    (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.
%   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 can never see the 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 can potentially see the plane.
%
%	 (Calculated by EDB1edgeo, but only for diffraction order >= 2:)
%   edgeseespartialedge     A sparse matrix, [nedges,nedges], with the
%                           values 0, pos. integer or neg. integer:
%                           0           Two edges can not see each other (or one
%                                       of the edges is inactive)
%                           pos.int.    A value from 1 to 2^nedgesubs-1
%                                       indicating how much of the edges
%                                       see each other. A pos. value
%                                       indicates that the edge-to-edge
%                                       path runs along a plane.
%                                       NB! A complete visibility test is
%                                       not made; subsegment 1 on edge 1 is only
%                                       checked against subsegment 1 on
%                                       edge 2, not against all other
%                                       subsegments of edge 2.
%                           neg.int.    Same as for the pos.int. except
%                                       that the edge-to-edge path does not
%                                       run along a plane.
%   
%   reftoshortlistE
%   re1sho, re2sho
%   thetae1sho, thetae2sho
%	ze1sho, ze2sho
%   examplecombE	
%   edgealignedwithedge     A sparse matrix, [nedges,nedges], with the value 1
%                           or 0 stating whether an (active) edge is aligned
%                           with another (active) edge or not.
%   edgeperptoplane         A sparse matrix, [nplanes,nedges], with the value 1
%                           or 0 stating whether an (active) edge is perpendicular
%                           to an active plane or not.
%   edgeplaneperptoplane1   A sparse matrix, [nplanes,nedges], with the value 1
%                           or 0 stating whether an (active) edge has one
%                           of its two defining planes perpendicular to
%                           another plane with which it shares an edge.
%   edgeplaneperptoplane2   A sparse matrix, [nplanes,nedges], with the value 1
%                           or 0 stating whether an (active) edge has one
%                           of its two defining planes perpendicular to
%                           another plane which has a flat edge (i.e., a 180 degree edge).
%
% Uses the functions EDB1coordtrans1, EDB1coordtrans2
%                    EDB1infrontofplane, EDB1compress7, EDB1strpend, EDB1cross
%                    EDB1checkobstr_edgetoedge
%
% ----------------------------------------------------------------------------------------------
%   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) 20110822
%
% outputfile = EDB1ed2geox(eddatafile,sdatafile,rdatafile,doSbatches,doRbatches,specorder,difforder,nedgesubs,ndiff2batches);

% 20110822 Bug found and fixed, giving the wrong amplitude for second- and higher-order
% 		   diffraction, but only for diffraction paths that do not travel along a surface. 

global SHOWTEXT 

geomacc = 1e-10;

eval(['load ',eddatafile])
edgestartcoordsnudge =   edgedata.edgestartcoordsnudge;
edgeendcoordsnudge =    edgedata.edgeendcoordsnudge;
edgenormvecs =    edgedata.edgenormvecs;    

[sfilepath,sfilestem,temp1] = fileparts(sdatafile);
sdatafile = [sfilepath,filesep,sfilestem];
if doSbatches == 0
    eval(['load ',sdatafile,'.mat'])
else
    eval(['load ',sdatafile,'_1.mat'])    
end
[rfilepath,rfilestem,temp1] = fileparts(rdatafile);
rdatafile = [rfilepath,filesep,rfilestem];
if doRbatches == 0
    eval(['load ',rdatafile,'.mat'])
else
    eval(['load ',rdatafile,'_1.mat'])
end

outputfile = eddatafile;

nplanes = size(planecorners,1);
nedges = size(edgecorners,1);

zerosvec1 = zeros(nplanes,1);

%---------------------------------------------------------------
% We check right away which edges that can not be seen by either the source
% or the receiver. We don't need to check edgeseesedge for combinations
% that involve any of these.
%
% Note that this is true only for difforder = 2! If difforder >= 3, then we
% can not exclude any edges based on this! (It could be possible to exclude
% edges, if we managed to "unwind" the
% edge-sees-an-edge-which-sees-an-edge-which-can-be-seen-by-the-source.

if difforder == 2
    edgesthatareseen = any( [vispartedgesfroms vispartedgesfromr].' );
    edgesnottoworryabout = [1:nedges];
    edgesnottoworryabout(edgesthatareseen) = [];
    clear vispartedgesfroms vispartedgesfromr visplanesfroms visplanesfromr
else
    edgesnottoworryabout = [];    
end

%---------------------------------------------------------------
%
% edgeperptoplane
%
%---------------------------------------------------------------
% Make matrices over:
%   1. Which edges are perpendicular to planes ("edgeperptoplane")
%   2. Which edges can give combinations with edge1-plane-edge1
%      where plane is perpendicular to one of the two planes that
%      defines the edge ("edgeplaneperptoplane1")
%   3. Which edges can give combinations with edge1-plane-edge1
%      where plane is a part of a (unnecessarily) divided plane
%      ("edgeplaneperptoplane2")
%
% NB! Only active edges and active planes can get
% such combinations.
%
% These are needed for diff-spec-diff combos, so only when
% difforder >= 2 (and specorder >= 3).

edgeperptoplane = [];
edgeplaneperptoplane1 = [];
edgeplaneperptoplane2 = [];
if specorder >= 3
    edgeperptoplane = sparse(zeros(nplanes,nedges));
    edgeplaneperptoplane1 = sparse(zeros(nplanes,nedges));
    edgeplaneperptoplane2 = sparse(zeros(nplanes,nedges));    
    if SHOWTEXT >= 4
		disp('         checking which edges are perpendicular to planes ...')
	end
    
    iv = find(edgeseesplane == 1);
    if ~isempty(iv)
        edgelist = floor((iv-1)/nplanes)+1;
        planelist = iv - (edgelist-1)*nplanes;

        % Check first which edges have vectors that are parallel to plane
        % normal vectors. These are "edgeperptoplane" cases. Clear them
        % afterwards from the edgelist and planelist.
        
        vec1 = edgenormvecs(edgelist,:);
        vec2 = planenvecs(planelist,:);
        ivperp = find( abs(  sum( (vec1.*vec2).' )  )==1);        
        if ~isempty(ivperp)
            edgeperptoplane(iv(ivperp)) = ones(size(ivperp));    
            clear iv
            edgelist(ivperp) = [];
            planelist(ivperp) = [];
            vec2(ivperp,:) = [];
        else
            clear iv    
        end

        % Now check which edges have plane1 or plane2 perpendicular to
        % other planes.

        if ~isempty(planelist)
            
            vec1 = planenvecs(planesatedge(edgelist,1),:);
            ivplaneperp1 = find( abs(  sum( (vec1.*vec2).' )  )==0);
            vec1 = planenvecs(planesatedge(edgelist,2),:);
            ivplaneperp2 = find( abs(  sum( (vec1.*vec2).' )  )==0);
            ivplaneperp = [ivplaneperp1 ivplaneperp2];
            if ~isempty(ivplaneperp)
                ivplaneperp = unique(ivplaneperp);
                edgelist = edgelist(ivplaneperp);
                planelist = planelist(ivplaneperp);
                edgeplane1 = planesatedge(edgelist,1);
                edgeplane2 = planesatedge(edgelist,2);

                % Check which of these other perpendicular planes that connnect 
                % to one of the edge planes via a 90 deg. corner.                

                ivtoclear = [];
                [ivfoundcombos,loc] = ismember([planelist edgeplane1],planesatedge,'rows');
                ivperp = find(ivfoundcombos);
                if ~isempty(ivperp)
                    edgenumb = loc(ivperp);
                    ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
                    if ~isempty(ivfinalcomb)
                        ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
                        edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));
                        ivtoclear = ivfinalcomb;
                    end
                end
                [ivfoundcombos,loc] = ismember([edgeplane1 planelist],planesatedge,'rows');
                ivperp = find(ivfoundcombos);
                if ~isempty(ivperp)
                    edgenumb = loc(ivperp);
                    ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
                    if ~isempty(ivfinalcomb)
                        ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
                        edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));

                        ivtoclear = ivtoclear.';
                        ivfinalcomb = ivfinalcomb.';
                        ivtoclear = [ivtoclear ivfinalcomb];
                    end
                end
                
                [ivfoundcombos,loc] = ismember([planelist edgeplane2],planesatedge,'rows');
                ivperp = find(ivfoundcombos);
                if ~isempty(ivperp)
                    edgenumb = loc(ivperp);
                    ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
                    if ~isempty(ivfinalcomb)
                        ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
                        edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));

                        if size(ivtoclear,2) == size(ivfinalcomb,2)
                            ivtoclear = [ivtoclear;ivfinalcomb];
                        elseif size(ivtoclear,1) == size(ivfinalcomb,1),                    
                            ivfinalcomb = ivfinalcomb.';
                            ivtoclear = [ivtoclear ivfinalcomb];
                        elseif isempty(ivtoclear)
                            error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')    
                        end
                    end
                end
                
                
                [ivfoundcombos,loc] = ismember([edgeplane2 planelist],planesatedge,'rows');
                ivperp = find(ivfoundcombos);
                if ~isempty(ivperp)
                    edgenumb = loc(ivperp);
                    ivfinalcomb = ivperp(find(closwedangvec(edgenumb)==3*pi/2));
                    if ~isempty(ivfinalcomb)
                        ivreftomatrix = planelist(ivfinalcomb) + (edgelist(ivfinalcomb)-1)*nplanes;
                        edgeplaneperptoplane1(ivreftomatrix) = ones(size(ivreftomatrix));

                        if size(ivtoclear,2) == size(ivfinalcomb,2)
                            ivtoclear = [ivtoclear;ivfinalcomb];
                        elseif size(ivtoclear,1) == size(ivfinalcomb,1),                    
                            ivfinalcomb = ivfinalcomb.';
                            ivtoclear = [ivtoclear ivfinalcomb];
                        elseif isempty(ivtoclear)
                            error('ERROR: Unexpected error involving sizes of ivtoclear and ivfinalcomb')    
                        end
                    end
                end
                
                % Finally check which edges have plane1 or plane2 perpendicular to
                % two other planes that connnect via a 180 deg. corner.

                if ~isempty(ivtoclear)
                    ivtoclear = unique(ivtoclear);
                    edgelist(ivtoclear) = [];
                    planelist(ivtoclear) = [];
                    if ~isempty(edgelist)
                        
                        ivflatedges = find(closwedangvec==pi);
                        if ~isempty(ivflatedges)
                            planecombos = planesatedge(ivflatedges,:);
                            for ii = 1:length(ivflatedges)
                                iv1 = find(planelist == planecombos(ii,1));                                
                                if ~isempty(iv1)
                                    ivreftomatrix = double(planecombos(ii,1)) + (edgelist(iv1)-1)*nplanes;
                                    edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));                                  
                                end
                                iv2 = find(planelist == planecombos(ii,2));
                                if ~isempty(iv2)
                                    ivreftomatrix = planecombos(ii,2) + (edgelist(iv2)-1)*nplanes;
                                    edgeplaneperptoplane2(ivreftomatrix) = ones(size(ivreftomatrix));                                  
                                end
                                
                            end
                        end
                        
                    end
                end
                  
            end
        end    
        
    end
       
end

%---------------------------------------------------------------
%
%			edgealignedwithedge
%
%---------------------------------------------------------------

edgealignedwithedge = [];

if specorder >= 2
    edgealignedwithedge = sparse(eye(nedges));
	if SHOWTEXT >= 4
		disp('         checking which edges are aligned with other edges ...')
	end
    
    listofactiveedges = [1:nedges].';
    listofactiveedges(offedges) = [];
    vecstocheck = edgenormvecs(listofactiveedges,:);
    iv = find(vecstocheck(:,1)<0);
    vecstocheck(iv,:) = -vecstocheck(iv,:);
    iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)<0);
    vecstocheck(iv,:) = -vecstocheck(iv,:);
    iv = find(vecstocheck(:,1)==0 & vecstocheck(:,2)==0 & vecstocheck(:,3)<0);
    vecstocheck(iv,:) = -vecstocheck(iv,:);
    if ~isempty(listofactiveedges)
        [uniquevecs,II,JJ] = unique(vecstocheck,'rows');
        N = histc(JJ,[1:length(listofactiveedges)]);
        iv = find(N>1);
        for ii = 1:length(iv)
            iv2 = find(JJ==iv(ii));
            edgestocheck = listofactiveedges(iv2);
            ntocheck = length(edgestocheck);
            cornernumberstocheck = edgecorners(edgestocheck,:);
            expandedstartcornermatrix = [reshape(repmat(cornernumberstocheck(:,1).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,1),[ntocheck 1])];
            expandedendcornermatrix   = [reshape(repmat(cornernumberstocheck(:,2).',[ntocheck 1]),ntocheck^2,1) repmat(cornernumberstocheck(:,2),[ntocheck 1])];
            toexpandededgestocheck   =         repmat(edgestocheck,  [ntocheck 1]);
            fromexpandededgestocheck = reshape(repmat(edgestocheck.',[ntocheck 1]),ntocheck^2,1);
            
            % Now we don't need to check edges against themselves, or edge1
            % vs edge2 if edge2 vs edge1 has already been tested.
            
            indmat = reshape([1:ntocheck^2],ntocheck,ntocheck);
            indmat = triu(indmat,1);
            indmat = indmat(find(indmat));
            
            expandedstartcornermatrix = expandedstartcornermatrix(indmat,:);
            expandedendcornermatrix   = expandedendcornermatrix(indmat,:);
            fromexpandededgestocheck = fromexpandededgestocheck(indmat,:);
            toexpandededgestocheck = toexpandededgestocheck(indmat,:);
            
            % We make the easiest check first: if two edges with the same
            % direction vector share one corner, they must be aligned.

            A = [expandedstartcornermatrix expandedendcornermatrix];
            A = sort(A.').';
            A = diff(A.').';
            iv3 = find(sum(A.'==0).');            
            validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
            if ~isempty(validcombs)
                indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
                edgealignedwithedge(indexvals) =  edgealignedwithedge(indexvals) + 1;     
                indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
                edgealignedwithedge(indexvals) =  edgealignedwithedge(indexvals) + 1;
                fromexpandededgestocheck(iv3) = [];
                toexpandededgestocheck(iv3) = [];
                expandedstartcornermatrix(iv3,:) = [];
                expandedendcornermatrix(iv3,:) = [];
            end
            
            % For the other edge-pair combinations we must check if the two
            % edges are really aligned. We do this by checking if two
            % vectors are aligned: edge1-vector and the vector from
            % edge1start to edge2start.            
                        
            direcvec = corners(expandedstartcornermatrix(:,2),:) - corners(expandedstartcornermatrix(:,1),:);
            direcveclengths = sqrt( sum(direcvec.'.^2) ).';
            direcvec = direcvec./(direcveclengths(:,ones(1,3)) + eps*10);
            test = abs(sum( (direcvec.*edgenormvecs(fromexpandededgestocheck,:)).' ));
            iv3 = find(abs(test-1)< 1e-8);
            validcombs = [fromexpandededgestocheck(iv3) toexpandededgestocheck(iv3)];
            if ~isempty(validcombs)
                indexvals = (validcombs(:,2)-1)*nedges + validcombs(:,1);
                edgealignedwithedge(indexvals) =  edgealignedwithedge(indexvals) + 1;     
                indexvals = (validcombs(:,1)-1)*nedges + validcombs(:,2);
                edgealignedwithedge(indexvals) =  edgealignedwithedge(indexvals) + 1;
                fromexpandededgestocheck(iv3) = [];
                toexpandededgestocheck(iv3) = [];
                expandedstartcornermatrix(iv3,:) = [];
                expandedendcornermatrix(iv3,:) = [];
            end            
        end
    end
    
end

%---------------------------------------------------------------
%
%			edgeseesedge
%
%---------------------------------------------------------------
%
% Make a matrix of which edges can see which edges
%
% Planes that are in front of the edge-planes might have visible edges.
% If closwedang = 0, then all other edges are potentially visible.
% If closwedang < pi, then all planes in front of plane1 or plane2 are OK
% If closwedang > pi, then all planes in front of plane1 and plane2 are OK

if SHOWTEXT >= 4
	disp('         checking which edges see which edges...')
end

edgeseesedge = int8(zeros(nedges,nedges));
for ii = 1:nedges
	closwed = closwedangvec(ii);
	if closwed == 0
		edgeseesedge(:,ii) = ones(nedges,1);
	else
		plane1 = planesatedge(ii,1);
		plane2 = planesatedge(ii,2);	
		okplanelist1 = find(planeseesplane(:,plane1)==1);
		okplanelist2 = find(planeseesplane(:,plane2)==1);		
		
		if closwed < pi
			okplanes = [okplanelist1;okplanelist2];
			okplanes = sort(okplanes);
		else
			okplanes = zerosvec1;
			okplanes(okplanelist1) = okplanes(okplanelist1)+1;
			okplanes(okplanelist2) = okplanes(okplanelist2)+1;
			okplanes = find(okplanes==2);
		end
		
		okedges = edgesatplane(okplanes,:);
		[n1,n2] = size(okedges);
		okedges = reshape(okedges,n1*n2,1);
		
		if ~isempty(okedges)
			okedges = okedges(find(okedges~=0));			
			edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1); 
		end

		edgesatplane1 = edgesatplane(plane1,:);
		edgesatplane2 = edgesatplane(plane2,:);
		okedges = [edgesatplane1 edgesatplane2].';
		okedges = okedges(find(okedges~=0));
		edgeseesedge(okedges,ii) = int8(double(edgeseesedge(okedges,ii)) + 1); 
	
	end	   	% else ..... if closwed == 0

end % for ii = 1:nedges

% Mask out all edge pairs that are aligned
edgeseesedge = double(edgeseesedge) - double(edgealignedwithedge);

% Make sure the edge-to-same-edge combos are masked out
% 	edgeseesedge = sign(edgeseesedge.*(edgeseesedge.')).*(1-eye(nedges));
edgeseesedge = int8( ( (edgeseesedge & (edgeseesedge.'))>0).*(1-eye(nedges)));

% Mask out all edges that should be switched off

mask = ones(nedges,nedges);
mask(offedges,:) = mask(offedges,:)*0;
mask(:,offedges.') = mask(:,offedges.')*0;

edgeseesedge = edgeseesedge & mask;
clear mask	


%-----------------------------------------------------------
% Go through all edgeseesedge combinations. If two edges see
% each other without belonging to the same plane, there should
% be a gain factor of 2.

if SHOWTEXT >= 4
	disp('         checking which edge-to-edge paths that run along planes')
end

edgetoedgestrength = 2*edgeseesedge;
for ii = 1:nedges
	if sum(edgetoedgestrength(:,ii)) ~= 0
		plane1 = planesatedge(ii,1);
		plane2 = planesatedge(ii,2);	
		edgesatplane1 = edgesatplane(plane1,:);
		edgesatplane2 = edgesatplane(plane2,:);
		sameplaneedges = [edgesatplane1 edgesatplane2].';
		sameplaneedges = sameplaneedges(find(sameplaneedges~=0));
		edgetoedgestrength(sameplaneedges,ii) = edgetoedgestrength(sameplaneedges,ii)*0 + 1;  
	end % % if sum(edge
end

% an edge can not contribute to itself

edgetoedgestrength = edgetoedgestrength.*(1-eye(nedges));

% The edges belonging to the same planes, but that are inactive (90 degrees etc)
% must be switched off here too.
% 011021 This is redundant.
 edgetoedgestrength = edgetoedgestrength.*(edgeseesedge>0);
edgetoedgestrength = int8(edgetoedgestrength);

%----------------------------------------------------------------------------
%
%		CHECK OBSTRUCTION OF EDGE-TO-EDGE PATHS
%
%----------------------------------------------------------------------------
%
% We construct two long lists: 'from-edges' and 'to-edges' for which
% obstruction should be tested. These lists are made up of all the
% combinations in edgeseesedge that have a value 1. 
% NB! We use symmetry - if edge 1 can see edge 2, then edge 2 can also see
% edge 1.
%
% Both edges should be subdivided into nedgesubs segments. So, the long
% lists must be expanded so that each from-edge/to-edge combination is
% replaced by nedgesubs^2 as many.
% For efficiency, we make a shortlist of the actual edge subdivisions and
% we then expand long lists with pointers to this shortlist.
%
% The final matrix, edgeseespartialedge, will have an integer value which
% tells whether two edges see each other:
%
%   0               Edges can not see each other, or one of the edges is inactive
%                   (an inactive edge has an integer wedge index, or has one plane
%                   which is TOTABS)
%   2^nedgesubs-1   Two edges see each other completely.
%   Other integers  Two edges see each other partly.
%   Neg. integer    As above, but in addition, the edges do not belong
%                   to the same plane.
%
% We check in the visedgesfroms and visedgesfromr lists to check which
% combinations we don't need to check.

if SHOWTEXT >= 3
	disp(['   Checking for obstructing planes between edges and edges'])
    disp(' ')
    disp('NB!!! The value of nedgesubs is temporarily set to 1 for the edge-to-edge vis. test')
end

obstructtestneeded = (sum(canplaneobstruct)~=0);
maxedgetoedgevisvalue = 2^(2*nedgesubs)-1;

if obstructtestneeded
    
    nedgesubsorig = nedgesubs;
    nedgesubs = 1;
    edgeseesedge = triu(edgeseesedge);
    
    iv = full(find(edgeseesedge~=0));
    
    if ~isempty(iv)
        fromedges = floor(iv/nedges)+1;
        toedges = iv - (fromedges-1)*nedges;
        ncombs = length(fromedges);

        ivcancel = find(ismember(fromedges,edgesnottoworryabout));
        fromedges(ivcancel) = [];
        toedges(ivcancel) = [];

        ivcancel = find(ismember(toedges,edgesnottoworryabout));        
        toedges(ivcancel) = [];
        fromedges(ivcancel) = [];
        clear ivcancel

        % Some of these fromedges-toedges combos might involve two
        % neighboring edges that form an indent. Those should be
        % removed before we start checking fro obstructions.
        %
        % We check the [fromedges toedges] matrix against the
        % indentingedgepairs matrix.
        
        [tf,loc]= ismember(indentingedgepairs,sort([fromedges toedges].').','rows');
        if ~isempty(loc)
            ivmatch = find(loc);
           fromedges(loc(ivmatch)) = [];
           toedges(loc(ivmatch)) = [];
        end
        ncombs = length(fromedges);

        [uniqueedges,iunique,junique] = unique([fromedges toedges]);
        clear fromedges toedges
        
        % The lists edgesubcoords are the shortlists.
        
        [edgesubcoords,edgeweightlist,edgenumberlist] = EDB1getedgepoints(edgestartcoords(uniqueedges,:),edgeendcoords(uniqueedges,:),edgelengthvec(uniqueedges,:),nedgesubs,1);

        % The two lists below contain pointers to the first edge segment in the
        % shortlist.
	
        fromcoords_refp1 = nedgesubs*(junique(1:ncombs)-1)+1;
        tocoords_refp1    = nedgesubs*(junique(ncombs+1:2*ncombs)-1)+1;
        clear junique
        
        % Now the original [fromedges toedges] could be recovered by:
        % uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)])
        % reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),6,2)
%         uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]);
%         npairs = length(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]));
%         reshape(uniqueedges([edgenumberlist(fromcoords_refp1) edgenumberlist(tocoords_refp1)]),npairs/2,2)

        addmask = [0:nedgesubs-1];
        onesvec1 = ones(1,nedgesubs);
        ntot1 = ncombs*nedgesubs;
	
        % First, the from-list gets repetitions of the same values
        expandfrom_ref = fromcoords_refp1(:,onesvec1);
        clear fromcoords_refp1;
        expandfrom_ref = reshape(expandfrom_ref.',ntot1,1);
        
        % Second, the expanded from-list gets repetitions that are
        % incremented in steps of 1 (because the shortlists have the
        % edge segments placed after each other).
        expandfrom_ref = expandfrom_ref(:,onesvec1);
        expandfrom_ref = expandfrom_ref + addmask(ones(ntot1,1),:);
        expandfrom_ref = reshape(expandfrom_ref.',ncombs*nedgesubs^2,1);
        
        % Same thing for the to-list except that the first expansion
        % gives repetitions that are incremented in steps of 1.
        expandto_ref   = tocoords_refp1(:,onesvec1);
        clear tocoords_refp1;
        expandto_ref = expandto_ref + addmask(ones(ncombs,1),:);
        expandto_ref = reshape(expandto_ref.',ntot1,1);
	
        % Second expansion for the to-list: repetitions
        expandto_ref = expandto_ref(:,onesvec1);
        expandto_ref = reshape(expandto_ref.',ncombs*nedgesubs^2,1);
	
        % Now we have expanded lists of all combinations, and we can 
        % make lists with the coordinates, the edge numbers, and the weights of
        % the segments.
        
        % edgesubcoords(expandfrom_ref,:) will give all the starting points
        % "fromcoords"
        % edgesubcoords(expandto_ref,:) will give all the ending points
        % "tocoords"
        % uniqueedges(edgenumberlist(expandfrom_ref)) will give all
        % the starting edges, "fromedges"
        % uniqueedges(edgenumberlist(expandto_ref)) will give all
        % the ending edges, "toedges"
        % We can make a simple test of which planes that could be excluded
        % from the obstruction test by checking the edgeseesplane for the
        % relevant edges. It can be multiplied by the canplaneobstruct list
        % when calling the findobstructedpaths function.
        isplaneactiveforobstruction = (sum(edgeseesplane(:,uniqueedges).'==1)>0);
        
        global STARTPLANES ENDPLANES
        global BIGFROMCOORDS BIGTOCOORDS BIGSTARTPLANES BIGENDPLANES
        global REFTOFROMCOSHO REFTOTOCOSHO
        
        FROMCOORDSSHORTLIST = edgesubcoords;
        REFTOFROMCOSHO = uint32(expandfrom_ref);
        TOCOORDSSHORTLIST   = edgesubcoords;
        REFTOTOCOSHO = uint32(expandto_ref);
        if nedges < 256
            bigfromedge   = uint8(uniqueedges(edgenumberlist(expandfrom_ref)).');
        else
            bigfromedge   = uint16(uniqueedges(edgenumberlist(expandfrom_ref)).');
        end
        clear expandfrom_ref
        bigfromedge = bigfromedge(:);
        if nedges < 256
            bigtoedge     = uint8(uniqueedges(edgenumberlist(expandto_ref)).');
        else
            bigtoedge     = uint16(uniqueedges(edgenumberlist(expandto_ref)).');
        end
        clear expandto_ref
        bigtoedge = bigtoedge(:);
        STARTPLANES = [planesatedge(bigfromedge,1) planesatedge(bigfromedge,2)];
        ENDPLANES   = [planesatedge(bigtoedge,1) planesatedge(bigtoedge,2)];        
        shouldplanebechecked = canplaneobstruct.*isplaneactiveforobstruction;
        nplanesperbatch = ceil(sum(shouldplanebechecked)/ndiff2batches);
        if nplanesperbatch > 0
            temp = cumsum(shouldplanebechecked);
            diff2batchlist = zeros(ndiff2batches,2);
            diff2batchlist(1,1) = 1;
            loopisfinished = 0;
            for ii = 1:ndiff2batches-1
                ivtemp = find(temp==nplanesperbatch*ii);
                if ~isempty(ivtemp)
                    diff2batchlist(ii,2) = ivtemp(1);
                    diff2batchlist(ii+1,1) = ivtemp(1)+1;
               else
                   loopisfinished = 1;
               end               
            end
            diff2batchlist(ii,2) = nplanes;
            ivtemp = find(diff2batchlist(:,1)>0);
            diff2batchlist = diff2batchlist(ivtemp,:);
            ndiff2batches = size(diff2batchlist,1);

            npathstocheck = size(REFTOFROMCOSHO,1);
            for ii = 1:ndiff2batches
                if npathstocheck > 0
                    maskedchkpla = shouldplanebechecked;
                    if ii > 1
                        maskedchkpla(1:diff2batchlist(ii,1)-1) = maskedchkpla(1:diff2batchlist(ii,1)-1)*0;    
                    end
                    if ii < ndiff2batches
                        maskedchkpla(diff2batchlist(ii,2)+1:end) = maskedchkpla(diff2batchlist(ii,2)+1:end)*0;                    
                    end
                    nplanestocheck = sum(maskedchkpla);
                    if SHOWTEXT >= 3, 
                        if round(ii/1)*1==ii
                            disp(['   Batch ',int2str(ii),' of ',int2str(ndiff2batches),': '])
                            disp(['   ',int2str(npathstocheck),' paths to check, ',int2str(nplanestocheck),' planes to check obstruction for'])
                        end
                    end
                
                    % We create the BIGFROMCOORDS matrices and BIGTOCOORDS matrices
                    % outside since they need to be exactly identical from batch to
                    % batch as long as nplanestocheck is the same.
        
                    npathstocheck = size(REFTOFROMCOSHO,1);
                
                    ntot = nplanestocheck*npathstocheck;
                
                    BIGFROMCOORDS = reshape(repmat(FROMCOORDSSHORTLIST(REFTOFROMCOSHO,:).',[nplanestocheck,1]),3,ntot).';
                    BIGTOCOORDS = reshape(repmat(TOCOORDSSHORTLIST(REFTOTOCOSHO,:).',[nplanestocheck,1]),3,ntot).';
                    BIGSTARTPLANES = reshape(repmat(STARTPLANES.',[nplanestocheck,1]),2,ntot).';
                    BIGENDPLANES = reshape(repmat(ENDPLANES.',[nplanestocheck,1]),2,ntot).';
                
                    [nonobstructedpaths,nobstructions,edgehits,cornerhits] = EDB1checkobstr_edgetoedge(maskedchkpla,planeseesplane,...
                        planeeqs,planenvecs,minvals,maxvals,planecorners,corners,ncornersperplanevec,rearsideplane);
                
                    
                    if ~isempty(edgehits)
                           % Look for the edge-to-edge combos that share a
                           % plane. When such combos have hit an edge, it
                           % must mean that they are planes with indents.
                          wasedgehitnotindent = prod(double(diff(sort([STARTPLANES(edgehits,:) ENDPLANES(edgehits,:)].')))).';
                        indentedgehits = find(wasedgehitnotindent==0);
                        if ~isempty(indentedgehits)
                          nonobstructedpaths = setdiff(nonobstructedpaths,edgehits(indentedgehits));
                        end
                        
                    end
                    
                    if ~isempty(nonobstructedpaths)
                        REFTOFROMCOSHO = REFTOFROMCOSHO(nonobstructedpaths);
                        REFTOTOCOSHO = REFTOTOCOSHO(nonobstructedpaths);
                        STARTPLANES = STARTPLANES(nonobstructedpaths,:);
                        ENDPLANES = ENDPLANES(nonobstructedpaths,:);
                        bigfromedge = bigfromedge(nonobstructedpaths);
                        bigtoedge = bigtoedge(nonobstructedpaths);
                        npathstocheck = size(REFTOFROMCOSHO,1);
                    else
                        REFTOFROMCOSHO = [];
                        REFTOTOCOSHO = [];
                        STARTPLANES = [];
                        ENDPLANES = [];
                        bigfromedge = [];
                        bigtoedge = [];
                        npathstocheck = 0;
                    end
                   
                end   %     if npathstocheck > 0
            
            end           %  for ii = 1:ndiff2batches
        else
            nonobstructedpaths = 1;
            
        end         %if nplanesperbatch > 0

            clear REFTOFROMCOSHO REFTOTOCOSHO STARTPLANES ENDPLANES


        if ~isempty(nonobstructedpaths)
            
            % Add the segment pieces together
            % NB! We don't try to implement the correct way to do it since
            % we would need nedgesubs^2 bits to represent all edge-segment
            % to edge-segment visibility possibilities.
            % Instead we just establish whether the two edges see each
            % other fully or partly.

            visiblesegmentscounter = ones(size(bigfromedge));
            test = [bigfromedge bigtoedge];
            
            ncombs = length(bigfromedge);
            dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
            ivremove = find(dtest==1);
            
            while ~isempty(ivremove)
                visiblesegmentscounter(ivremove+1) = visiblesegmentscounter(ivremove+1) + visiblesegmentscounter(ivremove);
                visiblesegmentscounter(ivremove) = [];
                bigfromedge(ivremove) = [];
                bigtoedge(ivremove,:) = [];
                
                test = [bigfromedge bigtoedge];
                ncombs = length(bigfromedge);
                dtest = diff([0;prod(   double(test(1:ncombs-1,:)==test(2:ncombs,:)).'  ).']);
                ivremove = find(dtest==1);
	
            end
            
            indexvec = uint32((double(bigfromedge)-1)*nedges + double(bigtoedge));
            if maxedgetoedgevisvalue < 128
                edgeseespartialedge = int8(zeros(nedges,nedges));            
            elseif maxedgetoedgevisvalue < 32768
                edgeseespartialedge = int16(zeros(nedges,nedges));            
            else
                edgeseespartialedge = int32(zeros(nedges,nedges));            
            end                
                
            iv = find(visiblesegmentscounter==nedgesubs^2);
            edgeseespartialedge(indexvec(iv)) = maxedgetoedgevisvalue*ones(size(iv));
            iv = find(visiblesegmentscounter>0 & visiblesegmentscounter<nedgesubs^2);
            edgeseespartialedge(indexvec(iv)) = maxedgetoedgevisvalue/2*ones(size(iv));
            
            if maxedgetoedgevisvalue < 128
               edgeseespartialedge = int8(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');
           elseif maxedgetoedgevisvalue < 32768
               edgeseespartialedge = int16(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');
           else
               edgeseespartialedge = int32(double(edgeseespartialedge) + double(triu(edgeseespartialedge)).');               
           end
        else
            edgeseespartialedge = sparse(zeros(nedges,nedges));            
        end
    else
        edgeseespartialedge = [];
    end
else
    if maxedgetoedgevisvalue < 128
        edgeseespartialedge = int8((edgeseesedge*maxedgetoedgevisvalue));            
    elseif maxedgetoedgevisvalue < 32768
        edgeseespartialedge = int16((edgeseesedge*maxedgetoedgevisvalue));            
    else
        edgeseespartialedge = int32((edgeseesedge*maxedgetoedgevisvalue));            
    end                
    disp('      No obstruction tests needed!')
    
end   %(obstructtestneeded)

% For the special case of thin planes that are all in the same plane, only
% edge pairs that belong to the same plane should be included. That is
% when edgeseespartialedge is negative, then we should set
% edgeseespartialedge to zero.

% Bug detected 20110822 PS
% Old: allthinplanes = sign(1-sum(~planeisthin));
allthinplanes = sign(1-sign(sum(~planeisthin)));

% Bug detected 20110822 PS
% Old: allplanesinsameplane = 1 - sum(sum(planeseesplane==1))
allplanesinsameplane = 1 - sign(sum(sum(planeseesplane==1)));

multfactor = 1 - allthinplanes*allplanesinsameplane;

iv = find(edgetoedgestrength==2);
if ~isempty(iv)
    if maxedgetoedgevisvalue < 128
        edgeseespartialedge(iv) = int8(-double(edgeseespartialedge(iv))*multfactor);            
    elseif maxedgetoedgevisvalue < 32768
        edgeseespartialedge(iv) = int16(-double(edgeseespartialedge(iv))*multfactor);            
    else
        edgeseespartialedge(iv) = int32(-double(edgeseespartialedge(iv))*multfactor);            
    end                
end

clear edgeseesedge edgetoedgestrength

%----------------------------------------------------------------------------
%
%		CYLINDRICAL EDGE-TO-EDGE PARAMETERS
%
%----------------------------------------------------------------------------
%
% For each edge, calculate the cylindrical coordinates of the starting
% and ending points of all other edges relative to the first edge.

if difforder >= 2 & ~isempty(edgeseespartialedge)
	
	if SHOWTEXT >= 4
		disp('         E2E')
	end
	
    zerosvec4 = zeros(nedges,nedges);
	Bigre1 = (zerosvec4);
	Bigthetae1 = (zerosvec4);
	Bigze1 = (zerosvec4);
	Bigre2 = (zerosvec4);
	Bigthetae2 = (zerosvec4);
	Bigze2 = (zerosvec4);
    clear zerosvec4
    
	for edge1 = 1:nedges
		if SHOWTEXT >= 4
            if round(edge1/1)*1 == edge1
    			disp(['            Edge no. ',int2str(edge1)])
            end
        end
		edge1coords = [edgestartcoords(edge1,:);edgeendcoords(edge1,:)];
		iv = find(edgeseespartialedge(edge1,:)~=0).';        

		if ~isempty(iv),		
			% First find the subset of edges which belong to the same plane as the
			% edge itself. These should be treated separately for higher accuracy

			% iv1 will be the edges on the reference plane. They should have theta = 0.
			
            
			refplane = planesatedge(edge1,1);
            ncornersperplanevec = double(ncornersperplanevec);
			iv1 = edgesatplane( refplane,1:ncornersperplanevec(  refplane ));
			iv1 = iv1( find(iv1~= edge1)).';

			edgestart = edgestartcoordsnudge(iv1,:);
			edgeend   = edgeendcoordsnudge(iv1,:);

            [rs,slask,zs,rr,slask,zr] = ...
			EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
			Bigre1(iv1,edge1) = rs;
			Bigthetae1(iv1,edge1) = 0*iv1;
			Bigze1(iv1,edge1) = zs;
			Bigre2(iv1,edge1) = rr;
			Bigthetae2(iv1,edge1) = 0*iv1;
			Bigze2(iv1,edge1) = zr;

            [samevalues,iivec,jjvec] = intersect(iv,iv1);
% Bug found 20050421!!  PS
% % % Old wrong version:
% % %             %             sameedges = [iivec jjvec];
% % % 			if sum(sum(sameedges)) ~= 0
% % %             	iv( sameedges(1,:).' ) = [];
% % %             end
            if ~isempty(samevalues)
                iv(iivec) = [];
            end

            % iv2 will be the edges on the non-reference plane.
			% They should have theta = 2*pi - closthetavec.

			if ~isempty(iv)

				if planesatedge(edge1,2) > 0
					secplane = planesatedge(edge1,2);

					iv2 = edgesatplane( secplane,1:ncornersperplanevec(secplane));
					iv2 = iv2( find(iv2~= edge1)).';

					if ~isempty(iv2)
						edgestart = edgestartcoordsnudge(iv2,:);
						edgeend   = edgeendcoordsnudge(iv2,:);

                        [rs,slask,zs,rr,slask,zr] = ...
						EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));
						Bigre1(iv2,edge1) = rs;
						Bigthetae1(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
						Bigze1(iv2,edge1) = zs;
						Bigre2(iv2,edge1) = rr;
						Bigthetae2(iv2,edge1) = (2*pi-closwedangvec(edge1))*ones(size(iv2));
						Bigze2(iv2,edge1) = zr;

%						sameedges = EDB1findsame(iv,iv2);
% Bug found 20050421!!! PS
% % % % Old wrong version:
% % % %                         [samevalues,iivec,jjvec] = intersect(iv,iv1);
% % % %                         sameedges = [iivec;jjvec];
% % % % 						if sum(sum(sameedges)) ~= 0
% % % % 							iv( sameedges(1,:).' ) = [];
% % % %                         end

                        [samevalues,iivec,jjvec] = intersect(iv,iv2);                        
                        if ~isempty(samevalues)
                            iv(iivec) = [];
                        end
					end
				end
			end
			if ~isempty(iv)

				% Move the edge coordinates to be checked a short distance away
	
				edgestart = edgestartcoordsnudge(iv,:);
				edgeend   = edgeendcoordsnudge(iv,:);
	
                [rs,thetas,zs,rr,thetar,zr] = ...
				EDB1coordtrans2(edgestart,edgeend,edge1coords,edgenvecs(edge1,:));							
				Bigre1(iv,edge1) = rs;
				Bigthetae1(iv,edge1) = thetas;
				Bigze1(iv,edge1) = zs;
				Bigre2(iv,edge1) = rr;
				Bigthetae2(iv,edge1) = thetar;
				Bigze2(iv,edge1) = zr;

			end
		end
    end

    %-----------------------------------------------------------
    % Go through all edges that are in-plane with each other, that is, two
    % edges have at least one plane each that are in-plane with each other.

    % First we identify all planes that have at least one more co-planar
    % plane
    
    planehascoplanar = any(planeseesplane == -1);
    
    % Then we go through the list of edges and every edge with a connected
    % plane that has another co-planar plane must also have potential
    % in-plane edges.

    % For each edge, we check if there are other edges (that don't belong to the same plane!!)
    % with thetaangle = 0 or thetaangle = 2*pi. If there are other such
    % edges, those edge-to-edge combinations should be shut off.
    %
    % After all such edge-to-edge combinations have been shut off, we still
    % need another pass to cancel edge-to-edge paths that pass entirely
    % across other edges.
    
    ivedges = 1:nedges;
    ivedges(offedges) = [];
    
    for ii = ivedges
       if any( planehascoplanar(planesatedge(1,:)) )
          
           ivcancelcombs = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ( (Bigthetae1(:,ii) == 0) | (Bigthetae1(:,ii) == 2*pi) ) );           
           if ~isempty(ivcancelcombs)
               edgeseespartialedge(ivcancelcombs,ii) = 0;
               edgeseespartialedge(ii,ivcancelcombs.') = 0;
           end
           
       end
        
    end
    
    % The only edge-to-edge combinations that we can easily discard are
    % those for which edges are aligned with each other (same zstart and
    % zend): only the closest of those should be kept!
    
    for ii = ivedges
       if any( planehascoplanar(planesatedge(1,:)) )

           
           ivotheredges = find( (edgeseespartialedge(:,ii) == -maxedgetoedgevisvalue) & ((Bigthetae1(:,ii) == pi)) );
           if ~isempty(ivotheredges)
               
               nudgeval = 1e-10*edgelengthvec(ii)*100;
                ivselectedges = find( abs(Bigze1(ivotheredges,ii))<nudgeval | abs(Bigze2(ivotheredges,ii))<nudgeval);
                if length(ivselectedges) > 1
                   disp(['Edge no. ',int2str(ii)])
                    
                    lengthok = abs([Bigze1(ivotheredges(ivselectedges),ii)  Bigze2(ivotheredges(ivselectedges),ii)] - edgelengthvec(ii)) < nudgeval;
                    ivstillok = find(any(lengthok.'));
                    
                    if ~isempty(ivstillok)
                        ivotheredges = ivotheredges(ivselectedges(ivstillok));
                    else
                       ivotheredges = []; 
                    end

                    meanradialdist = mean( [Bigre1(ivotheredges,ii) Bigre2(ivotheredges,ii)].'    );
                    
                    ivshutoff = ivotheredges;
                    noshutoff = find(meanradialdist == min(meanradialdist));
                    ivshutoff(noshutoff) = [];

                    if ~isempty(ivshutoff)
                       edgeseespartialedge(ivshutoff,ii) = 0; 
                       edgeseespartialedge(ii,ivshutoff) = 0; 
                    end
                    
                end           
            
           end
           
       end
        
    end    
    
	%-----------------------------------------------------------------------
	% Find identical combinations

    if SHOWTEXT >= 3
        disp('   Looking for identical combinations')    
    end
	[reftoshortlistE,re1sho,re2sho,thetae1sho,thetae2sho,...
	ze1sho,ze2sho,edgelengthsho,examplecombE] = EDB1compress7(Bigre1,Bigre2,...
	Bigthetae1,Bigthetae2,Bigze1,Bigze2, edgelengthvec(:,ones(1,nedges)).');    

else
	reftoshortlistE = [];
	re1sho = [];
	re2sho = [];
	thetae1sho = [];
	thetae2sho = [];
	ze1sho = [];
	ze2sho = [];
	examplecombE = [];
end

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

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

ncorners = size(corners,1);
if ncorners < 256
    planecorners = uint8(planecorners);
elseif ncorners < 65536
    planecorners = uint16(planecorners);    
end   
planeisthin = uint8(planeisthin);

Varlist =        [ ' edgecorners planesatedge closwedangvec planehasindents'];
Varlist = [Varlist,' planeisthin planeseesplane rearsideplane edgeseesplane canplaneobstruct'];
Varlist = [Varlist,' corners planecorners  planeeqs planenvecs ncornersperplanevec'];
Varlist = [Varlist,' edgestartcoords edgeendcoords edgenvecs reflfactors edgesatplane edgelengthvec'];
Varlist = [Varlist,' minvals maxvals offedges edgestartcoordsnudge edgeendcoordsnudge'];
Varlist = [Varlist,' reftoshortlistE re1sho re2sho thetae1sho thetae2sho ze1sho ze2sho examplecombE'];
Varlist = [Varlist,' edgeseespartialedge int_or_ext_model'];
Varlist = [Varlist,' edgealignedwithedge edgeperptoplane edgeplaneperptoplane1 edgeplaneperptoplane2'];

eval(['save ',outputfile,Varlist])
