【CAD算法】【计算机图形学】Bezier贝塞尔曲线生成程序(python/numpy实现)[1]

整个项目，从Bezier曲线的创建，到Coons Patch曲面的实现，再到网格的实现和优化，还有最后对表面的光顺，链接如下：
【CAD算法】【计算机图形学】Bezier贝塞尔曲线生成程序(python/numpy实现)[1]
【CAD算法】【计算机图形学】Coons Patch曲面生成程序(python/numpy实现)[2]
【CAD核心算法】基于Laplacian Smoothing的网格mesh优化算法(python实现)[3]
【CAD算法】【遗传算法】采用遗传算法优化Coons Patch曲面程序(python实现)[4]
【CAD算法】基于最速下降法(梯度下降)的最短路径实现(python实现)[5]

1.内容回顾

首先回顾一下，表示一个曲线，有两种方式，implicit表示法和参数表示法。
implicit表示法：通过XYZ坐标表示曲线。
优点：

• Compact
• Easy to check if a point belongs to the curve

缺点：

• Difficult for curve evaluation
• Difficult for partial curve definition

参数表示法：Represent the X, Y, Z coordinates as a function of a single parameter
优点：

• Easy for curve evaluation.
• Convenient for partial curve definition.
• Many others such as easy for manipulation, intersection.

很明显，参数表示法在工程领域中被广泛运用。

参数表示法有6种，分别是：

• Conic sections
• Polynomial Curves
• Hermite Curves
• Bezier Curves
• B-Splines
• Non-Uniform Rational B-Splines (NURBS)

Conic sections：适用于描述圆锥曲线(圆/椭圆/双曲线/抛物线)。
Polynomial Curves：P(u) = [x(u) y(u) z(u)]T = a0 + a1u + a2u2 + a3u3 (0<=u<=1),这种表示法不能直观地通过改变a0到a3的值预测到曲线形状，一般不被运用。

Hermite Curves：P(u) = f0(u)P0 + f1(u)P1 + f2(u)P0’+ f3(u)P1’ (0<=u<=1)，通过确定一条曲线两端点的坐标和斜率来确定一条曲线的坐标。同样不能直观地通过两端点的信息确定曲线的形状。
B spine：也就是B样条曲线，这种表示法得到最广泛的运用，它克服了Bezier曲线某一个控制点坐标变换会整体上地影响曲线形状这个缺点。
Bezier Curves：通过n+1个控制点的坐标，得到曲线的坐标。它的定义如下：

一定要注意，参数u的取值范围只能在0-1之间，包括端点。根据公式可知，只要知道了bezier曲线的维度n，也就是控制点数减一，就能计算出每一项的blending function的值，分别和控制点坐标相乘再求和，就能得到Bezier曲线的表达式。一旦确定了一个u的值，自然能得到对应曲线的坐标。如果想通过pyplot画出来，可以通过采样的方式，用很多个点来描述这个曲线。

2.题目描述

现在开始我们的程序。问题如下，用户通过txt或者控制台输入一组坐标，该程序能自动计算出Bezier曲线并通过画图软件画出来。
其中，输入信息格式如下：

Line 1           m          /*  Bezier 曲线的数量 */
Line 2           n1          /* 第一条Bezier曲线有多少个控制点 */
Line 3           0 or 1      /*0意味着该曲线首尾闭合, 1代表着首尾不相连 */
Line 4           第一个控制点坐标  /* 坐标之间通过“,”或空格分开 */
Line 5           第二个控制点坐标
.
.
.
Line n1+3        第n个控制点坐标

可能这么描述还不是很清楚，举一个例子：

2                           #一共需要绘制2条Bezier曲线
4                           #第一条曲线有4个控制点
1                           #第一条曲线首尾不相连
1.5		0.4	-3              #控制点的三维坐标
1.7   	-2.56   	4.78    #控制点的三维坐标
0   		245.6 	10      #控制点的三维坐标
0		0	0               #控制点的三维坐标
3                           #第二条曲线有4个控制点
0                           #第二条曲线首尾相连
3		5	-6              #控制点的三维坐标
-4		3.5	7               #控制点的三维坐标
1.5		-5	-10.5           #控制点的三维坐标

3.解决方案

3.1头文件引用和变量定义

import os
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd

np.set_printoptions(suppress = True)
#database = []
m = np.array([],dtype = np.uint8) # the number of bazier curves
pointsNumber = np.array([],dtype = np.uint8) #for the ith curve, the points number of it
closedOrOpen = np.array([],dtype = np.uint8) #for the ith curve, 0 means it is open, 1 means it is closed
point = np.array([],dtype = np.float64) #the control points for all the curves

解释一下，以上四个变量均采用array数据类型，方便在内存中存放很多数据点并可视化。
其中，m是Bezier曲线数量；pointsNumber依次存放着每一条Bezier曲线的控制点数；closedOrOpen依次存放着每一条Bezier曲线是闭合的还是开放的bool数据；point依次存放着每一条Bezier曲线的全部特征点数，配合pointsNumber[i]指针可以很快定位到每一条Bezier曲线的第一个控制点在point数组中的位置。

3.2实现从txt文件的读取

def inputFromText():
  global filepath
  filepath = os.path.abspath('.') + "/" + raw_input("please input the file name(eg:data.txt):")
  
  with open (filepath) as data:
    line_zero = data.readline()
    line_zero_list = list(line_zero.split())
    global m
    m = np.array(line_zero_list, dtype=np.uint8)

    for i in range(m[0]): 
      pointsNumber_str = data.readline()
      pointsNumber_list = list(pointsNumber_str.split())
      temp1 = np.array(pointsNumber_list,dtype=np.uint8)
      global pointsNumber
      pointsNumber = np.append(pointsNumber,temp1) # the number of control points of the ith Bezier curve
      closedOrOpen_str = data.readline()
      closedOrOpen_list = list(closedOrOpen_str.split())
      temp2 = np.array(closedOrOpen_list,dtype=np.uint8)
      global closedOrOpen
      closedOrOpen = np.append(closedOrOpen,temp2) #0 means the ith curve is closed, 1 means it is open
      
      global point
      for j in range(temp1[0]):
      #u += 1
        point_str = data.readline()
        point_list = list(point_str.split())
        temp3 = np.array(point_list,dtype=np.float64)
        point = np.append(point,temp3) # control points for the ith Bezier Curve
  point = np.reshape(point, (-1,3))

读取txt采用python内置的文件操作函数。为了便于理解，采用循环方式依次读取txt的每一行，代码中没有介绍过的变量都是中间变量。

3.3 实现从控制台读取数据并保存数据到txt

def inputFromKeyboard(): 
  
  m_int = input("the number of Beziercurve:") #the number of Beziercurve
  global m
  m = np.array([m_int],dtype = np.uint8)

  filename = raw_input("please input filename:(example: data.txt):")  #input can only read 1 word!
  
  with open(filename,mode="w") as file_txt:
    file_txt.write(str(m_int)+'\n') #change to string!!

    for i in range(m_int): 
      pointsNumber_str = str(input("numbers of control points of the {}th Bezier curve:".format(i)))
      pointsNumber_list = list(pointsNumber_str.split())
      temp1 = np.array(pointsNumber_list,dtype=np.uint8)
      global pointsNumber
      pointsNumber = np.append(pointsNumber,temp1) # the number of control points of the ith Bezier curve
    #u += 1
      file_txt.write(pointsNumber_str + '\n')

      closedOrOpen_str = str(input("0 or 1: 0 means the ith curve is closed, 1 means it is open:")) #0 means the ith curve is closed, 1 means it is open
      closedOrOpen_list = list(closedOrOpen_str.split())
      temp2 = np.array(closedOrOpen_list,dtype=np.uint8)
      global closedOrOpen
      closedOrOpen = np.append(closedOrOpen,temp2) #0 means the ith curve is closed, 1 means it is open
      file_txt.write(closedOrOpen_str + '\n')

      global point
      for j in range(temp1[0]):
        point_str = str(raw_input("the coordinates of control points for the {}th Bezier Curve (eg: '0 0 0'):".format(i))) # control points for the ith Bezier Curve
        point_list = list(point_str.split())
        temp3 = np.array(point_list,dtype=np.float64)
        point = np.append(point,temp3) # control points for the ith Bezier Curve
        file_txt.write(point_str + '\n')

需要注意的是，保存到txt里的是string类型。

3.4实现Bezier曲线的核心程序

