function [re_ipr, dx_ipr, im_ipr] = reconstruct(ipr, img_meas, re_fg, dx_fg)
	% [re_ipr, D_ipr, im_ipr] = reconstruct(ipr, img_ary, re_fg, dx_fg)
    %
	% apply RGB Iterative Particle Reconstruction to reconstruct
	% a particle field, given a first guess of the particle field
	%
	% Input:
	%	img_ary		set of images to perform reconstruction of
    %               1 x n_c cell array of images
    %
    %   re_part 3 x n_p     positions of particles at t-0
    %   dx_part 3 x n_p     displacements of particles at t=0
	%
	% Output:
    %   ipr         iterative particle reconstruction structure containing
    %      
    %   re_part     3 x n_p     positions of particles at t=0
    %   v_part      3 x n_p     velocities of particles at t=0
    %   A_part      1 x n_p     amplitudes (brightness) of particles
	%   img_recon   1 x n_c     cell array of reconstructed images
    
    %% load variables from calibration
	setup				= ipr.setup;
	n_cameras			= setup.n_cameras;
	camera				= setup.camera;
    delta_t             = setup.delta_t;
    opts                = ipr.opts;

	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;
    mopts.indent		= '        ';
	mopts.opts.b_verbose= ipr.opts.b_verbose;
    
	%% first guess
	ipr.n_part_start    = zeros(1, ipr.opts.n_ipr);
	ipr.n_part_found    = zeros(1, ipr.opts.n_ipr);
	ipr.n_part_remove   = zeros(5, ipr.opts.n_ipr);
	ipr.n_part_merge    = zeros(1, ipr.opts.n_ipr);
    ipr.initialise(img_meas, re_fg, dx_fg);
    
	%% iterate
	for riter = 1 : ipr.opts.n_ipr
		%% iteration parameters
		ipr.riter		= riter; % iteration counter
		ipr.n_part_start(riter) = ipr.count();

		%% override brightness
		%warning('override brightness');
		%ipr.d_part		= ones(1, ipr.count(), ipr.opts.precision) * sqrt(500);

		%% shake particles
		if ipr.opts.b_verbose
			fprintf('        shaking...\n');
		end
		t_shake=tic;
        ipr.shake();
		t_shake=toc(t_shake);
		fprintf('       took %0.1f seconds\n', t_shake);
        
		%% eliminate bad matches
		% retain only those particles where the signal to noise ratio is
		% better than all the alternatives, and if it is above a minimum
		% threshold, and is consistent with its neighbours
		b_reason			= ipr.validate();
		b_reject			= any(b_reason, 1);
		b_accept			= ~b_reject;
		ipr.n_part_remove(1:end-1,riter) = sum(b_reason,2);
		ipr.n_part_remove(    end,riter) = sum(b_reject,2);
		ipr.select(b_accept);
		
        if ipr.opts.b_verbose
			n_bad           = nnz(~b_accept);
			fprintf('I%02d:  %d particles with low SNR\n', ipr.riter, nnz(b_reason(1,:)));
			fprintf('I%02d:  %d particles too dim\n', ipr.riter, nnz(b_reason(2,:)));
			fprintf('I%02d:  %d particles outside velocity limits\n', ipr.riter, nnz(b_reason(3,:)));
			fprintf('I%02d:  %d particles inconsistent with neighbours\n', ipr.riter, nnz(b_reason(4,:)));
			fprintf('I%02d:  eliminated %5d bad fit / dim / too fast particles\n', ipr.riter, n_bad);
		end
		
		%% merge particles
		% merge particles which are too close together to discriminate
		n_before_merge		= ipr.count();
		ipr.merge(opts.matching.R_tolerance_obj);
		n_after_merge		= ipr.count();
		n_merged			= n_before_merge-n_after_merge;
		if ipr.opts.b_verbose
			fprintf('I%02d:  merged %5d -> %5d particles (%5d)\n', ipr.riter, n_before_merge, n_after_merge, n_merged);
		end
		
		if ipr.riter == ipr.opts.n_ipr
            break;
        end
		
        %% reproject particles and calculate residual image
        if ipr.opts.b_verbose
			fprintf('I%02d:  reprojecting %d particles ...\n', ipr.riter, ipr.count());
        end
        % reproject and measure residual
        ipr.project();
        
        %% retriangulate
        % identify particles from residual images in frame 2
        % then estimate a particle velocity from predictor field

		% build predictor field
		interpolant		= buildScatteredInterpolant3(ipr.re_part, ipr.dx_part / ipr.setup.delta_t);
		v_pred			= @(x) scatteredInterpolant3(interpolant, x);
		ipr.opts.matching.v_pred = v_pred;
		
        % perform particle tracking with predictor field
        [tracks,tracks2f,re_unmatched]= rgblpt(ipr.setup, ipr.img_res, ipr.opts);
        re_part_new     = cat(2, [tracks.x], [tracks2f.x]);
        dx_part_new     = cat(2, [tracks.u], [tracks2f.u])*delta_t;
		% add also unmatched particles predicted to time t=0
		re_pred_t0		= cell(1, 3);
		dx_pred_t0		= cell(1, 3);
		for rgb = 1 : 3
			dx_pred_t0{rgb}= t(rgb)*v_pred(re_unmatched{rgb});
			re_pred_t0{rgb}= re_unmatched{rgb} - dx_pred_t0{rgb};
		end
		re_part_new		= cat(2, re_part_new, cell2mat(re_pred_t0));
		dx_part_new		= cat(2, dx_part_new, cell2mat(dx_pred_t0));
		n_new			= size(re_part_new,2);
        
		% show residual at this step
		if ipr.opts.b_debug_plot
			ipr.plotresiduals(tracks);
		end
        
        % build new reconstruction and merge
        ipr1            = copy(ipr);
		% initialise new particles based on original images
        ipr1.initialise(ipr.img_meas, re_part_new, dx_part_new);
		% add to current model
		ipr.n_part_found(riter) = n_new;
		ipr.n_part_merge(riter) = ipr1.count();
        ipr.append(ipr1);
        ipr.merge(opts.matching.R_tolerance_obj);

		% report on this step
		if ipr.opts.b_verbose
			fprintf('I%02d:  found  %6d triplets\n', ipr.riter, length(tracks));
			fprintf('I%02d:  found  %6d pairs\n', ipr.riter, length(tracks2f));
			fprintf('I%02d:  found  %6d particles\n', ipr.riter, size(re_unmatched{2},2));
			fprintf('I%02d:  found  %6d total\n', ipr.riter, n_new);
			fprintf('I%02d:  accept %6d\n', ipr.riter, ipr1.count());
			fprintf('I%02d:  merged to %d particles ...\n', ipr.riter, ipr.count());
		end
    end
end