We show how to apply Diffusion Spectrum Imaging [Wedeen08] to diffusion MRI datasets of Cartesian keyhole diffusion gradients.

First import the necessary modules:

```
from dipy.data import fetch_taiwan_ntu_dsi, read_taiwan_ntu_dsi, get_sphere
from dipy.reconst.dsi import DiffusionSpectrumModel
```

Download and read the data for this tutorial.

```
fetch_taiwan_ntu_dsi()
img, gtab = read_taiwan_ntu_dsi()
```

img contains a nibabel Nifti1Image object (data) and gtab contains a GradientTable object (gradient information e.g. b-values). For example to read the b-values it is possible to write print(gtab.bvals).

Load the raw diffusion data and the affine.

```
data = img.get_data()
print('data.shape (%d, %d, %d, %d)' % data.shape)
```

data.shape `(96, 96, 60, 203)`

This dataset has anisotropic voxel sizes, therefore reslicing is necessary.

```
affine = img.affine
```

Read the voxel size from the image header.

```
voxel_size = img.header.get_zooms()[:3]
```

Instantiate the Model and apply it to the data.

```
dsmodel = DiffusionSpectrumModel(gtab)
```

Lets just use one slice only from the data.

```
dataslice = data[:, :, data.shape[2] // 2]
dsfit = dsmodel.fit(dataslice)
```

Load an odf reconstruction sphere

```
sphere = get_sphere('repulsion724')
```

Calculate the ODFs with this specific sphere

```
ODF = dsfit.odf(sphere)
print('ODF.shape (%d, %d, %d)' % ODF.shape)
```

ODF.shape `(96, 96, 724)`

In a similar fashion it is possible to calculate the PDFs of all voxels in one call with the following way

```
PDF = dsfit.pdf()
print('PDF.shape (%d, %d, %d, %d, %d)' % PDF.shape)
```

PDF.shape `(96, 96, 17, 17, 17)`

We see that even for a single slice this PDF array is close to 345 MBytes so we really have to be careful with memory usage when use this function with a full dataset.

The simple solution is to generate/analyze the ODFs/PDFs by iterating through each voxel and not store them in memory if that is not necessary.

```
from dipy.core.ndindex import ndindex
for index in ndindex(dataslice.shape[:2]):
pdf = dsmodel.fit(dataslice[index]).pdf()
```

If you really want to save the PDFs of a full dataset on the disc we recommend
using memory maps (`numpy.memmap`

) but still have in mind that even if you do
that for example for a dataset of volume size `(96, 96, 60)`

you will need about
2.5 GBytes which can take less space when reasonable spheres (with < 1000 vertices)
are used.

Let’s now calculate a map of Generalized Fractional Anisotropy (GFA) [Tuch04] using the DSI ODFs.

```
from dipy.reconst.odf import gfa
GFA = gfa(ODF)
import matplotlib.pyplot as plt
fig_hist, ax = plt.subplots(1)
ax.set_axis_off()
plt.imshow(GFA.T)
plt.savefig('dsi_gfa.png', bbox_inches='tight', origin='lower', cmap='gray')
```

See also Calculate DSI-based scalar maps for calculating different types of DSI maps.

- Wedeen08
Wedeen et al., Diffusion spectrum magnetic resonance imaging (DSI) tractography of crossing fibers, Neuroimage, vol 41, no 4, 1267-1277, 2008.

- Tuch04
Tuch, D.S, Q-ball imaging, MRM, vol 52, no 6, 1358-1372, 2004.

Example source code

You can download `the full source code of this example`

. This same script is also included in the dipy source distribution under the `doc/examples/`

directory.