Table of Contents

• 1  Import images
• 2  Setup
• 3  Partial volume example
• 4  Geometric transformation matrix (GTM)
• 5  Partial volume correction (PVC)
• 6  Save
• 7  Visualization
• 8  Discussion
• 9  Citation

from tqdm.notebook import tqdm

import numpy as np
import pandas as pd
import ants

import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter

!tree . -L 1

.
├── Algorithm_GTM_PVC.html
├── Algorithm_GTM_PVC.ipynb
├── img_MR_cropped.nii
├── img_PET_avg.nii
├── img_PET_cropped.nii
├── img_PET_pvc.nii
├── img_ROI_cropped.nii
└── original_3D_data

1 directory, 7 files

Import images

img_MR = ants.image_read('img_MR_cropped.nii')
img_PET = ants.image_read('img_PET_cropped.nii')
img_ROI = ants.image_read('img_ROI_cropped.nii')

arr_MR = np.squeeze(img_MR.numpy())
arr_PET = np.squeeze(img_PET.numpy())
arr_ROI = np.squeeze(img_ROI.numpy())

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, sharey=True, figsize = (15, 40))
ax1.imshow(arr_MR, cmap='gray', interpolation='nearest')
ax2.imshow(arr_PET, cmap='Blues', interpolation='nearest', vmin = -0, vmax = 15)
ax3.imshow(arr_ROI, cmap='jet', interpolation='nearest')
ax1.title.set_text('MR')
ax2.title.set_text('PET')
ax3.title.set_text('ROI')
plt.show()

ls_ROI = np.unique(arr_ROI.flatten())
print(f'{len(ls_ROI)} ROI indice was found: {ls_ROI}')

106 ROI indice was found: [  0.   1.   3.   4.   9.  12.  16.  17.  18.  20.  21.  23.  25.  26.
  27.  28.  29.  34.  38.  39.  40.  42.  43.  45.  47.  48.  49.  50.
  51.  54.  57.  59.  62.  63.  64.  65.  66.  67.  68.  69.  71.  72.
  73.  74.  75.  78.  80.  81.  84.  87.  88.  89.  94.  95.  97.  98.
  99. 100. 101. 102. 103. 104. 106. 108. 110. 115. 116. 119. 123. 124.
 127. 129. 130. 132. 133. 134. 135. 136. 137. 138. 139. 143. 144. 145.
 150. 151. 154. 157. 158. 163. 164. 165. 166. 167. 168. 169. 170. 171.
 172. 173. 175. 177. 179. 184. 185. 188.]

Setup

# setup 
FWHM = 5
sigma = 5/2.355
print(f'sigma = {sigma}')

sigma = 2.1231422505307855

Partial volume example

# smoothing example using a ROI index 
i = 23
arr_ROI_tmp = arr_ROI.copy()
arr_ROI_tmp[arr_ROI_tmp!=ls_ROI[i]]=0
arr_ROI_tmp[arr_ROI_tmp==ls_ROI[i]]=1
arr_ROI_tmp_smoothed = gaussian_filter(arr_ROI_tmp, sigma = sigma)

fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True, figsize = (10, 15))
ax1.imshow(arr_ROI_tmp, cmap='gray', interpolation='nearest')
ax2.imshow(arr_ROI_tmp_smoothed, cmap='gray', interpolation='nearest')
ax1.title.set_text('Original ROI')
ax2.title.set_text('Original ROI with partial volume effect')
plt.show()

Geometric transformation matrix (GTM)

# generate GTM
GTM = np.zeros((len(ls_ROI), len(ls_ROI)))
print(f'GTM with shape {GTM.shape} was initialized.')

arr_ROI_flatten = arr_ROI.flatten()
for i in tqdm(np.arange(0,len(ls_ROI)), desc = 'calculating weight according to ROI indice'):
    # smoothing
    arr_ROI_tmp = arr_ROI.copy()
    arr_ROI_tmp[arr_ROI_tmp!=ls_ROI[i]]=0
    arr_ROI_tmp[arr_ROI_tmp==ls_ROI[i]]=1
    arr_ROI_tmp_smoothed = gaussian_filter(arr_ROI_tmp, sigma = sigma)
    arr_ROI_tmp_smoothed_flatten = arr_ROI_tmp_smoothed.flatten()/arr_ROI_tmp_smoothed.sum()
    # collect weight
    for j in np.arange(0,len(ls_ROI)):
        w = arr_ROI_tmp_smoothed_flatten[arr_ROI_flatten==ls_ROI[j]].sum()
        GTM[i, j] = w

