|
20 | 20 | from tqdm.auto import tqdm |
21 | 21 |
|
22 | 22 |
|
23 | | - |
24 | 23 | # %% |
25 | 24 | import logging |
| 25 | + |
26 | 26 | logging.getLogger("snake").setLevel(level=logging.DEBUG) |
27 | 27 |
|
28 | 28 | # %% |
|
39 | 39 |
|
40 | 40 | sim_conf.hardware.n_coils = 8 # Update to get multi coil results. |
41 | 41 | sim_conf.hardware.field_strength = 7 |
42 | | -sim_conf.fov.res_mm = (1,1,1) # 1 mm isotropic resolution |
43 | | -sim_conf.fov.size = (181,217,1) # Single slice (in mm) |
44 | | -sim_conf.fov.offset = (-90,-125, -10) # Located in the center of the cortex |
45 | | -sim_conf.fov.angles = (0,0,0) |
| 42 | +sim_conf.fov.res_mm = (1, 1, 1) # 1 mm isotropic resolution |
| 43 | +sim_conf.fov.size = (181, 217, 1) # Single slice (in mm) |
| 44 | +sim_conf.fov.offset = (-90, -125, -10) # Located in the center of the cortex |
| 45 | +sim_conf.fov.angles = (0, 0, 0) |
46 | 46 | # - |
47 | 47 |
|
48 | | -phantom = Phantom.from_brainweb(sub_id=4, sim_conf=sim_conf, tissue_file="tissue_7T", output_res=1) |
49 | | -phantom = phantom.resample(new_affine=sim_conf.fov.affine, new_shape=sim_conf.fov.shape, use_gpu=True) # already select the phantom here, to reduce computational burden. |
| 48 | +phantom = Phantom.from_brainweb( |
| 49 | + sub_id=4, sim_conf=sim_conf, tissue_file="tissue_7T", output_res=1 |
| 50 | +) |
| 51 | +phantom = phantom.resample( |
| 52 | + new_affine=sim_conf.fov.affine, new_shape=sim_conf.fov.shape, use_gpu=True |
| 53 | +) # already select the phantom here, to reduce computational burden. |
50 | 54 |
|
51 | 55 | contrast = phantom.contrast(sim_conf=sim_conf) |
52 | 56 |
|
|
55 | 59 |
|
56 | 60 | activation_handler = BlockActivationHandler( |
57 | 61 | block_off=20, |
58 | | - block_on=20, |
59 | | - duration=300, |
| 62 | + block_on=20, |
| 63 | + duration=300, |
60 | 64 | atlas="hardvard-oxford__cort-maxprob-thr0-1mm", |
61 | | - atlas_label = 48 # occipital cortex |
| 65 | + atlas_label=48, # occipital cortex |
62 | 66 | ) |
63 | 67 |
|
64 | 68 | example_phantom = activation_handler.get_static(phantom.copy(), sim_conf) |
65 | 69 |
|
66 | 70 | roi = example_phantom.masks[example_phantom.labels_idx["ROI"]] |
67 | 71 |
|
68 | | -roi_resampled = apply_affine(roi,new_affine=sim_conf.fov.affine, old_affine=phantom.affine,new_shape= sim_conf.fov.shape) |
| 72 | +roi_resampled = apply_affine( |
| 73 | + roi, |
| 74 | + new_affine=sim_conf.fov.affine, |
| 75 | + old_affine=phantom.affine, |
| 76 | + new_shape=sim_conf.fov.shape, |
| 77 | +) |
69 | 78 |
|
70 | 79 | plt.imshow(contrast.squeeze().T, origin="lower", cmap="gray") |
71 | 80 | plt.axis("off") |
72 | | -plt.contour(roi_resampled.T.squeeze(), levels=(0.01, 0.5),origin="lower", colors=["tab:blue", "tab:orange"]) |
| 81 | +plt.contour( |
| 82 | + roi_resampled.T.squeeze(), |
| 83 | + levels=(0.01, 0.5), |
| 84 | + origin="lower", |
| 85 | + colors=["tab:blue", "tab:orange"], |
| 86 | +) |
73 | 87 |
|
74 | 88 | # + |
75 | 89 | # %% |
|
93 | 107 |
|
94 | 108 | # %% |
95 | 109 | # Defining Acquisition model |
96 | | -#-------------------------- |
| 110 | +# -------------------------- |
97 | 111 |
|
98 | 112 |
|
99 | 113 | # Here we assume an almost infinite SNR and a T2* decay signal model |
100 | 114 | # in the modeling of the fMRI signal |
101 | 115 |
|
102 | 116 | engine_t2s = NufftAcquisitionEngine( |
103 | | - model="T2s", |
104 | | - snr=30000, |
105 | | - slice_2d=True, # Does not matter, we are in 2D anyway (: |
| 117 | + model="T2s", |
| 118 | + snr=30000, |
| 119 | + slice_2d=True, # Does not matter, we are in 2D anyway (: |
106 | 120 | ) |
107 | 121 |
|
108 | 122 | engine_t2s( |
109 | | - "example_2D+T.mrd", |
110 | | - sampler, |
111 | | - phantom, |
112 | | - sim_conf, |
113 | | - handlers=[activation_handler], |
114 | | - worker_chunk_size=20, |
115 | | - n_workers=1, |
116 | | - nufft_backend="finufft", |
| 123 | + "example_2D+T.mrd", |
| 124 | + sampler, |
| 125 | + phantom, |
| 126 | + sim_conf, |
| 127 | + handlers=[activation_handler], |
| 128 | + worker_chunk_size=20, |
| 129 | + n_workers=1, |
| 130 | + nufft_backend="finufft", |
117 | 131 | ) |
118 | 132 |
|
119 | 133 |
|
120 | | - |
121 | | - |
122 | 134 | # %% |
123 | 135 |
|
124 | | -from snake.toolkit.reconstructors import ZeroFilledReconstructor, SequentialReconstructor |
125 | | -with NonCartesianFrameDataLoader("example_2D+T.mrd") as data_loader: |
126 | | - rec_zer = ZeroFilledReconstructor(nufft_backend="finufft", density_compensation=None).reconstruct(data_loader) |
| 136 | +from snake.toolkit.reconstructors import ( |
| 137 | + ZeroFilledReconstructor, |
| 138 | + SequentialReconstructor, |
| 139 | +) |
| 140 | + |
| 141 | +with NonCartesianFrameDataLoader("example_2D+T.mrd", squeeze_dims=False) as data_loader: |
| 142 | + rec_zer = ZeroFilledReconstructor( |
| 143 | + nufft_backend="finufft", density_compensation=None |
| 144 | + ).reconstruct(data_loader) |
127 | 145 | dyn_datas = data_loader.get_all_dynamic() |
128 | 146 | phantom_dl = data_loader.get_phantom() |
129 | 147 | sim_conf = data_loader.get_sim_conf() |
130 | 148 |
|
131 | 149 |
|
132 | 150 | # %% |
133 | | -plt.imshow(abs(rec_zer[0, ...]), cmap="gray", origin='lower') |
| 151 | +plt.imshow(abs(rec_zer[0, ...]), cmap="gray", origin="lower") |
134 | 152 |
|
135 | 153 | # %% |
136 | 154 |
|
137 | 155 | from snake.toolkit.analysis.stats import contrast_zscore, get_scores |
138 | 156 |
|
139 | | -+ |
| 157 | + |
140 | 158 | waveform_name = "activation-block_on" |
141 | 159 | good_d = None |
142 | 160 | for d in dyn_datas: |
|
151 | 169 | del phantom |
152 | 170 | del dyn_datas |
153 | 171 |
|
154 | | -# - |
155 | | - |
156 | 172 | TR_vol = sim_conf.seq.TR |
157 | | -z_score = contrast_zscore( |
158 | | - rec_zer, TR_vol, bold_signal, bold_sample_time, "block-on" |
159 | | -) |
| 173 | +z_score = contrast_zscore(rec_zer, TR_vol, bold_signal, bold_sample_time, "block-on") |
160 | 174 | stats_results = get_scores(z_score, roi_mask, 0.5) |
161 | | -np.save(data_zscore_file, z_score) |
162 | | - |
163 | 175 |
|
164 | | -rec_zer.shape |
165 | | - |
166 | | -# %% |
167 | 176 |
|
168 | 177 | # %% |
0 commit comments