function A = event_intersect(E1, E2)

% event_intersect - compute intersection of events
% ------------------------------------------------
%
% A = event_intersect(E1)
%
%   = event_intersect(E1, E2)
%
% Input:
% ------
%  E1, E2 - event arrays to intersect
%
% Output:
% -------
%  A - intersection area

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Revision: 539 $
% $Date: 2005-02-16 19:58:53 -0500 (Wed, 16 Feb 2005) $
%--------------------------------

% TODO: rename to 'intersect_events'

% TODO: extract component event computation from test code

%----------------------------
% TEST CODE
%----------------------------

% NOTE: there are many possible variations and uses for this kind of processing

if ((nargin < 1) || ~isstruct(E1))
	
	n = E1;
	
	for k = 1:length(n)
		test_code(n(k));
	end
	
	return;
	
end

%--
% compute depending on number of inputs
%--

switch (nargin)
	
	%--
	% self intersection
	%--
	
	case (1)
		
		%--
		% convert events to position representation
		%--
		
		pos = event_to_pos(E1);
						
		%--
		% compute intersection area of rectangles
		%--
		
		A = rectint(pos,pos);
				
		A = sparse(A);
		
	%--
	% cross intersection
	%--
	
	case (2)
		
		%--
		% convert events to position representation
		%--
		
		pos1 = event_to_pos(E1);
		
		pos2 = event_to_pos(E2);
				
		%--
		% compute intersection area of rectangles
		%--
		
		A = rectint(pos1,pos2);
				
		A = sparse(A);

end


%----------------------------------------------------------------------
% EVENT_TO_POS
%----------------------------------------------------------------------

function pos = event_to_pos(event)

% event_to_pos - get position bounds from event array
% ---------------------------------------------------
%
% pos = event_to_pos(event)
%
% Input:
% ------
%  event - event array
%
% Output:
% -------
%  pos - position rows for events

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Revision: 539 $
% $Date: 2005-02-16 19:58:53 -0500 (Wed, 16 Feb 2005) $
%--------------------------------

%--
% get time and frequency bounds from events
%--

% NOTE: we essentially convert from array of structures to structure of arrays

event = struct_field(event,'time','freq');

%--
% create position matrix
%--

% NOTE: each row is of the form [left, bottom, width, height]

pos = [event.time(:,1), event.freq(:,1), diff(event.time,1,2), diff(event.freq,1,2)];


%----------------------------------------------------------------------
% RECTINT
%----------------------------------------------------------------------

% NOTE: this is a modified copy of the MATLAB rectangle intersection

function out = rectint(A,B)

%  RECTINT Rectangle intersection area.
%
%  AREA = RECTINT(A, B) returns the area of intersection of the
%  rectangles specified by position vectors A and B.  
%
%  If A and B each specify one rectangle, the output AREA is a scalar.
%
%  A and B can also be matrices, where each row is a position vector.
%  AREA is then a matrix giving the intersection of all rectangles
%  specified by A with all the rectangles specified by B.  That is, if A
%  is M-by-4 and B is N-by-4, then AREA is an M-by-N matrix where
%  AREA(P,Q) is the intersection area of the rectangles specified by the
%  Pth row of A and the Qth row of B.
%
%  Note:  A position vector is a four-element vector [X,Y,WIDTH,HEIGHT],
%  where the point defined by X and Y specifies one corner of the
%  rectangle, and WIDTH and HEIGHT define the size in units along the x-
%  and y-axes respectively.

%  Copyright 1984-2002 The MathWorks, Inc.
%  $Revision: 539 $  $Date: 2005-02-16 19:58:53 -0500 (Wed, 16 Feb 2005) $

%-----------------------------
% HANDLE INPUT
%-----------------------------

%--
% represent A more conveniently
%--

numRectA = size(A,1);

leftA = A(:,1); bottomA = A(:,2); rightA = leftA + A(:,3); topA = bottomA + A(:,4);

%--
% represent B more conveniently
%--

numRectB = size(B,1);

% NOTE: consider transposition

leftB = B(:,1)'; bottomB = B(:,2)'; rightB = leftB + B(:,3)'; topB = bottomB + B(:,4)';

%--
% expand A representation for comparison
%--

