手把手教你用Python画一个绝美土星环

迈不过友情╰ 2023-03-12 12:27 51阅读 0赞

![format_png][]

**导读：**本文带你了解 Python 图形工具能够做哪些事情。

作者：B. J. 科里特斯

来源：大数据DT（ID：hzdashuju）

![format_png 1][]

土星的行星环非常出名。虽然木星、土星、天王星和海王星也有环，**但土星环是我们太阳系中最大、最亮、最广为人知的行星环。**

它由小到灰尘的颗粒，大到巨石的物体组成。这些物体的成分主要是冰，一般认为是彗星或较大的小行星与土星的一颗卫星相撞时产生的，两者都撞成了小碎块。在远古时代，土星就已为人所知，但直到1610年，伽利略才首次用望远镜对它进行观测。

这个行星以罗马农业之神土星Saturn命名为，也就是我们每个星期的第六天Saturday。

**图1至图5中的图像是由本文文末的代码生成的。每张图都呈现不同的方向角**，图标题中有相应的说明。

图标题中还列出了入射光线的单位矢量分量，例如lx=+0.707，ly=+0.707，lz=0 表示左上象限中的光源；lx=-1，ly=0，lz=0表示来自右侧的光源。在图像中请注意行星在环上投射的阴影，尤其是在图5中能够看到行星轮廓的曲率。

![format_png 2][]

▲图1 包含土星环和阴影的土星1：Rx=-20°, Ry=0, Rz=-10°, lx=1, ly=0, lz=0

![format_png 3][]

▲图2 包含土星环和阴影的土星2：Rx=-8°, Ry=0, Rz=-30°, lx=0.707, ly=.707, lz=0

![format_png 4][]

▲图3 包含土星环和阴影的土星3：Rx=20°, Ry=0, Rz=25°, lx=0.707, ly=0.707, lz=0

![format_png 5][]

▲图4 包含土星环和阴影的土星4：Rx=-10°, Ry=0, Rz=25°, lx=0.707, ly=-0.707, lz=0

![format_png 6][]

▲图5 包含土星环和阴影的土星5：Rx=20°, Ry=0, Rz=30°, lx=-1, ly=0, lz=0

为了进行比较，我们再看一张土星的摄像图：

![format_png 7][]

©NASA/JPL-Caltech/Space Science Institute/G. Ugarkovic (ISS), NASA/JPL-Caltech/University of Arizona/CNRS/LPG-Nantes (VIMS)

你可以在下面网址找到更多土星的摄像图：

https://www.jpl.nasa.gov/spaceimages/?search=saturn&category=\#submit

图6所示为构建土星环所用的数学模型。这里介绍一种实现**球体着色**的算法。首先创建一个直立球体，也就是说，经度是垂直的，纬度是水平的（即平行于XZ平面），然后从初始方向开始，围绕x，y和z轴对球体进行旋转。

![format_png 8][]

▲图6 土星环模型：行星俯视图和从XZ平面上Rx=0, Ry=0, Rz=0向下看环

我们对土星环也进行同样的操作。我们可以创建平行于XZ平面的水平环，然后将它和土星一起旋转相同的角度。土星环所处的平面穿过土星的球心，因此土星和环具有相同的旋转中心。

土星环绘制为一系列相邻的同心圆，每个同心圆由短线段组成。参考图6和文末代码，程序第42和43行设置了土星环的内半径和外半径，第44行设置同心圆的间距。土星环被分成七个同心环形带（图6中未画出）且具有不同的颜色，第45行的deltar是它们的宽度。

构成同心圆的每个线段都单独绘制。第48行从r1向r2进行绘制，通过径向循环绘制圆弧段。第49行是绕圆周方向绘制的循环。第50-61行执行旋转操作产生第62和63行中的全局绘图坐标xpg和ypg，旋转函数与先前程序中的相同。

接下来在第66-75行中设置线段的颜色。**土星环是由不同颜色的条带构成的，这和NASA观测图像中看到的物理组成结果一致。**从r = r1到r1 + deltar的第一个条带具有颜色clr=(.63，.54，.18)，剩余的条带也是如此。

第五个条带省略掉了，因为它是空的，背景颜色能显示出来。第六个条带的宽度是其他条带的两倍，并且为第七个条带提供了颜色。