# drop the first row and column if non-brain region is at the first place in ls_ROI and GTM
ls_ROI = ls_ROI[ls_ROI!=0]
GTM = GTM[1:, 1:]
inverse_GTM = np.linalg.inv(GTM)

# visualized GTM and inverse_GTM
fig, (ax1, ax2) = plt.subplots(1, 2, sharex=True, sharey=True, figsize = (10, 5))
img1 = ax1.imshow(GTM, cmap='gray', interpolation='nearest')
img2 = ax2.imshow(inverse_GTM, cmap='gray', interpolation='nearest')
fig.colorbar(img1, ax=ax1)
fig.colorbar(img2, ax=ax2)
ax1.title.set_text('GTM')
ax2.title.set_text('inverse_GTM')
plt.show()

GTM with shape (106, 106) was initialized.

Partial volume correction (PVC)

# Sampling PET
signal_matrix = np.zeros(len(ls_ROI))
for i in tqdm(np.arange(0,len(ls_ROI)), desc = 'calculating raw avg signal according to ROI indice'):
    signal_matrix[i] = arr_PET[arr_ROI==ls_ROI[i]].mean()
    
# regenerate the PET with avg
arr_PET_avg = np.zeros(arr_PET.shape)
for i in tqdm(np.arange(0,len(ls_ROI)), desc = 'generate avg PET image according to ROI indice'):
    arr_PET_avg[arr_ROI==ls_ROI[i]] = signal_matrix[i]

# regenerate the PET with PVC
arr_PET_pvc = np.zeros(arr_PET.shape)
signal_matrix_pvc = np.matmul(inverse_GTM, signal_matrix)
for i in tqdm(np.arange(0,len(ls_ROI)), desc = 'generate pvc PET image according to ROI indice'):
    arr_PET_pvc[arr_ROI==ls_ROI[i]] = signal_matrix_pvc[i]

# compare histgrame of signal_matrix and signal_matrix_pvc
fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, sharey=True, figsize = (14, 4))
ax1.hist(arr_PET[arr_ROI!=0], bins=50, density=True, range = (-10,30))
ax2.hist(arr_PET_avg[arr_ROI!=0], bins=50, density=True, range = (-10,30))
ax3.hist(arr_PET_pvc[arr_ROI!=0], bins=50, density=True, range = (-10,30))
ax1.title.set_text('arr_PET')
ax2.title.set_text('arr_PET_avg')
ax3.title.set_text('arr_PET_pvc')
plt.show()

Save

# save result
img_PET_avg = img_PET.new_image_like(arr_PET_avg[:, :, np.newaxis])
img_PET_pvc = img_PET.new_image_like(arr_PET_pvc[:, :, np.newaxis])
ants.image_write(img_PET_avg, 'img_PET_avg.nii')
ants.image_write(img_PET_pvc, 'img_PET_pvc.nii')

Visualization

arr_PET[arr_ROI==0] = np.nan
arr_PET_avg[arr_ROI==0] = np.nan
arr_PET_pvc[arr_ROI==0] = np.nan

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, sharex=True, sharey=True, figsize = (15, 5))
img1 = ax1.imshow(arr_PET, cmap='Blues', interpolation='nearest', vmin = -0, vmax = 15)
img2 = ax2.imshow(arr_PET_avg, cmap='Blues', interpolation='nearest', vmin = -0, vmax = 15)
img3 = ax3.imshow(arr_PET_pvc, cmap='Blues', interpolation='nearest', vmin = -0, vmax = 15)
ax1.title.set_text('PET')
ax2.title.set_text('PET_avg')
ax3.title.set_text('PET_pvc')
# add space for colour bar
fig.subplots_adjust(right=0.85)
cbar_ax = fig.add_axes([0.88, 0.23, 0.01, 0.55])
fig.colorbar(img3, cax=cbar_ax)
plt.show()

Discussion

• This is a GTM_PVC algorithm for 2D slice only. Due to small volume for some ROI in this 2D case, the values of these ROI are very vulnerable, which result in extremely large or small value after PVC.
• In 3D case where their complete ROI has more volume that are similar to the other, their value are more reasonable.

Citation

• Rousset OG, Ma Y, Evans AC. Correction for partial volume effects in PET: principle and validation. J Nucl Med. 1998 May;39(5):904-11. PMID: 9591599.

!python3 -m jupyter nbconvert Algorithm_GTM_PVC.ipynb --to html

[NbConvertApp] Converting notebook Algorithm_GTM_PVC.ipynb to html
[NbConvertApp] Writing 834751 bytes to Algorithm_GTM_PVC.html