PSF Fitting Test

Development of StackAnalyzer and related classes.


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


Populating the interactive namespace from numpy and matplotlib

In [4]:
PSF = tif.imread('testPSF.tif')

In [5]:
mip(PSF)


Out[5]:
(<matplotlib.figure.Figure at 0x102b217b8>,
 array([<matplotlib.axes._subplots.AxesSubplot object at 0x112ea5550>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1136706d8>,
        <matplotlib.axes._subplots.AxesSubplot object at 0x1136c1278>], dtype=object))

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]:
<matplotlib.text.Text at 0x11471ec88>
<matplotlib.figure.Figure at 0x114199240>

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)


1 loops, best of 3: 6.79 s per loop

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]:
amp offset sigma_x sigma_y x0 y0
slice
0 2.566925 109.684633 5.499640 9.904805 109.350162 99.786462
1 2.100198 108.576603 5.046779 16.541716 111.903155 104.934294
2 2.538099 110.009754 4.536449 11.272036 112.260666 104.351709
3 3.781665 110.804557 3.101141 10.217138 108.519851 105.472446
4 NaN NaN NaN NaN NaN NaN
5 NaN NaN NaN NaN NaN NaN
6 NaN NaN NaN NaN NaN NaN
7 NaN NaN NaN NaN NaN NaN
8 NaN NaN NaN NaN NaN NaN
9 NaN NaN NaN NaN NaN NaN
10 NaN NaN NaN NaN NaN NaN
11 NaN NaN NaN NaN NaN NaN
12 NaN NaN NaN NaN NaN NaN
13 NaN NaN NaN NaN NaN NaN
14 NaN NaN NaN NaN NaN NaN
15 NaN NaN NaN NaN NaN NaN
16 NaN NaN NaN NaN NaN NaN
17 NaN NaN NaN NaN NaN NaN
18 NaN NaN NaN NaN NaN NaN
19 NaN NaN NaN NaN NaN NaN
20 NaN NaN NaN NaN NaN NaN
21 NaN NaN NaN NaN NaN NaN
22 75.528906 57.480162 28.574622 25.454387 101.815951 99.021356
23 74.872473 60.282711 27.072357 22.836968 102.139563 99.125892
24 48.282252 86.557258 18.422779 17.286842 101.690110 99.192835
25 48.249555 94.770688 14.681794 13.633270 101.713145 98.879193
26 51.308968 96.965847 13.423922 12.391056 101.353158 98.855991
27 54.381971 99.594230 11.887906 11.152488 101.640588 98.814782
28 59.103236 102.023964 10.209443 9.804970 101.686206 98.856330
29 64.365063 102.268938 9.788064 9.227288 101.545753 98.877285
... ... ... ... ... ... ...
71 NaN NaN NaN NaN NaN NaN
72 NaN NaN NaN NaN NaN NaN
73 NaN NaN NaN NaN NaN NaN
74 NaN NaN NaN NaN NaN NaN
75 NaN NaN NaN NaN NaN NaN
76 57.889935 124.619516 2.888623 3.068994 99.531918 98.097346
77 46.794625 123.509526 3.253013 3.436115 99.281517 98.410855
78 36.251184 121.209680 3.580423 4.080789 99.091352 98.707792
79 30.137126 119.042338 3.954852 4.228009 99.322679 98.753047
80 30.162420 118.071642 4.279784 4.372407 98.901477 99.237506
81 28.301628 117.263369 4.169762 4.700976 98.578117 99.546949
82 23.242725 115.939692 4.849059 4.747162 98.975125 99.110425
83 23.286831 115.351641 4.811594 5.331294 99.472096 99.272161
84 23.410485 114.269689 5.059800 5.465686 99.398872 99.147351
85 21.375599 111.712063 5.128856 5.704025 98.456786 99.464749
86 19.687888 110.784739 5.992234 6.338373 98.892676 100.342437
87 19.145818 110.711646 5.666380 6.614547 97.888602 99.297559
88 16.714178 108.885284 6.524853 7.298967 98.096269 99.710828
89 17.439841 109.061699 6.280373 7.136874 98.328293 99.949196
90 16.930242 109.564516 6.946419 7.746323 98.229957 98.898101
91 16.706650 108.747874 7.523896 8.590451 98.682297 99.403406
92 12.619664 106.980207 9.970288 8.467009 98.540394 99.205013
93 12.366962 106.244699 8.935597 10.759295 98.476755 99.263104
94 11.190261 106.022839 10.635715 13.426878 98.036888 99.389104
95 11.689555 106.680741 10.163404 12.583265 97.937501 99.789088
96 12.317252 108.959820 5.636378 6.633117 98.303771 99.488755
97 12.511903 109.683583 4.402282 5.221743 98.483018 98.730089
98 11.710488 109.833204 4.540440 4.839100 98.911276 99.790293
99 9.369913 108.991228 4.655492 6.001445 98.827127 98.877037
100 9.396854 108.253157 4.795620 6.264721 98.764325 99.481612

