Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 48 additions & 33 deletions brukerapi/dataset.py
Original file line number Diff line number Diff line change
Expand Up @@ -245,11 +245,7 @@ def __init__(self, path, **state):
(name for name in content if re.fullmatch(r"rawdata\.job\d+", name)),
key=lambda name: int(name.rsplit("job", 1)[1]),
)
named_rawdata_jobs = sorted(
name
for name in content
if name not in rawdata_jobs and self.is_supported_path(self.path / name, {"rawdata"})
)
named_rawdata_jobs = sorted(name for name in content if name not in rawdata_jobs and self.is_supported_path(self.path / name, {"rawdata"}))
if "fid" in content:
self.path = self.path / "fid"
elif "2dseq" in content:
Expand All @@ -273,10 +269,7 @@ def __init__(self, path, **state):
if self.type not in DEFAULT_STATES:
raise UnsupportedDatasetType(self.type)
if self.type == "fid" and self.subtype in FID_COMPANION_SUBTYPES and not state.get("_auxiliary"):
raise UnsupportedDatasetType(
f"{self.path.name} is a fid companion (spec §3.5); load it as "
f"Dataset({with_suffix(self.path, '')!s}).fid_companions[{self.subtype!r}]"
)
raise UnsupportedDatasetType(f"{self.path.name} is a fid companion (spec §3.5); load it as Dataset({with_suffix(self.path, '')!s}).fid_companions[{self.subtype!r}]")
if not self._is_supported_subtype() and not (state.get("_auxiliary") and self.type == "fid" and self.subtype in FID_COMPANION_SUBTYPES):
raise UnsupportedDatasetType(self.path.name)

Expand Down Expand Up @@ -596,9 +589,9 @@ def get(self, name, default=None):
A name no configuration declares still raises, so a typo does not
quietly become `default`::

dataset.get("TE") # 4.0, or None where there is no echo time
dataset.get("TE", "n/a") # 4.0, or "n/a"
dataset.get("TEE") # AttributeError
dataset.get("TE") # 4.0, or None where there is no echo time
dataset.get("TE", "n/a") # 4.0, or "n/a"
dataset.get("TEE") # AttributeError

Parameters are a separate layer, read with ``dataset[key]`` or
``dataset.parameters``.
Expand Down Expand Up @@ -727,28 +720,56 @@ def rawdata_stored_scans(self):
settings_stored = int(settings[9])
if settings_stored != job_stored and not getattr(self, "_warned_rawdata_stored_scans", False):
warnings.warn(
f"ACQ_ScanPipeJobSettings nStoredScans={settings_stored} disagrees with "
f"ACQ_jobs nStoredScans={job_stored} for {self.path}; using the settings value",
f"ACQ_ScanPipeJobSettings nStoredScans={settings_stored} disagrees with ACQ_jobs nStoredScans={job_stored} for {self.path}; using the settings value",
RuntimeWarning,
stacklevel=2,
)
self._warned_rawdata_stored_scans = True
return settings_stored

@property
def rawdata_job_channel(self):
"""The RF channel this rawdata job is acquired on -- ``chanNum``, spec 3.3.

Only the nine-field PV360 form of ``ACQ_jobs`` carries it; the eight-field
PV6/PV7 form has no channel field, and there is one channel to select from.
"""
index = self._rawdata_job_index()
try:
jobs = self["ACQ_jobs"].nested
except KeyError:
return 1
if index is None or index >= len(jobs) or len(jobs[index]) != 9:
return 1
return int(jobs[index][7])

@property
def rawdata_channels(self):
"""Number of receivers recorded for this job (spec 3.3)."""
"""Number of receivers stored for this job (spec 3.3, 14.4).

Spec 3.3 keys the count on the job's own channel: look up
``c = ACQ_jobs[n].chanNum``, then count the ``Yes`` entries of
``ACQ_ReceiverSelectPerChan`` row ``c-1``. The flat ``ACQ_ReceiverSelect``
gives the same count for the common single-channel case, and
``PVM_EncNReceivers`` is the method-side mirror -- a last resort, not the
arbiter.
"""
selected = self._parameter_value("ACQ_ReceiverSelectPerChan")
fallback = int(self._parameter_value("PVM_EncNReceivers", 1))
if selected is not None:
values = np.atleast_1d(selected)
channels = sum(str(value).casefold() in {"yes", "on", "1"} for value in values)
# Some systems expose more physical receiver paths than the job
# writes. Only promote the per-channel declaration when it agrees
# with the method's logical receiver count.
if channels == fallback:
return channels
return fallback
selected = np.atleast_2d(selected)
channel = self.rawdata_job_channel
row = channel - 1 if 0 < channel <= selected.shape[0] else 0
return self._count_selected(selected[row])

