Skip to content

deeplabcut.refine_training_dataset.stitch

Classes:

Name Description
Tracklet

Functions:

Name Description
stitch_tracklets

Stitch sparse tracklets into full tracks via a graph-based, minimum-cost flow

Tracklet

Methods:

Name Description
__add__

Join this tracklet to another one.
__contains__

Test whether tracklets temporally overlap.
__gt__

Test whether this tracklet follows the other one.
__init__

Create a Tracklet object.
__lt__

Test whether this tracklet precedes the other one.
box_overlap_with

Calculate the overlap between each Tracklet's bounding box.
calc_rate_of_turn

Calculate the rate of turn (or angular velocity) of either the head or
calc_velocity

Calculate the linear velocity of either the head or tail of the Tracklet,
contains_duplicates

Evaluate whether the Tracklet contains duplicate time indices.
distance_to

Calculate the Euclidean distance between this Tracklet and another.
dynamic_dissimilarity_with

Compute a dissimilarity score between Hankelets. This metric efficiently
dynamic_similarity_with

Evaluate the complexity of the tracklets' underlying dynamics from the rank
estimate_rank

Estimate the (low) rank of a noisy matrix via hard thresholding of singular
immediately_follows

Test whether this Tracklet follows another within a tolerance ofmax_gap
motion_affinity_with

Evaluate the motion affinity of this Tracklet' with another one.
shape_dissimilarity_with

Calculate the dissimilarity in shape between this Tracklet and another.
time_gap_to

Return the time gap separating this Tracklet to another.

Attributes:

Name Type Description
centroid

Return the instantaneous 2D position of the Tracklet centroid.
end

Return the time at which the tracklet ends.
identity

Return the average predicted identity of all Tracklet detections.
is_continuous

Test whether there are gaps in the time indices.
likelihood

Return the average likelihood of all Tracklet detections.
start

Return the time at which the tracklet starts.
xy

