{"id":56,"date":"2006-05-02T08:00:48","date_gmt":"2006-05-02T12:00:48","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/?p=56"},"modified":"2023-01-09T16:08:22","modified_gmt":"2023-01-09T21:08:22","slug":"fast-local-sums","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2006\/05\/02\/fast-local-sums\/","title":{"rendered":"Fast local sums"},"content":{"rendered":"
\r\n\r\n
\r\n<\/span>\r\n

Note<\/strong><\/p>\r\n

See the 09-Jan-2023 post for updated information about this topic, including a discussion of integral images and integral box filtering.<\/p>\r\n<\/div>\r\n\r\n

The Image Processing Toolbox function normxcorr2<\/tt> computes the two-dimensional normalized cross correlation. One of the denominator terms in this computation is a matrix\r\n containing local sums of the input image.\r\n <\/p>\r\n

To see what this is, let's look at an example:<\/p>

A = magic(7)<\/pre>
\r\nA =\r\n\r\n    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\r\n<\/pre>
The 3-by-3 local sum of the (2,3) element of A can be computed this way:<\/p>
submatrix = A(1:3, 2:4);\r\nlocal_sum_2_3 = sum(submatrix(:))<\/pre>
\r\nlocal_sum_2_3 =\r\n\r\n   182\r\n\r\n<\/pre>
Computing an N-by-N local sum would seem to require N^2 - 1 additions for each location in the matrix.  Function normxcorr2<\/tt>, though, contains code contributed by Eli Horn that uses two tricks to speed up the computation.\r\n   <\/p>\r\n   
The first trick is that an N-by-N local sum can be computed by first summing along one dimension, giving a vector, and then\r\n      summing the vector.  Here's what the 3-by-3 local sum example from above looks like using this method:\r\n   <\/p>
column_sums = sum(submatrix, 1)<\/pre>
\r\ncolumn_sums =\r\n\r\n    92    63    27\r\n\r\n<\/pre>
local_sum_2_3 = sum(column_sums)<\/pre>
\r\nlocal_sum_2_3 =\r\n\r\n   182\r\n\r\n<\/pre>
To compute N-by-N local sums for the entire matrix, you could first compute N-point column sums for each matrix location,\r\n      and then you could compute N-point row sums of the result.  That gives you an average of 2N additions per matrix element,\r\n      instead of N^2-1.\r\n   <\/p>\r\n   
The second trick uses cumsum<\/tt><\/a> to reduce the computational cost of computing the column and row sums.  Let's explore this trick by looking at a single row\r\n      of A, zero-padded on the left and the right:\r\n   <\/p>
row = [0, 0, 0, A(1,:), 0, 0]<\/pre>
\r\nrow =\r\n\r\n     0     0     0    30    39    48     1    10    19    28     0     0\r\n\r\n<\/pre>
Now use cumsum<\/tt> to compute its cumulative sum:\r\n   <\/p>
c = cumsum(row)<\/pre>
\r\nc =\r\n\r\n     0     0     0    30    69   117   118   128   147   175   175   175\r\n\r\n<\/pre>
Next, form a new vector that is [c(4)-c(1), c(5)-c(2), ..., c(end)-c(end-3)].<\/p>
c_local_sums = c(4:end) - c(1:end-3)<\/pre>
\r\nc_local_sums =\r\n\r\n    30    69   117    88    59    30    57    47    28\r\n\r\n<\/pre>
Notice that c_local_sum(3) is the sum of the first three elements of A(1,:). Similarly, c_local_sum(4) is the sum of second,\r\n      third, and fourth elements of A(1,:). So this procedure gives us the local sums of a vector.\r\n   <\/p>\r\n   
The significance of this result is that it required on average only a single addition (the call to cumsum<\/tt>) and a single subtraction per vector element.  Most importantly, the number of floating-point operations per element in this\r\n      procedure is independent of N, the size of the local sum window<\/b>.\r\n   <\/p>\r\n   
Putting these two tricks together, we can compute two-dimensional local sums by first computing column sums and then computing\r\n      row sums.  The column and row sums can be computed using the cumsum<\/tt> trick.  The overall average computational cost per matrix element is thereby reduced to 2 additions and 2 subtractions. \r\n      Remarkably, this cost is independent of N.\r\n   <\/p>\r\n   
If you want to see the particular implementation of this idea in normxcorr2<\/tt>, take a look at the subfunction local_sum<\/tt>.\r\n   <\/p>\r\n