对于给定的光方向，从大多数角度上，行星体本身都会在环上投下阴影。参考图7，我们的目标是确定点p到底位于行星阴影区域内部还是外部。

球状的行星将产生圆形的阴影，阴影的直径与行星的尺寸相等，或者更准确地说，是球体的“大圆”。它是用**通过圆心的平面切割球体而得到的最大圆**，就像把橙子切成两半，你看到的是一个橙子的最大圆。

在图7中，这种阴影可能是由相同大小的圆盘投影产生的，也可能是由球状行星投影产生的，两种情况下，阴影的大小都是一样的。在土星的侧面图中，大圆显示为是一条通过平面中心的加粗线。

从图7的几何图形中可以看出，如果p位于|B| > rs的位置，则它位于阴影区之外，其中rs是土星的半径；如果|B| < rs，则p位于在阴影区之中。在绘制条带的时候，如果我们确定了p的位置在阴影区中，我们就把这个点涂成灰色，如果它在阴影区之外，我们就用第66-75行设置的条带颜色给它着色。

![format_png 9][]

▲图7 阴影模型

我们的目标是求出给定位置p时的|B|值，由图7可看出：

|B|=|V|sin(φ)

并且我们知道：

V×û=|V||û|sin(φ)

其中û=-î 。将上述方程与|û|=1合并得：

B=V×û

|B|=|V×û|

在代码第78行中确定了入射光矢量î 的长度为1，但如果在第23-25行中输入的分量不计算为1（即![format_png 10][]），则入射光矢量可能不等于1。有必要的话，需要在第79-81行进行重构。第82-84行建立矢量V的分量。第85-87行计算B的分量。第88行给出其幅度magB = |B|。

第89行确定p是否位于阴影区内，如果是，则执行第90行V与î的点积。这是用来确定p是否位于行星朝向光源的一侧，在这种情况下，它与行星的暗侧相对，而不在阴影区。因为在第78-89行的阴影算法中并未区分p的位置，所以此处必须进行明确。如果p确实位于阴影区域内的暗侧，则在第91行中将p设置为中等灰色。

相信你一定注意到，图6的土星环上有一个暗色条带，**这是因为土星环在该位置上是空的：那里没有固体颗粒物，不能反射阳光。**因此我们透过条带看到了背景颜色'midnightblue'。但这会产生一个问题，即阴影颜色的绘制会覆盖该空白处的背景颜色，因此在第93和94行将其重置为'midnightblue'。

既然条带的颜色都建立好了，我们就可以通过一个个短线段来绘制土星环了。第97-100行计算第一个线段的起始位置。参考图6，第100-101行确定该线段与土星的遮挡关系，线段在前就会被绘制。

103-108行确定线段是否在土星后面，位于后面就不会被绘制。这种遮挡关系是通过计算点的全局坐标与土星中心之间的距离c来完成的。

第107行的意思是，如果c大于球体半径的1.075倍，则绘制该段。因子1.075的作用是防止线段与球体的表面重合，这是有必要的，不然小于球体半径的可见区段将不会被绘制。

本文代码生成的图像有两点需要注意：首先是颜色。美国宇航局的摄影图像呈现的是一种几乎没有颜色灰色，但是**许多土星观察者都将它描述为金黄色，因此我们选择了金色。**

所有摄影师都知道，在摄影图像中呈现物体的真实颜色是十分困难的，颜色取决于入射光和物体本身的颜色，或许最好的方法是依靠肉眼观察。

**如果你不赞同本文代码所生成的图像中的颜色，可以通过更改程序中clr的定义来修改它们。**需要注意的第二点，是图5中土星表面阴影的曲率，它表示着色算法是否按预期工作。

在程序的使用上，你可以自行更改第24-26行中的入射光的方向和第32-34行中的旋转角度。

本文篇幅有限，更多更详细的讲解请参阅《Python图形编程：2D和3D图像的创建》一书。

*  **土星代码**

运行代码也需要一段时间，请耐心等待。

