Contrib

Experimental

The contrib module contains experimental sequence definitions that are not covered by the same stability expectations as the core API.

contrib

Module containing experimental contributions, not tested!

gen_4DFlow_sequence

gen_4DFlow_sequence(
    system_specs: SystemSpec,
    LUT: ndarray,
    prof_per_phase: int,
    matrix: ndarray,
    resolution: ndarray,
    venc_list: tuple,
    venc_dir: tuple,
    spoil_moments: tuple,
    slice_thickness: Quantity,
    adc_duration: Quantity,
    flip_angle: Quantity,
    pulse_duration: Quantity,
    slice_position_offset: Quantity = Quantity(0.0, "m"),
    time_bandwidth_product: float = 4.0,
    venc_duration: Quantity = Quantity(0.0, "ms"),
    rf_spoiling: bool = False,
    balanced: bool = False,
    bal_norewind: bool = False,
    rampup_shots: int = 0,
)

Source code in cmrseq/contrib/_4DFlow.py

def gen_4DFlow_sequence(system_specs: cmrseq.SystemSpec,
                        LUT: np.ndarray,
                        prof_per_phase: int,
                        matrix: np.ndarray,
                        resolution: np.ndarray,
                        venc_list:tuple,
                        venc_dir:tuple,
                        spoil_moments: tuple,
                        slice_thickness: Quantity,
                        adc_duration: Quantity,
                        flip_angle: Quantity,
                        pulse_duration: Quantity,
                        slice_position_offset: Quantity = Quantity(0., "m"),
                        time_bandwidth_product: float = 4.,
                        venc_duration:Quantity = Quantity(0.,'ms'),
                        rf_spoiling:bool=False,
                        balanced:bool=False,
                        bal_norewind:bool=False,
                        rampup_shots:int=0):

    # kspace extent
    k_ext = 1 / resolution.to('m')

    num_samples = matrix[0]

    # kspace step size
    dk = k_ext / matrix

    # --- Generate reference sequence blocks ---

    if balanced:
        spoil_moments = (0,0,0)
    # First calculate spoiling area
    spoiler_area = (Quantity(spoil_moments,'rad')/system_specs.gamma_rad/resolution).to('mT/m*ms')


    # Get reference sequence that will have P and S prewind/rewind gradients scaled later
    fastseq = spoiled_3D_cartesian_line(system_specs=system_specs,
                                         num_samples=num_samples,
                                         k_M_total=k_ext[0],
                                         k_P=k_ext[1] / 2,
                                         k_S=k_ext[2] / 2,
                                         adc_duration=adc_duration,
                                         spoiler_area=spoiler_area)
    # Scale P and S to unit area
    fastseq.get_block('prephaser_P_0').scale_gradients(1/fastseq.get_block('prephaser_P_0').area[1].m_as('mT/m*ms'))
    fastseq.get_block('prephaser_S_0').scale_gradients(1/fastseq.get_block('prephaser_S_0').area[2].m_as('mT/m*ms'))

    fastseq.get_block('prephaser_P_rewind_0').scale_gradients(1/fastseq.get_block('prephaser_P_rewind_0').area[1].m_as('mT/m*ms'))
    fastseq.get_block('prephaser_S_rewind_0').scale_gradients(1/fastseq.get_block('prephaser_S_rewind_0').area[2].m_as('mT/m*ms'))

    if slice_thickness is not None:
        # generate RF pulse
        rf_seq = cmrseq.seqdefs.excitation.slice_selective_sinc_pulse(system_specs=system_specs,
                                                                      slice_thickness=slice_thickness,
                                                                      flip_angle=flip_angle,
                                                                      pulse_duration=pulse_duration,
                                                                      time_bandwidth_product=time_bandwidth_product,
                                                                      slice_position_offset=slice_position_offset,
                                                                      slice_normal=np.array([0., 0., 1.]))

        #Add prewinder to RF
        rf_prewind = deepcopy(rf_seq.get_block('slice_select_rewind_0'))
        rf_prewind.shift(-rf_prewind.tmin)

        rf_seq.shift_in_time(rf_prewind.duration)
        rf_seq.add_block(rf_prewind)

    else:
        rf_pulse = cmrseq.bausteine.HardRFPulse(system_specs=system_specs,
                                                flip_angle=flip_angle,
                                                duration=pulse_duration,
                                                name='rf_excitation')
        rf_seq = cmrseq.Sequence([rf_pulse],system_specs=system_specs)

    # generate bipolar gradient
    # First generate fastest bipolar for the strongest case (minimum venc)
    vms = np.array([_.m_as('m/s') for _ in venc_list])
    vms = vms[np.nonzero(vms)]
    # If there are no non-zero values, set min venc to zero
    if len(vms)==0:
        venc_min = Quantity(0., 'm/s')
    else: # Otherwise we calculate the fastest gradient
        venc_min = Quantity(np.min(vms), 'm/s')
        bip_fastest = cmrseq.parametric_definitions.velocity.bipolar(system_specs=system_specs,
                                                                     venc=venc_min,
                                                                     direction=np.array([0.,0.,1.]))
        if bip_fastest.duration > venc_duration:
            venc_duration = bip_fastest.duration


    bipolars = []
    for venc,dir in zip(venc_list,venc_dir):
        if venc_duration>0:
            bip = cmrseq.parametric_definitions.velocity.bipolar(system_specs=system_specs,
                                                                venc=venc,
                                                                duration=venc_duration,
                                                                direction=dir)
            bipolars.append(bip)


    # Scaling area by this factor results in a traverse of dk
    dk_area_scale = (dk/system_specs.gamma).m_as('mT/m*ms')

    # number of expected blocks of all segs
    expected_blocks = np.ceil(LUT.shape[1] / prof_per_phase).astype(int)

    seq_list = []

    # Array containing empty lists
    LUT_return = np.zeros([4, LUT.shape[1]*LUT.shape[2]*LUT.shape[3] + rampup_shots])
    trcount = 0
    lutcount = 0
    rf_phase_init = Quantity(117,'degree').to('rad')

    # Ramp up shots
    rampup_list = []
    pbar_ramp = tqdm(range(rampup_shots), desc="Loop - Rampup",leave=False)
    #Extract P and S from first TR of sequence
    kP = LUT[0, 0, 0, 0]
    kS = LUT[1, 0, 0, 0]
    for tr in pbar_ramp:
         # Store some info for sorting during recon

        readout = deepcopy(fastseq)
        readout.get_block('prephaser_P_0').scale_gradients(kP * dk_area_scale[1])
        readout.get_block('prephaser_S_0').scale_gradients(kS * dk_area_scale[2])

        readout.get_block('prephaser_P_rewind_0').scale_gradients(kP * dk_area_scale[1])
        readout.get_block('prephaser_S_rewind_0').scale_gradients(kS * dk_area_scale[2])

        readout.remove_block('adc_0')

        seq = deepcopy(rf_seq)
        if venc_duration>0:
            seq.append(bipolars[seg],copy=True)

        seq.append(readout,copy=False)

        # RF quadratic RF spoiling formula based on:
        # 1. Zur Y, Wood ML, Neuringer LJ.
        # Spoiling of transverse magnetization in steady‐state sequences.
        # Magn. Reson. Med. 1991;21:251–263 doi: 10.1002/mrm.1910210210.
        if rf_spoiling and not balanced:
            rf_offset = rf_phase_init/2*(trcount**2+trcount+2)
            seq.get_block('rf_excitation_0').phase_offset = rf_offset

        # If balanced, we adjust RF phase and add rewinder for bipolar at same time
        elif balanced:

            rf_offset = Quantity(np.mod(trcount,2)*np.pi,'rad')
            seq.get_block('rf_excitation_0').phase_offset = rf_offset

            if venc_duration>0:
                # if no rewind flag set, only add delay
                if bal_norewind:
                    bip_rew = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                    duration=venc_duration,
                                                    name='rewind_delay')
                    bip_rew = cmrseq.Sequence([bip_rew],system_specs=system_specs)
                else:
                    bip_rew = deepcopy(bipolars[seg])
                    bip_rew.invert_gradients()

                seq.append(bip_rew)



        # Append this sequence to the list of the current heartbeat
        rampup_list.append(seq)
        trcount += 1

    seq_list.append(rampup_list)

    pbar = tqdm(range(expected_blocks), desc="Loop - Blocks")
    for block in pbar:
        pbar2 = tqdm(range(LUT.shape[3]), desc="Loop - Encoding directions",leave=False)
        for seg in pbar2:

            # Every segment represents a new heartbeat
            seq_beat_list = []
            # Loop over heart phases
            for phase in range(LUT.shape[2]):

                # Loop over profiles per phase
                for prof in range(prof_per_phase):

                    if prof+prof_per_phase*block >= LUT.shape[1]:
                        break

                    # Extract P and S locations and make one TR
                    kP = LUT[0, prof + prof_per_phase * block, phase, seg]
                    kS = LUT[1, prof + prof_per_phase * block, phase, seg]

                    # Store some info for sorting during recon
                    LUT_return[0,lutcount] = kP
                    LUT_return[1, lutcount] = kS
                    LUT_return[2, lutcount] = seg
                    LUT_return[3, lutcount] = phase

                    readout = deepcopy(fastseq)
                    readout.get_block('prephaser_P_0').scale_gradients(kP * dk_area_scale[1])
                    readout.get_block('prephaser_S_0').scale_gradients(kS * dk_area_scale[2])

                    readout.get_block('prephaser_P_rewind_0').scale_gradients(kP * dk_area_scale[1])
                    readout.get_block('prephaser_S_rewind_0').scale_gradients(kS * dk_area_scale[2])

                    seq = deepcopy(rf_seq)
                    if venc_duration>0:
                        seq.append(bipolars[seg],copy=True)

                    seq.append(readout,copy=False)

                    # RF quadratic RF spoiling formula based on:
                    # 1. Zur Y, Wood ML, Neuringer LJ.
                    # Spoiling of transverse magnetization in steady‐state sequences.
                    # Magn. Reson. Med. 1991;21:251–263 doi: 10.1002/mrm.1910210210.
                    if rf_spoiling and not balanced:
                        rf_offset = rf_phase_init/2*(trcount**2+trcount+2)
                        seq.get_block('rf_excitation_0').phase_offset = rf_offset
                        seq.get_block('adc_0').phase_offset = rf_offset

                    # If balanced, we adjust RF phase and add rewinder for bipolar at same time
                    elif balanced:

                        rf_offset = Quantity(np.mod(trcount,2)*np.pi,'rad')
                        seq.get_block('rf_excitation_0').phase_offset = rf_offset
                        seq.get_block('adc_0').phase_offset = rf_offset

                        if venc_duration>0:
                            # if no rewind flag set, only add delay
                            if bal_norewind:
                                bip_rew = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                                duration=venc_duration,
                                                                name='rewind_delay')
                                bip_rew = cmrseq.Sequence([bip_rew],system_specs=system_specs)
                            else:
                                bip_rew = deepcopy(bipolars[seg])
                                bip_rew.invert_gradients()

                            seq.append(bip_rew)



                    # Append this sequence to the list of the current heartbeat
                    seq_beat_list.append(seq)

                    trcount += 1
                    lutcount += 1

            # Once we have gone through all the phases, we are done with one heartbeat
            # This list, representing one heartbeat, is appended to the overall result
            seq_list.append(seq_beat_list)

    return seq_list, LUT_return

