''' Collection of fitting functions '''
"""
Author : Thomas Haslwanter
Date : Oct-2013
Version: 1.0
"""
from skimage import filter
import numpy as np
import os
import matplotlib.pyplot as plt
import scipy.stats as stats
from collections import namedtuple
[docs]def demo_ransac():
'''Find the best-fit circle in an image, using the RANSAC algorithm '''
debug_flag = 1
# Since this function is only used in demo_ransac, define it here
def drawCircle(center,r):
'''Draw a circle'''
nPts = 100
phi = np.linspace(0,2*np.pi,nPts)
x = center[0] + r*np.cos(phi)
y = center[1] + r*np.sin(phi)
plt.hold(True)
plt.plot(y,x,'r')
# Get the data
os.chdir(r'C:\Users\p20529\Coding\Matlab\imProc\ransac_fitCircle')
data = plt.imread('0_5.bmp')
# Eliminate edge artefacts
rim = 20
data = data[rim:-rim, rim:-rim]
imSize = data.shape
# Find edges
edges = filter.sobel(data)
edgePnts = edges>0.15
np.sum(edgePnts)
[x,y] = np.where(edgePnts==True)
# set RANSAC parameters
par={'eps': 3,
'ransac_threshold': 0.1,
'nIter': 500,
'lowerRadius': 5,
'upperRadius': 200
}
# Allocate memory, for center(2D), radius, numPoints (structured array)
fitted = np.zeros(par['nIter'], \
dtype={'names':['center', 'radius', 'nPts'], 'formats':['2f8', 'f8', 'i4']})
for ii in range(par['nIter']):
# Takes 3 random points, and find the corresponding circle
randEdges = np.random.permutation(len(x))[:3]
(center, radius) = fitCircle(x[randEdges], y[randEdges])
# Eliminate very small and very large circles, and centers outside the image
if not (par['lowerRadius'] < radius < par['upperRadius'] and \
0 <= center[0] < imSize[0] and \
0 <= center[1] < imSize[1]):
continue
# Make sure a reasonable number of points lies near that circle
centerDistance = np.sqrt((x-center[0])**2 + (y-center[1])**2)
inCircle = np.where(np.abs(centerDistance-radius)<par['eps'])[0]
inPts = len(inCircle)
if inPts < par['ransac_threshold'] *4*np.pi*radius*par['eps'] or inPts < 3:
continue
# Fit a circle to all good points, and save the corresponding parameters
(center, radius) = fitCircle(x[inCircle], y[inCircle])
fitted[ii] = (center, radius, inPts)
# If you want to see these points:
if debug_flag == 1:
plt.plot(y,x,'.')
plt.hold(True)
plt.plot(y[inCircle], x[inCircle],'r.')
plt.plot(y[randEdges], x[randEdges], 'g+', ms=15)
plt.axis('equal')
plt.show()
# Sort the circles, according to number of points included
fitted = np.sort(fitted,order='nPts')
# Show the best-fitting circle
plt.imshow(data, cmap='gray', origin='lower')
drawCircle(fitted[-1]['center'], fitted[-1]['radius'])
plt.show()
[docs]def fit_circle(x,y):
'''
Determine the best-fit circle to given datapoints.
Parameters
----------
x : array (N,)
x-values.
y : array (N,)
corresponding y-values.
Returns
-------
center : array (2,)
x/y coordinates of center of the circle
radius : float
Circle radius.
'''
M = np.vstack((2*x,2*y,np.ones(len(x)))).T
(par,_,_,_) = np.linalg.lstsq(M,x**2+y**2)
center = par[:2]
radius = np.sqrt(par[2]+np.sum(center**2))
return(center, radius)
[docs]def fit_exp(tFit, yFit):
'''
Calculates best-fit parameters for the exponential decay to an offset.
Parameters
----------
tFit : array (N,)
Time values.
yFit : array (N,)
Function values
Returns
-------
offset : float
Function offset/bias.
amp : float
Amplitude of exponential function
tau : float
Decay time.
'''
from scipy import optimize
# Define the fit-function and the error-function
fitfunc = lambda p, x: p[0] + p[1]*np.exp(-x/p[2])
errfunc = lambda p,x,y: fitfunc(p,x) - y
pInit = np.r_[0, 1, 1] # Initial values
# Make the fit
pFit, success = optimize.leastsq(errfunc, pInit, args=(tFit, yFit))
# Plot the data and the fit
plt.plot(tFit, yFit, label='rawdata')
plt.hold(True)
plt.plot(tFit, fitfunc(pFit,tFit), label='fit')
plt.legend()
plt.show()
ExpFit = namedtuple('ExpFit', ['offset', 'amplitude', 'tau'])
return ExpFit(pFit)
[docs]def fit_line(x, y, alpha=0.05, newx=[], plotFlag=1):
"""
Linear regression fit.
Parameters
----------
x : ndarray
Input / Predictor.
y : ndarray
Input / Estimator.
alpha : float
Confidence limit [default=0.05]
newx : float or ndarray
Values for which the fit and the prediction limits are calculated (optional)
plotFlag: int, optional
1 = plot, 0 = no_plot [default]
Returns
-------
a : float
Intercept
b : float
Slope
ci : ndarray
Lower and upper confidence interval for the slope
info : dictionary
contains return information on
- residuals
- var_res
- sd_res
- alpha
- tval
- df
newy : list(ndarray)
Predictions for (newx, newx-ciPrediction, newx+ciPrediction)
Examples
--------
>>> import numpy as np
>>> from fitLine import fitLine
>>> x = np.r_[0:10:11j]
>>> y = x**2
>>> (a,b,(ci_a, ci_b),_)=fitLine(x,y)
Notes
-----
Example data and formulas are taken from
D. Altman, "Practical Statistics for Medicine"
"""
# Summary data
n = len(x) # number of samples
Sxx = np.sum(x**2) - np.sum(x)**2/n
# Syy = np.sum(y**2) - np.sum(y)**2/n # not needed here
Sxy = np.sum(x*y) - np.sum(x)*np.sum(y)/n
mean_x = np.mean(x)
mean_y = np.mean(y)
# Linefit
b = Sxy/Sxx
a = mean_y - b*mean_x
# Residuals
fit = lambda xx: a + b*xx
residuals = y - fit(x)
var_res = np.sum(residuals**2)/(n-2)
sd_res = np.sqrt(var_res)
# Confidence intervals
se_b = sd_res/np.sqrt(Sxx)
se_a = sd_res*np.sqrt(np.sum(x**2)/(n*Sxx))
df = n-2 # degrees of freedom
tval = stats.t.isf(alpha/2., df) # appropriate t value
ci_a = a + tval*se_a*np.array([-1,1])
ci_b = b + tval*se_b*np.array([-1,1])
# create series of new test x-values to predict for
npts = 100
px = np.linspace(np.min(x),np.max(x),num=npts)
se_fit = lambda x: sd_res * np.sqrt( 1./n + (x-mean_x)**2/Sxx)
se_predict = lambda x: sd_res * np.sqrt(1+1./n + (x-mean_x)**2/Sxx)
print('Summary: a={0:5.4f}+/-{1:5.4f}, b={2:5.4f}+/-{3:5.4f}'.format(a,tval*se_a,b,tval*se_b))
print('Confidence intervals: ci_a=({0:5.4f} - {1:5.4f}), ci_b=({2:5.4f} - {3:5.4f})'.format(ci_a[0], ci_a[1], ci_b[0], ci_b[1]))
print('Residuals: variance = {0:5.4f}, standard deviation = {1:5.4f}'.format(var_res, sd_res))
print('alpha = {0:.3f}, tval = {1:5.4f}, df={2:d}'.format(alpha, tval, df))
# Return info
ri = {'residuals': residuals,
'var_res': var_res,
'sd_res': sd_res,
'alpha': alpha,
'tval': tval,
'df': df}
if plotFlag == 1:
# Plot the data
plt.figure()
plt.plot(px, fit(px),'k', label='Regression line')
#plt.plot(x,y,'k.', label='Sample observations', ms=10)
plt.plot(x,y,'k.')
x.sort()
limit = (1-alpha)*100
plt.plot(x, fit(x)+tval*se_fit(x), 'r--', lw=2, label='Confidence limit ({0:.1f}%)'.format(limit))
plt.plot(x, fit(x)-tval*se_fit(x), 'r--', lw=2 )
plt.plot(x, fit(x)+tval*se_predict(x), '--', lw=2, color=(0.2,1,0.2), label='Prediction limit ({0:.1f}%)'.format(limit))
plt.plot(x, fit(x)-tval*se_predict(x), '--', lw=2, color=(0.2,1,0.2))
plt.xlabel('X values')
plt.ylabel('Y values')
plt.title('Linear regression and confidence limits')
# configure legend
plt.legend(loc=0)
leg = plt.gca().get_legend()
ltext = leg.get_texts()
plt.setp(ltext, fontsize=10)
# show the plot
plt.show()
if newx != []:
try:
newx.size
except AttributeError:
newx = np.array([newx])
print('Example: x = {0}+/-{1} => se_fit = {2:5.4f}, se_predict = {3:6.5f}'\
.format(newx[0], tval*se_predict(newx[0]), se_fit(newx[0]), se_predict(newx[0])))
newy = (fit(newx), fit(newx)-se_predict(newx), fit(newx)+se_predict(newx))
return (a,b,(ci_a, ci_b), ri, newy)
else:
return (a,b,(ci_a, ci_b), ri)
[docs]def fit_sin(tList,yList,freq):
'''
Fit a sine wave with a known frequency to a given set of data.
y = amplitude * sin(2*pi*freq * tList + phase*pi/180) + bias
Parameters
----------
yList : array
datapoints
tList : float
time base, in sec
freq : float
in Hz
Returns
-------
phase : float
in degrees
amplitude : float
bias : float
Examples
--------
>>> (phase, amp, offset) = thLib.fitSin.fitSine(t, data, freq)
'''
# Get input data
b = yList
rows = [ [np.sin(freq*2*np.pi*t), np.cos(freq*2*np.pi*t), 1] for t in tList]
A = np.array(rows)
# Make the fit
(w,residuals,rank,sing_vals) = np.linalg.lstsq(A,b)
# Extract desired parameters ...
phase = math.atan2(w[1],w[0])*180/np.pi
amplitude = np.linalg.norm([w[0],w[1]],2)
bias = w[2]
# ... and return them
return (phase,amplitude,bias)
[docs]def regress(x, y, alpha=0.05):
'''
Linear regression and confidence intervals, for a linear regression y = k * x + d
(kd, ci) = regress(x,y,[alpha=0.05])
Parameters
----------
x : ndarray (N,)
predictors
y : ndarray(N,)
responses
alpha: float
Defines the 100*(1-alpha)% confidence level in ci.
Returns
-------
k : float
estimated slope
d : float
estimated intercept
ci : ndarray (2,)
confidence intervals for the slope
Note
----
Similar to the MATLAB command "regress"
Examples
--------
>>> x = arange(100)
>>> y = 0.3*x + 10 + randn(len(x))
>>> (kd, ci) = thLib.fits.regress(x,y)
See also
--------
fit_line
'''
regressed = stats.linregress(x,y)
kd = regressed[0:2]
se_k = regressed[4]
level = (1.-alpha/2.)
tVal = stats.t.ppf(level,len(x)-2)
ci = kd[0] + se_k*tVal*np.array([-1, 1])
RegressionLine = namedtuple('Regression', ['slope', 'intercept', 'ci_slope'])
return RegressionLine(kd[0], kd[1], ci)
if __name__=='__main__':
print ('Done')