重庆分公司,新征程启航
为企业提供网站建设、域名注册、服务器等服务
#!/usr/bin/env python
网页设计是网站建设的前奏,好的网页设计更深度的剖析产品和设计风格定位,结合最新的网页设计流行趋势,与WVI应用标准,设计出具企业表现力,大器而深稳的网站界面设。成都创新互联于2013年开始,是成都网站建设公司:提供企业网站设计,成都品牌网站建设,营销型企业网站建设方案,响应式网站设计,小程序设计,专业建站公司做网站。
# -*- coding: utf-8 -*-
# File name: parabolic
# Project name: parabolic_equation
"""
.. moduleauthor::
.. Module.. name parabolic of procjet parabolic_equation
"""
from sympy import *
import matplotlib.pyplot as plt
import numpy as np
def _filterComplex(inputvalue, description='inputvalue'):
try:
str(inputvalue).index('I')
except ValueError:
return False
else:
return True
def _checkBool(inputvalue, description='inputvalue'):
"""
:param inputvalue:
:param description:
:return:
"""
if not isinstance(inputvalue, bool):
raise TypeError(
'The {0} must be boolean. Given: {1!r}'.format(description, inputvalue))
def _checkNumerical(inputvalue, description='inputvalue'):
"""
:param inputvalue:
:param description:
:return:
"""
try:
inputvalue + 1
except TypeError:
raise TypeError(
'The {0} must be numerical. Given: {1!r}'.format(description, inputvalue))
def _drawTowPara(expr_1, expr_2, inputmin, inputmax ,step=0.1):
"""
:param expr_1:
:param expr_2:
:param inputmin:
:param inputmax:
:param step:
:param expr_1_evalwithY:
:param expr_2_evalwithY:
:return:
"""
_checkNumerical(inputmin, 'xmin')
_checkNumerical(inputmax, 'xmax')
_checkNumerical(step, 'step')
y1List = []
x1List = []
y2List = []
x2List = []
if expr_1.vertical is True:
x1List = np.arange(inputmin, inputmax, step)
for x in x1List:
y1List.append(expr_1.evaluates_Y(x))
else:
y1List = np.arange(inputmin, inputmax, step)
for y in y1List:
x1List.append(expr_1.evaluates_X(y))
if expr_2.vertical is True:
x2List = np.arange(inputmin, inputmax, step)
for x in x2List:
y2List.append(expr_2.evaluates_Y(x))
else:
y2List = np.arange(inputmin, inputmax, step)
for y in y2List:
x2List.append(expr_2.evaluates_X(y))
plt.plot(x1List, y1List, '+')
plt.plot(x2List, y2List, '-')
plt.show()
def _solveCrossing(expr_1, expr_2):
"""
:param expr_1:
:param expr_2:
:return:
"""
x = Symbol('x')
y = Symbol('y')
print "Given the first expression: {0!r}".format(expr_1.expr)
print "Given the first expression: {0!r}".format(expr_2.expr)
ResultList = solve([expr_1.expr, expr_2.expr], [x, y])
Complex = False
ResultListTrue = []
for i in range(0, (len(ResultList)),1):
if _filterComplex(ResultList[i][0], 'x') or _filterComplex(ResultList[i][1], 'y'):
Complex = True
else:
ResultListTrue.append(ResultList[i])
if len(ResultListTrue) == 0 and Complex:
print "Two hyperbolic do not intersect, and there is imaginary value."
elif len(ResultListTrue) == 1:
print "Two hyperbolic tangent.:"
print ResultListTrue
else:
print "Two hyperbolic intersection, and Points are:"
for iterm in ResultListTrue:
print iterm
class Parabolic():
"""
"""
def __init__(self, a, b, c, vertical=True):
"""
:return:
"""
_checkNumerical(a, 'a')
_checkNumerical(b, 'b')
_checkNumerical(c, 'c')
_checkBool(vertical, 'vertical')
self.a = a
self.b = b
self.c = c
self.vertical = vertical
self.y = Symbol('y')
self.x = Symbol('x')
self.xarray = []
self.yarray = []
if vertical is True:
self.expr = (self.x**2)*self.a + self.x*self.b + self.c
else:
self.expr = (self.y**2)*self.a + self.y*self.b + self.c
def __repr__(self):
"""
:return:
"""
if self.vertical is True:
return "The Equation look like: {0!r}".format(self.expr)
else:
return "The Equation look like: {0!r}".format(self.expr)
def evaluates_X(self, inputvalue):
"""
:param inputvalue:
:return:
"""
_checkNumerical(inputvalue, 'y')
return self.expr.subs(self.y, inputvalue)
def evaluates_Y(self, inputvalue):
"""
:param inputvalue:
:return:
"""
_checkNumerical(inputvalue, 'x')
return self.expr.subs(self.x, inputvalue)
def getArrays(self, inputmin, inputmax, step=1):
"""
:param inputmin:
:param inputmax:
:param step:
:return:
"""
_checkNumerical(inputmin, 'xmin')
_checkNumerical(inputmax, 'xmax')
_checkNumerical(step, 'step')
if self.vertical is True:
for x in range(inputmin, inputmax, step):
self.xarray.append(x)
self.yarray.append(self.evaluates_Y(x))
else:
for y in range(inputmin, inputmax, step):
self.yarray.append(y)
self.xarray.append(self.evaluates_X(y))
def drawPara(self, inputmin, inputmax, step=1):
"""
:param inputmin:
:param inputmax:
:param step:
:return:
"""
_checkNumerical(inputmin, 'xmin')
_checkNumerical(inputmax, 'xmax')
_checkNumerical(step, 'step')
yList = []
xList = []
if self.vertical is True:
xList = np.arange(inputmin, inputmax, step)
for x in xList:
yList.append(self.evaluates_Y(x))
else:
yList = np.arange(inputmin, inputmax, step)
for y in yList:
xList.append(self.evaluates_X(y))
plt.plot(xList, yList, '+')
plt.show()
if __name__ == '__main__':
pa1 = Parabolic(-5,3,6)
pa2 = Parabolic(-5,2,5, False)
print pa1
print pa2
_solveCrossing(pa1, pa2)
_drawTowPara(pa1, pa2, -10, 10, 0.1)
# 这就是你想要的,代码解决了你的大部分问题,可以求两条双曲线交点,或者直线与双曲线交#点,或者两直线交点. 不过定义双曲线时候使用的是一般式.也也尽可能做了测试,如果有#问题的话,追问吧
最优化
为什么要做最优化呢?因为在生活中,人们总是希望幸福值或其它达到一个极值,比如做生意时希望成本最小,收入最大,所以在很多商业情境中,都会遇到求极值的情况。
函数求根
这里「函数的根」也称「方程的根」,或「函数的零点」。
先把我们需要的包加载进来。import numpy as npimport scipy as spimport scipy.optimize as optimport matplotlib.pyplot as plt%matplotlib inline
函数求根和最优化的关系?什么时候函数是最小值或最大值?
两个问题一起回答:最优化就是求函数的最小值或最大值,同时也是极值,在求一个函数最小值或最大值时,它所在的位置肯定是导数为 0 的位置,所以要求一个函数的极值,必然要先求导,使其为 0,所以函数求根就是为了得到最大值最小值。
scipy.optimize 有什么方法可以求根?
可以用 scipy.optimize 中的 bisect 或 brentq 求根。f = lambda x: np.cos(x) - x # 定义一个匿名函数x = np.linspace(-5, 5, 1000) # 先生成 1000 个 xy = f(x) # 对应生成 1000 个 f(x)plt.plot(x, y); # 看一下这个函数长什么样子plt.axhline(0, color='k'); # 画一根横线,位置在 y=0
opt.bisect(f, -5, 5) # 求取函数的根0.7390851332155535plt.plot(x, y)plt.axhline(0, color='k')plt.scatter([_], [0], c='r', s=100); # 这里的 [_] 表示上一个 Cell 中的结果,这里是 x 轴上的位置,0 是 y 上的位置
求根有两种方法,除了上面介绍的 bisect,还有 brentq,后者比前者快很多。%timeit opt.bisect(f, -5, 5)%timeit opt.brentq(f, -5, 5)10000 loops, best of 3: 157 s per loopThe slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.10000 loops, best of 3: 35.9 s per loop
函数求最小化
求最小值就是一个最优化问题。求最大值时只需对函数做一个转换,比如加一个负号,或者取倒数,就可转成求最小值问题。所以两者是同一问题。
初始值对最优化的影响是什么?
举例来说,先定义个函数。f = lambda x: 1-np.sin(x)/xx = np.linspace(-20., 20., 1000)y = f(x)
当初始值为 3 值,使用 minimize 函数找到最小值。minimize 函数是在新版的 scipy 里,取代了以前的很多最优化函数,是个通用的接口,背后是很多方法在支撑。x0 = 3xmin = opt.minimize(f, x0).x # x0 是起始点,起始点最好离真正的最小值点不要太远plt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300); # 起始点画出来,用圆圈表示plt.scatter(xmin, f(xmin), marker='v', s=300); # 最小值点画出来,用三角表示plt.xlim(-20, 20);
初始值为 3 时,成功找到最小值。
现在来看看初始值为 10 时,找到的最小值点。x0 = 10xmin = opt.minimize(f, x0).xplt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300)plt.scatter(xmin, f(xmin), marker='v', s=300)plt.xlim(-20, 20);
由上图可见,当初始值为 10 时,函数找到的是局部最小值点,可见 minimize 的默认算法对起始点的依赖性。
那么怎么才能不管初始值在哪个位置,都能找到全局最小值点呢?
如何找到全局最优点?
可以使用 basinhopping 函数找到全局最优点,相关背后算法,可以看帮助文件,有提供论文的索引和出处。
我们设初始值为 10 看是否能找到全局最小值点。x0 = 10from scipy.optimize import basinhoppingxmin = basinhopping(f,x0,stepsize = 5).xplt.plot(x, y);plt.scatter(x0, f(x0), marker='o', s=300);plt.scatter(xmin, f(xmin), marker='v', s=300);plt.xlim(-20, 20);
当起始点在比较远的位置,依然成功找到了全局最小值点。
如何求多元函数最小值?
以二元函数为例,使用 minimize 求对应的最小值。def g(X): x,y = X return (x-1)**4 + 5 * (y-1)**2 - 2*x*yX_opt = opt.minimize(g, (8, 3)).x # (8,3) 是起始点print X_opt[ 1.88292611 1.37658521]fig, ax = plt.subplots(figsize=(6, 4)) # 定义画布和图形x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, g((X, Y)), 50) # 等高线图ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15) # 最小点的位置是个元组ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax) # colorbar 表示颜色越深,高度越高fig.tight_layout()
画3D 图。from mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmfig = plt.figure()ax = fig.gca(projection='3d')x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)surf = ax.plot_surface(X, Y, g((X,Y)), rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)cset = ax.contour(X, Y, g((X,Y)), zdir='z',offset=-5, cmap=cm.coolwarm)fig.colorbar(surf, shrink=0.5, aspect=5);
曲线拟合
曲线拟合和最优化有什么关系?
曲线拟合的问题是,给定一组数据,它可能是沿着一条线散布的,这时要找到一条最优的曲线来拟合这些数据,也就是要找到最好的线来代表这些点,这里的最优是指这些点和线之间的距离是最小的,这就是为什么要用最优化问题来解决曲线拟合问题。
举例说明,给一些点,找到一条线,来拟合这些点。
先给定一些点:N = 50 # 点的个数m_true = 2 # 斜率b_true = -1 # 截距dy = 2.0 # 误差np.random.seed(0)xdata = 10 * np.random.random(N) # 50 个 x,服从均匀分布ydata = np.random.normal(b_true + m_true * xdata, dy) # dy 是标准差plt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');
上面的点整体上呈现一个线性关系,要找到一条斜线来代表这些点,这就是经典的一元线性回归。目标就是找到最好的线,使点和线的距离最短。要优化的函数是点和线之间的距离,使其最小。点是确定的,而线是可变的,线是由参数值,斜率和截距决定的,这里就是要通过优化距离找到最优的斜率和截距。
点和线的距离定义如下:def chi2(theta, x, y): return np.sum(((y - theta[0] - theta[1] * x)) ** 2)
上式就是误差平方和。
误差平方和是什么?有什么作用?
误差平方和公式为:
误差平方和大,表示真实的点和预测的线之间距离太远,说明拟合得不好,最好的线,应该是使误差平方和最小,即最优的拟合线,这里是条直线。
误差平方和就是要最小化的目标函数。
找到最优的函数,即斜率和截距。theta_guess = [0, 1] # 初始值theta_best = opt.minimize(chi2, theta_guess, args=(xdata, ydata)).xprint(theta_best)[-1.01442005 1.93854656]
上面两个输出即是预测的直线斜率和截距,我们是根据点来反推直线的斜率和截距,那么真实的斜率和截距是多少呢?-1 和 2,很接近了,差的一点是因为有噪音的引入。xfit = np.linspace(0, 10)yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
最小二乘(Least Square)是什么?
上面用的是 minimize 方法,这个问题的目标函数是误差平方和,这就又有一个特定的解法,即最小二乘。
最小二乘的思想就是要使得观测点和估计点的距离的平方和达到最小,这里的“二乘”指的是用平方来度量观测点与估计点的远近(在古汉语中“平方”称为“二乘”),“最小”指的是参数的估计值要保证各个观测点与估计点的距离的平方和达到最小。
关于最小二乘估计的计算,涉及更多的数学知识,这里不想详述,其一般的过程是用目标函数对各参数求偏导数,并令其等于 0,得到一个线性方程组。具体推导过程可参考斯坦福机器学习讲义 第 7 页。def deviations(theta, x, y): return (y - theta[0] - theta[1] * x)theta_best, ier = opt.leastsq(deviations, theta_guess, args=(xdata, ydata))print(theta_best)[-1.01442016 1.93854659]
最小二乘 leastsq 的结果跟 minimize 结果一样。注意 leastsq 的第一个参数不再是误差平方和 chi2,而是误差本身 deviations,即没有平方,也没有和。yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
非线性最小二乘
上面是给一些点,拟合一条直线,拟合一条曲线也是一样的。def f(x, beta0, beta1, beta2): # 首先定义一个非线性函数,有 3 个参数 return beta0 + beta1 * np.exp(-beta2 * x**2)beta = (0.25, 0.75, 0.5) # 先猜 3 个 betaxdata = np.linspace(0, 5, 50)y = f(xdata, *beta)ydata = y + 0.05 * np.random.randn(len(xdata)) # 给 y 加噪音def g(beta): return ydata - f(xdata, *beta) # 真实 y 和 预测值的差,求最优曲线时要用到beta_start = (1, 1, 1)beta_opt, beta_cov = opt.leastsq(g, beta_start)print beta_opt # 求到的 3 个最优的 beta 值[ 0.25525709 0.74270226 0.54966466]
拿估计的 beta_opt 值跟真实的 beta = (0.25, 0.75, 0.5) 值比较,差不多。fig, ax = plt.subplots()ax.scatter(xdata, ydata) # 画点ax.plot(xdata, y, 'r', lw=2) # 真实值的线ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2) # 拟合的线ax.set_xlim(0, 5)ax.set_xlabel(r"$x$", fontsize=18)ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)fig.tight_layout()
除了使用最小二乘,还可以使用曲线拟合的方法,得到的结果是一样的。beta_opt, beta_cov = opt.curve_fit(f, xdata, ydata)print beta_opt[ 0.25525709 0.74270226 0.54966466]
有约束的最小化
有约束的最小化是指,要求函数最小化之外,还要满足约束条件,举例说明。
边界约束def f(X): x, y = X return (x-1)**2 + (y-1)**2 # 这是一个碗状的函数x_opt = opt.minimize(f, (0, 0), method='BFGS').x # 无约束最优化
假设有约束条件,x 和 y 要在一定的范围内,如 x 在 2 到 3 之间,y 在 0 和 2 之间。bnd_x1, bnd_x2 = (2, 3), (0, 2) # 对自变量的约束x_cons_opt = opt.minimize(f, np.array([0, 0]), method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x # bounds 矩形约束fig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X,Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 没有约束下的最小值,蓝色五角星ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 有约束下的最小值,红色星星bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]), bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0], facecolor="grey")ax.add_patch(bound_rect)ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
不等式约束
介绍下相关理论,先来看下存在等式约束的极值问题求法,比如下面的优化问题。
目标函数是 f(w),下面是等式约束,通常解法是引入拉格朗日算子,这里使用 ββ 来表示算子,得到拉格朗日公式为
l 是等式约束的个数。
然后分别对 w 和ββ 求偏导,使得偏导数等于 0,然后解出 w 和βiβi,至于为什么引入拉格朗日算子可以求出极值,原因是 f(w) 的 dw 变化方向受其他不等式的约束,dw的变化方向与f(w)的梯度垂直时才能获得极值,而且在极值处,f(w) 的梯度与其他等式梯度的线性组合平行,因此他们之间存在线性关系。(参考《最优化与KKT条件》)
对于不等式约束的极值问题
常常利用拉格朗日对偶性将原始问题转换为对偶问题,通过解对偶问题而得到原始问题的解。该方法应用在许多统计学习方法中。有兴趣的可以参阅相关资料,这里不再赘述。def f(X): return (X[0] - 1)**2 + (X[1] - 1)**2def g(X): return X[1] - 1.75 - (X[0] - 0.75)**4x_opt = opt.minimize(f, (0, 0), method='BFGS').xconstraints = [dict(type='ineq', fun=g)] # 约束采用字典定义,约束方式为不等式约束,边界用 g 表示x_cons_opt = opt.minimize(f, (0, 0), method='SLSQP', constraints=constraints).xfig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X, Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 蓝色星星,没有约束下的最小值ax.plot(x_, 1.75 + (x_-0.75)**4, '', markersize=15)ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color="grey")ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 在区域约束下的最小值ax.set_ylim(-1, 3)ax.set_xlabel(r"$x_0$", fontsize=18)ax.set_ylabel(r"$x_1$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
scipy.optimize.minimize 中包括了多种最优化算法,每种算法使用范围不同,详细参考官方文档。
# -*- coding: UTF-8 -*-
import matplotlib.pyplot as plt
import numpy as np
#生成x的等差数列0-10之间取100个数
x = np.linspace(0,10, 100)
#生成每个x对应的y
y = 0.5*x+3
#画直线
plt.plot(x, y, c='orange')
#画标题
plt.title("y=0.5x+3")
#显示
plt.show()
你好,既然你知道怎么pylab
画图的话,
那么画斜率的不是也一样的吗?
用斜率的公式,先计算出来,
然后传进函数里面,
你可以里面subplot,显示在同一个界面上。
霍夫变换(Hough Transform)是图像处理领域中,从图像中识别几何形状的基本方法之一。主要识别具有某些相同特征的几何形状,例如直线,圆形,本篇博客的目标就是从黑白图像中识别出直线。
翻阅霍夫直线变换的原理时候,橡皮擦觉得原理部分需要先略过,否则很容易在这个地方陷进去,但是问题来了,这个原理略过了,直接应用函数,里面有些参数竟然看不懂。例如极坐标,角度扫描范围,这种函数就属于绕不过去的知识点了,所以本文转移方向,死磕原理,下面的博文将语无伦次的为你展示如何学习原理知识。
因为数学知识的贫乏,所以在学习阶段会涉及到很多基础概念的学习,一起来吧。
首先找到相对官方的资料,打开该 地址
下面是一个数学小白对原理的学习经验。
教材说:众所周知,一条直线在图像二维空间可由两个变量表示。
抱歉,小白还真不知道……即使学习过,这些年也早已经还给老师了。
一开始难道要学习笛卡尔坐标系,不,你低估小白的能力了,我第一个查询的是 θ 读作 西塔 ,是一个希腊字母。
什么是笛卡尔坐标系?
这个比较简单,直角坐标系。
斜率和截距
斜率,亦称“角系数”,表示一条直线相对于横坐标轴的倾斜程度。
一条直线与某平面直角坐标系横坐标轴正半轴方向的夹角的正切值即该直线相对于该坐标系的斜率。
如果直线与 x 轴互相垂直,直角的正切直无穷大,故此直线不存在斜率。
对于一次函数 y=kx+b , k 就是该函数图像的斜率。
在学习的时候,也学到如下内容:
截距:对 x 的截距就是 y=0 时, x 的值,对 y 的截距就是 x=0 时, y 的值,
截距就是直线与坐标轴的交点的横(纵)坐标。 x 截距为 a , y 截距 b ,截距式就是: x/a+y/b=1(a≠0且b≠0) 。
斜率:对于任意函数上任意一点,其斜率等于其切线与 x 轴正方向所成的角,即 k=tanα 。 ax+by+c=0中,k=-a/b 。
什么是极坐标系?
关于极坐标系,打开 百度百科 学习一下即可。
重点学到下面这个结论就行:
找资料的时候,发现一个解释的比较清楚的 博客 ,后续可以继续学习使用。
继续阅读资料,看到如下所示的图,这个图也出现在了很多解释原理的博客里面,但是图下面写了一句话
在这里直接蒙掉了,怎么就表示成极坐标系了?上面这个公式依旧是笛卡尔坐标系表示直线的方式呀,只是把 k 和 b 的值给替换掉了。
为何是这样的,具体原因可以参照下图。
centerchou 图/center
继续寻找关于霍夫变换的资料,找到一个新的概念 霍夫空间 。
在笛卡尔坐标系中,一条直线可以用公式 表示,其中 k 和 b 是参数,表示的是斜率和截距。
接下来将方程改写为 ,这时就建立了一个基于 k - b 的笛卡尔坐标系。
此时这个新的方程在 k - b 坐标系也有一个新的直线。
你可以在纸上画出这两个方程对应的线和点,如下图所示即可。
centerchou 图/center
新的 k - b 坐标系就叫做霍夫空间,这时得到一个结论,图像空间 x - y 中的点 对应了 霍夫空间 k - b 中的一条直线 ,即图像空间的点与霍夫空间的直线发生了对应关系。
如果在图像空间 x - y 中在增加一个点 ,那相应的该点在霍夫空间也会产生相同的点与线的对应关系,并且 A 点与 B 点产生的直线会在霍夫空间相交于一个点。而这个点的坐标值 就是直线 AB 的参数。
如果到这里你掌握了,这个性质就为我们解决直线检测提供了方法,只需要把图像空间的直线对应到霍夫空间的点,然后统计交点就可以达到目的,例如图像空间中有 3 条直线,那对应到霍夫空间就会有 3 个峰值点。
遍历图像空间中的所有点,将点转换到霍夫空间,形成大量直线,然后统计出直线交会的点,每个点的坐标都是图像空间直线方程参数,这时就能得到图像空间的直线了。
上述的内容没有问题,但是存在一种情况是,当直线趋近于垂直时,斜率 k 会趋近于无穷大,这时就没有办法转换了,解决办法是使用法线来表示直线。
上文提及的斜截式如下:
通过第二个公式,可以得到下述公式:
此时,我们可以带入一些数值进行转换。
图像空间有如下的几个点:
转换后的函数,都可以在霍夫空间 θ - ρ (横坐标是 θ ,纵坐标是 ρ )进行表示。
原理这时就比较清晰了:
除了一些数学知识以外,经典的博客我们也有必要记录一下,方便后面学习的时候,进行复盘。
本部分用于记录本文中提及的相关数学原理,后续还要逐步埋坑。
今天涉及了一点点数学知识,能力限制,大家一起学习,有错误的地方,可以在评论区指出,不胜感激。
希望今天的 1 个小时(今天内容有点多,不一定可以看完),你有所收获,我们下篇博客见~
相关阅读
技术专栏
逗趣程序员
对于气象绘图来讲,第一步是对数据的处理,通过各类公式,或者统计方法将原始数据处理为目标数据。
按照气象统计课程的内容,我给出了一些常用到的统计方法的对应函数:
在计算气候态,区域平均时均要使用到求均值函数,对应NCL中的dim_average函数,在python中通常使用np.mean()函数
numpy.mean(a, axis, dtype)
假设a为[time,lat,lon]的数据,那么
需要特别注意的是,气象数据中常有缺测,在NCL中,使用求均值函数会自动略过,而在python中,当任意一数与缺测(np.nan)计算的结果均为np.nan,比如求[1,2,3,4,np.nan]的平均值,结果为np.nan
因此,当数据存在缺测数据时,通常使用np.nanmean()函数,用法同上,此时[1,2,3,4,np.nan]的平均值为(1+2+3+4)/4 = 2.5
同样的,求某数组最大最小值时也有np.nanmax(), np.nanmin()函数来补充np.max(), np.min()的不足。
其他很多np的计算函数也可以通过在前边加‘nan’来使用。
另外,
也可以直接将a中缺失值全部填充为0。
np.std(a, axis, dtype)
用法同np.mean()
在NCL中有直接求数据标准化的函数dim_standardize()
其实也就是一行的事,根据需要指定维度即可。
皮尔逊相关系数:
相关可以说是气象科研中最常用的方法之一了,numpy函数中的np.corrcoef(x, y)就可以实现相关计算。但是在这里我推荐scipy.stats中的函数来计算相关系数:
这个函数缺点和有点都很明显,优点是可以直接返回相关系数R及其P值,这避免了我们进一步计算置信度。而缺点则是该函数只支持两个一维数组的计算,也就是说当我们需要计算一个场和一个序列的相关时,我们需要循环来实现。
其中a[time,lat,lon],b[time]
(NCL中为regcoef()函数)
同样推荐Scipy库中的stats.linregress(x,y)函数:
slop: 回归斜率
intercept:回归截距
r_value: 相关系数
p_value: P值
std_err: 估计标准误差
直接可以输出P值,同样省去了做置信度检验的过程,遗憾的是仍需同相关系数一样循环计算。