function t_puzzle ( )

%*****************************************************************************80
%
%% t_puzzle() sets up the T puzzle.
%
%  Discussion:
%
%    Rearrange the four pieces to form a T, an arrow, a rhombus, or a fat T.
%
%    Click near the center of one of the four pieces to drag it.
%    Click near a vertex of a piece to rotate it.
%    Right click or control click on a piece to change its orientation.
%    Click on the subplot to cycle through the puzzle shapes.
%
%  Licensing:
%
%    Copyright 2014 The MathWorks, Inc.
%
%  Modified:
%
%    10 December 2021
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Martin Gardner,
%    Mathematical Games:
%    Advertising premiums to beguile the mind: classics by Sam Loyd,
%    master puzzle poser,
%    Scientific American,
%    Volume 225, Number 5, pages 114-121, November 1971.
%
%    Cleve Moler,
%    Experiments with MATLAB,
%    ebook: https://www.mathworks.com/moler/exm/index.html
%

%
%  Define the four playing pieces.
%
  s = sqrt(2);
  T{1} = [0 1 1+2*i 3*i 0];
  T{2} = [0 s/2+s/2*i 1+(s-1)*i 1+(2*s-1)*i 2*s*i 0];
  T{3} = [0 1 i 0];
  T{4} = [0 1 1+(2-s)*i (3-s)*i 0];
%
%  Define the four puzzle shapes for subplot.
%
  P{1} = [2 2 1 1 0 0 -1 -1 2]-1/2 + [2 1 1 -2 -2 1 1 2 2]*i;
  P{2} = [0 2 2-s 2-s 1-s 1-s -2 0] + [2 0 0 -2 -2 0 0 2]*i;
  P{3} = [1 3 0 -2 1]-1/2 + [1 -1 -1 1 1]*i;
  P{4} = [0 2 1 1 -1 -1 -2 0]+[2 0 0 -1 -1 0 0 2]*i;

  init_plot ( T );
  init_subplot ( P );
  init_buttons ( );
%
%  In order for "P" to be accessible by button_down(), the following
%  functions must be included in the body of the main function.
%
%  Hence, the "end" statement for this main function comes at the 
%  very end of this file!
%
%  There must be a better way.
%
%  return
%end
function button_motion ( varargin )

%*****************************************************************************80
%
%% button_motion() drags or rotates a piece, using complex arithmetic.
%      
  point = get(gca,'currentpoint');
  v = point(1,1) + point(1,2)*i;
  hit = get(gca,'userdata');
  if ( ~isempty(hit) & ishandle(hit) & hit ~= 1 )
    u = get(hit,'xdata') + get(hit,'ydata')*i;
    z = get(hit,'userdata');
%
%  Check if closer to center or vertex.
%
    w  = mean(u(1:end-1));
%
%  Drag.
%
    if ( abs(w-v) < min(abs(u-v)) )
      u = u - (z - v);
%
%  Rotate about the center by an integer
%
    else
%          
%  multiple of pi/20, which is 9 degrees.
%
      dtheta = pi/20;
      theta = angle(v-w) - angle(z-w);
      theta = round(theta/dtheta)*dtheta;
      omega = exp(i*theta);
      u = omega*(u-w) + w;
      v = omega*(z-w) + w;
    end

    set(hit,'xdata',real(u),'ydata',imag(u),'userdata',v);

  end

  return
end
function init_plot ( T )

%*****************************************************************************80
%
%% init_plot() initializes the primary plot.
%
  clf ( );
  shg ( );
  set(gcf,'numbertitle','off','menubar','none','name','T_puzzle');
  ax(1) = axes('pos',[0 0 1 1],'color',get(gcf,'color'), ...
         'xlim',[-4 4],'ylim',[-3 3],'userdata',[]);
  hold ( 'on' );
  axis ( 'off' );
  bluegreen = [0 1/2 1/2];
  init_shift = [-2-i -1/2-i 1-i 1+i/4];
  for k = 1 : 4
    t = T{k} + init_shift(k);
    fill(real(t),imag(t),bluegreen, ...
      'markeredgecolor','black','linewidth',2);
  end

  set(gcf,'userdata',ax);

  return
end  
function init_subplot ( P )

%*****************************************************************************80
%
%% init_subplot() initializes the puzzle shapes subplot.
%
  ax = get(gcf,'userdata');
  ax(2) = axes('units','normal','position',[.02 .02 .24 .24]);
  shape = 1;
  t = P{shape};
  buttoncolor = [.8314 .8157 .7843];
  bluegreen = [0 1/2 1/2];
  fill(real(t),imag(t),bluegreen, ...
         'markeredgecolor','black','linewidth',1);
  set(ax(2),'xtick',[],'ytick',[],'xlim',[-4 4],'ylim',[-3 3], ...
         'color',buttoncolor,'userdata',shape);
  set(gcf,'userdata',ax);

  return