generate_4Dflow_LUT

generate_4Dflow_LUT(
    matrix,
    total_prof_per_phase,
    card_phases,
    encoding_segments,
    prof_per_phase,
    spiral_inout: bool = False,
    self_gate: bool = True,
    r_max_search: int = 10,
)

Source code in cmrseq/contrib/_4DFlow.py

def generate_4Dflow_LUT(matrix, total_prof_per_phase, card_phases, encoding_segments,
                        prof_per_phase, spiral_inout:bool=False, self_gate:bool=True,
                        r_max_search: int = 10):

    # Initialize some parameters
    golden_increment = 1.8416

    spiral_twist = 1

    theta = 0
    radius = 0. if spiral_inout else 1.
    sampled_center = False
    prof_counter = 0

    d_r = 1 / (prof_per_phase - 1)

    # Define LUT
    LUT = np.zeros([2, total_prof_per_phase, card_phases, encoding_segments])

    # define sampling matrix
    sampling_mtx = np.zeros([matrix[1], matrix[2], card_phases, encoding_segments])

    mtx_center = (np.floor(matrix[1] / 2).astype(int), np.floor(matrix[2] / 2).astype(int))
    mtx_max = (np.floor((matrix[1] - 1) / 2), np.floor((matrix[2] - 1) / 2))

    # generate map for optimal nearest search
    mi, mj = np.meshgrid(np.arange(-r_max_search, r_max_search + 1), np.arange(-r_max_search, r_max_search + 1))
    r_map = np.sqrt(mi ** 2 + mj ** 2)
    rsearch_map = list(np.unravel_index(np.argsort(r_map.flatten()), np.shape(r_map)))
    rsearch_map[0] = rsearch_map[0] - rsearch_map[0][0]
    rsearch_map[1] = rsearch_map[1] - rsearch_map[1][0]

    for seg in range(encoding_segments):
        for card in range(card_phases):
            for prof in range(total_prof_per_phase):
                rad = radius

                # Jitter radius
                if np.abs(rad)>d_r/2:
                    rad = rad - d_r*np.random.rand()

                # Get profile positions
                py = np.round(np.cos(theta-spiral_twist*np.pi*rad/2)*rad*mtx_max[0]).astype(int)
                pz = np.round(np.sin(theta-spiral_twist*np.pi*rad/2)*rad*mtx_max[1]).astype(int)

                # Convert to sampling matrix coordianates
                iy = py + mtx_center[0]
                iz = pz + mtx_center[1]

                # Keep center sample if not sampled already
                if py==0 and pz==0 and not sampled_center:
                    sampling_mtx[iy,iz,card,seg] = sampling_mtx[iy,iz,card,seg] + 1
                    sampled_center=True
                else: # Search for nearest available point to fill
                    # find closest zero
                    # We loop over the shifts from the point we want
                    for si, sj in zip(rsearch_map[0], rsearch_map[1]):

                        if sampling_mtx[np.clip(iy + si, a_min=0, a_max=sampling_mtx.shape[0] - 1),
                                        np.clip(iz + sj, a_min=0, a_max=sampling_mtx.shape[1] - 1),
                                        card, seg] == 0:
                            # If we find an in-bounds zero, we update that point
                            iy = iy + si
                            iz = iz + sj
                            sampling_mtx[iy, iz, card, seg] = sampling_mtx[iy, iz, card, seg] + 1
                            break
                        # otherwise we will just stay at the original point

                LUT[0, prof, card, seg] = iy - mtx_center[0]
                LUT[1, prof, card, seg] = iz - mtx_center[1]

                # Increment profile counter
                prof_counter += 1

                # If we reach the number of profiles per cardiac interval,
                # the spiral angle is incremented by the golden angle and radius is reset
                if prof_counter == prof_per_phase:
                    prof_counter = 0
                    theta = theta + golden_increment
                    radius = 0. if spiral_inout else 1.
                    # If we are self-gating, also we resample the center point on the next interval
                    if self_gate: sampled_center = False
                else:
                    # Otherwise we just update the radius
                    radius = radius + d_r if spiral_inout else radius - d_r

            # Cardiac phase finished, reset
            prof_counter = 0
            radius = 0. if spiral_inout else 1.
            sampled_center = False

    return LUT, sampling_mtx

gen_WASP

gen_WASP(
    nspokes,
    nframes,
    polartilt=1,
    spiral_density=1,
    rangefactor=1,
    random_seed=None,
)

Source code in cmrseq/contrib/_4DRadial.py

