[图形学]smallpt代码详解（下）

一、简介

本文紧接在[图形学]smallpt代码详解（上）和[图形学]smallpt代码详解（中）之后，继续详细讲解smallpt中的main函数部分，包括相机、屏幕参数设置、遍历像素、Tent滤波、计算radiance和保存结果为ppm文件几部分。

二、smallpt代码详解

0.main函数（第75到99行）

int main(int argc, char *argv[]){ 
int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples 
Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir 
Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h]; 
#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP 
for (int y=0; y<h; y++){                       // Loop over image rows 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols 
    for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows 
        for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols 
        for (int s=0; s<samps; s++){ 
            double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
            double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
            Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + 
                    cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d; 
            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
        } // Camera rays are pushed ^^^^^ forward to start in interior 
        c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
        } 
} 
FILE *f = fopen("image.ppm", "w");         // Write image to PPM file. 
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
for (int i=0; i<w*h; i++) 
    fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
} 

smallpt源代码的main()函数在第75到99行，下面我们一行一行解释各部分代码的功能。

1.相机、屏幕参数设置

int w=1024, h=768, samps = argc==2 ? atoi(argv[1])/4 : 1; // # samples 
Ray cam(Vec(50,52,295.6), Vec(0,-0.042612,-1).norm()); // cam pos, dir 
Vec cx=Vec(w*.5135/h), cy=(cx%cam.d).norm()*.5135, r, *c=new Vec[w*h]; 

代码第一行设置渲染结果的像素尺寸，宽高分别为1024*768。同时设置子像素内的采样光线数为sample。因为smallpt将一个像素设为2*2的子像素，因此每个子像素内的采样光线数为像素内的光线数除以4。
第二行设置相机的位置，位置为(50,52,295.6)，相机的朝向向量设为(0,-0.042612,-1)。
第三行cx向量为屏幕水平方向在世界坐标系下的宽度向量，设置为(w*.5135/h,0,0)。cy向量为屏幕垂直方向的宽度向量，cy向量的单位向量使用cx和cam.d叉乘得到，然后再乘以屏幕在世界坐标系下的高度，得到最终的cy。并且申请一个大小为w*h的Vec数组，用于存储渲染的结果。
下图展示了相机、屏幕和渲染场景的示意图。

2.光线跟踪

1).遍历像素

#pragma omp parallel for schedule(dynamic, 1) private(r)       // OpenMP 
for (int y=0; y<h; y++){                       // Loop over image rows 
    fprintf(stderr,"\rRendering (%d spp) %5.2f%%",samps*4,100.*y/(h-1)); 
    for (unsigned short x=0, Xi[3]={0,0,y*y*y}; x<w; x++)   // Loop cols 

该部分代码第一行启用OpenMP，用于代码并行加速。
第二行遍历屏幕空间Y方向的各像素。
第三行用于在程序运行时动态显示程序的渲染结果的进度。
第四行遍历屏幕空间X方向的各像素，同时将y*y*y作赋值给Xi，用于在光线跟踪中作为erand()函数的随机数种子值。

2).遍历子像素

    for (int sy=0, i=(h-y-1)*w+x; sy<2; sy++)     // 2x2 subpixel rows 
        for (int sx=0; sx<2; sx++, r=Vec()){        // 2x2 subpixel cols 
        for (int s=0; s<samps; s++){ 

因为smallpt将屏幕各个像素又细分为2*2的子像素，该部分代码即为在像素[x,y]内，遍历子像素。像素[x,y]在数组c中的下标为i，最后渲染得到的像素[x,y]的颜色写入到数组c[i]的内。每个子像素内部使用sample根光线。

3).Tent滤波

            double r1=2*erand48(Xi), dx=r1<1 ? sqrt(r1)-1: 1-sqrt(2-r1); 
            double r2=2*erand48(Xi), dy=r2<1 ? sqrt(r2)-1: 1-sqrt(2-r2); 
            Vec d = cx*( ( (sx+.5 + dx)/2 + x)/w - .5) + 
                    cy*( ( (sy+.5 + dy)/2 + y)/h - .5) + cam.d; 

该部分代码通过在像素[x,y]的四个子像素内部使用Tent filter 采样smaple个点，确定采样光线的方向d。
Tent filter即在一个矩形区域内，在X和Y方向都使用下图的PDF概率密度函数采样点。使用Tent filter采样是为了减少渲染结果的噪音。

一个使用Tent filter方法得到的采样点分布如下图所示：

4).计算radiance

            r = r + radiance(Ray(cam.o+d*140,d.norm()),0,Xi)*(1./samps); 
        } // Camera rays are pushed ^^^^^ forward to start in interior 
        c[i] = c[i] + Vec(clamp(r.x),clamp(r.y),clamp(r.z))*.25; 
        } 

这部分代码通过调用radiance()函数，计算该子像素上的颜色值。
在调用radiance()函数时，使用前面计算得到的光线方向d，计算从相机发出的光线在屏幕上的交点，cam.o+d*140，作为光线的出发点。光线的方向依旧为d。并将计算得到的radiance除以采样光线数smaple，累加到变量r上。
最后一行为将各个子像素上颜色值r累加到数组c[i]上，由于每个像素被划分为4个子像素，因此累加r时需要将r除以4（即乘以0.25）。另外，在累加r时需要保证r的各个分量在[0,1]范围内，因此调用clamp()函数约束r各个分量值的范围。

3.保存结果为ppm文件

FILE *f = fopen("image.ppm", "w");         // Write image to PPM file. 
fprintf(f, "P3\n%d %d\n%d\n", w, h, 255); 
for (int i=0; i<w*h; i++) 
    fprintf(f,"%d %d %d ", toInt(c[i].x), toInt(c[i].y), toInt(c[i].z)); 
} 

这部分代码将存在数组c中的渲染结果输出到一个image.ppm文件中，用于可视化渲染结果。

三、参考

[1].smallpt: Global Illumination in 99 lines of C++
[2].smallpt: Global Illumination in 99 lines of C+±Presentation slides
[3].光线跟踪smallpt详解 (二)