end   
function init_buttons ( )

%*****************************************************************************80
%
%% init_buttons() initializes the buttons and callbacks.
%
  ax = get(gcf,'userdata');
  axis(ax(1));

  set(gcf,'windowbuttondownfcn',@button_down, ...
          'windowbuttonmotionfcn',@button_motion, ...
          'windowbuttonupfcn',@button_up);

  uicontrol('units','normal','position',[.83 .12 .12 .05],'style','toggle', ...
            'callback','t_puzzle','string','restart');

  uicontrol('units','normal','position',[.83 .06 .12 .05],'style','push', ...
            'callback','close(gcf)','string','exit');

  return
end
function button_down ( varargin )

%*****************************************************************************80
%
%% button_down() handles a button that has been pressed.
%
  ax = get(gcf,'userdata');

  if gca == ax(1)
    select_piece ( );
  else
    cycle_subplot ( P );
  end

  return
end  
function select_piece ( )

%*****************************************************************************80
%
%% select_piece() selects a piece and remembers it in ax(1) userdata.
%
  point = get(gca,'currentpoint');
  z = point(1,1) + point(1,2)*i;
  delete(findobj(gca,'type','line'));
  h = flipud(get(gca,'children'));
  h = h(1:4);

  hit  = [];
  for k = 1:length(h)
    x = get(h(k),'xdata');
    y = get(h(k),'ydata');
    if inregion(real(z),imag(z),x,y)
      hit = h(k);
      set(hit,'userdata',z);
      break
    end
  end

  set(gca,'userdata',hit);
%
%  Right click reflects piece about center horizontally.
%
  if ~isempty(hit) & isequal(get(gcf,'selectiontype'),'alt')
    x = 2*mean(x(1:end-1))-x;
    set(hit,'xdata',x);
  end

  return
end
function cycle_subplot ( P )

%*****************************************************************************80
%
%% cycle_subplot() cycles through puzzle shapes in subplot.
%
%  Discussion:
%
%    A fourth shape was added.
%
%  Modified:
%
%    10 December 2021
%
  ax = get(gcf,'userdata');
  shape = get(ax(2),'userdata');
  shape = mod(shape,4)+1;
  t = P{shape};
  f = get(ax(2),'child');
  set(f,'xdata',real(t),'ydata',imag(t));
  set(ax(1),'userdata',[]);
  set(ax(2),'userdata',shape);
  axes(ax(1));

  return
end  
function button_up ( varargin )

%*****************************************************************************80
%
%% button_up() snaps to any nearby piece and deselects.
%  
  delete(findobj(gca,'type','line'));
  hit = get(gca,'userdata');
  set(gca,'userdata',[]);
%     
%  Compute distance to nearest vertex of other pieces.
%  
  z = get(hit,'xdata') + get(hit,'ydata')*i;
  h = get(gca,'children');

  w = [];
  for k = 1:length(h)
    if h(k) ~= hit
      w = [w; get(h(k),'xdata')+get(h(k),'ydata')*i];
    end
  end

  dz = 1;
  for k = 1:length(z);
    d = z(k) - w;
    dw = d(find(abs(d)==min(abs(d)),1));
    if abs(dw) < abs(dz)
      dz = dw;
    end
  end
% 
%  If close enough, snap to nearby piece.
% 
  tol = 1/8;
  if abs(dz) < tol
    set(hit,'xdata',real(z-dz),'ydata',imag(z-dz));
    bluegreen = [0 1/2 1/2];
    set(hit,'facecolor',1.25*bluegreen);
    pause(.25);
    set(hit,'facecolor',bluegreen);
  end

  return
end
function [ in, on ] = inregion ( x, y, xv, yv )

%*****************************************************************************80
%
%% inregion() is true for points inside or on a polygonal region.
%
%  Discussion:
%
%    IN = INREGION(X, Y, XV, YV) returns a matrix IN the size of X and Y.
%    IN(p,q) = 1 if the point (X(p,q), Y(p,q)) is either strictly inside or
%    on the edge of the polygonal region whose vertices are specified by the
%    vectors XV and YV; otherwise IN(p,q) = 0.
%
%    [IN ON] = INREGION returns a second matrix, ON, which is the size of X
%    and Y.  ON(p,q) = 1 if the point (X(p,q), Y(p,q)) is on the edge of the
%    polygonal region; otherwise ON(p,q) = 0.
%
%    INREGION is a modification of INPOLYGON that uses a roundoff error
%    compensating tolerance in the cross product sign test.
%
%  Example:
%
%    xv = [-3 -3 1 1 3 1 -1 -1 -3];
%    yv = [-3 -1 3 1 1 -1 -1 -3 -3];
%    [x,y] = meshgrid(-3:1/2:3);
%    [in,on] = inregion(x,y,xv,yv);
%    p = find(in-on);
%    q = find(on);
%    plot(xv,yv,'-',x(p),y(p),'ko',x(q),y(q),'ro')
%    axis([-5 5 -4 4])
%
%  Licensing:
%
%   Copyright 2014 Cleve Moler
%   Copyright 2014 The MathWorks, Inc.
%
%  Modified:
%
%    09 December 2021
%
%  Author:
%
%    Cleve Moler
%
%  Reference:
%
%    Cleve Moler,
%    Experiments with MATLAB,
%    ebook: https://www.mathworks.com/moler/exm/index.html
%

