【光线追踪系列十八】混合密度函数采样_完结(附源码)

本文主要参照 Ray Tracing: The Rest of Your Life,其中只是主要精炼光追相关理论,具体实现可参照原文。
在这里插入图片描述

我们已经有2个pdf,一个和 cosine(θ) 相关,一个和跟光源采样相关。接下来,我们将他们结合起来生成一个混合密度。加权平均的pdf依然是一个pdf,例如:

在这里插入图片描述
伪代码示意:

if (random_double() < 0.5)
    Pick direction according to pdf_reflection
else
    Pick direction according to pdf_light

要估算pdf_mixture,就要估算 pdf_reflectionpdf_light,这需要知道两个值:

  • 当前方向的pdf值
  • 生成对应该pdf分布的随机数

接下来将pdf封装成类:

class pdf  {
    public:
        virtual float value(const vec3& direction) const = 0;	//当前方向的pdf值
        virtual vec3 generate() const = 0;	//生成对应该pdf分布的随机数
  	};

一、封装随机方向的pdf

先实现pdf_reflection,即之前和cosine(θ)相关的pdf。

class cosine_pdf : public pdf {
    public:
        cosine_pdf(const vec3& w) { uvw.build_from_w(w); }
        virtual float value(const vec3& direction) const {
            float cosine = dot(unit_vector(direction), uvw.w());
            if (cosine > 0)
                return cosine/M_PI;
            else
                return 0;
        }
        virtual vec3 generate() const  {
            return uvw.local(random_cosine_direction());
        }
        onb uvw;
};

修改color函数:

vec3 color(const ray& r, hittable *world, int depth) {
	hit_record rec;
	if (world->hit(r, 0.001, FLT_MAX, rec)) {
		ray scattered;
		vec3 attenuation;
		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
		float pdf;
		vec3 albedo;
		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {

			//new start
			float sizeX = lightEndX - lightStartX;
			float sizeZ = lightEndZ - lightStartZ;
			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
				lightY,
				lightStartZ + random_double()*sizeZ);
			vec3 to_light = on_light - rec.p;
			float distance_squared = to_light.squared_length();
			to_light.make_unit_vector();
			if (dot(to_light, rec.normal) < 0) return emitted;
			float light_cosine = fabs(to_light.y());
			if (light_cosine < 0.000001) return emitted;

			// origin
			//pdf = distance_squared / (light_cosine * light_area);
			//scattered = ray(rec.p, to_light, r.time());

			//cosine
			cosine_pdf p(rec.normal);
			scattered = ray(rec.p, p.generate(), r.time());
			pdf = p.value(scattered.direction());
			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
				*color(scattered, world, depth + 1) / pdf;
		}
		else
			return emitted;
	}
	else
		return vec3(0, 0, 0);
}

渲染效果和 【光线追踪系列十六】基于着色点的正向半球随机方向生成 中的一样:
在这里插入图片描述

二、封装对hittable采样的pdf

上一节我们实现了直接朝着光源方向采样的pdf,光源也是一种hittable类型。

接下来,封装hittable_pdf:

class hittable_pdf : public pdf {
    public:
        hittable_pdf(hittable *p, const vec3& origin) : ptr(p), o(origin) {}
        virtual float value(const vec3& direction) const {
            return ptr->pdf_value(o, direction);
        }
        virtual vec3 generate() const {
            return ptr->random(o);
        }
        vec3 o;
        hittable *ptr;
};

hittable_pdf类中声明了一个指向某个hittable的指针。
hittable基类中新增代码如下:

//物体的虚基类
class hittable
{
public:
	virtual bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const = 0;
	virtual bool bounding_box(float t0, float t1, aabb& box) const = 0;
	virtual float pdf_value(const vec3& o, const vec3& v) const { return 0.0; }
	virtual vec3 random(const vec3& o) const { return vec3(1, 0, 0); }
};

在灯光的hittable类xz_rect中实现pdf_value和random这两个方法:

class xz_rect : public hittable {
public:
	xz_rect() {}
	xz_rect(float _x0, float _x1, float _z0, float _z1, float _k, material *mat) :x0(_x0), x1(_x1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
	virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const override;

	virtual bool bounding_box(float time0, float time1, aabb& output_box) const override {
		// The bounding box must have non-zero width in each dimension, so pad the Y
		// dimension a small amount.
		output_box = aabb(vec3(x0, k - 0.0001, z0), vec3(x1, k + 0.0001, z1));
		return true;
	}
	virtual float  pdf_value(const vec3& o, const vec3& v) const {
		hit_record rec;
		if (this->hit(ray(o, v), 0.001, FLT_MAX, rec)) {
			float area = (x1 - x0)*(z1 - z0);
			float distance_squared = rec.t * rec.t * v.squared_length();
			float cosine = fabs(dot(v, rec.normal) / v.length());
			return  distance_squared / (cosine * area);
		}
		else
			return 0;
	}
	virtual vec3 random(const vec3& o) const {
		vec3 random_point = vec3(x0 + random_double()*(x1 - x0), k,
			z0 + random_double()*(z1 - z0));
		return random_point - o;
	}

	material  *mp;
	double x0, x1, z0, z1, k;
};

color函数修改如下:

vec3 color(const ray& r, hittable *world, int depth) {
	hit_record rec;
	if (world->hit(r, 0.001, FLT_MAX, rec)) {
		ray scattered;
		vec3 attenuation;
		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
		float pdf;
		vec3 albedo;
		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {

			//new start
			float sizeX = lightEndX - lightStartX;
			float sizeZ = lightEndZ - lightStartZ;
			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
				lightY,
				lightStartZ + random_double()*sizeZ);
			vec3 to_light = on_light - rec.p;
			float distance_squared = to_light.squared_length();
			to_light.make_unit_vector();
			if (dot(to_light, rec.normal) < 0) return emitted;
			float light_cosine = fabs(to_light.y());
			if (light_cosine < 0.000001) return emitted;

			// origin
			//pdf = distance_squared / (light_cosine * light_area);
			//scattered = ray(rec.p, to_light, r.time());

			// 1. cosine
			//cosine_pdf p(rec.normal);

			// 2. direction to light
			hittable *light_shape = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, 0);
			hittable_pdf p(light_shape, rec.p);			
			
			scattered = ray(rec.p, p.generate(), r.time());
			pdf = p.value(scattered.direction());
			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
				*color(scattered, world, depth + 1) / pdf;
		}
		else
			return emitted;
	}
	else
		return vec3(0, 0, 0);
}

渲染效果和 【光线追踪系列十七】直接光源采样 中的一样:

效果如下:
在这里插入图片描述

三、混合两种密度函数

正如开头所提到的,加权平均的pdf依然是一个pdf,将上述两种pdf的权值分别设置为0.5:

class mixture_pdf : public pdf {
    public:
        mixture_pdf(pdf *p0, pdf *p1 ) { p[0] = p0; p[1] = p1; }
        virtual float value(const vec3& direction) const {
            return 0.5 * p[0]->value(direction) + 0.5 *p[1]->value(direction);
        }
        virtual vec3 generate() const {
            if (random_double() < 0.5)
                return p[0]->generate();
            else
                return p[1]->generate();
        }
        pdf *p[2];
};

color函数修改如下:

vec3 color(const ray& r, hittable *world, int depth) {
	hit_record rec;
	if (world->hit(r, 0.001, FLT_MAX, rec)) {
		ray scattered;
		vec3 attenuation;
		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
		float pdf;
		vec3 albedo;
		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {

			//new start
			float sizeX = lightEndX - lightStartX;
			float sizeZ = lightEndZ - lightStartZ;
			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
				lightY,
				lightStartZ + random_double()*sizeZ);
			vec3 to_light = on_light - rec.p;
			float distance_squared = to_light.squared_length();
			to_light.make_unit_vector();
			if (dot(to_light, rec.normal) < 0) return emitted;
			float light_cosine = fabs(to_light.y());
			if (light_cosine < 0.000001) return emitted;

			// origin
			//pdf = distance_squared / (light_cosine * light_area);
			//scattered = ray(rec.p, to_light, r.time());

			// 1. cosine
			//cosine_pdf p(rec.normal);

			// 2. direction to light
			//hittable *light_shape = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, 0);
			//hittable_pdf p(light_shape, rec.p);			

			// 3.mixture density
			hittable *light_shape = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, 0);
			hittable_pdf p0(light_shape, rec.p);
			cosine_pdf p1(rec.normal);
			mixture_pdf p(&p0, &p1);
			
			scattered = ray(rec.p, p.generate(), r.time());
			pdf = p.value(scattered.direction());
			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
				*color(scattered, world, depth + 1) / pdf;
		}
		else
			return emitted;
	}
	else
		return vec3(0, 0, 0);
}

