started rework, points now use avx
Anton Lydike
4 years ago
parent
266f243703
commit
4131f7e2d1
@ -1,4 +1,6 @@
|
||||
*.bmp
|
||||
*.png
|
||||
*.out
|
||||
*.old.c
|
||||
test.c
|
||||
out/
|
||||
@ -0,0 +1,13 @@
|
||||
/*
|
||||
This handles code for generating rays
|
||||
*/
|
||||
|
||||
#ifndef MARCHER_CAMERA_H
|
||||
#define MARCHER_CAMERA_H
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
#endif
|
||||
@ -1,180 +1,64 @@
|
||||
#include <math.h>
|
||||
#include <stdio.h>
|
||||
|
||||
// basically a vector3
|
||||
inline Point pt_new(double x, double y, double z) {
|
||||
Point pt;
|
||||
pt.x = x;
|
||||
pt.y = y;
|
||||
pt.z = z;
|
||||
return pt;
|
||||
inline void pt_print(point x) {
|
||||
printf("(%f,%f,%f)\n", x[0], x[1], x[2]);
|
||||
}
|
||||
|
||||
// scale vector to length
|
||||
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 void pt_print(point x, const char* name) {
|
||||
printf("%s = (%f,%f,%f)\n", name, x[0], x[1], x[2]);
|
||||
}
|
||||
|
||||
inline Point pt_mult(Point pt, double scalar) {
|
||||
return pt_new(
|
||||
pt.x * scalar,
|
||||
pt.y * scalar,
|
||||
pt.z * scalar
|
||||
);
|
||||
// length related functions
|
||||
inline point pt_scale(point x, double length) {
|
||||
length = length / pt_length(x);
|
||||
return x * length;
|
||||
}
|
||||
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
|
||||
inline double pt_angle(Point a, Point b) {
|
||||
// angle related functions
|
||||
inline double pt_angle(point x, point y) {
|
||||
return acos(pt_dot(
|
||||
pt_normalize(a),
|
||||
pt_normalize(b)
|
||||
));
|
||||
}
|
||||
|
||||
// get the length of vector
|
||||
inline double pt_length(Point pt) {
|
||||
return sqrt((pt.x * pt.x) + (pt.y * pt.y) + (pt.z * pt.z));
|
||||
// miscelaneous vector operations:
|
||||
inline double pt_dot(point x, point y) {
|
||||
return x[0]*x[0] + x[1]*x[1] + x[2]*x[2];
|
||||
}
|
||||
|
||||
// add the vector add to the vector pt
|
||||
inline void pt_add(Point* pt, Point add) {
|
||||
pt->x = pt->x + add.x;
|
||||
pt->y = pt->y + add.y;
|
||||
pt->z = pt->z + add.z;
|
||||
inline point pt_cross(point x, point y) {
|
||||
return (point){
|
||||
x[1]*y[2] - x[2]*y[2],
|
||||
x[2]*y[0] - x[0]*y[2],
|
||||
x[0]*y[2] - x[2]*y[0]
|
||||
};
|
||||
}
|
||||
|
||||
// add the vector add to the vector pt
|
||||
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);
|
||||
inline point pt_mod(point x, double mod) {
|
||||
return (point){fmod(x[0], mod), fmod(x[1], mod), fmod(x[2], mod)};
|
||||
}
|
||||
|
||||
// dot product of two vectors
|
||||
inline double pt_dot(Point a, Point b) {
|
||||
return a.x*b.x + a.y*b.y + a.z*b.z;
|
||||
}
|
||||
|
||||
// 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);
|
||||
// call fabs on each entry
|
||||
inline point pt_abs(point x) {
|
||||
return (point){fabs(x[0]), fabs(x[1]), fabs(x[2])}
|
||||
}
|
||||
|
||||
// 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
|
||||
);
|
||||
}
|
||||
@ -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
|
||||
@ -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);
|
||||
}
|
||||
}
|
||||
|
||||
@ -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
|
||||
Loading…
Reference in New Issue