function [img_ary, pimg] = rgbipr_project(camera, otf, 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
	n_cameras			= length(camera);
    
    %% 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
    img_ary             = 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);
    
    for ic = 1 : n_cameras
        %% get image size
        n_i             = length(camera(ic).i_axis);
        n_j             = length(camera(ic).j_axis);
        
        %% 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);
        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(camera(ic), otf(ic), ...
                re_part_t{rgb}(:,b_visible), ...
                ipr.d_part(1, b_visible).^2, rgb);
        end
        %% accumulate
        img_ary{ic}     = sum(img_rgb, 3);
    end
    
    %% create pimg structure
    pimg        = struct('ind', pimg_ind, 'pos', pimg_pos, 'img', pimg_ary, 'map', pimg_map, 'jac', pimg_jac);
end