最终效果如下:

在这里插入图片描述

三、最终源码

/*
代码改动:
1.新增函数random_in_unit_disk()
2.修改相机类
3.修改入口函数
*/


#include "imgui.h"
#include "imgui_impl_glfw.h"
#include "imgui_impl_opengl2.h"
#include <stdio.h>
//#include <sys/time.h> // for gettimeofday()
#include <ctime>
//#include <unistd.h>
#include "vec3.h"
#include "draw.h"
#include <cstdlib>
#include <thread>
#include <fstream>
#include <vector>
#include <ctime>
#define STB_IMAGE_IMPLEMENTATION
#include "stb_image.h"

//#include <sys/time.h>
#include<string>
#ifdef __APPLE__
#define GL_SILENCE_DEPRECATION
#endif
#include <GLFW/glfw3.h>

#if defined(_MSC_VER) && (_MSC_VER >= 1900) && !defined(IMGUI_DISABLE_WIN32_FUNCTIONS)
#pragma comment(lib, "legacy_stdio_definitions")
#endif



static void glfw_error_callback(int error, const char *description)
{
	fprintf(stderr, "Glfw Error %d: %s\n", error, description);
}

using namespace std;

long curTime;
long totTime;
float timeRemaining;

//全局变量

ImVec4 clear_color = ImVec4(0.0f, 0.0f, 0.00f, 1.00f);
float progressDone = 0.0f, progressDir = 1.0f;
bool gRayTracingBegin = false;


int numPixelTotal;
int numThread = 5;
int numPixelRendered = 0;
int doneRecord = 0;
int framebufferInited = false;
int gFrameFinished = false;

int *lastFrameBuffer;

int gFov = 40;
float lightStartX = 213;
float lightEndX = 343;
float lightStartZ = 227;
float lightEndZ = 332;
float lightY = 554;

//vec3 lookfrom(22, 2, 3);
//vec3 lookat(0, 2.5, 0);
vec3 lookfrom(278, 278, -800);
vec3 lookat(278, 278, 0);
float dist_to_focus = 10.0;
float aperture = 0.0;
float M_PI = 3.1415926;

static int speedFactor = 1;

#pragma region Ray

class ray
{
public:
	vec3 A; //起点
	vec3 B; //方向
	float _time;
	ray() {}
	ray(const vec3 &a, const vec3 &b, float ti = 0.0f)
	{
		A = a;
		B = b;
		_time = ti;
	}
	float time() const { return _time; }
	vec3 origin() const { return A; }
	vec3 direction() const { return B; }
	vec3 point_at_parameter(float t) const { return A + t * B; } //终点的坐标
};

#pragma endregion

#pragma region BVH

inline float ffmin(float a, float b) { return a < b ? a : b; }
inline float ffmax(float a, float b) { return a > b ? a : b; }

class aabb {
public:
	vec3 _min; //左下角顶点
	vec3 _max; //右上角顶点
	aabb() {}
	aabb(const vec3& a, const vec3& b) { _min = a; _max = b; }

	vec3 min() const { return _min; }
	vec3 max() const { return _max; }

	//判断射线是否在t区间
	bool hit(const ray &r, float tmin, float tmax) const
	{
		for (int a = 0; a < 3; a++)
		{
			float invD = 1.0f / r.direction()[a];
			float t0 = (min()[a] - r.origin()[a]) * invD;
			float t1 = (max()[a] - r.origin()[a]) * invD;
			if (invD < 0.0f)
				std::swap(t0, t1);
			tmax = t1 < tmax ? t1 : tmax;   //F为两者终点的最小值
			tmin = t0 > tmin ? t0 : tmin;   //f为两者起点的最大值
			if (tmax <= tmin)       //F <= f
				return false;
		}
		return true;
	}
};

//计算2个box的包围盒
aabb surrounding_box(aabb box0, aabb box1)
{
	vec3 small(ffmin(box0.min().x(), box1.min().x()),
		ffmin(box0.min().y(), box1.min().y()),
		ffmin(box0.min().z(), box1.min().z()));
	vec3 big(ffmax(box0.max().x(), box1.max().x()),
		ffmax(box0.max().y(), box1.max().y()),
		ffmax(box0.max().z(), box1.max().z()));
	return aabb(small, big);
}

#pragma endregion

#pragma region obj

struct hit_record;

//class material
//{
//public:
//	//r_in为入射光线, scattered为散射光线, attenuation 意思为衰减量,实际为各通道的反射率
//	virtual bool scatter(const ray &r_in, const hit_record &rec, vec3 &attenuation, ray &scattered) const = 0;
//	virtual vec3 emitted(float u, float v, const vec3& p) const {
//		return vec3(0, 0, 0);
//	}
//};

class material {
public:
	virtual bool scatter(const ray& r_in,
		const hit_record& rec, vec3& albedo, ray& scattered, float& pdf) const {
		return false;
	}
	virtual float scattering_pdf(const ray& r_in, const hit_record& rec,
		const ray& scattered) const {
		return 0;
	}
	/*virtual vec3 emitted(float u, float v, const vec3& p) const {
		return vec3(0, 0, 0);
	}*/
	virtual vec3 emitted(const ray& r_in, const hit_record& rec, float u, float v, const vec3& p) const {
		return vec3(0, 0, 0);
	}
};


class texture {
public:
	//返回一个三通道的颜色值,p为命中终点坐标
	virtual vec3 value(float u, float v, const vec3& p) const = 0;
};

class constant_texture : public texture {
public:
	vec3 color;

	constant_texture() { }
	constant_texture(vec3 c) : color(c) { }

	//返回一个三通道的color
	virtual vec3 value(float u, float v, const vec3& p) const {
		return color;
	}
};

//棋盘纹理
class checker_texture : public texture {
public:
	texture *odd;
	texture *even;

	checker_texture() { }
	checker_texture(texture *t0, texture *t1) : even(t0), odd(t1) { }

	virtual vec3 value(float u, float v, const vec3& p) const {
		float sines = sin(10 * p.x())*sin(10 * p.y())*sin(10 * p.z());
		if (sines < 0)
			return odd->value(u, v, p);
		else
			return even->value(u, v, p);
	}
};


//记录命中信息
struct hit_record
{
	float t;     //命中射线的长度
	vec3 p;      //命中终点坐标
	vec3 normal; //命中点的法向量
	material *mat_ptr; //new
	float u = 1;
	float v = 1;
};

#define pi 3.1415926

//物体的虚基类
class hittable
{
public:
	virtual bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const = 0;
	virtual bool bounding_box(float t0, float t1, aabb& box) const = 0;
	virtual float pdf_value(const vec3& o, const vec3& v) const { return 0.0; }
	virtual vec3 random(const vec3& o) const { return vec3(1, 0, 0); }
};

#pragma endregion

#pragma region uv

class image_texture : public texture
{
public:
	unsigned char* data;
	int nx, ny;

	image_texture() {}
	image_texture(unsigned char* pixels, int A, int B) : data(pixels), nx(A), ny(B) {}

	//输入u和v,输出对应图片像素的rgb值
	virtual vec3 value(float u, float v, const vec3 &p) const
	{
		int i = int((u)* nx);//求出像素索引
		int j = int((1 - v)*ny - 0.001f);
		if (i < 0) i = 0;
		if (j < 0) j = 0;
		if (i > nx - 1) i = nx - 1;
		if (j > ny - 1) j = ny - 1;
		float r = int(data[3 * i + 3 * nx*j]) / 255.0f;
		float g = int(data[3 * i + 3 * nx*j + 1]) / 255.0f;
		float b = int(data[3 * i + 3 * nx*j + 2]) / 255.0f;
		return vec3(r, g, b);
	}
};

#pragma endregion
#pragma region bvH

int box_x_compare(const void * a, const void * b) {
	aabb box_left, box_right;
	hittable *ah = *(hittable**)a;
	hittable *bh = *(hittable**)b;

	if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right))
		std::cerr << "no bounding box in bvh_node constructor\n";

	if (box_left.min().x() - box_right.min().x() < 0.0)
		return -1;
	else
		return 1;
}