101 rows × 6 columns


In [14]:
peakfits_df.amp.plot()


Out[14]:
<matplotlib.axes._subplots.AxesSubplot at 0x11472c6d8>

In [15]:
peakfits_df.plot(y=['sigma_x','sigma_y'])
peakfits_df.plot(y=['y0','x0'])


Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x1138634e0>

In [77]:
peakfits_df


Out[77]:
amp offset sigma_x sigma_y x0 y0
slice
0 2.566925 109.684633 5.499640 9.904805 -109.350162 -99.786462
1 2.100198 108.576603 5.046779 16.541716 -111.903155 -104.934294
2 2.538099 110.009754 4.536449 11.272036 -112.260666 -104.351709
3 3.781665 110.804557 3.101141 10.217138 -108.519851 -105.472446
4 NaN NaN NaN NaN NaN NaN
5 NaN NaN NaN NaN NaN NaN
6 NaN NaN NaN NaN NaN NaN
7 NaN NaN NaN NaN NaN NaN
8 NaN NaN NaN NaN NaN NaN
9 NaN NaN NaN NaN NaN NaN
10 NaN NaN NaN NaN NaN NaN
11 NaN NaN NaN NaN NaN NaN
12 NaN NaN NaN NaN NaN NaN
13 NaN NaN NaN NaN NaN NaN
14 NaN NaN NaN NaN NaN NaN
15 NaN NaN NaN NaN NaN NaN
16 NaN NaN NaN NaN NaN NaN
17 NaN NaN NaN NaN NaN NaN
18 NaN NaN NaN NaN NaN NaN
19 NaN NaN NaN NaN NaN NaN
20 NaN NaN NaN NaN NaN NaN
21 NaN NaN NaN NaN NaN NaN
22 75.528882 57.480186 -28.574617 -25.454382 -101.815952 -99.021356
23 74.872519 60.282664 -27.072368 -22.836977 -102.139563 -99.125892
24 48.282250 86.557260 -18.422778 -17.286842 -101.690110 -99.192835
25 48.249555 94.770688 -14.681794 -13.633270 -101.713145 -98.879193
26 51.308969 96.965846 -13.423922 -12.391057 -101.353158 -98.855991
27 54.381971 99.594229 -11.887906 -11.152488 -101.640588 -98.814782
28 59.103236 102.023964 -10.209443 -9.804970 -101.686206 -98.856330
29 64.365063 102.268938 -9.788064 -9.227288 -101.545753 -98.877285
... ... ... ... ... ... ...
71 NaN NaN NaN NaN NaN NaN
72 NaN NaN NaN NaN NaN NaN
73 NaN NaN NaN NaN NaN NaN
74 NaN NaN NaN NaN NaN NaN
75 NaN NaN NaN NaN NaN NaN
76 57.889935 124.619516 2.888623 3.068994 -99.531918 -98.097346
77 46.794625 123.509526 3.253013 3.436115 -99.281517 -98.410855
78 36.251184 121.209680 3.580423 4.080789 -99.091352 -98.707792
79 30.137126 119.042338 3.954852 4.228009 -99.322679 -98.753047
80 30.162420 118.071642 4.279784 4.372407 -98.901477 -99.237506
81 28.301628 117.263369 4.169762 4.700976 -98.578117 -99.546949
82 23.242725 115.939692 4.849059 4.747162 -98.975125 -99.110425
83 23.286831 115.351641 4.811594 5.331294 -99.472096 -99.272161
84 23.410485 114.269689 5.059800 5.465686 -99.398872 -99.147351
85 21.375599 111.712063 5.128856 5.704025 -98.456786 -99.464749
86 19.687888 110.784739 5.992234 6.338373 -98.892676 -100.342437
87 19.145818 110.711646 5.666380 6.614547 -97.888602 -99.297559
88 16.714178 108.885284 6.524853 7.298967 -98.096269 -99.710828
89 17.439841 109.061699 6.280373 7.136874 -98.328293 -99.949196
90 16.930242 109.564516 6.946419 7.746323 -98.229957 -98.898101
91 16.706650 108.747874 7.523896 8.590451 -98.682297 -99.403406
92 12.619664 106.980207 9.970288 8.467009 -98.540394 -99.205013
93 12.366962 106.244699 8.935597 10.759295 -98.476755 -99.263104
94 11.190261 106.022839 10.635715 13.426878 -98.036888 -99.389104
95 11.689555 106.680741 10.163404 12.583265 -97.937501 -99.789088
96 12.317252 108.959820 5.636378 6.633117 -98.303771 -99.488755
97 12.511903 109.683583 4.402282 5.221743 -98.483018 -98.730089
98 11.710488 109.833204 4.540440 4.839100 -98.911276 -99.790293
99 9.369913 108.991228 4.655492 6.001445 -98.827127 -98.877037
100 9.396854 108.253157 4.795620 6.264721 -98.764325 -99.481612

