Initial commit
This commit is contained in:
frederic committed
commit
2ae821883c
10 files changed
+1841
No files matched your search
@@ -0,0 +1,11 @@
|
||||
from importlib.metadata import PackageNotFoundError, version
|
||||
|
||||
from .spexical import SpeXICAL
|
||||
from .spexiraw import SpeXIRAW
|
||||
|
||||
try:
|
||||
__version__ = version("spexipy")
|
||||
except PackageNotFoundError: # Running from source tree?
|
||||
__version__ = "0.0.0"
|
||||
|
||||
__all__ = ["SpeXICAL", "SpeXIRAW", "__version__"]
|
||||
@@ -0,0 +1,123 @@
|
||||
import h5py
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
||||
from .spexical import SpeXICAL
|
||||
|
||||
|
||||
def cube_inspect(file: str | h5py.File, e_range: None | tuple[float, float] = None,
|
||||
pixels: None | list[tuple[int, int]] = [(10, 10), (69, 69), (10, 69), (69, 10)]):
|
||||
if isinstance(file, str):
|
||||
data = h5py.File(file, 'r')
|
||||
elif isinstance(file, h5py.File):
|
||||
data = file
|
||||
else:
|
||||
raise ValueError(f"Invalid type for file parameter: {type(file)}")
|
||||
|
||||
all_px = np.sum(data['cube/cube'][:, :, :], axis=(1, 2))
|
||||
|
||||
energies = data['cube/energies'][:]
|
||||
|
||||
fig, ax = plt.subplots()
|
||||
ax.plot(energies, all_px, label='All pixels')
|
||||
for px in pixels:
|
||||
ax.plot(energies, data['cube/cube'][:, px[1], px[0]], label=f"({px[0]},{px[1]})")
|
||||
ax.set_yscale('log')
|
||||
if e_range:
|
||||
ax.set_xlim(*e_range)
|
||||
ax.legend()
|
||||
plt.show()
|
||||
|
||||
if isinstance(file, str):
|
||||
data.close()
|
||||
del data
|
||||
|
||||
|
||||
def cube_range_to_bins(data: h5py.File, e_lower: float, e_upper: float) -> tuple[int, int]:
|
||||
bin_count = len(data['cube/energies'])
|
||||
bin_size = data['cube/energies'][1] - data['cube/energies'][0]
|
||||
offset = data['cube/energies'][0] - (bin_size / 2.)
|
||||
|
||||
return (
|
||||
min(max(0, int(round((e_lower - offset) / bin_size))), bin_count),
|
||||
min(max(0, int(round((e_upper - offset) / bin_size))), bin_count)
|
||||
)
|
||||
|
||||
|
||||
def cube_mean_in_range(file: str | h5py.File, e_range: tuple[float, float]) -> float:
|
||||
if isinstance(file, str):
|
||||
data = h5py.File(file, 'r')
|
||||
elif isinstance(file, h5py.File):
|
||||
data = file
|
||||
else:
|
||||
raise ValueError(f"Invalid type for file parameter: {type(file)}")
|
||||
|
||||
bins = cube_range_to_bins(data, *e_range)
|
||||
energies = data['cube/energies'][bins[0]:bins[1]]
|
||||
weighted_sum = np.tensordot(energies, data['cube/cube'][bins[0]:bins[1], :, :], axes=(0, 0))
|
||||
total_counts = np.sum(data['cube/cube'][bins[0]:bins[1], :, :], axis=0)
|
||||
|
||||
with np.errstate(divide='ignore', invalid='ignore'):
|
||||
mean_energy = np.where(total_counts > 0,
|
||||
weighted_sum / total_counts,
|
||||
np.nan)
|
||||
|
||||
mean_energy[total_counts == 0] = np.nanmean(mean_energy)
|
||||
|
||||
if isinstance(file, str):
|
||||
data.close()
|
||||
del data
|
||||
|
||||
return mean_energy
|
||||
|
||||
|
||||
def cal_adjust(cal: SpeXICAL, lower: np.array, upper: np.array):
|
||||
old_gain = cal.pixel_gain.copy()
|
||||
old_offset = cal.pixel_offset.copy()
|
||||
|
||||
lower_mean = np.nanmean(lower)
|
||||
upper_mean = np.nanmean(upper)
|
||||
|
||||
pixel_gain = (lower - upper) / (lower_mean - upper_mean)
|
||||
pixel_offset = (lower_mean * upper - upper_mean * lower) / (lower_mean - upper_mean)
|
||||
|
||||
pixel_gain[np.where(np.isnan(pixel_gain))] = np.nanmean(pixel_gain)
|
||||
pixel_offset[np.where(np.isnan(pixel_offset))] = np.nanmean(pixel_offset)
|
||||
|
||||
# Merge with the old gains and offsets, because the new data has been obtained with those already in place
|
||||
cal.pixel_gain = (pixel_gain * old_gain).astype(np.float32)
|
||||
cal.pixel_offset = ((pixel_offset * old_gain) + old_offset).astype(np.float32)
|
||||
|
||||
|
||||
def cal_adjust_scale(cal: SpeXICAL, e_real: tuple[float, float], e_measured: tuple[float, float]):
|
||||
old_gain = cal.global_gain
|
||||
old_offset = cal.global_offset
|
||||
|
||||
global_gain = (e_measured[0] - e_measured[1]) / (e_real[0] - e_real[1]) * old_gain
|
||||
global_offset = (e_real[0] * e_measured[1] - e_real[1] * e_measured[0]) / (
|
||||
e_real[0] - e_real[1]) * old_gain + old_offset
|
||||
|
||||
cal.global_gain = global_gain
|
||||
cal.global_offset = global_offset
|
||||
|
||||
|
||||
def cal_adjust_per_size(cal: SpeXICAL, peaks: list[tuple[float, float]]):
|
||||
base = peaks[0]
|
||||
|
||||
out = []
|
||||
|
||||
for peak in peaks[1:]:
|
||||
gain = (peak[0] - peak[1]) / (base[0] - base[1])
|
||||
offset = (base[0] * peak[1] - base[1] * peak[0]) / (base[0] - base[1])
|
||||
|
||||
i = len(out)
|
||||
if len(cal.cluster_calibrations) > i:
|
||||
old_gain = cal.cluster_calibrations[i][1]
|
||||
old_offset = cal.cluster_calibrations[i][0]
|
||||
gain *= old_gain
|
||||
offset *= old_gain
|
||||
offset += old_offset
|
||||
|
||||
out.append((offset, gain))
|
||||
|
||||
cal.cluster_calibrations = out
|
||||
@@ -0,0 +1,130 @@
|
||||
import numpy as np
|
||||
import struct
|
||||
|
||||
class SpeXICAL:
|
||||
def __init__(self, path):
|
||||
self.is_open = False
|
||||
self._supported_versions = [2, 3]
|
||||
|
||||
self._tag = b'SPEXICAL'
|
||||
self._header_struct = struct.Struct('<8sBI36sxHHBffB')
|
||||
self._cluster_cal_struct = struct.Struct('<ff')
|
||||
|
||||
self._path = None
|
||||
self._file = None
|
||||
|
||||
self._open(path)
|
||||
|
||||
def _open(self, path):
|
||||
self._path = path
|
||||
self._file = open(path, mode='rb')
|
||||
|
||||
self._file.seek(0, 2) # move cursor to end
|
||||
self.file_size = self._file.tell()
|
||||
self._file.seek(0, 0) # move back to start
|
||||
|
||||
tag = self._file.read(len(self._tag))
|
||||
if tag != self._tag:
|
||||
raise RuntimeError("File does not look like a SpeXIDAQ calibration file")
|
||||
|
||||
version = struct.unpack('<B', self._file.read(1))[0]
|
||||
if version not in self._supported_versions:
|
||||
raise RuntimeError(f"File uses format version {version} but only versions {self._supported_versions} are supported")
|
||||
|
||||
if self.file_size < self._header_struct.size:
|
||||
raise RuntimeError("File too short to be a SpeXIDAQ calibration file")
|
||||
|
||||
self._file.seek(0, 0) # back to the start
|
||||
tag, self.format_version, self.version, self.uuid, self.width, \
|
||||
self.height, self.detector_gain, self.global_gain, \
|
||||
self.global_offset, self._cluster_cal_count = self._header_struct.unpack_from(self._file.read(self._header_struct.size))
|
||||
|
||||
self._px_count = self.width * self.height
|
||||
|
||||
self._calc_min_size()
|
||||
if self.file_size < self._min_size:
|
||||
raise RuntimeError("File too short to be a SpeXIDAQ calibration file")
|
||||
|
||||
self.cluster_calibrations = []
|
||||
if self.format_version >= 3:
|
||||
for i in range(self._cluster_cal_count):
|
||||
offset, gain = self._cluster_cal_struct.unpack_from(self._file.read(self._cluster_cal_struct.size))
|
||||
self.cluster_calibrations.append((offset, gain))
|
||||
else:
|
||||
cluster_cal_struct = struct.Struct('<f')
|
||||
for i in range(self._cluster_cal_count):
|
||||
offset = cluster_cal_struct.unpack_from(self._file.read(cluster_cal_struct.size))
|
||||
self.cluster_calibrations.append(offset)
|
||||
|
||||
self.dark_offset = self._load_frame(np.float32)
|
||||
self.thresholds = self._load_frame(np.float32)
|
||||
|
||||
self.pixel_gain = self._load_frame(np.float32)
|
||||
self.pixel_offset = self._load_frame(np.float32)
|
||||
|
||||
if self.format_version < 3:
|
||||
self.fine_gain = self._load_frame(np.float32)
|
||||
self.fine_offset = self._load_frame(np.float32)
|
||||
|
||||
self.pixel_mask = self._load_frame(np.uint16)
|
||||
|
||||
self.comment = self._file.read(self.file_size - self._min_size)
|
||||
|
||||
self._file.close()
|
||||
self._file = None
|
||||
self.is_open = True
|
||||
|
||||
def _calc_min_size(self):
|
||||
self._min_size = self._header_struct.size + self._cluster_cal_struct.size * self._cluster_cal_count
|
||||
self._min_size += self._px_count * (4 + 4 + 4 + 4 + 2) + 4 * 5
|
||||
self._min_size += 1 # null-termination of comment string
|
||||
|
||||
def _load_frame(self, dtype):
|
||||
w, h = struct.unpack('<HH', self._file.read(4))
|
||||
if w != self.width or h != self.height:
|
||||
raise RuntimeError(f"Unexpected frame size {w}x{h} in calibration, expected {self.width}x{self.height}")
|
||||
frame = np.frombuffer(self._file.read(self._px_count * np.dtype(dtype).itemsize), dtype=dtype).copy()
|
||||
frame.shape = (self.width, self.height)
|
||||
|
||||
return frame
|
||||
|
||||
def _store_frame(self, frame):
|
||||
if frame.shape != (self.width, self.height):
|
||||
raise RuntimeError("Trying to store unsupported frame shape")
|
||||
self._file.write(struct.pack('<HH', self.width, self.height))
|
||||
self._file.write(frame.tobytes())
|
||||
|
||||
def save_as(self, path):
|
||||
self._path = path
|
||||
self.save()
|
||||
|
||||
def save(self):
|
||||
self._file = open(self._path, mode='wb')
|
||||
|
||||
header_buf = self._header_struct.pack(self._tag, self.format_version, self.version, self.uuid, self.width, self.height, self.detector_gain, self.global_gain, self.global_offset, len(self.cluster_calibrations))
|
||||
|
||||
self._file.write(header_buf)
|
||||
|
||||
if self.format_version >= 3:
|
||||
for cal in self.cluster_calibrations:
|
||||
self._file.write(self._cluster_cal_struct.pack(cal[0], cal[1]))
|
||||
else:
|
||||
cluster_cal_struct = struct.Struct('<f')
|
||||
for cal in self.cluster_calibrations:
|
||||
self._file.write(cluster_cal_struct.pack(cal))
|
||||
|
||||
self._store_frame(self.dark_offset)
|
||||
self._store_frame(self.thresholds)
|
||||
self._store_frame(self.pixel_gain)
|
||||
self._store_frame(self.pixel_offset)
|
||||
|
||||
if self.format_version < 3:
|
||||
self._store_frame(self.fine_gain)
|
||||
self._store_frame(self.fine_offset)
|
||||
|
||||
self._store_frame(self.pixel_mask)
|
||||
|
||||
self._file.write(self.comment)
|
||||
self._file.write(b'\x00')
|
||||
self._file.close()
|
||||
self._file = None
|
||||
@@ -0,0 +1,83 @@
|
||||
import numpy as np
|
||||
import struct
|
||||
|
||||
class SpeXIRAW:
|
||||
def __init__(self, path):
|
||||
# initialise some defaults
|
||||
self.is_open = False
|
||||
self._supported_version = 1
|
||||
|
||||
# define some helper values and structs for reading
|
||||
self._tag = b'SPEXIRAW'
|
||||
self._header_struct = struct.Struct('<8sBHHfIB36s8s')
|
||||
self._frame_id_struct = struct.Struct('<Id')
|
||||
self._types = {1: 'uint8',
|
||||
2: 'uint16',
|
||||
3: 'uint32',
|
||||
4: 'uint64'}
|
||||
|
||||
self._open(path)
|
||||
|
||||
def _open(self, path):
|
||||
self._file = open(path, mode='rb')
|
||||
|
||||
self._file.seek(0, 2) # move cursor to end
|
||||
self.file_size = self._file.tell()
|
||||
self._file.seek(0, 0) # move back to start
|
||||
|
||||
tag1 = self._file.read(8)
|
||||
if tag1 != self._tag:
|
||||
raise RuntimeError("File does not look like a SpeXIRAW file")
|
||||
|
||||
version = struct.unpack('<B', self._file.read(1))[0]
|
||||
if version != self._supported_version:
|
||||
raise RuntimeError(f"File uses format version {version} but only version {self._supported_version} is supported")
|
||||
|
||||
if self.file_size < self._header_struct.size:
|
||||
raise RuntimeError("File too short to be a SpeXIRAW file")
|
||||
|
||||
self._file.seek(0, 0) # back to the start
|
||||
tag1, self.format_version, self.width, self.height, self.pixel_pitch, \
|
||||
self.fps, self.type, self.detector_uuid, tag2 = self._header_struct.unpack_from(self._file.read(self._header_struct.size))
|
||||
|
||||
if tag2 != self._tag:
|
||||
raise RuntimeError("File is probably corrupted")
|
||||
|
||||
self._dtype = np.dtype(self._types[self.type])
|
||||
self._frame_size = self.width * self.height * self._dtype.itemsize
|
||||
self.frame_count = int(np.floor((self.file_size - self._header_struct.size) / (self._frame_size + self._frame_id_struct.size)))
|
||||
|
||||
self.is_open = True
|
||||
|
||||
def frames(self):
|
||||
if not self.is_open:
|
||||
return
|
||||
|
||||
self._file.seek(self._header_struct.size, 0) # move to start of frame data
|
||||
for n in range(self.frame_count):
|
||||
#self._file.seek(self._frame_id_struct.size, 1) # skip header, unused in iterator
|
||||
index, timestamp = self._frame_id_struct.unpack(self._file.read(self._frame_id_struct.size))
|
||||
frame = np.fromfile(self._file, dtype=self._dtype, count=self.width * self.height)
|
||||
frame.shape = (self.width, self.height)
|
||||
yield index, timestamp, frame
|
||||
|
||||
def frame(self, index):
|
||||
if not self.is_open:
|
||||
raise RuntimeError("No open file")
|
||||
|
||||
if not isinstance(index, int):
|
||||
raise ValueError("Numeric frame index required")
|
||||
|
||||
if index >= self.frame_count:
|
||||
raise ValueError("Requested frame beyond end of file")
|
||||
|
||||
if index < 0:
|
||||
index += self.frame_count
|
||||
|
||||
if index < 0:
|
||||
raise ValueError("Requested frame before start of file")
|
||||
|
||||
self._file.seek(self._header_struct.size + (self._frame_id_struct.size + self._frame_size) * index + self._frame_id_struct.size)
|
||||
frame = np.fromfile(self._file, dtype=self._dtype, count=self.width * self.height)
|
||||
frame.shape = (self.width, self.height)
|
||||
return frame
|
||||
@@ -0,0 +1,45 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from mpl_toolkits.axes_grid1 import make_axes_locatable
|
||||
|
||||
|
||||
def plot_frame(frame, title=None, out=None):
|
||||
if out is None:
|
||||
out = plt.figure()
|
||||
|
||||
if isinstance(out, plt.Figure):
|
||||
ax = out.add_subplot(1, 1, 1)
|
||||
else:
|
||||
ax = out
|
||||
|
||||
im = ax.imshow(frame)
|
||||
ax.set_xlabel('X')
|
||||
ax.set_ylabel('Y')
|
||||
|
||||
divider = make_axes_locatable(ax)
|
||||
cax = divider.append_axes("right", size="5%", pad=0.05)
|
||||
plt.colorbar(im, cax=cax)
|
||||
|
||||
if title is not None:
|
||||
if out is plt:
|
||||
ax.title(title)
|
||||
else:
|
||||
ax.set_title(title)
|
||||
|
||||
|
||||
def plot_frames(frames):
|
||||
if isinstance(frames, np.ndarray):
|
||||
plot_frame(frames)
|
||||
elif len(frames) == 1:
|
||||
plot_frame(frames[0])
|
||||
else:
|
||||
rows = 1
|
||||
cols = len(frames)
|
||||
if len(frames) > 3:
|
||||
rows = ceil(len(frames) / 2)
|
||||
cols = 2
|
||||
fig, sub = plt.subplots(rows, cols, sharey=True, sharex=True, figsize=(14, 5 * rows))
|
||||
for i in range(len(frames)):
|
||||
plot_frame(frames[i], out=sub[i])
|
||||
|
||||
plt.tight_layout()
|
||||
Reference in new issue
Block a user