home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 2.ddi / SPARFUN.DI$ / UNMESH.M < prev   
Encoding:
Text File  |  1993-03-07  |  1.9 KB  |  81 lines
  1. function [A,xy] = unmesh(M)
  2. %UNMESH    Convert a list of mesh edges to a graph or matrix.
  3. %    [A,xy] = UNMESH(M)
  4. %
  5. %    Input:
  6. %      Each row of M is a mesh edge in two dimensions, [x1 y1 x2 y2].
  7. %
  8. %    Output:
  9. %      A is the Laplacian matrix of the mesh (the symmetric adjacency matrix 
  10. %           with -1 for edges and degrees on diagonal).
  11. %      Each row of xy is a coordinate [x y] of a mesh point.
  12. %
  13. %    See also MPLOT, GPLOT, SUBMESH.
  14.  
  15. %    John Gilbert, 1990.
  16. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  17.  
  18. % Discretize x and y with "range" steps, 
  19. % equating coordinates that round to the same step.
  20.  
  21.  
  22. range = round(eps^(-1/3));
  23.  
  24. [m,k] = size(M);
  25. if k ~= 4, error ('Mesh must have rows of the form [x1 y1 x2 y2].'), end;
  26.  
  27. x = [ M(:,1) ; M(:,3) ];
  28. y = [ M(:,2) ; M(:,4) ];
  29. xmax = max(x);
  30. ymax = max(y);
  31. xmin = min(x);
  32. ymin = min(y);
  33. xscale = (range-1) / (xmax-xmin);
  34. yscale = (range-1) / (ymax-ymin);
  35.  
  36. % The "name" of each (x,y) coordinate (i.e. vertex)
  37. % is scaledx + scaledy/range .
  38.  
  39. xnames = round( (x - xmin*ones(2*m,1)) * xscale );
  40. ynames = round( (y - ymin*ones(2*m,1)) * yscale );
  41. xynames = xnames+1 + ynames/range;
  42.  
  43. % vnames = the sorted list of vertex names, duplicates removed.
  44.  
  45. vnames = sort(xynames);
  46. f = find(diff( [-Inf; vnames] ));
  47. vnames = vnames(f);
  48. n = length(vnames);
  49. disp ([int2str(n) ' vertices:']);
  50.  
  51. % x and y are the rounded coordinates, un-scaled.
  52.  
  53. x = (floor(vnames) / xscale) + xmin;
  54. y = ((vnames-floor(vnames)) / yscale) * range + ymin;
  55. xy = [x y];
  56.  
  57. % Fill in the edge list one vertex at a time.
  58.  
  59. ij = zeros(2*m,1);
  60. for v = 1:n,
  61.     if ~rem(v,10),
  62.         disp ([int2str(v) '/' int2str(n)]);
  63.     end;
  64.     f = find( xynames == vnames(v) );
  65.     ij(f) = v*ones(length(f),1);
  66. end;
  67. if rem(n,10), disp ([int2str(n) '/' int2str(n)]); end;
  68.  
  69. % Fill in the edges of A.
  70.  
  71. i = ij(1:m);
  72. j = ij(m+1:2*m);
  73. A = sparse(i,j,1,n,n);
  74.  
  75. % Make A the symmetric Laplacian.
  76.  
  77. A = -spones(A+A');
  78. A = A - diag(sum(A));
  79.  
  80.  
  81.