[Bug]: mirror_odd_lines incorrectly applied to non-EPI sequences with RECO_inp_order='REV_ALT_ROWS'

[schuenke]

Before submitting

• I searched existing issues and did not find a duplicate
• If reporting dataset incompatiblity, I am willing to provide the data
• I am using the latest available version

Bug description

mirror_odd_lines reverses every other phase-encode line's readout when RECO_inp_order='REV_ALT_ROWS' is set. This is correct for symmetric EPI (where adjacent phase lines within one echo train have alternating readout polarity), but wrong for bipolar multi-echo sequences like MGE.

For the MGE example sequence in the PV7 Bruker API test data set (https://zenodo.org/records/4522220, exp number 26), the current code warns about a disagreement between the declared and inferred value, but uses the declared value anyway leading to wrong kspace data.

For PV360 MGE data, everything is working fine as it is, because the SchemaRawdata has no mirror_odd_lines logic

Steps to reproduce

When using this simple code snippet:

import numpy as np
from brukerapi.dataset import Dataset

data = Dataset(<path_to_MGE_data>)
img = np.fft.ifftshift(np.fft.ifftn(np.fft.fftshift(data.kspace, axes=(0, 1)), axes=(0, 1)), axes=(0, 1))

plt.imshow(np.abs(img[:, :, 16, 0, 0]))
plt.show()

the current output looks like this:

When adjusting the mirror_odd_lines() function to ignore the declared value for non-EPI sequences like this:

@property
def mirror_odd_lines(self):
    """Whether §6.2 odd-line reversal applies to phase-encode lines."""
    value = self._dataset._parameter_value("RECO_inp_order")
    declared = None if value is None else str(value).upper() == "REV_ALT_ROWS"
    inferred = "EPI" in self._dataset.scheme_id
    if not inferred and declared:
        return False
    return self._declared_or_inferred(
        "RECO_inp_order", declared, inferred, "_warned_reco_input_order_disagreement"
    )

the image looks correct:

I am not sure if you prefer to fix this in the mirror_odd_lines function or somewhere else, but ofc I can create a PR if you agree.

Project version

0.4.6

Paravision version

PV7

[headmeister]

Dear Patrick,

Let me think this through, and please correct me if I am wrong.

With a unipolar readout, it is clear: we do not want any flipping.

With a bipolar readout, however, consecutive echoes are acquired with opposite readout gradient polarity. So rather than blindly flipping every other k-space phase-encode step (e.g. every odd ky line), we should flip the readout direction of every other echo (e.g. every odd echo).

In other words, for bipolar MGE we should not omit the flipping altogether; we should apply it along the echo dimension (as in EPI case) rather than the phase-encode dimension.

So, if we already have this phase flipping functionality, then we would likely want this to work as I described, rather than omit phase flips altogether. Therefore, we would likely rather adjust the flipping in the _mirror_odd_lines, where instead of blindly flipping every ky line, we need to flip every other echo...

Jiri

[gdevenyi]

I can look into this later, but I'd like any comments you have on this versus the FILE_FORMAT.md document I constructed from manuals/example data.

https://github.com/CoBrALab/pvraw/blob/main/FILE_FORMAT.md

[headmeister]

Hi Gabriel,

yes, I believe my explanation is in line with Section 6.2. It states that RECO_inp_order should reverse every second scan—in other words, every other readout in acquisition order.

The main question from an implementation point of view is where to apply this correction. For MGE, or any other acquisition where the readout is separated into individual echoes, we could simply apply it at the beginning of the import and flip every other acquired readout/echo.

For EPI, however, this is a bit more complicated because a single shot usually contains multiple k-space lines, or even the entire k-space. In that case, the alternating readout polarity exists between the individual EPI lines within the shot, so the correction needs to be applied after the EPI data have been reordered into Cartesian k-space.

I still have not entirely decided how far we should go with BrukerAPI with regard to RECO. While we likely reorder EPI data correctly into k-space, we do not perform any ghost correction or regridding based on the measured trajectory. Therefore, anyone trying to use BrukerAPI's k-space output directly, especially for EPI, will see degraded images anyway.

There are also likely some aspects that we simply cannot infer from the reco files. Although we have information about the reconstruction graph, the internal implementation of essentially all ParaVision reconstruction filters is a black box, so in many cases we can only guess what they actually do internally.

Best,
Jiri

[headmeister]

@gdevenyi
Speaking of the FILE_FORMAT.md document, I would really like to reference it in BrukerAPI. I originally referenced the version that was part of brkraw-legacy, but I see that this now redirects to the version in pvraw.

Would it be possible to establish one maintained source of truth for this document? I don't particularly mind where it lives—one of your repositories, pvraw, or even the BrukerAPI repository would all be fine with me. I would just like to have a single canonical version of the file-format documentation that we can reference and contribute to, rather than having several versions that may diverge over time.

What do you think?

Best,
Jiri

[gdevenyi]

Sure, absolutely.

The one in pvraw is the one from brkraw-legacy, the redirect is correct, I re-branded in prep for proper release as I think its ready for usage. No other version exists.

It is 100% my intention to continue to maintain that copy, I have plans for a "call for more data" and to reach out to bruker directly for more manuals and examples to further refine it.

[headmeister]

Ok, I will reference the one in pvraw then.

Fingers crossed with Bruker. Our experience with them regarding manuals has usually been a bit rough. The manuals have been getting progressively less detailed since PV5, with less and less information available, and they are usually reluctant to provide any internal documentation (if it exists at all).

But you may have better luck!

[gdevenyi]

I reviewed this against FILE_FORMAT.md and the local test corpus (PV6.0.1 exp 6, PV7 exp 26, PV360 3.6 T2star_map_MGE_mod_all). Short version: the declared RECO_inp_order=REV_ALT_ROWS is correct for bipolar MGE. The defect is the axis on which brukerapi applies it. @headmeister's reading is correct: reverse every second echo, not every second k-space line.

What the spec says

1. §6.2 Input Reordering: REV_ALT_ROWS = "Reverse every second scan prior to processing (first direction only)". The unit is one scan of the raw stream, in acquisition order. It is not one sorted k-space line.
2. §5.2 Scan Loop Structure: the default nesting is NS > ACQ_phase_factor > NSLICES > NI or NSLICES > NA > ACQ_size[1]/ACQ_phase_factor > .... For MGE, ACQ_phase_factor=1, ACQ_scan_size=One_scan and NI = NSLICES * ACQ_n_echo_images (PV7 exp 26: 40 = 5 * 8). ACQ_obj_order = 0 1 2 3 4 5 6 7 16 17 ... 23 32 ... 39 8 ... 15 24 ... 31 shows that the objects of one slice are its 8 echoes, stored one after the other. So "every second scan" in fid is every odd echo.
3. §10.1 Standard Reconstruction: the first pass runs "4. Reversal of alternating rows (only if REV_ALT_ROWS)" before "6. Sorting of the profiles to resolve phase encoding and slice ordering". The on-disk graph (§10.2, RecoReverseFilter) shows the same order in the PV7 MGE reco and in the PV7 EPI reco: CAST0 -> REV0 -> FTS0 -> FT0 -> REORD0. The node is RecoReverseFilter REV0{segment=<RecoReverseSegment>}, with RecoReverseSegment=0 for MGE and 24 for the 8-shot EPI (exp 13).
4. §14.6 Acquisition schemes and §5.3: the EPI scheme is the one with ACQ_scan_size = ACQ_phase_factor_scans. MGE is One_scan. Scheme inference is therefore not the correct switch for this reversal. RECO_inp_order tracks the readout polarity directly.

What the data says

Every MGE reco in the corpus:

ParaVision	EchoAcqMode	RECO_inp_order	RecoReverseSegment
5.1 (0.2H2/3)	positiveReadOutEchoes	NO_REORDERING	absent
6.0.1 (exp 6)	allEchoes	REV_ALT_ROWS	0
7.0.0 (exp 26)	allEchoes	REV_ALT_ROWS	0
360 3.6 (T2star_map_MGE, _mod_pos)	positiveReadOutEchoes	NO_REORDERING	0
360 3.6 (T2star_map_MGE_mod_all)	allEchoes	REV_ALT_ROWS	0

Bipolar (allEchoes) always declares REV_ALT_ROWS; unipolar never does. The declaration is trustworthy. The "disagrees with scheme_id" warning is a false alarm for this case.

Numeric check against 2dseq (correlation of a plain 2-D iFFT magnitude with the 2dseq frame; orientation locked from echo 0; mean over even / odd echoes):

Data	current (odd-ky mirror)	no reversal (proposed patch)	reverse readout of odd echoes
PV6.0.1 exp 6, 6 echoes	0.58 / 0.50	1.000 / 0.526	1.000 / 1.000
PV7 exp 26, 8 echoes	0.42 / 0.53	0.673 / 0.497	0.673 / 0.673
PV360 3.6 _mod_all, rawdata.job0	n/a	0.772 / 0.235	0.772 / 0.758

So the proposed patch (return False when the scheme is not EPI) leaves each odd echo mirrored along the readout axis. Object 16 in the snippet is slice 2, echo 0. That is an even echo, so the picture looks correct with and without the patch. Object 17 (echo 1) shows the difference.

PV360 is not "working fine" either. SchemaRawdata returns the (sample, channel, scan) stream and applies no reversal, so the odd echoes of a bipolar PV360 MGE are still reversed in the stream. That is by design for a raw view, but a user who reshapes it must apply §6.2 themselves.

One more finding on the EPI path. _mirror_odd_lines runs after _reorder_fid_lines (the ky sort). For a multi-shot EPI with an even NSegments, PVM_EncSteps1 of PV7 exp 13 is -48 -40 -32 ... 40 | -47 -39 -31 ...: shot s acquires ky = s + 8*j. So ky parity equals shot parity, not echo-in-train parity, and the current code mirrors the whole odd shots. Only single-shot EPI gets the correct lines today. Applying the reversal before the sort (the §10.1 order) removes this problem at no extra cost.

Options

A. Ignore the declaration for non-EPI schemes (the patch in the report). Rejected. It leaves the odd echoes mirrored (table above) and it contradicts §6.2.

B. Apply §6.2 in acquisition order, before the ky sort. Two rules, keyed on ACQ_scan_size (§5.3):

• One_scan (MGE, any Cartesian multi-echo): reverse the readout of every second scan of the decoded stream. Verified on PV7 exp 26: a two-line patch in SchemaFid.deserialize after _decode_raw_stream (data[:, 1::2] = data[::-1, 1::2]) gives exactly the odd-echo reversal (np.array_equal is True). The physically meaningful form of the same rule is "echo index object_id % ACQ_n_echo_images is odd". The two forms agree on all corpus data; they differ only if NI is odd (odd echo count and odd slice count), and I have no such data.
• ACQ_phase_factor_scans (EPI): reverse the odd lines within each train, in acquisition order, before _reorder_fid_lines. This also repairs the multi-shot parity.
  This is my recommendation. It follows the spec order, one code path per scan type, and the existing mirror_odd_lines property keeps its role as the switch.

C. Port RecoReverseFilter semantics. Use RecoReverseSegment (PV6+, PV360) as the parity restart length (0 = whole stream, N = restart every N lines) and fall back to the rules in B when it is absent (PV5.1 FT_MODE recos have no RecoStage*). Most faithful to ParaVision. More code, and RecoReverseSegment is not documented yet in FILE_FORMAT.md (only the RecoSort* group is, in §10.2). I will add it.

D. Remove the reversal from the reader. Drop _mirror_odd_lines, keep RECO_inp_order as metadata, and document that Dataset.kspace keeps the acquisition polarity, as SchemaRawdata and ramp regridding already do. Smallest code, consistent between fid and rawdata. But it is a behaviour change for EPI users, and the result is then not a k-space.

E. Keep the current EPI mirror and add an echo-axis flip when ACQ_n_echo_images > 1 (@headmeister's minimal proposal). Smallest diff that fixes MGE. It keeps the multi-shot EPI parity problem and two code paths.

Whichever option we pick: the synthetic fixture in test/test_raw_layouts.py (FLASH with REV_ALT_ROWS) asserts the mirror of odd ky columns, so it encodes the wrong axis and must change. Good regression data are PV6.0.1 exp 6 and PV7 exp 26 (odd echoes against 2dseq) and PV7 exp 13 (8-shot EPI).

On the document side I will update FILE_FORMAT.md in pvraw: §6.2 to define "scan" with the MGE and EPI examples and the EchoAcqMode <-> RECO_inp_order table above, §10.2 to add RecoReverseSegment, and §5.2 to add the MGE ACQ_obj_order example.

@schuenke a PR along option B is welcome. If you prefer, I can do it.

[schuenke]

Hey @headmeister and @gdevenyi,
First of all, thanks for your work on this project!

So, if we already have this phase flipping functionality, then we would likely want this to work as I described, rather than omit phase flips altogether. Therefore, we would likely rather adjust the flipping in the _mirror_odd_lines, where instead of blindly flipping every ky line, we need to flip every other echo...

You are absolutely right. My code from above only removes the artefacts in the individual images, but does not align the images of different echos. Fixing this would actually require adjusting the _mirror_odd_lines() itself and not just the property.

A simple fix that seems to work for me would be the following, but this is probably more like a hot fix instead of a general approach

def _mirror_odd_lines(self, data):
    if "EPI" in self._dataset.scheme_id:
        # reverse readout of odd phase-encode lines
        for index in np.ndindex(data.shape[2:]):
            idx = (slice(None), slice(1, data.shape[1], 2)) + tuple(index)
            data[idx] = data[idx][::-1, :]
    else:
        # reverse readout of odd "objects"
        for index in np.ndindex(data.shape[3:]):
            idx = (slice(None), slice(None), slice(1, data.shape[2], 2)) + tuple(index)
            data[idx] = data[idx][::-1, :, :]
    return data

I still have not entirely decided how far we should go with BrukerAPI with regard to RECO. While we likely reorder EPI data correctly into k-space, we do not perform any ghost correction or regridding based on the measured trajectory. Therefore, anyone trying to use BrukerAPI's k-space output directly, especially for EPI, will see degraded images anyway.

I think this is completely fine. A ghost correction is part of the RECON and IMO at least the Dataset.kspace attribute should not go beyond reordering. However, since the kspace array does not include any header infos about the individual acquisitions / k-space lines, the flipping of reverse data should be included. To be consistent, it should not only be done for the EPI, but also for all other bruker methods like the MGE mentioned above.

As a general information why I came across this project and the MGE issue in general:
I started working on a new version of a bruker_to_ismrmrd converter, which, thanks to your latest efforts, is actually quite straight forward to realize. It's basically already working for all PV versions and methods that you support in brukerapi and even a little beyond: non-cartesian PV360 methods do work as well in the converter. So I think only PV360 EPI is missing at the moment.

Before starting the general bruker_to_ismrmrd converter, we implemented a converter for Bruker Pulseq (the interpreter / driver is still work in progress) data, which uses the Pulseq file to create an ISMRMRD skeleton and fills the missing acquisition data using Dataset.raw. This is also working fine already.

@schuenke a PR along option B is welcome. If you prefer, I can do it.

I will have a look at your. Your reply just came in while I was typing.

[schuenke]

Hey again,
I had a look at the suggested option B and the following SchemaFid.deserialize method seems to fix the PV7 MGE and PV7 multi-shot EPI issue:

def deserialize(self, data, layouts):
    data = self._decode_raw_stream(data, layouts)

    if self.mirror_odd_lines and not self.continuous_train:
        # One_scan (MGE): reverse every second scan in raw stream
        data[:, 1::2] = data[::-1, 1::2]

    data = self._acquisitions_to_encode(data, layouts)
    data = self._encode_to_permute(data, layouts)
    data = self._permute_to_kspace(data, layouts)

    if self.mirror_odd_lines and self.continuous_train:
        # EPI: reverse odd phase lines in acq order, BEFORE ky sort
        self._mirror_odd_lines(data)

    data = self._reorder_fid_lines(data, dir="FW")
    data = self._reorder_objects(data, dir="FW")
    return data

However, the issue in the PV360 bipolar MGE is not fixed by this, because SchemaRawdata currently used there does not use a _mirror_odd_lines method.

Should I create a PR with the fix above and we tackle the PV360 issue separately or do you want to include the PV360 fix in the same PR?

A "hacky" fix would be to apply something like

if str(self._dataset._parameter_value("EchoAcqMode", "")).lower() == "allechoes":
    data[:, :, 1::2] = data[::-1, :, 1::2]  # reverse readout of odd echoes

in SchemaRawdata, after reshaping to kspace, but this is not the idea of SchemaRawdata, is it?

[gdevenyi]

Should I create a PR with the fix above and we tackle the PV360 issue separately or do you want to include the PV360 fix in the same PR?

Seperately for sure please.

[gdevenyi]

@schuenke thanks, I ran your deserialize against the corpus. It is the correct shape. Three things to change before the PR, plus two facts from the data.

1. The One_scan flip must be per receiver. After _decode_raw_stream the first axis is (sample, receiver) in Fortran order (see SchemaFid.raw, which reshapes it to (samples, receivers, scans) with order="F"). data[::-1, 1::2] therefore also reverses the receiver order when PVM_EncNReceivers > 1. This was also wrong in my option B text. Use:

receivers = int(self._dataset.channels)
scans = data.reshape((-1, receivers, data.shape[1]), order="F")   # (sample, receiver, scan)
scans[:, :, 1::2] = scans[::-1, :, 1::2]
data = scans.reshape((-1, data.shape[1]), order="F")

On PV7 exp 26 this is bit-identical to your version (one receiver). No PV5.1/6/7 fid in the corpus has more than one receiver, so the PR needs a synthetic two-receiver fixture; the helpers in test/test_raw_layouts.py can build one. All PV360 EPI/DTI data in the corpus have 4 receivers, so the same form is mandatory for the later PV360 fix.

2. serialize must change with it. The write path still applies the old odd-ky mirror after the ky sort. Round trip stored -> deserialize -> serialize on PV7 exp 26: current code 0 mismatches; your patch 654,272 of 1,310,720 stored samples mismatch (the odd-echo scans). Move the reversal to the same positions in serialize (after _reorder_fid_lines(dir="BW") for the train case, after the reshape to the stream for the One_scan case). There is no fid round-trip test today, only test_trajectory.py::test_traj_deserialize_serialize_preserves_stored_values; please add one for fid on exp 26 and one EPI.

3. The EPI branch: parity restarts per train. Every EPI-family reco in the corpus (23 of 23, PV6.0.1, PV7, PV360) has RecoReverseSegment == PVM_EpiNEchoes; every MGE reco has 0. So ParaVision's RecoReverseFilter restarts the alternation at each echo train, and MGE alternates over the whole stream. Your _mirror_odd_lines along axis 1 in acquisition order gives the same result only when the lines per train are even. All corpus trains are even (12 to 80 lines), so it passes today. Either keep it with a comment, or use (line % PVM_EpiNEchoes) % 2 (fallback PVM_EncMatrix[1] // NSegments) and it is correct in general. I will document RecoReverseSegment in FILE_FORMAT.md §10.2 with this rule.

Two facts for the regression data.

• PV7 exp 13 and exp 10 are double-sampled EPI: PVM_EpiCombine=Yes and PVM_EpiNEchoes = 2 * PVM_EncMatrix[1] / NSegments (24 echoes, 12 lines). The reco graph splits odd and even echoes, sorts each (R10, R20), sums and divides by 2. They are not clean tests for the multi-shot branch, and whether the EPI layout handles them is a separate question. Clean multi-shot data: PV6.0.1 exp 20 (T1_EPI, 8 shots), 17 (FAIR_EPI, 4), 16 (EPI, 3); PV7 exp 16 (FAIR_EPI, 2). Clean single-shot: PV7 exp 14, 39, 41 and PV6.0.1 exp 45.
• The test_raw_layouts.py fixture (FLASH, ACQ_scan_size=ACQ_phase_factor_scans, REV_ALT_ROWS) asserts the odd-ky mirror; replace its expected array with the in-train form.

Random access is not affected: mmap=True is rejected for fid, so SchemaFid.ra is not reachable.

For PV360 later: keep RECO_inp_order as the switch there too. EchoAcqMode is an MGE-local parameter name, and the reco declaration is the generic signal (table in my previous comment).