Merge remote-tracking branch 'origin/master' into JA/patchBrukerIMT

.github/workflows/CompatHelper.yml

@@ -2,6 +2,7 @@ name: CompatHelper
 on:
   schedule:
     - cron: 0 0 * * *
+  workflow_dispatch:
 jobs:
   build:
     runs-on: ubuntu-latest

.github/workflows/breakage.yml

@@ -0,0 +1,76 @@
+name: Breakage 
+# Based on: https://github.com/JuliaSmoothOptimizers/LinearOperators.jl/blob/main/.github/workflows/Breakage.yml
+on:
+    pull_request:
+        branches:
+          - master
+    push:
+        branches:
+          - master
+        tags: '*'
+jobs:
+  break:
+    runs-on: ubuntu-latest
+    strategy:
+      fail-fast: false
+      matrix:
+        pkg: [
+          "MagneticParticleImaging/MPIReco.jl",
+        ]
+        pkgversion: [latest, stable]
+
+    steps:
+      - uses: actions/checkout@v4
+
+      # Install Julia
+      - uses: julia-actions/setup-julia@v2
+        with:
+          version: 1
+          arch: x64
+      - uses: actions/cache@v4
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+
+      # Breakage test
+      - name: 'Breakage of ${{ matrix.pkg }}, ${{ matrix.pkgversion }} version'
+        env:
+          URL: ${{ matrix.pkg }}
+          VERSION: ${{ matrix.pkgversion }}
+        run: |
+          set -v
+          mkdir -p ./pr
+          echo "${{ github.event.number }}" > ./pr/NR
+          git clone https://github.com/$URL
+          export PKG=$(echo $URL | cut -f2 -d/)
+          cd $PKG
+          if [ $VERSION == "stable" ]; then
+            TAG=$(git tag -l "v*" --sort=-creatordate | head -n1)
+            if [ -z "$TAG" ]; then
+              TAG="no_tag"
+            else
+              git checkout $TAG
+            fi
+          else
+            TAG=$VERSION
+          fi
+          export TAG
+          julia -e 'using Pkg;
+            PKG, TAG, VERSION = ENV["PKG"], ENV["TAG"], ENV["VERSION"]
+            joburl = joinpath(ENV["GITHUB_SERVER_URL"], ENV["GITHUB_REPOSITORY"], "actions/runs", ENV["GITHUB_RUN_ID"])
+            TAG == "no_tag" && error("Not tag for $VERSION")
+            pkg"activate .";
+            pkg"instantiate";
+            pkg"dev ../";
+            if TAG == "latest"
+              global TAG = chomp(read(`git rev-parse --short HEAD`, String))
+            end
+            pkg"build";
+            pkg"test";'

.github/workflows/ci.yml

@@ -44,6 +44,8 @@ jobs:
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
+        env:
+          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
         with:
           file: lcov.info
   docs:

Project.toml

@@ -1,7 +1,7 @@
 name = "MPIFiles"
 uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f"
 authors = ["Tobias Knopp <[email protected]>"]
-version = "0.15.3"
+version = "0.16.0"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"

README.md