def gen_WASP(nspokes, nframes, polartilt=1, spiral_density=1 , rangefactor=1, random_seed=None):

    rang = rangefactor
    total_spokes = nspokes
    total_frames = nframes
    fmod = spiral_density
    scale = polartilt

    # Golden angle for spiral rotation
    golden = np.pi*(3-np.sqrt(5))


    proj_per_frame = int(np.ceil(nspokes/nframes))
    if np.mod(proj_per_frame,2)==1:
        proj_per_frame += 1
    # First calculate spokes for single spiral
    lspokes = int(proj_per_frame/2)


    #proj_per_frame = int(total_spokes/total_frames) # Spokes per spiral up-down (ie per pole tile)
    #lspokes = int(proj_per_frame/2) # Spokes per spiral
    prev_angle=0

    fMODpar = fmod

    array_mp = np.zeros(lspokes+1) # z position of spoke
    array_ma = np.zeros(lspokes+1) # xy angle of spoke

    # Generate spiral template
    for i in range(lspokes+1):
        dh = -1 + 2 * i / lspokes / rang # goes from -1 to 1 (full sphere) or -1 to 0 (half sphere)
        array_mp[i] = np.arccos(dh)

        if i==0 or (i==lspokes and rangefactor==1):
            array_ma[i] = 0 # set angle to zero at start and end
        else:
            array_ma[i] = np.mod(prev_angle + fMODpar / np.sqrt(lspokes * rang * (1 - dh ** 2)), 2 * np.pi)

        prev_angle = array_ma[i]

    # Generate arrays for full spirals (up and down)
    full_mp = np.zeros(proj_per_frame)
    full_ma = np.zeros(proj_per_frame)


    # Write out full spiral up down
    for i in range(lspokes):

        full_ma[i] = array_ma[i]
        full_mp[i] = array_mp[i]

        full_ma[i + lspokes] = array_ma[lspokes - i] + np.pi
        full_mp[i + lspokes] = array_mp[lspokes - i]


    # Rescale?
    lspokes = lspokes*2

    #polar rotation angle
    polar = 2*np.pi/np.sqrt(lspokes)*scale
    if random_seed != None:
        np.random.seed(random_seed)
    ran = np.random.uniform(0,1,(2,total_frames))

    # Generate polar tilting angles
    pol_x = np.zeros(total_frames)
    pol_y = np.zeros(total_frames)
    azi_z = np.zeros(total_frames)

    for i in range(total_frames):

        flip_polar = polar*ran[0,i]
        flip_phi = 2*np.pi*ran[1,i]
        pol_x[i] = flip_polar*np.sin(flip_phi)
        pol_y[i] = flip_polar*np.cos(flip_phi)
        azi_z[i] = (i+1)*golden


    # Transform spiral into x,y,z, coordinates
    x_b = np.zeros(proj_per_frame)
    y_b = np.zeros(proj_per_frame)
    z_b = np.zeros(proj_per_frame)

    for i in range(proj_per_frame):
        x_b[i] = np.sin(full_mp[i]) * np.cos(full_ma[i])
        y_b[i] = np.sin(full_mp[i]) * np.sin(full_ma[i])
        z_b[i] = np.cos(full_mp[i])


    # Rotate all each frame by polar tilting
    x = np.zeros((total_frames,proj_per_frame))
    y = np.zeros((total_frames,proj_per_frame))
    z = np.zeros((total_frames,proj_per_frame))

    for j in range(total_frames):
        for i in range(proj_per_frame):
            x_t = x_b[i]
            y_t = np.cos(pol_x[j])*y_b[i] + np.sin(pol_x[j])*z_b[i]
            z_t = -np.sin(pol_x[j]) * y_b[i] + np.cos(pol_x[j]) * z_b[i]

            x_i = np.cos(pol_y[j]) * x_t + np.sin(pol_y[j]) * z_t
            y_i = y_t
            z_i = -np.sin(pol_y[j]) * x_t + np.cos(pol_y[j]) * z_t

            x_t = np.cos(azi_z[j]) * x_i + np.sin(azi_z[j]) * y_i
            y_t = -np.sin(azi_z[j]) * x_i + np.cos(azi_z[j]) * y_i
            z_t = z_i

            x[j,i] = x_t
            y[j,i] = y_t
            z[j,i] = z_t

    traj = (x,y,z)
    pole = (pol_x,pol_y,azi_z)
    spirals = (x_b,y_b,z_b)

    print(str(proj_per_frame)+' projections per interleave')
    print(str(int(proj_per_frame/2))+' projections per spiral direction')
    print(str(total_frames*proj_per_frame)+ ' / ' + str(total_spokes) + ' spokes populated')
    print(f"Max polar tilt angle of {polar:.2f} radians")
    wraps = fmod*(np.sqrt(proj_per_frame)*0.25-0.233)
    print(f"{(2*wraps):.2f} wraps per spiral up-down")
    return traj,pole,spirals

radial_bSSFP_3D_WASP

radial_bSSFP_3D_WASP(
    system_specs: SystemSpec,
    samples_per_spoke: int,
    resolution: Quantity,
    adc_duration: Quantity,
    flip_angle: Quantity,
    pulse_duration: Quantity,
    num_interleaves=int,
    spokes_per_interleave=int,
    dummy_shots: int = 0,
    half_sphere=False,
    polar_tilt=1.0,
    spiral_density=1.0,
    add_bipolar=False,
    venc_duration=None,
    venc_strength=None,
    venc_direction=np.array([0.0, 1.0, 0.0]),
    m1_compensation=None,
    random_seed=None,
    TR=None,
    TE=None,
    reverse_order=False,
    readout_ramp_sampling: bool = False,
    sample_prephaser: bool = False,
    sample_balanced: bool = False,
    prephaser_ramp_sampling: bool = False,
    balanced_ramp_sampling: bool = False,
    match_prephaser_amp_to_readout: bool = False,
    make_GRE: bool = False,
    gre_dephase_phase: Quantity = Quantity(
        2 * np.pi, "rad"
    ),
    gre_spoil_along_readout: bool = True,
    disable_bSSFP_phase_alternation_in_GRE: bool = True,
    delaycal: bool = False,
    num_delaycal: int = 16,
    delaycal_reps=10,
) -> List[cmrseq.Sequence]

Source code in cmrseq/contrib/_4DRadial.py