Return the x and y coordinates.
Source code in deeplabcut/refine_training_dataset/stitch.py 
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
class Tracklet:
    def __init__(self, data, inds):
        """Create a Tracklet object.

        Parameters
        ----------
        data : ndarray
            3D array of shape (nframes, nbodyparts, 3 or 4), where the last
            dimension is for x, y, likelihood and, optionally, identity.
        inds : array-like
            Corresponding time frame indices.
        """
        if data.ndim != 3 or data.shape[-1] not in (3, 4):
            raise ValueError("Data must of shape (nframes, nbodyparts, 3 or 4)")

        if data.shape[0] != len(inds):
            raise ValueError("Data and corresponding indices must have the same length.")

        self.data = data.astype(np.float64)
        self.inds = np.array(inds)
        monotonically_increasing = all(a < b for a, b in zip(inds, inds[1:], strict=False))
        if not monotonically_increasing:
            idx = np.argsort(inds, kind="mergesort")  # For stable sort with duplicates
            self.inds = self.inds[idx]
            self.data = self.data[idx]
        self._centroid = None

    def __len__(self):
        return self.inds.size

    def __add__(self, other):
        """Join this tracklet to another one."""
        data = np.concatenate((self.data, other.data))
        inds = np.concatenate((self.inds, other.inds))
        return Tracklet(data, inds)

    def __radd__(self, other):
        if other == 0:
            return self
        return self.__add__(other)

    def __sub__(self, other):
        mask = np.isin(self.inds, other.inds, assume_unique=True)
        if mask.all():
            return None
        return Tracklet(self.data[~mask], self.inds[~mask])

    def __lt__(self, other):
        """Test whether this tracklet precedes the other one."""
        return self.end < other.start

    def __gt__(self, other):
        """Test whether this tracklet follows the other one."""
        return self.start > other.end

    def __contains__(self, other_tracklet):
        """Test whether tracklets temporally overlap."""
        return np.isin(self.inds, other_tracklet.inds, assume_unique=True).any()

    def __repr__(self):
        return f"Tracklet of length {len(self)} from {self.start} to {self.end} with reliability {self.likelihood:.3f}"

    @property
    def xy(self):
        """Return the x and y coordinates."""
        return self.data[..., :2]

    @property
    def centroid(self):
        """Return the instantaneous 2D position of the Tracklet centroid.

        For Tracklets longer than 10 frames, the centroid is automatically smoothed
        using an exponential moving average. The result is cached for efficiency.
        """
        if self._centroid is None:
            self._update_centroid()
        return self._centroid

    def _update_centroid(self):
        like = self.data[..., 2:3] + 1e-10  # Avoid division by zero in very uncertain tracklets
        self._centroid = np.nansum(self.xy * like, axis=1) / np.nansum(like, axis=1)

    @property
    def likelihood(self):
        """Return the average likelihood of all Tracklet detections."""
        return np.nanmean(self.data[..., 2])

    @property
    def identity(self):
        """Return the average predicted identity of all Tracklet detections."""
        try:
            return mode(
                self.data[..., 3],
                axis=None,
                nan_policy="omit",
                keepdims=False,
            )[0]
        except IndexError:
            return -1

    @property
    def start(self):
        """Return the time at which the tracklet starts."""
        return self.inds[0]

    @property
    def end(self):
        """Return the time at which the tracklet ends."""
        return self.inds[-1]

    @property
    def flat_data(self):
        return self.data[..., :3].reshape((len(self)), -1)

    def get_data_at(self, ind):
        return self.data[np.searchsorted(self.inds, ind)]

    def set_data_at(self, ind, data):
        self.data[np.searchsorted(self.inds, ind)] = data

    def del_data_at(self, ind):
        idx = np.searchsorted(self.inds, ind)
        self.inds = np.delete(self.inds, idx)
        self.data = np.delete(self.data, idx, axis=0)
        self._update_centroid()

    def interpolate(self, max_gap=1):
        if max_gap < 1:
            raise ValueError("Gap should be a strictly positive integer.")

        gaps = np.diff(self.inds) - 1
        valid_gaps = (0 < gaps) & (gaps <= max_gap)
        fills = []
        for i in np.flatnonzero(valid_gaps):
            s, e = self.inds[[i, i + 1]]
            data1, data2 = self.data[[i, i + 1]]
            diff = (data2 - data1) / (e - s)
            diff[np.isnan(diff)] = 0
            interp = diff[..., np.newaxis] * np.arange(1, e - s)
            data = data1 + np.rollaxis(interp, axis=2)
            data[..., 2] = 0.5  # Chance detections
            if data.shape[1] == 4:
                data[:, 3] = self.identity
            fills.append(Tracklet(data, np.arange(s + 1, e)))
        if not fills:
            return self
        return self + sum(fills)

    def contains_duplicates(self, return_indices=False):
        """Evaluate whether the Tracklet contains duplicate time indices.

        If `return_indices`, also return the indices of the duplicates.
        """
        has_duplicates = len(set(self.inds)) != len(self.inds)
        if not return_indices:
            return has_duplicates
        return has_duplicates, np.flatnonzero(np.diff(self.inds) == 0)

    def calc_velocity(self, where="head", norm=True):
        """Calculate the linear velocity of either the `head` or `tail` of the Tracklet,
        computed over the last or first three frames, respectively.

        If `norm`, return the absolute
        speed rather than a 2D vector.
        """
        if where == "tail":
            vel = np.diff(self.centroid[:3], axis=0) / np.diff(self.inds[:3])[:, np.newaxis]
        elif where == "head":
            vel = np.diff(self.centroid[-3:], axis=0) / np.diff(self.inds[-3:])[:, np.newaxis]
        else:
            raise ValueError(f"Unknown where={where}")
        if norm:
            return np.sqrt(np.sum(vel**2, axis=1)).mean()
        return vel.mean(axis=0)

    @property
    def maximal_velocity(self):
        vel = np.diff(self.centroid, axis=0) / np.diff(self.inds)[:, np.newaxis]
        return np.sqrt(np.max(np.sum(vel**2, axis=1)))

    def calc_rate_of_turn(self, where="head"):
        """Calculate the rate of turn (or angular velocity) of either the `head` or
        `tail` of the Tracklet, computed over the last or first three frames,
        respectively."""
        if where == "tail":
            v = np.diff(self.centroid[:3], axis=0)
        else:
            v = np.diff(self.centroid[-3:], axis=0)
        theta = np.arctan2(v[:, 1], v[:, 0])
        return (theta[1] - theta[0]) / (self.inds[1] - self.inds[0])

    @property
    def is_continuous(self):
        """Test whether there are gaps in the time indices."""
        return self.end - self.start + 1 == len(self)

    def immediately_follows(self, other_tracklet, max_gap=1):
        """Test whether this Tracklet follows another within a tolerance of`max_gap`
        frames."""
        return 0 < self.start - other_tracklet.end <= max_gap

    def distance_to(self, other_tracklet):
        """Calculate the Euclidean distance between this Tracklet and another.

        If the Tracklets overlap in time, this is the mean distance over those frames.
        Otherwise, it is the distance between the head/tail of one to the tail/head of
        the other.
        """
        if self in other_tracklet:
            dist = (
                self.centroid[np.isin(self.inds, other_tracklet.inds)]
                - other_tracklet.centroid[np.isin(other_tracklet.inds, self.inds)]
            )
            return np.sqrt(np.sum(dist**2, axis=1)).mean()
        elif self < other_tracklet:
            return np.sqrt(np.sum((self.centroid[-1] - other_tracklet.centroid[0]) ** 2))
        else:
            return np.sqrt(np.sum((self.centroid[0] - other_tracklet.centroid[-1]) ** 2))

    def motion_affinity_with(self, other_tracklet):
        """Evaluate the motion affinity of this Tracklet' with another one.

        This evaluates whether the Tracklets could realistically be reached by one
        another, knowing the time separating them and their velocities. Return 0 if the
        Tracklets overlap.
        """
        time_gap = self.time_gap_to(other_tracklet)
        if time_gap > 0:
            if self < other_tracklet:
                d1 = self.centroid[-1] + time_gap * self.calc_velocity(norm=False)
                d2 = other_tracklet.centroid[0] - time_gap * other_tracklet.calc_velocity("tail", False)
                delta1 = other_tracklet.centroid[0] - d1
                delta2 = self.centroid[-1] - d2
            else:
                d1 = other_tracklet.centroid[-1] + time_gap * other_tracklet.calc_velocity(norm=False)
                d2 = self.centroid[0] - time_gap * self.calc_velocity("tail", False)
                delta1 = self.centroid[0] - d1
                delta2 = other_tracklet.centroid[-1] - d2
            return (np.sqrt(np.sum(delta1**2)) + np.sqrt(np.sum(delta2**2))) / 2
        return 0

    def time_gap_to(self, other_tracklet):
        """Return the time gap separating this Tracklet to another."""
        if self in other_tracklet:
            t = 0
        elif self < other_tracklet:
            t = other_tracklet.start - self.end
        else:
            t = self.start - other_tracklet.end
        return t

    def shape_dissimilarity_with(self, other_tracklet):
        """Calculate the dissimilarity in shape between this Tracklet and another."""
        if self in other_tracklet:
            dist = np.inf
        elif self < other_tracklet:
            dist = self.undirected_hausdorff(self.xy[-1], other_tracklet.xy[0])
        else:
            dist = self.undirected_hausdorff(self.xy[0], other_tracklet.xy[-1])
        return dist

    def box_overlap_with(self, other_tracklet):
        """Calculate the overlap between each Tracklet's bounding box."""
        if self in other_tracklet:
            overlap = 0
        else:
            if self < other_tracklet:
                bbox1 = self.calc_bbox(-1)
                bbox2 = other_tracklet.calc_bbox(0)
            else:
                bbox1 = self.calc_bbox(0)
                bbox2 = other_tracklet.calc_bbox(-1)
            overlap = calc_iou(bbox1, bbox2)
        return overlap

    @staticmethod
    def undirected_hausdorff(u, v):
        return max(directed_hausdorff(u, v)[0], directed_hausdorff(v, u)[0])

    def calc_bbox(self, ind):
        xy = self.xy[ind]
        bbox = np.empty(4)
        bbox[:2] = np.nanmin(xy, axis=0)
        bbox[2:] = np.nanmax(xy, axis=0)
        return bbox

    @staticmethod
    def hankelize(xy):
        ncols = int(np.ceil(len(xy) * 2 / 3))
        nrows = len(xy) - ncols + 1
        mat = np.empty((2 * nrows, ncols))
        mat[::2] = hankel(xy[:nrows, 0], xy[-ncols:, 0])
        mat[1::2] = hankel(xy[:nrows, 1], xy[-ncols:, 1])
        return mat

    def to_hankelet(self):
        # See Li et al., 2012. Cross-view Activity Recognition using Hankelets.
        # As proposed in the paper, the Hankel matrix can either be formed from
        # the tracklet's centroid or its normalized velocity.
        # vel = np.diff(self.centroid, axis=0)
        # vel /= np.linalg.norm(vel, axis=1, keepdims=True)
        # return self.hankelize(vel)
        return self.hankelize(self.centroid)

    def dynamic_dissimilarity_with(self, other_tracklet):
        """Compute a dissimilarity score between Hankelets. This metric efficiently
        captures the degree of alignment of the subspaces spanned by the columns of both
        matrices.

        See Li et al., 2012.     Cross-view Activity Recognition using Hankelets.
        """
        hk1 = self.to_hankelet()
        hk1 /= np.linalg.norm(hk1)
        hk2 = other_tracklet.to_hankelet()
        hk2 /= np.linalg.norm(hk2)
        min_shape = min(hk1.shape + hk2.shape)
        temp1 = (hk1 @ hk1.T)[:min_shape, :min_shape]
        temp2 = (hk2 @ hk2.T)[:min_shape, :min_shape]
        return 2 - np.linalg.norm(temp1 + temp2)

    def dynamic_similarity_with(self, other_tracklet, tol=0.01):
        """Evaluate the complexity of the tracklets' underlying dynamics from the rank
        of their Hankel matrices, and assess whether they originate from the same track.
        The idea is that if two tracklets are part of the same track, they can be
        approximated by a low order regressor. Conversely, tracklets belonging to
        different tracks will require a higher order regressor.

        See Dicle et al., 2013.
            The Way They Move: Tracking Multiple Targets with Similar Appearance.
        """
        # TODO Add missing data imputation
        joint_tracklet = self + other_tracklet
        joint_rank = joint_tracklet.estimate_rank(tol)
        rank1 = self.estimate_rank(tol)
        rank2 = other_tracklet.estimate_rank(tol)
        return (rank1 + rank2) / joint_rank - 1

    def estimate_rank(self, tol):
        """Estimate the (low) rank of a noisy matrix via hard thresholding of singular
        values.

        See Gavish & Donoho, 2013.     The optimal hard threshold for singular values is
        4/sqrt(3)
        """
        mat = self.to_hankelet()
        if np.any(mat):  # check that the matrix contains non-zero entries
            # nrows, ncols = mat.shape
            # beta = nrows / ncols
            # omega = 0.56 * beta ** 3 - 0.95 * beta ** 2 + 1.82 * beta + 1.43
            _, s, _ = sli.svd(mat, min(10, min(mat.shape)))
        else:
            s = np.zeros(min(10, min(mat.shape)))

        # return np.argmin(s > omega * np.median(s))
        eigen = s**2
        diff = np.abs(np.diff(eigen / eigen[0]))
        return np.argmin(diff > tol)

    def plot(self, centroid_only=True, color=None, ax=None, interactive=False):
        if ax is None:
            fig, ax = plt.subplots()
        centroid = np.full((self.end + 1, 2), np.nan)
        centroid[self.inds] = self.centroid
        lines = ax.plot(centroid, c=color, lw=2, picker=interactive)
        if not centroid_only:
            xy = np.full((self.end + 1, self.xy.shape[1], 2), np.nan)
            xy[self.inds] = self.xy
            ax.plot(xy[..., 0], c=color, lw=1)
            ax.plot(xy[..., 1], c=color, lw=1)
        return lines

