An error of 'DivideError: integer division' occurs in the internal functions of simulate when using a sequence with spiral acquisition

[jigna405]

What happened?

We are simulating various Pulseq sequences through KomaMRI, all of which work perfectly except for the sequence with spiral acquisition. The error occurs when simulation starts, indicating an attempt to perform an integer division where the divisor is 0. The division causing the error is happening in the function rem(x::Int64, y::Int64) in the int.jl file.

The test sequence being used is spiral.seq from Pulseq without modification. The first simulation was conducted using an AMD Ryzen 5 5500U with Radeon Graphics 2.10 GHz. Then, to verify if the error was due to hardware incompatibility, a second simulation was performed using an Nvidia RTX 3060 TI, which presented the same error.

Environment

OS x86_64-w64-mingw32
Julia 1.10.4
KomaMRIPlots 0.8.0
KomaMRIFiles 0.8.2
KomaMRI 0.8.2
KomaMRICore 0.8.0
KomaMRIBase 0.8.4

[cncastillo]

Hi @jigna405, thanks for reporting! I can confirm that I was able to reproduce this error.

It is throwing an error when generating the raw-data (ISMRMRD.jl). As we technically do not know how to reconstruct an image just with the sequence, we need some additional information (like the size Nx, Ny, Nz).

If that information is not available we guess it. This is a case in which the estimation is incorrect:

julia> seq.DEF["Nx"]
12000

julia> seq.DEF["Ny"]
0  #<-----------The problem

julia> seq.DEF["Nz"]
2 #<-----------Also the FatSat is confused with a new slice

You can avoid this guessing by specifying the dimension in the .seq. Specifically in the [DEFINITIONS] section. Try to specify

...
[DEFINITIONS]
...
Nx 100
Ny 100
Nz 1
...

We will definitely need to fix this. Or do the max(estimation, 1). As this also gave problems in #324.

[curtcorum]

It is interesting. The spiral.seq file used to work:

#311 (comment)

[cncastillo]

Yeah I think that it is because I changed KomaMRIFiles/src/Sequence/Pulseq.jl, for performance reasons, from

if !haskey(seq.DEF,"Nx")
    Nx = maximum(seq.ADC.N)
    RF_ex = (get_flip_angles(seq) .<= 90.01) .* is_RF_on.(seq)
    Nz = max(length(unique(seq.RF[RF_ex].Δf)), 1)
    Ny = sum(is_ADC_on.(seq)) / Nz |> x->floor(Int,x)

    seq.DEF["Nx"] = Nx  #Number of samples per ADC
    seq.DEF["Ny"] = Ny  #Number of ADC events
    seq.DEF["Nz"] = Nz  #Number of unique RF frequencies, in a 3D acquisition this should not work
end

to

# Guessing recon dimensions
seq.DEF["Nx"] = get(seq.DEF, "Nx", maximum(adc.N for adc = seq.ADC))
seq.DEF["Nz"] = get(seq.DEF, "Nz", length(unique(seq.RF.Δf)))
seq.DEF["Ny"] = get(seq.DEF, "Ny", sum(map(is_ADC_on, seq)) ÷ seq.DEF["Nz"])

6176d03#diff-d10352b8f4adb4ab29f5145c81811c3439646352a93598ae4bcd207a3e39e3b5L487-R475

That can generate some 0's.

[curtcorum]

It is a good one to add to a test!

[cncastillo]

Agreed, all of the provided sequences should work. We started doing that in #284. Related #367.