int box_y_compare(const void * a, const void * b) {
	aabb box_left, box_right;
	hittable *ah = *(hittable**)a;
	hittable *bh = *(hittable**)b;

	if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right))
		std::cerr << "no bounding box in bvh_node constructor\n";

	if (box_left.min().y() - box_right.min().y() < 0.0)
		return -1;
	else
		return 1;
}

int box_z_compare(const void * a, const void * b) {
	aabb box_left, box_right;
	hittable *ah = *(hittable**)a;
	hittable *bh = *(hittable**)b;

	if (!ah->bounding_box(0, 0, box_left) || !bh->bounding_box(0, 0, box_right))
		std::cerr << "no bounding box in bvh_node constructor\n";

	if (box_left.min().z() - box_right.min().z() < 0.0)
		return -1;
	else
		return 1;
}


class bvh_node : public hittable
{
public:
	bvh_node() {}

	hittable *left;
	hittable *right;
	aabb box;

	bool bounding_box(float t0, float t1, aabb &b) const
	{
		b = box;
		return true;
	}

	bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const
	{
		if (box.hit(r, t_min, t_max))
		{
			hit_record left_rec, right_rec;
			bool hit_left = left->hit(r, t_min, t_max, left_rec);
			bool hit_right = right->hit(r, t_min, t_max, right_rec);
			if (hit_left && hit_right)
			{
				if (left_rec.t < right_rec.t)
					rec = left_rec;
				else
					rec = right_rec;
				return true;
			}
			else if (hit_left)
			{
				rec = left_rec;
				return true;
			}
			else if (hit_right)
			{
				rec = right_rec;
				return true;
			}
			else
				return false;
		}
		else
			return false;
	}

	bvh_node(hittable **l, int n, float time0, float time1)
	{
		int axis = int(3 * random_double());

		if (axis == 0)
			qsort(l, n, sizeof(hittable *), box_x_compare);
		else if (axis == 1)
			qsort(l, n, sizeof(hittable *), box_y_compare);
		else
			qsort(l, n, sizeof(hittable *), box_z_compare);

		if (n == 1)
		{
			left = right = l[0];
		}
		else if (n == 2)
		{
			left = l[0];
			right = l[1];
		}
		else
		{
			left = new bvh_node(l, n / 2, time0, time1);
			right = new bvh_node(l + n / 2, n - n / 2, time0, time1);
		}

		aabb box_left, box_right;

		if (!left->bounding_box(time0, time1, box_left) ||
			!right->bounding_box(time0, time1, box_right))
		{

			std::cerr << "no bounding box in bvh_node constructor\n";
		}

		box = surrounding_box(box_left, box_right);
	}
};

#pragma endregion


#pragma region ONE

class onb
{
public:
	onb() {}
	inline vec3 operator[](int i) const { return axis[i]; }
	vec3 u() const { return axis[0]; }
	vec3 v() const { return axis[1]; }
	vec3 w() const { return axis[2]; }
	vec3 local(float a, float b, float c) const { return a * u() + b * v() + c * w(); }
	vec3 local(const vec3& a) const { return a.x()*u() + a.y()*v() + a.z()*w(); }
	void build_from_w(const vec3&);
	vec3 axis[3];
};


void onb::build_from_w(const vec3& n) {
	axis[2] = unit_vector(n);
	vec3 a;
	if (fabs(w().x()) > 0.9)
		a = vec3(0, 1, 0);
	else
		a = vec3(1, 0, 0);
	axis[1] = unit_vector(cross(w(), a));
	axis[0] = cross(w(), v());
}

#pragma endregion

#pragma region PDF

inline vec3 random_cosine_direction() {
	float r1 = random_double();
	float r2 = random_double();
	float z = sqrt(1 - r2);
	float phi = 2 * M_PI*r1;
	float x = cos(phi)*sqrt(r2);
	float y = sin(phi)*sqrt(r2);
	return vec3(x, y, z);
}

class pdf {
public:
	virtual float value(const vec3& direction) const = 0;	//当前方向的pdf值
	virtual vec3 generate() const = 0;	//生成对应该pdf分布的随机数
};

class cosine_pdf : public pdf {
public:
	cosine_pdf(const vec3& w) { uvw.build_from_w(w); }
	virtual float value(const vec3& direction) const {
		float cosine = dot(unit_vector(direction), uvw.w());
		if (cosine > 0)
			return cosine / M_PI;
		else
			return 0;
	}
	virtual vec3 generate() const {
		return uvw.local(random_cosine_direction());
	}
	onb uvw;
};

class hittable_pdf : public pdf {
public:
	hittable_pdf(hittable *p, const vec3& origin) : ptr(p), o(origin) {}
	virtual float value(const vec3& direction) const {
		return ptr->pdf_value(o, direction);
	}
	virtual vec3 generate() const {
		return ptr->random(o);
	}
	vec3 o;
	hittable *ptr;
};

class mixture_pdf : public pdf {
public:
	mixture_pdf(pdf *p0, pdf *p1) { p[0] = p0; p[1] = p1; }
	virtual float value(const vec3& direction) const {
		return 0.5 * p[0]->value(direction) + 0.5 *p[1]->value(direction);
	}
	virtual vec3 generate() const {
		if (random_double() < 0.5)
			return p[0]->generate();
		else
			return p[1]->generate();
	}
	pdf *p[2];
};

#pragma endregion


vec3 random_in_unit_sphere()
{
	vec3 p;
	do
	{
		p = 2.0 * vec3(random_double(), random_double(), random_double()) - vec3(1, 1, 1);
	} while (p.squared_length() >= 1.0);
	return unit_vector(p);
}

void DrawFrame(int *fb)
{

	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0.0, display_w, 0.0, display_h, 0.0, 1.0);
	glBegin(GL_POINTS);

	for (int i = 0; i < display_h; i++)
	{

		int r, g, b;
		for (int j = 0; j < display_w; j++)
		{
			Color2RGB(fb[i * display_w + j], r, g, b);
			glColor3f((float)r / 255, (float)g / 255, (float)b / 255);
			glVertex3f(j + 20, display_h - i - 300, 0);
		}
	}

	glEnd();
}
//记录进度和时间
void RecordProgressAndTime()
{
	progressDone = float(numPixelRendered) / (numPixelTotal);
	// totTime = (GetCurrentTimeMs() - curTime);
	// timeRemaining = totTime * (1 - progressDone) / progressDone / 1000; //根据过去的平均时间来计算剩余时间 x/t=pTodo/pDone -> x= t*pTodo/pDone = t*(1-pDone) /pDone
}

class camera
{
public:
	vec3 origin;
	vec3 lower_left_corner;
	vec3 horizontal;
	vec3 vertical;
	vec3 u, v, w;
	float lens_radius;
	float time0, time1;
	//lookfrom为相机位置,lookat为观察位置,vup传(0,1,0),vfov为视野角度,aspect为屏幕宽高比
	//aperture为光圈大小,focus_dist为相机到观察点的距离
	camera(vec3 lookfrom, vec3 lookat, vec3 vup, float vfov, float aspect, float aperture, float focus_dist, float t0, float t1)
	{
		time0 = t0;
		time1 = t1;
		lens_radius = aperture / 2;
		float theta = vfov * M_PI / 180;
		float half_height = tan(theta / 2);
		float half_width = aspect * half_height;
		origin = lookfrom;
		w = unit_vector(lookfrom - lookat);
		u = unit_vector(cross(vup, w));
		v = cross(w, u);
		lower_left_corner = origin - half_width * focus_dist * u - half_height * focus_dist * v - focus_dist * w;
		horizontal = 2 * half_width * focus_dist * u;
		vertical = 2 * half_height * focus_dist * v;
	}

	ray get_ray(float s, float t)
	{
		vec3 rd = lens_radius * random_in_unit_disk();
		vec3 offset = u * rd.x() + v * rd.y();
		float time = time0 + random_double() * (time1 - time0);
		return ray(origin + offset, lower_left_corner + s * horizontal + t * vertical - origin - offset, time);
	}

};

//输入命中点p的坐标,输出纹理坐标u,v
void get_sphere_uv(const vec3& p, float& u, float& v) {
	auto phi = atan2(p.z(), p.x());
	auto theta = asin(p.y());
	u = 1 - (phi + pi) / (2 * pi);
	v = (theta + pi / 2) / pi;
}