def drawCurve(closedOrOpen,pointsNumber,point,point_index):

  C = []
  n = pointsNumber - 1 # n is fewer in numbers than the total control points number. According to definition.
  point_show = np.array([],dtype=np.float64)
  for i in range(n+1):
    point_show = np.append(point_show,point[point_index + i])  

  if (closedOrOpen == 0): # if it is closed, means the oth and nth control points are the same.
    n += 1
    point_show = np.append(point_show,point[point_index])
  elif (closedOrOpen == 1):
    pass
  point_show = point_show.reshape((-1,3))
  #根据Bizier曲线定义
  if ((n+1) % 2 == 0):
    for i in range((n+1) / 2):
      up = 1
      down = 1
      j = n
      while (j > i):
        up *= j
        j = j - 1
      j = n - i
      while (j > 0):
        down *= j
        j = j - 1
      C.append(up / down)
  elif ((n+1) % 2 == 1):
    for i in range(n / 2):
      up = 1
      down = 1
      j = n
      while (j > i):
        up *= j
        j = j - 1
      j = n - i
      while (j > 0):
        down *= j
        j = j - 1
      C.append(up / down)
    up = 1
    down = 1
    j = n
    while (j > n/2):
      up *= j
      j = j - 1
    j = n/2
    while (j > 0):
      down *= j
      j = j - 1
    C.append(up / down)
  if (n%2 == 1):
    for i in range(int((n+1)/2)):
      C.append(C[int(n/2-i)])
  if (n%2 == 0):
    for i in range(int((n+1)/2)):
      C.append(C[int(n/2-i-1)])
  print("C",C) 
  f = np.array([],dtype=np.float64)

  u = 0
  fx = 0
  fy = 0 #not this place!!
  fz = 0
  #print(f_unit.shape)
  print("point_show\n",point_show)
  print("n+1",n+1)
  for j in range(100): 
    fx = 0
    fy = 0  #do not forget return 0!
    fz = 0
    for i in range(n+1):#这一条语句和Bezier曲线的表达式对应上了
      fx += C[i] * u**i * (1-u)**(n-i) * point_show[i][0]
      fy += C[i] * u**i * (1-u)**(n-i) * point_show[i][1]
      fz += C[i] * u**i * (1-u)**(n-i) * point_show[i][2]
    list = []
    list.append(fx)
    list.append(fy) 
    list.append(fz)
    array_list = np.array([list],dtype=np.float64) 
    u += 0.01
    f = np.append(f,array_list)
  f_show = f.reshape((-1,3))
  print(f_show.shape)
  #print(f_show)

  #show the control lines!
 ax.plot(f_show[:,0],f_show[:,1],f_show[:,2],'b.',markersize=1,label='open bazier curve')#bezier曲线，把'b.'改为'-'可以得到连续的曲线
  ax.plot(point_show[:,0],point_show[:,1],point_show[:,2],'-',markersize=1, label='control line')#控制点的连线
  ax.plot(point_show[:,0],point_show[:,1],point_show[:,2],'r.',markersize=8, label='control points')#控制点
  print(point_show.shape[0])

注意，当closedOrOpen为0时，也就是说这是一条封闭的Bezier，那么需要在末尾再把这个曲线的第一个控制点push进array里， 这样就能首尾相连。

3.5 main()函数

inputType = input("Please choose data from text/0 or keyboard/1 ?:")
#print(inputType.type())

if ( inputType == 0): # not '0'!!
  inputFromText()#从txt获取数据
elif ( inputType == 1):
  inputFromKeyboard()#从键盘获取数据

point_index = 0
fig = plt.figure()
ax = fig.gca(projection='3d')

for i in range(m[0]):#绘制m条Bezier曲线
  drawCurve(closedOrOpen[i],pointsNumber[i],point,point_index)
  point_index += pointsNumber[i]#获取每一条Bezier曲线第一个控制点在数组中的位置

ax.legend()
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

4.试验结果

这个结果是根据样例数据绘制，红色点为各个Bezier曲线的控制点，中间的点为Bezier曲线，如果在plt.plot()中增加"-"参数或者增加点数，j这条曲线就能更好看一些。

采用Ubuntu16.03进行调试，交互界面如上。下图为一个更加残暴的数据所绘制出来的图像。

所有代码的链接如下：https://github.com/iwander-all/CAD.git

参考文献
K. TANG, Fundamental Theories and Algorithms of CAD/CAE/CAM