centroid property 

centroid

Return the instantaneous 2D position of the Tracklet centroid.

For Tracklets longer than 10 frames, the centroid is automatically smoothed using an exponential moving average. The result is cached for efficiency.

end property 

end

Return the time at which the tracklet ends.

identity property 

identity

Return the average predicted identity of all Tracklet detections.

is_continuous property 

is_continuous

Test whether there are gaps in the time indices.

likelihood property 

likelihood

Return the average likelihood of all Tracklet detections.

start property 

start

Return the time at which the tracklet starts.

xy property 

xy

Return the x and y coordinates.

__add__ 

__add__(other)

Join this tracklet to another one.

Source code in deeplabcut/refine_training_dataset/stitch.py 
73
74
75
76
77
def __add__(self, other):
    """Join this tracklet to another one."""
    data = np.concatenate((self.data, other.data))
    inds = np.concatenate((self.inds, other.inds))
    return Tracklet(data, inds)

__contains__ 

__contains__(other_tracklet)

Test whether tracklets temporally overlap.

Source code in deeplabcut/refine_training_dataset/stitch.py 
 98
 99
100
def __contains__(self, other_tracklet):
    """Test whether tracklets temporally overlap."""
    return np.isin(self.inds, other_tracklet.inds, assume_unique=True).any()

__gt__ 

__gt__(other)

Test whether this tracklet follows the other one.

Source code in deeplabcut/refine_training_dataset/stitch.py 
94
95
96
def __gt__(self, other):
    """Test whether this tracklet follows the other one."""
    return self.start > other.end

__init__ 

__init__(data, inds)

Create a Tracklet object.

Parameters

data : ndarray 3D array of shape (nframes, nbodyparts, 3 or 4), where the last dimension is for x, y, likelihood and, optionally, identity. inds : array-like Corresponding time frame indices.

Source code in deeplabcut/refine_training_dataset/stitch.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def __init__(self, data, inds):
    """Create a Tracklet object.

    Parameters
    ----------
    data : ndarray
        3D array of shape (nframes, nbodyparts, 3 or 4), where the last
        dimension is for x, y, likelihood and, optionally, identity.
    inds : array-like
        Corresponding time frame indices.
    """
    if data.ndim != 3 or data.shape[-1] not in (3, 4):
        raise ValueError("Data must of shape (nframes, nbodyparts, 3 or 4)")

    if data.shape[0] != len(inds):
        raise ValueError("Data and corresponding indices must have the same length.")

    self.data = data.astype(np.float64)
    self.inds = np.array(inds)
    monotonically_increasing = all(a < b for a, b in zip(inds, inds[1:], strict=False))
    if not monotonically_increasing:
        idx = np.argsort(inds, kind="mergesort")  # For stable sort with duplicates
        self.inds = self.inds[idx]
        self.data = self.data[idx]
    self._centroid = None

__lt__ 

__lt__(other)

Test whether this tracklet precedes the other one.

Source code in deeplabcut/refine_training_dataset/stitch.py 
90
91
92
def __lt__(self, other):
    """Test whether this tracklet precedes the other one."""
    return self.end < other.start

box_overlap_with 

box_overlap_with(other_tracklet)

Calculate the overlap between each Tracklet's bounding box.

Source code in deeplabcut/refine_training_dataset/stitch.py 
304
305
306
307
308
309
310
311
312
313
314
315
316
def box_overlap_with(self, other_tracklet):
    """Calculate the overlap between each Tracklet's bounding box."""
    if self in other_tracklet:
        overlap = 0
    else:
        if self < other_tracklet:
            bbox1 = self.calc_bbox(-1)
            bbox2 = other_tracklet.calc_bbox(0)
        else:
            bbox1 = self.calc_bbox(0)
            bbox2 = other_tracklet.calc_bbox(-1)
        overlap = calc_iou(bbox1, bbox2)
    return overlap

calc_rate_of_turn 

calc_rate_of_turn(where='head')

Calculate the rate of turn (or angular velocity) of either the head or tail of the Tracklet, computed over the last or first three frames, respectively.

Source code in deeplabcut/refine_training_dataset/stitch.py 
223
224
225
226
227
228
229
230
231
232
def calc_rate_of_turn(self, where="head"):
    """Calculate the rate of turn (or angular velocity) of either the `head` or
    `tail` of the Tracklet, computed over the last or first three frames,
    respectively."""
    if where == "tail":
        v = np.diff(self.centroid[:3], axis=0)
    else:
        v = np.diff(self.centroid[-3:], axis=0)
    theta = np.arctan2(v[:, 1], v[:, 0])
    return (theta[1] - theta[0]) / (self.inds[1] - self.inds[0])

calc_velocity 

calc_velocity(where='head', norm=True)

Calculate the linear velocity of either the head or tail of the Tracklet, computed over the last or first three frames, respectively.

If norm, return the absolute speed rather than a 2D vector.