selected = self._parameter_value("ACQ_ReceiverSelect")
if selected is not None:
return self._count_selected(np.atleast_1d(selected))

return int(self._parameter_value("PVM_EncNReceivers", 1))

@staticmethod
def _count_selected(values):
return sum(str(value).casefold() in {"yes", "on", "1"} for value in values)

def _infer_scheme_id(self):
# Source of truth for format inference:
Expand Down Expand Up @@ -779,9 +800,7 @@ def _infer_scheme_id(self):
if descriptions and descriptions[0] == "Spectroscopic":
return "SPECTROSCOPY" if int(dim) == 1 else "CSI"

if n_projections is not None and int(n_projections) > 0 and (
(self.path.parent / "traj").exists() or "RADIAL" in family or "UTE" in family
):
if n_projections is not None and int(n_projections) > 0 and ((self.path.parent / "traj").exists() or "RADIAL" in family or "UTE" in family):
return "RADIAL"

acq_size = np.atleast_1d(self._parameter_value("ACQ_size", []))
Expand Down Expand Up @@ -940,9 +959,7 @@ def _read_binary_file(self, path, dtype, shape):
f"expected {expected_bytes} bytes for shape {tuple(shape)} and dtype {dtype}"
)

raise InvalidDataset(
f"Invalid dataset size at {path}: expected {expected_bytes} bytes for shape {tuple(shape)} and dtype {dtype}, got {actual_bytes} bytes"
)
raise InvalidDataset(f"Invalid dataset size at {path}: expected {expected_bytes} bytes for shape {tuple(shape)} and dtype {dtype}, got {actual_bytes} bytes")

return read_array(path, dtype, shape)

Expand Down Expand Up @@ -1343,9 +1360,7 @@ def affine(self):
"""
descriptors = np.atleast_1d(np.asarray(self._parameter_value("VisuCoreDimDesc", []))).astype(str)
if descriptors.size and any(descriptor != "spatial" for descriptor in descriptors):
raise UnsupportedDatasetType(
f"an image affine for {self.path}, whose frames are {sorted(set(descriptors))} rather than purely spatial (spec 7.2),"
)
raise UnsupportedDatasetType(f"an image affine for {self.path}, whose frames are {sorted(set(descriptors))} rather than purely spatial (spec 7.2),")
if len(self.slice_packages_index()) > 1:
warnings.warn(
f"{self.path} has multiple slice packages; a single affine cannot describe them -- use get_slice_packages() / affine_of_package(i)",
Expand Down
25 changes: 24 additions & 1 deletion brukerapi/jcampdx.py
Original file line number Diff line number Diff line change
Expand Up @@ -267,9 +267,32 @@ def value(self):
# declared dimensions, and a mismatch is a diagnosable
# condition, not an internal error
raise InvalidJcampdxFile(f"{self.key}: {value.size} values do not fill the declared size {self.size}") from error
return value
return np.reshape(value, self._text_shape(value.size), order="C")
return value

def _text_shape(self, count):
"""The shape `count` parsed text values take under the declared size.

Spec 2.3 writes a **string** array with one length indicator per array
dimension plus a final one for the string length -- ``( 6, 65 )`` is six
strings, not 390 -- while an **enum** array carries no string-length
dimension at all, so ``ACQ_ReceiverSelectPerChan=( 2, 7 )`` really is two
channels of seven receivers. Both parse to a text dtype, and the element
count is what separates them.

