function initialise(ipr, img_meas, re_fg, dx_fg)
    % initialise reconstruction with a first guess of particle positions
    % and displacements
    %
    % Input:
    %   img_ary     1 x n_cameras   particle images
    %   re_fg       3 x n_part      first guess particle position
    %   dx_fg       3 x n_part      first guess particle displacement
    
    %% sanity check input
    n_part          = size(re_fg, 2);
    assert(size(re_fg, 1) == 3);
    assert(size(re_fg, 1) == 3);
    assert(size(dx_fg, 2) == n_part);
    assert(numel(img_meas) == ipr.n_cameras);
    dtype			= ipr.opts.precision;
	n_cameras		= ipr.setup.n_cameras;

    %% assign image array
    ipr.img_meas    = img_meas(:).';
    ipr.re_part     = re_fg;
    ipr.dx_part     = dx_fg;
    ipr.d_part      = ones(1,n_part,dtype);
    ipr.lambda      = ipr.opts.lambda0 * ones(1, n_part, dtype);

	%% sort for jacobian sparsity
	re_lim          = ipr.setup.re_bb_lim;
    re_fov          = abs(re_lim(2,:) - re_lim(1,:));
	[~,i_long]      = max(re_fov);
    e_sort          = zeros(3,1);
    e_sort(i_long)  = 1;
	ipr.sort(e_sort);
    
    %% project with unit brightness
    pimg            = ipr.project();
    
    %% evaluate brightness based on middle frame
	JtJ_sum			= sparse(n_part, n_part);
	Jty_sum			= zeros(n_part, 1);
	for ic = 1 : n_cameras
		%% size of this image
		n_i					= size(ipr.img_meas{ic},1);
		n_j					= size(ipr.img_meas{ic},2);

        %% evaluate transfer function from particle brightness to image
		% jacobian is of size n_px x n_visible
		% we construct a model
		% img_meas{ic}(:) = jac_inten * I_part(:)
		% then solve in the least squares sense for the particle brightness
		
		% indices of particles 
		ind                 = pimg(2,ic).ind;
		n_visible			= size(ind, 2);
		% row and column coordinates of particle images
        row_map				= pimg(2,ic).map(:,:,1);
        col_map				= pimg(2,ic).map(:,:,2);
        % take care of out of bounds elements in map
        % dep_map is always in bounds by construction
        map_ok				= row_map >= 1 & row_map <= n_i ...
                            & col_map >= 1 & col_map <= n_j;
        % indices of pixels in image (rows of jacobian)
		n_support_px		= size(row_map,1);
        rows				= zeros([n_support_px n_visible]);
        rows(map_ok)		= sub2ind([n_i n_j], row_map(map_ok), col_map(map_ok));
		% indicies of visible particles (columns of jacobian)
		cols     			= repmat(ind, n_support_px, 1);
        % OTF values for each pixel/particle
		vals				= double(pimg(2,ic).img);
		jac_inten			= sparse(rows(map_ok), cols(map_ok), vals(map_ok), n_i*n_j, n_part);
		
		%% accumulate jacobian for linear least squares problem
		JtJ					= jac_inten.' * jac_inten;
		Jty					= jac_inten.' * double(ipr.img_meas{ic}(:));
		JtJ_sum				= JtJ_sum + JtJ;
		Jty_sum				= Jty_sum + Jty;
	end

	%% solve linear least squares problem
	% arg min_x (y - J*x)^2 + gamma*x^2
	% 
	% where y is set of observations and x is the set of brightnesses
	% we use Tikhonov regularisation with ridge parameter gamma to avoid 
	% ill-posed problems, which occur when particles are ambiguous or are
	% duplicated in the initialisation
	
	% determine the ridge parameter lambda
	% there is an optimal way of doing this:
	% https://en.wikipedia.org/wiki/Tikhonov_regularization
	% but this estimate seems to work ok...
	lambda					= 1/n_support_px;
	JtJ_regularised			= JtJ_sum + speye(n_part)*lambda;

	JtJ_bw					= bandwidth(JtJ_regularised);
	JtJ_bandden				= JtJ_bw / size(JtJ_regularised,1);
	fprintf('I%02d      JtJ sparsity = %0.4f%% (%d elements)\n', ipr.riter, 100*nnz(JtJ_regularised)/numel(JtJ_regularised), nnz(JtJ_regularised));
	fprintf('I%02d      JtJ bandwidth= %5d (den=%0.3f%%)\n', ipr.riter, JtJ_bw, 100*JtJ_bandden);
	% solve linear least squares problem
	I_part					= (JtJ_regularised \ Jty_sum).';
	% cast to working precision
	I_part                  = cast(I_part, ipr.opts.precision);
	I_part					= max(0, I_part);
	ipr.d_part				= sqrt(I_part);
	assert(isreal(ipr.d_part));
	
	%% remove dim particles
	b_accept		= ipr.d_part > sqrt(ipr.opts.I_min) ...
					& ipr.d_part < sqrt(ipr.opts.I_max);
	disp(prctile(ipr.d_part, [5 25 50 75 95]));
	ipr.select(b_accept);
	fprintf('I%02d     reject %d new particles because they are too dim/too bright\n', ipr.riter, nnz(~b_accept));
	%% merge close particles
	ipr.merge(ipr.opts.matching.R_tolerance_obj);
end