First pass at the data loading + simple tests, still need to do data saving tools

This commit is contained in:
Abe Levitan
2019-03-24 18:00:24 -04:00
parent f08a77d6da
commit b3b953fc8a
3 changed files with 739 additions and 1 deletions
+344
View File
@@ -0,0 +1,344 @@
from __future__ import division, print_function, absolute_import
import h5py
import numpy as np
__all__ = ['get_entry_info',
'get_sample_info',
'get_wavelength',
'get_detector_geometry',
'get_mask',
'get_data',
'get_ptycho_translations']
#
#
# I will put here some thoughts about how to load data into this program.
#
#
# The reconstructions should have the ability to generate datasets.
# So you could write a reconstruction engine and then it would be
# able to simulate data directly in the engine for you to use as a
# reconstruction
#
# I don't even think there needs to be a loading tool for loading cxi files
# because there isn't really a better method beyond just loading the
# file into an h5py object. This file could host the simple cxi file
# browser, perhaps. But I think the reality is that we need individual
# loaders for each kind of experiment. Perhaps we could put some basic
# reuseable tools for inspecting cxi-type h5 files in this file.
#
#
# Then, there can be some more sophisticated tools that load data for
# specific use cases that are common - loading data for a 2D CDI experiment,
# loading data for a 2D Ptycho experiment, loading data for Bragg Ptycho in
# 3D, loading data for a 3D CDI experiment, etc.
#
#
# Perhaps one good way to package this is for the kind of data associated
# with any particular experiment to have it's own kind of dataset or view.
# So there would be a "2D Ptychography" data viewer, which would contain
# all the measured data that comes from a 2D ptychography experiment.
# The specialized functions would plop out these data viewers, and the
# reconstruction classes could be designed around a particular kind of
# viewer with the most general kind just requiring a generic data viewer.
#
# Data viewers could have simple tools like the ability to send themselves
# to the GPU, CPU, change the datatype, etc. I think the most generic thing
# is as a subclass of the torch Data objects, where they would for each slice
# return the index, a set of defining parameters (translation, angle, energy,
# whatever), and a diffraction pattern. They would also have a "setup"
# attribute, or "metadata", or whatever you'd want to call it, that contain
# the various fixed experimental parameters (energy, distance, etc.)
#
# And I think the cxi visualizer should really go into it's own script,
# because it's not a reuseable component.
#
#
# Functions to inspect the basic attributes of a cxi file represented as an
# h5 file object
#
def get_entry_info(cxi_file):
"""Returns a dictionary with the basic metadata from the cxi file's entry_1 attribute
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
dict : A dictionary with basic metadata defined in the cxi file
"""
e1 = cxi_file['entry_1']
metadata_attrs = ['title',
'experiment_identifier',
'experiment_description',
'program_name',
'start_time',
'end_time']
metadata = {attr: str(e1[attr][()].decode()) for attr in metadata_attrs
if attr in e1}
return metadata
def get_sample_info(cxi_file):
"""Returns a dictionary with the basic metadata from the cxi file's entry_1/sample_1 attribute
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
dict : A dictionary with basic metadata from the sample defined in the cxi file
"""
if 'entry_1/sample_1' not in cxi_file:
return None
s1 = cxi_file['entry_1/sample_1']
metadata_attrs = ['name','description','unit_cell_group']
metadata = {attr: str(s1[attr][()].decode()) for attr in metadata_attrs
if attr in s1}
float_attrs = ['concentration',
'mass',
'temperature',
'thickness',
'unit_cell_volume']
for attr in float_attrs:
if attr in s1:
metadata[attr] = np.float32(s1[attr][()])
if 'unit_cell' in s1:
metadata['unit_cell'] = np.array(s1['unit_cell']).astype(np.float32)
# TODO: Add my nonstandard "surface normal" attribute here
# TODO: I should add the sample geometry as a valid metadata that can
# be copied over
return metadata
def get_wavelength(cxi_file):
"""Returns the wavelength of the source defined in the cxi file object, in m
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
np.float32 : The wavelength of the source defined in the cxi file
"""
i1 = cxi_file['entry_1/instrument_1']
if 'source_1/wavelength' in i1:
wavelength = np.float32(i1['source_1/wavelength'])
elif 'source_1/energy' in i1:
energy = np.float32(i1['source_1/energy'])
wavelength = 1.9864459e-25 / energy
else:
raise KeyError('Neither Wavelength or Energy Defined in provided .cxi File')
return wavelength
def get_detector_geometry(cxi_file):
"""Returns a standardized description of the detector geometry defined in the cxi file object
It makes intelligent assumptions based on the definitions in the cxi
file definition. The standardized description of the geometry that it
outputs includes the sample to detector distance, the corner location
of the detector, and the basis vectors defining the detector. It can
only handle detectors defined as rectangular grids of pixels.
The distance and corner_location values are technically overdetermining
the detector location, but for many experiments (particularly
transmission experiments), the distance is needed and the exact
corner location is not. If the corner location is not reported in
the cxi file, no attempt will be made to calculate it.
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
distance (np.float32) : The sample to detector distance, in m
basis_vectors (np.array) : The basis vectors for the detector
corner_location (np.array) : The location of the (0,0) pixel in the detector
"""
i1 = cxi_file['entry_1/instrument_1']
d1 = i1['detector_1']
if 'detector_1/basis_vectors' in i1:
basis_vectors = np.array(d1['basis_vectors'])
else:
# This whole thing just to account for all the ways people can
# implicitly define the x or y pixel size for a detector. I've
# seen too many of these in the wild, unfortunately...
try:
x_pixel_size = np.float32(d1['x_pixel_size'])
except:
x_pixel_size = None
try:
y_pixel_size = np.float32(d1['y_pixel_size'])
except:
y_pixel_size = None
if x_pixel_size is None and y_pixel_size is not None:
x_pixel_size = y_pixel_size
elif x_pixel_size is not None and y_pixel_size is None:
y_pixel_size = x_pixel_size
if x_pixel_size is None and y_pixel_size is None:
raise KeyError('Detector pixel size not defined in file.')
basis_vectors = np.array([[0,-y_pixel_size,0],
[-x_pixel_size,0,0]]).transpose()
try:
distance = np.float32(d1['distance'])
except:
distance = None
try:
corner_position = np.array(d1['corner_position'])
except:
corner_position = None
# Don't pretend to calculate corner position from distance if it's
# if it's not defined, but do calculate distance from corner position
# if distance is not defined. If neither is defined, then raise
# an error.
if distance is None and corner_position is not None:
detector_normal = np.cross(basis_vectors[:,0],
basis_vectors[:,1])
detector_normal /= np.linalg.norm(detector_normal)
distance = np.linalg.norm(np.dot(corner_position, detector_normal))
if distance is None and corner_position is not None:
raise KeyError('Neither sample to detector distance or corner position is defined in file.')
return distance, basis_vectors, corner_position
def get_mask(cxi_file):
"""Returns the detector mask defined in the cxi file object
This function converts from the format specified in the cxi file
definition to a simple on/off mask, where a value of 1 defines a
good pixel (on) and a value of 0 defines a bad pixel (off).
If any bit is set in the mask at all, it will be defined as a bad
pixel, with the exception of pixels marked exactly as 0x00001000,
which is defined to mean that the pixel has signal above the
background. These pixels are treated as on pixels
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
np.array : An array storing the mask from the cxi file
"""
i1 = cxi_file['entry_1/instrument_1']
if 'detector_1/mask' in i1:
mask = np.array(i1['detector_1/mask']).astype(np.uint32)
mask_on = np.equal(mask,np.uint32(0))
mask_has_signal = np.equal(mask,np.uint32(0x00001000))
return np.logical_or(mask_on,mask_has_signal).astype(np.uint8)
else:
return None
def get_data(cxi_file):
"""Returns an array with the full stack of detector data defined in the cxi file object
This function will make sure to check all the various places that it's
okay to store the data in, to ensure that it can find the data regardless
of whether the creator of the .cxi file has remembered to link the data
to all the required locations.
It will return the data array in whatever shape it's defined in.
It will also read out the axes attribute of the data into a list
of strings
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
np.array : An array storing the data defined in the cxi file
list : A list of the axes defined in the axes attribute, if any
"""
# Possible locations for the data
#
# entry_1/detector_1/data
if 'entry_1/data_1/data' in cxi_file:
pull_from = 'entry_1/data_1/data'
elif 'entry_1/instrument_1/detector_1/data' in cxi_file:
pull_from = 'entry_1/instrument_1/detector_1/data'
else:
raise KeyError('Data is not defined within cxi file')
data = np.array(cxi_file[pull_from]).astype(np.float32)
if 'axes' in cxi_file[pull_from].attrs:
axes = str(cxi_file[pull_from].attrs['axes'].decode()).split(':')
axes = [axis.strip().lower() for axis in axes]
else:
axes = None
return data, axes
def get_ptycho_translations(cxi_file):
"""Gets an array of x,y,z translations, if such an array has been defined in the file
It applies two operations to the translations. First, it negates them,
because the CXI file format is designed to specify translations of the
samples and the CDTools code specifies translations of the optics.
Second, it transposes the array so that the first axis is translation
ID and the second axis is the (x,y,z) components of the translation
Args:
cxi_file (h5py.File) : a file object to be read
Returns:
np.array : An array storing the translations defined in the cxi file
list : A list of the axes defined in the axes attribute, if any
"""
if 'entry_1/data_1/translation' in cxi_file:
pull_from = 'entry_1/data_1/translation'
elif 'entry_1/sample_1/geometry_1/translation' in cxi_file:
pull_from = 'entry_1/sample_1/geometry_1/translation'
elif 'entry_1/instrument_1/detector_1/translation' in cxi_file:
pull_from = 'entry_1/instrument_1/detector_1/translation'
else:
raise KeyError('Translations are not defined within cxi file')
translations = -np.array(cxi_file[pull_from]).astype(np.float32).transpose()
return translations
#
# It might be useful to make some helper functions to help write cxi files
#
#
# A function to place the skeleton of a cxi file down
#
#
# A function to define the source attributes
#
#
# A function to define the detector geometry
#
#
# A function to save out a mask, converting it to the correct format
#
#
# Perhaps a function to store the data and link it correctly? But this might
# have to change too much situation to situation
#
-1
View File
@@ -1 +0,0 @@
abe@paulsimon.2556:1553086360
+395
View File
@@ -0,0 +1,395 @@
from __future__ import division, print_function, absolute_import
from CDTools.tools import data
import numpy as np
import torch as t
import h5py
import pytest
import os
import datetime
from pathlib import Path
#
# First, we write a few fixtures to generate specific data files we
# want with deliberate pathologies.
#
# * Energy defined but no wavelength
# * Wavelength defined but no energy
# * Distance and pixel pitch defined but no corner position or basis
# * Corner position and basis defined but no distance or pixel pitch
# * No mask defined
# * Mask defined
#
# Each file will get a fixture that loads the cxi file but also loads
# a dictionary with the relevant information for the cxi file
#
#
# This just grabs the directory whose name matches the file we're running
@pytest.fixture(scope='module')
def datadir(request):
filename = request.module.__file__
test_dir, _ = os.path.splitext(filename)
return Path(test_dir)
@pytest.fixture(scope='module')
def ptycho_cxi_1():
"""Creates an example file for CXI ptychography. This file is defined
to have everything done as correctly as possible with lots of attributes
defined. It will return both a dictionary describing what is expected
to be loaded and a file with the data stored in it.
"""
expected = {}
f = h5py.File('ptycho_cxi_1',driver='core',backing_store=False)
# Start by defining the basic structure
f.create_dataset('cxi_version', data=150)
f.create_dataset('number_of_entries',data=1)
# Then define a bunch of metadata for entry_1
e1f = f.create_group('entry_1')
expected['entry metadata'] = {}
e1e = expected['entry metadata']
e1e['start_time'] = datetime.datetime.now().isoformat()
e1f['start_time'] = np.string_(e1e['start_time'])
e1e['end_time'] = datetime.datetime.now().isoformat()
e1f['end_time'] = np.string_(e1e['end_time'])
e1e['experiment_identifier'] = 'Fake Experiment 1'
e1f['experiment_identifier'] = np.string_(e1e['experiment_identifier'])
e1e['experiment_description'] = 'A fully defined ptychography experiment to test the data loading'
e1f['experiment_description'] = np.string_(e1e['experiment_description'])
e1e['program_name'] = 'CDTools'
e1f['program_name'] = np.string_(e1e['program_name'])
e1e['title'] = 'The one experiment we did'
e1f['title'] = np.string_(e1e['title'])
# Set up the sample info
s1f = e1f.create_group('sample_1')
expected['sample info'] = {}
s1e = expected['sample info']
s1e['name'] = 'Fake Sample'
s1f['name'] = np.string_(s1e['name'])
s1e['description'] = 'A sample that isn\'t real'
s1f['description'] = np.string_(s1e['description'])
s1e['unit_cell_group'] = 'P1'
s1f['unit_cell_group'] = np.string_(s1e['unit_cell_group'])
s1e['concentration'] = np.float32(np.random.rand())
s1f['concentration'] = s1e['concentration']
s1e['mass'] = np.float32(np.random.rand())
s1f['mass'] = s1e['mass']
s1e['temperature'] = np.float32(np.random.rand()*100)
s1f['temperature'] = s1e['temperature']
s1e['thickness'] = np.float32(np.random.rand()*1e-7)
s1f['thickness'] = s1e['thickness']
s1e['unit_cell_volume'] = np.float32(np.random.rand() * 1e-27)
s1f['unit_cell_volume'] = s1e['unit_cell_volume']
s1e['unit_cell'] = np.array([1,1,1,90,90,90]).astype(np.float32)
s1f.create_dataset('unit_cell',data = s1e['unit_cell'])
i1f = e1f.create_group('instrument_1')
source1f = i1f.create_group('source_1')
energy = np.float32(1.3618e-16) #Joules, = 850 eV
source1f['energy'] = energy
expected['wavelength'] = np.float32(1.9864459e-25) / energy
source1f['wavelength'] = expected['wavelength']
d1f = i1f.create_group('detector_1')
expected['detector'] = {}
d1e = expected['detector']
d1e['distance'] = np.float32(0.3)
d1f['distance'] = d1e['distance']
d1e['basis'] = np.array([[0,-30e-6,0],
[-20e-6,0,0]]).astype(np.float32).transpose()
d1f.create_dataset('basis_vectors',data=d1e['basis'])
d1f['x_pixel_size'] = np.float32(20e-6)
d1f['y_pixel_size'] = np.float32(30e-6)
d1e['corner'] = np.array((2550e-6,3825e-6,0.3)).astype(np.float32)
d1f.create_dataset('corner_position', data=d1e['corner'])
# Remember the format for the CXI file differs from the format used
# internally
mask = np.zeros((100,256,256)).astype(np.uint32)
expected['mask'] = np.ones((100,256,256)).astype(np.uint8)
d1f.create_dataset('mask',data=mask)
data1f = e1f.create_group('data_1')
data = np.random.rand(100,256,256).astype(np.float32)
expected['data'] = data
d1f.create_dataset('data',data=data)
data1f['data'] = h5py.SoftLink('/entry_1/instrument_1/detector_1/data')
d1f['data'].attrs['axes'] = np.string_('translation:y:x')
expected['axes'] = ['translation','y','x']
g1f = s1f.create_group('geometry_1')
translations = np.arange(300).reshape((100,3)).astype(np.float32)
g1f.create_dataset('translation',data=translations)
data1f['translation'] = h5py.SoftLink('/entry_1/sample_1/geometry_1/translation')
d1f['translation'] = h5py.SoftLink('/entry_1/sample_1/geometry_1/translation')
expected['translations'] = -translations.transpose()
return f, expected
@pytest.fixture(scope='module')
def ptycho_cxi_2():
"""Creates an example file for CXI ptychography. This file is defined
to have a subset of things missing. In particular, it:
* Defines the wavelength but not the energy
* Defines the corner position but not the sample-detector distance
* Defines pixel sizes but no basis vectors
* Doesn't define a mask
* Only defines data in the relevant places, not under data_1
* Doesn't explicitly define axes for the data arrays
* Is missing many allowed metadata attributes
"""
expected = {}
f = h5py.File('ptycho_cxi_2',driver='core',backing_store=False)
# Start by defining the basic structure
f.create_dataset('cxi_version', data=150)
f.create_dataset('number_of_entries',data=1)
# Then define a bunch of metadata for entry_1
e1f = f.create_group('entry_1')
expected['entry metadata'] = {}
e1e = expected['entry metadata']
e1e['title'] = 'The one experiment we did'
e1f['title'] = np.string_(e1e['title'])
# Set up the sample info
s1f = e1f.create_group('sample_1')
expected['sample info'] = {}
s1e = expected['sample info']
s1e['temperature'] = np.float32(np.random.rand()*100)
s1f['temperature'] = s1e['temperature']
i1f = e1f.create_group('instrument_1')
source1f = i1f.create_group('source_1')
energy = np.float32(1.3618e-16) #Joules, = 850 eV
expected['wavelength'] = np.float32(1.9864459e-25) / energy
source1f['wavelength'] = expected['wavelength']
d1f = i1f.create_group('detector_1')
expected['detector'] = {}
d1e = expected['detector']
d1e['distance'] = np.float32(0.3)
d1e['basis'] = np.array([[0,-30e-6,0],
[-20e-6,0,0]]).astype(np.float32).transpose()
d1f['x_pixel_size'] = np.float32(20e-6)
d1f['y_pixel_size'] = np.float32(30e-6)
d1e['corner'] = np.array((2550e-6,3825e-6,0.3)).astype(np.float32)
d1f.create_dataset('corner_position', data=d1e['corner'])
# Remember the format for the CXI file differs from the format used
# internally
expected['mask'] = None
data1f = e1f.create_group('data_1')
data = np.random.rand(100,256,256).astype(np.float32)
expected['data'] = data
d1f.create_dataset('data',data=data)
expected['axes'] = None
g1f = s1f.create_group('geometry_1')
translations = np.arange(300).reshape((100,3)).astype(np.float32)
g1f.create_dataset('translation',data=translations)
expected['translations'] = -translations.transpose()
return f, expected
@pytest.fixture(scope='module')
def ptycho_cxi_3():
"""Creates an example file for CXI ptychography. This file is defined
to have a different subset of information missing. In particular, it:
* Has no sample_1 group
* Defines the energy but not the wavelength of light
* Has the data only defined under the data_1 group, not in the relevant places
* Defines the detector basis but no pixel sizes
* Defines a mask as all pixels flagged as "above the background"
* Defines the sample to detector distance but no corner location
* Is missing some of the allowed metadata
"""
expected = {}
f = h5py.File('ptycho_cxi_3',driver='core',backing_store=False)
# Start by defining the basic structure
f.create_dataset('cxi_version', data=150)
f.create_dataset('number_of_entries',data=1)
# Then define a bunch of metadata for entry_1
e1f = f.create_group('entry_1')
expected['entry metadata'] = {}
e1e = expected['entry metadata']
e1e['start_time'] = datetime.datetime.now().isoformat()
e1f['start_time'] = np.string_(e1e['start_time'])
e1e['end_time'] = datetime.datetime.now().isoformat()
e1f['end_time'] = np.string_(e1e['end_time'])
# Set up the sample info
expected['sample info'] = None
i1f = e1f.create_group('instrument_1')
source1f = i1f.create_group('source_1')
energy = np.float32(1.3618e-16) #Joules, = 850 eV
source1f['energy'] = energy
expected['wavelength'] = np.float32(1.9864459e-25) / energy
d1f = i1f.create_group('detector_1')
expected['detector'] = {}
d1e = expected['detector']
d1e['distance'] = np.float32(0.3)
d1f['distance'] = d1e['distance']
d1e['basis'] = np.array([[0,-30e-6,0],
[-20e-6,0,0]]).astype(np.float32).transpose()
d1f.create_dataset('basis_vectors',data=d1e['basis'])
d1e['corner'] = None
# Remember the format for the CXI file differs from the format used
# internally
mask = np.ones((100,256,256)).astype(np.uint32) * 0x00001000
expected['mask'] = np.ones((100,256,256)).astype(np.uint8)
d1f.create_dataset('mask',data=mask)
data1f = e1f.create_group('data_1')
data = np.random.rand(100,256,256).astype(np.float32)
expected['data'] = data
data1f.create_dataset('data',data=data)
data1f['data'].attrs['axes'] = np.string_('translation:y:x')
expected['axes'] = ['translation','y','x']
translations = np.arange(300).reshape((100,3)).astype(np.float32)
data1f.create_dataset('translation',data=translations)
expected['translations'] = -translations.transpose()
return f, expected
@pytest.fixture(scope='module')
def real_ptycho_CSX(datadir):
"""Loads a real ptychography file from CSX @ NSLS-II along with a
dictionary describing what is expected to be loaded from it.
"""
#file_obj = h5py.File(datadir/'example_CSX_Bragg.cxi','r')
#expected = {}
#return file_obj, expected
pass
@pytest.fixture(scope='module')
def real_ptycho_HXN(datadir):
"""Loads a real ptychography file from HXN @ NSLS-II along with a
dictionary describing what is expected to be loaded from it.
"""
pass
@pytest.fixture(scope='module')
def real_ptycho_COSMIC(datadir):
"""Loads a real ptychography file from COSMIC @ ALS along with a
dictionary describing what is expected to be loaded from it.
"""
pass
@pytest.fixture(scope='module')
def test_ptycho_cxis(ptycho_cxi_1, ptycho_cxi_2, ptycho_cxi_3):
"""Loads a list of tuples of ptychography CXI files and dictionaries,
describing the expected output from various functions on being called
on the cxi files.
"""
return [ptycho_cxi_1, ptycho_cxi_2, ptycho_cxi_3]
#
# Now we have a bunch of tests of the data loading capabilities
#
def test_get_entry_info(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
entry_info = data.get_entry_info(cxi)
for key in expected['entry metadata']:
assert entry_info[key] == expected['entry metadata'][key]
def test_get_sample_info(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
sample_info = data.get_sample_info(cxi)
if sample_info is None and \
('sample info' not in expected or
expected['sample info'] is None):
# Valid if no sample info is defined at all
continue
for key in expected['sample info']:
if isinstance(expected['sample info'][key],np.ndarray):
assert np.allclose(sample_info[key],
expected['sample info'][key])
else:
assert sample_info[key] == expected['sample info'][key]
def test_get_wavelength(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
assert np.isclose(expected['wavelength'],data.get_wavelength(cxi))
def test_get_detector_geometry(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
distance, basis, corner = data.get_detector_geometry(cxi)
assert np.isclose(distance,expected['detector']['distance'])
assert np.allclose(basis,expected['detector']['basis'])
if isinstance(expected['detector']['corner'], np.ndarray):
assert np.allclose(corner, expected['detector']['corner'])
else:
assert corner == expected['detector']['corner']
def test_get_mask(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
mask = data.get_mask(cxi)
if expected['mask'] is None and mask is None:
continue
assert np.all(data.get_mask(cxi) == expected['mask'])
def test_get_data(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
patterns, axes = data.get_data(cxi)
assert np.allclose(patterns, expected['data'])
assert axes == expected['axes']
def test_get_ptycho_translations(test_ptycho_cxis):
for cxi, expected in test_ptycho_cxis:
assert np.allclose(data.get_ptycho_translations(cxi),
expected['translations'])
#
# Then, write a test for the data saving. It should create a .cxi file
# using the data seving tools, and then check that when read with the
# .cxi reading tools that it gets the same things that were written.
#