Improve SM SNR calculation

[jusack]

Currently the code regarding the SNR calculation of SM measurements is a bit disorganized and features multiple calculation options (maximum instead of median, masking) just commented out.

What I think is especially problematic is the if clause in L116 where only for the frequency components where the SNR is above 3.5 the SNR is recalculated, but selecting only the positions with maximum signal. This causes the SNR for each component that reaches this threshhold of 3.5 to be even higher, causing a gap in the SNR values at 3.5

MPIFiles.jl/src/SystemMatrix.jl

Lines 81 to 134 in 206b708

	function calculateSystemMatrixSNRInner(S, SNR, J, R, K, N, gridSize::NTuple{D,Int}) where D
	for j=1:J # Periods
	for r=1:R # Channels
	# precalc on entire FOV
	for k=1:K # Frequencies
	SBG = S[(N+1):end, k, r, j]
	SFG = S[1:N, k, r, j]
	meanBG = mean(SBG)
	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] = !iszero(noise) ? signal / noise : 0
	end
	# generate mask representing signal region
	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])
	maxSFG = maximum(SFG)
	mask[ SFG .> maxSFG*0.90 ] .= true
	end
	#@info sum(mask)/length(mask)
	# calc SNR on mask
	for k=1:K # Frequencies
	SBG = S[(N+1):end, k, r, j]
	SFG = S[1:N, k, r, j]
	meanBG = mean(SBG)
	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] )
	κmax = maximum(κ)
	signal = median( κ[(κ .> κmax*0.9 ) ] )
	# signal = mean( κ[ mask .&& phaseMask ] )
	end
	SNR[k, r, j] = signal / noise
	end
	end
	end
	SNR[:, :, :] .= median(SNR,dims=3)
	return
	end

[tknopp]

yes, this discontinuity is really an issue. But I have not found a better idea here.
Remarks:

• in general taking the maximum is nice because it is independent of the SM-Fov, discretization etc. It is thus most robust, when the SNR is high enough.
• taking the maximum has, however, the issue when there are outliers. This is the reason, I take a ROI (90% of max) and then take the medium.
• at low SNR the maximum method does not work, since even regular noise has some regular "off pixels". In practice this means that one gets a bias line of an SNR of (I think) 3-10 even though there is no noise.