In [2]:
%load_ext autoreload
%autoreload 1
#get our environment ready for data analysis times
%pylab inline
#import some os functionality so that we can be platform independent
import os
#import skimage components
from skimage.external import tifffile as tif #note that this can be achieved by using the
#skimage.io with the tifffile plugin.
#better filtering than skimage
from scipy import ndimage
#import pandas
import pandas as pd
from scipy.ndimage import measurements
from scipy.ndimage.filters import median_filter
import matplotlib.gridspec as gridspec #fancy subplot layout
from matplotlib.path import Path #Needed to create shapes
import matplotlib.patches as patches #needed so show shapes on top of graphs
%aimport peaks.gauss2d
Gauss2D = peaks.gauss2d.Gauss2D
%aimport peaks.utils
detrend = peaks.utils.detrend
%aimport peaks.peakfinder
PeakFinder = peaks.peakfinder.PeakFinder
%aimport dphplotting.mip
mip = dphplotting.mip.mip
In [4]:
PSF = tif.imread('testPSF.tif')
In [5]:
mip(PSF)
Out[5]:
In [6]:
def makeSlice(blob, width):
y,x,w,m = blob
return [slice(y-width//2,y+width//2), slice(x-width//2, x+width//2)]
In [7]:
set_cmap('gnuplot2')
width = 35
datapeaks = PeakFinder(PSF.max(0),1.68,'norot')
datapeaks.find_blobs()
blob = datapeaks.remove_edge_blobs(width//2)[-1]
myslice = makeSlice(blob,width)
#print(myslice)
fig,axes = mip(log(PSF[:,myslice[0],myslice[1]]), allaxes = True)
axes[3].matshow(log(PSF.max(0)))
y,x,_,_ = blob
rect = plt.Rectangle((x-width/2,y-width/2),width,width, color='r', linewidth=1, fill=False)
circ = plt.Circle((x, y), radius=width/2, color='r', linewidth=1, fill=False)
axes[3].add_patch(circ)
axes[3].add_patch(rect)
for ax in axes:
ax.axis('off')
fig.suptitle('PSF',y=1.05,fontsize=16)
#fig.savefig('mip_'+k.replace('.tif','.png'),dpi=300, transparent = True,bbox_inches='tight')
Out[7]:
In [8]:
#forward fit
def sliceMaker(y0, x0, width, ymax = np.inf, xmax = np.inf):
ystart = y0-width//2
xstart = x0-width//2
yend = ystart+width
xend = ystart+width
if ystart < 0:
ystart = 0
if xstart < 0 :
xstart = 0
if yend > ymax:
yend = ymax
if xend > xmax:
xend = xmax
return [slice(ystart,yend), slice(xstart, xend)]
def fitPeak(slices, stack, width, startingfit, **kwargs):
'''
Method to fit a peak in a stack.
The method will track the peak through the stack, assuming that moves are relatively small
from one slice to the next
Parameters
----------
slices : iterator
an iterator which dictates which slices to fit, should yeild integers only
stack : ndarray (3 dimensions)
width : integert
width of fitting window
startingfit : dict
fit coefficients
'''
#set up our variable to return
toreturn = []
#grab the starting fit parameters
popt_d = startingfit.copy()
y0 = int(round(popt_d['y0']))
x0 = int(round(popt_d['x0']))
if len(popt_d) == 6:
modeltype = 'norot'
elif len(popt_d) == 5:
modeltype = 'sym'
else:
modeltype = 'full'
for slic in slices:
#try to update the y0/x0 values
#if the fit has failed, these will be nan and the operation will raise a value error
#doing nothing leaves the values intact
#make the slice
myslice = sliceMaker(y0, x0, width, stack.shape[1],stack.shape[2])
#pull the starting values from it
ystart = myslice[0].start
xstart = myslice[1].start
#insert the z-slice number
myslice.insert(0,slic)
#set up the fit and perform it using last best params
fit = Gauss2D((PSF[myslice]))
#move our guess coefs back into the window
popt_d['x0']-=xstart
popt_d['y0']-=ystart
fit.optimize_params(popt_d, **kwargs)
#if there was an error performing the fit, try again without a guess
if fit.error:
fit.optimize_params(modeltype = modeltype, **kwargs)
#if there's still an error move on to the next fit
if not fit.error:
popt_d = fit.opt_params_dict()
popt_d['x0']+=xstart
popt_d['y0']+=ystart
popt_d['slice']=slic
toreturn.append(popt_d.copy())
y0 = int(round(popt_d['y0']))
x0 = int(round(popt_d['x0']))
else:
bad_fit = fit.opt_params_dict()
bad_fit['slice']=slic
toreturn.append(bad_fit.copy())
return toreturn
In [66]:
#find location of max intensity
my_max = unravel_index(PSF.argmax(),PSF.shape)
#find position
y0, x0 = my_max[1:3]
#make a slice
myslice = makeSlice([y0,x0,0,0],40)
ystart = myslice[0].start
xstart = myslice[1].start
myslice.insert(0,my_max[0])
#set up dataframe for looking for stuff
peakfits = []
#initial fit
max_z = Gauss2D(PSF[myslice])
max_z.optimize_params()
d = max_z.opt_params_dict()
d['slice']=my_max[0]
d['x0']+=xstart
d['y0']+=ystart
peakfits.append(d.copy())
d.pop('slice')
lastparams = max_z.opt_params
In [73]:
%timeit fitPeak(range(my_max[0]+1,PSF.shape[0]),PSF,40,d.copy(), quiet = True) + fitPeak(reversed(range(0, my_max[0])),PSF,40,d.copy(),quiet = True)
In [74]:
peakfits+=fitPeak(range(my_max[0]+1,PSF.shape[0]),PSF,40,d.copy(), quiet = True)\
+ fitPeak(reversed(range(0, my_max[0])),PSF,40,d.copy(),quiet = True)
In [68]:
peakfits_df = pd.DataFrame(peakfits).set_index('slice')
In [69]:
peakfits_df.sort(inplace=True)
In [13]:
peakfits_df
Out[13]:
In [14]:
peakfits_df.amp.plot()
Out[14]:
In [15]:
peakfits_df.plot(y=['sigma_x','sigma_y'])
peakfits_df.plot(y=['y0','x0'])
Out[15]:
In [77]:
peakfits_df
Out[77]:
In [76]:
peakfits_df[['x0','y0']] *= -1
In [16]:
peakfits_df[['x0','y0']].std()
Out[16]:
In [17]:
test = (((peakfits_df[['x0','y0']] - peakfits_df[['x0','y0']].mean()) / peakfits_df[['x0','y0']].std()).abs() < 3)
In [18]:
peakfits_df_filtered = peakfits_df[test.all(axis=1)]
In [19]:
peakfits_df_filtered.plot(y=['sigma_x','sigma_y'])
peakfits_df_filtered.plot(y=['x0','y0'])
Out[19]:
In [20]:
peakfits_df.assign(angle = lambda x : 90-rad2deg(arccos(x.rho))).dropna()
In [ ]:
peakfits_df.sort()
In [ ]:
%aimport peaks.stackanalysis
PSFStackAnalyzer = peaks.stackanalysis.PSFStackAnalyzer
In [25]:
myPSF = PSFStackAnalyzer(PSF)
In [75]:
%timeit myPSF.fitPeaks(40)
In [70]:
test = myPSF.fitPeaks(40)
In [71]:
test[0].amp.plot()
peakfits_df.amp.plot()
Out[71]:
In [ ]: