图论近似算法

结点的连通性

结点的连通性描述了结点直接连接的路劲数量。是对连通性的衡量指标

结点分组运算

给图中所有结点分组,按照连通和独立的条件给结点分组

聚类系数 clustering

图论中,集聚系数(也称群聚系数集群系数)是用来描述一个中的顶点之间结集成团的程度的系数。

图的匹配

对于一个给定的图G=(V,E),这幅图的一个匹配M是图G的一个子图(由原来的图的一部分顶点和一部分边构成的图),其中每两条边都不相邻(没有公共顶点)。在匹配图中,一个顶点连出的边数至多是一条。如果这个顶点连出一条边,就称这个顶点是已匹配的

顶点覆盖问题 Vertex cover

In graph theory, a vertex cover (sometimes node cover) of a graph is a set of vertices that includes at least one endpoint of every edge of the graph.

牛顿法求解n阶非线性方程

在数值计算的应用中,我们有时候会经常遇见需要求解非线性方程的情况。这里介绍牛顿求解法,来解决n次非线型方程求解的问题。

假设方程f(x) = a0 + a1*x + a2*x^2 + ... + an*x^n = 0 (其他非此形式的函数也可以泰勒展开成此形式,当然要能满足泰勒展开的条件。你不用担心做泰勒展开时会有太多系数,因为在迭代过程中,我们只需要用前2个系数即可)

牛顿求解法,也称为牛顿迭代法(有的地方只在牛顿求解法用在方程sqrt(x) - n = 0的求解时才称为牛顿迭代法)。 求解过程按如下步骤完成:

  1. 取一个初始值x0, 近似取f(x0 + d0) = f(x0) + d0*f'(x0)
  2. 令f(x0) + d0*f'(x0) = 0,求得d0 = -f(x0)/f'(x0)
  3. 再取x1 = x0 + d0, 并近似取f(x1 + d1) = f(x1) + d1*f'(x1)
  4. 不断重复2,3步骤,根据xn的值来求xn+1的值
  5. 直到|xn+1 - xn |足够小,或者满足其他指定条件时,停止迭代,并取xn+1为方程f(x)=0的解

初始值x0是比较影响求解答案的。一个不好的x0可能会导致迭代发散,而一个切当的x0可能会让问题很快收敛。在实际应用中如何取x0和如何判断收敛条件是算法应用的关键。

b样条曲线各参数数量的关系小结

参照关于b样条曲线的上一篇文献,这里可以给出一个关于knots, 阶数,poles数量的小总结。

  1. 结点数大于2倍的曲线阶数加1
  2. poles数量等于结点数减阶数再减1

如果,设结点组为[knot0,knot1,knot2, ..knotk],曲线阶数为j,极点数为p,则:

  1. k > 2*j (注意:k = 结点数 - 1)
  2. p = k - j (注意:k = 结点数 - 1)

另外,需要注意,在有的应用中,把重复的结点值算作multiplicity。例如它可能结点是(0,1), 但时multiplicity是(2,2),那么此处实际上指的结点组为(0,0,1,1)

从以上结论可以知道,每当我们给b样条曲线增加几个结点,就必须相应增加几个极点。

结点向量与B样条曲线形态的关系

1.定义

在文章开头对广义B样条多项式函数做一个简要介绍。广义的B样条多项式是由一个函数族定义的,这个函数族主要由3个基本要素确定:1,一组递增的结点向量;2,一个能有“只在局部起影响”的第0阶函数;3,一个从n到n+1阶递推的多项式变化公式。

1.1广义B样条多项式的3要素

