In [ ]:

    

#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 [ ]:

    

#prefix defines the filepath to your experiment of choice


prefix ='/Volumes/J/3d/170309/s1/fov1/SUM_flz4883.tif'


    

In [ ]:

    

#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 [ ]:

    

#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 [ ]:

    

#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 [ ]:

    

#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 [ ]:

    

#background subtraction: 
data=data-background


    

In [ ]:

    

#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 [ ]:

    

#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 =1,stop =2)


    

In [ ]:

    

#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()


    

In [ ]:

    

#define the filtered data as data
data=filtered


    

In [ ]:

    

#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 [ ]:

    

#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 [ ]: