Trouble while importing Pulseq files

[Anais-A]

Hello,

I have been trying to simulate sequences importing Pulseq files.

I tested the files provided on the Koma GitHub:

• All were imported correctly
• All simulated
• The JEMRIS ones (except epi) and the ge one had unexpected results.

Regardless of the odd results, I tested some examples from the Pulseq repository (https://github.com/pulseq/pulseq/tree/master/tests):

• Most of them did not load into Koma
• Only gre.seq loaded and the simulated image was disappointing.

I pursued my test trying to import my own work-in-progress Pulseq file. It did not work directly. As I was using arbitrary shapes for gradients and not trapezoids, I thought that it might be a small issue linked to it. So I investigated and found that, when the gradient key was 0, an error was thrown because no gradient id was 0.
I performed a quick fix by excluding this case in Pulseq.jl (l.458):
elseif gradLibrary[i]["type"] == 'g' && i!=0 #Arbitrary gradient waveform
Now I can import my sequence and nothing seems to be missing.
However, this quick fix does nothing for the other examples.

Do you encounter similar issues?
Is there anything specific to do while generating the Pulseq file to make it compatible with Koma?

Best,
Anais

[cncastillo]

Hi Anais,

Thanks for all the information. This is very helpful.

I performed a quick fix by excluding this case in Pulseq.jl (l.458):
elseif gradLibrary[i]["type"] == 'g' && i!=0 #Arbitrary gradient waveform

Regardless of the odd results, I tested some examples from the Pulseq repository (https://github.com/pulseq/pulseq/tree/master/tests):
Most of them did not load into Koma

The reason could be the same for these two points, and it should be pretty easy to fix (basically what you proposed). I have experienced a similar problem with the fid.seq file, as it has no gradients and the reader assumed that there is at least one gradient/adc/rf in the file. I was sure this issue was already fixed, so that was our mistake. @beorostica, do you know what happened there?

I will try all of the files at https://github.com/pulseq/pulseq/tree/master/tests. To see if that fixes it. But if you experienced problems with a file that is not fid*.seq please let me know.

Only gre.seq loaded, and the simulated image was disappointing.

Regarding the simulation errors, I know the cause. It is related to #52. When gradient spoiling is applied for a Phantom with spins' positions in a finite grid (like brain_phantom2D), it generates replicas in k-space. So that is why you get blurring in the spoiler directions. A quick solution would be to increase the number of spins (spin density per voxel). This can be done by either:

1. Increasing the voxel size (should also affect the gradient spoiler area)
2. Increasing the number of spins. For brain_phantom2D, you need to reduce the ss parameter to 1 (refer to the docs )

The best solution would be to implement an EPG simulation method, but we will do this after speeding up, even more, the simulations (if you are curious read this)

For the spoiler-related problems, tell me if putting ss=1 helps. I am very limited with time until mid-March, but I will try to take a look during the weekend. Sorry for the bad experience 😞 and thanks for the patience!

Cheers!

[Anais-A]

Hi Carlos,

Thank you for the detailed feedback.

Of all the files in https://github.com/pulseq/pulseq/tree/master/tests, I can only load gre.seq in Koma.

I tried to create a new phantom to avoid the spoiler issue. I was able to launch the brain_phantom2D command, but I cannot find how to output it as an HDF5 file. Did I miss something?

Best,
Anais

[cncastillo]

Hi! If you are using the UI, you don't necessary need to write the phantom to use it, you can do the following in the terminal

obj_ui[] = brain_phantom2D(...)

That should plot the phantom in the UI with the updated number of spins. The same is true for the sequence with the variable seq_ui.

If you want to write the phantom in HDF5 format, right now we do not provide a function to do that as we are developing it in #184. Using the HDF5 you could write the field sample/resolution ([ $\Delta x$, $\Delta y$, $\Delta z$ ] in mm), sample/offset ([ $x$, $y$, $z$ ] in mm) and sample/data:

ρ =   data[1,:,:,:]
T1 =  1e-3 ./ data[2,:,:,:]  # R1 [ms] => T1 [s]
T2 =  1e-3 ./ data[3,:,:,:]  # R2 [ms] => T2* [s]
T2s = 1e-3 ./ data[4,:,:,:] # R2* [ms] => T2* [s]
Δw = data[5,:,:,:]              # off-resonance field map [rad/s]

Basically the opposite of

KomaMRI.jl/KomaMRIFiles/src/Phantom/JEMRIS.jl

Lines 21 to 52 in 30042f9

	function read_phantom_jemris(filename)
	# A(:,:,:,1)=Sample.M0;
	# I=find(Sample.T1); R1 =zeros(size(Sample.T1)); R1(I) =1./Sample.T1(I);
	# I=find(Sample.T2); R2 =zeros(size(Sample.T2)); R2(I) =1./Sample.T2(I);
	# I=find(Sample.T2S);R2S=zeros(size(Sample.T2S)); R2S(I)=1./Sample.T2S(I);
	# A(:,:,:,2)=R1; #1/T1
	# A(:,:,:,3)=R2; #1/T2
	# A(:,:,:,4)=R2S; #1/T2s
	# A(:,:,:,5)=Sample.DB;
	fid = HDF5.h5open(filename)
	data = read(fid["sample/data"])
	Δx = read(fid["sample/resolution"]) * 1e-3 #[m]
	offset = read(fid["sample/offset"]) * 1e-3 #[m]
	mask = data[1,:,:,:] .!= 0
	#Maps
	ρ = data[1,:,:,:]
	T1 = 1e-3 ./ data[2,:,:,:]
	T2 = 1e-3 ./ data[3,:,:,:]
	T2s = 1e-3 ./ data[4,:,:,:]
	Δw = data[5,:,:,:]
	#Positions
	X, Y, Z = size(ρ)
	FOVx = (X-1)*Δx[1] #[m]
	FOVy = (Y-1)*Δx[2] #[m]
	FOVz = (Z-1)*Δx[3] #[m]
	xx = reshape((-FOVx/2:Δx[1]:FOVx/2),:,1,1) #[(-FOVx/2:Δx[1]:FOVx/2)...;]
	yy = reshape((-FOVy/2:Δx[2]:FOVy/2),1,:,1) #[(-FOVy/2:Δx[2]:FOVy/2)...;;]
	zz = reshape((-FOVz/2:Δx[3]:FOVz/2),1,1,:) #[(-FOVz/2:Δx[3]:FOVz/2)...;;;]
	x = xx*1 .+ yy*0 .+ zz*0 .+ offset[1] #spin x coordinates
	y = xx*0 .+ yy*1 .+ zz*0 .+ offset[2] #spin y coordinates
	z = xx*0 .+ yy*0 .+ zz*1 .+ offset[3] #spin z coordinates
	v = 0 # m/s

Sorry for the delays with the fix, do you have a particular sequence file that really want to work so I can prioritize it?

Cheers!

[curtcorum]

@Anais-A @cncastillo

Anais and Carlos,

I can confirm that I see the same behavior for the tests at:

https://github.com/pulseq/pulseq/tree/master/tests

i.e. most do not load, and the gre.seq does not reconstruct well. I did not attempt the work around.

I see the same as well for the KomaMRI examples, i.e. band on the ge.seq (possibly this is correct for the seq?) and gre_JEMRIS.seq and radial_JEMRIS.seq do not turn out well for reconstruction.

The above was discovered while trying to debug the simulation of the example from:

https://github.com/jfnielsen/TOPPEpsdSourceCode/blob/UserGuide/v6/examples/write2DGRE.m

which is in a private repository from Jon Nielsen's @jfnielsen TOPPE GR scanner interpreter.

At first it seemed like it could be sys parameters were mismatched, but it now seems like it might be a more general issue, given the above discussion. I believe the gre2d.seq has been tested successfully on GE scanners with TOPPE, so in theory should simulate as well.

I will certainly share any results or insight I find on further seq files that run and/or seem problematic with KomaMRI.

We will figure it out!

Curt

[curtcorum]

It looks like the phantom sampling is definitely interacting with the sequence...

Here is https://github.com/pulseq/pulseq/blob/master/tests/gre.seq run with the default

obj_ui[] = brain_phantom2D()

Based on the periodicity direction it looks like the K-space is transposed compared to the reconstruction:

Using...

obj_ui[] = brain_phantom2D(;ss=2)

the reconstruction still has ghosts...

but the K-space shows fewer peaks of artifactual or unspoiled signal, suggesting the phantom interaction:

[curtcorum]

Here is what the ge.seq from KomaMRI/examples/1.sequences/Pulseq_create_seq at 128x128 looks like with:

obj_ui[] = brain_phantom2D(;ss=1)

and the K-space clearly has periodic artifact/echo:

Definitely there is some sort of Cartesian sampling issue...probably why the non-Cartesian sequences seem to work well.

Why the pulseq gre.seq test sequence does not work at 64x64 is still puzzling and probably a clue.

[curtcorum]

Read #52. It looks like a Z spoiler should be used with any 2d gre, until jittering the sample points can be implemented?

[curtcorum]

Please see pull request #314 for proposed fixes when the seq file has no gradients, or a shaped gradient.

[beorostica]

Thanks @curtcorum, I'll take a look into that!

[aTrotier]

I have a working 2D MP2RAGE pulseq sequence (validated on our 3T prisma) :

When I perform the simulation, I get the same kind of artefacts (my spoiler are along the readout) but I also get the same artefacts if I put the spoiler along the Z axis.

ss = 1 does not seems to work.

[aTrotier]

Nice !

[curtcorum]

@Anais-A is this the fix, or equivalent?

Many Thanks!!!

function get_adc_phase_compensation(seq)
  phase = ComplexF32[]
  for i = 1:length(seq)
      if is_ADC_on(seq[i])
          N = seq.ADC[i].N
          ϕ = seq.ADC[i].ϕ
          #aux = ones(N) .* exp(-1im*ϕ)
          aux = ones(N) .* exp(1im*ϕ)
          append!(phase, aux)
      end
  end
  return phase
end

[Anais-A]

@curtcorum Yes, it works on my end at least! (I recompiled my local Koma version)

[beorostica]

Thanks guys for all you efforts!
Please refer to pull request #315, which effectively solves the issue with the minus sign.
We'll soon release a new version of KomaMRI!

[cncastillo]

Hi! We realized that the phase compensation was correct!! We added a test to verify it.

The problem was the RF phase, that was read with a minus sign from the Pulseq (KomaMRIFiles), and then plotted with another minus sign🤦(KomaMRIPlots), so it looked fine.

The fix is in #315

There was a sneaky A' instead of a transpose(A) when plotting that conjugated the RF complex amplitude and therefore the -phase. We will upload a patch ASAP.

[beorostica]

@curtcorum, are you still experiencing issues even after adding the minus sign in your tests? Maybe this is another problem.

[curtcorum]

@beorostica I will try the #315 to make sure I did the correct thing. Thanks.

[curtcorum]

@beorostica

It seems like things are improved...

Here is gre.seq from pulseq/tests and ss=1 which looks reasonable. ss=4 also works with some ghosting:

pulseq/tests/ge.seq seems improved, but still has a lot of the ghosting from Cartesian issues, this is ss=1:

The gre2d.seq from TOPPE is improved, but still has issues. This is what I am wanting to work on:

It is looking like the rf spoiling is fixed! The remainder is the Cartesian sampling interacting with the Cartesian phantom? It looks like 64x64 is mostly clean, 100x100 or 128x128 messier and the 200x200 is messiest.

I will try gre2d at 64x64 and see how it looks.

[curtcorum]

Here is the gre2d with 64x64 resolution. Much better than before, but still the Cartesian issues or ?

% write2DGRE.m
%
% 2D spoiled GRE demo sequence for Pulseq on GE v1.0 User Guide.
% Similar to writeGradientEcho.m from Pulseq demoSeq folder.

% System/design parameters.
% Reduce gradients by 1/sqrt(3) to allow for oblique scans.
% Reduce slew a bit further to reduce PNS.
% sys = mr.opts('maxGrad', 50/sqrt(3), 'gradUnit','mT/m', ...
%               'maxSlew', 120/sqrt(3), 'slewUnit', 'T/m/s', ...
%               'rfDeadTime', 100e-6, ...
%               'rfRingdownTime', 60e-6, ...
%               'adcDeadTime', 40e-6, ...
%               'adcRasterTime', 2e-6, ...
%               'blockDurationRaster', 10e-6, ...
%               'B0', 3.0);

sys = mr.opts('MaxGrad', 28, 'GradUnit', 'mT/m', ...
    'MaxSlew', 150, 'SlewUnit', 'T/m/s', 'rfRingdownTime', 20e-6, ...
    'rfDeadTime', 100e-6, 'adcDeadTime', 10e-6);

% Create a new sequence object
seq = mr.Sequence(sys);             

% % Acquisition parameters
% fov = 240e-3; Nx = 200; Ny = 200;   % FOV and resolution
% sliceThickness = 3e-3;              % slice thickness (m)
% alpha = 6;                          % flip angle (degrees)
% TR = 12e-3;                         % repetition time TR (s)
% TE = 5e-3;                          % echo time TE (s) 
% dwell = 10e-6;                      % ADC sample time (s)
% rfSpoilingInc = 117;                % RF spoiling increment

% Acquisition parameters
fov=230e-3; Nx=64; Ny=64;           % Define FOV and resolution
sliceThickness=3e-3;                % slice
%TE=[7.38 9.84]*1e-3;               % give a vector here to have multiple TEs (e.g. for field mapping)
TE=4.3e-3;
TR=20.0e-3;                           % only a single value for now
T1=0.833;
alpha=acos(exp(-TR/T1)) / pi * 180; % flip angle (Ernst angle)
rfSpoilingInc=117;              % RF spoiling increment

% Create alpha-degree slice selection pulse and gradient
[rf, gz] = mr.makeSincPulse(alpha*pi/180, 'Duration', 3e-3, ...
    'SliceThickness', sliceThickness, 'apodization', 0.42, ...
    'timeBwProduct', 4, 'system', sys);

% Define other gradients and ADC events
deltak = 1/fov;
%gx = mr.makeTrapezoid('x', 'FlatArea', Nx*deltak, 'FlatTime', Nx*dwell, 'system', sys);
gx = mr.makeTrapezoid('x', 'FlatArea', Nx*deltak, 'FlatTime', 3.2e-3, 'system', sys); % 3.2e-3 for FlatTIme, CAC 240304 ***
adc = mr.makeAdc(Nx, 'Duration', gx.flatTime, 'Delay', gx.riseTime, 'system', sys);
gxPre = mr.makeTrapezoid('x', 'Area', -gx.area/2, 'Duration', 1e-3, 'system',sys);
gzReph = mr.makeTrapezoid('z', 'Area', -gz.area/2, 'Duration', 1e-3, 'system', sys);
phaseAreas = ((0:Ny-1)-Ny/2)*deltak;
gyPre = mr.makeTrapezoid('y', 'Area', max(abs(phaseAreas)), ...
    'Duration', mr.calcDuration(gxPre), 'system', sys);
peScales = phaseAreas/gyPre.area;

% gradient spoiling
gxSpoil = mr.makeTrapezoid('x', 'Area', 2*Nx*deltak, 'system', sys);
gzSpoil = mr.makeTrapezoid('z', 'Area', 4/sliceThickness, 'system', sys);

% Calculate timing
delayTE = ceil((TE - mr.calcDuration(gx)/2 - mr.calcDuration(gxPre) - ...
    gz.fallTime - gz.flatTime/2)/seq.gradRasterTime)*seq.gradRasterTime;
delayTR = ceil((TR - mr.calcDuration(gxSpoil, gzSpoil) ...
    - mr.calcDuration(gx)/2 - TE ...
    - gz.flatTime/2 - gz.riseTime)/seq.gradRasterTime)*seq.gradRasterTime;
assert(all(delayTE >= 0));
assert(all(delayTR >= mr.calcDuration(gxSpoil, gzSpoil)));

% Loop over phase encodes and define sequence blocks
% iY <= -10        Dummy shots to reach steady state
% -10 < iY <= 0    ADC is turned on and used for receive gain calibration on GE
% iY > 0           Image acquisition

nDummyShots = 32;
pislquant = 0;    % number of shots/ADC events used for receive gain calibration

rf_phase = 0;
rf_inc = 0;

for iY = (-nDummyShots-pislquant+1):Ny
    isDummyTR = iY <= -pislquant;

    % RF spoiling
    rf.phaseOffset = rf_phase/180*pi;
    adc.phaseOffset = rf_phase/180*pi;
    rf_inc = mod(rf_inc+rfSpoilingInc, 360.0);
    rf_phase = mod(rf_phase+rf_inc, 360.0);
    
    % Excitation block
    % Mark start of segment (block group) by adding label.
    % Subsequent blocks in block group are NOT labelled.
    seq.addBlock(rf, gz, mr.makeLabel('SET', 'TRID', 2-isDummyTR));

    % Slice-select refocus and readout prephasing
    % Set phase-encode gradients to zero while iY < 1
    pesc = (iY>0) * peScales(max(iY,1));  % phase-encode gradient scaling
    seq.addBlock(gxPre, mr.scaleGrad(gyPre, pesc), gzReph);
    seq.addBlock(mr.makeDelay(delayTE));

    % ADC
    if isDummyTR
        seq.addBlock(gx);
    else
        seq.addBlock(gx, adc);
    end

    % Spoil and PE rephasing, and TR delay
    seq.addBlock(gxSpoil, mr.scaleGrad(gyPre, -pesc), gzSpoil);
    seq.addBlock(mr.makeDelay(delayTR));
end

%% Check sequence timing
[ok, error_report] = seq.checkTiming;
if (ok)
    fprintf('Timing check passed successfully\n');
else
    fprintf('Timing check failed! Error listing follows:\n');
    fprintf([error_report{:}]);
    fprintf('\n');
end

%% Output for execution and plot
seq.setDefinition('FOV', [fov fov sliceThickness]);
seq.setDefinition('Name', 'gre2d');
seq.write('gre2d_sim.seq')       % Write to pulseq file

seq.plot('timeRange', [0 3]*TR);

%% Optional slow step, but useful for testing during development,
%% e.g., for the real TE, TR or for staying within slewrate limits  
%rep = seq.testReport;
%fprintf([rep{:}]);

[aTrotier]

With bloch simulation we need more 'spins' in each voxels. I suppose you can increase the phantom resolution by a factor at least 10.

I don't know how much spins are required for each voxels in order to get something clean
, it will depend on your spatial resolution.

[curtcorum]

@aTrotier

With bloch simulation we need more 'spins' in each voxels. I suppose you can increase the phantom resolution by a factor at least 10.

I don't know how much spins are required for each voxels in order to get something clean , it will depend on your spatial resolution.

It is essentially a Nyquist condition, I believe.

I am working on up-sampling the existing 2d phantom using Interpolations.jl for testing and will share what I find.

There may be more to it in that it interacts with spoiling, etc.

This would not be too difficult to add or extend the existing phantom options in KomaMRI.

There is also the:
BigBrain-MR: a new digital phantom with anatomically-realistic magnetic resonance properties at 100-µm resolution
https://zenodo.org/records/7432527

[curtcorum]

@aTrotier

Here are some results with "upsampling" the phantom2D.

Using gre.m from pulseq/tests and ss=1 looks pretty good (and slightly different with the adc phase fixes?)

There is definitely still an echo in the K-space data:

With obj_ui[] = brain_phantom2D(;us=2) is looks a fair bit cleaner:

and the K-space does not have a strong aliasing echo (although some of the tail could still be in).

Using obj_ui[] = brain_phantom2D(;us=3) does not seem to change things much so us=2-3 is probably good enough for this sequence.

The changes are on:

https://github.com/curtcorum/KomaMRI.jl/tree/upsample_phantoms

And I'm putting in a pull request: #319

[cncastillo]

We registered new versions of KomaMRIBase, KomaMRIFiles, KomaMRIPlots, and KomaMRICore to fix the phase issue. It took us a little bit to register because we were waiting for JuliaRegistries/General#102440. We will soon also merge #319 and register a KomaMRIBase 0.8.3 😄 .

[aTrotier]

hmm still have some issues for gre sequence with RF spoiling and obj_ui[] = brain_phantom2D(;us=2)

I don't have the ghosting with RF phase advance = 0

gre.zip

Can you reproduce the artefacts ?

[beorostica]

We need to go more in depth about this issue. I can replicate the same artifacts in my local machine:

[curtcorum]

I get it with us=2, and us=3 below, reduced somewhat.

[aTrotier]

In term of MR physics, RFSpoiling is not linked to the phantom sampling no ?

[curtcorum]

The phantom sampling is the cause of non-physical (for an actual continuous object) spurious echoes in the simulation, especially for spoiled steady state sequences.

So yes, it is "linked" in the sense that RF-spoiling can reduce (or possibly increase) the image artifact caused by spurious echoes created by the phantom being under-sampled relative to the readout...

Looking at the raw data and the log(|k-data|+1) you can see the spurious echoes at the beginning and end of the readout. (yellow arrows)

RF-spoiling will change the phase of the individual echoes, relative to the rest of the readout signal from the current TR excitation. When readouts are combined in reconstruction, there will be the possibility that the rf spoiling or other phase cycle will help to cancel out the effects of the spurious echoes, and move the position of the ghosts or reduce the relative intensity in the image. The phase of each spurious echo is changed by the RF-spoiling, but it is not eliminated.

[curtcorum]

@aTrotier

Here is a "thin" 3d phantom, with us=2: obj_ui[] = brain_phantom3D(; us=2, start_end=[179,181]);

The image is substantially lower artifact and raw data cleaner. The spurious echoes are reduced by the de-phasing over the z-direction, either from any spoiling or just T2* being closer to reality.

The current brain_phantom2D is really 2D so will not capture any Z-gradient spoiling, and the T2* of a "voxel" is nonphysically long. My suggestion is to use this type of thin 3D phantom with any spoiled gradient echo sequence.

@cncastillo @beorostica Probably the 2D phantom should be given some thickness of one voxel?

[curtcorum]

@aTrotier It is slow to display (actually not display) in KomaUI, but the following gives nice results:

julia> KomaUI()
[ Info: Listening on: 127.0.0.1:6188, thread id: 4
[ Info: Loaded `Scanner` to `sys_ui[]`
[ Info: Loaded `Sequence` to `seq_ui[]`
[ Info: Loaded `Phantom` to `obj_ui[]`
[ Info: Loaded `RawAcquisitionData` to `raw_ui[]`
[ Info: Loaded image to `img_ui[]`
┌ Info: 2 CUDA capable device(s).
│   (0*) = "Quadro RTX 8000"
└   (1 ) = "Quadro RTX 8000"
┌ Info: Currently using package versions
│   KomaMRI = "0.8.0"
│   KomaMRICore = "0.8.1"
│   KomaMRIFiles = "0.8.2"
└   KomaMRIPlots = "0.8.1"

julia> [ Info: Loading sequence gre_label_pypulseq.seq ...

julia> obj_ui[] = brain_phantom3D(; us=3, start_end=[179,181]); # this takes a while on my system

julia> ┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 8440443
│   time_points = 10651
└   adc_points = 4096
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:02:20
  simulated_blocks:  403
  rf_blocks:         198
  acq_samples:       4096
141.036470 seconds (50.74 M allocations: 3.249 GiB, 1.69% gc time, 28.23% compilation time: <1% of which was recompilation)
[ Info: Exporting to ISMRMRD file: /tmp/Koma_signal.mrd
[ Info: Running reconstruction ...

[aTrotier]

I try to mimic one of my SIEMENS custom sequence and I think the spoiler was along X :)

I will run the 3D on the cpu during the night my GPU (mem=16go) cannot handle that :)

2D_MP2RAGE_spoilX_RF117.seq.zip

I have to label the profile to extract the 2 Inversion Time images

[curtcorum]

See PR #342 for anisotropic phantom sampling.

[curtcorum]

@aTrotier

Using: #342 and obj_ui[]=brain_phantom2D(; us=[10,1]);

julia> [ Info: Loading sequence 2D_MP2RAGE_spoilX_RF117.seq ...

julia> obj_ui[]=brain_phantom2D(; us=[10,1]);
[ Info: setting ss=1 since us=[10, 1] defined

julia> ┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 1039640
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:01:10
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
 71.082497 seconds (57.23 M allocations: 3.596 GiB, 1.66% gc time, 57.32% compilation time: <1% of which was recompilation)
[ Info: Exporting to ISMRMRD file: /tmp/Koma_signal.mrd
[ Info: Running reconstruction ...

Does not seem very clean.

[curtcorum]

Here seems to replicate the obj_ui[]=brain_phantom2D(; us=10); result, which looks pretty clean:

julia> [ Info: Loading sequence 2D_MP2RAGE_spoilX_RF117.seq ...
┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 10396400
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:04:39
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
280.360841 seconds (58.06 M allocations: 3.629 GiB, 1.65% gc time, 15.42% compilation time: <1% of which was recompilation)

[curtcorum]

Already obj_ui[]=brain_phantom3D(; us=2, start_end=[179, 180]); is looking not too bad:

julia> ┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 1669160
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:01:24
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
 84.770149 seconds (58.04 M allocations: 3.628 GiB, 1.62% gc time, 50.54% compilation time: <1% of which was recompilation)
[ Info: Exporting to ISMRMRD file: /tmp/Koma_signal.mrd
[ Info: Running reconstruction ...

[curtcorum]

obj_ui[]=brain_phantom3D(; us=3, start_end=[179, 180]); looks as clean or slightly cleaner than obj_ui[]=brain_phantom2D(; us=10); and takes less time (only ~56% as many spins):

[curtcorum]

@aTrotier

obj_ui[]=brain_phantom2D(; us=[1, 10]); does very good, as good, or almost identical to obj_ui[]=brain_phantom2D(; us=10); and takes much less time (10% of spins):

julia> ┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 1039640
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:25
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
 25.931142 seconds (13.55 M allocations: 730.193 MiB, 2.15% gc time, 0.00% compilation time)

[curtcorum]

The critical 2D parameter seems to be usy >= 5, obj_ui[]=brain_phantom2D(; us=[1, 5]); for a modestly clean image:

julia> ┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 519820
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:15
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
 15.481724 seconds (13.50 M allocations: 728.550 MiB, 2.38% gc time, 0.00% compilation time)

us=[1, 7] looks even cleaner than us=[1, 10] ???

julia> obj_ui[]=brain_phantom2D(; us=[1, 7]);
[ Info: setting ss=1 since us=[1, 7] defined

julia> [ Info: Loading sequence 2D_MP2RAGE_spoilX_RF117.seq ...
┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│   koma_version = v"0.8.1"
│   sim_method = Bloch()
│   spins = 727748
│   time_points = 21203
└   adc_points = 8192
Progress: 100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:19
  simulated_blocks:  787
  rf_blocks:         391
  acq_samples:       8192
 19.052634 seconds (13.54 M allocations: 729.511 MiB, 2.60% gc time, 0.00% compilation time)