Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Improve median filter #150

Merged
merged 2 commits into from
Nov 29, 2024
Merged

Improve median filter #150

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
01a8e0f
Improve median filter
cgohlke Nov 28, 2024
2794806
Revert changes to introductory tutorial
cgohlke Nov 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
137 changes: 66 additions & 71 deletions src/phasorpy/_phasorpy.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2041,72 +2041,66 @@ def _signal_denoise_vector(
###############################################################################
# Filtering functions

cdef float_t _quickselect(
float_t *arr,
ssize_t left,
ssize_t right,
const ssize_t k
) noexcept nogil:
"""Return k-th smallest array value using Quickselect algorithm."""

cdef float_t _median(float_t *values, const ssize_t size) noexcept nogil:
"""Return median of array values using Quickselect algorithm."""
cdef:
ssize_t i, pivot_index, pivot_new_index
ssize_t i, pivot_index, pivot_index_new
ssize_t left = 0
ssize_t right = size - 1
ssize_t middle = size // 2
float_t pivot_value, temp

if size % 2 == 0:
middle -= 1 # Quickselect sorts on right

while left <= right:
pivot_index = left + (right - left) // 2
pivot_value = arr[pivot_index]
temp = arr[pivot_index]
arr[pivot_index] = arr[right]
arr[right] = temp
pivot_new_index = left
pivot_value = values[pivot_index]
temp = values[pivot_index]
values[pivot_index] = values[right]
values[right] = temp
pivot_index_new = left
for i in range(left, right):
if arr[i] < pivot_value:
temp = arr[i]
arr[i] = arr[pivot_new_index]
arr[pivot_new_index] = temp
pivot_new_index += 1
temp = arr[right]
arr[right] = arr[pivot_new_index]
arr[pivot_new_index] = temp

if pivot_new_index == k:
return arr[k]
if pivot_new_index < k:
left = pivot_new_index + 1
if values[i] < pivot_value:
temp = values[i]
values[i] = values[pivot_index_new]
values[pivot_index_new] = temp
pivot_index_new += 1
temp = values[right]
values[right] = values[pivot_index_new]
values[pivot_index_new] = temp

if pivot_index_new == middle:
if size % 2 == 0:
return (values[middle] + values[middle + 1]) / <float_t> 2.0
return values[middle]
if pivot_index_new < middle:
left = pivot_index_new + 1
else:
right = pivot_new_index - 1
right = pivot_index_new - 1

return arr[k]
return values[middle] # unreachable code?
cgohlke marked this conversation as resolved.
Show resolved Hide resolved


cdef inline float_t _median(float_t *values, const ssize_t n) noexcept nogil:
"""Return median of array of values."""
if n % 2 == 0:
return (
(_quickselect(values, 0, n - 1, n // 2 - 1) +
_quickselect(values, 0, n - 1, n // 2)) / <float_t> 2.0
)

return _quickselect(values, 0, n - 1, n // 2)


def _apply_2d_median_filter(
def _median_filter_2d(
float_t[:, :] image,
float_t[:, ::1] filtered_image,
const ssize_t kernel_size,
const int repeat=1,
const int num_threads=1,
):
"""Apply a 2D median filter ignoring NaN."""
"""Apply 2D median filter ignoring NaN."""
cdef:
ssize_t rows = image.shape[0]
ssize_t cols = image.shape[1]
ssize_t k = kernel_size // 2
ssize_t i, j, di, dj, ki, kj, valid_count
ssize_t i, j, r, di, dj, ki, kj, valid_count
float_t element
float_t *kernel

if kernel_size <= 0:
raise ValueError("kernel_size must be greater than 0")
raise ValueError('kernel_size must be greater than 0')

with nogil, parallel(num_threads=num_threads):

Expand All @@ -2115,34 +2109,35 @@ def _apply_2d_median_filter(
)
if kernel == NULL:
with gil:
raise MemoryError("Unable to allocate memory for kernel")
raise MemoryError('failed to allocate kernel')

for i in prange(rows):
for j in range(cols):
if isnan(image[i, j]):
filtered_image[i, j] = <float_t> NAN
continue
valid_count = 0
for di in range(kernel_size):
ki = i - k + di
if ki < 0:
ki = 0
elif ki >= rows:
ki = rows - 1
for dj in range(kernel_size):
kj = j - k + dj
if kj < 0:
kj = 0
elif kj >= cols:
kj = cols - 1
element = image[ki, kj]
if not isnan(element):
kernel[valid_count] = element
valid_count = valid_count + 1
filtered_image[i, j] = _median(kernel, valid_count)
for r in range(repeat):
for i in prange(rows):
for j in range(cols):
if isnan(image[i, j]):
filtered_image[i, j] = <float_t> NAN
continue
valid_count = 0
for di in range(kernel_size):
ki = i - k + di
if ki < 0:
ki = 0
elif ki >= rows:
ki = rows - 1
for dj in range(kernel_size):
kj = j - k + dj
if kj < 0:
kj = 0
elif kj >= cols:
kj = cols - 1
element = image[ki, kj]
if not isnan(element):
kernel[valid_count] = element
valid_count = valid_count + 1
filtered_image[i, j] = _median(kernel, valid_count)

for i in prange(rows):
for j in range(cols):
image[i, j] = filtered_image[i, j]

free(kernel)

for i in prange(rows):
for j in range(cols):
image[i, j] = filtered_image[i, j]
Loading