首先是一组递增的结点向量[t0,t1,t2,...tk],这是一组大于0的单调递增的有理数,相邻两个数可以相等。这里用符号Nji(t)来定义第j阶的第i个关于t的多项式,那么所谓的第0阶函数族就可以用N00,N01,...N0n,N0n+1来表示。B样条多项式的模型设定为Nji(t)函数只在[ti,ti+j+1)区间起作用,那么对于N0i而言,将仅仅在t属于[ti,ti+1)区间时起作用,等于1.确切的说是函数N0i(t)在[ti,ti+1)区间取值为1,在这个区间以外都为0。接下来再定义一个从j-1阶到j阶的递推:Nji= (t-ti)*Nj-1i/(ti+j-ti)+ (ti+j+1 – t)Nj-1i+1/(tj+i+1– ti+1)。这就是由Cox和DeBoor提出的著名的B样条多项式递推公式。

在这个递推公式和N0i(t)函数的定义下,可以顺势给出函数N1i(t)的公式

N1i(t)= (t-ti)N0i(t)/(ti+1-ti)+ (ti+1-t)N0i+1(t)/(ti+2-ti+1)

在结点向量不存在相邻结点相等的情况下N1i(t)的函数图像如下:

图1

图1很清楚的描绘出N1i(t)函数族在结点向量中的作用方式。每个1阶的函数也都只在局部范围起作用,而且比0阶作用范围多了一段。每两个相邻的函数互相由一段作用范围发送重叠。除了首尾两端的结点段以外,其他的每个结点段都有2个函数起作用。这个性质可以推广到2阶,3阶到n阶。

1.2基于广义B样条多项式的B样条曲线

基于广义B样条的多项式创建B样条曲线的时候,为了使B样条曲线的开头与结尾的性质于由Riesenfeld多项式确定的均匀B样条曲线所保持一致性,做了一个这样的设定:由结点向量[t0,t1,t2,...tk]确定的广义B样条多项式的作用范围对于j阶曲线而言是从tj开始到tk-j结束。区间[tj,tk-j)被称为j阶B样条曲线的参数域。k和j的大小关系必须保证这个参数域的定义有效,也就是说k>2j。参数域的规定,保证了对于j阶B样条曲线,在每一小段[ti,ti+1)参数区间都有j+1个点起作用(否则在曲线的首尾两端部分将只有不足j+1个点起作用,将无法保证曲线的连续性)。

上面已经确定了多项式Nji(t)和参数域,下面给出B样条的参数点。结点向量[t0,t1,t2,...tk]对于j阶B样条曲线而言,将伴随着k-j个参数点的定义。因为对于j阶曲线而言,每个参数点在结点向量中的作用域都是j+1个小段。那么每个参数点的作用域为[t0,tj+1),[t1,tj+2),[t2,tj+3)... [tk-j-1,tk],一起正好k-j个作用域。设这些参数点为P0,P1,P2,...Pk-j-1,j阶的B样条曲线的定义公式为:

从这个公式还可以看出,因为k>2j所以imax>=j,参数域的规定使得这个公式中的i能够满足至少从0到j的完整性。值的一提的是,如果结点向量都是均匀分布(ti+1- ti等于常数)或者结点向量除了首尾j+1个结点相等外其他都是均匀分布,这个B样条曲线将退化成为均匀B样条曲线。如果结点向量数量刚刚等于2j+2(即k=2j+1>2j),而且前j+1个结点都等于0,后j+1个结点都等于1,那么这个B样条曲线将退化成Bezier曲线(中文译名有的叫贝塞尔曲线,有的叫贝济埃曲线,还有的叫贝兹曲线,总之它姓贝,所以也叫贝氏曲线)。最后,如果给每个控制点都再增加一个权重[p0,p1,...pk-j-1]的定义,让Nji(t)函数由多项式变成有理式:

这就是著名的非均匀有理B样条曲线的定义公式。

2结点向量的结点间隔对曲线的影响

在定义了B样条曲线后,我们再重新按照控制点Pi的作用范围来对1阶曲线的结点与多项式图做一个分析。

图2

