用C语言画颗心（一）——心形曲面-Charlie小站

最近在知乎上看到一篇几年前帖子，有大神分享如何用C语言画心形，图就是上面那张。大神的程序狂拽炫酷，仔细研读，记录于此。

在wolfram上搜索心形曲面(heart surface)，可以得到心形曲面的笛卡尔方程为$f(x,y,z)=0$，其中
$$f(x,y,z)=\left( x^2 + \frac{9y^2}{4}+z^2-1\right)^3-x^2z^3-\frac{9y^2z^3}{80}$$
这个曲面长成以下的样子

上图给出了$x,y,z$的取值范围，省了不少事。自己去找边界可就费死劲了。

我尝试通过纯粹的数学手段去计算边界，试了一天，结果屁都没搞出来。原来想的通过法向量的方向确定一个边界，就是下面这个3个方程：
$$
\frac{\partial f}{\partial y}(x,y,z) = 0 ,
\frac{\partial f}{\partial z}(x,y,z) = 0 ,
f(x,y,z) = 0
$$
三个未知数，三个方程，高次，不太好搞。于是撸了个牛顿法（代码在这里）。解出来发现不对，怎么查都感觉没错，但是结果始终在（1,0,0）附近，后来发现$\nabla f(1,0,0) = (0,0,0)$！果然，迭代法初值很重要。但我实在是无力去试各种初值了，就此作罢。

接下来还要算一下曲面的法向量，因为在后期计算亮度时需要用到入射方向与曲面表面法向量的正弦值。实际上，曲面的法向量非常容易计算：
$$ \vec{n} = \left[\frac{\partial f}{\partial x}, \frac{\partial f}{\partial y},\frac{\partial f}{\partial z}\right]$$
还记得等势面与电场线垂直的性质以及如何通过电势计算电场强度吗，跟上面这个式子是一样的。不过，这里坑也是有的，就像上面那样，有些地方的梯度为零。

接下来可以开始考虑画图的步骤，先看一下二维图

心形曲面是$f(x,y,z)$的0等值面，若$f(x,y,z)\le 0$，则表示该点位于等值面内部，反之亦然。定义函数$f$：

double f(double x, double y, double z){
	double a = x*x + 9.0f/4.0f*y*y + z*z - 1;
	return a*a*a - x*x*z*z*z - 9.0/80.0*y*y*z;
}

着色或者是光照的问题下次再讨论，这里先直接用字符画个平面的心形。逻辑很简单：设定宽和高之后，遍历每个“像素”，判断该像素是否在心内，若是，打印字符，若否，打印另一种字符。程序如下：

int main(int argc, char const *argv[])
{
	for (float z = 1.5; z > -1.5; z-=0.1)
	{
		for (float x = -1.5; x < 1.5; x+=0.05)
		{
			if (f(x,0,z) < 0)
			{
				putchar('*');
			}
			else
				putchar(' ');
		}
		putchar('\n');
	}
	return 0;
}

完整代码路径戳这里，最后再丢张效果图