Skip to content

Commit

Permalink
Merge branch 'imr-framework:dev' into dev
Browse files Browse the repository at this point in the history
  • Loading branch information
@tonggehua
tonggehua authored Jun 23, 2022
2 parents 564e715 + 5ef964c commit 281155f
Show file tree
Hide file tree
Showing 77 changed files with 7,153 additions and 2,731 deletions.
43 changes: 31 additions & 12 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@

# PyPulseq: A Python Package for MRI Pulse Sequence Design

`Compatible with Pulseq 1.3.1`
`Compatible with Pulseq 1.4.0`

## Table of contents 🧾
1. [📚 Relevant literature][section-relevant-literature]
Expand All @@ -28,12 +28,12 @@ exported as a `.seq` file to be run on Siemens/[GE]/[Bruker] hardware by levera
Pulseq interpreters. This tool is targeted at MRI pulse sequence designers, researchers, students and other interested
users. It is a translation of the Pulseq framework originally written in Matlab [[3]][section-references].

👉 Currently, PyPulseq is compatible with Pulseq 1.3.1. 👈
👉 Currently, PyPulseq is compatible with Pulseq 1.4.0. 👈

It is strongly recommended to first read the [Pulseq specification] before proceeding. The specification
document defines the concepts required for pulse sequence design using PyPulseq.