1"""
      2SATURN
      3"""
      4
      5import numpy as np
      6import matplotlib.pyplot as plt
      7from math import sin, cos, radians, sqrt
      8
      9plt.axis([0,150,100,0])
     10plt.axis('off')
     11plt.grid(False)
     12
     13print('running')
     14#—————————————————parameters
     15g=[0]*3
     16
     17xc=80 #———sphere center
     18yc=50
     19zc=0
     20
     21rs=25 #———sphere radius
     22
     23lx=-1 #———light ray unit vector components
     24ly=0
     25lz=0
     26
     27IA=0
     28IB=.8
     29+n=2
     30
     31Rx=radians(-20)
     32Ry=radians(0)
     33Rz=radians(30)
     34
     35#————————same as SHADESPHERE—————–
     36
     37#———————————————————rings
     38alpha1=radians(-10)
     39alpha2=radians(370)
     40dalpha=radians(.5)
     41
     42r1=rs*1.5
     43r2=rs*2.2
     44dr=rs*.02
     45deltar=(r2-r1)/7 #———ring band width
     46
     47#—————————————rotate ring point p which is at r, alpha
     48for r in np.arange(r1,r2,dr):
     49    for alpha in np.arange(alpha1,alpha2,dalpha):
     50        xp=r*cos(alpha)
     51        yp=0
     52        zp=-r*sin(alpha)
     53        rotx(xc,yc,zc,xp,yp,zp,Rx)
     54        xp=g[0]-xc
     55        yp=g[1]-yc
     56        zp=g[2]-zc
     57        roty(xc,yc,zc,xp,yp,zp,Ry)
     58        xp=g[0]-xc
     59        yp=g[1]-yc
     60        zp=g[2]-zc
     61        rotz(xc,yc,zc,xp,yp,zp,Rz)
     62        xpg=g[0]
     63        ypg=g[1]
     64
     65#—————————————————select ring band color
     66    if r1 <= r < r1+1*deltar:
     67        clr=(.63,.54,.18)
     68    if r1+1*deltar <= r <= r1+2*deltar:
     69        clr=(.78,.7,.1)
     70    if r1+2*deltar <= r <= r1+3*deltar:
     71        clr=(.95,.85,.1)
     72    if r1+3*deltar <= r <= r1+4*deltar:
     73        clr=(.87,.8,.1)
     74    if r1+5*deltar <= r <= r1+7*deltar:
     75        clr=(.7,.6,.2)
     76
     77#———————————————————————shadow
     78    magu=sqrt(lx*lx+ly*ly+lz*lz)
     79    ux=-lx/magu
     80    uy=-ly/magu
     81    uz=-lz/magu
     82    vx=xc-xpg
     83    vy=yc-ypg
     84    vz=zc-zpg
     85    Bx=uy*vz-uz*vy
     86    By=uz*vx-ux*vz
     87    Bz=ux*vy-uy*vx
     88    magB=sqrt(Bx*Bx+By*By+Bz*Bz)
     89    if magB < rs: #—————————if in the shadow region
     90        if vx*lx+vy*ly+vz*lz <= 0: #———if v points toward light source
     91            clr=(.5,.5,.2) #———shadow color
     92
     93    if r1+4*deltar <= r <= r1+5*deltar: #———overplot empty band
     94        clr='midnightblue' #———with background color
     95
     96#——————————————————–plot line segment
     97    if alpha == alpha1:
     98        xstart=xpg
     99        ystart=ypg
    100    if zpg <= zc: #–front (z axis points into the screen)
    101        plt.plot([xstart,xpg],[ystart,ypg],linewidth=2,color=clr)
    102
    103    if zpg >= zc: #–back
    104        a=xpg-xc
    105        b=ypg-yc
    106        c=sqrt(a*a+b*b)
    107        if c > rs*1.075: #——plot only the visible portion of rings
    108            plt.plot([xstart,xpg],[ystart,ypg],linewidth=2,color=clr)
    109        xstart=xpg
    110        ystart=ypg
    111
    112plt.show()

关于作者：B. J. 科里特斯，博士，一直从事计算机工程和科学应用。他是一名教育家、顾问，著有十多本关于几何建模、计算机图形学、人工智能、物理过程模拟、结构分析、计算机在科学和工程中的应用的书籍。

本文摘编自《Python图形编程：2D和3D图像的创建》，经出版方授权发布。

延伸阅读《Python图形编程》

点击上图了解及购买

