started rework, points now use avx

.gitignore

@@ -1,4 +1,6 @@
 *.bmp
 *.png
 *.out
+*.old.c
+test.c
 out/

.gitmodules

@@ -1,3 +1,3 @@
 [submodule "images"]
 	path = images
-	url = https://git.datenvorr.at/anton/images.h.git
+	url = gitlab@git.datenvorr.at:anton/images.h.git

images

@@ -1 +1 @@
-Subproject commit 329520da739f07aded44841eb1b5de6e5903a425
+Subproject commit e20f4cb891ae1f6ac34f7572ea6847957b7195db

marcher.h

@@ -3,42 +3,18 @@
 
 #include "images/images.h"
 
+#include "src/point.h"
+#include "src/rays.h"
+
 // define pi if not available
 #ifndef M_PI
 #define M_PI 3.14159265358979323846
 #endif
 
-struct __myvec;
-struct __mymtrx;
 struct __mycam;
 struct __myobject;
 struct __myscene;
 
-typedef struct __myvec {
-    double x;
-    double y;
-    double z;
-} Point;
-
-typedef struct __mymtrx {
-    double entries[9];
-} Matrix;
-
-inline Point pt_new(double x, double y, double z);
-inline Point pt_scale(Point pt, double length);
-inline Point pt_normalize(Point pt);
-inline double pt_length(Point pt);
-inline void pt_add(Point* pt, Point add);
-inline void pt_sub(Point* pt, Point sub);
-inline double pt_dist(Point p1, Point p2);
-inline Point pt_mod(Point pt, double mod);
-inline double pt_dot(Point a, Point b);
-inline Point pt_cross(Point a, Point b);
-inline double pt_angle(Point a, Point b);
-inline void pt_print(Point pt);
-inline void pt_print_n(const char* name, Point pt);
-
-
 typedef struct __mycam {
     Point location;
     Point direction;
@@ -90,7 +66,6 @@ void scene_add_obj(Scene* scene, SceneObject object);
 
 void scene_destroy(Scene scene);
 
-#include "src/point.c"
 #include "src/camera.c"
 #include "src/scene.c"
 

src/camera.h

@@ -0,0 +1,13 @@
+/*
+This handles code for generating rays
+*/
+
+#ifndef MARCHER_CAMERA_H
+#define MARCHER_CAMERA_H
+
+
+
+
+
+
+#endif 

src/point.c

@@ -1,180 +1,64 @@
 #include <math.h>
 #include <stdio.h>
 
-// basically a vector3
+inline void   pt_print(point x) {
-inline Point pt_new(double x, double y, double z) {
+    printf("(%f,%f,%f)\n", x[0], x[1], x[2]);
-    Point pt;
-    pt.x = x;
-    pt.y = y;
-    pt.z = z;
-    return pt;
 }
-
+inline void   pt_print(point x, const char* name) {
-// scale vector to length
+    printf("%s = (%f,%f,%f)\n", name, x[0], x[1], x[2]);
-inline Point pt_scale(Point pt, double length) {
-    double f = length / pt_length(pt);
-    return pt_new(
-        pt.x * f,
-        pt.y * f,
-        pt.z * f
-    );
 }
 
-inline Point pt_mult(Point pt, double scalar) {
+// length related functions
-    return pt_new(
+inline point  pt_scale(point x, double length) {
-        pt.x * scalar,
+    length = length / pt_length(x);
-        pt.y * scalar,
+    return x * length;
-        pt.z * scalar
+}
-    );
+inline double pt_length(point x) {
+    return sqrt((x[0] * x[0]) + (x[1] * x[1]) + (x[2] * x[2]));
+}
+inline point  pt_normalize(point x) {
+    return pt_scale(x, 1);
 }
 
-// return internal angle between a and b
+// angle related functions
-inline double pt_angle(Point a, Point b) {
+inline double pt_angle(point x, point y) {
     return acos(pt_dot(
         pt_normalize(a),
         pt_normalize(b)
     ));
 }
 
-// get the length of vector
+// miscelaneous vector operations:
-inline double pt_length(Point pt) {
+inline double pt_dot(point x, point y) {
-    return sqrt((pt.x * pt.x) + (pt.y * pt.y) + (pt.z * pt.z));
+        return x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
 }
-
+inline point  pt_cross(point x, point y) {
-// add the vector add to the vector pt
+    return (point){
-inline void pt_add(Point* pt, Point add) {
+        x[1]*y[2] - x[2]*y[2], 
-    pt->x = pt->x + add.x;
+        x[2]*y[0] - x[0]*y[2],
-    pt->y = pt->y + add.y;
+        x[0]*y[2] - x[2]*y[0]
-    pt->z = pt->z + add.z;
+    };
 }
-
+inline point  pt_mod(point x, double mod) {
-// add the vector add to the vector pt
+    return (point){fmod(x[0], mod), fmod(x[1], mod), fmod(x[2], mod)};
-inline void pt_sub(Point* pt, Point sub) {
-    pt->x -= sub.x;
-    pt->y -= sub.y;
-    pt->z -= sub.z;
-}
-
-inline double pt_dist(Point p1, Point p2) {
-    pt_sub(&p1, p2);
-    return pt_length(p1);
-}
-
-// normalize a vector
-inline Point pt_normalize(Point pt) {
-    return pt_scale(pt, 1);
 }
 
-// dot product of two vectors
+// call fabs on each entry
-inline double pt_dot(Point a, Point b) {
+inline point  pt_abs(point x) {
-    return a.x*b.x + a.y*b.y + a.z*b.z;
+    return (point){fabs(x[0]), fabs(x[1]), fabs(x[2])}
-}
-
-// cross product of two vectors
-inline Point pt_cross(Point a, Point b) {
-    return pt_new(
-        a.y*b.z - a.z*b.y, 
-        a.z*b.x - a.x*b.z,
-        a.x*b.y - a.y*b.x
-    );
-}
-
-inline void pt_print(Point pt) {
-    printf("(%f, %f, %f)\n", pt.x, pt.y, pt.z);
-}
-
-inline void pt_print_n(const char* name, Point pt) {
-    printf("%s: (%f, %f, %f)\n", name, pt.x, pt.y, pt.z);
 }
 
 // find two vectors that span the orthogonal plane, where 
 // span_xy is a vector lying on the xy-plane (and pointing left)
 // and span_z is orthogonal to span_xy pointing "upwards"
-void pt_orthogonal_plane(Point pt, Point *span_z, Point *span_xy) {
+void pt_orthogonal_plane(point pt, point *span_z, point *span_xy) {
     pt = pt_normalize(pt);
+    point left = (point){0,0,1};
 
     // get the vector lying on the xy axis
     // this is done by 
-    *span_xy = pt_normalize(pt_cross(pt_new(0,0,1), pt)); // points to the "left" (of the viewing direction)
+    *span_xy = pt_normalize(pt_cross(left, pt)); // points to the "left" (of the viewing direction)
 
     // now use this, to find the vector 
     *span_z = pt_normalize(pt_cross(pt, *span_xy));
 }
-
-inline Point pt_mod(Point pt, double mod) {
-    return pt_new(
-        fabs(fmod(pt.x, mod)),
-        fabs(fmod(pt.y, mod)),
-        fabs(fmod(pt.z, mod))
-    );
-}
-
-
-///////////////////////////////
-////// Matrix operations //////
-///////////////////////////////
-
-
-/* create a new matrix with entries:
-    x1 x2 x3
-    y1 y2 y3
-    z1 z2 z3
-*/
-inline Matrix mtrx_new(double x1, double x2, double x3,
-                       double y1, double y2, double y3,
-                       double z1, double z2, double z3) 
-{
-    Matrix m;
-    m.entries[0] = x1;
-    m.entries[1] = y1;
-    m.entries[2] = z1;
-    m.entries[3] = x2;
-    m.entries[4] = y2;
-    m.entries[5] = z2;
-    m.entries[6] = x3;
-    m.entries[7] = y3;
-    m.entries[8] = z3;
-    return m;
-}
-
-inline Point mtrx_mult(Matrix mtrx, Point pt) {
-    Point result;
-
-    double *m = mtrx.entries;
-
-    result.x = m[0] * pt.x + m[3] * pt.y + m[6] * pt.z;
-    result.y = m[1] * pt.x + m[4] * pt.y + m[7] * pt.z;
-    result.z = m[2] * pt.x + m[5] * pt.y + m[8] * pt.z;
-
-    return result;
-}
-
-// create a rotation matrix around an axis given by the normalized axis vector (u)
-// taken from https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle
-Matrix mtrx_rotation(Point u, double theta) {
-    double theta_rad = theta * (M_PI / 180);
-    double cost = cos(theta_rad);
-    double sint = sin(theta_rad);
-
-    return mtrx_new(
-        cost+u.x*u.x*(1-cost),     u.x*u.y*(1-cost)-u.z*sint, u.x*u.z*(1-cost)+u.y*sint,
-        u.y*u.x*(1-cost)+u.z*sint, cost+u.y*u.y*(1-cost),     u.y*u.z*(1-cost)-u.x*sint,
-        u.z*u.x*(1-cost)-u.y*sint, u.z*u.y*(1-cost)+u.x*sint, cost+u.z*u.z*(1-cost)
-    );
-}
-
-void mtrx_print(Matrix mtrx) {
-    printf("  %8.2f %8.2f %8.2f\n  %8.2f %8.2f %8.2f\n  %8.2f %8.2f %8.2f\n",
-        mtrx.entries[0], mtrx.entries[3], mtrx.entries[6],
-        mtrx.entries[1], mtrx.entries[4], mtrx.entries[7],
-        mtrx.entries[2], mtrx.entries[5], mtrx.entries[8]
-    );
-}
-
-inline Matrix mtrx_outer_prod(Point a, Point b) {
-    return mtrx_new(
-        a.x*b.x, a.x*b.y, a.x*b.z,
-        a.y*b.x, a.y*b.y, a.y*b.z,
-        a.z*b.x, a.z*b.y, a.z*b.z
-    );
-}

src/point.h

@@ -0,0 +1,48 @@
+/*
+This library for working with 3D Points uses AVX instructions to speed up computation.
+
+It declares the type <point>, which should be treated as immutable once initialized.
+
+All operations conducted on points are prefixed with pt_
+
+In addition to the functions declared here, the following operators are supported:
+    point [+|-|*|/] [point|double]          (2.0 is treated as {2.0, 2.0, 2.0})
+    point [=|<=|>=] [point|double]
+
+Constructing a point is done like this: point x = (point) {1.0, 2.0, 3.0}
+
+*/
+
+#ifndef MARCHER_POINT_H
+#define MARCHER_POINT_H
+
+// vector of type double with length 3
+typedef double point __attribute__ ((vector_size (sizeof(double) * 4)));
+
+inline void   pt_print(point x);
+inline void   pt_print(point x, const char* name);
+
+// length related functions
+inline point  pt_scale(point x, double length);
+inline double pt_length(point x);
+inline point  pt_normalize(point x);
+
+// angle related functions
+inline double pt_angle(point x, point y);
+
+// miscelaneous vector operations:
+inline double pt_dot(point x, point y);
+inline point  pt_cross(point x, point y);
+inline point  pt_mod(point x, double mod);
+
+// call fabs on each entry
+inline point  pt_abs(point x);
+
+// find two vectors that span the orthogonal plane, where 
+// span_xy is a vector lying on the xy-plane (and pointing left)
+// and span_z is orthogonal to span_xy pointing "upwards"
+void pt_orthogonal_plane(point pt, point *span_z, point *span_xy)
+
+#include "point.c"
+
+#endif

src/rays.c

@@ -0,0 +1,39 @@
+#import <math.h>
+#include <stdlib.h>
+
+void ray_advance(ray *r, double distance) {
+    ray->position += pt_scale(ray->direction, distance);
+}
+
+
+// try sqrt(size) random rays from the bunch and calulate the max distance for every one of them
+// take the smallest max distance, this is our guess.
+double bundle_dispersion(ray_bundle bundle) {
+    int samples = sqrt(bundle.size); // this should be enough points
+
+    // large number
+    double guess = 10000000;
+
+    for (int i = 0; i < samples; i++) {
+        searched[i] = index;
+        double dist = 0;
+        point x = bundle.rays[rand() % bundle.size].position;
+        for (int j = 0; j < bundle.size; j++) {
+            double curr_dist = pt_dist(x, bundle.rays[j].position);
+            if (curr_dist > dist) {
+                dist = curr_dist;
+            }
+        }
+
+        if (guess < dist) {
+            guess = dist;
+        }
+    }
+}
+
+void bundle_advance(ray_bundle bundle, double distance) {
+    for (int i = 0; i < bundle.size; i++) {
+        ray_advance(bundle.rays + i, distance);
+    }
+}
+

src/rays.h

@@ -0,0 +1,32 @@
+/*
+This file defines rays and ray bundles. A ray is at a location and has a direction.
+
+RAYS ARE MUTABLE!
+
+*/
+#ifndef MARCHER_RAYS_H
+#define MARCHER_RAYS_H
+
+#include "point.h"
+
+// 
+typedef struct _ray {
+    point direction;
+    point position;
+} ray;
+
+void ray_advance(ray r, double distance);
+
+typedef struct _ray_bundle {
+    ray[] rays;
+    int count;
+} ray_bundle;
+
+// upper bound for distances between points (not the least upper bound)
+double bundle_dispersion(ray_bundle bundle);
+
+void bundle_advance(ray_bundle bundle, double distance);
+
+#include "rays.c"
+
+#endif