polyedit.m

computation of conformal maps to polygonally bounded regions

  case 'return'
    % Point may not be outside axes box.
    P(1) = min(max(P(1),axlim(1)),axlim(2));
    P(2) = min(max(P(2),axlim(3)),axlim(4));
    ang = scangle([pts(m-2:m,1);P(1)]+i*[pts(m-2:m,2);P(2)]);
    ang = ang(2:3);
    ang(ang>0) = ang(ang>0) - 2;
    ang = sum(ang);
    mode = 2;
end

% Modify point to meet angle, length, or grid constraints.
if any(mode==[0,2]) & qang & (m > 0)	% quantized angle
  % Find arg of new side which meets quantization requirements. 
  if m==1
    % No reference side, so refer to positive real axis
    ang = angle(P(1)-pts(m,1)+i*(P(2)-pts(m,2)))/pi;
    theta = qang*round(ang/qang)*pi;
  elseif mode==0
    ang = scangle([pts(m-1:m,1);P(1)]+i*[pts(m-1:m,2);P(2)]);
    ang = qang*round(ang(2)/qang);
    theta = atan2(pts(m,2)-pts(m-1,2),pts(m,1)-pts(m-1,1))-pi*ang;
  elseif mode==2
    % ang was computed above
    ang = qang*round(ang/qang);
    theta = atan2(pts(m-1,2)-pts(m-2,2),pts(m-1,1)-pts(m-2,1))-pi*ang;
  end
  % Project P to correct angle.
  A = pts(m,:);
  BA = [cos(theta),sin(theta)];
  P = A + ((BA)*(P-A)')*(BA);
  grid = 0;
end
if (mode==0) & qlen & (m > 0)		% quantized length
  A = pts(m,:);
  len = norm(P-A);
  fixlen = qlen*(round(len/qlen));
  P = A + fixlen/len*(P-A);
  grid = 0;
end
if any(mode==[0,1,2]) & grid	        % snap to grid
  grid = [grid,grid];
  minxy = axlim([1,3]);
  P = minxy + grid.*(round((P-minxy)./grid));
end

% If returning from infinity, don't allow angle > -1.
if mode==2
  if ang > -1+2*eps
    % Would be illegal.  Project to make ang=-1.
    A = pts(m,:);
    B = pts(m-2,:) + A - pts(m-1,:);
    P = A + ((B-A)*(P-A)')/((B-A)*(B-A)')*(B-A);
    ang = -1;
    qang = 0;				% override other restrictions
    grid = 0;
  elseif ang < -3
    % It's illegal.  Is it even possible?
    P = [NaN,NaN];
    ang = NaN;
  end
end

% Update.
if m > 0 & (mode~=1)
  set(data.PreviewLine,'xdata',[pts(m,1),P(1)], 'ydata',[pts(m,2),P(2)]);
end
data.currentpoint = P;
guidata(obj,data)
drawnow
  
function PEAddFinish(obj,varargin)
data = guidata(obj);
P = data.currentpoint;
if ~isnan(P(1))
  set(data.Figure,'windowbuttonmotion','')
  set(data.Figure,'windowbuttonup','')
  set(data.PreviewLine,'xdata',[P(1) NaN],'ydata',[P(2) NaN])
  pts = PE_get_clicks(data);
  x = pts(:,1);  y = pts(:,2);
  %%x = get(data.Vertices,'xdata');
  %%y = get(data.Vertices,'ydata');
  m = length(x);
  n = m - sum(data.atinf)/2;
  axlim = axis;
  reflen = mean([diff(axlim(1:2)),diff(axlim(3:4))])/60;
  if (m > 0) & norm([P(1)-x(1), P(2)-y(1)]) < reflen
    %%x = x(1:m);
    %%y = y(1:m);
    %%set(data.Vertices,'xdata',x,'ydata',y)
    guidata(obj,data)
    PEFinish(obj)
    return
  end
  x(m+1,1) = P(1);  
  y(m+1,1) = P(2);
  
  % Beta (angle) computation.
  if n==1
    data.orient = sign( diff( x(1:2)+i*y(1:2) ) );
  elseif n > 1 
    neworient = sign( diff( x(m:m+1)+i*y(m:m+1) ) );
    if strcmp(data.addmode,'normal')
      % Last vertex was finite.
      data.polybeta(n) = angle(data.orient/neworient) / pi;
    elseif strcmp(data.addmode,'return')
      % Last vertex was infinite.
      b = scangle( x(m-2:m+1)+i*y(m-2:m+1) );
      b = b(2:3);
      b(b>0) = b(b>0) - 2;
      data.polybeta(n) = sum(b);
    end
    data.orient = neworient;
  end
 
  switch data.addmode
    case {'normal','first'}
      if P(1)>=axlim(1) & P(1)<=axlim(2) & P(2)>=axlim(3) & P(2)<=axlim(4)
	% Finite vertex.
	data.atinf(m+1) = 0;
	data.addmode = 'normal';
	data.polyvertex(n+1) = P(1)+i*P(2);
	data.polybeta(n+1) = NaN;
	data.Vertices(n+1) = plot(P(1),P(2),'o','user',[P(1) P(2)]);
      else
	data.atinf(m+1) = 1;
	data.addmode = 'infinite';
	% Set up for re-entry point next time.
	set(data.Figure,'pointer','cross')
	h = [data.InsertButton,data.MoveButton,data.DeleteButton,...
	      data.FinishButton];
	set(h,'enable','off')
	data.polyvertex(n+1) = Inf;
	data.polybeta(n+1) = NaN;
	data.Vertices(n+1) = plot(Inf,Inf,'o','user',[P(1) P(2)]);
      end
      if n > 0
	data.Edges(n) = plot(x(m:m+1),y(m:m+1),'-','erasemode','norm');
	set(data.Edges(n),'user',n)
      end
    case 'infinite'
      data.atinf(m+1) = 1;
      set(data.Figure,'pointer','crosshair')
      data.addmode = 'return';
      pt = get(data.Vertices(n+0.5),'user');
      set(data.Vertices(n+0.5),'user',[pt;[P(1) P(2)]])
    case 'return'
      data.Edges(n) = plot(x(m:m+1),y(m:m+1),'-','erasemode','norm');
      set(data.Edges(n),'user',n)
      data.atinf(m+1) = 0;
      data.polyvertex(n+1) = P(1)+i*P(2);
      data.polybeta(n+1) = NaN;
      data.Vertices(n+1) = plot(P(1),P(2),'o','user',[P(1) P(2)]);
      % Restore mode buttons.
      h = [data.InsertButton,data.MoveButton,data.DeleteButton,...
	    data.FinishButton];
      set(h,'enable','on')
      data.addmode = 'normal';
  end     
  
  set(data.Vertices(1:end-1),'markerfacecolor','b','color','b')
  set(data.Vertices(end),'markerfacecolor','r','color','r')
  data.currentpoint = [NaN NaN];
  guidata(obj,data)
  drawnow
  
end

function PEMoveStart(obj,varargin)
data = guidata(obj);
% Figure out which vertex is selected.
%%P = get(gca,'currentpoint');
%%[m,idx] = min( abs( P(1)+i*P(2) - data.polyvertex ) );
%%data.moveselected = idx(1);
idx = find( get(gcf,'currentobj')==data.Vertices );
data.moveselected = idx;
adj = PE_adjacent_edges(data,data.moveselected);
data.moveadj = adj;
guidata(obj,data)
set(data.Edges(adj(~isnan(adj))),'linest','--')
set(data.Figure,'windowbuttonmotion',@PEMoveMove)
set(data.Figure,'windowbuttonup',@PEMoveFinish)

function PEMoveMove(obj,varargin)
data = guidata(obj);
xy = get(gca,'currentpoint');
xy = xy(1,1:2);
w = data.polyvertex;
n = length(w);
% Keep inside axes box.
axlim = axis;
xy(1) = min( max(xy(1),axlim(1)), axlim(2) );
xy(2) = min( max(xy(2),axlim(3)), axlim(4) );
% Snap to grid if requested.
if get(data.SnapToggle,'value');
  grid = str2num(get(data.GridSpaceEdit,'string'))*[1 1];
  minxy = axlim([1 3]);
  xy = minxy + grid.*(round((xy-minxy)./grid));
end
k = data.moveselected;
set(data.Vertices(k),'xdata',xy(1),'ydata',xy(2),'user',xy)

adj = data.moveadj;
edges = data.Edges;
beta = data.polybeta;

scale = max( axlim(2)-axlim(1), axlim(4)-axlim(3) );
if ~isnan(adj(1))
  % Adjust the "predecessor" side.
  j = rem(k-2+n,n)+1;
  if isinf(w(j))
    xd = get(edges(j),'xd');
    yd = get(edges(j),'yd');
    phi = atan2(-diff(yd),-diff(xd));
    % Note: draw lines to infinity long enough
    r = 2*scale;
    %r = sqrt(diff(xd)^2+diff(yd)^2);
    y = xy(1)+i*xy(2) + [r*exp(i*phi),0];
    % The "clickpoint" for the inf vertex has changed.
    pts = get(data.Vertices(j),'user');
    pts(2,:) = [real(y(1)) imag(y(1))];
    set(data.Vertices(j),'user',pts)
  else
    y = [w(j),xy(1)+i*xy(2)];
    phi = angle(diff(y)/diff(w([j,k])));
    beta(k) = beta(k)+phi/pi;
    beta(j) = beta(j)-phi/pi;
  end
  set(edges(adj(1)),'xd',real(y),'yd',imag(y))
end

if ~isnan(adj(2))
  % Adjust the "successor" side.
  j = rem(k,n)+1;
  if isinf(w(j))
    xd = get(edges(k),'xd');
    yd = get(edges(k),'yd');
    phi = atan2(diff(yd),diff(xd));
    r = 2*scale;
    %%r = sqrt(diff(xd)^2+diff(yd)^2);
    y = xy(1)+i*xy(2) + [0,r*exp(i*phi)];
    % The "clickpoint" for the inf vertex has changed.
    pts = get(data.Vertices(j),'user');
    pts(1,:) = [real(y(2)) imag(y(2))];
    set(data.Vertices(j),'user',pts)
  else
    y = [xy(1)+i*xy(2),w(j)];
    phi = angle(-diff(y)/diff(w([j,k])));
    beta(k) = beta(k)-phi/pi;
    beta(j) = beta(j)+phi/pi;
  end
  set(edges(adj(2)),'xd',real(y),'yd',imag(y))
end
% Make changes effective
drawnow
data.polyvertex(k) = xy(1)+i*xy(2);
data.polybeta = beta;
guidata(obj,data)

function PEMoveFinish(obj,varargin)
data = guidata(obj);
set(data.Figure,'windowbuttonmotion','')
set(data.Figure,'windowbuttonup','')
set(data.Edges,'linesty','-')
% Update the orientation of the last side, if polygon is not closed.
if ~data.isclosed
  pts = PE_get_clicks(data);
  w = pts(:,1) + 1i*pts(:,2);
  data.orient = sign( diff( w(end-1:end) ) );
  guidata(obj,data)
end
    

function PEInsert(obj,varargin)
data = guidata(obj);
idx = find( get(gcf,'currentobj')==data.Edges );
[wn,bn] = scaddvtx(data.polyvertex,data.polybeta,idx,axis);
P = wn(idx+1);
data.polyvertex = wn;
data.polybeta = bn;
%%new = zeros(1,length(data.atinf)+1);
%%new([1:idx idx+2:end]) = data.atinf(:);
%%data.atinf = new;
data.atinf = [data.atinf(1:idx); 0; data.atinf(idx+1:end)];
newVertex = plot(real(P),imag(P),'bo',...
    'markerfacecolor','b','user',[real(P) imag(P)]);
data.Vertices = [data.Vertices(1:idx) newVertex data.Vertices(idx+1:end)];
pred = data.Edges(idx);
x = get(pred,'xdata');  y = get(pred,'ydata');
set(pred,'xdata',[x(1) real(P)],'ydata',[y(1) imag(P)])
succ = plot([real(P) x(2)],[imag(P) y(2)],'-','erase','norm');
data.Edges = [data.Edges(1:idx) succ data.Edges(idx+1:end)];
guidata(obj,data)


function PEDelete(obj,varargin)
data = guidata(obj);
k = find( get(gcf,'currentobj')==data.Vertices );
n = length(data.polyvertex);
if n==1
  PEClear(obj)
  return
end
% Special case for an unfinished polygon at an endpoint.
if ~data.isclosed & ( k==1 | k==n ) 
  if (k==1 & isinf(data.polyvertex(2))) | ...
	(k==n & isinf(data.polyvertex(n-1))) 
    errordlg('To be deleted, a vertex must have at least one finite neighbor.','PolyEdit Error')
    return
  end
  delete(data.Vertices(k))
  data.Vertices(k) = [];
  data.polyvertex(k) = [];
  data.polybeta(k) = [];
  data.atinf(k) = [];
  if k==1
    data.polybeta(1) = NaN;
    delete(data.Edges(1))
    data.Edges(1) = [];
  else
    data.orient = data.orient*exp(i*pi*data.polybeta(end));
    data.polybeta(end) = NaN;
    delete(data.Edges(n-1))
    data.Edges(n-1) = [];
    set(data.Vertices(n-1),'color','r','markerfacecolor','r')
  end
else
  w = data.polyvertex;
  beta = data.polybeta;
  % Neighbors
  pred = rem( k-2+n, n ) + 1;
  succ = rem(k,n) + 1;
  % There must be at least one finite neighbor.
  infp = isinf(w(pred));
  infs = isinf(w(succ));
  if ~infp & ~infs
    wp = w(pred);
    ws = w(succ);
    beta(pred) = angle( exp(i*pi*beta(pred))*(w(k)-wp)/(ws-wp) ) / pi;
    beta(succ) = angle( exp(i*pi*beta(succ))*(wp-ws)/(w(k)-ws) ) / pi;
    side = [wp ws];
  elseif infp & ~infs
    axlim = axis;
    scale = 2*max( axlim(2)-axlim(1), axlim(4)-axlim(3) );
    beta(pred) = beta(pred) + beta(k);
    side = w(succ) + [scale*sign(w(k)-w(succ)) 0];
    pts = get(data.Vertices(pred),'user');
    pts(2,:) = [real(side(1)) imag(side(1))];
    set(data.Vertices(pred),'user',pts)
  elseif ~infp & infs
    axlim = axis;
    scale = 2*max( axlim(2)-axlim(1), axlim(4)-axlim(3) );
    beta(succ) = beta(succ) + beta(k);
    side = w(pred) + [0 scale*sign(w(k)-w(pred))];
    pts = get(data.Vertices(succ),'user');
    pts(1,:) = [real(side(2)) imag(side(2))];
    set(data.Vertices(succ),'user',pts)
  else
    errordlg('To be deleted, a vertex must have at least one finite neighbor.','PolyEdit Error')
    return
  end

  data.polyvertex(k) = [];
  data.atinf(k) = [];
  data.polybeta = beta;
  data.polybeta(k) = [];
  delete(data.Vertices(k))
  data.Vertices(k) = [];
  set(data.Edges(pred),'xdata',real(side),'ydata',imag(side))