function pnew = seglink(pold,width)
%
% seglink - determines 3 best segment neighbors both above and
%           below each segment [puts in properties 16-21]
% On input:
%     pold: segment properties
%     width: known width of band (or use avg)
% On output:
%     pnew: segment properties with links to neighbors added
%           [segment link with self number means nil]
% Call:
%     pnew = seglink(pold,width)
%
[n,dummy] = size(pold); pnew = pold;
for s1 = 1:n                              % set link to self
  pnew(s1,16:21) = [s1,s1,s1,s1,s1,s1];   % segdist returns big num 
end;                                      % for self distance

CLEFT = 1; CRIGHT = 2; CROW = 10; CCOL = 9;
for s1 = 1:n
  s1row = pnew(s1,CROW); s1col = pnew(s1,CCOL);
  s1left = pnew(s1,CLEFT); s1right = pnew(s1,CRIGHT);
  for s2 = 1:n
    if s1 ~= s2
      sup1 = pnew(s1,16); sup2 = pnew(s1,17); sup3 = pnew(s1,18);
      sdn1 = pnew(s1,19); sdn2 = pnew(s1,20); sdn3 = pnew(s1,21);
          % distance up
      dup1 = segdist(pnew(sup1,CROW),pnew(sup1,CCOL),s1row,s1col); 
      dup2 = segdist(pnew(sup2,CROW),pnew(sup2,CCOL),s1row,s1col);
      dup3 = segdist(pnew(sup3,CROW),pnew(sup3,CCOL),s1row,s1col);
         % distance dn
      ddn1 = segdist(pnew(sdn1,CROW),pnew(sdn1,CCOL),s1row,s1col);
      ddn2 = segdist(pnew(sdn2,CROW),pnew(sdn2,CCOL),s1row,s1col);
      ddn3 = segdist(pnew(sdn3,CROW),pnew(sdn3,CCOL),s1row,s1col);
	  s2row = pnew(s2,CROW); s2col = pnew(s2,CCOL);
	  s2left = pnew(s2,CLEFT); s2right = pnew(s2,CRIGHT);
	  if (abs(s1left-s2left)<width/2) | (abs(s1right-s2right)<width/2) 
               % if within column range
	    d12 = segdist(s1row,s1col,s2row,s2col); % get seg distance
		if s1row<s2row   % if seg2 below seg1
		  segs = [sdn1,sdn2,sdn3,s2]; dist = [ddn1,ddn2,ddn3,d12];
		  [dists,disti] = sort(dist);
		  pnew(s1,19) = segs(disti(1));
		  pnew(s1,20) = segs(disti(2));
		  pnew(s1,21) = segs(disti(3));
		elseif s1row>s2row    % if seg2 above seg1
		  segs = [sup1,sup2,sup3,s2]; dist = [dup1,dup2,dup3,d12];
		  [dists,disti] = sort(dist);
		  pnew(s1,16) = segs(disti(3));
		  pnew(s1,17) = segs(disti(2));
		  pnew(s1,18) = segs(disti(1));
		end
	  end
	end
  end
end
