How to integrate a 2D function over a grid without for loops?
17 views (last 30 days)
Show older comments
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.
0 Comments
Accepted Answer
Walter Roberson
on 7 Mar 2016
[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);
More Answers (1)
Torsten
on 8 Mar 2016
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 Comments
Torsten
on 9 Mar 2016
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.
See Also
Categories
Find more on Annotations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!