101 rows × 6 columns


In [76]:
peakfits_df[['x0','y0']] *= -1

In [16]:
peakfits_df[['x0','y0']].std()


Out[16]:
x0    2.592903
y0    1.222642
dtype: float64

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]:
<matplotlib.axes._subplots.AxesSubplot at 0x112e5e4e0>

In [20]:
peakfits_df.assign(angle = lambda x : 90-rad2deg(arccos(x.rho))).dropna()


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-20-e6a17ce5eab2> in <module>()
----> 1 peakfits_df.assign(angle = lambda x : 90-rad2deg(arccos(x.rho))).dropna()

//anaconda/lib/python3.4/site-packages/pandas/core/frame.py in assign(self, **kwargs)
   2303 
   2304             if callable(v):
-> 2305                 results[k] = v(data)
   2306             else:
   2307                 results[k] = v

<ipython-input-20-e6a17ce5eab2> in <lambda>(x)
----> 1 peakfits_df.assign(angle = lambda x : 90-rad2deg(arccos(x.rho))).dropna()

//anaconda/lib/python3.4/site-packages/pandas/core/generic.py in __getattr__(self, name)
   2148                 return self[name]
   2149             raise AttributeError("'%s' object has no attribute '%s'" %
-> 2150                                  (type(self).__name__, name))
   2151 
   2152     def __setattr__(self, name, value):

AttributeError: 'DataFrame' object has no attribute 'rho'

In [ ]:
peakfits_df.sort()

In [ ]:
%aimport peaks.stackanalysis
PSFStackAnalyzer = peaks.stackanalysis.PSFStackAnalyzer

In [25]:
myPSF = PSFStackAnalyzer(PSF)

In [75]:
%timeit myPSF.fitPeaks(40)


1 loops, best of 3: 6.81 s per loop

In [70]:
test = myPSF.fitPeaks(40)

In [71]:
test[0].amp.plot()
peakfits_df.amp.plot()


Out[71]:
<matplotlib.axes._subplots.AxesSubplot at 0x115f72208>

In [ ]: