在“Q79”和“Q80”中用三角形网格细分曲面时，都是将每一个三角形的三个顶点的坐标都保存在内存中。这句话有两个重点：
其一，每个三角形的三个顶点的坐标都计算了一次。但是，每个顶点都是被好几个三角形公用的，所以每个顶点的坐标被重复计算了好几次。
其二，所有顶点的坐标都是临时计算的，然后保存到内存中的。如果顶点数量特别大（百万级、千万级），这样临时计算则会导致程序运行缓慢。是否可先将顶点的坐标都算好保存在文件中，程序运行时直接读这个文件就可以了呢？
这种保存“将图形用三角形网格细分后得到的三角形信息”的文件被称为“PLY文件”。
PLY文件由两部分组成：文件头、文件体。
文件头是由以“回车（carriage-return）”为结束符的ASCII字符串组成的。
文件体是由“元素”组成。一般情况，会包含两类“元素”：第一类，顶点信息（含顶点的坐标值、法向量、纹理映射值等等）；第二类，多边形各个顶点在第一类元素中的索引值。
如下方示意：
PLY文件的内容是非常灵活的。第一类元素“顶点vertex”，除了可以包含顶点的坐标值，还可以包含顶点的“法向量”、“纹理坐标”、“颜色值”等等；第二类元素“face”，其实是可以表示任意多边形的，而且不同顶点数的各种多边形可以同时被包含在同一个PLY文件中。
参考：
http://blog.csdn.net/libing_zeng/article/details/60878097
转换后的PLY文件截图如下：
这些文件的路径：http://download.csdn.net/detail/libing_zeng/9775786
Greg Turk发明PLY文件格式之后，就有提供通用的可拓展的读取PLY文件的程序。相关链接如下：http://www.cc.gatech.edu/projects/large_models/ply.html
但是，本人是直接移植了《Ray Tracing from the Ground Up》的相关代码。在移植过程中遇到两个问题：
《Ray Tracing from the Ground Up》中提供的调用读取PLY文件接口的函数如下：
// ----------------------------------------------------------------------------- read_ply_file
// Most of this function was written by Greg Turk and is released under the licence agreement
// at the start of the PLY.h and PLY.c files
// The PLY.h file is #included at the start of this file
// It still has some of his printf statements for debugging
// I've made changes to construct mesh triangles and store them in the grid
// mesh_ptr is a data member of Grid
// objects is a data member of Compound
// triangle_type is either flat or smooth
// Using the one function construct to flat and smooth triangles saves a lot of repeated code
// The ply file is the same for flat and smooth triangles
void
Grid::read_ply_file(char* file_name, const int triangle_type) {
// Vertex definition
typedef struct Vertex {
float x,y,z; // space coordinates
} Vertex;
// Face definition. This is the same for all files but is placed here to keep all the definitions together
typedef struct Face {
unsigned char nverts; // number of vertex indices in list
int* verts; // vertex index list
} Face;
// list of property information for a vertex
// this varies depending on what you are reading from the file
PlyProperty vert_props[] = {
{"x", PLY_FLOAT, PLY_FLOAT, offsetof(Vertex,x), 0, 0, 0, 0},
{"y", PLY_FLOAT, PLY_FLOAT, offsetof(Vertex,y), 0, 0, 0, 0},
{"z", PLY_FLOAT, PLY_FLOAT, offsetof(Vertex,z), 0, 0, 0, 0}
};
// list of property information for a face.
// there is a single property, which is a list
// this is the same for all files
PlyProperty face_props[] = {
{"vertex_indices", PLY_INT, PLY_INT, offsetof(Face,verts),
1, PLY_UCHAR, PLY_UCHAR, offsetof(Face,nverts)}
};
// local variables
int i,j;
PlyFile* ply;
int nelems; // number of element types: 2 in our case - vertices and faces
char** elist;
int file_type;
float version;
int nprops; // number of properties each element has
int num_elems; // number of each type of element: number of vertices or number of faces
PlyProperty** plist;
Vertex** vlist;
Face** flist;
char* elem_name;
int num_comments;
char** comments;
int num_obj_info;
char** obj_info;
// open a ply file for reading
ply = ply_open_for_reading(file_name, &nelems, &elist, &file_type, &version);
// print what we found out about the file
printf ("version %f\n", version);
printf ("type %d\n", file_type);
// go through each kind of element that we learned is in the file and read them
for (i = 0; i < nelems; i++) { // there are only two elements in our files: vertices and faces
// get the description of the first element
elem_name = elist[i];
plist = ply_get_element_description (ply, elem_name, &num_elems, &nprops);
// print the name of the element, for debugging
cout << "element name " << elem_name << " num elements = " << num_elems << " num properties = " << nprops << endl;
// if we're on vertex elements, read in the properties
if (equal_strings ("vertex", elem_name)) {
// set up for getting vertex elements
// the three properties are the vertex coordinates
ply_get_property (ply, elem_name, &vert_props[0]);
ply_get_property (ply, elem_name, &vert_props[1]);
ply_get_property (ply, elem_name, &vert_props[2]);
// reserve mesh elements
mesh_ptr->num_vertices = num_elems;
mesh_ptr->vertices.reserve(num_elems);
// grab all the vertex elements
for (j = 0; j < num_elems; j++) {
Vertex* vertex_ptr = new Vertex;
// grab an element from the file
ply_get_element (ply, (void *) vertex_ptr);
mesh_ptr->vertices.push_back(Point3D(vertex_ptr->x, vertex_ptr->y, vertex_ptr->z));
delete vertex_ptr;
}
}
// if we're on face elements, read them in
if (equal_strings ("face", elem_name)) {
// set up for getting face elements
ply_get_property (ply, elem_name, &face_props[0]); // only one property - a list
mesh_ptr->num_triangles = num_elems;
objects.reserve(num_elems); // triangles will be stored in Compound::objects
// the following code stores the face numbers that are shared by each vertex
mesh_ptr->vertex_faces.reserve(mesh_ptr->num_vertices);
vector<int> faceList;
for (j = 0; j < mesh_ptr->num_vertices; j++)
mesh_ptr->vertex_faces.push_back(faceList); // store empty lists so that we can use the [] notation below
// grab all the face elements
int count = 0; // the number of faces read
for (j = 0; j < num_elems; j++) {
// grab an element from the file
Face* face_ptr = new Face;
ply_get_element (ply, (void *) face_ptr);
// construct a mesh triangle of the specified type
if (triangle_type == flat) {
FlatMeshTriangle* triangle_ptr = new FlatMeshTriangle(mesh_ptr, face_ptr->verts[0], face_ptr->verts[1], face_ptr->verts[2]);
triangle_ptr->compute_normal(reverse_normal);
objects.push_back(triangle_ptr);
}
if (triangle_type == smooth) {
SmoothMeshTriangle* triangle_ptr = new SmoothMeshTriangle(mesh_ptr, face_ptr->verts[0], face_ptr->verts[1], face_ptr->verts[2]);
triangle_ptr->compute_normal(reverse_normal); // the "flat triangle" normal is used to compute the average normal at each mesh vertex
objects.push_back(triangle_ptr); // it's quicker to do it once here, than have to do it on average 6 times in compute_mesh_normals
// the following code stores a list of all faces that share a vertex
// it's used for computing the average normal at each vertex in order(num_vertices) time
mesh_ptr->vertex_faces[face_ptr->verts[0]].push_back(count);
mesh_ptr->vertex_faces[face_ptr->verts[1]].push_back(count);
mesh_ptr->vertex_faces[face_ptr->verts[2]].push_back(count);
count++;
}
}
if (triangle_type == flat)
mesh_ptr->vertex_faces.erase(mesh_ptr->vertex_faces.begin(), mesh_ptr->vertex_faces.end());
}
// print out the properties we got, for debugging
for (j = 0; j < nprops; j++)
printf ("property %s\n", plist[j]->name);
} // end of for (i = 0; i < nelems; i++)
// grab and print out the comments in the file
comments = ply_get_comments (ply, &num_comments);
for (i = 0; i < num_comments; i++)
printf ("comment = '%s'\n", comments[i]);
// grab and print out the object information
obj_info = ply_get_obj_info (ply, &num_obj_info);
for (i = 0; i < num_obj_info; i++)
printf ("obj_info = '%s'\n", obj_info[i]);
// close the ply file
ply_close (ply);
}
简单来看，分如下几个步骤：
1，打开文件；
2，读取头文件；
3，读取顶点元素；
4，读取三角形元素，依次：根据索引在顶点元素中找到对应的顶点值，创建一个新的FlatMeshTriangle/SmoothMeshTriangle对象。
5，关闭文件；
其中第4步的截图如下：
判断某个撞击点p是否在阴影之中，只需要判断从撞击点p到光源位置的阴影光线是否撞击到其他物体。这里我们考虑的是Grid图形的阴影，所以只需要判断阴影光线是否撞击到Grid图形。所以，我们需要实现Grid类的shadow_hit()函数：
bool
Grid::shadow_hit(const Ray& ray, double& t) const {
double ox = ray.o.x;
double oy = ray.o.y;
double oz = ray.o.z;
double dx = ray.d.x;
double dy = ray.d.y;
double dz = ray.d.z;
double x0 = bbox.x0;
double y0 = bbox.y0;
double z0 = bbox.z0;
double x1 = bbox.x1;
double y1 = bbox.y1;
double z1 = bbox.z1;
double tx_min, ty_min, tz_min;
double tx_max, ty_max, tz_max;
// the following code includes modifications from Shirley and Morley (2003)
double a = 1.0 / dx;
if (a >= 0) {
tx_min = (x0 - ox) * a;
tx_max = (x1 - ox) * a;
}
else {
tx_min = (x1 - ox) * a;
tx_max = (x0 - ox) * a;
}
double b = 1.0 / dy;
if (b >= 0) {
ty_min = (y0 - oy) * b;
ty_max = (y1 - oy) * b;
}
else {
ty_min = (y1 - oy) * b;
ty_max = (y0 - oy) * b;
}
double c = 1.0 / dz;
if (c >= 0) {
tz_min = (z0 - oz) * c;
tz_max = (z1 - oz) * c;
}
else {
tz_min = (z1 - oz) * c;
tz_max = (z0 - oz) * c;
}
double t0, t1;
if (tx_min > ty_min)
t0 = tx_min;
else
t0 = ty_min;
if (tz_min > t0)
t0 = tz_min;
if (tx_max < ty_max)
t1 = tx_max;
else
t1 = ty_max;
if (tz_max < t1)
t1 = tz_max;
if (t0 > t1)
return(false);
// initial cell coordinates
int ix, iy, iz;
if (bbox.inside(ray.o)) { // does the ray start inside the grid?
ix = clamp((ox - x0) * nx / (x1 - x0), 0, nx - 1);
iy = clamp((oy - y0) * ny / (y1 - y0), 0, ny - 1);
iz = clamp((oz - z0) * nz / (z1 - z0), 0, nz - 1);
}
else {
Point3D p = ray.o + t0 * ray.d; // initial hit point with grid's bounding box
ix = clamp((p.x - x0) * nx / (x1 - x0), 0, nx - 1);
iy = clamp((p.y - y0) * ny / (y1 - y0), 0, ny - 1);
iz = clamp((p.z - z0) * nz / (z1 - z0), 0, nz - 1);
}
// ray parameter increments per cell in the x, y, and z directions
double dtx = (tx_max - tx_min) / nx;
double dty = (ty_max - ty_min) / ny;
double dtz = (tz_max - tz_min) / nz;
double tx_next, ty_next, tz_next;
int ix_step, iy_step, iz_step;
int ix_stop, iy_stop, iz_stop;
if (dx > 0) {
tx_next = tx_min + (ix + 1) * dtx;
ix_step = +1;
ix_stop = nx;
}
else {
tx_next = tx_min + (nx - ix) * dtx;
ix_step = -1;
ix_stop = -1;
}
if (dx == 0.0) {
tx_next = kHugeValue;
ix_step = -1;
ix_stop = -1;
}
if (dy > 0) {
ty_next = ty_min + (iy + 1) * dty;
iy_step = +1;
iy_stop = ny;
}
else {
ty_next = ty_min + (ny - iy) * dty;
iy_step = -1;
iy_stop = -1;
}
if (dy == 0.0) {
ty_next = kHugeValue;
iy_step = -1;
iy_stop = -1;
}
if (dz > 0) {
tz_next = tz_min + (iz + 1) * dtz;
iz_step = +1;
iz_stop = nz;
}
else {
tz_next = tz_min + (nz - iz) * dtz;
iz_step = -1;
iz_stop = -1;
}
if (dz == 0.0) {
tz_next = kHugeValue;
iz_step = -1;
iz_stop = -1;
}
// traverse the grid
while (true) {
GeometricObject* object_ptr = cells[ix + nx * iy + nx * ny * iz];
if (tx_next < ty_next && tx_next < tz_next) {
if (object_ptr && object_ptr->shadow_hit(ray, t) && t < tx_next) {
return (true);
}
tx_next += dtx;
ix += ix_step;
if (ix == ix_stop)
return (false);
}
else {
if (ty_next < tz_next) {
if (object_ptr && object_ptr->shadow_hit(ray, t) && t < ty_next) {
return (true);
}
ty_next += dty;
iy += iy_step;
if (iy == iy_stop)
return (false);
}
else {
if (object_ptr && object_ptr->shadow_hit(ray, t) && t < tz_next) {
return (true);
}
tz_next += dtz;
iz += iz_step;
if (iz == iz_stop)
return (false);
}
}
}
} // end of hit
由于Grid的每个cell对应的物体的指针可能有两种情况：
1，单个物体的指针（这里即是单个FlatMeshTriangle对象的指针）。
2，多个物体的指针，即Compound对象的指针。
所以，我们需要实现FlatMeshTriangle类和Compound类的shadow_hit()。当然，由于我们这里的图形是由FlatMeshTriangle组成，Compound类的shadow_hit()最终也会调用到FlatMeshTriangle类的shadow_hit()。
这就有个问题了：是不是可以只实现“FlatMeshTriangle类的shadow_hit()”而不用实现“Compound类的shadow_hit()”呢？
当然不行。如果这样，也就是只考虑了之前两种情况的“单个物体的指针”的情况，而忽略了“多个物体的指针”的情况。如果这样，最后看到的图形基本没有阴影，因为cell中对应的一般是多个物体的指针。
bool
Compound::shadow_hit(const Ray& ray, double& tmin) const {
double t;
Normal normal;
Point3D local_hit_point;
bool hit = false;
tmin = kHugeValue;
int num_objects = objects.size();
for (int j = 0; j < num_objects; j++)
if (objects[j]->shadow_hit(ray, t) && (t < tmin)) {
hit = true;
tmin = t;
}
return (hit);
}
bool
FlatMeshTriangle::shadow_hit(const Ray& ray, double& tmin) const {
Point3D v0(mesh_ptr->vertices[index0]);
Point3D v1(mesh_ptr->vertices[index1]);
Point3D v2(mesh_ptr->vertices[index2]);
double a = v0.x - v1.x, b = v0.x - v2.x, c = ray.d.x, d = v0.x - ray.o.x;
double e = v0.y - v1.y, f = v0.y - v2.y, g = ray.d.y, h = v0.y - ray.o.y;
double i = v0.z - v1.z, j = v0.z - v2.z, k = ray.d.z, l = v0.z - ray.o.z;
double m = f * k - g * j, n = h * k - g * l, p = f * l - h * j;
double q = g * i - e * k, s = e * j - f * i;
double inv_denom = 1.0 / (a * m + b * q + c * s);
double e1 = d * m - b * n - c * p;
double beta = e1 * inv_denom;
if (beta < 0.0)
return (false);
double r = e * l - h * i;
double e2 = a * n + d * q + c * r;
double gamma = e2 * inv_denom;
if (gamma < 0.0)
return (false);
if (beta + gamma > 1.0)
return (false);
double e3 = a * p - b * r + d * s;
double t = e3 * inv_denom;
if (t < kEpsilon)
return (false);
tmin = t;
return (true);
}
World;;build()函数如下：
void
World::build(void) {
int num_samples = 100;
vp.set_hres(400);
vp.set_vres(400);
vp.set_samples(1);
tracer_ptr = new RayCast(this);
background_color = black;
Pinhole* pinhole_ptr = new Pinhole;
#if 0 // 300*300: happy, Bunny69K, dragon
pinhole_ptr->set_eye(0, 0.1, 2.0);
pinhole_ptr->set_lookat(0, 0.1, 0);
pinhole_ptr->set_view_distance(1600);
pinhole_ptr->compute_uvw();
set_camera(pinhole_ptr);
PointLight* light_ptr1 = new PointLight;
light_ptr1->set_location(13, 10, 10);
light_ptr1->scale_radiance(3.0);
light_ptr1->set_cast_shadow(true);
add_light(light_ptr1);
#endif // 0
#if 0 // 300*300: Horse97K
pinhole_ptr->set_eye(0, 0.1, 6.0);
pinhole_ptr->set_lookat(0, 0.1, 0);
pinhole_ptr->set_view_distance(1600);
pinhole_ptr->compute_uvw();
set_camera(pinhole_ptr);
PointLight* light_ptr1 = new PointLight;
light_ptr1->set_location(13, 10, 10);
light_ptr1->scale_radiance(3.0);
light_ptr1->set_cast_shadow(true);
add_light(light_ptr1);
#endif // 0
#if 1 // 400*400: hand
pinhole_ptr->set_eye(25, 10, 25);
pinhole_ptr->set_lookat(1, 1, 0);
pinhole_ptr->set_view_distance(1600);
pinhole_ptr->compute_uvw();
set_camera(pinhole_ptr);
PointLight* light_ptr1 = new PointLight;
light_ptr1->set_location(13, 10, 10);
light_ptr1->scale_radiance(3.0);
light_ptr1->set_cast_shadow(true);
add_light(light_ptr1);
#endif // 0
#if 0 // 400*400: goldfish_high_res,
pinhole_ptr->set_eye(75, 20, 80);
pinhole_ptr->set_lookat(-0.05, -0.5, 0);
pinhole_ptr->set_view_distance(1600);
pinhole_ptr->compute_uvw();
set_camera(pinhole_ptr);
PointLight* light_ptr1 = new PointLight;
light_ptr1->set_location(13, 10, 10);
light_ptr1->scale_radiance(3.0);
light_ptr1->set_cast_shadow(true);
add_light(light_ptr1);
#endif // 1
Phong* phong_ptr1 = new Phong;
phong_ptr1->set_ka(0.4);
phong_ptr1->set_kd(0.8);
phong_ptr1->set_cd(1.0, 0.2, 0.0);
phong_ptr1->set_ks(0.5);
// phong_ptr1->set_cs(1.0, 1.0, 0.0);
phong_ptr1->set_exp(50.0);
// char* file_name = ".\\PLYFiles\\goldfish_low_res.ply";
// char* file_name = ".\\PLYFiles\\goldfish_high_res.ply";
// char* file_name = ".\\PLYFiles\\Horse2K.ply";
// char* file_name = ".\\PLYFiles\\Horse97K.ply";
// char* file_name = ".\\PLYFiles\\Bunny16K.ply";
// char* file_name = ".\\PLYFiles\\Bunny69K.ply";
// char* file_name = ".\\PLYFiles\\dragon.ply";
// char* file_name = ".\\PLYFiles\\happy.ply";
char* file_name = ".\\PLYFiles\\hand.ply";
Grid* grid_ptr = new Grid(new Mesh);
grid_ptr->read_flat_triangles(file_name);
// grid_ptr->read_smooth_triangles(file_name);
grid_ptr->set_material(phong_ptr1);
grid_ptr->setup_cells();
add_object(grid_ptr);
Phong* phong_ptr2 = new Phong;
phong_ptr2->set_ka(0.4);
phong_ptr2->set_kd(0.8);
phong_ptr2->set_cd(0.5, 1.0, 0.5);
phong_ptr2->set_ks(0.5);
// phong_ptr2->set_cs(1.0, 1.0, 0.0);
phong_ptr2->set_exp(50.0);
Plane* plane_ptr1 = new Plane(Point3D(0, -2.0, 0), Normal(0, 1, 0));
plane_ptr1->set_material(phong_ptr2);
add_object(plane_ptr1);
}
该函数当前生成的图形是hand skeleton：
如果要生成其他图形：horse、bunny、goldfish、dragon、happy buddha。则需要调节各自对应的：下方绿色平面的位置、点光源的位置、lookfrom的位置、lookat的位置，以便从各自最佳的角度观测到图形及其阴影。
下方贴图是没有下方绿色平面，没有阴影时的图形。
顺序如下：
完整的代码，参考：http://download.csdn.net/detail/libing_zeng/9776662
各种经典图形的PLY文件（包含Unix/Mac、Windows系统下的两种格式）：
http://download.csdn.net/detail/libing_zeng/9775786
Referrance
[1]. Kevin Suffern, Ray Tracing from theGround Up, A K Peters Ltd, 2007.
[2]. http://www.cc.gatech.edu/projects/large_models/
