function [tracks, tracks2f, re_unmatched] = rgblpt(setup, src_file_ary, opts)
	% tracks = rgblpt(calib, src_file_fun, opts)
	%
	% perform Lagrangian particle tracking on triple-exposed RGB particle images
	%
	% The algorithm works as follows:
	% 1)  Identify particles in red, green & blue channels of particle
	% images
	% 2)  Triangulate to find a set of candidate particles in 3D 
	% 3)  Match particle triplets 
	%
	% input:
	%	setup			a calibration structure, such as that saved by
	%					rgb_sig.m
	%	src_file_ary	filenames of images for each camera
    %                   alternatively, a cell array of images
	%	opts			a structure containing options for the tracking algorithm
	%
	% output:
	%	tracks			an array of track structures, which contain information
	%					on particle position and velocity
	%
	% mopts structure:
	%	img_threshold	brightness threshod for identifying particle images
	%					one threshold per channel (R,G,B)
	%
	%	R_search_img	radius of search for next particle in track in
	%					image space (px)
	%	R_search_obj	radius of search for tracking particles in object
	%					space (mm)
	%	K_feature		feature size of particle sets
	%	n_search		number of candidate particles to consider when matching
	%	R_tolerance_img	tolerance for matching particle projections to
	%					their observations in image space (px)
	%	R_tolerance_obj	tolerance for matching particles to their back
	%					projected rays in object space (mm)
	%	v_pred			predictor field function
	%					v_pred(x) returns velocity of predictor field
	%					where x is 3xN array of position
	%
	%  b_debug_plot		boolean flag, if true create plots for debugging purposes
	%
	% 

	%% load variables calibration
	n_cameras			= setup.n_cameras;
	camera				= setup.camera;

	% image space coordinate system
	i_axis				= camera(1).i_axis;
	j_axis				= camera(1).j_axis;
	i_lim				= sort(i_axis([1 end]));
	j_lim				= sort(j_axis([1 end]));
	[imat, jmat]		= ndgrid(i_axis, j_axis);
	n_i					= length(i_axis);
	n_j					= length(j_axis);
	bayer				= camera(1).bayer;
	delta_t				= setup.delta_t;
	t					= [-1 0 1]*delta_t;
	
	% recover triangulation and matching options structures
	mopts				= opts.matching;
	% add limits to reconstructino
	mopts.re_lim		= setup.re_bb_lim;
	% add pre-built adjacency matrix
	mopts.adjacency		= setup.adjacency;
	% add time axis for tracking
	mopts.t				= t;
	
	%% load images
	fprintf('loading images...\n');
	img_ary				= cell(1, n_cameras);
	for ic = 1 : n_cameras
        if ischar(src_file_ary{ic}) || isstring(src_file_ary{ic})
            fprintf('[C%3d] load %s\n', ic, src_file_ary{ic});
            img_ary{ic}		= cast(imread(src_file_ary{ic}), opts.precision);
        elseif all(size(src_file_ary{ic},1,2) == [n_i n_j])
            % load image from array
            img_ary{ic}     = cast(src_file_ary{ic}, opts.precision);
        else
            error('invalid src_file_ary');
        end 
    end
	
	%% search particles in image
	im_part_ary			= cell(3, n_cameras);
	im_bright_ary		= cell(3, n_cameras);
	n_2D_part			= zeros(3, n_cameras);
	for ic = 1 : n_cameras
		%% identify particle images
		[im_part, im_bright]= rgb_find_particles(camera(ic), img_ary{ic}, mopts.img_threshold);
		im_part_ary(:,ic)	= im_part;
		im_bright_ary(:,ic) = im_bright;
		n_2D_part(:,ic)		= cellfun(@(x) size(x, 2), im_part);
		fprintf('[C%3d] %d R, %d G, %d B\n', ic, n_2D_part(:,ic));
	end
	
	%% triangulate particles
	frame				= cell(1, 3);
	for rgb	= 1 : 3
		fprintf('triangulating at time step %d ...\n', rgb);
		R				= struct();
		[R.re_part, R.im_matched, R.im_res, R.im_idx, R.c_idx]	= ...
			match3d3c(camera, im_part_ary(rgb,:), mopts);
		R.n_part		= size(R.re_part,2);
		frame{rgb}		= R;
	end
	
	%% track particles
	[tracks,tracks2f,re_unmatched]	= track_frame3(camera, frame, mopts);

	%% debug plot
	if opts.b_debug_plot
		figure(100);
		clf;
		cmap	= {'r','g','b'};
		x0		= [tracks.x];
		u		= [tracks.u];
		x0_2f   = [tracks2f.x];
		u_2f	= [tracks2f.u];
		for rgb = 1 : 3
			R	= frame{rgb};
			ip	= 1:size(R.re_part,2);% match123(:, rgb);
			plot3(R.re_part(1,ip), R.re_part(2,ip), R.re_part(3,ip), '.', 'color', cmap{rgb});
			hold on;
		end
		scale	= 5*delta_t;
		% 3 frame tracks
		quiver3(x0(1,:), x0(2,:), x0(3,:), ...
				u(1,:)*scale, u(2,:)*scale, u(3,:)*scale, 0, 'k-');
		% 2 frame tracks
		quiver3(x0_2f(1,:), x0_2f(2,:), x0_2f(3,:), ...
				u_2f(1,:)*scale, u_2f(2,:)*scale, u_2f(3,:)*scale, 0, 'r-');
		axis equal tight xy;
		xlabel('$x$ [mm]', 'interpreter', 'latex', 'fontsize', 20);
		ylabel('$y$ [mm]', 'interpreter', 'latex', 'fontsize', 20);
		zlabel('$z$ [mm]', 'interpreter', 'latex', 'fontsize', 20);
	end
end
	