function path=shortest_path(mesh,nfrom,nto)
%function path=shortest_path(mesh,nfrom,nto)
%
% returns a list of u,v points that makes the shortest path 
% that goes from  u,v point  nfrom to u,v point nto
%
% uses Dijkstra's algorythem from Aho, Hopcroft and Ullman
%

%make some shorthand
 edges=mesh.surround;
 x=mesh.uvnode(:,1);
 y=mesh.uvnode(:,2);
 N=length(x);
  
 
% set up Dpath(N), a list of the quickest way from and u,v point nto to  
% nfrom.  Dpath(n) is number of the node that is next on the journey from 
% m to nfrom;
%
% Also setup distpath, the distance to nto from any point in the mesh.
%  and setup freeedge, which is true when there is still an edge from
%  a u,v point to another u,v point which is not yet on a shortest path.
%  freeedge is used to make the search run faster, by preventing it from
%  checking nodes that have no free edges.
%
  Dpath=nan*zeros(N,1);
  distpath=nan*zeros(N,1);
  freepath=ones(N,1);
  Dpath(nfrom)=nfrom; %to go home from home, stay at home.....
  distpath(nfrom)=0;  %at home, you are 0 away from home....
 
% create the set "on_edge" which contains the points that are on a
% shortest path to nfrom, AND still has connected points which are 
% not on a shortest path.  Only the members of on_edge are candidates 
% for the next edge to be added to a shortest path
  on_edge=[nfrom];
  
%compute distance between nodes, and store in dist.
%mesh.surround (here stored in edges)
%contains the numbers of the u,v points surrounding a given u,v point.  
%if an ellement in edges is zero, it means that it is on a boundary.  
%those boundary points will be assumed to be infinitely far away, and thus
%not on any shortest path. The index in edge will be set to some arbitary
%to keep from getting bad array references
  dist=zeros(size(edges))*inf;
  for n=1:3
    edgevec=edges(:,n);
    
    %if an element of edgevec is zero, that means we are on a boundary--
    %make that edge a link back to itself, of infinite length
    indx=find(edgevec==0);
    edgevec(indx)=indx;
    edges(indx,n)=indx;
    
    %compute distance
    dist(:,n)=sqrt((x(edgevec)-x).^2+(y(edgevec)-y).^2);
    
    %make infinite the boundary points
    dist(indx,n)=inf;
  end

  
  %remove the starting point from posible routes on the shortest path
  %to the starting point
  newnode=nfrom;
  for n=1:3
    indx=find(edges(edges(newnode,n),:)==newnode);
    dist(edges(newnode,n),indx)=inf;
  end

%now do Dijkstra's algorythem.... Find an edge to add to the collection of 
%shortest paths which makes the shortest new path.  
%
%Once that is done, make the distance between the the edge added 
%to the list of shortest paths infinite
%
  done_nodes=0;
  while isnan(Dpath(nto))
    done_nodes=done_nodes+1;
    
    %find all points already on a shortest path
    
    % now find the edge between a point already on a shortest path, and 
    % any point not on a shortest path, which makes a the shortest possible
    % new path. 
    %
    % The shortes new edge will be
    % from the u,v point already on the path "on_edge(ndistsmall)" and will be
    % from that point along the edge "candidates(ndistsmall)"

    [distcand,candidates]=...
	min(dist(on_edge,:)+repmat(distpath(on_edge),1,3),[],2);
    [distsmall, ndistsmall]=min(distcand);
    
    %this node is the new node which is on a shortest path, via 
    %on_edge(ndistsmall)
    newnode=edges(on_edge(ndistsmall),candidates(ndistsmall));
    on_edge = [on_edge newnode];
    if (isnan(Dpath(newnode)))
      Dpath(newnode)=on_edge(ndistsmall);
      distpath(newnode)=distpath(on_edge(ndistsmall))+distcand(ndistsmall);
    else
      keyboard
      error('something wrong -- tried to add same node onto path twice')
    end

    %can no longer take this path again
    dist(on_edge(ndistsmall),candidates(ndistsmall))=inf;

    %check if there are any free edges from node on_edge(ndistsmall).
    %if there are not, take it off the search list.
    if isempty(find(isfinite(dist(on_edge(ndistsmall),:))))
      on_edge(ndistsmall)=[];
    end

    
    %now make the length FROM any other node to the new node 
    %infinite, so we don't get any loops.  This relies on the 
    %fact that our graph is not really directed, i.e. no paths are 
    %one way. 
    for n=1:3
      indx=find(edges(edges(newnode,n),:)==newnode);
      dist(edges(newnode,n),indx)=inf;    
    end

    %if newnode==9503
    %  keyboard
    %end
    
% $$$     plot([x(on_edge(ndistsmall)) ...
% $$$ 	  x(edges(on_edge(ndistsmall),candidates(ndistsmall)))],...
% $$$ 	[y(on_edge(ndistsmall)) ...
% $$$ 	  y(edges(on_edge(ndistsmall),candidates(ndistsmall)))],'r-');
% $$$     hold on

    if rem(done_nodes,1000)==0
      fprintf('done with %2.0f nodes in search for short path\n',done_nodes)
    end
  
  end

  fprintf('done\n')
  %pause
  
  %now extracte the path from node nto to node nfrom
  path=nto;
  while (path(end)~=nfrom)
    path=[path Dpath(path(end))];
  end