1. analysis-clinic
  2. nipy_image_slicing.ipynb

Slicing images with nipy

Sometimes it is useful to slice out image edges or subsample images. For example, your functional images might have a lot of empty space around the brain that will take up memory during processing, or your structural image might contain too much of the subject's neck, or the resolution may be to high.

You can always load the image data and slice away the parts of the image data array that you don't want, but then you will have to work out the effect this should have on the image affine (see http://nipy.org/nibabel/coordinate_systems.html).

Luckily nipy has the ability to slice images while working out the effect on the affine matrix for you.

In [1]:
# Import the numerical and plotting libraries
import numpy as np
import matplotlib.pyplot as plt
In [2]:
# Show plots inside the notebook
%matplotlib inline
In [3]:
# - set gray colormap and nearest neighbor interpolation by default
plt.rcParams['image.cmap'] = 'gray'
plt.rcParams['image.interpolation'] = 'nearest'
In [4]:
# - make numpy print numbers to 4dp only
np.set_printoptions(precision=4, suppress=True)

We load our example functional image again:

In [5]:
import nipy
In [6]:
img = nipy.load_image('ds114_sub009_t2r1.nii')
img.coordmap
Out[6]:
AffineTransform(
   function_domain=CoordinateSystem(coord_names=('i', 'j', 'k', 't'), name='voxels', coord_dtype=float64),
   function_range=CoordinateSystem(coord_names=('scanner-x=L->R', 'scanner-y=P->A', 'scanner-z=I->S', 't'), name='scanner', coord_dtype=float64),
   affine=array([[  -4.    ,    0.    ,    0.    ,    0.    ,  124.244 ],
                 [   0.    ,    3.9345,    0.7207,    0.    , -103.4497],
                 [   0.    ,   -0.7207,    3.9346,    0.    ,  -33.4929],
                 [   0.    ,    0.    ,    0.    ,    2.5   ,    0.    ],
                 [   0.    ,    0.    ,    0.    ,    0.    ,    1.    ]])
)

The functional data has a lot of empty space around the brain:

In [7]:
data = img.get_data()
data.shape
Out[7]:
(64, 64, 30, 173)
In [8]:
plt.imshow(data[:, :, 15, 0])
Out[8]:
<matplotlib.image.AxesImage at 0x108bd6b10>
In [9]:
plt.imshow(data[32, :, :, 0])  # A sagittal slice
Out[9]:
<matplotlib.image.AxesImage at 0x108e3cc90>

Maybe we would look to slice off the empty space to left and right and front and back of the brain.

Maybe we also want to halve the resolution in the bottom-top direction. We do that by slicing the image:

In [10]:
small_img = img[15:50, 8:55, ::2]

The slicing we did on the image shows up when we get the data from the image:

In [11]:
small_data = small_img.get_data()
small_data.shape
Out[11]:
(35, 47, 15, 173)
In [12]:
plt.imshow(small_data[:, :, 7, 0])
Out[12]:
<matplotlib.image.AxesImage at 0x108f0d250>

The sagittal view now looks more squashed up because we've removed every second slice with the image slicing:

In [13]:
plt.imshow(small_data[17, :, :, 0])
Out[13]:
<matplotlib.image.AxesImage at 0x109047f10>

Nipy adjusted the image affine to the slicing. Notice that the last column of the affine, reflecting the different translations to go from voxel 0, 0, 0 to the position in millimeters (voxel 0, 0, 0 is a different voxel after the slicing). Also, the second to last column contains the size of the voxel in the z (bottom to top) direction, and this is now around 8 instead of 4 previously. That is because we threw away every second slice, so the remaining slices are twice the distance apart from each other.

In [14]:
small_img.affine
Out[14]:
array([[ -4.    ,   0.    ,   0.    ,   0.    ,  64.244 ],
       [  0.    ,   3.9345,   1.4414,   0.    , -71.9734],
       [  0.    ,  -0.7207,   7.8691,   0.    , -39.2586],
       [  0.    ,   0.    ,   0.    ,   2.5   ,   0.    ],
       [  0.    ,   0.    ,   0.    ,   0.    ,   1.    ]])

In the notebook on diagnostics, we saw that the first few scans were different from the rest. This is the T1 effect, where the signal is largest when the radiofrequency stimulation first starts. We can drop the first four scans with slicing, to get rid of this effect.

In [15]:
img = img[..., 4:]
In [16]:
img.shape
Out[16]:
(64, 64, 30, 169)

We can save this image back out to disk for further processing by other packages:

In [17]:
nipy.save_image(img, 'fixed.nii')

/Users/mb312/Library/Python/2.7/lib/python/site-packages/nipy/io/files.py:145: FutureWarning: Default `strict` currently False; this will change to True in a future version of nipy
  ni_img = nipy2nifti(img, data_dtype = io_dtype)
Out[17]:
Image(
  data=memmap([[[[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 8, 10, ..., 10,  7],
                [11, 10, ...,  8, 10],
                ..., 
                [19, 11, ..., 20, 21],
                [28, 23, ..., 19, 35]],
       
               ..., 
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]]],
       
       
              [[[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[10, 10, ..., 10,  7],
                [11, 11, ...,  9, 11],
                ..., 
                [17, 19, ..., 24, 14],
                [24, 41, ..., 21, 22]],
       
               ..., 
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]]],
       
       
              ..., 
              [[[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 1,  1, ...,  2,  1],
                [ 1,  1, ...,  1,  1],
                ..., 
                [22, 19, ..., 12, 19],
                [14, 14, ..., 28, 23]],
       
               ..., 
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]]],
       
       
              [[[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 1,  1, ...,  2,  1],
                [ 1,  1, ...,  0,  1],
                ..., 
                [14, 15, ..., 17, 18],
                [20, 12, ..., 16, 13]],
       
               ..., 
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]],
       
               [[ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0],
                ..., 
                [ 0,  0, ...,  0,  0],
                [ 0,  0, ...,  0,  0]]]], dtype=int16),
  coordmap=AffineTransform(
            function_domain=CoordinateSystem(coord_names=('i', 'j', 'k', 't'), name='input-slice', coord_dtype=float64),
            function_range=CoordinateSystem(coord_names=('scanner-x=L->R', 'scanner-y=P->A', 'scanner-z=I->S', 't'), name='scanner', coord_dtype=float64),
            affine=array([[  -4.      ,    0.      ,    0.      ,    0.      ,  124.244003],
                          [   0.      ,    3.934535,    0.720723,    0.      , -103.449684],
                          [   0.      ,   -0.72072 ,    3.934551,    0.      ,  -33.492859],
                          [   0.      ,    0.      ,    0.      ,    2.5     ,   10.      ],
                          [   0.      ,    0.      ,    0.      ,    0.      ,    1.      ]])
         ))