function rte = tspspfillcur(XY,k)
%TSPSPFILLCUR Spacefilling curve algorithm for TSP route construction.
%   rte = tspspfillcur(XY,k)
%    XY = n x 2 matrix of vertex cooridinates
%     k = (optional) no. of binary digits to use in SFCPOS (= 8, default)
%   rte = (n + 1)-vector of vertices
%
% (Based on algorithm in J.J. Bartholdi and L.K. Platzman, Management Sci.,
%  34(3):291-305, 1988)
%
% Example:
% vrpnc1
% rte = tspspfillcur(XY);
% TD = rteTC(rte,dists(XY,XY,2))
% h = pplot(XY,'r.');
% pplot(XY,num2cell(1:size(XY,1)))
% pplot({rte},XY,'m')

% Copyright (c) 1994-2006 by Michael G. Kay
% Matlog Version 9 13-Jan-2006 (http://www.ie.ncsu.edu/kay/matlog)

% Input Error Checking ****************************************************
error(nargchk(1,2,nargin));
if nargin < 2, k = []; end
if size(XY,2) ~= 2, error('XY must be a 2-column matrix.'), end
% End (Input Error Checking) **********************************************

[stheta,rte] = sort(sfcpos(XY(:,1),XY(:,2),k));
i1 = find(rte == 1);
rte = [rte([i1:end 1:i1-1])' 1];  % Make vertex 1 first in route
if nargout == 2
   TD = sum(dists(XY(rte,:),[],2));
end