function [sel, range, i] = see(x, y, presel, iprev)

% SEE.M Browse and mark data graphically.
%  [SEL,RANGE,I]=SEE(X,Y,PRESEL,IPREV) invites
%   point selection and zooming on a plot of
%   Y(X).  The logical vector PRESEL marks
%   pre-selected points.  The result is the
%   logical vector SEL showing which points
%   remain selected when the routine exits.
%   PRESEL defaults to no preselections.
%   The vector IPREV is a list of indices
%   previously output as I by SEE, that
%   are used to restore the previous state
%   of zoom.  IPREV defaults to no previous
%   zooming.  RANGE is a list of ranges in X.
%  SEE(Z,PRESEL) derives X and Y from the
%   first two columns of array Z.  PRESEL
%   behaves as above.
%  SEE(no arguments) performs a demonstration
%   of the routine with random data.

%	Originally as see3.m by Charles R. Denham
%	Revised by Marinna Martini 
%	(mmartini@usgs.gov) for the
%	U.S. Geological Survey, 
%	Atlantic Marine Geology, Woods Hole, MA
%	9/18/93		to work with MATLAB version 4.0
%	10/14/94	to improve screen update time
%				to add range marking capability
%	Note: i is the point stack
%	This m-file file does not depend on any other m-files.

invisible = 'w';%'c0';
black = 'k';%'c1';
red = 'r';%'c2';
green = 'g';%'c3';
blue = 'b';%'c4';
cyan = 'c';%'c5';
magenta = 'm';%'c6';
yellow = 'y'; %'c7';

if nargin == 0   % Demonstration.
   help see
   m = 10;
   x = (1:m).';
   y = rand(m, 1);
   seltemp = see(x, y);
   disp(' ');
   disp('   Index Selected');
   disp([x seltemp]);
   return
end

if nargin == 1
   y = x(:, 2); x = x(:, 1);
   presel = 0 .* x;
end

if nargin == 2
   [m, n] = size(x);
   if m > 1 & n > 1
      presel = y;
      y = x(:, 2); x = x(:, 1);
     else
      x = x(:); y = y(:);
      presel = 0 .* x;
   end
   iprev = [];
end

% Setup.

x = x(:); y = y(:); sel = presel(:);
[m, n] = size(x);
if nargin == 4
   i = iprev;
  else
   i = [1 m];
end

% Monotonic behaviour in the x-direction is
%  controlled by the latter of the next two
%  lines.  When "ismonotonic" is true (non-zero),
%  points are selected by their x-position only.
%  Otherwise, selection is based on the nearest
%  point, assuming the plot is square visually.

ismonotonic = 0;                   % Disabled.
ismonotonic = all(diff(x) >= 0);   % Enabled.

xsave = []; ysave = [];
newplot = 1;
inited = 0;

preselected = sum(presel ~= 0);
disp([' Number preselected: ' int2str(preselected)]);

while 1

% Refresh the plot.
if newplot
   hold off
   cla
   i(1:2) = sort(i(1:2));
   imin = i(1); imax = i(2);
   ind = imin:imax;
   xi = x(ind); yi = y(ind);
   if i(1) ~= i(2)
      plot(xi, yi, [black '-'])   % Many points.
	  %hold on
	  % force axes to the selected zoom range
	  set(gca,'xlim',[min(xi) max(xi)]);
      if length(ind) < 50
         %plot(xi, yi,[black '+'])
	     line('xdata',xi,'ydata',yi,'linestyle','+','color','white',...
			'erasemode','none');
      end
      ii = (sel(ind) == 1);
      if any(ii), 
		 %plot(xi(ii), yi(ii), [red 'o']), 
    	  line('xdata',xi(ii),'ydata',yi(ii),'linestyle','o',...
			'color','red','erasemode','none');
	  end
	  if exist('range') == 1, 
		% plot range markers
		%plot(x(range), y(range), [green '*']), 
        %line('xdata',x(range),'ydata',y(range),'linestyle','*',...
		%	'color','green','erasemode','none');
		% change color of marked range
		for j=1:length(range(:,1)),
			line('xdata',x(range(j,1):range(j,2)),'ydata',y(range(j,1):range(j,2)),...
				'linestyle','-','color','green','erasemode','none');
		end
	  end
	  %hold off
     else
      i(1:2) = [];   % Never plot just one point.
   end
   title('[> Zoom <]            [> Unzoom <]')
   xlabel('[> Select Range <]                  [> Select Point <]')
   ylabel('[> Return <]')
   newplot = 0;
   disp([' Range: ' int2str(i(1)) ' ' int2str(i(2))])
end

