How to integrate a 2D function over a grid without for loops?

I have a CCD focal plane array with 2048x2048 pixels. I need to integrate over every pixel. That is, I need to integrate over every individual pixel and not over the whole FPA. The output should be a 2048x2048 matrix and not a scalar. I'm using for loops and its taking too long, especially when I iterate other CCD parameters.

x1 = 1;
x2 = 2048;
y1 = 1;
y2 = 2048;

mux   = 1000;
muy   = 1000;
sigma = 1;

q = zeros(2048);

fun = @(x,y) (1/2*pi*sigma^2)*exp(-((x - mux).^2 +...
                                    (y - muy).^2)/(2*sigma^2));

for j = y1:y2
   for i = x1:x2
      q(j,i) = integral2(fun,i-0.5,i+0.5,j-0.5,j+0.5);
   end
end

I've seen 2 answers for 1D functions using meshgrid and arrayfun. I tried them but I must be making a mistake because I keep getting "not enough input arguments" error.

[xGrid,yGrid] = meshgrid(x1:x2,y1:y2);

q = arrayfun(@(X,Y) integral2(fun,X,Y),xGrid,yGrid);

Any help is greatly appreaciated.

Answer 1

[Walter Roberson]

[xGrid,yGrid] = meshgrid(x1:x2,y1:y2);

q = arrayfun(@(X,Y) integral2(fun, X-0.5, X+0.5, Y-0.5, Y+0.5), xGrid, yGrid);

[Joe]

You sir are a genius! But not the speed improvement I was hoping for. On large arrays, it is actually a bit slower than using for loops.

Answer 2

[Torsten]

Why not using the analytical integral ?

Solution is

q(j,i)=(pi/2*sigma^2)^2*(erf((mux-(i+0.5))/(sqrt(2)*sigma))-erf((mux-(i-0.5))/(sqrt(2)*sigma)))*(erf((muy-(j+0.5))/(sqrt(2)*sigma))-erf((muy-(j-0.5))/(sqrt(2)*sigma))))

Best wishes

Torsten.

[Joe]

Wow!! Even more genius than Walter Roberson. In all fairness, and technically, he answered my original question. But I was under the impression a vectorized approach would be faster. But integral2 is a slow function to start with. YOU actually answered my real problem of speed. The error function approach is like 50x faster than either the for loop or arrayfun approaches, without loss of precession. The only thing is that the first parameter, (pi/2*sigma^2)^2, should be just 1/4 for the numbers to match.

[Torsten]

The only thing is that the first parameter, (pi/2*sigma^2)^2, should be just 1/4 for the numbers to match.

For the function you defined above, the factor (pi/2*sigma^2)^2 is correct.

Most probably, you meant f to be

fun = @(x,y) 1/(2*pi*sigma^2)*exp(-((x - mux).^2 +...
                                    (y - muy).^2)/(2*sigma^2));

Best wishes

Torsten.

[Joe]

My bad. Typo. That's what I meant. Thanks.