在“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/