% Current axes.
% This section changed by M. Martini
%[status, ax] = inquire('axis');
%xmin = ax(1); xmax = ax(2);
ax = get(gca, 'Xlim');
xmin=ax(1); xmax=ax(2);
%ymin = ax(3); ymax = ax(4);
ax = get(gca, 'Ylim');
ymin=ax(1); ymax=ax(2);
xrange = xmax - xmin;
yrange = ymax - ymin;
xmean = mean([xmin xmax]);
ymean = mean([ymin ymax]);

% Invitation.

if ~inited
   disp(' Ready to click!'), inited = 1;
end

% Get the click.
%gxprev = gx;
%gyprev = gy;
if exist('gx') ~=1, 
	gxprev=0;
else,
	gxprev = gx;
end
if exist('gy') ~=1, 
	gyprev=0;
else,
	gyprev = gy;
end
[gx, gy] = ginput(1);
gxy = gx + sqrt(-1) .* gy;


% Locate the click.
% add padding for 0 axis
pad = (xmax-xmin)*0.05;
isreturn = (gx < xmin-pad & gy > ymin & gy < ymax) | (gx > xmax & gy < ymin);
iszoom = gy > ymax & gx < xmean;
isunzoom = gy > ymax & gx > xmean;
isselect = gy < ymin & gx > xmean & gx < xmax;
isrange = gy < ymin & gx < xmean;
isgraph = ~any([isreturn isselect iszoom isunzoom isrange]);

% Return.

if isreturn
   selected = sum(sel ~= 0);
   changed = sum((presel ~= 0) ~= (sel ~= 0));
   if exist('range') == 1, 
		bracketed = length(range(:,1));
   else, bracketed=0;
   end
   disp([' Number preselected:       ' int2str(preselected)]);
   disp([' Number selected:          ' int2str(selected)]);
   disp([' Number changed:           ' int2str(changed)]);
   disp([' Number of marked ranges:  ' int2str(bracketed)]);
   return
end

% Zoom.

if iszoom & length(i) > 2
   if ~rem(length(i), 2)
      newplot = 1;
     else
      disp(' Choose another point.')
   end
end

% Unzoom.

if isunzoom
   if rem(length(i), 2)
      i(1) = [];
   end
   if length(i) > 2
      i(1:2) = [];
   end
   newplot = 1;
end

% Select Range
if isrange & length(i) > 2
	if ~rem(length(i),2)
		newplot=0;		% don't replot for now
		% push latest two points onto the range stack
		if exist('range') ~= 1,
			range = [min(i(1:2)) max(i(1:2))];
		else
			range = [min(i(1:2)) max(i(1:2)); range];
		end
		% mark the range on the plot
        line('xdata',x(range(1,1):range(1,2)),'ydata',y(range(1,1):range(1,2)),...
			'linestyle','-','color','green','erasemode','none');
		disp(['Marked range ' int2str(range(1,1)) ' to ' int2str(range(1,2))])
		% remove the range from the point stack
		i(1:2) = [];
	else
		disp('Choose another point.')
	end
end

% Mark or unmark the selected point.
% the following three lines added by M. Martini
if exist('f') ~= 1,
	f=0;
end
if isselect & any(f)
   if i(1) == f, i(1) = []; end
   sel(f) = ~sel(f);
   %hold on
   if sel(f)
      %plot(x(f), y(f), [red 'o']);
	  line('xdata',x(f),'ydata',y(f),'linestyle','o','color','red',...
			'erasemode','none');
      disp('        Selected')
     else
      %plot(x(f), y(f), [cyan 'o']);
	  line('xdata',x(f),'ydata',y(f),'linestyle','o','color','cyan',...
			'erasemode','none');
      disp('        Deselected')
   end
   %hold off
end

% Find the nearest data point.

if isgraph
   if ismonotonic   % Use x-distance only.
      temp = abs(xi - gx);
     else           % Use Euclidean distance.
      xxi = xi .* yrange; yyi = yi .* xrange;
      xytemp = xxi + sqrt(-1) .* yyi;
      gxtemp = gx .* yrange; gytemp = gy .* xrange;
      gxytemp = gxtemp + sqrt(-1) .* gytemp;
      temp = abs(xytemp - gxytemp);
   end
   f = find(temp == min(temp));
   if any(f), f = f(1) + imin - 1; end
   i = [f i];
   %hold on
   %plot(x(f), y(f), [black 'o']);
   line('xdata',x(f),'ydata',y(f),'marker','o','color','white',...
		'erasemode','none');
   %hold off
   disp([' Point: ' int2str(f)])
end

% diagnostics, display i
%disp(i)

end