function rods = findrods(u_axis, v_axis, im_proc, opts)
	% find rod-shaped objects in image with ridge detection using steerable filter
	%
	% Input:
	%	u_axis		1 x M vector image space u coordinte
	%	v_axis		1 x N vector image space v coordinate
	%	im_proc		M x N image array to search for rod shaped particles
	%	opts		optional parameters to the detection algorithm,
	%				with the following fields
	%
	%		kernel			K x K separable kernel for blob detection
	%		im_threshold	binarisation threshold to detect particle shadow
	%		b_multiple		allow multiple particles per shadow region?
	%		l_max_px		maximum length of rod, in pixels
	%
	% following
	%	Evaluation of steerable filter for detection of fibers in flowing
	%	suspensions
	%	Allan Carlsson • Karl Ha°kansson • Mathias Kvick • Fredrik Lundell • L. Daniel So¨derberg
	%	Exp Fluids (2011) 51:987–996 DOI 10.1007/s00348-011-1115-1

	%% extract important parameters
	im_size				= size(im_proc);
	% exclude edges from computation
	bw_mask				= padarray(true(im_size-2),[1 1], false, 'both');	
	% pixel pitch
	du					= mean(diff(u_axis));
	dv					= mean(diff(v_axis));
	% grid for line search
	s_line				= 0 : +opts.l_max_px;
	theta_tol			= pi / opts.n_theta;

	%% ridge detection using steerable filter
	% following
	%	Evaluation of steerable filter for detection of fibers in flowing
	%	suspensions
	%	Allan Carlsson • Karl Ha°kansson • Mathias Kvick • Fredrik Lundell • L. Daniel So¨derberg
	%	Exp Fluids (2011) 51:987–996 DOI 10.1007/s00348-011-1115-1

	% convolve with separable kernel
    kernel              = gpuArray(opts.kernel(:)); 
    % convolve in u direction
	im_filt				= convn(gpuArray(im_proc), kernel, 'same');
    % convolve in v direction
    im_filt             = convn(im_filt, kernel.', 'same');
	im_filt				= gather(double(im_filt));

	% evaluate gradient on non-shadow regions
	bw_shad				= im_proc < opts.im_threshold & bw_mask;
	[i_shad,j_shad]		= find(bw_shad);
	k_0p				= sub2ind(im_size, i_shad  , j_shad+1);
	k_p0				= sub2ind(im_size, i_shad+1, j_shad  );
	k_m0				= sub2ind(im_size, i_shad-1, j_shad  );
	k_0m				= sub2ind(im_size, i_shad  , j_shad-1);
	k_00				= sub2ind(im_size, i_shad  , j_shad  );
	k_pp				= sub2ind(im_size, i_shad+1, j_shad+1);
	k_pm				= sub2ind(im_size, i_shad+1, j_shad-1);
	k_mp				= sub2ind(im_size, i_shad-1, j_shad+1);
	k_mm				= sub2ind(im_size, i_shad-1, j_shad-1);
	% second order finite differences, see:
	%  https://en.wikipedia.org/wiki/Finite_difference
	L_x_shad			= 1/2 * (im_filt(k_p0) - im_filt(k_m0));
	L_y_shad			= 1/2 * (im_filt(k_0p) - im_filt(k_0m));
	L_xx_shad			= im_filt(k_p0) - 2*im_filt(k_00) + im_filt(k_m0);
	L_yy_shad			= im_filt(k_0p) - 2*im_filt(k_00) + im_filt(k_0m);
	L_xy_shad			= 1/4 * (im_filt(k_pp) - im_filt(k_pm) - im_filt(k_mp) + im_filt(k_mm));
	% convert to sparse format
	L_x					= sparse(i_shad, j_shad, L_x_shad, im_size(1), im_size(2));
	L_y					= sparse(i_shad, j_shad, L_y_shad, im_size(1), im_size(2));
	L_xx				= sparse(i_shad, j_shad, L_xx_shad, im_size(1), im_size(2));
	L_yy				= sparse(i_shad, j_shad, L_yy_shad, im_size(1), im_size(2));
	L_xy				= sparse(i_shad, j_shad, L_xy_shad, im_size(1), im_size(2));
	
	% optimise the objective function
	% I(theta) =  (L_xx - L_yy/3)*cos(theta)^2 
	%			+  8/3*Lxy*cos(theta)*sin(theta) 
	%			+ (L_yy - L_xx/3)*sin(theta)^2
	% this is solved as an eigenvalue problem	
	% only identify extrema in shadowed regions
	hessian				= cat(3, cat(2, L_xx_shad - L_yy_shad/3, 4/3*L_xy_shad), ...
								 cat(2, 4/3*L_xy_shad, L_yy_shad - L_xx_shad/3));
	hessian				= permute(hessian, [2 3 1]);
	% eigenvalues of hessian matrix 
	[L_pq, e_pq]		= ndsyev(single(hessian), 's+'); 
	L_pq				= permute(L_pq, [2 3 1]);
	% orientations of the maximum/minimum 
	e_p					= e_pq(:,1,:);
	e_q					= e_pq(:,2,:);
	% directional derivatives
	L_q_val				= squeeze(e_pq(2,1,:)).*L_x_shad + squeeze(e_pq(2,2,:)).*L_y_shad;
	L_p_val				= squeeze(e_pq(1,1,:)).*L_x_shad + squeeze(e_pq(1,2,:)).*L_y_shad;
	% convert to sparse
	theta_p				= atan2(e_pq(2,1,:), e_pq(1,1,:));
	theta_q				= atan2(e_pq(2,2,:), e_pq(1,2,:)); %should differ by pi/2
	L_p					= sparse(i_shad, j_shad, double(L_p_val), im_size(1), im_size(2));
	L_q					= sparse(i_shad, j_shad, double(L_q_val), im_size(1), im_size(2));
	L_pp				= sparse(i_shad, j_shad, double(L_pq(:,:,1)), im_size(1), im_size(2));
	L_qq				= sparse(i_shad, j_shad, double(L_pq(:,:,2)), im_size(1), im_size(2));
	theta_p				= sparse(i_shad, j_shad, double(theta_p(:)), im_size(1), im_size(2));
	theta_q				= sparse(i_shad, j_shad, double(theta_q(:)), im_size(1), im_size(2));
	% identify ridge pixels, which correspond to 
	% zero crossings in L_q where L_pp < 0
	bw_crossing			= padarray(L_q(2:end,:) .* L_q(1:end-1,:) < 0, [1 0], false, 'post') ...
						| padarray(L_q(:,2:end) .* L_q(:,1:end-1) < 0, [0 1], false, 'post');
	bw_ridge			= full(L_pp < 0 & L_qq > 0 & bw_crossing);
	bw_ridge			= bw_ridge & bw_shad;
	% dilate to capture near misses
	SE					= logical([0 1 0; 1 1 1; 0 1 0]);
	bw_ridge			= imdilate(bw_ridge, SE);

	%% find maxima on thresholded regions
	% get sparse representation of points on ridges
	[i_max,j_max,A_max]	= find_sparse_maxima(+L_qq, 0);
	k_max				= sub2ind(im_size, i_max, j_max);
	% restrict to extrema on ridge pixels
	b_on_ridge			= bw_ridge(k_max);
	i_pk				= i_max(b_on_ridge, :);
	j_pk				= j_max(b_on_ridge, :);
	A_pk				= A_max(:, b_on_ridge)';
	k_pk				= k_max(b_on_ridge, :);
	
	% label shadow regions
	[label,n_label]		= bwlabel(bw_shad);
	label_pk			= label(k_pk);
	
	if opts.b_multiple
		% handle multiple maxima in same shadow region
		% where there are multiple maxima, choose the largest amongst the set
		% which has similar orientation
		theta_pk			= full(theta_p(k_pk));
		% quantize into bins
		itheta_pk			= floor( mod(theta_pk,pi) / pi * opts.n_theta ) + 1; 
		% indices of adjacent bins
		itheta_pk_plus		= mod(itheta_pk  , opts.n_theta) + 1;	
		itheta_pk_minus		= mod(itheta_pk-2, opts.n_theta) + 1;
		n_bins				= [n_label opts.n_theta];
		% accumulate contributions into this bin and adjacent bins
		subs				= [label_pk(:) itheta_pk(:); label_pk(:) itheta_pk_plus(:); label_pk(:) itheta_pk_minus(:)];
		regional_maximum	= accumarray(subs, [A_pk(:); A_pk(:); A_pk(:)], n_bins, @max);
		b_is_region_max		= regional_maximum(sub2ind(n_bins, label_pk, itheta_pk)) == A_pk(:);
	else
		% identify only the largest peak in each shadow region
		regional_maximum	= accumarray(label_pk(:), A_pk(:), [n_label 1], @max);
		b_is_region_max		= regional_maximum(label_pk,1) == A_pk(:);
	end
	i_pk				= i_pk(b_is_region_max,1);
	j_pk				= j_pk(b_is_region_max,1);
	k_pk				= k_pk(b_is_region_max,1);

	%% identify ends of rods 
	n_pk				= length(i_pk);
	L_qq_full			= full(L_qq);
	theta_p_full		= full(theta_p);
	theta_pk			= theta_p_full(k_pk);
	% image coordinates of lines
	im_line_p			= cat(3, +cos(theta_pk) .* s_line * du + u_axis(1,i_pk).', ...
								 +sin(theta_pk) .* s_line * dv + v_axis(1,j_pk).');
	im_line_m			= cat(3, -cos(theta_pk) .* s_line * du + u_axis(1,i_pk).', ...
								 -sin(theta_pk) .* s_line * dv + v_axis(1,j_pk).');
	% local filter orientation sampled along line
	theta_line_p		= interpn(u_axis, v_axis, theta_p_full, im_line_p(:,:,1), im_line_p(:,:,2), 'linear');
	theta_line_m		= interpn(u_axis, v_axis, theta_p_full, im_line_m(:,:,1), im_line_m(:,:,2), 'linear');
	Lqq_line_p			= interpn(u_axis, v_axis, L_qq_full, im_line_p(:,:,1), im_line_m(:,:,2), 'linear');
	Lqq_line_m			= interpn(u_axis, v_axis, L_qq_full, im_line_m(:,:,1), im_line_m(:,:,2), 'linear');
	% find index of first pixel from center that is above threshold
	im_endpoints		= nan(n_pk, 2, 2);
	b_delete			= false(n_pk, 1);
	for n = 1 : n_pk
		% evaluate the objective function for this particle orientation
		% I(theta) =  (L_xx - L_yy/3)*cos(theta)^2 
		%			+  8/3*Lxy*cos(theta)*sin(theta)
		%			+ (L_yy - L_xx/3)*sin(theta)^2
% 		I_theta		= (L_xx - L_yy/3)*cos(theta_pk(n))^2 ...
% 					+ (L_yy - L_xx/3)*sin(theta_pk(n))^2 ...
% 					+ 8/3*L_xy*cos(theta_pk(n))*sin(theta_pk(n));
% 		I_theta			= full(-I_theta);
% 		I_theta_p		= interpn(u_axis, v_axis, I_theta, im_line_p(n,:,1), im_line_p(n,:,2), 'linear');
% 		I_theta_m		= interpn(u_axis, v_axis, I_theta, im_line_p(n,:,1), im_line_p(n,:,2), 'linear');
% 
% 		% difference in objective function relative to maximum
% 		kk_first_p		= find(I_theta_p < I_theta_p(1)/2, 1, 'first');
% 		kk_first_m		= find(I_theta_m < I_theta_m(1)/2, 1, 'first');

		% difference between angle 
		theta_diff_p	= mod(theta_line_p(n,:) - theta_pk(n) + pi/2, pi) - pi/2;
		theta_diff_m	= mod(theta_line_m(n,:) - theta_pk(n) + pi/2, pi) - pi/2;
		kk_first_p		= find(abs(theta_diff_p) > theta_tol | theta_line_p(n,:) == 0, 1, 'first');
		kk_first_m		= find(abs(theta_diff_m) > theta_tol | theta_line_m(n,:) == 0, 1, 'first');
		%Lqq_thresh		= 0.5 * Lqq_line_p(n,1);
		%kk_first_p		= find(Lqq_line_p(n,:) < Lqq_thresh, 1, 'first');
		%kk_first_m		= find(Lqq_line_m(n,:) < Lqq_thresh, 1, 'first');
		if isempty(kk_first_p) || isempty(kk_first_m)
			b_delete(n) = true;
			continue;
		end
		im_endpoints(n,1,:) = im_line_p(n,kk_first_p,:);
		im_endpoints(n,2,:) = im_line_m(n,kk_first_m,:);

		%% debug plot
		continue;
		figure(10);
		clf;
		imagesc(u_axis,v_axis,I_theta.');
		axis equal tight xy;
		hold on;
		plot(im_endpoints(n,:,1), im_endpoints(n,:,2), 'ro-');
		%plot(im_line_p(n,:,1), im_line_p(n,:,2), 'r.');
		%plot(im_line_m(n,:,1), im_line_m(n,:,2), 'g.');
		colorbar;
		xlim(im_line_p(n,1,1) + [-100 100]);
		ylim(im_line_p(n,1,2) + [-100 100]);
	end
	% apply maximum length criterion
	im_dx				= squeeze(im_endpoints(:,2,:) - im_endpoints(:,1,:));
	im_length			= sqrt(sum(im_dx.^2, 2));
	b_delete			= b_delete | im_length > opts.l_max_px;
	im_endpoints(b_delete,:,:) = [];
	n_rods				= size(im_endpoints, 1);
	im_endpoints		= permute(im_endpoints, [3 2 1]);

	% recalculate all neceessary things
	im_dx				= im_endpoints(:,2,:) - im_endpoints(:,1,:);
	im_centroid			= (im_endpoints(:,2,:) + im_endpoints(:,1,:))/2;
	im_length			= sqrt(sum(im_dx.^2, 1));

	%% create structure array
	rods				= struct('x0', [], 'x1', [], ...
							'centroid', [], ...
							'dx', [], ...
							'length', []);
	rods(1,1:n_rods)    = rods;
	
	for n = 1 : n_rods
		rods(n).x0		= im_endpoints(:,1,n);
		rods(n).x1		= im_endpoints(:,2,n);
		rods(n).centroid= im_centroid(:,1,n);
		rods(n).dx      = im_dx(:,1,n);
		rods(n).length  = im_length(n);
	end
end