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

16 次查看（过去 30 天）
Joe
Joe 2016-3-7
评论： Joe 2016-3-10
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.

请先登录，再回答此问题。

采纳的回答

Walter Roberson
Walter Roberson 2016-3-7
[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);
  1 个评论
Joe
Joe 2016-3-8
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.

请先登录，再进行评论。

更多回答（1 个）

Torsten
Torsten 2016-3-8
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.
  3 个评论
显示 1更早的评论
Torsten
Torsten 2016-3-9
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Annotations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by