1 /*

2 Simple introductory ray tracer

3 Copyright (C) 2012 John Tsiombikas <nuclear@member.fsf.org>

4

5 This program is free software: you can redistribute it and/or modify

6 it under the terms of the GNU General Public License as published by

7 the Free Software Foundation, either version 3 of the License, or

8 (at your option) any later version.

9

10 This program is distributed in the hope that it will be useful,

11 but WITHOUT ANY WARRANTY; without even the implied warranty of

12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

13 GNU General Public License for more details.

14

15 You should have received a copy of the GNU General Public License

16 along with this program. If not, see <http://www.gnu.org/licenses/>.

17 */

18 #include "cone.h"

19

20 Cone::Cone()

21 {

22 angle = M_PI / 4.0;

23 ymin = 0.0f;

24 ymax = 1.0f;

25 }

26

27 Cone::Cone(float angle, float ymin, float ymax)

28 {

29 this->angle = angle;

30 this->ymin = ymin;

31 this->ymax = ymax;

32 }

33

34 #define SQ(x) ((x) * (x))

35 bool Cone::intersect(const Ray &inray, HitPoint *pt) const

36 {

37 Ray ray = inray.transformed(inv_xform);

38

39 float slope = 2.0 * angle / M_PI;

40

41 float a = SQ(ray.dir.x) - SQ(slope * ray.dir.y) + SQ(ray.dir.z);

42 float b = 2.0 * ray.origin.x * ray.dir.x - 2.0 * SQ(slope) * ray.origin.y * ray.dir.y +

43 2.0 * ray.origin.z * ray.dir.z;

44 float c = SQ(ray.origin.x) - SQ(slope * ray.origin.y) + SQ(ray.origin.z);

45

46 float discr = b * b - 4.0 * a * c;

47 if(discr < 1e-4)

48 return false;

49

50 float sqrt_discr = sqrt(discr);

51 float t0 = (-b + sqrt_discr) / (2.0 * a);

52 float t1 = (-b - sqrt_discr) / (2.0 * a);

53

54 if(t0 < 1e-4)

55 t0 = t1;

56 if(t1 < 1e-4)

57 t1 = t0;

58

59 float t = t0 < t1 ? t0 : t1;

60 if(t < 1e-4)

61 return false;

62

63 // fill the HitPoint structure

64 pt->obj = this;

65 pt->dist = t;

66 pt->pos = ray.origin + ray.dir * t;

67

68 if(pt->pos.y < ymin || pt->pos.y >= ymax) {

69 return false;

70 }

71

72 float radius = pt->pos.y * slope;

73 pt->normal = Vector3(pt->pos.x, 0.0, pt->pos.z) / radius;

74

75 float ny = -sin(angle);

76 float xzlen = sqrt(1.0 - ny * ny);

77 pt->normal.x *= xzlen;

78 pt->normal.y = ny;

79 pt->normal.z *= xzlen;

80

81 pt->pos.transform(xform);

82 pt->normal.transform(dir_xform);

83 return true;

84 }