0
# Simulations and Phantoms
1
2
Signal simulation tools for generating synthetic diffusion MRI data, phantoms, and testing algorithms with known ground truth. DIPY provides comprehensive simulation capabilities for method validation and algorithm development.
3
4
## Capabilities
5
6
### Single Tensor Simulation
7
8
Generate diffusion signals from single tensor models with specified diffusion properties.
9
10
```python { .api }
11
def single_tensor(gtab, S0=1, evals=None, evecs=None, snr=None):
12
"""
13
Generate single tensor diffusion signal.
14
15
Parameters:
16
gtab (GradientTable): gradient table for simulation
17
S0 (float): baseline signal intensity
18
evals (array): tensor eigenvalues [lambda1, lambda2, lambda3]
19
evecs (array): tensor eigenvectors (3x3 matrix)
20
snr (float): signal-to-noise ratio (None for no noise)
21
22
Returns:
23
array: simulated diffusion signal
24
"""
25
26
def tensor(gtab, evals, evecs, S0=1, snr=None):
27
"""
28
Simulate tensor signal with specified parameters.
29
30
Parameters:
31
gtab (GradientTable): diffusion gradient table
32
evals (array): eigenvalues [AD, RD, RD] in mm²/s
33
evecs (array): eigenvector matrix
34
S0 (float): b=0 signal intensity
35
snr (float): signal-to-noise ratio
36
37
Returns:
38
array: tensor-based signal
39
"""
40
41
def design_matrix(gtab):
42
"""
43
Create design matrix for tensor fitting.
44
45
Parameters:
46
gtab (GradientTable): gradient table
47
48
Returns:
49
array: design matrix for linear fitting
50
"""
51
```
52
53
### Multi-Tensor Simulation
54
55
Simulate complex tissue with multiple fiber populations and crossing configurations.
56
57
```python { .api }
58
def multi_tensor(gtab, mevals, angles, fractions, S0=100, snr=None):
59
"""
60
Multi-tensor signal simulation for crossing fibers.
61
62
Parameters:
63
gtab (GradientTable): gradient table
64
mevals (list): eigenvalues for each tensor
65
angles (list): fiber orientations [(theta1, phi1), (theta2, phi2), ...]
66
fractions (list): volume fractions for each tensor (sum to 100)
67
S0 (float): baseline signal
68
snr (float): signal-to-noise ratio
69
70
Returns:
71
array: multi-tensor signal
72
"""
73
74
def two_tensor(gtab, evals1, evals2, angles, fractions, S0=100, snr=None):
75
"""
76
Two-tensor crossing fiber simulation.
77
78
Parameters:
79
gtab (GradientTable): gradient information
80
evals1 (array): first tensor eigenvalues
81
evals2 (array): second tensor eigenvalues
82
angles (tuple): crossing angles (theta1, phi1, theta2, phi2)
83
fractions (list): volume fractions [f1, f2]
84
S0 (float): b=0 signal
85
snr (float): noise level
86
87
Returns:
88
array: two-tensor signal
89
"""
90
91
def multi_tensor_odf(mevals, angles, fractions, sphere, tau=1/4*pi):
92
"""
93
Generate multi-tensor ODF for visualization.
94
95
Parameters:
96
mevals (list): tensor eigenvalues
97
angles (list): fiber angles
98
fractions (list): volume fractions
99
sphere (Sphere): sampling sphere
100
tau (float): diffusion time parameter
101
102
Returns:
103
array: ODF values on sphere
104
"""
105
```
106
107
### Sticks and Ball Model
108
109
Simplified model with stick-like fibers in isotropic background for fast simulation.
110
111
```python { .api }
112
def sticks_and_ball(gtab, d=0.0015, S0=1, angles=[(0, 0)], fractions=[100], snr=None):
113
"""
114
Sticks and ball model simulation.
115
116
Parameters:
117
gtab (GradientTable): gradient table
118
d (float): stick diffusivity (mm²/s)
119
S0 (float): baseline signal
120
angles (list): stick orientations [(theta, phi), ...]
121
fractions (list): volume fractions for each stick
122
snr (float): signal-to-noise ratio
123
124
Returns:
125
tuple: (signal, sticks) simulated signal and stick parameters
126
"""
127
128
def ball_and_stick_signal(gtab, d=0.0015, f=0.5, S0=1, angles=[(0, 0)], snr=None):
129
"""
130
Ball and stick model with isotropic and anisotropic components.
131
132
Parameters:
133
gtab (GradientTable): diffusion gradients
134
d (float): stick diffusivity
135
f (float): stick volume fraction (0-1)
136
S0 (float): baseline signal
137
angles (list): stick directions
138
snr (float): noise level
139
140
Returns:
141
array: ball and stick signal
142
"""
143
```
144
145
### Diffusion Kurtosis Simulation
146
147
Generate signals with non-Gaussian diffusion characteristics for DKI validation.
148
149
```python { .api }
150
def kurtosis_signal(gtab, dt, kt, S0=1, snr=None):
151
"""
152
Simulate diffusion kurtosis signal.
153
154
Parameters:
155
gtab (GradientTable): multi-shell gradient table
156
dt (array): diffusion tensor (6-element vector)
157
kt (array): kurtosis tensor (15-element vector)
158
S0 (float): baseline signal
159
snr (float): signal-to-noise ratio
160
161
Returns:
162
array: kurtosis-based signal
163
"""
164
165
def gaussian_signal(gtab, dt, S0=1, snr=None):
166
"""
167
Generate Gaussian (tensor) signal for comparison.
168
169
Parameters:
170
gtab (GradientTable): gradient table
171
dt (array): diffusion tensor
172
S0 (float): b=0 signal
173
snr (float): noise level
174
175
Returns:
176
array: Gaussian diffusion signal
177
"""
178
```
179
180
### Phantom Generation
181
182
Create realistic phantom datasets with known ground truth for validation.
183
184
```python { .api }
185
class SingleShellPhantom:
186
"""Single-shell diffusion phantom generator."""
187
def __init__(self, gtab, snr=None):
188
"""
189
Initialize phantom generator.
190
191
Parameters:
192
gtab (GradientTable): single-shell gradient table
193
snr (float): signal-to-noise ratio
194
"""
195
196
def create_voxel(self, tensor_params):
197
"""
198
Create single voxel with specified parameters.
199
200
Parameters:
201
tensor_params (dict): tensor parameters
202
203
Returns:
204
array: simulated voxel signal
205
"""
206
207
def create_volume(self, shape, tissue_map):
208
"""
209
Create full volume phantom.
210
211
Parameters:
212
shape (tuple): volume dimensions
213
tissue_map (array): tissue type labels
214
215
Returns:
216
array: 4D phantom volume
217
"""
218
219
class MultiShellPhantom:
220
"""Multi-shell phantom for advanced models."""
221
def __init__(self, gtab, snr=None):
222
"""Initialize multi-shell phantom."""
223
224
def add_crossing_region(self, center, radius, angles, fractions):
225
"""Add crossing fiber region to phantom."""
226
227
def add_single_fiber_region(self, center, radius, tensor_params):
228
"""Add single fiber region."""
229
230
def generate(self):
231
"""Generate complete phantom dataset."""
232
233
def spherical_phantom(gtab, sphere, angles, S0=100, snr=None):
234
"""
235
Create spherical phantom with specified fiber orientations.
236
237
Parameters:
238
gtab (GradientTable): gradient table
239
sphere (Sphere): sampling sphere
240
angles (list): fiber directions
241
S0 (float): baseline signal
242
snr (float): noise level
243
244
Returns:
245
array: spherical phantom data
246
"""
247
```
248
249
### Noise Models
250
251
Realistic noise models for diffusion MRI simulation including Rician and Gaussian noise.
252
253
```python { .api }
254
def add_noise(signal, snr, S0=None, noise_type='rician'):
255
"""
256
Add noise to simulated signal.
257
258
Parameters:
259
signal (array): clean signal
260
snr (float): signal-to-noise ratio
261
S0 (float): reference signal for SNR calculation
262
noise_type (str): noise model ('rician', 'gaussian')
263
264
Returns:
265
array: noisy signal
266
"""
267
268
def rician_noise(signal, sigma):
269
"""
270
Add Rician noise to magnitude signal.
271
272
Parameters:
273
signal (array): magnitude signal
274
sigma (float): noise standard deviation
275
276
Returns:
277
array: Rician-distributed noisy signal
278
"""
279
280
def gaussian_noise(signal, sigma):
281
"""
282
Add Gaussian noise to signal.
283
284
Parameters:
285
signal (array): input signal
286
sigma (float): noise standard deviation
287
288
Returns:
289
array: Gaussian noisy signal
290
"""
291
292
def estimate_sigma(data, mask=None):
293
"""
294
Estimate noise level from background regions.
295
296
Parameters:
297
data (array): diffusion data
298
mask (array): brain mask
299
300
Returns:
301
float: estimated noise standard deviation
302
"""
303
```
304
305
### Ground Truth Generation
306
307
Tools for creating datasets with known ground truth for algorithm validation.
308
309
```python { .api }
310
class GroundTruthDataset:
311
"""Generate datasets with known ground truth."""
312
def __init__(self, shape=(64, 64, 40), gtab=None):
313
"""
314
Initialize ground truth generator.
315
316
Parameters:
317
shape (tuple): volume dimensions
318
gtab (GradientTable): gradient table
319
"""
320
321
def add_bundle(self, streamlines, tensor_params):
322
"""
323
Add fiber bundle with known properties.
324
325
Parameters:
326
streamlines (list): bundle streamlines
327
tensor_params (dict): diffusion parameters
328
"""
329
330
def add_partial_volume(self, region_mask, tissue_fractions):
331
"""Add partial volume effects."""
332
333
def generate_data(self, snr=30):
334
"""
335
Generate complete dataset.
336
337
Returns:
338
tuple: (data, ground_truth_params)
339
"""
340
341
def create_qa_phantom(gtab, num_bundles=3, crossing_angle=60):
342
"""
343
Create quality assurance phantom with crossing bundles.
344
345
Parameters:
346
gtab (GradientTable): gradient table
347
num_bundles (int): number of fiber bundles
348
crossing_angle (float): angle between crossing fibers
349
350
Returns:
351
tuple: (phantom_data, bundle_masks, tensor_params)
352
"""
353
```
354
355
### Usage Examples
356
357
```python
358
# Single tensor simulation
359
from dipy.sims.voxel import single_tensor, multi_tensor
360
from dipy.core.gradients import gradient_table
361
from dipy.core.sphere import get_sphere
362
import numpy as np
363
364
# Create gradient table
365
bvals = np.array([0, 1000, 1000, 2000, 2000])
366
bvecs = np.array([[0, 0, 0], [1, 0, 0], [0, 1, 0], [0, 0, 1], [1, 1, 0]])
367
bvecs = bvecs / np.linalg.norm(bvecs, axis=1, keepdims=True)
368
bvecs[0] = 0 # b=0 has no direction
369
gtab = gradient_table(bvals, bvecs)
370
371
# Simulate single tensor
372
evals = np.array([1.7e-3, 0.3e-3, 0.3e-3]) # AD, RD, RD
373
evecs = np.eye(3) # aligned with x-axis
374
signal = single_tensor(gtab, S0=100, evals=evals, evecs=evecs, snr=30)
375
376
print(f"Simulated signal shape: {signal.shape}")
377
print(f"b=0 signal: {signal[0]:.1f}")
378
print(f"DW signal mean: {signal[1:].mean():.1f}")
379
380
# Multi-tensor crossing simulation
381
from dipy.sims.voxel import sticks_and_ball
382
383
# 90-degree crossing
384
angles = [(0, 0), (90, 0)] # two perpendicular fibers
385
fractions = [50, 50] # equal volume fractions
386
crossing_signal, sticks = sticks_and_ball(
387
gtab, d=0.0015, S0=100, angles=angles, fractions=fractions, snr=25
388
)
389
390
# Create phantom volume
391
from dipy.sims.phantom import SingleShellPhantom
392
393
phantom = SingleShellPhantom(gtab, snr=30)
394
395
# Define tissue regions
396
shape = (32, 32, 20)
397
phantom_data = np.zeros(shape + (len(gtab.bvals),))
398
399
# Add single fiber region
400
center = (16, 16, 10)
401
for i in range(shape[0]):
402
for j in range(shape[1]):
403
for k in range(shape[2]):
404
# Distance from center
405
dist = np.sqrt((i-center[0])**2 + (j-center[1])**2 + (k-center[2])**2)
406
if dist < 8: # Within sphere
407
if dist < 4: # Inner core - single fiber
408
signal = single_tensor(gtab, S0=100, evals=evals, evecs=evecs, snr=30)
409
else: # Outer shell - partial volume
410
mixed_signal = 0.7 * single_tensor(gtab, S0=100, evals=evals, evecs=evecs) + \
411
0.3 * single_tensor(gtab, S0=100, evals=[0.7e-3]*3, evecs=evecs)
412
mixed_signal = add_noise(mixed_signal, snr=30)
413
phantom_data[i, j, k] = signal if dist < 4 else mixed_signal
414
415
print(f"Phantom shape: {phantom_data.shape}")
416
417
# Validation with known ground truth
418
from dipy.reconst.dti import TensorModel
419
420
# Fit tensor model to simulated data
421
tensor_model = TensorModel(gtab)
422
tensor_fit = tensor_model.fit(phantom_data)
423
424
# Compare with ground truth
425
estimated_fa = tensor_fit.fa
426
ground_truth_fa = 0.8 # Expected FA for our simulation
427
428
print(f"Ground truth FA: {ground_truth_fa:.3f}")
429
print(f"Estimated FA (center): {estimated_fa[16, 16, 10]:.3f}")
430
print(f"Estimation error: {abs(estimated_fa[16, 16, 10] - ground_truth_fa):.3f}")
431
432
# Noise analysis
433
from dipy.sims.voxel import add_noise, estimate_sigma
434
435
# Test different noise levels
436
snr_levels = [10, 20, 30, 50, 100]
437
clean_signal = single_tensor(gtab, S0=100, evals=evals, evecs=evecs)
438
439
for snr in snr_levels:
440
noisy_signal = add_noise(clean_signal, snr, noise_type='rician')
441
sigma_estimated = estimate_sigma(noisy_signal.reshape(1, 1, 1, -1))
442
sigma_true = 100 / snr # True noise level
443
print(f"SNR {snr}: True σ={sigma_true:.2f}, Estimated σ={sigma_estimated:.2f}")
444
```