Source code in deeplabcut/refine_training_dataset/stitch.py 
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
def calc_velocity(self, where="head", norm=True):
    """Calculate the linear velocity of either the `head` or `tail` of the Tracklet,
    computed over the last or first three frames, respectively.

    If `norm`, return the absolute
    speed rather than a 2D vector.
    """
    if where == "tail":
        vel = np.diff(self.centroid[:3], axis=0) / np.diff(self.inds[:3])[:, np.newaxis]
    elif where == "head":
        vel = np.diff(self.centroid[-3:], axis=0) / np.diff(self.inds[-3:])[:, np.newaxis]
    else:
        raise ValueError(f"Unknown where={where}")
    if norm:
        return np.sqrt(np.sum(vel**2, axis=1)).mean()
    return vel.mean(axis=0)

contains_duplicates 

contains_duplicates(return_indices=False)

Evaluate whether the Tracklet contains duplicate time indices.

If return_indices, also return the indices of the duplicates.

Source code in deeplabcut/refine_training_dataset/stitch.py 
191
192
193
194
195
196
197
198
199
def contains_duplicates(self, return_indices=False):
    """Evaluate whether the Tracklet contains duplicate time indices.

    If `return_indices`, also return the indices of the duplicates.
    """
    has_duplicates = len(set(self.inds)) != len(self.inds)
    if not return_indices:
        return has_duplicates
    return has_duplicates, np.flatnonzero(np.diff(self.inds) == 0)

distance_to 

distance_to(other_tracklet)

Calculate the Euclidean distance between this Tracklet and another.

If the Tracklets overlap in time, this is the mean distance over those frames. Otherwise, it is the distance between the head/tail of one to the tail/head of the other.

Source code in deeplabcut/refine_training_dataset/stitch.py 
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def distance_to(self, other_tracklet):
    """Calculate the Euclidean distance between this Tracklet and another.

    If the Tracklets overlap in time, this is the mean distance over those frames.
    Otherwise, it is the distance between the head/tail of one to the tail/head of
    the other.
    """
    if self in other_tracklet:
        dist = (
            self.centroid[np.isin(self.inds, other_tracklet.inds)]
            - other_tracklet.centroid[np.isin(other_tracklet.inds, self.inds)]
        )
        return np.sqrt(np.sum(dist**2, axis=1)).mean()
    elif self < other_tracklet:
        return np.sqrt(np.sum((self.centroid[-1] - other_tracklet.centroid[0]) ** 2))
    else:
        return np.sqrt(np.sum((self.centroid[0] - other_tracklet.centroid[-1]) ** 2))

dynamic_dissimilarity_with 

dynamic_dissimilarity_with(other_tracklet)

Compute a dissimilarity score between Hankelets. This metric efficiently captures the degree of alignment of the subspaces spanned by the columns of both matrices.

See Li et al., 2012. Cross-view Activity Recognition using Hankelets.

Source code in deeplabcut/refine_training_dataset/stitch.py 
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def dynamic_dissimilarity_with(self, other_tracklet):
    """Compute a dissimilarity score between Hankelets. This metric efficiently
    captures the degree of alignment of the subspaces spanned by the columns of both
    matrices.

    See Li et al., 2012.     Cross-view Activity Recognition using Hankelets.
    """
    hk1 = self.to_hankelet()
    hk1 /= np.linalg.norm(hk1)
    hk2 = other_tracklet.to_hankelet()
    hk2 /= np.linalg.norm(hk2)
    min_shape = min(hk1.shape + hk2.shape)
    temp1 = (hk1 @ hk1.T)[:min_shape, :min_shape]
    temp2 = (hk2 @ hk2.T)[:min_shape, :min_shape]
    return 2 - np.linalg.norm(temp1 + temp2)

dynamic_similarity_with 

dynamic_similarity_with(other_tracklet, tol=0.01)

Evaluate the complexity of the tracklets' underlying dynamics from the rank of their Hankel matrices, and assess whether they originate from the same track. The idea is that if two tracklets are part of the same track, they can be approximated by a low order regressor. Conversely, tracklets belonging to different tracks will require a higher order regressor.

See Dicle et al., 2013. The Way They Move: Tracking Multiple Targets with Similar Appearance.

Source code in deeplabcut/refine_training_dataset/stitch.py 
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
def dynamic_similarity_with(self, other_tracklet, tol=0.01):
    """Evaluate the complexity of the tracklets' underlying dynamics from the rank
    of their Hankel matrices, and assess whether they originate from the same track.
    The idea is that if two tracklets are part of the same track, they can be
    approximated by a low order regressor. Conversely, tracklets belonging to
    different tracks will require a higher order regressor.

    See Dicle et al., 2013.
        The Way They Move: Tracking Multiple Targets with Similar Appearance.
    """
    # TODO Add missing data imputation
    joint_tracklet = self + other_tracklet
    joint_rank = joint_tracklet.estimate_rank(tol)
    rank1 = self.estimate_rank(tol)
    rank2 = other_tracklet.estimate_rank(tol)
    return (rank1 + rank2) / joint_rank - 1

estimate_rank 

estimate_rank(tol)

Estimate the (low) rank of a noisy matrix via hard thresholding of singular values.

See Gavish & Donoho, 2013. The optimal hard threshold for singular values is 4/sqrt(3)

Source code in deeplabcut/refine_training_dataset/stitch.py 
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
def estimate_rank(self, tol):
    """Estimate the (low) rank of a noisy matrix via hard thresholding of singular
    values.

    See Gavish & Donoho, 2013.     The optimal hard threshold for singular values is
    4/sqrt(3)
    """
    mat = self.to_hankelet()
    if np.any(mat):  # check that the matrix contains non-zero entries
        # nrows, ncols = mat.shape
        # beta = nrows / ncols
        # omega = 0.56 * beta ** 3 - 0.95 * beta ** 2 + 1.82 * beta + 1.43
        _, s, _ = sli.svd(mat, min(10, min(mat.shape)))
    else:
        s = np.zeros(min(10, min(mat.shape)))

    # return np.argmin(s > omega * np.median(s))
    eigen = s**2
    diff = np.abs(np.diff(eigen / eigen[0]))
    return np.argmin(diff > tol)

