Skip to content

Commit d742956

Browse files
evgueni-ovtchinnikovevgueni-ovtchinnikov
authored
Array views for sirf.Gadgetron data containers. (#1332)
Array views for sirf.Gadgetron data containers implemented, delivering significant acceleration of data algebra. Also addressed numpy 1.* compatibility
1 parent 306f0e7 commit d742956

File tree

21 files changed

+1256
-38
lines changed

21 files changed

+1256
-38
lines changed

‎CHANGES.md‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
# ChangeLog
2-
## vx.x.x
2+
3+
## v3.9.0
34

45
* CI
56
- made tests return value handling compatible with a future version of pytest.
@@ -12,7 +13,10 @@
1213
- Performance of acquisitions and images data algebra improved, acquisitions algebra running up to 3 times faster and images algebra up to 15 times faster.
1314
- `DataContainer.supports_array_view` to test for zero-copy compatibility.
1415
- SIRF interfaces (C++ and Python) for STIR Poisson noise generation utilities provided.
15-
- `ImageData` and `AcquisitionData` have `.asarray(copy=None)` (NumPy-like behaviour: default zero-copy if contiguous, fallback to deepcopy otherwise) via `__array_interface__`.
16+
- PET/Registration `ImageData` and PET `AcquisitionData` have `.asarray(copy=False)` (NumPy-like behaviour: default zero-copy if contiguous, fallback to deepcopy otherwise) via `__array_interface__`.
17+
18+
* SIRF/Gadgetron
19+
- `ImageDataView` and `AcquisitionDataView` classes implemented that encapsulate arrays of NumPy views of `ISMRMRD_ImageData` and `ISMRMRD_AcquisitionData` objects respectively, significantly accelerating the algebraic operations (up to a factor of about 10 for images).
1620

1721
## v3.8.1
1822

‎NOTICE.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
11
CCP SyneRBI (formerly PETMR) Synergistic Image Reconstruction Framework (SIRF).
2-
Copyright 2015-2024 Rutherford Appleton Laboratory, STFC, UK Research and Innovation
2+
Copyright 2015-2025 Rutherford Appleton Laboratory, STFC, UK Research and Innovation
33
Copyright 2015-2024 University College London (UCL)
44
Copyright 2017-2023 Physikalisch-Technische Bundesanstalt (PTB)
55
Copyright 2019 University of Hull

‎examples/Python/MR/data_views.py‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
# -*- coding: utf-8 -*-
2+
"""sirf.Gadgetron Views Test.
3+
v{version}
4+
5+
DataContainer algebra tests
6+
7+
Usage:
8+
data_views.py
9+
10+
{author}
11+
12+
{licence}
13+
"""
14+
from sirf.SIRF import norm, dot, copyto
15+
from sirf.Gadgetron import AcquisitionData, \
16+
AcquisitionDataView, ImageDataView, FullySampledReconstructor
17+
from sirf.Utilities import examples_data_path, error
18+
__version__ = "3.1.0"
19+
__author__ = "Evgueni Ovtchinnikov, Casper da Costa-Luis"
20+
21+
22+
def main():
23+
data_path = examples_data_path('MR')
24+
AcquisitionData.set_storage_scheme('memory')
25+
input_data = AcquisitionData(data_path + '/simulated_MR_2D_cartesian.h5')
26+
27+
prep_gadgets = ['RemoveROOversamplingGadget']
28+
processed_data = input_data.process(prep_gadgets)
29+
x = processed_data
30+
y = x + 0
31+
z = x + 0
32+
33+
x_view = AcquisitionDataView(x)
34+
y_view = AcquisitionDataView(y)
35+
z_view = AcquisitionDataView(z)
36+
37+
y = 2*x
38+
norm_y = y.norm()
39+
copyto(y_view, x_view)
40+
y_view *= 2
41+
norm_y_view = norm(y_view)
42+
print(norm_y, norm_y_view)
43+
z = x + y
44+
norm_z = z.norm()
45+
copyto(z_view, x_view)
46+
z_view += y_view
47+
norm_z_view = norm(z_view)
48+
print(norm_z, norm_z_view)
49+
s = norm(x)
50+
t = norm(x_view)
51+
print(s, t)
52+
s = x.sum()
53+
t = x_view.sum()
54+
print(s, t)
55+
s = x.dot(y)
56+
t = dot(x_view, y_view)
57+
print(s, t)
58+
59+
recon = FullySampledReconstructor()
60+
recon.set_input(processed_data)
61+
recon.process()
62+
complex_images = recon.get_output()
63+
64+
x = complex_images
65+
y = x + 0
66+
z = x + 0
67+
68+
x_view = ImageDataView(x)
69+
y_view = ImageDataView(y)
70+
z_view = ImageDataView(z)
71+
72+
y = 2*x
73+
norm_y = y.norm()
74+
copyto(y_view, x_view)
75+
y_view *= 2
76+
norm_y_view = norm(y_view)
77+
print(norm_y, norm_y_view)
78+
z = x + y
79+
norm_z = z.norm()
80+
copyto(z_view, x_view)
81+
z_view += y_view
82+
norm_z_view = norm(z_view)
83+
print(norm_z, norm_z_view)
84+
s = norm(x)
85+
t = norm(x_view)
86+
print(s, t)
87+
s = x.sum()
88+
t = x_view.sum()
89+
print(s, t)
90+
s = x.dot(y)
91+
t = dot(x_view, y_view)
92+
print(s, t)
93+
94+
try:
95+
main()
96+
print('\n=== done with %s' % __file__)
97+
98+
except error as err:
99+
# display error information
100+
print('??? %s' % err.value)
101+
exit(1)
102+

‎examples/Python/MR/mr_timings.py‎

Lines changed: 253 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,253 @@
1+
# -*- coding: utf-8 -*-
2+
3+
"""
4+
sirf.Gadgetron Views Test.
5+
Compares the performance of two approaches that avoid copying data between C++ and Python:
6+
- C++ algebra that is optimised by using templated loops
7+
- Python algebra that directly accesses C++ data via NumPy Array API views
8+
9+
Usage:
10+
mr_timings [--help | options]
11+
12+
Options:
13+
-f <file>, --file=<file> raw data file
14+
[default: simulated_MR_2D_cartesian.h5]
15+
-p <path>, --path=<path> path to data files, defaults to data/examples/MR
16+
subfolder of SIRF root folder
17+
-t <tsts>, --tests=<tsts> number of tests [default: 1]
18+
--non-interactive do not show plots
19+
"""
20+
21+
import numpy
22+
import timeit
23+
import importlib
24+
25+
from sirf.SIRF import copyto
26+
from sirf.Utilities import examples_data_path, existing_filepath
27+
28+
__version__ = "0.1.0"
29+
from docopt import docopt
30+
args = docopt(__doc__, version=__version__)
31+
#print(args)
32+
33+
# process command-line options
34+
data_file = args['--file']
35+
data_path = args['--path']
36+
if data_path is None:
37+
data_path = examples_data_path('MR')
38+
input_file = existing_filepath(data_path, data_file)
39+
ntests = int(args['--tests'])
40+
41+
# import engine module
42+
mr = importlib.import_module('sirf.Gadgetron')
43+
44+
# obtain acquisitions and image data
45+
acq = mr.AcquisitionData(input_file)
46+
processed_data = mr.preprocess_acquisition_data(acq)
47+
recon = mr.FullySampledReconstructor()
48+
recon.set_input(processed_data)
49+
print('---\n reconstructing...')
50+
recon.process()
51+
img = recon.get_output()
52+
53+
tests = 5
54+
'''
55+
tests:
56+
x * 2
57+
x + 2
58+
x + y
59+
x * y
60+
x / y
61+
'''
62+
63+
print('\n=== Comparing acquisitions algebra norms and timings...\n')
64+
65+
ax = acq
66+
ay = ax + 0
67+
az = ay + 0
68+
vx = mr.AcquisitionDataView(ax)
69+
vy = mr.AcquisitionDataView(ay)
70+
vz = mr.AcquisitionDataView(az)
71+
72+
temp_t = numpy.zeros(tests)
73+
view_t = numpy.zeros(tests)
74+
75+
for test in range(ntests):
76+
77+
start = timeit.default_timer()
78+
ay = 2*ax
79+
elapsed = timeit.default_timer() - start
80+
temp_t[0] += elapsed
81+
norm_y = ay.norm()
82+
start = timeit.default_timer()
83+
#vy.copy(vx)
84+
copyto(vy, vx)
85+
vy *= 2
86+
elapsed = timeit.default_timer() - start
87+
view_t[0] += elapsed
88+
if test == 0:
89+
print(f'norm(x * 2): {norm_y} {vy.norm()}')
90+
91+
start = timeit.default_timer()
92+
ay = ax + 2
93+
elapsed = timeit.default_timer() - start
94+
temp_t[1] += elapsed
95+
norm_y = ay.norm()
96+
start = timeit.default_timer()
97+
#vy.copy(vx)
98+
copyto(vy, vx)
99+
vy += 2
100+
elapsed = timeit.default_timer() - start
101+
view_t[1] += elapsed
102+
if test == 0:
103+
print(f'norm(x + 2): {norm_y} {vy.norm()}')
104+
105+
start = timeit.default_timer()
106+
az = ax + ay
107+
elapsed = timeit.default_timer() - start
108+
temp_t[2] += elapsed
109+
norm_z = az.norm()
110+
start = timeit.default_timer()
111+
#vz.copy(vx)
112+
copyto(vz, vx)
113+
vz += vy
114+
elapsed = timeit.default_timer() - start
115+
view_t[2] += elapsed
116+
if test == 0:
117+
print(f'norm(x + y): {norm_z} {vz.norm()}')
118+
119+
start = timeit.default_timer()
120+
az = ax * ay
121+
elapsed = timeit.default_timer() - start
122+
temp_t[3] += elapsed
123+
norm_z = az.norm()
124+
start = timeit.default_timer()
125+
#vz.copy(vx)
126+
copyto(vz, vx)
127+
vz *= vy
128+
elapsed = timeit.default_timer() - start
129+
view_t[3] += elapsed
130+
if test == 0:
131+
print(f'norm(x * y): {norm_z} {vz.norm()}')
132+
133+
ay += 1e-10
134+
start = timeit.default_timer()
135+
az = ax / ay
136+
elapsed = timeit.default_timer() - start
137+
temp_t[4] += elapsed
138+
norm_z = az.norm()
139+
start = timeit.default_timer()
140+
#vz.copy(vx)
141+
copyto(vz, vx)
142+
vz /= vy
143+
elapsed = timeit.default_timer() - start
144+
view_t[4] += elapsed
145+
if test == 0:
146+
print(f'norm(x / y): {norm_z} {vz.norm()}')
147+
148+
view_t /= ntests
149+
temp_t /= ntests
150+
151+
print('\ntest templates views')
152+
print(f'x * 2 {temp_t[0]:.2e} {view_t[0]:.2e}')
153+
print(f'x + 2 {temp_t[1]:.2e} {view_t[1]:.2e}')
154+
print(f'x + y {temp_t[2]:.2e} {view_t[2]:.2e}')
155+
print(f'x * y {temp_t[3]:.2e} {view_t[3]:.2e}')
156+
print(f'x / y {temp_t[4]:.2e} {view_t[4]:.2e}')
157+
158+
print('\n=== Comparing images algebra norms and timings...\n')
159+
160+
ix = img
161+
iy = ix + 0
162+
iz = iy + 0
163+
vx = mr.ImageDataView(ix)
164+
vy = mr.ImageDataView(iy)
165+
vz = mr.ImageDataView(iz)
166+
167+
temp_t = numpy.zeros(tests)
168+
view_t = numpy.zeros(tests)
169+
170+
for test in range(ntests):
171+
172+
start = timeit.default_timer()
173+
iy = 2*ix
174+
elapsed = timeit.default_timer() - start
175+
temp_t[0] += elapsed
176+
norm_y = iy.norm()
177+
start = timeit.default_timer()
178+
#vy.copy(vx)
179+
copyto(vy, vx)
180+
vy *= 2
181+
elapsed = timeit.default_timer() - start
182+
view_t[0] += elapsed
183+
if test == 0:
184+
print(f'norm(x * 2): {norm_y} {vy.norm()}')
185+
186+
start = timeit.default_timer()
187+
iy = ix + 2
188+
elapsed = timeit.default_timer() - start
189+
temp_t[1] += elapsed
190+
norm_y = iy.norm()
191+
start = timeit.default_timer()
192+
#vy.copy(vx)
193+
copyto(vy, vx)
194+
vy += 2
195+
elapsed = timeit.default_timer() - start
196+
view_t[1] += elapsed
197+
if test == 0:
198+
print(f'norm(x + 2): {norm_y} {vy.norm()}')
199+
200+
start = timeit.default_timer()
201+
iz = ix + iy
202+
elapsed = timeit.default_timer() - start
203+
temp_t[2] += elapsed
204+
norm_z = iz.norm()
205+
start = timeit.default_timer()
206+
#vz.copy(vx)
207+
copyto(vz, vx)
208+
vz += vy
209+
elapsed = timeit.default_timer() - start
210+
view_t[2] += elapsed
211+
if test == 0:
212+
print(f'norm(x + y): {norm_z} {vz.norm()}')
213+
214+
start = timeit.default_timer()
215+
iz = ix * iy
216+
elapsed = timeit.default_timer() - start
217+
temp_t[3] += elapsed
218+
norm_z = iz.norm()
219+
start = timeit.default_timer()
220+
#vz.copy(vx)
221+
copyto(vz, vx)
222+
vz *= vy
223+
elapsed = timeit.default_timer() - start
224+
view_t[3] += elapsed
225+
if test == 0:
226+
print(f'norm(x * y): {norm_z} {vz.norm()}')
227+
228+
ay += 1e-10
229+
start = timeit.default_timer()
230+
iz = ix / iy
231+
elapsed = timeit.default_timer() - start
232+
temp_t[4] += elapsed
233+
norm_z = iz.norm()
234+
start = timeit.default_timer()
235+
#vz.copy(vx)
236+
copyto(vz, vx)
237+
vz /= vy
238+
elapsed = timeit.default_timer() - start
239+
view_t[4] += elapsed
240+
if test == 0:
241+
print(f'norm(x / y): {norm_z} {vz.norm()}')
242+
243+
view_t /= ntests
244+
temp_t /= ntests
245+
246+
print('\ntest templates views')
247+
print(f'x * 2 {temp_t[0]:.2e} {view_t[0]:.2e}')
248+
print(f'x + 2 {temp_t[1]:.2e} {view_t[1]:.2e}')
249+
print(f'x + y {temp_t[2]:.2e} {view_t[2]:.2e}')
250+
print(f'x * y {temp_t[3]:.2e} {view_t[3]:.2e}')
251+
print(f'x / y {temp_t[4]:.2e} {view_t[4]:.2e}')
252+
253+
print('\n=== done with %s' % __file__)

‎examples/Python/PET/asarray.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -113,3 +113,4 @@ def main():
113113
print('\n=== done with %s' % __file__)
114114
except error as err:
115115
print('%s' % err.value)
116+

0 commit comments

Comments
 (0)