Here's a related problem I encountered recently. Suppose, instead of adjusting each array value as a weighted sum of its neighbors, you want to use each array value to adjust all its neighbors. Specifically, you've got an array of mostly 1's, interspersed with the occasional 0, and you want to take each zero and surround it with 8 other 0's. One way to do this is to have a nested for-loop that checks if each element is equal to 0, and if so, sets its neighbors to 0. This is the brute force method. One annoying feature of this approach, however, is that you have to include special exceptions for any values along the edges: for example, if there's a 0 in the first row of your array, then it only has 6 neighbors to set equal to 0. Also, if the array is of any reasonably large size, nested for-loops take forever and a day to run. Anyway, here's a code snippet I wrote that avoids some of these issues:
A = ones(5); A(1,1) = 0; A(2,4) = 0; pA = padarray(A,[1 1],1); A1 = padarray(A,[2 2],1,'post'); A2 = padarray(padarray(A,[0 1],1),[2 0],1,'post'); A3 = padarray(padarray(A,[0 2],1,'pre'),[2 0],1,'post'); A4 = padarray(padarray(A,[1 0],1),[0 2],1,'post'); A5 = padarray(padarray(A,[1 0],1),[0 2],1,'pre'); A6 = padarray(padarray(A,[2 0],1,'pre'),[0 2],1,'post'); A7 = padarray(padarray(A,[0 1],1),[2 0],1,'pre'); A8 = padarray(A,[2 2],1,'pre'); zA = pA.*A1.*A2.*A3.*A4.*A5.*A6.*A7.*A8; zA(1,:) = []; zA(:,1) = []; zA(6,:) = []; zA(:,6) = [];
So, what this does is sort of like a convolution. It first 'pads' the array with 1's along all its edges, then creates 8 arrays that are shifted to every position that needs to be zeroed out. Then it multiplies all the arrays together (element-wise multiplication, of course, not matrix multiplication!). The final four lines just chop off the 1's that the array was padded with.
While this works fine, and is much faster than looping through the indices explicitly, I have a nagging question of whether this can in fact be done using a standard 2-D convolution/deconvolution operation... Hmm. In any case, in the off chance that someone will find this code useful, it's here!
UPDATE (2/20/11): As KSP points out in the comments, this is really just an interpolation. With this in mind, it is straightforward to rewrite this as a convolution:
A = ones(5); A(1,1) = 0; A(2,4) = 0; F = ones(3); zA = ~conv2(double(~A),F,'same');
(F is the convolution filter.) This does the same thing as the previous code snippet, but it runs about 3.5 times faster. Not sure how I didn't see that you could do it this way to begin with! D'oh.
I'm not sure what it is, but its darn useful. It reminds me most of what we call "moving window analysis" in geostats. Imagine that you are walking in an evenly spaced forest in which you know the location of particular nesting trees for owls. You want to make a map showing where owls might be sited. You know that owls aren't just in their one tree, and you hypothesize that the owl probably uses an equal radius around its tree for space. Moving window analysis is a form of spatial interpolation that maintains local values (that is, the tree itself is not a "kernel" that gets filtered out) while also creating a full extent map. It is similar on some scales to IDW, although IDW can be extrapolated further.
ReplyDeleteDoes that help clarify what your array does?
That makes sense. I hadn't thought of it as an interpolation scheme, but I guess it is just a really primitive form of interpolation, isn't it? That also answers the question: it is in fact a convolution. I've added an update at the end of the post with the code to do this as a convolution. Thanks for your insight :)
ReplyDelete