* aclib #20 - kray.c の紹介ページ -(by [[K]], 2020.09.23) ** (1) 基本情報 -krayは、[[K]]がaclライブラリのデモのために用意したプログラムです。 -オリジナルは https://gist.github.com/yhara/ea0e66e0d8bdd114d2401dd133539fa3 です。 ** (2) 実行画面 -http://k.osask.jp/files/pic20200923a.png ** (3) ソースコード #include <acl.c> #define EPS 1.0e-4 typedef struct Vec_ { double x, y, z; } Vec; Vec Vec_new(double x, double y, double z) { Vec v; v.x = x; v.y = y; v.z = z; Vec v = { x, y, z }; return v; } double Util_clamp(double t, double min, double max) { if (t < min) t = min; if (t > max) t = max; return t; } Vec Vec_add(Vec a, Vec b) { return Vec_new(a.x + b.x, a.y + b.y, a.z + b.z); } Vec Vec_sub(Vec a, Vec b) { return Vec_new(a.x - b.x, a.y - b.y, a.z - b.z); } Vec Vec_mul(double t, Vec v) { return Vec_new(t * v.x, t * v.y, t * v.z); } Vec Vec_multi(Vec a, Vec b) { return Vec_new(a.x * b.x, a.y * b.y, a.z * b.z); } double Vec_dot(Vec a, Vec b) { return a.x * b.x + a.y * b.y + a.z * b.z; } double Vec_length(Vec v) { return sqrt(Vec_dot(v, v)); } Vec Vec_reflect(Vec v, Vec normal) { return Vec_add(v, Vec_mul(-2.0 * Vec_dot(normal, v), normal)); } AInt32 Util_color(double t) { return (AInt32) (255.99999 * Util_clamp(t, 0.0, 1.0)); } AInt32 Util_rgb(Vec v) { return aRgb8(Util_color(v.x), Util_color(v.y), Util_color(v.z)); } Vec Vec_normalize(Vec v) { double len = Vec_length(v); if (len > 1.0e-17) v = Vec_mul(1.0 / len, v); return v; } typedef struct Ray_ { Vec origin, dir; } Ray; Ray Ray_new(Vec origin, Vec dir) { Ray r; r.origin = origin; r.dir = dir; return r; } typedef struct Isect_ { Vec hit_point, normal, color, ray_dir; Vec hit_point, normal, color; double distance; } Isect; Isect Isect_new(Vec hit_point, Vec normal, Vec color, double distance, Vec ray_dir) { Isect i; i.hit_point = hit_point; i.normal = normal; i.color = color; i.distance = distance; i.ray_dir = ray_dir; return i; } double Util_mod2(double t) { t -= (AInt32) (t * (1.0 / 2.0)) * 2.0; if (t < 0.0) t += 2.0; return t; } typedef struct Sphere_ { double radius; Vec position, color; } Sphere; typedef struct Object_ { Vec pos, col, nor; double rad; // radius } Object; // Sphere(pos, col, rad) or Plane(pos, col, nor). Sphere Sphere_new(double radius, Vec position, Vec color) void Sphere_intersect(Object s, Vec ray_origin, Vec ray_dir, Vec light, Isect *i) { Sphere s; s.radius = radius; s.position = position; s.color = color; return s; } void Sphere_intersect(Sphere s, Ray ray, Vec light, Isect *i) { Vec rs = Vec_sub(ray.origin, s.position); double b = Vec_dot(rs, ray.dir); double c = Vec_dot(rs, rs) - s.radius * s.radius; Vec rs = Vec_sub(ray_origin, s.pos); double b = Vec_dot(rs, ray_dir); double c = Vec_dot(rs, rs) - s.rad * s.rad; double d = b * b - c; if (d < 0.0) return; double t = - b - sqrt(d); if (t < EPS || t > i->distance) return; i->hit_point = Vec_add(ray.origin, Vec_mul(t, ray.dir)); i->normal = Vec_normalize(Vec_sub(i->hit_point, s.position)); i->color = Vec_mul(Util_clamp(Vec_dot(light, i->normal), 0.1, 1.0), s.color); i->hit_point = Vec_add(ray_origin, Vec_mul(t, ray_dir)); i->normal = Vec_normalize(Vec_sub(i->hit_point, s.pos)); i->color = Vec_mul(Util_clamp(Vec_dot(light, i->normal), 0.1, 1.0), s.col); i->distance = t; i->ray_dir = ray.dir; } typedef struct Plane_ { Vec position, normal, color; } Plane; Plane Plane_new(Vec position, Vec normal, Vec color) void Plane_intersect(Object p, Vec ray_origin, Vec ray_dir, Vec light, Isect *i) { Plane p; p.position = position; p.normal = normal; p.color = color; return p; } void Plane_intersect(Plane p, Ray ray, Vec light, Isect *i) { double d = - Vec_dot(p.position, p.normal); double v = Vec_dot(ray.dir, p.normal); double d = - Vec_dot(p.pos, p.nor); double v = Vec_dot(ray_dir, p.nor); if (v * v < 1.0e-30) return; double t = - (Vec_dot(ray.origin, p.normal) + d) / v; double t = - (Vec_dot(ray_origin, p.nor) + d) / v; if (t < EPS || t > i->distance) return; i->hit_point = Vec_add(ray.origin, Vec_mul(t, ray.dir)); i->normal = p.normal; i->hit_point = Vec_add(ray_origin, Vec_mul(t, ray_dir)); i->normal = p.nor; double d2 = Util_clamp(Vec_dot(light, i->normal), 0.1, 1.0); if ((Util_mod2(i->hit_point.x) - 1.0) * (Util_mod2(i->hit_point.z) - 1.0) > 0.0) d2 *= 0.5; i->color = Vec_mul(d2 * (1.0 - Util_clamp(fabs(i->hit_point.z) * 0.04, 0.0, 1.0)), p.color); i->color = Vec_mul(d2 * (1.0 - Util_clamp(fabs(i->hit_point.z) * 0.04, 0.0, 1.0)), p.col); i->distance = t; i->ray_dir = ray.dir; } typedef struct Util_ { Vec light; Sphere s1, s2, s3; Plane p; Object s1, s2, s3, p; } Util; void Util_intersect(Util u, Ray ray, Isect *i) void Util_intersect(Util u, Vec ray_origin, Vec ray_dir, Isect *i) { i->distance = 1.0e+30; Sphere_intersect(u.s1, ray, u.light, i); Sphere_intersect(u.s2, ray, u.light, i); Sphere_intersect(u.s3, ray, u.light, i); Plane_intersect (u.p, ray, u.light, i); Sphere_intersect(u.s1, ray_origin, ray_dir, u.light, i); Sphere_intersect(u.s2, ray_origin, ray_dir, u.light, i); Sphere_intersect(u.s3, ray_origin, ray_dir, u.light, i); Plane_intersect (u.p, ray_origin, ray_dir, u.light, i); } void aMain() { Util u; Isect i; u.light = Vec_new(0.577, 0.577, 0.577); u.s1 = Sphere_new(0.5, Vec_new( 0.0, -0.5, 0.0), Vec_new(1.0, 0.0, 0.0)); u.s2 = Sphere_new(1.0, Vec_new( 2.0, 0.0, cos(6.66)), Vec_new(0.0, 1.0, 0.0)); u.s3 = Sphere_new(1.5, Vec_new(-2.0, 0.5, cos(3.33)), Vec_new(0.0, 0.0, 1.0)); u.p = Plane_new(Vec_new(0.0, -1.0, 0.0), Vec_new(0.0, 1.0, 0.0), Vec_new(1.0, 1.0, 1.0)); u.s1.rad = 0.5; u.s1.pos = Vec_new( 0.0, -0.5, 0.0); u.s1.col = Vec_new(1.0, 0.0, 0.0); u.s2.rad = 1.0; u.s2.pos = Vec_new( 2.0, 0.0, cos(6.66)); u.s2.col = Vec_new(0.0, 1.0, 0.0); u.s3.rad = 1.5; u.s3.pos = Vec_new(-2.0, 0.5, cos(3.33)); u.s3.col = Vec_new(0.0, 0.0, 1.0); u.p.pos = Vec_new(0.0, -1.0, 0.0); u.p.nor = Vec_new(0.0, 1.0, 0.0); u.p.col = Vec_new(1.0, 1.0, 1.0); AWindow *win = aOpenWin(512, 384, "kray", 1); AInt16 ix, iy; AInt16 ix, iy, j; for (iy = 0; iy < 384; iy++) { for (ix = 0; ix < 512; ix++) { double x = ix * (1.0 / 256.0) - 1.0; double y = (384 - iy) * (1.0 / 256.0) - 1.0; Ray ray = Ray_new(Vec_new(0.0, 2.0, 6.0), Vec_normalize(Vec_new(x, y, -1.0))); Util_intersect(u, ray, &i); Vec dest_col = Vec_mul(ray.dir.y, Vec_new(1.0, 1.0, 1.0)); Vec ray_dir = Vec_normalize(Vec_new(x, y, -1.0)); Util_intersect(u, Vec_new(0.0, 2.0, 6.0), ray_dir, &i); Vec dest_col = Vec_mul(ray_dir.y, Vec_new(1.0, 1.0, 1.0)); if (i.distance < 1.0e+30) { Vec temp_col = dest_col = i.color; AInt8 j; for (j = 1; j < 4; j++) { Ray q = Ray_new(i.hit_point, Vec_reflect(i.ray_dir, i.normal)); Util_intersect(u, q, &i); ray_dir = Vec_reflect(ray_dir, i.normal); Util_intersect(u, i.hit_point, ray_dir, &i); if (i.distance >= 1.0e+30) break; dest_col = Vec_add(dest_col, Vec_multi(temp_col, i.color)); temp_col = Vec_multi(temp_col, i.color); dest_col = Vec_add(dest_col, temp_col); } } aSetPix(win, ix, iy, Util_rgb(dest_col)); } aLeapFlushAll(win, 100); } aWait(-1); } ** (4) 備考 -構造体を引数に渡すとき、ポインタ渡しにしていないので、おそらく無駄なコピーが発生し、速度的には不利だろうと思われます。 -しかしその辺をきちんと書くと、行数が数十行増えてしまうので、あえてこの書き方にしています。 -mandel.cと比べればそれほど遅くはないので、まあこれでもいいかなと思っています。 * こめんと欄 #comment