转载请联系微信：DoctorData

**推荐语：**一本书教你学会Python图形图像编程及可视化。本书使用Python的内置特征创建用于数据可视化和技术说明的创新图形。从构建绘图空间和生成2维和3维对象的基本Python函数开始，你将学习如何构造更复杂的对象、转换和旋转、删除隐藏线、引入阴影以添加真实感以及将图像投影到任何数据集上。

![音符][format_png 11]

**划重点????**

**干货直达????**

*  [大厂和初创公司都在用！Redis好在哪？终于有人讲明白了][Redis]  
    
 *  [5种不会被机器智能替代的能力][5]  
    
 *  [世界一流大学的计算机专业，在用哪些书当教材？][Link 1]  
    
 *  [什么是机器阅读理解？跟自然语言处理有什么关系？][Link 2]

**更多精彩????**

在公众号对话框输入以下**关键词**

查看更多优质内容！

**PPT** | **读书** | **书单** | **硬核** | **干货**** **

**大数据** | **揭秘** | **Python** | **可视化**

**AI** | **人工智能** | **5G** | **中台**

**机器学习** | **深度学习** | **神经网络**

**合伙人**** **| **1024** | **大神** | **数学**

据统计，99%的大咖都完成了这个神操作

**????**

![format_png 12][]

![format_png 13][]

[format_png]: /images/20230312/72461cb58405447ca44bf472b3fc239f.png
[format_png 1]: /images/20230312/ab933ee632df4d1db09dfd1f05dad5c2.png
[format_png 2]: /images/20230312/6112c5c53cf846d597f44a7913cddd35.png
[format_png 3]: /images/20230312/af89bc96cd8948ecaece954c18e3f2ce.png
[format_png 4]: /images/20230312/921523d2de98463fa5738fbcc6a69c6d.png
[format_png 5]: /images/20230312/f1b4397e50fa41539d041590fc261b83.png
[format_png 6]: /images/20230312/fbe02653526c4d798f8dd515962d9191.png
[format_png 7]: /images/20230312/7e43dacaf0534887ae5b39e93a4b80b4.png
[format_png 8]: /images/20230312/121e2b42fbff4abdbde386d5bd84941d.png
[format_png 9]: /images/20230312/d85fee069b654fa8b6e9f2915064b986.png
[format_png 10]: /images/20230312/10c5f196bd314ae8b0867e01f8d87b87.png
[format_png 11]: /images/20230312/111279c646944ce298a84c48be08e5c1.png
[Redis]: http://mp.weixin.qq.com/s?__biz=MzI5OTk5OTM2Mw%3D%3D&chksm=ec8c8541dbfb0c5705308c740d81e2a258c46f4d829f786a97f5ed341f852374f623eb1335a6&idx=2&mid=2247503021&scene=21&sn=93f76cf827fd7c8b730e3df73a24364b#wechat_redirect
[5]: http://mp.weixin.qq.com/s?__biz=MzI5OTk5OTM2Mw%3D%3D&chksm=ec8c856bdbfb0c7d2f3bfdb2198d142f20684d9d81ef79e789a17013b673ae569f765e07e219&idx=2&mid=2247502983&scene=21&sn=37842690429f448692642f60daeb8c77#wechat_redirect
[Link 1]: http://mp.weixin.qq.com/s?__biz=MzI5OTk5OTM2Mw%3D%3D&chksm=ec8c8568dbfb0c7e8fbc8368d448018fe5a64c05f038e7d3e6f4cb3cd2eb78d78975e6b17f48&idx=1&mid=2247502980&scene=21&sn=7980d5a1bf337f12eda1c49f1782eda7#wechat_redirect
[Link 2]: http://mp.weixin.qq.com/s?__biz=MzI5OTk5OTM2Mw%3D%3D&chksm=ec8c85a0dbfb0cb6a1ae5ed37cbab88078a510e5852a0ba89c22312bfa8aecf45e093a00c2e9&idx=2&mid=2247502924&scene=21&sn=757324c685fe02c84a2e91fd988629af#wechat_redirect
[format_png 12]: /images/20230312/b474e2abf033498c9d41abb7e5038671.png
[format_png 13]: /images/20230312/82559bb7c3384b7794cc1af535a1629e.png