function newrange=sift(range);

%	function newrange = sift(range);
%
%	This function takes a matrix series of range indeces and
%	'sifts' them to remove overlaps, and sorts them in order.
%
%	range = a matrix of range indeces to time series
%			series a  series b 
%			ra1 ra2   rb1  rb2
%			ra1 ra2   rb1  rb2
%			where each pair of columns are the starting and ending
%				indeces for marked ranges in each time series
%			range should be padded with NaN's if the number of ranges
%				for each time series is not the same
%	newrange = sifted and sorted

%	Written by Marinna Martini
%	(mmartini@usgs.gov) for the
%	U.S. Geological Survey
%	Branch of Atlantic Marine Geology, Woods Hole, MA
%	Originally written to augment markaxes.m.  
%	Depends on the m-file(s)
%		intersect.m, obtained from the Mathworks contributions ftp site
%		union.m, obtained from the Mathworks contributions ftp site

% initialize
[lr, wr] = size(range);
newrange = ones(size(range)).*NaN;	
nj=1;  ni=1;						% index into newrange

% main loop
% if this works properly, range should be all NaN's
% the loop works down the rows of each pair (set) of range indexes,
% compares each set to the ones below.  The set is either absorbed 
% into another set below it or moved to newrange if there is no
% intersection.  Finally, the last remaining set is also moved to
% newrange.
for i=1:2:wr,	% work across each time series
	if ~all(isnan(range(1:lr,i:i+1))),	% are there any ranges for this time series?
		% determine the number of range sets in this time series
		nsets = lr-sum(isnan(range(:,i)));
		for j=1:nsets-1,	% work down the set of ranges for each time series
			if ~any(any(isnan(range(j:j+1,i:i+1)))),	% are there two sets of ranges to compare?
				x = range(j,i):range(j,i+1);	% expand index into a time series
				for k=j+1:nsets, % set we will compare to 'x'
					y = range(k,i):range(k,i+1);
					xy = intersect(x,y);
					if ~isempty(xy); % there's a match, absorb this range into the intersecting one
						if length(xy)==length(y) % y is a subset of x, replace y with x
							range(k,i:i+1) = [x(1) x(length(x))]; % or range(j,i:i+1);
						elseif length(xy)==length(x), % if x is a subset of y	
							% do nothing
						else % they overlap
							range(k,i) = min([min(x) min(y) min(xy)]);
							range(k,i+1) = max([max(x) max(y) max(xy)]);
						end
						range(j,i:i+1) = [NaN NaN];	% remove the absorbed range
						break;	% stop looking for more matches
					end
				end
				if ~any(isnan(range(j,i:i+1))),	% if not absorbed, save this unique set to newrange
					newrange(nj,i:i+1) = range(j,i:i+1);
					range(j,i:i+1) = [NaN NaN];	% remove the saved range
					nj=nj+1;
				end
			else
				% this is the very last set.  Save it.
				newrange(nj,i:i+1) = range(j,i:i+1);
				range(j,i:i+1) = [NaN NaN];	% remove the saved range
				nj=1;  % reset index into newrange for new axis
			end
		end
		% this is the very last set left after absorption.  Save it.
		% this is also the only set.  Save it.
		if nsets==1, j=0; end	% for the first and only set case
		newrange(nj,i:i+1) = range(j+1,i:i+1);
		range(j+1,i:i+1) = [NaN NaN];	% remove the saved range
		nj=1;  % reset index into newrange for new axis
	end
end

% sort the ranges into ascending order
[r,c]=size(newrange);
if r>1,
	newrange=sort(newrange);
end