immediately_follows 

immediately_follows(other_tracklet, max_gap=1)

Test whether this Tracklet follows another within a tolerance ofmax_gap frames.

Source code in deeplabcut/refine_training_dataset/stitch.py 
239
240
241
242
def immediately_follows(self, other_tracklet, max_gap=1):
    """Test whether this Tracklet follows another within a tolerance of`max_gap`
    frames."""
    return 0 < self.start - other_tracklet.end <= max_gap

motion_affinity_with 

motion_affinity_with(other_tracklet)

Evaluate the motion affinity of this Tracklet' with another one.

This evaluates whether the Tracklets could realistically be reached by one another, knowing the time separating them and their velocities. Return 0 if the Tracklets overlap.

Source code in deeplabcut/refine_training_dataset/stitch.py 
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
def motion_affinity_with(self, other_tracklet):
    """Evaluate the motion affinity of this Tracklet' with another one.

    This evaluates whether the Tracklets could realistically be reached by one
    another, knowing the time separating them and their velocities. Return 0 if the
    Tracklets overlap.
    """
    time_gap = self.time_gap_to(other_tracklet)
    if time_gap > 0:
        if self < other_tracklet:
            d1 = self.centroid[-1] + time_gap * self.calc_velocity(norm=False)
            d2 = other_tracklet.centroid[0] - time_gap * other_tracklet.calc_velocity("tail", False)
            delta1 = other_tracklet.centroid[0] - d1
            delta2 = self.centroid[-1] - d2
        else:
            d1 = other_tracklet.centroid[-1] + time_gap * other_tracklet.calc_velocity(norm=False)
            d2 = self.centroid[0] - time_gap * self.calc_velocity("tail", False)
            delta1 = self.centroid[0] - d1
            delta2 = other_tracklet.centroid[-1] - d2
        return (np.sqrt(np.sum(delta1**2)) + np.sqrt(np.sum(delta2**2))) / 2
    return 0

shape_dissimilarity_with 

shape_dissimilarity_with(other_tracklet)

Calculate the dissimilarity in shape between this Tracklet and another.

Source code in deeplabcut/refine_training_dataset/stitch.py 
294
295
296
297
298
299
300
301
302
def shape_dissimilarity_with(self, other_tracklet):
    """Calculate the dissimilarity in shape between this Tracklet and another."""
    if self in other_tracklet:
        dist = np.inf
    elif self < other_tracklet:
        dist = self.undirected_hausdorff(self.xy[-1], other_tracklet.xy[0])
    else:
        dist = self.undirected_hausdorff(self.xy[0], other_tracklet.xy[-1])
    return dist

time_gap_to 

time_gap_to(other_tracklet)

Return the time gap separating this Tracklet to another.

Source code in deeplabcut/refine_training_dataset/stitch.py 
284
285
286
287
288
289
290
291
292
def time_gap_to(self, other_tracklet):
    """Return the time gap separating this Tracklet to another."""
    if self in other_tracklet:
        t = 0
    elif self < other_tracklet:
        t = other_tracklet.start - self.end
    else:
        t = self.start - other_tracklet.end
    return t

stitch_tracklets 

stitch_tracklets(
    config_path,
    videos,
    video_extensions: str | Sequence[str] | None = None,
    shuffle=1,
    trainingsetindex=0,
    n_tracks=None,
    animal_names: list[str] | None = None,
    min_length=10,
    split_tracklets=True,
    prestitch_residuals=True,
    max_gap=None,
    weight_func=None,
    destfolder=None,
    modelprefix="",
    track_method="",
    output_name="",
    transformer_checkpoint="",
    save_as_csv=False,
    **kwargs
)

Stitch sparse tracklets into full tracks via a graph-based, minimum-cost flow optimization problem.

Parameters

config_path : str Path to the main project config.yaml file.

list

A list of strings containing the full paths to videos for analysis or a path to the directory, where all the videos with same extension are stored.

str | Sequence[str] | None, optional, default=None

Controls how videos are filtered, based on file extension. File paths and directory contents are treated differently: - None (default): file paths are accepted as-is; directories are scanned for files with a recognized video extension. - str or Sequence[str] (e.g. "mp4" or ["mp4", "avi"]): both file paths and directory contents are filtered by the given extension(s).

int, optional

An integer specifying the shuffle index of the training dataset used for training the network. The default is 1.

int, optional

Integer specifying which TrainingsetFraction to use. By default the first (note that TrainingFraction is a list in config.yaml).

int, optional

Number of tracks to reconstruct. By default, taken as the number of individuals defined in the config.yaml. Another number can be passed if the number of animals in the video is different from the number of animals the model was trained on.

list, optional

If you want the names given to individuals in the labeled data file, you can specify those names as a list here. If given and n_tracks is None, n_tracks will be set to len(animal_names). If n_tracks is not None, then it must be equal to len(animal_names). If it is not given, then animal_names will be loaded from the individuals in the project config.yaml file.

