Merge branch 'master' into JS/apapt-tf-loading

Project.toml

@@ -1,64 +1,66 @@
 name = "MPIFiles"
 uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f"
 authors = ["Tobias Knopp <[email protected]>"]
-version = "0.13.0"
+version = "0.14.1"
 
 [deps]
 AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9"
 CodecZlib = "944b1d66-785c-5afd-91f1-9de20f533193"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
 DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FileIO = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
 Graphics = "a2bd30eb-e257-5431-a919-1863eab51364"
 HDF5 = "f67ccb44-e63f-5c2f-98bd-6dc0ccc4ba2f"
-ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
 ImageAxes = "2803e5a7-5153-5ecf-9a86-9b4c37f5f5ac"
+ImageMetadata = "bc367c6b-8a6b-528e-b4bd-a4b897500b49"
 Inflate = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575"
 Mmap = "a63ad114-7e13-5084-954f-fe012c677804"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 REPL = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 Tar = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulAngles = "6fb2a4bd-7999-5318-a3b2-8ad61056cd98"
+UnitfulParsableString = "06c00241-927a-4d5b-bb5e-6b5a2ada3567"
 
 [compat]
+Aqua = "0.6"
 AxisArrays = "0.3, 0.4"
 CodecZlib = "0.7"
-FileIO = "1.0"
+DelimitedFiles = "1"
+DocStringExtensions = "0.8, 0.9"
 FFTW = "1.3"
+FileIO = "1.0"
 Graphics = "0.4, 1.0"
 HDF5 = "0.14, 0.15, 0.16"
-ImageMetadata = "0.9"
 ImageAxes = "0.6"
+ImageMetadata = "0.9"
 Inflate = "0.1.2"
 Interpolations = "0.12, 0.13, 0.14"
-Reexport = "1.0"
 LinearOperatorCollection = "1.1"
 StableRNGs = "1.0.0"
-Unitful = "0.17, 1.2"
+Unitful = "1.13, 1.14, 1.15, 1.16, 1.17"
 UnitfulAngles = "0.6.1"
-DocStringExtensions = "0.8, 0.9"
+UnitfulParsableString = "0.1.6"
 julia = "1.5"
 
 [extras]
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
 HTTP = "cd3eb016-35fb-5094-929b-558a96fad6f3"
 LazyArtifacts = "4af54fe1-eca0-43a8-85a7-787d91b784e3"
 Scratch = "6c6a2e73-6563-6170-7368-637461726353"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "HTTP", "LazyArtifacts", "Scratch"]
+test = ["Aqua", "Test", "HTTP", "LazyArtifacts", "Scratch"]

README.md

