--- /dev/null
+*.o
+*.d
+*.swp
+ray24
--- /dev/null
+obj = src/main.o src/geom.o src/rend.o
+dep = $(obj:.o=.d)
+bin = ray24
+
+CFLAGS = -pedantic -Wall -g -MMD `sdl-config --cflags`
+LDFLAGS = `sdl-config --libs` -lm
+
+$(bin): $(obj)
+ $(CC) -o $@ $(obj) $(LDFLAGS)
+
+-include $(dep)
+
+.c.o:
+ $(CC) -o $@ $(CFLAGS) -c $<
+
+.PHONY: clean
+clean:
+ rm -f $(obj) $(bin)
+
+.PHONY: cleandep
+cleandep:
+ rm -f $(dep)
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ *
+ * Function prefixes signify the data type of their operand(s):
+ * - cgm_v... functions are operations on cgm_vec3 vectors
+ * - cgm_w... functions are operations on cgm_vec4 vectors
+ * - cgm_q... functions are operations on cgm_quat quaternions (w + xi + yj + zk)
+ * - cgm_m... functions are operations on 4x4 matrices (stored as linear 16 float arrays)
+ * - cgm_r... functions are operations on cgm_ray rays
+ *
+ * NOTE: *ALL* matrix arguments are pointers to 16 floats. Even the functions
+ * which operate on 3x3 matrices, actually use the upper 3x3 of a 4x4 matrix,
+ * and still expect an array of 16 floats.
+ *
+ * NOTE: matrices are treated by all operations as column-major, to match OpenGL
+ * conventions, so everything is pretty much transposed.
+*/
+#ifndef CGMATH_H_
+#define CGMATH_H_
+
+#include <math.h>
+#include <string.h>
+
+#define CGM_PI 3.141592653589793
+
+typedef struct {
+ float x, y;
+} cgm_vec2;
+
+typedef struct {
+ float x, y, z;
+} cgm_vec3;
+
+typedef struct {
+ float x, y, z, w;
+} cgm_vec4, cgm_quat;
+
+typedef struct {
+ cgm_vec3 origin, dir;
+} cgm_ray;
+
+typedef enum cgm_euler_mode {
+ CGM_EULER_XYZ,
+ CGM_EULER_XZY,
+ CGM_EULER_YXZ,
+ CGM_EULER_YZX,
+ CGM_EULER_ZXY,
+ CGM_EULER_ZYX,
+ CGM_EULER_ZXZ,
+ CGM_EULER_ZYZ,
+ CGM_EULER_YXY,
+ CGM_EULER_YZY,
+ CGM_EULER_XYX,
+ CGM_EULER_XZX
+} cgm_euler_mode;
+
+#ifdef __cplusplus
+#define CGM_INLINE inline
+
+extern "C" {
+#else
+
+#if (__STDC_VERSION__ >= 199901) || defined(__GNUC__)
+#define CGM_INLINE inline
+#else
+#define CGM_INLINE __inline
+#endif
+
+#endif
+
+/* --- operations on cgm_vec3 --- */
+static CGM_INLINE void cgm_vcons(cgm_vec3 *v, float x, float y, float z);
+static CGM_INLINE cgm_vec3 cgm_vvec(float x, float y, float z);
+
+static CGM_INLINE void cgm_vadd(cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE void cgm_vadd_scaled(cgm_vec3 *a, const cgm_vec3 *b, float s); /* a+b*s */
+static CGM_INLINE void cgm_vsub(cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE void cgm_vsub_scaled(cgm_vec3 *a, const cgm_vec3 *b, float s); /* a-b*s */
+static CGM_INLINE void cgm_vmul(cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE void cgm_vscale(cgm_vec3 *v, float s);
+static CGM_INLINE void cgm_vmul_m4v3(cgm_vec3 *v, const float *m); /* m4x4 * v */
+static CGM_INLINE void cgm_vmul_v3m4(cgm_vec3 *v, const float *m); /* v * m4x4 */
+static CGM_INLINE void cgm_vmul_m3v3(cgm_vec3 *v, const float *m); /* m3x3 * v (m still 16 floats) */
+static CGM_INLINE void cgm_vmul_v3m3(cgm_vec3 *v, const float *m); /* v * m3x3 (m still 16 floats) */
+
+static CGM_INLINE float cgm_vdot(const cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE void cgm_vcross(cgm_vec3 *res, const cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE float cgm_vlength(const cgm_vec3 *v);
+static CGM_INLINE float cgm_vlength_sq(const cgm_vec3 *v);
+static CGM_INLINE float cgm_vdist(const cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE float cgm_vdist_sq(const cgm_vec3 *a, const cgm_vec3 *b);
+static CGM_INLINE void cgm_vnormalize(cgm_vec3 *v);
+
+static CGM_INLINE void cgm_vreflect(cgm_vec3 *v, const cgm_vec3 *n);
+static CGM_INLINE void cgm_vrefract(cgm_vec3 *v, const cgm_vec3 *n, float ior);
+
+static CGM_INLINE void cgm_vrotate_quat(cgm_vec3 *v, const cgm_quat *q);
+static CGM_INLINE void cgm_vrotate_axis(cgm_vec3 *v, int axis, float angle);
+static CGM_INLINE void cgm_vrotate(cgm_vec3 *v, float angle, float x, float y, float z);
+static CGM_INLINE void cgm_vrotate_euler(cgm_vec3 *v, float a, float b, float c, enum cgm_euler_mode mode);
+
+static CGM_INLINE void cgm_vlerp(cgm_vec3 *res, const cgm_vec3 *a, const cgm_vec3 *b, float t);
+
+#define cgm_velem(vptr, idx) ((&(vptr)->x)[idx])
+
+/* --- operations on cgm_vec4 --- */
+static CGM_INLINE void cgm_wcons(cgm_vec4 *v, float x, float y, float z, float w);
+static CGM_INLINE cgm_vec4 cgm_wvec(float x, float y, float z, float w);
+
+static CGM_INLINE void cgm_wadd(cgm_vec4 *a, const cgm_vec4 *b);
+static CGM_INLINE void cgm_wsub(cgm_vec4 *a, const cgm_vec4 *b);
+static CGM_INLINE void cgm_wmul(cgm_vec4 *a, const cgm_vec4 *b);
+static CGM_INLINE void cgm_wscale(cgm_vec4 *v, float s);
+
+static CGM_INLINE void cgm_wmul_m4v4(cgm_vec4 *v, const float *m);
+static CGM_INLINE void cgm_wmul_v4m4(cgm_vec4 *v, const float *m);
+static CGM_INLINE void cgm_wmul_m34v4(cgm_vec4 *v, const float *m); /* doesn't affect w */
+static CGM_INLINE void cgm_wmul_v4m43(cgm_vec4 *v, const float *m); /* doesn't affect w */
+static CGM_INLINE void cgm_wmul_m3v4(cgm_vec4 *v, const float *m); /* (m still 16 floats) */
+static CGM_INLINE void cgm_wmul_v4m3(cgm_vec4 *v, const float *m); /* (m still 16 floats) */
+
+static CGM_INLINE float cgm_wdot(const cgm_vec4 *a, const cgm_vec4 *b);
+
+static CGM_INLINE float cgm_wlength(const cgm_vec4 *v);
+static CGM_INLINE float cgm_wlength_sq(const cgm_vec4 *v);
+static CGM_INLINE float cgm_wdist(const cgm_vec4 *a, const cgm_vec4 *b);
+static CGM_INLINE float cgm_wdist_sq(const cgm_vec4 *a, const cgm_vec4 *b);
+static CGM_INLINE void cgm_wnormalize(cgm_vec4 *v);
+
+static CGM_INLINE void cgm_wlerp(cgm_vec4 *res, const cgm_vec4 *a, const cgm_vec4 *b, float t);
+
+#define cgm_welem(vptr, idx) ((&(vptr)->x)[idx])
+
+/* --- operations on quaternions --- */
+static CGM_INLINE void cgm_qcons(cgm_quat *q, float x, float y, float z, float w);
+
+static CGM_INLINE void cgm_qneg(cgm_quat *q);
+static CGM_INLINE void cgm_qadd(cgm_quat *a, const cgm_quat *b);
+static CGM_INLINE void cgm_qsub(cgm_quat *a, const cgm_quat *b);
+static CGM_INLINE void cgm_qmul(cgm_quat *a, const cgm_quat *b);
+
+static CGM_INLINE float cgm_qlength(const cgm_quat *q);
+static CGM_INLINE float cgm_qlength_sq(const cgm_quat *q);
+static CGM_INLINE void cgm_qnormalize(cgm_quat *q);
+static CGM_INLINE void cgm_qconjugate(cgm_quat *q);
+static CGM_INLINE void cgm_qinvert(cgm_quat *q);
+
+static CGM_INLINE void cgm_qrotation(cgm_quat *q, float angle, float x, float y, float z);
+static CGM_INLINE void cgm_qrotate(cgm_quat *q, float angle, float x, float y, float z);
+
+static CGM_INLINE void cgm_qslerp(cgm_quat *res, const cgm_quat *a, const cgm_quat *b, float t);
+static CGM_INLINE void cgm_qlerp(cgm_quat *res, const cgm_quat *a, const cgm_quat *b, float t);
+
+#define cgm_qelem(qptr, idx) ((&(qptr)->x)[idx])
+
+/* --- operations on matrices --- */
+static CGM_INLINE void cgm_mcopy(float *dest, const float *src);
+static CGM_INLINE void cgm_mzero(float *m);
+static CGM_INLINE void cgm_midentity(float *m);
+
+static CGM_INLINE void cgm_mmul(float *a, const float *b);
+static CGM_INLINE void cgm_mpremul(float *a, const float *b);
+
+static CGM_INLINE void cgm_msubmatrix(float *m, int row, int col);
+static CGM_INLINE void cgm_mupper3(float *m);
+static CGM_INLINE float cgm_msubdet(const float *m, int row, int col);
+static CGM_INLINE float cgm_mcofactor(const float *m, int row, int col);
+static CGM_INLINE float cgm_mdet(const float *m);
+static CGM_INLINE void cgm_mtranspose(float *m);
+static CGM_INLINE void cgm_mcofmatrix(float *m);
+static CGM_INLINE int cgm_minverse(float *m); /* returns 0 on success, -1 for singular */
+
+static CGM_INLINE void cgm_mtranslation(float *m, float x, float y, float z);
+static CGM_INLINE void cgm_mscaling(float *m, float sx, float sy, float sz);
+static CGM_INLINE void cgm_mrotation_x(float *m, float angle);
+static CGM_INLINE void cgm_mrotation_y(float *m, float angle);
+static CGM_INLINE void cgm_mrotation_z(float *m, float angle);
+static CGM_INLINE void cgm_mrotation_axis(float *m, int idx, float angle);
+static CGM_INLINE void cgm_mrotation(float *m, float angle, float x, float y, float z);
+static CGM_INLINE void cgm_mrotation_euler(float *m, float a, float b, float c, int mode);
+static CGM_INLINE void cgm_mrotation_quat(float *m, const cgm_quat *q);
+
+static CGM_INLINE void cgm_mtranslate(float *m, float x, float y, float z);
+static CGM_INLINE void cgm_mscale(float *m, float sx, float sy, float sz);
+static CGM_INLINE void cgm_mrotate_x(float *m, float angle);
+static CGM_INLINE void cgm_mrotate_y(float *m, float angle);
+static CGM_INLINE void cgm_mrotate_z(float *m, float angle);
+static CGM_INLINE void cgm_mrotate_axis(float *m, int idx, float angle);
+static CGM_INLINE void cgm_mrotate(float *m, float angle, float x, float y, float z);
+static CGM_INLINE void cgm_mrotate_euler(float *m, float a, float b, float c, int mode);
+static CGM_INLINE void cgm_mrotate_quat(float *m, const cgm_quat *q);
+
+static CGM_INLINE void cgm_mpretranslate(float *m, float x, float y, float z);
+static CGM_INLINE void cgm_mprescale(float *m, float sx, float sy, float sz);
+static CGM_INLINE void cgm_mprerotate_x(float *m, float angle);
+static CGM_INLINE void cgm_mprerotate_y(float *m, float angle);
+static CGM_INLINE void cgm_mprerotate_z(float *m, float angle);
+static CGM_INLINE void cgm_mprerotate_axis(float *m, int idx, float angle);
+static CGM_INLINE void cgm_mprerotate(float *m, float angle, float x, float y, float z);
+static CGM_INLINE void cgm_mprerotate_euler(float *m, float a, float b, float c, int mode);
+static CGM_INLINE void cgm_mprerotate_quat(float *m, const cgm_quat *q);
+
+static CGM_INLINE void cgm_mget_translation(const float *m, cgm_vec3 *res);
+static CGM_INLINE void cgm_mget_rotation(const float *m, cgm_quat *res);
+static CGM_INLINE void cgm_mget_scaling(const float *m, cgm_vec3 *res);
+static CGM_INLINE void cgm_mget_frustum_plane(const float *m, int p, cgm_vec4 *res);
+
+static CGM_INLINE void cgm_normalize_plane(cgm_vec4 *p);
+
+static CGM_INLINE void cgm_mlookat(float *m, const cgm_vec3 *pos, const cgm_vec3 *targ,
+ const cgm_vec3 *up);
+static CGM_INLINE void cgm_minv_lookat(float *m, const cgm_vec3 *pos, const cgm_vec3 *targ,
+ const cgm_vec3 *up);
+static CGM_INLINE void cgm_mortho(float *m, float left, float right, float bot, float top,
+ float znear, float zfar);
+static CGM_INLINE void cgm_mfrustum(float *m, float left, float right, float bot, float top,
+ float znear, float zfar);
+static CGM_INLINE void cgm_mperspective(float *m, float vfov, float aspect, float znear, float zfar);
+
+static CGM_INLINE void cgm_mmirror(float *m, float a, float b, float c, float d);
+
+/* --- operations on rays --- */
+static CGM_INLINE void cgm_rcons(cgm_ray *r, float x, float y, float z, float dx, float dy, float dz);
+
+static CGM_INLINE void cgm_rmul_mr(cgm_ray *ray, const float *m); /* m4x4 * ray */
+static CGM_INLINE void cgm_rmul_rm(cgm_ray *ray, const float *m); /* ray * m4x4 */
+
+static CGM_INLINE void cgm_rreflect(cgm_ray *ray, const cgm_vec3 *n);
+static CGM_INLINE void cgm_rrefract(cgm_ray *ray, const cgm_vec3 *n, float ior);
+
+/* --- miscellaneous utility functions --- */
+static CGM_INLINE float cgm_deg_to_rad(float deg);
+static CGM_INLINE float cgm_rad_to_deg(float rad);
+
+static CGM_INLINE float cgm_smoothstep(float a, float b, float x);
+static CGM_INLINE float cgm_lerp(float a, float b, float t);
+static CGM_INLINE float cgm_logerp(float a, float b, float t);
+static CGM_INLINE float cgm_bezier(float a, float b, float c, float d, float t);
+static CGM_INLINE float cgm_bspline(float a, float b, float c, float d, float t);
+static CGM_INLINE float cgm_spline(float a, float b, float c, float d, float t);
+
+static CGM_INLINE void cgm_discrand(cgm_vec3 *v, float rad);
+static CGM_INLINE void cgm_sphrand(cgm_vec3 *v, float rad);
+
+static CGM_INLINE void cgm_unproject(cgm_vec3 *res, const cgm_vec3 *norm_scrpos,
+ const float *inv_viewproj);
+static CGM_INLINE void cgm_glu_unproject(float winx, float winy, float winz,
+ const float *view, const float *proj, const int *vp,
+ float *objx, float *objy, float *objz);
+
+static CGM_INLINE void cgm_pick_ray(cgm_ray *ray, float nx, float ny,
+ const float *viewmat, const float *projmat);
+
+static CGM_INLINE void cgm_raypos(cgm_vec3 *p, const cgm_ray *ray, float t);
+
+/* calculate barycentric coordinates of point pt in triangle (a, b, c) */
+static CGM_INLINE void cgm_bary(cgm_vec3 *bary, const cgm_vec3 *a,
+ const cgm_vec3 *b, const cgm_vec3 *c, const cgm_vec3 *pt);
+
+/* convert between unit vectors and spherical coordinates */
+static CGM_INLINE void cgm_uvec_to_sph(float *theta, float *phi, const cgm_vec3 *v);
+static CGM_INLINE void cgm_sph_to_uvec(cgm_vec3 *v, float theta, float phi);
+
+#include "cgmvec3.inl"
+#include "cgmvec4.inl"
+#include "cgmquat.inl"
+#include "cgmmat.inl"
+#include "cgmray.inl"
+#include "cgmmisc.inl"
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* CGMATH_H_ */
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+static CGM_INLINE void cgm_mcopy(float *dest, const float *src)
+{
+ memcpy(dest, src, 16 * sizeof(float));
+}
+
+static CGM_INLINE void cgm_mzero(float *m)
+{
+ static float z[16];
+ cgm_mcopy(m, z);
+}
+
+static CGM_INLINE void cgm_midentity(float *m)
+{
+ static float id[16] = {1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1};
+ cgm_mcopy(m, id);
+}
+
+static CGM_INLINE void cgm_mmul(float *a, const float *b)
+{
+ int i, j;
+ float res[16];
+ float *resptr = res;
+ float *arow = a;
+
+ for(i=0; i<4; i++) {
+ for(j=0; j<4; j++) {
+ *resptr++ = arow[0] * b[j] + arow[1] * b[4 + j] +
+ arow[2] * b[8 + j] + arow[3] * b[12 + j];
+ }
+ arow += 4;
+ }
+ cgm_mcopy(a, res);
+}
+
+static CGM_INLINE void cgm_mpremul(float *a, const float *b)
+{
+ int i, j;
+ float res[16];
+ float *resptr = res;
+ const float *brow = b;
+
+ for(i=0; i<4; i++) {
+ for(j=0; j<4; j++) {
+ *resptr++ = brow[0] * a[j] + brow[1] * a[4 + j] +
+ brow[2] * a[8 + j] + brow[3] * a[12 + j];
+ }
+ brow += 4;
+ }
+ cgm_mcopy(a, res);
+}
+
+static CGM_INLINE void cgm_msubmatrix(float *m, int row, int col)
+{
+ float orig[16];
+ int i, j, subi, subj;
+
+ cgm_mcopy(orig, m);
+
+ subi = 0;
+ for(i=0; i<4; i++) {
+ if(i == row) continue;
+
+ subj = 0;
+ for(j=0; j<4; j++) {
+ if(j == col) continue;
+
+ m[subi * 4 + subj++] = orig[i * 4 + j];
+ }
+ subi++;
+ }
+
+ cgm_mupper3(m);
+}
+
+static CGM_INLINE void cgm_mupper3(float *m)
+{
+ m[3] = m[7] = m[11] = m[12] = m[13] = m[14] = 0.0f;
+ m[15] = 1.0f;
+}
+
+static CGM_INLINE float cgm_msubdet(const float *m, int row, int col)
+{
+ float tmp[16];
+ float subdet00, subdet01, subdet02;
+
+ cgm_mcopy(tmp, m);
+ cgm_msubmatrix(tmp, row, col);
+
+ subdet00 = tmp[5] * tmp[10] - tmp[6] * tmp[9];
+ subdet01 = tmp[4] * tmp[10] - tmp[6] * tmp[8];
+ subdet02 = tmp[4] * tmp[9] - tmp[5] * tmp[8];
+
+ return tmp[0] * subdet00 - tmp[1] * subdet01 + tmp[2] * subdet02;
+}
+
+static CGM_INLINE float cgm_mcofactor(const float *m, int row, int col)
+{
+ float min = cgm_msubdet(m, row, col);
+ return (row + col) & 1 ? -min : min;
+}
+
+static CGM_INLINE float cgm_mdet(const float *m)
+{
+ return m[0] * cgm_msubdet(m, 0, 0) - m[1] * cgm_msubdet(m, 0, 1) +
+ m[2] * cgm_msubdet(m, 0, 2) - m[3] * cgm_msubdet(m, 0, 3);
+}
+
+static CGM_INLINE void cgm_mtranspose(float *m)
+{
+ int i, j;
+ for(i=0; i<4; i++) {
+ for(j=0; j<i; j++) {
+ int a = i * 4 + j;
+ int b = j * 4 + i;
+ float tmp = m[a];
+ m[a] = m[b];
+ m[b] = tmp;
+ }
+ }
+}
+
+static CGM_INLINE void cgm_mcofmatrix(float *m)
+{
+ float tmp[16];
+ int i, j;
+
+ cgm_mcopy(tmp, m);
+
+ for(i=0; i<4; i++) {
+ for(j=0; j<4; j++) {
+ m[i * 4 + j] = cgm_mcofactor(tmp, i, j);
+ }
+ }
+}
+
+static CGM_INLINE int cgm_minverse(float *m)
+{
+ int i, j;
+ float tmp[16];
+ float inv_det;
+ float det = cgm_mdet(m);
+ if(det == 0.0f) return -1;
+ inv_det = 1.0f / det;
+
+ cgm_mcopy(tmp, m);
+
+ for(i=0; i<4; i++) {
+ for(j=0; j<4; j++) {
+ m[i * 4 + j] = cgm_mcofactor(tmp, j, i) * inv_det; /* transposed */
+ }
+ }
+ return 0;
+}
+
+static CGM_INLINE void cgm_mtranslation(float *m, float x, float y, float z)
+{
+ cgm_midentity(m);
+ m[12] = x;
+ m[13] = y;
+ m[14] = z;
+}
+
+static CGM_INLINE void cgm_mscaling(float *m, float sx, float sy, float sz)
+{
+ cgm_mzero(m);
+ m[0] = sx;
+ m[5] = sy;
+ m[10] = sz;
+ m[15] = 1.0f;
+}
+
+static CGM_INLINE void cgm_mrotation_x(float *m, float angle)
+{
+ float sa = sin(angle);
+ float ca = cos(angle);
+
+ cgm_midentity(m);
+ m[5] = ca;
+ m[6] = sa;
+ m[9] = -sa;
+ m[10] = ca;
+}
+
+static CGM_INLINE void cgm_mrotation_y(float *m, float angle)
+{
+ float sa = sin(angle);
+ float ca = cos(angle);
+
+ cgm_midentity(m);
+ m[0] = ca;
+ m[2] = -sa;
+ m[8] = sa;
+ m[10] = ca;
+}
+
+static CGM_INLINE void cgm_mrotation_z(float *m, float angle)
+{
+ float sa = sin(angle);
+ float ca = cos(angle);
+
+ cgm_midentity(m);
+ m[0] = ca;
+ m[1] = sa;
+ m[4] = -sa;
+ m[5] = ca;
+}
+
+static CGM_INLINE void cgm_mrotation_axis(float *m, int idx, float angle)
+{
+ switch(idx) {
+ case 0:
+ cgm_mrotation_x(m, angle);
+ break;
+ case 1:
+ cgm_mrotation_y(m, angle);
+ break;
+ case 2:
+ cgm_mrotation_z(m, angle);
+ break;
+ }
+}
+
+static CGM_INLINE void cgm_mrotation(float *m, float angle, float x, float y, float z)
+{
+ float sa = sin(angle);
+ float ca = cos(angle);
+ float invca = 1.0f - ca;
+ float xsq = x * x;
+ float ysq = y * y;
+ float zsq = z * z;
+
+ cgm_mzero(m);
+ m[15] = 1.0f;
+
+ m[0] = xsq + (1.0f - xsq) * ca;
+ m[4] = x * y * invca - z * sa;
+ m[8] = x * z * invca + y * sa;
+
+ m[1] = x * y * invca + z * sa;
+ m[5] = ysq + (1.0f - ysq) * ca;
+ m[9] = y * z * invca - x * sa;
+
+ m[2] = x * z * invca - y * sa;
+ m[6] = y * z * invca + x * sa;
+ m[10] = zsq + (1.0f - zsq) * ca;
+}
+
+static CGM_INLINE void cgm_mrotation_euler(float *m, float a, float b, float c, int mode)
+{
+ /* this array must match the EulerMode enum */
+ static const int axis[][3] = {
+ {0, 1, 2}, {0, 2, 1},
+ {1, 0, 2}, {1, 2, 0},
+ {2, 0, 1}, {2, 1, 0},
+ {2, 0, 2}, {2, 1, 2},
+ {1, 0, 1}, {1, 2, 1},
+ {0, 1, 0}, {0, 2, 0}
+ };
+
+ float ma[16], mb[16];
+ cgm_mrotation_axis(ma, axis[mode][0], a);
+ cgm_mrotation_axis(mb, axis[mode][1], b);
+ cgm_mrotation_axis(m, axis[mode][2], c);
+ cgm_mmul(m, mb);
+ cgm_mmul(m, ma);
+}
+
+static CGM_INLINE void cgm_mrotation_quat(float *m, const cgm_quat *q)
+{
+ float xsq2 = 2.0f * q->x * q->x;
+ float ysq2 = 2.0f * q->y * q->y;
+ float zsq2 = 2.0f * q->z * q->z;
+ float sx = 1.0f - ysq2 - zsq2;
+ float sy = 1.0f - xsq2 - zsq2;
+ float sz = 1.0f - xsq2 - ysq2;
+
+ m[3] = m[7] = m[11] = m[12] = m[13] = m[14] = 0.0f;
+ m[15] = 1.0f;
+
+ m[0] = sx;
+ m[1] = 2.0f * q->x * q->y + 2.0f * q->w * q->z;
+ m[2] = 2.0f * q->z * q->x - 2.0f * q->w * q->y;
+ m[4] = 2.0f * q->x * q->y - 2.0f * q->w * q->z;
+ m[5] = sy;
+ m[6] = 2.0f * q->y * q->z + 2.0f * q->w * q->x;
+ m[8] = 2.0f * q->z * q->x + 2.0f * q->w * q->y;
+ m[9] = 2.0f * q->y * q->z - 2.0f * q->w * q->x;
+ m[10] = sz;
+}
+
+static CGM_INLINE void cgm_mtranslate(float *m, float x, float y, float z)
+{
+ float tm[16];
+ cgm_mtranslation(tm, x, y, z);
+ cgm_mmul(m, tm);
+}
+
+static CGM_INLINE void cgm_mscale(float *m, float sx, float sy, float sz)
+{
+ float sm[16];
+ cgm_mscaling(sm, sx, sy, sz);
+ cgm_mmul(m, sm);
+}
+
+static CGM_INLINE void cgm_mrotate_x(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_x(rm, angle);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate_y(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_y(rm, angle);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate_z(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_z(rm, angle);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate_axis(float *m, int idx, float angle)
+{
+ float rm[16];
+ cgm_mrotation_axis(rm, idx, angle);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate(float *m, float angle, float x, float y, float z)
+{
+ float rm[16];
+ cgm_mrotation(rm, angle, x, y, z);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate_euler(float *m, float a, float b, float c, int mode)
+{
+ float rm[16];
+ cgm_mrotation_euler(rm, a, b, c, mode);
+ cgm_mmul(m, rm);
+}
+
+static CGM_INLINE void cgm_mrotate_quat(float *m, const cgm_quat *q)
+{
+ float rm[16];
+ cgm_mrotation_quat(rm, q);
+ cgm_mmul(m, rm);
+}
+
+
+static CGM_INLINE void cgm_mpretranslate(float *m, float x, float y, float z)
+{
+ float tm[16];
+ cgm_mtranslation(tm, x, y, z);
+ cgm_mpremul(m, tm);
+}
+
+static CGM_INLINE void cgm_mprescale(float *m, float sx, float sy, float sz)
+{
+ float sm[16];
+ cgm_mscaling(sm, sx, sy, sz);
+ cgm_mpremul(m, sm);
+}
+
+static CGM_INLINE void cgm_mprerotate_x(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_x(rm, angle);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate_y(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_y(rm, angle);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate_z(float *m, float angle)
+{
+ float rm[16];
+ cgm_mrotation_z(rm, angle);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate_axis(float *m, int idx, float angle)
+{
+ float rm[16];
+ cgm_mrotation_axis(rm, idx, angle);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate(float *m, float angle, float x, float y, float z)
+{
+ float rm[16];
+ cgm_mrotation(rm, angle, x, y, z);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate_euler(float *m, float a, float b, float c, int mode)
+{
+ float rm[16];
+ cgm_mrotation_euler(rm, a, b, c, mode);
+ cgm_mpremul(m, rm);
+}
+
+static CGM_INLINE void cgm_mprerotate_quat(float *m, const cgm_quat *q)
+{
+ float rm[16];
+ cgm_mrotation_quat(rm, q);
+ cgm_mpremul(m, rm);
+}
+
+
+static CGM_INLINE void cgm_mget_translation(const float *m, cgm_vec3 *res)
+{
+ res->x = m[12];
+ res->y = m[13];
+ res->z = m[14];
+}
+
+/* Algorithm in Ken Shoemake's article in 1987 SIGGRAPH course notes
+ * article "Quaternion Calculus and Fast Animation".
+ * adapted from: http://www.geometrictools.com/LibMathematics/Algebra/Wm5Quaternion.inl
+ */
+static CGM_INLINE void cgm_mget_rotation(const float *m, cgm_quat *res)
+{
+ static const int next[3] = {1, 2, 0};
+ float quat[4];
+ int i, j, k;
+
+ float trace = m[0] + m[5] + m[10];
+ float root;
+
+ if(trace > 0.0f) {
+ /* |w| > 1/2 */
+ root = sqrt(trace + 1.0f); /* 2w */
+ res->w = 0.5f * root;
+ root = 0.5f / root; /* 1 / 4w */
+ res->x = (m[6] - m[9]) * root;
+ res->y = (m[8] - m[2]) * root;
+ res->z = (m[1] - m[4]) * root;
+ } else {
+ /* |w| <= 1/2 */
+ i = 0;
+ if(m[5] > m[0]) {
+ i = 1;
+ }
+ if(m[10] > m[i * 4 + i]) {
+ i = 2;
+ }
+ j = next[i];
+ k = next[j];
+
+ root = sqrt(m[i * 4 + i] - m[j * 4 + j] - m[k * 4 + k] + 1.0f);
+ quat[i + 1] = 0.5f * root;
+ root = 0.5f / root;
+ quat[0] = (m[j + 4 + k] - m[k * 4 + j]) * root;
+ quat[j + 1] = (m[i * 4 + j] - m[j * 4 + i]) * root;
+ quat[k + 1] = (m[i * 4 + k] - m[k * 4 + i]) * root;
+ res->w = quat[0];
+ res->x = quat[1];
+ res->y = quat[2];
+ res->z = quat[3];
+ }
+}
+
+static CGM_INLINE void cgm_mget_scaling(const float *m, cgm_vec3 *res)
+{
+ res->x = sqrt(m[0] * m[0] + m[4] * m[4] + m[8] * m[8]);
+ res->y = sqrt(m[1] * m[1] + m[5] * m[5] + m[9] * m[9]);
+ res->z = sqrt(m[2] * m[2] + m[6] * m[6] + m[10] * m[10]);
+}
+
+static CGM_INLINE void cgm_mget_frustum_plane(const float *m, int p, cgm_vec4 *res)
+{
+ switch(p) {
+ case 0:
+ res->x = m[3] + m[0];
+ res->y = m[7] + m[4];
+ res->z = m[11] + m[8];
+ res->w = m[15] + m[12];
+ break;
+
+ case 1:
+ res->x = m[3] - m[0];
+ res->y = m[7] - m[4];
+ res->z = m[11] - m[8];
+ res->w = m[15] - m[12];
+ break;
+
+ case 2:
+ res->x = m[3] + m[1];
+ res->y = m[7] + m[5];
+ res->z = m[11] + m[9];
+ res->w = m[15] + m[13];
+ break;
+
+ case 3:
+ res->x = m[3] - m[1];
+ res->y = m[7] - m[5];
+ res->z = m[11] - m[9];
+ res->w = m[15] - m[13];
+ break;
+
+ case 4:
+ res->x = m[3] + m[2];
+ res->y = m[7] + m[6];
+ res->z = m[11] + m[10];
+ res->w = m[15] + m[14];
+ break;
+
+ case 5:
+ res->x = m[3] - m[2];
+ res->y = m[7] - m[6];
+ res->z = m[11] - m[10];
+ res->w = m[15] - m[14];
+ break;
+
+ default:
+ break;
+ }
+}
+
+static CGM_INLINE void cgm_normalize_plane(cgm_vec4 *p)
+{
+ float len = cgm_vlength((cgm_vec3*)p);
+ if(len != 0.0f) {
+ float s = 1.0f / len;
+ p->x *= s;
+ p->y *= s;
+ p->z *= s;
+ p->w *= s;
+ }
+}
+
+static CGM_INLINE void cgm_mlookat(float *m, const cgm_vec3 *pos, const cgm_vec3 *targ,
+ const cgm_vec3 *up)
+{
+ float trans[16];
+ cgm_vec3 dir = *targ, right, vup;
+
+ cgm_vsub(&dir, pos);
+ cgm_vnormalize(&dir);
+ cgm_vcross(&right, &dir, up);
+ cgm_vnormalize(&right);
+ cgm_vcross(&vup, &right, &dir);
+ cgm_vnormalize(&vup);
+
+ cgm_midentity(m);
+ m[0] = right.x;
+ m[1] = right.y;
+ m[2] = right.z;
+ m[4] = vup.x;
+ m[5] = vup.y;
+ m[6] = vup.z;
+ m[8] = -dir.x;
+ m[9] = -dir.y;
+ m[10] = -dir.z;
+
+ cgm_mtranslation(trans, pos->x, pos->y, pos->z);
+ cgm_mmul(m, trans);
+}
+
+static CGM_INLINE void cgm_minv_lookat(float *m, const cgm_vec3 *pos, const cgm_vec3 *targ,
+ const cgm_vec3 *up)
+{
+ float rot[16];
+ cgm_vec3 dir = *targ, right, vup;
+
+ cgm_vsub(&dir, pos);
+ cgm_vnormalize(&dir);
+ cgm_vcross(&right, &dir, up);
+ cgm_vnormalize(&right);
+ cgm_vcross(&vup, &right, &dir);
+ cgm_vnormalize(&vup);
+
+ cgm_midentity(rot);
+ rot[0] = right.x;
+ rot[4] = right.y;
+ rot[8] = right.z;
+ rot[1] = vup.x;
+ rot[5] = vup.y;
+ rot[9] = vup.z;
+ rot[2] = -dir.x;
+ rot[6] = -dir.y;
+ rot[10] = -dir.z;
+
+ cgm_mtranslation(m, -pos->x, -pos->y, -pos->z);
+ cgm_mmul(m, rot);
+}
+
+static CGM_INLINE void cgm_mortho(float *m, float left, float right, float bot, float top,
+ float znear, float zfar)
+{
+ float dx = right - left;
+ float dy = top - bot;
+ float dz = zfar - znear;
+
+ cgm_midentity(m);
+ m[0] = 2.0f / dx;
+ m[5] = 2.0f / dy;
+ m[10] = -2.0f / dz;
+ m[12] = -(right + left) / dx;
+ m[13] = -(top + bot) / dy;
+ m[14] = -(zfar + znear) / dz;
+}
+
+static CGM_INLINE void cgm_mfrustum(float *m, float left, float right, float bot, float top,
+ float znear, float zfar)
+{
+ float dx = right - left;
+ float dy = top - bot;
+ float dz = zfar - znear;
+
+ cgm_mzero(m);
+ m[0] = 2.0f * znear / dx;
+ m[5] = 2.0f * znear / dy;
+ m[8] = (right + left) / dx;
+ m[9] = (top + bot) / dy;
+ m[10] = -(zfar + znear) / dz;
+ m[14] = -2.0f * zfar * znear / dz;
+ m[11] = -1.0f;
+}
+
+static CGM_INLINE void cgm_mperspective(float *m, float vfov, float aspect, float znear, float zfar)
+{
+ float s = 1.0f / (float)tan(vfov / 2.0f);
+ float range = znear - zfar;
+
+ cgm_mzero(m);
+ m[0] = s / aspect;
+ m[5] = s;
+ m[10] = (znear + zfar) / range;
+ m[14] = 2.0f * znear * zfar / range;
+ m[11] = -1.0f;
+}
+
+static CGM_INLINE void cgm_mmirror(float *m, float a, float b, float c, float d)
+{
+ m[0] = 1.0f - 2.0f * a * a;
+ m[5] = 1.0f - 2.0f * b * b;
+ m[10] = 1.0f - 2.0f * c * c;
+ m[15] = 1.0f;
+
+ m[1] = m[4] = -2.0f * a * b;
+ m[2] = m[8] = -2.0f * a * c;
+ m[6] = m[9] = -2.0f * b * c;
+
+ m[12] = -2.0f * a * d;
+ m[13] = -2.0f * b * d;
+ m[14] = -2.0f * c * d;
+
+ m[3] = m[7] = m[11] = 0.0f;
+}
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+#include <stdlib.h>
+
+static CGM_INLINE float cgm_deg_to_rad(float deg)
+{
+ return M_PI * deg / 180.0f;
+}
+
+static CGM_INLINE float cgm_rad_to_deg(float rad)
+{
+ return 180.0f * rad / M_PI;
+}
+
+static CGM_INLINE float cgm_smoothstep(float a, float b, float x)
+{
+ if(x < a) return 0.0f;
+ if(x >= b) return 1.0f;
+
+ x = (x - a) / (b - a);
+ return x * x * (3.0f - 2.0f * x);
+}
+
+static CGM_INLINE float cgm_lerp(float a, float b, float t)
+{
+ return a + (b - a) * t;
+}
+
+static CGM_INLINE float cgm_logerp(float a, float b, float t)
+{
+ if(a == 0.0f) return 0.0f;
+ return a * pow(b / a, t);
+}
+
+static CGM_INLINE float cgm_bezier(float a, float b, float c, float d, float t)
+{
+ float omt, omt3, t3, f;
+ t3 = t * t * t;
+ omt = 1.0f - t;
+ omt3 = omt * omt * omt;
+ f = 3.0f * t * omt;
+
+ return (a * omt3) + (b * f * omt) + (c * f * t) + (d * t3);
+}
+
+static CGM_INLINE float cgm_bspline(float a, float b, float c, float d, float t)
+{
+ static const float mat[] = {
+ -1, 3, -3, 1,
+ 3, -6, 0, 4,
+ -3, 3, 3, 1,
+ 1, 0, 0, 0
+ };
+ cgm_vec4 tmp, qfact;
+ float tsq = t * t;
+
+ cgm_wcons(&qfact, tsq * t, tsq, t, 1.0f);
+ cgm_wcons(&tmp, a, b, c, d);
+ cgm_wmul_m4v4(&tmp, mat);
+ cgm_wscale(&tmp, 1.0f / 6.0f);
+ return cgm_wdot(&tmp, &qfact);
+}
+
+static CGM_INLINE float cgm_spline(float a, float b, float c, float d, float t)
+{
+ static const float mat[] = {
+ -1, 2, -1, 0,
+ 3, -5, 0, 2,
+ -3, 4, 1, 0,
+ 1, -1, 0, 0
+ };
+ cgm_vec4 tmp, qfact;
+ float tsq = t * t;
+
+ cgm_wcons(&qfact, tsq * t, tsq, t, 1.0f);
+ cgm_wcons(&tmp, a, b, c, d);
+ cgm_wmul_m4v4(&tmp, mat);
+ cgm_wscale(&tmp, 1.0f / 6.0f);
+ return cgm_wdot(&tmp, &qfact);
+}
+
+static CGM_INLINE void cgm_discrand(cgm_vec3 *pt, float rad)
+{
+ float theta = 2.0f * M_PI * (float)rand() / RAND_MAX;
+ float r = sqrt((float)rand() / RAND_MAX) * rad;
+ pt->x = cos(theta) * r;
+ pt->y = sin(theta) * r;
+ pt->z = 0.0f;
+}
+
+static CGM_INLINE void cgm_sphrand(cgm_vec3 *pt, float rad)
+{
+ float u, v, theta, phi;
+
+ u = (float)rand() / RAND_MAX;
+ v = (float)rand() / RAND_MAX;
+
+ theta = 2.0f * M_PI * u;
+ phi = acos(2.0f * v - 1.0f);
+
+ pt->x = cos(theta) * sin(phi) * rad;
+ pt->y = sin(theta) * sin(phi) * rad;
+ pt->z = cos(phi) * rad;
+}
+
+static CGM_INLINE void cgm_unproject(cgm_vec3 *res, const cgm_vec3 *norm_scrpos,
+ const float *inv_viewproj)
+{
+ cgm_vec4 pos;
+
+ pos.x = 2.0f * norm_scrpos->x - 1.0f;
+ pos.y = 2.0f * norm_scrpos->y - 1.0f;
+ pos.z = 2.0f * norm_scrpos->z - 1.0f;
+ pos.w = 1.0f;
+
+ cgm_wmul_m4v4(&pos, inv_viewproj);
+
+ res->x = pos.x / pos.w;
+ res->y = pos.y / pos.w;
+ res->z = pos.z / pos.w;
+}
+
+static CGM_INLINE void cgm_glu_unproject(float winx, float winy, float winz,
+ const float *view, const float *proj, const int *vp,
+ float *objx, float *objy, float *objz)
+{
+ cgm_vec3 npos, res;
+ float inv_pv[16];
+
+ cgm_mcopy(inv_pv, view);
+ cgm_mmul(inv_pv, proj);
+ cgm_minverse(inv_pv);
+
+ npos.x = (winx - vp[0]) / vp[2];
+ npos.y = (winy - vp[1]) / vp[4];
+ npos.z = winz;
+
+ cgm_unproject(&res, &npos, inv_pv);
+
+ *objx = res.x;
+ *objy = res.y;
+ *objz = res.z;
+}
+
+static CGM_INLINE void cgm_pick_ray(cgm_ray *ray, float nx, float ny,
+ const float *viewmat, const float *projmat)
+{
+ cgm_vec3 npos, farpt;
+ float inv_pv[16];
+
+ cgm_mcopy(inv_pv, viewmat);
+ cgm_mmul(inv_pv, projmat);
+ cgm_minverse(inv_pv);
+
+ cgm_vcons(&npos, nx, ny, 0.0f);
+ cgm_unproject(&ray->origin, &npos, inv_pv);
+ npos.z = 1.0f;
+ cgm_unproject(&farpt, &npos, inv_pv);
+
+ ray->dir.x = farpt.x - ray->origin.x;
+ ray->dir.y = farpt.y - ray->origin.y;
+ ray->dir.z = farpt.z - ray->origin.z;
+}
+
+static CGM_INLINE void cgm_raypos(cgm_vec3 *p, const cgm_ray *ray, float t)
+{
+ p->x = ray->origin.x + ray->dir.x * t;
+ p->y = ray->origin.y + ray->dir.y * t;
+ p->z = ray->origin.z + ray->dir.z * t;
+}
+
+static CGM_INLINE void cgm_bary(cgm_vec3 *bary, const cgm_vec3 *a,
+ const cgm_vec3 *b, const cgm_vec3 *c, const cgm_vec3 *pt)
+{
+ float d00, d01, d11, d20, d21, denom;
+ cgm_vec3 v0 = *b, v1 = *c, v2 = *pt;
+
+ cgm_vsub(&v0, a);
+ cgm_vsub(&v1, a);
+ cgm_vsub(&v2, a);
+
+ d00 = cgm_vdot(&v0, &v0);
+ d01 = cgm_vdot(&v0, &v1);
+ d11 = cgm_vdot(&v1, &v1);
+ d20 = cgm_vdot(&v2, &v0);
+ d21 = cgm_vdot(&v2, &v1);
+ denom = d00 * d11 - d01 * d01;
+
+ bary->y = (d11 * d20 - d01 * d21) / denom;
+ bary->z = (d00 * d21 - d01 * d20) / denom;
+ bary->x = 1.0f - bary->y - bary->z;
+}
+
+static CGM_INLINE void cgm_uvec_to_sph(float *theta, float *phi, const cgm_vec3 *v)
+{
+ *theta = atan2(v->z, v->x);
+ *phi = acos(v->y);
+}
+
+static CGM_INLINE void cgm_sph_to_uvec(cgm_vec3 *v, float theta, float phi)
+{
+ v->x = sin(theta) * cos(phi);
+ v->y = sin(phi);
+ v->z = cos(theta) * cos(phi);
+}
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+static CGM_INLINE void cgm_qcons(cgm_quat *q, float x, float y, float z, float w)
+{
+ q->x = x;
+ q->y = y;
+ q->z = z;
+ q->w = w;
+}
+
+
+static CGM_INLINE void cgm_qneg(cgm_quat *q)
+{
+ q->x = -q->x;
+ q->y = -q->y;
+ q->z = -q->z;
+ q->w = -q->w;
+}
+
+static CGM_INLINE void cgm_qadd(cgm_quat *a, const cgm_quat *b)
+{
+ a->x += b->x;
+ a->y += b->y;
+ a->z += b->z;
+ a->w += b->w;
+}
+
+static CGM_INLINE void cgm_qsub(cgm_quat *a, const cgm_quat *b)
+{
+ a->x -= b->x;
+ a->y -= b->y;
+ a->z -= b->z;
+ a->w -= b->w;
+}
+
+static CGM_INLINE void cgm_qmul(cgm_quat *a, const cgm_quat *b)
+{
+ float x, y, z, dot;
+ cgm_vec3 cross;
+
+ dot = a->x * b->x + a->y * b->y + a->z * b->z;
+ cgm_vcross(&cross, (cgm_vec3*)a, (cgm_vec3*)b);
+
+ x = a->w * b->x + b->w * a->x + cross.x;
+ y = a->w * b->y + b->w * a->y + cross.y;
+ z = a->w * b->z + b->w * a->z + cross.z;
+ a->w = a->w * b->w - dot;
+ a->x = x;
+ a->y = y;
+ a->z = z;
+}
+
+static CGM_INLINE float cgm_qlength(const cgm_quat *q)
+{
+ return sqrt(q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w);
+}
+
+static CGM_INLINE float cgm_qlength_sq(const cgm_quat *q)
+{
+ return q->x * q->x + q->y * q->y + q->z * q->z + q->w * q->w;
+}
+
+static CGM_INLINE void cgm_qnormalize(cgm_quat *q)
+{
+ float len = cgm_qlength(q);
+ if(len != 0.0f) {
+ float s = 1.0f / len;
+ q->x *= s;
+ q->y *= s;
+ q->z *= s;
+ q->w *= s;
+ }
+}
+
+static CGM_INLINE void cgm_qconjugate(cgm_quat *q)
+{
+ q->x = -q->x;
+ q->y = -q->y;
+ q->z = -q->z;
+}
+
+static CGM_INLINE void cgm_qinvert(cgm_quat *q)
+{
+ float len_sq = cgm_qlength_sq(q);
+ cgm_qconjugate(q);
+ if(len_sq != 0.0f) {
+ float s = 1.0f / len_sq;
+ q->x *= s;
+ q->y *= s;
+ q->z *= s;
+ q->w *= s;
+ }
+}
+
+static CGM_INLINE void cgm_qrotation(cgm_quat *q, float angle, float x, float y, float z)
+{
+ float hangle = angle * 0.5f;
+ float sin_ha = sin(hangle);
+ q->w = cos(hangle);
+ q->x = x * sin_ha;
+ q->y = y * sin_ha;
+ q->z = z * sin_ha;
+}
+
+static CGM_INLINE void cgm_qrotate(cgm_quat *q, float angle, float x, float y, float z)
+{
+ cgm_quat qrot;
+ cgm_qrotation(&qrot, angle, x, y, z);
+ cgm_qmul(q, &qrot);
+}
+
+static CGM_INLINE void cgm_qslerp(cgm_quat *res, const cgm_quat *quat1, const cgm_quat *q2, float t)
+{
+ float angle, dot, a, b, sin_angle;
+ cgm_quat q1 = *quat1;
+
+ dot = quat1->x * q2->x + quat1->y * q2->y + quat1->z * q2->z + quat1->w * q2->w;
+ if(dot < 0.0f) {
+ /* make sure we inteprolate across the shortest arc */
+ cgm_qneg(&q1);
+ dot = -dot;
+ }
+
+ /* clamp dot to [-1, 1] in order to avoid domain errors in acos due to
+ * floating point imprecisions
+ */
+ if(dot < -1.0f) dot = -1.0f;
+ if(dot > 1.0f) dot = 1.0f;
+ angle = acos(dot);
+
+ sin_angle = sin(angle);
+ if(sin_angle == 0.0f) {
+ /* use linear interpolation to avoid div/zero */
+ a = 1.0f;
+ b = t;
+ } else {
+ a = sin((1.0f - t) * angle) / sin_angle;
+ b = sin(t * angle) / sin_angle;
+ }
+
+ res->x = q1.x * a + q2->x * b;
+ res->y = q1.y * a + q2->y * b;
+ res->z = q1.z * a + q2->z * b;
+ res->w = q1.w * a + q2->w * b;
+}
+
+static CGM_INLINE void cgm_qlerp(cgm_quat *res, const cgm_quat *a, const cgm_quat *b, float t)
+{
+ res->x = a->x + (b->x - a->x) * t;
+ res->y = a->y + (b->y - a->y) * t;
+ res->z = a->z + (b->z - a->z) * t;
+ res->w = a->w + (b->w - a->w) * t;
+}
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+static CGM_INLINE void cgm_rcons(cgm_ray *r, float x, float y, float z, float dx, float dy, float dz)
+{
+ r->origin.x = x;
+ r->origin.y = y;
+ r->origin.z = z;
+ r->dir.x = dx;
+ r->dir.y = dy;
+ r->dir.z = dz;
+}
+
+static CGM_INLINE void cgm_rmul_mr(cgm_ray *ray, const float *m)
+{
+ cgm_vmul_m4v3(&ray->origin, m);
+ cgm_vmul_m3v3(&ray->dir, m);
+}
+
+static CGM_INLINE void cgm_rmul_rm(cgm_ray *ray, const float *m)
+{
+ cgm_vmul_v3m4(&ray->origin, m);
+ cgm_vmul_v3m3(&ray->dir, m);
+}
+
+static CGM_INLINE void cgm_rreflect(cgm_ray *ray, const cgm_vec3 *n)
+{
+ cgm_vreflect(&ray->dir, n);
+}
+
+static CGM_INLINE void cgm_rrefract(cgm_ray *ray, const cgm_vec3 *n, float ior)
+{
+ cgm_vrefract(&ray->dir, n, ior);
+}
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+static CGM_INLINE void cgm_vcons(cgm_vec3 *v, float x, float y, float z)
+{
+ v->x = x;
+ v->y = y;
+ v->z = z;
+}
+
+static CGM_INLINE cgm_vec3 cgm_vvec(float x, float y, float z)
+{
+ cgm_vec3 v;
+ v.x = x;
+ v.y = y;
+ v.z = z;
+ return v;
+}
+
+static CGM_INLINE void cgm_vadd(cgm_vec3 *a, const cgm_vec3 *b)
+{
+ a->x += b->x;
+ a->y += b->y;
+ a->z += b->z;
+}
+
+static CGM_INLINE void cgm_vadd_scaled(cgm_vec3 *a, const cgm_vec3 *b, float s)
+{
+ a->x += b->x * s;
+ a->y += b->y * s;
+ a->z += b->z * s;
+}
+
+static CGM_INLINE void cgm_vsub(cgm_vec3 *a, const cgm_vec3 *b)
+{
+ a->x -= b->x;
+ a->y -= b->y;
+ a->z -= b->z;
+}
+
+static CGM_INLINE void cgm_vsub_scaled(cgm_vec3 *a, const cgm_vec3 *b, float s)
+{
+ a->x -= b->x * s;
+ a->y -= b->y * s;
+ a->z -= b->z * s;
+}
+
+static CGM_INLINE void cgm_vmul(cgm_vec3 *a, const cgm_vec3 *b)
+{
+ a->x *= b->x;
+ a->y *= b->y;
+ a->z *= b->z;
+}
+
+static CGM_INLINE void cgm_vscale(cgm_vec3 *v, float s)
+{
+ v->x *= s;
+ v->y *= s;
+ v->z *= s;
+}
+
+static CGM_INLINE void cgm_vmul_m4v3(cgm_vec3 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[4] + v->z * m[8] + m[12];
+ float y = v->x * m[1] + v->y * m[5] + v->z * m[9] + m[13];
+ v->z = v->x * m[2] + v->y * m[6] + v->z * m[10] + m[14];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_vmul_v3m4(cgm_vec3 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[1] + v->z * m[2] + m[3];
+ float y = v->x * m[4] + v->y * m[5] + v->z * m[6] + m[7];
+ v->z = v->x * m[8] + v->y * m[9] + v->z * m[10] + m[11];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_vmul_m3v3(cgm_vec3 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[4] + v->z * m[8];
+ float y = v->x * m[1] + v->y * m[5] + v->z * m[9];
+ v->z = v->x * m[2] + v->y * m[6] + v->z * m[10];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_vmul_v3m3(cgm_vec3 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[1] + v->z * m[2];
+ float y = v->x * m[4] + v->y * m[5] + v->z * m[6];
+ v->z = v->x * m[8] + v->y * m[9] + v->z * m[10];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE float cgm_vdot(const cgm_vec3 *a, const cgm_vec3 *b)
+{
+ return a->x * b->x + a->y * b->y + a->z * b->z;
+}
+
+static CGM_INLINE void cgm_vcross(cgm_vec3 *res, const cgm_vec3 *a, const cgm_vec3 *b)
+{
+ res->x = a->y * b->z - a->z * b->y;
+ res->y = a->z * b->x - a->x * b->z;
+ res->z = a->x * b->y - a->y * b->x;
+}
+
+static CGM_INLINE float cgm_vlength(const cgm_vec3 *v)
+{
+ return sqrt(v->x * v->x + v->y * v->y + v->z * v->z);
+}
+
+static CGM_INLINE float cgm_vlength_sq(const cgm_vec3 *v)
+{
+ return v->x * v->x + v->y * v->y + v->z * v->z;
+}
+
+static CGM_INLINE float cgm_vdist(const cgm_vec3 *a, const cgm_vec3 *b)
+{
+ float dx = a->x - b->x;
+ float dy = a->y - b->y;
+ float dz = a->z - b->z;
+ return sqrt(dx * dx + dy * dy + dz * dz);
+}
+
+static CGM_INLINE float cgm_vdist_sq(const cgm_vec3 *a, const cgm_vec3 *b)
+{
+ float dx = a->x - b->x;
+ float dy = a->y - b->y;
+ float dz = a->z - b->z;
+ return dx * dx + dy * dy + dz * dz;
+}
+
+static CGM_INLINE void cgm_vnormalize(cgm_vec3 *v)
+{
+ float len = cgm_vlength(v);
+ if(len != 0.0f) {
+ float s = 1.0f / len;
+ v->x *= s;
+ v->y *= s;
+ v->z *= s;
+ }
+}
+
+static CGM_INLINE void cgm_vreflect(cgm_vec3 *v, const cgm_vec3 *n)
+{
+ float ndotv2 = cgm_vdot(v, n) * 2.0f;
+ v->x -= n->x * ndotv2;
+ v->y -= n->y * ndotv2;
+ v->z -= n->z * ndotv2;
+}
+
+static CGM_INLINE void cgm_vrefract(cgm_vec3 *v, const cgm_vec3 *n, float ior)
+{
+ float ndotv = cgm_vdot(v, n);
+ float k = 1.0f - ior * ior * (1.0f - ndotv * ndotv);
+ if(k < 0.0f) {
+ cgm_vreflect(v, n); /* TIR */
+ } else {
+ float sqrt_k = sqrt(k);
+ v->x = ior * v->x - (ior * ndotv + sqrt_k) * n->x;
+ v->y = ior * v->y - (ior * ndotv + sqrt_k) * n->y;
+ v->z = ior * v->z - (ior * ndotv + sqrt_k) * n->z;
+ }
+}
+
+static CGM_INLINE void cgm_vrotate_quat(cgm_vec3 *v, const cgm_quat *q)
+{
+ cgm_quat vq, inv_q = *q, tmp_q = *q;
+
+ cgm_qcons(&vq, v->x, v->y, v->z, 0.0f);
+ cgm_qinvert(&inv_q);
+ cgm_qmul(&tmp_q, &vq);
+ cgm_qmul(&tmp_q, &inv_q);
+ cgm_vcons(v, tmp_q.x, tmp_q.y, tmp_q.z);
+}
+
+static CGM_INLINE void cgm_vrotate_axis(cgm_vec3 *v, int axis, float angle)
+{
+ float m[16];
+ cgm_mrotation_axis(m, axis, angle);
+ cgm_vmul_m3v3(v, m);
+}
+
+static CGM_INLINE void cgm_vrotate(cgm_vec3 *v, float angle, float x, float y, float z)
+{
+ float m[16];
+ cgm_mrotation(m, angle, x, y, z);
+ cgm_vmul_m3v3(v, m);
+}
+
+static CGM_INLINE void cgm_vrotate_euler(cgm_vec3 *v, float a, float b, float c, enum cgm_euler_mode mode)
+{
+ float m[16];
+ cgm_mrotation_euler(m, a, b, c, mode);
+ cgm_vmul_m3v3(v, m);
+}
+
+static CGM_INLINE void cgm_vlerp(cgm_vec3 *res, const cgm_vec3 *a, const cgm_vec3 *b, float t)
+{
+ res->x = a->x + (b->x - a->x) * t;
+ res->y = a->y + (b->y - a->y) * t;
+ res->z = a->z + (b->z - a->z) * t;
+}
--- /dev/null
+/* gph-cmath - C graphics math library
+ * Copyright (C) 2018-2023 John Tsiombikas <nuclear@member.fsf.org>
+ *
+ * This program is free software. Feel free to use, modify, and/or redistribute
+ * it under the terms of the MIT/X11 license. See LICENSE for details.
+ * If you intend to redistribute parts of the code without the LICENSE file
+ * replace this paragraph with the full contents of the LICENSE file.
+ */
+static CGM_INLINE void cgm_wcons(cgm_vec4 *v, float x, float y, float z, float w)
+{
+ v->x = x;
+ v->y = y;
+ v->z = z;
+ v->w = w;
+}
+
+static CGM_INLINE cgm_vec4 cgm_wvec(float x, float y, float z, float w)
+{
+ cgm_vec4 v;
+ v.x = x;
+ v.y = y;
+ v.z = z;
+ v.w = w;
+ return v;
+}
+
+static CGM_INLINE void cgm_wadd(cgm_vec4 *a, const cgm_vec4 *b)
+{
+ a->x += b->x;
+ a->y += b->y;
+ a->z += b->z;
+ a->w += b->w;
+}
+
+static CGM_INLINE void cgm_wsub(cgm_vec4 *a, const cgm_vec4 *b)
+{
+ a->x -= b->x;
+ a->y -= b->y;
+ a->z -= b->z;
+ a->w -= b->w;
+}
+
+static CGM_INLINE void cgm_wmul(cgm_vec4 *a, const cgm_vec4 *b)
+{
+ a->x *= b->x;
+ a->y *= b->y;
+ a->z *= b->z;
+ a->w *= b->w;
+}
+
+static CGM_INLINE void cgm_wscale(cgm_vec4 *v, float s)
+{
+ v->x *= s;
+ v->y *= s;
+ v->z *= s;
+ v->w *= s;
+}
+
+static CGM_INLINE void cgm_wmul_m4v4(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[4] + v->z * m[8] + v->w * m[12];
+ float y = v->x * m[1] + v->y * m[5] + v->z * m[9] + v->w * m[13];
+ float z = v->x * m[2] + v->y * m[6] + v->z * m[10] + v->w * m[14];
+ v->w = v->x * m[3] + v->y * m[7] + v->z * m[11] + v->w * m[15];
+ v->x = x;
+ v->y = y;
+ v->z = z;
+}
+
+static CGM_INLINE void cgm_wmul_v4m4(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[1] + v->z * m[2] + v->w * m[3];
+ float y = v->x * m[4] + v->y * m[5] + v->z * m[6] + v->w * m[7];
+ float z = v->x * m[8] + v->y * m[9] + v->z * m[10] + v->w * m[11];
+ v->w = v->x * m[12] + v->y * m[13] + v->z * m[14] + v->w * m[15];
+ v->x = x;
+ v->y = y;
+ v->z = z;
+}
+
+static CGM_INLINE void cgm_wmul_m34v4(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[4] + v->z * m[8] + v->w * m[12];
+ float y = v->x * m[1] + v->y * m[5] + v->z * m[9] + v->w * m[13];
+ v->z = v->x * m[2] + v->y * m[6] + v->z * m[10] + v->w * m[14];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_wmul_v4m43(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[1] + v->z * m[2] + v->w * m[3];
+ float y = v->x * m[4] + v->y * m[5] + v->z * m[6] + v->w * m[7];
+ v->z = v->x * m[8] + v->y * m[9] + v->z * m[10] + v->w * m[11];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_wmul_m3v4(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[4] + v->z * m[8];
+ float y = v->x * m[1] + v->y * m[5] + v->z * m[9];
+ v->z = v->x * m[2] + v->y * m[6] + v->z * m[10];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE void cgm_wmul_v4m3(cgm_vec4 *v, const float *m)
+{
+ float x = v->x * m[0] + v->y * m[1] + v->z * m[2];
+ float y = v->x * m[4] + v->y * m[5] + v->z * m[6];
+ v->z = v->x * m[8] + v->y * m[9] + v->z * m[10];
+ v->x = x;
+ v->y = y;
+}
+
+static CGM_INLINE float cgm_wdot(const cgm_vec4 *a, const cgm_vec4 *b)
+{
+ return a->x * b->x + a->y * b->y + a->z * b->z + a->w * b->w;
+}
+
+static CGM_INLINE float cgm_wlength(const cgm_vec4 *v)
+{
+ return sqrt(v->x * v->x + v->y * v->y + v->z * v->z + v->w * v->w);
+}
+
+static CGM_INLINE float cgm_wlength_sq(const cgm_vec4 *v)
+{
+ return v->x * v->x + v->y * v->y + v->z * v->z + v->w * v->w;
+}
+
+static CGM_INLINE float cgm_wdist(const cgm_vec4 *a, const cgm_vec4 *b)
+{
+ float dx = a->x - b->x;
+ float dy = a->y - b->y;
+ float dz = a->z - b->z;
+ float dw = a->w - b->w;
+ return sqrt(dx * dx + dy * dy + dz * dz + dw * dw);
+}
+
+static CGM_INLINE float cgm_wdist_sq(const cgm_vec4 *a, const cgm_vec4 *b)
+{
+ float dx = a->x - b->x;
+ float dy = a->y - b->y;
+ float dz = a->z - b->z;
+ float dw = a->w - b->w;
+ return dx * dx + dy * dy + dz * dz + dw * dw;
+}
+
+static CGM_INLINE void cgm_wnormalize(cgm_vec4 *v)
+{
+ float len = cgm_wlength(v);
+ if(len != 0.0f) {
+ float s = 1.0f / len;
+ v->x *= s;
+ v->y *= s;
+ v->z *= s;
+ v->w *= s;
+ }
+}
+
+static CGM_INLINE void cgm_wlerp(cgm_vec4 *res, const cgm_vec4 *a, const cgm_vec4 *b, float t)
+{
+ res->x = a->x + (b->x - a->x) * t;
+ res->y = a->y + (b->y - a->y) * t;
+ res->z = a->z + (b->z - a->z) * t;
+ res->w = a->w + (b->w - a->w) * t;
+}
--- /dev/null
+#include "geom.h"
+
+
+int ray_sphere(cgm_ray *ray, union gobject *obj, struct rayhit *hit)
+{
+}
+
+int ray_plane(cgm_ray *ray, union gobject *obj, struct rayhit *hit);
+{
+}
+
+void texlookup(cgm_vec3 *col, struct texture *tex, float u, float v)
+{
+}
+
--- /dev/null
+#ifndef GEOM_H_
+#define GEOM_H_
+
+struct texture {
+ int width, height;
+ int xshift;
+ uint32_t *pixels;
+};
+
+struct material {
+ cgm_vec3 kd, ks;
+ float shininess;
+ struct texture *tex;
+};
+
+#define GEOM_COMMON \
+ int type; \
+ float xform[16], inv_xform[16]; \
+ struct material *mtl
+
+struct gsphere {
+ GEOM_COMMON;
+ float rad;
+};
+
+struct gplane {
+ GEOM_COMMON;
+ float dist;
+};
+
+union gobject {
+ int type;
+ struct gsphere sph;
+ struct gplane plane;
+};
+
+
+int ray_sphere(cgm_ray *ray, union gobject *obj, struct rayhit *hit);
+int ray_plane(cgm_ray *ray, union gobject *obj, struct rayhit *hit);
+
+void texlookup(cgm_vec3 *col, struct texture *tex, float u, float v);
+
+#endif /* GEOM_H_ */
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <SDL.h>
+#include "rend.h"
+
+static void handle_event(SDL_Event *ev);
+static int parse_args(int argc, char **argv);
+
+static int xres = 1280, yres = 800;
+static const char *scnfile;
+static SDL_Surface *fbsurf;
+static int quit;
+
+
+int main(int argc, char **argv)
+{
+ int i;
+ SDL_Event ev;
+
+ if(parse_args(argc, argv) == -1) {
+ return 1;
+ }
+ if(SDL_Init(SDL_INIT_VIDEO | SDL_INIT_TIMER | SDL_INIT_NOPARACHUTE) == -1) {
+ fprintf(stderr, "failed to initialize SDL\n");
+ return 1;
+ }
+
+ if(!(fbsurf = SDL_SetVideoMode(xres, yres, 32, SDL_SWSURFACE))) {
+ fprintf(stderr, "failed to initialize %dx%d video\n", xres, yres);
+ SDL_Quit();
+ return 1;
+ }
+
+ if(rend_init(xres, yres) == -1) {
+ SDL_Quit();
+ return 1;
+ }
+ if(rend_loadscene(scnfile) == -1) {
+ SDL_Quit();
+ return 1;
+ }
+
+ for(;;) {
+ struct rendrect dirty;
+ void *sptr, *dptr;
+
+ while(SDL_PollEvent(&ev)) {
+ handle_event(&ev);
+ if(quit) goto end;
+ }
+ if(rend_step(&dirty) == -1) {
+ break; /* done rendering */
+ }
+
+ SDL_LockSurface(fbsurf);
+ sptr = ((char*)rendimg.pixels + dirty.y * rendimg.pitch) + (dirty.x << 2);
+ dptr = ((char*)fbsurf->pixels + dirty.y * fbsurf->pitch) + (dirty.x << 2);
+ for(i=0; i<dirty.height; i++) {
+ memcpy(dptr, sptr, dirty.width * 4);
+ dptr = (char*)dptr + fbsurf->pitch;
+ sptr = (char*)sptr + rendimg.pitch;
+ }
+ SDL_UnlockSurface(fbsurf);
+ }
+
+ for(;;) {
+ if(SDL_WaitEvent(&ev)) {
+ handle_event(&ev);
+ if(quit) goto end;
+ }
+ }
+
+end:
+ rend_cleanup();
+ SDL_Quit();
+ return 0;
+}
+
+static void handle_event(SDL_Event *ev)
+{
+ switch(ev->type) {
+ case SDL_KEYDOWN:
+ if(ev->key.keysym.sym == 27) {
+ quit = 1;
+ }
+ break;
+ }
+}
+
+static int parse_args(int argc, char **argv)
+{
+ int i;
+
+ for(i=1; i<argc; i++) {
+ if(argv[i][0] == '-') {
+ if(strcmp(argv[i], "-size") == 0) {
+ if(!argv[++i] || sscanf(argv[i], "%dx%d", &xres, &yres) != 2) {
+ fprintf(stderr, "-size must be followed by WxH\n");
+ return -1;
+ }
+ } else {
+ fprintf(stderr, "invalid option: %s\n", argv[i]);
+ return -1;
+ }
+ } else {
+ if(!scnfile) {
+ scnfile = argv[i];
+ } else {
+ fprintf(stderr, "unexpected argument: %s\n", argv[i]);
+ return -1;
+ }
+ }
+ }
+
+ return 0;
+}
--- /dev/null
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <ctype.h>
+#include <float.h>
+#include <errno.h>
+#include "rend.h"
+#include "cgmath/cgmath.h"
+
+struct rendimg rendimg;
+
+static void raytrace(cgm_vec3 *col, cgm_ray *ray, int raydepth);
+static void skycolor(cgm_vec3 *col, cgm_ray *ray);
+static void shade(cgm_vec3 *col, cgm_ray *ray, struct rayhit *hit);
+static void albedo(cgm_vec3 *col, struct rayhit *hit);
+static int find_isect(cgm_ray *ray, struct rayhit *hit);
+static int intersect(cgm_ray *ray, union gobject *obj, struct rayhit *hit);
+static void primary_ray(cgm_ray *ray, int x, int y);
+static void add_sphere(char *args);
+static void add_plane(char *args);
+static void add_light(char *args);
+static void setup_view(char *args);
+static char *clean_line(char *s);
+
+static union gobject **objects;
+static int num_obj;
+static struct light **lights;
+static int num_lights;
+
+static float cam_fov, fov_dist;
+static float cam_xform[16];
+
+static float aspect;
+
+
+int rend_init(int xsz, int ysz)
+{
+ cgm_mtranslation(cam_xform, 0, 0, -10);
+ cam_fov = cgm_deg_to_rad(50.0f);
+
+ if(!(rendimg.pixels = malloc(xsz * ysz * 4))) {
+ fprintf(stderr, "failed to allocate %dx%d framebuffer\n", scanpix, ysz);
+ return -1;
+ }
+ rendimg.width = xsz;
+ rendimg.height = ysz;
+ rendimg.pitch = xsz * 4;
+ aspect = (float)xsz / (float)ysz;
+
+ num_obj = num_lights = 0;
+ return 0;
+}
+
+int rend_loadscene(const char *fname)
+{
+ FILE *fp;
+ char buf[256], *line;
+ int nlines = 0;
+
+ if(!(fp = fopen(fname, "rb"))) {
+ fprintf(stderr, "failed to load scene: %s: %s\n", fname, strerror(errno));
+ return -1;
+ }
+
+ while(fgets(buf, sizeof buf, fp)) {
+ nlines++;
+ line = clean_line(buf);
+ if(!line) continue;
+
+ if(isspace(line[1])) {
+ if(!(arg = clean_line(line + 2))) {
+ fprintf(stderr, "%d: ignoring invalid primitive: %s\n", nlines, line);
+ continue;
+ }
+ switch(line[0]) {
+ case 's':
+ add_sphere(arg);
+ break;
+ case 'p':
+ add_plane(arg);
+ break;
+ case 't':
+ add_triangle(arg);
+ break;
+ case 'l':
+ add_light(arg);
+ break;
+ case 'v':
+ setup_view(arg);
+ break;
+ default:
+ fprintf(stderr, "%d: ignoring unknown primitive '%c'\n", nlines, line[0]);
+ }
+ }
+ }
+
+ printf("added %d objects and %d lights\n", num_obj, num_lights);
+ fclose(fp);
+ return 0;
+}
+
+void rend_cleanup(void)
+{
+ free(rendimg.pixels);
+ rendimg.pixels = 0;
+}
+
+void rend_begin(void)
+{
+ cur_line = 0;
+
+ fov_dist = 1.0f / tan(cam_fov * 0.5f);
+}
+
+#define LINES_PER_STEP 1
+
+int rend_step(struct rendrect *dirty)
+{
+ int i, j;
+ unsigned int r, g, b;
+ uint32_t *pixptr;
+ cgm_ray ray;
+ cgm_vec3 color;
+
+ if(cur_line >= rendimg.height) {
+ return -1;
+ }
+
+ if(dirty) {
+ dirty->x = 0;
+ dirty->y = cur_line;
+ dirty->width = rendimg.width;
+ dirty->height = LINES_PER_STEP;
+ }
+
+ pixptr = (uint32_t*)((char*)rendimg.pixels + rendimg.pitch * cur_line) + rendimg.width;
+ for(i=0; i<LINES_PER_STEP; i++) {
+ for(j=0; j<rendimg.width; j++) {
+ primary_ray(&ray, j, cur_line);
+ raytrace(&color, &ray, 0);
+ r = color.x * 255.0f;
+ if(r > 255) r = 255;
+ g = color.y * 255.0f;
+ if(g > 255) g = 255;
+ b = color.z * 255.0f;
+ if(b > 255) b = 255;
+ *pixptr++ = r | (g << 8) | (b << 16);
+ }
+ cur_line++;
+ pixptr = (uint32_t*)((char*)pixptr - (rendimg.width << 2) + rendimg.pitch);
+ }
+ return 0;
+}
+
+static void raytrace(cgm_vec3 *col, cgm_ray *ray, int raydepth)
+{
+ struct rayhit hit;
+
+ if(raydepth >= MAX_RAY_DEPTH) {
+ return;
+ }
+
+ if(find_isect(ray, &hit)) {
+ shade(&col, ray, &best, raydepth);
+ } else {
+ skycolor(&col, ray);
+ }
+}
+
+static void skycolor(cgm_vec3 *col, cgm_ray *ray)
+{
+ col->x = col->y = col->z = 0.0f; /* TODO */
+}
+
+static void shade(cgm_vec3 *col, cgm_ray *ray, struct rayhit *hit)
+{
+ int i;
+ cgm_vec3 n, l, v, h, alb;
+ cgm_ray sray;
+ float diff, spec;
+ struct light *lt;
+ struct material *mtl = hit->obj->mtl;
+
+ n = hit->n;
+ v = ray->origin; cgm_vsub(&v, &hit->p);
+
+ cgm_vnormalize(&n);
+ cgm_vnormalize(&v);
+
+ albedo(&alb, hit);
+ *col = ambient; cgm_vmul(col, &alb);
+
+ lt = lights;
+ for(i=0; i<num_lights; i++) {
+ shadow_ray(&sray, lt, &hit->p);
+ if(!find_isect(&sray, 0)) {
+ l = sray.dir;
+ cgm_vnormalize(&l);
+ h = l; cgm_vadd(&h, &v); cgm_vnormalize(&h);
+
+ diff = cgm_vdot(&n, &l);
+ if(diff < 0.0f) diff = 0.0f;
+ spec = cgm_vdot(&n, &h);
+ if(spec < 0.0f) spec = 0.0f;
+ spec = pow(spec, mtl->shininess);
+
+ col->x += alb.x * diff + mtl->ks.x * spec;
+ col->y += alb.y * diff + mtl->ks.y * spec;
+ col->z += alb.z * diff + mtl->ks.z * spec;
+ }
+ lt++;
+ }
+}
+
+static void albedo(cgm_vec3 *col, struct rayhit *hit)
+{
+ cgm_vec3 texel;
+ struct material *mtl = hit->obj->mtl;
+
+ *col = mtl->kd;
+
+ if(mtl->texmap) {
+ texlookup(texel, mtl->texmap, hit->uv.x, hit->uv.y);
+ cgm_vmul(col, &texel);
+ }
+}
+
+static int find_isect(cgm_ray *ray, struct rayhit *hit)
+{
+ int i;
+ struct rayhit curhit;
+
+ if(hit) {
+ hit->t = FLT_MAX;
+ for(i=0; i<num_obj; i++) {
+ if(intersect(ray, objects[i], &curhit) && curhit.t < hit->t) {
+ *hit = curhit;
+ }
+ }
+
+ return hit->t == FLT_MAX ? 0 : 1;
+ }
+
+ for(i=0; i<num_obj; i++) {
+ if(intersect(ray, objects[i], 0)) {
+ return 1;
+ }
+ }
+ return 0;
+}
+
+static int intersect(cgm_ray *ray, union gobject *obj, struct rayhit *hit)
+{
+ switch(obj->type) {
+ case GOBJ_SPHERE:
+ return ray_sphere(ray, &obj->sph, hit);
+ case GOBJ_PLANE:
+ return ray_plane(ray, &obj->plane, hit);
+ default:
+ break;
+ }
+ return 0;
+}
+
+static void primary_ray(cgm_ray *ray, int x, int y)
+{
+ cgm_vcons(&ray->origin, 0, 0, 0);
+ ray->dir.x = ((float)(x << 1) / (float)rendimg.width - 1.0f) * aspect;
+ ray->dir.y = 1.0f - (float)(y << 1) / (float)rendimg.height;
+ ray->dir.z = fov_dist;
+
+ cgm_rmul_mr(ray, cam_xform);
+}
+
+static void add_sphere(char *args)
+{
+}
+
+static void add_plane(char *args);
+static void add_light(char *args);
+static void setup_view(char *args);
+
+static char *clean_line(char *s)
+{
+ char *end;
+
+ while(*s && isspace(*s)) s++;
+ if(!*s) return 0;
+
+ if((end = strchr(s, '#'))) {
+ *end-- = 0;
+ } else {
+ end = s + strlen(s) - 1;
+ }
+ while(end > s && isspace(*end)) *end-- = 0;
+
+ return *s ? s : 0;
+}
--- /dev/null
+#ifndef REND_H_
+#define REND_H_
+
+#ifdef __sgi
+#include <sys/types.h>
+#else
+#include <stdint.h>
+#endif
+
+#define MAX_RAY_DEPTH 6
+
+struct rendrect {
+ int x, y;
+ int width, height;
+};
+
+struct rendimg {
+ int width, height, pitch;
+ uint32_t *pixels;
+};
+
+struct rayhit {
+ float t;
+ cgm_vec3 p;
+ cgm_vec3 n;
+ cgm_vec2 uv;
+ void *obj, *subobj;
+};
+
+extern struct rendimg rendimg;
+
+int rend_init(int xsz, int ysz);
+void rend_cleanup(void);
+
+void rend_begin(void);
+int rend_step(struct rendrect *dirty);
+
+#endif /* REND_H_ */
projects / ray24 / commitdiff