光线追踪算法综述
Problem Formulation
Ray Tracing的目标是生成一张包含场景内物体，具有真实感的图像，因此实现一个简单的Ray Tracing算法并不需要显式地构建一个可视的三维场景，只需要隐式地构建三维空间就可以了（也就是说这个三维场景只要存在你的脑袋里就可以了）。生成包含酒杯的渲染图像并不是一件很简单的事情，但是只生成包含几个简单几何体的渲染图，只需要两三百行代码。不需要图形库函数，只需要最基本的STL库。
Ray Tracing能够实现一些使得画面更具真实感的效果，包括阴影、折射和反射，这些效果的本质是这张图片中颜色的变化，接下来我们将讨论如何量化这些效果，也就是量化这些颜色的变化。为了这些效果，我们也要对object的性质量化：表面颜色，反射性质，透射性质，然后利用公式计算得到每一个像素点的颜色。
首先确定目标，我们生成的是一张图片，图片上的每一个点就是一个pixel，我们要计算的是每一个pixel的RGB值。
Basic Knowledge
主要用到的数学工具：线性代数、几何知识，其中向量和向量的运算在整个光线追踪中非常重要。
线性代数函数：实现加减、数乘、向量长度、向量正则化、点乘（用于判断两个向量的方向，计算投影长度）、叉乘（计算和两个向量构成的平面垂直的向量）。
涉及到的三维物体：球体、圆柱体、圆环、立方体和任意形状的模型。这里要求我们对这些三维物体有一个参数化的定义，这样我们才能通过数学工具精确地定义这些物体的位置，任意形状的模型可以使用三角面片来表示。。
我们这里将光源设计成球体，也可以设计成perfect points。
物理学知识：在正常场景上，光线是沿直线传播的，当然假如你觉得特别无聊的话，可以为物体赋予质量，利用质能方程可以让光线进行偏移，实现引力红移效果:)。
Fresnel定律：
光线照射到透明物体上时，一部分发生反射，一部分会在介质交界处发生折射，被反射和被折射的光通量存在一定的比率关系，这个关系可以根据Fresnel定律量化。根据Fresnel定律计算得出的数据是Fresnel系数，分为反射系数和折射系数，它们的和为1。一个完整的Fresnel公式依赖于折射率，消光率和入射角度等因素。
F=f_0+(1-f_0)(1-V*H)^5
f_0是入射角度接近0时的反射系数，V是指向视点的观察方向，H是半角向量。随着入射角趋近于直角，反射系数趋近1，则所有入射光都会被反射。F是反射系数，将用于折射光线和反射光线产生的效果的混合。
程序结构
我们给出了一个demo程序，其中只定义了球体，这意味着你可以参考这个代码，然后扩展:)
我们在mac OS和Windows都进行了编译：
使用clang编译：
cl -o a.exe RayTracer.cpp
a.exe
使用g++编译：
g++ -o a RayTracer.cpp
./a
首先我们定义了一个向量class，可以用来处理RGB信息、点的坐标和射线方向，随后重载了向量的运算，这里要注意的几个地方：要定义向量长度计算，定义向量正则化计算。
接下来我们定义了球体类，这个类有两个用处，第一个用处帮助我们创建场景内的球体，第二可以帮我们定义一个光源。
定义一个球体需要使用球心，半径，但是我们还需要定义这个球体的材质，能否折射（透明），能否发光（光源），能否反射（反射系数），表面颜色。这里要注意漫反射和反射材质有所区别。
这个类内最重要的函数intersect：
输入射线的起点和方向，计算射线和该球体是否有交点，如果有，交点位置。
我们可以通过几何的方法简单得到这些点，计算起点和球心的向量，计算这个向量向射线的投影，这样我们就可以构建一个直角三角形，然后根据勾股定理，观察是否相交，如果相交，再了解球心和交点，构造第二个直角三角形，然后计算得到起点沿射线方向到交点的距离，这样知道起点、方向和距离，我们就可以知道交点在哪里了。
定义一个球体有明显的缺陷，这意味着我们的场景里只会存在球类物体，我们需要添加更多种类的物体，但是万变不离其宗，除了定义物体的性质，这个类最重要的工作是计算任意一条射线和这个物体的相交情况，我们可以设计一个父类，可以处理表面材质这些性质，然后其他物体形状类继承这个类，分别实现定义和相交函数。
接下来就是具体的trace函数实现，最重要的观点是：这是一个递归函数，我们必须让它能停下来，这就是设置最大递归深度的原因。另一个是从人眼出来的反向光线，并不能在交点留下任何信息，那么交点的颜色信息来自哪里呢？是来自光源的照射（包括直接照射，经过反射、折射的间接照射），没有光源的照射，就没有颜色信息，因此最关键的是判断这个交点和光源的关系，以及反射光线和折射光线从光源带来的颜色信息。
因此trace函数的输入是：从观察点到图像中的像素点的射线的起点和方向，整个场景中所有的物体，还有终止条件：递归深度。
这条射线可能和场景中很多物体相交，但是起作用的只有第一个交点，因此整个函数的第一步是遍历场景中所有的物体，找到第一个交点。这里就用到了我们定义的球体类中的intersect函数。这时会有不同的情况：
如果不存在交点，返回背景颜色，这个颜色是自己定义的。
如果存在交点，定义一个surface color，用于记录这个像素的颜色。我们可以通过向量加减乘除，简单计算出和球体的交点，交点的法线方向（当然不同的几何体法线计算不同，你可以考虑自己在物体类内实现一个法线计算函数）。记得判断一下射线起点是否在第一个交点所在物体内，进行一下符号变换，这个在折射光线计算的时候很重要。
要定义一小段距离bias，可以在场景中给定点上投影阴影的很小的一段距离，这样可以避免对象对自己投射阴影。
接下来，我们就要量化反射、折射和阴影的效果了！反射，折射和阴影效果和物体材质有关，和物体位置与光源的关系有关，和入射夹角和法线方向有关。
首先要查看物体能否折射或反射，这取决于transparency和reflection两个参数，这两个参数代表了物体的反射和折射性质，如果存在反射折射性质，我们就要递归跟踪反射折射光线了，当然记住要求递归深度没有越界。
计算Fresnel Effect，它和入射光线方向和法线方向有关，这里我们用三次方来算，为了更快。f0=0.1和物体的材质有关系。
float facingratio = - raydirection.dot(nhit);
float fresneleffect = mix (pow(1 - facingratio, 3), 1, 0.1);
float mix(const float &a, const float &b, const float &mix)
{
return b * mix + a * (1 - mix);
}
计算反射射线：构造一个等腰三角形，底是两倍的入射光线向法线方向的投影。随后追踪这条射线，计算这条射线对这个点产生的颜色，记得depth加一。
多层递归反射在镜面反射中很重要，构建更真实的镜面效果。
计算折射射线：
$$
n_isin\theta_i=n_Tsin\theta_T
cos\theta_T=\sqrt{1-\frac{n_i2(1-cos2\theta_i)}{n^2_T}}IN=cos\theta_i|I||N|
(-N)T=cos\theta_T|N||T|
T=-\frac{n_i}{n_T}I+(\frac{n_i}{n_T}(IN)-cos\theta_T)*N
$$
（这部分是Latex公式，在简书中支持不是很好，可以直接查看夏老师的PPT）
如果球体是透明的，就可以产生折射了，在这里我们假设跟物体材质有关的参数ior，也就是n_i/n_T，当然要时刻记住我们可能在物体外，也可能在物体内。随后跟踪这条折射光线。
这样我们计算得到了反射光线和折射光线在这一点产生的颜色，如何让这里的颜色具备真实感呢，我们要引入物理学知识，最后的颜色还和物体表面颜色有关，因此利用公式计算得到surface color。
surfaceColor = (next_reflection * fresneleffect + next_refraction * (1 - fresneleffect) * sphere->transparency) * sphere->surfaceColor;
还有阴影效果没有考虑，考虑阴影效果时，我们必须考虑光线和光源的相交情况，根据前面的公式，我们在else条件里处理的是反射折射参数为0的交点，当然要定义光源表面的反射折射参数都为0。
还有一种情况是物体表面非常粗糙，这样就会产生一种漫反射的现象，此时这一点的颜色直接根据表面颜色、交点与光源的位置关系决定，如果交点和光源之间有阻挡，我们就要产生阴影效果了。
遍历整个场景，找到光源，计算光源到交点的射线会不会在中途和其他物体相交，如果相交，令此点rgb=0。当然为了让此点的颜色具备真实感，我们仍然要借用物理学的知识，这点的颜色和物体表面颜色，光源颜色，光线方向有关。为什么用+=，因为可能不止一个光源。
surfaceColor += sphere->surfaceColor * transmission * std::max(float(0), nhit.dot(lightDirection)) * spheres[i].emissionColor;
最后，除了光源，其他球体的emissionColor均为0，因此返回surfaceColor，假如此时相交的球体是光源，那么我们直接赋值为光源的光，因为光源的surfaceColor我们定义为0。
Render程序：此时就要生成一个图像了，因此我们的输入是整个场景，我们从起点到每一个像素坐标做一条射线，然后执行trace程序，返回这个像素的颜色信息，最后保存成图片格式就可以了。
你可以设置自己的图像大小，然后创建一个数组来保存每个图像的像素值，我们设置起点坐标为（0，0，0），你要计算每一个像素值的空间位置，这样你才能定义这条射线。因此我们需要定义我们的视角范围，这样我们才能用几何方法确认坐标，进而计算surfaceColor。
在这里为了方便输出文件格式的方便，选用了ppm文件格式，这样按照要求输入文件头、宽度、高度以及各点的颜色信息，就可以保存为ppm文件，这是在Linux下的图片文件格式，可以使用https://www.coolutils.com/online/PPM-to-PNG 来转成合适的图片格式，但是我们欢迎直接生成png、jpeg这种文件格式，你可以直接使用openCV库来输出合适的文件格式。
最后在main函数中我们添加场景中的物体，添加光源，渲染整个场景，就能得到渲染后的结果了。
关键结构
定义Vector：主要实现：加减，点乘，长度以及正则化，这个类主要用于定义光线的起点和方向；
定义球体：主要包含球心、半径、表面颜色，能否发光，还需要定义球体的材质，包含反射系数，折射系数。 最重要的是计算从某点发出的射线与这个球体相交时的情况，并返回距离，这样就可以根据起点，方向和大小得到交点的位置。
定义光源：光源可以看作是特殊的球体，发光系数为正值，表面颜色为0。从物体到光源仍然计算的是能否相交。
定义多面体：与定义球体类似，要包含它的性质，同时也要包含它与某点发来的射线的交点的计算。
定义光线追踪函数：
首先这是一个递归函数，因此这个函数最关键的是对边界条件的处理，以及不能让这个函数无限度地继续下去，因此需要定义最大递归深度。输入是：射线起点和方向，以及场景内的所有物体和递归深度。首先计算此射线是否与场景内的物体相交，如果不相交，返回背景颜色。如果相交，返回所有相交点中与之距离最近的那个点和这个点所在的物体。接着计算点的坐标和法线方向，法线方向将用于计算反射和折射效果。还要判断此时射线的起点是在物体的内部还是外部。如果物体表面是diffuse，也就是漫反射的，我们就不关心反射和折射的问题了，此时我们关心这个点处在阴影内还是阴影外，也就是添加阴影效果，是否处在阴影内取决于这个点和光源之间是否有障碍物，此时判断点和光源是否有其他交点，如果没有的话，就意味着没有阴影效果，添加物体表面颜色，如果有的话就要添加黑色阴影效果了。
如果这个物体表面是可以反射或者折射的，那么就要根据法线方向和入射方向计算反射方向和折射方向。而反射或者折射后的方向需要继续进行追踪，得到正向过程中的光线到达这里时产生的颜色。在这里，需要使用fresnel function来计算反射和折射的综合效果。
定义Render函数：
上述的trace对应某一条射线，但是从camera到成像平面有很多条射线，因此需要对每一条射线都计算一次光线追踪，并得到此射线对应的像素值，最后把结果写入图像文件中。
定义main函数：
对于一个场景中，需要包含一个光源，通常光源是球体，还需要定义多个多面体，然后调用render函数进行渲染。
程序展示
x, y, z轴的右手守则
Vector Class：定义各种向量运算。
Light Source Class：可以看成球体的一个特例。但是如果不加光源，只会是黑色的。
Object Class：定义Object的运动，旋转、平移，和某一射线的交点（根据不同的形状）。如果想做视频的话，你可以考虑定义物体的运动或者视角的运动，这样你可以逐帧渲染，最后生成一个好的视频。
有的时候我们用Triangle Mesh表示三维模型，这个时候我们计算到的交点的性质需要包含这个点的三角面片的三个顶点决定。
Trace Ray Function：just return a color。
继续优化
1.构建一个Scene（包括背景颜色），向Scene中添加物体或者Set Union。
2.构建多样化的多面体，包括圆环、cube等。可以构造一个SetUnion作为允许装载多种object的容器。光线与平面、三角形、多边形、长方体等求交。
3.添加纹理效果（面片顶点指定纹理，对面片和光线的交点进行插值）。
4.加速算法：包围盒加速、层次结构加速。
程序
这个程序来自于refrence的网站中，我们做了一些修改，按照上述方法可直接编译运行。
//Compile using clang under Windows: cl -o RayTracer.exe RayTracer.cpp
#include <cstdlib>
#include <cstdio>
#include <cmath>
#include <fstream>
#include <vector>
#include <iostream>
#include <cassert>
#include <algorithm>
#define M_PI 3.141592653589
#define INFINITY 1e8
//define class Vec_3, used in ray direction
template<typename T>
class Vec_3
{
public:
T x, y, z;
Vec_3(): x(T(0)), y(T(0)), z(T(0)) {}
Vec_3(T xx): x(xx), y(xx), z(xx) {}
Vec_3(T xx, T yy, T zz): x(xx), y(yy), z(zz){}
Vec_3<T> operator * (const T &f) const
{ return Vec_3<T>(x * f, y * f, z * f);}
Vec_3<T> operator * (const Vec_3<T> &v) const
{ return Vec_3<T>(x * v.x, y * v.y, z * v.z);}
T dot(const Vec_3<T> &v) const
{ return x * v.x + y * v.y + z * v.z;}
Vec_3<T> operator - (const Vec_3<T> &v) const
{ return Vec_3<T>( x - v.x, y - v.y, z - v.z);}
Vec_3<T> operator + (const Vec_3<T> &v) const
{ return Vec_3<T>( x + v.x, y + v.y, z + v.z);}
Vec_3<T>& operator += (const Vec_3<T> &v)
{
x += v.x;
y += v.y;
z += v.z;
return *this;
}
Vec_3<T>& operator *= (const Vec_3<T> &v)
{
x *= v.x;
y *= v.y;
z *= v.z;
return *this;
}
Vec_3<T> operator - () const
{
return Vec_3<T>(-x, -y, -z);
}
T length2() const
{
return x * x + y * y + z * z;
}
T length() const
{
return sqrt(length2());
}
Vec_3& normal()
{
T nor2= length2();
if (nor2 > 0)
{
T nor2_inv= 1 / sqrt(nor2);
x *= nor2_inv;
y *= nor2_inv;
z *= nor2_inv;
}
return *this;
}
friend std::ostream & operator << (std::ostream &os, const Vec_3<T> &v)
{
os<< "[" << v.x << " " << v.y << " " << v.z << "]";
return os;
}
};
typedef Vec_3<float> Vec_3f;
//Define Sphere Class
class Sphere
{
public:
Vec_3f center;
float radius, radius2;
Vec_3f surfaceColor, emissionColor;
float transparency, reflection;
Sphere(
const Vec_3f &c,
const float &r,
const Vec_3f &sc,
const float &refl = 0,
const float &transp = 0,
const Vec_3f &ec = 0):
center(c), radius(r), radius2(r * r), surfaceColor(sc), emissionColor(ec),
transparency(transp), reflection(refl)
{}
//Use geometric solution to solve a ray-sphere intersection
bool intersect(const Vec_3f &rayorigin, const Vec_3f & raydirection, float &t0, float &t1) const
{
Vec_3f l = center - rayorigin;
//Determine whether reverse direction
float tca = l.dot(raydirection);
if (tca < 0) return false;
//a^2=b^2+c^2
float dist = l.dot(l) - tca * tca;
if (dist > radius2) return false;
float thc = sqrt(radius2 - dist);
//t0: first intersection distance, t1: second intersection distance
t0 = tca - thc;
t1 = tca + thc;
return true;
}
};
//Define the maximum recursion depth
#define MAX_DEPTH 5
//Calculate the mix value for reflection and refraction
float mix(const float &a, const float &b, const float &mix)
{
return b * mix + a * (1 - mix);
}
//Ray Tracing Function: takes a ray (defined by its origin and direction) as argument.
//Through the function, we can know if the ray intersects any of the geometry in the scene.
//If the ray intersects an object, calculate the intersection point and its normal, then shade the point.
//Shading depends on the surface (transparent, reflective, diffuse)
//If the ray intersects an object, then return the color of the object at the intersection point, otherwise return the backgroud color.
Vec_3f trace(
const Vec_3f &rayorigin,
const Vec_3f &raydirection,
const std::vector<Sphere> &spheres,
const int &depth
)
{
float tnear= INFINITY;
const Sphere* sphere=NULL;
//calculate intersection of this ray with the sphere in the scene
for(unsigned i=0; i < spheres.size(); i++)
{
float t0=INFINITY;
float t1=INFINITY;
if(spheres[i].intersect(rayorigin, raydirection, t0, t1))
{
//If the point in the sphere
if(t0 < 0) t0= t1;
if(t0 < tnear)
{
tnear = t0;
sphere = &spheres[i];
}
}
}
//If there is no intersection, then return backgroud color
if(!sphere) return Vec_3f(0);
//Color of ray
Vec_3f surfaceColor = 0;
//point of intersect
Vec_3f phit = rayorigin + raydirection * tnear;
//normal of the intersection point
Vec_3f nhit = phit - sphere->center;
//normalize the normal direction
nhit.normal();
//If the normal and the view direction's dot is positive, means the view point inside sphere
float bias = 1e-4;
bool inside = false;
if(raydirection.dot(nhit) > 0)
{
nhit = -nhit;
inside = true;
}
//Tackle with relection and refraction
if((sphere->transparency > 0 || sphere->reflection > 0) && depth < MAX_DEPTH)
{
//Compute fresnel effect
float facingratio = - raydirection.dot(nhit);
float fresneleffect = mix (pow(1 - facingratio, 3), 1, 0.1);
//Compute reflection direction
Vec_3f reflect_direction = raydirection - nhit * 2 * raydirection.dot(nhit);
reflect_direction.normal();
Vec_3f next_reflection = trace(phit + nhit * bias, reflect_direction, spheres, depth + 1);
//Vec_3f next_reflection = trace(phit, reflect_direction, spheres, depth + 1);
Vec_3f next_refraction = 0;
//Only if the sphere is transparent, then compute refraction ray
if(sphere->transparency)
{
//judge whether we are inside or outside? ior is the index of two materials
float ior = 1.1, eta = (inside) ? ior : 1 / ior;
float cosi = -nhit.dot(raydirection);
float k = 1 - eta * eta * (1 - cosi * cosi);
Vec_3f refraction_direction = raydirection * eta + nhit * (eta * cosi - sqrt(k));
refraction_direction.normal();
next_refraction = trace(phit - nhit * bias, refraction_direction, spheres, depth+1);
//next_refraction = trace(phit, refraction_direction, spheres, depth+1);
}
//The surface is a mix of reflection and refraction (if the sphere is transparent)
surfaceColor = (next_reflection * fresneleffect + next_refraction * (1 - fresneleffect) * sphere->transparency) * sphere->surfaceColor;
}
//If it is a diffuse object, no need to ray tracing.
else
{
for(unsigned i = 0; i < spheres.size(); i++)
{
//This is a light
if(spheres[i].emissionColor.x > 0)
{
Vec_3f transmission = 1;
Vec_3f lightDirection = spheres[i].center - phit;
lightDirection.normal();
//Check whether have an obstacle between light and object, add shadow
for(unsigned j = 0; j < spheres.size(); ++j)
{
if(i != j)
{
float t0, t1;
if(spheres[j].intersect(phit + nhit * bias, lightDirection, t0, t1))
//if(spheres[j].intersect(phit, lightDirection, t0, t1))
{
transmission = 0;
break;
}
}
}
//If nhit and lightDirection's dot is less than 0, then no light.
surfaceColor += sphere->surfaceColor * transmission * std::max(float(0), nhit.dot(lightDirection)) * spheres[i].emissionColor;
}
}
}
return surfaceColor + sphere->emissionColor;
}
//Render function, compute each pixel of the image.
void render(const std::vector<Sphere> &spheres)
{
unsigned width = 640, height = 480;
Vec_3f *img = new Vec_3f[width * height], *pixel = img;
float invWidth = 1 / float(width), invHeight = 1 / float(height);
float fov = 30;
float aspectratio = width / float(height);
float angle = tan(M_PI * 0.5 * fov / 180.);
//Trace all ray
for(unsigned y = 0; y < height; y++)
{
for(unsigned x = 0; x < width; x++, pixel++)
{
float xx = (2 * ((x + 0.5) * invWidth) - 1) * angle * aspectratio;
float yy = (1 - 2 * ((y + 0.5) * invHeight)) * angle;
Vec_3f raydir(xx, yy, -1);
raydir.normal();
*pixel = trace(Vec_3f(0), raydir, spheres, 0);
}
}
//Save the result
std::ofstream ofs("./1.ppm", std::ios::out | std::ios::binary);
ofs << "P6\n" << width << " " << height << "\n255\n";
for(unsigned i = 0; i < width * height; i++)
{
//0,255
ofs << (unsigned char)(std::min(float(1), img[i].x) * 255) <<
(unsigned char)(std::min(float(1), img[i].y) * 255) <<
(unsigned char)(std::min(float(1), img[i].z) * 255);
}
ofs.close();
delete [] img;
}
//Create a sign including 5 spheres and 1 light (which is also a sphere), then render it.
int main()
{
std::vector<Sphere> spheres;
//argument: position, radius, surfaceColor, reflectivity, transparency, emissionColor
spheres.push_back(Sphere(Vec_3f( 0.0, 0, -20), 4, Vec_3f(1.00, 0.00, 0.00), 1, 0.5));
spheres.push_back(Sphere(Vec_3f( 5.0, -1, -15), 2, Vec_3f(0.00, 1.00, 0.00), 1, 0.0));
spheres.push_back(Sphere(Vec_3f( 5.0, 0, -25), 3, Vec_3f(0.00, 0.00, 1.00), 1, 0.0));
spheres.push_back(Sphere(Vec_3f(-5.5, 0, -15), 3, Vec_3f(0.00, 1.00, 0.00), 1, 0.0));
//Light
spheres.push_back(Sphere(Vec_3f(0.0, 20, -30), 3, Vec_3f(0.0, 0.0, 0.0), 0, 0.0, Vec_3f(3)));
render(spheres);
return 0;
}
结果
正常结果：
1 (2).png
不包含bias结果：
1 (1).png
包含背景光，不包含光源结果：
1 (3).png
包含光源，不包含背景光结果：
1 (4).png
不包含光源、背景光结果：
1 (5).png
最短的光线追踪程序
#include <stdlib.h> // card > aek.ppm
#include <stdio.h>
#include <math.h>
typedef int i;typedef float f;struct v{f x,y,z;v operator+(v r){return v(x+r.x,y+r.y,z+r.z);}v operator*(f r){return v(x*r,y*r,z*r);}f operator%(v r){return x*r.x+y*r.y+z*r.z;}v(){}v operator^(v r){return v(y*r.z-z*r.y,z*r.x-x*r.z,x*r.y-y*r.x);}v(f a,f b,f c){x=a;y=b;z=c;}v operator!(){return*this*(1/sqrt(*this%*this));}};i G[]={247570,280596,280600,249748,18578,18577,231184,16,16};f R(){return(f)rand()/RAND_MAX;}i T(v o,v d,f&t,v&n){t=1e9;i m=0;f p=-o.z/d.z;if(.01<p)t=p,n=v(0,0,1),m=1;for(i k=19;k--;)for(i j=9;j--;)if(G[j]&1<<k){v p=o+v(-k,0,-j-4);f b=p%d,c=p%p-1,q=b*b-c;if(q>0){f s=-b-sqrt(q);if(s<t&&s>.01)t=s,n=!(p+d*t),m=2;}}return m;}v S(v o,v d){f t;v n;i m=T(o,d,t,n);if(!m)return v(.7,.6,1)*pow(1-d.z,4);v h=o+d*t,l=!(v(9+R(),9+R(),16)+h*-1),r=d+n*(n%d*-2);f b=l%n;if(b<0||T(h,l,t,n))b=0;f p=pow(l%r*(b>0),99);if(m&1){h=h*.2;return((i)(ceil(h.x)+ceil(h.y))&1?v(3,1,1):v(3,3,3))*(b*.2+.1);}return v(p,p,p)+S(h,r)*.5;}i main(){printf("P6 512 512 255 ");v g=!v(-6,-16,0),a=!(v(0,0,1)^g)*.002,b=!(g^a)*.002,c=(a+b)*-256+g;for(i y=512;y--;)for(i x=512;x--;){v p(13,13,13);for(i r=64;r--;){v t=a*(R()-.5)*99+b*(R()-.5)*99;p=S(v(17,16,8)+t,!(t*-1+(a*(R()+x)+b*(y+R())+c)*16))*3.5+p;}printf("%c%c%c",(i)p.x,(i)p.y,(i)p.z);}}
