view novelty.py @ 18:b4bf37f94e92

prepared to add another annotation

#!/usr/bin/env python
# encoding: utf-8
"""
novelty.py

Created by mi tian on 2015-04-02.
Copyright (c) 2015 __MyCompanyName__. All rights reserved.
"""

import sys, os
import numpy as np
from scipy.signal import correlate2d, convolve2d
import matplotlib.pyplot as plt
from utils.PeakPickerUtil import PeakPicker

	
peak_picker = PeakPicker()
peak_picker.params.alpha = 9.0 # Alpha norm
peak_picker.params.delta = 0.0 # Adaptive thresholding delta
peak_picker.params.QuadThresh_a = (100 - 20.0) / 1000.0
peak_picker.params.QuadThresh_b = 0.0
peak_picker.params.QuadThresh_c = (100 - 20.0) / 1500.0
peak_picker.params.rawSensitivity = 20
peak_picker.params.aCoeffs = [1.0000, -0.5949, 0.2348]
peak_picker.params.bCoeffs = [0.1600,	 0.3200, 0.1600]
peak_picker.params.preWin = 5
peak_picker.params.postWin = 5 + 1
peak_picker.params.LP_on = True
peak_picker.params.Medfilt_on = True
peak_picker.params.Polyfit_on = True
peak_picker.params.isMedianPositive = False
peak_picker.params.gt = None

def getNoveltyCurve(ssm, kernel_size, normalise=False):
	'''Return novelty score from ssm.'''

	kernel_size = int(np.floor(kernel_size/2.0) + 1)
	stripe = getDiagonalSlice(ssm, kernel_size)
	kernel = gaussian_kernel(kernel_size)
	xc = convolve2d(stripe,kernel,mode='same')
	xc[abs(xc)>1e+10]=0.00001			
	
	novelty = xc[int(np.floor(xc.shape[0]/2.0)),:]	
	novelty = [0.0 if (np.isnan(x) or np.isinf(x) or x > 1e+100) else x for x in novelty]
	
	if normalise:
		novelty = (novelty - np.min(novelty)) / (np.max(novelty) - np.min(novelty))
	return novelty	

def getNoveltyDiff(novelty, N=1, relative=False):
	'''Return the second order differece in the novelty curve.'''
	
	diff = np.zeros_like(novelty)
	diff[:-N] = np.diff(novelty, n=N)
	
	if relative:
		diff /= novelty 
		
	return diff
	
def getDiagonalSlice(ssm, width):
	''' Return a diagonal stripe of the ssm given its width, with 45 degrees rotation. 
	Note: requres 45 degrees rotated kernel also.'''
	w = int(np.floor(width/2.0))
	length = len(np.diagonal(ssm))
	stripe = np.zeros((2*w+1,length))
	# print 'diagonal', length, w, stripe.shape
	for i in xrange(-w, w+1) :		
		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))) ))
	return stripe

def gaussian_kernel(size):
	'''Create a gaussian tapered 45 degrees rotated checkerboard kernel. 
	TODO: Unit testing: Should produce this with kernel size 3:
	0.1353	 -0.3679	0.1353
	0.3679	  1.0000	0.3679
	0.1353	 -0.3679	0.1353
	'''
	n = float(np.ceil(size / 2.0))
	kernel = np.zeros((size,size))
	for i in xrange(1,size+1) :
		for j in xrange(1,size+1) :
			gauss = np.exp( -4.0 * (np.square( (i-n)/n ) + np.square( (j-n)/n )) )
			# gauss = 1			
			if np.logical_xor( j - n > np.floor((i-n) / 2.0), j - n > np.floor((n-i) / 2.0) ) :
				kernel[i-1,j-1] = -gauss
			else:
				kernel[i-1,j-1] = gauss
				
	return kernel	
	
def process(ssm, peak_picker=peak_picker, kernel_size=48, normalise=False, plot=False):
	'''Detect segment boundaries in the ssm.'''
	# peak_picker for the 1st round boudary detection

	
	novelty = getNoveltyCurve(ssm, kernel_size, normalise=False)
	smoothed_novelty, novelty_peaks = peak_picker.process(novelty)
	
	if plot:
		plot_detection(smoothed_novelty, novelty_peaks)
	return novelty, smoothed_novelty, novelty_peaks
	
def plot_detection(smoothed_novelty, novelty_peaks):
	pass