If you use PyPulseq in your work, please include both citations:
If you use PyPulseq in your work, cite the following publications:
```
Ravi, Keerthi, Sairam Geethanath, and John Vaughan. "PyPulseq: A Python Package for MRI Pulse Sequence Design." Journal
of Open Source Software 4.42 (2019): 1725.
Expand All @@ -48,23 +48,42 @@ Design pulse sequences using `pypulseq` in your browser! Check out the [⚡ Ligh
learn how!

## 📚 [Relevant literature][scholar-citations] (reverse chronological)
1. Ravi, Keerthi Sravan, and Sairam Geethanath. "Autonomous Magnetic Resonance Imaging." medRxiv (2020).
2. Nunes, Rita G., et al. "Implementation of a Diffusion-Weighted Echo Planar Imaging sequence using the Open Source
1. Niso, G., Botvinik-Nezer, R., Appelhoff, S., De La Vega, A., Esteban, O., Etzel, J.A., Finc, K., Ganz, M., Gau, R.,
Halchenko, Y.O. and Herholz, P., 2022. Open and reproducible neuroimaging: from study inception to publication.
2. Tong, G., Gaspar, A.S., Qian, E., Ravi, K.S., Vaughan, J.T., Nunes, R.G. and Geethanath, S., 2022. Open-source
magnetic resonance imaging acquisition: Data and documentation for two validated pulse sequences. Data in Brief, 42,
p.108105.
3. Tong, G., Gaspar, A.S., Qian, E., Ravi, K.S., Vaughan Jr, J.T., Nunes, R.G. and Geethanath, S., 2022. A framework
for validating open-source pulse sequences. Magnetic resonance imaging, 87, pp.7-18.
4. Karakuzu, A., Appelhoff, S., Auer, T., Boudreau, M., Feingold, F., Khan, A.R., Lazari, A., Markiewicz, C., Mulder,
M.J., Phillips, C. and Salo, T., 2021. qMRI-BIDS: an extension to the brain imaging data structure for quantitative
magnetic resonance imaging data. medRxiv.
5. Karakuzu, A., Biswas, L., Cohen‐Adad, J. and Stikov, N., 2021. Vendor‐neutral sequences and fully transparent
workflows improve inter‐vendor reproducibility of quantitative MRI. Magnetic Resonance in Medicine.
6. Geethanath, S., Single echo reconstruction for rapid and silent MRI. (ISMRM) (2021).
7. Qian, E. and Geethanath, S., Open source Magnetic rEsonance fingerprinting pAckage (OMEGA). (ISMRM) (2021).
8. Ravi, K.S., O'Reilly, T., Vaughan Jr, J.T., Webb, A. and Geethanath, S., Seq2prospa: translating PyPulseq for
low-field imaging. (ISMRM) (2021).
9. Ravi, K.S., Vaughan Jr, J.T. and Geethanath, S., PyPulseq in a web browser: a zero footprint tool for collaborative
and vendor-neutral pulse sequence development. (ISMRM) (2021).
10. Ravi, K.S. and Geethanath, S., 2020. Autonomous magnetic resonance imaging. Magnetic Resonance Imaging, 73,
pp.177-185.
11. Nunes, Rita G., et al. "Implementation of a Diffusion-Weighted Echo Planar Imaging sequence using the Open Source
Hardware-Independent PyPulseq Tool." ISMRM & SMRT Virtual Conference & Exhibition, International Society for Magnetic
Resonance in Medicine (ISMRM) (2020).
3. Loktyushin, Alexander, et al. "MRzero--Fully automated invention of MRI sequences using supervised learning." arXiv
12. Loktyushin, Alexander, et al. "MRzero--Fully automated invention of MRI sequences using supervised learning." arXiv
preprint arXiv:2002.04265 (2020).
4. Jimeno, Marina Manso, et al. "Cross-vendor implementation of a Stack-of-spirals PRESTO BOLD fMRI sequence using
13. Jimeno, Marina Manso, et al. "Cross-vendor implementation of a Stack-of-spirals PRESTO BOLD fMRI sequence using
TOPPE and Pulseq." ISMRM & SMRT Virtual Conference & Exhibition, International Society for Magnetic Resonance in
Medicine (ISMRM) (2020).
5. Clarke, William T., et al. "Multi-site harmonization of 7 tesla MRI neuroimaging protocols." NeuroImage 206 (2020): 116335.
6. Geethanath, Sairam, and John Thomas Vaughan Jr. "Accessible magnetic resonance imaging: a review." Journal of
14. Clarke, William T., et al. "Multi-site harmonization of 7 tesla MRI neuroimaging protocols." NeuroImage 206 (2020): 116335.
15. Geethanath, Sairam, and John Thomas Vaughan Jr. "Accessible magnetic resonance imaging: a review." Journal of
Magnetic Resonance Imaging 49.7 (2019): e65-e77.
7. Tong, Gehua, et al. "Virtual Scanner: MRI on a Browser." Journal of Open Source Software 4.43 (2019): 1637.
8. Archipovas, Saulius, et al. "A prototype of a fully integrated environment for a collaborative work in MR sequence
16. Tong, Gehua, et al. "Virtual Scanner: MRI on a Browser." Journal of Open Source Software 4.43 (2019): 1637.
17. Archipovas, Saulius, et al. "A prototype of a fully integrated environment for a collaborative work in MR sequence
development for a reproducible research." ISMRM 27th Annual Meeting & Exhibition, International Society for
Magnetic Resonance in Medicine (ISMRM) (2019).
9. Pizetta, Daniel Cosmo. PyMR: a framework for programming magnetic resonance systems. Diss. Universidade de São
18. Pizetta, Daniel Cosmo. PyMR: a framework for programming magnetic resonance systems. Diss. Universidade de São
Paulo (2018).
---

Expand Down
2 changes: 1 addition & 1 deletion VERSION
View file
Original file line number Diff line number Diff line change
@@ -1 +1 @@
1.3.1post1
1.4.0
126 changes: 84 additions & 42 deletions pypulseq/SAR/SAR_calc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,17 +55,19 @@ def _load_Q() -> Tuple[np.ndarray, np.ndarray]:
Contains the Q-matrix of global SAR values for body-mass and head-mass respectively.
"""
# Load relevant Q matrices computed from the model - this code will be integrated later - starting from E fields
path_Q = str(Path(__file__).parent / 'QGlobal.mat')
path_Q = str(Path(__file__).parent / "QGlobal.mat")
Q = sio.loadmat(path_Q)
Q = Q['Q']
Q = Q["Q"]
val = Q[0, 0]

Qtmf = val['Qtmf']
Qhmf = val['Qhmf']
Qtmf = val["Qtmf"]
Qhmf = val["Qhmf"]
return Qtmf, Qhmf


def _SAR_from_seq(seq: Sequence, Qtmf: np.ndarray, Qhmf: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
def _SAR_from_seq(
seq: Sequence, Qtmf: np.ndarray, Qhmf: np.ndarray
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Compute global whole body and head only SAR values for the given `seq` object.
Expand All @@ -89,7 +91,7 @@ def _SAR_from_seq(seq: Sequence, Qtmf: np.ndarray, Qhmf: np.ndarray) -> Tuple[np
"""
# Identify RF blocks and compute SAR - 10 seconds must be less than twice and 6 minutes must be less than
# 4 (WB) and 3.2 (head-20)
block_events = seq.dict_block_events
block_events = seq.block_events
num_events = len(block_events)
t = np.zeros(num_events)
SAR_wbg = np.zeros(t.shape)
Expand All @@ -101,7 +103,7 @@ def _SAR_from_seq(seq: Sequence, Qtmf: np.ndarray, Qhmf: np.ndarray) -> Tuple[np
block_dur = calc_duration(block)
t[block_counter - 1] = t_prev + block_dur
t_prev = t[block_counter - 1]
if hasattr(block, 'rf'): # has rf
if hasattr(block, "rf"): # has rf
rf = block.rf
signal = rf.signal
# This rf could be parallel transmit as well
Expand Down Expand Up @@ -135,8 +137,18 @@ def _SAR_interp(SAR: np.ndarray, t: np.ndarray) -> Tuple[np.ndarray, np.ndarray]
return SAR_interp, t_sec


def _SAR_lims_check(SARwbg_lim_s, SARhg_lim_s, tsec) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray,
np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
def _SAR_lims_check(
SARwbg_lim_s, SARhg_lim_s, tsec
) -> Tuple[
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
np.ndarray,
]:
"""
Check for SAR violations as compared to IEC 10 second and 6 minute averages;
returns SAR values that are interpolated for the fixed IEC time intervals.
Expand Down Expand Up @@ -165,45 +177,63 @@ def _SAR_lims_check(SARwbg_lim_s, SARhg_lim_s, tsec) -> Tuple[np.ndarray, np.nda
six_min_threshold_hg = 3.2
ten_sec_threshold_hg = 6.4

SAR_wbg_lim_app = np.concatenate((np.zeros(5), SARwbg_lim_s, np.zeros(5)), axis=0)
SAR_wbg_lim_app = np.concatenate(
(np.zeros(5), SARwbg_lim_s, np.zeros(5)), axis=0
)
SAR_hg_lim_app = np.concatenate((np.zeros(5), SARhg_lim_s, np.zeros(5)), axis=0)

SAR_wbg_tensec = _do_sw_sar(SAR_wbg_lim_app, tsec, 10) # < 2 SARmax
SAR_hg_tensec = _do_sw_sar(SAR_hg_lim_app, tsec, 10) # < 2 SARmax
SAR_wbg_tensec_peak = np.round(np.max(SAR_wbg_tensec), 2)
SAR_hg_tensec_peak = np.round(np.max(SAR_hg_tensec), 2)

if (np.max(SAR_wbg_tensec) > ten_sec_threshold_wbg) or (np.max(SAR_hg_tensec) > ten_sec_threshold_hg):
print('Pulse exceeding 10 second Global SAR limits, increase TR')
SAR_wbg_sixmin = 'NA'
SAR_hg_sixmin = 'NA'
SAR_wbg_sixmin_peak = 'NA'
SAR_hg_sixmin_peak = 'NA'
if (np.max(SAR_wbg_tensec) > ten_sec_threshold_wbg) or (
np.max(SAR_hg_tensec) > ten_sec_threshold_hg
):
print("Pulse exceeding 10 second Global SAR limits, increase TR")
SAR_wbg_sixmin = "NA"
SAR_hg_sixmin = "NA"
SAR_wbg_sixmin_peak = "NA"
SAR_hg_sixmin_peak = "NA"

if tsec[-1] > 600:
SAR_wbg_lim_app = np.concatenate((np.zeros(300), SARwbg_lim_s, np.zeros(300)), axis=0)
SAR_hg_lim_app = np.concatenate((np.zeros(300), SARhg_lim_s, np.zeros(300)), axis=0)
SAR_wbg_lim_app = np.concatenate(
(np.zeros(300), SARwbg_lim_s, np.zeros(300)), axis=0
)
SAR_hg_lim_app = np.concatenate(
(np.zeros(300), SARhg_lim_s, np.zeros(300)), axis=0
)

SAR_hg_sixmin = _do_sw_sar(SAR_hg_lim_app, tsec, 600)
SAR_wbg_sixmin = _do_sw_sar(SAR_wbg_lim_app, tsec, 600)
SAR_wbg_sixmin_peak = np.round(np.max(SAR_wbg_sixmin), 2)
SAR_hg_sixmin_peak = np.round(np.max(SAR_hg_sixmin), 2)

if (np.max(SAR_hg_sixmin) > six_min_threshold_wbg) or (np.max(SAR_hg_sixmin) > six_min_threshold_hg):
print('Pulse exceeding 10 second Global SAR limits, increase TR')
if (np.max(SAR_hg_sixmin) > six_min_threshold_wbg) or (
np.max(SAR_hg_sixmin) > six_min_threshold_hg
):
print("Pulse exceeding 10 second Global SAR limits, increase TR")
else:
print('Need at least 10 seconds worth of sequence to calculate SAR')
SAR_wbg_tensec = 'NA'
SAR_wbg_sixmin = 'NA'
print("Need at least 10 seconds worth of sequence to calculate SAR")
SAR_wbg_tensec = "NA"
SAR_wbg_sixmin = "NA"
SAR_hg_tensec = "NA"
SAR_hg_sixmin = "NA"
SAR_wbg_sixmin_peak = 'NA'
SAR_hg_sixmin_peak = 'NA'
SAR_wbg_tensec_peak = 'NA'
SAR_hg_tensec_peak = 'NA'

return SAR_wbg_tensec, SAR_wbg_sixmin, SAR_hg_tensec, SAR_hg_sixmin, \
SAR_wbg_sixmin_peak, SAR_hg_sixmin_peak, SAR_wbg_tensec_peak, SAR_hg_tensec_peak
SAR_wbg_sixmin_peak = "NA"
SAR_hg_sixmin_peak = "NA"
SAR_wbg_tensec_peak = "NA"
SAR_hg_tensec_peak = "NA"

return (
SAR_wbg_tensec,
SAR_wbg_sixmin,
SAR_hg_tensec,
SAR_hg_sixmin,
SAR_wbg_sixmin_peak,
SAR_hg_sixmin_peak,
SAR_wbg_tensec_peak,
SAR_hg_tensec_peak,
)


def _do_sw_sar(SAR: np.ndarray, tsec: np.ndarray, t: np.ndarray) -> np.ndarray:
Expand All @@ -225,9 +255,13 @@ def _do_sw_sar(SAR: np.ndarray, tsec: np.ndarray, t: np.ndarray) -> np.ndarray:
Sliding window time average of SAR values.
"""
SAR_time_avg = np.zeros(len(tsec) + int(t))
for instant in range(int(t / 2), int(t / 2) + (int(tsec[-1]))): # better to go from -sw / 2: sw / 2
SAR_time_avg[instant] = sum(SAR[range(instant - int(t / 2), instant + int(t / 2) - 1)]) / t
SAR_time_avg = SAR_time_avg[int(t / 2):int(t / 2) + (int(tsec[-1]))]
for instant in range(
int(t / 2), int(t / 2) + (int(tsec[-1]))
): # better to go from -sw / 2: sw / 2
SAR_time_avg[instant] = (
sum(SAR[range(instant - int(t / 2), instant + int(t / 2) - 1)]) / t
)
SAR_time_avg = SAR_time_avg[int(t / 2) : int(t / 2) + (int(tsec[-1]))]
return SAR_time_avg


Expand Down Expand Up @@ -255,28 +289,36 @@ def calc_SAR(file: Union[str, Path, Sequence]) -> None:
seq_obj.read(str(file))
seq_obj = seq_obj
else:
raise ValueError('Seq file does not exist.')
raise ValueError("Seq file does not exist.")
else:
seq_obj = file

Q_tmf, Q_hmf = _load_Q()
SAR_wbg, SAR_hg, t = _SAR_from_seq(seq_obj, Q_tmf, Q_hmf)
SARwbg_lim, tsec = _SAR_interp(SAR_wbg, t)
SARhg_lim, tsec = _SAR_interp(SAR_hg, t)
SAR_wbg_tensec, SAR_wbg_sixmin, SAR_hg_tensec, SAR_hg_sixmin, SAR_wbg_sixmin_peak, SAR_hg_sixmin_peak, \
SAR_wbg_tensec_peak, SAR_hg_tensec_peak = _SAR_lims_check(SARwbg_lim, SARhg_lim, tsec)
(
SAR_wbg_tensec,
SAR_wbg_sixmin,
SAR_hg_tensec,
SAR_hg_sixmin,
SAR_wbg_sixmin_peak,
SAR_hg_sixmin_peak,
SAR_wbg_tensec_peak,
SAR_hg_tensec_peak,
) = _SAR_lims_check(SARwbg_lim, SARhg_lim, tsec)

# Plot 10 sec average SAR
if (tsec[-1] > 10):
plt.plot(tsec, SAR_wbg_tensec, 'x-', label='Whole Body: 10sec')
plt.plot(tsec, SAR_hg_tensec, '.-', label='Head only: 10sec')
if tsec[-1] > 10:
plt.plot(tsec, SAR_wbg_tensec, "x-", label="Whole Body: 10sec")
plt.plot(tsec, SAR_hg_tensec, ".-", label="Head only: 10sec")

# plt.plot(t, SARwbg, label='Whole Body - instant')
# plt.plot(t, SARhg, label='Whole Body - instant')

plt.xlabel('Time (s)')
plt.ylabel('SAR (W/kg)')
plt.title('Global SAR - Mass Normalized - Whole body and head only')
plt.xlabel("Time (s)")
plt.ylabel("SAR (W/kg)")
plt.title("Global SAR - Mass Normalized - Whole body and head only")

plt.legend()
plt.grid(True)
Expand Down
Loading

0 comments on commit 281155f

Please sign in to comment.