From 73c989303475da758f27990b58d1e0397864c252 Mon Sep 17 00:00:00 2001 From: Matteo Cencini Date: Wed, 30 Oct 2024 18:36:08 +0100 Subject: [PATCH] update wrapper for SPGR --- src/pulserver/_apps/SPGR2D.py | 203 +-------------------- src/pulserver/_apps/SPGR3D.py | 184 ------------------- src/pulserver/_server/_server.py | 4 + src/pulserver/parsing/_cartesian_params.py | 26 ++- 4 files changed, 36 insertions(+), 381 deletions(-) delete mode 100644 src/pulserver/_apps/SPGR3D.py diff --git a/src/pulserver/_apps/SPGR2D.py b/src/pulserver/_apps/SPGR2D.py index 0fe59b1..003c9ec 100644 --- a/src/pulserver/_apps/SPGR2D.py +++ b/src/pulserver/_apps/SPGR2D.py @@ -2,205 +2,22 @@ __all__ = ["SPGR2D"] -from collections.abc import Iterable +from pulserver.parsing import Cartesian2DParams +from pulserver.sequences import design_2D_spgr -import numpy as np -import pypulseq as pp - - -from pulserver import Sequence -from pulserver.plan import RfPhaseCycle - - -def SPGR2D( - fov: Iterable[float], - slice_thickness: float, - npix: Iterable[int], - nslices: int, - alpha: float, - rectime: float, - max_grad: float, - max_slew: float, - raster_time: float, - slice_spacing: float = 0.0, - seqformat: str | bool = "bytes", -): +def SPGR2D(kwargs): """ Generate a 2D Spoiled Gradient Recalled Echo (SPGR) pulse sequence. - This function designs a 2D SPGR sequence based on the provided field of view (FOV), matrix size, - number of slices, slice thickness and spacing, flip angle, recovery time, - and hardware constraints such as maximum gradient amplitude and slew rate. - The output can be formatted in different sequence file formats if specified. - - Parameters - ---------- - fov : Iterable[float] - Field of view along each spatial dimension [fov_x, fov_y] in mm. - If scalar, assume squared fov. - slice_thickness : float) - Slice thickness in mm. - npix : Iterable[int] - Number of voxels along each spatial dimension [nx, ny] (matrix size). - If scalar, assume squared matrix size. - nslices : int - Number of slices. - alpha : float - Flip angle in degrees. - rectime : float - Recovery time after spoiling in ms. - max_grad : float - Maximum gradient amplitude in mT/m. - max_slew : float - Maximum gradient slew rate in T/m/s. - raster_time : float - Waveform raster time in seconds (the time between successive gradient samples). - slice_spacing : float, optional - Additional slice spacing in mm. The default is 0.0 (contiguous slices). - seqformat : str or bool, optional - Output sequence format. If a string is provided, it specifies the desired output format (e.g., 'pulseq', 'bytes'). - If False, the sequence is returned as an internal object. Default is False. - - Returns - ------- - seq : object or dict - The generated SPGR sequence. If `seqformat` is a string, the sequence is returned in the specified format. - If `seqformat` is False, the sequence is returned as an internal representation. - - Notes - ----- - - This function is designed to work within the constraints of MRI scanners, taking into account the physical limits - on gradient amplitude and slew rates. - - The flip angle (`alpha`) controls the excitation of spins and directly impacts the signal-to-noise ratio (SNR) and contrast. - - Examples - -------- - Generate a 2D SPGR sequence for a single 5 mm thick slice and 240x240 mm FOV, 256x256 matrix size, - 15-degree flip angle, 5s recovery time and hardware limits 4mT/m, 150T/m/s, 4e-6 s raster time as: - - >>> from pulseforge import SPGR2D - >>> SPGR2D(240.0, 5.0, 256, 1, 15.0, 5.0, 40, 150, 4e-6) - - Generate the same sequence and export it in bytes format: - - >>> SPGR2D(240.0, 5.0, 256, 1, 15.0, 5.0, 40, 150, 4e-6, format='bytes') + This function wraps ``pulserver.sequences.design_2D_spgr`` - see + the corresponding docstrings for more information. """ - # RF specs - rf_spoiling_inc = 117.0 # RF spoiling increment - - # initialize system limits - system = pp.Opts( - max_grad=max_grad, - grad_unit="mT/m", - max_slew=max_slew, - slew_unit="T/m/s", - grad_raster_time=raster_time, - rf_raster_time=raster_time, - ) - - # initialize sequence - seq = Sequence(system=system, format=seqformat) - - # initialize prescription - fov, slice_thickness, slice_spacing = ( - fov * 1e-3, - slice_thickness * 1e-3, - slice_spacing * 1e-3, - ) - slice_spacing += slice_thickness - - if np.isscalar(npix): - Nx, Ny, Nz = npix, npix, nslices # in-plane resolution, slice thickness - else: - Nx, Ny, Nz = npix[0], npix[1], nslices # in-plane resolution, slice thickness - - # initialize event events - # RF pulse - rf, gss, _ = pp.make_sinc_pulse( - flip_angle=np.deg2rad(alpha), - duration=3e-3, - slice_thickness=slice_thickness, - apodization=0.42, - time_bw_product=4, - system=system, - return_gz=True, - ) - gss_reph = pp.make_trapezoid( - channel="z", area=-gss.area / 2, duration=1e-3, system=system - ) - seq.register_block(name="excitation", rf=rf, gz=gss) - seq.register_block(name="slab_rephasing", gz=gss_reph) - - # readout - delta_kx, delta_ky = 1 / fov, 1 / fov - g_read = pp.make_trapezoid( - channel="x", flat_area=Nx * delta_kx, flat_time=3.2e-3, system=system - ) - adc = pp.make_adc( - num_samples=Nx, duration=g_read.flat_time, delay=g_read.rise_time, system=system - ) - seq.register_block("dummy_readout", gx=g_read) - seq.register_block("readout", gx=g_read, adc=adc) - - # phase encoding - gx_phase = pp.make_trapezoid( - channel="x", area=-g_read.area / 2, duration=1e-3, system=system - ) - gy_phase = pp.make_trapezoid(channel="y", area=delta_ky * Ny, system=system) - seq.register_block("g_phase", gx=gx_phase, gy=gy_phase) - - # crusher gradient - gz_spoil = pp.make_trapezoid(channel="z", area=32 / slice_thickness, system=system) - delay = pp.make_delay(rectime) - seq.register_block("g_spoil", gz=gz_spoil, delay=delay) - - # phase encoding plan TODO: helper routine - pey_steps = ((np.arange(Ny)) - (Ny / 2)) / Ny - FOVz = Nz * slice_spacing - foff_steps = np.linspace(-FOVz / 2, FOVz / 2, Nz) * gss.amplitude - encoding_plan = np.meshgrid(pey_steps, foff_steps, indexing="xy") - encoding_plan = [enc.ravel() for enc in encoding_plan] - - # scan duration - dummy_scans = 10 - calib_scans = 10 - imaging_scans = Ny * nslices - - # generate rf phases - rf_phases = RfPhaseCycle( - num_pulses=dummy_scans + imaging_scans + calib_scans, - phase_increment=rf_spoiling_inc, - ) - - # construct sequence - seq.section(name="ss_prep") - for n in range(dummy_scans): - rf_phase = rf_phases() - seq.add_block("excitation") - seq.add_block("slab_rephasing") - seq.add_block("g_phase", gy_amp=0.0) - seq.add_block("dummy_readout", rf_phase=rf_phase) - seq.add_block("g_phase", gy_amp=0.0) - seq.add_block("g_spoil") - - seq.section(name="scan_loop") - for n in range(imaging_scans + calib_scans): - rf_phase = rf_phases() - - if n < calib_scans: - rf_freq = 0.0 - y_amp = 0.0 - else: - rf_freq = encoding_plan[1][n - calib_scans] - y_amp = encoding_plan[0][n - calib_scans] + # parse parameters + params = Cartesian2DParams(**kwargs, fudge_factor=0.9) - seq.add_block("excitation", rf_freq=rf_freq) - seq.add_block("slab_rephasing") - seq.add_block("g_phase", gy_amp=y_amp) - seq.add_block("readout", rf_phase=rf_phase, adc_phase=rf_phase) - seq.add_block("g_phase", gy_amp=-y_amp) - seq.add_block("g_spoil") + # call design function (exclude header for now) + seq, _ = design_2D_spgr(**params.asdict(), platform="gehc") - return seq.export() + return seq.export(), None diff --git a/src/pulserver/_apps/SPGR3D.py b/src/pulserver/_apps/SPGR3D.py deleted file mode 100644 index 5afc6dd..0000000 --- a/src/pulserver/_apps/SPGR3D.py +++ /dev/null @@ -1,184 +0,0 @@ -"""3D Spoiled Gradient Echo sequence.""" - -__all__ = ["SPGR3D"] - -from collections.abc import Iterable - - -import numpy as np -import pypulseq as pp - - -from pulserver import Sequence -from pulserver.plan import RfPhaseCycle - - -def SPGR3D( - fov: Iterable[float], - npix: Iterable[int], - alpha: float, - max_grad: float, - max_slew: float, - raster_time: float, - seqformat: str | bool = "bytes", -): - """ - Generate a 3D Spoiled Gradient Recalled Echo (SPGR) pulse sequence. - - This function designs a 3D SPGR sequence based on the provided field of view (FOV), matrix size, - flip angle, and hardware constraints such as maximum gradient amplitude and slew rate. The output - can be formatted in different sequence file formats if specified. - - Parameters - ---------- - fov : Iterable[float] - Field of view along each spatial dimension [fov_plane, fov_z] in mm. - If scalar, assume cubic fov. - npix : Iterable[int] - Number of voxels along each spatial dimension [plane_mtx, nz] (matrix size). - If scalar, assume cubic matrix size. - alpha : float - Flip angle in degrees. - max_grad : float - Maximum gradient amplitude in mT/m. - max_slew : float - Maximum gradient slew rate in T/m/s. - raster_time : float - Waveform raster time in seconds (the time between successive gradient samples). - seqformat : str or bool, optional - Output sequence format. If a string is provided, it specifies the desired output format (e.g., 'pulseq', 'bytes'). - If False, the sequence is returned as an internal object. Default is False. - - Returns - ------- - seq : object or dict - The generated SPGR sequence. If `seqformat` is a string, the sequence is returned in the specified format. - If `seqformat` is False, the sequence is returned as an internal representation. - - Notes - ----- - - This function is designed to work within the constraints of MRI scanners, taking into account the physical limits - on gradient amplitude and slew rates. - - The flip angle (`alpha`) controls the excitation of spins and directly impacts the signal-to-noise ratio (SNR) and contrast. - - Examples - -------- - Generate a 3D SPGR sequence for a 256x256x128 matrix with a 240x240x120 mm FOV, a 15-degree flip angle, and hardware limits: - - >>> from pulseforge import SPGR3D - >>> SPGR3D([240, 120], [256, 128], 15, 0.04, 150, 4e-6) - - Generate the same sequence and export it in bytes format: - - >>> SPGR3D([240, 120], [256, 128], 15, 0.04, 150, 4e-6, seqformat='bytes') - - """ - # RF specs - rf_spoiling_inc = 117.0 # RF spoiling increment - - # initialize system limits - system = pp.Opts( - max_grad=max_grad, - grad_unit="mT/m", - max_slew=max_slew, - slew_unit="T/m/s", - grad_raster_time=raster_time, - rf_raster_time=raster_time, - ) - - # initialize sequence - seq = Sequence(system=system, format=seqformat) - - # initialize prescription - if np.isscalar(fov): - fov, slab_thickness = fov * 1e-3, fov * 1e-3 # isotropic - else: - fov, slab_thickness = ( - fov[0] * 1e-3, - fov[1] * 1e-3, - ) # in-plane FOV, slab thickness - - if np.isscalar(npix): - Nx, Ny, Nz = npix, npix, npix # in-plane resolution, slice thickness - else: - Nx, Ny, Nz = npix[0], npix[0], npix[1] # in-plane resolution, slice thickness - - # initialize event events - # RF pulse - rf, gss, _ = pp.make_sinc_pulse( - flip_angle=np.deg2rad(alpha), - duration=3e-3, - slice_thickness=slab_thickness, - apodization=0.42, - time_bw_product=4, - system=system, - return_gz=True, - ) - gss_reph = pp.make_trapezoid( - channel="z", area=-gss.area / 2, duration=1e-3, system=system - ) - seq.register_block(name="excitation", rf=rf, gz=gss) - seq.register_block(name="slab_rephasing", gz=gss_reph) - - # readout - delta_kx, delta_ky, delta_kz = 1 / fov, 1 / fov, 1 / slab_thickness - g_read = pp.make_trapezoid( - channel="x", flat_area=Nx * delta_kx, flat_time=3.2e-3, system=system - ) - adc = pp.make_adc( - num_samples=Nx, duration=g_read.flat_time, delay=g_read.rise_time, system=system - ) - seq.register_block("dummy_readout", gx=g_read) - seq.register_block("readout", gx=g_read, adc=adc) - - # phase encoding - gx_phase = pp.make_trapezoid( - channel="x", area=-g_read.area / 2, duration=1e-3, system=system - ) - gy_phase = pp.make_trapezoid(channel="y", area=delta_ky * Ny, system=system) - gz_phase = pp.make_trapezoid(channel="z", area=delta_kz * Nz, system=system) - seq.register_block("g_phase", gx=gx_phase, gy=gy_phase, gz=gz_phase) - - # crusher gradient - gz_spoil = pp.make_trapezoid(channel="z", area=32 / slab_thickness, system=system) - seq.register_block("g_spoil", gz=gz_spoil) - - # phase encoding plan TODO: helper routine - pey_steps = ((np.arange(Ny)) - (Ny / 2)) / Ny - pez_steps = ((np.arange(Nz)) - (Nz / 2)) / Nz - encoding_plan = np.meshgrid(pey_steps, pez_steps, indexing="xy") - encoding_plan = [enc.ravel() for enc in encoding_plan] - - # scan duration - dummy_scans = Ny - imaging_scans = Ny * Nz - - # generate rf phases - rf_phases = RfPhaseCycle( - num_pulses=dummy_scans + imaging_scans, phase_increment=rf_spoiling_inc - ) - - # construct sequence - seq.section(name="ss_prep") - for n in range(dummy_scans): - rf_phase = rf_phases() - seq.add_block("excitation") - seq.add_block("slab_rephasing") - seq.add_block("g_phase", gy_amp=0.0, gz_amp=0.0) - seq.add_block("dummy_readout", rf_phase=rf_phase) - seq.add_block("g_phase", gy_amp=0.0, gz_amp=0.0) - seq.add_block("g_spoil") - - seq.section(name="scan_loop") - for n in range(imaging_scans): - rf_phase = rf_phases() - seq.add_block("excitation") - seq.add_block("slab_rephasing") - seq.add_block("g_phase", gy_amp=encoding_plan[0][n], gz_amp=encoding_plan[1][n]) - seq.add_block("readout", rf_phase=rf_phase, adc_phase=rf_phase) - seq.add_block( - "g_phase", gy_amp=-encoding_plan[0][n], gz_amp=-encoding_plan[1][n] - ) - seq.add_block("g_spoil") - - return seq.export() diff --git a/src/pulserver/_server/_server.py b/src/pulserver/_server/_server.py index b195bd5..74cec45 100644 --- a/src/pulserver/_server/_server.py +++ b/src/pulserver/_server/_server.py @@ -222,6 +222,7 @@ def send_to_recon_server(optional_buffer, config): with socket.socket(socket.AF_INET, socket.SOCK_STREAM) as s: s.connect((RECON_SERVER_ADDRESS, RECON_SERVER_PORT)) s.sendall(optional_buffer) + s.shutdown(socket.SHUT_WR) def handle_client_connection(config, client_socket, plugins, logger): @@ -244,6 +245,9 @@ def handle_client_connection(config, client_socket, plugins, logger): # Send the result buffer to the client client_socket.sendall(result_buffer) + # Signal that no more data will be sent + client_socket.shutdown(socket.SHUT_WR) + # Optionally send the reconstruction info to the secondary server if optional_buffer is not None: send_to_recon_server(optional_buffer, config) diff --git a/src/pulserver/parsing/_cartesian_params.py b/src/pulserver/parsing/_cartesian_params.py index 98710fa..4070ea4 100644 --- a/src/pulserver/parsing/_cartesian_params.py +++ b/src/pulserver/parsing/_cartesian_params.py @@ -37,6 +37,7 @@ def __init__( rf_ringdown_time: float | None = 60e-6, adc_dead_time: float | None = 40e-6, psd_rf_wait: float | None = 0.0, + fudge_factor: float | None = None, *args, **kwargs, ): # noqa @@ -73,13 +74,21 @@ def __init__( self.flip_angle = flip # TE / TR - self.TE = TE * 1e-3 - self.TR = TR * 1e-3 + if self.TE is not None: + self.TE = TE * 1e-3 + if self.TR is not None: + self.TR = TR * 1e-3 # Accelerations self.R = R self.PF = PF + # apply fudge + if fudge_factor is not None and gmax is not None: + gmax = fudge_factor * gmax + if fudge_factor is not None and smax is not None: + smax = fudge_factor * smax + # Build opts super().__init__( dwell, @@ -128,6 +137,7 @@ def __init__( rf_ringdown_time: float | None = 60e-6, adc_dead_time: float | None = 40e-6, psd_rf_wait: float | None = 0.0, + fudge_factor: float | None = None, *args, **kwargs, ): # noqa @@ -161,8 +171,10 @@ def __init__( self.flip_angle = flip # TE / TR - self.TE = TE * 1e-3 - self.TR = TR * 1e-3 + if self.TE is not None: + self.TE = TE * 1e-3 + if self.TR is not None: + self.TR = TR * 1e-3 # Accelerations self.R = R @@ -172,6 +184,12 @@ def __init__( self.Rshift = Rshift self.PF = PF + # apply fudge + if fudge_factor is not None and gmax is not None: + gmax = fudge_factor * gmax + if fudge_factor is not None and smax is not None: + smax = fudge_factor * smax + # Build opts super().__init__( dwell,