'get_full_kspace()' Underestimates Time Duration

[mzaiss]

Description

get_full_kspace() underestimates the total sequence time because rep_traj[:, 3] spans event_count - 1 intervals instead of the full event_count durations from rep.event_time. This causes a ~50% time shortfall in downstream applications (e.g., plotting).

Steps to Reproduce

1. Create a Sequence with rep.event_time.sum() = 0.010s per repetition.
2. Call seq.get_full_kspace().
3. Check rep_traj[-1, 3] - rep_traj[0, 3] vs. rep.event_time.sum().

Expected vs. Actual

• Expected: rep_traj[-1, 3] - rep_traj[0, 3] = rep.event_time.sum() (e.g., 0.010s); total matches get_duration() (e.g., 0.142829s).
• Actual: rep_traj[-1, 3] - rep_traj[0, 3] = 0.009s; total ~0.072829s (ratio: 0.51).

Evidence

Diagnostics from a 15-rep sequence:

• Rep 1: event_time sum = 0.010000s, plotted = 0.009000s.
• Total: get_duration() = 0.142829s, plotted = 0.072829s.

Suggested Fix

Modify get_full_kspace() to include the starting point, yielding event_count + 1 points:

def get_full_kspace(self):
    k_pos = torch.zeros(4, device=self.device)
    trajectory = []
    stored = torch.zeros(4, device=self.device)

    for rep in self:
        if rep.pulse.usage == PulseUsage.EXCIT:
            k_pos = stored
        elif rep.pulse.usage == PulseUsage.REFOC:
            k_pos = -k_pos
        elif rep.pulse.usage == PulseUsage.STORE:
            stored = k_pos

        grad_time = torch.cat([rep.gradm, rep.event_time[:, None]], 1)
        cum_traj = torch.cumsum(grad_time, dim=0)
        rep_traj = torch.cat([k_pos.unsqueeze(0), k_pos + cum_traj], dim=0)
        k_pos = rep_traj[-1, :]
        trajectory.append(rep_traj)

    return trajectory

Here the time difference:

[mzaiss]

to plot the k-space I also updated the function with an own timing calculation (can be updated when issue is resolved), but alos with a time_range option which is very helpful for debiggung:

def plot_kspace_trajectory(seq0,
                          figsize: tuple[float, float] = (5, 5),
                          plotting_dims: str = 'xy',
                          plot_timeline: bool = True,
                          time_range: tuple[float, float] | None = None):
   Plot the kspace trajectory produced by seq0 within an optional time range.

    Parameters
    ----------
    seq0 : Sequence
        Sequence object that provides get_full_kspace() method and iterable repetitions
    figsize : tuple[float, float], optional
        The size of the plotted matplotlib figure. Default: (5, 5)
    plotting_dims : str, optional
        String defining what is plotted on the x and y axis ('xy', 'zy', etc). Default: 'xy'
    plot_timeline : bool, optional
        Plot a second subfigure with kx, ky, kz components over time. Default: True
    time_range : tuple[float, float] | None, optional
        Time range (start, end) in seconds to plot. If None, plots the full timeline. Default: None
    """
    assert len(plotting_dims) == 2
    assert plotting_dims[0] in ['x', 'y', 'z']
    assert plotting_dims[1] in ['x', 'y', 'z']
    dim_map = {'x': 0, 'y': 1, 'z': 2}

    if time_range is not None:
        assert len(time_range) == 2 and time_range[0] <= time_range[1], \
            "time_range must be a tuple (start, end) with start <= end"

    kspace = seq0.get_full_kspace()
    adc_mask = [rep.adc_usage > 0 for rep in seq0]

    cmap = plt.get_cmap('rainbow')
    plt.figure(figsize=figsize)
    
    if plot_timeline:
        plt.subplot(211)
    
    # K-space trajectory plot with time filtering
    for i, (rep_traj, mask) in enumerate(zip(kspace, adc_mask)):
        times = torch.cat([torch.tensor([0.0], device=seq0[i].device), 
                          torch.cumsum(seq0[i].event_time, dim=0)])
        if i > 0:
            time_offset = sum(seq0[j].event_time.sum() for j in range(i))
            times = times + time_offset

        # Filter by time_range if specified
        if time_range is not None:
            mask_time = (times >= time_range[0]) & (times <= time_range[1])
            if not mask_time.any():  # Skip if no points in range
                continue
            rep_traj = rep_traj[mask_time[:-1]]  # Adjust for event endpoints
            mask = mask[mask_time[:-1]]
            if len(rep_traj) == 0:  # Skip if filtered to empty
                continue

        kx = rep_traj[:, dim_map[plotting_dims[0]]]
        ky = rep_traj[:, dim_map[plotting_dims[1]]]
        plt.plot(kx, ky, c=cmap(i / len(kspace)))
        plt.plot(kx[mask], ky[mask], 'r.')
        plt.plot(kx[~mask], ky[~mask], 'k.')
    plt.xlabel(f"$k_{plotting_dims[0]}$")
    plt.ylabel(f"$k_{plotting_dims[1]}$")
    plt.grid()

    if plot_timeline:
        plt.subplot(212)
        time_offset = 0.0
        
        for i, rep in enumerate(seq0):
            rep_traj = kspace[i]
            times = torch.cat([torch.tensor([0.0], device=rep.device), 
                              torch.cumsum(rep.event_time, dim=0)])
            absolute_times = times + time_offset
            kx = torch.cat([rep_traj[0:1, 0], rep_traj[:, 0]])
            ky = torch.cat([rep_traj[0:1, 1], rep_traj[:, 1]])
            kz = torch.cat([rep_traj[0:1, 2], rep_traj[:, 2]])

            # Filter by time_range if specified
            if time_range is not None:
                mask_time = (absolute_times >= time_range[0]) & (absolute_times <= time_range[1])
                if not mask_time.any():  # Skip if no points in range
                    time_offset = absolute_times[-1]
                    continue
                absolute_times = absolute_times[mask_time]
                kx = kx[mask_time]
                ky = ky[mask_time]
                kz = kz[mask_time]
                if len(absolute_times) == 0:  # Skip if filtered to empty
                    time_offset = absolute_times[-1]
                    continue

            if i == 0:
                plt.plot(absolute_times, kx, c='r', label="$k_x$")
                plt.plot(absolute_times, ky, c='g', label="$k_y$")
                plt.plot(absolute_times, kz, c='b', label="$k_z$")
            else:
                plt.plot(absolute_times, kx, c='r', label="_")
                plt.plot(absolute_times, ky, c='g', label="_")
                plt.plot(absolute_times, kz, c='b', label="_")
            
            time_offset = times[-1] + time_offset  # Update offset to full rep duration

        plt.xlabel("Time (s)")
        plt.ylabel("k-space Position")
        plt.legend()
        plt.grid()

    plt.tight_layout()
    plt.show()

[mzaiss]

@J-Endres , could this affect the diffusion weighting?

[J-Endres]

No, the diffusion calculations are not based on get_full_kspace()

I agree that for plotting, the trajectory should include the last time point as well (the same problem should be true for gradient moments).

Before changing this, I would like to collect all places that are influenced by this - comments are welcome!