"Dense Amygdala": Extensive Complex-valued Functional MRI of the Ventral and Medial Temporal Lobe during Passive Movie Watching in Three Individuals.

1. """
2. Class for estimating the respiratory waveform from the phase channel of an EPI timeseries
3. - Developed specifically for the Dense Amygdala project, but more widely useful for complex-valued fMRI
4. - Required temporally unwrapped complex phase difference BOLD timeseries with critical or better temporal sampling of respiration
5. - Do NOT use Laplacian phase unwrapping which removes most of the respiration phase modulation
7. AUTHOR : Mike Tyszka
8. PLACE : Caltech Brain Imaging Center
9. DATES : 2024-12-05 JMT Build from ideas in Calhoun complex phase fMRI
10. """
12. import json
13. import os.path as op
14. import numpy as np
15. import pandas as pd
16. import nibabel as nib
17. import scipy as sp
18. from sklearn.decomposition import PCA
19. from .bidsphysio import save_bids_physio
22. class PhaseRespWave:
24. def __init__(self, dphi_pname, wmask_pname, method='weighted_mean', n_pca=1):
25. """
27. :param dphi_pname: str, Pathlike
28. Path to the temporally unwrapped complex phase difference 3D x time Nifti image
29. :param wmask_pname: str, Pathlike
30. Path to weighting mask Nifti image. Contains floats [0, 1]
31. :param method: str, optional
32. Method for respiration waveform estimation ['weighted_mean' or 'pca']
33. :param n_pca: int, optional
34. Number of PCA components to retain [1]
35. """
37. # Load the 4D phase difference image
38. self.dphi_nii = nib.load(dphi_pname)
39. self.dphi = self.dphi_nii.get_fdata()
41. # Create a time vector from the TR and number of timepoints
42. self.tr_s = self.dphi_nii.header['pixdim'][4]
43. self.nt = self.dphi.shape[3]
45. # Time vector for mid-TR of each volume
46. self.t_s = (np.arange(0, self.nt) + 0.5) * self.tr_s
48. # Load the 3D probabilistic computation mask
49. self.wmask_nii = nib.load(wmask_pname)
50. self.w = self.wmask_nii.get_fdata()
52. # Clamp weights to [0, 1]
53. self.w = np.clip(self.w, 0, 1)
55. # PCA components to retain
56. self.n_pca = n_pca
58. assert method in ['weighted_mean', 'pca'], f'Unknown method {method}'
59. self.method = method
61. def estimate(self):
62. """
63. Respiratory waveform estimation from phase difference 3D x time image
64. """
66. print('\nRespiration Waveform Reconstruction')
67. print('-------------------------------------')
68. print(f' Method: {self.method}')
69. if self.method == 'pca':
70. print(f' Number of PCA components: {self.n_pca}')
71. print(f' TR: {self.tr_s:.3f} s')
72. print(f' Number of timepoints: {self.nt}')
73. print(f' Voxel dimensions: {self.dphi.shape[0:3]}')
75. # Flatten the spatial dimensions and transpose to create temporal-spatial data matrix X
76. # X is Nt x Nvox
77. X = self.dphi.reshape(-1, self.dphi.shape[3]).T
79. match self.method:
81. case 'weighted_mean':
82. # Weighted mean of signal
83. # Calculate the mask weighted spatial mean phase difference for all volumes
84. mean_dphi = np.average(X, axis=1, weights=self.w.flatten())
86. # High pass filter by subtracting savgol baseline
87. # For similarity to typical first PCA component
88. k = int(10.0 / self.tr_s)
89. bline = sp.signal.savgol_filter(mean_dphi, window_length=k, polyorder=3)
90. mean_dphi_hpf = mean_dphi - bline
92. respwave = np.zeros([self.nt, 2])
93. respwave_dict = {
94. 'Time': self.t_s,
95. 'Resp': mean_dphi_hpf,
96. }
97. self.resp_df = pd.DataFrame(respwave_dict)
99. # Add peak inhalation/exhalation labels
100. self.resp_peaks()
102. case 'pca':
103. # Weighted temporal PCA of signal
104. # Returns temporal components and spatial weights
105. temporal_components, spatial_weights, _ = self.weighted_pca(X, self.w.flatten(), self.n_pca)
107. # Retain first component as respiration waveform
108. respwave = np.zeros([self.nt, 1+self.n_pca])
109. respwave[:, 0] = self.t_s
110. respwave[:, 1:] = temporal_components
112. # Construct respiration waveform dataframe
113. resp_headers = [f'Resp{pc+1}' for pc in range(self.n_pca)]
114. self.resp_df = pd.DataFrame(respwave, columns=['Time'] + resp_headers)
116. # Add peak inhalation/exhalation labels
117. self.resp_peaks()
119. return self.resp_df
121. def resp_peaks(self):
122. """
123. Find peak inhalation and exhalation in the respiration waveform
124. """
126. # Find positive peaks in respiration waveform
127. min_period = 1.0 / self.tr_s
128. pos_peaks, _ = sp.signal.find_peaks(self.resp_df['Resp'], distance=min_period)
129. neg_peaks, _ = sp.signal.find_peaks(-self.resp_df['Resp'], distance=min_period)
131. # Add peak locations to respiration waveform dataframe
132. self.resp_df['Peak'] = 'n/a'
133. self.resp_df.loc[pos_peaks, 'Peak'] = 'Inhalation'
134. self.resp_df.loc[neg_peaks, 'Peak'] = 'Exhalation'
136. return self.resp_df
138. def to_bids(self, bids_outdir, bids_stub):
140. # Save the respiration waveform to BIDS physio format (tsv.gz + json)
141. save_bids_physio(bids_outdir, bids_stub, self.resp_df)
143. @staticmethod
144. def weighted_pca(X, w, n_pca=1):
145. """
146. Weighted PCA of temporal (Nt rows) x spatial (Nvox columns) data matrix X
147. - spatial weights are applied to each column of X
149. :param X: np.ndarray, shape (Nt, Nvox)
150. Temporal spatial data matrix
151. :param w: np.ndarray, shape (Nvox,)
152. Voxel Weights for each column of X
153. :param n_pca: int, optional (default = 1)
154. Number of PCA components to retain
156. """
158. # Step 1: Apply spatial weights
159. Xw = X * w # Broadcasting
161. # Step 2: Standardize the data (optional, depending on your needs)
162. mean = np.mean(Xw, axis=0)
163. sd = np.std(Xw, axis=0)
164. Xstd = np.zeros_like(Xw)
165. nonzero = sd > 0.0
166. Xstd[:, nonzero] = (Xw[:, nonzero] - mean[nonzero]) / sd[nonzero]
168. # Step 3: Perform PCA
169. pca = PCA(n_components=n_pca)
170. pca.fit(Xstd)
172. # Report key PCA results
173. print(f' PCA singular values : {pca.singular_values_}')
174. print(f' PCA exp var ratio : {pca.explained_variance_ratio_}')
176. # Step 4: Retrieve results
177. temporal_components = pca.transform(Xstd) # Principal components
178. spatial_patterns = pca.components_ # Eigenvectors for spatial patterns
180. return temporal_components, spatial_patterns, pca

phaseresp.py at commit 86a8854, under MIT · at the source