精华内容
下载资源
问答
  • 高数立体图形

    2014-06-21 21:18:23
    高等数学微积分中常用的立体曲线,制作讲义幻灯式可用
  • 使用函数创建 # 如果生成一定规则的数据,可以使用NumPy提供的专门函数 # arange函数类似于python的range函数:指定起始值、终止值和步长来创建数组 # 和Python的range类似,arange同样不包括终值;但arange可以生成...

    !/usr/bin/python

    -- coding:utf-8 --

    import numpy as np

    import matplotlib as mpl

    from mpl_toolkits.mplot3d import Axes3D

    from matplotlib import cm

    import time

    from scipy.optimize import leastsq

    from scipy import stats

    import scipy.optimize as opt

    import scipy

    import matplotlib.pyplot as plt

    from scipy.stats import norm, poisson

    from scipy.interpolate import BarycentricInterpolator

    from scipy.interpolate import CubicSpline

    import math

    import seaborn

    def residual(t, x, y):

    return y - (t[0] * x ** 2 + t[1] * x + t[2])

    def residual2(t, x, y):

    print(t[0], t[1])

    return y - t[0]np.sin(t[1]x)

    x ** x x > 0

    (-x) ** (-x) x < 0

    def f(x):

    y = np.ones_like(x)

    i = x > 0

    y[i] = np.power(x[i], x[i])

    i = x < 0

    y[i] = np.power(-x[i], -x[i])

    return y

    if name == "main":

    # # 开场白:

    # numpy是非常好用的数据包,如:可以这样得到这个二维数组

    # [[ 0 1 2 3 4 5]

    # [10 11 12 13 14 15]

    # [20 21 22 23 24 25]

    # [30 31 32 33 34 35]

    # [40 41 42 43 44 45]

    # [50 51 52 53 54 55]]

    # a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)

    # print a

    # 正式开始 -:)

    # 标准Python的列表(list)中,元素本质是对象。

    # 如:L = [1, 2, 3],需要3个指针和三个整数对象,对于数值运算比较浪费内存和CPU。

    # 因此,Numpy提供了ndarray(N-dimensional array object)对象:存储单一数据类型的多维数组。

    # # 1.使用array创建

    # 通过array函数传递list对象

    # L = [1, 2, 3, 4, 5, 6]

    # print("L = ", L)

    # a = np.array(L)

    # print("a = ", a)

    # print(type(a))

    # # # 若传递的是多层嵌套的list,将创建多维数组

    # b = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]])

    # print(b)

    # # #

    # # # # 数组大小可以通过其shape属性获得

    # print(a.shape)

    # print(b.shape)

    # # #

    # # # 也可以强制修改shape

    # b.shape = 4, 3

    # print(b)

    # # # # 注:从(3,4)改为(4,3)并不是对数组进行转置,而只是改变每个轴的大小,数组元素在内存中的位置并没有改变

    # # #

    # # # 当某个轴为-1时,将根据数组元素的个数自动计算此轴的长度

    # b.shape = 2, -1

    # print(b)

    # print(b.shape)

    # # #

    # b.shape = 3, 4

    # # # 使用reshape方法,可以创建改变了尺寸的新数组,原数组的shape保持不变

    # c = b.reshape((4, -1))

    # print("b = \n", b

    # print('c = \n', c)

    # #

    # # # 数组b和c共享内存,修改任意一个将影响另外一个

    # b[0][1] = 20

    # print("b = \n", b)

    # print("c = \n", c)

    # # #

    # # # 数组的元素类型可以通过dtype属性获得

    # print(a.dtype)

    # print(b.dtype)

    # # #

    # # # # 可以通过dtype参数在创建时指定元素类型

    # d = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.float)

    # f = np.array([[1, 2, 3, 4], [5, 6, 7, 8], [9, 10, 11, 12]], dtype=np.complex)

    # print(d)

    # print(f)

    # # #

    # # # 如果更改元素类型,可以使用astype安全的转换

    # f = d.astype(np.int)

    # print(f)

    # #

    # # # 但不要强制仅修改元素类型,如下面这句,将会以int来解释单精度float类型

    # d.dtype = np.int

    # print(d)

    #

    # 2.使用函数创建

    # 如果生成一定规则的数据,可以使用NumPy提供的专门函数

    # arange函数类似于python的range函数:指定起始值、终止值和步长来创建数组

    # 和Python的range类似,arange同样不包括终值;但arange可以生成浮点类型,而range只能是整数类型

    # a = np.arange(1, 10, 0.5)

    # print(a)

    # # #

    # # # linspace函数通过指定起始值、终止值和元素个数来创建数组,缺省包括终止值

    # b = np.linspace(1, 10, 10)

    # print('b = ', b)

    # # #

    # # 可以通过endpoint关键字指定是否包括终值

    # c = np.linspace(1, 10, 10, endpoint=False)

    # print('c = ', c)

    # # #

    # # 和linspace类似,logspace可以创建等比数列

    # # 下面函数创建起始值为10^1,终止值为10^2,有20个数的等比数列

    # d = np.logspace(1, 2, 10, endpoint=True)

    # print(d)

    # # #

    # # # 下面创建起始值为2^0,终止值为2^10(包括),有10个数的等比数列

    # f = np.logspace(0, 10, 11, endpoint=True, base=2)

    # print(f)

    # # #

    # # # 使用 frombuffer, fromstring, fromfile等函数可以从字节序列创建数组

    # s = 'abcd'

    # g = np.fromstring(s, dtype=np.int8)

    # print(g)

    # #

    # 3.存取

    # 3.1常规办法:数组元素的存取方法和Python的标准方法相同

    # a = np.arange(10)

    # print(a)

    # # # 获取某个元素

    # print(a[3])

    # # # # 切片[3,6),左闭右开

    # print(a[3:6])

    # # # # 省略开始下标,表示从0开始

    # print(a[:5])

    # # # # 下标为负表示从后向前数

    # print(a[3:])

    # # # # 步长为2

    # print(a[1:9:2])

    # # # # 步长为-1,即翻转

    # print(a[::-1])

    # # # # 切片数据是原数组的一个视图,与原数组共享内容空间,可以直接修改元素值

    # a[1:4] = 10, 20, 30

    # print(a)

    # # # # 因此,在实践中,切实注意原始数据是否被破坏,如:

    # b = a[2:5]

    # b[0] = 200

    # print(a)

    #

    # # 3.2 整数/布尔数组存取

    # # 3.2.1

    # 根据整数数组存取:当使用整数序列对数组元素进行存取时,

    # 将使用整数序列中的每个元素作为下标,整数序列可以是列表(list)或者数组(ndarray)。

    # 使用整数序列作为下标获得的数组不和原始数组共享数据空间。

    # a = np.logspace(0, 9, 10, base=2)

    # print(a)

    # i = np.arange(0, 10, 2)

    # print i

    # # # # 利用i取a中的元素

    # b = a[i]

    # print b

    # # # b的元素更改,a中元素不受影响

    # b[2] = 1.6

    # print b

    # print a

    # # 3.2.2

    # 使用布尔数组i作为下标存取数组a中的元素:返回数组a中所有在数组b中对应下标为True的元素

    # # 生成10个满足[0,1)中均匀分布的随机数

    # a = np.random.rand(10)

    # print a

    # # # 大于0.5的元素索引

    # print a > 0.5

    # # # # 大于0.5的元素

    # b = a[a > 0.5]

    # print b

    # # # # 将原数组中大于0.5的元素截取成0.5

    # a[a > 0.5] = 0.5

    # print a

    # # # # b不受影响

    # print b

    # 3.3 二维数组的切片

    # [[ 0 1 2 3 4 5]

    # [10 11 12 13 14 15]

    # [20 21 22 23 24 25]

    # [30 31 32 33 34 35]

    # [40 41 42 43 44 45]

    # [50 51 52 53 54 55]]

    # a = np.arange(0, 60, 10) # 行向量

    # print 'a = ', a

    # b = a.reshape((-1, 1)) # 转换成列向量

    # print b

    # c = np.arange(6)

    # print c

    # f = b + c # 行 + 列

    # print f

    # # 合并上述代码:

    # a = np.arange(0, 60, 10).reshape((-1, 1)) + np.arange(6)

    # print a

    # # 二维数组的切片

    # print a[[0, 1, 2], [2 ,3, 4]]

    # print a[4, [2, 3, 4]]

    # print a[4:, [2, 3, 4]]

    # i = np.array([True, False, True, False, False, True])

    # print a[i]

    # print a[i, 3]

    # # 4.1 numpy与Python数学库的时间比较

    # for j in np.logspace(0, 7, 10):

    # j = int(j)

    # x = np.linspace(0, 10, j)

    # start = time.clock()

    # y = np.sin(x)

    # t1 = time.clock() - start

    #

    # x = x.tolist()

    # start = time.clock()

    # for i, t in enumerate(x):

    # x[i] = math.sin(t)

    # t2 = time.clock() - start

    # print j, ": ", t1, t2, t2/t1

    # 4.2 元素去重

    # 4.2.1直接使用库函数

    # a = np.array((1, 2, 3, 4, 5, 5, 7, 3, 2, 2, 8, 8))

    # print '原始数组:', a

    # # 使用库函数unique

    # b = np.unique(a)

    # print '去重后:', b

    # # 4.2.2 二维数组的去重,结果会是预期的么?

    # c = np.array(((1, 2), (3, 4), (5, 6), (1, 3), (3, 4), (7, 6)))

    # print '二维数组', c

    # print '去重后:', np.unique(c)

    # # 4.2.3 方案1:转换为虚数

    # # r, i = np.split(c, (1, ), axis=1)

    # # x = r + i * 1j

    # x = c[:, 0] + c[:, 1] * 1j

    # print '转换成虚数:', x

    # print '虚数去重后:', np.unique(x)

    # print np.unique(x, return_index=True) # 思考return_index的意义

    # idx = np.unique(x, return_index=True)[1]

    # print '二维数组去重:\n', c[idx]

    # # 4.2.3 方案2:利用set

    # print '去重方案2:\n', np.array(list(set([tuple(t) for t in c])))

    # # 4.3 stack and axis

    # a = np.arange(1, 10).reshape((3, 3))

    # b = np.arange(11, 20).reshape((3, 3))

    # c = np.arange(101, 110).reshape((3, 3))

    # print 'a = \n', a

    # print 'b = \n', b

    # print 'c = \n', c

    # print 'axis = 0 \n', np.stack((a, b, c), axis=0)

    # print 'axis = 1 \n', np.stack((a, b, c), axis=1)

    # print 'axis = 2 \n', np.stack((a, b, c), axis=2)

    # 5.绘图

    # 5.1 绘制正态分布概率密度函数

    # mu = 0

    # sigma = 1

    # x = np.linspace(mu - 3 * sigma, mu + 3 * sigma, 50)

    # y = np.exp(-(x - mu) ** 2 / (2 * sigma ** 2)) / (math.sqrt(2 * math.pi) * sigma)

    # print x.shape

    # print 'x = \n', x

    # print y.shape

    # print 'y = \n', y

    # # plt.plot(x, y, 'ro-', linewidth=2)

    # plt.plot(x, y, 'r-', x, y, 'go', linewidth=2, markersize=8)

    # plt.grid(True)

    # plt.show()

    mpl.rcParams['font.sans-serif'] = [u'Bold'] #FangSong/黑体 FangSong/KaiTi

    mpl.rcParams['axes.unicode_minus'] = False

    # # 5.2 损失函数:Logistic损失(-1,1)/SVM Hinge损失/ 0/1损失

    =============================================================================

    x = np.array(np.linspace(start=-2, stop=3, num=1001, dtype=np.float))

    y_logit = np.log(1 + np.exp(-x)) / math.log(2)

    y_boost = np.exp(-x)

    y_01 = x < 0

    y_hinge = 1.0 - x

    y_hinge[y_hinge < 0] = 0

    plt.plot(x, y_logit, 'r-', label='Logistic Loss', linewidth=2)

    plt.plot(x, y_01, 'g-', label='0/1 Loss', linewidth=2)

    plt.plot(x, y_hinge, 'b-', label='Hinge Loss', linewidth=2)

    plt.plot(x, y_boost, 'm--', label='Adaboost Loss', linewidth=2)

    plt.grid()

    plt.legend(loc='upper right')

    plt.title(u"SVM loss 函数图像")

    plt.savefig('1.png')

    plt.show()

    =============================================================================

    # # 5.3 x^x

    =============================================================================

    x = np.linspace(-1.3, 1.3, 101)

    y = f(x)

    plt.plot(x, y, 'g-', label='x^x', linewidth=2)

    plt.grid()

    plt.legend(loc='upper left')

    plt.show()

    =============================================================================

    # # 5.4 胸型线

    =============================================================================

    x = np.arange(1, 0, -0.001)

    y = (-3 * x * np.log(x) + np.exp(-(40 * (x - 1 / np.e)) ** 4) / 25) / 2

    plt.figure(figsize=(5,7))

    plt.plot(y, x, 'r-', linewidth=2)

    plt.grid(True)

    plt.show()

    =============================================================================

    # 5.5 心形线

    =============================================================================

    t = np.linspace(0, 7, 100)

    x = 16 * np.sin(t) ** 3

    y = 13 * np.cos(t) - 5 * np.cos(2t) - 2 * np.cos(3t) - np.cos(4*t)

    plt.plot(x, y, 'r-', linewidth=2)

    plt.grid(True)

    plt.show()

    =============================================================================

    # # 5.6 渐开线

    =============================================================================

    t = np.linspace(0, 50, num=1000)

    x = t*np.sin(t) + np.cos(t)

    y = np.sin(t) - t*np.cos(t)

    plt.plot(x, y, 'y+', linewidth=2)

    plt.grid()

    plt.show()

    =============================================================================

    # # Bar

    =============================================================================

    mpl.rcParams['font.sans-serif'] = [u'SimHei'] #黑体 FangSong/KaiTi

    mpl.rcParams['axes.unicode_minus'] = False

    x = np.arange(0, 10, 0.05)

    y = np.sin(x)

    plt.bar(x, y, width=0.04, linewidth=0.2)

    plt.plot(x, y, 'r--', linewidth=2)

    plt.title(u'Sin曲线')

    plt.xticks(rotation=-60)

    plt.xlabel('X')

    plt.ylabel('Y')

    plt.grid()

    plt.show()

    =============================================================================

    # # 6. 概率分布

    # # 6.1 均匀分布

    =============================================================================

    x = 100*np.random.rand(10000)

    t = np.arange(len(x))

    plt.hist(x, 30, color='m', alpha=0.5, label=u'均匀分布')

    plt.plot(t, x, 'r-', label=u'均匀分布')

    plt.legend(loc='upper left')

    plt.grid()

    plt.show()

    =============================================================================

    # # 6.2 验证中心极限定理

    =============================================================================

    t = 1000

    a = np.zeros(10000)

    for i in range(t):

    a += np.random.uniform(-5, 5, 10000)

    a /= t

    plt.hist(a, bins=30, color='g', alpha=0.5, normed=True, label=u'均匀分布叠加')

    plt.legend(loc='upper left')

    plt.grid()

    plt.show()

    =============================================================================

    # 6.21 其他分布的中心极限定理

    =============================================================================

    lamda = 10

    p = stats.poisson(lamda)

    y = p.rvs(size=1000)

    mx = 30

    r = (0, mx)

    bins = r[1] - r[0]

    plt.figure(figsize=(10, 8), facecolor='w')

    plt.subplot(121)

    plt.hist(y, bins=bins, range=r, color='g', alpha=0.8, normed=True)

    t = np.arange(0, mx+1)

    plt.plot(t, p.pmf(t), 'ro-', lw=2)

    plt.grid(True)

    =============================================================================

    #

    =============================================================================

    lamda=10

    N = 1000

    M = 10000

    plt.subplot(111)

    a = np.zeros(M, dtype=np.float)

    p = stats.poisson(lamda)

    for i in np.arange(N):

    y = p.rvs(size=M)

    a += y

    a /= N

    plt.hist(a, bins=20, color='g', alpha=0.8, normed=True)

    plt.grid(b=True)

    plt.show()

    =============================================================================

    # # 6.3 Poisson分布

    =============================================================================

    x = np.random.poisson(lam=5, size=10000)

    print(x)

    pillar = 15

    a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)

    plt.grid()

    plt.show()

    print(a)

    print(a[0].sum())

    =============================================================================

    # # 6.4 直方图的使用

    =============================================================================

    mu = 2

    sigma = 3

    data = mu + sigma * np.random.randn(1000)

    h = plt.hist(data, 30, normed=1, color='#a0a0ff')

    x = h[1]

    y = norm.pdf(x, loc=mu, scale=sigma)

    plt.plot(x, y, 'r--', x, y, 'ro', linewidth=2, markersize=4)

    plt.grid()

    plt.show()

    =============================================================================

    # # 6.5 插值

    =============================================================================

    x = np.random.poisson(lam=5, size=10000)

    print(x)

    pillar = 15

    a = plt.hist(x, bins=pillar, normed=True, range=[0, pillar], color='g', alpha=0.5)

    rv = poisson(5)

    x1 = a[1]

    y1 = rv.pmf(x1)

    itp = BarycentricInterpolator(x1, y1) # 重心插值

    x2 = np.linspace(x1.min(), x1.max(), 50)

    y2 = itp(x2)

    cs = scipy.interpolate.CubicSpline(x1, y1) # 三次样条插值

    plt.plot(x2, cs(x2), 'm--', linewidth=5, label='CubicSpine') # 三次样条插值

    plt.plot(x2, y2, 'g-', linewidth=3, label='BarycentricInterpolator') # 重心插值

    plt.plot(x1, y1, 'r-', linewidth=1, label='Actural Value') # 原始值

    plt.legend(loc='upper right')

    plt.grid()

    plt.show()

    =============================================================================

    # 7. 绘制三维图像

    =============================================================================

    x, y = np.ogrid[-3:3:100j, -3:3:100j]

    u = np.linspace(-3, 3, 101)

    x, y = np.meshgrid(u, u)

    z = xynp.exp(-(x2 + y2)/2) / math.sqrt(2*math.pi)

    z = xynp.exp(-(x2 + y2)/2) / math.sqrt(2*math.pi)

    fig = plt.figure()

    ax = fig.add_subplot(111, projection='3d')

    ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.coolwarm, linewidth=0.1)

    ax.plot_surface(x, y, z, rstride=5, cstride=5, cmap=cm.Accent, linewidth=0.5)

    plt.show()

    cmaps = [('Perceptually Uniform Sequential',

    ['viridis', 'inferno', 'plasma', 'magma']),

    ('Sequential', ['Blues', 'BuGn', 'BuPu',

    'GnBu', 'Greens', 'Greys', 'Oranges', 'OrRd',

    'PuBu', 'PuBuGn', 'PuRd', 'Purples', 'RdPu',

    'Reds', 'YlGn', 'YlGnBu', 'YlOrBr', 'YlOrRd']),

    ('Sequential (2)', ['afmhot', 'autumn', 'bone', 'cool',

    'copper', 'gist_heat', 'gray', 'hot',

    'pink', 'spring', 'summer', 'winter']),

    ('Diverging', ['BrBG', 'bwr', 'coolwarm', 'PiYG', 'PRGn', 'PuOr',

    'RdBu', 'RdGy', 'RdYlBu', 'RdYlGn', 'Spectral',

    'seismic']),

    ('Qualitative', ['Accent', 'Dark2', 'Paired', 'Pastel1',

    'Pastel2', 'Set1', 'Set2', 'Set3']),

    ('Miscellaneous', ['gist_earth', 'terrain', 'ocean', 'gist_stern',

    'brg', 'CMRmap', 'cubehelix',

    'gnuplot', 'gnuplot2', 'gist_ncar',

    'nipy_spectral', 'jet', 'rainbow',

    'gist_rainbow', 'hsv', 'flag', 'prism'])]

    =============================================================================

    # 8.1 scipy

    # 线性回归例1

    =============================================================================

    x = np.linspace(-2, 2, 50)

    A, B, C = 2, 3, -1

    y = (A * x ** 2 + B * x + C) + np.random.rand(len(x))*0.75

    t = leastsq(residual, [0, 0, 0], args=(x, y))

    theta = t[0]

    print('真实值:', A, B, C)

    print('预测值:', theta)

    y_hat = theta[0] * x ** 2 + theta[1] * x + theta[2]

    plt.plot(x, y, 'r-', linewidth=2, label=u'Actual')

    plt.plot(x, y_hat, 'g-', linewidth=2, label=u'Predict')

    plt.legend(loc='upper left')

    plt.grid()

    plt.show()

    =============================================================================

    # # 线性回归例2

    x = np.linspace(0, 5, 100)

    A = 5

    w = 1.5

    y = A * np.sin(w*x) + np.random.rand(len(x)) - 0.5

    t = leastsq(residual2, [3, 1], args=(x, y))

    theta = t[0]

    print('真实值:', A, w)

    print('预测值:', theta)

    y_hat = theta[0] * np.sin(theta[1] * x)

    plt.plot(x, y, 'r-', linewidth=2, label='Actual')

    plt.plot(x, y_hat, 'g-', linewidth=2, label='Predict')

    plt.legend(loc='lower left')

    plt.grid()

    plt.show()

    # # 8.2 使用scipy计算函数极值

    # a = opt.fmin(f, 1)

    # b = opt.fmin_cg(f, 1)

    # c = opt.fmin_bfgs(f, 1)

    # print a, 1/a, math.e

    # print b

    # print c

    # markerdescription

    # ”.”point

    # ”,”pixel

    # “o”circle

    # “v”triangle_down

    # “^”triangle_up

    # “<”triangle_left

    # “>”triangle_right

    # “1”tri_down

    # “2”tri_up

    # “3”tri_left

    # “4”tri_right

    # “8”octagon

    # “s”square

    # “p”pentagon

    # “*”star

    # “h”hexagon1

    # “H”hexagon2

    # “+”plus

    # “x”x

    # “D”diamond

    # “d”thin_diamond

    # “|”vline

    # “_”hline

    # TICKLEFTtickleft

    # TICKRIGHTtickright

    # TICKUPtickup

    # TICKDOWNtickdown

    # CARETLEFTcaretleft

    # CARETRIGHTcaretright

    # CARETUPcaretup

    # CARETDOWNcaretdown

    展开全文
  • cao congyong - Lecture Notes * 4 三维立体图形绘制 4.3 三维曲面图绘图函数 为了绘制定义在平面区域 D =[x0,xm][y0,yn ]上的三维曲面z=f(x,y) 首先将[x0,xm]在 x 方向分成 m 份将[y0,yn]在 y 方向分成 n 份由各划...
  • 绘制三维图形的基本函数

    万次阅读 2018-10-29 17:17:16
    1.绘制三维图形的基本函数   1 2 3 4 5 6 最基本的三维绘图函数为plot3; plot3与plot用法十分相似,调用格式:   plot(x1,y1,z1,选项1,x2,y2,z2,选项2,...,...

    1.绘制三维图形的基本函数

        

    1

    2

    3

    4

    5

    6

    最基本的三维绘图函数为plot3;

    plot3与plot用法十分相似,调用格式:

     

    plot(x1,y1,z1,选项1,x2,y2,z2,选项2,...,xn,yn,zn,选项n)

    当x,y,z是同维向量时,则x,y,z,对应元素构成一条三维曲线;

    当x,y,z是同维矩阵时,则以x,y,z对应列元素绘制三维曲线,曲线条数等于矩阵列数。

      例:

    程序如下:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    t=0:pi/50:2*pi;

    x=8*cos(t);

    y=4*sqrt(2)*sin(t);

    z=-4*sqrt(2)*sin(t);

    plot3(x,y,z,'p');

    title('Line in 3-D Space');

    text(0,0,0,'origin');

    xlabel('x'),ylabel('y'),zlabel('z');

    grid;

      运行结果:

    2.三维曲面

    2.1平面网格坐标矩阵的生成

      绘制z=f(x,y)所代表的三维曲面图,先要在xy平面选定一个矩形区域,假定矩形区域D=[a,b]*[c,d],然后将[a,b]在x方向分成m份,将[c,d]在y方向分成n份,由各划分点分别作平行于两坐标轴的直线,将区域D分成m*n个小矩形,生成代表每一个小矩形顶点坐标的平面网格坐标矩阵,最后利用有关函数绘图。

      产生平面区域内的网格坐标矩阵有两种方法:

      1.利用矩阵运算生成、

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    x=a:dx:b;

    y=(c:dy:d)';

    X=ones(size(y))*x;

    Y=y*ones(size(x));

    语句执行后,

    矩阵X的每一行都是向量x,行数等于向量y的元素个数,

    矩阵Y的每一列都是向量y,列数等于向量x的元素个数。

    于是对于矩阵X,Y来说,它们位置(i,j)上的元素值(X(i,j),Y(i,j))就是所要形成的平面网格

    在位置(i,j)上的X,Y坐标。可根据每一个网格点上的x,y坐标求这个点对应的z,则得到Z矩阵。

    显然,X,Y,Z各列或各行所对应坐标,对应于一条空间曲线,空间曲线的集合将可组成空间曲面。

      2.利用meshgrid函数生成。

     

    1

    2

    3

    4

    5

    6

    调用格式:

    x=a:dx:b;

    y=c:dy:d;

    [x,y]=meshgrid(x,y);

    语句执行后得到与方法1相同的矩阵X,Y。

    当向量x=y时,函数可写成meshgrid(x);

      例:利用法网格坐标阵巧解不定方程:

      已知6<x<30,15<y<36,求不定方程2x+5y=126的整数解。程序如下:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    25

    26

    27

    28

    29

    30

    31

    32

    33

    34

    35

    36

    37

    38

    39

    40

    41

    42

    43

    44

    45

    46

    47

    48

    49

    50

    x=7:29;

    y=16:35;

    [x,y]=meshgrid(x,y); %在[7,29]*[16,35]区域生成网格坐标

    z=2*x+5*y;

    k=find(z==126);%找出解的位置,即k为z中元素等于126的元素的位置

    x(k)',y(k)'%输出对应位置的x,y即方程的解

     

    输出:

    ans =

     

         8    13    18    23

     

     

    ans =

     

        22    20    18    16

    %即方程有4组解:(8,22),(13,20),(18,18)(23,16).

     

    输出:

    >> k

     

    k =

     

        27

       125

       223

       321

     

    输出(关于find函数):

    >> [a,b]=find(z==126)

     

    a =

     

         7

         5

         3

         1

     

     

    b =

     

         2

         7

        12

        17

    >> x(7,2)

     

    ans =

     

         8

    2.2 绘制三维曲面的函数

    两个函数:

    1

    2

    3

    4

    5

    6

    7

    8

    mesh(x,y,z,c)%用于绘制三维网格图,在不需要绘制特别精细三维曲面时使用。

    surf(x,y,z,c)%用于绘制三维曲面,各线条之间的补面用颜色填充。

     

    关于x,y,z,c:

     

    one:通常x,y,z是同维矩阵,x,y是网格坐标矩阵,z是网格点的高度矩阵,c用于指定在不同高度下的颜色范围。

    two:c省略时,MATLAB认为c=z,即颜色的高度正比于图形高度,以得到层次分明的三维图形。当x,y省略时,把z矩阵的列下标当做x轴坐标,把z矩阵的行下标当做y轴坐标,然后绘制三维曲面图。

    three:当x,y是向量时,要求x的长度必须等于z矩阵的列数,y的长度等于z矩阵的行数,x,y向量元素的组合构成网格点的x,y坐标,z坐标则取自z矩阵,然后绘制三维曲面图。

    例5.15:用三维曲面图表现函数z=sinycosx。

    program1:用meshgrid+mesh

    1

    2

    3

    4

    5

    6

    x=0:0.1:2*pi;

    [x,y]=meshgrid(x);

    z=sin(y).*cos(x);

    mesh(x,y,z);

    xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis');

    title('mesh');<br><br>效果同:

    x=0:0.1:2*pi;
    y=0:0.1:2*pi;
    z=sin(y')*cos(x);
    mesh(x,y,z);
    xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis');

      运行结果:

    program2:用meshgrid+surf

    1

    2

    3

    4

    5

    6

    x=0:0.1:2*pi;

    [x,y]=meshgrid(x);

    z=sin(y).*cos(x);

    surf(x,y,z);

    xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis');

    title('meshgrid+surf');

      

    program3:用一般绘图函数plot3

    1

    2

    3

    4

    5

    6

    7

    x=0:0.1:2*pi;

    [x,y]=meshgrid(x);

    z=sin(y).*cos(x);

    plot3(x,y,z);

    xlabel('x-axis'),ylabel('y-axis'),zlabel('z-axis');

    title('meshgrid+plot3-1f');

    grid;

      

    例5.16:绘制两个直径相等的圆管的相交图形。

    程序如下:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    21

    22

    23

    24

    25

    26

    27

    28

    %两个等直径圆管的交线

     

    m=60;%m是圆圈的密集度,表示画60个圆圈

    z=1.2*(0:m)/m;%1.2是圆柱高

    r=ones(size(z));

    theta=(0:m)/m*2*pi;

     

    x1=r'*cos(theta);%每行都是一个cos(theta)

    y1=r'*sin(theta);%每行都是一个sin(theta)

    %y1=y1';

    z1=z'*ones(1,m+1);%每行的z相同

     

    surf(x1,y1,z1);%绘图,立起的圆柱

     

    %axis equal,axis off

    hold on

     

    x=(-m:2:m)/m;

    x2=x'*ones(1,m+1);%m+1个x列

    y2=r'*cos(theta);%以y和z为底画圆

    %y2=y2';

    z2=r'*sin(theta);

     

    surf(x2,y2,z2);

     

    axis equal,axis off

    title('两个等直径圆管的交线');

    hold off

      运行结果:

     将上述例5.16中程序的%备注取消,即将第一图的y阵第二图的z阵转置,这样在底层面就不再是圆线了,效果如下:

    例5.17 分析由函数z=x^2-2y^2构成的曲面形状及与平面z=a的交线。

     

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    [x,y]=meshgrid(-10:0.2:10);

    z1=(x.^2-2*y.^2)+eps;%第一个曲面

    a=input('a=?');

    z2=a*ones(size(x));%第二个曲面(本质是一个数乘)

    subplot(1,2,1);

    mesh(x,y,z1);hold on;mesh(x,y,z2);%分别画出两个曲面

    v=[-10,10,-10,10,-100,100];axis(v);grid;%第一个子图的坐标设置

    hold off;

    r0=(abs(z1-z2)<=1);%求两曲面z坐标差小于1的点,r0只有0、1值

    xx=r0.*x;yy=r0.*y;zz=r0.*z2;%求这些点上的x,y,z坐标,即交线坐标

    subplot(1,2,2);

    plot3(xx(r0~=0),yy(r0~=0),zz(r0~=0),'*');%在第2子图画出交线

    axis(v);grid;%第2子图的坐标设置

    a=?8
    size(x)

    ans =

    101 101

    1

    <br>

      

      此外,还有两个和mesh函数相似的函数,即带等高线的三维网格曲面函数meshc带底座的三维网格曲面函数meshz。其用法与mesh类似,不同的是meshc还在xy平面上绘制曲面在z轴方向的等高线,meshz还在xy平面上绘制曲面的底座。

          例5.18 在xy平面内选择区域[-8,8]*[-8,8],绘制函数的4种三维曲面图。

    程序如下:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    [x,y]=meshgrid(-8:0.5:8);

    z=sin(sqrt(x.^2)+y.^2)./sqrt(x.^2+y.^2+eps);

    subplot(2,2,1);

    meshc(x,y,z);

    title('meshz(x,y,z)')

    subplot(2,2,2);

    meshz(x,y,z);

    title('meshz(x,y,z)')

    subplot(2,2,3);

    surfc(x,y,z)

    title('surfc(x,y,z)')

    subplot(2,2,4);

    surfl(x,y,z)

    title('surf1(x,y,z)')

      

      3.标准三维曲面

      MATLAB提供了一些函数用于绘制标准三维曲面,这些函数可以产生相应的绘图数据,常用于三维图形的演示。例如:

    sphere函数和cylinder函数分别用于绘制三维球面和柱面。其调用格式为:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    <strong>sphere函数的调用格式为:</strong>

    [x,y,z]=sphere(n)

    该函数将产生(n+1)*(n+1)矩阵x,y,z,采用这3个矩阵可以绘制出圆心位于原点、半径为1的单位球体。

    若在调用该函数时不带输出参数,则直接绘制所需球面。n决定了球面的圆滑程度,其默认值为20.若n值取得较小,则将绘制出多面体表面图。

     

     

     

     

    <strong>cylinder函数调用格式为:</strong>

    [x,y,z]=cylinder(R,n)

    其中,R是一个向量,存放柱面各个等间隔高度上的半径,n表示在圆柱圆周上有n个间隔点,默认有20个间隔点。例如,cylinder(3)生成一个圆柱,cylinder([10,1])生成一个圆锥,

    t=0:pi/100:4*pi;

    R=sin(t);

    cylinder(R,30)

    生成一个正弦柱面。

    另外,生成矩阵的大小与R向量的长度及n有关。其余与sphere函数相同。

      

     

    MATLAB还有一个peaks函数,称为多峰函数,,常用于三维曲面的演示。该函数可以用来生成绘图数据矩阵,矩阵元素由函数:

    在矩形区域[-3,3]*[-3,3]的等分网格点的函数值确定。例如:

    z=peaks(30);

    将生成一个30*30矩阵z,即分别沿x和y方向将区间[-3,3]等分成29份,并计算这些网格点上的函数值。默认的等分数是48,即p=peaks将生成一个49*49矩阵p。也可以根据风格坐标矩阵x、y重新计算函数值矩阵。例如:

    [x,y]=meshgrid(-5:0.1:5);

    z=peaks(x,y);

    生成的数值矩阵可以作为mesh、surf等函数参数而绘制出发多峰函数曲面图。另处,若在调用peaks函数时不带输出参数,则直接绘制出多峰函数曲面图。

    例5.19 绘制标准三维曲面图形。

    程序如下:

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    t=0:pi/20:2*pi;

    [x,y,z]=cylinder(2+sin(t),30);

    subplot(1,3,1);

    surf(x,y,z);

    subplot(1,3,2);

    [x,y,z]=sphere;

    surf(x,y,z);

    subplot(1,3,3);

    [x,y,z]=peaks(30);

    meshz(x,y,z);

      

    4.其它三维图形

    4.1 三维条形图

    1

    2

    3

    4

    5

    bar3函数绘制三维条形图,调用格式为:

    bar3(y)

    bar3(x,y)

    在第一种格式中,y的每个元素对应一个条形。

    第二种格式在x指定 的位置上绘制y中元素的条形图,X为向量,当y为向量时,x元素个数与y列数相同,当y为矩阵时,x元素与y的行数相同。

    例:1.bar3(y)

    (1)当y为矩阵时,以元素下标为坐标,以元素值为高度,绘制条形图。

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    20

    >> y=magic(5)

     

    y =

     

        17    24     1     8    15

        23     5     7    14    16

         4     6    13    20    22

        10    12    19    21     3

        11    18    25     2     9

     

    >> y(5,:)=[]%删除第五行

     

    y =

     

        17    24     1     8    15

        23     5     7    14    16

         4     6    13    20    22

        10    12    19    21     3

     

    >> bar3(y)

      

    (2)当y为向量时,也是以下标为坐标,为值为高度。

    1

    2

    3

    4

    5

    6

    7

    8

    >> y=[1 3 5 7 2]

     

    y =

     

         1     3     5     7     2

     

    >> bar3(y)

    >>

      

    2.bar3(x,y)

    (1)x为向量,y为向量

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    >> y

     

    y =

     

         1     3     5     7     9    11

     

    >> x

     

    x =

     

         1     3     5     4     8    11

     

    >> bar3(x,y)

    >>

      

    (1)x为向量,y为矩阵(x元素改变y矩阵的x坐标)

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    y =

     

        17    24     1     8    15

        23     5     7    14    16

         4     6    13    20    22

        10    12    19    21     3

     

    >> x=[1 3 5 9]

     

    x =

     

         1     3     5     9

     

    >> bar3(x,y)

      

     

    1

    <strong style="font-size: 1.5em; background-color: rgb(225, 230, 215); font-family: verdana, Arial, Helvetica, sans-serif; line-height: 1.5;">4.2 三维杆图</strong>

    1

    2

    3

    4

    stem3(z)

    stem3(x,y,z)

    第一种格式将 数据序列z 表示为xy平面向上延伸的杆图,x和y自动生成。

    第二种格式在x和y指定位置上绘制 数据序列z的杆图,x,y,z的维数必须相同。

    1.stem3(z)

    (1)z为矩阵,以下标为坐标,值为杆值

    1

    2

    3

    4

    5

    6

    7

    8

    z =

     

        17    24     1     8    15

        23     5     7    14    16

         4     6    13    20    22

        10    12    19    21     3

     

    >> stem3(z)

      

    (2)z为向量,以下标为坐标,值为杆值

    1

    2

    3

    4

    5

    6

    7

    8

    z=y(1,:)

     

    z =

     

        17    24     1     8    15

     

    >> stem3(z)

    >> stem(z)

      

      

    2.stem3(x,y,z)

    (1)x,y,z均为向量,以(x,y)为对应坐标z为值

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    x =

     

         1     2     5     9     6

     

    >> y=x

     

    y =

     

         1     2     5     9     6

     

    >> z=x

     

    z =

     

         1     2     5     9     6

     

    >> stem3(x,y,z)

    >>

      

    4.3 三维饼图

    1

    2

    pie3(x)

    其中x为向量,用x中的数据绘制一个三维饼图。

      

    1

    2

    pie3([2347 1827 2043 3025])

    pie([2347 1827 2043 3025])

      

    4.4 填充多边形

    1

    2

    fill3(x,y,z,c)

    使用x,y,z作为多边形的顶点,用c指定了填充的颜色。

      

    1

    2

    3

    4

    5

    6

    7

    8

    9

    10

    11

    12

    13

    14

    15

    16

    17

    18

    19

    fill(x,y,c)

    以(x,y)为点,c为颜色图,连点并填充点间面。

    x =

     

         1     5     6     8

     

    >> y

     

    y =

     

         2     6     4     6

     

    >> z=[1 2 3 4 5 ]

     

    z =

     

         1     2     3     4     5

     

    >> fill(x,y,z)

      

    1

    fill3(rand(3,5),rand(3,5),rand(3,5),'y')

      

    分类: matlab

    展开全文
  • plot 模块中绘制基础图的常用函数 plt.plot(x,y,fmt,…):绘制一个坐标图 plt.boxplot(data,notch,position):绘制一个箱形图 plt.bar(left,height,bottom):绘制一个条形图 plt.barh(width,bottom,left,height):...

    plot 模块中绘制基础图的常用函数

    plt.plot(x,y,fmt,…):绘制一个坐标图
    plt.boxplot(data,notch,position):绘制一个箱形图
    plt.bar(left,height,bottom):绘制一个条形图
    plt.barh(width,bottom,left,height):绘制一个横向条形图
    plt.polar(theta,r):绘制极坐标图
    plt.pie(data,explode):绘制饼图
    plt.psd(x,NFFT = 256,pad_to,Fs):绘制功率谱密度图
    plt.specgram(x,NFFT = 256,pad_to,F):绘制谱图
    plt.cohere(x,y,NFFT = 256,Fs):绘制X-Y的相关性函数
    plt.scatter(x,y):绘制散点图,其中,x和y长度相同
    plt.step(x,y,where):绘制步阶图
    plt.hist(x,bins,normed):绘制直方图
    plt.contour(X,Y,Z,N):绘制等值图
    plt.vlines():绘制垂直图
    plt.stem(x,y,linefmt,markerfmt):绘制柴火图
    plt.plot_date():绘制数据日期
    
    

    1.折线图(plt.plot()函数)

    import matplotlib.pyplot as plt
    import numpy as np
    
    x=[2016,2017,2018,2019,2020]
    y=[12,15,18,20,45]
    #x与y一一对应
    
    plt.plot(x,y)
    #生成图表
    
    #设置x,y显示的刻度
    plt.yticks([10,20,30,40,50])
    plt.xticks([2016,2017,2018,2019,2020])
    
    plt.show()
    #显示图表
    

    在这里插入图片描述
    2.条形图(又名柱状图)(plt.bar()函数)

    import matplotlib.pyplot as plt
    import numpy as np
    
    n=12
    X=np.arange(n)#arange函数,类似于range()函数,用来创建一维数组
    Y1=(1-X/float(n))*np.random.uniform(0.5,1.0,n)
    #random.uniform()从一个均匀分布[low,high)中随机采样,注意定义域是左闭右开,即包含low,不包含high.
    
    Y2=(1-X/float(n))*np.random.uniform(0.5,1.0,n)
    
    plt.bar(X,+Y1,facecolor='#9999ff',edgecolor='white')
    #facecolor填充颜色,edgecolor边框颜色
    plt.bar(X,-Y2,facecolor='#ff9999',edgecolor='white')
    
    
    for x,y in zip(X,Y1):
        #zip(a,b)函数分别从a和b中取一个元素组成元组,再次将组成的元组组合成一个新的迭代器。
        plt.text(x+0.2,y+0.05,'%.2f'%y,ha='center',va='bottom')
        #ha:horizontal alignment 水平对齐方式
        #va:纵向对齐方式
    
    for x,y in zip(X,Y2):
        #ha:horizontal alignment 水平对齐方式
        #va:纵向对齐方式
        plt.text(x+0.2,-y-0.05,'%.2f'%y,ha='center',va='top')
    
    
    plt.xlim(-5,n) #x显示的范围
    plt.xticks(())#隐藏x轴上的刻度
    plt.ylim(-1.25,1.25)
    plt.yticks(())
    
    plt.show()#显示图像
    
    

    在这里插入图片描述
    3.直方图(plt.hist()函数)

    import matplotlib.pyplot as plt
    import numpy as np
    
    #本例是标准正态分布
    mu=100  #设置均值,中心所在点
    sigma=20  #用于将每个点都扩大的倍数
    
    x=mu+sigma*np.random.randn(20000)
    #x中的点分布在mu旁边
    
    plt.hist(x,bins=100,color='green',density=True)
    #x:这个参数是arrays,指定每个x的分布位置
    #bins:指定箱子的个数,也就是有几个条形图
    #color:每个箱子的颜色
    #原来有的是normed=True,但是python3里面已经移除了,变成了density
    #density:是否对y轴数据进行标准化,若为Ture,则是在本区间的点占所有点的概率
    
    plt.show()
    
    

    这个每次运行的时候这个直方图的图形都是不同的,因为x是随机数

    4.散点图(plt.scatter()函数)

    import matplotlib.pyplot as plt
    import numpy as np
    
    x=np.random.normal(0,1,1024)
    #np.random.normal()的意思是一个正态分布,这里是随机分布一些点
    y=np.random.normal(0,1,1024)
    T=np.arctan2(y,x)#设置散点颜色,可以变成不同的颜色,无需深入研究
    
    plt.scatter(x,y,s=50,alpha=0.5)#s是指散点大小,c是指散点颜色,alpha是指散点的透明度
    
    plt.xlim(-1.5,1.5)#x轴的显示范围
    plt.ylim(-1.5,1.5)
    plt.xticks(())#同下
    plt.yticks(())#隐藏y坐标上的值
    
    plt.show()
    
    

    在这里插入图片描述
    5.等高线图(plt.scatter()函数)

    import matplotlib.pyplot as plt
    import numpy as np
    
    def f(x,y):
        #等高线公式
        return (1-x/2+x**5+y**3)*np.exp(-x**2-y**2)
    n=256
    #np.linspace在指定的间隔内返回均匀间隔的数字,下面是-3到3之间,且元素个数为256个。
    x=np.linspace(-3,3,n)#x轴
    y=np.linspace(-3,3,n)#y轴
    
    
    X,Y=np.meshgrid(x,y)#底部是网格形式
    #输入的x,y,就是网格点的横纵坐标列向量(非矩阵)
    #输出的X,Y,就是坐标矩阵。
    
    
    #等高线的颜色设置
    plt.contourf(X,Y,f(X,Y),8,alpha=0.75,cmap=plt.cm.hot)#8就是分成10部分
    #登高线的线
    C=plt.contour(X,Y,f(X,Y),8,colors='black',linewidths=0.5)
    
    #数字描述
    plt.clabel(C,inline=True,fontsize=10)
    #inline=Ture是围着数字,等于false的话是线盖着数字,fontsize是数字大小
    
    
    
    plt.xticks(())#隐藏x轴的刻度
    plt.yticks(())
    plt.show()
    
    

    在这里插入图片描述
    6.饼状图(plt.pie()函数)

    import matplotlib.pyplot as plt
    
    
    plt.rcParams['font.sans-serif']=['SimHei']
    # 正常显示中文标签
    plt.rcParams['axes.unicode_minus']=False
    # 用来正常显示负号
    
    labels='玫瑰花','月季','菊花','百合花'
    sizes=15,20,45,10
    colors='yellowgreen','gold','lightskyblue','lightcoral'
    explode=0,0.1,0,0
    plt.pie(sizes,explode=explode,labels=labels,colors=colors,autopct='%1.1f%%',shadow=True,startangle=50)
    #sizes指不同扇区的面积(数值),labels是指标签,
    #explode是指英文意思是爆炸,表示把饼状图的各个部分像爆炸一样,从中心向四周分散出去。,哪个是0.1哪个往外出来,
    #autopct表示自动进行百分比运算,autopct等号后面的值,.2就代表小数点后精确到2位,如果是3就是小数点后精确到3位,f代表浮点数。
    #shadow=Ture是指阴影,使饼图更立体,false是没阴影,
    #startangle=50表示开始第一块是正时针50度的角度开始画,且逆时针画的
    plt.axis('equal')
    #设置x,y的刻度一样,使其饼图为正圆
    plt.show()
    
    

    在这里插入图片描述
    7.3D图(Axes3D 对象的 plot_surface() 函数)

    import numpy as np
    import matplotlib.pyplot as plt
    from mpl_toolkits.mplot3d import Axes3D
    
    fig=plt.figure()
    ax=Axes3D(fig)
    
    X=np.arange(-4,4,0.25)
    Y=np.arange(-4,4,0.25)
    X,Y=np.meshgrid(X,Y)#变成矩阵,也就是网格
    R=np.sqrt(X**2+Y**2)#计算,得出z轴
    Z=np.sin(R)  #z轴
    
    ax.plot_surface(X,Y,Z,rstride=1,cstride=1,cmap=plt.get_cmap('rainbow'))
    #rstide,线与线之间的行距,cstride线与线之间的列跨,数值越大,线与线之间的距离越大
    #cmap,表示映射颜色,这里的 rainbow表示彩虹色
    
    #等高线的设置,这里是底下那个面
    ax.contourf(X,Y,Z,zdir='z',offset=-2,cmap='rainbow')
    #zdir表示从哪个轴压下去,
    #offset表示把等高线图放到距离x=0轴平面多远的地方,整数表在上面,负数表示在下面
    
    ax.set_zlim(-2,2)#z轴让他显示的高度范围
    plt.show()
    
    

    在这里插入图片描述

    展开全文
  • 本文由@浅墨_毛星云出品,首发于...作为基于物理的渲染(PBR)技术中材质高光质感的决定因素,更先进的法线分布函数(Normal Distribution Function,NDF)的问世和发展,是PBR能够在游戏和电影工业日益普及的重要...

            本文由@浅墨_毛星云 出品,首发于知乎专栏,转载请注明出处  

            文章链接: https://zhuanlan.zhihu.com/p/69380665

     

    作为基于物理的渲染(PBR)技术中材质高光质感的决定因素,更先进的法线分布函数(Normal Distribution Function,NDF)的问世和发展,是PBR能够在游戏和电影工业日益普及的重要推动力之一。

    从法线分布函数和微平面理论(Microfacet Theory)的视角来看,基于物理的渲染代表着一种从宏观表现到微观细节的渲染理念的进化,从而对材质有了亚像素级更精细的把控和更科学的定量表示,从而推动了游戏和电影业界渲染品质的提升,以及工作流的升级。

    图 基于PBR渲染的《雷神3:诸神黄昏》 @Arnold Renderer

    图 基于PBR渲染的《银河护卫队2》@Arnold Renderer

    图 基于PBR渲染的真实感枪械 @ArtSation,linus scheffel

    图 基于PBR渲染的Agent 327@ Blender Animation Studio

    图 基于PBR渲染的场景 @BattleField V

    历史上主流的法线分布函数,按提出时间进行排序,可以总结为:

    • Berry [1923]

    • Beckmann [1963]

    • Phong [1973]

    • Blinn-Phong [1977]

    • ABC [1989]

    • GGX [2007] / Trowbridge-Reitz [1975]

    • Shifted Gamma Distribution,SGD [2012]

    • Trowbridge-Reitz(GTR)[2012]

    • Student’s T-Distribution , STD [2017]

    • Exponential Power Distribution , EPD [2017]

    而在如今的PBR时代,业界主流的法线分布函数已从传统的Blinn-Phong等分布,迁移到更接近真实材质外观表现,具备能量守恒,具有更宽尾部和更高峰值的GGX分布,而且在朝着多高光波瓣(multiple specular lobes)的方向发展。

    本文将从以下几个方面,对基于物理的渲染中法线分布函数(Normal Distribution Function,NDF)相关内容进行一个系统的总结:

    • 基于物理的渲染理念:从宏观表现到微观细节

    • 法线分布函数与微平面理论

    • 法线分布函数的基本性质

    • 各项同性NDF相关总结

    • 法线分布函数的形状不变性

    • 各向异性NDF相关总结

    • NDF的性能优化

    • 多高光波瓣(multiple specular lobes)

    • PBR中的高光抗锯齿

    • 主流NDF的局限性和发展趋势

    按照惯例,开始正文前,先放出本文内容总结的思维导图:

    OK,让我们从第一部分开始。



     

    一、基于物理的渲染理念:从宏观表现到微观细节

    关于基于物理的渲染理念,其实跟图形学中对几何体的建模尺度有一定关联。图形学中,对几何体外观的建模,总会假设一定的建模尺度和观察尺度:

    • 宏观尺度(Macroscale), 几何体通过三角形网格进行建模, 由顶点法线(Vertex Normal)提供每顶点法线信息

    • 中尺度(Mesoscale), 几何体通过纹理进行建模,由法线贴图(Normal Map)提供每像素法线信息

    • 微观尺度(Microscale), 几何体通过BRDF进行建模,由粗糙度贴图(Roughness Map)配合法线分布函数,提供每亚像素(subpixel)法线信息

    图 PBR物体的渲染建模。从左到右:渲染后物体,几何体三角形网格,法线贴图,粗糙度贴图。(图片来自Sébastien Lagarde ,Physically-Based Material, where are we,SIGGRAPH 2017)

    传统光照模型中,一般只将几何体建模到中尺度的法线贴图(Normal Map)层面。虽说Blinn-Phong等分布也是基于微平面理论推导而来,但并没有配套粗糙度贴图(Roughness Map)为其提供亚像素级精度的细节,而且传统的NDF一般都没有经过归一化,不满足能量守恒,容易出现失真。

    而在基于物理的渲染工作流中,通过将粗糙度贴图(Roughness Map)与微平面归一化的法线分布函数结合使用,将需渲染的几何体的建模尺度细化到了微观尺度(Microscale)的亚像素层面,对材质的微观表现更加定量,所以能够带来更加接近真实的渲染质量和更全面的材质外观质感把控。如下图。

    图 多尺度的几何建模表示(图片来自Sébastien Lagarde ,Physically-Based Material, where are we,SIGGRAPH 2017)



     

    二、法线分布函数与微平面理论

    • 微平面理论(microfacet theory)作为一种研究微观几何(microgeometry)对反射率影响的数学分析方法,基于将微观几何(microgeometry)建模为微平面(microfacets)的集合的思想。其最初由光学物理领域开发,用于研究统计表面上的散射[Beckmann Spizzichino 1963]。在图形社区中,我们使用它来推导基于物理的BRDF。

    • 微平面模型的一个重要特性是微平面法线m的统计分布(statistical distribution)。 此分布由曲面的法线分布函数(Normal Distribution Function,NDF)定义。

    图 微观表面(Micro Surface)与宏观表面(Macro Surface)法线方向(图片来自[Walter 2007])

    • 法线分布函数(Normal Distribution Function,NDF)在一些文献中也用Specular D进行表示。

    • 微平面的法线分布函数D(m)描述了微观表面上的表面法线m的统计分布。给定以m为中心的无穷小立体角dωm和无穷小宏观表面区域dA,则是相应微表面部分的总面积,其法线位于指定的立体角内。因此NDF的本质是一个密度函数,单位为1/球面度(1/steradians)。

    • 从直觉上来说,NDF就像是微平面法线分布的直方图(histogram)。 它在微平面法线更可能指向的方向上具有更高的值。大多数表面都具有在宏观表面法线n处显示出很强的法线分布峰值。

    • 若以函数输入和输出的角度来看NDF,则其输入为微表面粗糙度(微表面法线集中程度)和宏观法线与视线的中间矢量(微表面法线方向),输出为此方向上的微表面法线强度。

    • 一般我们用宏观表面的半矢量h来表示微观表面法线m,因为仅m = h的表面点的朝向才会将光线l反射到视线v的方向,其他朝向的表面点对BRDF没有贡献(正负相互抵消)。

    图 仅m = h的表面点的朝向才会将光线l反射到视线v的方向,其他表面点对BRDF没有贡献(图片来自《Real-Time Rendering 4th》)

    • 目前业界广泛采用的Microfacet Cook-Torrance BRDF形式如下:

    • 其中 D(h) 即法线分布函数 (Normal Distribution Function),描述了平面法线分布的概率,即具有正确朝向的微表面法线浓度。即具有正确朝向,能够将来自l的光反射到v的表面点的相对于表面面积的浓度。D(h)常被直接写作D(m)。

    • 可以将法线分布函数 D(m) 理解为微观几何表面区域上的微平面表面法线的统计分布。 对 D(m) 在整个微平面法线上积分,会得到微表面的面积。更有用的是对D(m)(n · m)进行积分,即将 D(m) 的投影到宏观表面平面上,会得到宏观表面片元(patch)的面积,其被约定等于1,如下图所示:

    图 微表面的侧视图。积分D(m)(n · m),得到投影到宏观表面平面上的微平面区域,等于宏观表面的面积,被约定为1。(图片来自《Real-Time Rendering 4th》)

    • 换句话说,投影D(m)(n · m)是被归一化的:

    • 上式在积分时用到了Θ符号,表示在在整个球体上积分。而在以n为中心的半球上积分时,一般用Ω表示。实际上,图形学中使用的大多数微结构模型都是高度场(heightfields),这意味着对于Ω外的所有方向,D(m) = 0。 但是,上式也适用于非高度场微观结构。

    • 更一般地,微观表面(microsurface)和宏观表面(macrosurface)在垂直于任何视图方向v的平面上的投影是相等的:

    • 上面两个积分公式定义了成为一个合法的基于物理的法线分布函数必须服从的约束。

    • 另外,上面两个积分公式中的点积不用被约束为大于等于0,因为投影会产生正负抵消,如下图所示。

    图 对D(m)(v · m)进行积分,微平面区域投影到垂直于v的平面(图中黑色的cos θo线段),产生宏观表面到该平面的投影,即cos θo或(v·n)。当多个微平面的投影重叠时,背向(backfacing)微平面的负投影区域抵消了“额外的”前向(frontfacing)微平面。(图片来自《Real-Time Rendering 4th》)

    可以发现,尽管存在许多具有重叠(overlapping)投影的微平面,但用于最终渲染而言,我们仅关注可见的微平面,即在每个重叠集合中最接近相机的微平面。 这一事实表明了将投影的微观区域与投影的宏观几何区域相关联的另一种方法:可见微平面的投影面积之和等于宏观表面的投影面积。 我们可以通过定义遮蔽函数(masking function)G1(m,v)来对其进行数学表达,其给出了沿着视图向量v可见的具有法线m的微平面的比率。G1(m,v)D(m) 在球体上的积分给出投影到垂直于v的平面上的宏观表面的面积(其中x+ 表示为被约束为大于等于0):

    图 对可见微平面的投影区域(亮红色)进行积分,得到宏观表面在垂直于v的平面上的投影面积

    • 上式中,通过x+的表示方法表达将v · m限制为大于等于0。 背面微平面不可见,因此在这种情况下不计算它们。 乘积G1(m,v)D(m)则表示了可见法线的分布

    • 上式对G1(m,v)施加约束,但并不能唯一地确定它。有无数个函数满足给定微平面法线分布D(m)的约束。 这是因为D(m)没有完全指定微表面(microsurface)。 它仅告诉我们有多少百分比的微平面(microfacets)的法线指向了某些方向,而没有告诉我们这些法线是如何进行排列。

    本文的重点是法线分布函数D,遮蔽阴影函数G相关的讨论暂时就不继续展开,下篇文章将对遮蔽阴影函数G做更详细的探讨。



     

    三、法线分布函数的基本性质

    一个基于物理的的微平面法线分布的基本性质,可以总结如下:

    1. 微平面法线密度始终为非负值:

    2. 微表面的总面积始终不小于宏观表面总面积:

    3. 任何方向上微观表面投影面积始终与宏观表面投影面积相同:

    4. 若观察方向为法线方向,则其积分可以归一化。即v = n时,有



     

    四、各项同性NDF相关总结

    最常见的法线分布函数是各向同性(isotropic)的,它们围绕由宏观表面法线n定义的轴旋转对称(rotationally symmetrical)。常见的各项同性法线分布函数按出现时间进行排序,可以总结如下:

    • Berry [1923]

    • Beckmann [1963]

    • Phong [1973]

    • Blinn-Phong [1977]

    • ABC [1989]

    • GGX [2007] / Trowbridge-Reitz [1975]

    • Shifted Gamma Distribution,SGD [2012]

    • Trowbridge-Reitz(GTR)[2012]

    • Student’s T-Distribution , STD [2017]

    • Exponential Power Distribution , EPD [2017]

    本文将对其中在历史上使用和讨论较为广泛的几种法线分布函数,Blinn-Phong、Beckmann、GGX(Trowbridge-Reitz)和GTR进行总结。

    4.1 Blinn-Phong分布

    • Blinn-Phong分布的前身Phong分布是计算机图形学文献中提出的最早的着色方程之一。

    • Blinn-Phong法线分布函数由Blinn推导出,作为(非基于物理的)Phong着色模型的改进,以更好地拟合微平面BRDF的结构。

    • Blinn-Phong 分布不具备形状不变性(shape-invariant)。

    • 虽然Blinn在提出Blinn-Phong分布时没有指定归一化因子,但很容易计算,关于Blinn-Phong的归一化的讨论,可以参见(http://www.thetenthplanet.de/archives/255)。下式是较为主流的归一化的Blinn-Phong(Normalized Blinn-Phong)的形式:

    • 其中,幂αp是Blinn-Phong NDF的“粗糙度参数”;高值表示光滑表面,低值表示粗糙表面。对于非常光滑的曲面,值可以任意高(一个完美的镜面αp=∞),并且通过将αp设置为0可以实现最大随机曲面(均匀NDF)。

    • αp参数不便于艺术家操纵或直接绘制,因为它带来的视觉变化非常不均匀。出于这个原因,经常让美术师们操纵“界面值”,即通过非线性函数从中导出αp。例如:αp=ms,其中s是0到1之间的艺术家操纵值,m是给定的电影或游戏中αp的上限。这种映射被多款游戏使用,包括《使命召唤:黑色行动(Call of Duty: Black Ops)》,其中m被设置为值8192。

    • UE4中,则采用映射 ,那么得到的Blinn-Phong的形式为:

    • UE4中对Blinn-Phong的实现代码如下:
    // [Blinn 1977, "Models of light reflection for computer synthesized pictures"]
    float D_Blinn( float a2, float NoH )
    {
            float n = 2 / a2 - 2;
            return (n+2) / (2*PI) * PhongShadingPow( NoH, n );		// 1 mad, 1 exp, 1 mul, 1 log
    }
    

    4.2 Beckmann分布

    • Beckmann分布是光学业界开发的第一批微平面模型中使用的法线分布。[Beckmann 1963],也是Cook-Torrance BRDF在提出时选择的NDF [Cook 1981] [Cook 1982]。

    • Beckmann NDF具备形状不变性(shape-invariant)

    • 正确归一化后,Beckmann分布具有以下形式:

    • Beckmann分布在某些方面与Phong分布非常相似。 两种法线分布的参数可以使用关系式 进行等效。

    图 Blinn-Phong(蓝色虚线)和Beckmann(绿色)分布,αb值在0.025到0.2之间(使用参数关系)。(图片来自《Real-Time Rendering 4th》)

    UE4中对Beckmann分布的实现代码如下:

    // [Beckmann 1963, "The scattering of electromagnetic waves from rough surfaces"]
    float D_Beckmann( float a2, float NoH )
    {
    	float NoH2 = NoH * NoH;
    	return exp( (NoH2 - 1) / (a2 * NoH2) ) / ( PI * a2 * NoH2 * NoH2 );
    }
    

    4.3 GGX(Trowbridge-Reitz)分布

    • GGX即Trowbridge-Reitz分布,最初由Trowbridge和Reitz [Trowbridge 1975]推导出,在Blinn 1977年的论文 [Blinn 1977]中也有推荐此分布函数,但一直没有受到图形学界的太多关注。30多年后,Trowbridge-Reitz分布被Walter等人独立重新发现[Walter 2007],并将其命名为GGX分布。

    • 在[Walter 2007]重新发现并提出GGX分布之后,GGX分布采用风潮开始在电影 [Burley 2012]和游戏 [Karis 2013][Lagarde 2014]行业中广泛传播,成为了如今游戏行业和电影行业中最常用的法线分布函数。

    • GGX分布的公式为:

    • 在流行的模型中,GGX拥有最长的尾部。这是GGX能够日益普及的主要原因:

    图 主流法线分布函数高光长尾对比

    • GGX分布具备形状不变性(shape-invariant),而与其对标的GTR等分布不具备形状不变性,这是GGX能普及的深层次原因。

    • 在迪士尼原理着色模型(Disney principled shading model)中,Burley推荐将粗糙度控制以α= r2暴露给用户,其中r是0到1之间的用户界面粗糙度参数值,以让分布以更线性的方式变化。这种方式实用性较好,不少使用GGX分布的引擎与游戏都采用了这种映射,如UE4和Unity。

    // GGX / Trowbridge-Reitz
    // [Walter et al. 2007, "Microfacet models for refraction through rough surfaces"]
    float D_GGX( float a2, float NoH )
    {
    	float d = ( NoH * a2 - NoH ) * NoH + 1;	// 2 mad
    	return a2 / ( PI*d*d );			// 4 mul, 1 rcp
    }
    

    4.4 Generalized-Trowbridge-Reitz(GTR)分布

    • Burley [ Burley 2012]根据对Berry,GGX等分布的观察,提出了广义的Trowbridge-Reitz(Generalized-Trowbridge-Reitz,GTR)法线分布函数,其目标是允许更多地控制NDF的形状,特别是分布的尾部:

    • 其中,γ参数用于控制尾部形状。 当γ= 2时,GTR等同于GGX。 随着γ的值减小,分布的尾部变得更长。而随着γ值的增加,分布的尾部变得更短。上式中:

    • γ=1时,GTR即Berry分布

    • γ=2时,GTR即GGX(Trowbridge-Reitz)分布

    以下为各种γ值的GTR分布曲线与θh的关系图示:

    图 各种γ值的GTR分布曲线与θh的关系

    • GTR分布不具备形状不变性(shape-invariant),导致其发布以来,无法被广泛使用。

    以下是γ= 1和γ= 2时GTR分布的Shader实现代码:

    // Generalized-Trowbridge-Reitz distribution
    float D_GTR1(float alpha, float dotNH)
    {
        float a2 = alpha * alpha;
        float cos2th = dotNH * dotNH;
        float den = (1.0 + (a2 - 1.0) * cos2th);
    
        return (a2 - 1.0) / (PI * log(a2) * den);
    }
    
    float D_GTR2(float alpha, float dotNH)
    {
        float a2 = alpha * alpha;
        float cos2th = dotNH * dotNH;
        float den = (1.0 + (a2 - 1.0) * cos2th);
    
        return a2 / (PI * den * den);
    }
    

    4.5 其他分布

    2017年以来新发布的学生T-分布(Student’s T-Distribution , STD) [ Ribardière 2017]和指数幂分布(exponential power distribution ,EPD)[ Holzschuch 2017] NDF包括形状控制参数。两者都具有形状不变性,但由于发布时间较新,实用性尚不明朗,目前很少有听说有实际使用。

    五、法线分布函数的形状不变性

    • 形状不变性(shape-invariant)是一个合格的法线分布函数需要具备的重要性质。具有形状不变性(shape-invariant)的法线分布函数,可以用于推导该函数的归一化的各向异性版本,并且可以很方便地推导出对应的遮蔽阴影项G。

    • 若一个各向同性的NDF可以改写成以下形式,则这个NDF具有形状不变性(shape-invariant):

    还有一种等价的写法是:

    • 对于形状不变的NDF,缩放粗糙度参数相当于通过倒数拉伸微观几何,如下图所示。

    图 对于形状不变的NDF,缩放粗糙度参数相当于通过倒数拉伸微观几何(图片来自Naty Hoffman, Recent Advances in Physically Based Shading, SIGGRAPH 2016)

    • 为了更容易理解形状不变性,可以将NDF视为P22,即一个2D斜率的分布。原始的NDF是一个3D的矢量分布。然后我们可以发现,对于具有形状不变形式的分布,线性缩放粗糙度α会导致斜率空间中的分布线性拉伸。

    • 2D斜率的分布P22和法线分布函数的关系以及变换式可写作:

    • 关于形状不变性的好处,可以总结为:

      • 方便推导出该NDF归一化的各向异性版本

      • 方便推导出遮蔽阴影项 Smith G()

      • 方便基于NDF或可见法线分布推导其重要性采样

        • 对于Smith G(),可用低维函数或表格处理所有粗糙度和各向异性
    • 对于常见的法线分布函数的旋转不变性的分类:

      • 具备形状不变性的常用法线分布函数为:

        • GGX

        • Beckmann

      • 不具备形状不变性的常用法线分布函数为:

        • Phong

        • Blinn-Phong

        • GTR

    • 下图显示了如何通过拉伸表面(stretching the surface)将各向同性的形状不变分布转换为各向异性分布。相反,任何具有各向异性分布的配置都可以转换回具有各向同性分布的配置。

    图 通过拉伸表面,可以将各向同性具备形状不变性的分布,转换为各向异性分布。(图片来自[Heitz 2014])

    通过使用这种方法,我们可以推导出Beckmann和GGX分布的各向异性形式。

    六、各向异性NDF相关总结

    现实世界中,大多数材质具有各向同性的表面外观,但有些特殊材质的微观结构具有显著的各向异性(Anisotropy),从而显著影响其外观。

    图 各向异性(Anisotropy)材质渲染表现 @Arnold Renderer

    图 各向异性(Anisotropy)材质渲染表现 @Arnold Renderer

    创建各向异性NDF的常用方法是基于现有各向同性NDF进行推导。而推导所使用的方法是通用的,可以应用于任何形状不变的(shape-invariant)各向同性NDF,这便是GGX等形状不变的NDF能更加普及的另一个原因。

    如上文所述,若一个各向同性(isotropic)的NDF具备形状不变性(shape-invariant),则其可以用以下形式写出:

    其中g代表一个表示了NDF形状的一维函数。而通过此形式,可得到各向异性的(anisotropic)版本:

    • 其中,参数αxαy分别表示沿切线(tangent)方向t和副法线(binormal)方向b的粗糙度。若αx = αy,则上式缩减回各向同性形式。

    • 需要注意的是,一般的shader写法,会将切线方向t写作X,副法线(binormal)b方向写作Y。

    6.1 Anisotropic Beckmann Distribution

    各项异性的Beckmann分布形式如下:

    • 其中,m为微表面法线(可以理解为half半矢量),n为宏观表面法线,t为切线方向,b为副法线方向。

    以下为UE4中Anisotropic Beckmann分布的Shader实现代码:

    // Anisotropic Beckmann
    float D_Beckmann_aniso( float ax, float ay, float NoH, float3 H, float3 X, float3 Y )
    {
    	float XoH = dot( X, H );
    	float YoH = dot( Y, H );
    	float d = - (XoH*XoH / (ax*ax) + YoH*YoH / (ay*ay)) / NoH*NoH;
    	return exp(d) / ( PI * ax*ay * NoH * NoH * NoH * NoH );
    }
    

    6.2 Anisotropic GGX Distribution

    各项异性的GGX分布形式如下:

    以下为UE4中Anisotropic GGX分布的Shader实现代码:

    // Anisotropic GGX
    // [Burley 2012, "Physically-Based Shading at Disney"]
    float D_GGXaniso( float ax, float ay, float NoH, float3 H, float3 X, float3 Y )
    {
    	float XoH = dot( X, H );
    	float YoH = dot( Y, H );
    	float d = XoH*XoH / (ax*ax) + YoH*YoH / (ay*ay) + NoH*NoH;
    	return 1 / ( PI * ax*ay * d*d );
    }
    
    • 其中,X为tangent,t切线方向,Y为binormal,b,副法线方向

    • 需要注意的是,将法线贴图与各向异性BRDF组合时,重要的是要确保法线贴图扰动(perturbs)切线和副切线矢量以及法线。

    6.3 其他各项异性参数化方法

    • 虽然参数化各向异性NDF的最直接的方法是使用各向同性粗糙度进行两次参数化,一次用于αx,一次用于αy,有时也使用其他参数化。 在迪士尼原理着色模型中,各向同性粗糙度参数r与第二标量参数kaniso组合,范围为[0,1]。 因此,从这些参数计算αxαy的值:

    • 其中,上式中的因子0.9将纵横比限制为10:1。

    • Sony Imageworks 则使用了一种不同的参数化,允许任意程度的各向异性:



     

    七、NDF的性能优化

    7.1 Blinn-Phong不一定比GGX更省性能

    需要注意的是,Normalized Blinn-Phong不一定比GGX更省,要具体看GPU架构。让我们从指令数对两者的计算量进行量化。两者的计算Shader代码如下:

    // [Blinn 1977, "Models of light reflection for computer synthesized pictures"]
    float D_Blinn( float a2, float NoH )
    {
    	float n = 2 / a2 - 2;
    	return (n+2) / (2*PI) * PhongShadingPow( NoH, n );	// 1 mad, 1 exp, 1 mul, 1 log
    }
    
    // GGX / Trowbridge-Reitz
    // [Walter et al. 2007, "Microfacet models for refraction through rough surfaces"]
    float D_GGX( float a2, float NoH )
    {
    	float d = ( NoH * a2 - NoH ) * NoH + 1;	// 2 mad
    	return a2 / ( PI*d*d );			// 4 mul, 1 rcp
    }
    

    若统计一下两种分布函数的计算指令,可以得到:

    • 标准的GGX计算指令为两次mad,4次mul,一次rcp。共7次运算。

    • 标准的Blinn-Phong计算指令为1次mad, 1次exp, 1次mul, 1次log。共4次运算。

    不难发现,Blinn-Phong虽然比GGX的总运算次数少3次,但具有exp、log等稍复杂的运算。两者的性能差异主要还是要看GPU架构,对某些架构的GPU而言, GGX可能会更快。

    7.2 GGX分布的移动端性能优化

    标准的GGX的公式和一般Shader实现如下:

    // GGX / Trowbridge-Reitz
    // [Walter et al. 2007, "Microfacet models for refraction through rough surfaces"]
    float D_GGX( float a2, float NoH )
    {
    	float d = ( NoH * a2 - NoH ) * NoH + 1;	// 2 mad
    	return a2 / ( PI*d*d );			// 4 mul, 1 rcp
    }
    

    在上述实现中,用float进行数据的存储与计算。其实我们可以通过使用半精度浮点数(half precision floats)来对此实现进行改进。这种优化需要改变原始方程,因为在半浮点数half(即mediump)中计算1-(n·h)^2时存在两个问题:

    由于n和h都是单位矢量,那么| n×h |^ 2 = 1-(n·h)^2。

    于是,我们可以通过使用简单的叉积| n×h |^2来直接计算半精度浮点数下的1-(n·h)^2。

    总的来说,此优化方案会带来更好的性能,并保持所有计算都在half(mediump)内进行。

    UE4对GGX的移动端优化的Shader实现代码如下:

    #ifndef MOBILE_GGX_USE_FP16
    	#define MOBILE_GGX_USE_FP16 1
    #endif
    
    #define MEDIUMP_FLT_MAX    65504.0
    #define MEDIUMP_FLT_MIN    0.00006103515625
    
    #if MOBILE_GGX_USE_FP16
    	#define saturateMediump(x) min(x, MEDIUMP_FLT_MAX)
    #else
    	#define saturateMediump(x) (x)
    #endif
    
    half GGX_Mobile(half Roughness, half NoH, half3 H, half3 N)
    {
    
    #if MOBILE_GGX_USE_FP16
    	float3 NxH = cross(N, H);
    	float OneMinusNoHSqr = dot(NxH, NxH);
    #else
        float OneMinusNoHSqr = 1.0 - NoH * NoH;
    #endif
    
    	half a = Roughness * Roughness;
    	float n = NoH * a;
    	float p = a / (OneMinusNoHSqr + n * n);
    	float d = p * p;
    	return saturateMediump(d);
    }
    

    更多关于NDF性能优化的讨论,由于篇幅原因,这边就不一一展开了。这里列举一些其他的法线分布函数优化相关的材料供大家参考:



     

    八、多高光波瓣(multiple specular lobes)

    • 即便是GGX的高光长尾,仍然不足以媲美真实世界材质中的高光表现。在不增加NDF本身的复杂性的前提下,更好地匹配测量材质的替代解决方案是使用多个高光波瓣(multiple specular lobes)。

    • Cook 和Torrance [Cook Torrance1981]首先提出了这个想法。之后 [Ngan 2005]进行了实验测试,他发现对于大部分材质来说,添加第二个高光波瓣确实显着改善了贴合性。

    • 皮克斯PxrSurface材质具有“粗糙镜面反射(rough specular)”波瓣,旨在用于此目的(与主镜面波瓣一起使用)。附加波瓣是一个完整的镜面微平面BRDF,包含所有相关的参数与项。

    • Sony Imageworks使用更外部的方法,使用两个GGX NDF的混合作为扩展的NDF暴露给用户,而不是整个单独的镜面BRDF项。在这种情况下,所需的唯一附加参数是第二粗糙度值和混合量。

    • Disney Principled BRDF也使用了两个固定的高光波瓣(specular lobe),且都使用GTR分布。 主波瓣(primary lobe)使用γ=2的GTR(即GGX分布),代表基础底层材质(Base Material)的反射,可为各项异性(anisotropic)或各项同性(isotropic)的金属或非金属。次级波瓣(secondary lobe)使用γ=1的GTR(即Berry分布),代表基础材质上的清漆层(ClearCoat Layer)的反射,一般为各项同性(isotropic)的非金属材质,即清漆层(ClearCoat Layer)

    • 另外,多高光波瓣(multiple specular lobes)常与Layered mixture model结合使用。

    图 离线渲染器的多层混合建模(Layered mixture model)。 @ Autodesk Standard Surface

    图:Autodesk Standard Surface多层建模的multiple specular lobes材质。整体框架为“原子”闭包的加权和,每个闭包的重量是沿着从其相应叶节点到根节点的路径的边缘重量的乘积。shader参数权重以粗体显示。@ Autodesk Standard Surface

    图 通过组合clearcoat层和标准层的多高光波瓣(multiple specular lobes)架构,可以实现各种效果。 从左到右:铝花(flakes),表面水滴(raindrops),碳纤维(carbon fiber)。@Arnold Renderer



     

    九、PBR中的高光抗锯齿

    • 锯齿(Aliasing)是实时渲染和图形学中经常会面对的问题。而PBR由于使用了标准化的法线分布函数(normalized NDF),以及无处不在的反射现象,加上实时渲染中较少的采样率,让其高光的锯齿问题更加明显。这导致了基于物理的渲染中,高光锯齿是实践中经常会遇到的问题。

    • 模型精度越高、工作流越倾向于全PBR方式、光照计算精确程度越高,则反射的高光锯齿问题就越明显。

    • 以下是带高光锯齿的PBR渲染图和经过TAA处理锯齿后的对比图:

    图 PBR赛车渲染(未经过抗锯齿处理) @CTAA , Unity AssetStore

    图 PBR赛车渲染(基于TAA(Temporal Anti-Aliasing,时域抗锯齿)处理)@CTAA , Unity AssetStore

    • 出现高光锯齿的原因可以总结为: 法线分布函数作为亚像素表面结构(subpixel surface structure)的统计描述。 当相机和表面之间的距离增加时,先前覆盖多个像素的表面结构会减小到亚像素(subpixel)大小,从而从法线贴图的领域移动到法线分布函数的亚像素领域。在亚像素领域,纹理的mipmap一般以平均的方式进行处理,会丢失原有的细节,从而让该像素处的法线分布过于狭窄和集中,于是便会出现高光在像素级别的不连续性,以闪烁高光(flickering highlights)形式引起锯齿。

    • 关于高光锯齿,业界的解决方案分为两大流派:屏幕空间抗锯齿(Anti-Aliasing)和预过滤(Pre-Filtering),下面分别进行总结:

    • 屏幕空间抗锯齿(Anti-Aliasing)。 MSAA(MultiSample Anti-Aliasing,多采样抗锯齿),SSAA(SuperSample Anti-Aliasing,超采样抗锯齿), FXAA(Fast Approximate Anti-Aliasing,快速近似抗锯齿)和TAA(Temporal Anti-Aliasing,时域抗锯齿)等抗锯齿技术可以求解单像素上多个点的入射光,找出其中小的变化点,从而减少可见的锯齿。其中目前较为有效的PBR高光抗锯齿的技术,当属TAA(Temporal Anti-Aliasing,时域抗锯齿)。有关TAA技术的更多详细信息,可以参考UE4在SIGGPRACH 2014中的分享:http://advances.realtimerendering.com/s2014/epic/TemporalAA.pptx

    • 预过滤(Pre-Filtering)。 类似Toksvig,LEAN Mapping、CLEAN Mapping和LEADR Mapping等技术方案,按照像素覆盖区域,将宏观几何(macro- geometric)(曲率)和中观几何(meso-geometric)(法线贴图,位移贴图)的变化,转移到材质属性的微观几何变化,来保持采样数量较少。这种变换可以更容易和更快地求解像素覆盖区域内发生的所有交互。其中,Toksvig和LEAN Mapping专注于法线贴图的过滤,LEAN映射的简单变体CLEAN Mapping需要较少的存储,代价是失去各向异性支持,而LEADR Mapping则用于位移贴图的过滤。而其他技术则通过将曲率转换为材质属性来近似宏观的几何过滤。Stephen Hill的blog文章(https://blog.selfshadow.com/2011/07/22/specular-showdown/)对此内方案进行了结合改进,而这边是一个对应的WebGL实现demo:http://www.selfshadow.com/sandbox/gloss.html



     

    十、主流NDF的局限性和发展趋势

    当前主流NDF的局限性可以主要总结为如下三点:

    • 缺少更好的形状控制NDF

    • 无法表示粗粒度微观结构

    • 单次散射建模的局限性

    下面分别进行说明。

    9.1 缺少更好的形状控制NDF

    • 现有的主流NDF缺少更好的形状不变性(Shape Invariance)+形状控制(Shape Control)的结合。对此,Naty Hoffman在SIGGRAPH 2016上提出,广义Beckmann分布(Generalized Beckmann)和超柯西分布(Hyper-Cauchy)可以作为实践中的选择。

    图 广义Beckmann分布(Generalized Beckmann)(图片来自Naty Hoffman, Recent Advances in Physically Based Shading, SIGGRAPH 2016)

    图 超柯西分布(Hyper-Cauchy)(图片来自Naty Hoffman, Recent Advances in Physically Based Shading, SIGGRAPH 2016)

    9.2 现有NDF无法表示粗粒度微观结构

    • 当今使用的NDF从外观而言都很平滑,如下图中左边的NDF。这种NDF每个像素覆盖了数万个表面细节,是对细粒度的微观几何的一种良好表示形式。

    图 细粒度的NDF vs 粗粒度的NDF(图片来自Naty Hoffman, Recent Advances in Physically Based Shading, SIGGRAPH 2016)

    • 但其实真实世界中的许多表面材质,具有粗粒度的微观结构,像素仅覆盖了几十个表面元素。在这种情况下,法线分布的表现更像是如如上图的右边所示,表面有一个复杂而闪烁外观,而不仅仅的各项异性这么简单。目前提出的模型都无法表示出这种类型的法线分布。期待未来有更多能解决此问题的法线分布函数的问世。

    图 真实世界中的法线分布(图片来自[Yan 2014])

    9.3 单次散射建模的局限性和发展趋势

    • 目前渲染领域广泛采用的Cook-Torrance microfacet BRDF微平面模型,实际上是人们可以想到的最简单的模型,它仅对几何光学系统中的单层微表面上的单次散射进行建模。没有考虑多次散射,没有考虑衍射,也没有考虑波动光学。其假设所有遮挡的光线都被丢失,会导致与现实行为相比的能量损失。
    • 对此,一些论文提出,可以采用非基于物理的修正因子(corrective factors)来尝试对缺失的能量进行补偿。例如迪士尼模型中的的“Sheen”光泽项。
    • 图 现有Microfacet建模未考虑图中蓝色部分的multiple surface bounce反射(图片来自Naty Hoffman, Recent Advances in Physically Based Shading, SIGGRAPH 2016)

       

    • 另外,近年不少渲染器也开始结合multiple-scattering BSDF使用多次散射进行建模,如cycles renderer的 Multiscatter GGX(https://developer.blender.org/D2002)。但目前的multiple-scattering BSDF方案主要为随机求解,所以不适用于实时渲染和游戏领域。

    图 cycles renderer Multiscatter GGX @ cycles renderer

    虽然Multiscatter GGX等方案离实时渲染还有很长一段路。这边也列举一些相关资料,感兴趣的同学不妨根据需要深入研究。



     

    十一、本文内容要点总结

    正文到这里已经结束。不妨使用本文主要内容提炼出的思维导图作为全文的内容总结:



     

    Reference

    [1] Akenine-Moller T, Haines E, Hoffman N. Real-Time Rendering 4th[M]. AK Peters CRC Press, 2018

    [2] Blinn, James F., "Models of Light Reflection for Computer Synthesized Pictures," ACM Computer Graphics (SIGGRAPH '77 Proceedings), vol. 11, no. 2, pp. 192-198, July 1977. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.131.7741&rep=rep1&type=pdf

    [3] Walter, Bruce, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance, "Microfacet Models for Refraction through Rough Surfaces," Rendering Techniques 2007, Eurographics Association, pp. 195-206, June 2007.https://www.cs.cornell.edu/\~srm/publications/EGSR07-btdf.pdf

    [4] https://lesterbanks.com/2014/01/v-ray-new-ggx-shader-offers-better-microfacet-distribution-over-blinn-or-phong/

    [5] Ribardière, Mickaël, Benjamin Bringier, Daniel Meneveaux, and Lionel Simonot, "STD: Student's t-Distribution of Slopes for Microfacet Based BSDFs," Computer Graphics Forum, vol. 36, no. 2, pp. 421-429, 2017. https://mribar03.bitbucket.io/projects/eg_2017/

    [6] Holzschuch, Nicolas, and Romain Pacanowski, "A Two-Scale Microfacet Reflectance Model Combining Reflection and Diffraction," ACM Transactions on Graphics (SIGGRAPH 2017), vol. 36, no. 4, pp. 66:1-66:12, July 2017. https://hal.inria.fr/hal-01515948

    [7] Hoffman N. Background: physics and math of shading[J]. Physically Based Shading in Theory and Practice, 2013

    [8] Cook, Robert L., and Kenneth E. Torrance, "A Reflectance Model for Computer Graphics," Computer Graphics (SIGGRAPH '81 Proceedings),vol. 15, no. 3, pp. 307-316, Aug. 1981. http://www.irisa.fr/prive/kadi/Lopez/p307-cook.pdf

    [9] Sébastien Lagarde , Physically Based Material Where Are We , SIGGRAPH 2017

    [10] Heitz E. Understanding the masking-shadowing function in microfacet-based BRDFs[J]. Journal of Computer Graphics Techniques, 2014

    [11] https://stackoverflow.com/questions/4414041/what-is-the-precision-of-highp-floats-in-glsl-es-2-0-for-iphone-ipod-touch-ipad

    [12] https://en.wikipedia.org/wiki/Floating_point\#Addition_and_subtraction

    [13] Dupuy, Jonathan, "Antialiasing Physically Based Shading with LEADR Mapping," SIGGRAPH Physically Based Shading in Theory and Practice course, Aug. 2014.

    https://blog.selfshadow.com/publications/s2014-shading-course/dupuy/s2014_pbs_leadr_slides.pdf

    [14] Beckmann, Petr, and André Spizzichino, "The Scattering of Electromagnetic Waves from Rough Surfaces", Pergamon Press, 1963. http://www.gbv.de/dms/ilmenau/toc/017290422.PDF

    [15] Blinn, James F., "Models of Light Reflection for Computer Synthesized Pictures," ACM Computer Graphics (SIGGRAPH '77 Proceedings), vol. 11, no. 2, pp. 192-198, July 1977.http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.131.7741&rep=rep1&type=pdf

    [16] Phong,BuiTuong, “Illumination for Computer Generated Pictures,” Communications of the ACM, vol. 18, no. 6, pp. 311–317, June 1975. Cited on p. 118, 340, 416 http://users.cs.northwestern.edu/\~ago820/cs395/Papers/Phong_1975.pdf

    [17] Trowbridge, T. S., and K. P. Reitz, "Average Irregularity Representation of a Roughened Surface for Ray Reflection," Journal of the Optical Society of America, vol. 65, no. 5, pp. 531-536, May 1975.

    [18] Walter, Bruce, Stephen R. Marschner, Hongsong Li, and Kenneth E. Torrance, "Microfacet Models for Refraction through Rough Surfaces," Rendering Techniques 2007, Eurographics Association, pp. 195-206, June 2007. https://www.cs.cornell.edu/\~srm/publications/EGSR07-btdf.pdf

    [19] Church, Eugene L., Peter Z. Takacs, and Thomas A. Leonard. "The prediction of BRDFs from surface profile measurements." Scatter from Optical Components. Vol. 1165. International Society for Optics and Photonics, 1990.

    [20] Bagher, Mahdi M., Cyril Soler, Nicolas Holzschuch, “Accurate Fitting of Measured Reflectances using a Shifted Gamma Microfacet Distribution,” Eurographics Symposium on Rendering (2012), 1509–1518, June 2012. http://hal.inria.fr/hal-00702304/en

    [21] Burley, Brent, "Physically Based Shading at Disney," SIGGRAPH Practical Physically Based Shading in Film and Game Production course, Aug. 2012.https://blog.selfshadow.com/publications/s2012-shading-course/burley/s2012_pbs_disney_brdf_notes_v3.pdf

    [22] https://www.arnoldrenderer.com/gallery/

    [23] https://scicomp.stackexchange.com/questions/5253/cancellation-problem-in-float-point-numbers

    [24] https://github.com/EpicGames/UnrealEngine

    [25] https://blender.stackexchange.com/questions/93790/is-the-roughness-setting-is-based-on-a-standard-measuring-unit

    [26] https://learnopengl.com/PBR/Theory

    [27] https://80.lv/articles/the-future-of-real-time-rendering-with-lumberyard/

    [28] https://assetstore.unity.com/packages/vfx/shaders/ctaa-cinematic-temporal-anti-aliasing-pc-vr-106995

    [29] http://advances.realtimerendering.com/s2014/epic/TemporalAA.pptx

    [30] M. Toksvig. “Mipmapping Normal Maps”. In: Journal of Graphics, GPU, and Game Tools 10.3 (2005), pp. 65–71. doi: 10.1080/2151237X.2005.10129203. http://www.nvidia.com/object/mipmapping_normal_maps.html

    [31] Olano, Marc, and Dan Baker, “LEAN Mapping,” in Proceedings of the 2010 ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games, ACM, pp. 181–188, 2010. Cited on p. 370 http://www.csee.umbc.edu/\~olano/papers/lean/

    [32] Baker, Dan, "Spectacular Specular--LEAN and CLEAN Specular Highlights," Game Developers Conference, Mar. 2011.http://twvideo01.ubm-us.net/o1/vault/gdc2011/slides/Dan_Baker_SpectacularSpecular.ppt

    [33] 题图来自《刺客信条:奥德赛》

     

    展开全文
  • 图形图像库集合

    千次阅读 2014-08-10 15:14:42
    2D图形库 AGG Google 图形处理引擎 skia 三维图形渲染引擎 OGRE 开源图形库 FreeImage 3D引擎 Irrlicht Engine 移动设备上的OpenGL OpenGL ES 高质量图形图表库 MathGL 开源图形库 CxImage
  • 对傅里叶函数以及级数的理解

    千次阅读 多人点赞 2018-04-09 21:29:50
    在完整的立体图中,我们将投影得到的时间差依次除以所在频率的周期,就得到了最下面的相位谱。所以,频谱是从侧面看,相位谱是从下面看。下次偷看女生裙底被发现的话,可以告诉她:“对不起,我只是想看看你的相位谱...
  • 1.图形与图像的区别图形是用数学方法来描述一幅图,强调图形的几何表示,图形是由场景的几何模型和景物的物理属性共同组成的。图像是把彩色图分成许许多多的象素,每个象素用若干个二进制位来指定该象素的颜色、亮度...
  • 数学 三角函数

    2012-04-23 22:34:00
    角 θ的所有三角函数 三角函数(Trigonometric)是数学中属于初等函数中的超越函数的一类函数。它们的本质是任意角的集合与一个比值的集合的变量之间的映射。通常的三角函数是在平面直角坐标系中定义的,其定义...
  • 【Matlab】图形用户界面设计

    万次阅读 多人点赞 2018-08-06 14:02:30
    专题八 MATLAB图形用户界面设计 一 图形窗口与坐标轴 1. 图形对象的句柄 1.1 句柄的概念 在MATLAB中,每一个具体的图形都是由若干个不同的图形对象组成的。 在MATLAB中,用句柄来标识对象,通过句柄来访问相应...
  • 5.2 二维纹理域的映射 在纹理映射技术中,最常见的纹理是二维纹理,通过映射将二维纹理变换到三维物体的表面,形 成最终的图像。 二维纹理可以通过函数和图像的方式表示。当使用函数的方式来表示的时候,使用将参数...
  • 2.1 常见的数学函数 398 2.1.1 求整数的绝对值 398 范例2-1 求整数的绝对值 398 ∷相关函数:abs函数 2.1.2 求长整型整数的绝对值 399 范例2-2 求长整型整数的绝对值 399 ∷相关函数:labs函数 2.1.2 求...
  • 立体渲染

    千次阅读 2015-10-13 19:49:56
    立体渲染 外文名 volume rendering 又 称 立体绘制 实 质 二维投影的技术 目录 1 简介 2 直接立体渲染 ▪ 立体光线投射 ▪ Splatting ▪ Shear Warp ▪ 纹理映射 ...
  • 计算机图形学大作业 课程设计 实验报告 消隐算法安徽建筑工业学院计算机图形学 大作业大作业名称: 消隐算法 演示院(系)名称:专 业:班 级:姓 名:学 号:指 导 老 师:2011 ~ 2012年度第 一 学期计算机图形学——...
  • 计算机图形学——OpenGL学习系列之Graphics3D 一、OpenGl中的坐标系 跟数学中常见的坐标系有点不同,Z轴垂直纸面,反正我刚开始是不习惯的 二、OpenGl中的几何变换 在OpenGl中,无论2D还是3D都可以进行几何变换,...
  • 本文主要介绍对极几何(Epipolar Geometry)与立体视觉(Stereo Vision)的相关知识。对极几何简单点来说,其目的就是描述是两幅视图之间的内部对应关系,用来对立体视觉进行建模,实际上就是一种约束条件,这样可以确定...
  • 图形图像处理

    千次阅读 2011-07-14 13:58:58
    形图像处理起源于20世纪20年代,当时通过海底电缆从英国...此后由于遥感等领域的应用,使得图形图像处理技术逐步得到发展。一直到20世纪50年代,随着大型数字计算机和太空科学研究计划的出现,人们才注意到图像处理的潜
  • 一元、二元函数图像绘制

    千次阅读 2015-09-09 13:42:00
    本篇博客主要是在上一篇《每个人都该懂点函数式编程》的基础上,进一步说明“函数”在函数式编程中的重要作用。强调了函数和普通类型一样,可以赋值、存储、传参以及作为另外函数的返回值。 本文附带了一个Demo,该...
  • Mathematica函数大全

    千次阅读 2013-10-31 10:34:25
    Normal[expr] 化简并给出最常见的表达式 SeriesCoefficient[series, n] 给出级数中第n 次项的系数 SeriesCoefficient[series, {n1,n2...}] '或Derivative[n1,n2...][f] 一阶导数 InverseSeries[s, x] 给出逆...
  • 高中数学在高考过程中以立体几何为背景的新颖问题常见的有折叠问题,与函数图像相结合问题、最值问题、探索性问题等,对探索性、开放性、存在型问题的考察,探索性试题使问题具有不确定性、探究性和开放性,对学生的...
  • 安装各种库东奔西走......  GLUT(英文全写:OpenGLUtilityToolkit)是一个处理OpenGL程式的...以一个函数呼叫绘制某些常用的立体图形,例如长方体、球、以及犹他茶壶(实心或只有骨架,如glutWireTeapot()) 提...
  • 什么是计算机图形学?

    千次阅读 2016-06-06 12:49:54
    什么是计算机图形学? 刘利刚 中国科技大学 http://staff.ustc.edu.cn/~lgliu   【注】 由于时常有本科学生来向笔者询问计算机图形学是做什么的,为了使得学生能够快速了解计算机图形学,有利于他们在选择...
  • 清华大学计算机图形学课程

    千次阅读 2016-06-06 13:00:35
    1.1 计算机图形学的研究内容 1.2 发展的历史回顾 1.3 应用及研究前沿 1.4 图形设备 2学时 第二章 颜色模型、图像基本知识、Phong光照模型 2.1 颜色模型  2.1.1 颜色模型的视觉基础  ...
  • 第一讲 计算机图形学概论第一周测验题1、显示颜色64K,分辨率为1024*1024的显示器,至少需要的帧缓存容量为A、2MBB、1MBC、3MBD、512KB2、在下列有关显示器的叙述中,正确的论述为A、CRT显示器已经被淘汰B、笔记本...
  • 图形学数学基础之采样分布映射

    千次阅读 2018-02-03 17:39:57
    接下来我们就尝试推导几个常见的PDF,以后也会慢慢添加其他用到的。 这里向大家推荐一个求解积分函数的工具,我下面的积分求解都是通过这个工具来实现的: WolframAlpha Perfect Diffuse 1. p ...
  • 图形文件的格式

    千次阅读 2007-07-20 15:18:00
    近年来,个人计算机和工作站上的图形工具比几年前的巨型机上的图形工具还要多,计算机图形学的领域也随之扩展。过去,当人们编出越来越多的图形应用程序后,需要把图像文件存储下来以作日后的处理或显示之用。在缺乏...
  • 图形加速卡技术介绍

    千次阅读 2012-03-16 12:34:08
    图形加速卡技术论坛:1.入门篇--图形加速之 图形 (发表于GZeasy.com: Jul 20 2003, 04:14 PM) 来了这么久,也对这里的朋友有了一个大概的了解。 恕我罗索一两句,这里的名字是“图形加速卡技术论坛”,...
  • 计算机图形学——光线追踪

    千次阅读 2017-06-15 14:19:52
    计算机图形学——光线追踪  参考资料 :用JavaScript玩转计算机图形学  实验内容:使用光线追踪进行场景渲染。  实验效果:   1、简介  光线跟踪(ray tracing)是一个在二维(2D)屏幕上呈现三维(3D)图像的方法。...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 4,000
精华内容 1,600
关键字:

常见立体图形函数