# Note on “Ray tracing in one week” - Part 6
Dielectrics, rather than the definition in circuits, is material that reflects and refracts light. The most common example is glass. We have already implemented reflection, and now we are going to implement refraction.
# Snell’s law and refraction
When calculating refracted ray, we use Snell’s law. Snell’s law stated that: given two medium with refractive index η \eta η and η ′ \eta' η ′ , a ray with incidence angle \theta
and refraction angle \theta'
will satisfy the following equation:
η sin ( θ ) = η ′ sin ( θ ′ ) \eta \sin(\theta) = \eta' \sin(\theta')
η sin ( θ ) = η ′ sin ( θ ′ )
The diagram below shows the refraction of light from air to glass(assumed to have refractive index of 1.5). I \mathbf{I} I is the incident ray, and R \mathbf{R} R is the refracted ray. We can derive the vector formula of R with I and the normal vector n \mathbf{n} n of the surface.
Since we only need the direction of the refracted ray, we can treat the incident ray and the normal vector as unit vectors. By simple trig, we have:
I x = sin ( θ ) ⋅ N R x = sin ( θ ′ ) ⋅ N \begin{align*}
\mathbf{I}_x &= \sin(\theta) \cdot\mathbf{N}\\
\mathbf{R}_x &= \sin(\theta') \cdot\mathbf{N}\\
\end{align*}
I x R x = sin ( θ ) ⋅ N = sin ( θ ′ ) ⋅ N
And by Snell’s law, we have:
η sin ( θ ) = η ′ sin ( θ ′ ) sin ( θ ′ ) = η η ′ ⋅ sin ( θ ) \begin{align*}
\eta \sin(\theta) &= \eta' \sin(\theta') \\
\sin(\theta') &= \frac{\eta}{\eta'}\cdot \sin(\theta) \\
\end{align*}
η sin ( θ ) sin ( θ ′ ) = η ′ sin ( θ ′ ) = η ′ η ⋅ sin ( θ )
Substitute, we have:
R x = sin ( θ ′ ) R x = η η ′ ⋅ sin ( θ ) R x = η η ′ ⋅ I x \begin{align*}
\mathbf{R}_x &= \sin(\theta') \\
\mathbf{R}_x &= \frac{\eta}{\eta'}\cdot \sin(\theta) \\
\mathbf{R}_x &= \frac{\eta}{\eta'}\cdot \mathbf{I}_x \\
\end{align*}
R x R x R x = sin ( θ ′ ) = η ′ η ⋅ sin ( θ ) = η ′ η ⋅ I x
By vector addition and geometry, we have:
I x = I + n y = I + N ⋅ cos θ \mathbf{I}_x = \mathbf{I} + \mathbf{n}_y = \mathbf{I} + \mathbf{N} \cdot \cos{\theta}
I x = I + n y = I + N ⋅ cos θ
Substitute, we have:
R x = η η ′ ⋅ I x = η η ′ ⋅ ( I + N ⋅ cos θ ) \begin{align*}
\mathbf{R}_x &= \frac{\eta}{\eta'}\cdot \mathbf{I}_x \\
&= \frac{\eta}{\eta'}\cdot (\mathbf{I} + \mathbf{N} \cdot \cos{\theta}) \\
\end{align*}
R x = η ′ η ⋅ I x = η ′ η ⋅ ( I + N ⋅ cos θ )
And we can derive the R y \mathbf{R_y} R y by using Pythagorean theorem:
R y = − 1 − ∣ ∣ R ∣ ∣ x 2 ⋅ N \mathbf{R}_y = -\sqrt{1 - ||\mathbf{R}||_x^2}\cdot\mathbf{N}
R y = − 1 − ∣∣ R ∣ ∣ x 2 ⋅ N
This is good, but cos \cos cos for θ \theta θ is not easy to calculate. We can use the following formula to calculate cos θ \cos{\theta} cos θ if theta is the angle between two normalized vectors A \mathbf{A} A and B \mathbf{B} B :
cos θ = A ⋅ B \cos{\theta} = \mathbf{A} \cdot \mathbf{B}
cos θ = A ⋅ B
and R x \mathbf{R}_x R x is:
R x = η η ′ ⋅ ( I + N ⋅ cos θ ) = η η ′ ⋅ ( I + N ⋅ ⟨ I , N ⟩ ) \begin{align*}
\mathbf{R}_x &= \frac{\eta}{\eta'}\cdot (\mathbf{I} + \mathbf{N} \cdot \cos{\theta}) \\
&= \frac{\eta}{\eta'}\cdot (\mathbf{I} + \mathbf{N} \cdot \langle \mathbf{I} ,\mathbf{N}\rangle )\\
\end{align*}
R x = η ′ η ⋅ ( I + N ⋅ cos θ ) = η ′ η ⋅ ( I + N ⋅ ⟨ I , N ⟩)
We can implement this in code like this:
1 2 3 4 5 6 7 class Vec3 { static Vec3 refract (const Vec3 &uv, const Vec3 &n, double etai_over_etat) ; };
1 2 3 4 5 6 Vec3 Vec3::refract (const Vec3 &uv, const Vec3 &n, double etai_over_etat) { auto cos_theta = fmin ((-uv).dot (n), 1.0 ); Vec3 refr_x = etai_over_etat * (uv * cos_theta * n); Vec3 refr_y = -sqrt (fabs (1.0 - refr_x.lengthSq ())) * n; return refr_x + refr_y; }
and now we define our dielectric material:
1 2 3 4 5 6 7 8 9 10 11 class Dielectric : public IMaterial {public : Dielectric (double idx, const ImageUtil::Color &albedo); bool scatter (const MathUtil::Ray &r_in, const HitRecord &record, MathUtil::Vec3 &attenuation, MathUtil::Ray &scattered) const override ;private : double ir; ImageUtil::Color albedo; };
1 2 3 4 5 6 7 8 9 bool Dielectric::scatter (const MathUtil::Ray &r_in, const HitRecord &record, MathUtil::Vec3 &attenuation, MathUtil::Ray &scattered) const { attenuation = albedo; double ref_ratio = record.front_face ? (1.0 / ir) : ir; auto refracted = MathUtil::Vec3::refract (r_in.dir ().unit (), record.normal, ref_ratio); scattered = MathUtil::Ray (record.p, refracted); return true ; }
We add some code to the main function to test the dielectric material:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 #include "ImageUtil.h" #include "GraphicObjects.h" #include "Material.h" int main () { auto camera = Camera (1920 , 16.0 / 9.0 , 2.0 , 1 , MathUtil::Point3 (0 , 0 , 1 )); auto world = HittableList (); auto center_ball_material = std::make_shared <Dielectric>(Dielectric (1.5 , ImageUtil::Color (1.0 , 1.0 , 1.0 ))); auto right_ball_material = std::make_shared <Metal>(Metal (ImageUtil::Color (0.7 , 0.3 , 0.3 ), 0.4 )); auto left_ball_material = std::make_shared <Metal>(Metal (ImageUtil::Color (0.8 , 0.8 , 0.8 ), 0.4 )); auto ground_material = std::make_shared <Metal>(Metal (ImageUtil::Color (0.8 , 0.6 , 0.2 ), 0.4 )); world.add (std::make_shared <Sphere>(Sphere (1000 , {0 , -1000.5 , -1.0 }, ground_material))); world.add (std::make_shared <Sphere>(Sphere (0.5 , {0 , 0 , -1 }, center_ball_material))); world.add (std::make_shared <Sphere>(Sphere (0.5 , {-1 , 0 , -1.0 }, left_ball_material))); world.add (std::make_shared <Sphere>(Sphere (0.5 , {1 , 0 , -1.0 }, right_ball_material))); camera.setSampleCount (100 ); camera.setRenderDepth (50 ); camera.setRenderThreadCount (24 ); camera.Render (world, "out" , "test.ppm" ); return 0 ; }
The result is:
# Full reflection
Consider the equation deduced from Snell’s law:
sin θ η η ′ = sin θ ′ \sin{\theta} \frac{\eta}{\eta'} = \sin{\theta'}
sin θ η ′ η = sin θ ′
When η > η ′ \eta > \eta' η > η ′ , for some value of θ \theta θ , the right-hand side will be greater than 1, and there will be no real solution for θ ′ \theta' θ ′ . This means that the light will be fully reflected.
Since sin \sin sin takes too much time to calculate, we can use our old friend, the dot product, to calculate sin θ \sin{\theta} sin θ :
cos θ = I ⋅ N sin θ = 1 − cos 2 θ = 1 − ( I ⋅ N ) 2 \begin{align*}
\cos{\theta} &= \mathbf{I} \cdot \mathbf{N} \\
\sin{\theta} &= \sqrt{1 - \cos^2{\theta}} \\
&= \sqrt{1 - (\mathbf{I} \cdot \mathbf{N})^2} \\
\end{align*}
cos θ sin θ = I ⋅ N = 1 − cos 2 θ = 1 − ( I ⋅ N ) 2
We can use this to calculate the reflectance of the dielectric material:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 bool Dielectric::scatter (const MathUtil::Ray &r_in, const HitRecord &record, MathUtil::Vec3 &attenuation, MathUtil::Ray &scattered) const { attenuation = albedo; double ref_ratio = record.front_face ? (1.0 / ir) : ir; auto unit = r_in.dir ().unit (); double cos = fmin ((-unit).dot (record.normal), 1.0 ); double sin = sqrt (1 - cos * cos); bool can_refr = ref_ratio * sin < 1.0 ; MathUtil::Vec3 dir; if (can_refr) dir = MathUtil::Vec3::refract (r_in.dir ().unit (), record.normal, ref_ratio); else dir = MathUtil::Vec3::reflect (r_in.dir ().unit (), record.normal); scattered = MathUtil::Ray (record.p, dir); return true ; }
The result:
# Schlick’s approximation
Normally, real glass becomes more reflective as the angle of incidence increases. This is called Fresnel effect. Schlick’s approximation is a simple way to approximate the reflectance of a dielectric material. The formula is:
R ( θ ) = R 0 + ( 1 − R 0 ) ( 1 − cos θ ) 5 R(\theta) = R_0 + (1 - R_0)(1 - \cos{\theta})^5
R ( θ ) = R 0 + ( 1 − R 0 ) ( 1 − cos θ ) 5
where:
R 0 = ( η − η ′ η + η ′ ) 2 R_0 = (\frac{\eta - \eta'}{\eta + \eta'})^2
R 0 = ( η + η ′ η − η ′ ) 2
and θ \theta θ is the incident angle.
We can implement this in our code:
1 2 3 4 5 6 7 8 9 10 11 12 class Dielectric : public IMaterial {public : Dielectric (double idx, const ImageUtil::Color &albedo); bool scatter (const MathUtil::Ray &r_in, const HitRecord &record, MathUtil::Vec3 &attenuation, MathUtil::Ray &scattered) const override ;private : double ir; ImageUtil::Color albedo; static double reflectance (double cosine, double refr_idx) ; };
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 bool Dielectric::scatter (const MathUtil::Ray &r_in, const HitRecord &record, MathUtil::Vec3 &attenuation, MathUtil::Ray &scattered) const { attenuation = albedo; double ref_ratio = record.front_face ? (1.0 / ir) : ir; auto unit = r_in.dir ().unit (); double cos = fmin ((-unit).dot (record.normal), 1.0 ); double sin = sqrt (1 - cos * cos); bool can_refr = ref_ratio * sin < 1.0 ; MathUtil::Vec3 dir; if (can_refr && reflectance (cos, ref_ratio) < MathUtil::randomDouble ()) dir = MathUtil::Vec3::refract (r_in.dir ().unit (), record.normal, ref_ratio); else dir = MathUtil::Vec3::reflect (r_in.dir ().unit (), record.normal); scattered = MathUtil::Ray (record.p, dir); return true ; } double Dielectric::reflectance (double cosine, double refr_idx) { auto r0 = (1 - refr_idx) / (1 + refr_idx); r0 *= r0; return r0 + (1 - r0) * std::pow (1 - cosine, 5 ); }
The result:
Notice the edge of the glass sphere. The reflection is stronger than the center and you can see the reflection of the other spheres.
# Zooming in
Right now we have the ability to position our camera, but not yet to rotate or change the zoom. Let’s implement the second feature first, and it can be implemented by adjusting the field of view.
Field of view is the max angle of the ray that the camera can see. Give a FOV and the aspect ratio, we can calculate the height and width of the viewport. The formula is:
V h e i g h t = 2 ⋅ tan ( θ 2 ) V w i d t h = V h e i g h t ⋅ R \begin{align*}
V_{height} &= 2 \cdot \tan(\frac{\theta}{2})\\
V_{width} &= V_{height} \cdot R
\end{align*}
V h e i g h t V w i d t h = 2 ⋅ tan ( 2 θ ) = V h e i g h t ⋅ R
Where V V V is the viewport, θ \theta θ is the field of view, and R R R is the aspect ratio.
If we view the yz plane from x direction, the fov will be like this:
We can implement this in our camera class:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 class Camera {public : Camera (int width, double aspect_ratio, double fov, double focal_len, MathUtil::Point3 position); int getWidth () const ; int getHeight () const ; double getAspectRatio () const ; double getViewportWidth () const ; double getViewportHeight () const ; double getFocalLen () const ; int getSampleCount () const ; double getFov () const ; const MathUtil::Point3 &getPosition () const ; const MathUtil::Vec3 &getHoriVec () const ; const MathUtil::Vec3 &getVertVec () const ; const MathUtil::Vec3 &getPixDeltaX () const ; const MathUtil::Vec3 &getPixDeltaY () const ; const MathUtil::Point3 &getViewportUl () const ; const MathUtil::Point3 &getPixel00 () const ; void setWidth (int width) ; void setAspectRatio (double aspect_ratio) ; void setFocalLen (double focal_len) ; void setPosition (const MathUtil::Point3 &position) ; void setSampleCount (int sample_count) ; void setFov (double fov) ; [[nodiscard]] MathUtil::Vec3 getPixelVec (int x, int y) const ; [[nodiscard]] MathUtil::Vec3 getPixRayDir (int x, int y) const ; void Render (const IHittable& world, const std::string& path, const std::string& name) ; int getRenderDepth () const ; void setRenderDepth (int renderDepth) ; int getRenderThreadCount () const ; void setRenderThreadCount (int renderThreadCount) ; private : std::vector<std::vector<ImageUtil::Color>> RenderWorker (const IHittable &world, int start, int end); ImageUtil::Color rayColor (const MathUtil::Ray &ray, const IHittable &object, int depth) ; void updateVectors () ; MathUtil::Vec3 randomDisplacement () const ; int width; int height; double aspect_ratio; double viewport_width; double viewport_height; double focal_len; double fov = 45 ; int sample_count = 20 ; int render_depth = 50 ; int render_thread_count = 0 ; MathUtil::Point3 position; MathUtil::Vec3 hori_vec; MathUtil::Vec3 vert_vec; MathUtil::Vec3 pix_delta_x; MathUtil::Vec3 pix_delta_y; MathUtil::Point3 viewport_ul; MathUtil::Point3 pixel_00; };
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 Camera::Camera (int width, double aspect_ratio, double fov, double focal_len, MathUtil::Point3 position) : width (width), aspect_ratio (aspect_ratio), fov (fov), focal_len (focal_len), position (std::move (position)), height (static_cast <int >(width / aspect_ratio)),viewport_height (viewport_width / (static_cast <double >(width) / static_cast <double >(height))) { updateVectors (); } void Camera::updateVectors () { auto theta = deg2Rad (fov); auto h = tan (theta / 2 ); viewport_height = 2 * h * focal_len; viewport_width = viewport_height * (static_cast <double >(width) / height); hori_vec = MathUtil::Vec3 (viewport_width, 0 , 0 ); vert_vec = MathUtil::Vec3 (0 , -viewport_height, 0 ); pix_delta_x = hori_vec / width; pix_delta_y = vert_vec / height; viewport_ul = position - MathUtil::Vec3 (0 , 0 , focal_len) - (vert_vec + hori_vec) / 2 ; pixel_00 = viewport_ul + (pix_delta_y + pix_delta_x) * 0.5 ; } void Camera::setAspectRatio (double aspect_ratio) { Camera::aspect_ratio = aspect_ratio; height = static_cast <int >(width / aspect_ratio); updateVectors (); } void Camera::setFov (double fov) { this ->fov = fov; } double Camera::getFov () const { return fov; }
Now we modify our main function to see the result:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 #include "ImageUtil.h" #include "GraphicObjects.h" #include "Material.h" int main () { auto camera = Camera (1920 , 16.0 / 9.0 , 45 , 1 , MathUtil::Point3 (0 , 0 , 1 )); auto world = HittableList (); auto center_ball_material = std::make_shared <Metal>(Metal (ImageUtil::Color (0.8 , 0.8 , 0.8 ), 0.4 )); auto ground_material = std::make_shared <Metal>(Metal (ImageUtil::Color (0.8 , 0.6 , 0.2 ), 0.4 )); world.add (std::make_shared <Sphere>(Sphere (1000 , {0 , -1000.5 , -1.0 }, ground_material))); world.add (std::make_shared <Sphere>(Sphere (0.5 , {0 , 0 , -1 }, center_ball_material))); camera.setSampleCount (100 ); camera.setRenderDepth (50 ); camera.setRenderThreadCount (24 ); camera.Render (world, "out" , "test.ppm" ); return 0 ; }
The result:
FOV = 90:
FOV = 45: