A = <\/span>7\u00d77<\/span><\/div><\/div>
30 39 48 1 10 19 28\r\n 38 47 7 9 18 27 29\r\n 46 6 8 17 26 35 37\r\n 5 14 16 25 34 36 45\r\n 13 15 24 33 42 44 4\r\n 21 23 32 41 43 3 12\r\n 22 31 40 49 2 11 20\r\n<\/div>
<\/div>
<\/div><\/div><\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div>
The local sum of the <\/span>$3 \\times 3$<\/span> neighborhood around (2,3) element of <\/span>A<\/span> can be computed this way:<\/span><\/div>submatrix = A(1:3, 2:4)<\/span><\/span><\/div>submatrix = <\/span>3\u00d73<\/span><\/div><\/div>    39    48     1\r\n    47     7     9\r\n     6     8    17\r\n<\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>local_sum_2_3 = sum(submatrix,<\/span>\"all\"<\/span>)<\/span><\/span><\/div>local_sum_2_3 = 182<\/div><\/div><\/div><\/div>In image processing applications, we are often interested in the local sum divided by the number of elements in the neighborhood, which is 9 in this case:<\/span><\/div>local_sum_2_3\/9<\/span><\/span><\/div>ans = 20.2222<\/div><\/div><\/div><\/div>Computing all the local sums for the matrix, divided by the neighborhood size, is the same as performing 2-D filtering with a local mean filter, like this:<\/span><\/div>B = imfilter(A,ones(3,3)\/9);<\/span><\/span><\/div><\/div>B(2,3)<\/span><\/span><\/div>ans = 20.2222<\/div><\/div><\/div><\/div>This is often called a <\/span>box filter<\/span>. Computing an <\/span>$N \\times N$<\/span> box filter for an entire <\/span>$P \\times Q$<\/span> image would seem to require approximately <\/span>$PQN^2$<\/span> operations. In other words, for a fixed image size, we might expect the box filter computational time to increase in proportion to <\/span>$N^2$<\/span>, where <\/span>N<\/span> is the local sum window size. <\/span><\/div>Let's see how long <\/span>imfilter<\/span> takes with respect to a changing local sum neighborhood size.<\/span><\/div>I = rand(3000,4000);<\/span><\/span><\/div><\/div><\/div>nn = 15:10:205;<\/span><\/span><\/div><\/div>imfilter_times = zeros(size(nn));<\/span><\/span><\/div><\/div>for <\/span>k = 1:length(nn)<\/span><\/span><\/div><\/div>    n = nn(k);<\/span><\/span><\/div><\/div>    f = @() imfilter(I,ones(n,n)\/n^2);<\/span><\/span><\/div><\/div>    imfilter_times(k) = timeit(f);<\/span><\/span><\/div><\/div>end<\/span><\/span><\/div><\/div><\/div>plot(nn,imfilter_times)<\/span><\/span><\/div><\/div>title(<\/span>\"Box filter execution time using imfilter\"<\/span>)<\/span><\/span><\/div><\/div>xlabel(<\/span>\"window size\"<\/span>)<\/span><\/span><\/div><\/div>ylabel(<\/span>\"time (s)\"<\/span>)<\/span><\/span><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>The function <\/span>imfilter<\/span> is doing better than expected, with an execution time that increases only linearly with the window size. (That's because <\/span>imfilter<\/span> is taking advantage of filter separability, <\/span>which I have written about previously<\/span><\/a>.)<\/span><\/div>But we can compute local sums even faster, by using something called an <\/span>integral image<\/span>, also known as a <\/span>summed-area table<\/span><\/a>. Let's go back to the magic square example and use the function <\/span>integralImage<\/span>.<\/span><\/div>Ai = integralImage(A)<\/span><\/span><\/div>Ai = <\/span>8\u00d78<\/span><\/div><\/div>           0           0           0           0           0           0           0           0\r\n           0          30          69         117         118         128         147         175\r\n           0          68         154         209         219         247         293         350\r\n           0         114         206         269         296         350         431         525\r\n           0         119         225         304         356         444         561         700\r\n           0         132         253         356         441         571         732         875\r\n           0         153         297         432         558         731         895        1050\r\n           0         175         350         525         700         875        1050        1225\r\n<\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>Each element of <\/span>Ai<\/span> is the cumulative sum of all the elements in A that are above and to the left. <\/span>Ai(4,3)<\/span>, which is 206, is the sum of all the elements of <\/span>A(1:3,1:2)<\/span>. <\/span><\/div>A(1:3,1:2)<\/span><\/span><\/div>ans = <\/span>3\u00d72<\/span><\/div><\/div>    30    39\r\n    38    47\r\n    46     6\r\n<\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>sum(ans,<\/span>\"all\"<\/span>)<\/span><\/span><\/div>ans = 206<\/div><\/div><\/div>Ai(4,3)<\/span><\/span><\/div>ans = 206<\/div><\/div><\/div><\/div>Expressed as a general formula, the integral image is <\/span>$A_i(m,n) = \\sum_{p=1}^{m-1} \\sum_{q=1}^{n-1} A(p,q)$<\/span>. <\/span><\/div>The really interesting thing about having the integral image, <\/span>Ai<\/span>, is that it lets you compute the sum of <\/span>any<\/span> rectangular submatrix of <\/span>A<\/span> by an addition and subtraction of just four numbers. Specifically, the sum of all the elements in <\/span>A(m1:m2,n1:n2)<\/span> is <\/span>Ai(m2+1,n2+1) - Ai(m1,n2+1) - Ai(m2+1,n1) + Ai(m1,n1)<\/span>. Consider computing the sum of all elements in <\/span>A(3:5,2:6)<\/span>.<\/span><\/div>m1 = 3;<\/span><\/span><\/div><\/div>m2 = 5;<\/span><\/span><\/div><\/div>n1 = 2;<\/span><\/span><\/div><\/div>n2 = 6;<\/span><\/span><\/div><\/div>A(m1:m2,n1:n2)<\/span><\/span><\/div>ans = <\/span>3\u00d75<\/span><\/div><\/div>     6     8    17    26    35\r\n    14    16    25    34    36\r\n    15    24    33    42    44\r\n<\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>sum(ans,<\/span>\"all\"<\/span>)<\/span><\/span><\/div>ans = 375<\/div><\/div><\/div>Ai(m2+1,n2+1) - Ai(m1,n2+1) - Ai(m2+1,n1) + Ai(m1,n1)<\/span><\/span><\/div>ans = 375<\/div><\/div><\/div><\/div>If we have the integral image, then, we can compute every element of the box filter output by adding just four terms together. In other words, the box filter computation can be independent of the box filter window size!<\/span><\/div>The Image Processing Toolbox function <\/span>integralBoxFilter<\/span> does this for you. You give it the integral image, and it returns the result of the box filter applied to the original image.<\/span><\/div>Let's see how long that takes as we vary the box filter window size.<\/span><\/div>integral_box_filter_times = zeros(1,length(nn));<\/span><\/span><\/div><\/div>for <\/span>k = 1:length(nn)<\/span><\/span><\/div><\/div>    n = nn(k);<\/span><\/span><\/div><\/div>    h = @() integralBoxFilter(Ai,n);<\/span><\/span><\/div><\/div>    integral_box_filter_times(k) = timeit(h);<\/span><\/span><\/div><\/div>end<\/span><\/span><\/div><\/div><\/div>plot(nn,integral_box_filter_times, <\/span>...<\/span><\/span><\/div><\/div>    nn,imfilter_times)<\/span><\/span><\/div><\/div>legend([<\/span>\"integralBoxFilter\" \"imfilter\"<\/span>], <\/span>...<\/span><\/span><\/div><\/div>    Location = <\/span>\"northwest\"<\/span>)<\/span><\/span><\/div><\/div>title(<\/span>\"Box filter execution time\"<\/span>)<\/span><\/span><\/div><\/div>xlabel(<\/span>\"window size\"<\/span>)<\/span><\/span><\/div><\/div>ylabel(<\/span>\"time (s)\"<\/span>)<\/span><\/span><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>That looks great! As expected, <\/span>integralBoxFilter<\/span> is fast, with execution time that is independent of the window size. To be completely fair, though, we should include the time required to compute the integral image.<\/span><\/div>integral_box_filter_total_times = zeros(1,length(nn));<\/span><\/span><\/div><\/div>for <\/span>k = 1:length(nn)<\/span><\/span><\/div><\/div>    n = nn(k);<\/span><\/span><\/div><\/div>    r = @() integralBoxFilter(integralImage(I),n);<\/span><\/span><\/div><\/div>    integral_box_filter_total_times(k) = timeit(r);<\/span><\/span><\/div><\/div>end<\/span><\/span><\/div><\/div><\/div>plot(nn,imfilter_times,<\/span>...<\/span><\/span><\/div><\/div>    nn,integral_box_filter_times,<\/span>...<\/span><\/span><\/div><\/div>    nn,integral_box_filter_total_times)<\/span><\/span><\/div><\/div>title(<\/span>\"Box filter execution time\"<\/span>)<\/span><\/span><\/div><\/div>xlabel(<\/span>\"window size\"<\/span>)<\/span><\/span><\/div><\/div>ylabel(<\/span>\"time (s)\"<\/span>)<\/span><\/span><\/div><\/div>legend([<\/span>\"imfilter\" \"integralBoxFilter\" \"integralImage + integralBoxFilter\"<\/span>],<\/span>...<\/span><\/span><\/div><\/div>    Location = <\/span>\"northwest\"<\/span>)<\/span><\/span><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div><\/div>From the plot, we can see that using the integral image for box filtering is the way to go, even when you take into consideration the time needed to compute the integral image.<\/span><\/div>\r\n