By example:
import numpy as np
k=2
kern=np.ones(2*k+1)/(2*k+1)
arr=np.random.random((10))
out=np.convolve(arr,kern, mode='same')
By example:
import numpy as np
k=2
kern=np.ones(2*k+1)/(2*k+1)
arr=np.random.random((10))
out=np.convolve(arr,kern, mode='same')
I found a solution that is inspired from a matlab function I stumbled across while trying to implement a different edge behavior with no padding in mean filtering: https://www.mathworks.com/matlabcentral/fileexchange/23287-smooth2a . It works with matrix multiplication, no loop or convolve is involved. With 1D array, it would give something like this:
import scipy.sparse
import numpy as np
def mean_filter(arr, k):
p = len(arr)
diag_offset = np.linspace(-(k//2), k//2, k, dtype=int)
eL = scipy.sparse.diags(np.ones((k, p)), offsets=diag_offset, shape=(p, p))
nrmlize = eL @ np.ones_like(arr)
return (eL @ arr) / nrmlize
How this works is you create a matrix with as many ones diagonals as your kernel size, then do the dot product with your array. The result is an array where each element is the sum of the elements of the original array over the kernel. Basically, this is the equivalent of a convolution. Finally you normalize this sum by the kernel size.
One thing to note is that at the edges of the array where you cannot build a full kernel, as much of the kernel that fits on your array is used (this is why you cannot just normalize by a scalar afterwards).
This solution can also be extended to work on multidimensional arrays.
Videos
Code review
Peilonrays points out a mixup with the out-of-bounds testing that is valid. The statement if j ... must be within the loop for k .... One of the results is that you add a different number of elements to temp depending on which boundary you're at. But there are better ways to avoid out-of-bounds indexing, see below.
Your biggest bug, however, is that you write the result of the filter into the image you are processing. Median filtering cannot be done in-place. When you update data[i][j], you'll be reading the updated value to compute data[i][j+1]. You need to allocate a new image, and write the result there.
I would suggest not adding zeros for out-of-bounds pixels at all, because it introduces a bias to the output pixels near the boundary. The clearest example is for the pixels close to any of the corners. At the corner pixel, with a 3x3 kernel, you'll have 4 image pixels covered by the kernel. Adding 5 zeros for the out-of-bounds pixels guarantees that the output will be 0. For larger kernels this happens in more pixels of course. Instead, it is easy to simply remove the temp.append(0) statements, leading to a non-biased result. Other options are to read values from elsewhere in the image, for example mirroring the image at the boundary or extending the image by extrapolation. For median filtering this has little effect, IMO.
You set temp = [] at the very beginning of your function, then reset it when you're done using it, in preparation for the next loop. Instead, initialize it once inside the main double-loop over the pixels:
for i in range(len(data)):
for j in range(len(data[0])):
temp = []
# ...
You're looping over i and j as image indices, then over z and c or k for filter kernel indices. c and k have the same function in two different loops, I would suggest using the same variable for that. z doesn't really fit in with either c or k. I would pick two names that are related in the way that i and j are, such as m and n. The choice of variable names is always very limited if it's just one letter. Using longer names would make this code clearer: for example img_column, img_row, kernel_column, kernel_row.
Out-of-bounds checking
This concludes my comments on your code. Now I'd like to offer some alternatives for out-of-bounds checking. These tests are rather expensive when performed for every pixel -- it's a test that is done \$n k\$ times (with \$n\$ pixels in the image and \$k\$ pixels in the kernel). Maybe in Python the added cost is relatively small, it's an interpreted language after all, but for a compiled language these tests can easily amount to doubling processing time. There are 3 common alternatives that I know of. I will use border = filter_size // 2, and presume filter_size is odd. It is possible to adjust all 3 methods to even-sized filters.
Separate loops for image border pixels
The idea here is that the loop over the first and last border pixels along each dimension are handled separately from the loop over the core of the image. This avoids all tests. But it does require some code duplication (all in the name of speed!).
for i in range(border):
# here we loop over the kernel from -i to border+1
for i in range(border, len(data)-border):
# here we loop over the full kernel
for i in range(len(data)-border, len(data)):
# here we loop over the kernel from -border to len(data)-i
Of course, within each of those loops, a similar set of 3 loops is necessary to loop over j. The filter logic is thus repeated 9 times. In a compiled language, where this is the most efficient method, code duplication can be avoided with inlined functions or macros. I don't know how a Python function call compares to a bunch of tests for out-of-bounds access, so can't comment on the usefulness of this method in Python.
A separate code path for border pixels
The idea here is to do out-of-bounds checking only for those pixels that are close to the image boundary. For pixels within the border, you use a version of the filtering logic with out-of-bounds checking. For the pixels in the core of the image (which is the big majority of pixels), you use a second version of the logic without out-of-bounds checking.
for i in range(len(data)):
i_border = i < border or i >= len(data)-border
for j in range(len(data[0])):
j_border = j < border or j >= len(data)-border
if i_border or j_border:
# filtering with bounds checking
else:
# filtering without bounds checking
Padding the image
The simplest solution, and also the most flexible one, is to create a temporary image that is larger than the input image by 2*border along each dimension, and copy the input image into it. The "new" pixels can be filled with zeros (to replicate what OP intended to do), or with values taken from the input image (for example by mirroring the image at the boundary or extrapolating in some other way).
The filter now never needs to check for out-of-bounds reads. When the filtering kernel is placed over any of the input image pixels, all samples fall within the padded image.
Since for this type of filtering it is necessary to create a new output image anyway (it is not possible to compute it in-place, as I mentioned before), this is not a huge cost: the original input image can now be re-used as output image.
This solution leads to the simplest code, allows for all sorts of boundary extension methods without complicating the filtering code, and often results in the fastest code too.
You seem to have a few bugs.
if i + z - indexer < 0 or i + z - indexer > len(data) - 1:If
iandzare0, whereindexeris 1, then you'll have0 + 0 - 1 < 0. This would mean that you'd replace the data in(-1, j),(0, j)and(1, j)to 0. Since 0 and 1 probably do contain data this is just plain wrong.if j + z - indexer < 0 or j + indexer > len(data[0]) - 1: temp.append(0)This removes some data, meaning that the median is shifted. Say you should have
(0, 0, 0, 1, 2, 3), however you removed the first three because of this you'd have(0, 1, 2, 3). Now the median is1rather than0.
Your code would be simpler if you:
- Made a window list, that contained all the indexes that you want to move to.
- Have an if to check if the data in that index is out of bounds.
- If it's out of bounds default to 0.
- If it's not out of bounds use the data.
This could become:
def median_filter(data, filter_size):
temp = []
indexer = filter_size // 2
window = [
(i, j)
for i in range(-indexer, filter_size-indexer)
for j in range(-indexer, filter_size-indexer)
]
index = len(window) // 2
for i in range(len(data)):
for j in range(len(data[0])):
data[i][j] = sorted(
0 if (
min(i+a, j+b) < 0
or len(data) <= i+a
or len(data[0]) <= j+b
) else data[i+a][j+b]
for a, b in window
)[index]
return data