GCC Code Coverage Report

Line	Branch	Exec	Source
1			// Copyright Contributors to the OpenVDB Project
2			// SPDX-License-Identifier: MPL-2.0
3
4			/// @file RayTracer.h
5			///
6			/// @author Ken Museth
7			///
8			/// @brief Defines two simple but multithreaded renders, a level-set
9			/// ray tracer and a volume render. To support these renders we also define
10			/// perspective and orthographic cameras (both designed to mimic a Houdini camera),
11			/// a Film class and some rather naive shaders.
12			///
13			/// @note These classes are included mainly as reference implementations for
14			/// ray-tracing of OpenVDB volumes. In other words they are not intended for
15			/// production-quality rendering, but could be used for fast pre-visualization
16			/// or as a starting point for a more serious render.
17
18			#ifndef OPENVDB_TOOLS_RAYTRACER_HAS_BEEN_INCLUDED
19			#define OPENVDB_TOOLS_RAYTRACER_HAS_BEEN_INCLUDED
20
21			#include <openvdb/Types.h>
22			#include <openvdb/math/BBox.h>
23			#include <openvdb/math/Ray.h>
24			#include <openvdb/math/Math.h>
25			#include <openvdb/tools/RayIntersector.h>
26			#include <openvdb/tools/Interpolation.h>
27			#include <openvdb/openvdb.h>
28			#include <deque>
29			#include <iostream>
30			#include <fstream>
31			#include <limits>
32			#include <memory>
33			#include <string>
34			#include <type_traits>
35			#include <vector>
36
37			namespace openvdb {
38			OPENVDB_USE_VERSION_NAMESPACE
39			namespace OPENVDB_VERSION_NAME {
40			namespace tools {
41
42			// Forward declarations
43			class BaseCamera;
44			class BaseShader;
45
46			/// @brief Ray-trace a volume.
47			template<typename GridT>
48			void rayTrace(const GridT&,
49			const BaseShader&,
50			BaseCamera&,
51			size_t pixelSamples = 1,
52			unsigned int seed = 0,
53			bool threaded = true);
54
55			/// @brief Ray-trace a volume using a given ray intersector.
56			template<typename GridT, typename IntersectorT>
57			void rayTrace(const GridT&,
58			const IntersectorT&,
59			const BaseShader&,
60			BaseCamera&,
61			size_t pixelSamples = 1,
62			unsigned int seed = 0,
63			bool threaded = true);
64
65
66			///////////////////////////////LEVEL SET RAY TRACER ///////////////////////////////////////
67
68			/// @brief A (very) simple multithreaded ray tracer specifically for narrow-band level sets.
69			/// @details Included primarily as a reference implementation.
70			template<typename GridT, typename IntersectorT = tools::LevelSetRayIntersector<GridT> >
71			class LevelSetRayTracer
72			{
73			public:
74			using GridType = GridT;
75			using Vec3Type = typename IntersectorT::Vec3Type;
76			using RayType = typename IntersectorT::RayType;
77
78			/// @brief Constructor based on an instance of the grid to be rendered.
79			LevelSetRayTracer(const GridT& grid,
80			const BaseShader& shader,
81			BaseCamera& camera,
82			size_t pixelSamples = 1,
83			unsigned int seed = 0);
84
85			/// @brief Constructor based on an instance of the intersector
86			/// performing the ray-intersections.
87			LevelSetRayTracer(const IntersectorT& inter,
88			const BaseShader& shader,
89			BaseCamera& camera,
90			size_t pixelSamples = 1,
91			unsigned int seed = 0);
92
93			/// @brief Copy constructor
94			LevelSetRayTracer(const LevelSetRayTracer& other);
95
96			/// @brief Destructor
97			~LevelSetRayTracer();
98
99			/// @brief Set the level set grid to be ray-traced
100			void setGrid(const GridT& grid);
101
102			/// @brief Set the intersector that performs the actual
103			/// intersection of the rays against the narrow-band level set.
104			void setIntersector(const IntersectorT& inter);
105
106			/// @brief Set the shader derived from the abstract BaseShader class.
107			///
108			/// @note The shader is not assumed to be thread-safe so each
109			/// thread will get its only deep copy. For instance it could
110			/// contains a ValueAccessor into another grid with auxiliary
111			/// shading information. Thus, make sure it is relatively
112			/// light-weight and efficient to copy (which is the case for ValueAccesors).
113			void setShader(const BaseShader& shader);
114
115			/// @brief Set the camera derived from the abstract BaseCamera class.
116			void setCamera(BaseCamera& camera);
117
118			/// @brief Set the number of pixel samples and the seed for
119			/// jittered sub-rays. A value larger than one implies
120			/// anti-aliasing by jittered super-sampling.
121			/// @throw ValueError if pixelSamples is equal to zero.
122			void setPixelSamples(size_t pixelSamples, unsigned int seed = 0);
123
124			/// @brief Perform the actual (potentially multithreaded) ray-tracing.
125			void render(bool threaded = true) const;
126
127			/// @brief Public method required by tbb::parallel_for.
128			/// @warning Never call it directly.
129			void operator()(const tbb::blocked_range<size_t>& range) const;
130
131			private:
132			const bool mIsMaster;
133			double* mRand;
134			IntersectorT mInter;
135			std::unique_ptr<const BaseShader> mShader;
136			BaseCamera* mCamera;
137			size_t mSubPixels;
138			};// LevelSetRayTracer
139
140
141			///////////////////////////////VOLUME RENDER ///////////////////////////////////////
142
143			/// @brief A (very) simple multithreaded volume render specifically for scalar density.
144			/// @details Included primarily as a reference implementation.
145			/// @note It will only compile if the IntersectorT is templated on a Grid with a
146			/// floating-point voxel type.
147			template <typename IntersectorT, typename SamplerT = tools::BoxSampler>
148			class VolumeRender
149			{
150			public:
151
152			using GridType = typename IntersectorT::GridType;
153			using RayType = typename IntersectorT::RayType;
154			using ValueType = typename GridType::ValueType;
155			using AccessorType = typename GridType::ConstAccessor;
156			using SamplerType = tools::GridSampler<AccessorType, SamplerT>;
157			static_assert(std::is_floating_point<ValueType>::value,
158			"VolumeRender requires a floating-point-valued grid");
159
160			/// @brief Constructor taking an intersector and a base camera.
161			VolumeRender(const IntersectorT& inter, BaseCamera& camera);
162
163			/// @brief Copy constructor which creates a thread-safe clone
164			VolumeRender(const VolumeRender& other);
165
166			/// @brief Perform the actual (potentially multithreaded) volume rendering.
167			void render(bool threaded=true) const;
168
169			/// @brief Set the camera derived from the abstract BaseCamera class.
170		✗	void setCamera(BaseCamera& camera) { mCamera = &camera; }
171
172			/// @brief Set the intersector that performs the actual
173			/// intersection of the rays against the volume.
174			void setIntersector(const IntersectorT& inter);
175
176			/// @brief Set the vector components of a directional light source
177			/// @throw ArithmeticError if input is a null vector.
178		✗	void setLightDir(Real x, Real y, Real z) { mLightDir = Vec3R(x,y,z).unit(); }
179
180			/// @brief Set the color of the directional light source.
181		✗	void setLightColor(Real r, Real g, Real b) { mLightColor = Vec3R(r,g,b); }
182
183			/// @brief Set the integration step-size in voxel units for the primay ray.
184		✗	void setPrimaryStep(Real primaryStep) { mPrimaryStep = primaryStep; }
185
186			/// @brief Set the integration step-size in voxel units for the primay ray.
187		✗	void setShadowStep(Real shadowStep) { mShadowStep = shadowStep; }
188
189			/// @brief Set Scattering coefficients.
190		✗	void setScattering(Real x, Real y, Real z) { mScattering = Vec3R(x,y,z); }
191
192			/// @brief Set absorption coefficients.
193		✗	void setAbsorption(Real x, Real y, Real z) { mAbsorption = Vec3R(x,y,z); }
194
195			/// @brief Set parameter that imitates multi-scattering. A value
196			/// of zero implies no multi-scattering.
197		✗	void setLightGain(Real gain) { mLightGain = gain; }
198
199			/// @brief Set the cut-off value for density and transmittance.
200		✗	void setCutOff(Real cutOff) { mCutOff = cutOff; }
201
202			/// @brief Print parameters, statistics, memory usage and other information.
203			/// @param os a stream to which to write textual information
204			/// @param verboseLevel 1: print parameters only; 2: include grid
205			/// statistics; 3: include memory usage
206			void print(std::ostream& os = std::cout, int verboseLevel = 1);
207
208			/// @brief Public method required by tbb::parallel_for.
209			/// @warning Never call it directly.
210			void operator()(const tbb::blocked_range<size_t>& range) const;
211
212			private:
213
214			AccessorType mAccessor;
215			BaseCamera* mCamera;
216			std::unique_ptr<IntersectorT> mPrimary, mShadow;
217			Real mPrimaryStep, mShadowStep, mCutOff, mLightGain;
218			Vec3R mLightDir, mLightColor, mAbsorption, mScattering;
219			};//VolumeRender
220
221			//////////////////////////////////////// FILM ////////////////////////////////////////
222
223			/// @brief A simple class that allows for concurrent writes to pixels in an image,
224			/// background initialization of the image, and PPM file output.
225			class Film
226			{
227			public:
228			/// @brief Floating-point RGBA components in the range [0, 1].
229			/// @details This is our preferred representation for color processing.
230			struct RGBA
231			{
232			using ValueT = float;
233
234			RGBA() : r(0), g(0), b(0), a(1) {}
235			explicit RGBA(ValueT intensity) : r(intensity), g(intensity), b(intensity), a(1) {}
236		✗	RGBA(ValueT _r, ValueT _g, ValueT _b, ValueT _a = static_cast<ValueT>(1.0)):
237		✗	r(_r), g(_g), b(_b), a(_a)
238			{}
239			RGBA(double _r, double _g, double _b, double _a = 1.0)
240			: r(static_cast<ValueT>(_r))
241			, g(static_cast<ValueT>(_g))
242			, b(static_cast<ValueT>(_b))
243			, a(static_cast<ValueT>(_a))
244			{}
245
246		✗	RGBA operator* (ValueT scale) const { return RGBA(r*scale, g*scale, b*scale);}
247			RGBA operator+ (const RGBA& rhs) const { return RGBA(r+rhs.r, g+rhs.g, b+rhs.b);}
248			RGBA operator* (const RGBA& rhs) const { return RGBA(r*rhs.r, g*rhs.g, b*rhs.b);}
249		✗	RGBA& operator+=(const RGBA& rhs) { r+=rhs.r; g+=rhs.g; b+=rhs.b; a+=rhs.a; return *this;}
250
251			void over(const RGBA& rhs)
252			{
253			const float s = rhs.a*(1.0f-a);
254			r = a*r+s*rhs.r;
255			g = a*g+s*rhs.g;
256			b = a*b+s*rhs.b;
257			a = a + s;
258			}
259
260			ValueT r, g, b, a;
261			};
262
263
264			Film(size_t width, size_t height)
265			: mWidth(width), mHeight(height), mSize(width*height), mPixels(new RGBA[mSize])
266			{
267			}
268			Film(size_t width, size_t height, const RGBA& bg)
269			: mWidth(width), mHeight(height), mSize(width*height), mPixels(new RGBA[mSize])
270			{
271			this->fill(bg);
272			}
273
274			const RGBA& pixel(size_t w, size_t h) const
275			{
276			assert(w < mWidth);
277			assert(h < mHeight);
278			return mPixels[w + h*mWidth];
279			}
280
281		✗	RGBA& pixel(size_t w, size_t h)
282			{
283		✗	assert(w < mWidth);
284		✗	assert(h < mHeight);
285		✗	return mPixels[w + h*mWidth];
286			}
287
288			void fill(const RGBA& rgb=RGBA(0)) { for (size_t i=0; i<mSize; ++i) mPixels[i] = rgb; }
289			void checkerboard(const RGBA& c1=RGBA(0.3f), const RGBA& c2=RGBA(0.6f), size_t size=32)
290			{
291			RGBA *p = mPixels.get();
292			for (size_t j = 0; j < mHeight; ++j) {
293			for (size_t i = 0; i < mWidth; ++i, ++p) {
294			*p = ((i & size) ^ (j & size)) ? c1 : c2;
295			}
296			}
297			}
298
299			template <typename Type = unsigned char>
300			std::unique_ptr<Type[]> convertToBitBuffer(const bool alpha = true) const
301			{
302			const size_t totalSize = mSize * (alpha ? 4 : 3);
303			std::unique_ptr<Type[]> buffer(new Type[totalSize]);
304			Type *q = buffer.get();
305			const RGBA* p = this->pixels();
306			size_t n = mSize;
307			while (n--) {
308			*q++ = static_cast<Type>(255.0f*(*p).r);
309			*q++ = static_cast<Type>(255.0f*(*p).g);
310			*q++ = static_cast<Type>(255.0f*(*p).b);
311			if(alpha)
312			*q++ = static_cast<Type>(255.0f*(*p).a);
313			++p;
314			}
315			return buffer;
316			}
317
318			void savePPM(const std::string& fileName)
319			{
320			std::string name(fileName);
321			if (name.find_last_of(".") == std::string::npos) name.append(".ppm");
322
323			std::ofstream os(name.c_str(), std::ios_base::binary);
324			if (!os.is_open()) {
325			std::cerr << "Error opening PPM file \"" << name << "\"" << std::endl;
326			return;
327			}
328
329			auto buf = this->convertToBitBuffer<unsigned char>(/*alpha=*/false);
330			unsigned char* tmp = buf.get();
331
332			os << "P6\n" << mWidth << " " << mHeight << "\n255\n";
333			os.write(reinterpret_cast<const char*>(&(*tmp)), 3 * mSize * sizeof(unsigned char));
334			}
335
336		✗	size_t width() const { return mWidth; }
337		✗	size_t height() const { return mHeight; }
338			size_t numPixels() const { return mSize; }
339			const RGBA* pixels() const { return mPixels.get(); }
340
341			private:
342			size_t mWidth, mHeight, mSize;
343			std::unique_ptr<RGBA[]> mPixels;
344			};// Film
345
346
347			//////////////////////////////////////// CAMERAS ////////////////////////////////////////
348
349			/// Abstract base class for the perspective and orthographic cameras
350			class BaseCamera
351			{
352			public:
353			BaseCamera(Film& film, const Vec3R& rotation, const Vec3R& translation,
354			double frameWidth, double nearPlane, double farPlane)
355			: mFilm(&film)
356			, mScaleWidth(frameWidth)
357			, mScaleHeight(frameWidth * double(film.height()) / double(film.width()))
358			{
359			assert(nearPlane > 0 && farPlane > nearPlane);
360			mScreenToWorld.accumPostRotation(math::X_AXIS, rotation[0] * M_PI / 180.0);
361			mScreenToWorld.accumPostRotation(math::Y_AXIS, rotation[1] * M_PI / 180.0);
362			mScreenToWorld.accumPostRotation(math::Z_AXIS, rotation[2] * M_PI / 180.0);
363			mScreenToWorld.accumPostTranslation(translation);
364			this->initRay(nearPlane, farPlane);
365			}
366
367			virtual ~BaseCamera() {}
368
369		✗	Film::RGBA& pixel(size_t i, size_t j) { return mFilm->pixel(i, j); }
370
371		✗	size_t width() const { return mFilm->width(); }
372		✗	size_t height() const { return mFilm->height(); }
373
374			/// Rotate the camera so its negative z-axis points at xyz and its
375			/// y axis is in the plane of the xyz and up vectors. In other
376			/// words the camera will look at xyz and use up as the
377			/// horizontal direction.
378			void lookAt(const Vec3R& xyz, const Vec3R& up = Vec3R(0.0, 1.0, 0.0))
379			{
380			const Vec3R orig = mScreenToWorld.applyMap(Vec3R(0.0));
381			const Vec3R dir = orig - xyz;
382			try {
383			Mat4d xform = math::aim<Mat4d>(dir, up);
384			xform.postTranslate(orig);
385			mScreenToWorld = math::AffineMap(xform);
386			this->initRay(mRay.t0(), mRay.t1());
387			} catch (...) {}
388			}
389
390			Vec3R rasterToScreen(double i, double j, double z) const
391			{
392			return Vec3R( (2 * i / double(mFilm->width()) - 1) * mScaleWidth,
393			(1 - 2 * j / double(mFilm->height())) * mScaleHeight, z );
394			}
395
396			/// @brief Return a Ray in world space given the pixel indices and
397			/// optional offsets in the range [0, 1]. An offset of 0.5 corresponds
398			/// to the center of the pixel.
399			virtual math::Ray<double> getRay(
400			size_t i, size_t j, double iOffset = 0.5, double jOffset = 0.5) const = 0;
401
402			protected:
403			void initRay(double t0, double t1)
404			{
405			mRay.setTimes(t0, t1);
406			mRay.setEye(mScreenToWorld.applyMap(Vec3R(0.0)));
407			mRay.setDir(mScreenToWorld.applyJacobian(Vec3R(0.0, 0.0, -1.0)));
408			}
409
410			Film* mFilm;
411			double mScaleWidth, mScaleHeight;
412			math::Ray<double> mRay;
413			math::AffineMap mScreenToWorld;
414			};// BaseCamera
415
416
417			class PerspectiveCamera: public BaseCamera
418			{
419			public:
420			/// @brief Constructor
421			/// @param film film (i.e. image) defining the pixel resolution
422			/// @param rotation rotation in degrees of the camera in world space
423			/// (applied in x, y, z order)
424			/// @param translation translation of the camera in world-space units,
425			/// applied after rotation
426			/// @param focalLength focal length of the camera in mm
427			/// (the default of 50mm corresponds to Houdini's default camera)
428			/// @param aperture width in mm of the frame, i.e., the visible field
429			/// (the default 41.2136 mm corresponds to Houdini's default camera)
430			/// @param nearPlane depth of the near clipping plane in world-space units
431			/// @param farPlane depth of the far clipping plane in world-space units
432			///
433			/// @details If no rotation or translation is provided, the camera is placed
434			/// at (0,0,0) in world space and points in the direction of the negative z axis.
435			PerspectiveCamera(Film& film,
436			const Vec3R& rotation = Vec3R(0.0),
437			const Vec3R& translation = Vec3R(0.0),
438			double focalLength = 50.0,
439			double aperture = 41.2136,
440			double nearPlane = 1e-3,
441			double farPlane = std::numeric_limits<double>::max())
442			: BaseCamera(film, rotation, translation, 0.5*aperture/focalLength, nearPlane, farPlane)
443			{
444			}
445
446			~PerspectiveCamera() override = default;
447
448			/// @brief Return a Ray in world space given the pixel indices and
449			/// optional offsets in the range [0,1]. An offset of 0.5 corresponds
450			/// to the center of the pixel.
451			math::Ray<double> getRay(
452			size_t i, size_t j, double iOffset = 0.5, double jOffset = 0.5) const override
453			{
454			math::Ray<double> ray(mRay);
455			Vec3R dir = BaseCamera::rasterToScreen(Real(i) + iOffset, Real(j) + jOffset, -1.0);
456			dir = BaseCamera::mScreenToWorld.applyJacobian(dir);
457			dir.normalize();
458			ray.scaleTimes(1.0/dir.dot(ray.dir()));
459			ray.setDir(dir);
460			return ray;
461			}
462
463			/// @brief Return the horizontal field of view in degrees given a
464			/// focal lenth in mm and the specified aperture in mm.
465			static double focalLengthToFieldOfView(double length, double aperture)
466			{
467			return 360.0 / M_PI * atan(aperture/(2.0*length));
468			}
469			/// @brief Return the focal length in mm given a horizontal field of
470			/// view in degrees and the specified aperture in mm.
471			static double fieldOfViewToFocalLength(double fov, double aperture)
472			{
473			return aperture/(2.0*(tan(fov * M_PI / 360.0)));
474			}
475			};// PerspectiveCamera
476
477
478			class OrthographicCamera: public BaseCamera
479			{
480			public:
481			/// @brief Constructor
482			/// @param film film (i.e. image) defining the pixel resolution
483			/// @param rotation rotation in degrees of the camera in world space
484			/// (applied in x, y, z order)
485			/// @param translation translation of the camera in world-space units,
486			/// applied after rotation
487			/// @param frameWidth width in of the frame in world-space units
488			/// @param nearPlane depth of the near clipping plane in world-space units
489			/// @param farPlane depth of the far clipping plane in world-space units
490			///
491			/// @details If no rotation or translation is provided, the camera is placed
492			/// at (0,0,0) in world space and points in the direction of the negative z axis.
493			OrthographicCamera(Film& film,
494			const Vec3R& rotation = Vec3R(0.0),
495			const Vec3R& translation = Vec3R(0.0),
496			double frameWidth = 1.0,
497			double nearPlane = 1e-3,
498			double farPlane = std::numeric_limits<double>::max())
499			: BaseCamera(film, rotation, translation, 0.5*frameWidth, nearPlane, farPlane)
500			{
501			}
502			~OrthographicCamera() override = default;
503
504			math::Ray<double> getRay(
505			size_t i, size_t j, double iOffset = 0.5, double jOffset = 0.5) const override
506			{
507			math::Ray<double> ray(mRay);
508			Vec3R eye = BaseCamera::rasterToScreen(Real(i) + iOffset, Real(j) + jOffset, 0.0);
509			ray.setEye(BaseCamera::mScreenToWorld.applyMap(eye));
510			return ray;
511			}
512			};// OrthographicCamera
513
514
515			//////////////////////////////////////// SHADERS ////////////////////////////////////////
516
517
518			/// Abstract base class for the shaders
519			class BaseShader
520			{
521			public:
522			using RayT = math::Ray<Real>;
523			BaseShader() {}
524			BaseShader(const BaseShader&) = default;
525			virtual ~BaseShader() = default;
526			/// @brief Defines the interface of the virtual function that returns a RGB color.
527			/// @param xyz World position of the intersection point.
528			/// @param nml Normal in world space at the intersection point.
529			/// @param dir Direction of the ray in world space.
530			virtual Film::RGBA operator()(const Vec3R& xyz, const Vec3R& nml, const Vec3R& dir) const = 0;
531			virtual BaseShader* copy() const = 0;
532			};
533
534
535			/// @brief Shader that produces a simple matte.
536			///
537			/// @details The color can either be constant (if GridT =
538			/// Film::RGBA which is the default) or defined in a separate Vec3
539			/// color grid. Use SamplerType to define the order of interpolation
540			/// (default is zero order, i.e. closes-point).
541			template<typename GridT = Film::RGBA,
542			typename SamplerType = tools::PointSampler>
543			class MatteShader: public BaseShader
544			{
545			public:
546			MatteShader(const GridT& grid) : mAcc(grid.getAccessor()), mXform(&grid.transform()) {}
547			MatteShader(const MatteShader&) = default;
548			~MatteShader() override = default;
549			Film::RGBA operator()(const Vec3R& xyz, const Vec3R&, const Vec3R&) const override
550			{
551			typename GridT::ValueType v = zeroVal<typename GridT::ValueType>();
552			SamplerType::sample(mAcc, mXform->worldToIndex(xyz), v);
553			return Film::RGBA(v[0], v[1], v[2]);
554			}
555			BaseShader* copy() const override { return new MatteShader<GridT, SamplerType>(*this); }
556
557			private:
558			typename GridT::ConstAccessor mAcc;
559			const math::Transform* mXform;
560			};
561
562			// Template specialization using a constant color of the material.
563			template<typename SamplerType>
564			class MatteShader<Film::RGBA, SamplerType>: public BaseShader
565			{
566			public:
567			MatteShader(const Film::RGBA& c = Film::RGBA(1.0f)): mRGBA(c) {}
568			MatteShader(const MatteShader&) = default;
569			~MatteShader() override = default;
570			Film::RGBA operator()(const Vec3R&, const Vec3R&, const Vec3R&) const override
571			{
572			return mRGBA;
573			}
574			BaseShader* copy() const override { return new MatteShader<Film::RGBA, SamplerType>(*this); }
575
576			private:
577			const Film::RGBA mRGBA;
578			};
579
580
581			/// @brief Color shader that treats the surface normal (x, y, z) as an
582			/// RGB color.
583			///
584			/// @details The color can either be constant (if GridT =
585			/// Film::RGBA which is the default) or defined in a separate Vec3
586			/// color grid. Use SamplerType to define the order of interpolation
587			/// (default is zero order, i.e. closes-point).
588			template<typename GridT = Film::RGBA,
589			typename SamplerType = tools::PointSampler>
590			class NormalShader: public BaseShader
591			{
592			public:
593			NormalShader(const GridT& grid) : mAcc(grid.getAccessor()), mXform(&grid.transform()) {}
594			NormalShader(const NormalShader&) = default;
595			~NormalShader() override = default;
596			Film::RGBA operator()(const Vec3R& xyz, const Vec3R& normal, const Vec3R&) const override
597			{
598			typename GridT::ValueType v = zeroVal<typename GridT::ValueType>();
599			SamplerType::sample(mAcc, mXform->worldToIndex(xyz), v);
600			return Film::RGBA(v[0]*(normal[0]+1.0), v[1]*(normal[1]+1.0), v[2]*(normal[2]+1.0));
601			}
602			BaseShader* copy() const override { return new NormalShader<GridT, SamplerType>(*this); }
603
604			private:
605			typename GridT::ConstAccessor mAcc;
606			const math::Transform* mXform;
607			};
608
609			// Template specialization using a constant color of the material.
610			template<typename SamplerType>
611			class NormalShader<Film::RGBA, SamplerType>: public BaseShader
612			{
613			public:
614			NormalShader(const Film::RGBA& c = Film::RGBA(1.0f)) : mRGBA(c*0.5f) {}
615			NormalShader(const NormalShader&) = default;
616			~NormalShader() override = default;
617			Film::RGBA operator()(const Vec3R&, const Vec3R& normal, const Vec3R&) const override
618			{
619			return mRGBA * Film::RGBA(normal[0] + 1.0, normal[1] + 1.0, normal[2] + 1.0);
620			}
621			BaseShader* copy() const override { return new NormalShader<Film::RGBA, SamplerType>(*this); }
622
623			private:
624			const Film::RGBA mRGBA;
625			};
626
627
628			/// @brief Color shader that treats position (x, y, z) as an RGB color in a
629			/// cube defined from an axis-aligned bounding box in world space.
630			///
631			/// @details The color can either be constant (if GridT =
632			/// Film::RGBA which is the default) or defined in a separate Vec3
633			/// color grid. Use SamplerType to define the order of interpolation
634			/// (default is zero order, i.e. closes-point).
635			template<typename GridT = Film::RGBA,
636			typename SamplerType = tools::PointSampler>
637			class PositionShader: public BaseShader
638			{
639			public:
640			PositionShader(const math::BBox<Vec3R>& bbox, const GridT& grid)
641			: mMin(bbox.min())
642			, mInvDim(1.0/bbox.extents())
643			, mAcc(grid.getAccessor())
644			, mXform(&grid.transform())
645			{
646			}
647			PositionShader(const PositionShader&) = default;
648			~PositionShader() override = default;
649			Film::RGBA operator()(const Vec3R& xyz, const Vec3R&, const Vec3R&) const override
650			{
651			typename GridT::ValueType v = zeroVal<typename GridT::ValueType>();
652			SamplerType::sample(mAcc, mXform->worldToIndex(xyz), v);
653			const Vec3R rgb = (xyz - mMin) * mInvDim;
654			return Film::RGBA(v[0],v[1],v[2]) * Film::RGBA(rgb[0], rgb[1], rgb[2]);
655			}
656			BaseShader* copy() const override { return new PositionShader<GridT, SamplerType>(*this); }
657
658			private:
659			const Vec3R mMin, mInvDim;
660			typename GridT::ConstAccessor mAcc;
661			const math::Transform* mXform;
662			};
663
664			// Template specialization using a constant color of the material.
665			template<typename SamplerType>
666			class PositionShader<Film::RGBA, SamplerType>: public BaseShader
667			{
668			public:
669			PositionShader(const math::BBox<Vec3R>& bbox, const Film::RGBA& c = Film::RGBA(1.0f))
670			: mMin(bbox.min()), mInvDim(1.0/bbox.extents()), mRGBA(c) {}
671			PositionShader(const PositionShader&) = default;
672			~PositionShader() override = default;
673			Film::RGBA operator()(const Vec3R& xyz, const Vec3R&, const Vec3R&) const override
674			{
675			const Vec3R rgb = (xyz - mMin)*mInvDim;
676			return mRGBA*Film::RGBA(rgb[0], rgb[1], rgb[2]);
677			}
678			BaseShader* copy() const override { return new PositionShader<Film::RGBA, SamplerType>(*this); }
679
680			private:
681			const Vec3R mMin, mInvDim;
682			const Film::RGBA mRGBA;
683			};
684
685
686			/// @brief Simple diffuse Lambertian surface shader.
687			///
688			/// @details The diffuse color can either be constant (if GridT =
689			/// Film::RGBA which is the default) or defined in a separate Vec3
690			/// color grid. Lambertian implies that the (radiant) intensity is
691			/// directly proportional to the cosine of the angle between the
692			/// surface normal and the direction of the light source. Use
693			/// SamplerType to define the order of interpolation (default is
694			/// zero order, i.e. closes-point).
695			template<typename GridT = Film::RGBA,
696			typename SamplerType = tools::PointSampler>
697			class DiffuseShader: public BaseShader
698			{
699			public:
700			DiffuseShader(const GridT& grid): mAcc(grid.getAccessor()), mXform(&grid.transform()) {}
701			DiffuseShader(const DiffuseShader&) = default;
702			~DiffuseShader() override = default;
703			Film::RGBA operator()(const Vec3R& xyz, const Vec3R& normal, const Vec3R& rayDir) const override
704			{
705			typename GridT::ValueType v = zeroVal<typename GridT::ValueType>();
706			SamplerType::sample(mAcc, mXform->worldToIndex(xyz), v);
707			// We take the abs of the dot product corresponding to having
708			// light sources at +/- rayDir, i.e., two-sided shading.
709			return Film::RGBA(v[0],v[1],v[2])
710			* static_cast<Film::RGBA::ValueT>(math::Abs(normal.dot(rayDir)));
711			}
712			BaseShader* copy() const override { return new DiffuseShader<GridT, SamplerType>(*this); }
713
714			private:
715			typename GridT::ConstAccessor mAcc;
716			const math::Transform* mXform;
717			};
718
719			// Template specialization using a constant color of the material.
720			template <typename SamplerType>
721			class DiffuseShader<Film::RGBA, SamplerType>: public BaseShader
722			{
723			public:
724			DiffuseShader(const Film::RGBA& d = Film::RGBA(1.0f)): mRGBA(d) {}
725			DiffuseShader(const DiffuseShader&) = default;
726			~DiffuseShader() override = default;
727			Film::RGBA operator()(const Vec3R&, const Vec3R& normal, const Vec3R& rayDir) const override
728			{
729			// We assume a single directional light source at the camera,
730			// so the cosine of the angle between the surface normal and the
731			// direction of the light source becomes the dot product of the
732			// surface normal and inverse direction of the ray. We also ignore
733			// negative dot products, corresponding to strict one-sided shading.
734			//return mRGBA * math::Max(0.0, normal.dot(-rayDir));
735
736			// We take the abs of the dot product corresponding to having
737			// light sources at +/- rayDir, i.e., two-sided shading.
738			return mRGBA * static_cast<Film::RGBA::ValueT>(math::Abs(normal.dot(rayDir)));
739			}
740			BaseShader* copy() const override { return new DiffuseShader<Film::RGBA, SamplerType>(*this); }
741
742			private:
743			const Film::RGBA mRGBA;
744			};
745
746
747			//////////////////////////////////////// RAYTRACER ////////////////////////////////////////
748
749			template<typename GridT>
750		✗	void rayTrace(const GridT& grid,
751			const BaseShader& shader,
752			BaseCamera& camera,
753			size_t pixelSamples,
754			unsigned int seed,
755			bool threaded)
756			{
757			LevelSetRayTracer<GridT, tools::LevelSetRayIntersector<GridT> >
758		✗	tracer(grid, shader, camera, pixelSamples, seed);
759		✗	tracer.render(threaded);
760			}
761
762
763			template<typename GridT, typename IntersectorT>
764		✗	void rayTrace(const GridT&,
765			const IntersectorT& inter,
766			const BaseShader& shader,
767			BaseCamera& camera,
768			size_t pixelSamples,
769			unsigned int seed,
770			bool threaded)
771			{
772		✗	LevelSetRayTracer<GridT, IntersectorT> tracer(inter, shader, camera, pixelSamples, seed);
773		✗	tracer.render(threaded);
774			}
775
776
777			//////////////////////////////////////// LevelSetRayTracer ////////////////////////////////////////
778
779
780			template<typename GridT, typename IntersectorT>
781		✗	inline LevelSetRayTracer<GridT, IntersectorT>::
782			LevelSetRayTracer(const GridT& grid,
783			const BaseShader& shader,
784			BaseCamera& camera,
785			size_t pixelSamples,
786			unsigned int seed)
787			: mIsMaster(true),
788			mRand(nullptr),
789			mInter(grid),
790		✗	mShader(shader.copy()),
791		✗	mCamera(&camera)
792			{
793		✗	this->setPixelSamples(pixelSamples, seed);
794			}
795
796			template<typename GridT, typename IntersectorT>
797		✗	inline LevelSetRayTracer<GridT, IntersectorT>::
798			LevelSetRayTracer(const IntersectorT& inter,
799			const BaseShader& shader,
800			BaseCamera& camera,
801			size_t pixelSamples,
802			unsigned int seed)
803			: mIsMaster(true),
804			mRand(nullptr),
805			mInter(inter),
806		✗	mShader(shader.copy()),
807		✗	mCamera(&camera)
808			{
809		✗	this->setPixelSamples(pixelSamples, seed);
810			}
811
812			template<typename GridT, typename IntersectorT>
813		✗	inline LevelSetRayTracer<GridT, IntersectorT>::
814			LevelSetRayTracer(const LevelSetRayTracer& other) :
815			mIsMaster(false),
816		✗	mRand(other.mRand),
817			mInter(other.mInter),
818		✗	mShader(other.mShader->copy()),
819		✗	mCamera(other.mCamera),
820		✗	mSubPixels(other.mSubPixels)
821			{
822			}
823
824			template<typename GridT, typename IntersectorT>
825		✗	inline LevelSetRayTracer<GridT, IntersectorT>::
826			~LevelSetRayTracer()
827			{
828		✗	if (mIsMaster) delete [] mRand;
829			}
830
831			template<typename GridT, typename IntersectorT>
832			inline void LevelSetRayTracer<GridT, IntersectorT>::
833			setGrid(const GridT& grid)
834			{
835			assert(mIsMaster);
836			mInter = IntersectorT(grid);
837			}
838
839			template<typename GridT, typename IntersectorT>
840			inline void LevelSetRayTracer<GridT, IntersectorT>::
841			setIntersector(const IntersectorT& inter)
842			{
843			assert(mIsMaster);
844			mInter = inter;
845			}
846
847			template<typename GridT, typename IntersectorT>
848			inline void LevelSetRayTracer<GridT, IntersectorT>::
849			setShader(const BaseShader& shader)
850			{
851			assert(mIsMaster);
852			mShader.reset(shader.copy());
853			}
854
855			template<typename GridT, typename IntersectorT>
856			inline void LevelSetRayTracer<GridT, IntersectorT>::
857			setCamera(BaseCamera& camera)
858			{
859			assert(mIsMaster);
860			mCamera = &camera;
861			}
862
863			template<typename GridT, typename IntersectorT>
864		✗	inline void LevelSetRayTracer<GridT, IntersectorT>::
865			setPixelSamples(size_t pixelSamples, unsigned int seed)
866			{
867		✗	assert(mIsMaster);
868		✗	if (pixelSamples == 0) {
869		✗	OPENVDB_THROW(ValueError, "pixelSamples must be larger than zero!");
870			}
871		✗	mSubPixels = pixelSamples - 1;
872		✗	delete [] mRand;
873		✗	if (mSubPixels > 0) {
874		✗	mRand = new double[16];
875			math::Rand01<double> rand(seed);//offsets for anti-aliaing by jittered super-sampling
876		✗	for (size_t i=0; i<16; ++i) mRand[i] = rand();
877			} else {
878		✗	mRand = nullptr;
879			}
880			}
881
882			template<typename GridT, typename IntersectorT>
883		✗	inline void LevelSetRayTracer<GridT, IntersectorT>::
884			render(bool threaded) const
885			{
886		✗	tbb::blocked_range<size_t> range(0, mCamera->height());
887		✗	threaded ? tbb::parallel_for(range, *this) : (*this)(range);
888			}
889
890			template<typename GridT, typename IntersectorT>
891		✗	inline void LevelSetRayTracer<GridT, IntersectorT>::
892			operator()(const tbb::blocked_range<size_t>& range) const
893			{
894			const BaseShader& shader = *mShader;
895			Vec3Type xyz, nml;
896		✗	const float frac = 1.0f / (1.0f + float(mSubPixels));
897		✗	for (size_t j=range.begin(), n=0, je = range.end(); j<je; ++j) {
898		✗	for (size_t i=0, ie = mCamera->width(); i<ie; ++i) {
899		✗	Film::RGBA& bg = mCamera->pixel(i,j);
900		✗	RayType ray = mCamera->getRay(i, j);//primary ray
901		✗	Film::RGBA c = mInter.intersectsWS(ray, xyz, nml) ? shader(xyz, nml, ray.dir()) : bg;
902		✗	for (size_t k=0; k<mSubPixels; ++k, n +=2 ) {
903		✗	ray = mCamera->getRay(i, j, mRand[n & 15], mRand[(n+1) & 15]);
904		✗	c += mInter.intersectsWS(ray, xyz, nml) ? shader(xyz, nml, ray.dir()) : bg;
905			}//loop over sub-pixels
906		✗	bg = c*frac;
907			}//loop over image height
908			}//loop over image width
909			}
910
911			//////////////////////////////////////// VolumeRender ////////////////////////////////////////
912
913			template<typename IntersectorT, typename SampleT>
914		✗	inline VolumeRender<IntersectorT, SampleT>::
915			VolumeRender(const IntersectorT& inter, BaseCamera& camera)
916			: mAccessor(inter.grid().getConstAccessor())
917			, mCamera(&camera)
918		✗	, mPrimary(new IntersectorT(inter))
919		✗	, mShadow(new IntersectorT(inter))
920			, mPrimaryStep(1.0)
921			, mShadowStep(3.0)
922			, mCutOff(0.005)
923			, mLightGain(0.2)
924		✗	, mLightDir(Vec3R(0.3, 0.3, 0).unit())
925			, mLightColor(0.7, 0.7, 0.7)
926			, mAbsorption(0.1)
927		✗	, mScattering(1.5)
928			{
929			}
930
931			template<typename IntersectorT, typename SampleT>
932		✗	inline VolumeRender<IntersectorT, SampleT>::
933			VolumeRender(const VolumeRender& other)
934			: mAccessor(other.mAccessor)
935		✗	, mCamera(other.mCamera)
936		✗	, mPrimary(new IntersectorT(*(other.mPrimary)))
937		✗	, mShadow(new IntersectorT(*(other.mShadow)))
938		✗	, mPrimaryStep(other.mPrimaryStep)
939		✗	, mShadowStep(other.mShadowStep)
940		✗	, mCutOff(other.mCutOff)
941		✗	, mLightGain(other.mLightGain)
942			, mLightDir(other.mLightDir)
943			, mLightColor(other.mLightColor)
944			, mAbsorption(other.mAbsorption)
945		✗	, mScattering(other.mScattering)
946			{
947			}
948
949			template<typename IntersectorT, typename SampleT>
950		✗	inline void VolumeRender<IntersectorT, SampleT>::
951			print(std::ostream& os, int verboseLevel)
952			{
953		✗	if (verboseLevel>0) {
954		✗	os << "\nPrimary step: " << mPrimaryStep
955		✗	<< "\nShadow step: " << mShadowStep
956		✗	<< "\nCutoff: " << mCutOff
957		✗	<< "\nLightGain: " << mLightGain
958		✗	<< "\nLightDir: " << mLightDir
959		✗	<< "\nLightColor: " << mLightColor
960		✗	<< "\nAbsorption: " << mAbsorption
961		✗	<< "\nScattering: " << mScattering << std::endl;
962			}
963		✗	mPrimary->print(os, verboseLevel);
964			}
965
966			template<typename IntersectorT, typename SampleT>
967		✗	inline void VolumeRender<IntersectorT, SampleT>::
968			setIntersector(const IntersectorT& inter)
969			{
970		✗	mPrimary.reset(new IntersectorT(inter));
971		✗	mShadow.reset( new IntersectorT(inter));
972			}
973
974			template<typename IntersectorT, typename SampleT>
975		✗	inline void VolumeRender<IntersectorT, SampleT>::
976			render(bool threaded) const
977			{
978		✗	tbb::blocked_range<size_t> range(0, mCamera->height());
979		✗	threaded ? tbb::parallel_for(range, *this) : (*this)(range);
980			}
981
982			template<typename IntersectorT, typename SampleT>
983		✗	inline void VolumeRender<IntersectorT, SampleT>::
984			operator()(const tbb::blocked_range<size_t>& range) const
985			{
986		✗	SamplerType sampler(mAccessor, mShadow->grid().transform());//light-weight wrapper
987
988			// Any variable prefixed with p (or s) means it's associated with a primary (or shadow) ray
989			const Vec3R extinction = -mScattering-mAbsorption, One(1.0);
990			const Vec3R albedo = mLightColor*mScattering/(mScattering+mAbsorption);//single scattering
991		✗	const Real sGain = mLightGain;//in-scattering along shadow ray
992		✗	const Real pStep = mPrimaryStep;//Integration step along primary ray in voxel units
993		✗	const Real sStep = mShadowStep;//Integration step along shadow ray in voxel units
994		✗	const Real cutoff = mCutOff;//Cutoff for density and transmittance
995
996			// For the sake of completeness we show how to use two different
997			// methods (hits/march) in VolumeRayIntersector that produce
998			// segments along the ray that intersects active values. Comment out
999			// the line below to use VolumeRayIntersector::march instead of
1000			// VolumeRayIntersector::hits.
1001			#define USE_HITS
1002			#ifdef USE_HITS
1003			std::vector<typename RayType::TimeSpan> pTS, sTS;
1004			//std::deque<typename RayType::TimeSpan> pTS, sTS;
1005			#endif
1006
1007		✗	RayType sRay(Vec3R(0), mLightDir);//Shadow ray
1008		✗	for (size_t j=range.begin(), je = range.end(); j<je; ++j) {
1009		✗	for (size_t i=0, ie = mCamera->width(); i<ie; ++i) {
1010		✗	Film::RGBA& bg = mCamera->pixel(i, j);
1011		✗	bg.a = bg.r = bg.g = bg.b = 0;
1012		✗	RayType pRay = mCamera->getRay(i, j);// Primary ray
1013		✗	if( !mPrimary->setWorldRay(pRay)) continue;
1014			Vec3R pTrans(1.0), pLumi(0.0);
1015			#ifndef USE_HITS
1016			Real pT0, pT1;
1017			while (mPrimary->march(pT0, pT1)) {
1018			for (Real pT = pStep*ceil(pT0/pStep); pT <= pT1; pT += pStep) {
1019			#else
1020			mPrimary->hits(pTS);
1021		✗	for (size_t k=0; k<pTS.size(); ++k) {
1022		✗	Real pT = pStep*ceil(pTS[k].t0/pStep), pT1=pTS[k].t1;
1023		✗	for (; pT <= pT1; pT += pStep) {
1024			#endif
1025		✗	Vec3R pPos = mPrimary->getWorldPos(pT);
1026		✗	const Real density = sampler.wsSample(pPos);
1027		✗	if (density < cutoff) continue;
1028		✗	const Vec3R dT = math::Exp(extinction * density * pStep);
1029			Vec3R sTrans(1.0);
1030			sRay.setEye(pPos);
1031		✗	if( !mShadow->setWorldRay(sRay)) continue;
1032			#ifndef USE_HITS
1033			Real sT0, sT1;
1034			while (mShadow->march(sT0, sT1)) {
1035			for (Real sT = sStep*ceil(sT0/sStep); sT <= sT1; sT+= sStep) {
1036			#else
1037			mShadow->hits(sTS);
1038		✗	for (size_t l=0; l<sTS.size(); ++l) {
1039		✗	Real sT = sStep*ceil(sTS[l].t0/sStep), sT1=sTS[l].t1;
1040		✗	for (; sT <= sT1; sT+= sStep) {
1041			#endif
1042		✗	const Real d = sampler.wsSample(mShadow->getWorldPos(sT));
1043		✗	if (d < cutoff) continue;
1044		✗	sTrans *= math::Exp(extinction * d * sStep/(1.0+sT*sGain));
1045		✗	if (sTrans.lengthSqr()<cutoff) goto Luminance;//Terminate sRay
1046			}//Integration over shadow segment
1047			}// Shadow ray march
1048		✗	Luminance:
1049		✗	pLumi += albedo * sTrans * pTrans * (One-dT);
1050			pTrans *= dT;
1051		✗	if (pTrans.lengthSqr()<cutoff) goto Pixel; // Terminate Ray
1052			}//Integration over primary segment
1053			}// Primary ray march
1054		✗	Pixel:
1055		✗	bg.r = static_cast<Film::RGBA::ValueT>(pLumi[0]);
1056		✗	bg.g = static_cast<Film::RGBA::ValueT>(pLumi[1]);
1057		✗	bg.b = static_cast<Film::RGBA::ValueT>(pLumi[2]);
1058		✗	bg.a = static_cast<Film::RGBA::ValueT>(1.0f - pTrans.sum()/3.0f);
1059			}//Horizontal pixel scan
1060			}//Vertical pixel scan
1061			}
1062
1063
1064			////////////////////////////////////////
1065
1066
1067			// Explicit Template Instantiation
1068
1069			#ifdef OPENVDB_USE_EXPLICIT_INSTANTIATION
1070
1071			#ifdef OPENVDB_INSTANTIATE_RAYTRACER
1072			#include <openvdb/util/ExplicitInstantiation.h>
1073			#endif
1074
1075			#define _FUNCTION(TreeT) \
1076			void rayTrace(const Grid<TreeT>&, const BaseShader&, BaseCamera&, size_t, unsigned int, bool)
1077			OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
1078			#undef _FUNCTION
1079
1080			#define _FUNCTION(TreeT) \
1081			void rayTrace(const Grid<TreeT>&, const tools::LevelSetRayIntersector<Grid<TreeT>>&, const BaseShader&, BaseCamera&, size_t, unsigned int, bool)
1082			OPENVDB_REAL_TREE_INSTANTIATE(_FUNCTION)
1083			#undef _FUNCTION
1084
1085			OPENVDB_INSTANTIATE_CLASS VolumeRender<tools::VolumeRayIntersector<FloatGrid>, tools::BoxSampler>;
1086			OPENVDB_INSTANTIATE_CLASS VolumeRender<tools::VolumeRayIntersector<DoubleGrid>, tools::BoxSampler>;
1087
1088			#endif // OPENVDB_USE_EXPLICIT_INSTANTIATION
1089
1090
1091			} // namespace tools
1092			} // namespace OPENVDB_VERSION_NAME
1093			} // namespace openvdb
1094
1095			#endif // OPENVDB_TOOLS_RAYTRACER_HAS_BEEN_INCLUDED
1096