fork download
copy 
  1. #include <Eigen\Dense>
  2. #include <Eigen\StdVector>
  3.  
  4. #include <time.h>		
  5. #include <iostream>
  6.  
  7. using Vec3f = Eigen::Vector3f;
  8.  
  9. class RandomLCG
  10. {
  11. private:
  12. 	unsigned int m_seed;
  13.  
  14. public:
  15. 	RandomLCG(unsigned seed = 0) :
  16. 		m_seed(seed)
  17. 	{}
  18.  
  19. 	float operator()()
  20. 	{
  21. 		m_seed = 214013u * m_seed + 2531011u;
  22. 		return m_seed * (1.0f / 4294967296.0f);
  23. 	}
  24. };
  25.  
  26. class Ray
  27. {
  28. public:
  29. 	Vec3f origin;
  30. 	Vec3f direction;
  31.  
  32. 	EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  33.  
  34. 		Ray(Vec3f const & ori, Vec3f const & dir) :
  35. 		origin(ori), direction(dir.normalized()) {}
  36. };
  37.  
  38. class Material
  39. {
  40. public:
  41. 	enum class Type { DIFFUSE, SPECULAR, REFRACTION };
  42.  
  43. 	Type type;
  44. 	float fresnel_ratio;
  45. };
  46.  
  47. class Sphere
  48. {
  49. public:
  50. 	Vec3f m_center;
  51. 	float m_radius;
  52.  
  53. 	Vec3f m_color;
  54. 	Vec3f m_emission;
  55. 	Material m_mateiral;
  56.  
  57. 	float m_sqr_radius;
  58. 	float m_max_color_component;
  59. 	Vec3f m_color_ratio;
  60.  
  61. 	EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  62.  
  63. public:
  64. 	Sphere(Vec3f && center, float rad,
  65. 		Vec3f && emission, Vec3f && color, Material && material) :
  66. 		m_center(center), m_color(color), m_emission(emission), m_radius(rad), m_mateiral(material)
  67. 	{
  68. 		m_sqr_radius = m_radius * m_radius;
  69. 		m_max_color_component = m_color.maxCoeff();
  70. 		m_color_ratio = (1.0f / m_max_color_component) * m_color;
  71. 	}
  72.  
  73. 	bool intersect(Ray const & ray, float & root, float threshold = 1e-20) const
  74. 	{
  75. 		root = 0.0f;
  76. 		// Solve t^2*d.d + 2*t*(o-p).d + (o-p).(o-p)-R^2 = 0
  77. 		Vec3f op = m_center - ray.origin;
  78. 		float b = op.transpose() * ray.direction;
  79. 		float det = b * b - op.squaredNorm() + m_sqr_radius;
  80. 		float const eps = 1e-4;
  81.  
  82. 		if (det < 0) return false;
  83.  
  84. 		float dets = (std::sqrt)(det);
  85.  
  86. 		if (b - dets > eps)
  87. 		{
  88. 			root = b - dets;
  89. 			return true;
  90. 		}
  91. 		else if (b + dets > eps)
  92. 		{
  93. 			root = b + dets;
  94. 			return true;
  95. 		}
  96.  
  97. 		return false;
  98. 	}
  99. };
  100.  
  101. class Intersection
  102. {
  103. public:
  104. 	enum class Type { GO_INWARD, GO_OUTWARD };
  105.  
  106. 	float root;
  107. 	Vec3f position;
  108. 	Vec3f normal;
  109. 	Vec3f normal_oriented;
  110. 	Type type;
  111. 	Vec3f color;
  112. 	Material material;
  113. 	Sphere const * obj;
  114.  
  115. 	EIGEN_MAKE_ALIGNED_OPERATOR_NEW
  116. };
  117.  
  118. class Scene
  119. {
  120. private:
  121. 	std::vector<Sphere, Eigen::aligned_allocator<Sphere> > m_components;
  122.  
  123. public:
  124. 	void add(Sphere const & sphere)
  125. 	{
  126. 		m_components.push_back(std::move(sphere));
  127. 	}
  128.  
  129. 	bool const & intersect(Ray const & ray, Intersection & intersection) const
  130. 	{
  131. 		bool is_intersect = false;
  132. 		float root = (std::numeric_limits<float>::max)();
  133. 		Sphere const * pivot = nullptr;
  134.  
  135. 		for (Sphere const & com : m_components)
  136. 		{
  137. 			float cur_root = 0.0f;
  138. 			if (!com.intersect(ray, cur_root)) continue;
  139. 			if (root <= cur_root) continue;
  140.  
  141. 			pivot = &com;
  142. 			root = cur_root;
  143. 			is_intersect = true;
  144. 		}
  145.  
  146. 		if (is_intersect)
  147. 		{
  148. 			//std::cout << "intersected" << std::endl;
  149. 			intersection.root = root;
  150. 			intersection.position = ray.origin + root * ray.direction;
  151. 			intersection.normal = (intersection.position - pivot->m_center).normalized();
  152. 			intersection.color = pivot->m_color;
  153. 			intersection.type = (intersection.normal.transpose() * ray.direction > 0) ?
  154. 				Intersection::Type::GO_INWARD : Intersection::Type::GO_OUTWARD;
  155. 			intersection.normal_oriented = intersection.normal * (intersection.type == Intersection::Type::GO_INWARD ? 1.0f : -1.0f);
  156. 			intersection.obj = pivot;
  157. 			intersection.material = pivot->m_mateiral;
  158.  
  159. 		}
  160.  
  161. 		//std::cout << "return value " << is_intersect << std::endl;
  162. 		return is_intersect;
  163. 	}
  164. };
  165.  
  166. Vec3f radiance(Scene const & scene, Ray const & ray, int const depth,
  167. 	int const min_depth, int const max_depth, RandomLCG & rand)
  168. {
  169. 	//std::cout << "depth " << depth << std::endl;
  170. 	Intersection inte;
  171. 	bool is_intersect = scene.intersect(ray, inte);
  172. 	//std::cout << "intersect test " << is_intersect << std::endl;
  173.  
  174. 	if (!is_intersect)
  175. 		return Vec3f::Zero();
  176.  
  177. 	Sphere const * sphere = inte.obj;
  178.  
  179. 	int const new_depth = depth + 1;
  180. 	bool reach_max_depth = (new_depth >= max_depth);
  181. 	bool reach_min_depth = (new_depth >= min_depth);
  182.  
  183. 	/* use Russian roulette for path termination */
  184. 	bool isRR = rand() < sphere->m_max_color_component;
  185.  
  186. 	if (reach_max_depth || (reach_min_depth && !isRR))
  187. 		return sphere->m_color;
  188.  
  189. 	/* WARNING: remain consideration */
  190.  
  191. 	Vec3f emission = inte.obj->m_emission;
  192.  
  193. 	switch (inte.material.type)
  194. 	{
  195. 	case Material::Type::DIFFUSE:
  196. 	{
  197. 		//std::cout << "DIFFUSE" << std::endl;
  198. 		/* local coordinate */
  199. 		Vec3f w = inte.normal_oriented;
  200. 		Vec3f wo = w.x() < -0.1f || w.x() > 0.1f ? Vec3f(0.0f, 1.0f, 0.0f) : Vec3f(1.0f, 0.0f, 0.0f);
  201. 		Vec3f u = (wo.cross(w)).normalized();
  202. 		Vec3f v = w.cross(u);
  203.  
  204. 		/* semi-sphere sampling */
  205. 		float r = 2.0f * M_PI * rand();
  206. 		float d = rand();
  207. 		float dsqt = std::sqrt(d);
  208.  
  209. 		Vec3f diffuse_dir = (u * std::cos(r) * dsqt + v * std::sin(r) * dsqt + w * std::sqrt(1.0f - d)).normalized();
  210. 		Vec3f diffuse = inte.color.cwiseProduct(radiance(scene,
  211. 			Ray{ inte.position, diffuse_dir }, new_depth, min_depth, max_depth, rand));
  212.  
  213. 		return emission + diffuse;
  214. 	}
  215. 	break;
  216. 	case Material::Type::SPECULAR:
  217. 	{
  218. 		//std::cout << "SPECULAR" << std::endl;
  219. 		Vec3f reflect_dir = ray.direction - inte.normal * (2.0f * inte.normal.transpose() * ray.direction);
  220. 		Vec3f spectular = inte.color.cwiseProduct(radiance(scene,
  221. 			Ray{ inte.position, reflect_dir }, new_depth, min_depth, max_depth, rand));
  222. 		return emission + spectular;
  223. 	}
  224. 	break;
  225. 	case Material::Type::REFRACTION:
  226. 	{
  227. 		//std::cout << "REFRACTION" << std::endl;
  228. 		Vec3f reflect_dir = ray.direction - inte.normal * (2.0f * inte.normal.transpose() * ray.direction);
  229.  
  230. 		bool go_inward = (inte.type == Intersection::Type::GO_INWARD);
  231. 		float refr_ratio = (go_inward ? inte.material.fresnel_ratio : 1.0f / inte.material.fresnel_ratio);
  232.  
  233. 		float ddn = ray.direction.transpose() * inte.normal_oriented;
  234. 		float cos2t = 1.0f - refr_ratio * refr_ratio * (1.0f - ddn * ddn);
  235.  
  236. 		/* total internal reflection */
  237. 		if (cos2t < 0.0f)
  238. 		{
  239. 			Vec3f reflection = inte.color.cwiseProduct(radiance(scene,
  240. 				Ray{ inte.position, reflect_dir }, new_depth, min_depth, max_depth, rand));
  241. 			return emission + reflection;
  242. 		}
  243.  
  244. 		/* compute refraction direction */
  245. 		Vec3f refraction_dir =
  246. 			(ray.direction * refr_ratio - inte.normal * (go_inward ? 1.0f : -1.0f) * (ddn * refr_ratio + std::sqrt(cos2t))).normalized();
  247.  
  248. 		float R0 = std::abs(refr_ratio - 1.0f) / (refr_ratio + 1.0f);
  249. 		float c = 1 - (go_inward ? -ddn : refraction_dir.transpose() * inte.normal);
  250. 		float fresnel_reflectance = R0 + (1.0f - R0) * c * c * c * c;
  251. 		float P = 0.25f + 0.5f * fresnel_reflectance;
  252.  
  253. 		Vec3f res;
  254. 		if (new_depth > 2)
  255. 		{
  256. 			if (rand() < P)
  257. 				res = (fresnel_reflectance / P) *
  258. 				radiance(scene, Ray{ inte.position, reflect_dir }, new_depth, min_depth, max_depth, rand);
  259. 			else
  260. 				res = ((1.0f - fresnel_reflectance) / (1.0f - P)) *
  261. 				radiance(scene, Ray{ inte.position, refraction_dir }, new_depth, min_depth, max_depth, rand);
  262. 		}
  263. 		else
  264. 			res =
  265. 			fresnel_reflectance * radiance(scene, Ray{ inte.position, reflect_dir }, new_depth, min_depth, max_depth, rand) +
  266. 			(1.0f - fresnel_reflectance) * radiance(scene, Ray{ inte.position, refraction_dir }, new_depth, min_depth, max_depth, rand);
  267.  
  268. 		return emission + inte.color.cwiseProduct(res);
  269. 	}
  270. 	break;
  271. 	default:
  272. 		std::cout << "ERROR : reach invalid code section" << std::endl;
  273. 		return Vec3f::Zero();
  274. 		break;
  275. 	}
  276. }
  277.  
  278. inline float clamp(float x) {
  279. 	if (x < 0.0f)
  280. 		return 0.0f;
  281. 	else if (x > 1.0f)
  282. 		return 1.0f;
  283. 	else
  284. 		return x;
  285. }
  286.  
  287. inline int toInt(float x) {
  288. 	return int(pow(clamp(x), 1 / 2.2) * 255 + .5);
  289. }
  290.  
  291. int main(int argc, char *argv[]) {
  292. 	clock_t start = clock();
  293.  
  294. 	Scene scene;
  295.  
  296. 	Sphere spheres[] = {
  297. 		Sphere(Vec3f(1e5f + 1.0f, 40.8f, 81.6f),	1e5f,	Vec3f::Zero(),	Vec3f(.75f, .25f, .25f),	{ Material::Type::DIFFUSE, 0.0f }),//Left
  298. 		Sphere(Vec3f(-1e5f + 99.0f, 40.8f, 81.6f),	1e5f,	Vec3f::Zero(),	Vec3f(.25f, .25f, .75f),	{ Material::Type::DIFFUSE, 0.0f }),//Rght
  299. 		Sphere(Vec3f(50.0f, 40.8f, 1e5f),			1e5f,	Vec3f::Zero(),	Vec3f(.75f, .75f, .75f),	{ Material::Type::DIFFUSE, 0.0f }),//Back
  300. 		Sphere(Vec3f(50.0f, 40.8f, -1e5f + 170.0f),	1e5f,	Vec3f::Zero(),	Vec3f::Zero(),				{ Material::Type::DIFFUSE, 0.0f }),//Frnt
  301. 		Sphere(Vec3f(50.0f, 1e5f, 81.6f),			1e5f,	Vec3f::Zero(),	Vec3f(.75f, .75f, .75f),	{ Material::Type::DIFFUSE, 0.0f }),//Botm
  302. 		Sphere(Vec3f(50.0f, -1e5f + 81.6f, 81.6f),	1e5f,	Vec3f::Zero(),	Vec3f(.75f, .75f, .75f),	{ Material::Type::DIFFUSE, 0.0f }),//Top
  303. 		Sphere(Vec3f(27.0f, 16.5f, 47.0f),			16.5f,	Vec3f::Zero(),	Vec3f(1.0f, 1.0f, 1.0f) * .999f,	{ Material::Type::SPECULAR, 0.0f }),//Mirr
  304. 		Sphere(Vec3f(73.0f, 16.5f, 78.0f),			16.5f,	Vec3f::Zero(),	Vec3f(1.0f, 1.0f, 1.0f) * .999f,	{ Material::Type::REFRACTION, 2.6f }),//Glas
  305. 		Sphere(Vec3f(50.0f, 681.6f - 0.27f, 81.6f),	600.0f,	Vec3f(12.0f, 12.0f, 12.0f), Vec3f::Zero(),	{ Material::Type::DIFFUSE, 0.0f }) //Lite
  306. 	};
  307.  
  308. 	for (auto sph : spheres)
  309. 		scene.add(sph);
  310.  
  311. 	const int w = 256;
  312. 	const int h = 256;
  313.  
  314. 	const int samps = argc == 2 ? atoi(argv[1]) / 4 : 100; // # samples
  315.  
  316. 	const Ray cam(Vec3f(50.0f, 52.0f, 295.6f), Vec3f(0.0f, -0.042612f, -1.0f)); // cam pos, dir
  317. 	const Vec3f cx(w * .5135f / h, 0.0f, 0.0f);
  318. 	const Vec3f cy = cx.cross(cam.direction).normalized() * .5135f;
  319.  
  320. 	Vec3f * buffer = new Vec3f[w * h];
  321.  
  322. //#pragma omp parallel for schedule(dynamic, 1)       // OpenMP
  323. 	// Loop over image rows
  324. 	for (int y = 0; y < h; y++) 
  325. 	{
  326. 		fprintf(stderr, "\rRendering (%d spp) %5.2f%%", samps * 4, 100.*y / (h - 1));
  327. 		RandomLCG rand(y);
  328.  
  329. 		// Loop cols
  330. 		for (int x = 0; x < w; x++) 
  331. 		{
  332. 			int const i = (h - y - 1) * w + x;
  333. 			//int const i = y * w + x;
  334. 			buffer[i] = Vec3f::Zero();
  335.  
  336. 			// 2x2 subpixel 
  337. 			for (int sy = 0; sy < 2; sy++) for (int sx = 0; sx < 2; sx++)
  338. 			{
  339. 				Vec3f r = Vec3f::Zero();
  340.  
  341. 				for (int s = 0; s < samps; s++) 
  342. 				{
  343. 					float r1 = 2.0f * rand();
  344. 					float r2 = 2.0f * rand();
  345. 					float dx = r1 < 1.0f ? sqrt(r1) - 1.0f : 1.0f - sqrt(2.0f - r1);
  346. 					float dy = r2 < 1.0f ? sqrt(r2) - 1.0f : 1.0f - sqrt(2.0f - r2);
  347.  
  348. 					Vec3f d = cx * (((sx + 0.5f + dx) / 2.0f + x) / w - 0.5f) +
  349. 						cy * (((sy + 0.5f + dy) / 2.0f + y) / h - 0.5f) + cam.direction;
  350.  
  351. 					r = r + radiance(scene, Ray(cam.origin + d * 140.0f, d), 0, 5, 30, rand) * (1.0f / samps);
  352. 				}
  353. 				buffer[i] = buffer[i] + Vec3f(clamp(r.x()), clamp(r.y()), clamp(r.z()));
  354. 			}
  355. 			buffer[i] = buffer[i] * 0.25f;
  356.  
  357. 			if (x % 16 == 0 && y % 16 == 0)
  358. 				std::cout << buffer[i] << std::endl;
  359.  
  360. 		}
  361. 	}
  362.  
  363. 	printf("\n%f sec\n", (float)(clock() - start) / CLOCKS_PER_SEC);
  364.  
  365. 	FILE *f = fopen("myimage4.ppm", "w"); // Write image to PPM file.
  366. 	fprintf(f, "P3\n%d %d\n%d\n", w, h, 255);
  367. 	for (int i = 0; i < w * h; i++)
  368. 		fprintf(f, "%d %d %d\n", toInt(buffer[i][0]), toInt(buffer[i][1]), toInt(buffer[i][2]));
  369. 	fclose(f);
  370.  
  371. 	return 0;
  372. }
  373.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin
copy 
Standard input is empty
compilation info 
prog.cpp:1:23: fatal error: Eigen\Dense: No such file or directory
compilation terminated.
stdout
copy 
Standard output is empty
https://ideone.com/C3t5Db
language:
C++14 (gcc 8.3)
created:
9 years ago
visibility:
secret

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!