<jupyter_start><jupyter_text>Similar to fiducial drift correction, 3D imaging based on astigmatism is implemented in B-Store in separate parts:1. the `CalibrateAstigmatism` processor that is used to launch the interactive calibration, and2. a `ComputeTrajectories` class that describes the algorithm for fitting smoothed curves to the beads' x- and y-widths.After fitting the calibration curves, CalibrateAstigmatism provides a function known as `calibrationCurve` that takes a set of localizations as inputs and computes their axial positions as a result.<jupyter_code># Be sure not to use the %pylab inline option
%pylab
from bstore import processors as proc
from pathlib import Path
import pandas as pd<jupyter_output>Using matplotlib backend: Qt5Agg
Populating the interactive namespace from numpy and matplotlib<jupyter_text>Load the test dataThe test data for this example is in the [B-Store test files repository](https://github.com/kmdouglass/bstore_test_files). Download or clone this repository, and set the variable below to point to */processor_test_files/sequence-as-stack-Beads-AS-Exp_Localizations.csv*<jupyter_code>pathToData = Path('../../bstore_test_files/processor_test_files/sequence-as-stack-Beads-AS-Exp_Localizations.csv')

# Load the test data
with open(str(pathToData), 'r') as f:
    df = pd.read_csv(f)
df.describe()<jupyter_output><empty_output><jupyter_text>This dataset contains localizations from six fluorescent beads that are scanned axially and with a cylindrical lens in the imaging path. The dataset is synthetic and comes from the [2016 SMLM Challenge](http://bigwww.epfl.ch/smlm/challenge2016/index.html); for convenience, the localizations have already been computed from the z-stacks.The known z-position is in the column named **z [nm]**. The PSF widths in x and y are in the columns named **sigma_x [nm]** and **sigma_y [nm]** respectively.Let's plot the localizations' x- and y-positions:<jupyter_code>plt.scatter(df['x [nm]'], -df['y [nm]'], s=2)
plt.xlabel('x-position')
plt.ylabel('y-position')
plt.axis('equal')
plt.grid(True)
plt.show()<jupyter_output><empty_output><jupyter_text>Axial calibrationsIf you have already worked through the Fiducial Drift Correction notebook, then this part will look familiar. We start by defining a `CalibrateAstigmatism` processor with a few custom parameters. Then, we feed it with our localization file. A window appears that shows a 2D histogram of the density of localizations, allowing us to manually select the beads. We can select any number of regions we like by clicking and dragging a rectangle around the regions.After a region is drawn, **press the space bar to add it to the processor**. You may then select another region in the same manner. To finish searching for beads, simply close the window.Try selecting the bead at \\( x = 4.3 \, \mu m \\), \\( y = 7 \, \mu m \\) and closing the window afterward.<jupyter_code># coordCols = ['x', 'y'] by default
ca = proc.CalibrateAstigmatism(coordCols=['x [nm]', 'y [nm]'],
                               sigmaCols=['sigma_x [nm]', 'sigma_y [nm]'],
                               zCol='z [nm]')
ca.astigmatismComputer.smoothingWindowSize=20
ca.astigmatismComputer.smoothingFilterSize=3

processed_df = ca(df)
processed_df.describe()<jupyter_output><empty_output><jupyter_text>The processed DataFrame is actually the same as the input; no changes were made. However, if all went well, the CalibrateAstigmatism processor has fit splines to the PSF widths as a function of z and computed the calibration curve.Before discussing what the `CalibrateAstigmatism` processor did, let's go over how it was used. First, we create the processor by setting some of its optional parameters:```ca = proc.CalibrateAstigmatism(coordCols=['x [nm]', 'y [nm]'],                               sigmaCols=['sigma_x [nm]', 'sigma_y [nm]'],                               zCol='z [nm]')````coordCols` contains the names of the columns of the x- and y-positions of the localizations and `sigmaCols` contains the names of the columns containing the widths of the PSF in the x- and y-directions. Finally, `zCol` is the name of the known z-position of the beads. (This is known because it is controlled during the acquisition of the image stack of the fluorescent beads.) Some of the other optional parameters we could set are1. `startz` : where the z-fitting should begin2. `stopz` : where the z-fitting should end3. `astigmatismComputer` : this is a `ComputeTrajectories` object for calculating the spline fits to the individual bead trajectories in z. You can write your own algorithm for fitting z-trajectories and feed it to the processor using this parameter.4. `interactiveSearch` : a True/False value that determines whether a window is displayed to allow the user to select regions containing beads. You would set this to False if you already found beads but want to refit them using some different spline fitting parameters of the astigmatismComputer.5. `wobbleComputer` : this is another `ComputeTrajectories` object that calculates the PSF's centroid's position as a function of z, also known as wobble.Next, we adjust some of the smoothing spline parameters. These parameters are not part of the `CalibrateAstigmatism` processor; rather they belong to the `DefaultAstigmatismComputer` which belongs to the processor. The `DefaultAstigmatismComputer` is simply a type of `ComputeTrajectories` object for computing astigmatism calibration curves.```ca.astigmatismComputer.smoothingWindowSize=20ca.astigmatismComputer.smoothingFilterSize=3````smoothingWindowSize` is the size of the moving window that is used to weight the points in the trajectory during the spline fitting; `smoothingFilterSize` is the standard deviation of the Gaussian weighting function.Finally, we perform the calibration by calling the processor on the DataFrame:```processed_df = ca(df)```  Plot the calibration curves and bead localizationsJust like with the FiducialDriftCorrect processor, we can plot the individual localizations belonging to the bead we selected as a function of z, as well as the average smoothing spline.<jupyter_code>ca.astigmatismComputer.plotBeads()<jupyter_output><empty_output><jupyter_text>You should see a single window appear containing two plots. The top is a plot of the PSF width in x vs. z, and the bottom is a plot of the PSF width in y vs. z. The data points are the individual localizations and the curves are the two splines that fit to each trajectory.  Modifying the bead fits Changing which beads are used in the average trajectoryLet's now rerun the processor. This time, select at least two regions containing beads. (I selected the same one as before and another at \\(x = 9.3 \, \mu m \\) and \\(y = 2.7 \, \mu m\\).<jupyter_code>processed_df = ca(df)
ca.astigmatismComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>By selecting multiple beads, we tell B-Store to compute the average of the individually-fitted splines. This average spline is displayed in the plots as the solid, continuous curve plotted over the data ponts. If, for some reason, one of the beads was noisy or simply not good, then the average spline may not accurately represent the astigmatism present in the system. We can request that the `DefaultAstigmatismComputer` use only certain beads by setting its `useTrajectories` parameter.<jupyter_code># Recompute the average spline without selecting beads first
ca.interactiveSearch = False
ca.astigmatismComputer.useTrajectories = [1] # Use only bead number 1

_ = ca(df) # underscore means don't bother capturing the output
ca.astigmatismComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>Now the points belonging to bead number 0 will appear in gray; this indicates that they were not used in the fit. If you look closely, you will also see that the spline has changed very slightly and fits only the localizations belonging to bead number 1.If you decide that you really do want to use all the beads, we can indicate this by setting `useTrajectories` to the empty list (`[]`).<jupyter_code>ca.astigmatismComputer.useTrajectories = [] # Use all beads

_ = ca(df)
ca.astigmatismComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>Changing the fit rangeYou may also find that the full axial range in the data contains regions that are noisy or not well fit. We can select a smaller axial region to fit using the `startz` and `stopz` parameters of the `DefaultAstigmatismComputer`.<jupyter_code>ca.astigmatismComputer.startz = -300
ca.astigmatismComputer.stopz  = 300

_ = ca(df)
ca.astigmatismComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>You should now see gray x's corresponding to data points that are outside the fitting range. You should also see that the average spline now only covers the range \\( \left[ -300 \, \mu m, 300 \, \mu m \right] \\).You will also see a notice that startz and stopz parameters of the wobble computer were updated as well. Its startz and stopz parameters are always synchronized with the astigmatism computer to ensure that all fits are performed on the same range.  Changing the spline smoothing parametersSimilarly, we can change the smoothing parameters of the cubic spline after we have already selected beads.<jupyter_code>ca.astigmatismComputer.reset()
ca.astigmatismComputer.smoothingWindowSize = 50
ca.astigmatismComputer.smoothingFilterSize = 25

_ = ca(df)
ca.astigmatismComputer.plotBeads()
ca.astigmatismComputer.reset()

_ = ca(df)
ca.astigmatismComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>Adjust the wobble curvesWobble is the x- and y-position of the PSF centroid as a function of the axial position. The trajectories of the beads' centroid in x and y as a function z is a wobble curve and may be used to correct errors made by false assumptions about the aberrations present in the PSF. (See [Carlini et al., PLoS One 2015](http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0142949) for more information.)The wobble computer is tuned in much the same way as the astigmatism computer. If you don't want to correct for wobble, then simply ignore this feature.<jupyter_code>ca.wobbleComputer.plotBeads()
ca.wobbleComputer.smoothingWindowSize = 10
ca.wobbleComputer.smoothingFilterSize = 2

# The following are locked to the value of startz and stopz from the
# astigmatism computer and therefore do not do anything.
ca.wobbleComputer.startz = -300 # Does nothing!
ca.wobbleComputer.stopz  = 300  # Does nothing!

_ = ca(df)
ca.wobbleComputer.plotBeads()<jupyter_output>Setting wobble fiting range to the match the astigmatism fit range. startz and stopz are set in the astigmatism computer.
Performing spline fits...
Performing spline fits...<jupyter_text>Using the calibrations to axially localize a datasetOnce calibrated, the `CalibrateAstigmatism` processor contains a property called `calibrationCurves` that holds the spline fits to \\( W_x \\) vs. \\(z\\) and \\( W_y \\) vs. z, where \\( W_x \\) and \\( W_y \\) are the PSF widths in x and y, respectively. These fits are functions, which means they that they accept a single number (the z-coordinate) as an input and produce the width in x and y as outputs.<jupyter_code>ca.calibrationCurves<jupyter_output><empty_output><jupyter_text>We can use these functions in B-Store's `ComputeZPosition` processor. To initialize the processor, we need to specify which functions to use and, optionally, the names of the columns containing the z-positions and PSF widths.<jupyter_code>cz = proc.ComputeZPosition(ca.calibrationCurves,
                           zCol='z [nm]',
                           sigmaCols=['sigma_x [nm]', 'sigma_y [nm]'])<jupyter_output><empty_output><jupyter_text>Though we didn't specify it above, the ComputeZPosition also accepts a parameter called `fittype` that takes one of two values: `diff` (the default value) and `huang`. `diff` computes the calibration curve by first sampling the two spline fits, taking their difference, and then reinterpolating to produce a monotonic calibration curve that transforms \\(W_x - W_y\\) into \\( z \\). In general, it is very fast. The `huang` method computes the z-position by minimizing an objective function related to the distance between the experimental and calibrated PSF widths. This method was used in the first astigmatic STORM paper [Huang, et al. Science 319, 810-813 (2008)](http://science.sciencemag.org/content/319/5864/810). Because each localization requires a call to an optimiztion method, it is much smaller.Having created the processor, let's now load a test dataset and localize it in z.<jupyter_code>pathToData = Path('../../bstore_test_files/processor_test_files/MT0.N1.LD-AS-Exp_Localizations.csv')

with open(str(pathToData), 'r') as f:
    locs = pd.read_csv(f)
    
locs.head()<jupyter_output><empty_output><jupyter_text>You will notice that we already specified the correct column names when creating the `ComputeZPosition` processor.Now, we simply pass these localizations to the processor and the z-position is computed automatically.<jupyter_code>locs_z = cz(locs)
locs_z.head()<jupyter_output><empty_output><jupyter_text>Correcting wobbleTo also correct for wobble, we can specify a few extra parameters to the `ComputeZPosition` processor, including the wobble curves calculated by the `CalibrateAstigmatism` processor.<jupyter_code>ca.wobbleCurves
cz = proc.ComputeZPosition(ca.calibrationCurves,
                           zCol='z [nm]',
                           coordCols=['x [nm]', 'y [nm]'],
                           sigmaCols=['sigma_x [nm]', 'sigma_y [nm]'],
                           wobbleFunc=ca.wobbleCurves)
locs_z_wobble = cz(locs)
locs_z_wobble.head()<jupyter_output>/home/douglass/anaconda3/envs/bstore/lib/python3.5/site-packages/scipy/interpolate/interpolate.py:612: RuntimeWarning: invalid value encountered in less
  below_bounds = x_new < self.x[0]
/home/douglass/anaconda3/envs/bstore/lib/python3.5/site-packages/scipy/interpolate/interpolate.py:613: RuntimeWarning: invalid value encountered in greater
  above_bounds = x_new > self.x[-1]<jupyter_text>Now you see that a small offset has been applied to the **x [nm]** and **y [nm]** columns to correct these localizations for wobble. The value of the offset has been saved in the **dx** and **dy** columns.  Modifying the trajectory-fitting algorithm*You may skip this section if you do not want to program your own astigmatism computer.*By default, B-Store uses a curve fitting algorithm based on a cubic smoothing spline with weights determined by a Gaussian filter. The algorithm is implemented in a class called `DefaultAstigmatismComputer` which uses the `ComputeTrajectories` interface. You can write your own astigmatism computer by inheriting this interface.<jupyter_code>import inspect
print(inspect.getsource(proc.ComputeTrajectories))<jupyter_output>class ComputeTrajectories(metaclass=ABCMeta):
    """Basic functionality for computing trajectories from localizations.
    
    This is used to compute trajectories from regions of a dataset containing
    localizations, such as fiducial drift trajectories (position vs. frame
    number)or astigmatic calibration curves (PSF width vs. z).

    Attributes
    ----------
    regionLocs : Pandas DataFrame
        The localizations for individual regions.

    """

    def __init__(self):
        """Initializes the trajectory computer.

        """
        self._regionData = None
        
    def _plotCurves(self, curveNumber=None, coordCols=['x', 'y'],
                    horizontalLabels=['', 'time'], verticalLabels=['x', 'y'],
                    title='trajectories', splineCols=['t','x','y'],
                    offsets=[0,0], ylims=[-100, 500, -100, 500]):
        """Make a plot of each region's trajectory and the average spline fit.

        plotCurves allows the user to check the tr[...]<jupyter_text>The `ComputeTrajectories` interface provides a property and four methods:1. `regionLocs` contains a DataFrame with all of the localizations. It must have at least one index with the label 'region_id' that identifies which region the localizations came from.2. `clearRegionLocs()` removes the localization information that is held by the computer.3. `_plotCurves()` is the code used to plot the trajectories.4. `_movingAverage()` is the sliding window Gaussian filter used to weight the datapoints for the cubic smoothing spline.5. `reset()` resets the computer to its initial state.In addition, there is one abstract method called `computeTrajectory`. Any class that implements this interface must define a function with this name.As an example, the actual implementation of this interface by the `DefaultAstigmatismComputer` is printed below:<jupyter_code>print(inspect.getsource(proc.DefaultAstigmatismComputer.computeTrajectory))<jupyter_output>def computeTrajectory(self, locs):
        """Computes the final drift trajectory from fiducial localizations.

        Parameters
        ----------
        locs        : Pandas DataFrame
            DataFrame containing the localizations belonging to beads.

        Returns
        -------
        avgSpline : Pandas DataFrame
            A dataframe containing z-positions and PSF widths in x- and y- for
            calibrating an astigmatic imaging measurement.

        """
        z = self.zCol
        if self.startz:
            startz = self.startz
        else:
            startz = locs[z].min()
        
        if self.stopz:
            stopz = self.stopz
        else:
            stopz  = locs[z].max()
        
        self.clearRegionLocs()
        self.regionLocs = locs
        self._removeOutliers(startz, stopz)
        self.fitCurves()
        self.combineCurves(startz, stopz)

        return self.avgSpline<jupyter_text>The method returns the averaged splines for the PSF widths in each direction. This is a Pandas DataFrame with  columns named `z`, `xS` and `yS`.<jupyter_code># Print the first five values of the DataFrame returned by the drift computer
ca.astigmatismComputer.avgSpline.head()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup