diff --git a/.gitignore b/.gitignore index 877906e..bf48f6c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,4 +1,5 @@ *.bmp *.png *.out -out/ \ No newline at end of file +obj +out diff --git a/.gitmodules b/.gitmodules index 3c9ded2..bb3729c 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,3 +1,3 @@ [submodule "images"] path = images - url = gitlab@git.datenvorr.at:anton/images.h.git + url = https://git.datenvorr.at/anton/images.h.git diff --git a/.idea/.gitignore b/.idea/.gitignore new file mode 100644 index 0000000..13566b8 --- /dev/null +++ b/.idea/.gitignore @@ -0,0 +1,8 @@ +# Default ignored files +/shelf/ +/workspace.xml +# Editor-based HTTP Client requests +/httpRequests/ +# Datasource local storage ignored files +/dataSources/ +/dataSources.local.xml diff --git a/.idea/misc.xml b/.idea/misc.xml new file mode 100644 index 0000000..79b3c94 --- /dev/null +++ b/.idea/misc.xml @@ -0,0 +1,4 @@ + + + + \ No newline at end of file diff --git a/.idea/modules.xml b/.idea/modules.xml new file mode 100644 index 0000000..7a8cdab --- /dev/null +++ b/.idea/modules.xml @@ -0,0 +1,8 @@ + + + + + + + + \ No newline at end of file diff --git a/.idea/vcs.xml b/.idea/vcs.xml new file mode 100644 index 0000000..46bff8d --- /dev/null +++ b/.idea/vcs.xml @@ -0,0 +1,7 @@ + + + + + + + \ No newline at end of file diff --git a/Makefile b/Makefile new file mode 100644 index 0000000..fd0b93f --- /dev/null +++ b/Makefile @@ -0,0 +1,31 @@ + + +OPTIMIZATION=-O3 -flto +CC=gcc +CFLAGS=-Isrc/ -lm -lpthread -Wall -Wextra -pedantic-errors $(OPTIMIZATION) -DPOINT_DTYPE=float + +.PHONY: directories + +directories: + mkdir -p obj out + +clean: + rm -rf obj/* out/* + +obj/scene.o: src/scene.c src/scene.h src/point.h + $(CC) $(CFLAGS) -c -o $@ src/scene.c + +obj/camera.o: src/camera.c src/camera.h src/point.h + $(CC) $(CFLAGS) -c -o $@ src/camera.c + +obj/images.o: images/src/images.c images/src/images.h src/point.h + $(CC) $(CFLAGS) -c -o $@ images/src/images.c + +march: obj/camera.o obj/scene.o obj/images.o src/point.h + $(CC) $(CFLAGS) -o out/march $^ marcher.c + +bench: obj/camera.o obj/scene.o obj/images.o src/point.h + $(CC) $(CFLAGS) -o out/bench $^ bench.c + +gpu: obj/camera.o obj/images.o src/point.h + $(CC) -fopenacc $(CFLAGS) -o out/gpu $^ gpu.c diff --git a/bench.c b/bench.c new file mode 100644 index 0000000..138e588 --- /dev/null +++ b/bench.c @@ -0,0 +1,181 @@ +#include +#include +#include +#include +#include +#include +#include +#include "images/src/images.h" +#include "src/scene.h" +#include "src/camera.h" +#include "src/point.h" + + +#define BENCH_VERSION "1.0" + + +/* + + Mandelbulb scene object + + Currently cannot be set at a specific location, always resides at origin (0,0,0) + + Color function is just a flat shader, detail is displayed with ambient occlusion + +*/ +double mandelbulb_dist(struct point pt, struct scene_object *self) { + int iters = (int) self->args[0]; + double power = self->args[1]; + + struct point z = pt; + double dr = 1.0; + double r = 0.0; + for (int i = 0; i < iters ; i++) { + r = pt_length(z); + + if (r>2) { + break; + } + + // convert to polar coordinates + double theta = acos(z.z/r); + double phi = atan2(z.y,z.x); + dr = pow(r, power-1.0)*power*dr + 1.0; + + // scale and rotate the point + double zr = pow(r, power); + theta = theta*power; + phi = phi*power; + + // convert back to cartesian coordinates + z = (struct point) { + .x = sin(theta)*cos(phi) * zr + pt.x, + .y = sin(phi)*sin(theta) * zr + pt.y, + .z = cos(theta) * zr + pt.z + }; + } + + return 0.5*log(r)*r/dr; + +} + +Color mandelbulb_color(struct point hit, struct point direction, struct scene_object *self) { + return self->color; +} + +// constructs the scene object +struct scene_object mandelbulb_new(struct point location, int iters, double power) { + struct scene_object so; + so.location = location; + so.args = malloc(sizeof(double) * 3); + so.args[0] = iters; // iterations + so.args[1] = power; // power + so.args[2] = -1; // reserved for color calculations + so.distance = mandelbulb_dist; + so.get_color = mandelbulb_color; + so.color = color_new(255,255,255); + return so; +} + +int run_bench(int size, double pow, int threads, const char path[], int save) { + double cam_position = 1.15; + int steps = 2000; + int iters = 1000; + double threshold = 0.0001; + + struct camera cam; + cam.fov = 90; + + camera_set_looking_at(&cam, (struct point){.x = cam_position, .y = cam_position, .z = cam_position}, PT_ZERO); + + // create basic scene with up to 10 objects + struct scene scene = scene_new(size, size, 1); + scene.perf_opts.max_steps = steps; + scene.perf_opts.threshold = threshold; + scene.perf_opts.speed_cutoff = 10; + scene.background = color_new(0,0,0); + + scene_add_obj(&scene, mandelbulb_new(PT_ZERO, iters, pow)); + + Image *img = render_scene(&scene, &cam, threads); + + if (save) { + image_save_bmp(*img, path); + } + + image_destroy(*img); + scene_destroy(scene); + + return 0; +} + +struct timer { + struct timeval start; + struct timeval end; + char *name; +}; + +void timer_start(struct timer *t) { + printf("\nStarting bench %s\n", t->name); + gettimeofday(&t->start, NULL); +} +void timer_end(struct timer *t) { + gettimeofday(&t->end, NULL); +} +void timer_print(struct timer t) { + long time, secs, usecs; + secs = t.end.tv_sec - t.start.tv_sec; + usecs = t.end.tv_usec - t.start.tv_usec; + time = ((secs) * 1000 + usecs/1000.0) + 0.5; + printf("\nBenchmark %s took %ldms (%.2fs)\n", t.name, time, time / 1000.0f); +} + +int main() +{ + int threads = get_nprocs() / 2; + struct timer bench; + + printf("Mandelbulb Benchmark v%s\n\nDetected %d threads...\n", BENCH_VERSION, threads); + + sleep(2); + + bench.name = "1080px render with saving"; + timer_start(&bench); + run_bench(1080, 3.0, threads, "bench-pow3-1080p.bmp", 1); + timer_end(&bench); + timer_print(bench); + + sleep(2); + + bench.name = "1080px render without saving"; + timer_start(&bench); + run_bench(1080, 3.0, threads, "", 0); + timer_end(&bench); + timer_print(bench); + + sleep(2); + + bench.name = "10 megapixel render with saving"; + timer_start(&bench); + run_bench(3162, 3.0, threads, "bench-pow3-10mpx.bmp", 1); + timer_end(&bench); + timer_print(bench); + + sleep(2); + + bench.name = "40 megapixel render with saving"; + timer_start(&bench); + run_bench(6324, 3.0, threads, "bench-pow3-40mpx.bmp", 1); + timer_end(&bench); + timer_print(bench); + + sleep(2); + + bench.name = "1080px render single threaded without saving"; + timer_start(&bench); + run_bench(1080, 3.0, 1, "", 0); + timer_end(&bench); + timer_print(bench); + + return 0; +} diff --git a/gpu.c b/gpu.c new file mode 100644 index 0000000..44470b4 --- /dev/null +++ b/gpu.c @@ -0,0 +1,186 @@ +// use floats instead of doubles +#define DTYPE float +// #define POINT_DTYPE DTYPE + +#include +#include + +#include "src/point.h" +#include "src/camera.h" +#include "images/src/images.h" + +#define ITERS 100 +#define POWER 3 + +#define SIZE 1000 +#define STEP_WIDTH (2 / ((DTYPE) (SIZE - 1))) + +#define CAM_POSITION 1.35 + +#define STEPS 100 +#define THRESHOLD 0.001 + +#include + + +float mandelbulb_dist(struct point pt) +{ + DTYPE power = POWER; + + struct point z = pt; + DTYPE dr = 1.0; + DTYPE r = 0.0; + for (int i = 0; i < ITERS ; i++) { + r = pt_length(z); + + if (r>2) { + break; + } + + // convert to polar coordinates + DTYPE theta = acos(z.z/r); + DTYPE phi = atan2(z.y,z.x); + dr = pow(r, power-1.0)*power*dr + 1.0; + + // scale and rotate the point + DTYPE zr = pow(r, power); + theta = theta*power; + phi = phi*power; + + // convert back to cartesian coordinates + z = (struct point) { + .x = sin(theta)*cos(phi) * zr + pt.x, + .y = sin(phi)*sin(theta) * zr + pt.y, + .z = cos(theta) * zr + pt.z + }; + } + + return 0.5*log(r)*r/dr; +} + +struct setup { + struct point p0; // origin + struct point direction; // ray direction + struct point x; // ray movement in col + struct point y; // ray movement in row +}; + +struct setup make_setup() +{ + // set up camera + struct camera cam; + cam.fov = 90; + + camera_set_looking_at(&cam, (struct point){.x=CAM_POSITION, .y= CAM_POSITION, .z = CAM_POSITION}, PT_ZERO); + + struct point span_z, span_xy; + + // get rotation axis + pt_orthogonal_plane(cam.direction, &span_z, &span_xy); + + printf("rendering %ix%ipx\n", SIZE, SIZE); + + // distance each ray has from anothe on the ortogonal plane + //DTYPE step_dist = 2 / (DTYPE) (width - 1); + + // vectors to move on the projection plane + struct point move_right = pt_scale(span_xy, STEP_WIDTH); + struct point move_up = pt_scale(span_z, STEP_WIDTH);; + + // set starting point + struct point starting_point = pt_normalize(cam.direction); + + // rotate starting point to (0,0) + starting_point = pt_add(starting_point, pt_mult(move_right, - SIZE / (DTYPE) 2)); + starting_point = pt_add(starting_point, pt_mult(move_up, - SIZE / (DTYPE) 2)); + + return (struct setup) { + cam.location, starting_point, span_xy, span_z + }; +} + +int main() +{ + struct setup setup = make_setup(); + struct point start = setup.p0; + + //printf("device num acc_device_current: %d\n", acc_get_num_devices(acc_device_current)); + //printf("device num acc_device_none: %d\n", acc_get_num_devices(acc_device_none)); + //printf("device num acc_device_default: %d\n", acc_get_num_devices(acc_device_default)); + //printf("device num acc_device_host: %d\n", acc_get_num_devices(acc_device_host)); + //printf("device num acc_device_not_host: %d\n", acc_get_num_devices(acc_device_not_host)); + //printf("device num acc_device_nvidia: %d\n", acc_get_num_devices(acc_device_nvidia)); + //printf("device num acc_device_radeon: %d\n", acc_get_num_devices(acc_device_radeon)); + + // get backing buff + int* buff = calloc(sizeof(int), SIZE * SIZE); + + // indicate malloc failure + if (buff == NULL) + return -255; + + printf("Before kernel\n"); + +#pragma acc data copy(buff) +{ + // kernel goes brr + #pragma acc kernels + for (int x = 0; x < SIZE; x++) { + for (int y = 0; y < SIZE; y++) { + // get direction + struct point offset = pt_add(pt_mult(setup.x, STEP_WIDTH * x), pt_mult(setup.y, STEP_WIDTH * y)); + struct point direction = pt_add(setup.direction, offset); + // get start + struct point loc = start; + // march! + DTYPE dist; + int res = -1; + + for (int i = 0; i < STEPS; i++) { + dist = mandelbulb_dist(loc); + if (dist < THRESHOLD) { + res = i; + break; + } + if (dist > 100) { + break; + } + loc = pt_add(loc, pt_scale(direction, dist)); + } + buff[y * SIZE + x] = res; + } + } +} + + printf("after kernel\n"); + + // convert distance field into image + + Image img; + + // initialize shared pixel buffer + image_new(SIZE, SIZE, &img); + + #pragma acc parallel + for (unsigned int x = 0; x < SIZE; x++) { + #pragma acc loop + for (unsigned int y = 0; y < SIZE; y++) { + if (buff[y * SIZE + x] < 0) { + image_set_px(img, x, y, 0, 0, 0); + } else { + // float in range of [0,1] + DTYPE fac = buff[y * SIZE + x] / (DTYPE) STEPS; + // calc shade + int shade = ((1-fac) * 255); + image_set_px(img, x, y, shade, shade, shade); + } + } + } + + printf("after encoding\n"); + + + image_save_bmp(img, "gpu-goes-brrrrrrr.bmp"); + + return 0; +} \ No newline at end of file diff --git a/images b/images index e20f4cb..f503f18 160000 --- a/images +++ b/images @@ -1 +1 @@ -Subproject commit e20f4cb891ae1f6ac34f7572ea6847957b7195db +Subproject commit f503f1835fd47de323b19d8f3f1ea0fc542dd6f6 diff --git a/main.c b/marcher.c similarity index 55% rename from main.c rename to marcher.c index 685501d..4a786ae 100644 --- a/main.c +++ b/marcher.c @@ -1,6 +1,10 @@ #include -#include "images/images.h" -#include "marcher.h" +#include +#include +#include "images/src/images.h" +#include "src/scene.h" +#include "src/camera.h" +#include "src/point.h" #define SCENE_MOD 2 @@ -15,15 +19,15 @@ */ -double circle_dist(Point x, SceneObject *self) { +double circle_dist(struct point x, struct scene_object *self) { double r = self->args[0]; return pt_dist(pt_mod(x, SCENE_MOD), self->location) - r; } -Color circle_color(Point hit, Point direction, SceneObject *self) { - Point obj_direction = self->location; +Color circle_color(struct point hit, struct point direction, struct scene_object *self) { + struct point obj_direction = self->location; - pt_sub(&obj_direction, pt_mod(hit, SCENE_MOD)); + obj_direction = pt_sub(obj_direction, pt_mod(hit, SCENE_MOD)); double angle = pt_angle(direction, obj_direction) / M_PI * 180; Color color = self->color; @@ -34,15 +38,16 @@ Color circle_color(Point hit, Point direction, SceneObject *self) { return color; } // constructs the scene object -SceneObject circle_new(Point loc, double radius) { - SceneObject so; - so.location = loc; - so.args = malloc(sizeof(double) * 2); - so.args[0] = radius; - so.distance = circle_dist; - so.get_color = circle_color; - so.color = color_new(255,255,255); - return so; +struct scene_object circle_new(struct point loc, double radius) { + double * args = malloc(sizeof (double) * 2); + args[0] = radius; + return (struct scene_object) { + .location = loc, + .args = args, + .distance = circle_dist, + .get_color = circle_color, + .color = color_new(255, 255, 255), + }; } @@ -55,11 +60,11 @@ SceneObject circle_new(Point loc, double radius) { Color function is just a flat shader, detail is displayed with ambient occlusion */ -double mandelbulb_dist(Point pt, SceneObject *self) { +double mandelbulb_dist(struct point pt, struct scene_object *self) { int iters = self->args[0]; double power = self->args[1]; - Point z = pt; + struct point z = pt; float dr = 1.0; float r = 0.0; for (int i = 0; i < iters ; i++) { @@ -79,39 +84,43 @@ double mandelbulb_dist(Point pt, SceneObject *self) { theta = theta*power; phi = phi*power; - // convert back to cartesian coordinates - z = pt_mult(pt_new(sin(theta)*cos(phi), sin(phi)*sin(theta), cos(theta)), zr); - pt_add(&z, pt); + // convert back to cartesian coordinates, add zr and the old pt + z = (struct point) { + .x = sin(theta)*cos(phi) * zr + pt.x, + .y = sin(phi)*sin(theta) * zr + pt.y, + .z = cos(theta) * zr + pt.z + }; } return 0.5*log(r)*r/dr; - } -Color mandelbulb_color(Point hit, Point direction, SceneObject *self) { +Color mandelbulb_color(struct point hit, struct point direction, struct scene_object *self) { return self->color; } // constructs the scene object -SceneObject mandelbulb_new(Point location, int iters, double power) { - SceneObject so; - so.location = location; - so.args = malloc(sizeof(double) * 3); - so.args[0] = iters; // iterations - so.args[1] = power; // power - so.args[2] = -1; // reserved for color calculations - so.distance = mandelbulb_dist; - so.get_color = mandelbulb_color; - so.color = color_new(255,255,255); - return so; +struct scene_object mandelbulb_new(struct point location, int iters, double power) { + double * args = malloc(sizeof(double) * 3); + args[0] = iters; + args[1] = power; + args[2] = -1; + + return (struct scene_object) { + .location= location, + .args= args, + .distance= mandelbulb_dist, + .get_color= mandelbulb_color, + .color= color_new(255,255,255), + }; } int main(int argc, char* argv[]) { - float dpi = 800; + float dpi = 200; int threads = 32; int size = dpi * 27.56f; // 400dpi by 70cm size float pow = 3; - float cam_position = 1.15; + float cam_position = 1.35; int steps = 1000; int iters = 500; float threshold = 0.001; @@ -131,15 +140,13 @@ int main(int argc, char* argv[]) { printf("Rendering to %s\n", path); - Camera cam; + struct camera cam; cam.fov = 90; - camera_set_looking_at(&cam, pt_new(cam_position, cam_position, cam_position), pt_new(0,0,0)); - - + camera_set_looking_at(&cam, (struct point){.x=cam_position, .y= 0, .z = cam_position}, PT_ZERO); // create basic scene with up to 10 objects - Scene scene = scene_new(size, size, 1); + struct scene scene = scene_new(size, size, 1); scene.perf_opts.max_steps = steps; scene.perf_opts.threshold = threshold; scene.perf_opts.speed_cutoff = 10; @@ -147,13 +154,13 @@ int main(int argc, char* argv[]) { //scene_add_obj(&scene, circle_new(pt_new(SCENE_MOD / 2.0, SCENE_MOD/ 2.0, SCENE_MOD / 2.0), .2)); - scene_add_obj(&scene, mandelbulb_new(pt_new(0,0,0), iters, pow)); + scene_add_obj(&scene, mandelbulb_new(PT_ZERO, iters, pow)); Image *img = render_scene(&scene, &cam, threads); image_save_bmp(*img, path); - image_destroy_shared(*img); + image_destroy(*img); scene_destroy(scene); return 0; diff --git a/marcher.h b/marcher.h deleted file mode 100644 index acc539a..0000000 --- a/marcher.h +++ /dev/null @@ -1,97 +0,0 @@ -#ifndef __MARCHER_H__ -#define __MARCHER_H__ - -#include "images/images.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; - unsigned int fov; -} Camera; - -Camera camera_new(Point direction, unsigned int fov); -void camera_set_looking_at(Camera *cam, Point origin, Point thing); - -// Scene objects have a position, some args, and a distance calculation function -// the distance calc function has the following signature: -// double distanceTo(Point myLocation, double * myArgs, Point externalPoint) -// where myLocation is this.location, myArgs is this.args and externalPoint is the point from wich we want to know the distance -// the get_color function takes args: point_hit, direction_hit, myArgs, MyLocation, MyColor -typedef struct __myobject { - Point location; - double * args; - double (*distance)(Point, struct __myobject *); - Color (*get_color)(Point, Point, struct __myobject *); - Color color; - struct __myscene* scene; -} SceneObject; - - -typedef struct __perfopts { - int speed_cutoff; - int max_steps; - double threshold; -} PerformanceOptimizations; - - -typedef struct __myscene { - unsigned int width; - unsigned int height; - SceneObject * objects; - int object_count; - int allocated_space; - // performance opts - PerformanceOptimizations perf_opts; - // colors etc - Color background; -} Scene; - -Image* render_scene(Scene *scene, Camera *camera, unsigned int threads); - -Scene scene_new(unsigned int width, unsigned int height, int obj_count); - -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" - -#endif \ No newline at end of file diff --git a/src/camera.c b/src/camera.c index fd205e4..7bcfe6e 100644 --- a/src/camera.c +++ b/src/camera.c @@ -1,12 +1,32 @@ -#include "../marcher.h" -#include #include #include - - -Camera camera_new(Point direction, unsigned int fov) { - Camera camera; - camera.location = pt_new(0,0,0); +#include +#include +#include +#include "point.h" +#include "camera.h" + +struct thread_args { + struct point start; + int thread_id; + int thread_count; + int height; + int width; + struct point move_up; + struct point move_right; + void (*callback)(struct point, int, int); +}; + +static void * camera_iterate_rays_const_dist_thread(void* args); + + +struct camera camera_new(struct point direction, unsigned int fov) { + struct camera camera; + camera.location = (struct point) { + .x = 0, + .y = 0, + .z = 0 + }; camera.fov = fov; // normalize camera direction @@ -15,160 +35,118 @@ Camera camera_new(Point direction, unsigned int fov) { return camera; } -void camera_set_looking_at(Camera *cam, Point origin, Point thing) { +void camera_set_looking_at(struct camera *cam, struct point origin, struct point thing) { cam->location = origin; - pt_sub(&thing, origin); - cam->direction = pt_normalize(thing); + cam->direction = pt_normalize(pt_sub(thing, origin)); } -void camera_iterate_rays_const_angle(Camera camera, int width, int height, int threads, void (*callback)(Point, int, int)) { - // negative threads => single threaded. - if (threads < 0) threads = 0; +void camera_iterate_rays_const_dist(struct camera camera, int width, int height, int thread_count, void (*callback)(struct point, int, int)) { + // negative thread_count => single threaded. + if (thread_count < 0) thread_count = 0; - Point span_z, span_xy; + struct point span_z, span_xy; // get rotation axis pt_orthogonal_plane(camera.direction, &span_z, &span_xy); - printf("rendering %ix%i px", width, height); - - pt_print_n("span_xy", span_xy); - pt_print_n("span_z", span_z); - - // angle between rays - double angle_step = camera.fov / (double) (width - 1); - - // rotation applied to reach the outmost end of the view - double angle_start_h = - (camera.fov / 2.0); - double angle_start_v = ((angle_step * (height - 1)) / 2) ; - - printf("step: %f\nstart_h: %f\nstart_v: %f\n", angle_step, angle_start_h, angle_start_v); + printf("rendering %ix%ipx\n", width, height); - // calculate both rotation matrices (expensive!) - Matrix rot_z = mtrx_rotation(span_z, angle_step); - Matrix rot_xy = mtrx_rotation(span_xy, -angle_step); + // distance each ray has from anothe on the ortogonal plane + double step_dist = 2 / (double) (width - 1); - // rotate vector to starting location (bot left of screen) - // (very expensive!) - Point starting_point = mtrx_mult( - mtrx_rotation(span_xy, angle_start_v), - mtrx_mult( - mtrx_rotation(span_z, angle_start_h), - camera.direction - ) - ); + // vectors to move on the projection plane + struct point move_right = pt_scale(span_xy, step_dist); + struct point move_up = pt_scale(span_z, step_dist);; - // initialize threads - int thread_id = 0; - for (int i = 0; i < threads - 1; i++) { - if (fork() == 0) { - thread_id = i + 1; - break; - } + // set starting point + struct point starting_point = pt_normalize(camera.direction); + + // rotate starting point to (0,0) + starting_point = pt_add(starting_point, pt_mult(move_right, - width / (double) 2)); + starting_point = pt_add(starting_point, pt_mult(move_up, - height / (double) 2)); + + if (thread_count < 2) { + //TODO implement single threaded work here + struct thread_args arg = { + .start = starting_point, + .thread_id = 0, + .thread_count = 1, + .height = height, + .width = width, + .move_up = move_up, + .move_right = move_right, + .callback = callback + }; + + camera_iterate_rays_const_dist_thread(&arg); + + return; } - // this point is rotated for every pixel - Point curr_pt = starting_point; - // (0,0) screenspace is bottom left corner - for (int y = 0; y < height; y++) { - curr_pt = mtrx_mult(rot_xy, starting_point); - // move starting point one row down - starting_point = curr_pt; - - if (y % threads != thread_id) continue; - - for (int x = 0; x < width; x++) { - callback(curr_pt, x, y); - curr_pt = mtrx_mult(rot_z, curr_pt); // rotate point - } + // initialize thread_count + pthread_t * threads = malloc(sizeof(pthread_t) * thread_count); + struct thread_args* args = malloc(sizeof(struct thread_args) * thread_count); + + for (int i = 0; i < thread_count; i++) { + args[i] = (struct thread_args) { + .start = starting_point, + .thread_id = i, + .thread_count = thread_count, + .height = height, + .width = width, + .move_up = move_up, + .move_right = move_right, + .callback = callback + }; + + pthread_create(threads + i, NULL, camera_iterate_rays_const_dist_thread, (void*) (args + i)); } - if (thread_id != 0) { - printf("Thread %i is finished\n", thread_id); - exit(0); - } + //struct thread_args { + // struct point start; + // int thread_id; + // int thread_count; + // int height; + // int width; + // struct point move_up; + // struct point move_right; + // void (*callback)(struct point, int, int); + //}; - int status; - for (int i = 0; i < threads - 1; i++) { - while(wait(&status) > 0) {} - } - printf("got threads\n"); + for (int i = 0; i < thread_count; i++) { + pthread_join(threads[i], NULL); + } } -void camera_iterate_rays_const_dist(Camera camera, int width, int height, int threads, void (*callback)(Point, int, int)) { - // negative threads => single threaded. - if (threads < 0) threads = 0; - - Point span_z, span_xy; - - // get rotation axis - pt_orthogonal_plane(camera.direction, &span_z, &span_xy); - - printf("rendering %ix%i px\n", width, height); - - pt_print_n("span_xy", span_xy); - pt_print_n("span_z", span_z); - - - // distance each ray has from anothe on the ortogonal plane - double step_dist = 2 / (double) (width - 1); - - // vectors to move on the projection plane - Point move_right = pt_scale(span_xy, step_dist); - Point move_up = pt_scale(span_z, step_dist);; - - printf("step: %f\n", step_dist); - // set starting point - Point starting_point = pt_normalize(camera.direction); - - // rotate starting point to (0,0) - pt_add(&starting_point, pt_mult(move_right, - width / (double) 2)); - pt_add(&starting_point, pt_mult(move_up, - height / (double) 2)); - - // initialize threads - int thread_id = 0; - for (int i = 0; i < threads - 1; i++) { - if (fork() == 0) { - thread_id = i + 1; - break; - } - } +static void * camera_iterate_rays_const_dist_thread(void* arg_ptr) { + // explicit cast to make gcc happy + struct thread_args* args = (struct thread_args*) arg_ptr; // this point is moved for every pixel - Point curr_pt = starting_point; + struct point curr_pt = args->start; // (0,0) screenspace is bottom left corner - for (int y = 0; y < height; y++) { + for (int y = 0; y < args->height; y++) { // move one row up (this has to be done in every thread!) - pt_add(&starting_point, move_up); + args->start = pt_add(args->start, args->move_up); // only render the lines this thread is responsible for - if (y % threads != thread_id) continue; - + if (y % args->thread_count != args->thread_id) continue; + // display progress in percent - if (height > 200 && y % (height / 100) == 0 && y != 0) { - printf("\r%02i%%", (y * 100) / height); + if (args->height > 200 && y % (args->height / 100) == 0 && y != 0) { + printf("\r%02i%%", (y * 100) / args->height); fflush(stdout); } // actually iterate this line - curr_pt = starting_point; - for (int x = 0; x < width; x++) { - callback(curr_pt, x, y); - pt_add(&curr_pt, move_right); // move pt right to next pt + curr_pt = args->start; + for (int x = 0; x < args->width; x++) { + args->callback(curr_pt, x, y); + curr_pt = pt_add(curr_pt, args->move_right); // move pt right to next pt } } - if (thread_id != 0) { - printf("Thread %i is finished\n", thread_id); - exit(0); - } - - int status; - for (int i = 0; i < threads - 1; i++) { - while(wait(&status) > 0) {} - } - - printf("got threads\n"); -} \ No newline at end of file + return NULL; +} diff --git a/src/camera.h b/src/camera.h new file mode 100644 index 0000000..2fb6318 --- /dev/null +++ b/src/camera.h @@ -0,0 +1,11 @@ +#pragma once + +struct camera { + struct point location; + struct point direction; + unsigned int fov; +}; + +struct camera camera_new(struct point direction, unsigned int fov); +void camera_set_looking_at(struct camera *cam, struct point origin, struct point thing); +void camera_iterate_rays_const_dist(struct camera camera, int width, int height, int thread_count, void (*callback)(struct point, int, int)); \ No newline at end of file diff --git a/src/point.c b/src/point.c deleted file mode 100644 index 1324312..0000000 --- a/src/point.c +++ /dev/null @@ -1,180 +0,0 @@ -#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; -} - -// 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 Point pt_mult(Point pt, double scalar) { - return pt_new( - pt.x * scalar, - pt.y * scalar, - pt.z * scalar - ); -} - -// return internal angle between a and b -inline double pt_angle(Point a, Point b) { - 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)); -} - -// 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; -} - -// 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); -} - -// 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); -} - -// 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) { - pt = pt_normalize(pt); - - // 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) - - // 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..a1f4b7d --- /dev/null +++ b/src/point.h @@ -0,0 +1,149 @@ +#pragma once + +#ifndef POINT_DTYPE +#define POINT_DTYPE double +#endif + + +struct point { + POINT_DTYPE x; + POINT_DTYPE y; + POINT_DTYPE z; +}; + +#define PT_ZERO ((struct point) {.x=0, .y=0, .z=0}) +#define PT_NEW(x,y,z) ((struct point){(x), (y), (z)}) + +static struct point pt_scale(struct point pt, POINT_DTYPE length); +static struct point pt_normalize(struct point pt); +static struct point pt_mult(struct point pt, POINT_DTYPE scalar); +static POINT_DTYPE pt_length(struct point pt); +static struct point pt_add(struct point pt, struct point add); +static struct point pt_sub(struct point pt, struct point sub); +static POINT_DTYPE pt_dist(struct point p1, struct point p2); +static struct point pt_mod(struct point pt, POINT_DTYPE mod); +static POINT_DTYPE pt_dot(struct point a, struct point b); +static struct point pt_cross(struct point a, struct point b); +static POINT_DTYPE pt_angle(struct point a, struct point b); +static void pt_print(struct point pt); +static void pt_print_n(const char* name, struct point pt); +static void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy); + +#include +#include +#include "point.h" + +// get the length of vector +static inline POINT_DTYPE pt_length_inline (struct point pt) __attribute__((always_inline)); + +static inline struct point pt_mult(struct point pt, POINT_DTYPE scalar) { + return (struct point) { + .x = pt.x * scalar, + .y = pt.y * scalar, + .z = pt.z * scalar + }; +} + +// return internal angle between a and b +static inline POINT_DTYPE pt_angle(struct point a, struct point b) { + return acos(pt_dot( + pt_normalize(a), + pt_normalize(b) + )); +} + +// get the length of vector +static inline POINT_DTYPE pt_length_inline (struct point pt) { + return sqrt((pt.x * pt.x) + (pt.y * pt.y) + (pt.z * pt.z)); +} + +static inline POINT_DTYPE pt_length(struct point pt) { + return pt_length_inline(pt); +} + +// add the vector add to the vector pt +static inline struct point pt_add(struct point pt, struct point add) { + return (struct point) { + .x = pt.x + add.x, + .y = pt.y + add.y, + .z = pt.z + add.z, + }; +} + +// add the vector add to the vector pt +static inline struct point pt_sub(struct point pt, struct point sub) { + return (struct point) { + .x = pt.x - sub.x, + .y = pt.y - sub.y, + .z = pt.z - sub.z, + }; +} + +static inline POINT_DTYPE pt_dist(struct point p1, struct point p2) { + return pt_length_inline(pt_sub(p1, p2)); +} + +// normalize a vector +static inline struct point pt_normalize(struct point pt) { + POINT_DTYPE length = pt_length(pt); + + return (struct point) { + .x = pt.x / length, + .y = pt.y / length, + .z = pt.z / length + }; +} + +// scale vector to length +static inline struct point pt_scale(struct point pt, POINT_DTYPE length) { + POINT_DTYPE f = length / pt_length_inline(pt); + return (struct point) { + .x = pt.x * f, + .y = pt.y * f, + .z = pt.z * f + }; +} + +// dot product of two vectors +static inline POINT_DTYPE pt_dot(struct point a, struct point b) { + return a.x*b.x + a.y*b.y + a.z*b.z; +} + +// cross product of two vectors +static inline struct point pt_cross(struct point a, struct point b) { + return (struct point) { + .x = a.y*b.z - a.z*b.y, + .y = a.z*b.x - a.x*b.z, + .z = a.x*b.y - a.y*b.x + }; +} + +static inline void pt_print(struct point pt) { + printf("(%f, %f, %f)\n", pt.x, pt.y, pt.z); +} + +static inline void pt_print_n(const char* name, struct 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" +static inline void pt_orthogonal_plane(struct point pt, struct point *span_z, struct point *span_xy) { + pt = pt_normalize(pt); + + // 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) + + // now use this, to find the vector + *span_z = pt_normalize(pt_cross(pt, *span_xy)); +} + +static inline struct point pt_mod(struct point pt, POINT_DTYPE mod) { + return (struct point) { + .x = fabs(fmod(pt.x, mod)), + .y = fabs(fmod(pt.y, mod)), + .z = fabs(fmod(pt.z, mod)) + }; +} \ No newline at end of file diff --git a/src/scene.c b/src/scene.c index 6b08202..b62fab6 100644 --- a/src/scene.c +++ b/src/scene.c @@ -1,36 +1,37 @@ -#include "../marcher.h" -#include "../images/images.h" - +#include "../images/src/images.h" +#include "scene.h" +#include "camera.h" +#include "point.h" #include +#include static Image* current_image; -static Scene* current_scene; -static Camera* current_camera; +static struct scene* current_scene; +static struct camera* current_camera; -Color march_ray(Point origin, Point direction, Scene* scene); -void camera_iter_callback(Point direction, int x, int y); +Color march_ray(struct point origin, struct point direction, struct scene* scene); +void camera_iter_callback(struct point direction, int x, int y); -Scene scene_new(unsigned int width, unsigned int height, int obj_count) { - Scene scene; +struct scene scene_new(unsigned int width, unsigned int height, int obj_count) { + struct scene scene; scene.height = height; scene.width = width; scene.object_count = 0; - scene.objects = malloc(obj_count * sizeof(SceneObject)); + scene.objects = malloc(obj_count * sizeof(struct scene_object)); scene.allocated_space = obj_count; scene.background = color_new(0,0,0); - PerformanceOptimizations perf_opts; - perf_opts.speed_cutoff = 0; - perf_opts.max_steps = 32; - perf_opts.threshold = 0.02; - - scene.perf_opts = perf_opts; + scene.perf_opts = (struct performance_optimization) { + .speed_cutoff = 0, + .max_steps = 32, + .threshold = 0.02 + }; return scene; } -void scene_add_obj(Scene* scene, SceneObject object) { +void scene_add_obj(struct scene * scene, struct scene_object object) { if (scene->object_count >= scene->allocated_space) return; // limit reached // TODO realloc @@ -44,13 +45,13 @@ void scene_add_obj(Scene* scene, SceneObject object) { // render out the scene with threads // creates a shared image, so destroy with image_destroy_shared then free struct with free_shared_memory -Image* render_scene(Scene *scene, Camera *camera, unsigned int threads) { +Image* render_scene(struct scene *scene, struct camera *camera, unsigned int threads) { current_image = malloc(sizeof(Image)); current_scene = scene; current_camera= camera; // initialize shared pixel buffer - image_new_shared(scene->width, scene->height, current_image); + image_new(scene->width, scene->height, current_image); // iterate over the rays camera_iterate_rays_const_dist(*camera, scene->width, scene->height, threads, camera_iter_callback); @@ -61,20 +62,20 @@ Image* render_scene(Scene *scene, Camera *camera, unsigned int threads) { } // march the ray, set the color. repeated for each direction generated by the camera -void camera_iter_callback(Point direction, int x, int y) { +void camera_iter_callback(struct point direction, int x, int y) { Color c = march_ray(current_camera->location, direction, current_scene); image_set_px_c(*current_image, x, y, c); } -Color march_ray(Point origin, Point direction, Scene* scene) { +Color march_ray(struct point origin, struct point direction, struct scene* scene) { // some local variables - Point pos = origin; + struct point pos = origin; double closest_encounter = DBL_MAX; double dist = closest_encounter; // the closest object we have - SceneObject* closest_obj = scene->objects; + struct scene_object* closest_obj = scene->objects; // get steps, threshold from scene int steps = scene->perf_opts.max_steps; @@ -88,7 +89,7 @@ Color march_ray(Point origin, Point direction, Scene* scene) { // find distance to closest object for(int i = 0; i < scene->object_count; i++) { // get pointer to scene obj - SceneObject* obj = scene->objects + i; + struct scene_object* obj = scene->objects + i; double curr_dist = scene->objects[i].distance(pos, obj); // if we are close @@ -106,8 +107,8 @@ Color march_ray(Point origin, Point direction, Scene* scene) { if (dist < closest_encounter) closest_encounter = dist; // scale direction vector to distance, then add it to our position - Point step_vector = pt_scale(direction, dist); - pt_add(&pos, step_vector); + struct point step_vector = pt_scale(direction, dist); + pos = pt_add(pos, step_vector); // one step taken... steps--; @@ -137,7 +138,7 @@ Color march_ray(Point origin, Point direction, Scene* scene) { } } -void scene_destroy(Scene scene) { +void scene_destroy(struct scene scene) { for (int i = 0; i < scene.object_count; i++) { // free args memory free(scene.objects[i].args); diff --git a/src/scene.h b/src/scene.h new file mode 100644 index 0000000..eacad6c --- /dev/null +++ b/src/scene.h @@ -0,0 +1,52 @@ +#pragma once + +#include "point.h" +#include "camera.h" + +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + +struct scene; + +// Scene objects have a position, some args, and a distance calculation function +// the distance calc function has the following signature: +// double distanceTo(struct point myLocation, double * myArgs, struct point externalPoint) +// where myLocation is this.location, myArgs is this.args and externalPoint is the point from wich we want to know the distance +// the get_color function takes args: point_hit, direction_hit, myArgs, MyLocation, MyColor +struct scene_object { + struct point location; + double * args; + double (*distance)(struct point, struct scene_object *); + Color (*get_color)(struct point, struct point, struct scene_object *); + Color color; + struct scene* scene; +}; + + +struct performance_optimization { + int speed_cutoff; + int max_steps; + double threshold; +}; + + +struct scene { + unsigned int width; + unsigned int height; + struct scene_object * objects; + int object_count; + int allocated_space; + // performance opts + struct performance_optimization perf_opts; + // colors etc + Color background; +}; + +Image* render_scene(struct scene *scene, struct camera *camera, unsigned int threads); + +struct scene scene_new(unsigned int width, unsigned int height, int obj_count); + +void scene_add_obj(struct scene* scene, struct scene_object object); + +void scene_destroy(struct scene scene); \ No newline at end of file