mergeseg(seglist, angtol, linkrad, linedevtol)

% MERGESEG - Line segment merging function.
%
% Usage: newseglist = mergeseg(seglist, angtol, linkrad)
%
% Arguments:  seglist - an Nx4 array storing line segments in the form
%                       [x1 y1 x2 y2
%                        x1 y1 x2 y2
%                           . . .   ] etc
%
%             angtol  - Angle tolerance used when attempting to merge line
%                       segements (radians).
%             linkrad - Maximum distance between end points of line segments
%                       for segments to be elegible for linking (pixels).
%
% Function scans through the list of line segments seeking to merge segments
% together.  Segments are merged if the orientation difference is less than
% angtol and if the ends of the segments are within linkrad of each other.
% The comparison is performed by first sorting line segments by angle, then
% each line segment is only tested against other line segments that satisfy
% the orientation constraint.
%
% See also:  EDGELINK, LINESEG, MAXLINEDEV, DRAWSEG
%

% Copyright (c) 2000-2005 Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% March    2000 - Original version
% November 2005 - Bug fix for case Nseg < 4 (thanks to Ming Luo)
% May      2006 - Bug fix for termination of while loop (thanks to Martin Labrecque)

function newseglist = mergeseg(seglist, angtol, linkrad, linedevtol)
    
    Nseg = size(seglist,1);
    
    if Nseg <= 1    % Nothing to merge
	newseglist = seglist;
	return;
    end
    
    cosAngTol = cos(angtol);
 
    % Build up some data in the following form
    % [x1 y1 x2 y2 dx dy angle merged]
    
    seg = zeros(Nseg,8);
    seg(:,1:4) = seglist;
    seg(:,5:6) = seg(:,3:4)-seg(:,1:2);     % delta x and delta y
    seg(:,7) = atan2(seg(:,6), seg(:,5));   % =segment angle [-pi ... pi]
    neg = seg(:,7) < 0;                     % compress angle range [0...pi]
    seg(:,7) = seg(:,7) + neg*pi;
    [Y,I] = sort(seg(:,7));                 % sort by angle
    seg = seg(I,:);                      
    seg(:,7) = 2*seg(:,7);                  % double the angles
    
    TotalSegs = Nseg;
    
    for s1 = 1:Nseg
%	fprintf('Merging Segment %d/%d\r',s1,Nseg);	    
	if ~seg(s1,8)  % if segment s1 has not already been merged
	    s2 = s1+1; if s2>Nseg, s2=s2-Nseg; end  % `modulo' arithetic
            ang1 = seg(s1,7); ang2 = seg(s2,7);
	    deltaSin = sin(ang1)*cos(ang2) - cos(ang1)*sin(ang2);
	    deltaCos = cos(ang1)*cos(ang2) + sin(ang1)*sin(ang2);
	    deltaAng = 0.5*abs(atan2(deltaSin,deltaCos));

	    while (s2 ~= s1) && (deltaAng < angtol)
		if ~seg(s2,8)   % if segment s2 has not already been merged
		    [newseg, merged] = compareseg(seg(s1,:), seg(s2,:),linkrad, ...
						  linedevtol);
		    if merged
			seg(s1,1:4) = newseg;
			%  update delta x, delta y, and angle ? (I think not)
			seg(s2,8) = merged;  % eliminate s2;
			TotalSegs = TotalSegs - 1;
		    end
		end
		s2 = s2+1; if s2>Nseg, s2=s2-Nseg; end		
		ang2 = seg(s2,7);
		deltaSin = sin(ang1)*cos(ang2) - cos(ang1)*sin(ang2);
		deltaCos = cos(ang1)*cos(ang2) + sin(ang1)*sin(ang2); 
		deltaAng = 0.5*abs(atan2(deltaSin,deltaCos));
	    end
	end
    end

%    fprintf('\n');
    
    % Now clean up list by extracting the unmerged entries
    
    newseglist = zeros(TotalSegs,4);
    newNseg = 0;
    for s = 1:Nseg
	if ~seg(s,8)                 % if this segment has not been merged
	    newNseg = newNseg+1;     % store it in the final list.
	    newseglist(newNseg,:) =  seg(s,1:4);  
	end
    end
    
    
%----------------------------------------------------------------------    
% Internal function that does all the work.  
% Two segments, that already meet the angle constraint, are tested for end
% point proximity.  If they meet this constraint they then tested to see
% what the maximum line deviation would be if they were merged.  If this
% is within tolerance they are merged.
%----------------------------------------------------------------------    
    
function [newseglist, merged] = compareseg(seg1, seg2, linkrad, linedevtol)
    
    s1p1 = seg1(1:2);   % point1 of segment1
    s1p2 = seg1(3:4);   % point2 of segment1
    
    s2p1 = seg2(1:2);   % point1 of segment2
    s2p2 = seg2(3:4);   % point2 of segment2
    
    dir1 = seg1(5:6);
    dir2 = seg2(5:6);
    
    cosAng = dot(dir1,dir2);
    
    if cosAng < 0          % Reverse order of points on seg2 so that both
                           % segments are in the `same' direction
        tmp  = s2p1;
	s2p1 = s2p2;
	s2p2 = tmp;
    end

    merged = 0;                        % Default result: not merged
    newseglist = [seg1; seg2];         % and segments left unchanged
    
    if norm(s1p2-s2p1) < linkrad       % If we satisfy endpoint proximity
	rc = [s1p1; s1p2; s2p1; s2p2];
	[maxdev, index, D] = maxlinedev(rc(:,1),rc(:,2));  % test line deviation
	if maxdev < linedevtol
	    newseglist = [s1p1 s2p2];      % Merge the segments one way
	    merged = 1;
	end
    elseif norm(s2p2-s1p1) < linkrad 
	rc = [s2p1; s2p2; s1p1; s1p2];
	[maxdev, index, D] = maxlinedev(rc(:,1),rc(:,2));
	if maxdev < linedevtol	
	    newseglist = [s2p1 s1p2];      % ... or merge the other way.
	    merged = 1;	
	end
    end