% interactive NaN boundary marker                           Steve Mann
%
% mouse clicks... to define polygonal mask
%
% Examples: B = inan(A);      % returns image B set to NaN inside mask

function [B,pts] = interactiveinonewindow(A)

if nargout == 0
  error('inan: you must have an output argument to capture the inan image')
end%if

[M,N]=size(A);

k = 1; % point index
xs=[];
ys=[];
% ps=[];   % p=point;  ps=points
while k ~= 0  % point index
  fprintf('inan: click left for verticex %g of polygon, middle or right inside\r',k)
  [y,x,button]=ginput(1);  % Matlab x,y reversed from m,n
  fprintf('')
  if button==1
    xs=[xs;x];
    ys=[ys;y];
    % ps=[ps [cross(x,y)]]
    k = k+1;
  else%if
    k = 0;  % fall out of loop if middle or right button pressed
  end%if
end%while

K=length(xs);  % K=number of pts
% disp(sprintf('inan: %g points have been selected to define %g-gon \r',K,K))

[Y,X]=meshgrid(1:N,1:M);
polynomial_basis=0;
polynomial_basis0=0;
for k=1:K
  kplus1=k+1;
  if kplus1 > K % e.g. if kplus1 = K+1
    kplus1 = 1; % e.g. modulo K
  end%if
  thispoint = [xs(k);ys(k);1];
  nextpoint = [xs(kplus1);ys(kplus1);1];
  abc=cross(thispoint,nextpoint);
  % cross product: [x1;y1;1] X [x2,y2,1] is normal to plane containing line
  % raised to ``z'' (height) = 1 and the origin at (x,y,z) = 0
  a=abc(1);
  b=abc(2);
  c=abc(3);
  %a=ps(1,k);
  %b=ps(2,k);
  %c=ps(3,k);
  equation  = a*X+b*Y+c;  % equation should be 0 on the line
  %figure;tvs(equation);title(sprintf('equation %g',k));
  equation0 = a*x+b*y+c;  % last point selected by middle or right button
  %figure;tvs(sign(equation));title(sprintf('sign(equation %g)',k));
  polynomial_basis  = polynomial_basis  + sign(equation)*10 .^(k-1);  % matrix
  polynomial_basis0 = polynomial_basis0 + sign(equation0)*10 .^(k-1); % scalar
  %figure;tvs(polynomial_basis);title(sprintf('polynomial_basis %g',k));
end%for

absdiff = abs(polynomial_basis - polynomial_basis0);
%figure;tvs(diff);title(sprintf('diff = polynomial_basis - polynomial_basis0'));
if button==2
  mask=sign(absdiff); % 0 where inside, 1 where outside
  disp(sprintf(...
   'inan: done; region inside %g-gon set to NaN                            ',K))
end%if
if button==3
  mask=1-sign(absdiff); % 0 where inside, 1 where outside
  disp(sprintf(...
    'inan: done; only image inside %g-gon is saved; the rest=NaN           ',K))
end%if

B=A.*chval(mask,0,NaN);  % 1 means identity, 0 means change to NaN