Houdini学习笔记035_Poincaré Disc（庞加莱圆盘）<上>

先说说我为什么会学这个，起因是下面这幅画。这是我崇拜的荷兰艺术家Escher的代表作之一——Circle Limit IV。Escher的画作中充满了数学计算和空间几何的概念，很难想象不借助计算机单靠尺规作图能完成这样的作品。

于是我开始搜索关于这幅作品的几何知识，发现这个圆叫做Poincaré Disc（庞加莱圆盘），它属于一种罗氏几何模型。我们熟知的几何是欧式几何，而罗氏几何是非欧几何。在这一几何体系中，“第五公设”的表述为：过直线外一点至少有两条直线和已知直线平行。

是不是有点懵？

与罗氏几何对应的还有黎曼几何，它直接认为过直线外一点不存在已知直线的平行线。

其实，这三种几何都是成立的，只是适应的体系不同。我们熟悉的欧式几何反映的是平直空间的几何关系，而罗氏几何和黎曼几何分别适用于负曲率空间和正曲率空间。爱因斯坦的广义相对论就是受到黎曼几何的启发而提出的。

简单来说，我们可以把庞加莱圆盘看作一个异空间，在这个空间中，圆心处距离为0，圆周上的点相当于无穷远。而连接圆内两点（非圆心）的直线用圆弧来表示，称之为hyperbolic line。

花了好几天的时间终于厘清了其中的数学关系，甚至找到了几篇相关的文献，比如下面这篇Dunham教授1986年发表的文章。文献中给出了类似上述Circle Limit IV图案的画法，分为两个基本步骤，先画出一个单元图案，然后利用在空间中的变换关系拷贝得到其他的图案。

但是文献中给出的是变换矩阵迭代的方法，本来就是非欧几何空间，还要不断计算更新变换矩阵，我果断放弃了。然后搜了下Horikawa小哥的教程，居然有这方面的介绍（点这里），然后就找到了这个网站http://www.malinc.se。

又花了一个周末的时间，费了一大堆稿纸，确定文献中给的步骤是合理的。但要在Houdini中实现，我们需要借助一个工具——Mobius Inverse（莫比乌斯变换）。具体见下图，对于圆内一点A，过A点作垂直于OA的弦，并与圆相交于H点。过H点作圆的切线，与OA的延长线交于M点。由相似三角形关系可知：OA · OM = OH^2 = r^2。M就是A点相对于圆O的莫比乌斯变换点。

今天我们先学习第一步，初始单元图案的制作。这里需要定义两个参数：p和q。p是多边形的边数，q是交点的连线数。比如下图中，p = 4意为每个图形都是四边形，q = 5意为每个交点处连接有五条边（不经过圆心的边在欧式几何中<即我们在屏幕上看到的>都是弧线）。

如果初始图案的中心点位于圆心，绘制起来相对简单一些。但是为了赋予其更多变化，我们让它可以是圆内的任意一点。问题是，当中心点移动时，初始图案并不是在欧式几何平面内那样简单的移动旋转，我们需要在新的中心点和圆心之间找到垂直平分线（也是一条弧线），然后进行莫比乌斯变换。

中垂线可由下图得到，对于圆内任一点B（除圆心O外），M为OB延长线上一点（圆外），A是OB在庞加莱圆盘中的中点，注意这里看上去两端不相等，因为越靠近圆周，同样长度的线段表示的距离就越大（www.malinc.se这里有具体的计算方法）。对于圆M来说，O和B满足莫比乌斯反演变换关系。圆M在圆O内的弧线就是OB的垂直平分线。

创建一个null节点（取名“CONTROLLER”），自定义如下参数：p和q为整数，需满足(p-2)(q-2) > 4才有解，取值范围可分别设为3~10。r是初始图案中心点B距圆心O的距离，范围为0~1。theta是B点绕圆心的角度，其坐标可表示为(r*cos(theta), r*sin(theta), 0)。rotation是初始图案绕中心点B旋转的角度，默认为0。最后，segments是画弧线取的分段数，默认设为10。

创建一个Attribute Wrangle节点，用于初始图案的绘制，Run Over方式改为Detail (only once)。全部代码如下，我们逐一来解释：

首先定义了两个函数，一个是mobiusInverse（莫比乌斯变换），给定圆心o和半径r，对于点m求其变换后的坐标，返回的是一个矢量。因为所有的图案都是在xy平面内，z方向上的坐标值为0，计算时可忽略。

第二个自定义函数getCircle是根据三个点坐标求外接圆圆心，网上套用的计算方法，最后返回圆心坐标center和半径radius。

单元图案的半径d用如下公式计算，跟p和q的值有关。由

如果初始图案的中心位于圆心，假设A点坐标为(d,0,0)，第二个点B点坐标为(d*cos(θ), d*sin(θ), 0)。这里的θ = 2π/p。另外一个点A'的坐标为(1/d,0,0)。点A、A'、B三点在一个圆上，用getCircle函数找出这个圆的圆心和半径。AB两点之间的弧线就是初始图案的其中一条边。绕O点旋转p-1次就可以得到完整的图案。

如果初始图案的中心不在O点，而是在C点（由自定义参数 r 和 theta确定）。下面的代码就是计算OC的中垂线（弧）对应的圆心cm和半径rm。注意是r != 0时，如果r=0，表示初始图案中心就是圆心O，无需进行莫比乌斯变换。

再后面是创建弧线的数学计算部分，用了两层for循环，内层循环是在一段弧上挨个生成点，最后连成弧线；外层循环是通过旋转分别得到其他的弧线。所以外层循环共p次，内层循环segments+1次。

这里大家不熟悉的一个是ident()，一个是rotate函数。前者定义一个标准矩阵，后者是对矩阵进行旋转的函数（得到的是一个旋转矩阵）。变量4 = ZAXIS，也可以用{0,0,1}表示。

坐标乘上旋转矩阵就会得到旋转后的坐标：

看下最终的结果——

可以看到，当图案远离圆心时会越来越小。但是在罗氏空间中，其半径和每边的边长都是不变的。后面我们会对这个初始形状不断进行莫比乌斯变换，充满整个庞加莱圆盘。今天的内容到此为止，不指望有几个人能完全理解，看看就好。