raymarcher.h/gpu.c

// use floats instead of doubles
#define DTYPE float
// #define POINT_DTYPE DTYPE

#include <stdlib.h>
#include <stdio.h>

#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 <openacc.h>


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;
}