// key:223ee46f765a1e0eee2fdb7ca1e02850 // cvoronoi.sq // Copyright 1999, Yves Piguet. // 26.10.1999 variable x, y init (x, y) = init figure "Color Voronoi" draw drawColVoronoi(x, y) mousedrag (x, y) = dragTriangulation(x, y, _nb, _x1, _y1) mouseover _cursor = overTriangulation(_id) functions {@ function (x, y) = init n = 13; x = rand(n, 1); y = rand(n, 1); subplots('Color Voronoi'); function drawColVoronoi(x, y) scale('equal', [0,1,0,1]); tr = tri(x, y); voronoi(x, y, tr, 'mcgy'); links = [tr(:,1), tr(:,2); tr(:,2), tr(:,3); tr(:,3), tr(:,1)]; plot(x(links), y(links)); plot(x, y, '[', 1); cxy = ccircle(x, y, tr); function (x, y) = dragTriangulation(x, y, nb, x1, y1) if isempty(nb) cancel; end x(nb) = x1 + 1e-5 * rand; % rand to avoid having twice the same point y(nb) = y1 + 1e-5 * rand; function cursor = overTriangulation(id) cursor = ~isempty(id); % --- function tr = tri(x, y) % give tr = list of triangles based on x and y (indices, one triangle per row) if length(x) < 3 tr = []; else x = x - mean(x); y = y - mean(y); z = x.^2 + y.^2; c = cmb(length(x)); tr = []; for i = 1:size(c, 1) % find plane z=ab*[x;y;1] from points x(c(i,:)'), y(c(i,:)'), z(c(i,:)') ab = ([x(c(i,:)'), y(c(i,:)'), [1; 1; 1]] \ z(c(i,:)'))'; % add triangle if no other point is over plane ab if ~any(ab * [x(~c(i,:)')'; y(~c(i,:)')'; ones(1, length(x)-3)] > z(~c(i,:)')') tr(end+1, :) = find(c(i,:)); end end end function c = cmb(n) c = logical(zeros(n * (n - 1) * (n - 2) / 6, n)); l = 1; for i = 1:n-2 for j = i+1:n-1 for k = j+1:n c(l, [i, j, k]) = true; l = l + 1; end end end function voronoi(x, y, tr, colStr) % draw a voronoi diagram based on points x, filled is colStr % tr is the result of tri cxy = ccircle(x, y, tr); col = triCol(x, y, tr, colStr); tr0 = tr; for i = 1:length(x) l = find(any(tr0 == i, 2)); % list of triangles which include point i corners = tr(l, :); % check border corners = sort(corners(corners ~= i)); corners([corners(1:end-1)==corners(2:end);false] ... | [false;corners(1:end-1)==corners(2:end)]) = []; % orphews if ~isempty(corners) % add 3 dummy triangles to the end of tr and cxy (not included in the i loop) tr = [tr; ... i, corners(1), max(tr(:))+1; ... i, corners(2), max(tr(:))+2; ... i, max(tr(:))+1, max(tr(:))+2]; d = [y(i)-y(corners(1)), x(corners(1))-x(i)]; if [x(i)-mean(x),y(i)-mean(y)] * d' < 0 d = -d; % away end p1 = [x(i)+x(corners(1)), y(i)+y(corners(1))] / 2 + 2 * d / norm(d); d = [y(i)-y(corners(2)), x(corners(2))-x(i)]; if [x(i)-mean(x),y(i)-mean(y)] * d' < 0 d = -d; % away end p2 = [x(i)+x(corners(2)), y(i)+y(corners(2))] / 2 + 2 * d / norm(d); d = p1 + p2 - 2 * [x(i), y(i)]; p3 = [x(i), y(i)] + 2 * d / norm(d); cxy = [cxy; p1; p2; p3]; l = [l; size(cxy, 1)+(-2:0)']; end ix = 1; % start with one of them cix = find(tr(l(ix), :) ~= i); cix = cix(1); % index of a corner ~= point i c = cxy(l(ix), :); for j = 2:length(l) cix2 = find(tr(l(ix), :) ~= i & tr(l(ix), :) ~= tr(l(ix), cix)); % next corner trix = find(any(tr(l, :) == tr(l(ix), cix2), 2) & (1:length(l))' ~= ix); % next triangle (index in l) if length(trix) ~= 1 % none or too many (?) 'oooops!' break; % too bad (should never occur) end cix = find(tr(l(trix), :) == tr(l(ix), cix2)); ix = trix; c(end+1,:) = cxy(l(ix), :); end if isempty(col) plot(c([1:end,1], 1)', c([1:end,1], 2)'); else plot(c(:, 1)', c(:, 2)', [col(i), 'f']); end end function cxy = ccircle(x, y, tr) % center cxy=[cx,cy] of the circumscribed circle of triangles given by x,y,tr b1 = x(tr(:,1)).^2 - x(tr(:,2)).^2 + y(tr(:,1)).^2 - y(tr(:,2)).^2; b2 = x(tr(:,2)).^2 - x(tr(:,3)).^2 + y(tr(:,2)).^2 - y(tr(:,3)).^2; a11 = x(tr(:,1)) - x(tr(:,2)); a21 = x(tr(:,2)) - x(tr(:,3)); a12 = y(tr(:,1)) - y(tr(:,2)); a22 = y(tr(:,2)) - y(tr(:,3)); d = a11 .* a22 - a12 .* a21; cxy = [a22.*b1-a12.*b2, a11.*b2-a21.*b1] ./ (2 * [d,d]); function c = triCol(x, y, tr, cstr) % gives a string with color codes (taked from cstr) for each point in x,y % tr is the result of tri % special case if length(x) < 3 c = cstr(1:length(x)); return; end % build connect, where connect(i,j) is true if points i and j are connected sides = [tr(:,[1,2]); tr(:,[2,3]); tr(:,[1,3])]; connect = logical(zeros(length(x))); for i = 1:size(sides, 1) connect(sides(i, 1), sides(i, 2)) = true; end connect = connect | connect'; % build col, where col(i,j) is true if point i may have color j col = logical(ones(length(x), 4)); attempt = []; % in case of problem, can come back to attemptChoice = []; % attempt(:,end-3:end) and choose next attemptChoice(end) while any(sum(col, 2) > 1) | any(all(~col,2)) if any(all(~col,2)) % no more possible color if ~isempty(attemptChoice) % can try something else col = attempt(:,end-3:end); attempt(:,end-3:end) = []; jx = attemptChoice(end) + 1; attemptChoice(end) = []; (dmy,ix) = min(100*(sum(col, 2) <= 1) + sum(col, 2)); ixCol = find(col(ix,:)); if length(ixCol) > jx attempt = [attempt, col]; attemptChoice = [attemptChoice, jx]; end ixCol = ixCol(jx); col(ix, :) = (1:4) == ixCol; else error('Needs more than 4 colors.'); end else (dmy,ix) = min(100*(sum(col, 2) <= 1) + sum(col, 2)); ixCol = find(col(ix,:)); if length(ixCol) > 1 attempt = [attempt, col]; attemptChoice = [attemptChoice, 1]; end ixCol = ixCol(1); col(ix, :) = (1:4) == ixCol; end while any(sum(col,2)==1) ix = find(sum(col,2)==1); % points with one color ch = false; for i = 1:length(ix) ixCol = find(col(ix(i),:)); cnct = find(connect(ix(i),:)); ch = ch | any(col(cnct, ixCol)); col(cnct, ixCol) = false; end if ~ch break; end end end % take first possible color of each point for i = 1:length(x) ix = find(col(i,:)); c(i) = cstr(ix(1)); end @}