|
mi@0
|
1 #!/usr/bin/env python
|
|
mi@0
|
2 # encoding: utf-8
|
|
mi@0
|
3 """
|
|
mi@0
|
4 novelty.py
|
|
mi@0
|
5
|
|
mi@0
|
6 Created by mi tian on 2015-04-02.
|
|
mi@0
|
7 Copyright (c) 2015 __MyCompanyName__. All rights reserved.
|
|
mi@0
|
8 """
|
|
mi@0
|
9
|
|
mi@0
|
10 import sys, os
|
|
mi@0
|
11 import numpy as np
|
|
mi@0
|
12 from scipy.signal import correlate2d, convolve2d
|
|
mitian@1
|
13 import matplotlib.pyplot as plt
|
|
mitian@7
|
14 from utils.PeakPickerUtil import PeakPicker
|
|
mi@0
|
15
|
|
mitian@7
|
16
|
|
mitian@7
|
17 peak_picker = PeakPicker()
|
|
mitian@7
|
18 peak_picker.params.alpha = 9.0 # Alpha norm
|
|
mitian@7
|
19 peak_picker.params.delta = 0.0 # Adaptive thresholding delta
|
|
mitian@7
|
20 peak_picker.params.QuadThresh_a = (100 - 20.0) / 1000.0
|
|
mitian@7
|
21 peak_picker.params.QuadThresh_b = 0.0
|
|
mitian@7
|
22 peak_picker.params.QuadThresh_c = (100 - 20.0) / 1500.0
|
|
mitian@7
|
23 peak_picker.params.rawSensitivity = 20
|
|
mitian@7
|
24 peak_picker.params.aCoeffs = [1.0000, -0.5949, 0.2348]
|
|
mitian@7
|
25 peak_picker.params.bCoeffs = [0.1600, 0.3200, 0.1600]
|
|
mitian@7
|
26 peak_picker.params.preWin = 5
|
|
mitian@7
|
27 peak_picker.params.postWin = 5 + 1
|
|
mitian@7
|
28 peak_picker.params.LP_on = True
|
|
mitian@7
|
29 peak_picker.params.Medfilt_on = True
|
|
mitian@7
|
30 peak_picker.params.Polyfit_on = True
|
|
mitian@7
|
31 peak_picker.params.isMedianPositive = False
|
|
mitian@18
|
32 peak_picker.params.gt = None
|
|
mitian@18
|
33
|
|
mi@0
|
34 def getNoveltyCurve(ssm, kernel_size, normalise=False):
|
|
mi@0
|
35 '''Return novelty score from ssm.'''
|
|
mi@0
|
36
|
|
mi@0
|
37 kernel_size = int(np.floor(kernel_size/2.0) + 1)
|
|
mi@0
|
38 stripe = getDiagonalSlice(ssm, kernel_size)
|
|
mi@0
|
39 kernel = gaussian_kernel(kernel_size)
|
|
mi@0
|
40 xc = convolve2d(stripe,kernel,mode='same')
|
|
mi@0
|
41 xc[abs(xc)>1e+10]=0.00001
|
|
mi@0
|
42
|
|
mi@0
|
43 novelty = xc[int(np.floor(xc.shape[0]/2.0)),:]
|
|
mi@0
|
44 novelty = [0.0 if (np.isnan(x) or np.isinf(x) or x > 1e+100) else x for x in novelty]
|
|
mi@0
|
45
|
|
mi@0
|
46 if normalise:
|
|
mi@0
|
47 novelty = (novelty - np.min(novelty)) / (np.max(novelty) - np.min(novelty))
|
|
mi@0
|
48 return novelty
|
|
mi@0
|
49
|
|
mitian@8
|
50 def getNoveltyDiff(novelty, N=1, relative=False):
|
|
mitian@8
|
51 '''Return the second order differece in the novelty curve.'''
|
|
mitian@8
|
52
|
|
mitian@8
|
53 diff = np.zeros_like(novelty)
|
|
mitian@8
|
54 diff[:-N] = np.diff(novelty, n=N)
|
|
mitian@8
|
55
|
|
mitian@8
|
56 if relative:
|
|
mitian@8
|
57 diff /= novelty
|
|
mitian@8
|
58
|
|
mitian@8
|
59 return diff
|
|
mitian@8
|
60
|
|
mi@0
|
61 def getDiagonalSlice(ssm, width):
|
|
mi@0
|
62 ''' Return a diagonal stripe of the ssm given its width, with 45 degrees rotation.
|
|
mi@0
|
63 Note: requres 45 degrees rotated kernel also.'''
|
|
mi@0
|
64 w = int(np.floor(width/2.0))
|
|
mi@0
|
65 length = len(np.diagonal(ssm))
|
|
mi@0
|
66 stripe = np.zeros((2*w+1,length))
|
|
mi@0
|
67 # print 'diagonal', length, w, stripe.shape
|
|
mi@0
|
68 for i in xrange(-w, w+1) :
|
|
mi@0
|
69 stripe[w+i,:] = np.hstack(( np.zeros(int(np.floor(abs(i)/2.0))), np.diagonal(ssm,i), np.zeros(int(np.ceil(abs(i)/2.0))) ))
|
|
mi@0
|
70 return stripe
|
|
mi@0
|
71
|
|
mi@0
|
72 def gaussian_kernel(size):
|
|
mi@0
|
73 '''Create a gaussian tapered 45 degrees rotated checkerboard kernel.
|
|
mi@0
|
74 TODO: Unit testing: Should produce this with kernel size 3:
|
|
mi@0
|
75 0.1353 -0.3679 0.1353
|
|
mi@0
|
76 0.3679 1.0000 0.3679
|
|
mi@0
|
77 0.1353 -0.3679 0.1353
|
|
mi@0
|
78 '''
|
|
mi@0
|
79 n = float(np.ceil(size / 2.0))
|
|
mi@0
|
80 kernel = np.zeros((size,size))
|
|
mi@0
|
81 for i in xrange(1,size+1) :
|
|
mi@0
|
82 for j in xrange(1,size+1) :
|
|
mi@0
|
83 gauss = np.exp( -4.0 * (np.square( (i-n)/n ) + np.square( (j-n)/n )) )
|
|
mi@0
|
84 # gauss = 1
|
|
mi@0
|
85 if np.logical_xor( j - n > np.floor((i-n) / 2.0), j - n > np.floor((n-i) / 2.0) ) :
|
|
mi@0
|
86 kernel[i-1,j-1] = -gauss
|
|
mi@0
|
87 else:
|
|
mi@0
|
88 kernel[i-1,j-1] = gauss
|
|
mi@0
|
89
|
|
mi@0
|
90 return kernel
|
|
mi@0
|
91
|
|
mitian@17
|
92 def process(ssm, peak_picker=peak_picker, kernel_size=48, normalise=False, plot=False):
|
|
mi@0
|
93 '''Detect segment boundaries in the ssm.'''
|
|
mitian@7
|
94 # peak_picker for the 1st round boudary detection
|
|
mitian@7
|
95
|
|
mitian@7
|
96
|
|
mi@0
|
97 novelty = getNoveltyCurve(ssm, kernel_size, normalise=False)
|
|
mi@0
|
98 smoothed_novelty, novelty_peaks = peak_picker.process(novelty)
|
|
mi@0
|
99
|
|
mitian@1
|
100 if plot:
|
|
mitian@1
|
101 plot_detection(smoothed_novelty, novelty_peaks)
|
|
mitian@1
|
102 return novelty, smoothed_novelty, novelty_peaks
|
|
mitian@1
|
103
|
|
mitian@1
|
104 def plot_detection(smoothed_novelty, novelty_peaks):
|
|
mitian@17
|
105 pass
|
