function [im_part, f0] = rgb_find_particles(camera, img, threshold)
	% im_part, part_bright = find_particles(cmodel, img, threshold)
	%
	% find locations of particles in the image img, which is a 2D array
	% with dimensions n_i x n_j with a bayer filter
	%
	% particle locations are returned in the 2 x Np array part_ij
	% where Np is the number of particles that were found. The location is
	% determined to sub-pixel accuracy using Gaussian peak finding
	%
	% part_bright is a 1 x Np vector of the peak brightness of the particle image
	%
	% i_axis and j_axis are Ni- and Nj- vectors respectively
	% which define the coordinate system of the image
	%
	% particle images are identified by:
	%   pixel being above threshold
	%   pixel being a local maximum over its four nearest neighbours

	im_size			= size(img);
	if numel(threshold) == 1
		threshold	= [threshold threshold threshold];
	end
	
	%% get camera properties
	i_axis			= camera.i_axis;
	j_axis			= camera.j_axis;
	n_i				= length(i_axis);
	n_j				= length(j_axis);
	bayer			= camera.bayer;
	assert(size(img,1) == n_i && size(img, 2) == n_j, 'image is wrong size');
	
	%% identify particles as local maxima above a threshold
	%img				= single(img);
	[i_part,j_part,f0,b_color]= rgb_find_maxima(img, threshold);
	
	%% sub-pixel peak location
	f0				= f0(:).';
	fi				= [img(sub2ind(im_size, i_part-2, j_part  ))';
					   img(sub2ind(im_size, i_part+2, j_part  ))'];
	fj				= [img(sub2ind(im_size, i_part  , j_part-2))';
					   img(sub2ind(im_size, i_part  , j_part+2))'];
	% logarithm of brightness
	ln_fi			= log(max(fi, eps));
	ln_fj			= log(max(fj, eps));
	ln_f0			= log(f0);
	
	% subpixel shift
	% note factor of two is missing because stride is 2 rather than 1
	part_di			= (ln_fi(1,:) - ln_fi(2,:)) ./ (ln_fi(1,:) - 2*ln_f0 + ln_fi(2,:));
	part_dj			= (ln_fj(1,:) - ln_fj(2,:)) ./ (ln_fj(1,:) - 2*ln_f0 + ln_fj(2,:));
	
	%% translate to image coords
	du				= (i_axis(end) - i_axis(1)) / (im_size(1) - 1);
	dv				= (j_axis(end) - j_axis(1)) / (im_size(2) - 1);
	im_part			= [(i_part' + part_di - 1) * du + i_axis(1); 
					   (j_part' + part_dj - 1) * dv + j_axis(1)];
					
	%% decompose result in cell array
	f0				= { f0(b_color(:,1)), ...
						f0(b_color(:,2)), ...
						f0(b_color(:,3))};
	im_part			= { im_part(:,b_color(:,1)), ...
						im_part(:,b_color(:,2)), ...
						im_part(:,b_color(:,3))};
end

function [i,j,val,b_color] = rgb_find_maxima(img, threshold)
	% [i,j] = find_maxima(im, threshold)
	%
	% find indices i,j of entries in im which are local maxima
	% and above a given threshold
	
	%% find indices above lowest threshold, sufficiently far from edges
	im_size			= size(img);
	min_thresh		= min(threshold);
	[i,j]			= find(img >= min_thresh);
	ok				= i > 2 & i < im_size(1)-1 ...
					& j > 2 & j < im_size(2)-1;
	i				= i(ok);
	j				= j(ok);
	
	%% categories as red, blue or green
	% if the local maximum occurs on a red pixel, it is a "red" particle
	% likewise for green and blue
	b_blue			= mod(i,2) == 1 & mod(j, 2) == 1;
	b_red			= mod(i,2) == 0 & mod(j, 2) == 0;
	b_green			= ~b_blue & ~b_red;
	val				= img(sub2ind(im_size, i  ,j  ));
	% apply red/green/blue threshold independently
	ok				= false(size(i));
	ok(b_red)		= val(b_red ) > threshold(1);
	ok(b_green)		= val(b_green)> threshold(2);
	ok(b_blue)		= val(b_blue) > threshold(3);
	i				= i(ok);
	j				= j(ok);
	
	%% of these, find indices that are local maxima
	val				= img(sub2ind(im_size, i  ,j  ));
	is_max			= img(sub2ind(im_size, i  ,j-2)) < val ...
					& img(sub2ind(im_size, i  ,j+2)) < val ...
					& img(sub2ind(im_size, i-2,j  )) < val ...
					& img(sub2ind(im_size, i+2,j  )) < val;
	val				= val(is_max);
	i				= i(is_max);
	j				= j(is_max);
	b_blue			= mod(i,2) == 1 & mod(j, 2) == 1;
	b_red			= mod(i,2) == 0 & mod(j, 2) == 0;
	b_green			= ~b_blue & ~b_red;
	b_color			= [b_red b_green b_blue];
	
	%% filter close neighbours of green pixels
	% green pixels must be additionally maxima in a local neighbourhood
	% to avoid identifying the same particle twice close together 
	i_green			= i(b_green);
	j_green			= j(b_green);
	v_green			= val(b_green);
	i_green			= i_green(:);
	j_green			= j_green(:);
	v_green			= v_green(:);
	is_green_max	= img(sub2ind(im_size, i_green+1, j_green+1)) < v_green ...
					& img(sub2ind(im_size, i_green+1, j_green-1)) < v_green ...
					& img(sub2ind(im_size, i_green-1, j_green+1)) < v_green ...
					& img(sub2ind(im_size, i_green-1, j_green-1)) < v_green;
	n_ok			= nnz(is_green_max);
	% cut and replace
	i(b_green)		= [];
	j(b_green)		= [];
	val(b_green)	= [];
	b_color(b_green,:)=[];
	i				= [i; i_green(is_green_max)];
	j				= [j; j_green(is_green_max)];
	val				= [val; v_green(is_green_max)];
	b_color			= [b_color; false(n_ok,1), true(n_ok,1), false(n_ok,1)];
end