function xy = cyclops(a, xy0, wts)
% Place nodes in radial pattern
xy = cyclopsR(a, xy0);
dxy = round(wts*(xy-xy0)./sum(wts));
xy = round(bsxfun(@minus,xy,dxy));
end
function xy = cyclopsR(a, xy0)
nn = length(a);
deg = sum(a);
la = logical(a);
xy = zeros(size(xy0));
left = true(1,nn);
fill = zeros(1,nn);
seq = zeros(1,nn);
ind = 0;
ap = zeros(1,nn);
ap1 = zeros(1,nn);
ap2 = zeros(1,nn);
nr = zeros(1,nn);
misc = 0.0001;
% highest degree node is root
rt = find(deg==max(deg),1);
left(rt) = false;
% add level 1 nodes radially
lvl = find(a(rt,:));
lvl = pickOrder(lvl,la(lvl,lvl));
ring = ceil(deg(rt)/8);
ringn = 8*ring;
spot = round(linspace(0,ringn-1,numel(lvl)));
ap(lvl) = spot/ringn;
r = ring./cos(mod(2*pi*ap(lvl)+pi/4,pi/2)-pi/4);
nr(lvl) = ring;
ap1(lvl) = (ap(lvl)+[ap(lvl(end))-1 ap(lvl(1:end-1))])/2+misc;
ap2(lvl) = (ap(lvl)+[ap(lvl(2:end)) ap(lvl(1))+1])/2-misc;
xy(lvl,:) = round([r.*cos(2*pi*ap(lvl));r.*sin(2*pi*ap(lvl))])';
left(lvl) = false;
%gplot(a(~left,~left),xy(~left,:),'o-');
%axis equal
%waitforbuttonpress;
% add remaining nodes level by level
lvlb = lvl;
lvl = find(left&any(a(lvl,:),1));
while any(lvl)
% place children of previous level
for par = lvlb
ch = find(left&a(par,:));
if isempty(ch), continue; end;
ch = pickOrder(ch,la(ch,ch));
nch = numel(ch);
ring = nr(par)+1;
while (floor(8*ring*ap2(par))-ceil(8*ring*ap1(par))+1 < nch)
ring = ring+1;
end;
ringn = 8*ring;
nr(ch) = ring;
ap(ch) = round(linspace(ceil(8*ring*ap1(par)),floor(8*ring*ap2(par)),nch))/(8*ring);
ap1(ch) = [ap1(par) (ap(ch(1:end-1))+ap(ch(2:end)))/2]+misc*(ap2(par)-ap1(par));
ap2(ch) = [ap1(ch(2:end)) ap2(par)]-misc*(ap2(par)-ap1(par));
xy(ch,:) = ringPts(ring,ap(ch));
left(ch) = false;
%clf
%gplot(a(~left,~left),xy(~left,:),'o-');
%hold on
%gplot(a([par ch],[par ch]),xy([par ch],:),'ro-');
%hold off
%axis equal
%waitforbuttonpress;
end;
% prepare for next level
lvlb = lvl;
lvl = find(left&any(a(lvl,:),1));
end;
%gplot(a,xy,'o-');
%axis equal
if any(left)
rxy = cyclopsR(a(left,left), xy0(left,:));
rxy(:,1) = rxy(:,1)-min(rxy(:,1))+max(xy(~left,:))+1;
xy(left,:) = rxy;
end;
end
function xy = ringPts(ring,ap)
r = ring./cos(mod(2*pi*ap+pi/4,pi/2)-pi/4);
s = sin(2*pi*ap);
c = cos(2*pi*ap);
d = abs(mod(ap*8*ring+ring,2*ring)-ring);
%d = sqrt(r.^2-ring.^2);
xy = [ring*sign(c)' ring*sign(s)'];
sd = abs(s)>abs(c);
xy(sd,1) = d(sd).*sign(c(sd));
xy(~sd,2) = d(~sd).*sign(s(~sd));
end
function id = pickOrder(id0,la)
nid = numel(id0);
id = zeros(1,nid);
nc = sum(la);
nf = zeros(1,nid);
left = true(1,nid);
for i = 1:nid
s = nf-nc;
j = find(left&(s==max(s(left))),1);
id(i) = id0(j);
nf(la(j,:)) = nf(la(j,:))+1;
nc(la(j,:)) = nc(la(j,:))-1;
left(j) = false;
end;
end
function N = numKnots(XY,A)
% This function finds the number of knots in a network of points XY whose
% connectivity is given by A.
% Consider only unique lines (A is symmetric with ones on diagonal)
A = triu(A,1);
nLines = nnz(A);
% Find xy coordinates for first and second point on each line
[p1i,p2i] = find(A);
x1 = repmat(XY(p1i,1),1,nLines);
y1 = repmat(XY(p1i,2),1,nLines);
x2 = repmat(XY(p2i,1),1,nLines);
y2 = repmat(XY(p2i,2),1,nLines);
x3 = repmat(XY(p1i,1)',nLines,1);
y3 = repmat(XY(p1i,2)',nLines,1);
x4 = repmat(XY(p2i,1)',nLines,1);
y4 = repmat(XY(p2i,2)',nLines,1);
pick = tril(true(nLines),-1);
x1 = x1(pick);
y1 = y1(pick);
x2 = x2(pick);
y2 = y2(pick);
x3 = x3(pick);
y3 = y3(pick);
x4 = x4(pick);
y4 = y4(pick);
N = sum(areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4));
end
function bool = areIntersecting(x1,y1,x2,y2,x3,y3,x4,y4)
% Determine if two line segments intersect.
%
% line1 = [x1, y1, x2, y2]
% line2 = [x3, y3, x4, y4]
%
% This function is robust to floating point precision issues assuming all
% input coordinates are integer-valued. If all coordinates have absolute
% values that do not exceed x, then the arithmetic operations used in this
% function could, as an intermediate step, produce a value with magnitude
% as large as 8*x^4. For double-precision floating point, this implies
% that the following must be satisfied:
%
% 8*x^4 < 2^53
%
% or equivalently, x < 5792.
%
% http://en.wikipedia.org/wiki/Line-line_intersection
bool = false(1,numel(x1));
gen = true(numel(x1),1);
% Test to see if the lines share exactly one endpoint. If so, the lines
% intersect if either of the free points lies on the line segment formed by
% the other two points.
s13 = find((x1==x3)&(y1==y3));
%s13 = find((x1==x3)&(y1==y3)&((x2~=x4)|(y2~=y4)));
bool(s13) = (isPointOnSegment(x2(s13), y2(s13), x1(s13), y1(s13), x4(s13), y4(s13))|isPointOnSegment(x4(s13), y4(s13), x1(s13), y1(s13), x2(s13), y2(s13)));
gen(s13) = false;
s14 = find((x1==x4)&(y1==y4));
%s14 = find((x1==x4)&(y1==y4)&((x2~=x3)|(y2~=y3)));
bool(s14) = (isPointOnSegment(x2(s14), y2(s14), x1(s14), y1(s14), x3(s14), y3(s14))|isPointOnSegment(x3(s14), y3(s14), x1(s14), y1(s14), x2(s14), y2(s14)));
gen(s14) = false;
s23 = find((x2==x3)&(y2==y3));
%s23 = find((x2==x3)&(y2==y3)&((x1~=x4)|(y1~=y4)));
bool(s23) = (isPointOnSegment(x1(s23), y1(s23), x2(s23), y2(s23), x4(s23), y4(s23))|isPointOnSegment(x4(s23), y4(s23), x1(s23), y1(s23), x2(s23), y2(s23)));
gen(s23) = false;
s24 = find((x2==x4)&(y2==y4));
%s24 = find((x2==x4)&(y2==y4)&((x1~=x3)|(y1~=y3)));
bool(s24) = (isPointOnSegment(x1(s24), y1(s24), x2(s24), y2(s24), x3(s24), y3(s24))|isPointOnSegment(x3(s24), y3(s24), x1(s24), y1(s24), x2(s24), y2(s24)));
gen(s24) = false;
% Next we check for parallel and coincident lines. Parallel lines
% obviously don't intersect. Coincident lines intersect if an endpoint from
% one of the lines lies on the other segment.
ss = find(gen&haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4));
gen(ss) = false;
si = ss(haveSameIntercept(x1(ss), y1(ss), x2(ss), y2(ss), x3(ss), y3(ss), x4(ss), y4(ss)));
bool(si) = ((isPointBetween(x3(si), x4(si), x1(si), 1) & isPointBetween(y3(si), y4(si), y1(si), 1)) | ...
(isPointBetween(x3(si), x4(si), x2(si), 1) & isPointBetween(y3(si), y4(si), y2(si), 1)) | ...
(isPointBetween(x1(si), x2(si), x3(si), 1) & isPointBetween(y1(si), y2(si), y3(si), 1)) | ...
(isPointBetween(x1(si), x2(si), x4(si), 1) & isPointBetween(y1(si), y2(si), y4(si), 1)));
% If we get this far, we're in the general case of two well-defined
% non-parallel lines. The lines formed by the two segments must intersect
% somewhere. Find this intersection point and determine if it lies on the
% segments.
gen = find(gen);
Px_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(x3(gen) - x4(gen)) - (x1(gen) - x2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Px_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));
Py_n = (x1(gen).*y2(gen) - y1(gen).*x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen).*y4(gen) - y3(gen).*x4(gen));
Py_d = (x1(gen) - x2(gen)).*(y3(gen) - y4(gen)) - (y1(gen) - y2(gen)).*(x3(gen) - x4(gen));
bool(gen) = (isPointBetween(x1(gen), x2(gen), Px_n, Px_d) & ...
isPointBetween(y1(gen), y2(gen), Py_n, Py_d) & ...
isPointBetween(x3(gen), x4(gen), Px_n, Px_d) & ...
isPointBetween(y3(gen), y4(gen), Py_n, Py_d));
end
function bool = haveSameIntercept(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (x4 - x3).*(y1.*(x2 - x1) - x1.*(y2 - y1)) == ...
(x2 - x1).*(y3.*(x4 - x3) - x3.*(y4 - y3));
end
function bool = haveSameSlope(x1, y1, x2, y2, x3, y3, x4, y4)
bool = (y2 - y1).*(x4 - x3) == (x2 - x1).*(y4 - y3);
end
function bool = isPointOnSegment(x1, y1, x2, y2, x3, y3)
% Determine if point (x1, y1) is on the segment formed by (x2, y2) and (x3, y3).
bool = haveSameSlope(x2, y2, x1, y1, x1, y1, x3, y3) & ...
isPointBetween(x2, x3, x1, 1) & ...
isPointBetween(y2, y3, y1, 1);
end
function bool = isPointBetween(pt1, pt2, num, den)
% Determine if either of the following is satisfied:
% pt1 <= num/den <= pt2
% pt1 >= num/den >= pt2
bool = ((pt1.*den <= num) & (num <= pt2.*den)) | ...
((pt1.*den >= num) & (num >= pt2.*den));
end
|