//aarect.h
class xy_rect : public hittable {
public:
	xy_rect() {}
	xy_rect(float _x0, float _x1, float _y0, float _y1, float _k, material *mat) : x0(_x0), x1(_x1), y0(_y0), y1(_y1), k(_k), mp(mat) {};
	virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
	virtual bool bounding_box(float t0, float t1, aabb& box) const {
		box = aabb(vec3(x0, y0, k - 0.0001), vec3(x1, y1, k + 0.0001));
		return true;
	}
	material  *mp;
	float x0, x1, y0, y1, k;
};

bool xy_rect::hit(const ray& r, float t0, float t1, hit_record& rec) const {
	float t = (k - r.origin().z()) / r.direction().z();
	if (t < t0 || t > t1)        return false;
	float x = r.origin().x() + t * r.direction().x();
	float y = r.origin().y() + t * r.direction().y();
	if (x < x0 || x > x1 || y < y0 || y > y1)         return false;
	rec.u = (x - x0) / (x1 - x0);
	rec.v = (y - y0) / (y1 - y0);
	rec.t = t;
	rec.mat_ptr = mp;
	rec.p = r.point_at_parameter(t);
	rec.normal = vec3(0, 0, 1);	//默认法线为(0,0,1),但会有潜在的朝向问题
	return true;
}

class xz_rect : public hittable {
public:
	xz_rect() {}
	xz_rect(float _x0, float _x1, float _z0, float _z1, float _k, material *mat) :x0(_x0), x1(_x1), z0(_z0), z1(_z1), k(_k), mp(mat) {};
	virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const override;

	virtual bool bounding_box(float time0, float time1, aabb& output_box) const override {
		// The bounding box must have non-zero width in each dimension, so pad the Y
		// dimension a small amount.
		output_box = aabb(vec3(x0, k - 0.0001, z0), vec3(x1, k + 0.0001, z1));
		return true;
	}
	virtual float  pdf_value(const vec3& o, const vec3& v) const {
		hit_record rec;
		if (this->hit(ray(o, v), 0.001, FLT_MAX, rec)) {
			float area = (x1 - x0)*(z1 - z0);
			float distance_squared = rec.t * rec.t * v.squared_length();
			float cosine = fabs(dot(v, rec.normal) / v.length());
			return  distance_squared / (cosine * area);
		}
		else
			return 0;
	}
	virtual vec3 random(const vec3& o) const {
		vec3 random_point = vec3(x0 + random_double()*(x1 - x0), k,
			z0 + random_double()*(z1 - z0));
		return random_point - o;
	}

	material  *mp;
	double x0, x1, z0, z1, k;
};

bool xz_rect::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
	auto t = (k - r.origin().y()) / r.direction().y();
	if (t < t_min || t > t_max)
		return false;
	auto x = r.origin().x() + t * r.direction().x();
	auto z = r.origin().z() + t * r.direction().z();
	if (x < x0 || x > x1 || z < z0 || z > z1)
		return false;
	rec.u = (x - x0) / (x1 - x0);
	rec.v = (z - z0) / (z1 - z0);
	rec.t = t;
	rec.normal = vec3(0, 1, 0);
	rec.mat_ptr = mp;
	rec.p = r.point_at_parameter(t);
	return true;
}

class yz_rect : public hittable {
public:

	yz_rect(float _y0, float _y1, float _z0, float _z1, float _k,material *mat): y0(_y0), y1(_y1), z0(_z0), z1(_z1), k(_k), mp(mat) {};

	virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const override;

	virtual bool bounding_box(float time0, float time1, aabb& output_box) const override {
		// The bounding box must have non-zero width in each dimension, so pad the X
		// dimension a small amount.
		output_box = aabb(vec3(k - 0.0001, y0, z0), vec3(k + 0.0001, y1, z1));
		return true;
	}

	material  *mp;
	double y0, y1, z0, z1, k;
};

bool yz_rect::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
	auto t = (k - r.origin().x()) / r.direction().x();
	if (t < t_min || t > t_max)
		return false;
	auto y = r.origin().y() + t * r.direction().y();
	auto z = r.origin().z() + t * r.direction().z();
	if (y < y0 || y > y1 || z < z0 || z > z1)
		return false;
	rec.u = (y - y0) / (y1 - y0);
	rec.v = (z - z0) / (z1 - z0);
	rec.t = t;
	rec.normal = vec3(1, 0, 0);

	rec.mat_ptr = mp;
	rec.p = r.point_at_parameter(t);
	return true;
}

class sphere : public hittable
{
public:
	vec3 center;
	float radius;
	material *mat_ptr; /* NEW */
	sphere() {}
	sphere(vec3 cen, float r, material *m) : center(cen), radius(r), mat_ptr(m) {}; //new


	//如果命中了,命中记录保存到rec
	virtual bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const
	{
		vec3 oc = r.origin() - center;
		float a = dot(r.direction(), r.direction());
		float b = dot(oc, r.direction());
		float c = dot(oc, oc) - radius * radius;
		float discriminant = b * b - a * c;

		if (discriminant > 0)
		{
			float temp = (-b - sqrt(discriminant)) / a; //小实数根
			if (temp < t_max && temp > t_min)
			{
				rec.t = temp;
				rec.p = r.point_at_parameter(rec.t);
				rec.normal = (rec.p - center) / radius;
				rec.mat_ptr = mat_ptr; /* NEW */
				vec3 outward_normal = (rec.p - center) / radius;
				get_sphere_uv(outward_normal, rec.u, rec.v);
				return true;
			}
			temp = (-b + sqrt(discriminant)) / a; //大实数根
			if (temp < t_max && temp > t_min)
			{
				rec.t = temp;
				rec.p = r.point_at_parameter(rec.t);
				rec.normal = (rec.p - center) / radius;
				rec.mat_ptr = mat_ptr; /* NEW */
				vec3 outward_normal = (rec.p - center) / radius;
				get_sphere_uv(outward_normal, rec.u, rec.v);
				return true;
			}
		}
		return false;
	}

	bool bounding_box(float t0, float t1, aabb &box) const
	{
		box = aabb(center - vec3(radius, radius, radius), center + vec3(radius, radius, radius));
		return true;
	}
};

class moving_sphere : public hittable
{
public:
	vec3 center0, center1;
	float time0, time1;
	float radius;
	material *mat_ptr; /* NEW */
	vec3 center(float time) const
	{
		return center0 + ((time - time0) / (time1 - time0)) * (center1 - center0);
	}
	moving_sphere() {}
	moving_sphere(
		vec3 cen0, vec3 cen1, double t0, double t1, double r, material *m)
		: center0(cen0), center1(cen1), time0(t0), time1(t1), radius(r), mat_ptr(m)
	{};

	//如果命中了,命中记录保存到rec
	virtual bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const
	{
		vec3 oc = r.origin() - center(r.time());
		float a = dot(r.direction(), r.direction());
		float b = dot(oc, r.direction());
		float c = dot(oc, oc) - radius * radius;
		float discriminant = b * b - a * c;
		if (discriminant > 0)
		{
			float temp = (-b - sqrt(discriminant)) / a; //小实数根
			if (temp < t_max && temp > t_min)
			{
				rec.t = temp;
				rec.p = r.point_at_parameter(rec.t);
				rec.normal = (rec.p - center(r.time())) / radius;
				rec.mat_ptr = mat_ptr; /* NEW */
				return true;
			}
			temp = (-b + sqrt(discriminant)) / a; //大实数根
			if (temp < t_max && temp > t_min)
			{
				rec.t = temp;
				rec.p = r.point_at_parameter(rec.t);
				rec.normal = (rec.p - center(r.time())) / radius;
				rec.mat_ptr = mat_ptr; /* NEW */
				return true;
			}
		}
		return false;
	}

	bool bounding_box(float t0, float t1, aabb &box) const
	{
		aabb box0(center(t0) - vec3(radius, radius, radius), center(t0) + vec3(radius, radius, radius));
		aabb box1(center(t1) - vec3(radius, radius, radius), center(t1) + vec3(radius, radius, radius));
		box = surrounding_box(box0, box1);
		return true;
	}
};

class hittable_list : public hittable
{
public:
	hittable **list;
	int list_size;

