diff --git a/src/MDFCommon.jl b/src/MDFCommon.jl index d749b15..427a824 100644 --- a/src/MDFCommon.jl +++ b/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 \ No newline at end of file diff --git a/src/SystemMatrix.jl b/src/SystemMatrix.jl index 0ec9f4e..9bcb22d 100644 --- a/src/SystemMatrix.jl +++ b/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 #=