diff --git a/Project.toml b/Project.toml index b49a0607..163dab2d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MPIFiles" uuid = "371237a9-e6c1-5201-9adb-3d8cfa78fa9f" authors = ["Tobias Knopp "] -version = "0.12.8" +version = "0.13.0" [deps] AxisArrays = "39de3d68-74b9-583c-8d2d-e117c070f3a9" diff --git a/docs/src/lowlevel.md b/docs/src/lowlevel.md index 5399dee8..9af862f4 100644 --- a/docs/src/lowlevel.md +++ b/docs/src/lowlevel.md @@ -40,7 +40,7 @@ measIsBGFrame, measIsSpectralLeakageCorrected, measFramePermutation # calibrations calibSNR, calibFov, calibFovCenter, calibSize, calibOrder, calibPositions, -calibOffsetField, calibDeltaSampleSize, calibMethod, calibIsMeanderingGrid +calibOffsetFields, calibDeltaSampleSize, calibMethod, calibIsMeanderingGrid # reconstruction results recoData, recoFov, recoFovCenter, recoSize, recoOrder, recoPositions diff --git a/src/Brukerfile.jl b/src/Brukerfile.jl index e7cc1438..d01eda03 100644 --- a/src/Brukerfile.jl +++ b/src/Brukerfile.jl @@ -540,22 +540,17 @@ function systemMatrix(b::BrukerFileCalib, rows, bgCorrection=true) nFreq = div(rxNumSamplingPoints(b)*numSubPeriods(b),2)+1 if numSubPeriods(b) > 1 - rows_ = collect(rows) - NFreq = rxNumFrequencies(b) - NRx = rxNumChannels(b) stepsize = numSubPeriods(b) for k=1:length(rows) - freq = mod1(rows[k],NFreq) - rec = div(rows[k]-1,NFreq) - rows_[k] = (freq-1)*stepsize+1 + rec*nFreq + freq = rows[k][1] + rec = rows[k][2] + rows[k] = CartesianIndex((freq-1)*stepsize+1, rec) end - else - rows_ = rows end s = open(sfFilename) - data = Mmap.mmap(s, Array{ComplexF64,2}, (prod(calibSize(b)),nFreq*rxNumChannels(b))) - S = data[:,rows_] + data = Mmap.mmap(s, Array{ComplexF64,3}, (prod(calibSize(b)), nFreq, rxNumChannels(b))) + S = data[:,rows] close(s) rmul!(S, 1.0/acqNumAverages(b)) return S @@ -659,7 +654,7 @@ calibFovCenter(b::BrukerFile) = calibSize(b::BrukerFile) = [parse(Int64,s) for s in b["PVM_Matrix"]] calibOrder(b::BrukerFile) = "xyz" calibPositions(b::BrukerFile) = nothing -calibOffsetField(b::BrukerFile) = nothing +calibOffsetFields(b::BrukerFile) = nothing calibDeltaSampleSize(b::BrukerFile) = [0.0, 0.0, 0.0] #TODO calibMethod(b::BrukerFile) = "robot" calibIsMeanderingGrid(b::BrukerFile) = true diff --git a/src/Conversion.jl b/src/Conversion.jl index eb565a0c..bee5fee3 100644 --- a/src/Conversion.jl +++ b/src/Conversion.jl @@ -115,7 +115,7 @@ end function loadCalibParams(f, params = Dict{Symbol,Any}()) if experimentIsCalibration(f) for op in [:calibFov, :calibFovCenter, - :calibSize, :calibOrder, :calibPositions, :calibOffsetField, + :calibSize, :calibOrder, :calibPositions, :calibOffsetFields, :calibDeltaSampleSize, :calibMethod] setparam!(params, op, eval(op)(f)) end @@ -825,7 +825,7 @@ function saveasMDF(file::HDF5.File, params::Dict{Symbol,Any}) writeIfAvailable(file, "/calibration/size", params, :calibSize) writeIfAvailable(file, "/calibration/order", params, :calibOrder) writeIfAvailable(file, "/calibration/positions", params, :calibPositions) - writeIfAvailable(file, "/calibration/offsetField", params, :calibOffsetField) + writeIfAvailable(file, "/calibration/offsetFields", params, :calibOffsetFields) writeIfAvailable(file, "/calibration/deltaSampleSize", params, :calibDeltaSampleSize) writeIfAvailable(file, "/calibration/method", params, :calibMethod) if hasKeyAndValue(params, :calibIsMeanderingGrid) diff --git a/src/FrequencyFilter.jl b/src/FrequencyFilter.jl index ad715281..0bd1b4f8 100644 --- a/src/FrequencyFilter.jl +++ b/src/FrequencyFilter.jl @@ -1,4 +1,19 @@ -export filterFrequencies +export filterFrequencies, sortFrequencies + +function getCalibSNR(f::MPIFile; numPeriodGrouping = 1, stepsize = 1) + nFreq = rxNumFrequencies(f, numPeriodGrouping) + nReceivers = rxNumChannels(f) + SNR = zeros(nFreq, nReceivers) + idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq + + SNRAll = calibSNR(f) + if SNRAll != nothing + SNR[idx,:] = SNRAll[:,:,1] + end + + SNR = SNR[1:stepsize:nFreq,:,:] + return SNR +end """ filterFrequencies(f; kargs...) => Vector{Int64} @@ -15,66 +30,40 @@ Supported keyword arguments: * sortByMixFactors """ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0, - maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f), - sortBySNR=false, numUsedFreqs=-1, stepsize=1, maxMixingOrder=-1, - sortByMixFactors=false, numPeriodAverages=1, numPeriodGrouping=1) + maxFreq=rxBandwidth(f), recChannels=1:rxNumChannels(f),numUsedFreqs=-1, stepsize=1, + maxMixingOrder=-1, numPeriodAverages=1, numPeriodGrouping=1, numSidebandFreqs = -1) nFreq = rxNumFrequencies(f, numPeriodGrouping) nReceivers = rxNumChannels(f) nPeriods = 1 #acqNumPeriodsPerFrame(f) + freqIndices = collect(vec(CartesianIndices((nFreq, nReceivers)))) - minIdx = floor(Int, minFreq / rxBandwidth(f) * (nFreq-1) ) + 1 - maxIdx = ceil(Int, maxFreq / rxBandwidth(f) * (nFreq-1) ) + 1 - - freqMask = zeros(Bool,nFreq,nReceivers,nPeriods) - - freqMask[:,recChannels,:] .= true - - if measIsFrequencySelection(f) - freqMask[:,recChannels,:] .= false - idx = measFrequencySelection(f) - freqMask[idx,recChannels,:] .= true - else - freqMask[:,recChannels,:] .= true - end - + filterFrequenciesBySelection!(freqIndices, f) + filterFrequenciesByChannel!(freqIndices, recChannels) + if minFreq > 0 - freqMask[1:(minIdx),:,:] .= false + filterFrequenciesByMinFreq!(freqIndices, f, minFreq) end - - if maxIdx < nFreq - freqMask[(maxIdx):end,:,:] .= false + + if maxFreq < nFreq + filterFrequenciesByMaxFreq!(freqIndices, f, maxFreq) end if maxMixingOrder > 0 if numPeriodGrouping == 1 - mf = mixingFactors(f) - for l=1:size(mf,1) - if mf[l,4] > maxMixingOrder || mf[l,4] > maxMixingOrder - freqMask[(l-1)+1,recChannels] .= false - end - end - else # This is a hack until we have a general solution - - numPatches = div(acqNumPeriodsPerFrame(f), numPeriodAverages) - fBands = (collect(0:(rxNumFrequencies(f)-1)).-1) .* numPatches - freqMaskMO = zeros(Bool,nFreq,nReceivers,nPeriods) - - for i=1:length(fBands) - for y=-maxMixingOrder:maxMixingOrder - index = fBands[i]+y - if 1 <= index <= size(freqMask,1) - freqMaskMO[index,recChannels] .= true - end - end - end - freqMask .&= freqMaskMO + filterFrequenciesByMaxMixingOrder!(freqIndices, maxMixingOrder, f) + else + error("Can not apply max mixing order with a period grouping larger than one, found: $numPeriodGrouping") end end + if numSidebandFreqs > 0 && numPeriodGrouping > 1 + filterFrequenciesByNumSidebandFreqs!(freqIndices, numSidebandFreqs, f, numPeriodGrouping = numPeriodGrouping) + end + SNRAll = nothing - if SNRThresh > 0 || numUsedFreqs > 0 || sortBySNR + if SNRThresh > 0 || numUsedFreqs > 0 SNR = zeros(nFreq, nReceivers) idx = measIsFrequencySelection(f) ? measFrequencySelection(f) : idx = 1:nFreq @@ -84,15 +73,85 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0, end end - if SNRThresh > 0 && numUsedFreqs > 0 + if SNRThresh > 0 && numUsedFreqs <= 0 + filterFrequenciesBySNRThresh!(freqIndices, SNRThresh, SNR) + elseif numUsedFreqs > 0 && SNRThresh <= 0 + filterFrequenciesByNumUsedFrequencies!(freqIndices, numUsedFreqs) + elseif numUsedFreqs > 0 && SNRThresh > 0 error("It is not possible to use SNRThresh and SNRFactorUsedFreq similtaneously") end - if SNRThresh > 0 && SNRAll != nothing - freqMask[ findall(vec(SNR) .< SNRThresh) ] .= false + + if stepsize > 1 + filterFrequenciesByStepsize!(freqIndices, stepsize) end - if numUsedFreqs > 0 && SNRAll != nothing + return freqIndices +end + +export filterFrequenciesBySelection! +function filterFrequenciesBySelection!(indices::Vector{CartesianIndex{2}}, 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) + +export filterFrequenciesByChannel! +filterFrequenciesByChannel!(indices::Vector{CartesianIndex{2}}, channels) = filter!(x-> in(x[2], channels), indices) + +export filterFrequenciesByMinFreq! +function filterFrequenciesByMinFreq!(indices::Vector{CartesianIndex{2}}, 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 + +export filterFrequenciesByMaxFreq! +function filterFrequenciesByMaxFreq!(indices::Vector{CartesianIndex{2}}, 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) + +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) + +export filterFrequenciesByNumSidebandFreqs! +function filterFrequenciesByNumSidebandFreqs!(indices::Vector{CartesianIndex{2}}, numSidebandFreqs::Int64, f::MPIFile; numPeriodGrouping = 1) + # https://en.wikipedia.org/wiki/Sideband + # Because of period grouping we have more frequencies than in original data + + # Find "virtual" frequency indices of original frequencies + fBands = (collect(0:(rxNumFrequencies(f))).-1) .* numPeriodGrouping + + delete = Int64[] + for (i, cart) in enumerate(indices) + freq = cart[1] + # Check if there is no frequency that is at most numSideBandFreqs away from our original frequency + if !any(fBand -> abs(fBand - freq) <= numSidebandFreqs, fBands) + push!(delete, i) + end + end + deleteat!(indices, delete) +end + +export filterFrequenciesBySNRThresh! +function filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, f::MPIFile, snrthresh; numPeriodGrouping = 1) + SNR = getCalibSNR(f, numPeriodGrouping = numPeriodGrouping) + return filterFrequenciesBySNRThresh!(indices, snrthresh, SNR) +end +filterFrequenciesBySNRThresh!(indices::Vector{CartesianIndex{2}}, SNRThresh, SNR) = filter!(x-> SNR[x] >= SNRThresh , indices) + +export filterFrequenciesByNumUsedFrequencies! +function filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile, maxFreq) + error("Not implemented") +end +filterFrequenciesByNumUsedFrequencies!(indices::Vector{CartesianIndex{2}}, maxIdx) = error("not implemented") +#= numFreqsAlreadyFalse = sum(!freqMask) numFreqsFalse = round(Int,length(freqMask)* (1-numUsedFreqs)) S = sortperm(vec(SNR)) @@ -106,50 +165,62 @@ function filterFrequencies(f::MPIFile; SNRThresh=-1, minFreq=0, end l += 1 end +=# - end - - - if stepsize > 1 - freqStepsizeMask = zeros(Bool,nFreq, nReceivers, nPatches) - freqStepsizeMask[1:stepsize:nFreq,:,:] = freqMask[1:stepsize:nFreq,:,:] - freqMask = freqStepsizeMask - end - - freq = findall( vec(freqMask) ) - - if sortBySNR && !sortByMixFactors && SNRAll != nothing - SNR = vec(SNR[1:stepsize:nFreq,:,:]) +export filterFrequenciesByStepsize! +function filterFrequenciesByStepsize!(indices::Vector{CartesianIndex{2}}, stepsize) + stepIndices = 1:stepsize:maximum(map(x -> x[1], indices)) + filter!(x -> in(x[1], stepIndices), indices) +end - freq = freq[reverse(sortperm(SNR[freq]),dims=1)] +function sortFrequencies!(indices::Vector{CartesianIndex{2}}, f::MPIFile; numPeriodGrouping = 1, stepsize = 1, sortBySNR = false, sortByMixFactors = false) + if sortBySNR && !sortByMixFactors + indices = sortFrequenciesBySNR(indices, f, numPeriodGrouping = numPeriodGrouping, stepsize = stepsize) + elseif !sortBySNR && sortByMixFactors + indices = sortFrequenciesByMixFactors(indices, f, numPeriodGrouping = numPeriodGrouping) + else + error("Can not apply multiple sorting algorithms to frequencies") end + return indices +end - if !sortBySNR && sortByMixFactors - mfNorm = zeros(nFreq,nReceivers,nPeriods) - mf = mixingFactors(f) - for k=1:nFreq - mfNorm[k,:,:] = norm(mf[k,1:3]) - end +export sortFrequenciesBySNR +function sortFrequenciesBySNR(indices::Vector{CartesianIndex{2}}, 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)] - freq = freq[sortperm(mfNorm[freq])] +export sortFrequenciesByMixFactors +function sortFrequenciesByMixFactors(indices, f::MPIFile; numPeriodGrouping = 1) + nFreq = rxNumFrequencies(f, numPeriodGrouping) + nReceivers = rxNumChannels(f) + nPeriods = 1 #acqNumPeriodsPerFrame(f) + mfNorm = zeros(nFreq,nReceivers,nPeriods) + mf = mixingFactors(f) + for k=1:nFreq + mfNorm[k,:,:] = norm(mf[k,1:3]) end - - freq + sortFrequenciesByMixFactors(indices, mfNorm) end +sortFrequenciesByMixFactors(indices::Vector{CartesianIndex{2}}, mfNorm::AbstractArray) = indices[sortperm(mfNorm[indices])] function rowsToSubsampledRows(f::MPIFile, rows) if measIsFrequencySelection(f) # In this case we need to convert indices - tmp = zeros(Int64, rxNumFrequencies(f), rxNumChannels(f) ) + tmp = Array{Union{Nothing, CartesianIndex{2}}}(undef, rxNumFrequencies(f), rxNumChannels(f)) idxAvailable = measFrequencySelection(f) - for d=1:rxNumChannels(f) - tmp[idxAvailable, d] = (1:length(idxAvailable)) .+ (d-1)*length(idxAvailable) + for (i, idx) in enumerate(idxAvailable) + for d=1:rxNumChannels(f) + tmp[idx, d] = CartesianIndex(i, d) + end end - rows_ = vec(tmp)[rows] - if findfirst(x -> x == 0, rows_) != nothing + rows_ = tmp[rows] + if any(isnothing, rows_) @error "Indices applied to systemMatrix are not available in the file" end + return identity.(rows_) # Go from Vector{Union{Nothing, Index}} to Vector{Index} else rows_ = rows end diff --git a/src/IMT.jl b/src/IMT.jl index 2bacaabf..adb7e150 100644 --- a/src/IMT.jl +++ b/src/IMT.jl @@ -201,7 +201,7 @@ calibSize(f::IMTFile) = [div(size(f["/systemResponseFrequencies"],1),2), size(f["/systemResponseFrequencies"],2), size(f["/systemResponseFrequencies"],3)] calibOrder(f::IMTFile) = "xyz" -calibOffsetField(f::IMTFileCalib) = nothing +calibOffsetFields(f::IMTFileCalib) = nothing calibDeltaSampleSize(f::IMTFile) = [0.0,0.0,0.0] calibMethod(f::IMTFile) = "simulation" #calibIsMeanderingGrid(f::IMTFile) = Bool(f["/calibration/isMeanderingGrid", 0]) diff --git a/src/MDF.jl b/src/MDF.jl index f96e0032..642f11f7 100644 --- a/src/MDF.jl +++ b/src/MDF.jl @@ -419,7 +419,7 @@ calibFov(f::MDFFile) = @keyoptional f["/calibration/fieldOfView"] calibFovCenter(f::MDFFile) = @keyoptional f["/calibration/fieldOfViewCenter"] calibSize(f::MDFFile) = @keyoptional f["/calibration/size"] calibOrder(f::MDFFile) = @keyoptional f["/calibration/order"] -calibOffsetField(f::MDFFile) = @keyoptional f["/calibration/offsetField"] +calibOffsetFields(f::MDFFile) = @keyoptional f["/calibration/offsetFields"] calibDeltaSampleSize(f::MDFFile) = @keyoptional f["/calibration/deltaSampleSize",[0.0,0.0,0.0]] calibMethod(f::MDFFile) = @keyrequired f["/calibration/method"] calibIsMeanderingGrid(f::MDFFile) = @keyoptional Bool(f["/calibration/isMeanderingGrid", 0]) diff --git a/src/MDFCommon.jl b/src/MDFCommon.jl index 8fba3525..2c225356 100644 --- a/src/MDFCommon.jl +++ b/src/MDFCommon.jl @@ -34,9 +34,7 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru data_ = measDataRaw(f) - data_ = reshape(data_, size(data_,1), - size(data_,2)*size(data_,3), - size(data_,4))[:, rows_, :] + data_ = data_[:, rows_, :] data = reshape(data_, Val(2)) fgdata = data[measFGFrameIdx(f),:] @@ -47,9 +45,7 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru B = createLinearOperator(measSparsityTransformation(f), ComplexF32; shape=tuple(calibSize(f)...)) tmp = measSubsamplingIndices(f) - subsamplingIndices_ = reshape(tmp, size(tmp,1), - size(tmp,2)*size(tmp,3), - size(tmp,4))[:, rows_, :] + subsamplingIndices_ = tmp[:, rows_, :] subsamplingIndices = reshape(subsamplingIndices_, Val(2)) for l=1:size(fgdata,2) diff --git a/src/MDFInMemory.jl b/src/MDFInMemory.jl index c8c1b3a3..8eb7adb0 100644 --- a/src/MDFInMemory.jl +++ b/src/MDFInMemory.jl @@ -1303,10 +1303,6 @@ for (fieldname, fieldtype) in zip(fieldnames(MDFv2InMemory), fieldtypes(MDFv2InM end end - -# The MDF specification uses the plural -calibOffsetFields(f::MDFFile) = calibOffsetField(f) - # And some utility functions measIsCalibProcessed(mdf::MDFv2InMemory)::Union{Bool, Missing} = measIsFramePermutation(mdf) && measIsFourierTransformed(mdf) && diff --git a/src/MPIFiles.jl b/src/MPIFiles.jl index 3405d97c..9b0df446 100644 --- a/src/MPIFiles.jl +++ b/src/MPIFiles.jl @@ -90,7 +90,7 @@ export measData, measDataTDPeriods, measIsFourierTransformed, measIsTFCorrected, # calibrations export calibSNR, calibSnr, calibFov, calibFieldOfView, calibFovCenter, calibFieldOfViewCenter, calibSize, calibOrder, calibPositions, - calibOffsetField, calibDeltaSampleSize, + calibOffsetFields, calibDeltaSampleSize, calibMethod, calibIsMeanderingGrid # reconstruction results @@ -199,7 +199,7 @@ abstract type MPIFile end @mustimplement calibSize(f::MPIFile) @mustimplement calibOrder(f::MPIFile) @mustimplement calibPositions(f::MPIFile) -@mustimplement calibOffsetField(f::MPIFile) +@mustimplement calibOffsetFields(f::MPIFile) @mustimplement calibDeltaSampleSize(f::MPIFile) @mustimplement calibMethod(f::MPIFile) @mustimplement calibIsMeanderingGrid(f::MPIFile) diff --git a/src/Measurements.jl b/src/Measurements.jl index db49dd96..cd63dae4 100644 --- a/src/Measurements.jl +++ b/src/Measurements.jl @@ -344,7 +344,6 @@ function getMeasurementsFD(f::MPIFile, args...; if frequencies != nothing # here we merge frequencies and channels - data = reshape(data, size(data,1)*size(data,2), size(data,3), size(data,4)) data = data[frequencies, :, :] end diff --git a/src/MultiMPIFile.jl b/src/MultiMPIFile.jl index e682e6ca..acde4ab2 100644 --- a/src/MultiMPIFile.jl +++ b/src/MultiMPIFile.jl @@ -144,7 +144,7 @@ end # TODO: define functions for multi calibration data # if experimentIsCalibration(f) # for op in [:calibSNR, :calibFov, :calibFovCenter, -# :calibSize, :calibOrder, :calibPositions, :calibOffsetField, +# :calibSize, :calibOrder, :calibPositions, :calibOffsetFields, # :calibDeltaSampleSize, :calibMethod] # setparam!(params, string(op), eval(op)(f)) # end diff --git a/src/Positions/Positions.jl b/src/Positions/Positions.jl index 8c5af8f5..57b3ab6f 100644 --- a/src/Positions/Positions.jl +++ b/src/Positions/Positions.jl @@ -80,6 +80,8 @@ function range(grid::RegularGridPositions, dim::Int) return 1:1 end end +ndims(grid::RegularGridPositions) = length(grid.shape) +axes(grid::RegularGridPositions) = tuple([range(grid, i) for i in 1:ndims(grid)]...) RegularGridPositions(shape, fov, center) = RegularGridPositions(shape, fov, center, ones(Int,length(shape))) @@ -209,13 +211,18 @@ function getindex(grid::RegularGridPositions, idx::Vector{T}) where T<:Number return 0.5.*fieldOfView(grid) .* (-1 .+ (2 .* idx .- 1) ./ shape(grid)) .+ fieldOfViewCenter(grid) end -function posToIdxFloat(grid::RegularGridPositions,pos::Vector) +function getindex(grid::RegularGridPositions, idx::CartesianIndex) + return getindex(grid, LinearIndices(tuple(grid.shape...))[idx]) +end + +function posToIdxFloat(grid::RegularGridPositions, pos) idx = 0.5 .* (shape(grid) .* ((pos .- fieldOfViewCenter(grid)) ./ ( 0.5 .* fieldOfView(grid) ) .+ 1) .+ 1) + idx = [isnan(val) ? one(eltype(idx)) : val for val in idx] return idx end -function posToIdx(grid::RegularGridPositions,pos::Vector) +function posToIdx(grid::RegularGridPositions, pos) idx = round.(Int64, posToIdxFloat(grid,pos)) for d=1:length(idx) if grid.sign[d] == -1 @@ -225,7 +232,7 @@ function posToIdx(grid::RegularGridPositions,pos::Vector) return idx end -function posToLinIdx(grid::RegularGridPositions,pos::Vector) +function posToLinIdx(grid::RegularGridPositions, pos) return (LinearIndices(tuple(shape(grid)...)))[posToIdx(grid,pos)...] end diff --git a/test/Cartesian.jl b/test/Cartesian.jl index 05e26cb9..95636ee4 100644 --- a/test/Cartesian.jl +++ b/test/Cartesian.jl @@ -7,27 +7,27 @@ fnSMFDSP = joinpath(tmpdir,"mdf","systemMatrixCartesianSP.mdf") smTD = MPIFile(fnSMTD) @test typeof(smTD) <: MDFFileV2 - +freqs = collect(vec(CartesianIndices((10, 1)))) @test acqNumPeriodsPerFrame(smTD) == 6500 @test rxNumSamplingPoints(smTD) == 76 -@test size(getSystemMatrix(smTD,1:10)) == (81,6500*10) +@test size(getSystemMatrix(smTD,freqs)) == (81,6500*10) @test size(getMeasurements(smTD)) == (76, 2, 6500, 81) @test size(getMeasurements(smTD, numPeriodAverages=65)) == (76, 2, 100, 81) @test size(getMeasurements(smTD, numPeriodAverages=65, numPeriodGrouping=100)) == (7600, 2, 1, 81) -@test size(getMeasurementsFD(smTD, frequencies=1:10)) == (10, 6500, 81) +@test size(getMeasurementsFD(smTD, frequencies=freqs)) == (10, 6500, 81) saveasMDF(fnSMFDMP, smTD, numPeriodAverages=65, applyCalibPostprocessing=true) smFDMP = MPIFile(fnSMFDMP) @test acqNumPeriodsPerFrame(smFDMP) == 100 @test rxNumSamplingPoints(smFDMP) == 76 -@test size(getSystemMatrix(smFDMP,1:10)) == (81,10*100) +@test size(getSystemMatrix(smFDMP,freqs)) == (81,10*100) saveasMDF(fnSMFDSP, smTD, numPeriodAverages=65, applyCalibPostprocessing=true, numPeriodGrouping=100) smFDSP = MPIFile(fnSMFDSP) @test acqNumPeriodsPerFrame(smFDSP) == 1 @test rxNumSamplingPoints(smFDSP) == 7600 -@test size(getSystemMatrix(smFDSP,1:10)) == (81,10) +@test size(getSystemMatrix(smFDSP,freqs)) == (81,10) end \ No newline at end of file diff --git a/test/General.jl b/test/General.jl index 156adee5..9b24da47 100644 --- a/test/General.jl +++ b/test/General.jl @@ -110,7 +110,7 @@ @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10) - @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50) + @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50) end end @@ -141,7 +141,7 @@ @info "Test $sm" @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1) - @test size( systemMatrix(sm,1:10) ) == (1936,10) + @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (1936,10) @test size(rxTransferFunction(sm)) == (817, 3) @test rxHasTransferFunction(sm) == true @@ -164,9 +164,9 @@ @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,) #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working - @test size(getSystemMatrix(sm,1:10)) == (1936,10) - @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20) - @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),bgCorrection=true)) == (1936,10) # test on the data level if the conversion was successful SNRThresh = 2 freq = filterFrequencies(smBruker,SNRThresh=SNRThresh) @@ -174,7 +174,7 @@ S = getSystemMatrix(sm,frequencies=freq) relativeDeviation = zeros(Float32,length(freq)) - for f in 1:length(freq) + for f in 1:length(map(f->f[1], freq)) relativeDeviation[f] = norm(SBruker[:,f]-S[:,f])/norm(SBruker[:,f]) end # test if relative deviation for most of the frequency components is below 0.003 @@ -253,7 +253,7 @@ @info "Test $sm" @test size( systemMatrixWithBG(sm) ) == (67,52,3,1) - @test size( systemMatrix(sm,1:10) ) == (60,10) + @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (60,10) @test size( systemMatrix(sm) ) == (60,52,3,1) @test size(calibSNR(sm)) == (52, 3, 1) diff --git a/test/IMT.jl b/test/IMT.jl index 1395b129..2a150681 100644 --- a/test/IMT.jl +++ b/test/IMT.jl @@ -93,7 +93,7 @@ end @test calibSize(calibIMT) == [20; 20; 1] @test calibOrder(calibIMT) == "xyz" @test calibPositions(calibIMT) == nothing - @test calibOffsetField(calibIMT) == nothing + @test calibOffsetFields(calibIMT) == nothing @test calibDeltaSampleSize(calibIMT) == [0.0; 0.0; 0.0] @test calibMethod(calibIMT) == "simulation" diff --git a/test/MDFInMemory.jl b/test/MDFInMemory.jl index 8698cce8..e94f8d80 100644 --- a/test/MDFInMemory.jl +++ b/test/MDFInMemory.jl @@ -147,7 +147,7 @@ end @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:100, loadasreal=true)) == (1634,3,1,10) - @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50) + @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50) end @testset "Calibration" begin @@ -168,7 +168,7 @@ end @test typeof(sm) <: MDFFileV2 @test size( systemMatrixWithBG(sm) ) == (1959,817,3,1) - @test size( systemMatrix(sm,1:10) ) == (1936,10) + @test size( systemMatrix(sm,collect(vec(CartesianIndices((10, 1))))) ) == (1936,10) @test size(rxTransferFunction(sm)) == (817, 3) @test rxHasTransferFunction(sm) == true @@ -198,9 +198,9 @@ end @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,) #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working - @test size(getSystemMatrix(sm,1:10)) == (1936,10) - @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20) - @test size(getSystemMatrix(sm,1:10,bgCorrection=true)) == (1936,10) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),bgCorrection=true)) == (1936,10) # test on the data level if the conversion was successfull SNRThresh = 2 freq = filterFrequencies(smBruker,SNRThresh=SNRThresh) diff --git a/test/MDFv1.jl b/test/MDFv1.jl index 11f178ca..9afcc348 100644 --- a/test/MDFv1.jl +++ b/test/MDFv1.jl @@ -96,7 +96,7 @@ for mdf in (mdfv1,mdfv2) @test size(getMeasurementsFD(mdf, numAverages=10, frames=1:500, loadasreal=true)) == (1634,3,1,50) - @test size(getMeasurementsFD(mdf,frequencies=1:10, numAverages=10)) == (10,1,50) + @test size(getMeasurementsFD(mdf,frequencies=collect(vec(CartesianIndices((10, 1)))), numAverages=10)) == (10,1,50) end @@ -125,15 +125,15 @@ for sm in (smv1,smv2) @test calibSize(sm) == [44; 44; 1] @test calibOrder(sm) == "xyz" @test calibPositions(smv1) == nothing - @test calibOffsetField(smv1) == nothing + @test calibOffsetFields(smv1) == nothing @test calibDeltaSampleSize(sm) == [0.001; 0.001; 0.001] @test calibMethod(sm) == "robot" @test size(filterFrequencies(sm, SNRThresh = 5)) == (147,) #@test size(filterFrequencies(sm, numUsedFreqs = 100)) == (100,) # not working - @test size(getSystemMatrix(sm,1:10)) == (1936,10) - @test size(getSystemMatrix(sm,1:10,loadasreal=true)) == (1936,20) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))))) == (1936,10) + @test size(getSystemMatrix(sm,collect(vec(CartesianIndices((10, 1)))),loadasreal=true)) == (1936,20) end end diff --git a/test/Positions.jl b/test/Positions.jl index 53f97dcd..638a60cb 100644 --- a/test/Positions.jl +++ b/test/Positions.jl @@ -6,6 +6,8 @@ pospath = joinpath(tmpdir,"positions","Positions.h5") ctr = [0.0,0.0,0.0]Unitful.mm caG = RegularGridPositions(shp,fov,ctr) @test shape(caG) == shp + @test ndims(caG) == 3 + @test axes(caG) == ((-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm", (-1.0:1.0:1.0)u"mm") @test fieldOfView(caG) == fov @test fieldOfViewCenter(caG) == ctr @test_throws BoundsError caG[0] @@ -15,6 +17,8 @@ pospath = joinpath(tmpdir,"positions","Positions.h5") @test caG[3] == [1,-1,-1]Unitful.mm @test caG[4] == [-1,0,-1]Unitful.mm @test caG[27] == [1,1,1]Unitful.mm + @test getindex(caG, CartesianIndex(1,1,1)) == [-1,-1,-1]Unitful.mm + @test getindex(caG, CartesianIndex(3,3,3)) == [1,1,1]Unitful.mm h5open(pospath, "w") do file write(file, caG) end