@@ -4,7 +4,7 @@ MPIFiles.jl is a Julia package for loading and storing magnetic particle imaging
 
 [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://magneticparticleimaging.github.io/MPIFiles.jl/dev)
 [![DOI](http://joss.theoj.org/papers/10.21105/joss.01331/status.svg)](https://doi.org/10.21105/joss.01331)
-[![Build status](https://github.com/MagneticParticleImaging/MPIFiles.jl/workflows/CI/badge.svg)](https://github.com/MagneticParticleImaging/MPIFiles.jl/actions)
+[![Build status](https://github.com/MagneticParticleImaging/MPIMeasurements.jl/actions/workflows/ci.yml/badge.svg)](https://github.com/MagneticParticleImaging/MPIMeasurements.jl/actions/workflows/ci.yml)
 [![codecov.io](http://codecov.io/github/MagneticParticleImaging/MPIFiles.jl/coverage.svg?branch=master)](http://codecov.io/github/MagneticParticleImaging/MPIFiles.jl?branch=master)
 [![License](https://img.shields.io/github/license/MagneticParticleImaging/MPIFiles.jl?color=green&style=flat)](https://github.com/MagneticParticleImaging/MPIFiles.jl/blob/master/LICENSE)
-[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)
+[![Aqua QA](https://raw.githubusercontent.com/JuliaTesting/Aqua.jl/master/badge.svg)](https://github.com/JuliaTesting/Aqua.jl)

src/FrequencyFilter.jl

@@ -7,7 +7,7 @@ function getCalibSNR(f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
   idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq
 
   SNRAll = calibSNR(f)
-  if SNRAll != nothing
+  if !isnothing(SNRAll)
     SNR[idx,:] = SNRAll[:,:,1]
   end
 
@@ -36,11 +36,15 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   nReceivers = rxNumChannels(f)
   nPeriods = 1 #acqNumPeriodsPerFrame(f)
-  freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers))))
+  freqs = measIsFrequencySelection(f) ? measFrequencySelection(f) : 1:nFreq
 
-  filterFrequenciesBySelection!(freqIndices, f)
-  filterFrequenciesByChannel!(freqIndices, recChannels)
-
+  if numPeriodGrouping == 1
+    freqIndices = vec([CartesianIndex{2}(i, j) for i in freqs, j in recChannels])
+  else
+    freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers))))
+    filterFrequenciesBySelection!(freqIndices, f) 
+    filterFrequenciesByChannel!(freqIndices, recChannels)
+  end  
   if minFreq > 0
     filterFrequenciesByMinFreq!(freqIndices, f, minFreq)
   end
@@ -64,12 +68,11 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
   SNRAll = nothing
 
   if SNRThresh > 0 || numUsedFreqs > 0
-    SNR = zeros(nFreq, nReceivers)
-    idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq
+    SNR = zeros(Float64, nFreq, nReceivers)
 
     SNRAll = calibSNR(f)
     if !isnothing(SNRAll)
-      SNR[idx,:] = SNRAll[:,:,1]
+      SNR[freqs, :] = SNRAll[:, :, 1]
     end
   end
 
@@ -81,47 +84,46 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
     error("It is not possible to use SNRThresh and SNRFactorUsedFreq similtaneously")
   end
 
-
   if stepsize > 1
     filterFrequenciesByStepsize!(freqIndices, stepsize)
   end
 
-  return freqIndices
+  return collect(vec(freqIndices))
 end
 
 export filterFrequenciesBySelection!
-function filterFrequenciesBySelection!(indices::Vector{CartesianIndex{2}}, f::MPIFile)
+function filterFrequenciesBySelection!(indices, f::MPIFile)
   if measIsFrequencySelection(f)
     return filterFrequenciesBySelection!(indices, measFrequencySelection(f))
   end
 end
-filterFrequenciesBySelection!(indices::Vector{CartesianIndex{2}}, selection::Vector{Int64}) = filter!(x -> in(x[1], selection), indices)
+filterFrequenciesBySelection!(indices, selection::Vector{Int64}, sorted = issorted(selection)) = filter!(x -> sorted ? insorted(x[1], selection) : in(x[1], selection), indices)
 
 export filterFrequenciesByChannel!
-filterFrequenciesByChannel!(indices::Vector{CartesianIndex{2}}, channels) = filter!(x-> in(x[2], channels), indices)
+filterFrequenciesByChannel!(indices, channels, sorted = issorted(channels)) = filter!(x-> sorted ? insorted(x[2], channels) : in(x[2], channels), indices)
 
 export filterFrequenciesByMinFreq!
-function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, minFreq; numPeriodGrouping = 1)
+function filterFrequenciesByMinFreq!(indices, f::MPIFile, minFreq; numPeriodGrouping = 1)
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMinIdx!(indices, minIdx)
 end
-filterFrequenciesByMinIdx!(indices::Vector{CartesianIndex{2}}, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
+filterFrequenciesByMinIdx!(indices, minIdx) = minIdx > 0 ?  filter!(x -> x[1] > minIdx, indices) : indices 
 
 export filterFrequenciesByMaxFreq!
-function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq; numPeriodGrouping = 1)
+function filterFrequenciesByMaxFreq!(indices, f::MPIFile, maxFreq; numPeriodGrouping = 1)
   nFreq = rxNumFrequencies(f, numPeriodGrouping)
   maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1
   return filterFrequenciesByMaxIdx!(indices, maxIdx)
 end
-filterFrequenciesByMaxIdx!(indices::Vector{CartesianIndex{2}}, maxIdx) = filter!(x-> x[1] < maxIdx, indices)
+filterFrequenciesByMaxIdx!(indices, maxIdx) = filter!(x-> x[1] < maxIdx, indices)
 
 export filterFrequenciesByMaxMixingOrder!
-filterFrequenciesByMaxMixingOrder!(indices::Vector{CartesianIndex{2}}, maxMixingOrder, f::MPIFile) = filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, mixingFactors(f))
-filterFrequenciesByMaxMixingOrder!(indices::Vector{CartesianIndex{2}}, maxMixingOrder, mf::Matrix) = filter!(x-> mf[x[1], 4] <= maxMixingOrder, indices)
+filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, f::MPIFile) = filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, mixingFactors(f))
+filterFrequenciesByMaxMixingOrder!(indices, maxMixingOrder, mf::Matrix) = filter!(x-> mf[x[1], 4] <= maxMixingOrder, indices)
 
 export filterFrequenciesByNumSidebandFreqs!
-function filterFrequenciesByNumSidebandFreqs!(indices::Vector{CartesianIndex{2}}, numSidebandFreqs::Int64, f::MPIFile; numPeriodGrouping = 1)
+function filterFrequenciesByNumSidebandFreqs!(indices, numSidebandFreqs::Int64, f::MPIFile; numPeriodGrouping = 1)
   # https://en.wikipedia.org/wiki/Sideband
   # Because of period grouping we have more frequencies than in original data
 
@@ -140,17 +142,19 @@ function filterFrequenciesByNumSidebandFreqs!(indices::Vector{CartesianIndex{2}}
 end
 
 export filterFrequenciesBySNRThresh!
-function filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, f::MPIFile, snrthresh; numPeriodGrouping = 1)
+function filterFrequenciesBySNRThresh!(indices, f::MPIFile, snrthresh::T; numPeriodGrouping = 1) where T <: Real
   SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping)
   return filterFrequenciesBySNRThresh!(indices, snrthresh, SNR)
 end
-filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, SNRThresh, SNR) = filter!(x-> SNR[x] >= SNRThresh , indices)
+filterFrequenciesBySNRThresh!(indices, SNRThresh::T, SNR::Matrix) where T <: Real = filter!(x-> SNR[x] >= SNRThresh , indices)
+filterFrequenciesBySNRThresh!(indices, SNRThresh::T, SNR::Dict{CartesianIndex{2}, Float64}) where T <: Real = filter!(x-> get(SNR, x, 0.0) >= SNRThresh , indices)
+
 
 export filterFrequenciesByNumUsedFrequencies!