	hittable_list() {}
	hittable_list(hittable **l, int n)
	{
		list = l;
		list_size = n;
	}
	//如果命中了,命中记录保存到rec
	virtual bool hit(const ray &r, float t_min, float t_max, hit_record &rec) const
	{
		hit_record temp_rec;
		bool hit_anything = false;
		double closest_so_far = t_max; //记录目前最近的t值
		for (int i = 0; i < list_size; i++)
		{
			if (list[i]->hit(r, t_min, closest_so_far, temp_rec))
			{
				hit_anything = true;
				closest_so_far = temp_rec.t;
				rec = temp_rec; //只记录打到的最近的球
			}
		}
		return hit_anything;
	}

	bool bounding_box(float t0, float t1, aabb &box) const
	{
		if (list_size < 1)
			return false;
		aabb temp_box;
		bool first_true = list[0]->bounding_box(t0, t1, temp_box);
		if (!first_true)
			return false;
		else
			box = temp_box;
		for (int i = 1; i < list_size; i++)
		{
			if (list[i]->bounding_box(t0, t1, temp_box))
			{
				box = surrounding_box(box, temp_box);
			}
			else
				return false;
		}
		return true;
	}
};

class flip_normals : public hittable {
public:
	flip_normals(hittable *p) : ptr(p) {}
	virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
		if (ptr->hit(r, t_min, t_max, rec)) {
			rec.normal = -rec.normal;
			return true;
		}
		else
			return false;
	}
	virtual bool bounding_box(float t0, float t1, aabb& box) const {
		return ptr->bounding_box(t0, t1, box);
	}
	hittable *ptr;
};

class box : public hittable {
public:
	box() {}
	box(const vec3& p0, const vec3& p1, material *ptr);
	virtual bool hit(const ray& r, float t0, float t1, hit_record& rec) const;
	virtual bool bounding_box(float t0, float t1, aabb& box) const {
		box = aabb(pmin, pmax);
		return true;
	}
	vec3 pmin, pmax;
	hittable *list_ptr;
};

box::box(const vec3& p0, const vec3& p1, material *ptr) {
	pmin = p0;
	pmax = p1;
	hittable **list = new hittable*[6];
	list[0] = new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p1.z(), ptr);
	list[1] = new flip_normals(new xy_rect(p0.x(), p1.x(), p0.y(), p1.y(), p0.z(), ptr));
	list[2] = new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p1.y(), ptr);
	list[3] = new flip_normals(new xz_rect(p0.x(), p1.x(), p0.z(), p1.z(), p0.y(), ptr));
	list[4] = new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p1.x(), ptr);
	list[5] = new flip_normals(new yz_rect(p0.y(), p1.y(), p0.z(), p1.z(), p0.x(), ptr));
	list_ptr = new hittable_list(list, 6);
}

bool box::hit(const ray& r, float t0, float t1, hit_record& rec) const {
	return list_ptr->hit(r, t0, t1, rec);
}

class translate : public hittable {
public:
	translate(hittable *p, const vec3& displacement)
		: ptr(p), offset(displacement) {}
	virtual bool hit(
		const ray& r, float t_min, float t_max, hit_record& rec) const;
	virtual bool bounding_box(float t0, float t1, aabb& box) const;
	hittable *ptr;
	vec3 offset;
};

bool translate::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
	ray moved_r(r.origin() - offset, r.direction(), r.time());
	if (ptr->hit(moved_r, t_min, t_max, rec)) {
		rec.p += offset;
		return true;
	}
	else
		return false;
}

bool translate::bounding_box(float t0, float t1, aabb& box) const {
	if (ptr->bounding_box(t0, t1, box)) {
		box = aabb(box.min() + offset, box.max() + offset);
		return true;
	}
	else
		return false;
}

class rotate_y : public hittable {
public:
	rotate_y(hittable *p, float angle);
	virtual bool hit(const ray& r, float t_min, float t_max, hit_record& rec) const;
	virtual bool bounding_box(float t0, float t1, aabb& box) const {
		box = bbox; return hasbox;
	}
	hittable *ptr;
	float sin_theta;
	float cos_theta;
	bool hasbox;
	aabb bbox;
};

rotate_y::rotate_y(hittable *p, float angle) : ptr(p) {
	float radians = (M_PI / 180.) * angle;
	sin_theta = sin(radians);
	cos_theta = cos(radians);
	hasbox = ptr->bounding_box(0, 1, bbox);
	vec3 min(FLT_MAX, FLT_MAX, FLT_MAX);
	vec3 max(-FLT_MAX, -FLT_MAX, -FLT_MAX);
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 2; j++) {
			for (int k = 0; k < 2; k++) {
				float x = i * bbox.max().x() + (1 - i)*bbox.min().x();
				float y = j * bbox.max().y() + (1 - j)*bbox.min().y();
				float z = k * bbox.max().z() + (1 - k)*bbox.min().z();
				float newx = cos_theta * x + sin_theta * z;
				float newz = -sin_theta * x + cos_theta * z;
				vec3 tester(newx, y, newz);
				for (int c = 0; c < 3; c++)
				{
					if (tester[c] > max[c])
						max[c] = tester[c];
					if (tester[c] < min[c])
						min[c] = tester[c];
				}
			}
		}
	}
	bbox = aabb(min, max);
}

bool rotate_y::hit(const ray& r, float t_min, float t_max, hit_record& rec) const {
	vec3 origin = r.origin();
	vec3 direction = r.direction();
	origin[0] = cos_theta * r.origin()[0] - sin_theta * r.origin()[2];
	origin[2] = sin_theta * r.origin()[0] + cos_theta * r.origin()[2];
	direction[0] = cos_theta * r.direction()[0] - sin_theta * r.direction()[2];
	direction[2] = sin_theta * r.direction()[0] + cos_theta * r.direction()[2];
	ray rotated_r(origin, direction, r.time());
	if (ptr->hit(rotated_r, t_min, t_max, rec)) {
		vec3 p = rec.p;
		vec3 normal = rec.normal;
		p[0] = cos_theta * rec.p[0] + sin_theta * rec.p[2];	//命中点坐标需要旋转
		p[2] = -sin_theta * rec.p[0] + cos_theta * rec.p[2];
		normal[0] = cos_theta * rec.normal[0] + sin_theta * rec.normal[2];	//命中点法线也要旋转
		normal[2] = -sin_theta * rec.normal[0] + cos_theta * rec.normal[2];
		rec.p = p;
		rec.normal = normal;
		return true;
	}
	else
		return false;
}

//class lambertian : public material
//{
//public:
//	texture *albedo; //反射率
//	lambertian(texture *a) : albedo(a) {}
//	virtual bool scatter(const ray &r_in, const hit_record &rec, vec3 &attenuation, ray &scattered) const
//	{
//		vec3 s_world = rec.p + rec.normal + random_in_unit_sphere();
//		scattered = ray(rec.p, s_world - rec.p, r_in.time()); //scattered为散射光线
//		attenuation = albedo->value(rec.u, rec.v, rec.p);     //注意这是各通道的反射率!
//		//attenuation = vec3(1, 0, 0);     //注意这是各通道的反射率!
//		return true;
//	}
//};

class lambertian : public material {
    public:
        lambertian(texture *a) : albedo(a) {}
        float scattering_pdf(const ray& r_in,
            const hit_record& rec, const ray& scattered) const {
            float cosine = dot(rec.normal, unit_vector(scattered.direction()));
            if (cosine < 0)
                return 0;
            return cosine / M_PI;
        }

		//origin
        //bool scatter(const ray& r_in,
        //    const hit_record& rec, vec3& alb, ray& scattered, float& pdf) const {
        //    vec3 target = rec.p + rec.normal + random_in_unit_sphere();
        //    scattered = ray(rec.p, unit_vector(target-rec.p), r_in.time());
        //    alb = albedo->value(rec.u, rec.v, rec.p);
        //    pdf = dot(rec.normal, scattered.direction()) / M_PI;
        //    return true;
        //}

		//pdf
		//bool scatter(const ray& r_in,
		//	const hit_record& rec, vec3& alb, ray& scattered, float& pdf) const {
		//	vec3 direction;
		//	do {
		//		direction = random_in_unit_sphere();
		//	} while (dot(direction, rec.normal) < 0);
		//	scattered = ray(rec.p, unit_vector(direction), r_in.time());
		//	alb = albedo->value(rec.u, rec.v, rec.p);
		//	pdf = 0.5 / M_PI;
		//	return true;
		//}

		//random_cosine_direction
		bool scatter(const ray& r_in, const hit_record& rec, vec3& alb, ray& scattered, float& pdf) const
		{
			onb uvw;
			uvw.build_from_w(rec.normal);
			vec3 direction = uvw.local(random_cosine_direction());
			scattered = ray(rec.p, unit_vector(direction), r_in.time());
			alb = albedo->value(rec.u, rec.v, rec.p);
			pdf = dot(uvw.w(), scattered.direction()) / M_PI;
			return true;
		}


