format compact
global myim ilab pts probs;
load myim;
%load myim2
[N,M] = size(myim); out = 4*N*M+1; all = (1:N*M)'; 
myim = reshape(1./myim, 1,N*M);
phi = 2*(.6728:.00157:1.12182); theta = 4.0198:.00316:4.82244; % for myim
%phi = 2*(.5298:.00157:1.10913); theta = 1.6158:.00315:2.51355; % for myim2
pts = [reshape(sin(phi')*cos(theta),1,N*M); ...
       reshape(sin(phi')*sin(theta),1,N*M); ...
       reshape(cos(phi')*ones(1,M),1,N*M)];
 
% Vertices are numbered 1 to N*M going down column by column, as in 'all'.
% 'neigh' gives the 4 neighbors of each vertex, (or 'out' for outside grid)
% in the order 3 o'clock, 6 o'clock, 9 o'clock, 12 o'clock. 
% Oriented edges are numbered from 1 to 4*M*N, using first by their 
% beginning vertex, second the 4 directions and including edges which point 
% outside the grid. 'edge' enumerates the edges coming into each vertex 
% (numbering as 'out' edges outside the grid).
neigh = zeros(N*M,4);
neigh(all,1) = all+N; neigh((N*M-N+1):(N*M),1) = out;
neigh(all,2) = all+1; neigh((N:N:N*M),2) = out;
neigh(all,3) = all-N; neigh(1:N,3) = out;
neigh(all,4) = all-1; neigh(1:N:(N*M-N+1),4) = out;
edge(all,1:4) = min(out, neigh(all,1:4) + N*M*ones(N*M,1)*[2 3 0 1]);

% 'f' takes 1,2 or 3 as 1st arg, then an oriented edge and outputs  
% the three edges feeding into that edge in clockwise order. 
% 'beg' assigns to an edge its beginning vertex.
f = zeros(N*M,4,3);
for j=1:3
   f(all,1:j,j)     = (neigh(all,j+1)+ N*M*mod(j+2,4))*ones(1,j);
   f(all,(j+1):4,j) = (neigh(all,j) + N*M*mod(j+1,4))*ones(1,4-j);
end
f = reshape(min(f,out),4*M*N,3)'; 
clear neigh;