function w=getpts(bds)
% GETPTS plots a grid and prompts the user to
%        enter by mouse click points (x,y) 
%        with bds(1) <x <bds(2), bds(3)<y<bds(4)
% nrx = number of grid subdivisions in x-direction
% nry = number of grid subdivisions in y-direction

% First plot grid
nrx=10; nry=10;
figure
hold on;
axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
axis('equal')
for i=1:(nrx+1)
  x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
  line([ x x],[ bds(3 )  bds(4)]);
end;
for j=1:(nry+1)
  y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
  line([bds(1) bds(2)],[ y y]);
end
n=0;
% Get points with mouse:
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while (but == 1)
  [xi,yi,but] = ginput(1);
  plot(xi,yi,'go')
  n = n+1;
  x(n,1) = xi;
  y(n,1) = yi;
end
% Plot the points
pts=[x y];
plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)
% Prompt for correction:
[d1,d2]=size(pts);
yn=input('Enter 1/0 to edit/not to edit points: ')
while ((yn==1)| (isempty(yn)==1))
  disp('Click with mouse the point to be edited. ')
  pointed=ginput(1);
  %find index of pointed point:
  index=1; tst=10;
  for i=1:d1
    tmp=norm(pointed-pts(i,:));
    if (tmp<tst)
      index=i; tst=tmp;
    end;
  end
  disp('Enter a new point: ');
  new=ginput(1);
  pts(index,:)=new;
% plot the grid again
  hold off;
  clf;
  axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
  axis('equal')
  for i=1:(nrx+1)
    x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
    line([ x x],[ bds(3 )  bds(4)]);
  end;
  for j=1:(nry+1)
    y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
    line([bds(1) bds(2)],[ y y]);
  end
  hold on;
% plot the corrected points
  plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
  plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)  
  yn=input('Enter 1/0 to edit/not to edit points: ');
end;
% End of correction loop
hold on;
axis([bds(1)-0.1 bds(2)+0.1 bds(3)-0.1 bds(4)+0.1]);
axis('equal')
for i=1:(nrx+1)
  x=(bds(2)- bds(1))*(i-1)/nrx+ bds(1);
  line([ x x],[ bds(3 )  bds(4)]);
end;
for j=1:(nry+1)
  y=(bds(4)- bds(3))*(j-1)/nry+ bds(3);
  line([bds(1) bds(2)],[ y y]);
end
plot(pts(:,1),pts(:,2),'ko','LineWidth',1.5)
plot(pts(:,1),pts(:,2),'b-','LineWidth',1.5)
% Save points
ok=input('Enter 1/0 to save/not save points: ');
if (ok==1)
  fid=fopen('getpts.dat','w');
  fprintf(fid, '%12.6f %12.6f \n',pts');
  fclose(fid);
  disp('Points saved in the file getpts.dat');
end;
w=pts;
% FILE getpts.m ends.
