精华内容
下载资源
问答
  • vb绘制极坐标曲线图形
    2021-07-06 08:43:22

    第十一章 真实感图形技术

    1,简单光照明模型

    2,多边形绘制方法

    3,透明

    4,整体观照明模型

    5,光线跟踪算法

    第十章 真实感图形绘制

    光照模型 (Illumination Model):计算某

    一点的光强度的模型

    11.1 真实感图形的 特点

    ? 能反映物体表面颜色和亮度的细微变化

    ? 能表现物体表面的质感

    ? 能通过光照下的物体阴影,极大地改善场景的

    深度感和层次感,充分体现物体间的相互遮挡

    关系

    ? 能模拟透明物体的透明效果和镜面物体的镜面

    效果

    影响观察者看到的表面颜色的因素

    ① 物体的几何形状

    ②光源 位置、距离、颜色、数量、强度、种类

    ③环境 遮挡关系、光的反射与折射, 阴影

    ④视点位置

    ⑤物性 材料、颜色、透明度 折射性

    ⑥表面光洁度

    光源

    ① 几何性质

    – 点光源

    – 线光源

    – 面光源

    ② 光谱组成

    – 白色光等能量的各种波长可见光的组合

    – 彩色光

    – 单色光

    11.2 真实感图形学早期发展

    ? 1967年,Wylie等人第一次在显示物体时

    加进光照效果,认为光强与距离成反比。

    ? 1970年,Bouknight提出第一个光反射模

    型,Lambert漫反射+环境光

    ? 1971年,Gouraud提出漫反射模型加插值

    的思想

    ? 1975年,Phong提出图形学中第一个有影

    响的光照明模型

    相关物理知识

    ? 光的传播

    – 反射定律:入射角等于反射角,而且反射光

    线、入射光线与法向量在同一平面上

    光源 法向量

    入射光 反射光

    视线

    折射定律

    – 折射定律:折射线在入射线与法线构成的平

    面上,折射角与入射角满足

    入射光

    折射光

    ?

    ?

    2

    ?

    1

    ?

    1

    2

    s in

    s in

    ??

    ???

    能量关系

    – 在光的反射和折射现象中的能量分布:

    – 下标为 i,d,s,t,v的能量项分别表示 为入射光强,

    漫反射光强,镜面反射光强,透射光强,吸

    收光强

    – 能量是守恒的

    i d s t vI I I I I? ? ? ?

    11.3 简单光照明模型

    模拟物体表面的光照明物理现象的数学

    模型-光照明模型

    简单光照明模型 亦称局部光照明模型,

    其假定物体是不透明的,只考虑光源的

    直接照射,而将光在物体之间的传播效

    果笼统地模拟为环境光。

    可以处理物体之间光照的相互作用的模

    型称为整体光照明模型

    简单光照明模型

    光照射到物体表面,主要发生:

    反射

    透射(对透明物体)

    部分被吸收成热能

    反射光,透射光决定了物体所呈现的颜

    简单光照明模型 -环境光

    假定物体是不透明的(即无透射光)

    ? 环境光,在空间中近似均匀分布,即在任何位置、

    任何方向上强度一样,记为 Ia

    ? 环境光反射系数 Ka,在分布均匀的环境光照射

    下,不同物体表面所呈现的亮度未必相同,因为它们

    的环境光反射系数不同。

    ? 光照明方程(仅含环境光),Ie = KaIa

    Ie为物体表面所呈现的亮度。

    简单光照明模型 - 环境光例子

    ? 具有不同环境光反射系数的两个球

    0.1?aI

    4.0?aK 8.0?

    aK

    简单光照明模型 -环境光

    ? 缺点:虽然不同的物体具有不同的亮度,

    但是同一物体的表面的亮度是一个恒定

    的值,没有明暗的自然过度。

    简单光照明模型

    ? 考虑引入点光源。

    ? 点光源:几何形状为一个点,位于空间中的某

    个位置,向周围所有的方向上辐射等强度的光。

    记其亮度为 Ip

    ? 点光源的照射,在物体的不同部分其亮度也不

    同,亮度的大小依赖于物体的朝向及它与点光源之间

    的距离,

    简单光照明模型,-漫反射 角度

    余弦的推导

    ? 漫反射

    – 粗糙、无光泽物体(如粉笔)表面对光的反射

    – 光照明方程

    ? 漫反射的亮度

    ? 点光源的亮度

    ? 漫反射系数

    ? 入射角

    漫反射光的强度

    只与入射角有关

    ]2,0[c o s ??? ?? dpd KII

    pI

    dK

    ?

    dI

    简单光照明模型 -漫反射

    ? 将环境光与漫反射结合起来

    一般取 Ia= (0.02~0.2)Id

    ? 例子

    )( NLKIKIIII dpaade ?????

    简单光照明模型 -漫反射

    缺点:对于许多物体,使用上式计算其反

    射光是可行的,但对于大多数的物体,

    如擦亮的金属、光滑的塑料等是不适用

    的,原因是这些物体还会产生镜面发射。

    简单光照明模型 -镜面反射

    ? 镜面反射

    – 光滑物体(如金属或塑料)表面对光的反射

    ? 高光

    – 入射光在光滑物体表面形成的特别亮的区域

    简单光照明模型 -镜面反射

    ? 理想镜面反射

    ? 观察者只能在反射方向上才能看到反射

    光,偏离了该方向则看不到任何光。

    简单光照明模型 -镜面反射

    ? 非理想镜面反射

    ? P为物体表面上一点,L为从 P指向光源的单位

    矢量,N为单位法矢量,R为反射单位矢量,V

    为从 P指向视点的单位矢量

    光滑平面

    I = Ip K scosna

    镜面

    简单光照明模型 -镜面反射

    ? 镜面反射

    ? Is为 镜面反射光强。 点光源的亮度

    – Ks是与物体有关的镜面反射系数。 n为 镜面反射指数,n越

    大,则 Is随 a的增大衰减的越快。

    – n的取值与表面粗糙程度有关。

    – n越大,表面越平滑(散射现象少,稍一偏离,明

    暗亮度急剧下降)

    – n越小,表面越毛糙(散射现象严重)

    ansps KII c o s? nsps RVKII )( ??

    pI

    简单光照明模型 -镜面反射

    – 反射方向计算

    – L在 N上的投影矢量为 Ncosu,则 S+L= Ncosu

    记矢量 S= Ncosu -L

    则有 R= Ncosu +S

    ? ? LLNNLNR ????? 2c o s2 ?

    N

    R

    L V

    ??

    S S

    简单光照明模型 -Phong光照明模型

    ? 简单光照明模型模拟物体表面对光的反

    射作用,光源为点光源

    ? 反射作用分为

    – 物体间作用用环境光 (Ambient Light)

    – 漫反射 (Diffuse Reflection)

    – 镜面反射 (Specular Reflection)

    简单光照明模型 -Phong光照明模型

    ? Phong光照明模型的综合表述:由物体表

    面上一点 P反射到视点的光强 I为环境光

    的反射光强 Ie、理想漫反射光强 Id、和镜

    面反射光 Is的总和。

    ])()([ nsdpaa

    sde

    RVKNLKIKI

    IIII

    ?????

    ???

    简单光照明模型 -Phong光照明模

    型的实现

    ? 对物体表面上的每个点 P,均需计算光线

    的反射方向。为了减少计算量,假设:

    – 光源在无穷远处,L为常向量

    – 视点在无穷远处,V为常向量

    – ( H?N)近似( R?V),H为 L与 V的平分向量

    N HL R

    a

    Vb

    H----L和 V的角平分线

    ? 对所有的点总共

    只需计算一次 H的

    值,节省了计算

    时间

    简单光照明模型 -Phong光照明模型

    ? Phong模型几何

    P

    L N H

    R

    V

    简单光照明模型 -光的衰减

    ? 光的衰减

    两个阶段:

    1)从光源到物体表面的过程中的衰减

    2)从物体表面到人眼过程中的衰减

    总的效果:物体表面的亮度降低

    ? 光照明方程

    1)有效衰减函数的加入

    2)深度暗示技术的加入

    简单光照明模型 -光的衰减

    ? 光的衰减

    – 光在光源到物体表面过程中的衰减

    – 光强按 1/d2 进行衰减:

    缺点:当 d很大时,变化很小;当 d很小时,变

    化很大。

    ? 衰减函数

    ? 光照明方程

    )1,1m i n ()( 2

    210 dcdcc

    df ???

    ])()([)( nsdpaa RVKNLKIdfKII ?????

    简单光照明模型 -光的衰减

    – 光在物体表面到人眼过程中的衰减

    ? 深度暗示( Depth Cueing)技术:最初用于线框

    图形的显示,使距离远的点比近的点暗一些。经

    过改进,这种技术同样适用于真实感图形显示。

    ? 设前参考面 Z=Zf,后参考面 Z=Zb;其比例因子

    分别为 Sf和 Sb( Sf和 Sb e[0,1])。给定物体上一

    点的深度值 Z0,该点对

    应的比例因子 S0按如下

    方式确定

    前参考面

    后参考面

    简单光照明模型 -光的衰减

    ? 当 Z0 >Zf时,取 S0=Sf

    ? 当 Z0

    ? 当 Z0e [Zb,Zf]时,取

    ? 原亮度 I按比例 S0与融和亮度 Idc混合,目

    的是获得最终用于显示的亮度 I’, Idc由

    用户指定,

    )( 00 b

    bf

    bf

    b ZZZZ

    SSSS ?

    ?

    ???

    dcISISI )1( 00 ????

    简单光照明模型 -光的衰减

    ? 特例:

    ? 取 Sf=1,Sb=0,Idc=0,则当物体位于参考面

    之前时,S0= Sf= 1,I’ =I,即亮度没有

    被衰减。当物体位于后参考面之后时,

    S0 = Sb =0,I’ =Idc=0,即亮度衰减为 0。

    而当 Z0e [Zb,Zf]时,I’ =S0I,亮度被部分

    衰减。由此可以产生较好的效果。

    dcISISI )1( 00 ????

    简单光照明模型 -彩色场景的产生

    ? 产生彩色

    – 选择合适的颜色模型 ----RGB模型

    – 为颜色模型中的每一种基色建立光照明方程

    ?

    ?

    ?

    ?

    ?

    ?????

    ?????

    ?????

    ])()([)(

    ])()([)(

    ])()([)(

    n

    sBdBpBaBaBB

    n

    sGdGpGaGaGG

    n

    sRdRpRaRaRR

    RVKNLKIdfKII

    RVKNLKIdfKII

    RVKNLKIdfKII

    简单光照明模型 -彩色场景的产生

    – 系数分解

    ? 上述各等式中,右端的矢量用来控制表面的基本

    颜色,当选定了物体表面的颜色之后,它们就固

    定不变了。用户通过调节 Ka,Kd,Ks来改变表

    面的反射率。

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    dB

    dG

    dR

    a

    aB

    aG

    aR

    C

    C

    C

    K

    K

    K

    K

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    dB

    dG

    dR

    d

    dB

    dG

    dR

    C

    C

    C

    K

    K

    K

    K

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    sB

    sG

    sR

    s

    sB

    sG

    sR

    C

    C

    C

    K

    K

    K

    K

    简单光照明模型 -彩色场景的产生

    – 新的光照明方程

    – 统一表示

    ?

    ?

    ?

    ?

    ?

    ?????

    ?????

    ?????

    ])()([)(

    ])()([)(

    ])()([)(

    n

    sBsdBdpBaBdBaB

    n

    sGsdGdpGaGdGaG

    n

    sRsdRdpRaRdRaR

    RVCKNLCKIdfICKI

    RVCKNLCKIdfICKI

    RVCKNLCKIdfICKI

    ])()([)( nssddpada RVCKNLCKIdfICKI ????? ??????

    BGR,,??

    简单光照明模型 -多个光源

    ? 采用多个光源

    – 采用 m个光源的光照明方程

    ?

    ?

    ???

    ??

    m

    i

    n

    issiddpi

    ada

    RVCKNLCKIdf

    ICKI

    i

    1

    ])()([)( ???

    ???

    简单光照明模型 -多个光源

    – 例子:其中 a图:线框图 b图:环境光

    c图:增加漫反射 d图:增加镜面反射

    e图:增加光的衰减 f图:两个点光源

    Phong光照明模型的 不足

    ? Phong光照明模型是真实感图形学中提出

    的第一个有影响的光照明模型

    ? 经验模型,Phong模型存在不足:

    – 显示出的物体象塑料,无质感变化

    – 没有考虑物体间相互反射光

    – 镜面反射颜色与材质无关

    – 镜面反射大入射角失真现象

    11.4 多边形绘制方法

    ? 分类,均匀着色与光滑着色

    ? 均匀着色

    方法:任取多边形上一点,利用光照明方

    程计算出它的颜色,用这个颜色填充整

    个多边形

    适用场合,1)光源在无穷远处;

    2)视点在无穷远处;

    3)多边形是物体表面的精确表示;

    多边形绘制方法

    ? 缺点:产生的图形效果不好。

    ? 如左图:相邻两个多边形的法向

    不同,计算出来的颜色也不同,

    因此造成整个物体表面的颜色过

    渡不光滑。

    ? 如何解决?

    ? 光滑着色,亦称插值着色

    Gouraud着色方法

    Phong着色方法

    Gouraud着色方法

    ? Gouraud于 1971年提出,又被称 Gouraud

    明暗处理

    ? 基本思想:在每个多边形顶点处计算颜

    色,然后在各个多边形内部进行线性插

    值,得到多边形内部各点颜色。即它是

    一种 颜色插值着色方法。

    ? 注意,Gouraud着色方法并不是孤立的处

    理单个多边形,而是将构成一个物体表

    面的所有多边形(多边形网格)作为一

    个整体来处理。

    Gourand 着色方法

    ? 对多边形网格中的每一个多边形,

    Gourand 着色处理分为如下四个步骤:

    – 步骤

    1、计算多边形的单位法矢量

    2、计算多边形顶点的单位法矢量

    ? 与某个顶点相邻的所有多边形的法向平

    均值近似作为该顶点的近似法向量

    ? 计算出的平均法向一般与该多边形物体

    近似曲面的切平面比较接近

    Gouraud着色方法 -顶点法向计算

    ?

    ?

    ?

    ??

    n

    i

    n

    i

    i

    v

    Ni

    N

    N

    1

    1

    Gourand 着色方法

    3、利用光照明方程计算顶点光强(颜色)

    4、对多边形顶点光强(颜色)进行双线性

    插值,获得多边形内部各点的光强(颜色)

    Gourand 着色方法 -光强插值

    ? 双线性光强插值:假设待绘制的三角形投影为

    P1P2P3,Pi的坐标为 (xi,yi),i=1,2,3;一条扫描

    线与三角形的两条边分别交于 A(xA,yA),B(xB,yB)

    两点。 P(x,y)是 AB上的一点。 A点的颜色 IA由 P1、

    P2点的颜色 I1,I2线性插值得到

    B

    AB

    A

    A

    AB

    B

    P

    BB

    B

    AA

    A

    I

    xx

    xx

    I

    xx

    xx

    I

    I

    yy

    yy

    I

    yy

    yy

    I

    I

    yy

    yy

    I

    yy

    yy

    I

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    3

    31

    1

    1

    31

    3

    2

    21

    1

    1

    21

    2

    Gourand 着色方法 -增量算法

    ? 采用增量算法可以加速计算。

    ? 1)当扫描线 y递增一个单位变为 y+1时,

    IA,IB的增量分别为 DIA, DIB,即

    31

    31

    21

    21

    ,1,,1,

    yy

    II

    I

    yy

    II

    I

    IIIIII

    BA

    ByByBAyAyA

    ?

    ?

    ?D

    ?

    ?

    ?D

    D??D?? ??

    其中:

    Gourand 着色方法 -增量算法

    ? 2)当 x递增一个单位时,IP的增量为 DIP即

    AB

    AB

    P

    PxPxP

    xx

    II

    I

    III

    ?

    ?

    ?D

    D??

    ?

    其中:

    ,1,

    Gourand 着色方法

    ? 优点,能有效的显示漫反射曲面,计算量小

    ? 缺点:

    ? 1、高光有时会异常

    ? 2、当对曲面采用不同的多边形进行分割时会产生不同

    的效果。

    ? 3,Gouraud明暗处理会造成表面上出现过亮或过暗的

    条纹,称为马赫带( Mach_band)效应

    ? 改进- Phong提出双线性法向插值,以时间为代价,解

    决高光问题

    Phong着色方法

    ? 基本思想:通过对多边形顶点的法矢量进

    行插值,获得其内部各点的法矢量,又称

    为法向插值着色方法。

    – 步骤

    1、计算多边形单位法矢量

    2、计算多边形顶点单位法矢量

    3、对多边形顶点法矢量进行双线性插值,获

    得内部各点的法矢量

    4、利用光照明方程计算多边形内部各点颜色

    Phong着色方法 -法向插值

    NA由 N1,N2线性插值得到:

    B

    AB

    A

    A

    AB

    B

    P

    BB

    B

    AA

    A

    N

    xx

    xx

    N

    xx

    xx

    N

    N

    yy

    yy

    N

    yy

    yy

    N

    N

    yy

    yy

    N

    yy

    yy

    N

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    ?

    3

    31

    1

    1

    31

    3

    2

    21

    1

    1

    21

    2

    Phong着色方法 -增量算法

    ? 采用增量算法可以加速计算。

    ? 1)当扫描线 y递增一个单位变为 y+1时,

    NA,NB的增量分别为 DNA, DNB,即

    31

    31

    21

    21

    ,1,,1,

    yy

    NN

    N

    yy

    NN

    N

    NNNNNN

    BA

    ByByBAyAyA

    ?

    ?

    ?D

    ?

    ?

    ?D

    D??D?? ??

    Phong着色方法 -增量算法

    ? 2)当 x递增一个单位时,IP的增量为 DIP即

    AB

    AB

    P

    PxPxP

    xx

    NN

    N

    NNN

    ?

    ?

    ?D

    D??

    ?

    其中:

    ,1,

    Phong着色方法

    优点:

    Phong着色方法绘制的图形比 Gouraud方法更

    真实,体现在两个方面:高光区域的扩散,产

    生正确的高光区域

    缺点:

    1,Phong着色方法计算量远大于 Gouraud着

    色方法

    2、在处理某些多边形分割的曲面时,Phong

    算法还不如 Gouraud算法好。

    增量式模型示例

    ? 牛的三角网格模型

    ? 用简单光照明模型显示

    ? 用增量式光照明模型显示

    插值多边形绘制方法

    ? 着色方法存在的问题

    – 不光滑的物体轮廓:物体边缘轮廓是折线段

    而非光滑曲线

    插值着色多边形绘制方法

    – 透视变形

    – 方向依赖性

    插值着色多边形绘制方法

    – 公共顶点处颜色不连续

    – 顶点方向不具有代表性

    11.5 透明

    ? 现实世界中有许多透明物体,如玻璃等。

    透过透明物体,可以观察到其后面的景

    物。如何模拟这种透明效果呢?

    ? 模拟透明的最简单的方法是忽略光线在

    穿过透明体时所发生的折射。虽然这种

    模拟方法产生的结果不真实,但在许多

    场合往往非常有用。例如:我们有时希

    望能够看到透过某透明物体观察其后面

    的景物,而又不希望景物应为折射而发

    生变形。

    透明效果的简单模拟

    ? 不考虑透明体对光的折射以及透明物体

    本身的厚度

    ? 光通过物体表面不改变方向

    ? 产生简单透明效果的方法

    插值透明方法

    过滤透明方法

    ),( yx

    a

    I

    b

    I

    透明体

    不透明体

    简单透明 -插值透明

    ? 假设:多边形 1是透明的,它

    位于观察者与不透明的多边形

    2之间。像素的颜色 I?由 A,B

    两点的颜色 I?1和 I?2插值产生,

    ? 其中 Kt1是多边形 1的透射系数。

    ? Kt1范围( 0,1)

    ? Kt1 =0表示多边形完全不透明,

    所以 I? = I?1

    ? Kt1 =1表示多边形完全透明,

    所以 I? = I?2

    21 11 )1( ??? IKIKI tt ???

    简单透明 -插值透明

    ? 为了产生逼真的效果,

    通常只对两个多边形表

    面颜色的环境光分量和

    漫反射分量采用

    ?

    ? 进行计算,得到的结果

    再加上多边形 1的镜面反

    射分量作为像素的颜色

    值。

    21 11 )1( ??? IKIKI tt ???

    简单透明 -过滤透明

    ? 过滤透明方法将透明物体看

    作一个过滤器,有选择的允

    许某些光透过而屏蔽了其余

    的光。对右图有:

    ? 其中 Kt1仍是多边形 1的透射系

    数,但不再局限于 (0~1)。 Kt1

    越大,多边形 2的颜色透过来

    的越多。 Ct?对不同的颜色各

    不相同。 Ct? =0表示某种颜色

    的光不能透过多边形 1。

    21 1 ???? ICKII tt??

    简单透明

    ? 无论采用插值透明方法还是采用过滤透

    明方法,当多边形 1之前还有其它的透明

    多边形时,I?都要递归计算。

    ? 简单透明比较容易结合到多边形绘制算

    法中。

    考虑折射的透明

    ? 折射定律

    ? 其中,?i, ?t分别是入射

    光线在空气,物体中的

    折射率,?i, ?t分别是入

    射角和折射角

    i

    t

    t

    i

    ?

    ?

    ?

    ? ?

    s in

    s in

    考虑折射的透明 -透射矢量的计算

    – 设单位入射光矢量为 I

    (方向与光线的入射

    方向相反),单位法

    矢量为 N,单位透射

    光矢量为 T,则

    INT

    INM

    M

    IN

    t

    NMT

    ti

    i

    t

    i

    i

    ti

    i

    t

    t

    ????

    ??

    ?

    ?

    ?

    ?

    ??

    ?

    ?

    ?

    ?

    ?

    ???

    ????

    ?

    ?

    ??

    ??

    )c o sc o s(

    c o s

    c o s

    s i n

    i

    s i n

    c o s

    考虑折射的透明

    ? 当光线从高密度介质向低密度介质时,

    ?i>?t,即 ?t>?i。如果入射角不断增大,到

    一定的程度,折射角 ?t?90度,此时透射

    光线沿着平行于分界面的方向传播,称

    此时的 ?i为临界角度,记为 ?c 。当 ?i > ?c

    时,发生全反射,透射与反射光合二为

    一。

    ? 如何产生带有折射的透明效果呢?

    光透射模型的研究

    ? 早期简单透射现象的模拟

    ? 1980年,Whitted光透射模型,首次考虑

    了光线的折射现象

    ? 1983年,在 Whitted的基础上,Hall光透

    射模型,考虑了漫透射和规则透射光

    11.6 整体光照明模型

    简单光照模型(亦称局部光照模型)不考

    虑周围环境对当前景物表面的光照明影

    响,忽略了光在环境景物之间的传递,

    很难表现自然界复杂场景的高质量真实

    感图形。为了增加图形的真实感,必须

    考虑环境的漫射、镜面反射和规则投射

    对景物表面产生的整体照明效果。

    整体光照明模型

    ? 物体表面入射光的构成

    ( 1)光源直接照射

    ( 2)其它物体的反射光

    ( 3)透射光

    ? 局部光照明模型仅考虑了( 1)

    整体光照明模型

    ? 例如:从视点观察到的物体 A表

    面的亮度来源于三方面的贡献:

    ( 1)光源直接照射到 A的表面,然

    后被反射到人眼中的光产生的。

    ( 2)光源或其它物体的光经 A物体

    折射到人眼中的光产生的。

    ( 3)物体 B的表面将光反射到物体

    A的表面,再经物体 A的表面反

    射到人眼中产生的。

    ? 局部光照明模型仅考虑了( 1)

    Witted光照模型

    ? Whitted光照模型基于如下假设:

    ? 物体表面向视点方向 V辐射的光亮度 I?由三

    部分组成:

    ( 1)光源直接照射引起的反射光亮度 Il?。

    ( 2)来自 V的镜面反射方向 R的其它物体反

    射或折射来的光的亮度 Is?。

    ( 3)来自 V的透射方向 T的其它物体反射或

    折射来的光的亮度 It?

    Witted光照模型

    ? Witted光照模型,I?= Il? + Ks Is? + KtIt?

    ? 或

    – Is?为镜面反射方向的入射光强度; Ks为镜面

    反射系数,为 0~ 1之间的一个常数

    – It?为折射方向光强,Kt为透射系数,是 0 ~

    1之间的常数

    – Il?的计算可采用 Phong模型

    因此,关键是 Is和 It的计算。如何计算呢?

    ?????? tttsssl ICKICKII ???

    Witted光照模型 -反射、折射方向计算

    ? 已知视线方向 V,求其反射方向 R与折射

    方向 T( N是表面的法向方向 )

    ? 视线 V的反射方向 R

    ? 折射方向 T

    N

    V

    R

    L

    T

    i

    ?

    i

    ?

    t

    ?

    i

    ?

    t

    ?

    VVNNR ??? )(2

    INT ti

    t

    i

    ????

    ?

    ?

    ?

    ???

    ?

    )c o sc o s(

    光线跟踪算法的基本原理

    ? 自然界中光线的传播过程

    光源 物体表面 物体表面

    人眼

    ? 光线跟踪过程 ----光线传播的逆过程

    光线跟踪算法的基本原理

    ? 从视点向每个 象素发出一

    条光线,它与场景中的一

    些物体表面相交,最近的

    交点即为可见点,记为 P,

    像素的亮度即由 P点的亮

    度确定。由 Whitted光照

    模型可知,P点的亮度由

    三部分组成:其中 Il?可以

    直接由局部光照模型计算

    得到。

    光线跟踪算法的基本原理

    ? 为了求 Is?和 It?,从 P点发出反射光线和透射光线,它分别交场

    景中的物体表面于 Ps和 Pt,Ps和 Pt点的亮度即分别为 Is?和 It?,

    将它们求出代入 Whitted模型即可。但是,Is?和 It?同样由

    Whitted模型确定,即 Whitted模型是一个递归式,从而计算 Is?

    和 It?需要重复以上的计算过程,计算局部光亮度、发出反射光

    线与透射光线。。。可以用一棵光线树来表示

    光线跟踪算法的基本原理

    ? 递归终止条件:

    1、光线不与场景中的任何物体相交

    2、被跟踪的光线达到了给定的层次

    3、由于 Ks和 Kt都小于 0,当光线经过反射

    和折射后,其亮度会衰减。因此可以预

    先设置一个阈值,在进行光线跟踪时,

    若被跟踪光线对像素亮度的贡献小于这

    个阈值,便停止跟踪。

    光线跟踪算法 -算法描述

    设置视点,投影平面以及窗口的参数;

    For (窗口内的每一条扫描线)

    for (扫描线上的每一个像素)

    { 确定从视点指向像素中心的光线 ray;

    像素的颜色 =RayTracing(ray,1);

    }

    光线跟踪算法描述Color RayTracing(Ray ray,int depth)

    { 求 ray与物体表面最近的交点 P;

    if (有交点)

    { 用局部光照明模型计算 P点的 Ic;

    color = Ic;

    if ( depth

    { 计算 ray的反射光线;

    Is=RayTracing(反射光线,depth+1);

    if (物体是透明的)

    { 计算 ray的透射光线;

    It=RayTracing(透射光线,depth+1);

    }

    color = Ic + Is + It ;

    }

    } else color = black;

    return color ;

    }

    光线跟踪算法

    根据光线跟踪基本原理,假设在显示分辨率为

    N× M的屏幕上生成图形,必须从视点出发通过

    屏幕向景物发射 N× M条光线。设每根光线在场

    景中经过反射和折射平均派生出 d根光线,并

    设每一景物交点朝光源发射 m 条阴影探测光线,

    则总的光线数增加至 (m+1)dNM。当 m为 2,d为 5,

    N为 800,M为 600时,其光线数目是 720万条。

    这意味着需进行 720万次直线与景物的求交计

    算才能完成图形的绘制。庞大的求交量使图形

    绘制耗费大大增加。生成一幅中等复杂程度的

    真实感图形需要数分钟至数小时。

    光线跟踪算法

    优点:

    能够方便的产生阴影,模拟镜面反射与

    折射现象。

    缺点:

    计算量大,每一条光线都要与场景中的

    物体进行求交、计算光照模型等。

    更多相关内容
  • 极坐标函数 ρ=exp⁡(cos⁡(θ))−2cos⁡(4θ)+(sin⁡(θ/12))5 \rho=\exp(\cos(\theta))-2\cos(4\theta)+(\sin(\theta/12))^5 ρ=exp(cos(θ))−2cos(4θ)+(sin(θ/12))5 的曲线是何样的?估计用MATLAB作图最...

    先看MATLAB如何做

    极坐标函数
    ρ = exp ⁡ ( cos ⁡ ( θ ) ) − 2 cos ⁡ ( 4 θ ) + ( sin ⁡ ( θ / 12 ) ) 5 \rho=\exp(\cos(\theta))-2\cos(4\theta)+(\sin(\theta/12))^5 ρ=exp(cos(θ))2cos(4θ)+(sin(θ/12))5
    的曲线是何样的?估计用MATLAB作图最简单,三句话即可:

    theta=0:0.01:20*pi;
    rho=exp(cos(theta))-2*cos(4*theta)+(sin(theta/12)).^5;
    polarplot(theta,rho);
    

    结果:
    在这里插入图片描述

    如果不用MATLAB呢?

    如果不用MATLAB呢?估计想得到作图能力强大的软件还有R语言、Python语言。Gnuplot也是一个不错的选择,直接用svg也可做图,或采用百度开发的Echarts也能完成。总之,方法很多。但用Gnuplot是借助了第3方作图工具了,用svg,Echarts或其他javascript包也要借助浏览器来显示。

    假如在有些情况下,不允许用任何第3方工具。不允许用windows编程架构或MFC,仅用C/C++语言结合windows的基本绘图函数,也能作图吗?

    我们追求简单。在此,我们仅用26行代码,在VC基本编译环境下,在控制台上也作出了如MATLAB那样的曲线结果。

    单纯依靠C函数完成作图

    windows下作图是使用gdi32.lib, user32.lib中的函数,所以编译要链接相应的库文件。库文件可以通过编译器选项链接,如:

    cl.exe  simpledraw.cpp /link user32.lib gdi32.lib
    

    也可在Cpp代码中加入预编指令:

    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    

    为了引用GDI函数(GDI:图形设备接口),需在main函数前加入:

    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    

    在main函数中,采用

    HWND hwnd = GetConsoleWindow(); 
    HDC hdc = GetDC(hwnd); 
    

    来获取控制台窗口句柄,从而操纵控制台使之以图形模式而不是字符模式作图。

    【注】还可有第2种方法获得控制台窗口的绘图hdc:先用GetConsoleTitle得到控制台窗口标题,再通过FindWindow函数借助标题获取窗口句柄HWND,之后再利用GetDC函数借助HWND获取HDC。这种做法可避免使用GetConsoleWindow函数。局部代码如下:

    TCHAR title[256];//控制台程序标题
    ......
    int main(int argc,char *argv[]) 
    { 
    	GetConsoleTitle(title, 256);//获取控制台标题
    	HDC hdc = GetDC(FindWindow(0, title));
    ......
    

    如果以默认画笔作图,则直接调用作图函数即可。最简单的作图函数就是画直线(LineTo)。如

    MoveToEx(hdc,400,300,NULL); //笔移到坐标(400,300)点
    LineTo(hdc,(int)x,(int)y); //画线到(x,y)坐标
    

    函数中,hdc就是指向控制台的作图句柄。

    好了,完整代码如下:

    #include<windows.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double rho, theta, x, y;//极坐标和直角坐标
    	system("cls");//清屏
    	system("color f0");//白底黑字
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	MoveToEx(hdc,400,300,NULL); //笔移到坐标原点(屏幕坐标400,300)
    	for (int i=0; i<10000; i++)
    	{
    		theta=0.01*i;
    		rho=exp(cos(theta))-2*cos(4*theta)+pow(sin(theta/12.0),5);
    		x=400+50*rho*cos(theta);
    		y=300-50*rho*sin(theta);
    		LineTo(hdc,(int)x,(int)y); //画线作图
    	}
    	TextOut(hdc, 550, 45, "美丽的蝴蝶结", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    其中,LineTo用以画线,TextOut用以在指定坐标处写字。退出程序前,用ReleaseDC(hwnd,hdc);释放作图资源,不用也行,但这是良好的编程习惯。默认画笔是黑色的,所以画图前将背景设为白色,并清屏。程序中,以屏幕坐标点(400,300)作为数学坐标原点。以极坐标计算后转换为直角坐标(并适当放大)作图。
    编译:
    1、直接在VC++6.0中编译。或
    2、控制台编译:

    cl.exe  simpledraw.cpp
    

    3、或链接上库文件(代码中可不要#pragma)

    cl.exe  simpledraw.cpp /link user32.lib gdi32.lib
    

    运行结果:
    在这里插入图片描述

    再来几个例子

    三叶草

    theta=0:0.01:2*pi;
    rho=1+cos(3*theta)+1.5*(sin(3*theta)).^2;
    polarplot(theta,rho);
    

    在这里插入图片描述

    Cpp代码:

    #include<windows.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double rho, theta, x, y;//极坐标和直角坐标
    	system("cls");//清屏
    	system("color f0");//白底黑字
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<10000; i++)
    	{
    		theta=0.01*i;
    		rho=1+cos(3*theta)+pow(1.5*(sin(3*theta)),2);
    		x=400+50*rho*cos(theta);
    		y=300-50*rho*sin(theta);
    		if(i==0){
    			MoveToEx(hdc,(int)x,(int)y,NULL); //笔移到起点
    		}else
    		{
    			LineTo(hdc,(int)x,(int)y); //画线作图
    		}	
    	}
    	TextOut(hdc, 550, 45, "美丽的三叶草", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    程序中演示如何设定画笔起始位置。
    运行结果:
    在这里插入图片描述

    正弦波:

    #include<windows.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double x, y;//直角坐标
    	system("cls");//清屏
    	system("color f0");//白底黑字
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<50000; i++)
    	{
    		x=0.01*i;
    		y=300+200*sin(2*3.14159*0.01*x);
    		if(i==0){
    			MoveToEx(hdc,130+(int)x,(int)y,NULL); //笔移到起点
    		}else
    		{
    			LineTo(hdc,130+(int)x,(int)y); //画线作图
    		}	
    	}
    	TextOut(hdc, 550, 45, "美丽的正弦波", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    在这里插入图片描述

    改画笔画刷(颜色、粗细、透明)

    也很简单。在作图前设置一下就行。

    HPEN hPen = CreatePen(PS_SOLID, 2, RGB(0,255,0));//添加画笔属性。2代表宽度,RGB(0,255,0)用于调颜色为绿色
    SelectObject(hdc, hPen);//画笔绑定到句柄
    HBRUSH hBrush=(HBRUSH)GetStockObject(NULL_BRUSH) ;//透明的画刷
    SelectObject (hdc, hBrush);
    

    完整代码:

    #include<windows.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double x, y;//直角坐标
    	system("cls");//清屏
    	system("color 0f");//黑底白字
    	HPEN hPen = CreatePen(PS_SOLID, 2, RGB(0,255,0));//添加画笔属性。2代表宽度,RGB(0,255,0)用于调颜色为绿色
        SelectObject(hdc, hPen);//画笔绑定到句柄
    	HBRUSH hBrush=(HBRUSH)GetStockObject(NULL_BRUSH) ;//透明
        SelectObject (hdc, hBrush);
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<50000; i++)
    	{
    		x=0.01*i;
    		y=300+200*sin(2*3.14159*0.01*x);
    		if(i==0){
    			MoveToEx(hdc,130+(int)x,(int)y,NULL); //笔移到起点
    		}else
    		{
    			LineTo(hdc,130+(int)x,(int)y); //画线作图
    		}	
    	}
    	TextOut(hdc, 550, 45, "美丽的正弦波", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    结果:
    在这里插入图片描述

    可见字色还是默认的。如果要改字背景和前景色,可用:

    SetBkColor(hdc, RGB(0,0,0));
    SetTextColor(hdc, RGB(0,255,0));
    //然后再写字:
    TextOut(hdc, 550, 45, "美丽的正弦波", 12);//写字
    

    图就成了这样:
    在这里插入图片描述

    如果还想设置字体,代码就要麻烦点,设置字体后,用SelectObject(hdc, font)选择字体即可,如下:

    ......
    	HFONT font = CreateFont(
    		28, // nHeight
    		0, // nWidth
    		50, // nEscapement
    		0, // nOrientation
    		FW_BOLD, // nWeight
    		TRUE, // bItalic
    		FALSE, // bUnderline
    		0, // cStrikeOut
    		ANSI_CHARSET, // nCharSet
    		OUT_DEFAULT_PRECIS, // nOutPrecision
    		CLIP_DEFAULT_PRECIS, // nClipPrecision
    		DEFAULT_QUALITY, // nQuality
    		DEFAULT_PITCH | FF_SWISS,
    		"Arial" // nPitchAndFamily Arial
    		);
    	SelectObject(hdc, font);
    ......
    

    作图结果:
    在这里插入图片描述

    可见,用C语言从底层控制作图是极为灵活的。

    完整代码:

    #include<windows.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double x, y;//直角坐标
    	system("cls");//清屏
    	system("color 0f");//黑底白字
    	HPEN hPen = CreatePen(PS_SOLID, 2, RGB(0,255,0));//添加画笔属性。2代表宽度,RGB(0,255,0)用于调颜色为绿色
        SelectObject(hdc, hPen);//画笔绑定到句柄
    	HBRUSH hBrush=(HBRUSH)GetStockObject(NULL_BRUSH) ;//透明
        SelectObject (hdc, hBrush);
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<50000; i++)
    	{
    		x=0.01*i;
    		y=300+200*sin(2*3.14159*0.01*x);
    		if(i==0){
    			MoveToEx(hdc,130+(int)x,(int)y,NULL); //笔移到起点
    		}else
    		{
    			LineTo(hdc,130+(int)x,(int)y); //画线作图
    		}	
    	}
    
    	HFONT font = CreateFont(
    		28, // nHeight
    		0, // nWidth
    		50, // nEscapement
    		0, // nOrientation
    		FW_BOLD, // nWeight
    		TRUE, // bItalic
    		FALSE, // bUnderline
    		0, // cStrikeOut
    		ANSI_CHARSET, // nCharSet
    		OUT_DEFAULT_PRECIS, // nOutPrecision
    		CLIP_DEFAULT_PRECIS, // nClipPrecision
    		DEFAULT_QUALITY, // nQuality
    		DEFAULT_PITCH | FF_SWISS,
    		"Arial" // nPitchAndFamily Arial
    		);
    	SelectObject(hdc, font);
    
    	SetBkColor(hdc, RGB(0,0,0));
        SetTextColor(hdc, RGB(0,255,0));
    	TextOut(hdc, 520, 45, "美丽的正弦波", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    稍作变化

    #include<windows.h>
    #include<stdio.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double rho, theta, x, y;//极坐标和直角坐标
    	HPEN hPen1,hPen2,hPen3;
    	hPen1 = CreatePen(PS_SOLID, 50, RGB(255,155,0));//添加画笔属性
    	hPen2 = CreatePen(PS_SOLID, 0, RGB(0,255,0));//添加画笔属性
    	hPen3 = CreatePen(PS_SOLID, 0, RGB(0,0,255));//添加画笔属性
    
    	const double pi=3.14159265;
    	system("cls");//清屏
    	system("color f0");//白底黑字
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<20000; i++)
    	{
    		theta=-pi+pi/10000.0*i;
    		rho=1+cos(3*theta)+pow(1.5*(sin(3*theta)),2);
    		x=400+50*rho*cos(theta);
    		y=300-50*rho*sin(theta);
    
            if(theta<(-pi/3.0)){
    		    SelectObject(hdc, hPen1);//画笔绑定到句柄
    		}	
    		if(theta>-pi/3.0 )
    		{
    			SelectObject(hdc, hPen2);//画笔绑定到句柄
    		}	
    		if(theta>pi/3.0)
    		{
    			SelectObject(hdc, hPen3);//画笔绑定到句柄
    		}
    
    		MoveToEx(hdc,400,300,NULL); //笔移到起点
    		LineTo(hdc,(int)x,(int)y); //画线作图
    	}
    	TextOut(hdc, 550, 45, "美丽的三叶草", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    

    在这里插入图片描述

    #include<windows.h>
    #include<stdio.h>
    #include<math.h>
    #pragma comment(lib,"user32.lib") 
    #pragma comment(lib,"gdi32.lib")
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    int main(int argc,char *argv[]) 
    { 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	double rho, theta, x, y;//极坐标和直角坐标
    	HPEN hPen1,hPen2,hPen3;
    	hPen1 = CreatePen(PS_SOLID, 0, RGB(255,155,0));//添加画笔属性
    	hPen2 = CreatePen(PS_SOLID, 0, RGB(0,255,0));//添加画笔属性
    	hPen3 = CreatePen(PS_SOLID, 0, RGB(0,0,255));//添加画笔属性
    
    	const double pi=3.14159265;
    	system("cls");//清屏
    	system("color f0");//白底黑字
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//画个边框
    	for (int i=0; i<20000; i++)
    	{
    		theta=-pi+pi/10000.0*i;
    		rho=1+cos(7*theta)+pow(1.5*(sin(3*theta)),2);
    		x=400+50*rho*cos(theta);
    		y=300-50*rho*sin(theta);
    
            if(theta<(-pi/3.0)){
    		    SelectObject(hdc, hPen1);//画笔绑定到句柄
    		}	
    		if(theta>-pi/3.0 )
    		{
    			SelectObject(hdc, hPen2);//画笔绑定到句柄
    		}	
    		if(theta>pi/3.0)
    		{
    			SelectObject(hdc, hPen3);//画笔绑定到句柄
    		}
    
    		MoveToEx(hdc,400,300,NULL); //笔移到起点
    		LineTo(hdc,(int)x,(int)y); //画线作图
    	}
    	TextOut(hdc, 550, 45, "美丽的三叶草", 12);//写字
        ReleaseDC(hwnd,hdc);
    	return 0; 
    }
    
    

    在这里插入图片描述

    还能有比这更简单的C代码吗?

    也许如下,画个矩形,只有8行,核心代码3行–建立图形句柄调并用之。

    #include<windows.h>
    extern "C" WINBASEAPI HWND WINAPI GetConsoleWindow ();
    void main(int argc,char *argv[]){ 
    	HWND hwnd = GetConsoleWindow(); 
    	HDC hdc = GetDC(hwnd); 
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//用句柄hdc画矩形
    }
    
    

    或是

    #include<windows.h>
    TCHAR title[256];//控制台程序标题
    void main(int argc,char *argv[]){ 
    	GetConsoleTitle(title, 256);//获取控制台标题
    	HDC hdc = GetDC(FindWindow(0, title));
    	Rectangle(hdc, 120, 50, 120+600, 50+500);//用句柄hdc画矩形
    }
    
    

    编译要加link选项(在VC++6.0集成环境中默认已加了):

    cl.exe  simpledraw.cpp /link user32.lib gdi32.lib
    

    还可让图形动起来

    让图形动起来的办法就是不断定时清屏、重绘。像放电影那样。然而,直接向屏幕绘图是耗时的,在屏幕上画很多线,将多次调用绘图接口,这样,不但绘制效率低(慢),而且视觉上闪烁感明显。

    解决的办法是:先在内存中将多条线的绘制结果准备好,再一次性调用向屏幕绘图的接口。这称为双缓冲法。即先向内存中的缓冲区作图,再将作图结果一次性拷贝到屏幕的显示区。

    详細方法下篇再讲。

    附录: GDI绘图基本步骤总结

    一、获得绘图的窗口句柄

    方法:多种

    
    HWND FindWindow(LPCTSTR lpClassName, LPCTSTR lpWindowName)
    
    HWND FindWindowEx(HWND hwndParent, HWND hwndChildAfter,LPCTSTR lpClassName, LPCTSTR lpWindowName)
    
    HWND  GetConsoleWindow()   //直接找控制台窗口句柄
    
    HWND WindowFromPoint(POINT& Point)
    
    BOOL CALLBACK EnumChildProc(HWND hwnd,LPARAM lParam)   
    BOOL CALLBACK EnumChildWindows(HWND hWndParent, WNDENUMPROC lpEnumFunc,LPARAM lParam)   
    BOOL CALLBACK EnumWindows(WNDENUMPROC lpEnumFunc, LPARAM lParam)
    BOOL CALLBACK EnumWindowsProc(HWND hwnd, LPARAM lParam)
    

    二、由窗口句柄得到设备环境句柄HDC

    方法(3种):BeginPaint、GetWindowDC、GetDC。

    这些函数都需要步骤一中的HWND的句柄。调用这些函数后要释放句柄,相应的有EndPaint、ReleaseDC进行清理。

    1、 采用BeginPaint获取HDC

    HDC hdc;
    PAINTSTRUCT ps;// 保存
    hdc = ::BeginPaint( hwnd, &ps );
    // 此处添加绘图代码
    ::EndPaint( hwnd, &ps );
    

    说明:获得的hdc的有效区域仅限于客户区无效区域的设备环境句柄,不包括标题栏、边框等

    2、 采用GetWindowDC获取HDC

    HDC hdc = ::GetWindowDC( hwnd );
    // 此处添加绘图代码
    ::ReleaseDC( hwnd, hdc );
    

    说明:绘制区域是整个窗口**(边框、标题栏、客户区的总和)**。

    3、 采用GetDC获取HDC

    HDC hdc = ::GetDC( hwnd );
    // 此处添加绘图代码
    ::ReleaseDC( hwnd, hdc );
    

    说明:获得的hdc的有效区域仅限于客户区有效区域的设备环境句柄,不包括标题栏、边框等。

    三、图形绘制方法

    1、 画笔CreatePen

    绘画之前先选择画笔,画笔的功能主要是绘制边框,其函数原型如下:

    WINGDIAPI HPEN WINAPI CreatePen(
    __in int iStyle,           // 画笔的类型,比如是实线,还是虚线等等。
    __in int cWidth,          // 线的宽度。
    __in COLORREF color       // 线的颜色。
    );
    
    其中:
    // iStyle参数可选值:
    PS_SOLID       = 0;// 实线
    PS_DASH        = 1;// 段线; 要求笔宽<=1
    PS_DOT         = 2;// 点线; 要求笔宽<=1
    PS_DASHDOT     = 3;// 线、点; 要求笔宽<=1
    PS_DASHDOTDOT  = 4;// 线、点、点; 要求笔宽<=1
    PS_NULL        = 5;// 不可见
    PS_INSIDEFRAME = 6;// 实线; 但笔宽是向里扩展
    

    返回值为画笔类型,SelectObject函数选中。选中后,返回原来画刷的句柄用来恢复时使用。图形绘制完毕后使用DeleteObject函数将其释放。

    【SelectObject函数说明】:

    函数功能:该函数选择一对象到指定的设备上下文环境中,该新对象替换先前的相同类型的对象。

    函数原型:
    HGDIOBJ SelectObject(HDC hdc, HGDIOBJ hgdiobj);

    参数:

    hdc:设备上下文环境的句柄。
    hgdiobj:被选择的对象的句型,该指定对象必须由如下的函数创建。
    
    位图:CreateBitmap, CreateBitmapIndirect, CreateCompatible Bitmap, CreateDIBitmap, CreateDIBsection(只有内存设备上下文环境可选择位图,并且在同一时刻只能一个设备上下文环境选择位图)。
    
    画笔:CreateBrushIndirect, CreateDIBPatternBrush, CreateDIBPatternBrushPt, CreateHatchBrush, CreatePatternBrush, CreateSolidBrush。
    
    字体:CreateFont, CreateFontIndirect。
    笔:CreatePen, CreatePenIndirect。
    
    区域:CombineRgn, CreateEllipticRgn, CreateEllipticRgnIndirect, CreatePolygonRgn, CreateRectRgn, CreateRectRgnIndirect。
    

    返回值:如果选择对象不是区域并且函数执行成功,那么返回值是被取代的对象的句柄。

    【DeleteObject函数说明】:

    函数功能:该函数删除一个逻辑笔、画笔、字体、位图、区域或者调色板,释放所有与该对象有关的系统资源,在对象被删除之后,指定的句柄也就失效了。

    函数原型:BOOL DeleteObject(HGDIOBJ hObject);

    参数:

    hObject:逻辑笔、画笔、字体、位图、区域或者调色板的句柄。
    返回值:成功,返回非零值;如果指定的句柄无效或者它已被选入设备上下文环境,则返回值为零。

    2、 画刷

    画刷的功能主要是填充区域内的颜色,创建画刷的方法如下:

    A、CreateSolidBrush函数

    函数功能:该函数创建一个具有指定颜色的逻辑刷子。
    函数原理:HBRUSH CreateSolidBrush(COLORREF crColor);
    参数:
    crColor:指定刷子的颜色。
    返回值:如果该函数执行成功,那么返回值标识一个逻辑实心刷子;如果函数失败,那么返回值为NULL。

    B、GetStockObject函数

    函数功能:该函数检索预定义的备用笔、刷子、字体或者调色板的句柄。
    函数原型:HGDIOBJ GetStockObject(int fnObject);
    参数:

    fnObject:指定对象的类型,该参数可取如下值之一;
    BLACK_BRUSH:黑色画笔;
    DKGRAY_BRUSH:暗灰色画笔;
    DC_BRUSH:在Windows98,Windows NT 5.0和以后版本中为纯颜色画笔,缺省色为白色,可以用SetDCBrushColor函数改变颜色,更多的信息参见以下的注释部分。
    GRAY_BRUSH:灰色画笔;
    HOLLOW_BRUSH:空画笔(相当于HOLLOW_BRUSH);
    LTGRAY_BRUSH:亮灰色画笔;
    NULL_BRUSH:空画笔(相当于HOLLOW_BRUSH);
    WHITE_BRUSH:白色画笔;BLACK_PEN:黑色钢笔;
    DC_PEN:在Windows98、Windows NT 5.0和以后版本中为纯色钢笔,缺省色为白色,使用SetDCPenColor函数可以改变色彩,更多的信息,参见下面的注释部分。
    WHITE_PEN:白色钢笔;
    ANSI_FIXED_FONT:在Windows中为固定间距(等宽)系统字体;
    ANSI_VAR_FONT:在Windows中为变间距(比例间距)系统字体;
    DEVICE_DEFAUCT_FONT:在WindowsNT中为设备相关字体;
    DEFAULT_GUI_FONT:用户界面对象缺省字体,如菜单和对话框;
    OEM_FIXED_FONT:原始设备制造商(OEM)相关固定间距(等宽)字体;
    SYSTEM_FONT:系统字体,在缺省情况下,系统使用系统字体绘制菜单,对话框控制和文本;
    SYSTEM_FIXED_FONT:固定间距(等宽)系统字体,该对象仅提供给兼容16位Windows版本;
    DEFAULT_PALETTE:缺省调色板,该调色板由系统调色板中的静态色彩组成。
    返回值:如果成功,返回值标识声请的逻辑对象,如果失败,返回值为NULL。

    C、CreateHatchBrush函数

    函数功能:该函数可以创建一个具有指定阴影模式和颜色的逻辑刷子。
    函数原型:HBRUSH CreateHatchBrush(int fnStyle, COLORREF clrref);
    参数:
    fnStyle:指定刷子的阴影样式。该参数可以取下列值,这些值的含义为:
    HS_BDIAGONAL:表示45度向下,从左至右的阴影;
    HS_CROSS:水平和垂直交叉险影;
    HS_DIAGCROSS:45度交叉阴影;
    HS_FDIAGONAL:45度向上,自左至右阴影;
    HS_HORIZONTAL:水平阴影;
    HS_VERTICAL:垂直阴影。
    cirref:指定用于阴影的刷子的前景色。
    返回值:如果函数执行成功,那么返回值标识为逻辑刷子;如果函数执行失败,那么返回值为NULL。

    画刷的选中和释放,请参照画笔。

    3、 点SetPixel

    函数功能:该函数将指定坐标处的像素设为指定的颜色。
    函数原型:COLORREF SetPixel(HDC hdc, int X, int Y, COLORREF crColor);
    参数:
    hdc:设备环境句柄。
    X:指定要设置的点的X轴坐标,按逻辑单位表示坐标。
    Y:指定要设置的点的Y轴坐标,按逻辑单位表示坐标。
    crColor:指定要用来绘制该点的颜色。
    返回值:如果函数执行成功,那么返回值就是函数设置像素的RGB颜色值。这个值可能与crColor指定的颜色有不同,之所以有时发生这种情况是因为没有找到对指定颜色进行真正匹配造成的;如果函数失败,那么返回值是0。

    4、 直线MoveToEx、LineTo

    A、 MoveToEx
    函数功能:将当前位置指定为特定的某一点
    函数原型:BOOL MoveToEx( __in HDC hdc, __in int X, __in int Y, __out LPPoint lpPoint )
    参数:
    hdc:设备环境句柄。
    X:指定要设置的点的X轴坐标,按逻辑单位表示坐标。
    Y:指定要设置的点的Y轴坐标,按逻辑单位表示坐标。
    lpPoint:指向一个POINT结构,用来接收前一位置,为空时,当前位置不被返回。
    返回值:执行成功返回非零,否则返回值为零。

    B、 LineTo

    函数功能:从当前点到目标点进行画线。
    函数原型:BOOL LineTo( int x, int y )
    参数说明:
    X:目标点的横坐标。
    Y:目标点的纵坐标。
    返回值:成功非零,其它返回零。

    5、 矩形Rectangle

    函数功能:该函数画一个矩形,用当前的画笔画矩形轮廓,用当前画刷进行填充。
    函数原型:BOOL Rectangle(HDC hdc, int nLeftRect, int nTopRect, int nRightRect, int nBottomRect);
    参数:
    hdc:设备环境句柄。
    nLeftRect:指定矩形左上角的逻辑X坐标。
    nTopRect:指定矩形左上角的逻辑Y坐标。
    nRightRect:指定矩形右下角的逻辑X坐标。
    nBottomRect:指定矩形右下角的逻辑Y坐标。
    返回值:如果函数调用成功,返回值非零,否则返回值为0。

    6、 椭圆Ellipse

    函数功能:该函数画一个椭圆形,用当前的画笔画矩形轮廓,用当前画刷进行填充。
    函数原型:
    BOOL Ellipse( HDC hdc, int x1, int y1, int x2, int y2 )

    参数:
    hdc:设备环境句柄。
    x1:指定椭圆形左上角的逻辑X坐标。
    y1:指定椭圆形左上角的逻辑Y坐标。
    x2:指定椭圆形右下角的逻辑X坐标。
    y2:指定椭圆形右下角的逻辑Y坐标。
    返回值:如果函数调用成功,返回值非零,否则返回值为0。

    展开全文
  • Photon是功能齐全的图形计算器,能够进行数值计算以及图形功能,包括参数图和极坐标图。 除了常规图形绘制外,Photon还可以绘制积分,一阶和二阶导数,跟踪,动画,矩阵数学以及度/弧度模式。 光子可以找到根,...
  • VB程序设计补修课复习资料
  • 概述:本文为使用Iocomp工控图表工具绘制实时曲线探索及研究教程,为大家介绍了Iocomp控件、实时曲线绘制方法、Iocomp界面操作,属性分类等。帮助学习者更好的运用Iocomp。

    概述本文为使用Iocomp工控图表工具绘制实时曲线探索及研究教程,为大家介绍了Iocomp控件、实时曲线绘制方法、Iocomp界面操作,属性分类等。帮助学习者更好的运用Iocomp。

    [摘要]数据采集是控制系统最常见的任务,对于大量的实时采集数据采用曲线加以分析已成为很重要的一种手段。文章将介绍如何用 Iocomp控件实现控制软件中的实时曲线的设计与绘制,并结合实例程序加以具体说明。实践表明,该方法简单可靠,对工业实时控制应用有一定借鉴意义。


    Iocomp控件:


          0.引言

    在工业控制领域,需要进行大量的数据处理和可视化显示。实时监测软件中,常常需要将采集到的数据实时显示到界面上来,以便于工作人员观测,及时发现问题和解决问题,通常还要求曲线可伸缩、可漫游、可取值,可若干条曲线的比较,以增强其可分析性。传统的控制软件开发工具多用 DOS或 Windows下的 C语言开发,导致系统开发周期长,可维护性差,并且不具有标准的 Windows图形用户界面。VB作为一种 Windows软件开发工具,既具有效率高界面友好的功能,又可以使用 DLL来实现 I/O端口的输入功能,还可以通过 API函数或 Mscomm控件实现串口通信。而本文将详细介绍一种基于 VB的 Iocomp控件实现实时曲线的绘制方法,该方法简便易行,编程也比较简单,在实际应用中得到了良好的效果。

    1、Iocomp Software简介

    Iocomp Software是一个让 Iocomp公司引以为豪的全新的、100%托管的、领先的、高速的、易用的、能实时绘制的控件。拥有许多在其他同类图表控件产品中所不能找到的全新特征功能及性能。Iocomp Software特征如下:

    属性定制编辑器:每一个控件都带有一个将属性以逻辑形式分组的属性定制编辑器,它们使得用户能够轻松的设置控件属性。在属性窗口中可以无限的搜索,属性窗口也可以无限的层叠。

    值相关联:大多数的控件都拥有一些值,它们在设置前可能是相关联的,这就提供了一个相关联的机制以及程序灵活性。

    实时 -高效:该产品的工具控件是当今最高效的。具有实时显示、缩放、滚动、即使是在进行数据绘制时、不限制 X以及 Y坐标轴、不限制通道、曲线拟合、直角坐标轴、可视化的布局管理器、支持 EMF、BMP、 JPG以及 TXT格式的文件输出。

    2、实时曲线绘制方法

    实时曲线的绘制方法多种多样,根据对曲线的要求,我们可以采用不同的方法来绘制,从而达到最佳的曲线效果,以下列出了几种常用的绘制实时曲线的方法:

    方法一采用 TeeChart实现。TeeChart Pro是一款提供上百种 2D和 3D图形风格、40种数学和统计功能、加上无限制的轴和 22种调色板组件供选择的绘图控件。TeeChart还包括一个强大的、完整的编辑对话框,几乎可用于每个组件和子组件,允许你快速的设计复杂图表应用程序。图表编辑器通过 TeeCommander组件进一步得到增强,它提供一次点击访问图表编辑器和共同特征。

    方法二在 VB中绘制实时曲线是比较难的,一般要应用第三方控件或是 Windows API函数来完成,但是如果你对实时曲线的要求不是很高,只要能表示出当前的一般情况的话,我们可以直接应用 VB提供给我们的空间来完成。

    方法三可以采用工控组态软件来实现。组态软件具有动画显示、流程控制、数据采集、设备控制与输出、工程报表、数据与曲线等强大功能,在自动控制中占据主力军的位置,已逐渐成为工业自动化的灵魂。

    方法四采用 Iocomp控件来实现。结合 Iocomp Components图形仪表组件,易于实现 VB组态。

    结合以上四种方法的总结与比较,采用 Iocomp控件来实现实时曲线的绘制是最佳的,在工业控制软件当中,它将为更多的程序员提供更广的应用领域和发展空间。


          3、Iocomp控件介绍

    Iocomp控件主要有三种组件:⑴iPlot ⑵iXYPlot ⑶iScope iPlot是一款即时绘图组件,支持具有连续递增的 X坐标的数据序列的绘图。其典型的用途是图表记录或滚动图表类型的应用。该组件可
    用于所有的绘图应用。iPlot组件的应用使得绘图程序得到了最优化, iPlot组件还可提供高速接入和绘图方法。

    iXYPlot是一个即时绘图组件,它支持具备任意 X、Y坐标值的数据。

    iScope是一个真正的实时模拟和数字范围的组件。支持实时处理和显示的数据时,5MHz的信号输入率如果是连续的数据或使用更高的数据传输率采用间断数据。

    iPlot组件和 iXYPlot组件的等级是相同的,除了通道的对象。其层次结构图如图 1。


    图 1层次结构图常用属性简介如下表:


    Channel属性,可以访问指定的通道。可以使用 AddChannel, DeleteChannel和 RemoveAllChannel方法在运行时间中添加或者移除某个通道。在设计时,可以使用内置的属性编辑器来改变通道。
          Annotation属性,可以访问指定的注释。可以使用 AddAnnotation, DeleteAnnotation和 RemoveAllAnnotation方法在运行时间中添加或者移除某个注释。注释无法在设计时间中添加。
          X-axis属性,您可以访问指定的 X轴线。您可以使用 AddXAxis, DeleteXAxis和 RemoveAllXAxes方法在运行时间中添加或者移除某个 X轴。在设计时间,可以使用内置的属性编辑器来改变某个 X轴。
          YAxis属性,可以访问指定的 Y轴线。可以使用 AddYAxis, Delete YAxis和 RemoveAllYAxes方法在运行时间中添加或者移除某个 Y轴。在设计时间,可以使用内置的属性编辑器来改变某个 Y轴。
          DataView属性通过索引号,可以访问数据视图。该功能接口支持多个数据视图,但是在当前版本中未提供该功能支持,期待在下面版本中推出。
          ToolBar属性通过索引号,可以访问工具栏。该功能接口支持多个工具栏,但是在当前版本中未提供该功能支持,期待在下面版本中推出。
          Legend属性通过索引号,可以访问图标符号。该功能接口支持多个图标符号,但是在当前版本中未提供该功能支持,期待在下面版本中推出。

    4、编辑主界面

    对界面的设置是必不可少的,它可以更直观更简捷的为我们提供绘制实时曲线的各个属性功能,我们可以直接修改程序来进行对一些属性要求的设置,也可以直接在面板的分项属性中进行设置,对各属性

    不同的设置得到的结果显示也将不一样,可以根据设计员自己的要求来具体设置。下图 2是显示了本文所应用的 iPlot组建的编辑主界面。


    5、具体实例

    下面来看一下如何采用 Iocomp控件中的 iPlot组件来实现实时曲线的绘制,编写程序既简单又方便。
    采用 Iocomp控件中的 iPlot组件绘制实时曲线,其实时测量值曲线绘制显示界面如图 3。


    实时曲线反映的是现场数据的实时性和当前趋势,绘制实时采集数据曲线是为了实时观测,以便掌握实时采集数据变动的趋势,使曲线显示效果最佳,因此在实现时需显示曲线的动态变化,当前点在曲线的最右端显示,而整个曲线动态地向左移动。实时曲线在动态的移动时,测量值、峰值、谷值这三个通道同时相应的显示具体采样数值,一目了然。由于篇幅有限,程序只保留核心部分。编写程序如下:

    Private Sub Form_Load()
    iPlotX1.Channel(0).TitleText =  " 测量值"
    iPlotX1.Channel(1).TitleText =  " 峰值"
    iPlotX1.Channel(2).TitleText =  " 谷值"
    Call Comm_initial
    XValue = 0
    If Right(App.Path, 1) = "\" Then
    fpname = App.Path &"data\"
    Else
    fpname = App.Path &"\data\"
    End If
    ComD1.InitDir = fpname
    End Sub

    以上程序完成对 iPlotX控件的初始化,初始化具有三个通道,名字分别为测量值、峰值、谷值,并调用 Comm_initial函数完成对串口的初始化,设定 App.Path & "data\"为采集数据存放路径。

    Private Sub Timer1_Timer()
    Dim a
    Dim i%, j%, k%, l%
    Dim bjsta As String
    Dim otime As Long
    Dim delayt As Integer
    delayt = 60
    i=0 '命令重发次数计数变量
    j=0
    k=0
    l=0
    fs1:
    MSComm.InBufferCount = 0
    MSComm.OutBufferCount = 0 '清空输出缓冲区
     
    MSComm.Output =  "#01" & vbCr '测量值读取 otime = GetTickCount i=i+1 Do
    a = DoEvents() Loop Until MSComm.InBufferCount >= 10 Or GetTickCount >= otime +
    delayt clclz = MSComm.Input If Left(clclz, 1) =  "=" And IsNumeric(Mid(clclz, 2, 6)) Then
    lbclz.Caption = Mid(clclz, 2, 6) clclz = lbclz.Caption
    Else If i > 10 Then GoTo fs2 GoTo fs1
    End If ……………… //省略程序数据采集峰值、谷值,其采集方法同上 fs5:
    YValue = Val(lbclz.Caption) iPlotX1.Channel(0).AddYNow Yvalue YValue = Val(Lbfz.Caption) iPlotX1.Channel(1).AddYNow YValue YValue = Val(Lbgz.Caption) iPlotX1.Channel(2).AddYNow YValue YValue = yboutv2 XValue = XValue + 1 MSComm.InBufferCount = 0 MSComm.OutBufferCount = 0 ……………… //数据保存部分省略
    End Sub

    以上程序为实时采集函数,其功能:用 MSComm实现与串口的通信, MSComm.Output = "#01" & vbCr语句为读仪表命令。每次采样发出 3条仪表读取测量值命令,1条报警开关读取命令。采样频率在有线通信和无线通信模式下最高每秒可以达到 15次,在无线通信出错较大情况下可以保证每秒正确采样 8次,通信采用的强烈的纠错方法;采样数据测量值、峰值、谷值送入 Iocomp控件实时显示,程序会在每次采样过后将数据存入本软件的 data目录下。

    6.结束语

    以上程序采用 Iocomp控件实现实时曲线的绘制,曲线的采样频率可以每秒达到 50次,最高可达 100次,灵敏度极高,绘制出的曲线也比较流畅,它不仅可用在工业控制的历史数据处理,而且可以广泛地适用于商业、管理及很多有大量数据需曲线显示的应用。
    因此,Iocomp Components图形仪表组件对于构建逼真的人机界面、处理实时数据将非常有用,相信它们会被广泛地应用到越来越多的领域和行业当中。由于控件具有界面友好、操作简单、可扩展性强等特点,因此使用这些组件可以非常方便地扩展出风格多种多样的数据图形,整个设计过程都无需编程。


    作者:陕西理工学院物理系 翟世磊 李明波 李福 刘东

    文章转载自:慧都控件网 [http://www.evget.com]










    展开全文
  • 已知一点,另一点相对于该点的角(弧度)和轴长度,求另一点的位置 两点法绘制圆、三点法绘制绘制椭圆、椭圆弧的中心线 绘制面域中心线 交换两个数组变量 给任用一个实体绘制边框 将三个变量转换成一个点...
  • 曲线拟合数值方法求解,主要是利用最小二乘法及其线性化法、样条函数插值、三角多项式以及贝塞尔函数法进行拟合,并分析它们的原理以及应用场景。然后,分析各算法流程和核心计算代码模块组织结构,然后对各算法的...

    《数值计算方法》系列总目录

    第一章 误差序列实验
    第二章 非线性方程f(x)=0求根的数值方法
    第三章 CAD模型旋转和AX=B的数值方法
    第四章 插值与多项式逼近的数值计算方法
    第五章 曲线拟合的数值方法
    第六章 数值微分计算方法
    第七章 数值积分计算方法
    第八章 数值优化方法



    一、算法原理

    1、最小二乘拟合曲线

      在科学和工程试验中,经常产生一组数据( x 1 , y 1 x_1,y_1 x1,y1),…,( x N , y N x_N,y_N xN,yN),这里的横坐标{ x k x_k xk}是明确的的。数值方法的目标之一是确定一个将这些变量联系起来的函数 y = f ( x ) y = f(x) y=f(x),即通过算法形成一种数学刻画的公式来寻求其中的规律。通常会选择一类可用的函数并确定它们的系数,然而这类函数不是随便选取的,必须根据实际的实验数据进行判断,选取最符合的函数模型,因此选择函数的可能性是多种多样的。一般会根据物理情况采用一个基本数学模型来确定函数的形式。这一节主要分析线性函数,其形式为
    y = f ( x ) = A x + B y = f(x) = Ax + B y=f(x)=Ax+B  根据插值多项式逼近原理,如果所有的数值{ x k x_k xk},{ y k y_k yk}已知有多位有效数字精度,则能成功地使用多项式插值;否则,不能使用多项式插值。因为插值多项式是根据已有的精确的实验数据进行多项式拟合的,其拟合结果一定会经过插值点,如果插值点数据本身存在较大误差,则逼近结果就会出现更大误差。然而,许多试验数据可能并没有如此高的精度,而且通常在试验中还存在试验误差,所以尽管记录的{ x k x_k xk},{ y k y_k yk}有多位有效数字,但还是与其真实值存在范围外的误差,即真实值 f ( x k ) f(x_k) f(xk)满足 f ( x k ) = y k + e k f(x_k)=y_k+e_k f(xk)=yk+ek,其中, e k e_k ek表示测量误差,因此需要寻找一种经过测试点附近(不总是穿过)的最佳线性逼近表达式。下面首先对误差(又称偏差或残差)概念进行一个介绍:
    e k = f ( x k ) − y k , f o r 1 ≤ k ≤ N {e_k} = f({x_k}) - {y_k},{\rm{ }}for{\rm{ }}1 \le k \le N ek=f(xk)yk,for1kN  其中, y k y_k yk测量实验值, f ( x k ) f(x_k) f(xk)表示真实值,误差还有其他表示如下形式:
    E ∞ ( f ) = max ⁡ 1 ≤ k ≤ N { ∣ f ( x k ) − y k ∣ } {E_\infty }(f) = \mathop {\max }\limits_{1 \le k \le N} \{ |f({x_k}) - {y_k}|\} E(f)=1kNmax{f(xk)yk} E 1 ( f ) = 1 N ∑ k = 1 N ∣ f ( x k ) − y k ∣ {E_1}(f) = \frac{1}{N}\sum\limits_{k = 1}^N {\left| {f({x_k}) - {y_k}} \right|} E1(f)=N1k=1Nf(xk)yk E 2 ( f ) = ( 1 N ∑ k = 1 N ∣ f ( x k ) − y k ∣ 2 ) 1 / 2 {E_2}(f) = {\left( {\frac{1}{N}\sum\limits_{k = 1}^N {{{\left| {f({x_k}) - {y_k}} \right|}^2}} } \right)^{1/2}} E2(f)=(N1k=1Nf(xk)yk2)1/2  那么如何确定系数 A , B A,B AB使得曲线拟合达到我们想要的效果呢?其中,最小二乘法是在进行曲线拟合时常用的一种求解拟合曲线表达式的方法,无论是拟合何种曲线,其系数求解基本思想都是基于最小二乘法的。分析可知,设 是一个N个点的集合,其中横坐标是确定的。那么,理想的最小二乘拟合曲线 的拟合结果应满足是均方根误差最小的曲线。那么,所选定的 A , B A,B AB构成的拟合曲线应当使 E 2 ( f ) E_2(f) E2(f)的值当且仅当所选取的 A , B A,B AB值处最小。换言之,也就是说数据点到曲线的垂直距离平方和最小,即误差的平方和,其最小值为 N ( E 2 ( f ) 2 ) = ∑ k = 1 N ∣ A x k + B − y k ∣ 2 N({E_2}{(f)^2}) = \sum\limits_{k = 1}^N {{{\left| {A{x_k} + B - {y_k}} \right|}^2}} N(E2(f)2)=k=1NAxk+Byk2。那么,基于实验结果的所求出选取 A , B A,B AB的值的表达式到底是怎样的呢?因此,有如下推导过程:
      设 { ( x k , y k ) } k = 1 N \{ ({x_k},{y_k})\} _{k = 1}^N {(xk,yk)}k=1N有N个点,其中横坐标是确定的,当数据点到曲线的垂直距离平方和最小时, A , B A,B AB取值点应在误差表达式的极值处,即表达式对于 A , B A,B AB两变量的偏导数为0,即
    ∂ E ( A , B ) ∂ A = ∑ k = 1 N 2 ( A x k + B − y k ) x k = 0 ∂ E ( A , B ) ∂ B = ∑ k = 1 N 2 ( A x k + B − y k ) = 0 \frac{{\partial E(A,B)}}{{\partial A}} = \sum\limits_{k = 1}^N {2\left( {A{x_k} + B - {y_k}} \right){x_k}} = 0\\ \frac{{\partial E(A,B)}}{{\partial B}} = \sum\limits_{k = 1}^N {2\left( {A{x_k} + B - {y_k}} \right)} = 0 AE(A,B)=k=1N2(Axk+Byk)xk=0BE(A,B)=k=1N2(Axk+Byk)=0  通过移项变换即得到关于系数 A , B A,B AB的表达式如下
    ( ∑ k = 1 N x k 2 ) A + ( ∑ k = 1 N x k ) B = ∑ k = 1 N x k y k ( ∑ k = 1 N x k ) A + N B = ∑ k = 1 N y k \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)B = \sum\limits_{k = 1}^N {{x_k}{y_k}} \\ \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)A + NB = \sum\limits_{k = 1}^N {{y_k}} (k=1Nxk2)A+(k=1Nxk)B=k=1Nxkyk(k=1Nxk)A+NB=k=1Nyk  并称式(5)表示的方程组为正规方程。得到的曲线称为最小二乘拟合曲线。
      除此之外,对于形如 f ( x ) = A x M f(x) = A{x^M} f(x)=AxM幂函数拟合,其中M为已知常数,根据线性曲线的拟合推导过程,写出幂函数数据点到曲线的垂直距离平方和表达式,即
    E ( A ) = ∑ k = 1 N ( A x k M − y k ) 2 E(A) = \sum\limits_{k = 1}^N {{{\left( {Ax_k^M - {y_k}} \right)}^2}} E(A)=k=1N(AxkMyk)2  求出A为其极小值点时候的值,即表达式对于A的偏导数为0,即
    E ′ ( A ) = 2 ∑ k = 1 N ( A x k M − y k ) x k M = 2 ∑ k = 1 N ( A x k 2 M − x k M y k ) = 0 E'(A) = 2\sum\limits_{k = 1}^N {\left( {Ax_k^M - {y_k}} \right)} x_k^M = 2\sum\limits_{k = 1}^N {\left( {Ax_k^{2M} - x_k^M{y_k}} \right) = 0} E(A)=2k=1N(AxkMyk)xkM=2k=1N(Axk2MxkMyk)=0  最终可以得出幂函数拟合曲线的系数A的表达式为

    2、曲线拟合原理

      对于不是线性关系的曲线拟合,可以对其进行线性化操作再利用最小二乘法求出其系数的表达式即可。也就是说,对于一组不具有线性关系的数据点集合,可以对其进行诸如取对数、指数的方法将其变换为具有线性关系的表达式,即将数据点 x k , y k {xk,yk} xk,yk变换成 X Y XY XY平面上具有线性关系的点集 X k , Y k {Xk,Yk} Xk,Yk,再利用最小二乘法求出变换后的曲线表达式,最后再变换为原变量的表达式。这个过程称为数据线性化。例如,对于非线性函数[y = C{e^{Ax}}],可以对其进行取对数变化,即
    ln ⁡ ( y ) = A x + ln ⁡ ( C ) \ln (y) = Ax + \ln (C) ln(y)=Ax+ln(C)  并引入变量, ,得到变换后的线性关系式 ,根据最小二乘法的结果,系数 A , B A,B AB表达式为
    ( ∑ k = 1 N X k 2 ) A + ( ∑ k = 1 N X k ) B = ∑ k = 1 N X k Y k ( ∑ k = 1 N X k ) A + N B = ∑ k = 1 N Y k \left( {\sum\limits_{k = 1}^N {X_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {X_k^{}} } \right)B = \sum\limits_{k = 1}^N {{X_k}{Y_k}} \\ \left( {\sum\limits_{k = 1}^N {X_k^{}} } \right)A + NB = \sum\limits_{k = 1}^N {{Y_k}} (k=1NXk2)A+(k=1NXk)B=k=1NXkYk(k=1NXk)A+NB=k=1NYk C = e B C = {e^B} C=eB  则可以得到拟合曲线的表达式。然而,也可以直接写出数据点到曲线的垂直距离平方和表达式,并求出其极值点,即为拟合曲线的系数表达式,其推导过程如下:
    E ( A , C ) = ∑ k = 1 N ( C e A x k − y k ) 2 l e t { ∂ E ∂ A = 2 ∑ k = 1 N ( C e A x k − y k ) ( C x k e A x k ) = 0 ∂ E ∂ C = 2 ∑ k = 1 N ( C e A x k − y k ) ( e A x k ) s o { C ∑ k = 1 N x k e 2 A x k − ∑ k = 1 N y k x k e A x k = 0 C ∑ k = 1 N e A x k − ∑ k = 1 N y k e A x k = 0 E(A,C) = \sum\limits_{k = 1}^N {{{\left( {C{e^{A{x_k}}} - {y_k}} \right)}^2}} \\ {\rm{let}}\left\{ \begin{array}{l} \frac{{\partial E}}{{\partial A}} = 2\sum\limits_{k = 1}^N {\left( {C{e^{A{x_k}}} - {y_k}} \right)} \left( {C{x_k}{e^{A{x_k}}}} \right){\rm{ = 0}}\\ \frac{{\partial E}}{{\partial C}} = 2\sum\limits_{k = 1}^N {\left( {C{e^{A{x_k}}} - {y_k}} \right)} \left( {{e^{A{x_k}}}} \right) \end{array} \right.\\ {\rm{so}}\\ \left\{ \begin{array}{l} C\sum\limits_{k = 1}^N {{x_k}{e^{2A{x_k}}} - } \sum\limits_{k = 1}^N {{y_k}{x_k}{e^{A{x_k}}}} = 0\\ C\sum\limits_{k = 1}^N {{e^{A{x_k}}} - } \sum\limits_{k = 1}^N {{y_k}{e^{A{x_k}}}} = 0 \end{array} \right. E(A,C)=k=1N(CeAxkyk)2letAE=2k=1N(CeAxkyk)(CxkeAxk)=0CE=2k=1N(CeAxkyk)(eAxk)soCk=1Nxke2Axkk=1NykxkeAxk=0Ck=1NeAxkk=1NykeAxk=0  其中未知数 A , C A,C AC是非线性的,可用牛顿法求解,可以看出,直接利用最小二乘法求解更加复杂。除此之外,在求解其他函数的拟合曲线的时候,也可以进行类似的数据线性变换。具体如下图所示。

      除此之外,还可以使用线性最小二乘法来进行拟合。即设有N个数据点 { ( x k , y k ) } \{(x_k,y_k)\} {(xk,yk)},并给定了M个线性独立的函数 f i ( x ) {f_i(x)} fi(x)。利用这些线性独立的函数的线性组合表示函数 f ( x ) f(x) f(x),即
    f ( x ) = ∑ j = 1 M c j f j ( x ) f(x) = \sum\limits_{j = 1}^M {{c_j}{f_j}(x)} f(x)=j=1Mcjfj(x)  同样,利用最小二乘法求解系数的思想,使其数据点到曲线的垂直距离平方和最小,即其表达式对于每一个参数 c i c_i ci的偏导数为0,则可以得到其正规方程:
    ∑ j = 1 M ( ∑ k = 1 N f i ( x k ) f j ( x k ) ) c j = ∑ k = 1 N f i ( x k ) y k , i = 1 , 2 , . . . , M \sum\limits_{j = 1}^M {\left( {\sum\limits_{k = 1}^N {{f_i}({x_k}){f_j}({x_k})} } \right)} {c_j} = \sum\limits_{k = 1}^N {{f_i}({x_k})} {y_k},i = 1,2,...,M j=1M(k=1Nfi(xk)fj(xk))cj=k=1Nfi(xk)yk,i=1,2,...,M  不难看出,等式左边可以表达为矩阵乘积的结果,因此构造如下矩阵:

      则等式右边可以写成如下形式

      而方程左边未知参数c的系数可以写为矩阵F与其逆矩阵的乘积,因此,方程组可以表达为如下线性方程组的形式:
    F ′ F C = F ′ Y F'FC = F'Y FFC=FY  对于不太大的M,可以利用第三章求解线性方程组的方法求解系数c,从而得出函数的表达式。因此,当所取线性独立函数 f i ( x ) {f_i(x)} fi(x)为幂函数时,则函数 f ( x ) f(x) f(x)可以表示为M阶多项式的形式,然而事实上,这种多项式拟合通常具有多项式摆动的特性,因此会造成拟合结果的失真,而且阶数越高越容易发生,因此通常情况下只会取到二阶多项式进行拟合,除非已知使用的多项式是真实的。因此,进行最小二乘抛物线拟合的函数表达式为
    y = f ( x ) = C + B x + A x 2 y = f(x) = C + Bx + A{x^2} y=f(x)=C+Bx+Ax2  同样,利用最小二乘法求解系数的思想,使其数据点到曲线的垂直距离平方和最小,即其表达式对于每一个参数的偏导数为0,则可以得到其正规方程:
    ( ∑ k = 1 N x k 4 ) A + ( ∑ k = 1 N x k 3 ) B + ( ∑ k = 1 N x k 2 ) C = ∑ k = 1 N y k x k 2 ( ∑ k = 1 N x k 3 ) A + ( ∑ k = 1 N x k 2 ) B + ( ∑ k = 1 N x k ) C = ∑ k = 1 N y k x k ( ∑ k = 1 N x k 2 ) A + ( ∑ k = 1 N x k ) B + N C = ∑ k = 1 N y k \left( {\sum\limits_{k = 1}^N {x_k^4} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^3} } \right)B + \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)C = \sum\limits_{k = 1}^N {{y_k}x_k^2} \\ \left( {\sum\limits_{k = 1}^N {x_k^3} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)B + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)C = \sum\limits_{k = 1}^N {{y_k}x_k^{}} \\ \left( {\sum\limits_{k = 1}^N {x_k^2} } \right)A + \left( {\sum\limits_{k = 1}^N {x_k^{}} } \right)B + NC = \sum\limits_{k = 1}^N {{y_k}} (k=1Nxk4)A+(k=1Nxk3)B+(k=1Nxk2)C=k=1Nykxk2(k=1Nxk3)A+(k=1Nxk2)B+(k=1Nxk)C=k=1Nykxk(k=1Nxk2)A+(k=1Nxk)B+NC=k=1Nyk  以上系数表达式即为进行二阶多项式拟合的各阶项系数表达式。

    3、样条插值

      对N+1个点 的多项式插值经常不令人满意。一个N阶多项式可能有N-1个相对极大值和极小值,同时曲线可能会摆动,以经过这些点。另一个方法是将图形分段,每段为一个低阶多项式 S ( x ) S(x) S(x),并在相邻点 ( x k , y k ) (x_k,y_k) (xk,yk) ( x k + 1 , y k + 1 ) (x_{k+1},y_{k+1}) (xk+1,yk+1)之间进行插值.这样可以保证插值多项式没有高阶多项式插值的摆动性,而且这种插值方法一般是N段的分段函数。这样的分段拟合的函数称之为样条函数。其中一个重要的性质就是插值多项式经过结点。那么,接下来要思考的问题就是每一小段的插值多项式应该选择什么样的多项式。
      首先考虑一阶多项式,即一阶拉格朗日多项式。然而不难想到这种拟合方式在大多数情况下下是不适用的,因为拟合后的曲线为一段一段的折线段,显然不符合拟合的要求。然后再来考虑二阶多项式,二阶多项式能满足平滑的特点,但是会发现其在偶数结点处其曲率变化很大,这很有可能导致有非期望的弯曲或畸变。而且二次样条的导数在其结点处不连续。那么如果使用三次样条,则可以使其一阶二阶导数都连续,然而这个条件仅仅需要其多项式系数满足一定要求即可。不再继续考虑更高阶的样条函数是因为三次样条已经能满足需求,不必要牺牲产生震荡失真的代价去考虑高次样条,而且高次样条的系数增多,使求解变得极其困难。那么三次样条的多项式系数该如何确定呢?
      首先,我们需要给出三次样条需要满足的条件使其是一个光滑连续的函数,即一阶二阶导数在结点处连续。

      其中,性质Ⅰ描述了由分段三次多项式构成的 S ( x ) S(x) S(x)。性质 Ⅱ Ⅱ 描述了对给定数据点集的分段三次插值。性质 Ⅲ Ⅲ 和性质 Ⅳ Ⅳ 保证了分段三次多项式函数是一个光滑连续函数。性质 V V V保证了函数的二阶导数也是连续的。
      是否可能构造一个三次样条满足性质 Ⅰ Ⅰ 到性质 V V V?每个三次多项式 S ( x ) S(x) S(x)有4个未知常数,因此需要求解 4 N 4N 4N个系数。宽松的讲,要确定 4 N 4N 4N个自由度或条件。数据点提供了 N + 1 N+1 N+1个条件,性质 Ⅲ Ⅲ 、性质 Ⅳ Ⅳ 和性质 V V V都提供了 N − 1 N-1 N1个条件,总共确定了 N + 1 + 3 ( N − 1 ) = 4 N − 2 N+1+3(N - 1)=4N-2 N+1+3(N1)=4N2个条件。这样剩下了两个自由度。.可称之为端点约束( end-point constraints)涉及在点 x 0 x0 x0 x N xN xN处的导数 S ′ ( x ) S'(x) S(x) S " ( x ) S"(x) S"(x)。通过推导,由于二阶导数满足拉格朗日多项式,因此令其满足性质 V V V,再通过积分两次,可以求解出 S ( x ) S(x) S(x)的表达式,再利用性质 I I II II求出积分后不定系数的表达式,因此其表达式为
    S k ( x ) = − m k 6 h k ( x k + 1 − x ) 3 + m k + 1 6 h k ( x − x k ) 3 + ( y k h k − m k h k 6 ) ( x k + 1 − x ) + ( y k + 1 h k − m k + 1 h k 6 ) ( x − x k ) {S_k}(x) = - \frac{{{m_k}}}{{6{h_k}}}{\left( {{x_{k + 1}} - x} \right)^3} + \frac{{{m_{k + 1}}}}{{6{h_k}}}{\left( {x - {x_k}} \right)^3} + \left( {\frac{{{y_k}}}{{{h_k}}} - \frac{{{m_k}{h_k}}}{6}} \right)\left( {{x_{k + 1}} - x} \right) + \left( {\frac{{{y_{k + 1}}}}{{{h_k}}} - \frac{{{m_{k + 1}}{h_k}}}{6}} \right)\left( {x - {x_k}} \right) Sk(x)=6hkmk(xk+1x)3+6hkmk+1(xxk)3+(hkyk6mkhk)(xk+1x)+(hkyk+16mk+1hk)(xxk)   其中系数 h k h_k hk为第 k k k个和第 k + 1 k+1 k+1个结点的差分值, m k = S ′ ′ ( x k ) {m_k} = S''({x_k}) mk=S(xk)。那么只剩下一阶导数的性质没有用上。通过推到,利用一阶导数的性质,可以确定 m k m_k mk的递推关系式,即
    h k − 1 m k − 1 + 2 ( h k − 1 + h k ) m k + h k m k + 1 = u k , u k = 6 ( d k − d k − 1 ) , k = 1 , 2 , . . . , N − 1 {h_{k - 1}}{m_{k - 1}} + 2\left( {{h_{k - 1}} + {h_k}} \right){m_k} + {h_k}{m_{k + 1}} = {u_k},{u_k} = 6({d_k} - {d_{k - 1}}),k = 1,2,...,N - 1 hk1mk1+2(hk1+hk)mk+hkmk+1=uk,uk=6(dkdk1),k=1,2,...,N1  可以发现,可以构造一个有包含 N + 1 N+1 N+1个未知数具有 N − 1 N-1 N1个方程的线性方程组,即端点处的二阶导数值,因此有以下五种策略求解线性方程组。

      因此采用策略可以定出 m 0 m_0 m0的值,则第一个方程可以化简为
    2 ( h 0 + h 1 ) m 1 + h 1 m 2 = u 1 − h 0 m 0 , w h e n k = 1 2\left( {{h_0} + {h_1}} \right){m_1} + {h_1}{m_2} = {u_1} - {h_0}{m_0},{\rm{ }}when{\rm{ }}k = 1 2(h0+h1)m1+h1m2=u1h0m0,whenk=1  同理,已知 m N m_N mN的值也可以知道最后一个方程
    h N − 2 m N − 2 + 2 ( h N − 2 + h N − 1 ) m N − 1 = u N − 1 − h N − 1 m N , w h e n k = N − 1 {h_{N - 2}}{m_{N - 2}} + 2\left( {{h_{N - 2}} + {h_{N - 1}}} \right){m_{N - 1}} = {u_{N - 1}} - {h_{N - 1}}{m_N},{\rm{ }}when{\rm{ }}k = N - 1 hN2mN2+2(hN2+hN1)mN1=uN1hN1mN,whenk=N1  这样,剩下的系数就可以用一个带状的线性方程组来表示,即

      表示为 H M = V HM=V HM=V,利用第三章的几种求解线性方程组的方法,即可确定所有的结点二阶导数值,将这些值带回原插值函数,得到最终的插值函数多项式。到此,三次样条函数插值已经全部构造完成。
      对于五种不同的策略,其求解的 有所不同。具体为:

    1. 紧压样条。存在惟一的三次样条曲线,其一阶导数的边界条件是 S ′ ( a ) = d 0 S'(a)=d_0 S(a)=d0 S ′ ( b ) = d N S'(b)=d_N S(b)=dN,于是求解的方程组为

      Mark:紧压样条在端点有斜率。紧压样条可想像为用外力使柔软而有弹性的木杆经过数据点,并在端点处使其具有固定斜率。这样的样条对于画经过多个点的光滑曲线的绘图员相当有用。
    2. natural样条。存在惟一的三次样条曲线,其一阶导数的边界条件是 =0和 =0,于是求解的方程组为

      Mark:natural样条是柔软有弹性的木杆经过所有数据点形成的曲线,但让端点的斜率自由地在某一位置保持平衡,使得曲线的摇摆最小。它在对有多位有效数字精度的试验数据进行曲线拟合时很有用。
    3. 外推样条。存在惟一的三次样条曲线,其中通过对点 x 1 x_1 x1 x 2 x_2 x2进行外推得到 ,同时通过对点 x N − 1 x_{N-1} xN1 x N − 2 x_{N-2} xN2进行外推得到 。于是求解的方程组为

      Mark:外推样条等价于端点处的三次多项式曲线是相邻三次多项式曲线的扩展形成的样条,也就是说,样条曲线在区间 [ x 0 , x 2 ] [x_0,x_2] [x0,x2]内形成单个三次多项式曲线,同时在区间 [ x N − 2 , x N ] [x_{N-2},x_N] [xN2,xN]内形成另一个三次多项式曲线。
    4. 抛物线终结样条。存在唯一的三次样条,其中在区间 [ x 0 , x 1 ] [x_0,x_1] [x0,x1]内 ,而在 [ x N − 1 , x N ] [x_{N-1},x_N] [xN1,xN]内, .于是求解的方程组为

      Mark:在区间 [ x 0 , x 1 ] [x_0,x_1] [x0,x1] S " ( x ) = 0 S"(x)=0 S"(x)=0使得在区间 [ x 0 , x 1 ] [x_0,x_1] [x0,x1]内三次多项式曲线退化为二次多项式曲线,同时在区间 [ x N − 1 , x N ] [x_{N-1},x_N] [xN1,xN]内也发生同样的情况。
    5. 端点曲率调整样条。存在唯一的三次样条曲线,其中二阶导数的边界条件 S ’ ’ ( a ) S’’(a) S(a) S ’ ’ ( b ) S’’(b) S(b)使确定的。于是求解的方程组为

      Mark:直接对 S ’ ’ ( a ) S’’(a) S(a) S ’ ’ ( b ) S’’(b) S(b)赋值,可调整每个端点的曲率。

      综上所述,无论是哪种策略,其都给出了在端点处的一阶或二阶导数取值以及相应的线性方程组,以此来得出除端点外的所有结点二阶导数值,即 m k m_k mk。然后再通过五种策略的根据其他点二阶导数值来确定端点的二阶导数值,从而确定所有系数,最终求出三次样条插值函数。

    4、贝塞尔曲线

      20世纪70年代,雷诺汽车公司的Pierre Bezier和雪铁龙汽车公司的Paul de Castaliau各自独立推导出了CAD/CAM应用中的贝塞尔(Bezier)曲线。这些参数多项式是一类逼近样条。贝塞尔曲线是计算机图形学(CAD/CAM,计算机辅助几何设计)中曲线和曲面表示的基本方法。最初推导的贝塞尔曲线是用递归方法隐式定义的。将它们用伯恩斯坦多项式(Bemnstein polynomial)显式地定义,有助于推导贝塞尔曲线的性质。
      其中,N阶伯恩斯坦多项式的定义为
    B i , N ( t ) = ( N i ) t i ( 1 − t ) N − i i = 0 , 1 , 2 , . . . , N ( N i ) = N ! t ! ( N − i ) ! {B_{i,N}}(t) = \left( \begin{array}{l} N\\ i \end{array} \right){t^i}{(1 - t)^{N - i}}{\rm{ }}i = 0,1,2,...,N{\rm{ }}\left( \begin{array}{l} N\\ i \end{array} \right) = \frac{{N!}}{{t!(N - i)!}} Bi,N(t)=(Ni)ti(1t)Nii=0,1,2,...,N(Ni)=t!(Ni)!N!  一般,N阶伯恩斯坦多项式由N+1项。其中,伯恩斯坦多项式具有几个重要的性质。首先,伯恩斯坦多项式具有递归关系,表现为
    B 0 , 0 ( t ) = 1 , B i , N ( t ) = 0 , i < 0 o r i > N B i , N ( t ) = ( 1 − t ) B i , N − 1 ( t ) + t B i − 1 , N − 1 ( t ) = 0 , i = 1 , 2 , 3 , . . . , N − 1 {B_{0,0}}(t) = 1,{B_{i,N}}(t) = 0,{\rm{ }}i < 0{\rm{ }}or{\rm{ }}i > N\\ {B_{i,N}}(t) = \left( {1 - t} \right){B_{i,N - 1}}(t) + t{B_{i - 1,N - 1}}(t) = 0,{\rm{ }}i = 1,2,3,...,N - 1 B0,0(t)=1,Bi,N(t)=0,i<0ori>NBi,N(t)=(1t)Bi,N1(t)+tBi1,N1(t)=0,i=1,2,3,...,N1  因此,任意一个伯恩斯坦多项式可由其递归关系生成,并为后来的贝塞尔函数求解提供手段。其次,伯恩斯坦多项式还具有非负性,即在区间 [ 0 , 1 ] [0,1] [0,1]上是非负函数。因为伯恩斯坦多项式可以看作n次重复实验的伯努利分布的概率公式,其中t可以视为概率p,即只能在区间 [ 0 , 1 ] [0,1] [0,1]内取值,因此,其值也就是概率一定是处在 [ 0 , 1 ] [0,1] [0,1]内的非负数。第三个性质,伯恩斯坦多项式还具有规范性,即
    ∑ i = 0 N B i , N ( t ) = 1 \sum\limits_{i = 0}^N {{B_{i,N}}(t)} = 1 i=0NBi,N(t)=1  其规范性可以由二项式定理证明。第四个性质,也就是伯恩斯坦多项式的导数性质,即
    d d t B i , N ( t ) = N ( B i − 1 , N − 1 ( t ) − B i , N − 1 ( t ) ) \frac{d}{{dt}}{B_{i,N}}(t) = N\left( {{B_{i - 1,N - 1}}(t) - {B_{i,N - 1}}(t)} \right) dtdBi,N(t)=N(Bi1,N1(t)Bi,N1(t))  可以对伯恩斯坦多项式的定义是直接求导,即可证明其导数性质。最后一个性质,也是其可以推导出贝塞尔曲线的最重要的一个性质,就是其可以作为基,即N阶伯恩斯坦多项式组成阶数小于等于N的所有多项式的一个基空间。性质5说明,任何阶数小于等于N的多项式都可以惟一地表示为伯恩斯坦多项式的线性组合。
      给出贝塞尔曲线的定义,给定一个控制点集 ,一条N阶贝塞尔曲线定义为N阶伯恩斯坦多项式的加权和。即N阶贝塞尔曲线可以表示成
    P ( t ) = ∑ i = 0 N P i B i , N ( t ) , i = 0 , 1 , . . . , N , t ∈ [ 0 , 1 ] P(t) = \sum\limits_{i = 0}^N {{P_i}{B_{i,N}}(t)} ,{\rm{ }}i = 0,1,...,N,t \in [0,1] P(t)=i=0NPiBi,N(t),i=0,1,...,N,t[0,1]  其中 为N阶伯恩斯坦多项式。控制点是表示平面中x和y坐标的有序对。可将控制点作为向量,而对应的伯恩斯坦多项式作为标量处理,这样公式(28)可参数化表示为 P ( t ) = ( x ( t ) , y ( t ) ) P(t)=(x(t),y(t)) P(t)=(x(t),y(t)).其中
    x ( t ) = ∑ i = 0 N x i B i , N ( t ) , y ( t ) = ∑ i = 0 N y i B i , N ( t ) , x(t) = \sum\limits_{i = 0}^N {{x_i}{B_{i,N}}(t)} ,{\rm{ }}y(t) = \sum\limits_{i = 0}^N {{y_i}{B_{i,N}}(t)} ,{\rm{ }} x(t)=i=0NxiBi,N(t),y(t)=i=0NyiBi,N(t),  其中 0 < t < 1 0<t<1 0<t<1。函数 P ( t ) P(t) P(t)称为向量值函数,或等价地说,函数的值域是 x y xy xy平面上的一个点集。贝塞尔曲线因此也具有一些重要的性质,即 P 0 , P 1 P_0,P_1 P0P1在曲线的端点上,即
    P ( 0 ) = ∑ i = 0 N P i B i , N ( 0 ) = P 0 , P ( 1 ) = ∑ i = 0 N P i B i , N ( 1 ) = P N P(0) = \sum\limits_{i = 0}^N {{P_i}{B_{i,N}}(0)} = {P_0},P(1) = \sum\limits_{i = 0}^N {{P_i}{B_{i,N}}(1)} = {P_N} P(0)=i=0NPiBi,N(0)=P0,P(1)=i=0NPiBi,N(1)=PN  而其他点却不一定在曲线上,事实上大多数的点都不在曲线上。其次,贝塞尔曲线在 [ 0 , 1 ] [0,1] [0,1]上具有连续性,且具有任意阶导数。因此推导后有
    P ′ ( 0 ) = ∑ i = 0 N P i N ( B i − 1 , N − 1 ( 0 ) − B i , N − 1 ( 0 ) ) = N ( P 1 − P 0 ) P ′ ( 1 ) = N ( P N − P N − 1 ) P'(0) = \sum\limits_{i = 0}^N {{P_i}N\left( {{B_{i - 1,N - 1}}(0) - {B_{i,N - 1}}(0)} \right)} = N({P_1} - {P_0})\\ P'(1) = N({P_N} - {P_{N - 1}}) P(0)=i=0NPiN(Bi1,N1(0)Bi,N1(0))=N(P1P0)P(1)=N(PNPN1)  最后一个性质基于凸集的概念。即凸集的定义为,如果 x y xy xy平面上的子集C中任意两点连线上的所有点都是集合C中的元素,则称C为凸集。并且有如下定义,集合C的凸包是包含C的所有凸集的交。点的一个凸组合必定是该点集的凸包的子集。由伯恩斯坦多项式性质知,贝塞尔曲线是控制点的凸组合,因此它必定落在控制点的凸包内部。即有性质4贝塞尔曲线落在它的控制点集的凸包内。该性质说明,N阶贝塞尔曲线是一条连续曲线,并为其控制点集 的凸包所限,该曲线分别在点 P 0 P_0 P0 P 1 P_1 P1处开始和结束。贝塞尔观察到,曲线依次被"拉向"其余的合个不同。
      贝塞尔曲线的有效性在于,通过对控制点进行微小的调整,可以方便地调整(通过鼠标、键盘或其他图形界面)曲线的形状。改变任意一个控制点的坐标,将使整条曲线在参数区间 0 ≤ t ≤ 1 0≤t≤1 0t1的形状发生变化。然而这一变化在一定程度上是局部的,因为对应于控制点P的伯恩斯坦多项式在参数值 t = k / N t= k/N t=k/N处有最大值。因此,该贝塞尔曲线形状的最大变化应该在点 P ( k / N ) P(k/N) P(k/N)附近。这样,创建一条特定形状的曲线只需对初始控制点集做相对较少的改变。这样一条重要的特性在计算机辅助设计绘图中具有重要的作用,相较于样条插值,其没有高阶项,因其具有规范性,因此在变换时其可以有效地避免失真。

    5、三角多项式原理

       通常在工业生产中,会遇到一些具有周期性的信号。那么对于一些具有周期性的函数,怎样的拟合方式能够使得拟合更加完善呢?经过思考发现,一个具有周期性的函数,能够表示为由 a j c o s ( j x ) a_jcos(jx) ajcos(jx) b j s i n ( j x ) b_jsin(jx) bjsin(jx)线性组合的多项式表示。因此下面给出如下定义:如果存在值 t 0 , t 1 , … , t x t_0, t_1,…, t_x t0,t1,,tx满足 a = t 0 < t 1 < … < t x = b a = t_0 < t_1<…<t_x= b a=t0<t1<<tx=b,函数 f ( x ) f(x) f(x)在每个开区间 t i − 1 < x < t i , i = 1 , 2 , … , K t_{i-1}<x<t_i,i= 1,2,…,K ti1<x<ti,i=1,2,,K内是连续的,而且函数在每个端点 t i t_i ti有左极限和右极限,则称函数 f ( x ) f(x) f(x)在区间 [ a , b ] [a,b] [a,b]内分段连续。那么对于在一定区间内是分段连续且具有周期性的函数,可以将其展开为傅里叶级数的形式。设 f ( x ) f(x) f(x)是周期函数,周期为 2 Π 2Π 2Π而且 f ( x ) f(x) f(x)在区间 [ − Π , Π ] [-Π,Π] [Π,Π]内分段连续,则 f ( x ) f(x) f(x)的傅里叶级数 S ( x ) S(x) S(x)表示为
    S ( x ) = a 0 2 + ∑ j = 1 ∞ ( a j cos ⁡ ( j x ) + b j sin ⁡ ( j x ) ) S(x) = \frac{{{a_0}}}{2} + \sum\limits_{j = 1}^\infty {({a_j}} \cos (jx) + {b_j}\sin (jx)) S(x)=2a0+j=1(ajcos(jx)+bjsin(jx))   其中系数可以由欧拉公式得到,也可由三角函数系的正交性及其模得出,因此其中系数表达式为

       其中,设 S ( x ) S(x) S(x) f ( x ) f(x) f(x)在区间 [ − Π , Π ] [-Π,Π] [Π,Π]内的傅里叶级数。如果 f ( x ) f(x) f(x)在区间 [ − Π , Π ] [-Π,Π] [Π,Π]内是分段连续的,而且在区间内的每个端点有左导数和右导数,则 S ( x ) S(x) S(x)对于所有 x ∈ [ − Π , Π ] x∈[-Π,Π] x[Π,Π]是收敛的。对于所有的点 x ∈ [ − Π , Π ] x∈[-Π,Π] x[Π,Π],存在关系式 S ( x ) = f ( x ) S(x)=f(x) S(x)=f(x)其中 f ( x ) f(x) f(x)是连续的。如果 x = a x=a x=a是函数f的不连续点,则 ,这里分别表示左极限和右极限。这样,可得到傅里叶级数扩展表达式:
    f ( x ) = a 0 2 + ∑ j = 1 ∞ ( a j cos ⁡ ( j x ) + b j sin ⁡ ( j x ) ) f(x) = \frac{{{a_0}}}{2} + \sum\limits_{j = 1}^\infty {({a_j}} \cos (jx) + {b_j}\sin (jx)) f(x)=2a0+j=1(ajcos(jx)+bjsin(jx))   如果函数 f ( x ) f(x) f(x)为奇函数或者偶函数的时候,则由系数的确定公式可以知道,其中 a j a_j aj b j b_j bj会等于0,也就是说, f ( x ) f(x) f(x)的傅里叶级数仅含正弦函数或者余弦函数,通过公式可以确定此时的系数。最后,我们称形如
    T M ( x ) = a 0 2 + ∑ j = 1 M ( a j cos ⁡ ( j x ) + b j sin ⁡ ( j x ) ) {T_M}(x) = \frac{{{a_0}}}{2} + \sum\limits_{j = 1}^M {({a_j}} \cos (jx) + {b_j}\sin (jx)) TM(x)=2a0+j=1M(ajcos(jx)+bjsin(jx))   含有正余弦函数的多项式为三角多项式。因此,可以看见,三角多项式是傅里叶级数的近似解形式,当阶数越高时其误差越小,因此可以用来拟合函数。因为在实际工程中,当知道一个函数的某些点取值时,可以利用离散傅里叶级数进行三角多项式逼近,可以证明,三角多项式是使其均方误差最小的多项式。经过推导可知,其三角多项式的系数为;
    a j = 2 N ∑ k = 1 N f ( x k ) cos ⁡ ( j x k ) , j = 0 , 1 , . . . , M b j = 2 N ∑ k = 1 N f ( x k ) sin ⁡ ( j x k ) , j = 1 , 2 , . . . , M {a_j} = \frac{2}{N}\sum\limits_{k = 1}^N {f({x_k})\cos (j{x_k})} ,j = 0,1,...,M\\ {b_j} = \frac{2}{N}\sum\limits_{k = 1}^N {f({x_k})\sin (j{x_k})} ,j = 1,2,...,M aj=N2k=1Nf(xk)cos(jxk),j=0,1,...,Mbj=N2k=1Nf(xk)sin(jxk),j=1,2,...,M
       其中,需已知函数的 N + 1 N+1 N+1个点,并且横坐标之间等距且 x ∈ [ − Π , Π ] x∈[-Π,Π] x[Π,Π],则可展开为M阶三角多项式。尽管用最小二乘法定义,但可看成是连续傅里叶级数的数值近似值。欧拉公式给出了连续函数的傅里叶级数的系数,而给出了对数据点集进行曲线拟合的三角多项式系数。当使用更多的数据点时,三角多项式的系数更接近傅里叶级数的系数。

    二、实验内容及核心算法代码

    1、最小二乘拟合曲线原理实现

      由第一节原理部分可知,设 { ( x k , y k ) } k = 1 N \{ ({x_k},{y_k})\} _{k = 1}^N {(xk,yk)}k=1N是一个N个点的集合,其中横坐标是确定的。那么,理想的最小二乘拟合曲线 y = f ( x ) = A x + B y = f(x) = Ax + B y=f(x)=Ax+B的拟合结果应满足是均方根误差最小的曲线。那么,所选定的 A , B A,B AB构成的拟合曲线应当使 E 2 ( f ) E_2(f) E2(f)的值当且仅当所选取的 A , B A,B AB值处最小。也就是说,数据点到曲线的垂直距离平方和最小,其最小值为 。其系数 A , B A,B AB满足所示的正规方程。因此,仅需要知道点集合就可求出其拟合曲线的相关系数。
    其大致过程为:

    1. 获取所要拟合的点集合,并写出误差平方和表达式
    2. 对每个待定参数求偏导,并列出正规方程
    3. 求解线性方程组,得到最小二乘拟合曲线系数

      以下通过实例来进行最小线性二乘拟合的实验。
    (a) 根据数据 { ( x k , y k ) } k = 1 50 \{ ({x_k},{y_k})\} _{k = 1}^{50} {(xk,yk)}k=150,其中 x k = ( 0.1 ) k 且 y k = x k + c o s ( k 1 / 2 ) x_k=(0.1)k且y_k=x_k+cos(k_{1/2}) xk=(0.1)kyk=xk+cos(k1/2),求解 最小二乘线性拟合。
    (b) 计算 E 2 ( f ) E_2(f) E2(f).
    (c) 在同一坐标系下,画出点集和最小二乘线性拟合曲线。
      由题意可知,根据给出的 x , y x,y xy之间的关系可以生成由50个点所组成的集合。因此当进行最小线性二乘拟合的时候,可以列出其正规方程,并计算求解即得出其拟合曲线。然后,根据求出的拟合曲线,即可以计算其均方根误差,并绘制其曲线。
    其核心函数算法代码为:

    //the xpt,ypt are the point set coordinate,and the A,B are the coefficient of the curve
    //n is the size of the point set
    void CmpcurvebyLeastSquare(float* xpt, float* ypt, float& A, float& B, const int n, float** val = NULL)
    {
    	float xk2, xk, xkyk, yk;
    	float sumxk2, sumxk, sumxkyk, sumyk;
    	sumxk2 = 0, sumxk = 0, sumxkyk = 0, sumyk = 0;
    	//if user give the val argument,note down the step-value
    	for (int i = 0; i < n && val != NULL; ++i)
    	{
    		xk = xpt[i];
    		yk = ypt[i];
    		xk2 = xk * xk;
    		xkyk = xk * yk;
    		val[i][0] = xk;
    		val[i][1] = yk;
    		val[i][2] = xk2;
    		val[i][3] = xkyk;
    		sumxk += xk;
    		sumyk += yk;
    		sumxk2 += xk2;
    		sumxkyk += xkyk;
    	}
    	if (val != NULL)
    	{
    		val[n][0] = sumxk;
    		val[n][1] = sumyk;
    		val[n][2] = sumxk2;
    		val[n][3] = sumxkyk;
    	}
    	float** mat;
    	float* vB, * X;
    	mat = MallocArr(2, 2);
    	vB = MallocArr(2);
    	X = MallocArr(2);
    	mat[0][0] = sumxk2; mat[0][1] = sumxk; mat[1][0] = sumxk; mat[1][1] = n;
    	vB[0] = sumxkyk; vB[1] = sumyk;
    	//as solve the AB linear system,use gauss elem to solve out the AB
    	GaussElimmethodtoLinearsystem(mat, vB, X, 2);
    	A = X[0]; B = X[1];
    	FreeArr(mat, 2); FreeArr(vB); FreeArr(X);
    }
    

    2、曲线拟合原理实现

      由第一节原理部分可知,对于不是线性关系的曲线拟合,可以对其进行线性化操作再利用最小二乘法求出其系数的表达式即可。也就是说,对于一组不具有线性关系的数据点集合,可以对其进行诸如取对数、指数的方法将其变换为具有线性关系的表达式,即将数据点 x k , y k {xk,yk} xk,yk变换成 X Y XY XY平面上具有线性关系的点集 X k , Y k {Xk,Yk} Xk,Yk,再利用最小二乘法求出变换后的曲线表达式,最后再变换为原变量的表达式。
    因此,无论对于具有什么样关系的数据点集合,只要能够大致确定其函数类型,就可以通过相应的变换与换元进行最小二乘线性拟合,从而刻画出其函数关系。
    其大致过程如下;

    1. 确定数据点集合的函数关系,写出其通解与待定系数。
    2. 将通解进行如取指对数、三角换元等方法转化成具有线性关系的变量,并求解新变量的最小二乘拟合曲线的系数与表达式。
    3. 将求出系数的表达式变换回原变量的关系式,得到最终曲线拟合的表达式。

    以下通过实例来进行曲线拟合的实验。
    洛杉矶郊区在11月8日的温度记录如下表所示。共有24个数据点。
    (a) 使用fmins命令,对给定的数据集求解最小二乘曲线 f ( x ) = A c o s ( B x ) + C s i n ( D x ) + E f(x)=Acos(Bx)+Csin(Dx)+E f(x)=Acos(Bx)+Csin(Dx)+E
    (b) 求 E 2 ( f ) E_2(f) E2(f)
    (c) 在同一坐标系下绘制点集和(a)得出的最小二乘曲线。

    时间,p.m.温度时间,a.m.温度
    166158
    266258
    365358
    464458
    563557
    663657
    762757
    861858
    960960
    10601064
    11591167
    午夜58正午68

    由题意可知,使用matlab中的fminsearch命令求解最小化 E ( A , B , C , D , E ) E(A,B,C,D,E) E(ABCDE)后的A和C的近似值。首先在matlab中定义 E ( A , B , C , D , E ) E(A,B,C,D,E) E(A,B,C,D,E)为一个M文件,在另外一个源文件中使用fminsearch命令和初始值全设置为1.0,可以得到非线性最小二乘拟合的表达式的系数。
      其核心函数算法代码为:

    function z=E(u)
    A=u(1);
    B=u(2);
    C=u(3);
    D=u(4);
    E=u(5);
    z=(A.*cos(B.*1)+C.*sin(D.*1)+E-58).^2+(A.*cos(B.*2)+C.*sin(D.*2)+E-58).^2+...
    (A.*cos(B.*3)+C.*sin(D.*3)+E-58).^2+(A.*cos(B.*4)+C.*sin(D.*4)+E-58).^2+... (A.*cos(B.*5)+C.*sin(D.*5)+E-57).^2+(A.*cos(B.*6)+C.*sin(D.*6)+E-57).^2+... (A.*cos(B.*7)+C.*sin(D.*7)+E-57).^2+(A.*cos(B.*8)+C.*sin(D.*8)+E-58).^2+... (A.*cos(B.*9)+C.*sin(D.*9)+E-60).^2+(A.*cos(B.*10)+C.*sin(D.*10)+E-64).^2+... (A.*cos(B.*11)+C.*sin(D.*11)+E-67).^2+(A.*cos(B.*12)+C.*sin(D.*12)+E-68).^2+... (A.*cos(B.*13)+C.*sin(D.*13)+E-66).^2+(A.*cos(B.*14)+C.*sin(D.*14)+E-66).^2+... (A.*cos(B.*15)+C.*sin(D.*15)+E-65).^2+(A.*cos(B.*16)+C.*sin(D.*16)+E-64).^2+... (A.*cos(B.*17)+C.*sin(D.*17)+E-63).^2+(A.*cos(B.*18)+C.*sin(D.*18)+E-63).^2+... (A.*cos(B.*19)+C.*sin(D.*19)+E-62).^2+(A.*cos(B.*20)+C.*sin(D.*20)+E-61).^2+... (A.*cos(B.*21)+C.*sin(D.*21)+E-60).^2+(A.*cos(B.*22)+C.*sin(D.*22)+E-60).^2+... (A.*cos(B.*23)+C.*sin(D.*23)+E-59).^2+(A.*cos(B.*24)+C.*sin(D.*24)+E-58).^2;
    

    3、样条函数插值原理实现

      有第一节原理部分可知,将图形分段,每段为一个低阶多项式 S ( x ) S(x) S(x),并在相邻点 ( x k , y k ) (x_k,y_k) (xk,yk) ( x k + 1 , y k + 1 ) (x_{k+1},y_{k+1}) (xk+1,yk+1)之间进行插值.这样可以保证插值多项式没有高阶多项式插值的摆动性,而且这种插值方法一般是N段的分段函数。这样的分段拟合的函数称之为样条函数。使用三次样条,则可以使其一阶二阶导数都连续,然而这个条件仅仅需要其多项式系数满足一定要求即可。不再继续考虑更高阶的样条函数是因为三次样条已经能满足需求,不必要牺牲产生震荡失真的代价去考虑高次样条,而且高次样条的系数增多,使求解变得极其困难。
      因此,求解一个三次样条插值函数最重要的点就在于如何选取策略以及求解端点二阶导数值,其中求解除端点处二阶导数值的步骤大同小异,因此不需要放在重点考虑。
      因此,下面通过一个实例来展现natural样条的逼近效果。

    • 下面的表给出了在洛杉矶的郊区12小时内每个小时的温度(华氏温度)。根据这些数据求natural三次样条插值。在同–坐标系中,画出natural三次样条插值和这些数据。根据natural三次样条插值和习题12的(a)部分的结论求12小时内的平均温度近似值
    时间,a.m.温度时间,a.m.温度
    158757
    258858
    358960
    4581064
    5571167
    6571268

      根据前述步骤,其基本流程为:

    1. 根据数据算出已知的基本的参数,包括 h k , d k , u k h_k,d_k,u_k hk,dk,uk等。
    2. 根据所选择的策略列出线性方程组,并求解除端点外的二阶导数值。
    3. 根据求出的系数以及选定的策略,求解端点处的二阶导数值。
    4. 将系数带回原多项式,得到最终的三次样条插值多项式。

      其核心函数算法代码为:

    //cmp the hk coe,so must pass the point info,and save the data in the Hk,so have N+1-points,N-Hk;k=0,1,2...N-1
    //the formula is hk = x(k+1) - x(k),the N>=1
    void CmpHkcoeofSplFunc(float* xpt, const int N, float* Hk)
    {
    	if (N < 1)
    	{
    		cout << "arguments error!" << endl;
    		return;
    	}
    	for (int i = 0; i < N; ++i)
    	{
    		Hk[i] = xpt[i + 1] - xpt[i];
    	}
    }
    //cmp the dk coe,so must pass the point info,and save the data in the Hk,so have N+1-points,N-Dk;k=0,1,2..N-1 
    //the formula is dk = ( y(k+1)-y(k) ) / hk
    void CmpDkcoeofSplFunc(float* xpt, float* ypt, const int N, float* Dk)
    {
    	float* Hk;
    	float dy;
    	Hk = MallocArr(N);
    	CmpHkcoeofSplFunc(xpt, N, Hk);
    	for (int i = 0; i < N; ++i)
    	{
    		dy = ypt[i + 1] - ypt[i];
    		Dk[i] = dy / Hk[i];
    	}
    	FreeArr(Hk);
    }
    //cmp the uk coe,so must pass the point info,and save the data in the Hk,so have N+1-points,N-1-Uk;k=1,2,3...N-1
    // so the Uk only exist in the knot except the endpoints
    //the formula is uk = 6 * ( d(k) - d(k-1) )
    void CmpUkcoeofSplFunc(float* xpt, float* ypt, const int N, float* Uk)
    {
    	float* Dk;
    	Dk = MallocArr(N);
    	CmpDkcoeofSplFunc(xpt, ypt, N, Dk);
    	for (int i = 1; i < N; ++i)
    	{
    		Uk[i-1] = 6 * (Dk[i] - Dk[i - 1]);
    	}
    	FreeArr(Dk);
    }
    //as have the Mk,so cmp the Ski,which have N Hk,Dk,Ski,N+1 Mk
    void CmpSkiwithMk(float *ypt, float* Mk, float* Hk, float* Dk, float** Ski, const int N)
    {
    	float val1, val2;
    	for (int i = 0; i < N; ++i)
    	{
    		val1 = 2 * Mk[i] + Mk[i + 1];
    		val2 = Mk[i + 1] - Mk[i];
    		Ski[i][0] = ypt[i];
    		Ski[i][1] = Dk[i] - Hk[i] * val1 / 6;
    		Ski[i][2] = Mk[i] / 2;
    		Ski[i][3] = val2 / (6 * Hk[i]);
    	}
    }
    //Spline function(样条函数)
    //the natural spline function
    //pass the poitn coordinate as the arguments and the num of the spline function,equally how many points
    //save the coefficient in the Sk,so solve out the polynomial
    //the point coordinate xpt,ypt,and have N+1 point and Ski,so have N spline function
    //so should cmp the N<1,
    void CmpNaturalSplFunSki(float* xpt, float* ypt, const int N, float** Ski)
    {
    	//if N<=1,the spline function is 1-order Lagrange polynimial or does not exist
    	if (N <= 1)
    	{
    		cout << "the num of spline function is an error!" << endl;
    		return;
    	}
    	//construct the Mk
    	float* Mk, * Hk, * Uk, * Dk;
    	Mk = MallocArr(N + 1);
    	Hk = MallocArr(N);
    	Uk = MallocArr(N - 1);
    	Dk = MallocArr(N);
    	CmpHkcoeofSplFunc(xpt, N, Hk);
    	CmpUkcoeofSplFunc(xpt, ypt, N, Uk);
    	CmpDkcoeofSplFunc(xpt, ypt, N, Dk);
    	//the N=2
    	if (N == 2)
    	{
    		Mk[0] = 0; Mk[N] = 0;
    		Mk[1] = Uk[0] / (2 * (Hk[0] - Hk[1]));
    		CmpSkiwithMk(ypt, Mk, Hk, Dk, Ski, N);
    		return;
    	}
    	//the N=3
    	if (N == 3)
    	{
    		Mk[0] = 0; Mk[N] = 0;
    		float** H;
    		H = MallocArr(N - 1, N - 1);
    		H[0][0] = 2 * (Hk[0] + Hk[1]); H[0][1] = Hk[1];
    		H[1][0] = Hk[1]; H[1][1] = 2 * (Hk[1] + Hk[2]);
    		float* M;
    		M = MallocArr(N - 1);
    		GaussElimmethodtoLinearsystem(H, Uk, M, N - 1);
    		Mk[1] = M[0]; Mk[2] = M[1];
    		CmpSkiwithMk(ypt, Mk, Hk, Dk, Ski, N);
    		PrintMatrix(Ski,N,4);
    		return;
    	}
    	//the N>=4,so have M matrix
    	//construct the H matrix
    	float** H;
    	H = MallocArr(N - 1, N - 1);
    	H[0][0] = 2 * (Hk[0] - Hk[1]); H[0][1] = Hk[1];
    	H[N - 2][N - 3] = Hk[N - 2]; H[N - 2][N - 2] = 2 * (Hk[N - 2] + Hk[N - 1]);
    	for (int i = 1; i < N - 2; ++i)
    	{
    		H[i][i - 1] = Hk[i];
    		H[i][i] = 2 * (Hk[i] + Hk[i + 1]);
    		H[i][i + 1] = Hk[i + 1];
    	}
    	//solve Hk with the strategy natural spline function
    	Mk[0] = 0; Mk[N] = 0;
    	float* M;
    	M = MallocArr(N - 1);
    	GaussElimmethodtoLinearsystem(H, Uk, M, N - 1);
    	for (int i = 1; i <= N - 1; ++i)
    	{
    		Mk[i] = M[i - 1];
    	}
    	//solve the Ski coefficient
    	CmpSkiwithMk(ypt, Mk, Hk, Dk, Ski, N);
    	//free
    	FreeArr(Hk);
    	FreeArr(Dk);
    	FreeArr(Uk);
    	FreeArr(Mk);
    	FreeArr(H, N - 1);
    	FreeArr(M);
    }
    

    4、三角多项式原理实现

      由第一节原理部分可知,一个具有周期性的函数,能够表示为由 a j c o s ( j x ) a_jcos(jx) ajcos(jx) b j s i n ( j x ) b_jsin(jx) bjsin(jx)线性组合的多项式表示。因此下面给出如下定义:如果存在值 t 0 , t 1 , … , t x t_0, t_1,…, t_x t0,t1,,tx满足 a = t 0 < t 1 < … < t x = b a = t_0 < t_1<…<t_x= b a=t0<t1<<tx=b,函数 f ( x ) f(x) f(x)在每个开区间 t i − 1 < x < t i , i = 1 , 2 , … , K t_{i-1}<x<t_i,i= 1,2,…,K ti1<x<ti,i=1,2,,K内是连续的,而且函数在每个端点ti有左极限和右极限,则称函数 f ( x ) f(x) f(x)在区间 [ a , b ] [a,b] [a,b]内分段连续。那么对于在一定区间内是分段连续且具有周期性的函数,可以将其展开为傅里叶级数的形式。
      因此,对于一个具有周期性的函数,只需要对其一段周期上进行离散傅里叶级数展开,即可得到其三角多项式的逼近公式。在利用其周期性,即可求出其他区间上的点的值。
      通过以下实例可以展现如何求解并利用三角多项式拟合函数。
      与前述实例相同。
    洛杉矶郊区在11月8日的温度记录如下表所示。共有24个数据点。
    (a)求三角多项式 T 7 ( x ) T_7(x) T7(x).
    (b)在同一坐标系下绘制点集和(a)得出的三角多项式曲线。

    时间,p.m.温度时间,a.m.温度
    166158
    266258
    365358
    464458
    563557
    663657
    762757
    861858
    960960
    10601064
    11591167
    午夜58正午68

      根据已知条件可知,已知时间为1-24的各个温度值,当利用三角多项式进行拟合时,首先需要将其横坐标转换为 [ − Π , Π ] [-Π,Π] [Π,Π]上的等距节点,才能求解出其三角多项式,因此与最小二乘拟合不同,需要一个提前对横坐标进行预处理。
    利用已知点求解三角多项式的基本步骤为:

    1. 得到N+1个已知节点,并对其进行坐标变换 ,使其横坐标为位于区间 [ − Π , Π ] [-Π,Π] [Π,Π]等距已知节点。
    2. 利用系数公式,根据所需要求解多项式的阶数,求解其三角多项式。

      其核心函数算法代码为

    //the xk,yk is the data point set,the size of it is N,and the M is the order-number of the trigonometric polynomial
    //store the result in the aj,bj
    void CmpTrigopolyCoe(float* xk, float* yk, const int N, const int M, float* aj, float* bj)
    {
    	//cmp the aj,bj
    	aj[0] = 0;
    	for (int i = 0; i < N; ++i)
    	{
    		aj[0] += yk[i];
    	}
    	//cmp the first ele
    	aj[0] = 2 * aj[0] / float(N);
    	bj[0] = NULL;
    	for (int i = 1; i <= M; ++i)
    	{
    		float sumcos = 0, sumsin = 0;
    		//cmp the least-square coe
    		for (int j = 0; j < N; ++j)
    		{
    			sumcos += yk[j] * cos(i * xk[j]);
    			sumsin += yk[j] * sin(i * xk[j]);
    		}
    		//cmp the i-th aj,bj
    		aj[i] = 2 * sumcos / float(N);
    		bj[i] = 2 * sumsin / float(N);
    	}
    }
    //cmp the value of the trigonometric polynomial,the n is the order of the poly
    float CmpTrigonoPolyValue(float* aj, float* bj, const int n, float x)
    {
    	float y = aj[0] / 2;
    	float cosv, sinv;
    	for (int i = 1; i <= n; ++i)
    	{
    		cosv = aj[i] * cos(i * x);
    		sinv = bj[i] * sin(i * x);
    		y += (cosv + sinv);
    	}
    	return y;
    }
    

    5、贝塞尔曲线原理实现

      由第一节原理部分可知,贝塞尔曲线是计算机图形学(CAD/CAM,计算机辅助几何设计)中曲线和曲面表示的基本方法。推导的贝塞尔曲线是用递归方法隐式定义的。贝塞尔曲线是控制点的凸组合,因此它必定落在控制点的凸包内部。贝塞尔曲线的有效性在于,通过对控制点进行微小的调整,可以方便地调整(通过鼠标、键盘或其他图形界面)曲线的形状。改变任意一个控制点的坐标,将使整条曲线在参数区间 0 ≤ t ≤ 1 0≤t≤1 0t1的形状发生变化。
      因此,可以通过设计控制点的位置来设计贝塞尔曲线的样式,以此设计出更加复杂与大小位置可控的图形。由于贝塞尔曲线的定义及其性质使得其具有非常良好的不失真的特性,即其始终位于凸包范围内,从而便于计算机的几何设计。
      因此,通过以下实例来绘制并体会贝塞尔曲线的特性。

    • 编写程序,生成并绘制组合贝塞尔曲线。利用该程序生成和绘制过3个控制点集{(0,0),(1,2),(1,1),(3,0)},{(3,0),(4,-1),(5,-2),(6,1),(7,0)}和{(7,0),(4,- 3),(2,- 1),(0,0)}的贝塞尔曲线。

      由贝塞尔曲线的显式定义可知,既然给出了贝塞尔曲线的控制点集,就可以确定其阶数N。然后,利用伯恩斯坦多项式的递归的特性可以求出 t t t [ 0 , 1 ] [0,1] [0,1]范围内的函数表达式,将控制点集作为向量处理,由此可以得到不同t对应的点的向量,从而求出其贝塞尔曲线的 x , y x,y xy方向上的分量,从而绘制出控制点集下的N阶贝塞尔曲线。
      其基本流程为:

    1. 根据控制点集确定贝塞尔曲线的阶数N。
    2. 利用伯恩斯坦多项式的递归特性,将控制点集作为向量处理,伯恩斯坦多项式作为标量处理,求出其在 x , y x,y xy上的分量的表达式。
    3. t t t [ 0 , 1 ] [0,1] [0,1]范围内绘制贝塞尔曲线。

      其核心函数算法代码为:

    //recursion to cmp the value of the Bernstein polynomial value
    //the polynomial is B(i,N)(t)
    float BernsteinPolyValue(const int i, const int N, float t)
    {
    	if (i == 0)
    	{
    		return powf(1 - t, float(N));
    	}
    	if (i == N)
    	{
    		return powf(t, float(N));
    	}
    	return (1 - t) * BernsteinPolyValue(i, N - 1, t) + t * BernsteinPolyValue(i - 1, N - 1, t);
    }
    //cmp the Bezier curve value in each point given by tk
    //it's N-order Bezier curve,so have N+1 set of control points,and samply n point t in [0,1]
    //xk and yk are the set of control points,Xt and Yt are the value of each tk
    void CmpBezierCurveValueVec(float* xk, float* yk, float* t, float* Xt, float* Yt, const int N, const int n)
    {
    	float xval, yval, tk, Berval;
    	for (int i = 0; i < n; ++i)
    	{
    		xval = 0, yval = 0;
    		tk = t[i];
    		for (int j = 0; j <= N; ++j)
    		{
    			Berval = BernsteinPolyValue(j, N, tk);
    			xval += xk[j] * Berval;
    			yval += yk[j] * Berval;
    		}
    		Xt[i] = xval;
    		Yt[i] = yval;
    	}
    }
    

    三、结果分析

    1、最小二乘拟合曲线实例分析

      由实验数据结果可以看出,其线性拟合效果是非常好的,其中题目中所给的 x , y x,y xy关系式中含有余弦项,即可视为 x , y x,y xy线性关系的波动,因此进行线性拟合的时候其理性误差应当控制在1以内,而结果可以看出,其效果相当好。通过图2可以直观地看出其拟合效果是特别好的,因此最小二乘线性拟合法对于具有线性关系的数据具有良好的拟合效果。

    2、曲线拟合实例分析

      由图(3)可以直观看出,对于没有规律的数据点集合,进行非线性最小二乘曲线拟合后的结果的曲线走向大致与数据点相同,因此可以判定其拟合效果优良。除此之外,还可以发现数据点均匀的落在拟合曲线两侧,说明拟合曲线可以很好的刻画其变化规律且其数值对于原数据没有较大偏差,因此具有一定的参考价值。通过表(3)也可以看出,拟合曲线的误差在0附近摆动,且相对误差较小,可以更加确定判定拟合曲线对数据具有良好的拟合效果。综上所述,通过finsearch命令求解的拟合曲线具有良好的拟合特性。

    3、样条函数插值实例分析

      由表(4)可以看出三次样条插值拟合的取值情况,而从图(3)可以直接看出,利用三次样条的natural样条策略能够使得曲线拟合变得平滑,并且可以直观看出,在结点处满足样条插值的性质I和性质II,因此,三次样条具有良好的插值特性,可以推测使用其他样条策略同样能够具有如此的精度,并且相较于拉格朗日多项式插值而言,其拟合特性以及平滑度提高了很多,是一种非常不错的拟合方法。

    4、三角多项式实例分析

      由表(4)可以看出,在求解出来的三角多项式中,变换后的坐标对应的函数值如表所示,通过绘制曲线可以直观看出,将原数据点也进行坐标变换,使其位于区间[-Π,Π]上,可以看出,拟合曲线在端点处出现较大偏差,而在中间区域具有较好的拟合效果。其原因为由于三角多项式的特点,在端点处会有较大波动,因为三角函数的极大值与零点出现在端点处,因此会出现端点拟合失真的情况,但是从中间部分数据可以看出,其拟合效果以及曲线走势与原数据点相同,因此三角多项式具有良好的拟合特性,可以推测当多项式的阶数增加时,其拟合效果会更加良好。

    5、贝塞尔曲线实例分析

      由表(5)可以看出在不同控制点集下贝塞尔曲线的取值情况,可以看出,在每个控制点集的端点处,都位于贝塞尔曲线上,因此可以验证贝塞尔曲线的一条性质。由图(4)可以看出贝塞尔曲线的形状。可以看出,在控制点集发生明显突变的情况下,贝塞尔曲线并没有出现明显失真,曲线走向与控制点集的走向相同,且在各个分段的凸包内,因此贝塞尔曲线具有良好的图形表示性质。


    四、结论

      综上所述,通过分析各个曲线拟合方法的原理,可以发现各个拟合方法的优点。其中,线性或曲线最小二乘拟合方法适用于寻找数据点的规律,即通过最小二乘法来求取推测函数模型的系数,但是并不要求其拟合曲线一定要过各个数据点,因此适合于寻找一般数学物理规律的情况下。其次,样条插值拟合以及三角多项式拟合适用于希望无限逼近于数据点的情况下,即刻画数据点的曲线模型,直观的看出数据点的数学规律,并且三次样条插值具有非常好的曲线特征。最后,贝塞尔曲线适用于设计曲线形状,通过给定控制点集来求解位于凸包内的不失真的函数图像,以此来通过计算机辅助作图。因此,可以通过判定不同的场景,来选取最适当的曲线拟合方式。


    声明:本系列《数值计算方法》所有章节源代码均为作者编写,大家可以自行下载,仅供学习交流使用,欢迎大家评论留言或者私信交流,请勿用于任何商业行为。其中,各章实验数据结果也在资源链接中,大家可以自行实验对照,欢迎大家批评指正!文件中每个算法均由作者进行初步测试,部分算法可直接用于工程开发作为模块参考。

    各章算法源代码(C++)

    展开全文
  • 数字化测图详解

    千次阅读 2021-06-27 06:38:38
    数字测图系统主要由数据输入、数据处理和数据输出三部分组成流程如下: 地形图采集 数据处理与采集 成果与图形输出 数字化测图的基本思想: 数字化测图就是将采集的各种有关的地物和地貌信息转化为数字形式,通过数据...
  • 076_《OpenGL图形程序开发实务》

    千次阅读 2010-12-13 12:38:00
    主要介绍了计算机图形学的发展及开发现状、OpenGL的有关理论、OpenGL的特点及工作原理、OpenGL的基础知识及高级应用技巧,最后结合实例系统地讲述了如何综合利用OpenGL技术绘制复杂的三维真实感图形。本书深入浅出、...
  • 超过150种图形类型,包括隐式、等高图、复数、、向量场、密度、保角、常微分方程、偏微分方程、工程和统计图。  照明模式、标题、说明、光泽度、坐标轴控制、网格线、透明度等。  在图形标题、说明、刻度、标签...
  • CAD文件转换为SVG文件

    千次阅读 2021-06-09 15:28:33
    详细剖析了CAD的图形交换格式———DXF文件的结构,分析了SVG文件的框架元素和图形元素,建立了CAD中的对象、DXF文件中的实体和SVG中的元素三者之间的对应表,并对转换中遇到的难点问题提出了解决方案,最后,通过编程...
  • GDI+ 在 .NET 里有很多用途,包括向打印机输出文档、在一个 Windows 应用程序里显示图形、以及在网页里呈现图形。 你可以创建采用了用户指定信息的富图形,也可以基于数据库记录动态呈现图表或图形。 GDI+ 编程的...
  • ,把大的图形分割为更小的图形。 Dice ,把图形转换成微多边形网格,每个大概一个像素大小。 Shade ,计算每个微多边形网格顶点的灯光和颜色。 Bust ,把网格炸开成单个的微多边形,对每个计算边界并判断是否可见。 ...
  • matplotlib的原理或者说基础逻辑是,用Artist对象在画布(canvas)上绘制(Render)图形。这与人作画的步骤类似: 准备一块画布或画纸 准备好颜料、画笔等制图工具 作画 所以相对,matplotlib有三个层次的API: ...
  • 本控件提供直角坐标或极坐标模式,对任意数值与任意数学表达式组合,进行图形绘制。 C# (C# 2.0) Windows, .NET (.NET 2.0) Win32, VS, GDI+ * 下载演示项目...
  • Matplotlib之Python可视化

    千次阅读 2021-07-20 06:18:24
    学习绘制图形的API,熟悉各API的参数。 图形颜色和线条美化,选择适合所分析行业的颜色和线条,例如分析的行业是金融业就选择黑灰商务色,看起来严谨认真的线条和字体;分析的是教育行业就选择鲜活可爱的颜色主题和...
  • 离散点插值方法、等值线的绘制及平滑技巧吕勇平 戴景茹(广东省气候应用研究所 510080) 由于等值线图看起来非常直观、形象,因此在天气预报、气候预测分析等方面用得非常多,已成为预报员不可缺少的工具之一。...
  • 1不能应用修剪命令“trim”进行修剪的对象是(D、文字)。2.命令行(B.不能随意移动)3. 布尔运算中差集的热键为(A.SU)4.... 既可以绘制直线,又可以绘制曲线的命令是( C.多段线)9.在修改编辑中,只能采用交叉窗口或...
  • tilcon所有控件说明

    2014-07-05 11:15:23
    tilcon是VxWorks下强大的图形界面开发工具,tilcon控件说明告诉大家怎么使用tilcon里面的控件 非常不错
  • 音乐的频谱分析

    2014-06-01 18:03:45
    可以对以下情况进行帮助 1、利用MATLAB的声音频谱分析。 2、MATLAB降噪 3、MATLAB合成音乐并加谐波
  • 第十一章 AWT编程

    2022-04-03 20:48:58
    总体上,AWT是图形用户界面编程的基础,Swing组件替代了绝大部分AWT组件,对AWT图形用户界面编程有好的补充和加强。 Java9的AWT和Swing组件可以自适应高分辨率屏。在Java9之前,如果使用高分辨率屏,由于这种屏幕...
  • GDI+编程小结

    万次阅读 多人点赞 2010-10-28 21:13:00
    矢量图形包括坐标系统中的系列点指定的绘图基元(如直线、曲线图形)。例如,直线可通过它的两个端点来指定,而矩形可通过确定其左上角位置的点并给出其宽度和高度的一对数字来指定。简单路径可由通过直线连接的点...
  • MEDICI仿真NMOS器件晶体管语法笔记

    千次阅读 2019-09-29 01:16:04
    MEDICI 仿真 NMOS 器件晶体管   ... TMA MEDICI Example 1 - 1.5 Micron N-Channel MOSFET ...下面将要求解漏电压和漏电流的关系,因为是 NMOS 器件,所以设置载流子类型为电子 SOLVE  V(Drain)=...
  • WPF学习第十二集-绘图和动画

    千次阅读 2014-01-15 15:59:02
    Path:路径(闭合区域),基本图形中功能最强的一个,可由若干直线,圆弧,被塞尔曲线组成。 1 直线 直线是最简单的图形。使用X1,Y1两个属性值可以设置它的起点坐标,X2,Y2两个属性值可以设置它的终点坐标...
  • 数据分析总集

    2021-08-30 17:57:39
    数据分析 什么是数据分析? 数据分析是指用适当的统计分析方法对收集来的大量数据进行分析,提取有用信息和形成结论而对数据加以详细研究和概括总结的过程。 数据分析经典案例 (一)啤酒与尿布 ...
  • LaTeX语法集合 cmap颜色映射表 等高线图 热成像图 极坐标系 3D图像绘制 简单动画 numpy常用函数 加载文件 算数平均值 加权平均值 numpy常用函数 最值 中位数 标准差 案例应用 时间数据处理 数组的轴向汇总 移动均线 ...
  • 绘制二维直角坐标 plot(x,y,z) 绘制三维直角坐标 极坐标polar(x,y,‘-*’); 极坐标 [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IgBhUXRe-1651042249325)(C:\Users\31219\Documents\...
  • 自动控制理论

    2021-05-20 02:49:26
    掌握求线性定常连续系统输出响应的方法,运用连续系统时域响应函数(impulse,step,lsim),得到系统的时域响应曲线。 2.掌握使用MATLAB软件作出系统根轨迹;利用根轨迹图对控制系统进行分析;掌...