1D Convolution over a slice within 2D data

If I have a matrix R which tabulates some function "R(x,y)" and I wish to calculate

I(x,y) = int_{-inf}^{inf} dz P(z) R(x+z,y+z)

which essentially a convolution over a function of a different variable (X = x+y) what is the best way to go about this in Matlab? Currently I have x and y as not the same range, but they are uniformly spaced.

I have done this by making a function

FF = @(w,xq,yq) PP(w).*interpn(x,y,R,xq+w,yq+w);

and taking

norm_fct = integral(PP,wmin,wmax);
I(lp1,lp2) = integral(FF,wmin,wmax,'RelTol',1e-6)/norm_fct;

looping over lp1 and lp2 with xq=x(lp1) etc and wmin and wmax are set so I only consider the range of which I have tabulated the function. This *works* but is very slow, partly because it doesn't work vectorised and partly because I think this is quite an order N^2 way of doing an order N log(N) type problem.

Is there a function in matlab / clever method for doing this which I haven't considered?