        texture *albedo;
};


class metal : public material
{
public:
	vec3 albedo;
	float fuzz;

	metal(const vec3 &a, float f) : albedo(a)
	{
		fuzz = f < 1 ? f : 1;
	}

	virtual bool scatter(const ray &r_in, const hit_record &rec, vec3 &attenuation, ray &scattered) const
	{
		vec3 v = unit_vector(r_in.direction());
		vec3 n = rec.normal;
		vec3 p = rec.p;
		vec3 r = reflect(v, n);
		vec3 offset = fuzz * random_in_unit_sphere();
		scattered = ray(p, r + offset);

		attenuation = albedo;
		return (dot(scattered.direction(), rec.normal) > 0);
	}
};

class dielectric : public material
{
public:
	dielectric(float ri) : ref_idx(ri) {} //n2/n1
	virtual bool scatter(const ray &r_in, const hit_record &rec, vec3 &attenuation, ray &scattered) const
	{
		vec3 outward_normal;
		vec3 reflected = reflect(r_in.direction(), rec.normal);
		float ni_over_nt;
		attenuation = vec3(1.0, 1.0, 1.0);
		vec3 refracted;

		float reflect_prob; //反射概率
		float cosine;

		if (dot(r_in.direction(), rec.normal) > 0) //从里到外,即入射向量在法向量另外一侧的情况
		{
			outward_normal = -rec.normal;   //对法向量取反
			ni_over_nt = ref_idx;
			cosine = ref_idx * dot(r_in.direction(), rec.normal) / r_in.direction().length();
			//不知道为什么这里要乘一个ref_idx
		}
		else    //从外到里,即入射向量在法向量同一侧
		{
			outward_normal = rec.normal;    //法向量不变
			ni_over_nt = 1.0 / ref_idx;
			cosine = -dot(r_in.direction(), rec.normal) / r_in.direction().length();
		}

		if (refract(r_in.direction(), outward_normal, ni_over_nt, refracted))
		{
			reflect_prob = schlick(cosine, ref_idx);
		}
		else    //若无折射,则全反射
		{
			reflect_prob = 1.0;
		}

		if (random_double() < reflect_prob)
		{
			scattered = ray(rec.p, reflected);
		}
		else
		{
			scattered = ray(rec.p, refracted);
		}

		return true;
	}

	float ref_idx;
};

//material.h
class diffuse_light : public material {
public:
	texture *emit;
	diffuse_light(texture *a) : emit(a) {}
	virtual bool scatter(const ray& r_in, const hit_record& rec, vec3& attenuation, ray& scattered) const {
		return false;
	}
	virtual vec3 emitted(const ray& r_in, const hit_record& rec, float u, float v, const vec3& p) const {
		
		if (dot(rec.normal, r_in.direction()) >= 0.0)
			return emit->value(u, v, p);
		else
			return vec3(0, 0, 0);
	}
	/*virtual vec3 emitted(float u, float v, const vec3& p) const {
		return emit->value(u, v, p);
	}*/
};


//发射一条射线,并采样该射线最终输出到屏幕的颜色值值
//vec3 color(const ray &r, hittable *world, int depth) {
//	hit_record rec;
//	if (world->hit(r, 0.001, FLT_MAX, rec)) //射线命中物体
//	{
//		ray scattered; //散射光线
//		vec3 attenuation; //其实是反射率!
//		if (depth < 66 && rec.mat_ptr->scatter(r, rec, attenuation, scattered))
//		{
//			return attenuation * color(scattered, world, depth + 1);
//		}
//		else
//		{
//			return vec3(0, 0, 0);
//		}
//	}
//	else
//	{
//		vec3 unit_direction = unit_vector(r.direction());
//		float t = 0.5 * (unit_direction.y() + 1.0);
//		return (1.0 - t) * vec3(1.0, 1.0, 1.0) + t * vec3(0.5, 0.7, 1.0);
//	}
//}

//vec3 color(const ray& r, hittable *world, int depth) {
//	hit_record rec;
//	if (world->hit(r, 0.001, FLT_MAX, rec)) {
//		ray scattered;
//		vec3 attenuation;
//		vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
//		if (depth < 50 && rec.mat_ptr->scatter(r, rec, attenuation, scattered))
//			return emitted + attenuation * color(scattered, world, depth + 1);
//		else
//			return emitted;
//	}
//	else
//		return vec3(0, 0, 0);
//}

//pdf
//vec3 color(const ray& r, hittable *world, int depth) {
//	hit_record rec;
//	if (world->hit(r, 0.001, FLT_MAX, rec)) {
//		ray scattered;
//		vec3 attenuation;
//		vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
//		float pdf;
//		vec3 albedo;
//		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {
//			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
//				*color(scattered, world, depth + 1) / pdf;
//		}
//		else
//			return emitted;
//	}
//	else
//		return vec3(0, 0, 0);
//}
//light
//vec3 color(const ray& r, hittable *world, int depth) {
//	hit_record rec;
//	if (world->hit(r, 0.001, FLT_MAX, rec)) {
//		ray scattered;
//		vec3 attenuation;
//		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
//		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
//		float pdf;
//		vec3 albedo;
//		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {
//
//			//new start
//			float sizeX = lightEndX - lightStartX;
//			float sizeZ = lightEndZ - lightStartZ;
//			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
//			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
//				lightY,
//				lightStartZ + random_double()*sizeZ);
//			vec3 to_light = on_light - rec.p;
//			float distance_squared = to_light.squared_length();
//			to_light.make_unit_vector();
//			if (dot(to_light, rec.normal) < 0) return emitted;
//			float light_cosine = fabs(to_light.y());
//			if (light_cosine < 0.000001) return emitted;
//			pdf = distance_squared / (light_cosine * light_area);
//			scattered = ray(rec.p, to_light, r.time());
//			//new end
//
//			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
//				*color(scattered, world, depth + 1) / pdf;
//		}
//		else
//			return emitted;
//	}
//	else
//		return vec3(0, 0, 0);
//}

vec3 color(const ray& r, hittable *world, int depth) {
	hit_record rec;
	if (world->hit(r, 0.001, FLT_MAX, rec)) {
		ray scattered;
		vec3 attenuation;
		//vec3 emitted = rec.mat_ptr->emitted(rec.u, rec.v, rec.p);
		vec3 emitted = rec.mat_ptr->emitted(r, rec, rec.u, rec.v, rec.p);
		float pdf;
		vec3 albedo;
		if (depth < 50 && rec.mat_ptr->scatter(r, rec, albedo, scattered, pdf)) {

			//new start
			float sizeX = lightEndX - lightStartX;
			float sizeZ = lightEndZ - lightStartZ;
			float light_area = (lightEndX - lightStartX)*(lightEndZ - lightStartZ);
			vec3 on_light = vec3(lightStartX + random_double()*sizeX,
				lightY,
				lightStartZ + random_double()*sizeZ);
			vec3 to_light = on_light - rec.p;
			float distance_squared = to_light.squared_length();
			to_light.make_unit_vector();
			if (dot(to_light, rec.normal) < 0) return emitted;
			float light_cosine = fabs(to_light.y());
			if (light_cosine < 0.000001) return emitted;

			// origin
			//pdf = distance_squared / (light_cosine * light_area);
			//scattered = ray(rec.p, to_light, r.time());

			// 1. cosine
			//cosine_pdf p(rec.normal);

			// 2. direction to light
			//hittable *light_shape = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, 0);
			//hittable_pdf p(light_shape, rec.p);			

			// 3.mixture density
			hittable *light_shape = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, 0);
			hittable_pdf p0(light_shape, rec.p);
			cosine_pdf p1(rec.normal);
			mixture_pdf p(&p0, &p1);
			
			scattered = ray(rec.p, p.generate(), r.time());
			pdf = p.value(scattered.direction());
			return emitted + albedo * rec.mat_ptr->scattering_pdf(r, rec, scattered)
				*color(scattered, world, depth + 1) / pdf;
		
		}
		else
			return emitted;
	}
	else
		return vec3(0, 0, 0);
}