% NOTE: check that ones multiplication is not more efficient, observed elsewhere

leftA = repmat(leftA, 1, numRectB);

bottomA = repmat(bottomA, 1, numRectB);

rightA = repmat(rightA, 1, numRectB);

topA = repmat(topA, 1, numRectB);

%--
% expand B representation for comparison
%--

% NOTE: check that ones multiplication is not more efficient, observed elsewhere

leftB = repmat(leftB, numRectA, 1);

bottomB = repmat(bottomB, numRectA, 1);

rightB = repmat(rightB, numRectA, 1);

topB = repmat(topB, numRectA, 1);

%--
% intersection computation
%--

% NOTE: separate computation to get join and intersection rectangle plus area

out = ...
	(max(0, min(rightA, rightB) - max(leftA, leftB))) .* (max(0, min(topA, topB) - max(bottomA, bottomB)))    ...
;


%----------------------------------------------------------------------
% RECTINT
%----------------------------------------------------------------------

function test_code(n)

% test_code - test event intersection code
% ----------------------------------------
%
% test_code(n)
%
% Input:
% ------
%  n - number of events to intersect (def: 20)

%--------------------------------
% Author: Harold Figueroa
%--------------------------------
% $Revision: 539 $
% $Date: 2005-02-16 19:58:53 -0500 (Wed, 16 Feb 2005) $
%--------------------------------

%--
% set number of rectangles
%--

if (~nargin || isempty(n))
	n = 20;
end

%----------------------------------------------------
% CREATE RANDOM EVENTS
%----------------------------------------------------

%--
% create random times and frequency bounds
%--

time = sort(rand(n,2),2);

time(:,2) = time(:,1) + (time(:,2) / 5);

freq = sort(rand(n,2),2);

freq(:,2) = freq(:,1) + (freq(:,2) / 5);

%--
% create events from bounds
%--

base = event_create;

for k = 1:n

	base.time = time(k,:);
	
	base.freq = freq(k,:);

	event(k) = base;

end

%----------------------------------------------------
% INTERSECT EVENTS AND COMPUTE COMPONENTS
%----------------------------------------------------

% NOTE: generalize to apply intersection recursively

A = event_intersect(event);

%--
% produce graph representation of intersection
%--

% NOTE: discard intersection weights, weights are very useful information

A = double(A ~= 0);

% NOTE: where should graph tools go in the XBAT tree?

G.E = sparse_to_edge(A);

%--
% compute connected components of intersection graph
%--

% NOTE: this code can be improved

[step,comp] = graph_dfs(G);

label = unique(comp);

%--
% compute component rectangles
%--

j = 1;

% NOTE: pad is set for display, in XBAT this is part of 'event_view'

pad = 0.005; 

for k = 1:length(label)

	%--
	% select event using component labels
	%--

	ix = find(comp == label(k));

	%--
	% disregard the trivial components, single event components
	%--

	if (length(ix) == 1)
		continue;
	end

	%--
	% set time and frequency bounds
	%--

	time = struct_field(event(ix),'time');
	
	time = [min(time(:,1)) - pad, max(time(:,2)) + pad];

	freq = struct_field(event(ix),'freq');
	
	freq = [min(freq(:,1)) - pad, max(freq(:,2)) + pad];

	%--
	% create component event
	%--

	% NOTE: here we can set the level of the events as well

	base.time = time;
	
	base.freq = freq;

	comp_event(j) = base;

	j = j + 1;

end

%----------------------------------------------------
% DISPLAY RESULTS
%----------------------------------------------------

fig;

%--
% display intitial random rectangle events
%--

for k = 1:length(event)
	
	tmp = rectangle('position',event_to_pos(event(k)));
	
	set(tmp, ...
		'edgecolor',color_to_rgb('Cyan'), ...
		'linewidth',1 ...
	);

end	

%--
% display component events
%--

for k = 1:length(comp_event)
	
	tmp = rectangle( ...
		'position',event_to_pos(comp_event(k)) ...
	);

	set(tmp, ...
		'edgecolor',[1 1 0], ... 
		'linestyle',':' ...
	);

end	

%--
% set display properties
%--

set(gca,'visible','off');

set(gcf,'color',zeros(1,3));

