function mat_conv_rs = gen_conv_mat(h,w,hu,hfull,...
    K,coeff,N,func, t_exposure,t_readout, nsmp)
       
rx = coeff(1:N);
ry = coeff(N+1:2*N);
if numel(coeff) == 2*N
    rz = zeros(1,N);
else
    rz = coeff(2*N+1:3*N);
end

ind_obs = [];
ind_bpj = [];
ind_val = [];
for ih = 1:h
    % since boundaries are padded replicately, kernels in the padded
    % region should be constant...
    ts = max(0,(hu-1+ih-1)/(hfull-1) * t_readout);
    te = min(1,ts + t_exposure);
    tsmp = linspace(ts,te,nsmp);
    wt = 1/nsmp;
    for it = 1:nsmp
        t = tsmp(it);
        r = expm([0 -func(rz,t) func(ry,t); func(rz,t) 0 -func(rx,t); -func(ry,t) func(rx,t) 0]);
        H = K*(r)/K;
        [tmp_ind_obs,tmp_ind_bpj,tmp_ind_val] = gen_warp_mat(h,w,H,ih);
        ind_obs = [ind_obs, tmp_ind_obs];
        ind_bpj = [ind_bpj, tmp_ind_bpj];
        ind_val = [ind_val, tmp_ind_val.*wt];
    end
end
mat_conv_rs = sparse(ind_obs,ind_bpj,ind_val,h*w,h*w);
    
        
        