From 4131f7e2d110aab8e1805f3229cc38676922f4f4 Mon Sep 17 00:00:00 2001 From: Anton Lydike Date: Thu, 13 Aug 2020 16:00:08 +0200 Subject: [PATCH] started rework, points now use avx --- .gitignore | 2 + marcher.h | 31 +-------- src/camera.h | 13 ++++ src/point.c | 182 ++++++++++----------------------------------------- src/point.h | 48 ++++++++++++++ src/rays.c | 39 +++++++++++ src/rays.h | 32 +++++++++ test | Bin 0 -> 18416 bytes 8 files changed, 170 insertions(+), 177 deletions(-) create mode 100644 src/camera.h create mode 100644 src/point.h create mode 100644 src/rays.c create mode 100644 src/rays.h create mode 100755 test diff --git a/.gitignore b/.gitignore index 877906e..4f02e38 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,6 @@ *.bmp *.png *.out +*.old.c +test.c out/ \ No newline at end of file diff --git a/marcher.h b/marcher.h index acc539a..e19bf57 100644 --- a/marcher.h +++ b/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" diff --git a/src/camera.h b/src/camera.h new file mode 100644 index 0000000..fceb64b --- /dev/null +++ b/src/camera.h @@ -0,0 +1,13 @@ +/* +This handles code for generating rays +*/ + +#ifndef MARCHER_CAMERA_H +#define MARCHER_CAMERA_H + + + + + + +#endif \ No newline at end of file diff --git a/src/point.c b/src/point.c index 1324312..9d71545 100644 --- a/src/point.c +++ b/src/point.c @@ -1,180 +1,64 @@ #include #include -// 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 - ); } \ No newline at end of file diff --git a/src/point.h b/src/point.h new file mode 100644 index 0000000..3c477eb --- /dev/null +++ b/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 , 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 \ No newline at end of file diff --git a/src/rays.c b/src/rays.c new file mode 100644 index 0000000..68ef6d8 --- /dev/null +++ b/src/rays.c @@ -0,0 +1,39 @@ +#import +#include + +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); + } +} + diff --git a/src/rays.h b/src/rays.h new file mode 100644 index 0000000..3aeff22 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test b/test new file mode 100755 index 0000000000000000000000000000000000000000..bab91cebb7d21f3a1d9d759d7c231f960ac8df3e GIT binary patch literal 18416 zcmeHPdu$v>8K1qiW2a8eXD4Z#rX|}pq3Pr9{H~qgWY3P1bL%v29#Cl0^*#0;_1&Ji z+jATj6<47JtE3I6hzL{(sr(WCKq>{H@CVc(rA4ZQ(iVge1qDQ@fj-;^QAE`o-#0to zd3Sw5w9x-{;nCebNn(_IYzISbtVQJ_>&)2XPR=;n=HY67%E2U|dnTCDCrqF-DYd zcLEOU`QIiA`+ia0<>E3X>Q7XJmxYb&bhu;|ih6iD-J8wD`_uh1nQZ4wx+6K6?Vd?Y zcJ}oZGscur3Z`?}cp}&t>S_-|E}TmSb6KN29h~m#4fb}2N`+7-j}_Jl{qMQs5Wfs8 z8m~!#gW4y`&<6Qh33YVT(qB&228z}i${}#3)Q~q_#-8D|;wAMmO`9xc4Kt-BG84d{ zLMxfEqN(L$S%d8v9T^@Gg!9kxT0F=MK|qcfdpkHsc3 zebdH7XS}OtDwdhZ&D`3Vp9rQCiC}-IGt>@^0@_A44xvit3x)`l264uFPWM{^4(~6u zUvootC2DUHvhr2P`&cjVm=<@coMJmd9HomkxoF~CwIJjTz)_im@u(L8F9Kc!ya;#^ z@FL(vz>9zvf&X^|qQ_ouj6PETMwqeaL-VGuaxQx8*~YWgVk$jPgQ~LSM{sW5q(Tn) z4C?=OzEY_m{+;H)mS>w819}b8=<%NKS#14Jfj)!d;L|v`O*L=*H5@N2{1}ms)IUaA zK(70ig)9I`MVpOZg8qoe&vW#`ETB@&Z@q-&n4kh$j}*%52jVE_g@yfK{ciIE>+Z8G zSMRqhn=Q8fb{FE_mlhPORvSnIe*fgRYIBZP2~VOf?_H`?I_7z-QlaKeU&4SufLO;H zBXcj_8=d=e^w?h(#}1A>KYyZCWzpx)v~%*j3hx~%>!w1Cz1O^H7sMr?f_ebY=9<&! z1aHI2XHZZ%pE_|iWp}T#pK<#K-{fR&F?#0BozXKF<*4*r^uF`-zAGiHK9(9%<`OL`NZzG=`nY$1<6q)-Va&YwcmdwqJ z9f}?gB0VsAj%74fB z-F-N6Pvmgq-iQV@VA%3_n{>27jEk4N2zU|jBH%^9i+~paF9Kc!ya;#^@FL(v;D3n# ze&5=3@|x9Ituhj{u8=>>&(TiE|MbzJStYME|z0cnl@wZ3( zTL=BhpuY|52mOuw`{8#1u7Dr=cR_UNB`*SA1iT1%5%415MZk-I7XdE|ik`Pvq|tf+xH8E7tXi>N5BY%&%7w{@{Zu=SdIK zRL-{wy>R!j5n3niE})mD@NR{R-!n9{P1y=7qHk3CV5x9BvkCaJS(p7Sz+3mgiV*D^ zu+Poy5Ntm_hM%A*qINA~pW=FBNiXW%D$Jnvi*j}VrtZf~CchD#eLw?l`IHCF;IDBisgS5VB8{}8Yv~fVLtA|C} ztgZ{B3|p6vjlj84I)(yB9e`RY!Lqm;*P% z&Szb`^uB*JlOQ1!{z9Q}>1ZjbPnN>4&|d59W{j^s4;?@?guq9!ly^vSAT)UWj=h^V zm#5@_nej>B>3Z%PA7GN(H_IAXGn7C#ECz>r{dzGdYkp$^Z?fi~_I#FoT6Dy_m`sj(~7%WC+e+($^o%8nAsr2R5%4(|RzL z%LNk!1FCD7C6E>iWh0q2(omk6K`S6Zy;w8~z;NfW#6)niSTOYj?B4*PnJWeJg`}Pv zz}uRQ<%0UOo?zO@K3IeguC#JVPh#USvx#6)&lM7I>z4;i+}fH7yd6Bj^Oz@#gEKu( z1f@GqH>heSGHyKG-1P!axF^JO#f|?@ybM2#Q$BT>i40BeVokYR?Pu*E8J!`<0CPQK$ThYJcmLL+>u)=9I5uN=^A{rq-1E8QnLh zpJqmnS*JX}Xg`xv-oj2Qs$InMGHYct9=>|1hsu>$3LtP3~uTd;5;C zcemRga{2Q1{~Y(TynVg`@>aNBT(yXa!-1~#F{ZHP?e=vr2!Pz(ZvTpYe9M2Q@1mT0 zcAi?nKl*Y1E;jLUowmq@K;GiJk~ONTL%AL0Z0R{bfUW&lUSN@^uvf+Nihn}_uf}syjDTFd zqN$f~Tx8pR7w?)WO>}27xyfucKh-srO7-=|<2~`|bS9mgiT5vUfa9AZVLO}{t_91+ zgi$zRh|QAF_lNm_KMehBTu*0>P=O`&BIHV@UewIImdF(hy#)K`k_9cDE5u_tEeSn* zNsE=IS)!1i%;~0{4Bgt**IwI5OTkrXu_A0M)O5ov9%ZRwEU#vx!b z+}ErdNw{K7+r2-scUT*~W0$4@=cq=KX}dmqM`Z8FkW<5lav*~7nKm319-_PUGi}f4 zzQM?-wr}_D1H%WkgOS0}Ve|^y5EB{<^LRj*rYB=&jKtg--@>V?w1*S-Z9>?=%G$OF zjAIWY=t>1G6El)I-94_dm4+7s(0HjN-0+Q}_I;#FMdG$mO*^=ENQ8ogN=Nf%EDqEx zS~NpBnBeHeNft6-%omCbjs(qEnuVZ6g)*^HhJ}(x4FFrzELs{EdzP{V!@+2vgYjPs z0U|e(_kJm($Llgn>A7e%E#jIqCQ0&4wWen;h?O5Fa7>MY7TpFcw|AJWF;Cu z(mwhzAt!>i_JotH1~>>$<4oE&uGozsPtQ*h_do_`n1rYIK5epPu~`_>{2Jz$JdvLc z;4Z*B5yGpo%Jr(iebM^h?)Uf|~i zp58+VH$2vqM|wLPFej+ipT^6rqQX=j)t|~CdKcjFoe}Xgu2uw~9SecvuD>N^EDafq zCE;oOJ?6qw`BWde{?7~ipl~=V-b>920NEjn@7~<_`$2|n3hT!cjn4xXOYhr+c3JzG z0)-u~2)rWj3IPQbV(4?N^dJyBp88u=7~DyLap48azU;!&JmPUND53X%WKZpu_~Qa! zEk;o-m12=`i&W24ZxVb1GFX4Y)4XI{u`J-Ni6gnoe-gBIJk4KDDi&p~R z{BVfoUY+OA|6`#0P4>hS{Rt>wbXCREeCV|Je-&hRIs7j`fo+G%r+E|ozYF@m7Sz7T zAI0%kfW_-4Jk77p;R_*1?G~w?2~YH8(Ae?S@umP$IplzF#QzR5D5LUe{1y;+w>_3+ zPk#gkW$t