Leaving every text array flat, as this used to, made the 2-D enum arrays
unindexable: spec 3.3 asks for ``ACQ_ReceiverSelectPerChan[chanNum-1]``
and there were no rows to index.
"""
size = tuple(int(length) for length in self.size)
if count == int(np.prod(size)):
return size
if len(size) > 1 and count == int(np.prod(size[:-1])):
return size[:-1]
# An unrecognised shape stays flat rather than raising: unlike the numeric
# case above, a text array has always been delivered this way.
return (count,)

@value.setter
def value(self, value):
size = self.size
Expand Down
72 changes: 18 additions & 54 deletions brukerapi/schemas.py
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,7 @@
"k_space",
"encoded_dim",
"shape_storage",
"dim_type"
"dim_type",
],
"2dseq": [
"pv_version",
Expand Down Expand Up @@ -223,9 +223,7 @@ def layouts(self):
for name in ("encoding_space", "k_space"):
logical_samples = int(np.prod(layouts[name]))
if stored_samples % logical_samples:
raise InvalidDataset(
f"real-only AQ_mod=qf sample count {stored_samples} is incompatible with {name} layout {layouts[name]}"
)
raise InvalidDataset(f"real-only AQ_mod=qf sample count {stored_samples} is incompatible with {name} layout {layouts[name]}")
ratio = stored_samples // logical_samples
if ratio > 1:
layouts[name] = (layouts[name][0] * ratio,) + tuple(layouts[name][1:])
Expand Down Expand Up @@ -282,10 +280,7 @@ def raw(self):
data = self._decode_raw_stream(stored, self.layouts)
receivers = int(self._dataset.channels)
if data.shape[0] % receivers:
raise InvalidDataset(
f"decoded FID sample count {data.shape[0]} is not divisible by "
f"the receiver count {receivers} for {self._dataset.path}"
)
raise InvalidDataset(f"decoded FID sample count {data.shape[0]} is not divisible by the receiver count {receivers} for {self._dataset.path}")
samples = data.shape[0] // receivers
return np.transpose(
np.reshape(data, (samples, receivers, data.shape[1]), order="F"),
Expand All @@ -302,9 +297,7 @@ def to_kspace(self, data=None, *, bart=False):
try:
axes = tuple(BART_DIM_BY_TYPE[label] for label in self._dataset.dim_type)
except KeyError as error:
raise UnknownAcqSchemeException(
f"cannot map FID axis {error.args[0]!r} to BART for {self._dataset.path}"
) from error
raise UnknownAcqSchemeException(f"cannot map FID axis {error.args[0]!r} to BART for {self._dataset.path}") from error
return self._as_bart(data, axes)

def _reorder_objects(self, data, dir="FW"):
Expand Down Expand Up @@ -391,9 +384,7 @@ def _reorder_fid_lines(self, data, dir="FW"):
PVM_EncSteps1_sorted = self.permutation_inverse(PVM_EncSteps1_sorted)

if data.shape[1] != len(PVM_EncSteps1_sorted):
raise InvalidDataset(
f"phase-encode reorder length {len(PVM_EncSteps1_sorted)} does not match k-space axis length {data.shape[1]} for scheme {self._dataset.scheme_id}"
)
raise InvalidDataset(f"phase-encode reorder length {len(PVM_EncSteps1_sorted)} does not match k-space axis length {data.shape[1]} for scheme {self._dataset.scheme_id}")

if np.array_equal(PVM_EncSteps1_sorted, np.arange(len(PVM_EncSteps1_sorted))):
return data
Expand Down Expand Up @@ -712,17 +703,13 @@ def to_kspace(self, data=None, *, bart=False):
scheme_id = self._dataset._infer_scheme_id()
if scheme_id is not None:
raise UnknownAcqSchemeException(
f"rawdata-to-k-space is currently supported only for Cartesian PV-360 jobs, "
f"but {self._dataset.path} is {scheme_id}; use its acquisition-specific reader"
f"rawdata-to-k-space is currently supported only for Cartesian PV-360 jobs, but {self._dataset.path} is {scheme_id}; use its acquisition-specific reader"
)

encoded_dim = self._dataset._parameter_value("ACQ_dim")
matrix = np.atleast_1d(self._dataset._parameter_value("PVM_EncMatrix", []))
if encoded_dim not in (2, 3) or matrix.size < encoded_dim:
raise UnknownAcqSchemeException(
f"cannot establish a Cartesian rawdata layout for {self._dataset.path}; "
"pass data through an acquisition-specific reader"
)
raise UnknownAcqSchemeException(f"cannot establish a Cartesian rawdata layout for {self._dataset.path}; pass data through an acquisition-specific reader")

matrix = tuple(int(value) for value in matrix[:encoded_dim])
receivers = int(self._dataset.channels)
Expand All @@ -747,22 +734,15 @@ def to_kspace(self, data=None, *, bart=False):
permute = (0, 2, 3, 4, 1)
else:
if objects != 1:
raise UnknownAcqSchemeException(
f"cannot establish a 3-D Cartesian rawdata layout with NI={objects} for "
f"{self._dataset.path}"
)
raise UnknownAcqSchemeException(f"cannot establish a 3-D Cartesian rawdata layout with NI={objects} for {self._dataset.path}")
encoding_space = (readout, receivers, phase, matrix[2], repetitions)
permute = (0, 2, 3, 4, 1)

if data.ndim != 3 or data.shape[0] != readout or data.shape[1] != receivers:
raise InvalidDataset(
f"rawdata sample layout {data.shape} does not match Cartesian metadata "
f"(readout={readout}, receivers={receivers}) for {self._dataset.path}"
)
raise InvalidDataset(f"rawdata sample layout {data.shape} does not match Cartesian metadata (readout={readout}, receivers={receivers}) for {self._dataset.path}")
if data.size != int(np.prod(encoding_space)):
raise InvalidDataset(
f"rawdata contains {data.size} complex samples but Cartesian layout {encoding_space} "
f"requires {int(np.prod(encoding_space))} for {self._dataset.path}"
f"rawdata contains {data.size} complex samples but Cartesian layout {encoding_space} requires {int(np.prod(encoding_space))} for {self._dataset.path}"
)

k_space = np.transpose(np.reshape(data, encoding_space, order="F"), permute)
Expand All @@ -779,10 +759,7 @@ def _self_gated_k_space(self, data, readout, phase, *, receivers, objects, repet
"""Arrange retrospectively gated Cartesian data before cine binning."""
steps = np.atleast_1d(self._dataset._parameter_value("PVM_EncGenSteps1", []))
if steps.size == 0 or steps.size % phase:
raise InvalidDataset(
f"self-gated phase-encode sequence has {steps.size} steps, which is incompatible "
f"with phase size {phase} for {self._dataset.path}"
)
raise InvalidDataset(f"self-gated phase-encode sequence has {steps.size} steps, which is incompatible with phase size {phase} for {self._dataset.path}")
acquired_frames = steps.size // phase
output_frames = int(self._dataset._parameter_value("PVM_NMovieFrames", objects))
if output_frames != objects or acquired_frames % (objects * repetitions):
Expand All @@ -793,10 +770,7 @@ def _self_gated_k_space(self, data, readout, phase, *, receivers, objects, repet
acquisition_cycles = acquired_frames // (objects * repetitions)
layout = (readout, receivers, phase, acquired_frames)
if data.ndim != 3 or data.shape[0] != readout or data.shape[1] != receivers or data.size != int(np.prod(layout)):
raise InvalidDataset(
f"rawdata sample layout {data.shape} does not match self-gated Cartesian layout "
f"{layout} for {self._dataset.path}"
)
raise InvalidDataset(f"rawdata sample layout {data.shape} does not match self-gated Cartesian layout {layout} for {self._dataset.path}")

k_space = np.transpose(np.reshape(data, layout, order="F"), (0, 2, 3, 1))
order = np.argsort(np.reshape(steps, (phase, acquired_frames), order="F"), axis=0)
Expand All @@ -813,10 +787,7 @@ def _reorder_phase_lines(self, data):
return data
indices = np.argsort(np.atleast_1d(steps))
if indices.size != data.shape[1]:
raise InvalidDataset(
f"phase-encode reorder length {indices.size} does not match k-space axis length "
f"{data.shape[1]} for {self._dataset.path}"
)
raise InvalidDataset(f"phase-encode reorder length {indices.size} does not match k-space axis length {data.shape[1]} for {self._dataset.path}")
return np.take(data, indices, axis=1)

def _reorder_objects(self, data):
Expand All @@ -827,12 +798,10 @@ def _reorder_objects(self, data):
return data
indices = np.argsort(np.atleast_1d(order).astype(int))
if indices.size != data.shape[2]:
raise InvalidDataset(
f"object-order length {indices.size} does not match k-space axis length "
f"{data.shape[2]} for {self._dataset.path}"
)
raise InvalidDataset(f"object-order length {indices.size} does not match k-space axis length {data.shape[2]} for {self._dataset.path}")
return np.take(data, indices, axis=2)


class Schema2dseq(Schema):
"""
Schema2dseq class
Expand Down Expand Up @@ -948,8 +917,7 @@ def _apply_disk_slice_order(self, data):
axis = 2
if axis is None or axis >= data.ndim:
warnings.warn(
"VisuCoreDiskSliceOrder requests reversed slices but no slice axis is identifiable "
f"for {getattr(self._dataset, 'path', '<unknown>')}; leaving order unchanged",
f"VisuCoreDiskSliceOrder requests reversed slices but no slice axis is identifiable for {getattr(self._dataset, 'path', '<unknown>')}; leaving order unchanged",
RuntimeWarning,
stacklevel=2,
)
Expand All @@ -974,9 +942,7 @@ def _complex_frame_axis(self, data):
# component: there is nothing to combine, keep the axis as it is
return None
if data.shape[axis] != 2:
raise InvalidDataset(
f"complex 2dseq requires a two-element real/imag frame-group axis, got shape {data.shape} on axis {axis}"
)
raise InvalidDataset(f"complex 2dseq requires a two-element real/imag frame-group axis, got shape {data.shape} on axis {axis}")
return axis

def _combine_complex_frames(self, data):
Expand Down Expand Up @@ -1100,9 +1066,7 @@ def ra(self, slice_):
# Frame-group selection above does not alter encoded dimensions.
# Apply their requested selection after deserializing the selected
# frames, preserving singleton axes for the final squeeze below.
encoded_slice = tuple(
slice(index, index + 1) if isinstance(index, int) else index for index in slice_[: self._dataset.encoded_dim]
)
encoded_slice = tuple(slice(index, index + 1) if isinstance(index, int) else index for index in slice_[: self._dataset.encoded_dim])
array_ra = array_ra[encoded_slice + (slice(None),) * (array_ra.ndim - self._dataset.encoded_dim)]

singletons = tuple(i for i, v in enumerate(slice_) if isinstance(v, int))
Expand Down
Loading