int, optional

Tracklets less than min_length frames of length are considered to be residuals; i.e., they do not participate in building the graph and finding the solution to the optimization problem, but are rather added last after "almost-complete" tracks are formed. The higher the value, the lesser the computational cost, but the higher the chance of discarding relatively long and reliable tracklets that are essential to solving the stitching task. Default is 10, and must be 3 at least.

bool, optional

By default, tracklets whose time indices are not consecutive integers are split in shorter tracklets whose time continuity is guaranteed. This is for example very powerful to get rid of tracking errors (e.g., identity switches) which are often signaled by a missing time frame at the moment they occur. Note though that for long occlusions where tracker re-identification capability can be trusted, setting split_tracklets to False is preferable.

bool, optional

Residuals will by default be grouped together according to their temporal proximity prior to being added back to the tracks. This is done to improve robustness and simultaneously reduce complexity.

int, optional

Maximal temporal gap to allow between a pair of tracklets. This is automatically determined by the TrackletStitcher by default.

callable, optional

Function accepting two tracklets as arguments and returning a scalar that must be inversely proportional to the likelihood that the tracklets belong to the same track; i.e., the higher the confidence that the tracklets should be stitched together, the lower the returned value.

string, optional

Specifies the destination folder for analysis data (default is the path of the video). Note that for subsequent analysis this folder also needs to be passed.

string, optional

Specifies the tracker used to generate the pose estimation data. For multiple animals, must be either 'box', 'skeleton', or 'ellipse' and will be taken from the config.yaml file if none is given.

str, optional

Name of the output h5 file. By default, tracks are automatically stored into the same directory as the pickle file and with its name.

bool, optional

Whether to write the tracks to a CSV file too (False by default).

additional arguments.

For torch-based shuffles, can be used to specify: - snapshot_index - detector_snapshot_index

Returns

A TrackletStitcher object