图2清晰的给出了1阶曲线中每个点在结点区间中的作用范围。按照前文参数域的定义,1阶曲线的参数域从t1开始,因此,P0也是从t1开始起作用。每一段结点之间有两个控制点对其起作用。在1阶曲线中,结点间隔影响N1i(t)函数在每一段结点间的斜率,但这个斜率并不影响B样条曲线的形状。1阶B样条曲线是个多段线,无论t1到t2之间的间隔怎样变化,它都是从P0到P1之间的一条直线。但我们能隐约感觉得到,这个规律在更高的阶次曲线上就不成立了。事实上,单从B样条多项式来看,结点间隔所占一个多项式变化过程中的结点间隔的比例更像是对一个函数曲线在不同区间进行不同比例的“拉伸”和“压缩”。例如从图2中可以看出,N11(t)在[t1,t3]之间,[t1,t2]所占[t1,t3]的比例,描述了P1控制点在[t1,t2]区间的控制能力所占P1对整个控制领域的控制能力的权重。例如,t2-t1> t3-t2的意思是,P1把更多的控制能力放在了与P0的协同控制上而把相对更少的控制能力放在了与P2的协同上。

在1阶B样条曲线上,我们观察不到这种“控制能力”对曲线的影响,下面我们将观察它对2阶B样条曲线的影响。

我们考察两个2阶B样条曲线,它们的结点向量分别是[0,0, 0, 0.5, 1, 1, 1]和[0,0, 0, 0.9, 1, 1, 1]。这两条B样条曲线有着共同的控制点P0(0,0,0);P1(10,0,0);P2(10,10,0); P3(0,10,0)。接下来首先分别计算这两个结点向量所产生的B样条多项式。

首先,对于前者(下面的计算过程涉及到了对于重复点的一个约定,即在B样条多项式间,对重复点区间的单项做省略处理。):

对于后者:

分别做出两种情况的B样条曲线的控制点影响坐标图如下:

对于前者:

图3

对于后者:

图4

对比图3和图4中的[0,0.5)和[0,0.9)区间。P0,P1,P2的作用函数曲线变得更缓更长,对于[0,1]区间,P0,P1,P2在后一种情况中,对整个曲线的驾驭权重增加了,这直接影响了曲线的形状。前一种情况在t=0.5时,和后一种情况的t=0.9时P0,P1,P2的重心系数明显发生了变化(P0没有变化,P1减弱了,P2加强了)。

如果一个j阶B样条曲线(j>1)的某个结点区间[tm,tm+1)相对于[tm-j,tm+j]长度增加(tm在参数域中),这个区间对应起作用的控制点为Pv,Pv+1,...Pv+j这个区间的以下规律可以按曲线的阶次推广:

1.Pv,Pv+1,...Pv+j所控制的曲线段长度相对增加。

2.Pv,Pv+1,...Pv+j所控制的曲线段弯曲变化率降低。

3.对于Pv+h点,如果h<= j/2, 如果[tm,tm+1)的长度相对于[tm-j+h,tm+h)区间增加则点Pv+h的控制权重则减少,反之则增加。如果h> j/2, 如果[tm,tm+1)的长度相对于[tm-j+h,tm+h)区间增加则点Pv+h的控制权重则增加,反之则减少。

如果结点区间[tm,tm+1)相对于长度减少,则上诉1,2结论也取反。例如,图3到图4的变化,会导致B样条曲线向P2,及其相邻的点靠拢,顺序上离P2越远的点所发挥的控制力量都减弱。

下面给出实际绘制出来的B样条曲线,可以验证上述结论的正确性:

这两个图即是上文所讨论的结点发生变化的两个2阶B样条曲线,很明显右图中P2的控制权重增加了。

再看看多一个控制点的情况:


对比这两幅图中的样条曲线

1.结点在[t3,t4)区间长度相对减少,根据结论3,这会导致P0,P1权重增加,P2权重减少。

2.结点在[t4,t5)区间长度相对增加,根据结论3,这会导致P2权重减少,P1,P3权重增加。

3.结点在[t5,t6)区间长度相对增加,根据结论3,这会导致P3权重增加,P2,P4权重减少。

