There seems to be an error in the integral image technique. The filtering window should be
h = zeros(WS+2,WS+2,'double');
h(1,1) = 1; h(1,WS-1) = -1;
h(WS-1,1) = -1; h(WS-1,WS-1) = 1;
that is, for a 3x3 summing window, the filtering window must be 5x5 to implement the algorithm correclty (see http://en.wikipedia.org/wiki/Summed_area_table#The_algorithm and note that points A, B and C are outside the summing window)
I tested this with a simple matrix for validation. Otherwise, thanks for the useful code.
Comment only