KTS [Kernel Temporal Segmentation]
视频分割
python 2.X 版本
KTS 2.X Code
KTS 3.X Code
import numpy as np
from cpd_nonlin import cpd_nonlin
def cpd_auto(K, ncp, vmax, desc_rate=1, **kwargs):
"""Main interface
Detect change points automatically selecting their number
K - kernel between each pair of frames in video
ncp - maximum ncp
vmax - special parameter
Optional arguments:
lmin - minimum segment length
lmax - maximum segment length
desc_rate - rate of descriptor sampling (vmax always corresponds to 1x)
Note:
- cps are always calculated in subsampled coordinates irrespective to
desc_rate
- lmin and m should be in agreement
---
Returns: (cps, costs)
cps - best selected change-points
costs - costs for 0,1,2,...,m change-points
Memory requirement: ~ (3*N*N + N*ncp)*4 bytes ~= 16 * N^2 bytes
That is 1,6 Gb for the N=10000.
"""
m = ncp
(_, scores) = cpd_nonlin(K, m, backtrack=False, **kwargs)
N = K.shape[0]
N2 = N*desc_rate # length of the video before subsampling
penalties = np.zeros(m+1)
# Prevent division by zero (in case of 0 changes)
ncp = np.arange(1, m+1)
penalties[1:] = (vmax*ncp/(2.0*N2))*(np.log(float(N2)/ncp)+1)
costs = scores/float(N) + penalties
m_best = np.argmin(costs)
(cps, scores2) = cpd_nonlin(K, m_best, **kwargs)
return (cps, costs)
# ------------------------------------------------------------------------------
# Extra functions (currently not used)
def estimate_vmax(K_stable):
"""K_stable - kernel between all frames of a stable segment"""
n = K_stable.shape[0]
vmax = np.trace(centering(K_stable)/n)
return vmax
def centering(K):
"""Apply kernel centering"""
mean_rows = np.mean(K, 1)[:, np.newaxis]
return K - mean_rows - mean_rows.T + np.mean(mean_rows)
def eval_score(K, cps):
""" Evaluate unnormalized empirical score
(sum of kernelized scatters) for the given change-points """
N = K.shape[0]
cps = [0] + list(cps) + [N]
V1 = 0
V2 = 0
for i in range(len(cps)-1):
K_sub = K[cps[i]:cps[i+1], :][:, cps[i]:cps[i+1]]
V1 += np.sum(np.diag(K_sub))
V2 += np.sum(K_sub) / float(cps[i+1] - cps[i])
return (V1 - V2)
def eval_cost(K, cps, score, vmax):
""" Evaluate cost function for automatic number of change points selection
K - kernel between all frames
cps - selected change-points
score - unnormalized empirical score (sum of kernelized scatters)
vmax - vmax parameter"""
N = K.shape[0]
penalty = (vmax*len(cps)/(2.0*N))*(np.log(float(N)/len(cps))+1)
return score/float(N) + penalty
cpd_nolin.py
import numpy as np
import weave
def calc_scatters(K):
"""
Calculate scatter matrix:
scatters[i,j] = {scatter of the sequence with starting frame i and ending frame j}
"""
n = K.shape[0]
K1 = np.cumsum([0] + list(np.diag(K)))
K2 = np.zeros((n+1, n+1))
K2[1:, 1:] = np.cumsum(np.cumsum(K, 0), 1); # TODO: use the fact that K - symmetric
scatters = np.zeros((n, n));
code = r"""
for (int i = 0; i < n; i++) {
for (int j = i; j < n; j++) {
scatters(i,j) = K1(j+1)-K1(i) - (K2(j+1,j+1)+K2(i,i)-K2(j+1,i)-K2(i,j+1))/(j-i+1);
}
}
"""
weave.inline(code, ['K1','K2','scatters','n'], global_dict =
{'K1':K1, 'K2':K2, 'scatters':scatters, 'n':n}, type_converters=weave.converters.blitz)
return scatters
def cpd_nonlin(K, ncp, lmin=1, lmax=100000, backtrack=True, verbose=True,
out_scatters=None):
""" Change point detection with dynamic programming
K - square kernel matrix
ncp - number of change points to detect (ncp >= 0)
lmin - minimal length of a segment
lmax - maximal length of a segment
backtrack - when False - only evaluate objective scores (to save memory)
Returns: (cps, obj)
cps - detected array of change points: mean is thought to be constant on [ cps[i], cps[i+1] )
obj_vals - values of the objective function for 0..m changepoints
"""
m = int(ncp) # prevent numpy.int64
(n, n1) = K.shape
assert(n == n1), "Kernel matrix awaited."
assert(n >= (m + 1)*lmin)
assert(n <= (m + 1)*lmax)
assert(lmax >= lmin >= 1)
if verbose:
#print "n =", n
print "Precomputing scatters..."
J = calc_scatters(K)
if out_scatters != None:
out_scatters[0] = J
if verbose:
print "Inferring best change points..."
# I[k, l] - value of the objective for k change-points and l first frames
I = 1e101*np.ones((m+1, n+1))
I[0, lmin:lmax] = J[0, lmin-1:lmax-1]
if backtrack:
# p[k, l] --- "previous change" --- best t[k] when t[k+1] equals l
p = np.zeros((m+1, n+1), dtype=int)
else:
p = np.zeros((1,1), dtype=int)
code = r"""
#define max(x,y) ((x)>(y)?(x):(y))
for (int k=1; k'm':m, 'n':n, 'p':p, 'I':I, 'J':J,
'lmin':lmin, 'lmax':lmax, 'backtrack': int(1) if backtrack else int(0)},
type_converters=weave.converters.blitz)
# Collect change points
cps = np.zeros(m, dtype=int)
if backtrack:
cur = n
for k in range(m, 0, -1):
cps[k-1] = p[k, cur]
cur = cps[k-1]
scores = I[:, n].copy()
scores[scores > 1e99] = np.inf
return cps, scores
demo.py
import numpy as np
from cpd_nonlin import cpd_nonlin
from cpd_auto import cpd_auto
def gen_data(n, m, d=1):
"""Generates data with change points
n - number of samples
m - number of change-points
WARN: sigma is proportional to m
Returns:
X - data array (n X d)
cps - change-points array, including 0 and n"""
np.random.seed(1)
# Select changes at some distance from the boundaries
cps = np.random.permutation((n*3/4)-1)[0:m] + 1 + n/8
cps = np.sort(cps)
cps = [0] + list(cps) + [n]
mus = np.random.rand(m+1, d)*(m/2) # make sigma = m/2
X = np.zeros((n, d))
for k in range(m+1):
X[cps[k]:cps[k+1], :] = mus[k, :][np.newaxis, :] + np.random.rand(cps[k+1]-cps[k], d)
return (X, np.array(cps))
if __name__ == "__main__":
from matplotlib import pyplot as plt
plt.ioff()
print "Test 1: 1-dimensional signal"
plt.figure("Test 1: 1-dimensional signal")
n = 1000
m = 10
(X, cps_gt) = gen_data(n, m)
print "Ground truth:", cps_gt
plt.plot(X)
K = np.dot(X, X.T)
cps, scores = cpd_nonlin(K, m, lmin=1, lmax=10000)
print "Estimated:", cps
mi = np.min(X)
ma = np.max(X)
for cp in cps:
plt.plot([cp, cp], [mi, ma], 'r')
plt.show()
print "="*79
print "Test 2: multidimensional signal"
plt.figure("Test 2: multidimensional signal")
n = 1000
m = 20
(X, cps_gt) = gen_data(n, m, d=50)
print "Ground truth:", cps_gt
plt.plot(X)
K = np.dot(X, X.T)
cps, scores = cpd_nonlin(K, m, lmin=1, lmax=10000)
print "Estimated:", cps
mi = np.min(X)
ma = np.max(X)
for cp in cps:
plt.plot([cp, cp], [mi, ma], 'r')
plt.show()
print "="*79
print "Test 3: automatic selection of the number of change-points"
plt.figure("Test 3: automatic selection of the number of change-points")
(X, cps_gt) = gen_data(n, m)
print "Ground truth: (m=%d)" % m, cps_gt
plt.plot(X)
K = np.dot(X, X.T)
cps, scores = cpd_auto(K, 2*m, 1)
print "Estimated: (m=%d)" % len(cps), cps
mi = np.min(X)
ma = np.max(X)
for cp in cps:
plt.plot([cp, cp], [mi, ma], 'r')
plt.show()
print "="*79
