Skip to content

Calc sfnr #224

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mihaic merged 10 commits into from
May 12, 2017
Merged

Calc sfnr #224

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
10 commits
Select commit Hold shift + click to select a range
c002be1
Substantial edits to how noise is calculated. Some edits to how noise…
Mar 29, 2017
1734764
pep8 update
Mar 29, 2017
9120c62
Made _calc_fwhm less complex
Mar 29, 2017
9d24c09
Update for edits
Mar 29, 2017
aab55be
[reprsimil] Fix BRSA bugs (#197)
lcnature Mar 22, 2017
34a8b02
[dev] Fix LDFLAGS for LLVM on MacOS
mihaic Mar 23, 2017
9434611
Merge branch 'master' of https://github.com/IntelPNI/BrainIAK
Mar 29, 2017
d9c65cc
Merge branch 'master' of https://github.com/IntelPNI/BrainIAK
May 12, 2017
5c73f97
_calc_snr turned into _calc_sfnr with the necessary changes associate…
May 12, 2017
3b28226
Pep8 errors corrected
May 12, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
86 changes: 55 additions & 31 deletions brainiak/utils/fmrisim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -785,12 +785,13 @@ def _calc_fwhm(volume,
return fwhm


def _calc_snr(volume,
mask,
):
""" Calculate the SNR of a volume
def _calc_sfnr(volume,
mask,
):
""" Calculate the SFNR of a volume
This takes the middle of the volume and averages the signal within the
brain and compares to the std in the voxels outside the brain.
brain and compares to the temporal standard deviation in the voxels
outside the brain.

Parameters
----------
Expand All @@ -803,9 +804,15 @@ def _calc_snr(volume,
Returns
-------

float
The snr of the volume (how much greater in activity the brain is from
the background)
float

The sd of the temporal variability of brain voxels.



float

The sfnr of the volume (mean brain activity divided by
 temporal
variability in the average non brain voxels)



float

What is the max activity measured here, a point of reference for masking

"""

Expand All @@ -821,35 +828,45 @@ def _calc_snr(volume,
# then you would have a lower mean signal by default)
mean_signal = (slice_volume[slice_mask > 0]).mean()

# What is the noise in the non brain
background_noise = float(slice_volume[slice_mask == 0].std())
# What are the brain and non-brain voxels
brain_voxels = volume[mask[:, :, :, 0] > 0]
mask_voxels = volume[mask[:, :, :, 0] == 0]

# What is the max activation of this volume
max_activity = volume.max()
# What is the noise in the masked voxels
mask_noise = np.std(mask_voxels, 1).mean()

# What is the brain voxel by time activity.
brain_voxels = volume[mask[:, :, :, 0] > 0]
# Assume the mask noise is a combination of random + drift. Thus average
# all masked voxels and find the variation over time, this is the
# variance due to drift.
drift_noise = np.mean(mask_voxels, 0).std()

# Subtract the drift variance from total mask noise to find the system
# noise
background_noise = np.sqrt(mask_noise ** 2 - drift_noise ** 2)

# Find the noise to brain voxels
temporal_noise = np.std(brain_voxels, 1).mean()

# Convert temporal noise into percent signal change
temporal_noise = temporal_noise / mean_signal * 100

# Calculate snr
snr = mean_signal / background_noise
# Calculate sfnr
sfnr = mean_signal / background_noise

# Convert from memmap
snr = float(snr)
sfnr = float(sfnr)

# What is the max activation of this volume
max_activity = volume.max()

return temporal_noise, snr, max_activity
return temporal_noise, sfnr, max_activity


def _calc_autoregression(volume,
def _calc_temporal_noise(volume,
mask,
auto_reg_order=1,
):
""" Calculate the autoregressive sigma of the data.
""" Calculate the mix of autoregressive and drift noise.

Parameters
----------
Expand Down Expand Up @@ -928,9 +945,9 @@ def calc_noise(volume,
# Since you are deriving the 'true' values then you want your noise to
# be set to that level

# Calculate the temporal_noise noise and SNR of the volume
noise_dict['temporal_noise'], noise_dict['snr'], noise_dict[
'max_activity'] = _calc_snr(volume, mask)
# Calculate the temporal_noise noise and SFNR of the volume
noise_dict['temporal_noise'], noise_dict['sfnr'], noise_dict[
'max_activity'] = _calc_sfnr(volume, mask)

# Calculate the fwhm on a subset of volumes

Expand All @@ -954,7 +971,7 @@ def calc_noise(volume,
noise_dict['fwhm'] = np.mean(fwhm)

# Calculate the autoregressive and drift noise
auto_reg_sigma, drift_sigma = _calc_autoregression(volume, mask)
auto_reg_sigma, drift_sigma = _calc_temporal_noise(volume, mask)

# Calibrate for how sigma is originally calculated
auto_reg_sigma = auto_reg_sigma / noise_dict['temporal_noise']
Expand Down Expand Up @@ -1557,8 +1574,8 @@ def _noise_dict_update(noise_dict):
noise_dict['auto_reg_sigma'] = 0.45
if 'physiological_sigma' not in noise_dict:
noise_dict['physiological_sigma'] = 0.1
if 'snr' not in noise_dict:
noise_dict['snr'] = 30
if 'sfnr' not in noise_dict:
noise_dict['sfnr'] = 30
if 'max_activity' not in noise_dict:
noise_dict['max_activity'] = 1000
if 'voxel_size' not in noise_dict:
Expand Down Expand Up @@ -1645,7 +1662,7 @@ def generate_noise(dimensions,
# Create the base (this inverts the process to make the mask)
base = mask * noise_dict['max_activity']

# Set the amount of background based on the snr value
# Set the amount of background based on the SFNR value

# What are the midpoints to be extracted
mid_x_idx = int(np.ceil(base.shape[0] / 2))
Expand All @@ -1658,19 +1675,26 @@ def generate_noise(dimensions,
# What is the mean signal of the non masked voxels in this slice?
mean_signal = (slice_volume[slice_mask > 0]).mean()

# Set up the machine noise
noise_system = _generate_noise_system(dimensions_tr=dimensions_tr)

# What is the standard deviation of the background activity
system_sigma = mean_signal / noise_dict['snr']
# (N.B. You need to subtract the mean system noise from this number
# since system_sigma * noise_system.mean() will later be added to the base)
system_sigma = mean_signal / (noise_dict['sfnr'] - noise_system.mean())

# Set up the machine noise
noise_system = (_generate_noise_system(dimensions_tr=dimensions_tr) *
system_sigma)
# Increase the size of the system noise based on the SFNR
noise_system *= system_sigma

# Convert temporal noise (in percent) to real numbers
abs_change = noise_dict['temporal_noise'] * mean_signal / 100

# Sum up the noise of the brain
noise = base + (noise_temporal * abs_change) + noise_system

# Reject negative values (only happens outside of the brain
noise[noise < 0] = 0

return noise


Expand Down
8 changes: 4 additions & 4 deletions tests/utils/test_fmrisim.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -209,7 +209,7 @@ def test_generate_noise():
stimfunction_tr=stimfunction_tr,
tr_duration=tr_duration,
mask=mask,
noise_dict={'temporal_noise': 0, 'snr': 10000},
noise_dict={'temporal_noise': 0, 'sfnr': 10000},
)

temporal_noise = np.std(noise[mask[:,:,:,0]>0],1).mean()
Expand Down Expand Up @@ -274,7 +274,7 @@ def test_calc_noise():
nd_orig = {'auto_reg_sigma': 0.6,
'drift_sigma': 0.4,
'temporal_noise': 5,
'snr': 30,
'sfnr': 30,
'max_activity': 1000,
'fwhm': 4,
}
Expand Down Expand Up @@ -302,5 +302,5 @@ def test_calc_noise():
precision = abs(nd_calc['temporal_noise'] - nd_orig['temporal_noise'])
assert precision < 1, 'temporal_noise calculated incorrectly'

precision = abs(nd_calc['snr'] - nd_orig['snr'])
assert precision < 10, 'snr calculated incorrectly'
precision = abs(nd_calc['sfnr'] - nd_orig['sfnr'])
assert precision < 5, 'sfnr calculated incorrectly'