%
% If (xv,yv) is not closed, close it.
%
  xv = xv(:);
  yv = yv(:);
  Nv = length(xv);
  if ((xv(1) ~= xv(Nv)) | (yv(1) ~= yv(Nv)))
    xv = [xv ; xv(1)];
    yv = [yv ; yv(1)];
    Nv = Nv + 1;
  end

  inputSize = size(x);

  x = x(:).';
  y = y(:).';

  mask = (x >= min(xv)) & (x <= max(xv)) & (y>=min(yv)) & (y<=max(yv));

  if ~any(mask)
    in = zeros(inputSize);
    on = in;
    return
  end

  inbounds = find(mask);
  x = x(mask);
  y = y(mask);
%
% Choose block_length to keep memory usage of vec_inpolygon around 10 Megabytes.
%
  block_length = 1e5;

  M = prod(size(x));

  if M*Nv < block_length
    if nargout > 1
      [in on] = vec_inpolygon(Nv,x,y,xv,yv);
    else
      in = vec_inpolygon(Nv,x,y,xv,yv);
    end
  else
    % Process at most N elements at a time
    N = ceil(block_length/Nv);
    in = false(1,M);
    if nargout > 1
      on = false(1,M);
    end
    n1 = 0;
    n2 = 0;
    while n2 < M,
      n1 = n2+1;
      n2 = n1+N;
      if n2 > M,
        n2 = M;
      end
      if nargout > 1
        [in(n1:n2) on(n1:n2)] = vec_inpolygon(Nv,x(n1:n2),y(n1:n2),xv,yv);
      else
        in(n1:n2) = vec_inpolygon(Nv,x(n1:n2),y(n1:n2),xv,yv);
      end
    end
  end

  if nargout > 1
    onmask = mask;
    onmask(inbounds(~on)) = 0;
    on = reshape(onmask, inputSize);
  end

  mask(inbounds(~in)) = 0;
%
% Reshape output matrix.
%
  in = reshape(mask, inputSize);

  return
end
function [ in, on ] = vec_inpolygon ( Nv, x, y, xv, yv )

%*****************************************************************************80
%
%% vec_inpolygon() vectorizes the computation.
%

%
%  Translate the vertices so that the test points are at the origin.
%
  Np = length(x);
  x = x(ones(Nv,1),:);
  y = y(ones(Nv,1),:);
  xv = xv(:,ones(1,Np)) - x;
  yv = yv(:,ones(1,Np)) - y;
%
% Compute the quadrant number for the vertices relative to the test points.
%
  posX = xv > 0;
  posY = yv > 0;
  negX = ~posX;
  negY = ~posY;
  quad = (negX & posY) + 2*(negX & negY) + 3*(posX & negY);
%
% Compute the sign() of the cross product and dot product
% of adjacent vertices.
% Modified 09/17/03 to use a tolerance in the cross product sign test.
%
  m = 1:Nv-1;
  mp1 = 2:Nv;
  crossProduct = xv(m,:) .* yv(mp1,:) - xv(mp1,:) .* yv(m,:);
  tol = 10*Nv*(max(abs(xv(:)))+max(abs(yv(:))))*eps;
  crossProduct(abs(crossProduct)<tol) = 0;
  signCrossProduct = sign(crossProduct);
  dotProduct = xv(m,:) .* xv(mp1,:) + yv(m,:) .* yv(mp1,:);
%
% Compute the vertex quadrant changes for each test point.
%
  diffQuad = diff(quad);
%
% Fix up the quadrant differences.  Replace 3 by -1 and -3 by 1.
% Any quadrant difference with an absolute value of 2 should have
% the same sign as the cross product.
%
  idx = (abs(diffQuad) == 3);
  diffQuad(idx) = -diffQuad(idx)/3;
  idx = (abs(diffQuad) == 2);
  diffQuad(idx) = 2*signCrossProduct(idx);
%
% Find the inside points.
%
  in = (sum(diffQuad) ~= 0);
%
% Find the points on the polygon.  If the cross product is 0 and
% the dot product is nonpositive anywhere, then the corresponding
% point must be on the contour.
%
  on = any((signCrossProduct == 0) & (dotProduct <= 0));

  in = in | on;

  return
end

end