class perlin {
public:
	float noise(const vec3& p) const {
		int i = (int)floor(p.x()) & 255;
		int j = (int)floor(p.y()) & 255;
		int k = (int)floor(p.z()) & 255;
		// return ranfloat[perm_x[i] ^ perm_y[j] ^ perm_z[k]];
		return get_ranfloat(i, j, k);
	}

	float get_ranfloat(int x, int y, int z) const {
		return ranfloat[perm_x[x & 255] ^ perm_y[y & 255] ^ perm_z[z & 255]];
	}
	static float *ranfloat;
	static int *perm_x;
	static int *perm_y;
	static int *perm_z;
};

static float* perlin_generate() {
	float * p = new float[256];
	for (int i = 0; i < 256; ++i)
		p[i] = random_double();
	return p;
}

void permute(int *p, int n) {
	for (int i = n - 1; i > 0; i--) {
		int target = int(random_double()*(i + 1));
		int tmp = p[i];
		p[i] = p[target];
		p[target] = tmp;
	}
	return;
}

static int* perlin_generate_perm() {
	int * p = new int[256];
	for (int i = 0; i < 256; i++)
		p[i] = i;
	permute(p, 256);
	return p;
}

float *perlin::ranfloat = perlin_generate();
int *perlin::perm_x = perlin_generate_perm();
int *perlin::perm_y = perlin_generate_perm();
int *perlin::perm_z = perlin_generate_perm();

class noise_texture : public texture {
public:
	noise_texture() {}
	virtual vec3 value(float u, float v, const vec3& p) const {
		return vec3(1, 1, 1) * noise.noise(p);
	}
	perlin noise;
};

//hittable *random_scene() {
//	int n = 500;
//	texture *pertext = new noise_texture();
//	hittable **list = new hittable*[n + 1];
//	list[0] = new sphere(vec3(0, -1000, 0), 1000, new lambertian(new constant_texture(vec3(0.5, 0.5, 0.5))));
//	texture *checker = new checker_texture(new constant_texture(vec3(0.2, 0.3, 0.1)), new constant_texture(vec3(0.9, 0.9, 0.9)));
//	list[0] = new sphere(vec3(0, -1000, 0), 1000, new lambertian(checker));
//	int i = 1;
//	//list[0] = new sphere(vec3(0, -1000, 0), 1000, new lambertian(pertext));
//	//list[i++] = new sphere(vec3(5, 1, 0), 1.0, new lambertian(pertext));
//	for (int a = -11; a < 11; a++) {
//		for (int b = -11; b < 11; b++) {
//			float choose_mat = random_double();
//			vec3 center(a + 0.9*random_double(), 0.2, b + 0.9*random_double());
//			if ((center - vec3(4, 0.2, 0)).length() > 0.9) {
//				if (choose_mat < 0.8) {  // diffuse
//					if (b % 2 == 0) //动态模糊的球体
//					{
//						auto center2 = center + vec3(0, random_double(), 0);
//						list[i++] = new moving_sphere(center, center2, 0.0, 1.0, 0.2,
//							new lambertian(new constant_texture(vec3(random_double()*random_double(),
//								random_double()*random_double(),
//								random_double()*random_double()))
//							)
//						);
//					}
//					else
//					{
//						list[i++] = new sphere(center, 0.2,
//							new lambertian(new constant_texture(vec3(random_double()*random_double(),
//								random_double()*random_double(),
//								random_double()*random_double()))
//							)
//						);
//					}
//				}
//				else if (choose_mat < 0.95) { // metal
//					list[i++] = new sphere(center, 0.2,
//						new metal(vec3(0.5*(1 + random_double()),
//							0.5*(1 + random_double()),
//							0.5*(1 + random_double())),
//							0.5*random_double()));
//				}
//				else {  // glass
//					list[i++] = new sphere(center, 0.2, new dielectric(1.5));
//				}
//			}
//		}
//	}
//
//	int nx, ny, nn;
//	unsigned char *earthmapjpg = stbi_load("F://earthmap.jpg", &nx, &ny, &nn, 0);
//	material *earthmapJpg = new lambertian(new image_texture(earthmapjpg, nx, ny));
//	unsigned char *Cristiano = stbi_load("F://Cristiano.jpg", &nx, &ny, &nn, 0);
//	material *CristianoJpg = new lambertian(new image_texture(Cristiano, nx, ny));
//	unsigned char *earthmappng = stbi_load("F://earthmap.png", &nx, &ny, &nn, 0);
//	material *earthmapPng = new lambertian(new image_texture(earthmappng, nx, ny));
//
//	list[i++] = new sphere(vec3(5, 1, 0), 1.0, earthmapJpg);
//	list[i++] = new sphere(vec3(0, 1, 0), 1.0, CristianoJpg);
//	list[i++] = new sphere(vec3(-5, 1, 0), 1.0, earthmapPng);
//	list[i++] = new sphere(vec3(5, 0.35, 4), 0.35, new dielectric(1.2));
//	list[i++] = new sphere(vec3(0, 0.4, 3), 0.4, new lambertian(new constant_texture(vec3(0.0, 1.0, 1.0))));
//	list[i++] = new sphere(vec3(-6, 0.5, 2), 0.5, new metal(vec3(0.8, 0.8, 0.8), 0.1));
//
//	//return new hittable_list(list, i);
//	return new bvh_node(list, i, 0.0, 1.0);
//}

hittable *random_scene1() {
	texture *pertext = new noise_texture();
	hittable **list = new hittable*[4];
	int i = 0;
	list[i++] = new sphere(vec3(0, -1000, 0), 1000, new lambertian(new constant_texture(vec3(0.7, 0.5, 0.3))));
	
	int nx, ny, nn;
	unsigned char *Cristiano = stbi_load("F://Cristiano.jpg", &nx, &ny, &nn, 0);
	material *CristianoJpg = new lambertian(new image_texture(Cristiano, nx, ny));
	//list[i++] = new sphere(vec3(0, 2, 0), 2, new lambertian(pertext));
	list[i++] = new sphere(vec3(0, 2, 0), 2, CristianoJpg);
	list[i++] = new sphere(vec3(0, 6, 0), 1, new diffuse_light(new constant_texture(vec3(4, 4, 4))));
	list[i++] = new xy_rect(3, 5, 1, 3, -2, new diffuse_light(new constant_texture(vec3(4, 4, 4))));
	//return new hittable_list(list, i);
	return new bvh_node(list, i, 0.0, 1.0);
}

hittable *random_scene() {

	//cam.SetCamParams(vec3(278, 278, -800), vec3(278, 278, 0), vec3(0, 1, 0), 40, float(nx) / float(ny), 0.0, 10.0, 0.0, 1.0);
	hittable **list = new hittable*[8];
	int i = 0;
	material *red = new lambertian(new constant_texture(vec3(0.65, 0.05, 0.05)));
	material *white = new lambertian(new constant_texture(vec3(0.73, 0.73, 0.73)));
	material *green = new lambertian(new constant_texture(vec3(0.12, 0.45, 0.15)));
	material *light = new diffuse_light(new constant_texture(vec3(15, 15, 15)));

	list[i++] = new flip_normals(new yz_rect(0, 555, 0, 555, 555, green));
	list[i++] = new yz_rect(0, 555, 0, 555, 0, red);
	//list[i++] = new xz_rect(213, 343, 227, 332, 554, light);
	list[i++] = new xz_rect(lightStartX, lightEndX, lightStartZ, lightEndZ, lightY, light);
	list[i++] = new flip_normals(new xz_rect(0, 555, 0, 555, 555, white));
	list[i++] = new xz_rect(0, 555, 0, 555, 0, white);
	list[i++] = new flip_normals(new xy_rect(0, 555, 0, 555, 555, white));
	//list[i++] = new box(vec3(130, 0, 65), vec3(295, 165, 230), white);
	//list[i++] = new box(vec3(265, 0, 295), vec3(430, 330, 460), white);

	hittable *box1 = new box(vec3(0, 0, 0), vec3(165, 165, 165), white);
	hittable *box2 = new box(vec3(0, 0, 0), vec3(165, 330, 165), white);
	list[i++] = new translate(new rotate_y(box1, -18), vec3(130, 0, 65));
	list[i++] = new translate(new rotate_y(box2, 15), vec3(265, 0, 295));

	//int nn;
	//unsigned char *Cristiano = stbi_load("F://Cristiano.jpg", &nx, &ny, &nn, 0);
	//material *CristianoJpg = new lambertian(new image_texture(Cristiano, nx, ny));
	list[i++] = new sphere(vec3(0, 2, 0), 2, new lambertian(pertext));
	//list[i++] = new sphere(vec3(200, 130, 330), 130, CristianoJpg);
	return new bvh_node(list, i, 0.0, 1.0);
}





