drawpoly.m

computation of conformal maps to polygonally bounded regions

function [w,beta] = drawpoly(fig,cmd)
%DRAWPOLY Draw a polygon with the mouse.
%   P = DRAWPOLY allows the user to draw a polygon with the
%   mouse.  Use the mouse to position the crosshair and press the left
%   mouse button to create a vertex.  For use with other S-C Toolbox
%   functions, the vertices must be specified in a "positively oriented"
%   manner; i.e.  counterclockwise for interior polygons and clockwise
%   for exterior regions.  There are several GUI elements added to the
%   figure to help you snap vertices to a grid, get specfic angles, etc.
%   For the last vertex, use the middle or right mouse button, or double
%   click.  Upon return, P is a polygon object.
%
%   [W,BETA] = DRAWPOLY returns vectors of vertices and turning angles
%   instead of a polygon object.
%
%   See the user's guide for full details.
%	
%   See also POLYGON methods PLOT, EDIT, MODIFY.

%   Copyright 1998 by Toby Driscoll.
%   $Id: drawpoly.m 136 2001-05-11 19:49:16Z driscoll $

persistent youve_been_warned
if isempty(youve_been_warned)
  a = questdlg('DRAWPOLY is obsolete. Use POLYEDIT instead.','Obsolete function','Run DRAWPOLY','Run POLYEDIT','Run POLYEDIT');
  youve_been_warned = 1;
  if strcmp(a,'Run POLYEDIT')
    if nargin < 1
      w = polyedit;
    else
      w = polyedit(fig);
    end
    return
  end
end
  
% If called from the GUI, reroute to subfunction cback
if nargin > 1
  [w,beta] = cback(fig,cmd);
  return
end

if nargin < 1
  fig = gcf;
end

figure(fig);
turn_off_hold = ~ishold;

% Save current properties to be restored
figprops = {'unit';'pos';'resizefcn';'pointer';'windowbuttondown';...
    'windowbuttonup';'pointer'};
figstate = get(fig,figprops);
axprops = {'unit';'pos';'xgrid';'ygrid'};
axstate = get(gca,axprops);

% Get figure position and screen size
set(fig,'unit','point')
figpos = get(fig,'pos');
unit = get(0,'unit');
set(0,'unit','point')
maxpos = get(0,'screensize');
set(0,'unit',unit)

% Set up figure size
%figsize = [462 396];
%delta = figsize - figpos(3:4);
%figpos(1) = max(min(figpos(1),maxpos(3)-figsize(1)),0); % move left if needed
%figpos(2) = max(figpos(2)-delta(2),0); % move bottom down, not top up
%figpos(3:4) = figsize;
%set(fig,'resizefcn','')
%set(fig,'pos',figpos)
set(fig,'resizefcn','drawpoly(gcbo,''resize'');')

if ~ishold
  cla
  axis auto
end

% Reset axes limits?
if strcmp(get(gca,'xlimmode'),'auto') & strcmp(get(gca,'ylimmode'),'auto')
  axlim = [-4 4 -4 4];
else
  axlim = axis;
  axis(axlim)
end

view(2)
set(gca,'unit','point','box','on','xgrid','off','ygrid','off')
set(gca,'xlim',axlim(1:2),'ylim',axlim(3:4),'plotboxaspectratio',[1,1,1])
hold on

% Set up controls
dp_make_widgets(fig);
dp_resize_widgets(fig);

set(fig, 'pointer','crosshair');
preview = line(NaN,NaN,'linestyle','--','erasemode','xor',...
    'clipping','off','tag','preview');
vertices = line(NaN,NaN,'linesty','none','marker','.',...
    'markersize',12,'tag','vertices','erasemode','xor');
line(NaN,NaN,'visible','off','tag','DRP_PT');
  
% Preparation.
if ~strcmp(computer,'SUN4')
  % Kludge.  Draw preview line when button is pressed.
  set(fig,'windowbuttondownfcn','drawpoly(gcf,''move'');');
end
set(fig,'windowbuttonupfcn', 'drawpoly(gcf,''up'');');

set(preview,'user','first','xdata',[],'ydata',[]);
set(vertices,'xdata',[],'ydata',[],'zdata',[])
  