在nurbs曲线中,相比于控制点的权重值,结点之间的间隔大小影响的是一段控制点对曲线操控能力的相对权重。例如,在一段曲线中,我们希望P1,P2,P3这3个控制点相比于其他部分更起主导作用,而且控制得更细致一些,则应该在这3个点所对应的结点部分加长一些。值的注意的是,这样做,既增加了这段曲线的控制能力而又没有增加控制点和曲线的阶数,唯一损失的是另外一部分相对减少了的结点长度。

对拓扑空间的一种理解

拓扑空间是为了定向研究一类数学问题而给出来的一个抽象而且广义的一个数学定义。其定义的公理主要是想构造一个抽象的逻辑能够更好的描述连续的变化,或者说连续的映射。为了摆脱“连续”的概念对“距离”概念的依赖, 拓扑空间的一般定义,是通过集合的表述来表达的。

定义研究对象

描述拓扑空间时,首先需要确定的是研究对象。一开始我们需要有一个整体上研究对象的集合。例如,如果是研究3维空间的问题,那么定义拓扑空间时,通常就是对整个3维空间的点集所定义的。无论整么样,在定义一个拓扑空间时,我们需要有一个研究对象的总集合,这个集合中存在至少一个单位元素。我们假设这个集合为X。

定义邻域

对于集合X我们定义一种规则。对于X集合内的每一个单位元素x,按照这个规则能够从X集合中找出该单位元素所有的一组子集。而且这组子集都是这个点的邻域。

只要满足以下4个规则的子集,都叫做x的邻域:

  1. x在自己的每一个邻域里。
  2. x的任何两个领域的交集也在这一组邻域里。
  3. 如果N是x的邻域,U为X的子集,而N包含于U,则U也在x的这一组邻域里。
  4. 如果N是x的邻域,Ni表示一个这样的集合:{z属于N|N是z的邻域},则Ni也是x的邻域。(Ni也叫做N的内部)(这一条定义也可以理解为,如果Ni是N中排除所有边界点后剩下的点的集合,则Ni也是x的邻域。然后我们通过如下方式来定义N的边界点:该点虽然属于N,但N不是该点的邻域)

于是,拓扑空间就是,有研究对象的集合,而且给这个组研究对象中每一个单位元素定义了邻域。这样的一个结构,就叫做拓扑空间。

python中用4元数实现3D旋转

利用4元数实现3d旋转有非常多的优点,相比于4x4矩阵,欧拉角的方式具有更好的计算高效性。一般的,4元数(x,y,z,w)表达空间旋转的表达方式如下:

设空间单位向量(a,b,c), 和一个夹角theta。则表示绕空间向量(a,b,c)沿顺时针(沿该向量相同方向看)旋转theta角的四元数为q = (a*sin(theta/2), b*sin(theta/2), c*sin(theta/2), cos(theta/2))。

若,空间上某向量v,按照四元数q旋转后得到的值为u,则

u = qvq(-1) (这里的乘法方式为哈密尔顿积

其中q(-1)指的是四元数q的逆,对于单位4元数的逆,恰恰等于原四元数的虚数部分去反。即,

若q = xi + yj + zk + w, 则q(-1) = -xi -yj -zk + w

除此以外,空间的多个旋转还可以串连起来,例如:

q = q2q1 计算所的的q为先按照q1旋转,再按照q2旋转的4元数。

实际工程应用中,可以使用numpy的numpy-quaternion和scipy的scipy.spatial.transform.Rotation来配合实现4元数的运算,以及4元数与3x3或4x4矩阵的转化。

示例代码如下:(如下代码要求scipy >= 1.4.0)

from scipy.spatial.transform import Rotation as R
import numpy as np
def makeRotationMatrix(theta, axis):
     m4 = np.identity(4)
     s = np.sin(theta/2)
     c = np.cos(theta/2)
     quat = np.zeros([4])
     quat[0:3] = axis * s
     quat[3] = c
     r = R.from_quat(quat)
     m3 = r.as_matrix()
     m4[0:3, 0:3] = m3
     return m4

def makeTranslationMatrix(vec):
     m4 = np.identity(4)
     m4[0:3,3] = vec
     return m4