Getting Started¶

Import precon:

import precon as pr
from pathlib import Path

Reading Parameters¶

The Parameter class reads all scan parameters and scanner labels from a Philips raw file. Parameters are automatically sourced from the embedded ReconFrame patch data, or from a companion .sin file if no patch is present.

pars = pr.Parameter(Path('my_rawfile.raw'))

Access a commonly used parameter directly:

voxel_sizes = pars.get_voxel_sizes(mix=0)
recon_resolution = pars.get_recon_resolution(mix=0, xovs=False, yovs=True, zovs=True)

Retrieve any scanner GOAL parameter by name:

flip_angle = pars.get_value('EX_ACQ_flip_angle', default=[90.0])

get_value always returns a list. If the parameter is absent the default is returned; if no default is given, None is returned.

Access a GOAL object attribute directly:

sq_base = pars.goal.get_object('SQ', 'base')
sq_base_dur = sq_base.get_attribute('dur', cmp=0)

Reading and Sorting Data¶

Use Parameter2Read to specify which data to read. After reading, the data is in temporal order — use sort() to assemble it into a proper k-space.

parameter2read = pr.Parameter2Read(pars.labels)
for mix in parameter2read.mix:
    for stack in parameter2read.stack:
        parameter2read.stack = stack
        parameter2read.mix = mix

        with open('my_rawfile.raw', 'rb') as raw:
            data, labels = pr.read(raw, parameter2read, pars.labels, pars.coil_info)

        cur_recon_resolution = pars.get_recon_resolution(mix=mix, xovs=False, yovs=True, zovs=True)
        data, labels = pr.sort(data, labels, output_size=cur_recon_resolution)

The read function always returns data in the shape samples × coils × profiles. sort reorganises it into a 12-dimensional k-space: (kx, ky, kz, channels, dynamics, cardiac phases, echoes, locations, mixes, extra1, extra2, averages).

Note

read only handles data of a single size and geometry. Always loop over mixes and stacks.

Basic Reconstruction¶

A minimal Cartesian reconstruction after reading and sorting:

# Fourier transform kx, ky, kz → image space
data = pr.k2i(data, axis=(0, 1, 2))

# Sum-of-squares coil combination
data = pr.sos(data, axis=3)

Common Patterns¶

Partial Fourier (Homodyne)¶

kx_range = pars.get_range(enc=0, mix=mix, stack=stack, ovs=False)
ky_range = pars.get_range(enc=1, mix=mix, stack=stack)
kz_range = pars.get_range(enc=2, mix=mix, stack=stack)

if pr.is_partial_fourier(kx_range) or pr.is_partial_fourier(ky_range) or pr.is_partial_fourier(kz_range):
    data = pr.homodyne(data, kx_range, ky_range, kz_range)

data = pr.k2i(data, axis=(0, 1, 2))

SENSE Reconstruction¶

sens = pr.calculate_sensitivities(Path(refscan), Path(rawfile), match_target_size=False)
regularization_factor = pars.get_value(pars.SENSE_REGULARIZATION_FACTOR, at=0, default=2)
output_size = pars.get_recon_resolution(mix=mix, xovs=False, yovs=True, zovs=True, folded=False)
data = pr.sense_unfold(data, sens, output_size, regularization_factor=regularization_factor)

Array Compression¶

Array compression reduces the number of coil channels before reconstruction, cutting memory and compute cost:

ref_pars = pr.Parameter(Path(refscan))
qbc, coil = pr.reconstruct_refscan(ref_pars)
sens = pr.reformat_refscan(qbc, coil, ref_pars, pars, match_target_size=True)

A = pr.get_array_compression_matrix(sens.surfacecoil)
virtual_coils = 8
A = A[0:virtual_coils, :]

# Compress the sensitivity maps for SENSE
sens = pr.compress_sensitivity(sens, A)

# Pass A to read() so compression is applied on-the-fly during reading
with open(rawfile, 'rb') as raw:
    data, labels = pr.read(raw, parameter2read, pars.labels, pars.coil_info, array_compression=A)

Flow Reconstruction¶

data = pr.divide_flow_segments(data, pars.is_hadamard_encoding())
data = pr.format_flow(data, pars.get_coordinate_system(), pars.get_venc(), pars.is_hadamard_encoding())
data = pr.fit_flow_phase(data, order=3)

Radiological Formatting¶

Transform images into radiological convention as a final step:

transforms = pars.get_in_plane_transformation(mix=0, stack=0)
data = pr.format(data, transforms)

Non-Cartesian Reconstruction¶

See the Reference page for the full API. A radial example:

no_spokes = labels_sorted[0, 0, ...].shape[1]  # number of acquired spokes
traj = pr.noncart.get_radial_trajectory(no_samples, no_spokes)
weights = pr.noncart.get_radial_weights(traj)
images = pr.noncart.nufft(traj, data, weights, out_res=[matrix_size, matrix_size, 1])