@@ -6,4 +6,5 @@ MPIFiles.jl is a Julia package for loading and storing magnetic particle imaging
 [![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)
 [![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)
+[![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)

docs/Project.toml

@@ -1,2 +1,5 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+
+[compat]
+Documenter = "1"

docs/make.jl

@@ -13,9 +13,11 @@ makedocs(
         "Measurements" => "measurements.md",
         "System Matrices" => "systemmatrix.md",
         "Frequency Filter" => "frequencyFilter.md",
-        "Reconstruction Results" => "reconstruction.md"
+        "Reconstruction Results" => "reconstruction.md",
+        "Transfer Functions" => "transferfunction.md"
       #  "Positions" => "positions.md"
-    ]
+    ],
+    warnonly = [:missing_docs]
 )
 
 deploydocs(

docs/src/lowlevel.md

@@ -23,7 +23,7 @@ scannerFacility, scannerOperator, scannerManufacturer, scannerName, scannerTopol
 
 # acquisition parameters
 acqStartTime, acqNumFrames, acqNumAverages, acqGradient, acqOffsetField,
-acqNumPeriodsPerFrame, acqSize
+acqNumPeriodsPerFrame
 
 # drive-field parameters
 dfNumChannels, dfStrength, dfPhase, dfBaseFrequency, dfCustomWaveform, dfDivider,

docs/src/transferfunction.md

@@ -0,0 +1,69 @@
+# Transfer functions
+
+MPIFiles defines the type `TransferFunction` which represents the system-properties of a linear time-invariant system in frequency space ([see also](https://en.wikipedia.org/wiki/Transfer_function)). In the context of MPIFiles this usually includes the properties of the receive chain of an MPI system, but can also hold information about other parts of the signal chain. 
+
+## Basic construction and data access
+The `TransferFunction` object is constructed from samples in frequency space and offers two ways to access the underlying data. The first uses indexing to directly access the underlying data:
+
+```@setup tf
+using MPIFiles
+```
+```@repl tf
+    f = collect(range(0,1e6,step=1e3));
+    # TransferFunction of a simple lowpass filter
+    tf = TransferFunction(f, 1 ./ (1 .+ im*f/1e4 )) 
+    tf[1] # Value of tf at first frequency sample (f = 0 Hz)
+```
+
+The second method has to be called on the object itself and uses linear interpolation to access the tf at any frequency:
+```@repl tf
+    tf(0) # Value of tf at f = 0 Hz
+    tf(1e4) # Value of tf at f = 10000 Hz
+```
+
+A `TransferFunction` can have multiple channels, which adds a second index to both data access functions. Directly accessing multiple channels is also possible. The complex array used to construct the `TransferFunction` needs to have the shape [number of frequency samples, channels].
+
+```@repl tf
+tf = TransferFunction(f, [1 ./(1 .+im*f/1e4) 1 ./(1 .+im*f/1e3)])
+tf[11,1]
+tf[11,2]
+tf(1e4,1)
+tf(1e4,2)
+tf(1e4, [1,2])
+tf(1e4,:)
+```
+## Units 
+To attach units to the `TransferFunction` the keyword-parameter `units` can to be used to give a `Unitful` unit to every channel of the tf. This can be useful if the transfer function is not dimensionless but relates two physical quantities, e.g. voltage and current in the form of an impedance. All **interpolated** accesses to tf data then return a `Unitful.Quantity`.
+
+```@repl tf
+R = 1; # Ohm
+L = 10e-6; # Henry
+RL = TransferFunction(f, R .+ im*2pi*f*L, units=["V/A"])
+RL([0,100e3])
+```
+
+## Saving and loading
+A `TransferFunction` object can be saved to and loaded from a .h5 file.
+
+```@docs
+MPIFiles.save(filename::String, tf::TransferFunction)
+MPIFiles.TransferFunction(::String)
+```
+
+## Additional constructors
+
+In addition to the constructor taking a single (complex) array it is also possible to give two arrays representing amplitude and phase.
+
+It is also possible to construct a `TransferFunction` from the transfer function data included in an `MPIFile`.
+
+```@docs
+    MPIFiles.TransferFunction(freq::Vector{<:Real}, ampdata::Array{<:Real,N}, phasedata::Array{<:Real,N}) where N
+    MPIFiles.TransferFunction(::MPIFile)
+```
+
+## Other interesting functions
+```@docs
+  TransferFunction
+  combine(::TransferFunction, ::TransferFunction)
+```
+

src/Brukerfile.jl

@@ -540,17 +540,20 @@ function systemMatrix(b::BrukerFileCalib, rows, bgCorrection=true)
   nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1
 
   if numSubPeriods(b) > 1
+    rows_ = collect(rows)
     stepsize = numSubPeriods(b)
     for k=1:length(rows)
       freq = rows[k][1]
       rec = rows[k][2]
-      rows[k] = CartesianIndex((freq-1)*stepsize+1, rec)
+      rows_[k] = CartesianIndex((freq-1)*stepsize+1, rec)
     end
+  else
+    rows_ = rows
   end
 
   s = open(sfFilename)
   data = Mmap.mmap(s, Array{ComplexF64,3}, (prod(calibSize(b)), nFreq, rxNumChannels(b)))
-  S = data[:,rows]
+  S = data[:,rows_]
   close(s)
   rmul!(S, 1.0/acqNumAverages(b))
   return S

src/Conversion.jl

@@ -1,21 +1,21 @@
 # This file contains routines to generate MDF files
-export saveasMDF, loadDataset, loadMetadata, loadMetadataOnline, setparam!, compressCalibMDF
+export saveasMDF, loadDataset, loadMetadata, setparam!, compressCalibMDF
 export getFFdataPerPos, prepareAsMDFSingleMeasurement, convertCustomSF, blockAverage
 
 
 function setparam!(params::Dict, parameter, value)
-  if value != nothing
+  if !(isnothing(value) || ismissing(value))
     params[parameter] = value
   end
 end
 
 # we do not support all conversion possibilities
 function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocessing=false,
-                     numPeriodAverages=1, numPeriodGrouping=1, experimentNumber=nothing,
+                     numPeriodAverages=1, numPeriodGrouping=1, experimentNumber=missing,
                      fixDistortions=false, kargs...)
   params = loadMetadata(f)
 
-  if experimentNumber != nothing
+  if !ismissing(experimentNumber)
     params[:experimentNumber] = experimentNumber
   end
 
@@ -63,7 +63,7 @@ function loadDataset(f::MPIFile; frames=1:acqNumFrames(f), applyCalibPostprocess
           cat(zeros(Bool,acqNumFGFrames(f)),ones(Bool,acqNumBGFrames(f)), dims=1))
 
         snr = calibSNR(f)
-	      if snr == nothing
+	      if isnothing(snr)
           @info "calculate SNR"
           snr = calculateSystemMatrixSNR(f, data, numPeriodAverages=numPeriodAverages, numPeriodGrouping=numPeriodGrouping)
         end
@@ -151,8 +151,6 @@ function loadMeasParams(f, params = Dict{Symbol,Any}(); skipMeasData = false)
   return params
 end
 
-
-
 function appendBGDataset(params::Dict, filenameBG::String; kargs...)
   params = MPIFile(filenameBG) do fBG 
     appendBGDataset(params, fBG; kargs...)
@@ -171,7 +169,6 @@ function appendBGDataset(params::Dict, fBG::MPIFile; frames=1:acqNumFrames(fBG))
   return params
 end
 
-
 isConvertibleToMDF(f::MPIFile) = true
 function isConvertibleToMDF(f::BrukerFile)
   # check if raw data is consistent
@@ -196,16 +193,15 @@ function saveasMDF(filenameOut::String, filenameIn::String; kargs...)
 end
 
 function saveasMDF(filenameOut::String, f::MPIFile; filenameBG = nothing, enforceConversion=false, kargs...)
-
   # This is a hack. Needs to be fixed properly
-  if(haskey(kargs, :SNRThresh) || haskey(kargs, :sparsityTrafoRedFactor)) && calibSNR(f) != nothing
+  if(haskey(kargs, :SNRThresh) || haskey(kargs, :sparsityTrafoRedFactor)) && !isnothing(calibSNR(f))
     compressCalibMDF(filenameOut, f; kargs...)
     return
   end
 
   if enforceConversion || isConvertibleToMDF(f)
-    params = loadDataset(f;kargs...)
-    if filenameBG != nothing
+    params = loadDataset(f; kargs...)
+    if !isnothing(filenameBG)
       appendBGDataset(params, filenameBG)
     end
     saveasMDF(filenameOut, params)
@@ -746,7 +742,7 @@ function saveasMDF(file::HDF5.File, params::Dict{Symbol,Any})
   write(file, "/scanner/topology", get(params,:scannerTopology,"FFP"))
 
   # acquisition parameters
-  write(file, "/acquisition/numAverages",  params[:acqNumAverages])
+  write(file, "/acquisition/numAverages", params[:acqNumAverages])
   write(file, "/acquisition/numFrames", get(params,:acqNumFrames,1))
   write(file, "/acquisition/numPeriodsPerFrame", get(params,:acqNumPeriodsPerFrame,1))
   write(file, "/acquisition/startTime", "$( get(params,:acqStartTime, Dates.unix2datetime(time())) )")

src/DatasetStore/BrukerDatasetStore.jl

@@ -1,10 +1,11 @@
-export BrukerDatasetStore, BrukerStore, findBrukerFiles, getNewCalibNum
+export BrukerDatasetStore, findBrukerFiles, getNewCalibNum
 
 struct BrukerDatasetStore <: DatasetStore
   path::String
 end
 
-if ispath("/opt/mpidata")
+@static if ispath("/opt/mpidata")
+  export BrukerStore
   const BrukerStore = BrukerDatasetStore("/opt/mpidata")
 end
 

src/DatasetStore/MDFDatasetStore.jl

@@ -1,4 +1,4 @@
-export MDFDatasetStore, MDFStore, addStudy, getMDFStudyFolderName, calibdir, getMDFStore, getCalibStudy
+export MDFDatasetStore, addStudy, getMDFStudyFolderName, calibdir, getMDFStore, getCalibStudy
 
 struct MDFDatasetStore <: DatasetStore
   path::String
@@ -24,7 +24,8 @@ struct MDFDatasetStore <: DatasetStore
   end
 end
 
-if ispath("/opt/mpidata/Bruker")
+@static if ispath("/opt/mpidata/Bruker")
+  export MDFStore
   const MDFStore = MDFDatasetStore("/opt/data/Bruker")
 end
 path(e::Experiment{MDFDatasetStore}) = joinpath( path(e.study), string(e.num)*".mdf" )

src/FrequencyFilter.jl

@@ -45,7 +45,7 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0,
     filterFrequenciesByMinFreq!(freqIndices, f, minFreq)
   end
 
-  if maxFreq < nFreq
+  if maxFreq < rxBandwidth(f)
     filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq)
   end
 

src/MDF.jl

@@ -69,13 +69,7 @@ function h5haskey(filename, parameter)
   end
 end
 
-function getindex(f::MDFFile, parameter)
-  if haskey(f.file, parameter)
-    return read(f.file, parameter)
-  else
-    return nothing
-  end
-end
+getindex(f::MDFFile, parameter) = read(f.file, parameter)
 
 function getindex(f::MDFFile, parameter, default)
   #if !haskey(f.param_cache,parameter)
@@ -106,8 +100,8 @@ studyUuid(f::MDFFileV2) = @keyrequired UUID(f["/study/uuid"])
 studyDescription(f::MDFFileV1)::Union{String, Missing} = "n.a."
 studyDescription(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/study/description"]
 function studyTime(f::MDFFile)
-  t = f["/study/time"]
-  if typeof(t)==String
+  t = @keyoptional f["/study/time"]
+  if typeof(t) == String
    return DateTime(t)
   else
    return nothing
@@ -149,14 +143,15 @@ tracerSolute(f::MDFFileV2)::Union{Vector{String}, Missing} = @keyrequired _makeS
 tracerSolute(f::MDFFileV1)::Union{Vector{String}, Missing} = ["Fe"]
 function tracerInjectionTime(f::MDFFile)::Union{Vector{DateTime}, Nothing}
   p = typeof(f) <: MDFFileV1 ? "/tracer/time" : "/tracer/injectionTime"
-  if isnothing(f[p])
+  time = @keyoptional f[p]
+  if isnothing(time)
     return nothing
   end
 
-  if typeof(f[p]) == String
-    return [DateTime(f[p])]
+  if typeof(time) == String
+    return [DateTime(time)]
   else
-    return [DateTime(y) for y in f[p]]
+    return [DateTime(y) for y in time]
   end
 end
 #tracerInjectionTime(f::MDFFileV2) = DateTime( f["/tracer/injectionTime"] )
@@ -269,12 +264,12 @@ function rxHasTransferFunction(f::MDFFile)
   haskey(f.file, "/acquisition/receiver/transferFunction")
 end
 rxInductionFactor(f::MDFFileV1) = nothing
-rxInductionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/inductionFactor"]
+rxInductionFactor(f::MDFFileV2) = @keyoptional f["/acquisition/receiver/inductionFactor"]
 
 rxUnit(f::MDFFileV1)::Union{String, Missing} = "a.u."
 rxUnit(f::MDFFileV2)::Union{String, Missing} = @keyrequired f["/acquisition/receiver/unit"]
 rxDataConversionFactor(f::MDFFileV1) = repeat([1.0, 0.0], outer=(1,rxNumChannels(f)))
-rxDataConversionFactor(f::MDFFileV2) = @keyrequired f["/acquisition/receiver/dataConversionFactor"]
+rxDataConversionFactor(f::MDFFileV2) = @keyoptional f["/acquisition/receiver/dataConversionFactor"]
 
 # measurements
 function measData(f::MDFFileV1, frames=1:acqNumFrames(f), periods=1:acqNumPeriodsPerFrame(f),