Source code in deeplabcut/refine_training_dataset/stitch.py 
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
1203
@renamed_parameter(old="videotype", new="video_extensions", since="3.0.0")
def stitch_tracklets(
    config_path,
    videos,
    video_extensions: str | Sequence[str] | None = None,
    shuffle=1,
    trainingsetindex=0,
    n_tracks=None,
    animal_names: list[str] | None = None,
    min_length=10,
    split_tracklets=True,
    prestitch_residuals=True,
    max_gap=None,
    weight_func=None,
    destfolder=None,
    modelprefix="",
    track_method="",
    output_name="",
    transformer_checkpoint="",
    save_as_csv=False,
    **kwargs,
):
    """Stitch sparse tracklets into full tracks via a graph-based, minimum-cost flow
    optimization problem.

    Parameters
    ----------
    config_path : str
        Path to the main project config.yaml file.

    videos : list
        A list of strings containing the full paths to videos for analysis or a path to the directory, where all the
        videos with same extension are stored.

    video_extensions : str | Sequence[str] | None, optional, default=None
        Controls how ``videos`` are filtered, based on file extension.
        File paths and directory contents are treated differently:
        - ``None`` (default): file paths are accepted as-is; directories are
          scanned for files with a recognized video extension.
        - ``str`` or ``Sequence[str]`` (e.g. ``"mp4"`` or ``["mp4", "avi"]``):
          both file paths and directory contents are filtered by the given
          extension(s).

    shuffle: int, optional
        An integer specifying the shuffle index of the training dataset used for training the network. The default is 1.

    trainingsetindex: int, optional
        Integer specifying which TrainingsetFraction to use. By default the first (note that TrainingFraction is a list
        in config.yaml).

    n_tracks : int, optional
        Number of tracks to reconstruct. By default, taken as the number
        of individuals defined in the config.yaml. Another number can be
        passed if the number of animals in the video is different from
        the number of animals the model was trained on.

    animal_names: list, optional
        If you want the names given to individuals in the labeled data file, you can
        specify those names as a list here. If given and `n_tracks` is None, `n_tracks`
        will be set to `len(animal_names)`. If `n_tracks` is not None, then it must be
        equal to `len(animal_names)`. If it is not given, then `animal_names` will
        be loaded from the `individuals` in the project config.yaml file.

    min_length : int, optional
        Tracklets less than `min_length` frames of length
        are considered to be residuals; i.e., they do not participate
        in building the graph and finding the solution to the
        optimization problem, but are rather added last after
        "almost-complete" tracks are formed. The higher the value,
        the lesser the computational cost, but the higher the chance of
        discarding relatively long and reliable tracklets that are
        essential to solving the stitching task.
        Default is 10, and must be 3 at least.

    split_tracklets : bool, optional
        By default, tracklets whose time indices are not consecutive integers
        are split in shorter tracklets whose time continuity is guaranteed.
        This is for example very powerful to get rid of tracking errors
        (e.g., identity switches) which are often signaled by a missing
        time frame at the moment they occur. Note though that for long
        occlusions where tracker re-identification capability can be trusted,
        setting `split_tracklets` to False is preferable.

    prestitch_residuals : bool, optional
        Residuals will by default be grouped together according to their
        temporal proximity prior to being added back to the tracks.
        This is done to improve robustness and simultaneously reduce complexity.

    max_gap : int, optional
        Maximal temporal gap to allow between a pair of tracklets.
        This is automatically determined by the TrackletStitcher by default.

    weight_func : callable, optional
        Function accepting two tracklets as arguments and returning a scalar
        that must be inversely proportional to the likelihood that the tracklets
        belong to the same track; i.e., the higher the confidence that the
        tracklets should be stitched together, the lower the returned value.

    destfolder: string, optional
        Specifies the destination folder for analysis data (default is the path of the
        video). Note that for subsequent analysis this folder also needs to be passed.

    track_method: string, optional
         Specifies the tracker used to generate the pose estimation data.
         For multiple animals, must be either 'box', 'skeleton', or 'ellipse'
         and will be taken from the config.yaml file if none is given.

    output_name : str, optional
        Name of the output h5 file.
        By default, tracks are automatically stored into the same directory
        as the pickle file and with its name.

    save_as_csv: bool, optional
        Whether to write the tracks to a CSV file too (False by default).

    kwargs: additional arguments.
        For torch-based shuffles, can be used to specify:
            - snapshot_index
            - detector_snapshot_index

    Returns
    -------
    A TrackletStitcher object
    """
    vids = collect_video_paths(videos, extensions=video_extensions)
    if not vids:
        print("No video(s) found. Please check your path!")
        return

    cfg = auxiliaryfunctions.read_config(config_path)
    track_method = auxfun_multianimal.get_track_method(cfg, track_method=track_method)
    if track_method == "ctd":
        raise ValueError(
            "CTD tracking occurs directly during video analysis. No need to call "
            "`stitch_tracklets` with `track_method=='ctd'`."
        )

    if animal_names is None:
        animal_names = cfg["individuals"]
    elif n_tracks is not None and n_tracks != len(animal_names):
        raise ValueError(
            "When setting both `n_tracks` and `animal_names`, `n_tracks` must be equal "
            f"to len(animal_names)`. Found `n_tracks`={n_tracks} and `animal_names`="
            f"{animal_names} of length {len(animal_names)}.`"
        )

    if n_tracks is None:
        n_tracks = len(animal_names)

    DLCscorer, _ = deeplabcut.utils.auxiliaryfunctions.get_scorer_name(
        cfg,
        shuffle,
        cfg["TrainingFraction"][trainingsetindex],
        modelprefix=modelprefix,
        **kwargs,
    )

    if transformer_checkpoint:
        from deeplabcut.pose_tracking_pytorch import inference

        dlctrans = inference.DLCTrans(checkpoint=transformer_checkpoint)

    def trans_weight_func(tracklet1, tracklet2, nframe, feature_dict):
        zfill_width = int(np.ceil(np.log10(nframe)))
        if tracklet1 < tracklet2:
            ind_img1 = tracklet1.inds[-1]
            coord1 = tracklet1.data[-1][:, :2]
            ind_img2 = tracklet2.inds[0]
            coord2 = tracklet2.data[0][:, :2]
        else:
            ind_img2 = tracklet2.inds[-1]
            ind_img1 = tracklet1.inds[0]
            coord2 = tracklet2.data[-1][:, :2]
            coord1 = tracklet1.data[0][:, :2]
        t1 = (coord1, ind_img1)
        t2 = (coord2, ind_img2)

        dist = dlctrans(t1, t2, zfill_width, feature_dict)
        dist = (dist + 1) / 2

        return -dist

    base_weight_func = weight_func
    for video in vids:
        print("Processing... ", video)
        nframe = len(VideoWriter(video))
        videofolder = str(Path(video).parents[0])
        dest = destfolder or videofolder
        deeplabcut.utils.auxiliaryfunctions.attempt_to_make_folder(dest)
        vname = Path(video).stem

        feature_dict_path = os.path.join(dest, vname + DLCscorer + "_bpt_features.pickle")
        # should only exist one
        if transformer_checkpoint:
            import dbm

            try:
                feature_dict = shelve.open(feature_dict_path, flag="r")
            except dbm.error as err:
                raise FileNotFoundError(f"{feature_dict_path} does not exist. Did you run transformer_reID()?") from err

        dataname = os.path.join(dest, vname + DLCscorer + ".h5")

        method = TRACK_METHODS[track_method]
        pickle_file = dataname.split(".h5")[0] + f"{method}.pickle"
        try:
            stitcher = TrackletStitcher.from_pickle(
                pickle_file, n_tracks, min_length, split_tracklets, prestitch_residuals
            )
            current_weight_func = base_weight_func
            with_id = any(tracklet.identity != -1 for tracklet in stitcher)
            if with_id and weight_func is None:
                # Add in identity weighing before building the graph
                def current_weight_func(t1, t2, stitcher=stitcher):
                    w = 0.01 if t1.identity == t2.identity else 1
                    return w * stitcher.calculate_edge_weight(t1, t2)

            if transformer_checkpoint:
                stitcher.build_graph(
                    max_gap=max_gap,
                    weight_func=partial(trans_weight_func, nframe=nframe, feature_dict=feature_dict),
                )
            else:
                stitcher.build_graph(max_gap=max_gap, weight_func=current_weight_func)

            stitcher.stitch()
            if transformer_checkpoint:
                stitcher.write_tracks(
                    output_name=output_name,
                    animal_names=animal_names,
                    suffix="tr",
                    save_as_csv=save_as_csv,
                )
            else:
                stitcher.write_tracks(
                    output_name=output_name,
                    animal_names=animal_names,
                    suffix="",
                    save_as_csv=save_as_csv,
                )
        except FileNotFoundError as e:
            print(e, "\nSkipping...")