-function filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq)
+function filterFrequenciesByNumUsedFrequencies!(indices, f::MPIFile, maxFreq)
   error("Not implemented")
 end
-filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, maxIdx) = error("not implemented")
+filterFrequenciesByNumUsedFrequencies!(indices, maxIdx) = error("not implemented")
 #=
     numFreqsAlreadyFalse = sum(!freqMask)
     numFreqsFalse = round(Int,length(freqMask)* (1-numUsedFreqs))
@@ -168,12 +172,12 @@ filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, maxId
 =#
 
 export filterFrequenciesByStepsize!
-function filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsize)
+function filterFrequenciesByStepsize!(indices, stepsize)
   stepIndices = 1:stepsize:maximum(map(x -> x[1], indices))
-  filter!(x -> in(x[1], stepIndices), indices)
+  filter!(x -> insorted(x[1], stepIndices), indices)
 end
 
-function sortFrequencies(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false)
+function sortFrequencies(indices, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false)
   if sortBySNR && !sortByMixFactors
     indices = sortFrequenciesBySNR(indices, f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
   elseif !sortBySNR && sortByMixFactors
@@ -185,11 +189,11 @@ function sortFrequencies(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeri
 end
 
 export sortFrequenciesBySNR
-function sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
+function sortFrequenciesBySNR(indices, f::MPIFile; numPeriodGrouping = 1, stepsize = 1)
   SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize)
   sortFrequenciesBySNR(indices, SNR)
 end
-sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, SNR::AbstractArray) = indices[reverse(sortperm(SNR[indices]),dims=1)]
+sortFrequenciesBySNR(indices, SNR::AbstractArray) = indices[reverse(sortperm(SNR[indices]),dims=1)]
 
 export sortFrequenciesByMixFactors
 function sortFrequenciesByMixFactors(indices, f::MPIFile; numPeriodGrouping = 1)
@@ -203,7 +207,7 @@ function sortFrequenciesByMixFactors(indices, f::MPIFile; numPeriodGrouping = 1)
   end
   sortFrequenciesByMixFactors(indices, mfNorm)
 end
-sortFrequenciesByMixFactors(indices::Vector{CartesianIndex{2}}, mfNorm::AbstractArray) = indices[sortperm(mfNorm[indices])]
+sortFrequenciesByMixFactors(indices, mfNorm::AbstractArray) = indices[sortperm(mfNorm[indices])]
 
 
 function rowsToSubsampledRows(f::MPIFile, rows)
@@ -212,7 +216,7 @@ function rowsToSubsampledRows(f::MPIFile, rows)
     tmp = Array{Union{Nothing, CartesianIndex{2}}}(undef, rxNumFrequencies(f), rxNumChannels(f))
     idxAvailable = measFrequencySelection(f)
     for (i, idx) in enumerate(idxAvailable)
-      for d=1:rxNumChannels(f)
+      for d=1:size(tmp, 2)
         tmp[idx, d] = CartesianIndex(i, d)
       end
     end

src/MDFCommon.jl

