In [1]:
#all necessary imports to run SAMUROI
%matplotlib qt4
import numpy
import scipy.signal
import matplotlib.pyplot as plt
import samuroi
from samuroi.plugins.baseline import bandstop, power_spectrum
from samuroi.plugins.baseline import linbleeched_deltaF
from samuroi.plugins.baseline import stdv_deltaF
from samuroi.plugins.stabilize import Stabilization
from samuroi.plugins.swc import load_swc
from samuroi.plugins.tif import load_tif
In [2]:
#prefix defines the filepath to your experiment of choice
prefix ='/Volumes/J/3d/170309/s1/fov1/AVG_flz4883.tif'
In [3]:
#creates data, a numpy object out of the multitif file. Data is the image sequence.
data = load_tif(prefix)
#alternatively, one can upload data from a numpy file, which is also an image sequence that may have been processed before
##data=numpy.load(prefix+'.npy')
In [4]:
#stabilization is optional. if possible, data will be modified into a stabilized version of data."""
try:
#the create stabilization tool and run it
app = Stabilization(data = data)
#correct the data by the stabilization transformations
data = app.apply(data)
except Exception as e:
#Print the reason why stabilization failed.
print e
In [5]:
#percentile based evaluation of background value.
#q values indicate which percentiles are used for background and foreground. prints foreground and background.
background=numpy.percentile(data,q=5)
foreground=numpy.percentile(data,q=98)
print background
#foreground is just to indicate brightest pixels in image to illustrate range
print foreground
In [6]:
#background subtraction:
data=data-background
In [7]:
#creates a projection image of the data image sequence which will used as a morphology imag. Here, numpy.max() is used for maximal intensity projections.
#different types of numpy operations can be applied here.
morphology = numpy.max(data, axis=-1)
In [9]:
morphology.max()
Out[9]:
In [10]:
#functions needed for filter with spectral analysis and display
def power_spectrum(data, fs):
#Calculate the power spectral density (psd or power spectrum) and corresponding frequencies of the data.
#Arguments: data: the video data of the recording.
#fs: the sampling frequency of the data (used to scale the x axis).
#Returns: tuple holding df, avgpsd
#df: the frequency values defined by fs.
#avgpsd: the power spectrum of the data.
#First each pixels time series gets fourier transformed. Then, the power spectrum is
#calculated by averaging over all pixel fourier transforms and taking the absolute value.
N = data.shape[-1] / 2
_dfft = numpy.fft.fft(data,axis = -1)
avgpower = (_dfft*numpy.conjugate(_dfft)).mean(axis = (0,1))[0:N].real
df = numpy.fft.fftfreq(n = data.shape[-1], d = 1./fs)[0:N]
return df,avgpower
def bandstop(data, fs, start, stop):
#Apply band stop filter on the data.
# Arguments: data, the video data of a recording.
#fs, the sampling frequency of the data.
#start, the lower frequency, where to start filtering.
#stop, the upper frequency, where to stop filtering.
#Returns: filtered, A numpy array with the same shape as data.
#but with frequency filter applied.
# filter out frequencies from 12-13 Hz where data is sampled on 50 Hz basis.
assert(start < stop)
nyq = 0.5 * fs
high = stop/nyq
low = start/nyq
order = 3
b, a = scipy.signal.butter(order, [low, high], btype='bandstop')
#zi = scipy.signal.lfiltic(b, a, y=[0.])
dataf = scipy.signal.lfilter(b, a, data)
return dataf
In [11]:
#Bandstop filtering the data is especially useful when data was acquired with a spinning disc confocal. 'fs'
#corresponds to sampling frequency in Hz, 'start' and 'stop' define the bandstop filtering window.
#filter frequencies have to be manually read out from the plot of the power spectrum and the filtered data in cell below.
# filter out frequencies from 18.95-19.05 Hz where data is sampled on 51.13 Hz basis.
filtered = bandstop(data, fs =20.83 , start =2.35,stop =2.45)
In [12]:
#power spectrum of the unfiltered data
df,psd = power_spectrum(data,fs =20.83)
#power spectrum of the unfiltered data
df,psdfiltered = power_spectrum(filtered,fs =20.83)
#plot power spectra of unfiltered and filtered data to check the bandstop
plt.figure()
plt.plot(df,psd, label = 'non-filtered')
plt.plot(df,psdfiltered, label = 'filtered')
plt.yscale('log')
plt.xlabel('frequency (Hz)')
plt.ylabel('Power Spectral Density')
plt.legend()
Out[12]:
In [13]:
#define the filtered data as data
data=filtered
In [14]:
#calculate mean and deltaF for data
mean = numpy.mean(data,axis = -1)
#different modes function deltaF 'linear_bleech', 'stdv' , 'median'"""
#data = linbleeched_deltaF(data)
##data = median_deltaF(data)
data = stdv_deltaF(data)
In [15]:
#show the gui
app = samuroi.SamuROIWindow(data = data,morphology=morphology)
app.show()
In [ ]:
#loads swc file from e.g. neutube into the gui
swc = load_swc('/Users/friedrichjohenning/Dropbox/aprojects/SAMUROI/examples/subcellular_examples/exampledendrite.swc')
#this sets the width of the branch ROI generated from the swc. Can be larger than the dendritic diameter as it will
#only incorporate unmasked pixels.
swc.radius *= 2
#will apply the swc to the GUI as a branchmask
app.segmentation.load_swc(swc)
In [ ]: