-
Notifications
You must be signed in to change notification settings - Fork 17
/
perform_iterative_qPAI_reconstruction.py
253 lines (209 loc) · 11.1 KB
/
perform_iterative_qPAI_reconstruction.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
# SPDX-FileCopyrightText: 2021 Division of Intelligent Medical Systems, DKFZ
# SPDX-FileCopyrightText: 2021 Janek Groehl
# SPDX-License-Identifier: MIT
import os
import matplotlib.pyplot as plt
import numpy as np
from scipy.ndimage import zoom
import simpa as sp
from simpa import Tags
from simpa.utils.profiling import profile
from argparse import ArgumentParser
# FIXME temporary workaround for newest Intel architectures
os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE"
# TODO: Please make sure that a valid path_config.env file is located in your home directory, or that you
# point to the correct file in the PathManager().
def run_perform_iterative_qPAI_reconstruction(spacing: float | int = 0.2, path_manager=None,
visualise: bool = True):
"""
:param spacing: The simulation spacing between voxels
:param path_manager: the path manager to be used, typically sp.PathManager
:param visualise: If VISUALIZE is set to True, the reconstruction result will be plotted
:return: a run through of the example
"""
if path_manager is None:
path_manager = sp.PathManager()
VOLUME_TRANSDUCER_DIM_IN_MM = 30
VOLUME_PLANAR_DIM_IN_MM = 30
VOLUME_HEIGHT_IN_MM = 30
RANDOM_SEED = 471
VOLUME_NAME = "MyqPAIReconstruction_" + str(RANDOM_SEED)
# If VISUALIZE is set to True, the reconstruction result will be plotted
def create_example_tissue():
"""
This is a very simple example script of how to create a tissue definition.
It contains a muscular background, an epidermis layer on top of the muscles
and a blood vessel.
"""
background_dictionary = sp.Settings()
background_dictionary[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(0.05, 30, 0.9)
background_dictionary[Tags.STRUCTURE_TYPE] = Tags.BACKGROUND
epidermis_structure = sp.Settings()
epidermis_structure[Tags.PRIORITY] = 1
epidermis_structure[Tags.STRUCTURE_START_MM] = [0, 0, 2]
epidermis_structure[Tags.STRUCTURE_END_MM] = [0, 0, 2.5]
epidermis_structure[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(2.2, 100.0, 0.9)
epidermis_structure[Tags.CONSIDER_PARTIAL_VOLUME] = True
epidermis_structure[Tags.ADHERE_TO_DEFORMATION] = True
epidermis_structure[Tags.STRUCTURE_TYPE] = Tags.HORIZONTAL_LAYER_STRUCTURE
vessel_structure_1 = sp.Settings()
vessel_structure_1[Tags.PRIORITY] = 2
vessel_structure_1[Tags.STRUCTURE_START_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM / 2.5, 0,
VOLUME_HEIGHT_IN_MM / 2]
vessel_structure_1[Tags.STRUCTURE_END_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM / 2.5,
VOLUME_PLANAR_DIM_IN_MM, VOLUME_HEIGHT_IN_MM / 2]
vessel_structure_1[Tags.STRUCTURE_RADIUS_MM] = 1.75
vessel_structure_1[Tags.STRUCTURE_ECCENTRICITY] = 0.85
vessel_structure_1[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(5.2, 100.0, 0.9)
vessel_structure_1[Tags.CONSIDER_PARTIAL_VOLUME] = True
vessel_structure_1[Tags.ADHERE_TO_DEFORMATION] = True
vessel_structure_1[Tags.STRUCTURE_TYPE] = Tags.ELLIPTICAL_TUBULAR_STRUCTURE
vessel_structure_2 = sp.Settings()
vessel_structure_2[Tags.PRIORITY] = 3
vessel_structure_2[Tags.STRUCTURE_START_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM / 2, 0,
VOLUME_HEIGHT_IN_MM / 3]
vessel_structure_2[Tags.STRUCTURE_END_MM] = [VOLUME_TRANSDUCER_DIM_IN_MM / 2,
VOLUME_PLANAR_DIM_IN_MM, VOLUME_HEIGHT_IN_MM / 3]
vessel_structure_2[Tags.STRUCTURE_RADIUS_MM] = 0.75
vessel_structure_2[Tags.MOLECULE_COMPOSITION] = sp.TISSUE_LIBRARY.constant(3.0, 100.0, 0.9)
vessel_structure_2[Tags.CONSIDER_PARTIAL_VOLUME] = True
vessel_structure_2[Tags.STRUCTURE_TYPE] = Tags.CIRCULAR_TUBULAR_STRUCTURE
tissue_dict = sp.Settings()
tissue_dict[Tags.BACKGROUND] = background_dictionary
tissue_dict["epidermis"] = epidermis_structure
tissue_dict["vessel_1"] = vessel_structure_1
tissue_dict["vessel_2"] = vessel_structure_2
return tissue_dict
# set settings for volume creation, optical simulation and iterative qPAI method
np.random.seed(RANDOM_SEED)
general_settings = {
# These parameters set the general properties of the simulated volume
Tags.RANDOM_SEED: RANDOM_SEED,
Tags.VOLUME_NAME: VOLUME_NAME,
Tags.SIMULATION_PATH: path_manager.get_hdf5_file_save_path(),
Tags.SPACING_MM: spacing,
Tags.DIM_VOLUME_Z_MM: VOLUME_HEIGHT_IN_MM,
Tags.DIM_VOLUME_X_MM: VOLUME_TRANSDUCER_DIM_IN_MM,
Tags.DIM_VOLUME_Y_MM: VOLUME_PLANAR_DIM_IN_MM,
Tags.WAVELENGTHS: [700]
}
settings = sp.Settings(general_settings)
settings.set_volume_creation_settings({
# These parameters set the properties for the volume creation
Tags.SIMULATE_DEFORMED_LAYERS: True,
Tags.STRUCTURES: create_example_tissue()
})
settings.set_optical_settings({
# These parameters set the properties for the optical Monte Carlo simulation
Tags.OPTICAL_MODEL_NUMBER_PHOTONS: 1e7,
Tags.OPTICAL_MODEL_BINARY_PATH: path_manager.get_mcx_binary_path(),
Tags.OPTICAL_MODEL: Tags.OPTICAL_MODEL_MCX,
Tags.LASER_PULSE_ENERGY_IN_MILLIJOULE: 50
})
settings["noise_model"] = {
Tags.NOISE_MEAN: 1.0,
Tags.NOISE_STD: 0.01,
Tags.NOISE_MODE: Tags.NOISE_MODE_MULTIPLICATIVE,
Tags.DATA_FIELD: Tags.DATA_FIELD_INITIAL_PRESSURE,
Tags.NOISE_NON_NEGATIVITY_CONSTRAINT: True
}
settings["iterative_qpai_reconstruction"] = {
# These parameters set the properties of the iterative reconstruction
Tags.DOWNSCALE_FACTOR: 0.75,
Tags.ITERATIVE_RECONSTRUCTION_CONSTANT_REGULARIZATION: False,
# the following tag has no effect, since the regularization is chosen to be SNR dependent, not constant
Tags.ITERATIVE_RECONSTRUCTION_REGULARIZATION_SIGMA: 0.01,
Tags.ITERATIVE_RECONSTRUCTION_MAX_ITERATION_NUMBER: 20,
# for this example, we are not interested in all absorption updates
Tags.ITERATIVE_RECONSTRUCTION_SAVE_INTERMEDIATE_RESULTS: False,
Tags.ITERATIVE_RECONSTRUCTION_STOPPING_LEVEL: 1e-3
}
# run pipeline including iterative qPAI method
pipeline = [
sp.ModelBasedAdapter(settings),
sp.MCXAdapter(settings),
sp.GaussianNoise(settings, "noise_model"),
sp.IterativeqPAI(settings, "iterative_qpai_reconstruction")
]
class CustomDevice(sp.PhotoacousticDevice):
def __init__(self):
super(CustomDevice, self).__init__(device_position_mm=np.asarray([general_settings[Tags.DIM_VOLUME_X_MM] / 2,
general_settings[Tags.DIM_VOLUME_Y_MM] / 2,
0]))
self.add_illumination_geometry(sp.DiskIlluminationGeometry(beam_radius_mm=20))
device = CustomDevice()
device.update_settings_for_use_of_model_based_volume_creator(settings)
sp.simulate(pipeline, settings, device)
# visualize reconstruction results
if visualise:
# get simulation output
data_path = path_manager.get_hdf5_file_save_path() + "/" + VOLUME_NAME + ".hdf5"
settings = sp.load_data_field(data_path, Tags.SETTINGS)
wavelength = settings[Tags.WAVELENGTHS][0]
# get reconstruction result
absorption_reconstruction = sp.load_data_field(data_path, Tags.ITERATIVE_qPAI_RESULT, wavelength)
# get ground truth absorption coefficients
absorption_gt = sp.load_data_field(data_path, Tags.DATA_FIELD_ABSORPTION_PER_CM, wavelength)
# rescale ground truth to same dimension as reconstruction (necessary due to resampling in iterative algorithm)
scale = np.shape(absorption_reconstruction)[0] / np.shape(absorption_gt)[0] # same as Tags.DOWNSCALE_FACTOR
absorption_gt = zoom(absorption_gt, scale, order=1, mode="nearest")
# compute reconstruction error
difference = absorption_gt - absorption_reconstruction
median_error = np.median(difference)
q3, q1 = np.percentile(difference, [75, 25])
iqr = q3 - q1
# visualize results
x_pos = int(np.shape(absorption_gt)[0] / 2)
y_pos = int(np.shape(absorption_gt)[1] / 2)
if np.min(absorption_gt) > np.min(absorption_reconstruction):
cmin = np.min(absorption_reconstruction)
else:
cmin = np.min(absorption_gt)
if np.max(absorption_gt) > np.max(absorption_reconstruction):
cmax = np.max(absorption_gt)
else:
cmax = np.max(absorption_reconstruction)
results_x_z = [absorption_gt[:, y_pos, :], absorption_reconstruction[:, y_pos, :], difference[:, y_pos, :]]
results_y_z = [absorption_gt[x_pos, :, :], absorption_reconstruction[x_pos, :, :], difference[x_pos, :, :]]
label = ["Absorption coefficients: ${\mu_a}^{gt}$", "Reconstruction: ${\mu_a}^{reconstr.}$",
"Difference: ${\mu_a}^{gt} - {\mu_a}^{reconstr.}$"]
plt.figure(figsize=(20, 15))
plt.subplots_adjust(hspace=0.5)
plt.suptitle("Iterative qPAI Reconstruction \n median error = " + str(np.round(median_error, 4)) +
"\n IQR = " + str(np.round(iqr, 4)), fontsize=10)
for i, quantity in enumerate(results_x_z):
plt.subplot(2, len(results_x_z), i + 1)
if i == 0:
plt.ylabel("x-z", fontsize=10)
plt.title(label[i], fontsize=10)
plt.imshow(quantity.T)
plt.xticks(fontsize=6)
plt.yticks(fontsize=6)
plt.colorbar()
if i != 2:
plt.clim(cmin, cmax)
else:
plt.clim(np.min(difference), np.max(difference))
for i, quantity in enumerate(results_y_z):
plt.subplot(2, len(results_x_z), i + len(results_x_z) + 1)
if i == 0:
plt.ylabel("y-z", fontsize=10)
plt.title(label[i], fontsize=10)
plt.imshow(quantity.T)
plt.xticks(fontsize=6)
plt.yticks(fontsize=6)
plt.colorbar()
if i != 2:
plt.clim(cmin, cmax)
else:
plt.clim(np.min(difference), np.max(difference))
plt.show()
plt.close()
if __name__ == "__main__":
parser = ArgumentParser(description='Run the iterative qPAI reconstruction example')
parser.add_argument("--spacing", default=0.2, type=float, help='the voxel spacing in mm')
parser.add_argument("--path_manager", default=None, help='the path manager, None uses sp.PathManager')
parser.add_argument("--visualise", default=True, type=bool, help='whether to visualise the result')
config = parser.parse_args()
run_perform_iterative_qPAI_reconstruction(spacing=config.spacing, path_manager=config.path_manager,
visualise=config.visualise)