function pimg = project(ipr)
    % rgbipr_project
    % create synthetic particle images
    % 
    % Input:
    %   camera 1 x n_c  camera model structure
    %   otf    1 x n_c  optical transfer function model
    %   ipr    IPR structure containing the fields
    %
    %   re_part 3 x n   particle position
    %   dx_part 3 x n   particle displacement
    %   d_part  1 x n   sqrt of particle brightness (amplitude=d^2)
    %
    % Output
    %   img_ary     1 x n_c  array of projected images
    %
    %   pimg structure (3 x n_c) containing the following fields
    %   ind         index of particle
    %   pos         image space coordinate of projected particles
    %   img         particle images
    %   map         map of particle images to image indices
    %   jac         jacobian of particle image wrt. object space coord
    %     
    %% load variables from calibration
    camera              = ipr.setup.camera;
    otf                 = ipr.setup.otf;
	n_cameras			= length(camera);
    n_part              = ipr.count();
	dtype				= ipr.opts.precision;
    
    %% evaluate displacement
    re_part_t           = cell(1,3);
    re_part_t{1}        = ipr.re_part - ipr.dx_part;
    re_part_t{2}        = ipr.re_part;
    re_part_t{3}        = ipr.re_part + ipr.dx_part;

    %% reproject and measure residual
    ipr.img_proj        = cell(1, n_cameras);
    ipr.img_res         = cell(1, n_cameras);
    pimg_pos            = cell(3, n_cameras);
    pimg_ind            = cell(3, n_cameras);
    pimg_ary            = cell(3, n_cameras);
    pimg_jac            = cell(3, n_cameras);
    pimg_map            = cell(3, n_cameras);
    pimg_res            = cell(3, n_cameras);
    pimg_ref            = cell(3, n_cameras);
    pimg_res_norm       = cell(3, n_cameras);
    
    ipr.n_views         = zeros(3, n_part, dtype);
    ipr.pimg_res        = zeros(1, n_part, dtype);
    ipr.pimg_norm       = zeros(1, n_part, dtype);

    for ic = 1 : n_cameras
        %% get image size
        n_i             = length(camera(ic).i_axis);
        n_j             = length(camera(ic).j_axis);
        i_lim			= sort(camera(ic).i_axis([1 end]));
		j_lim			= sort(camera(ic).j_axis([1 end]));
        %% reproject and test visibility
        im_part_t0      = pinhole_model(ipr.re_part, camera(ic).calib, true);
        b_visible       = im_part_t0(1,:)>=i_lim(1) & im_part_t0(1,:) <= i_lim(2) ...
                        & im_part_t0(2,:)>=j_lim(1) & im_part_t0(2,:) <= j_lim(2);
        n_visible       = nnz(b_visible);
        pimg_ind(:,ic)  = {find(b_visible)};
        fprintf('[C%d] %d particles visible (%0.3f x 10^3 part/px)\n', ic, n_visible, n_visible/(n_i*n_j)*1e3);
            
        %% project at each time step
        img_rgb         = zeros(n_i, n_j, 3, dtype);
        for rgb = 1 : 3
            %% create image
            [pimg_pos{rgb,ic}, pimg_ary{rgb,ic}, pimg_map{rgb,ic}, img_rgb(:,:,rgb), pimg_jac{rgb,ic}] = ...
                rgb_project_gpu(camera(:,ic), otf(:,ic), ...
                re_part_t{rgb}(:,b_visible), ...
                ipr.d_part(1, b_visible).^2, rgb);
        end
        %% accumulate and evaluate residual
        ipr.img_proj{ic}    = sum(img_rgb, 3);
        ipr.img_res{ic}     = ipr.img_meas{ic} - ipr.img_proj{ic};
        
        %% map residual back to individual particles
        % perform mapping for each timestep
        for rgb = 1 : 3
            % row and column coordinates of particle images
            row_map				= pimg_map{rgb,ic}(:,:,1);
            col_map				= pimg_map{rgb,ic}(:,:,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;
            % create map of indices where map is valid
            n_support_px        = size(row_map, 1);
            n_visible           = size(row_map, 2);
            ind_map				= zeros([n_support_px n_visible]);
            ind_map(map_ok)     = sub2ind([n_i n_j], row_map(map_ok), col_map(map_ok));

            % look up residual on particle images, ignoring masked regions
			%proj				= pimg_ary{rgb,ic};
			ref					= zeros([n_support_px n_visible]);
			ref(map_ok)			= ipr.img_meas{ic}(ind_map(map_ok));
			res					= zeros([n_support_px n_visible]);
			res(map_ok)			= ipr.img_res{ic}(ind_map(map_ok));
            pimg_ref{rgb,ic}    = ref; % reference to measured image
            pimg_res{rgb,ic}    = res; % residual 
            % norm of residual
            pimg_res_norm{rgb,ic}= mean(res.^2, 1); % 1 x n_visible
            % accumulate
            ind                 = pimg_ind{rgb,ic};
            ipr.n_views(rgb,ind)= ipr.n_views(rgb,ind) + 1;
            ipr.pimg_res(1,ind) = ipr.pimg_res(1,ind) + pimg_res_norm{rgb,ic};
            ipr.pimg_norm(1,ind)= ipr.pimg_norm(1,ind)+ mean(ref.^2,1);
        end
    end
    
    % normalise residual by number of views particle is visible in
    ipr.pimg_res     = ipr.pimg_res ./ sum(ipr.n_views,1);
    ipr.pimg_norm    = ipr.pimg_norm./ sum(ipr.n_views,1);
    % compute signal to noise ratio
    ipr.pimg_snr     = 10*log10(ipr.pimg_norm ./ ipr.pimg_res);
    
    %% create pimg structure
    pimg        = struct('ind', pimg_ind, ...
                         'pos', pimg_pos, ...
                         'img', pimg_ary, ...
                         'map', pimg_map, ...
                         'jac', pimg_jac, ...
                         'res', pimg_res, ...
                         'ref', pimg_ref, ...
                         'res_norm', pimg_res_norm);
end
    
    