How to perform iterative 2D operation on 4D numpy array
image-processing, multidimensional-array, numpy, python, scipy
Solution
When you need to multiply element-wise, then reduce with addition, think `np.dot` or `np.einsum`:
from numpy.lib.stride_tricks import as_strided
arr = np.random.rand(256, 256)
mask = np.random.rand(3, 3)
arr_view = as_strided(arr, shape=(254, 254, 3, 3), strides=arr.strides*2)
arr[1:-1, 1:-1] = np.einsum('ijkl,kl->ij', arr_view, mask)
Problem
Let me preface this post by saying that I'm pretty new to Python and NumPy, so I'm sure I'm overlooking something simple. What I'm trying to do is image processing over a PGM (grayscale) file using a mask (a mask convolution operation); however, I don't want to do it using the SciPy all-in-one imaging processing libraries that are available—I'm trying to implement the masking and processing operations myself. What I want to do is the following: - Iterate a 3x3 sliding window over a 256x256 array - At each iteration, I want to perform an operation with a 3x3 image mask (array that consists of fractional values < 1 ) and the 3x3 window from my original array - The operation is that the image mask gets multiplied by the 3x3 window, and that the results get summed up into one number, which represents a weighted average of the original 3x3 area - This sum should get inserted back into the center of the 3x3 window, with the original surrounding values left untouched - However, the output of one of these operations shouldn't be the input of the next operation, so a new array should be created or the original 256x256 array shouldn't be updated until all operations have completed. The process is sort of like this, except I need to put the result of the convolved feature back into the center of the window it came from: (source: stanford.edu) So, in this above example, the `4` would go back into the center position of the 3x3 window it came from (after all operations had concluded), so it would look like `[[1, 1, 1], [0, 4, 1], [0, 0, 1]]` and so on for every other convolved feature obtained. A non-referential copy could also be made of the original and this new value inserted into that. So, this is what I've done so far: I have a 256x256 2D numpy array which is my source image. Using `as_strided,` I convert it into a 4D numpy array of 3x3 slices. The main problem I'm facing is that I want to execute the operation I've specified over each slice. I'm able to perform it on one slice, but in `npsum` operations I've tried, it adds up all the slices' results into one value. After this, I either want to create a new 256x256 array with the results, in the fashion that I've described, or iterate over the original, replacing the middle values of each 3x3 window as appropriate. I've tried using `ndenumerate` to change just the same value `(v, x, 1, 1)` of my 4D array each time, but since the index returned from my 4D array is of the form `(v, x, y, z),` I can't seem to figure out how to only iterate through `(v, x)` and leave the last two parts as constants that shouldn't change at all. Here's my code thus far: ``` import numpy as np from numpy.lib import stride_tricks # create 256x256 NumPy 2D array from image data and image size so we can manipulate the image data, then create a 4D array of strided windows # currently, it's only creating taking 10 slices to test with imageDataArray = np.array(parsedPGMFile.imageData, dtype=int).reshape(parsedPGMFile.numRows, parsedPGMFile.numColumns) xx = stride_tricks.as_strided(imageDataArray, shape=(1, 10, 3, 3), strides=imageDataArray.strides + imageDataArray.strides) # create the image mask to be used mask = [1,2,1,2,4,2,1,2,1] mask = np.array(mask, dtype=float).reshape(3, 3)/16 # this will execute the operation on just the first 3x3 element of xx, but need to figure out how to iterate through all elements and perform this operation individually on each element result = np.sum(mask * xx[0,0]) ``` Research from sources like http://wiki.scipy.org/Cookbook/GameOfLifeStrides, http://www.johnvinyard.com/blog/?p=268, and http://chintaksheth.wordpress.com/2013/07/31/numpy-the-tricks-of-the-trade-part-ii/ were very helpful (as well as SO), but they don't seem to address what I'm trying to do exactly (unless I'm missing something obvious). I could probably use a ton of `for` loops, but I'd rather learn how to do it using these awesome Python libraries we have. I also realize I'm combining a few questions together, but that's only because I have the sneaking suspicion that this can all be done very simply! Thanks in advance for any help!