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
	I_part_mean		= ones(1, n_part, ipr.opts.precision);
	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(1:n_visible, 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_visible);
		
		%% solve linear least squares problem
		% this identifies the brightness we should make each particle 
		JtJ					= jac_inten.' * jac_inten;
		JtDy				= jac_inten.' * double(ipr.img_meas{ic}(:));
		% add Tikhonov regularisation
		gamma				= 1 / n_support_px;
		JtJ_regular			= JtJ + gamma * speye(n_visible);
		JtJ_bw				= bandwidth(JtJ_regular);
		Jtj_bandden			= JtJ_bw / size(JtJ_regular,1);
		fprintf('JtJ sparsity = %0.4f%% (%d elements)\n', 100*nnz(JtJ)/numel(JtJ), nnz(JtJ));
		fprintf('JtJ bandwidth= %5d (den=%0.3f%%)\n', JtJ_bw, 100*Jtj_bandden);
		I_part				= (JtJ_regular \ JtDy).';
		%I_part				= lsqlin(jac_inten, double(ipr.img_meas{ic}(:))).';
		%[I_part, stats]		= lassoglm(jac_inten, double(ipr.img_meas{ic}(:)));
		% enforce positivity
		I_part				= max(0, min(I_part, ipr.opts.I_max));
		% accumulate geometric mean
		I_part_mean(1,ind)	= I_part_mean(1,ind) .* I_part;
	end
	% geometric mean across cameras
	I_part_mean			= I_part_mean.^(1./ipr.n_views(1,:)); 
	% intensity = diameter^2
	ipr.d_part			= sqrt(I_part_mean);
	assert(isreal(ipr.d_part));
	
	%% remove dim particles
	b_accept		= ipr.d_part > sqrt(ipr.opts.I_min);
	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\n', ipr.riter, nnz(~b_accept));
	%% merge close particles
	ipr.merge(ipr.opts.matching.R_tolerance_obj);
end