function [I,J,L] = bnd_shares(bnd,Z,eps) %BND_SHARES computes the boundary shares for shape vector files %Used as a component of BND_SHARE_MATRIX %Written by: T.E.SMITH, 1/4/02 %INPUTS: bnd = boundary file output from Arcview with row number % as the identifier % Z = polyform(bnd); % eps = (optional) value for error parameter [default defined by % by differences between adjacent boundary points] %OUTPUTS: This is a 'sparse' representation of the boundary-share % matrix in which each row [I(r),J(r),L(r)] of the three % vectors I,J,L corresponds to a positive element in the % boundary-share matrix, i.e., spatial units I(r) and J(r) % share a boundary of length L(r). %PROGRAMS CALLED: polyform.m %Define "small" differences cutoff [assumes that absolute differences %between the first two points are typical] if nargin == 2 del = abs(bnd(2,:) - bnd(3,:)); %assumes differences eps = min([sum(del)/1000 10^(-6)]); %allows for the possibility that %sum(del) is large (as in a regular %lattice) end n = length(Z); M = zeros(n,4); % Form matrix of maxs and mins (enclosing rectangles) for i = 1:n Vi = bnd(Z(i,1):Z(i,2)-1,:); %i-th poly M(i,1) = max(Vi(:,1)); M(i,2) = min(Vi(:,1)); M(i,3) = max(Vi(:,2)); M(i,4) = min(Vi(:,2)); end %Estimate size of sparse(W) by elimination of non-adjacent pairs count = 0; for i = 1:n-1 Vi1_max = M(i,1); Vi1_min = M(i,2); Vi2_max = M(i,3); Vi2_min = M(i,4); % mbscalar(Vi1_max); % mbscalar(Vi1_min); % mbscalar(Vi2_max); % mbscalar(Vi2_min); Li = Z(i,2) - Z(i,1); % Number of i-faces Ri = Z(i,1) - 1; %Position of i-th lable row in bnd matrix for j = i+1:n %Grand j-loop Vj1_max = M(j,1); Vj1_min = M(j,2); Vj2_max = M(j,3); Vj2_min = M(j,4); % mbscalar(Vj1_max); % mbscalar(Vj1_min); % mbscalar(Vj2_max); % mbscalar(Vj2_min); if (Vi1_max + eps < Vj1_min)|(Vj1_max + eps < Vi1_min)|... (Vi2_max + eps < Vj2_min)|(Vj2_max + eps < Vi2_min) count = count; else count = count + 1; end end end %End count procedure W = zeros(2*count,3); % Given symmetry, 2*count is an upper bound on nonzero elements % Now start the Grand Loop m = 1 ; %row counter for i = 1:n-1 %Repeat Adjaceny Test above and Record Results Vi1_max = M(i,1); Vi1_min = M(i,2); Vi2_max = M(i,3); Vi2_min = M(i,4); Li = Z(i,2) - Z(i,1); % Number of i-faces Ri = Z(i,1) - 1; %Position of i-th lable row in bnd matrix for j = i+1:n %Grand j-loop %******************************** % %CHECK: % % if (i == 353)&(j == 341) % % i = i; % % end %***************************** Vj1_max = M(j,1); Vj1_min = M(j,2); Vj2_max = M(j,3); Vj2_min = M(j,4); if (Vi1_max + eps < Vj1_min)|(Vj1_max + eps < Vi1_min)|... (Vi2_max + eps < Vj2_min)|(Vj2_max + eps < Vi2_min) m = m ; % i-polygon and j-polygon not adjacent else % Check i-polygon Bij = 0 ; % Accumulated Shared-Boundary Length for (i,j) for ii = 1:Li %Check faces of i-polygon P11 = bnd(Ri + ii,1); P12 = bnd(Ri + ii,2); P1 = [P11,P12]; %First point in i-face P21 = bnd(Ri + ii + 1,1); P22 = bnd(Ri + ii + 1,2); P2 = [P21,P22]; %Second point in i-face Max1 = max([P11,P21]); Min1 = min([P11,P21]); Max2 = max([P12,P22]); Min2 = min([P12,P22]); % mbscalar(Max1); % mbscalar(Min1); % mbscalar(Max2); % mbscalar(Min2); if (Max1 + eps < Vj1_min)|(Vj1_max + eps < Min1)|... (Max2 + eps < Vj2_min)|(Vj2_max + eps < Min2) Bij = Bij ; % i-face and j-polygon not adjacent else %Check j-polygon Lj = Z(j,2) - Z(j,1); % Number of j-faces Rj = Z(j,1) - 1; %Position of j-th lable row in bnd matrix for jj = 1:Lj %check faces of j-polygon Q11 = bnd(Rj + jj,1); Q12 = bnd(Rj + jj,2); Q1 = [Q11,Q12]; %First point in j-face Q21 = bnd(Rj + jj + 1,1); Q22 = bnd(Rj + jj + 1,2); Q2 = [Q21,Q22]; %Second point in j-face MaxQ1 = max([Q11,Q21]); MinQ1 = min([Q11,Q21]); MaxQ2 = max([Q12,Q22]); MinQ2 = min([Q12,Q22]); % mbscalar(MaxQ1); % mbscalar(MinQ1); % mbscalar(MaxQ2); % mbscalar(MinQ2); if (Max1 + eps < MinQ1)|(MaxQ1 + eps < Min1)|... (Max2 + eps < MinQ2)|(MaxQ2 + eps < Min2) Bij = Bij ; % [Q1,Q2] not adjacent to [P1,P2] else % Must check for alignment of faces X = P1 - P2; Y1 = Q1 - P1; Y2 = Q2 - P1; if abs(X(1)*Y1(2)-X(2)*Y1(1)) + abs(X(1)*Y2(2)-X(2)*Y2(1)) > eps Bij = Bij ; %At least one j-point not aligned with i-face else %have alignment and must check for overlap if abs(P12 - P22) > abs(P11 - P21) %Then abs(P12 - P22) >0 lam1 = (Q12 - P22)/(P12 - P22); lam2 = (Q22 - P22)/(P12 - P22); else %must have P11 - P21 != 0 lam1 = (Q11 - P21)/(P11 - P21); lam2 = (Q21 - P21)/(P11 - P21); end %Now check relative positions using lam1 and lam2 if lam1 > 1 + eps if lam2 > 1 + eps Bij = Bij ; % no overlap elseif lam2 < eps % must have [P1,P2] in [Q1,Q2] Bij = Bij + sum((P1 - P2).^2); else %must have [P1,Q2] in [Q1,P2] S = sum((P1 - Q2).^2); % mbscalar(S); if S < eps Bij = Bij; %Do not add zero distance else Bij = Bij + sqrt(sum((P1 - Q2).^2)); end end elseif lam1 > eps if lam2 > 1 + eps %must have [P1,Q1] in [Q2,P2] S = sum((P1 - Q1).^2); % mbscalar(S); if S < eps Bij = Bij ; %Do not add zero distance else Bij = Bij + sqrt(sum((P1 - Q1).^2)); end elseif lam2 < eps %must have [Q1,P2] in [P1,Q2] S = sum((Q1 - P2).^2); % mbscalar(S); if S < eps Bij = Bij; %Do not add zero distance else Bij = Bij + sqrt(sum((Q1 - P2).^2)); end else %must have [Q1,Q2] in [P1,P2] Bij = Bij + sqrt(sum((Q1 - Q2).^2)); end else %must have lam1 < eps if lam2 < eps Bij = Bij ; %no overlap elseif lam2 > 1 + eps %must have [P1,P2] in [Q1,Q2] Bij = Bij + sqrt(sum((P1 - P2).^2)); else %must have [Q2,P2] in [P1,Q1] Bij = Bij + sqrt(sum((Q2 - P2).^2)); end end %All positions of lam1 checked end %All cases of overlap checked end %All alignment possibilities checked end %All j-faces checked end %j-polygon checked end %faces of i-polygon checked % mbscalar(Bij); if Bij > 0 W(m,1) = i; W(m,2) = j; W(m,3) = Bij; W(m+1,1) = j; %Add symmetric relation W(m+1,2) = i; W(m+1,3) = Bij; m = m + 2; else m = m; % No sharing end end %i-polygon checked end %Go to next j-polygon percent_completed = 100*(i/n) end %Go to next i-polygon %Prepare Output Data %***************************************************** %NOTE: By these conventions, we can have max(I,J) < n since % the last few polygons could be ISLANDS and hence not % be positive. So a check needs to be made coming out of % this program to readjust the size of the matrix. %****************************************************** pos = find(W(:,1) > 0 ); I = W(pos,1); J = W(pos,2); L = W(pos,3);