Skip to content
This repository has been archived by the owner on Dec 20, 2024. It is now read-only.

Commit

Permalink
enh: create new filtering submodule
Browse files Browse the repository at this point in the history
  • Loading branch information
oesteban committed Apr 10, 2024
1 parent 6128b37 commit 9753f7d
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 41 deletions.
64 changes: 64 additions & 0 deletions src/eddymotion/data/filtering.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
# emacs: -*- mode: python; py-indent-offset: 4; indent-tabs-mode: nil -*-
# vi: set ft=python sts=4 ts=4 sw=4 et:
#
# Copyright 2022 The NiPreps Developers <[email protected]>
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#
# We support and encourage derived works from this project, please read
# about our expectations at
#
# https://www.nipreps.org/community/licensing/
#
"""Filtering data."""


def advanced_clip(
data, p_min=35, p_max=99.98, nonnegative=True, dtype="int16", invert=False
):
"""
Remove outliers at both ends of the intensity distribution and fit into a given dtype.
This interface tries to emulate ANTs workflows' massaging that truncate images into
the 0-255 range, and applies percentiles for clipping images.
For image registration, normalizing the intensity into a compact range (e.g., uint8)
is generally advised.
To more robustly determine the clipping thresholds, spikes are removed from data with
a median filter.
Once the thresholds are calculated, the denoised data are thrown away and the thresholds
are applied on the original image.
"""
import numpy as np
from scipy import ndimage
from skimage.morphology import ball

# Calculate stats on denoised version, to preempt outliers from biasing
denoised = ndimage.median_filter(data, footprint=ball(3))

a_min = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_min)
a_max = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_max)

# Clip and cast
data = np.clip(data, a_min=a_min, a_max=a_max)
data -= data.min()
data /= data.max()

if invert:
data = 1.0 - data

if dtype in ("uint8", "int16"):
data = np.round(255 * data).astype(dtype)

return data
45 changes: 4 additions & 41 deletions src/eddymotion/estimator.py
Original file line number Diff line number Diff line change
Expand Up @@ -188,49 +188,12 @@ def fit(
return dwdata.em_affines


def _advanced_clip(data, p_min=35, p_max=99.98, nonnegative=True, dtype="int16", invert=False):
"""
Remove outliers at both ends of the intensity distribution and fit into a given dtype.
This interface tries to emulate ANTs workflows' massaging that truncate images into
the 0-255 range, and applies percentiles for clipping images.
For image registration, normalizing the intensity into a compact range (e.g., uint8)
is generally advised.
To more robustly determine the clipping thresholds, spikes are removed from data with
a median filter.
Once the thresholds are calculated, the denoised data are thrown away and the thresholds
are applied on the original image.
"""
import numpy as np
from scipy import ndimage
from skimage.morphology import ball

# Calculate stats on denoised version, to preempt outliers from biasing
denoised = ndimage.median_filter(data, footprint=ball(3))

a_min = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_min)
a_max = np.percentile(denoised[denoised > 0] if nonnegative else denoised, p_max)

# Clip and cast
data = np.clip(data, a_min=a_min, a_max=a_max)
data -= data.min()
data /= data.max()

if invert:
data = 1.0 - data

if dtype in ("uint8", "int16"):
data = np.round(255 * data).astype(dtype)

return data


def _to_nifti(data, affine, filename, clip=True):
data = np.squeeze(data)
if clip:
data = _advanced_clip(data)
from eddymotion.data.filtering import advanced_clip

data = advanced_clip(data)
nii = nb.Nifti1Image(
data,
affine,
Expand Down Expand Up @@ -312,7 +275,7 @@ def _prepare_kwargs(dwdata, kwargs):
kwargs["mask"] = dwdata.brainmask

if hasattr(dwdata, "bzero") and dwdata.bzero is not None:
kwargs["S0"] = _advanced_clip(dwdata.bzero)
kwargs["S0"] = dwdata.bzero

if hasattr(dwdata, "gradients"):
kwargs["gtab"] = dwdata.gradients
Expand Down

0 comments on commit 9753f7d

Please sign in to comment.