<jupyter_start><jupyter_text>Tracking objects in 3D space---Here we are going to look at calculating the linkage probability between a new object observation and the existing measurements.Using BayesianTracker from https://github.com/quantumjot---Kalman filter:Excellent introduction to Kalman filters:http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdfPrediction:\begin{equation}\begin{aligned}\hat{x}_t &= A\hat{x}_{t-1} \\P_t &= AP_{t-1}A^T+Q\end{aligned}\end{equation}The process covariance can be estimated like this:\begin{equation}Q = G \cdot G^T \cdot \sigma_v^2\end{equation}Update:\begin{equation}\begin{aligned}K_t &= P_tH^T(HP_tH^T+R)^{-1} \\\hat{x}_t &= \hat{x}_t + K_t(z_k-H\hat{x}_t) \\P_k &= (I-K_tH)P_t\end{aligned}\end{equation}<jupyter_code>%matplotlib inline

# temporarily added BayesianTracker to the path
import sys
sys.path.append('../btrack/')



# math
import numpy as np
import scipy.stats as stats
import scipy.integrate as integrate

# plotting
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D


# btrack components
import core
#import btypes
import utils
from agents import AgentManager
from render import plot_tracks

# instantiate agents
num_chaff = 500
num_frames = 1500

manager = AgentManager(num_agents=1, noise=0.1, periodic_boundary=False)
manager.simulate(num_frames)

def chaff(n_frames):
    """ Random observation """
    position = np.random.randint(0,high=1000, size=(3,))
    obj = btypes.PyTrackObject()
    obj.ID = int(1e5)
    obj.x = position[0]
    obj.y = position[1]
    obj.z = position[2]
    obj.t = np.random.randint(0,high=n_frames)
    obj.dummy = False
    obj.label = 0
    obj.probability = np.zeros((5,),dtype='float')
    return obj<jupyter_output><empty_output><jupyter_text>Start a BayesianTracker session<jupyter_code>data = manager.get_data()
# shuffle the data (BT sorts it anyway) and dropout a fraction to simulate segmentation misses
dropout_rate = 0.01
import random
random.shuffle(data)
dropout = int(float(len(data))*dropout_rate)
for i in range(dropout):
    data.pop(0)
# count the number of metaphase observations
dl = [d.label for d in data]
dl.count(2)
print num_frames
# data = utils.crop_volume(data, ((0,1600), (0,1600), (0,1600), (0, num_frames)))
# start a session
with core.BayesianTracker() as trk:

    trk.append(data)

    trk.configure_from_file('MDCK_config.json')
    trk.return_kalman = True

    
    # trk.volume = ((0,1024),(0,1024),(0,1024))
    trk.track_interactive(step_size=500)
    
    # get the tracks
    tracks = trk.tracks

    
    # get the imaging volume (xyzt)
    volume = trk.volume

    # generate some hypotheses and optimise
    hypotheses = trk.hypotheses()
    trk.optimise()
    
    # get the optimized tracks
    tracks2 = trk.tracks


tracks == tracks2
print volume

# count how many divisions we see
mitotic_tracks = [trk for trk in tracks2 if trk.fate==4]
for t in mitotic_tracks: print t.ID,'->', t.children<jupyter_output>7 -> [13, 12]
12 -> [21, 20]
1 -> [2, 3]
2 -> [7, 6]
6 -> [11, 10]
3 -> [4, 5]
4 -> [8, 9]
8 -> [14, 15]
10 -> [16, 17]
13 -> [18, 19]<jupyter_text>Plot the trajectories and the original data---Plot the output of the tracking...<jupyter_code>plot_tracks(tracks, order='xyt', box=volume, cmap=plt.get_cmap('hsv'), labels=True)
plot_tracks(tracks2, order='xt', box=volume, cmap=plt.get_cmap('hsv'), labels=True, title='BayesianTracker output after optimization')
trk = tracks[0]


kal_mu_x = trk.kalman(1)
kal_mu_y = trk.kalman(2)
kal_pred_x = trk.kalman(13)
kal_pred_y = trk.kalman(14)

#a = agents[0]

plt.figure(figsize=(16,10))
plt.plot(trk.x,trk.y,'k.-')
# plt.plot(trk2.x,trk2.y,'k-')
plt.plot(kal_mu_x, kal_mu_y,'b:')
plt.plot(kal_pred_x, kal_pred_y,'r:')
plt.axis('image')
plt.legend(['raw', 'kalman prediction', 'tracker prediction'])
plt.show()<jupyter_output><empty_output><jupyter_text>inspect a tracklet<jupyter_code>def inspect(trk_ID):
    trk = tracks2[trk_ID]
    dims = ['x','y','z']
    plt.figure(figsize=(16,10))
    for n, d in enumerate(dims):
        plt.subplot(3,1,n+1)
        mu = trk.kalman(n+1)
        cov = trk.kalman(4+n*3+n) # argh! complex indexing
        plt.errorbar(trk.kalman(0), mu, yerr=3.*cov, color='lightgrey')
        plt.plot(trk.kalman(0), mu, color='b')
        plt.xlabel('Time (frames)')
        plt.ylabel(r'{0:s} (Pixels $\pm \sigma_{0:s}^2$)'.format(d))
        plt.autoscale(enable=True, axis='x', tight=True)
        plt.autoscale(enable=True, axis='y', tight=True)
    plt.show()
inspect(0)
from optimise import lineage
lineage_tree = lineage.LineageTree(tracks2)
trees = lineage_tree.create()
plt.figure(figsize=(16,10))
lineage_tree.plot()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1509995746490855, 'fraction_numerical': 0.033390046788600594, 'mean_word_length': 4.462253193960511, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 27, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '236227', 'license': 'CNRI-Python', 'n_tokens_mistral': 1899, 'n_tokens_neox': 1767, 'n_words': 388, 'path': 'Tracker_Development/Jupyter_Notebooks/TrackingObjectsIn3D.ipynb', 'repo_name': 'The-Kristina/CellComp', 'size': '4702'}