0
# Image Processing and Filtering
1

2
Denoising algorithms, filtering methods, and preprocessing tools for diffusion MRI data quality improvement. DIPY provides state-of-the-art methods for enhancing signal quality while preserving important diffusion information.
3

4
## Capabilities
5

6
### Denoising Algorithms
7

8
Advanced denoising methods specifically designed for diffusion MRI data.
9

10
```python { .api }
11
def nlmeans(arr, sigma=None, mask=None, patch_radius=1, block_radius=5, rician=True, num_threads=None):
12
    """
13
    Non-local means denoising for diffusion MRI.
14
    
15
    Parameters:
16
        arr (array): 4D diffusion data (x, y, z, diffusion_directions)
17
        sigma (float): noise standard deviation (estimated if None)
18
        mask (array): binary mask for processing region
19
        patch_radius (int): radius of comparison patches
20
        block_radius (int): radius of search neighborhood  
21
        rician (bool): assume Rician noise distribution
22
        num_threads (int): number of parallel threads
23
        
24
    Returns:
25
        array: denoised diffusion data
26
    """
27

28
def localpca(arr, sigma=None, mask=None, patch_radius=2, pca_method='eig', tau_factor=2.3, return_sigma=False, out_dtype=None):
29
    """
30
    Local PCA-based denoising using Marchenko-Pastur distribution.
31
    
32
    Parameters:
33
        arr (array): 4D diffusion data
34
        sigma (float): noise level (auto-estimated if None)
35
        mask (array): processing mask
36
        patch_radius (int): local patch radius
37
        pca_method (str): PCA computation method ('eig', 'svd')
38
        tau_factor (float): threshold factor for eigenvalue cutoff
39
        return_sigma (bool): return estimated noise level
40
        out_dtype (dtype): output data type
41
        
42
    Returns:
43
        array or tuple: denoised data, optionally with sigma estimate
44
    """
45

46
def mppca(arr, mask=None, patch_radius=2, pca_method='eig', return_sigma=False, out_dtype=None):
47
    """
48
    Marchenko-Pastur PCA denoising.
49
    
50
    Parameters:
51
        arr (array): 4D diffusion data
52
        mask (array): binary mask
53
        patch_radius (int): patch size for local PCA
54
        pca_method (str): eigenvalue computation method
55
        return_sigma (bool): return noise estimate
56
        out_dtype (dtype): output data type
57
        
58
    Returns:
59
        array or tuple: denoised data and optional sigma
60
    """
61

62
def patch2self(data, bvals, patch_radius=0, model='ols', b0_threshold=50, alpha=1.0, verbose=False):
63
    """
64
    Patch2Self self-supervised denoising.
65
    
66
    Parameters:
67
        data (array): 4D diffusion data
68
        bvals (array): b-values for diffusion encoding
69
        patch_radius (int): spatial patch radius
70
        model (str): regression model ('ols', 'ridge', 'lasso')
71
        b0_threshold (float): threshold for b=0 identification
72
        alpha (float): regularization parameter
73
        verbose (bool): print progress information
74
        
75
    Returns:
76
        array: self-supervised denoised data
77
    """
78
```
79

80
### Motion Correction
81

82
Algorithms for correcting subject motion during diffusion MRI acquisition.
83

84
```python { .api }
85
def motion_correction(data, gtab, *, affine=None, b0_ref=0, pipeline=None, static_mask=None):
86
    """
87
    Motion correction for DWI datasets (inter-volume motion).
88
    
89
    Parameters:
90
        data (array/str): 4D DWI data or file path
91
        gtab (GradientTable): gradient table
92
        affine (array): voxel-to-world transformation
93
        b0_ref (int): reference B0 volume index
94
        pipeline (list): registration transformation sequence
95
        static_mask (array): mask for reference volume
96
        
97
    Returns:
98
        tuple: (corrected_img, affine_array) motion-corrected data and transforms
99
    """
100

101
def register_dwi_series(data, gtab, *, affine=None, b0_ref=0, pipeline=None, static_mask=None):
102
    """
103
    Register DWI series to mean B0 for motion correction.
104
    
105
    Parameters:
106
        data (array/str): DWI data
107
        gtab (GradientTable): gradient information
108
        affine (array): affine transformation
109
        b0_ref (int): reference B0 index
110
        pipeline (list): registration steps (['center_of_mass', 'translation', 'rigid', 'affine'])
111
        static_mask (array): processing mask
112
        
113
    Returns:
114
        tuple: (registered_img, transformation_matrices)
115
    """
116
```
117

118
### Gibbs Ringing Removal
119

120
Methods for removing Gibbs ringing artifacts from MRI images.
121

122
```python { .api }
123
def gibbs_removal(vol, slice_axis=2, n_points=3, num_processes=1):
124
    """
125
    Remove Gibbs ringing artifacts from MRI data.
126
    
127
    Parameters:
128
        vol (array): 3D or 4D MRI volume
129
        slice_axis (int): axis along which to apply correction
130
        n_points (int): number of oscillations to correct
131
        num_processes (int): parallel processing threads
132
        
133
    Returns:
134
        array: corrected volume with reduced Gibbs ringing
135
    """
136
```
137

138
### Eddy Current Correction
139

140
Correction for eddy current-induced distortions in diffusion data.
141

142
```python { .api }
143
def eddy_correct(dwi, gtab, *, ref_b0=0, num_threads=1, slice_drop_correction=True):
144
    """
145
    Eddy current correction for diffusion data.
146
    
147
    Parameters:
148
        dwi (array): 4D diffusion-weighted images
149
        gtab (GradientTable): gradient table
150
        ref_b0 (int): reference B0 volume
151
        num_threads (int): number of threads
152
        slice_drop_correction (bool): correct for slice dropouts
153
        
154
    Returns:
155
        array: eddy current corrected data
156
    """
157
```
158

159
### Bias Field Correction
160

161
Methods for correcting intensity inhomogeneity (bias fields) in MRI.
162

163
```python { .api }
164
def n4_bias_correction(data, mask=None, num_iterations=50, convergence_threshold=0.001, bspline_grid_resolution=(1, 1, 1)):
165
    """
166
    N4 bias field correction algorithm.
167
    
168
    Parameters:
169
        data (array): 3D or 4D image data
170
        mask (array): binary processing mask
171
        num_iterations (int): maximum iterations
172
        convergence_threshold (float): convergence criteria
173
        bspline_grid_resolution (tuple): B-spline knot spacing
174
        
175
    Returns:
176
        tuple: (corrected_data, bias_field) corrected image and estimated bias
177
    """
178
```
179

180
### Noise Estimation
181

182
Methods for estimating noise characteristics in diffusion MRI data.
183

184
```python { .api }
185
def estimate_sigma(arr, disable_background_masking=False, N=0):
186
    """
187
    Estimate noise standard deviation in DWI data.
188
    
189
    Parameters:
190
        arr (array): 4D diffusion data or 3D single volume
191
        disable_background_masking (bool): disable automatic background detection
192
        N (int): number of coils for Rician correction
193
        
194
    Returns:
195
        float or array: estimated noise standard deviation
196
    """
197

198
def piesno(data, N=1, alpha=0.01, l=100, itermax=100, eps=1e-5, return_mask=False):
199
    """
200
    Probabilistic Identification and Estimation of Noise (PIESNO).
201
    
202
    Parameters:
203
        data (array): diffusion MRI data
204
        N (int): number of receiver coils
205
        alpha (float): probability of false detection
206
        l (int): number of background voxels to select
207
        itermax (int): maximum iterations
208
        eps (float): convergence tolerance
209
        return_mask (bool): return background mask
210
        
211
    Returns:
212
        float or tuple: noise standard deviation, optionally with mask
213
    """
214
```
215

216
### Gradient Non-linearity Correction
217

218
Correction for gradient non-linearities in diffusion encoding.
219

220
```python { .api }
221
def gradient_nonlinearity_correction(dwi, gtab, grad_dev=None, scanner='siemens'):
222
    """
223
    Correct for gradient non-linearity effects.
224
    
225
    Parameters:
226
        dwi (array): diffusion-weighted images
227
        gtab (GradientTable): gradient table
228
        grad_dev (array): gradient deviation tensor field
229
        scanner (str): scanner manufacturer
230
        
231
    Returns:
232
        array: corrected diffusion data
233
    """
234
```
235

236
### Signal Drift Correction
237

238
Methods for correcting signal drift during long acquisitions.
239

240
```python { .api }
241
def signal_drift_correction(data, gtab, *, b0_threshold=50, method='linear'):
242
    """
243
    Correct for signal drift over acquisition time.
244
    
245
    Parameters:
246
        data (array): 4D diffusion data
247
        gtab (GradientTable): gradient table with acquisition timing
248
        b0_threshold (float): threshold for b=0 identification
249
        method (str): drift correction method ('linear', 'polynomial')
250
        
251
    Returns:
252
        array: drift-corrected diffusion data
253
    """
254
```
255

256
### Usage Examples
257

258
```python
259
# Comprehensive preprocessing pipeline
260
from dipy.denoise import localpca, nlmeans, patch2self
261
from dipy.denoise.gibbs import gibbs_removal
262
from dipy.denoise.noise_estimate import estimate_sigma, piesno
263
from dipy.align import motion_correction
264
from dipy.data import read_stanford_hardi
265
import numpy as np
266

267
# Load example data
268
img, gtab = read_stanford_hardi()
269
data = img.get_fdata()
270

271
print(f"Original data shape: {data.shape}")
272

273
# Step 1: Gibbs ringing removal
274
print("Removing Gibbs ringing artifacts...")
275
data_gibbs = gibbs_removal(data, slice_axis=2)
276

277
# Step 2: Noise estimation
278
print("Estimating noise level...")
279
sigma = estimate_sigma(data_gibbs)
280
print(f"Estimated noise sigma: {sigma:.4f}")
281

282
# Alternative: PIESNO noise estimation
283
sigma_piesno, mask = piesno(data_gibbs, N=1, return_mask=True)
284
print(f"PIESNO sigma estimate: {sigma_piesno:.4f}")
285

286
# Step 3: Denoising with MP-PCA
287
print("Applying MP-PCA denoising...")
288
from dipy.denoise import mppca
289
data_denoised = mppca(data_gibbs, patch_radius=2)
290

291
# Alternative: Non-local means denoising
292
data_nlm = nlmeans(data_gibbs, sigma=sigma, patch_radius=1, block_radius=2)
293

294
# Alternative: Patch2Self denoising
295
data_p2s = patch2self(data_gibbs, gtab.bvals, patch_radius=0, model='ols')
296

297
# Step 4: Motion correction
298
print("Applying motion correction...")
299
data_corrected, transforms = motion_correction(data_denoised, gtab)
300

301
print(f"Applied {transforms.shape[-1]} motion correction transforms")
302

303
# Check improvement in data quality
304
original_std = np.std(data[data > 0])
305
denoised_std = np.std(data_corrected.get_fdata()[data_corrected.get_fdata() > 0])
306
print(f"Noise reduction: {original_std:.4f} -> {denoised_std:.4f}")
307

308
# Advanced denoising comparison
309
print("\nComparing denoising methods...")
310

311
# Local PCA with return sigma
312
data_lpca, sigma_lpca = localpca(data_gibbs, patch_radius=2, return_sigma=True)
313
print(f"Local PCA estimated sigma: {sigma_lpca:.4f}")
314

315
# Non-local means with mask
316
mask = data_gibbs[..., 0] > np.percentile(data_gibbs[..., 0], 50)
317
data_nlm_masked = nlmeans(data_gibbs, mask=mask, patch_radius=1, 
318
                         block_radius=3, rician=True)
319

320
# Save results for comparison
321
from dipy.io.image import save_nifti
322

323
save_nifti('data_original.nii.gz', data, img.affine)
324
save_nifti('data_denoised.nii.gz', data_denoised, img.affine)
325
save_nifti('data_motion_corrected.nii.gz', data_corrected.get_fdata(), 
326
           data_corrected.affine)
327

328
print("Preprocessing pipeline completed successfully!")
329

330
# Noise estimation in practice
331
print("\nNoise estimation methods comparison:")
332

333
# Background-based estimation
334
sigma_bg = estimate_sigma(data, disable_background_masking=False)
335

336
# Manual background masking
337
background_mask = data[..., 0] < np.percentile(data[..., 0], 10)
338
sigma_manual = np.std(data[background_mask, :])
339

340
print(f"Auto background sigma: {sigma_bg:.4f}")
341
print(f"Manual background sigma: {sigma_manual:.4f}")
342
print(f"PIESNO sigma: {sigma_piesno:.4f}")
343
```