hittable *world;

void RayTracingInOneThread(int k)
{
	float R = cos(M_PI / 4);
	// list[0] = new sphere(vec3(-R, 0, -1), R, new lambertian(vec3(0, 0, 1)));
	// list[1] = new sphere(vec3(R, 0, -1), R, new lambertian(vec3(1, 0, 0)));
	// world = new hittable_list(list, 2);

	// list[0] = new sphere(vec3(0, 0, -1), 0.5, new lambertian(vec3(0.8, 0.3, 0.3)));
	// list[1] = new sphere(vec3(0, -100.5, -1), 100, new lambertian(vec3(0.1/2, 0.2/2, 0.5/2)));
	// list[2] = new sphere(vec3(1, 0, -1), 0.5, new metal(vec3(0.8, 0.6, 0.2), 0.3));    
	// list[3] = new sphere(vec3(-1, 0, -1), 0.5, new dielectric(1.5));
	// world = new hittable_list(list, 4); 

	// camera cam(vec3(-2, 2, 1), vec3(-0.5, 0, -1), vec3(0, 1, 0), gFov, float(nx) / float(ny));

	camera cam(lookfrom, lookat, vec3(0, 1, 0), 40, float(nx) / float(ny), aperture, dist_to_focus, 0, 1.0);

	for (int j = ny - k; j >= 0; j -= numThread)
	{
		for (int i = 0; i < nx; i++)
		{
			RecordProgressAndTime();
			vec3 col(0, 0, 0);
			for (int s = 0; s < ns; s++)
			{
				float u = float(i + random_double()) / float(nx);
				float v = float(j + random_double()) / float(ny);
				ray r = cam.get_ray(u, v);
				col += color(r, world, 0);
			}
			col /= float(ns);
			col = vec3(sqrt(col[0]), sqrt(col[1]), sqrt(col[2]));

			int ir = int(255.99 * col[0]);
			int ig = int(255.99 * col[1]);
			int ib = int(255.99 * col[2]);
			DrawPixel(i, j, ir, ig, ib);
			numPixelRendered++;
		}
	}
}

void RayTracing()
{
	nx = 600;
	ny = 400;
	ns = 1000;
	// hittable *world = random_scene();
	world = random_scene();
	numPixelTotal = nx * ny;
	while (true)
	{
		while (!gRayTracingBegin)
		{
			//wait until begin
		}

		// curTime = GetCurrentTimeMs();
		gFrameFinished = false;
		if (!framebufferInited)
		{
			framebuffer = new int[display_w * display_h];
			lastFrameBuffer = new int[display_w * display_h];

			for (int i = 0; i < display_h * display_w; i++)
			{
				lastFrameBuffer[i] = 0;
			}
			framebufferInited = true;
		}
		numPixelRendered = 0;

		ofstream outFile("output_" + to_string(nx) + "x" + to_string(ny) + ".ppm");
		outFile << "P3\n"
			<< nx << " " << ny << "\n255\n";
		vector<thread> threads;

		for (int k = 0; k < numThread; k++)
		{
			threads.push_back(thread(RayTracingInOneThread, k));
		}

		for (auto &thread : threads)
		{
			thread.join();
		}
		gFrameFinished = true;
		Framebuffer2File(nx, ny, ns, framebuffer, outFile, progressDone);
		for (int i = 0; i < display_h*display_w; i++)
		{
			lastFrameBuffer[i] = framebuffer[i];
		}
		gRayTracingBegin = false;
	}

	return;
}

vec3 random_on_unit_sphere() {
	vec3 p;
	do {
		p = 2.0*vec3(random_double(), random_double(), random_double()) - vec3(1, 1, 1);
	} while (dot(p, p) >= 1.0);
	return unit_vector(p);
}


inline float pdf(const vec3& p) {
	return 1 / (4 * M_PI);
}

int main(int, char **)
{

	int N = 1000000;
	float sum=0;
	for (int i = 0; i < N; i++) {
		vec3 d = random_on_unit_sphere();
		float cosine_squared = d.z()*d.z();
		sum += cosine_squared / pdf(d);
}
	std::cout << "I =" << sum / N << "\n";
	totTime = 0;
#ifdef __APPLE__
	glfwWindowHint(GLFW_OPENGL_FORWARD_COMPAT, GL_TRUE);
#endif

	thread t(RayTracing);

	// Setup window
	glfwSetErrorCallback(glfw_error_callback);
	if (!glfwInit())
		return 1;
	GLFWwindow *window = glfwCreateWindow(1800, 900, "Ray Tracing", NULL, NULL);
	if (window == NULL)
		return 1;
	glfwMakeContextCurrent(window);
	glfwSwapInterval(1); // Enable vsync

	// Setup Dear ImGui context
	IMGUI_CHECKVERSION();
	ImGui::CreateContext();
	ImGuiIO &io = ImGui::GetIO();
	(void)io;

	// Setup Dear ImGui style
	ImGui::StyleColorsDark();

	// Setup Platform/Renderer bindings
	ImGui_ImplGlfw_InitForOpenGL(window, true);
	ImGui_ImplOpenGL2_Init();

	// Our state
	bool show_demo_window = true;
	bool show_another_window = false;

	// Main loop
	while (!glfwWindowShouldClose(window))
	{

		glfwPollEvents();

		// Start the Dear ImGui frame
		ImGui_ImplOpenGL2_NewFrame();
		ImGui_ImplGlfw_NewFrame();
		ImGui::NewFrame();

		// 1. Show the big demo window (Most of the sample code is in ImGui::ShowDemoWindow()! You can browse its code to learn more about Dear ImGui!).
		if (show_demo_window)
			ImGui::ShowDemoWindow(&show_demo_window);

		// 2. Show a simple window that we create ourselves. We use a Begin/End pair to created a named window.
		{
			static float f = 0.0f;
			static int counter = 0;

			ImGui::Begin("Ray Tracing"); // Create a window called "Hello, world!" and append into it.

			if (ImGui::Button("Start")) // Buttons return true when clicked (most widgets return true when edited/activated)
			{
				gRayTracingBegin = true;
			}

			ImGui::ProgressBar(progressDone, ImVec2(0.0f, 0.0f));
			// ImGui::SameLine(0.0f, ImGui::GetStyle().ItemInnerSpacing.x);
			ImGui::InputInt(" speed", &speedFactor);
			// ImGui::SameLine();
			ImGui::Text("Thread Num: %d ", numThread);
			ImGui::Text("Image Size:  %d x %d ", nx, ny);
			ImGui::Text("Camera fov:  %d ", gFov);
			ImGui::Text("Camera aperture: %.3f ", (float)aperture);
			ImGui::Text("Camera dist to focus: %.3f ", (float)dist_to_focus);

			// ImGui::Text("Progress Bar");
			ImGui::Text("Total time %.3f s", (float)totTime / 1000);
			ImGui::Text("timeRemaining time %.3f s", timeRemaining);
			ImGui::End();
		}

		// 3. Show another simple window.
		if (show_another_window)
		{
			ImGui::Begin("Another Window", &show_another_window); // Pass a pointer to our bool variable (the window will have a closing button that will clear the bool when clicked)
			ImGui::Text("Hello from another window!");
			if (ImGui::Button("Close Me"))
				show_another_window = false;
			ImGui::End();
		}

		ImGui::Render();

		glfwGetFramebufferSize(window, &display_w, &display_h);

		glViewport(0, 0, display_w, display_h);
		glClearColor(0.0f, 0.0f, 0.00f, 1.00f);
		glClear(GL_COLOR_BUFFER_BIT);

		if (framebufferInited)
		{
			if (gFrameFinished)
			{
				DrawFrame(framebuffer);
			}
			else
			{
				// DrawFrame(lastFrameBuffer);
				DrawFrame(framebuffer);
			}
			// DrawFrame(framebuffer);
		}


		ImGui_ImplOpenGL2_RenderDrawData(ImGui::GetDrawData());

		glfwMakeContextCurrent(window);

		glfwSwapBuffers(window);
	}

	// Cleanup
	ImGui_ImplOpenGL2_Shutdown();
	ImGui_ImplGlfw_Shutdown();
	ImGui::DestroyContext();

	glfwDestroyWindow(window);
	glfwTerminate();

	return 0;
}

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值