@@ -25,15 +25,18 @@ function measDataTDPeriods(f::Union{MDFFileV2, MDFv2InMemory}, periods=1:acqNumP
 end
 
 function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=true)
-  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) ||
-    !measIsFourierTransformed(f)
+  if !experimentHasMeasurement(f) || !measIsFourierTransformed(f)
     return nothing
   end
 
   rows_ = rowsToSubsampledRows(f, rows)
 
   data_ = measDataRaw(f)
 
+  if !measIsFastFrameAxis(f)
+    data_ = permutedims(data_, [4, 1, 2, 3])
+  end
+
   data_ = data_[:, rows_, :]
   data = reshape(data_, Val(2))
 
@@ -48,18 +51,25 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru
     subsamplingIndices_ = tmp[:, rows_, :]
     subsamplingIndices = reshape(subsamplingIndices_, Val(2))
 
-    for l=1:size(fgdata,2)
+    for l ∈ axes(fgdata, 2)
       dataBackTrafo[:,l] .= 0.0
       dataBackTrafo[subsamplingIndices[:,l],l] .= fgdata[:,l]
       dataBackTrafo[:,l] .= adjoint(B) * vec(dataBackTrafo[:,l])
     end
     fgdata = dataBackTrafo
   end
 
-  if bgCorrection # this assumes equidistant bg frames
+  if bgCorrection && length(measBGFrameIdx(f)) > 0 # this assumes equidistant bg frames
     @debug "Applying bg correction on system matrix (MDF)"
     bgdata = data[measBGFrameIdx(f),:]
-    blockLen = measBGFrameBlockLengths( invpermute!(deepcopy(measIsBGFrame(f)), measFramePermutation(f)) ) # Added deepcopy to be side-effect free in in-memory MDF
+
+    if measIsFramePermutation(f)
+      mask = invpermute!(deepcopy(measIsBGFrame(f)), measFramePermutation(f)) # Added deepcopy to be side-effect free in in-memory MDF
+    else
+      mask = deepcopy(measIsBGFrame(f)) # Added deepcopy to be side-effect free in in-memory MDF
+    end
+    blockLen = measBGFrameBlockLengths(mask)
+
     st = 1
     for j=1:length(blockLen)
       bgdata[st:st+blockLen[j]-1,:] .=
@@ -69,16 +79,18 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru
 
     bgdataInterp = interpolate(bgdata, (BSpline(Linear()), NoInterp()))
     # Cubic does not work for complex numbers
-    origIndex = measFramePermutation(f)
     M = size(fgdata,1)
     K = size(bgdata,1)
     N = M + K
+    origIndex = measIsFramePermutation(f) ? measFramePermutation(f) : 1:N
     for m=1:M
       alpha = (origIndex[m]-1)/(N-1)*(K-1)+1
       for k=1:size(fgdata,2)
         fgdata[m,k] -= bgdataInterp(alpha,k)
       end
     end
+  elseif bgCorrection && length(measBGFrameIdx(f)) == 0
+    @warn "Ignoring parameter `bgCorrection` since there are no background frames in the file."
   end
   return fgdata
 end

src/Measurements.jl

@@ -188,7 +188,7 @@ Supported keyword arguments:
 """
 function getMeasurements(f::MPIFile, neglectBGFrames=true;
       frames=neglectBGFrames ? (1:acqNumFGFrames(f)) : (1:acqNumFrames(f)),
-      bgCorrection=false, interpolateBG=false, tfCorrection=rxHasTransferFunction(f),
+      bgCorrection=false, bgFrames = 1:length(measBGFrameIdx(f)), interpolateBG=false, tfCorrection=rxHasTransferFunction(f),
       sortFrames=false, numAverages=1, numPeriodGrouping=1, kargs...)
 
   if neglectBGFrames
@@ -197,7 +197,7 @@ function getMeasurements(f::MPIFile, neglectBGFrames=true;
     data = getAveragedMeasurements(f; frames=idx[frames],
                                       numAverages=numAverages, kargs...)
 
-    idxBG = measBGFrameIdx(f)
+    idxBG = measBGFrameIdx(f)[bgFrames]
     hasBGFrames = length(idxBG) > 0
     if bgCorrection && !hasBGFrames
       @warn "Background correction was selected but there are no background frames in the file."

src/SystemMatrix.jl

@@ -19,7 +19,7 @@ function getSystemMatrix(f::MPIFile, frequencies=nothing;
     frequencies = filterFrequencies(f)
   end
 
-  if measIsFastFrameAxis(f) && measIsFourierTransformed(f)
+  if measIsFourierTransformed(f) && numPeriodGrouping == 1
     # This is the fast path if the SM lays on disk in the required format
     data = systemMatrix(f, frequencies, bgCorrection)
   else