Merge branch 'master' into JS/change-getindex

Project.toml

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

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

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

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)

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

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])

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])

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)