Fix SNR calculation for single BG frame files (#116)

Co-authored-by: jonschumacher <[email protected]>

src/MDFCommon.jl

@@ -96,22 +96,32 @@ function systemMatrix(f::Union{MDFFileV2, MDFv2InMemory}, rows, bgCorrection=tru
 end
 
 function systemMatrixWithBG(f::Union{MDFFileV2, MDFv2InMemory})
-  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) || !measIsFourierTransformed(f)
+  if !experimentHasMeasurement(f) || !measIsFourierTransformed(f)
     return nothing
   end
 
   data_ = measDataRaw(f)
+
+  if !measIsFastFrameAxis(f)
+    data_ = permutedims(data_, [4, 1, 2, 3])
+  end
+
   data = data_[:, :, :, :]
   return data
 end
 
 # This is a special variant used for matrix compression
 function systemMatrixWithBG(f::Union{MDFFileV2, MDFv2InMemory}, freq)
-  if !experimentHasMeasurement(f) || !measIsFastFrameAxis(f) || !measIsFourierTransformed(f)
+  if !experimentHasMeasurement(f) || !measIsFourierTransformed(f)
     return nothing
   end
 
   data_ = measDataRaw(f)
+
+  if !measIsFastFrameAxis(f)
+    data_ = permutedims(data_, [4, 1, 2, 3])
+  end
+
   data = data_[:, freq, :, :]
   return data
 end

src/SystemMatrix.jl

@@ -79,51 +79,57 @@ function calculateSystemMatrixSNR(f::MPIFile, S::Array; numPeriodAverages=1, num
 end
 
 function calculateSystemMatrixSNRInner(S, SNR, J, R, K, N, gridSize::NTuple{D,Int}) where D
-  for j=1:J
-    for r=1:R
+  for j=1:J # Periods
+    for r=1:R # Channels
       # precalc on entire FOV
-      for k=1:K
-        SBG = S[(N+1):end,k,r,j]
-        SFG = S[1:N,k,r,j]
+      for k=1:K # Frequencies
+        SBG = S[(N+1):end, k, r, j]
+        SFG = S[1:N, k, r, j]
         meanBG = mean(SBG)
-        noise = mean(abs.(SBG .- meanBG))
+        noise = length(SBG) > 1 ? noise = mean(abs.(SBG .- meanBG)) : abs.(meanBG) # Prevent Inf due to substracting the same value
         signal = median( abs.(SFG.-meanBG) ) 
-        SNR[k,r,j] = signal / noise
+        SNR[k, r, j] = !iszero(noise) ? signal / noise : 0
       end
+
       # generate mask representing signal region 
-      idx = sortperm(vec(SNR[:,r,j]), rev=true)
+      idx = sortperm(vec(SNR[:, r, j]), rev=true)
       mask = zeros(Bool, N)
       for q = 1:20
-        SFG = abs.(S[1:N,idx[q],r,j])
+        SFG = abs.(S[1:N, idx[q], r, j])
         maxSFG = maximum(SFG)
         mask[ SFG .> maxSFG*0.90 ] .= true
       end
+
       #@info sum(mask)/length(mask)
       # calc SNR on mask
-      for k=1:K
-        SBG = S[(N+1):end,k,r,j]
-        SFG = S[1:N,k,r,j]
+      for k=1:K # Frequencies
+        SBG = S[(N+1):end, k, r, j]
+        SFG = S[1:N, k, r, j]
         meanBG = mean(SBG)
-        noise = mean(abs.(SBG .- meanBG))
-        κ = abs.(SFG.-meanBG) 
+        noise = length(SBG) > 1 ? noise = mean(abs.(SBG .- meanBG)) : abs.(meanBG) # Prevent Inf due to substracting the same value
+        κ = abs.(SFG .- meanBG) 
         #κ_ = mapwindow(median!, reshape(κ, gridSize), ntuple(d->5,D))
         #signal = maximum( κ_ ) 
         signal = median( κ[mask .== true] ) 
         #signal = maximum( κ[mask .== true] ) 
+
         if signal > 3.5*noise
-          phaseMaskA = angle.(SFG.-meanBG) .> 0
-          phaseMaskB = angle.(SFG.-meanBG) .<= 0
-	  phaseMask = sum(phaseMaskA) > sum(phaseMaskB) ? phaseMaskA : phaseMaskB
-	  #signal = maximum( κ[mask .== true] )
+          # phaseMaskA = angle.(SFG.-meanBG) .> 0
+          # phaseMaskB = angle.(SFG.-meanBG) .<= 0
+          # phaseMask = sum(phaseMaskA) > sum(phaseMaskB) ? phaseMaskA : phaseMaskB
+          # signal = maximum( κ[mask .== true] )
           κmax = maximum(κ)
-          signal = median( κ[(κ .> κmax*0.9 ) ] ) 
-          #signal = mean( κ[ mask .&&  phaseMask ] ) 
+          signal = median( κ[(κ .> κmax*0.9 ) ] )
+          # signal = mean( κ[ mask .&&  phaseMask ] ) 
         end
-        SNR[k,r,j] = signal / noise
+
+        SNR[k, r, j] = signal / noise
       end
     end
   end
-  SNR[:,:,:] .= median(SNR,dims=3)
+
+  SNR[:, :, :] .= median(SNR,dims=3)
+
   return
 end
 #=