GCC Code Coverage Report

Line	Branch	Exec	Source
1			/*
2			* Copyright (c) 2000-2022 Inria
3			* All rights reserved.
4			*
5			* Redistribution and use in source and binary forms, with or without
6			* modification, are permitted provided that the following conditions are met:
7			*
8			* * Redistributions of source code must retain the above copyright notice,
9			* this list of conditions and the following disclaimer.
10			* * Redistributions in binary form must reproduce the above copyright notice,
11			* this list of conditions and the following disclaimer in the documentation
12			* and/or other materials provided with the distribution.
13			* * Neither the name of the ALICE Project-Team nor the names of its
14			* contributors may be used to endorse or promote products derived from this
15			* software without specific prior written permission.
16			*
17			* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
18			* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
19			* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
20			* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
21			* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22			* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23			* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24			* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25			* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26			* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27			* POSSIBILITY OF SUCH DAMAGE.
28			*
29			* Contact: Bruno Levy
30			*
31			* https://www.inria.fr/fr/bruno-levy
32			*
33			* Inria,
34			* Domaine de Voluceau,
35			* 78150 Le Chesnay - Rocquencourt
36			* FRANCE
37			*
38			*/
39
40			#ifndef RAYTRACING_H
41			#define RAYTRACING_H
42
43
44
45			#include <geogram/basic/geometry.h>
46			#include <geogram/mesh/mesh.h>
47			#include <geogram/mesh/mesh_io.h>
48			#include <geogram/mesh/mesh_geometry.h>
49			#include <geogram/mesh/mesh_AABB.h>
50
51			#ifdef RAYTRACE_GUI
52			#include <geogram_gfx/third_party/imgui/imgui.h>
53			#include <geogram_gfx/imgui_ext/imgui_ext.h>
54			#endif
55
56			// We have everything in the .h (a bit ugly but it is used
57			// by two demos, so it is simpler like that), this makes
58			// Clang complain about vtbls generation (make it ignore
59			// the warning).
60			#ifdef __clang__
61			#pragma GCC diagnostic ignored "-Wweak-vtables"
62			#endif
63
64			/**
65			* \file raytracing.h
66			* \brief Implementation of a simple raytracer, to demonstrate
67			* AABB usage. See MeshObject below.
68			*/
69
70			namespace GEO {
71
72			/*******************************************************************/
73
74			/**
75			* \brief Normalizes the coordinates of a mesh
76			* in the unit box.
77			* \details A uniform scaling is applied.
78			* \param[in,out] M the mesh to be normalized.
79			*/
80		✗	inline void normalize_mesh(Mesh& M) {
81			double mesh_xyz_min[3];
82			double mesh_xyz_max[3];
83		✗	get_bbox(M, mesh_xyz_min, mesh_xyz_max);
84			double dim[3];
85		✗	dim[0] = (mesh_xyz_max[0] - mesh_xyz_min[0]);
86		✗	dim[1] = (mesh_xyz_max[1] - mesh_xyz_min[1]);
87		✗	dim[2] = (mesh_xyz_max[2] - mesh_xyz_min[2]);
88			double max_dim = std::max(dim[0], std::max(dim[1], dim[2]));
89		✗	double s = 1.0 / max_dim;
90			double T[3];
91		✗	T[0] = 0.5 * (max_dim - dim[0]);
92		✗	T[1] = 0.5 * (max_dim - dim[1]);
93		✗	T[2] = 0.5 * (max_dim - dim[2]);
94		✗	for(index_t v=0; v<M.vertices.nb(); ++v) {
95			double* p = M.vertices.point_ptr(v);
96		✗	double x = s * (T[0] + p[0] - mesh_xyz_min[0]);
97		✗	double y = s * (T[1] + p[1] - mesh_xyz_min[1]);
98		✗	double z = s * (T[2] + p[2] - mesh_xyz_min[2]);
99		✗	p[0] = x;
100		✗	p[1] = z;
101		✗	p[2] = 1.0-y;
102			}
103		✗	}
104
105			/*******************************************************************/
106
107			class Object;
108
109			/**
110			* \brief The small constant to shift a little bit ray intersections
111			* in order to avoid false positives when detecting shadows.
112			*/
113			const double epsilon_t = 1e-6;
114
115			/**
116			* \brief Multiplies two vec3 componentwise.
117			* \param[in] U , V the two vec3 to be multiplied
118			* \return the component-wise product
119			*/
120			inline vec3 mul(const vec3& U, const vec3& V) {
121			return vec3(
122		✗	U.x*V.x,
123		✗	U.y*V.y,
124		✗	U.z*V.z
125			);
126			}
127
128			/*******************************************************************/
129
130			/**
131			* \brief The Camere.
132			* \details Launches rays and stores the resulting image.
133			*/
134		✗	class Camera {
135			public:
136
137			/**
138			* \brief Camera constructor.
139			* \param[in] position the position of the camera.
140			* \param[in] target a point that is looked at.
141			* \param[in] image_width , image_height dimension of the image.
142			* \param[in] zoom viewing angle of the camera, in degrees.
143			*/
144		✗	Camera(
145			const vec3& position,
146			const vec3& target,
147			index_t image_width,
148			index_t image_height,
149			double zoom = 40.0
150		✗	) :
151		✗	image_width_(image_width),
152		✗	image_height_(image_height),
153		✗	image_(image_width*image_height*3),
154		✗	bpp_(3)
155			{
156		✗	update(position, target, zoom);
157		✗	}
158
159			/**
160			* \brief Camera constructor.
161			* \details Viewing parameters are not initialized.
162			* \param[in] image_width , image_height dimension of the image.
163			* \param[in] bpp bytes per pixel, 3 or 4
164			*/
165		✗	Camera(
166			index_t image_width, index_t image_height, index_t bpp
167		✗	) : image_width_(image_width),
168		✗	image_height_(image_height),
169		✗	image_(image_width*image_height*bpp),
170		✗	bpp_(bpp)
171			{
172		✗	geo_assert(bpp == 3 || bpp == 4);
173		✗	}
174
175			/**
176			* \brief Resizes the image.
177			* \param[in] new_width , new_height new image size, in pixels.
178			*/
179			void resize(index_t new_width, index_t new_height) {
180		✗	image_width_ = new_width;
181		✗	image_height_ = new_height;
182		✗	image_.resize(image_width_*image_height_*bpp_);
183			}
184
185			/**
186			* \brief Updates the camera parameters.
187			* \param[in] position the position of the camera.
188			* \param[in] target a point that is looked at.
189			* \param[in] zoom viewing angle of the camera, in degrees.
190			*/
191		✗	void update(
192			const vec3& position,
193			const vec3& target,
194			double zoom = 20.0
195			) {
196		✗	position_ = position;
197		✗	target_ = target;
198
199			// Viewing vector
200		✗	Z_ = normalize(target_ - position_);
201
202			// Horizontal direction
203			// We construct it as a vector both orthogonal
204			// to Z_ and to the vertical direction (0 0 1).
205		✗	X_ = cross(Z_, vec3(0.0, 0.0, 1.0));
206
207			// Vertical direction
208		✗	Y_ = cross(Z_,X_);
209
210			// Coordinate of the viewing plane along viewing vector
211		✗	double zp =
212		✗	(double(image_height_) / 2.0) / tan(zoom * M_PI / 180.0);
213
214			// Center of the viewing plane
215		✗	center_ = position_+zp*Z_;
216		✗	}
217
218			/**
219			* \brief gets the image width.
220			* \return the image width, in pixels.
221			*/
222			index_t image_width() const {
223		✗	return image_width_;
224			}
225
226			/**
227			* \brief gets the image height.
228			* \return the image height, in pixels.
229			*/
230			index_t image_height() const {
231		✗	return image_height_;
232			}
233
234			/**
235			* \brief gets the image data.
236			* \return a raw pointer to the image data.
237			*/
238			const Memory::byte* image_data() const {
239			return image_.data();
240			}
241
242			/**
243			* \brief Launches a primary ray.
244			* \param[in] X , Y the pixel coordinates of the ray.
245			* \return the launched ray.
246			*/
247		✗	Ray launch_ray(index_t X, index_t Y) const {
248			vec3 pixel3d = center_ +
249		✗	(double(X) - double(image_width_/2 ))*X_ +
250		✗	(double(Y) - double(image_height_/2))*Y_ ;
251
252			return Ray(
253			position_,
254			pixel3d - position_
255		✗	);
256			}
257
258			/**
259			* \brief Sets a pixel of the image.
260			* \param[in] X , Y the pixel coordinates of the ray.
261			* \param[in] color the RGB color of the pixel, with components
262			* in [0.0,1.0]. Components larger than 1.0 are clamped to 1.0.
263			*/
264		✗	void set_pixel(index_t X, index_t Y, const vec3& color) {
265			geo_debug_assert(X < image_width_);
266			geo_debug_assert(Y < image_height_);
267		✗	Memory::byte* pixel_base = &image_[(Y*image_width_+X)*bpp_];
268		✗	pixel_base[0] = Memory::byte(std::min(color.x, 1.0)*255.0);
269		✗	pixel_base[1] = Memory::byte(std::min(color.y, 1.0)*255.0);
270		✗	pixel_base[2] = Memory::byte(std::min(color.z, 1.0)*255.0);
271		✗	if(bpp_ == 4) {
272		✗	pixel_base[3] = 255;
273			}
274		✗	}
275
276			/**
277			* \brief Saves the image in PPM file format.
278			* \param[in] filename the name of the file where to save the image.
279			*/
280		✗	void save_image(const std::string& filename) const {
281		✗	geo_assert(bpp_ == 3);
282		✗	FILE* f = fopen(filename.c_str(),"wb");
283		✗	if(f == nullptr) {
284			std::cerr << "Could not create file: " << filename << std::endl;
285		✗	return;
286			}
287		✗	fprintf(
288		✗	f,"P6 %d %d %d ", int(image_width_), int(image_height_), 255
289			);
290		✗	fwrite(image_.data(), 1, image_.size(), f);
291		✗	fclose(f);
292			}
293
294			private:
295			vec3 position_;
296			vec3 target_;
297			vec3 center_;
298			vec3 X_;
299			vec3 Y_;
300			vec3 Z_;
301			index_t image_width_;
302			index_t image_height_;
303			vector<Memory::byte> image_;
304			index_t bpp_;
305			};
306
307			/*******************************************************************/
308
309			/**
310			* \brief A material.
311			*/
312			struct Material {
313			vec3 Kd; /**< Diffuse */
314			vec3 Kr; /**< Reflection */
315			vec3 Ke; /**< Emmission */
316		✗	Material():
317			Kd(0.7, 0.7, 0.7),
318			Kr(0.0, 0.0, 0.0),
319			Ke(0.0, 0.0, 0.0) {
320			}
321			/**
322			* \brief Tests whether this material is reflective.
323			* \retval true if this material has non-zero Kr, false otherwise.
324			*/
325			bool reflective() const {
326		✗	return (Kr.x != 0.0 || Kr.y != 0.0 || Kr.z != 0.0);
327			}
328			/**
329			* \brief Tests whether this material is emissive.
330			* \retval true if this material has non-zero Ke, false otherwise.
331			*/
332			bool emissive() const {
333		✗	return (Ke.x != 0.0 || Ke.y != 0.0 || Ke.z != 0.0);
334			}
335			};
336
337			/*******************************************************************/
338
339			/**
340			* \brief A Ray-Object intersection.
341			*/
342			struct Intersection {
343			/**
344			* \brief Intersection default constructor.
345			*/
346			Intersection() :
347		✗	t(Numeric::max_float64()),
348		✗	object(nullptr),
349			K(0.1, 0.1, 0.1) {
350			}
351			vec3 position; /**< position of the intersection. */
352			vec3 normal; /**< normal to the object. */
353			double t; /**< ray parameter of the intersection. */
354			const Object* object; /**< intersected object. */
355			vec3 K; /**< current computed ray color. */
356			Material material; /**< current material. */
357			};
358
359			/*******************************************************************/
360
361			/**
362			* \brief An object that can be raytraced.
363			*/
364			class Object {
365			public:
366
367			/**
368			* \brief Object destructor.
369			*/
370		✗	virtual ~Object() {}
371
372			/**
373			* \brief Computes the intersection with a ray.
374			* \details If there is an intersection and if it is nearer
375			* than the previous one, then replace it.
376			* \param[in] R the ray
377			* \param[in,out] I the nearest intersection along the ray.
378			*/
379			virtual void get_nearest_intersection(
380			const Ray& R, Intersection& I
381			) const = 0;
382
383			/**
384			* \brief Tests whether this object shadows a ray.
385			* \details This object shadows the ray R if there is an
386			* intersection between R.origin and R.origin + R.direction.
387			* Intersections further away than R.origin + R.direction
388			* are ignored.
389			* \param[in] R the ray. R.origin corresponds to a point
390			* queried for shadow. R.origin + R.direction corresponds to
391			* the light-source.
392			*/
393			virtual bool in_shadow(const Ray& R) const = 0;
394
395			/**
396			* \brief Sets the diffuse coefficient.
397			* \param[in] K the diffuse coefficient.
398			* \return a pointer to the Object, to allow chaining operations.
399			*/
400			Object* set_diffuse_coefficient(const vec3& K) {
401		✗	material_.Kd = K;
402			return this;
403			}
404
405			/**
406			* \brief Sets the reflection coefficient.
407			* \param[in] K the reflection coefficient.
408			* \return a pointer to the Object, to allow chaining operations.
409			*/
410			Object* set_reflection_coefficient(const vec3& K) {
411		✗	material_.Kr = K;
412			return this;
413			}
414
415			/**
416			* \brief Sets the emission coefficient.
417			* \param[in] K the emission coefficient.
418			* \return a pointer to the Object, to allow chaining operations.
419			*/
420			Object* set_emission_coefficient(const vec3& K) {
421			material_.Ke = K;
422			return this;
423			}
424
425			/**
426			* \brief Gets the Material.
427			* \return a const reference to the Material.
428			*/
429			const Material& material() const {
430			return material_;
431			}
432
433			/**
434			* \brief Gets the Material.
435			* \return a modifiable reference to the Material.
436			*/
437			Material& material() {
438			return material_;
439			}
440
441			Object* rename(const std::string& name) {
442		✗	name_ = name;
443		✗	return this;
444			}
445
446			const std::string& name() const {
447		✗	return name_;
448			}
449
450			#ifdef RAYTRACE_GUI
451
452		✗	bool edit_color(const std::string& name, vec3& K) {
453			float Kf[3];
454		✗	Kf[0] = float(K.x);
455		✗	Kf[1] = float(K.y);
456		✗	Kf[2] = float(K.z);
457		✗	bool result = ImGui::ColorEdit3WithPalette(name.c_str(), Kf);
458		✗	if(result) {
459		✗	K.x = double(Kf[0]);
460		✗	K.y = double(Kf[1]);
461		✗	K.z = double(Kf[2]);
462			}
463		✗	return result;
464			}
465
466		✗	bool edit_vector(const std::string& name, vec3& V) {
467			float Vf[3];
468		✗	Vf[0] = float(V.x);
469		✗	Vf[1] = float(V.y);
470		✗	Vf[2] = float(V.z);
471		✗	ImGui::SetNextItemWidth(-ImGui::CalcTextSize(name.c_str()).x);
472		✗	bool result = ImGui::DragFloat3(
473			name.c_str(), Vf, 0.1f, 0.0f, 0.0f, "%.3f"
474			);
475		✗	if(result) {
476		✗	V.x = double(Vf[0]);
477		✗	V.y = double(Vf[1]);
478		✗	V.z = double(Vf[2]);
479			}
480		✗	return result;
481			}
482
483		✗	bool edit_scalar(const std::string& name, double& V) {
484		✗	double zero = 0.0;
485		✗	ImGui::SetNextItemWidth(-ImGui::CalcTextSize(name.c_str()).x);
486		✗	return ImGui::DragScalar(
487			name.c_str(),
488			ImGuiDataType_Double,
489			&V,
490			0.005f,
491			&zero,
492			nullptr,
493			"%.3f"
494		✗	);
495			}
496
497			/**
498			* \brief Draws and handle the GUI.
499			* \retval true if an element was changed.
500			* \retval false otherwise.
501			*/
502		✗	virtual bool draw_gui() {
503		✗	ImGui::PushID(this);
504			bool result = false;
505		✗	ImGui::Separator();
506		✗	if(ImGui::Button("X")) {
507		✗	to_delete_ = this;
508			}
509		✗	ImGui::SameLine();
510		✗	ImGui::Text("%s", name().c_str());
511		✗	if(edit_color("Diffuse", material().Kd)) {
512			result = true;
513			}
514		✗	if(edit_color("Reflect", material().Kr)) {
515			result = true;
516			}
517		✗	ImGui::PopID();
518		✗	return result;
519			}
520			#endif
521
522			protected:
523			std::string name_;
524			Material material_;
525			static Object* to_delete_;
526			};
527
528			Object* Object::to_delete_ = nullptr;
529
530			/*******************************************************************/
531
532			/**
533			* \brief A sphere object.
534			*/
535			class Sphere : public Object {
536			public:
537
538			/**
539			* \brief Sphere constructor.
540			* \param[in] center the center of the sphere.
541			* \param[in] radius the radius of the sphere.
542			*/
543		✗	Sphere(const vec3& center, double radius) :
544		✗	center_(center), radius_(radius) {
545			}
546
547			/**
548			* \brief Gets the center.
549			* \return the center of the sphere.
550			*/
551			const vec3& center() const {
552			return center_;
553			}
554
555			/**
556			* \brief Gets the radius.
557			* \return the radius of the sphere.
558			*/
559			double radius() const {
560			return radius_;
561			}
562
563			/**
564			* \copydoc Object::get_nearest_intersection()
565			*/
566		✗	void get_nearest_intersection(
567			const Ray& R, Intersection& I
568			) const override {
569		✗	double t = get_intersection_t(R);
570		✗	if(t > epsilon_t && t < I.t) {
571		✗	I.t = t;
572		✗	I.object = this;
573		✗	I.material = material_;
574		✗	I.position = R.origin + t * R.direction;
575		✗	I.normal = normalize(I.position - center_);
576			}
577		✗	}
578
579			/**
580			* \copydoc Object::in_shadow()
581			*/
582		✗	bool in_shadow(const Ray& R) const override {
583		✗	double t = get_intersection_t(R);
584		✗	return (t > 0.0 && t < 1.0);
585			}
586
587			#ifdef RAYTRACE_GUI
588			/**
589			* \copydoc Object::draw_gui()
590			*/
591		✗	bool draw_gui() override {
592		✗	ImGui::PushID(this);
593		✗	bool result = Object::draw_gui();
594		✗	if(edit_vector("C", center_)) {
595			result = true;
596			}
597		✗	if(edit_scalar("R", radius_)) {
598			result = true;
599			}
600		✗	ImGui::PopID();
601		✗	return result;
602			}
603			#endif
604
605
606			protected:
607
608			/**
609			* \brief Gets the coordinate of the intersection between
610			* this sphere and a ray.
611			* \return the coordinate of the intersection along \p R
612			* or a negative number if there is no intersection.
613			*/
614		✗	double get_intersection_t(const Ray& R) const {
615			// Detail of the computation:
616			// M = O + tD (1) (parametric ray eqn)
617			// (M-C)^2 = M^2 - 2M.C + C^2 = R^2 (2) (implicit sphere eqn)
618			//(O + tD)^2 - 2(O+tD).C + C^2 = R^2 (inject (1) into (2)
619			// O^2 + 2tO.D + t^2D^2 -2O.C -2tD.C + C^2 = R^2
620			// t^2 (D^2) + 2t D.(O-C) + O^2 - 2 O.C + C^2 - R^2 = 0
621			// t^2 (D^2) + 2t D.(O-C) + (O-C)^2 - R^2 = 0
622			// This is a quadratic equation in t (a t^2 + b t + c = 0),
623			// let us solve it for t now !
624
625			double t = -1.0;
626
627			vec3 CO = R.origin - center_;
628			double a = length2(R.direction);
629		✗	double b = 2.0*dot(R.direction,CO);
630		✗	double c = length2(CO) - radius_*radius_;
631		✗	double delta = b*b - 4.0 * a * c;
632
633		✗	if(delta < 0.0) {
634			return -1.0;
635			}
636		✗	double sqrt_delta = sqrt(delta);
637		✗	t = (-b-sqrt_delta) / (2.0 * a);
638		✗	if(t > 0) {
639			return t;
640			}
641		✗	t = (-b+sqrt_delta) / (2.0 * a);
642		✗	return t;
643			}
644
645
646			protected:
647			vec3 center_;
648			double radius_;
649			};
650
651			/*******************************************************************/
652
653			/**
654			* \brief Light object.
655			* \details A Light appears as a colored sphere. The radius does not play
656			* a role in the lighting, this is just a point light.
657			*/
658			class Light : public Sphere {
659			public:
660			/**
661			* \brief Light constructor.
662			* \param[in] center the position of the light.
663			* \param[in] R the radius.
664			* \param[in] K the color of the light.
665			*/
666		✗	Light(const vec3& center, double R, const vec3& K) : Sphere(center, R) {
667		✗	material_.Kd = vec3(0.0, 0.0, 0.0);
668		✗	material_.Kr = vec3(0.0, 0.0, 0.0);
669		✗	material_.Ke = K;
670		✗	on_ = true;
671			}
672
673			bool on() const {
674		✗	return on_;
675			}
676
677
678			/**
679			* \copydoc Object::get_nearest_intersection()
680			*/
681		✗	void get_nearest_intersection(
682			const Ray& R, Intersection& I
683			) const override {
684		✗	Sphere::get_nearest_intersection(R,I);
685		✗	if(!on()) {
686		✗	I.material.Ke = vec3(0.0, 0.0, 0.0);
687			}
688		✗	}
689
690			#ifdef RAYTRACE_GUI
691			/**
692			* \copydoc Object::draw_gui()
693			*/
694		✗	bool draw_gui() override {
695		✗	ImGui::PushID(this);
696			bool result = false;
697		✗	ImGui::Separator();
698		✗	if(ImGui::Button("X")) {
699		✗	to_delete_ = this;
700			}
701		✗	ImGui::SameLine();
702		✗	ImGui::Text("%s", name().c_str());
703		✗	if(ImGui::Checkbox("##On", &on_)) {
704			result = true;
705			}
706		✗	ImGui::SameLine();
707		✗	if(edit_color("Emit.", material().Ke)) {
708			result = true;
709			}
710		✗	if(edit_vector("C", center_)) {
711			result = true;
712			}
713		✗	if(edit_scalar("R", radius_)) {
714			result = true;
715			}
716		✗	ImGui::PopID();
717		✗	return result;
718			}
719			#endif
720
721			private:
722			bool on_;
723
724			};
725
726			/*******************************************************************/
727
728			/**
729			* \brief Mesh object.
730			* \details Optimized ray-mesh intersections computed using an axis-aligned
731			* bounding box tree (geogram's MeshAABB).
732			*/
733			class MeshObject : public Object {
734			public:
735			/**
736			* \brief MeshObject constructor.
737			* \param[in] filename the name of the file that contains the mesh.
738			* \param[in] normalize if set, the mesh is normalized in
739			* the unit box after loading.
740			*/
741		✗	MeshObject(const std::string& filename, bool normalize=true) {
742		✗	mesh_load(filename, mesh_);
743		✗	if(normalize) {
744		✗	normalize_mesh(mesh_);
745			}
746		✗	AABB_.initialize(mesh_);
747		✗	}
748
749			/**
750			* \copydoc Object::get_nearest_intersection()
751			*/
752		✗	void get_nearest_intersection(
753			const Ray& R, Intersection& I
754			) const override {
755			MeshFacetsAABB::Intersection cur_I;
756		✗	if(AABB_.ray_nearest_intersection(R, cur_I)) {
757		✗	if(cur_I.t > epsilon_t && cur_I.t < I.t) {
758		✗	I.t = cur_I.t;
759		✗	I.object = this;
760		✗	I.material = material_;
761		✗	I.position = cur_I.p;
762		✗	I.normal = normalize(cur_I.N);
763			}
764			}
765		✗	}
766
767			/**
768			* \copydoc Object::in_shadow()
769			*/
770		✗	bool in_shadow(const Ray& R) const override {
771			vec3 p2 = R.origin + R.direction;
772		✗	return AABB_.segment_intersection(R.origin, p2);
773			}
774
775			private:
776			Mesh mesh_;
777			MeshFacetsAABB AABB_;
778			};
779
780			/*******************************************************************/
781
782			/**
783			* \brief The traditional checkerboard.
784			* \details Cannot avoid to have this in a raytracer
785			* (this is the tradition).
786			*/
787			class HorizontalCheckerboardPlane : public Object {
788			public:
789
790			/**
791			* \brief HorizontalCheckerboardPlane constructor;
792			* \param[in] z altitude of the plane.
793			*/
794		✗	HorizontalCheckerboardPlane(double z) : Z_(z) {
795			}
796
797			/**
798			* \copydoc Object::get_nearest_intersection()
799			*/
800		✗	void get_nearest_intersection(
801			const Ray& R, Intersection& I
802			) const override {
803		✗	if(R.direction.z != 0.0) {
804		✗	double t = (Z_ - R.origin.z) / R.direction.z;
805		✗	if(t > epsilon_t && t < I.t) {
806		✗	I.t = t;
807		✗	I.position = R.origin + t * R.direction;
808		✗	int X = int((I.position.x + 1000.0)* 2.0);
809		✗	int Y = int((I.position.y + 1000.0)* 2.0);
810		✗	double color = (((X&1) ^ (Y&1)) == 0) ? 0.0 : 1.0;
811		✗	I.object = this;
812		✗	I.material = material_;
813		✗	I.material.Kd.x *= color;
814		✗	I.material.Kd.y *= color;
815		✗	I.material.Kd.z *= color;
816		✗	I.normal = vec3(0.0, 0.0, 1.0);
817			}
818			}
819		✗	}
820
821			/**
822			* \copydoc Object::in_shadow()
823			*/
824		✗	bool in_shadow(const Ray& R) const override {
825		✗	if(R.direction.z == 0.0) {
826			return false;
827			}
828		✗	double t = (Z_ - R.origin.z) / R.direction.z;
829		✗	return (t >= 0.0 && t <= 1.0);
830			}
831
832			private:
833			double Z_;
834			};
835
836			/*******************************************************************/
837
838			/**
839			* \brief A scene to be ray-traced.
840			*/
841		✗	class Scene : public Object {
842			public:
843
844			/**
845			* \brief Scene destructor.
846			*/
847		✗	~Scene() override {
848		✗	for(index_t i=0; i<objects_.size(); ++i) {
849		✗	delete objects_[i];
850			}
851		✗	}
852
853			/**
854			* \brief Adds an object to the scene.
855			* \param[in] O a pointer to the object to be added.
856			* Pointer ownership is transfered to this Scene.
857			*/
858		✗	Object* add_object(Object* O) {
859		✗	objects_.push_back(O);
860		✗	Light* L = dynamic_cast<Light*>(O);
861		✗	if(L == nullptr) {
862		✗	real_objects_.push_back(O);
863			} else {
864		✗	lights_.push_back(L);
865			}
866		✗	if(O->name() == "") {
867		✗	O->rename("object " + String::to_string(objects_.size()));
868			}
869		✗	return O;
870			}
871
872
873			/**
874			* \brief Removes an object from the scene.
875			* \param[in] o the object to be removed.
876			* \details this does not deallocates the object.
877			*/
878		✗	void remove_object(Object* o) {
879		✗	for(index_t i=0; i<objects_.size(); ++i) {
880		✗	if(objects_[i] == o) {
881		✗	objects_.erase(objects_.begin() + std::ptrdiff_t(i));
882		✗	return;
883			}
884			}
885		✗	geo_assert_not_reached;
886			}
887
888			/**
889			* \brief Gets the number of objects in the scene.
890			* \return the number of objects, comprising lights.
891			*/
892			index_t nb_objects() const {
893			return objects_.size();
894			}
895
896			/**
897			* \brief Gets an object by index.
898			* \param[in] i the index, in [0..nb_objects()-1]
899			* \return a pointer to the ith object.
900			*/
901			Object* ith_object(index_t i) {
902		✗	return objects_[i];
903			}
904
905			/**
906			* \copydoc Object::get_nearest_intersection()
907			*/
908		✗	void get_nearest_intersection(
909			const Ray& R, Intersection& I
910			) const override {
911		✗	for(index_t i=0; i<objects_.size(); ++i) {
912		✗	objects_[i]->get_nearest_intersection(R,I);
913			}
914		✗	}
915
916			/**
917			* \copydoc Object::in_shadow()
918			*/
919		✗	bool in_shadow(const Ray& R) const override {
920		✗	for(index_t i=0; i<real_objects_.size(); ++i) {
921		✗	if(real_objects_[i]->in_shadow(R)) {
922			return true;
923			}
924			}
925			return false;
926			}
927
928			/**
929			* \brief Computes the lighting at an intersection.
930			* \details Launches shadow rays to the light sources.
931			* \param[in,out] I a reference to the intersection.
932			*/
933		✗	void compute_lighting(Intersection& I) const {
934			if(I.material.emissive()) {
935		✗	I.K = I.material.Ke;
936			} else {
937		✗	for(index_t i=0; i<lights_.size(); ++i) {
938		✗	if(!lights_[i]->on()) {
939			continue;
940			}
941			vec3 L = lights_[i]->center() - I.position;
942		✗	if(!in_shadow(Ray(I.position + epsilon_t*L, L))) {
943			double Lambert = dot(I.normal, L);
944		✗	if(Lambert > 0.0) {
945		✗	Lambert /= length(L);
946			I.K += Lambert*mul(
947		✗	I.material.Kd,lights_[i]->material().Ke
948			);
949			}
950			}
951			}
952			}
953		✗	}
954
955			/**
956			* \brief Launches a ray and computes the color.
957			* \details Reflected rays are recursively computed.
958			* \param[in] R the ray to be launched.
959			* \return the computed color.
960			*/
961		✗	vec3 raytrace(const Ray& R, index_t level=0) const {
962			Intersection I;
963		✗	get_nearest_intersection(R,I);
964		✗	if(I.object != nullptr) {
965		✗	compute_lighting(I);
966		✗	if(I.material.reflective() && level < 3) {
967		✗	vec3 D = R.direction;
968			Ray Reflected(
969			I.position,
970		✗	D - 2.0*dot(D, I.normal)*I.normal
971			);
972		✗	vec3 Kreflect = raytrace(Reflected,level+1);
973			I.K += mul(I.material.Kr, Kreflect);
974			}
975			}
976		✗	return I.K;
977			}
978
979			#ifdef RAYTRACE_GUI
980			/**
981			* \copydoc Object::draw_gui()
982			*/
983		✗	virtual bool draw_gui() override {
984			bool result = false;
985
986		✗	for(index_t i=0; i<objects_.size(); ++i) {
987		✗	if(to_delete_ == objects_[i]) {
988		✗	delete objects_[i];
989		✗	objects_.erase(objects_.begin() + std::ptrdiff_t(i));
990			result = true;
991		✗	break;
992			}
993			}
994
995		✗	for(index_t i=0; i<real_objects_.size(); ++i) {
996		✗	if(to_delete_ == real_objects_[i]) {
997		✗	real_objects_.erase(
998			real_objects_.begin() + std::ptrdiff_t(i)
999			);
1000		✗	break;
1001			}
1002			}
1003
1004		✗	for(index_t i=0; i<lights_.size(); ++i) {
1005		✗	if(to_delete_ == lights_[i]) {
1006		✗	lights_.erase(lights_.begin() + std::ptrdiff_t(i));
1007		✗	break;
1008			}
1009			}
1010
1011		✗	to_delete_ = nullptr;
1012
1013		✗	for(index_t i=0; i<objects_.size(); ++i) {
1014		✗	if(objects_[i]->draw_gui()) {
1015			result = true;
1016			}
1017			}
1018		✗	return result;
1019			}
1020			#endif
1021
1022			private:
1023			vector<Object*> objects_; /**< all the objects. */
1024			vector<Object*> real_objects_; /**< all the objects but the lights. */
1025			vector<Light*> lights_; /**< the lights are here. */
1026			};
1027
1028			/*******************************************************************/
1029			/** Utilities */
1030			/*******************************************************************/
1031
1032			/**
1033			* \brief Sets a 4x4 homogeneous transform matrix from a translation
1034			* and a quaternion.
1035			* \param[out] M the matrix.
1036			* \param[in] Tx , Ty , Tz the translation.
1037			* \param[in] Qx , Qy , Qz , Qw the quaternion.
1038			*/
1039			inline void set_mat4_from_translation_and_quaternion(
1040			mat4& M,
1041			double Tx, double Ty, double Tz,
1042			double Qx, double Qy, double Qz, double Qw
1043			) {
1044			// for unit q, just set s = 2 or set xs = Qx + Qx, etc.
1045			double s = 2.0 / (Qx*Qx + Qy*Qy + Qz*Qz + Qw*Qw);
1046
1047			double xs = Qx * s;
1048			double ys = Qy * s;
1049			double zs = Qz * s;
1050
1051			double wx = Qw * xs;
1052			double wy = Qw * ys;
1053			double wz = Qw * zs;
1054
1055			double xx = Qx * xs;
1056			double xy = Qx * ys;
1057			double xz = Qx * zs;
1058
1059			double yy = Qy * ys;
1060			double yz = Qy * zs;
1061			double zz = Qz * zs;
1062
1063			M(0,0) = 1.0 - (yy + zz);
1064			M(0,1) = xy - wz;
1065			M(0,2) = xz + wy;
1066			M(0,3) = 0.0;
1067
1068			M(1,0) = xy + wz;
1069			M(1,1) = 1 - (xx + zz);
1070			M(1,2) = yz - wx;
1071			M(1,3) = 0.0;
1072
1073			M(2,0) = xz - wy;
1074			M(2,1) = yz + wx;
1075			M(2,2) = 1 - (xx + yy);
1076			M(2,3) = 0.0;
1077
1078			M(3,0) = Tx;
1079			M(3,1) = Ty;
1080			M(3,2) = Tz;
1081			M(3,3) = 1.0;
1082			}
1083
1084		✗	inline vec3 random_color() {
1085			vec3 result;
1086		✗	while(length2(result) < 0.1) {
1087		✗	result = vec3(
1088		✗	Numeric::random_float64(),
1089		✗	Numeric::random_float64(),
1090		✗	Numeric::random_float64()
1091			);
1092			}
1093		✗	return result;
1094			}
1095			}
1096
1097
1098			#endif
1099