% Get vertices.
done = 0;
while ~done
  [x0,y0] = drawpoly(fig,'getpoint');
  x = get(vertices,'xdata');
  y = get(vertices,'ydata');
  atinf = get(vertices,'zdata');
  m = length(x);
  n = m - sum(atinf)/2;
  mode = get(preview,'user');
  axlim = axis;
  reflen = mean([diff(axlim(1:2)),diff(axlim(3:4))])/60;
  if m > 0 
    if norm([x0-x(1), y0-y(1)]) < reflen
      mode = 'done';
    end
  end
  x(m+1) = x0;  
  y(m+1) = y0;
  edges = findobj(fig,'tag','edges');
  
  switch mode
  case {'normal','first'}
    if n > 0
      edges(n) = plot(x(m:m+1),y(m:m+1),'-','erasemode','back');
      set(edges(n),'user',n)
    end
    if x0>=axlim(1) & x0<=axlim(2) & y0>=axlim(3) & y0<=axlim(4)
      % Finite vertex.
      atinf(m+1) = 0;
      set(preview,'user','normal')
    else
      atinf(m+1) = 1;
      set(preview,'user','infinite')
      set(fig,'pointer','cross')
    end
  case 'infinite'
    % Re-entry point for next edge.
    atinf(m+1) = 1;
    set(fig,'pointer','crosshair')
    set(preview,'user','return')
  case 'return'
    edges(n) = plot(x(m:m+1),y(m:m+1),'-','erasemode','back');
    set(edges(n),'user',n)
    atinf(m+1) = 0;
    set(preview,'user','normal')
  case 'done'
    x = x(1:m);
    y = y(1:m);
    done = 1;
  end  
  
  % Check buttonpress.
  if ~strcmp(get(fig,'selectiontype'),'normal') & (m > 0)
    done = 1;
  end
  
  set(vertices,'xdata',x,'ydata',y,'zdata',atinf);
  set(edges,'tag','edges')
  drawnow
end
  
m = length(x);
n = length(x) - sum(atinf)/2;
if n > 0
  edges(n) = plot(x([m,1]),y([m,1]),'-');
end
drawnow

w = x(:)+i*y(:);
beta = scangle(w);
s = sum(atinf);
if s > 0
  b = beta(logical(atinf));
  b(b>0) = b(b>0) - 2;
  b = sum(reshape(b,2,s/2))';
  idx = find(atinf);
  w(idx(1:2:s)) = Inf*ones(s/2,1);
  beta(idx(1:2:s)) = b;
  w(idx(2:2:s)) = [];
  beta(idx(2:2:s)) = [];
end

% Clean up the mess.
set(fig,figprops,figstate)
set(gca,axprops,axstate)
set(vertices,'erasemode','normal')
set(edges,'erasemode','normal')
frame(1) = findobj(fig,'tag','dp_settings');
frame(2) = findobj(fig,'tag','dp_actions');
delete(get(frame(1),'userdata'))
delete(get(frame(2),'userdata'))
delete(frame)
delete(preview)
delete(vertices)
delete(edges)
  
% Construct a polygon object
if sum(beta) < 0
  poly = polygon(w,1+beta);
else
  poly = polygon(flipud(w),1-flipud(beta));
end
if nargout < 2
  w = poly;
end

axis auto
if turn_off_hold
  hold off
end
plot(poly)
set(gca,'xtickmode','auto','ytickmode','auto')
set(gca,'xticklabelmode','auto','yticklabelmode','auto')


function [out1,out2] = cback(fig,cmd)

% Invocations as a callback
if nargout > 0
  out1 = [];
  if nargout > 1
    out2 = [];
  end
end

switch cmd

case 'getpoint'
  DRP_PT = findobj(fig,'tag','DRP_PT');
  set(fig,'windowbuttonmotionfcn','drawpoly(gcf,''move'');');
  drawnow
  waitfor(DRP_PT,'user')
  data = get(DRP_PT,'user');
  out1 = data(1);
  out2 = data(2);
  
case 'move'
  ptrpos = get(gca,'currentpoint');
  P = ptrpos(1,1:2);
  axlim = axis;
  preview = findobj(fig,'tag','preview');
  vertices = findobj(fig,'tag','vertices');
  x = get(vertices,'xdata');
  y = get(vertices,'ydata');
  pts = [x(:),y(:)];
  m = size(pts,1);
  mode = get(preview,'user');
  grid = str2num(get(findobj(fig,'tag','grid_space'),'string'));
  grid = grid * get(findobj(fig,'tag','snap'),'value');
  qang = 1/str2num(get(findobj(fig,'tag','ang_quant'),'string'));
  qang = qang * get(findobj(fig,'tag','qang'),'value');
  qlen = str2num(get(findobj(fig,'tag','len_quant'),'string'));
  qlen = qlen * get(findobj(fig,'tag','qlen'),'value');
  
  % Modify point to meet mode constraints.
  if strcmp(mode,'normal')
    % No constraints.
    mode = 0;
  elseif strcmp(mode,'first')
    % 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));
    mode = 0;
  elseif strcmp(mode,'infinite')
    % Point may not be inside axes box.
    if all(P>axlim([1,3])) & all(P<axlim([2,4]))
      [junk,j] = min(abs([P(1)-axlim(1:2);P(2)-axlim(3:4)]'));
      P = axlim(j);
    end
    mode = 1;
  elseif strcmp(mode,'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]);