故宫下雪了！我花了45秒，用Python给它画了一组手绘图

欢迎关注天善智能，我们是专注于商业智能BI，人工智能AI，大数据分析与挖掘领域的垂直社区，学习，问答、求职一站式搞定！

对商业智能BI、大数据分析挖掘、机器学习，python，R等数据领域感兴趣的同学加微信：tstoutiao，邀请你进入数据爱好者交流群，数据爱好者们都在这儿。

作者：丁彦军 

知乎专栏|恋习Python：

https://zhuanlan.zhihu.com/p/35667053

个人公众号：恋习Python

历史文章推荐：

我用Python，3分钟快速实现，9种经典排序算法的可视化

这是恋习Python之手把手系列第10篇原创首发文章

这几天，许多城市，迎来了2019年的第一场雪

13日早晨，当北京市民拉开窗帘时发现，窗外雪花纷纷扬扬在空中飘落

而且越下越大，树上、草地、屋顶、道路上...都落满雪花

京城银装素裹，这是今冬以来北京迎来的第二场降雪

一下雪，北京就变成了北平，故宫就变成了紫禁城

八万张门票在雪花飘下来之前，便早已预订一空

一、概念与原理

我们都知道手绘图效果的特征主要有：

黑白灰色；边界线条较重；相同或相近色彩趋于白色；略有光源效果

核心原理：利用像素之间的梯度值和虚拟深度值对图像进行重构，根据灰度变化来模拟人类视觉的模拟程度

把图像看成二维离散函数，灰度梯度其实就是这个二维离散函数的求导，用差分代替微分，求取图像的灰度梯度。常用的一些灰度梯度模板有：Roberts 梯度、Sobel 梯度、Prewitt 梯度、Laplacian 梯度。

以Sobel 梯度计算来解释：

二、图像的数组形式与变换

其中，需要用到的方法：

Image.open( ): 打开图片

np.array( ) : 将图像转化为数组

convert("L"): 将图片转换成二维灰度图片

Image.fromarray( ): 将数组还原成图像uint8格式

代码如下：

from PIL import Image
import numpy as np

im = Image.open(r"C:\Users\Administrator\Desktop\gugong\微信图片_20190216152248.jpg").convert('L')
a=np.asarray(im).astype('float')
print(a.shape,a.dtype)
(1080, 608) float64
#(1080, 608)分别表示高度，宽度

三、图像的手绘效果处理

实现思路步骤：

1、梯度的重构

numpy的梯度函数的介绍

np.gradient(a) ： 计算数组a中元素的梯度，f为多维时，返回每个维度的梯度 

离散梯度： xy坐标轴连续三个x轴坐标对应的y轴值：a, b, c 其中b的梯度是（c-a）/2 

而c的梯度是： (c-b)/1

当为二维数组时，np.gradient(a) 得出两个数组，第一个数组对应最外层维度的梯度，第二个数组对应第二层维度的梯度。 

代码如下：

grad=np.gradient(a)
grad_x,grad_y=grad
grad_x = grad_x * depth / 100.#对grad_x值进行归一化
grad_y = grad_y * depth / 100.#对grad_y值进行归一化

2、构造guan光源效果

设计一个位于图像斜上方的虚拟光源

光源相对于图像的视角为Elevation,方位角为Azimuth

建立光源对各点梯度值的影响函数

运算出各点的新像素值

其中：

np.cos(evc.el) : 单位光线在地平面上的投射长度

dx,dy,dz ：光源对x，y，z三方向的影响程度

3、梯度归一化

构造x和y轴梯度的三维归一化单位坐标系；

梯度与光源相互作用，将梯度转化为灰度。

4、图像生成

具体详情代码如下：

from PIL import Image
import numpy as np
import os
import join
import time

def image(sta,end,depths=10):
    a = np.asarray(Image.open(sta).convert('L')).astype('float')
    depth = depths  # 深度的取值范围(0-100)，标准取10
    grad = np.gradient(a)  # 取图像灰度的梯度值
    grad_x, grad_y = grad  # 分别取横纵图像梯度值
    grad_x = grad_x * depth / 100.#对grad_x值进行归一化
    grad_y = grad_y * depth / 100.#对grad_y值进行归一化
    A = np.sqrt(grad_x ** 2 + grad_y ** 2 + 1.)
    uni_x = grad_x / A
    uni_y = grad_y / A
    uni_z = 1. / A
    vec_el = np.pi / 2.2  # 光源的俯视角度，弧度值
    vec_az = np.pi / 4.  # 光源的方位角度，弧度值
    dx = np.cos(vec_el) * np.cos(vec_az)  # 光源对x 轴的影响
    dy = np.cos(vec_el) * np.sin(vec_az)  # 光源对y 轴的影响
    dz = np.sin(vec_el)  # 光源对z 轴的影响
    b = 255 * (dx * uni_x + dy * uni_y + dz * uni_z)  # 光源归一化
    b = b.clip(0, 255)
    im = Image.fromarray(b.astype('uint8'))  # 重构图像
    im.save(end)

def main():
    xs=10
    start_time = time.clock()
    startss = os.listdir(r"C:\Users\Administrator\Desktop\gugong")
    time.sleep(2)
    for starts in startss:
        start = ''.join(starts)
        sta = 'C:/Users/Administrator/Desktop/gugong/' + start
        end = 'C:/Users/Administrator/Desktop/gugong/' + 'HD_' + start
        image(sta=sta,end=end,depths=xs)

    end_time = time.clock()
    print('程序运行了  ----' + str(end_time - start_time) + '   秒')
    time.sleep(3)

main()
程序运行了  ----43.01828205879955   秒  #一共35张图片

最终效果图对比：

其他图片就不一一列举，若需要获取更多图片素材，关注公众号后台回复故宫雪景即可获得；你也可以通过此代码为自己画一张手绘图；也可以为自己的家乡或母校画。

参考资料：北京理工大学的嵩天老师的网络课程

http://www.icourse163.org/learn/BIT-1001870002?tid=1001963001#/learn/announce

最后，恋习Python温馨提示：

人生苦短，我用Python；

除了生娃，啥都能干！！

欢迎大家在留言处，留言自己曾经用Python做过哪些有意思的事！

Python的爱好者社区历史文章大合集：

2018年Python爱好者社区历史文章合集（作者篇）

2018年Python爱好者社区历史文章合集（类型篇）

关注后在公众号内回复“ 课程 ”即可获取：

小编的转行入职数据科学（数据分析挖掘/机器学习方向）【最新免费】

小编的Python的入门免费视频课程！

小编的Python的快速上手matplotlib可视化库！

崔老师爬虫实战案例免费学习视频。

陈老师数据分析报告扩展制作免费学习视频。

玩转大数据分析！Spark2.X + Python精华实战课程免费学习视频。