Test Failed
Pull Request — master (#865)
by Daniil
03:26
created

BaseTomophantomLoader._populate_nexus_file()   B

Complexity

Conditions 6

Size

Total Lines 38
Code Lines 23

Duplication

Lines 0
Ratio 0 %

Importance

Changes 0
Metric Value
cc 6
eloc 23
nop 2
dl 0
loc 38
rs 8.3946
c 0
b 0
f 0
1
# Copyright 2014 Diamond Light Source Ltd.
2
#
3
# Licensed under the Apache License, Version 2.0 (the "License");
4
# you may not use this file except in compliance with the License.
5
# You may obtain a copy of the License at
6
#
7
#     http://www.apache.org/licenses/LICENSE-2.0
8
#
9
# Unless required by applicable law or agreed to in writing, software
10
# distributed under the License is distributed on an "AS IS" BASIS,
11
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12
# See the License for the specific language governing permissions and
13
# limitations under the License.
14
15
"""
16
.. module:: base_tomophantom_loader
17
   :platform: Unix
18
   :synopsis: A loader that generates synthetic 3D projection full-field tomo data\
19
        as hdf5 dataset of any size.
20
21
.. moduleauthor:: Daniil Kazantsev <[email protected]>
22
"""
23
24
import os
25
import h5py
26
import logging
27
import numpy as np
28
from mpi4py import MPI
29
30
from savu.data.chunking import Chunking
31
from savu.plugins.utils import register_plugin
32
from savu.plugins.loaders.base_loader import BaseLoader
33
from savu.plugins.savers.utils.hdf5_utils import Hdf5Utils
34
35
from savu.data.plugin_list import PluginList
36
37
import tomophantom
38
from tomophantom import TomoP2D, TomoP3D
39
40
@register_plugin
41
class BaseTomophantomLoader(BaseLoader):
42
    def __init__(self, name='BaseTomophantomLoader'):
43
        super(BaseTomophantomLoader, self).__init__(name)
44
45
    def setup(self):
46
        exp = self.exp
47
        data_obj = exp.create_data_object('in_data', 'synth_proj_data')
48
49
50
        data_obj.set_axis_labels(*self.parameters['axis_labels'])
51
        self.__convert_patterns(data_obj,'synth_proj_data')
52
        self.__parameter_checks(data_obj)
53
54
        self.tomo_model = self.parameters['tomo_model']
55
        # setting angles for parallel beam geometry
56
        self.angles = np.linspace(0.0,180.0-(1e-14), self.parameters['proj_data_dims'][0], dtype='float32')
57
        path = os.path.dirname(tomophantom.__file__)
58
        self.path_library3D = os.path.join(path, "Phantom3DLibrary.dat")
59
60
61
        data_obj.backing_file = self.__get_backing_file(data_obj, 'synth_proj_data')
62
        data_obj.data = data_obj.backing_file['/']['entry1']['tomo_entry']['data']['data']
63
        #data_obj.data.dtype # Need to do something to .data to keep the file open!
64
65
        # create a phantom file
66
        data_obj2 = exp.create_data_object('in_data', 'phantom')
67
        data_obj2.set_axis_labels(*['voxel_x.voxel', 'voxel_y.voxel', 'voxel_z.voxel'])
68
        self.__convert_patterns(data_obj2, 'phantom')
69
        self.__parameter_checks(data_obj2)
70
71
        data_obj2.backing_file = self.__get_backing_file(data_obj2, 'phantom')
72
        data_obj2.data = data_obj2.backing_file['/']['phantom']['data']
73
        #data_obj2.data.dtype # Need to do something to .data to keep the file open!
74
        data_obj.set_shape(data_obj.data.shape)
75
        group_name = '1-TomoPhantomLoader-phantom'
76
77
        self.n_entries = data_obj.get_shape()[0]
78
        cor_val=0.5*(self.parameters['proj_data_dims'][2])
79
        self.cor=np.linspace(cor_val, cor_val, self.parameters['proj_data_dims'][1], dtype='float32')
80
        self._set_metadata(data_obj, self._get_n_entries())
81
82
83
        return data_obj, data_obj2
84
85
    def __get_backing_file(self, data_obj, file_name):
86
        fname = '%s/%s.h5' % \
87
            (self.exp.get('out_path'), file_name)
88
89
        if os.path.exists(fname):
90
            return h5py.File(fname, 'r')
91
92
        self.hdf5 = Hdf5Utils(self.exp)
93
94
        dims_temp = self.parameters['proj_data_dims'].copy()
95
        proj_data_dims = tuple(dims_temp)
0 ignored issues
show
Comprehensibility Best Practice introduced by
The variable tuple does not seem to be defined.
Loading history...
96
        if (file_name == 'phantom'):
97
            dims_temp[0]=dims_temp[1]
98
            dims_temp[2]=dims_temp[1]
99
            proj_data_dims = tuple(dims_temp)
100
101
        patterns = data_obj.get_data_patterns()
102
        p_name = list(patterns.keys())[0]
103
        p_dict = patterns[p_name]
104
        p_dict['max_frames_transfer'] = 1
105
        nnext = {p_name: p_dict}
106
107
        pattern_idx = {'current': nnext, 'next': nnext}
108
        chunking = Chunking(self.exp, pattern_idx)
109
        chunks = chunking._calculate_chunking(proj_data_dims, np.int16)
110
111
        h5file = self.hdf5._open_backing_h5(fname, 'w')
112
113
        if file_name == 'phantom':
114
            group = h5file.create_group('/phantom', track_order=None)
115
        else:
116
            group = h5file.create_group('/entry1/tomo_entry/data', track_order=None)
117
118
        data_obj.dtype = np.dtype('<f4')
119
        dset = self.hdf5.create_dataset_nofill(group, "data", proj_data_dims, data_obj.dtype, chunks=chunks)
120
121
        self.exp._barrier()
122
123
124
        slice_dirs = list(nnext.values())[0]['slice_dims']
125
        nDims = len(dset.shape)
126
        total_frames = np.prod([dset.shape[i] for i in slice_dirs])
127
        sub_size = \
128
            [1 if i in slice_dirs else dset.shape[i] for i in range(nDims)]
129
130
        # need an mpi barrier after creating the file before populating it
131
        idx = 0
132
        sl, total_frames = \
133
            self.__get_start_slice_list(slice_dirs, dset.shape, total_frames)
134
        # calculate the first slice
135
        for i in range(total_frames):
136
            if (file_name == 'synth_proj_data'):
137
                #generate projection data
138
                gen_data = TomoP3D.ModelSinoSub(self.tomo_model, proj_data_dims[1], proj_data_dims[2], proj_data_dims[1], (i, i+1), -self.angles, self.path_library3D)
139
            else:
140
                #generate phantom data
141
                gen_data = TomoP3D.ModelSub(self.tomo_model, proj_data_dims[1], (i, i+1), self.path_library3D)
142
            dset[tuple(sl)] = np.swapaxes(gen_data,0,1)
143
            if sl[slice_dirs[idx]].stop == dset.shape[slice_dirs[idx]]:
144
                idx += 1
145
                if idx == len(slice_dirs):
146
                    break
147
            tmp = sl[slice_dirs[idx]]
148
            sl[slice_dirs[idx]] = slice(tmp.start+1, tmp.stop+1)
149
150
        self.exp._barrier()
151
152
153
154
        try:
155
            #nxsfile = NXdata(h5file)
156
            #nxsfile.save(file_name + ".nxs")
157
158
            h5file.close()
159
        except IOError as exc:
160
            logging.debug('There was a problem trying to close the file in random_hdf5_loader')
161
162
        return self.hdf5._open_backing_h5(fname, 'r')
163
164 View Code Duplication
    def __get_start_slice_list(self, slice_dirs, shape, n_frames):
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
165
        n_processes = len(self.exp.get('processes'))
166
        rank = self.exp.get('process')
167
        frames = np.array_split(np.arange(n_frames), n_processes)[rank]
168
        f_range = list(range(0, frames[0])) if len(frames) else []
169
        sl = [slice(0, 1) if i in slice_dirs else slice(None)
170
              for i in range(len(shape))]
171
        idx = 0
172
        for i in f_range:
173
            if sl[slice_dirs[idx]] == shape[slice_dirs[idx]]-1:
174
                idx += 1
175
            tmp = sl[slice_dirs[idx]]
176
            sl[slice_dirs[idx]] = slice(tmp.start+1, tmp.stop+1)
177
178
        return sl, len(frames)
179
180
    def __convert_patterns(self, data_obj, object_type):
181
        if (object_type == 'synth_proj_data'):
182
            pattern_list = self.parameters['patterns']
183
        else:
184
            pattern_list = self.parameters['patterns_tomo']
185
        for p in pattern_list:
186
            p_split = p.split('.')
187
            name = p_split[0]
188
            dims = p_split[1:]
189
            core_dims = tuple([int(i[0]) for i in [d.split('c') for d in dims]
190
                              if len(i) == 2])
191
            slice_dims = tuple([int(i[0]) for i in [d.split('s') for d in dims]
192
                               if len(i) == 2])
193
            data_obj.add_pattern(
194
                    name, core_dims=core_dims, slice_dims=slice_dims)
195
196
197
198
    def _set_metadata(self, data_obj, n_entries):
199
        n_angles = len(self.angles)
200
        data_angles = n_entries
201
        if data_angles != n_angles:
202
            raise Exception("The number of angles %s does not match the data "
203
                            "dimension length %s", n_angles, data_angles)
204
        data_obj.meta_data.set(['rotation_angle'], self.angles)
205
        data_obj.meta_data.set(['centre_of_rotation'], self.cor)
206
207
    def __parameter_checks(self, data_obj):
208
        if not self.parameters['proj_data_dims']:
209
            raise Exception(
210
                    'Please specifiy the dimensions of the dataset to create.')
211
212
    def _get_n_entries(self):
213
        return self.n_entries
214
215
216
    def post_process(self, data_obj, data_obj2):
217
218
        filename = self.exp.meta_data.get('nxs_filename')
219
        fsplit = filename.split('/')
220
        plugin_number = len(self.exp.meta_data.plugin_list.plugin_list)
221
        if plugin_number == 1:
222
            fsplit[-1] = 'synthetic_data.nxs'
223
        else:
224
            fsplit[-1] = 'synthetic_data_processed.nxs'
225
        filename = '/'.join(fsplit)
226
        self.exp.meta_data.set('nxs_filename', filename)
227
228
        plugin_list = PluginList()
229
        #plugin_list._save_plugin_list(filename)
230
        self.exp._finalise_setup(plugin_list)
231
        self._link_nexus_file(data_obj2, 'phantom', plugin_list)
232
        self._link_nexus_file(data_obj, 'synth_proj_data', plugin_list)
233
234
235
236
    def _link_nexus_file(self, data_obj, name, plugin_list):
237
        """Link phantom + synthetic projection data h5 files to a single nexus file containing both."""
238
239
240
        if name == 'phantom':
241
            data_obj.exp.meta_data.set(['group_name', 'phantom'], 'phantom')
242
            data_obj.exp.meta_data.set(['link_type', 'phantom'], 'final_result')
243
            data_obj.meta_data.set(["meta_data", "PLACEHOLDER", "VOLUME_XZ"], [10])
244
245
        else:
246
            data_obj.exp.meta_data.set(['group_name', 'synth_proj_data'], 'entry1/tomo_entry/data')
247
            data_obj.exp.meta_data.set(['link_type', 'synth_proj_data'], 'entry1')
248
249
        self._populate_nexus_file(data_obj)
250
        self._link_datafile_to_nexus_file(data_obj)
251
252
253
    def _populate_nexus_file(self, data):
254
        """"""
255
256
        filename = self.exp.meta_data.get('nxs_filename')
257
        name = data.data_info.get('name')
258
        #driver = "mpio", comm = MPI.COMM_WORLD
259
        with h5py.File(filename, 'a', driver="mpio", comm = MPI.COMM_WORLD) as nxs_file:
260
        #nxs_file = self.hdf5._open_backing_h5(filename, 'a', mpi=False)
261
262
            group_name = self.exp.meta_data.get(['group_name', name])
263
            link_type = self.exp.meta_data.get(['link_type', name])
264
265
            if name == 'phantom':
266
                if 'entry' not in list(nxs_file.keys()):
267
                    nxs_entry = nxs_file.create_group('entry')
268
                else:
269
                    nxs_entry = nxs_file['entry']
270
                if link_type == 'final_result':
271
                    group_name = 'final_result_' + data.get_name()
272
                else:
273
                    link = nxs_entry.require_group(link_type.encode("ascii"))
274
                    link.attrs['NX_class'] = 'NXcollection'
275
                    nxs_entry = link
276
277
                # delete the group if it already exists
278
                if group_name in nxs_entry:
279
                    del nxs_entry[group_name]
280
281
                plugin_entry = nxs_entry.require_group(group_name)
282
283
            else:
284
                plugin_entry = nxs_file.create_group(f'/{group_name}')
285
286
            self.__output_data_patterns(data, plugin_entry)
287
            self._output_metadata_dict(plugin_entry, data.meta_data.get_dictionary())
288
            self.__output_axis_labels(data, plugin_entry)
289
290
            plugin_entry.attrs['NX_class'] = 'NXdata'
291
292
293
    def __output_axis_labels(self, data, entry):
294
        axis_labels = data.data_info.get("axis_labels")
295
        ddict = data.meta_data.get_dictionary()
296
297
        axes = []
298
        count = 0
299
        dims_temp = self.parameters['proj_data_dims'].copy()
300
        if data.data_info.get('name') == 'phantom':
301
            dims_temp[0] = dims_temp[1]
302
            dims_temp[2] = dims_temp[1]
303
        dims = tuple(dims_temp)
304
305
        for labels in axis_labels:
306
            name = list(labels.keys())[0]
307
            axes.append(name)
308
            entry.attrs[name + '_indices'] = count
309
310
            mData = ddict[name] if name in list(ddict.keys()) \
311
                else np.arange(dims[count])
312
313
            if isinstance(mData, list):
314
                mData = np.array(mData)
315
316
            if 'U' in str(mData.dtype):
317
                mData = mData.astype(np.string_)
318
            if name not in list(entry.keys()):
319
                axis_entry = entry.require_dataset(name, mData.shape, mData.dtype)
320
                axis_entry[...] = mData[...]
321
                axis_entry.attrs['units'] = list(labels.values())[0]
322
            count += 1
323
        entry.attrs['axes'] = axes
324
325 View Code Duplication
    def __output_data_patterns(self, data, entry):
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
326
        data_patterns = data.data_info.get("data_patterns")
327
        entry = entry.require_group('patterns')
328
        entry.attrs['NX_class'] = 'NXcollection'
329
        for pattern in data_patterns:
330
            nx_data = entry.require_group(pattern)
331
            nx_data.attrs['NX_class'] = 'NXparameters'
332
            values = data_patterns[pattern]
333
            self.__output_data(nx_data, values['core_dims'], 'core_dims')
334
            self.__output_data(nx_data, values['slice_dims'], 'slice_dims')
335
336
    def _output_metadata_dict(self, entry, mData):
337
        entry.attrs['NX_class'] = 'NXcollection'
338
        for key, value in mData.items():
339
            if key != 'rotation_angle':
340
                nx_data = entry.require_group(key)
341
                if isinstance(value, dict):
342
                    self._output_metadata_dict(nx_data, value)
343
                else:
344
                    nx_data.attrs['NX_class'] = 'NXdata'
345
                    self.__output_data(nx_data, value, key)
346
347 View Code Duplication
    def __output_data(self, entry, data, name):
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
348
        if isinstance(data, dict):
349
            entry = entry.require_group(name)
350
            entry.attrs['NX_class'] = 'NXcollection'
351
            for key, value in data.items():
352
                self.__output_data(entry, value, key)
353
        else:
354
            try:
355
                self.__create_dataset(entry, name, data)
356
            except Exception:
357
                try:
358
                    import json
359
                    data = np.array([json.dumps(data).encode("ascii")])
360
                    self.__create_dataset(entry, name, data)
361
                except Exception:
362
                    try:
363
                        self.__create_dataset(entry, name, data)
364
                    except:
365
                        raise Exception('Unable to output %s to file.' % name)
366
367
    def __create_dataset(self, entry, name, data):
368
        if name not in list(entry.keys()):
369
            entry.create_dataset(name, data=data)
370
        else:
371
            entry[name][...] = data
372
373 View Code Duplication
    def _link_datafile_to_nexus_file(self, data):
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
374
        filename = self.exp.meta_data.get('nxs_filename')
375
376
        with h5py.File(filename, 'a', driver="mpio", comm = MPI.COMM_WORLD) as nxs_file:
377
        #nxs_file = self.hdf5._open_backing_h5(filename, 'a', mpi=False)
378
            # entry path in nexus file
379
            name = data.get_name()
380
            group_name = self.exp.meta_data.get(['group_name', name])
381
            link = self.exp.meta_data.get(['link_type', name])
382
            name = data.get_name(orig=True)
383
            nxs_entry = self.__add_nxs_entry(nxs_file, link, group_name, name)
384
            self.__add_nxs_data(nxs_file, nxs_entry, link, group_name, data)
385
386
    def __add_nxs_entry(self, nxs_file, link, group_name, name):
387
        if name == 'phantom':
388
            nxs_entry = '/entry/' + link
389
        else:
390
            nxs_entry = ''
391
        nxs_entry += '_' + name if link == 'final_result' else "/" + group_name
392
        nxs_entry = nxs_file[nxs_entry]
393
        nxs_entry.attrs['signal'] = 'data'
394
        return nxs_entry
395
396 View Code Duplication
    def __add_nxs_data(self, nxs_file, nxs_entry, link, group_name, data):
0 ignored issues
show
Duplication introduced by
This code seems to be duplicated in your project.
Loading history...
397
        data_entry = nxs_entry.name + '/data'
398
        # output file path
399
        h5file = data.backing_file.filename
400
401
        if link == 'input_data':
402
            dataset = self.__is_h5dataset(data)
403
            if dataset:
404
                nxs_file[data_entry] = \
405
                    h5py.ExternalLink(os.path.abspath(h5file), dataset.name)
406
        else:
407
            # entry path in output file path
408
            m_data = self.exp.meta_data.get
409
            if not (link == 'intermediate' and
410
                    m_data('inter_path') != m_data('out_path')):
411
                h5file = h5file.split(m_data('out_folder') + '/')[-1]
412
            nxs_file[data_entry] = \
413
                h5py.ExternalLink(h5file, group_name + '/data')
414