def radial_bSSFP_3D_WASP(system_specs: cmrseq.SystemSpec,
                     samples_per_spoke: int,
                     resolution: Quantity,
                     adc_duration: Quantity,
                     flip_angle: Quantity,
                     pulse_duration: Quantity,
                     num_interleaves = int,
                     spokes_per_interleave = int,
                     dummy_shots: int = 0,
                     half_sphere=False,
                     polar_tilt=1.,
                     spiral_density=1.,
                     add_bipolar=False,
                     venc_duration = None,
                     venc_strength = None,
                     venc_direction = np.array([0.,1.,0.]),
                     m1_compensation = None, # Defines the M1 balancing at TE. Either None, "seperate" (balancing seperate from VENC) or "combined_highgrad"
                     random_seed = None,
                     TR=None,
                     TE=None, # TE of the main readout
                     reverse_order=False,
                     readout_ramp_sampling: bool = False,
                     sample_prephaser: bool = False,
                     sample_balanced: bool = False,
                     prephaser_ramp_sampling: bool = False,
                     balanced_ramp_sampling: bool = False,
                     match_prephaser_amp_to_readout: bool = False,
                     make_GRE: bool = False,
                     gre_dephase_phase: Quantity = Quantity(2*np.pi, "rad"),   # “π across one voxel” by default
                     gre_spoil_along_readout: bool = True,
                     disable_bSSFP_phase_alternation_in_GRE: bool = True,
                     delaycal: bool = False,
                     num_delaycal: int = 16,
                     delaycal_reps = 10,
                     ) -> List[cmrseq.Sequence]:

    # Input handling
    if m1_compensation == "None":
        m1_compensation = None

    if venc_strength is None and venc_duration is None and add_bipolar:
        raise ValueError("Either venc_strength or venc_duration must be set to use bipolar gradients.")
    if venc_strength is None:
        venc_strength = Quantity(0., "mT/m")
    if venc_duration is None:
        venc_duration = Quantity(0., "ms")

    venc_direction = venc_direction/np.linalg.norm(venc_direction)

    if not add_bipolar and m1_compensation=="combined":
        # This is a reduntant case
        m1_compensation = "seperate"

    # Step 1: RF Hardpulses
    rf_block = cmrseq.bausteine.HardRFPulse(system_specs=system_specs, flip_angle=flip_angle,
                                            duration=pulse_duration,
                                            delay=Quantity(0.,'ms'),name="rf_excitation")

    rf_seq = cmrseq.Sequence([rf_block],system_specs=system_specs)

    # FA/2 RF hardpulse
    rf_block2 = cmrseq.bausteine.HardRFPulse(system_specs=system_specs, flip_angle=flip_angle/2,
                                             duration=pulse_duration,
                                             delay=Quantity(0.,'ms'),name="rf_catalyst")

    rf_seq2 = cmrseq.Sequence([rf_block2],system_specs=system_specs)


    # Step 2: Simple case for non-merged bipolars
    # Bipolar gradients, only if they are needed and not merged with M1 compensation

    bipolar = None
    bipolar_rewind = None
    if add_bipolar and m1_compensation != "combined":

        if venc_duration>0 and venc_strength==0:
            # Add delay only
            bipolar = cmrseq.bausteine.Delay(system_specs=system_specs, duration=venc_duration,
                                             name="velocity_encode_delay")
            bipolar = cmrseq.Sequence([bipolar],system_specs=system_specs)
            bipolar_rewind = cmrseq.bausteine.Delay(system_specs=system_specs, duration=venc_duration,
                                                    name="velocity_encode_rewinder_delay")
            bipolar_rewind = cmrseq.Sequence([bipolar_rewind],system_specs=system_specs)
        elif venc_strength>0:
            # Create bipolars
            bipolar = cmrseq.parametric_definitions.velocity.bipolar(system_specs=system_specs,
                                                                    venc=venc_strength,
                                                                    duration = venc_duration,
                                                                    direction=venc_direction)

            bipolar_rewind = cmrseq.parametric_definitions.velocity.bipolar(system_specs=system_specs,
                                                                            venc=venc_strength,
                                                                            duration = venc_duration,
                                                                            direction=-venc_direction)

        # Step 3: Generate 3D radial spokes
    kr_max = 1 / (2 * resolution)

    total_spokes = num_interleaves * spokes_per_interleave
    rangefactor = 2 if half_sphere else 1

    if total_spokes>0:
        traj_wasp, _, _ = gen_WASP(
            total_spokes,
            num_interleaves,
            polar_tilt,
            spiral_density,
            rangefactor=rangefactor,
            random_seed=random_seed,
        )
        traj = np.stack([d.flatten() for d in traj_wasp], axis=1)     # Flatten nominal WASP to (N, 3)

    else:
        traj = np.empty((0,3))




    # Prepend delay-calibration spokes, matching the C++ intent:
    # acquire 3 * num_delaycal orthogonal radial calibration spokes first,
    # then continue with nominal WASP ordering unchanged.
    if delaycal:
        delaycal_traj, num_delaycal = gen_delaycal_spokes(num_delaycal)
        delaycal_traj = np.repeat(delaycal_traj,delaycal_reps, axis = 0)
        traj = np.concatenate([delaycal_traj, traj], axis=0)
        print(
            f"Prepended delay calibration: {delaycal_traj.shape[0]} spokes "
            f"(3 x { num_delaycal } x {delaycal_reps})"
        )

    if traj.shape[0]==0:
        print('Warning: trajectory is empty. Adding a spoke in kz: [0,0,3]')
        traj = np.array([[0, 0, 1]])
    # Generate actual readout blocks
    ro_blocks = cmrseq.parametric_definitions.readout.radial_3D(system_specs=system_specs,
                                                                spoke_directions=traj,
                                                                samples_per_spoke=samples_per_spoke,
                                                                kr_max=kr_max,
                                                                adc_duration=adc_duration,
                                                                balanced=True,
                                                                readout_ramp_sampling=readout_ramp_sampling,
                                                                sample_prephaser=sample_prephaser,
                                                                sample_balanced=sample_balanced,
                                                                prephaser_ramp_sampling=prephaser_ramp_sampling,
                                                                balanced_ramp_sampling=balanced_ramp_sampling,
                                                                match_prephaser_amp_to_readout=match_prephaser_amp_to_readout)

    # Step 4: If compensated, re-calculate prephaser timings
    if m1_compensation == "seperate" or m1_compensation == "combined":
        # Get M1 and M0 of just readout gradient up to echo time
        # We only care about the magnitude
        ro = deepcopy(ro_blocks[0])
        if not reverse_order:
            ro.remove_block('radial_prephaser_0')
            if sample_prephaser:
                ro.remove_block('adc_prephaser_0')
            if sample_balanced or sample_prephaser:
                ro.remove_block('adc_dead_delay_0')
        ro.remove_block('radial_prephaser_balance_0')
        if sample_balanced:
            ro.remove_block('adc_balanced_0')
        if sample_balanced or sample_prephaser:
            ro.remove_block('adc_dead_delay_1')
        ro.shift_in_time(-ro.start_time)
        echo_time = ro['adc_0'].adc_center
        readout_M1 = Quantity(np.linalg.norm(ro.calculate_moment(1,end_time=echo_time).m_as('mT/m*ms**2')),'mT/m*ms**2')
        readout_M0 = Quantity(np.linalg.norm(ro.calculate_moment(0,end_time=echo_time).m_as('mT/m*ms')),'mT/m*ms')
        if reverse_order:
            readout_M0 = Quantity(0.,'mT/m*ms')
        desired_M1 = Quantity(0.,'mT/m*ms**2')
        if m1_compensation == "combined" and venc_strength > 0:
            # In this case, the compensation should account for the need to add the velocity encoding M1
            desired_M1 = Quantity(np.pi,'rad')/system_specs.gamma_rad/venc_strength


        # Calculate the prephaser timings
         # The worst case scenario is when the desired M1 is negative (encoding opposite to readout)
        T1, T2, d1, d2 = _optimize_prephaser_timing(system_specs=system_specs,
                                                    readout_M1=readout_M1,
                                                    readout_M0=readout_M0,
                                                    desired_M1=-desired_M1)
    # Step4b: Calculate the desired M1 for MPS
    if m1_compensation == "combined" and venc_strength > 0:
        desired_M1_MPS = desired_M1 * venc_direction
    else:
        desired_M1_MPS = Quantity(np.array([0.,0.,0.]),'mT/m*ms**2')



    # Step 5: Assemble sequence
    seq_list = []
    # for ro_idx, ro_b in enumerate(ro_blocks):
    #     # Start with RF pulse
    #     seq = deepcopy(rf_seq)

    #     # Adjust RF and phase
    #     seq['rf_excitation_0'].phase_offset = Quantity(np.mod(ro_idx,2)*np.pi,'rad')

    rf_dead_ms = system_specs.rf_dead_time.m_as("ms")
    rf_ringdown_ms = system_specs.rf_ringdown_time.m_as("ms")
    rf_pad_ms = float(np.max([rf_dead_ms, rf_ringdown_ms, 0.0]))

    for ro_idx, ro_b in enumerate(ro_blocks):
        # Build a subsequence starting with [delay → RF → delay]
        rf_delay_before = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                 duration=Quantity(rf_pad_ms, "ms"),
                                                 name="rf_deadtime_pre")
        rf_delay_after = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                duration=Quantity(rf_pad_ms, "ms"),
                                                name="rf_deadtime_post")

        # deepcopy RF to avoid reuse issues
        rf_exc = deepcopy(rf_seq['rf_excitation_0'])
        rf_exc.phase_offset = Quantity(np.mod(ro_idx, 2) * np.pi, 'rad')
        # Create sequence: [delay → RF → delay]
        rf_subseq = cmrseq.Sequence([], system_specs=system_specs)
        rf_subseq.append(rf_delay_before)
        rf_subseq.append(rf_exc)
        rf_subseq.append(rf_delay_after)

        # Make a working sequence from that
        seq = cmrseq.Sequence([], system_specs=system_specs)
        seq.append(rf_subseq)

        ro_b['adc_0'].phase_offset = Quantity(np.mod(ro_idx,2)*np.pi,'rad')
        if sample_prephaser:
            ro_b['adc_prephaser_0'].phase_offset = Quantity(np.mod(ro_idx,2)*np.pi,'rad')
        if sample_balanced:
            ro_b['adc_balanced_0'].phase_offset = Quantity(np.mod(ro_idx,2)*np.pi,'rad')

        # Add standalone bipolar if needed

        if reverse_order:
            prephaser = deepcopy(ro_b)
            prephaser.remove_block('radial_prephaser_balance_0')
            if sample_balanced:
                prephaser.remove_block('adc_balanced_0')
            if sample_balanced or sample_prephaser:
                prephaser.remove_block('adc_dead_delay_1')
            prephaser.remove_block('radial_readout_0')
            prephaser.remove_block('adc_0')
            prephaser.shift_in_time(-prephaser.start_time)
            seq.append(prephaser, copy=True)

            prephaser_balance = deepcopy(ro_b)
            prephaser_balance.remove_block('radial_prephaser_0')
            if sample_prephaser:
                prephaser_balance.remove_block('adc_prephaser_0')
            if sample_balanced or sample_prephaser:
                prephaser_balance.remove_block('adc_dead_delay_0')
            prephaser_balance.remove_block('radial_readout_0')
            prephaser_balance.remove_block('adc_0')
            prephaser_balance.shift_in_time(-prephaser_balance.start_time)

        if add_bipolar and m1_compensation != "combined":
            # the the first bipolar
            seq.append(bipolar, copy=True)

        # If compensating, need to create the new pre phasers and rewinders, per direction
        cur_readout_dir = traj[ro_idx]
        cur_readout_dir = cur_readout_dir/np.linalg.norm(cur_readout_dir)
        if m1_compensation == "seperate" or m1_compensation == "combined":
            p1_MPS = []
            p2_MPS = []
            for i in range(3):
                # First we need to find the component of the readout along the given direction
                cur_readout_M1 = cur_readout_dir[i]*readout_M1
                cur_readout_M0 = cur_readout_dir[i]*readout_M0
                G1,G2 = _calculate_optimized_prephasers(system_specs=system_specs,
                                                        readout_M1 = cur_readout_M1,
                                                        readout_M0 = cur_readout_M0,
                                                        desired_M1 = desired_M1_MPS[i],
                                                        T1=T1, T2=T2, d1=d1, d2=d2)
                p1_MPS.append(G1.m_as('mT/m'))
                p2_MPS.append(G2.m_as('mT/m'))
            p1_MPS = np.array(p1_MPS)
            p2_MPS = np.array(p2_MPS)

            # Create and append the prephasers
            p1_mag = np.linalg.norm(p1_MPS)
            p2_mag = np.linalg.norm(p2_MPS)
            if p1_mag>0:
                p1_dir = p1_MPS/p1_mag
                prephaser_1 = cmrseq.bausteine.TrapezoidalGradient(system_specs=system_specs,
                                                                   orientation=p1_dir,
                                                                   amplitude=Quantity(p1_mag,'mT/m'),
                                                                   flat_duration=T1,
                                                                   rise_time=d1,
                                                                   name='radial_prephaser_comp1')
            else:
                # need to add a delay
                prephaser_1 = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                    duration=T1 + 2*d1,
                                                    name='radial_prephaser_comp1')
            seq.append(prephaser_1)

            if p2_mag>0:
                p2_dir = p2_MPS/p2_mag
                prephaser_2 = cmrseq.bausteine.TrapezoidalGradient(system_specs=system_specs,
                                                                   orientation=p2_dir,
                                                                   amplitude=Quantity(p2_mag,'mT/m'),
                                                                   flat_duration=T2,
                                                                   rise_time=d2,
                                                                   name='radial_prephaser_comp2')
            else:
                # need to add a delay
                prephaser_1 = cmrseq.bausteine.Delay(system_specs=system_specs,
                                                    duration=T2 + 2*d2,
                                                    name='radial_prephaser_comp2')
            seq.append(prephaser_2)

            # Add the readout (remove prephaser/rewinders first)
            ro_b.remove_block('radial_prephaser_0')
            if sample_prephaser:
                ro_b.remove_block('adc_prephaser_0')
            ro_b.remove_block('radial_prephaser_balance_0')
            if sample_balanced:
                ro_b.remove_block('adc_balanced_0')
            if sample_balanced or sample_prephaser:
                ro_b.remove_block('adc_dead_delay_1')
                ro_b.remove_block('adc_dead_delay_0')
            ro_b.shift_in_time(-ro_b.start_time)
            seq.append(ro_b)

            # Add the rewinders
            if p2_mag>0:
                rewinder_2 = deepcopy(prephaser_2)
                rewinder_2.name = 'radial_rewinder_comp2'
                seq.append(rewinder_2)

            if p1_mag>0:
                rewinder_1 = deepcopy(prephaser_1)
                rewinder_1.name = 'radial_rewinder_comp1'
                seq.append(rewinder_1)

        else:
            # No compensation, we simply add the readout as is
            if reverse_order:
                ro_b.remove_block('radial_prephaser_0')
                if sample_prephaser:
                    ro_b.remove_block('adc_prephaser_0')
                ro_b.remove_block('radial_prephaser_balance_0')
                if sample_balanced:
                    ro_b.remove_block('adc_balanced_0')
                if sample_balanced or sample_prephaser:
                    ro_b.remove_block('adc_dead_delay_1')
                    ro_b.remove_block('adc_dead_delay_0')
                ro_b.shift_in_time(-ro_b.start_time)
                seq.append(ro_b)
            else:
                seq.append(ro_b)


        if add_bipolar and m1_compensation != "combined":
            # the the bipolar rewinder
            seq.append(bipolar_rewind, copy=True)

        if reverse_order:
            seq.append(prephaser_balance, copy=True)

        # Handle setting custom TE and TR (clunky code, maybe revisit)
        if TE is not None and TR is not None: # Add delay before readout to match TE and then after readout to match TR
            if ro_idx == 0:
                print("Careful, setting both TE and TR creates an assymetric TR --> TE is not at TR//2")
            current_te = seq.get_block('adc_0').anchor_time - seq.get_block('rf_excitation_0').tmin - seq.get_block('rf_excitation_0').duration/2
            delay_duration = ((TE - current_te))//system_specs.grad_raster_time* system_specs.grad_raster_time
            if delay_duration > 0:
                delay = cmrseq.bausteine.Delay(system_specs=system_specs,
                                               duration=delay_duration,
                                               name="repetition_delay")
                for block in seq[2:]: # shift everything after rf blocks
                    block.shift(time_shift=delay_duration)

            elif delay_duration < 0 and ro_idx == 0:
                print(f"Warning: Echo time {TE} is smaller than current echo time {current_te}. Requested echo time could not be handled. Setting to minimum echo time.")
            current_tr = seq.duration
            delay_duration = ((TR - current_tr))//system_specs.grad_raster_time* system_specs.grad_raster_time
            if delay_duration > 0:
                delay = cmrseq.bausteine.Delay(system_specs=system_specs,
                                               duration=delay_duration,
                                               name="repetition_delay")
                seq.append(delay)
            elif delay_duration < 0 and ro_idx == 0:
                print(f"Warning: Repetition time {TR} is smaller than sequence duration {current_tr}. Requested repetition time could not be handled. Setting to minimum repetition time")
        if TR is None and TE is not None: # Add delay before and after readout to match TE
            current_te = seq.get_block('adc_0').anchor_time - seq.get_block('rf_excitation_0').tmin - seq.get_block('rf_excitation_0').duration/2
            delay_duration = ((TE - current_te))//system_specs.grad_raster_time* system_specs.grad_raster_time
            if delay_duration > 0:
                delay = cmrseq.bausteine.Delay(system_specs=system_specs,
                                               duration=delay_duration,
                                               name="repetition_delay")
                for block in seq[2:]: # shift everything after rf blocks
                    block.shift(time_shift=delay_duration)

                seq.append(delay)
            elif delay_duration < 0 and ro_idx == 0:
                print(f"Warning: Echo time {TE} is smaller than current echo time {current_te}. Requested echo time could not be handled. Setting to minimum echo time.")
        if TR is not None and TE is None: # Add delay before and after readout to match TR
            current_tr = seq.duration
            delay_duration = ((TR - current_tr)/2)//system_specs.grad_raster_time* system_specs.grad_raster_time
            if delay_duration > 0:
                delay = cmrseq.bausteine.Delay(system_specs=system_specs,
                                               duration=delay_duration,
                                               name="repetition_delay")
                for block in seq[2:]: # shift everything after rf blocks
                    block.shift(time_shift=delay_duration)

                seq.append(delay)
            elif delay_duration < 0 and ro_idx == 0:
                print(f"Warning: Repetition time {TR} is smaller than sequence duration {current_tr}. Requested repetition time could not be handled. Setting to minimum repetition time")


        if make_GRE:
            # remove bSSFP balancing/rewinding
            # for b in list(seq):
            #     n = getattr(b, "name", "").lower()
            #     if ("balance" in n) or ("rewinder" in n):
            #         try:
            #             seq.remove_block(getattr(b, "name", ""))
            #         except Exception:
            #             pass
            # for bn in ["radial_prephaser_balance_0", "adc_balanced_0", "adc_dead_delay_1"]:
            #     try:
            #         seq.remove_block(bn)
            #     except Exception:
            #         pass

            # RF spoiling: quadratic phase progression with 117° increment
            if disable_bSSFP_phase_alternation_in_GRE:
                phi0 = np.deg2rad(117.0)
                rf_inc = np.mod(ro_idx * phi0, 2*np.pi)
                rf_phase = np.mod(ro_idx * (ro_idx - 1) / 2 * phi0, 2*np.pi)
                try:
                    seq.get_block("rf_excitation_0").phase_offset = Quantity(rf_phase, "rad")
                except Exception:
                    pass
                for bn in ["adc_0", "adc_prephaser_0", "adc_balanced_0"]:
                    try:
                        seq.get_block(bn).phase_offset = Quantity(rf_phase, "rad")
                    except Exception:
                        pass

            # spoiler gradient
            spoildir = traj[ro_idx] / np.linalg.norm(traj[ro_idx]) if gre_spoil_along_readout else np.array([0.0, 0.0, 1.0])
            dx = resolution.to("m").m_as("m")
            phi = gre_dephase_phase.to("rad").m_as("rad")
            gamma_rad = system_specs.gamma_rad.m_as("rad/T/s")
            crusher_area = Quantity(phi / (gamma_rad * dx) * 1e6, "mT/m*ms")

            spoiler = cmrseq.bausteine.TrapezoidalGradient.from_area(
                system_specs=system_specs,
                orientation=spoildir.astype(float),
                area=abs(crusher_area),
                delay=0,
                name="gre_spoiler",
            )

            spoiler.shift(seq.duration//system_specs.grad_raster_time* system_specs.grad_raster_time-spoiler.duration)
            seq += cmrseq.Sequence([spoiler], system_specs=system_specs)
        # --- end GRE post-pass ---

        seq_list.append(seq)

    # Set up dummy shots
    dummy_ref = deepcopy(seq_list[0])
    dummy_ref.remove_block('adc_0')
    if sample_prephaser:
        dummy_ref.remove_block('adc_prephaser_0')
    if sample_balanced:
        dummy_ref.remove_block('adc_balanced_0')
    for i in range(dummy_shots):
        cur_dummy = deepcopy(dummy_ref)
        cur_dummy['rf_excitation_0'].phase_offset = Quantity(np.mod(i+1,2)*np.pi,'rad')
        seq_list.insert(0, cur_dummy)

    # Add catalyst
    rf_catalyst = deepcopy(rf_seq2['rf_catalyst_0'])

    rf_delay_before_t =(Quantity(rf_pad_ms , "ms")+ seq_list[0].get_block('rf_excitation_0').duration/2 - rf_catalyst.duration/2)// system_specs.grad_raster_time * system_specs.grad_raster_time
    rf_delay_after_t = Quantity(rf_pad_ms, "ms")
    catalyst_delay_t = (seq_list[0].duration/2-rf_delay_before_t-rf_delay_after_t-rf_catalyst.duration) // system_specs.grad_raster_time * system_specs.grad_raster_time

    rf_delay_before = cmrseq.bausteine.Delay(system_specs=system_specs,
                                            duration=rf_delay_before_t,
                                            name="rf_deadtime_pre")
    rf_delay_after = cmrseq.bausteine.Delay(system_specs=system_specs,
                                            duration=rf_delay_after_t,
                                            name="rf_deadtime_post")
    catalyst_delay = cmrseq.bausteine.Delay(system_specs=system_specs,
                                            duration=catalyst_delay_t,
                                            name="catalyst_delay")


    rf_catalyst.phase_offset = Quantity(np.pi * np.mod(dummy_shots+1, 2), 'rad')
    rf_catalyst_subseq = cmrseq.Sequence([], system_specs=system_specs)
    rf_catalyst_subseq.append(rf_delay_before)
    rf_catalyst_subseq.append(rf_catalyst)
    rf_catalyst_subseq.append(rf_delay_after)
    rf_catalyst_subseq.append(catalyst_delay)
    # rf_seq2.append(catalyst_delay)
    seq_list.insert(0, rf_catalyst_subseq)

    return seq_list

pc_gre

pc_gre(
    system_specs: SystemSpec,
    matrix_size: ndarray,
    inplane_resolution: Quantity,
    slice_thickness: Quantity,
    adc_duration: Quantity,
    flip_angle: Quantity,
    pulse_duration: Quantity,
    repetition_time: Quantity,
    echo_time: Quantity,
    venc: Quantity,
    venc_direction: ndarray,
    venc_duration: Quantity = Quantity(0.0, "ms"),
    slice_position_offset: Quantity = Quantity(0.0, "m"),
    time_bandwidth_product: float = 4.0,
    dummy_shots: int = None,
    crusher_area: Quantity = Quantity(0.0, "mT/m*ms"),
    crusher_duration: Quantity = Quantity(0.0, "mT/m*ms"),
) -> List[cmrseq.Sequence]

Defines a 2D gradient echo sequence with bipolar velocity encoding.

Parameters:

Name	Type	Description	Default
system_specs	SystemSpec	SystemSpecifications	required
matrix_size	ndarray	array of shape (2, ) containing the resulting matrix dimensions	required
inplane_resolution	Quantity	Quantity[Length] of shape (2, ) containing the in-plane voxel dimensions	required
slice_thickness	Quantity	Quantity[Length] containing the required slice-thickness	required
adc_duration	Quantity	Quantity[time] Total duration of adc-sampling for a single TR	required
repetition_time	Quantity	Quantity[Time] containing the required repetition_time	required
echo_time	Quantity	Quantity[Time] containing the required echo-time. If too short for given system specifications, it is increased to minimum and a warning is raised.	required
venc	Quantity	Quantity[Velocity] strength of velocity encoding gradient	required
venc_duration	Quantity	Quantity[Time] denoting the duration of applied VENC-gradients. If 0. the resulting gradients will be the shortest for given system limits	Quantity(0.0, 'ms')
venc_direction	ndarray	Vector (3, ) denoting the direction of velocity encoding in MPS coordinates	required
flip_angle	Quantity	Quantity[Angle] containing the required flip_angle	required
pulse_duration	Quantity	Quantity[Time] Total pulse duration (corresponds to flat_duration of the slice selection gradient)	required
slice_position_offset	Quantity	Quantity[Length] positional offset in slice normal direction defining the frequency offset of the RF pulse	Quantity(0.0, 'm')
time_bandwidth_product	float	float used to calculate the rf bandwidth from duration	4.0
dummy_shots	int	number of dummy shots (TRs) without adc-events, with k-space center phase encoding	None
crusher_area	Quantity	Quantity[Gradient Area] crusher gradient area along slice direction. If set to zero no crusher will be applied and phase encoder will not be rewound	Quantity(0.0, 'mT/m*ms')
crusher_duration	Quantity	Quantity[Time] duration of crusher. If set too short will default to duration of phase encoder or shortest possible crusher	Quantity(0.0, 'mT/m*ms')

Returns:

Type	Description
List of sequence objects, that each represent a single TR

Source code in cmrseq/contrib/_2D_flow.py

def pc_gre(system_specs: cmrseq.SystemSpec,
           matrix_size: np.ndarray,
           inplane_resolution: Quantity,
           slice_thickness: Quantity,
           adc_duration: Quantity,
           flip_angle: Quantity,
           pulse_duration: Quantity,
           repetition_time: Quantity,
           echo_time: Quantity,
           venc: Quantity,
           venc_direction: np.ndarray,
           venc_duration: Quantity = Quantity(0., "ms"),
           slice_position_offset: Quantity = Quantity(0., "m"),
           time_bandwidth_product: float = 4.,
           dummy_shots: int = None,
           crusher_area: Quantity = Quantity(0.,'mT/m*ms'),
           crusher_duration:Quantity = Quantity(0.,'mT/m*ms')) -> List[cmrseq.Sequence]:
    r"""Defines a 2D gradient echo sequence with bipolar velocity encoding.

    Parameters
    ----------
    system_specs
        SystemSpecifications
    matrix_size
        array of shape (2, ) containing the resulting matrix dimensions
    inplane_resolution
        Quantity[Length] of shape (2, ) containing the in-plane voxel dimensions
    slice_thickness
        Quantity[Length] containing the required slice-thickness
    adc_duration
        Quantity[time] Total duration of adc-sampling for a single TR
    repetition_time
        Quantity[Time] containing the required repetition_time
    echo_time
        Quantity[Time] containing the required echo-time. If too short for given system specifications, it is increased to minimum and a warning is raised.
    venc
        Quantity[Velocity] strength of velocity encoding gradient
    venc_duration
        Quantity[Time] denoting the duration of applied VENC-gradients. If 0. the resulting gradients will be the shortest for given system limits
    venc_direction
        Vector (3, ) denoting the direction of velocity encoding in MPS coordinates
    flip_angle
        Quantity[Angle] containing the required flip_angle
    pulse_duration
        Quantity[Time] Total pulse duration (corresponds to flat_duration of the slice selection gradient)
    slice_position_offset
        Quantity[Length] positional offset in slice normal direction defining the frequency offset of the RF pulse
    time_bandwidth_product
        float used to calculate the rf bandwidth from duration
    dummy_shots
        number of dummy shots (TRs) without adc-events, with k-space center phase encoding
    crusher_area
        Quantity[Gradient Area] crusher gradient area along slice direction. If set to zero no crusher will be applied and phase encoder will not be rewound
    crusher_duration
        Quantity[Time] duration of crusher. If set too short will default to duration of phase encoder or shortest possible crusher

    Returns
    -------
    List of sequence objects, that each represent a single TR
    """
    rf_seq = cmrseq.seqdefs.excitation.slice_selective_sinc_pulse(system_specs=system_specs,
                                                                  slice_thickness=slice_thickness,
                                                                  flip_angle=flip_angle,
                                                                  pulse_duration=pulse_duration,
                                                                  time_bandwidth_product=time_bandwidth_product,
                                                                  slice_position_offset=slice_position_offset,
                                                                  slice_normal=np.array([0., 0., 1.]))

    ro_blocks = cmrseq.seqdefs.readout.multi_line_cartesian(system_specs=system_specs,
                                                            fnc=cmrseq.seqdefs.readout.gre_cartesian_line,
                                                            matrix_size=matrix_size,
                                                            inplane_resolution=inplane_resolution,
                                                            adc_duration=adc_duration,
                                                            dummy_shots=dummy_shots)

    venc_gradient = cmrseq.seqdefs.velocity.bipolar(system_specs=system_specs, venc=venc, duration=venc_duration,
                                                    direction=venc_direction)

    if crusher_area != 0:
        crusher_area = crusher_area.to("mT/m*ms") - rf_seq.get_block("slice_select_0").area[2] / 2
        if crusher_duration == 0:
            crusher_duration = ro_blocks[0].get_block("pe_prephaser_0").duration
        if crusher_duration < (2*system_specs.get_shortest_gradient(crusher_area)[1] + system_specs.get_shortest_gradient(crusher_area)[2]):
            crusher = cmrseq.bausteine.TrapezoidalGradient.from_area(system_specs, orientation=np.sign(crusher_area)*np.array([0., 0., 1.]),
                                                                     area=np.abs(crusher_area.to("mT/m*ms")),
                                                                     name="crusher")
        else:
            crusher = cmrseq.bausteine.TrapezoidalGradient.from_dur_area(system_specs,
                                                                         orientation=np.sign(crusher_area)*np.array([0., 0., 1.]),
                                                                         duration=crusher_duration,
                                                                         area=np.abs(crusher_area.to("mT/m*ms")),
                                                                         name="crusher")
        for seq in ro_blocks:
            pe_rewind = deepcopy(seq.get_block("pe_prephaser_0"))
            pe_rewind.name = "pe_rewind"
            pe_rewind.scale_gradients(-1)

            post = cmrseq.Sequence([crusher,pe_rewind],system_specs=system_specs)
            seq.append(post)

    adc_center = system_specs.time_to_raster(ro_blocks[0].get_block('adc_0').adc_center)

    minimal_tr = ro_blocks[0].duration + venc_gradient.duration + rf_seq.duration
    minimal_te = (rf_seq.duration - rf_seq.get_block("rf_excitation_0").rf_events[0] + adc_center + venc_gradient.duration)

    repetition_time = system_specs.time_to_raster(repetition_time)
    if repetition_time < minimal_tr:
        warn(f"PC GRE Sequence: Repetition time too short to be feasible, set TR to {minimal_tr}")
        repetition_time = minimal_tr

    maximum_te = repetition_time - (ro_blocks[0].duration - adc_center) \
                 - rf_seq.get_block("rf_excitation_0").rf_events[0]

    echo_time = system_specs.time_to_raster(echo_time)
    if echo_time < minimal_te:
        warn(f"PC GRE Sequence: Echo time too short to be feasible, set TE to {minimal_te}")
        echo_time = minimal_te

    if echo_time > maximum_te:
        warn(f"PC GRE Sequence: Echo time too long for given TR, set TE to {maximum_te}")
        echo_time = maximum_te

    te_delay = system_specs.time_to_raster(echo_time - minimal_te)
    tr_delay = system_specs.time_to_raster(repetition_time - minimal_tr - te_delay)

    # Concatenate readout blocks
    seq_list = []
    for ro_b in ro_blocks:
        seq = deepcopy(rf_seq)
        if te_delay > 0:
            seq.append(cmrseq.bausteine.Delay(system_specs=system_specs, duration=te_delay))
        seq.append(venc_gradient)
        seq.append(ro_b)
        if tr_delay > 0:
            seq.append(cmrseq.bausteine.Delay(system_specs=system_specs, duration=tr_delay))
        seq_list.append(seq)
    return seq_list

pc_gre_multivenc

pc_gre_multivenc(
    system_specs: SystemSpec,
    matrix_size: ndarray,
    inplane_resolution: Quantity,
    slice_thickness: Quantity,
    adc_duration: Quantity,
    flip_angle: Quantity,
    pulse_duration: Quantity,
    repetition_time: Quantity,
    echo_time: Quantity,
    venc_list: List[Quantity],
    venc_direction_list: List[ndarray],
    venc_duration: Quantity = Quantity(0.0, "ms"),
    slice_position_offset: Quantity = Quantity(0.0, "m"),
    time_bandwidth_product: float = 4.0,
    dummy_shots: int = None,
    crusher_area: Quantity = Quantity(0.0, "mT/m*ms"),
    crusher_duration: Quantity = Quantity(0.0, "mT/m*ms"),
) -> List[List[cmrseq.Sequence]]

Source code in cmrseq/contrib/_2D_flow.py

def pc_gre_multivenc(system_specs: cmrseq.SystemSpec,
                     matrix_size: np.ndarray,
                     inplane_resolution: Quantity,
                     slice_thickness: Quantity,
                     adc_duration: Quantity,
                     flip_angle: Quantity,
                     pulse_duration: Quantity,
                     repetition_time: Quantity,
                     echo_time: Quantity,
                     venc_list: List[Quantity],
                     venc_direction_list: List[np.ndarray],
                     venc_duration: Quantity = Quantity(0., "ms"),
                     slice_position_offset: Quantity = Quantity(0., "m"),
                     time_bandwidth_product: float = 4.,
                     dummy_shots: int = None,
                     crusher_area: Quantity = Quantity(0., 'mT/m*ms'),
                     crusher_duration: Quantity = Quantity(0., 'mT/m*ms')) -> List[List[cmrseq.Sequence]]:

    mvenc = np.asarray([v.m_as("m/s") for v in venc_list])
    mvenc = mvenc[mvenc != 0.]
    min_venc = Quantity(np.min(np.abs(mvenc)), "m/s")

    venc_duration = cmrseq.seqdefs.velocity.bipolar(system_specs=system_specs, venc=min_venc, duration=venc_duration,
                                                    direction=np.array([1,0,0])).duration

    seq_list = []
    for venc, dir in zip(venc_list, venc_direction_list):

        seq_list.append(pc_gre(system_specs=system_specs,
                               matrix_size=matrix_size,
                               inplane_resolution=inplane_resolution,
                               slice_thickness=slice_thickness,
                               adc_duration=adc_duration,
                               flip_angle=flip_angle,
                               pulse_duration=pulse_duration,
                               repetition_time=repetition_time,
                               echo_time=echo_time,
                               venc=venc,
                               venc_direction=dir,
                               venc_duration=venc_duration,
                               slice_position_offset=slice_position_offset,
                               time_bandwidth_product=time_bandwidth_product,
                               dummy_shots=dummy_shots,
                               crusher_area=crusher_area,
                               crusher_duration=crusher_duration))

    return seq_list

se_m012_ssepi

se_m012_ssepi(
    echo_time: Quantity,
    slice_thickness: Quantity,
    slice_pos_offset: Quantity,
    field_of_view: Quantity,
    matrix_size: ndarray,
    b_vectors: Quantity,
    max_bval: Quantity = Quantity(450, "s/mm^2"),
    water_fat_shift="minimum",
    diff_raster_time: Quantity = Quantity(0.1, "ms"),
    spoiler_duration=Quantity(1.5, "ms"),
) -> List[cmrseq.Sequence]

Defines a spin echo single shot EPI with second order motion compensated diffusion Weighting

!!!WIP!!!

Source code in cmrseq/contrib/_cardiac_diffusion.py

def se_m012_ssepi(echo_time: Quantity, slice_thickness: Quantity, slice_pos_offset: Quantity,
                  field_of_view: Quantity, matrix_size: np.ndarray, b_vectors: Quantity,
                  max_bval: Quantity = Quantity(450, "s/mm^2"), water_fat_shift = "minimum",
                  diff_raster_time: Quantity = Quantity(0.1, "ms"),
                  spoiler_duration = Quantity(1.5, "ms"),
                 ) -> List[cmrseq.Sequence]:
    r"""Defines a spin echo single shot EPI with second order motion compensated
    diffusion Weighting

    !!!WIP!!!
    """
    system_specs_diff = cmrseq.SystemSpec(max_grad=Quantity(80, "mT/m"),
                                          max_slew=Quantity(100., "mT/m/ms"),
                                          b0=Quantity(1.5, "T"),
                                          rf_peak_power=Quantity(30, "uT"),
                                          grad_raster_time=diff_raster_time,
                                          rf_raster_time=diff_raster_time/10,
                                          adc_raster_time=diff_raster_time/100)

    system_specs_epi = cmrseq.SystemSpec(max_grad=Quantity(45, "mT/m"),
                                         max_slew=Quantity(80., "mT/m/ms"),
                                         b0=Quantity(1.5, "T"),
                                         rf_peak_power=Quantity(30, "uT"),
                                         grad_raster_time=diff_raster_time,
                                         rf_raster_time=diff_raster_time/10,
                                         adc_raster_time=diff_raster_time/100)

    excitation_seq = cmrseq.seqdefs.excitation.spectral_spatial_excitation(
                                        system_specs_diff, binomial_degree=3,
                                        total_flip_angle=Quantity(90, "degree"),
                                        slice_thickness=slice_thickness,
                                        time_bandwidth_product=4.5, chemical_shift=3.4)
    excite_rf_center = np.stack([excitation_seq.get_block(bn).rf_events[0]
                                 for bn in excitation_seq.blocks if "rf" in bn]).mean()

    refocus_seq = cmrseq.seqdefs.excitation.slice_selective_sinc_pulse(
                                    system_specs_diff, slice_thickness=slice_thickness,
                                    flip_angle=Quantity(180, "degree"),
                                    time_bandwidth_product=4,
                                    slice_position_offset=slice_pos_offset)

    refocus_seq.remove_block("slice_select_rewind_0")
    refocus_seq.rename_blocks(["rf_excitation_0", "slice_select_0"],
                              ["rf_refocus", "slice_select_refocus"])
    refocus_rf_center = refocus_seq.get_block("rf_refocus_0").rf_events[0]
    refocus_seq.shift_in_time(excite_rf_center - refocus_rf_center + echo_time / 2)


    diffusion_seq = cmrseq.seqdefs.diffusion.shortest_m012(system_specs_diff,
                                                           np.array([1, 0, 0]),
                                                           bvalues=max_bval, flip_decoding=True)
    diff_decode_start = diffusion_seq.get_block("diffusion_decode_0").tmin

    diffusion_seq.shift_in_time(- diff_decode_start + refocus_seq.end_time + spoiler_duration)
    # diffusion_seq._system_specs = system_specs

    max_epi_duration = 2 * (echo_time + excite_rf_center - (diffusion_seq.end_time - refocus_rf_center)) - Quantity(5, "ms")

    epi_seq = cmrseq.seqdefs.readout.single_shot_epi(system_specs_epi,
                                                     field_of_view=field_of_view,
                                                     matrix_size=matrix_size,
                                                     slope_sampling=True,
                                                     water_fat_shift=water_fat_shift,
                                                     max_total_duration=max_epi_duration,
                                                     partial_fourier_lines=0
                                                     )

    k_center_time = epi_seq.get_block(f"adc_{np.floor(matrix_size[1]/2).astype(int)}").adc_center
    epi_seq.shift_in_time(-k_center_time + echo_time + excite_rf_center)
    epi_seq._system_specs = system_specs_diff
    base_seq = deepcopy(excitation_seq)
    base_seq += refocus_seq
    base_seq += epi_seq

    zeta = diffusion_seq.get_block("diffusion_encode_0").rise_time
    lambda_ = diffusion_seq.get_block("diffusion_encode_0").flat_duration
    sequences = []
    for b_vec in b_vectors:
        b_value = Quantity(np.linalg.norm(b_vec.m_as("(s/mm^2)^(1/2)")) ** 2, "s/mm^2")
        if b_value > 0:
            diff_dir = (b_vec / np.sqrt(b_value)).m
        else:
            diff_dir = np.array([1., 0., 0.])
        diff_seq =  cmrseq.seqdefs.diffusion.m012(system_specs_diff, zeta, lambda_, diff_dir,
                                                  bvalue=b_value, flip_decoding=True)
        diff_seq.shift_in_time(- diff_decode_start + refocus_seq.end_time + spoiler_duration)
        # diff_seq._system_specs = system_specs

        temp_seq = deepcopy(base_seq)
        temp_seq += diff_seq
        # temp_seq._system_specs = system_specs_diff
        sequences.append(temp_seq)

    sequences_flip = [deepcopy(seq) for seq in sequences]
    for seq in sequences_flip:
        for bn in ["diffusion_encode_0", "diffusion_encode_1",
                   "diffusion_decode_0", "diffusion_decode_1"]:
            seq.get_block(bn).scale_gradients(-1)
    sequences.extend(sequences_flip)

    return sequences