diff --git a/bin/defocus.c b/bin/defocus.c index 0641921..344e23c 100644 --- a/bin/defocus.c +++ b/bin/defocus.c @@ -1,259 +1,63 @@ -#include -#include -#include +/* Main command line executable. */ + #include -#include +#define STB_IMAGE_WRITE_IMPLEMENTATION +#include -int pinhole_fn(int argc, char **argv); -int thin_lense_fn(int argc, char **argv); -int test_fn(int argc, char **argv); +int main(int argc, char **argv) { + df_sphere spheres[2]; + spheres[0].center_x = 0.f; + spheres[0].center_y = 0.f; + spheres[0].center_z = -1.f; + spheres[0].radius = .5f; + spheres[1].center_x = 0.25f; + spheres[1].center_y = 0.25f; + spheres[1].center_z = -1.5f; + spheres[1].radius = .75f; -static const char *_model_names[] = {"pinhole", "thin_lense", "test"}; -typedef int (*model_fn)(int argc, char **argv); + int image_width = 1024; + int image_height = 512; -static model_fn _model_fns[] = { - pinhole_fn, - thin_lense_fn, - test_fn, -}; + df_image_handle input_image = df_load_image("../test_image.png", &image_width, &image_height); + float aspect = (float)image_width / (float)image_height; -void usage(const char *pname) -{ - printf("Usage: %s [model parameters]\n", pname); - printf("Available models:\n"); - for (int i = 0; i < DF_ARRAY_COUNT(_model_names); ++i) { - printf(" %s\n", _model_names[i]); - } -} + df_plane plane; + plane.base_x = 0.f; + plane.base_y = 0.f; + plane.base_z = -1.f; + plane.normal_x = 0.f; + plane.normal_y = 0.f; + plane.normal_z = 1.f; + plane.img_p0_x = -aspect; + plane.img_p0_y = -1.f; + plane.img_p0_z = plane.base_z; + plane.img_p1_x = aspect; + plane.img_p1_y = 1.f; + plane.img_p1_z = plane.base_z; + plane.img_ax0_x = 1.f; + plane.img_ax0_y = 0.f; + plane.img_ax0_z = 0.f; + plane.img_ax1_x = 0.f; + plane.img_ax1_y = 1.f; + plane.img_ax1_z = 0.f; + plane.image = input_image; -int main(int argc, char **argv) -{ - if (argc < 2) { - fprintf(stderr, "Missing model name!\n"); - usage(argv[0]); - return 1; - } - - model_fn fn = NULL; - const char *model = argv[1]; - for (int i = 0; i < DF_ARRAY_COUNT(_model_names); ++i) { - if (strcmp(model, _model_names[i]) == 0) { - fn = _model_fns[i]; - break; - } - } - if (!fn) { - fprintf(stderr, "Invalid model name!\n"); - usage(argv[0]); - return 1; - } - - return fn(argc - 2, &argv[2]); -} - -int pinhole_fn(int argc, char **argv) -{ - const char *in_path = NULL; - const char *out_path = "out.png"; - const char *focal_str = NULL; - const char *in_z_str = NULL; - const char *out_z_str = NULL; - - for (int i = 0; i < argc; ++i) { - if ((strcmp(argv[i], "-i") == 0 || strcmp(argv[i], "--input") == 0) && i < argc - 1) { - ++i; - in_path = argv[i]; - } else if ((strcmp(argv[i], "-o") == 0 || strcmp(argv[i], "--output") == 0) && i < argc - 1) { - ++i; - out_path = argv[i]; - } else if ((strcmp(argv[i], "-f") == 0 || strcmp(argv[i], "--focal-length") == 0) && i < argc - 1) { - ++i; - focal_str = argv[i]; - } else if ((strcmp(argv[i], "-iz") == 0 || strcmp(argv[i], "--input-z") == 0) && i < argc - 1) { - ++i; - in_z_str = argv[i]; - } else if ((strcmp(argv[i], "-oz") == 0 || strcmp(argv[i], "--output-z") == 0) && i < argc - 1) { - ++i; - out_z_str = argv[i]; - } else { - in_path = argv[i]; - } - } - - if (!in_path || !focal_str || !in_z_str || !out_z_str) { - fprintf(stderr, "Missing model parameters!\n"); - printf("Pinhole model parameters:\n"); - printf(" -i | --input \t\tPath to the input image.\n"); - printf(" -o | --output \t\tPath to the output image (Default: out.png).\n"); - printf(" -iz | --input-z \t\tDistance to the input image plane.\n"); - printf(" -oz | --output-z \t\tDistance to the output image plane.\n"); - printf(" -f | --focal-length \t\tThe focal length of the camera lens.\n"); - return 1; - } - - float focal, in_z, out_z; - focal = atof(focal_str); - in_z = atof(in_z_str); - out_z = atof(out_z_str); - - df_image *in_image, *out_image; - - int w, h; - df_result res = df_load_image(in_path, &w, &h, &in_image); - if (res != df_result_success) { - fprintf(stderr, "Failed to load %s! (%d)\n", in_path, (int)res); - return 1; - } - - res = df_create_image(w, h, &out_image); - if (res != df_result_success) { - fprintf(stderr, "Failed to create output image (%d)\n", (int)res); - df_release_image(in_image); - return 1; - } - - df_pinhole((df_pinhole_params){.orig_z = in_z, .new_z = out_z, .focal_length = focal}, in_image, out_image); - - int error_code = 0; - res = df_write_image(out_image, out_path); - if (res != df_result_success) { - fprintf(stderr, "Failed to write to output image %s (%d)\n", out_path, (int)res); - error_code = 1; - } - - df_release_image(in_image); - df_release_image(out_image); - - return error_code; -} - -int thin_lense_fn(int argc, char **argv) -{ - const char *in_path = NULL; - const char *out_path = "out.png"; - const char *focal_str = NULL; - const char *aperture_str = NULL; - const char *in_z_str = NULL; - const char *out_z_str = NULL; - const char *samples_str = "64"; - - for (int i = 0; i < argc; ++i) { - if ((strcmp(argv[i], "-i") == 0 || strcmp(argv[i], "--input") == 0) && i < argc - 1) { - ++i; - in_path = argv[i]; - } else if ((strcmp(argv[i], "-o") == 0 || strcmp(argv[i], "--output") == 0) && i < argc - 1) { - ++i; - out_path = argv[i]; - } else if ((strcmp(argv[i], "-f") == 0 || strcmp(argv[i], "--focal-length") == 0) && i < argc - 1) { - ++i; - focal_str = argv[i]; - } else if ((strcmp(argv[i], "-iz") == 0 || strcmp(argv[i], "--input-z") == 0) && i < argc - 1) { - ++i; - in_z_str = argv[i]; - } else if ((strcmp(argv[i], "-oz") == 0 || strcmp(argv[i], "--output-z") == 0) && i < argc - 1) { - ++i; - out_z_str = argv[i]; - } else if ((strcmp(argv[i], "-a") == 0 || strcmp(argv[i], "--aperture") == 0) && i < argc - 1) { - ++i; - aperture_str = argv[i]; - } else if ((strcmp(argv[i], "-s") == 0 || strcmp(argv[i], "--samples") == 0) && i < argc - 1) { - ++i; - samples_str = argv[i]; - } else { - in_path = argv[i]; - } - } - - if (!in_path || !focal_str || !in_z_str || !out_z_str) { - fprintf(stderr, "Missing model parameters!\n"); - printf("Pinhole model parameters:\n"); - printf(" -i | --input \t\tPath to the input image.\n"); - printf(" -o | --output \t\tPath to the output image (Default: out.png).\n"); - printf(" -iz | --input-z \t\tDistance to the input image plane.\n"); - printf(" -oz | --output-z \t\tDistance to the output image plane.\n"); - printf(" -f | --focal-length \t\tThe focal length of the camera lens.\n"); - printf(" -a | --aperture \t\tThe size of the camera lens.\n"); - printf(" -s | --samples \t\tThe number of samples per pixel (Default: 64).\n"); - return 1; - } - - float focal, in_z, out_z, aperture; - int samples; - focal = atof(focal_str); - in_z = atof(in_z_str); - out_z = atof(out_z_str); - aperture = atof(aperture_str); - samples = atoi(samples_str); - - df_image *in_image, *out_image; - - int w, h; - df_result res = df_load_image(in_path, &w, &h, &in_image); - if (res != df_result_success) { - fprintf(stderr, "Failed to load %s! (%d)\n", in_path, (int)res); - return 1; - } - - res = df_create_image(w, h, &out_image); - if (res != df_result_success) { - fprintf(stderr, "Failed to create output image (%d)\n", (int)res); - df_release_image(in_image); - return 1; - } - - df_thin_lense((df_thin_lense_params){.focal_distance = in_z, - .out_distance = out_z, - .focal_length = focal, - .aperture = aperture, - .sample_count = samples}, - in_image, - out_image); - - int error_code = 0; - res = df_write_image(out_image, out_path); - if (res != df_result_success) { - fprintf(stderr, "Failed to write to output image %s (%d)\n", out_path, (int)res); - error_code = 1; - } - - df_release_image(in_image); - df_release_image(out_image); - - return error_code; -} - -int test_fn( int argc, char **argv ) { - - df_v3 v = {1.f, 2.f, 3.f}; - df_m4 t = df_translate(-1.f, 2.f, 1.5f); - t = df_inverse_transform(t); - - df_v3 tv = df_transform_v3(t, v); - - df_camera_i cam = df_create_perspective_camera(1024.0, 1024.0, 1024.0, 1024.0, 2.e-4f, 1000.f, 350.0, 21.0); - df_ray_packet packet = cam.build_ray_packet(cam.o); - - int w = 1024; - int h = 1024; - df_image *out_image; - df_result res = df_create_image(w, h, &out_image); - if (res != df_result_success) { - fprintf(stderr, "Failed to create output image (%d)\n", (int)res); - return 1; - } + uint8_t *image = NULL; - df_evaluate_ray_packet(&packet, out_image); + df_trace_rays((df_trace_rays_settings){ + .focal_length = 1.f, + .image_width = image_width, + .image_height = image_height + + }, + spheres, 0, + &plane, 1, + &image); + + stbi_write_png("raytracer.png", image_width, image_height, 3, image, image_width * 3); - const char *out_path = "out.png"; - res = df_write_image(out_image, out_path); - if (res != df_result_success) { - fprintf(stderr, "Failed to write to output image %s (%d)\n", out_path, (int)res); - } - df_release_ray_packet(&packet); - - cam.release(cam.o); return 0; } \ No newline at end of file diff --git a/include/defocus/base.h b/include/defocus/base.h deleted file mode 100644 index dafec37..0000000 --- a/include/defocus/base.h +++ /dev/null @@ -1,278 +0,0 @@ -#ifndef DEFOCUS_BASE_H -#define DEFOCUS_BASE_H - -#include -#include - -#include - -/** @file base.h - * @brief basic utilities - */ - -/** @brief Result codes used throughout the codebase */ -typedef enum -{ - /** @brief operation successfull */ - df_result_success = 0, - - /** @brief memory allocation failed */ - df_result_out_of_memory = 1, - - /** @brief file i/o error */ - df_result_io_error = 2, -} df_result; - -/** @brief Used to silence warnings about unused variables */ -#define DF_UNUSED(x) ((void)sizeof((x))) - -/** @brief Evaluates to the number of elements in a static array. */ -#define DF_ARRAY_COUNT(a) (sizeof((a)) / sizeof((a)[0])) - -/** @brief Zero a memory range */ -#define DF_ZERO_MEMORY(ptr, n) memset((ptr), 0, (n)) - -/** @brief Zero a struct */ -#define DF_ZERO_STRUCT(ptr) DF_ZERO_MEMORY(ptr, sizeof(*(ptr))) - -/** @brief Zero an array */ -#define DF_ZERO_ARRAY(a, n) DF_ZERO_MEMORY(a, sizeof((a)[0]) * (n)) - -#ifdef _MSC_VER -#ifdef __cplusplus__ -#define DF_API extern "C" __declspec(dllexport) -#else -#define DF_API __declspec(dllexport) -#endif /* _MSC_VER && __cplusplus__ */ -#else -#ifdef __cplusplus__ -#define DF_API extern "C" -#else -#define DF_API -#endif -#endif - -/* Simple logging function */ - -#ifndef DF_ENABLE_LOGGING -#define DF_ENABLE_LOGGING 1 -#endif - -enum -{ - df_log_level_verbose = 0, - df_log_level_info = 1, - df_log_level_warn = 2, - df_log_level_error = 3, -}; - -#if DF_ENABLE_LOGGING - -/** @brief The actual log implementation. - * - * There should be no need to call this directly. - * - * Instead use one of: - * - df_log_verbose() - * - df_log_info() - * - df_log_warn() - * - df_log_error() - - * @param level the log level. Used to filter the messages. - * @param file file name of the log location - * @param line line number of the log location - * @param fmt format string of the message. - */ -DF_API void df_log_impl(int level, const char *file, int line, const char *fmt, ...); - -#define df_log_verbose(...) df_log_impl(df_log_level_verbose, __FILE__, __LINE__, __VA_ARGS__) -#define df_log_info(...) df_log_impl(df_log_level_info, __FILE__, __LINE__, __VA_ARGS__) -#define df_log_warn(...) df_log_impl(df_log_level_warn, __FILE__, __LINE__, __VA_ARGS__) -#define df_log_error(...) df_log_impl(df_log_level_error, __FILE__, __LINE__, __VA_ARGS__) - -#else /* DF_ENABLE_LOGGING = 0 */ - -#define df_log_verbose(fmt, ...) -#define df_log_info(fmt, ...) -#define df_log_warn(fmt, ...) -#define df_log_error(fmt, ...) - -#endif - -/* Math constants. - */ -#define DF_PI 3.1415926f -#define DF_PI_OVER_2 1.570796f -#define DF_PI_OVER_4 0.785398f - -/* Basic types */ - -/** @brief RGBA color */ -typedef union { - struct - { - uint8_t r; - uint8_t g; - uint8_t b; - uint8_t a; - }; - uint8_t e[4]; -} df_color; - -/** @brief Linear interpolation between two colors. - * - * @param a first color - * @param b second color - * @param t interpolation factor - * @return a + (b - a) * t - */ -DF_API df_color df_lerp_color(df_color a, df_color b, double t); - -/** @brief 2d vector */ -typedef union { - struct - { - float x; - float y; - }; - struct - { - float u; - float v; - }; - float e[2]; -} df_v2; - -/** @brief Add two 2d vectors */ -DF_API df_v2 df_add_v2(df_v2 a, df_v2 b); - -/** @brief Subtract two 2d vectors */ -DF_API df_v2 df_sub_v2(df_v2 a, df_v2 b); - -/** @brief Calculate the dot product of 2d vectors a and b. - */ -DF_API float df_dot_v2(df_v2 a, df_v2 b); - -/** @brief Multiply a 2d vector with a scalar. - */ -DF_API df_v2 df_mul_v2(float t, df_v2 v); - -/** @brief Returns the normalized version of a 2d vector v - */ -DF_API df_v2 df_normalize_v2(df_v2 v); - -/** @brief 3d vector */ -typedef union { - struct - { - float x; - float y; - float z; - }; - float e[3]; -} df_v3; - -/** @brief Add two 3d vectors */ -DF_API df_v3 df_add_v3(df_v3 a, df_v3 b); - -/** @brief Subtract two 3d vectors */ -DF_API df_v3 df_sub_v3(df_v3 a, df_v3 b); - -/** @brief Calculate the dot product of 3d vectors a and b. - */ -DF_API float df_dot_v3(df_v3 a, df_v3 b); - -/** @brief Multiply a 3d vector with a scalar. - */ -DF_API df_v3 df_mul_v3(float t, df_v3 v); - -/** @brief Returns the normalized version of a 3d vector v - */ -DF_API df_v3 df_normalize_v3(df_v3 v); - -/** @brief A plane in 3d space. - * - * Expressed as the set of all points p for which - * dot((p - base), normal) = 0 - */ -typedef struct -{ - df_v3 base; - df_v3 normal; -} df_plane; - -/** @brief A line in 3d space. - * - * Expressed as the vector equation p = base + t * direction. - */ -typedef struct -{ - df_v3 base; - df_v3 direction; -} df_line; - -/** @brief The three possible result types of a line-plane intersection test. - */ -typedef enum -{ - /** The line and plane intersect in one point */ - df_line_plane_intersection_vector, - - /** The line is contained in the plane */ - df_line_plane_intersection_contained, - - /** The line and plane don't intersect (-> they are parallel) */ - df_line_plane_intersection_none, -} df_line_plane_intersection_type; - -/** @brief Result of a line-plane intersection test. - * The contained vector is only valid if @c type = @c df_line_plane_intersection_vector. - */ -typedef struct -{ - df_line_plane_intersection_type type; - df_v3 vec; -} df_line_plane_intersection; - -/** @brief Calculate the intersection between a line and a plane in 3d space */ -DF_API df_line_plane_intersection df_calc_line_plane_intersection(df_line line, df_plane plane); - -/** @brief A 4x4 matrix. - * - * Stored in row major order. - */ -typedef struct df_m4 -{ - _Alignas(16) union { - float e[16]; - __m128 vec[4]; - }; -} df_m4; - -/** Accesses the element at the given row and column of a 4x4 matrix */ -#define DF_M4_AT(m, row, col) ((m).e[(row)*4 + (col)]) - -/** @brief Matrix multiply 4x4 matrix a and b */ -DF_API df_m4 df_mul_m4(df_m4 a, df_m4 b); - -/** @brief Get a scale matrix */ -DF_API df_m4 df_scale(float x, float y, float z); - -DF_API df_m4 df_translate(float x, float y, float z); - -/** @brief Transform (i.e. multiply) a 3d vector v by the transformation matrix T */ -DF_API df_v3 df_transform_v3(df_m4 T, df_v3 v); - -/** @brief Calculate the inverse of a non-scaling transform matrix. - * - * Special fast case. - */ -DF_API df_m4 df_inverse_transform_no_scale(df_m4 M); - -/** @brief Calculate the inverse of a transform matrix */ -DF_API df_m4 df_inverse_transform(df_m4 M); - -/** @brief Calculate the inverse of a general (invertible) 4x4 matrix */ -DF_API df_m4 df_inverse(df_m4 M); - -#endif \ No newline at end of file diff --git a/include/defocus/camera.h b/include/defocus/camera.h deleted file mode 100644 index 7980cd5..0000000 --- a/include/defocus/camera.h +++ /dev/null @@ -1,28 +0,0 @@ -#ifndef DEFOCUS_CAMERA_H -#define DEFOCUS_CAMERA_H - -#include "base.h" -#include "raytracing.h" - -/** @file camera.h - * @brief basic camera functions - */ - -/** @brief Interface for cameras. */ -typedef struct -{ - /** Release the camera object */ - void (*release)(void *o); - - /** Builds a ray packet of all rays that need to be evaluated to generate the image. - */ - df_ray_packet (*build_ray_packet)(void *o); - - /** Opaque pointer to the camera object */ - void *o; -} df_camera_i; - -DF_API df_camera_i df_create_perspective_camera( - float image_width, float image_height, float raster_width, float raster_height, float far_dist, float near_dist, float focal_dist, float lens_radius); - -#endif \ No newline at end of file diff --git a/include/defocus/defocus.h b/include/defocus/defocus.h index 2efba8d..ce4e6aa 100644 --- a/include/defocus/defocus.h +++ b/include/defocus/defocus.h @@ -1,11 +1,94 @@ -#ifndef DEFOCUS_DEFOCUS_H -#define DEFOCUS_DEFOCUS_H +#ifndef DEFOCUS_H +#define DEFOCUS_H -/** @file defocus.h - * @brief meta header +#include + +#ifdef __cplusplus + +#ifdef _MSC_VER +#define DF_API extern "C" __declspec(dllexport) +#else +#define DF_API extern "C" +#endif + +#else + +#ifdef _MSC_VER +#define DF_API __declspec(dllexport) +#else +#define DF_API +#endif + +#endif /* __cplusplus */ + +#define DF_EPSF32 0.00001f + +#define DF_ARRAY_COUNT(A) (sizeof((A)) / sizeof((A)[0])) + +/* Images */ +typedef unsigned int df_image_handle; + +DF_API df_image_handle df_load_image(const char *path, int *w, int *h); + +/* Settings for the basic ray tracing function. */ +typedef struct +{ + /* Width needs to be divisible by 4! */ + int image_width; + int image_height; + + /* Distance between lense and image plane */ + float focal_length; + float lens_radius; +} df_trace_rays_settings; + +/* Sphere shape */ +typedef struct +{ + float center_x; + float center_y; + float center_z; + float radius; +} df_sphere; + +typedef struct +{ + float base_x; + float base_y; + float base_z; + float normal_x; + float normal_y; + float normal_z; + + /* Coordinates on the plane */ + float img_p0_x; + float img_p0_y; + float img_p0_z; + float img_p1_x; + float img_p1_y; + float img_p1_z; + + /* TODO(Kevin): These could be calculated from p0 and p1... */ + float img_ax0_x; + float img_ax0_y; + float img_ax0_z; + float img_ax1_x; + float img_ax1_y; + float img_ax1_z; + df_image_handle image; +} df_plane; + + +DF_API float df_max_f32(const float *list, unsigned int count); + +/* Core function, implements raytracing. + * + * The result buffer receives RGB pixels, with 8 bits per channel stored in RGB order. + * It is allocated via malloc(). + * + * Returns 1 on success, 0 if allocating the result image failed. */ -#include "base.h" -#include "image.h" -#include "models.h" +DF_API int df_trace_rays(df_trace_rays_settings settings, const df_sphere *spheres, unsigned int sphere_count, + const df_plane *planes, unsigned int plane_count, uint8_t **result); #endif diff --git a/include/defocus/image.h b/include/defocus/image.h deleted file mode 100644 index 72ce747..0000000 --- a/include/defocus/image.h +++ /dev/null @@ -1,62 +0,0 @@ -#ifndef DEFOCUS_IMAGE_H -#define DEFOCUS_IMAGE_H - -/** @file image.h - * @brief image type and functions - */ - -#include "base.h" - -/** @brief Opaque rgba8 image object */ -typedef struct df_image df_image; - -/** @brief create an image - * - * The contents of the image will be undefined. - * @param w the image width - * @param h the image height - * @param out_image receives the image object - * @return error code - */ -DF_API df_result df_create_image(int w, int h, df_image **out_image); - -/** @brief load an image file - * - * The image data will be converted to rgba8 - * - * @param w the image width - * @param h the image height - * @param out_image receives the image object - * @return error code - */ -DF_API df_result df_load_image(const char *path, int *out_w, int *out_h, df_image **out_image); - -/** @brief Write an image to a PNG file - * @param img the image - * @param path the path - */ -DF_API df_result df_write_image(df_image *img, const char *path); - -/** @brief Free an image. - * - * Any pointer to the image will be invalid after this. - * @param img the image - */ -DF_API void df_release_image(df_image *img); - -/** @brief Returns the dimensions of the image */ -DF_API void df_get_image_size(const df_image *image, int *w, int *h); - -/** @brief Returns the color value at pixel coordinates x, y - * - * Returns black for coordinates outside the image. - */ -DF_API df_color df_get_image_pixel(const df_image *image, int x, int y); - -/** @brief Set the color value at pixel coordinates x, y. - * - * Does nothing for coordinates outside the image. - */ -DF_API void df_set_image_pixel(df_image *image, int x, int y, df_color c); - -#endif diff --git a/include/defocus/intrinsic_helper.h b/include/defocus/intrinsic_helper.h deleted file mode 100644 index 996354b..0000000 --- a/include/defocus/intrinsic_helper.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef DEFOCUS_INTRINSIC_HELPER_H -#define DEFOCUS_INTRINSIC_HELPER_H - -/** @file intrinsic_helper.h - * @brief Helper macros for intrinsics - * - * Assumes that at least SS3 is available. - */ - -#include -#include - -/* Shuffle mask for _mm_shuffle_epi32 */ -#define DF_MAKE_SHUFFLE_MASK(x, y, z, w) (x | (y << 2) | (z << 4) | (w << 6)) - -#define DF_VEC_SWIZZLE_MASK(vec, mask) _mm_castsi128_ps(_mm_shuffle_epi32(_mm_castps_si128(vec), mask)) - -#define DF_VEC_SWIZZLE(vec, x, y, z, w) DF_VEC_SWIZZLE_MASK(vec, DF_MAKE_SHUFFLE_MASK(x, y, z, w)) - -// "broadcasts" element x [0, 1, 2, 3] to all elements of the result */ -#define DF_VEC_SWIZZLE1(vec, x) DF_VEC_SWIZZLE_MASK(vec, DF_MAKE_SHUFFLE_MASK(x, x, x, x)) - -/* returns [vec0, vec0, vec2, vec2] */ -#define DF_VEC_SWIZZLE_0022(vec) _mm_moveldup_ps(vec) - -/* returns [vec1, vec1, vec3, vec3] */ -#define DF_VEC_SWIZZLE_1133(vec) _mm_movehdup_ps(vec) - -/* returns [vec1[x], vec1[y], vec2[z], vec2[w]] */ -#define DF_VEC_SHUFFLE(vec1, vec2, x, y, z, w) _mm_shuffle_ps(vec1, vec2, DF_MAKE_SHUFFLE_MASK(x, y, z, w)) - -/* returns [vec1[0], vec1[1], vec2[0], vec2[1]] */ -#define DF_VEC_SHUFFLE_0101(vec1, vec2) _mm_movelh_ps(vec1, vec2) - -/* returns [vec1[2], vec1[3], vec2[2], vec2[3]] */ -#define DF_VEC_SHUFFLE_2323(vec1, vec2) _mm_movehl_ps(vec2, vec1) - -/* returns the n-th float inside the vector */ -inline float extract_float(__m128 v, unsigned int n) -{ - /* _mm_shuffle_ps expects a constant */ - switch (n) { - case 0: - return _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 0))); - case 1: - return _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 1))); - case 2: - return _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 2))); - case 3: - return _mm_cvtss_f32(_mm_shuffle_ps(v, v, _MM_SHUFFLE(0, 0, 0, 3))); - } - int nan = 0x7F800001; - return *(float *)&nan; -} - -#endif diff --git a/include/defocus/models.h b/include/defocus/models.h deleted file mode 100644 index b172c73..0000000 --- a/include/defocus/models.h +++ /dev/null @@ -1,77 +0,0 @@ -#ifndef DF_MODELS_H -#define DF_MODELS_H - -/** @file models.h - * @brief Camera models for defocus. - * - * Camera models are usually implemented as a function taking an input image and parameters required for the model - * and produce an output image. - * - * Parameters are contained in a parameter struct, because this enables us to something similar to named parameters. - * We can also do default parameters, by leaving some values set to zero. - * - * @code{c} - * df_thin_lense((df_thin_lense){ - * .focal_length = 42.f, - * .aperture = 21.f, - * .focal_distance = 150.f, - * .out_distance = 200.f}, - * in_image, - * out_image); - * @endcode - * - * All camera models use the same coordinate space: The camera looks along the positive z-axis in a right-handed - * coordinate system. (-> positive y is down, positive x is right.) - - */ - -#include "image.h" - -/** Parameters for the pinhole model. */ -typedef struct df_pinhole_params -{ - float focal_length; /**< the cameras focal length in millimeters. */ - float orig_z; /**< the distance of the input image from the camera. */ - float new_z; /**< the distance of the generated output image. */ -} df_pinhole_params; - -/** @brief Simple pinhole camera model. - * - * The simplest possible(?) camera model. - * This function takes an input image at a known distance from the camera and moves it to a new distance. - * - * @param params model parameters. - * @param in_image input image. - * @param out_image the output image. - */ -DF_API void df_pinhole(df_pinhole_params params, const df_image *in_image, df_image *out_image); - -/** Parameters for the thin lense model. - */ -typedef struct df_thin_lense_params -{ - /** The cameras focal length in millimeters. */ - float focal_length; - - /** The aperture size in millimeters */ - float aperture; - - /** The distance between the plane of focus and the camera lens, in millimeters */ - float focal_distance; - - /** The distance between the output image and the camera lens, in millimeters */ - float out_distance; - - /** The number of samples per pixel. Leave at 0 for the default value (64). */ - int sample_count; -} df_thin_lense_params; - -/** @brief Thin-lense model. - * - * @param params model parameters - * @param in_image input image - * @param out_image the output image. - */ -DF_API void df_thin_lense(df_thin_lense_params params, const df_image *in_image, df_image *out_image); - -#endif diff --git a/include/defocus/raytracing.h b/include/defocus/raytracing.h deleted file mode 100644 index 3a10d95..0000000 --- a/include/defocus/raytracing.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef DF_RAYTRACING_H -#define DF_RAYTRACING_H - -#include "image.h" - -/** @brief Stores information for rays in a SIMD friendly way */ -typedef struct -{ - /** Packs all SIMD data into a single buffer */ - float *simd_mem; - - /** Source uvs for rays */ - df_v2 *ray_uvs; - - /* Points into simd_mem */ - float *base_x; - float *base_y; - float *base_z; - float *dir_x; - float *dir_y; - float *dir_z; - - size_t ray_count; -} df_ray_packet; - -/** @brief Free ray packet memory */ -DF_API void df_release_ray_packet(df_ray_packet *rays); - -DF_API void df_evaluate_ray_packet(const df_ray_packet *rays, df_image *output); - -#endif diff --git a/include/defocus/scene.h b/include/defocus/scene.h deleted file mode 100644 index 89cbe02..0000000 --- a/include/defocus/scene.h +++ /dev/null @@ -1,37 +0,0 @@ -#ifndef DEFOCUS_SCENE_H -#define DEFOCUS_SCENE_H - -#include "base.h" -#include "image.h" -#include "camera.h" - -/** @file scene.h - * @brief scene structure */ - -/** Opaque scene structure */ -typedef struct df_scene df_scene; - -df_scene *df_create_scene(); - -void df_release_scene(df_scene *scene); - -void df_scene_add_plane(df_scene *scene, df_plane plane, df_image *texture); - -df_image *df_scene_generate_image(df_scene *scene, df_camera_i camera); - -#if 0 -/* Usage example */ - -df_camera_i camera = df_create_perspective_camera(...); -df_scene* scene = df_create_scene(); - -df_scene_add_plane(..., in_image); - -df_image *out_image = df_scene_generate_image(scene, camera); - -df_release_scene(scene); -camera->release(camera->o); - -#endif - -#endif diff --git a/lib/camera.c b/lib/camera.c deleted file mode 100644 index 6ac9cd9..0000000 --- a/lib/camera.c +++ /dev/null @@ -1,156 +0,0 @@ -#include - -#include -#include -#include - -#ifdef _MSC_VER -#include -#endif - -/* ******************************************************** - * - * Ray Packet Functions - * - * ********************************************************/ - -DF_API void df_release_ray_packet(df_ray_packet *rays) -{ -#ifdef _MSC_VER - _aligned_free(rays->simd_mem); -#else - free(rays->simd_mem); -#endif - free(rays->ray_uvs); - - memset(rays, 0, sizeof(*rays)); -} - -/* ******************************************************** - * - * Perspective Camrea - * - * ********************************************************/ - -typedef struct -{ - float focal_dist; - float lens_radius; - - /* Raster space, i.e. the space for which we generate rays */ - float raster_width; - float raster_height; - - /* Image(=screen) space, i.e. the space in which we store pixels */ - float image_width; - float image_height; - - df_m4 screen_to_raster; - df_m4 raster_to_screen; - df_m4 raster_to_camera; -} df_perspective_camera; - -static void pc_release(void *o) { df_perspective_camera *camera = o; } - -static df_ray_packet pc_build_ray_packet(void *o) -{ - df_perspective_camera *camera = o; - df_ray_packet packet; - memset(&packet, 0, sizeof(packet)); - - /* Generate 1 ray per pixel (ignore lens for now) */ - size_t count = (size_t)camera->raster_width * (size_t)camera->raster_height; - - size_t alloc_count = count; - /* Round up to nearest multiple of 4, for SSE. - * Also satisfies 16 byte alignment (assuming simd_mem is 16 byte aligned) - * because 4 * sizeof(float) = 16. - */ - if ((alloc_count % 4) != 0) { - alloc_count = ((alloc_count + 3) / 4) * 4; - } - -#ifdef _MSC_VER - packet.simd_mem = _aligned_malloc(sizeof(float) * alloc_count * 6, 16); -#elif defined(_ISOC11_SOURCE) /* Feature test macro for GCC */ - packet.simd_mem = aligned_alloc(16, sizeof(float) * alloc_count * 6); -#else - /* Fall back to regular malloc and hope for the best */ - packet.simd_mem = malloc(sizeof(float) * alloc_count * 6); -#endif - - packet.base_x = packet.simd_mem; - packet.base_y = packet.base_x + alloc_count; - packet.base_z = packet.base_y + alloc_count; - packet.dir_x = packet.base_z + alloc_count; - packet.dir_y = packet.dir_x + alloc_count; - packet.dir_z = packet.dir_y + alloc_count; - - packet.ray_uvs = malloc(sizeof(df_v2) * count); - - packet.ray_count = count; - - size_t i = 0; - for (float y = 0; y < camera->raster_height; y += 1.f) { - for (float x = 0; x < camera->raster_width; x += 1.f) { - packet.base_x[i] = 0.f; - packet.base_y[i] = 0.f; - packet.base_z[i] = 0.f; - - df_v3 raster_p = {x, y, 0.f}; - df_v3 camera_p = df_transform_v3(camera->raster_to_camera, raster_p); - df_v3 dir = df_normalize_v3(camera_p); - - packet.dir_x[i] = dir.x; - packet.dir_y[i] = dir.y; - packet.dir_z[i] = dir.z; - - df_v3 img_p = df_transform_v3(camera->raster_to_screen, raster_p); - packet.ray_uvs[i].u = img_p.x / camera->image_width; - packet.ray_uvs[i].v = img_p.y / camera->image_height; - - ++i; - assert(i <= count); - } - } - - return packet; -} - -DF_API df_camera_i df_create_perspective_camera(float image_width, - float image_height, - float raster_width, - float raster_height, - float far_dist, - float near_dist, - float focal_dist, - float lens_radius) -{ - df_perspective_camera *camera = malloc(sizeof(*camera)); - - camera->screen_to_raster = df_scale(raster_width, raster_height, 1.f); - camera->screen_to_raster = - df_mul_m4(camera->screen_to_raster, df_scale(1.f / image_width, -1.f / image_height, 1.f)); - camera->screen_to_raster = df_mul_m4(camera->screen_to_raster, df_translate(0.f, -image_height, 0.f)); - camera->raster_to_screen = df_inverse_transform(camera->screen_to_raster); - - /* Perspective projection matrix */ - df_m4 persp = {0.f}; - DF_M4_AT(persp, 0, 0) = 1.f; - DF_M4_AT(persp, 1, 1) = 1.f; - DF_M4_AT(persp, 2, 2) = far_dist / (far_dist - near_dist); - DF_M4_AT(persp, 2, 3) = -1.f * far_dist * near_dist / (far_dist - near_dist); - DF_M4_AT(persp, 3, 2) = 1.f; - - camera->raster_to_camera = df_mul_m4(df_inverse(persp), camera->raster_to_screen); - - camera->focal_dist = focal_dist; - camera->lens_radius = lens_radius; - camera->raster_width = raster_width; - camera->raster_height = raster_height; - camera->image_width = image_width; - camera->image_height = image_height; - - df_camera_i iface = {.release = pc_release, .build_ray_packet = pc_build_ray_packet, .o = camera}; - return iface; -} diff --git a/lib/color.c b/lib/color.c deleted file mode 100644 index 0332bf6..0000000 --- a/lib/color.c +++ /dev/null @@ -1,14 +0,0 @@ -#include - -#include - -df_color df_lerp_color(df_color a, df_color b, double t) -{ - df_color r; - for (int i = 0; i < 4; ++i) { - double af = (double)a.e[i]; - double bf = (double)b.e[i]; - r.e[i] = (uint8_t)floor(af + (bf - af) * t); - } - return r; -} diff --git a/lib/image.c b/lib/image.c deleted file mode 100644 index 84f59b2..0000000 --- a/lib/image.c +++ /dev/null @@ -1,110 +0,0 @@ -#include - -#include -#include -#include - -#define STB_IMAGE_IMPLEMENTATION -#include - -#define STB_IMAGE_WRITE_IMPLEMENTATION -#include - -struct df_image -{ - int width; - int height; - uint8_t *pixels; -}; - -DF_API df_result df_create_image(int w, int h, df_image **out_image) -{ - df_image *img = malloc(sizeof(df_image)); - if (!img) - return df_result_out_of_memory; - DF_ZERO_STRUCT(img); - - img->width = w; - img->height = h; - img->pixels = malloc(4 * sizeof(uint8_t) * w * h); - if (!img->pixels) { - free(img); - return df_result_out_of_memory; - } - memset(img->pixels, 0, 4 * w * h); - - *out_image = img; - return df_result_success; -} - -DF_API df_result df_load_image(const char *path, int *out_w, int *out_h, df_image **out_image) -{ - int w, h, c; - stbi_uc *pixels = stbi_load(path, &w, &h, &c, 4); - if (!pixels) { - fprintf(stderr, "ERROR: failed to load %s: %s\n", path, stbi_failure_reason()); - return df_result_io_error; - } - - df_result res = df_result_success; - df_image *img = malloc(sizeof(df_image)); - if (!img) { - res = df_result_out_of_memory; - goto err; - } - DF_ZERO_STRUCT(img); - - img->width = w; - img->height = h; - img->pixels = (uint8_t *)pixels; - - *out_image = img; - *out_w = w; - *out_h = h; - - goto out; -err: - free(img); -out: - return res; -} - -DF_API void df_release_image(df_image *img) -{ - if (img) { - free(img->pixels); - free(img); - } -} - -DF_API df_result df_write_image(df_image *image, const char *path) -{ - df_result res = stbi_write_png(path, image->width, image->height, 4, image->pixels, image->width * 4) - ? df_result_success - : df_result_io_error; - return res; -} - -DF_API void df_get_image_size(const df_image *image, int *w, int *h) -{ - if (w) - *w = image->width; - if (h) - *h = image->height; -} - -DF_API df_color df_get_image_pixel(const df_image *image, int x, int y) -{ - df_color c = {0, 0, 0, 255}; - if (x >= 0 && x < image->width && y >= 0 && y < image->height) { - memcpy(&c.e[0], &image->pixels[4 * (y * image->width + x)], 4); - } - return c; -} - -DF_API void df_set_image_pixel(df_image *image, int x, int y, df_color c) -{ - if (x >= 0 && x < image->width && y >= 0 && y < image->height) { - memcpy(&image->pixels[4 * (y * image->width + x)], &c.e[0], 4); - } -} \ No newline at end of file diff --git a/lib/log.c b/lib/log.c deleted file mode 100644 index 7d54533..0000000 --- a/lib/log.c +++ /dev/null @@ -1,15 +0,0 @@ -#include -#include -#include - -static const char *log_level_names[] = {"VERBOSE", "INFO", "WARN", "ERROR"}; - -DF_API void df_log_impl(int level, const char *file, int line, const char *fmt, ...) -{ - va_list ap; - va_start(ap, fmt); - fprintf(stderr, "[%s] %s:%d: ", log_level_names[level], file, line); - vfprintf(stderr, fmt, ap); - fprintf(stderr, "\n"); - va_end(ap); -} diff --git a/lib/math.c b/lib/math.c deleted file mode 100644 index 5d1017d..0000000 --- a/lib/math.c +++ /dev/null @@ -1,326 +0,0 @@ -#include -#include - -#include -#include -#include - -DF_API df_v2 df_add_v2(df_v2 a, df_v2 b) -{ - df_v2 v = {a.x + b.x, a.y + b.y}; - return v; -} - -DF_API df_v2 df_sub_v2(df_v2 a, df_v2 b) -{ - df_v2 v = {a.x - b.x, a.y - b.y}; - return v; -} - -DF_API float df_dot_v2(df_v2 a, df_v2 b) { return a.x * b.x + a.y * b.y; } - -DF_API df_v2 df_mul_v2(float t, df_v2 v) -{ - df_v2 r = {t * v.x, t * v.y}; - return r; -} - -DF_API df_v2 df_normalize_v2(df_v2 v) -{ - float len_square = df_dot_v2(v, v); - float len = sqrtf(len_square); - df_v2 n = {v.x / len, v.y / len}; - return n; -} - -DF_API df_v3 df_add_v3(df_v3 a, df_v3 b) -{ - df_v3 v = {a.x + b.x, a.y + b.y, a.z + b.z}; - return v; -} - -DF_API df_v3 df_sub_v3(df_v3 a, df_v3 b) -{ - df_v3 v = {a.x - b.x, a.y - b.y, a.z - b.z}; - return v; -} - -DF_API float df_dot_v3(df_v3 a, df_v3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; } - -DF_API df_v3 df_mul_v3(float t, df_v3 v) -{ - df_v3 r = {t * v.x, t * v.y, t * v.z}; - return r; -} - -DF_API df_v3 df_normalize_v3(df_v3 v) -{ - float len_square = df_dot_v3(v, v); - float len = sqrtf(len_square); - df_v3 n = {v.x / len, v.y / len, v.z / len}; - return n; -} - -DF_API df_line_plane_intersection df_calc_line_plane_intersection(df_line line, df_plane plane) -{ - /* check case */ - float dot = df_dot_v3(line.direction, plane.normal); - df_v3 base_diff = df_sub_v3(plane.base, line.base); - float dot2 = df_dot_v3(base_diff, plane.normal); - if (dot != 0) { - float t = dot2 / dot; - df_v3 point = df_add_v3(line.base, df_mul_v3(t, line.direction)); - df_line_plane_intersection intersection = { - .type = df_line_plane_intersection_vector, - .vec = point, - }; - return intersection; - } else { - if (dot2 == 0) { - df_line_plane_intersection intersection = { - .type = df_line_plane_intersection_contained, - }; - return intersection; - } else { - df_line_plane_intersection intersection = { - .type = df_line_plane_intersection_none, - }; - return intersection; - } - } -} - -DF_API df_m4 df_mul_m4(df_m4 a, df_m4 b) -{ - /* Super simple, we could probably do it a lot better via SIMD. */ - df_m4 p; - for (int row = 0; row < 4; ++row) { - for (int col = 0; col < 4; ++col) { - DF_M4_AT(p, row, col) = 0.f; - for (int i = 0; i < 4; ++i) { - DF_M4_AT(p, row, col) += DF_M4_AT(a, row, i) * DF_M4_AT(b, i, col); - } - } - } - return p; -} - -DF_API df_m4 df_scale(float x, float y, float z) -{ - /* clang-format off */ - df_m4 s = {{ - x, 0.f, 0.f, 0.f, - 0.f, y, 0.f, 0.f, - 0.f, 0.f, z, 0.f, - 0.f, 0.f, 0.f, 1.f, - }}; - /* clang-format on */ - return s; -} - -DF_API df_m4 df_translate(float x, float y, float z) -{ - /* clang-format off */ - df_m4 t = {{ - 1.f, 0.f, 0.f, x, - 0.f, 1.f, 0.f, y, - 0.f, 0.f, 1.f, z, - 0.f, 0.f, 0.f, 1.f, - }}; - /* clang-format on */ - return t; -} - -/** SSE3 version of a horizontal sum: - * computes the sum of 4 32bit floats inside v - */ -static float hsum(__m128 v) -{ - /* v = [v0, v1, v2, v3] */ - - /* shuf = [v1, v1, v3, v3] */ - __m128 shuf = _mm_movehdup_ps(v); - - /* sum = [v0 + v1, 2*v1, v2 + v3, 2*v3] */ - __m128 sum = _mm_add_ps(v, shuf); - - /* shuf = [v2 + v3, 2*v3, v3, v3] */ - shuf = _mm_movehl_ps(shuf, sum); - - /* sum = [v0 + v1 + v2 + v3, ...] */ - sum = _mm_add_ss(sum, shuf); - - return _mm_cvtss_f32(sum); -} - -DF_API df_v3 df_transform_v3(df_m4 T, df_v3 v) -{ - df_v3 transf; - _Alignas(16) float tmp_v[4] = {v.x, v.y, v.z, 1.f}; - __m128 homov = _mm_load_ps(tmp_v); - - __m128 row0 = _mm_load_ps(&T.e[0]); - __m128 row1 = _mm_load_ps(&T.e[4]); - __m128 row2 = _mm_load_ps(&T.e[8]); - - __m128 prod = _mm_mul_ps(row0, homov); - transf.x = hsum(prod); - prod = _mm_mul_ps(row1, homov); - transf.y = hsum(prod); - prod = _mm_mul_ps(row2, homov); - transf.z = hsum(prod); - - return transf; -} - -/* Fast 4x4 matrix inverse via SIMD adapted from: - * https://lxjk.github.io/2017/09/03/Fast-4x4-Matrix-Inverse-with-SSE-SIMD-Explained.html - */ - -/* Only works if M is a transform matrix without scale: - * - * | R T | - * | 0 1 | - */ -DF_API df_m4 df_inverse_transform_no_scale(df_m4 M) -{ - df_m4 I = {0.f}; - - /* transpose 3x3 */ - __m128 t0 = DF_VEC_SHUFFLE_0101(M.vec[0], M.vec[1]); /* 00, 01, 10, 11 */ - __m128 t1 = DF_VEC_SHUFFLE_2323(M.vec[0], M.vec[1]); /* 02, 03, 12, 13 */ - I.vec[0] = DF_VEC_SHUFFLE(t0, M.vec[2], 0, 2, 0, 3); /* 00, 01, 20, 23(=0) */ - I.vec[1] = DF_VEC_SHUFFLE(t0, M.vec[2], 1, 3, 1, 3); /* 01, 11, 21, 23(=0) */ - I.vec[2] = DF_VEC_SHUFFLE(t1, M.vec[2], 0, 2, 2, 3); /* 02, 12, 22, 23(=0) */ - - __m128 t = _mm_setr_ps(DF_M4_AT(M, 0, 3), DF_M4_AT(M, 1, 3), DF_M4_AT(M, 2, 3), 0.f); - DF_M4_AT(I, 0, 3) = -1.f * hsum(_mm_mul_ps(I.vec[0], t)); - - DF_M4_AT(I, 1, 3) = -1.f * hsum(_mm_mul_ps(I.vec[1], t)); - DF_M4_AT(I, 2, 3) = -1.f * hsum(_mm_mul_ps(I.vec[2], t)); - I.vec[3] = _mm_setr_ps(0.f, 0.f, 0.f, 1.f); - - return I; -} - -DF_API df_m4 df_inverse_transform(df_m4 M) -{ - df_m4 I = {0.f}; - - /* transpose 3x3 */ - __m128 t0 = DF_VEC_SHUFFLE_0101(M.vec[0], M.vec[1]); /* 00, 01, 10, 11 */ - __m128 t1 = DF_VEC_SHUFFLE_2323(M.vec[0], M.vec[1]); /* 02, 03, 12, 13 */ - I.vec[0] = DF_VEC_SHUFFLE(t0, M.vec[2], 0, 2, 0, 3); /* 00, 01, 20, 23(=0) */ - I.vec[1] = DF_VEC_SHUFFLE(t0, M.vec[2], 1, 3, 1, 3); /* 01, 11, 21, 23(=0) */ - I.vec[2] = DF_VEC_SHUFFLE(t1, M.vec[2], 0, 2, 2, 3); /* 02, 12, 22, 23(=0) */ - - __m128 t = _mm_setr_ps(DF_M4_AT(M, 0, 3), DF_M4_AT(M, 1, 3), DF_M4_AT(M, 2, 3), 0.f); - DF_M4_AT(I, 0, 3) = -1.f * hsum(_mm_mul_ps(I.vec[0], t)); - - DF_M4_AT(I, 1, 3) = -1.f * hsum(_mm_mul_ps(I.vec[1], t)); - DF_M4_AT(I, 2, 3) = -1.f * hsum(_mm_mul_ps(I.vec[2], t)); - I.vec[3] = _mm_setr_ps(0.f, 0.f, 0.f, 1.f); - - return I; -} - - - -// for row major matrix -// we use __m128 to represent 2x2 matrix as A = | A0 A1 | -// | A2 A3 | -// 2x2 row major Matrix multiply A*B -static __forceinline __m128 mat2_mul(__m128 vec1, __m128 vec2) -{ - return _mm_add_ps(_mm_mul_ps(vec1, DF_VEC_SWIZZLE(vec2, 0, 3, 0, 3)), - _mm_mul_ps(DF_VEC_SWIZZLE(vec1, 1, 0, 3, 2), DF_VEC_SWIZZLE(vec2, 2, 1, 2, 1))); -} -// 2x2 row major Matrix adjugate multiply (A#)*B -static __forceinline __m128 mat2_adj_mul(__m128 vec1, __m128 vec2) -{ - return _mm_sub_ps(_mm_mul_ps(DF_VEC_SWIZZLE(vec1, 3, 3, 0, 0), vec2), - _mm_mul_ps(DF_VEC_SWIZZLE(vec1, 1, 1, 2, 2), DF_VEC_SWIZZLE(vec2, 2, 3, 0, 1))); -} -// 2x2 row major Matrix multiply adjugate A*(B#) -static __forceinline __m128 mat2_mul_adj(__m128 vec1, __m128 vec2) -{ - return _mm_sub_ps(_mm_mul_ps(vec1, DF_VEC_SWIZZLE(vec2, 3, 0, 3, 0)), - _mm_mul_ps(DF_VEC_SWIZZLE(vec1, 1, 0, 3, 2), DF_VEC_SWIZZLE(vec2, 2, 1, 2, 1))); -} - -DF_API df_m4 df_inverse( df_m4 M ) { - // use block matrix method - // A is a matrix, then i(A) or iA means inverse of A, A# (or A_ in code) means adjugate of A, |A| (or detA in code) - // is determinant, tr(A) is trace - - // sub matrices - __m128 A = DF_VEC_SHUFFLE_0101(M.vec[0], M.vec[1]); - __m128 B = DF_VEC_SHUFFLE_2323(M.vec[0], M.vec[1]); - __m128 C = DF_VEC_SHUFFLE_0101(M.vec[2], M.vec[3]); - __m128 D = DF_VEC_SHUFFLE_2323(M.vec[2], M.vec[3]); - -#if 0 - __m128 detA = _mm_set1_ps(M.m[0][0] * M.m[1][1] - M.m[0][1] * M.m[1][0]); - __m128 detB = _mm_set1_ps(M.m[0][2] * M.m[1][3] - M.m[0][3] * M.m[1][2]); - __m128 detC = _mm_set1_ps(M.m[2][0] * M.m[3][1] - M.m[2][1] * M.m[3][0]); - __m128 detD = _mm_set1_ps(M.m[2][2] * M.m[3][3] - M.m[2][3] * M.m[3][2]); -#else - // determinant as (|A| |B| |C| |D|) - __m128 detSub = _mm_sub_ps( - _mm_mul_ps(DF_VEC_SHUFFLE(M.vec[0], M.vec[2], 0, 2, 0, 2), DF_VEC_SHUFFLE(M.vec[1], M.vec[3], 1, 3, 1, 3)), - _mm_mul_ps(DF_VEC_SHUFFLE(M.vec[0], M.vec[2], 1, 3, 1, 3), DF_VEC_SHUFFLE(M.vec[1], M.vec[3], 0, 2, 0, 2))); - __m128 detA = DF_VEC_SWIZZLE1(detSub, 0); - __m128 detB = DF_VEC_SWIZZLE1(detSub, 1); - __m128 detC = DF_VEC_SWIZZLE1(detSub, 2); - __m128 detD = DF_VEC_SWIZZLE1(detSub, 3); -#endif - - // let iM = 1/|M| * | X Y | - // | Z W | - - // D#C - __m128 D_C = mat2_adj_mul(D, C); - // A#B - __m128 A_B = mat2_adj_mul(A, B); - // X# = |D|A - B(D#C) - __m128 X_ = _mm_sub_ps(_mm_mul_ps(detD, A), mat2_mul(B, D_C)); - // W# = |A|D - C(A#B) - __m128 W_ = _mm_sub_ps(_mm_mul_ps(detA, D), mat2_mul(C, A_B)); - - // |M| = |A|*|D| + ... (continue later) - __m128 detM = _mm_mul_ps(detA, detD); - - // Y# = |B|C - D(A#B)# - __m128 Y_ = _mm_sub_ps(_mm_mul_ps(detB, C), mat2_mul_adj(D, A_B)); - // Z# = |C|B - A(D#C)# - __m128 Z_ = _mm_sub_ps(_mm_mul_ps(detC, B), mat2_mul_adj(A, D_C)); - - // |M| = |A|*|D| + |B|*|C| ... (continue later) - detM = _mm_add_ps(detM, _mm_mul_ps(detB, detC)); - - // tr((A#B)(D#C)) - __m128 tr = _mm_mul_ps(A_B, DF_VEC_SWIZZLE(D_C, 0, 2, 1, 3)); - tr = _mm_hadd_ps(tr, tr); - tr = _mm_hadd_ps(tr, tr); - // |M| = |A|*|D| + |B|*|C| - tr((A#B)(D#C) - detM = _mm_sub_ps(detM, tr); - - const __m128 adjSignMask = _mm_setr_ps(1.f, -1.f, -1.f, 1.f); - // (1/|M|, -1/|M|, -1/|M|, 1/|M|) - __m128 rDetM = _mm_div_ps(adjSignMask, detM); - - X_ = _mm_mul_ps(X_, rDetM); - Y_ = _mm_mul_ps(Y_, rDetM); - Z_ = _mm_mul_ps(Z_, rDetM); - W_ = _mm_mul_ps(W_, rDetM); - - df_m4 r; - - // apply adjugate and store, here we combine adjugate shuffle and store shuffle - r.vec[0] = DF_VEC_SHUFFLE(X_, Y_, 3, 1, 3, 1); - r.vec[1] = DF_VEC_SHUFFLE(X_, Y_, 2, 0, 2, 0); - r.vec[2] = DF_VEC_SHUFFLE(Z_, W_, 3, 1, 3, 1); - r.vec[3] = DF_VEC_SHUFFLE(Z_, W_, 2, 0, 2, 0); - - return r; -} diff --git a/lib/pinhole.c b/lib/pinhole.c deleted file mode 100644 index 40e3b49..0000000 --- a/lib/pinhole.c +++ /dev/null @@ -1,76 +0,0 @@ -#include -#include -#include - -#include - -void df_pinhole(df_pinhole_params params, const df_image *in_image, df_image *out_image) -{ - /* orig_z is the z-coordinate of the original image plane (=> in_image is located there) - * new_z is the z-coordinate of the output image plane (=> out_image is located there) - * - * We can map from image coordinates to world coordinates, because we know the z coordinate of - * the input image plane and the focal length of the camera. - * - * Let x,y,z be the world coordinates and u, v the image coordintes. - * The pinhole camera model gives us: - * u = f/z * x; v = f/z * y - * => x = u * z / f - * => y = v * z / f - */ - - int w, h; - df_get_image_size(in_image, &w, &h); - - int out_w, out_h; - df_get_image_size(out_image, &out_w, &out_h); - - double in_z_over_f = params.orig_z / params.focal_length; - double out_f_over_z = params.focal_length / params.new_z; - - int prev_out_iy = -1; - for (int iy = 0; iy < h; ++iy) { - double v = (double)iy / (double)h; - double y = in_z_over_f * v; - double out_v = out_f_over_z * y; - - int out_iy = (int)floor(out_v * (double)out_h); - - int prev_out_ix = -1; - for (int ix = 0; ix < w; ++ix) { - double u = (double)ix / (double)w; - double x = in_z_over_f * u; - double out_u = out_f_over_z * x; - - int out_ix = (int)floor(out_u * (double)out_w); - - df_color px = df_get_image_pixel(in_image, ix, iy); - df_set_image_pixel(out_image, out_ix, out_iy, px); - - /* Go back and interpolate between this pixel and the previous one. */ - if ((out_ix - prev_out_ix) > 1) { - df_color prev_px = df_get_image_pixel(out_image, prev_out_ix, out_iy); - for (int d = prev_out_ix + 1; d < out_ix; ++d) { - double t = (double)(d - prev_out_ix) / (double)(out_ix - prev_out_ix); - df_color color = df_lerp_color(prev_px, px, t); - df_set_image_pixel(out_image, d, out_iy, color); - } - } - prev_out_ix = out_ix; - } - - /* Go back and interpolate between this row and the previous one */ - if ((out_iy - prev_out_iy) > 1) { - for (int y = prev_out_iy + 1; y < out_iy; ++y) { - double t = (double)(y - prev_out_iy) / (double)(out_iy - prev_out_iy); - for (int x = 0; x < out_w; ++x) { - df_color a = df_get_image_pixel(out_image, x, prev_out_iy); - df_color b = df_get_image_pixel(out_image, x, out_iy); - df_color color = df_lerp_color(a, b, t); - df_set_image_pixel(out_image, x, y, color); - } - } - } - prev_out_iy = out_iy; - } -} diff --git a/lib/raytracer.c b/lib/raytracer.c new file mode 100644 index 0000000..d0b9f91 --- /dev/null +++ b/lib/raytracer.c @@ -0,0 +1,347 @@ +#include + +#include +#include +#include + +#define STB_IMAGE_IMPLEMENTATION +#include "stb_image.h" + +/* *********** * + * * + * Vec 3 * + * * + * *********** */ + +/* We can later use this to store 4 vecs for simd */ +typedef struct +{ + float x; + float y; + float z; +} df_v3; + +static df_v3 normalize(df_v3 v) +{ + float l = sqrtf(v.x * v.x + v.y * v.y + v.z * v.z); + return (df_v3){v.x / l, v.y / l, v.z / l}; +} + +static float dot(df_v3 a, df_v3 b) { return a.x * b.x + a.y * b.y + a.z * b.z; } + +/* ****************** * + * * + * List of hits * + * * + * ****************** */ + +typedef struct +{ + df_v3 at; + df_v3 normal; + float ray_t; + int front_face; /* 1 if we hit the front face */ +} df_hit_record; + +static void set_face_normal(df_hit_record *record, df_v3 ray_dir, df_v3 outward_normal) +{ + record->front_face = dot(ray_dir, outward_normal) < 0; + record->normal = + record->front_face ? outward_normal : (df_v3){-outward_normal.x, -outward_normal.y, -outward_normal.z}; +} + +typedef struct +{ + df_hit_record *hits; + size_t count; + size_t capacity; +} df_hit_list; + +static void emit_hit(df_hit_list *list, df_hit_record record) +{ + if (list->count == list->capacity) { + size_t cap2 = (list->capacity > 0) ? 2 * list->capacity : 128; + df_hit_record *tmp = realloc(list->hits, sizeof(df_hit_record) * cap2); + if (!tmp) + return; + list->hits = tmp; + list->capacity = cap2; + } + list->hits[list->count++] = record; +} + +static df_hit_list merge_hit_lists(const df_hit_list **lists, unsigned int count) +{ + size_t total_size = 0; + for (unsigned int i = 0; i < count; ++i) + total_size += lists[i]->count; + + df_hit_list out; + out.count = total_size; + out.hits = malloc(sizeof(df_hit_record) * total_size); + if (!out.hits) + return out; + + df_hit_record *dst = out.hits; + for (unsigned int i = 0; i < count; ++i) { + memcpy(dst, lists[i]->hits, sizeof(df_hit_record) * lists[i]->count); + dst += lists[i]->count; + } + + return out; +} + +/* ********************* * + * * + * Image functions * + * * + * ********************* */ + +#define MAX_IMAGES 1024 + +struct +{ + struct + { + stbi_uc *pixels; + int w; + int h; + } images[MAX_IMAGES]; + unsigned int image_count; +} _image_table; + +DF_API df_image_handle df_load_image(const char *path, int *w, int *h) +{ + if (_image_table.image_count == MAX_IMAGES) + return 0; + + int width, height, c; + stbi_uc *pixels = stbi_load(path, &width, &height, &c, 4); + if (!pixels) { + return 0; + } + df_image_handle handle = (_image_table.image_count + 1); + _image_table.image_count++; + _image_table.images[handle].pixels = pixels; + _image_table.images[handle].w = width; + _image_table.images[handle].h = height; + + if (w) *w = width; + if (h) *h = height; + return handle; +} + +/* ********************************* * + * * + * Intersection test functions * + * + * ********************************* */ + +typedef struct +{ + float ray_t; + float img_u; + float img_v; + df_image_handle image; +} df_hit; + +/* -1 => does not hit, >= 0 => hit (sphere index) */ +static float sphere_test(float ray_origin_x, + float ray_origin_y, + float ray_origin_z, + float ray_dx, + float ray_dy, + float ray_dz, + const df_sphere *spheres, + unsigned int sphere_count) +{ + float result = -1.f; + for (unsigned int i = 0; i < sphere_count; ++i) { + float delta_x = ray_origin_x - spheres[i].center_x; + float delta_y = ray_origin_y - spheres[i].center_y; + float delta_z = ray_origin_z - spheres[i].center_z; + + float a = ray_dx * ray_dx + ray_dy * ray_dy + ray_dz * ray_dz; + float b = (delta_x * ray_dx + delta_y * ray_dy + delta_z * ray_dz); + float c = (delta_x * delta_x + delta_y * delta_y + delta_z * delta_z) - (spheres[i].radius * spheres[i].radius); + float discriminant = b * b - a * c; + + /* we can get the hit location t as: (-b - sqrtf(disciminant)) / (2.f * a) */ + if (discriminant > 0) { + float t = (-b - sqrtf(discriminant)) / (2.f * a); +#if 0 + df_hit_record hit; + hit.at.x = ray_origin_x + t * ray_dx; + hit.at.y = ray_origin_y + t * ray_dy; + hit.at.z = ray_origin_z + t * ray_dz; + hit.ray_t = t; + df_v3 outward_normal = { + hit.at.x - spheres[i].center_x, + hit.at.y - spheres[i].center_y, + hit.at.z - spheres[i].center_z + }; + set_face_normal(&hit, (df_v3){ray_dx, ray_dy, ray_dz}, outward_normal); +#endif + + if (t > result) + result = t; + + } + } + return result; +} + +static df_hit plane_test(float ray_origin_x, + float ray_origin_y, + float ray_origin_z, + float ray_dx, + float ray_dy, + float ray_dz, + const df_plane *planes, + unsigned int plane_count) +{ + df_hit result = { + .ray_t = -1.f, + }; + for (unsigned int i = 0; i < plane_count; ++i) { + float dot = (ray_dx * planes[i].normal_x + ray_dy * planes[i].normal_y + ray_dz * planes[i].normal_z); + if (dot > DF_EPSF32 || dot < -DF_EPSF32) { + float delta_x = planes[i].base_x - ray_origin_x; + float delta_y = planes[i].base_y - ray_origin_y; + float delta_z = planes[i].base_z - ray_origin_z; + + float num = delta_x * planes[i].normal_x + delta_y * planes[i].normal_y + delta_z * planes[i].normal_z; + float t = num / dot; + + if (t > result.ray_t) { + result.ray_t = t; + /* Project point on plane to image coordinate system */ + float px = ray_origin_x + t * ray_dx; + float py = ray_origin_y + t * ray_dy; + float pz = ray_origin_z + t * ray_dz; + + float img_p3_x = px - planes[i].img_p0_x; + float img_p3_y = py - planes[i].img_p0_y; + float img_p3_z = pz - planes[i].img_p0_z; + + /* FIXME(Kevin): We would need to take plane rotation into account. + * Alternatively, just pass w & h into the plane and let the user + * (i.e. higher level code) worry about that */ + float w = planes[i].img_p1_x - planes[i].img_p0_x; + float h = planes[i].img_p1_y - planes[i].img_p0_y; + + result.img_u = + img_p3_x * planes[i].img_ax0_x + img_p3_y * planes[i].img_ax0_y + img_p3_z * planes[i].img_ax0_z; + result.img_v = + img_p3_x * planes[i].img_ax1_x + img_p3_y * planes[i].img_ax1_y + img_p3_z * planes[i].img_ax1_z; + result.img_u /= w; + result.img_v /= h; + result.image = planes[i].image; + } + } + /* Line is parallel to the plane, either contained or not. Do we care? */ + } + return result; +} + +DF_API int df_trace_rays(df_trace_rays_settings settings, + const df_sphere *spheres, + unsigned int sphere_count, + const df_plane *planes, + unsigned int plane_count, + uint8_t **result) +{ + int image_width = settings.image_width; + int image_height = settings.image_height; + + float aspect_ratio = (float)image_width / (float)image_height; + + float viewport_height = 2.f; + float viewport_width = aspect_ratio * viewport_height; + + float focal_length = settings.focal_length; + + /* Simple perspective projection. + * The lens is placed at (0, 0, 0) + */ + + float lower_left_x = -viewport_width / 2.f; + float lower_left_y = -viewport_height / 2.f; + float lower_left_z = -focal_length; + + uint8_t *pixels = malloc(image_width * image_height * 3); + if (!pixels) + return 0; + + float max_img_u = 0.f; + float max_img_v = 0.f; + + for (int y = 0; y < image_height; ++y) { + /* TODO(Kevin): SIMD */ + uint8_t *row = pixels + y * image_width * 3; + for (int x = 0; x < image_width; ++x) { + float u = (float)x / (float)(image_width - 1); + float v = (float)y / (float)(image_height - 1); + + /* Target = Delta because origin is (0, 0, 0) */ + float target_x = lower_left_x + u * viewport_width; + float target_y = lower_left_y + v * viewport_height; + float target_z = lower_left_z; + + /* Raycast against all spheres */ + float sphere_hit_t = sphere_test(0, 0, 0, target_x, target_y, target_z, spheres, sphere_count); + df_hit plane_hit = plane_test(0, 0, 0, target_x, target_y, target_z, planes, plane_count); + + + #if 0 + float hits[] = {sphere_hit_t, plane_hit.ray_t}; + float hit_t = df_max_f32(hits, DF_ARRAY_COUNT(hits)); + + if (hit_t >= 0) { + df_v3 hit_p = {target_x * hit_t, target_y * hit_t, target_z * hit_t}; + hit_p.z -= -1.f; + df_v3 normal = normalize(hit_p); + + row[x * 3 + 0] = (uint8_t)((.5f * (normal.x + 1.f)) * 255); + row[x * 3 + 1] = (uint8_t)((.5f * (normal.y + 1.f)) * 255); + row[x * 3 + 2] = (uint8_t)((.5f * (normal.z + 1.f)) * 255); + } + #endif + if (plane_hit.ray_t >= 0) { + float img_u = plane_hit.img_u; + float img_v = plane_hit.img_v; + + // Temporary, handle arbitrary scale (esp. resulting from different orientation) + float view_aspect = viewport_width / viewport_height; + + if (img_u >= 0 && img_v >= 0 && img_u <= 1.f && img_v <= 1.f) { + + int pixelx = (int)floorf(img_u * (_image_table.images[plane_hit.image].w - 1)); + int pixely = (int)floorf(img_v * (_image_table.images[plane_hit.image].h - 1)); + + stbi_uc *pixel = _image_table.images[plane_hit.image].pixels + + 4 * (pixely * _image_table.images[plane_hit.image].w + pixelx); + + row[x * 3 + 0] = pixel[0]; + row[x * 3 + 1] = pixel[1]; + row[x * 3 + 2] = pixel[2]; + } + } + else { + /* Gradient background color */ + float len = sqrtf(target_x * target_x + target_y * target_y + target_z * target_z); + float t = .5f * (target_y / len + 1.f); + float r = (1.f - t) + t * 0.5f; + float g = (1.f - t) + t * 0.7f; + float b = (1.f - t) + t * 1.0f; + + row[x * 3 + 0] = (uint8_t)(r * 255); + row[x * 3 + 1] = (uint8_t)(g * 255); + row[x * 3 + 2] = (uint8_t)(b * 255); + } + } + } + + *result = pixels; + return 1; +} diff --git a/lib/raytracing.c b/lib/raytracing.c deleted file mode 100644 index 8208ce4..0000000 --- a/lib/raytracing.c +++ /dev/null @@ -1,65 +0,0 @@ -#include -#include - -#include - -#include - -DF_API void df_evaluate_ray_packet(const df_ray_packet *rays, df_image *output) -{ - const __m128 *base_x = (const __m128 *)rays->base_x; - const __m128 *base_y = (const __m128 *)rays->base_y; - const __m128 *base_z = (const __m128 *)rays->base_z; - - const __m128 *dir_x = (const __m128 *)rays->dir_x; - const __m128 *dir_y = (const __m128 *)rays->dir_y; - const __m128 *dir_z = (const __m128 *)rays->dir_z; - - /* Simple test: Let rays intersect with plane at z = 350.0 */ - - float PLANE_Z = 350.0f; - - __m128 plane_z = _mm_set1_ps(PLANE_Z); - - size_t ray_count = rays->ray_count; - - /* TODO(kevin): divide to multiple threads */ - for (size_t i = 0; i < ray_count; i += 4) { - size_t index = i / 4; - - __m128 rays_base_x = base_x[index]; - __m128 rays_base_y = base_y[index]; - __m128 rays_base_z = base_z[index]; - __m128 rays_dir_x = dir_x[index]; - __m128 rays_dir_y = dir_y[index]; - __m128 rays_dir_z = dir_z[index]; - - /* Solve for t: base.z + t * dir.z = plane_z - * t = (plane_z - base.z) / dir.z - */ - __m128 delta = _mm_sub_ps(plane_z, rays_base_z); - __m128 t = _mm_div_ps(delta, rays_dir_z); - - /* Sample p = base.z + t * dir */ - __m128 sample_p_x = _mm_mul_ps(t, rays_dir_x); - __m128 sample_p_y = _mm_mul_ps(t, rays_dir_y); - __m128 sample_p_z = _mm_mul_ps(t, rays_dir_z); - sample_p_x = _mm_add_ps(sample_p_x, rays_base_x); - sample_p_y = _mm_add_ps(sample_p_y, rays_base_y); - sample_p_z = _mm_add_ps(sample_p_z, rays_base_z); - - /* Map the sample position to a 2d (uv) point on the plane, so that we can sample a color */ - df_color color[4]; - for (int j = 0; j < 4; ++j) { - float r = extract_float(sample_p_x, j); - float g = extract_float(sample_p_y, j); - - printf("%f %f\n", r, g); - } - } - - /* Handle remaining (< 4) rays */ - if ((ray_count % 4) != 0) { - - } -} \ No newline at end of file diff --git a/lib/thin_lense.c b/lib/thin_lense.c deleted file mode 100644 index fc74ce3..0000000 --- a/lib/thin_lense.c +++ /dev/null @@ -1,129 +0,0 @@ -#include -#include -#include - -#include -#include - -/* The thin lense model, implemented as a quasi-raytracer. */ - -/** Map uv coordinates in [0, 1) to the unit circle centered at (0, 0). */ -static df_v2 concentric_map_disk(df_v2 in_p) -{ - df_v2 offset = df_sub_v2(df_mul_v2(2.f, in_p), (df_v2){{1.f, 1.f}}); - if (offset.u == 0.f && offset.v == 0.f) - return offset; - - float theta, r; - if (fabsf(offset.u) > fabsf(offset.v)) { - r = offset.u; - theta = DF_PI_OVER_4 * (offset.v / offset.u); - } else { - r = offset.v; - theta = DF_PI_OVER_2 - DF_PI_OVER_4 * (offset.u / offset.v); - } - return df_mul_v2(r, (df_v2){{cosf(theta), sinf(theta)}}); -} - -void df_thin_lense(df_thin_lense_params params, const df_image *in_image, df_image *out_image) -{ - /* We do the following: - * - Orthographic projection (i.e. don't do anything really except flip the y axis) - * - For each pixel: - * - Transform from raster space into camera space - * - Generate ray from pixel to plane of focus through lense center - * - For each sample: - * - Choose random point on lense - * - Trace ray from lense point through the above intersection point to the image. - * - Sum all samples - */ - - int sample_count = (params.sample_count > 0) ? params.sample_count : 64; - - int w, h; - df_get_image_size(out_image, &w, &h); - - /* TODO(Kevin): It would be pretty trivial to - * a) parallelize this, and - * b) use SIMD for multiple rays - */ - for (int y = 0; y < h; ++y) { - float cam_y = (float)(y - h) / 2.f; - for (int x = 0; x < w; ++x) { - float cam_x = (float)(x - w) / 2.f; - - df_v3 ray_dir = {{cam_x, cam_y, 1.f}}; - ray_dir = df_normalize_v3(ray_dir); - ray_dir.z = 1.f; - - /* Calculate the intersection with the plane of focus */ - df_v3 focus_p = df_mul_v3(params.focal_distance, ray_dir); - - uint32_t color[4] = {0, 0, 0, 0}; - for (int sample_idx = 0; sample_idx < sample_count; ++sample_idx) { - df_v2 lens_uv = {(float)(rand() % 1024) / 1024.f, (float)(rand() % 1024) / 1024.f}; - df_v2 lens_p = df_mul_v2(params.aperture, concentric_map_disk(lens_uv)); - - ray_dir.x = focus_p.x - lens_p.x; - ray_dir.y = focus_p.y - lens_p.y; - ray_dir.z = focus_p.z; - ray_dir = df_normalize_v3(ray_dir); - - df_v3 sample_p = df_mul_v3(params.focal_distance, ray_dir); - - int img_x = (int)sample_p.x + w / 2; - int img_y = (int)sample_p.y + h / 2; - - df_color sample_color = df_get_image_pixel(in_image, img_x, img_y); - - for (int i = 0; i < 4; ++i) - color[i] += (uint32_t)sample_color.e[i]; - } - - df_color pixel_color = { - {color[0] / sample_count, color[1] / sample_count, color[2] / sample_count, color[3] / sample_count}}; - - df_set_image_pixel(out_image, x, y, pixel_color); - } - } - -#if 0 - /* The guassian lens equation 1/z' - 1/z = 1/f - * with z' := distance between film and lens, - * z := focal distance, i.e. the distance between lens and the plane of focus - * f := focal length of the lens - * - * gives us: z' = fz / (f + z) */ - float film_dist = (params.focal_length * params.focal_distance) / (params.focal_length + params.focal_distance); - - int sample_count = (params.sample_count > 0) ? params.sample_count : 64; - - int w, h; - df_get_image_size(out_image, &w, &h); - - for (int y = 0; y < h; ++y) { - float v = (float)y / (float)h; - for (int x = 0; x < w; ++x) { - float u = (float)x / (float)w; - df_v2 img_p = {u, v}; - float color[4] = {0, 0, 0, 0}; - - df_v3 focus_ray = {}; - focus_ray.z = 1.f; - focus_ray.x = 0.f; - focus_ray.y = 0.f; - - for (int sample_idx = 0; sample_idx < sample_count; ++sample_idx) { - /* Generate a random sample on the lense */ - df_v2 sample_p = {(float)(rand() % 1024) / 1024.f, (float)(rand() % 1024) / 1024.f}; - df_v2 lens_p = df_mul_v2(params.aperture, concentric_map_disk(sample_p)); - - /* Compute the intersection point on the plane of focus. */ - /* TODO: Use a ray from img_p to the center of the lens, - * trace to focal_distance. - */ - } - } - } -#endif -} \ No newline at end of file diff --git a/lib/utils.c b/lib/utils.c new file mode 100644 index 0000000..20eafbd --- /dev/null +++ b/lib/utils.c @@ -0,0 +1,10 @@ +#include +#include + +DF_API float df_max_f32(const float *list, unsigned int count) +{ + float max = -INFINITY; + for (unsigned int i = 0; i < count; ++i) + max = (max < list[i]) ? list[i] : max; + return max; +} diff --git a/meson.build b/meson.build index c5dd12a..4502ba9 100644 --- a/meson.build +++ b/meson.build @@ -9,28 +9,20 @@ if cc.get_id() == 'gcc' or cc.get_id() == 'clang' add_project_arguments([ '-msse3', '-msse4.1', '-Wno-missing-braces' ], language: 'c') endif -lib = library('df', - 'lib/log.c', - 'lib/math.c', - 'lib/camera.c', - 'lib/pinhole.c', - 'lib/image.c', - 'lib/color.c', - 'lib/thin_lense.c', - 'lib/raytracing.c', - 'include/defocus/base.h', - 'include/defocus/camera.h', - 'include/defocus/defocus.h', - 'include/defocus/image.h', - 'include/defocus/intrinsic_helper.h', - 'include/defocus/models.h', - 'include/defocus/scene.h', - 'include/defocus/raytracing.h', +# Main library +lib = library('df', + 'lib/raytracer.c', + 'lib/utils.c', + 'include/defocus/defocus.h', + '3p/stb_image.h', include_directories: incdir, dependencies: m_dep) # Command Line Executable -executable('defocus', 'bin/defocus.c', include_directories: incdir, link_with: lib) +executable('defocus', + 'bin/defocus.c', + '3p/stb_image_write.h', + include_directories: incdir, link_with: lib) # Test driver #munit_dep = dependency('munit', fallback: ['munit', 'munit_dep']) diff --git a/test_image.png b/test_image.png new file mode 100644 index 0000000..80a41d1 --- /dev/null +++ b/test_image.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4ae7ea9a98c3ced9719a2597d02149ea3ec22abacd367a1cd8e0762828708b67 +size 359433