diff --git a/3p/pcg/LICENSE.txt b/3p/pcg/LICENSE.txt new file mode 100644 index 0000000..8dada3e --- /dev/null +++ b/3p/pcg/LICENSE.txt @@ -0,0 +1,201 @@ + Apache License + Version 2.0, January 2004 + http://www.apache.org/licenses/ + + TERMS AND CONDITIONS FOR USE, REPRODUCTION, AND DISTRIBUTION + + 1. Definitions. + + "License" shall mean the terms and conditions for use, reproduction, + and distribution as defined by Sections 1 through 9 of this document. + + "Licensor" shall mean the copyright owner or entity authorized by + the copyright owner that is granting the License. + + "Legal Entity" shall mean the union of the acting entity and all + other entities that control, are controlled by, or are under common + control with that entity. For the purposes of this definition, + "control" means (i) the power, direct or indirect, to cause the + direction or management of such entity, whether by contract or + otherwise, or (ii) ownership of fifty percent (50%) or more of the + outstanding shares, or (iii) beneficial ownership of such entity. + + "You" (or "Your") shall mean an individual or Legal Entity + exercising permissions granted by this License. + + "Source" form shall mean the preferred form for making modifications, + including but not limited to software source code, documentation + source, and configuration files. + + "Object" form shall mean any form resulting from mechanical + transformation or translation of a Source form, including but + not limited to compiled object code, generated documentation, + and conversions to other media types. + + "Work" shall mean the work of authorship, whether in Source or + Object form, made available under the License, as indicated by a + copyright notice that is included in or attached to the work + (an example is provided in the Appendix below). + + "Derivative Works" shall mean any work, whether in Source or Object + form, that is based on (or derived from) the Work and for which the + editorial revisions, annotations, elaborations, or other modifications + represent, as a whole, an original work of authorship. For the purposes + of this License, Derivative Works shall not include works that remain + separable from, or merely link (or bind by name) to the interfaces of, + the Work and Derivative Works thereof. + + "Contribution" shall mean any work of authorship, including + the original version of the Work and any modifications or additions + to that Work or Derivative Works thereof, that is intentionally + submitted to Licensor for inclusion in the Work by the copyright owner + or by an individual or Legal Entity authorized to submit on behalf of + the copyright owner. For the purposes of this definition, "submitted" + means any form of electronic, verbal, or written communication sent + to the Licensor or its representatives, including but not limited to + communication on electronic mailing lists, source code control systems, + and issue tracking systems that are managed by, or on behalf of, the + Licensor for the purpose of discussing and improving the Work, but + excluding communication that is conspicuously marked or otherwise + designated in writing by the copyright owner as "Not a Contribution." + + "Contributor" shall mean Licensor and any individual or Legal Entity + on behalf of whom a Contribution has been received by Licensor and + subsequently incorporated within the Work. + + 2. Grant of Copyright License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + copyright license to reproduce, prepare Derivative Works of, + publicly display, publicly perform, sublicense, and distribute the + Work and such Derivative Works in Source or Object form. + + 3. Grant of Patent License. Subject to the terms and conditions of + this License, each Contributor hereby grants to You a perpetual, + worldwide, non-exclusive, no-charge, royalty-free, irrevocable + (except as stated in this section) patent license to make, have made, + use, offer to sell, sell, import, and otherwise transfer the Work, + where such license applies only to those patent claims licensable + by such Contributor that are necessarily infringed by their + Contribution(s) alone or by combination of their Contribution(s) + with the Work to which such Contribution(s) was submitted. If You + institute patent litigation against any entity (including a + cross-claim or counterclaim in a lawsuit) alleging that the Work + or a Contribution incorporated within the Work constitutes direct + or contributory patent infringement, then any patent licenses + granted to You under this License for that Work shall terminate + as of the date such litigation is filed. + + 4. Redistribution. You may reproduce and distribute copies of the + Work or Derivative Works thereof in any medium, with or without + modifications, and in Source or Object form, provided that You + meet the following conditions: + + (a) You must give any other recipients of the Work or + Derivative Works a copy of this License; and + + (b) You must cause any modified files to carry prominent notices + stating that You changed the files; and + + (c) You must retain, in the Source form of any Derivative Works + that You distribute, all copyright, patent, trademark, and + attribution notices from the Source form of the Work, + excluding those notices that do not pertain to any part of + the Derivative Works; and + + (d) If the Work includes a "NOTICE" text file as part of its + distribution, then any Derivative Works that You distribute must + include a readable copy of the attribution notices contained + within such NOTICE file, excluding those notices that do not + pertain to any part of the Derivative Works, in at least one + of the following places: within a NOTICE text file distributed + as part of the Derivative Works; within the Source form or + documentation, if provided along with the Derivative Works; or, + within a display generated by the Derivative Works, if and + wherever such third-party notices normally appear. The contents + of the NOTICE file are for informational purposes only and + do not modify the License. You may add Your own attribution + notices within Derivative Works that You distribute, alongside + or as an addendum to the NOTICE text from the Work, provided + that such additional attribution notices cannot be construed + as modifying the License. + + You may add Your own copyright statement to Your modifications and + may provide additional or different license terms and conditions + for use, reproduction, or distribution of Your modifications, or + for any such Derivative Works as a whole, provided Your use, + reproduction, and distribution of the Work otherwise complies with + the conditions stated in this License. + + 5. Submission of Contributions. Unless You explicitly state otherwise, + any Contribution intentionally submitted for inclusion in the Work + by You to the Licensor shall be under the terms and conditions of + this License, without any additional terms or conditions. + Notwithstanding the above, nothing herein shall supersede or modify + the terms of any separate license agreement you may have executed + with Licensor regarding such Contributions. + + 6. Trademarks. This License does not grant permission to use the trade + names, trademarks, service marks, or product names of the Licensor, + except as required for reasonable and customary use in describing the + origin of the Work and reproducing the content of the NOTICE file. + + 7. Disclaimer of Warranty. Unless required by applicable law or + agreed to in writing, Licensor provides the Work (and each + Contributor provides its Contributions) on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or + implied, including, without limitation, any warranties or conditions + of TITLE, NON-INFRINGEMENT, MERCHANTABILITY, or FITNESS FOR A + PARTICULAR PURPOSE. You are solely responsible for determining the + appropriateness of using or redistributing the Work and assume any + risks associated with Your exercise of permissions under this License. + + 8. Limitation of Liability. In no event and under no legal theory, + whether in tort (including negligence), contract, or otherwise, + unless required by applicable law (such as deliberate and grossly + negligent acts) or agreed to in writing, shall any Contributor be + liable to You for damages, including any direct, indirect, special, + incidental, or consequential damages of any character arising as a + result of this License or out of the use or inability to use the + Work (including but not limited to damages for loss of goodwill, + work stoppage, computer failure or malfunction, or any and all + other commercial damages or losses), even if such Contributor + has been advised of the possibility of such damages. + + 9. Accepting Warranty or Additional Liability. While redistributing + the Work or Derivative Works thereof, You may choose to offer, + and charge a fee for, acceptance of support, warranty, indemnity, + or other liability obligations and/or rights consistent with this + License. However, in accepting such obligations, You may act only + on Your own behalf and on Your sole responsibility, not on behalf + of any other Contributor, and only if You agree to indemnify, + defend, and hold each Contributor harmless for any liability + incurred by, or claims asserted against, such Contributor by reason + of your accepting any such warranty or additional liability. + + END OF TERMS AND CONDITIONS + + APPENDIX: How to apply the Apache License to your work. + + To apply the Apache License to your work, attach the following + boilerplate notice, with the fields enclosed by brackets "{}" + replaced with your own identifying information. (Don't include + the brackets!) The text should be enclosed in the appropriate + comment syntax for the file format. We also recommend that a + file or class name and description of purpose be included on the + same "printed page" as the copyright notice for easier + identification within third-party archives. + + Copyright {yyyy} {name of copyright owner} + + Licensed under the Apache License, Version 2.0 (the "License"); + you may not use this file except in compliance with the License. + You may obtain a copy of the License at + + http://www.apache.org/licenses/LICENSE-2.0 + + Unless required by applicable law or agreed to in writing, software + distributed under the License is distributed on an "AS IS" BASIS, + WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + See the License for the specific language governing permissions and + limitations under the License. diff --git a/3p/pcg/pcg_basic.c b/3p/pcg/pcg_basic.c new file mode 100644 index 0000000..8c2fd0d --- /dev/null +++ b/3p/pcg/pcg_basic.c @@ -0,0 +1,116 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the full C implementation, which is in turn + * derived from the canonical C++ PCG implementation. The C++ version + * has many additional features and is preferable if you can use C++ in + * your project. + */ + +#include "pcg_basic.h" + +// state for global RNGs + +static pcg32_random_t pcg32_global = PCG32_INITIALIZER; + +// pcg32_srandom(initstate, initseq) +// pcg32_srandom_r(rng, initstate, initseq): +// Seed the rng. Specified in two parts, state initializer and a +// sequence selection constant (a.k.a. stream id) + +void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, uint64_t initseq) +{ + rng->state = 0U; + rng->inc = (initseq << 1u) | 1u; + pcg32_random_r(rng); + rng->state += initstate; + pcg32_random_r(rng); +} + +void pcg32_srandom(uint64_t seed, uint64_t seq) +{ + pcg32_srandom_r(&pcg32_global, seed, seq); +} + +// pcg32_random() +// pcg32_random_r(rng) +// Generate a uniformly distributed 32-bit random number + +uint32_t pcg32_random_r(pcg32_random_t* rng) +{ + uint64_t oldstate = rng->state; + rng->state = oldstate * 6364136223846793005ULL + rng->inc; + uint32_t xorshifted = ((oldstate >> 18u) ^ oldstate) >> 27u; + uint32_t rot = oldstate >> 59u; + return (xorshifted >> rot) | (xorshifted << ((-rot) & 31)); +} + +uint32_t pcg32_random() +{ + return pcg32_random_r(&pcg32_global); +} + + +// pcg32_boundedrand(bound): +// pcg32_boundedrand_r(rng, bound): +// Generate a uniformly distributed number, r, where 0 <= r < bound + +uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, uint32_t bound) +{ + // To avoid bias, we need to make the range of the RNG a multiple of + // bound, which we do by dropping output less than a threshold. + // A naive scheme to calculate the threshold would be to do + // + // uint32_t threshold = 0x100000000ull % bound; + // + // but 64-bit div/mod is slower than 32-bit div/mod (especially on + // 32-bit platforms). In essence, we do + // + // uint32_t threshold = (0x100000000ull-bound) % bound; + // + // because this version will calculate the same modulus, but the LHS + // value is less than 2^32. + + uint32_t threshold = -bound % bound; + + // Uniformity guarantees that this loop will terminate. In practice, it + // should usually terminate quickly; on average (assuming all bounds are + // equally likely), 82.25% of the time, we can expect it to require just + // one iteration. In the worst case, someone passes a bound of 2^31 + 1 + // (i.e., 2147483649), which invalidates almost 50% of the range. In + // practice, bounds are typically small and only a tiny amount of the range + // is eliminated. + for (;;) { + uint32_t r = pcg32_random_r(rng); + if (r >= threshold) + return r % bound; + } +} + + +uint32_t pcg32_boundedrand(uint32_t bound) +{ + return pcg32_boundedrand_r(&pcg32_global, bound); +} + diff --git a/3p/pcg/pcg_basic.h b/3p/pcg/pcg_basic.h new file mode 100644 index 0000000..e2b526a --- /dev/null +++ b/3p/pcg/pcg_basic.h @@ -0,0 +1,78 @@ +/* + * PCG Random Number Generation for C. + * + * Copyright 2014 Melissa O'Neill + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + * For additional information about the PCG random number generation scheme, + * including its license and other licensing options, visit + * + * http://www.pcg-random.org + */ + +/* + * This code is derived from the full C implementation, which is in turn + * derived from the canonical C++ PCG implementation. The C++ version + * has many additional features and is preferable if you can use C++ in + * your project. + */ + +#ifndef PCG_BASIC_H_INCLUDED +#define PCG_BASIC_H_INCLUDED 1 + +#include + +#if __cplusplus +extern "C" { +#endif + +struct pcg_state_setseq_64 { // Internals are *Private*. + uint64_t state; // RNG state. All values are possible. + uint64_t inc; // Controls which RNG sequence (stream) is + // selected. Must *always* be odd. +}; +typedef struct pcg_state_setseq_64 pcg32_random_t; + +// If you *must* statically initialize it, here's one. + +#define PCG32_INITIALIZER { 0x853c49e6748fea9bULL, 0xda3e39cb94b95bdbULL } + +// pcg32_srandom(initstate, initseq) +// pcg32_srandom_r(rng, initstate, initseq): +// Seed the rng. Specified in two parts, state initializer and a +// sequence selection constant (a.k.a. stream id) + +void pcg32_srandom(uint64_t initstate, uint64_t initseq); +void pcg32_srandom_r(pcg32_random_t* rng, uint64_t initstate, + uint64_t initseq); + +// pcg32_random() +// pcg32_random_r(rng) +// Generate a uniformly distributed 32-bit random number + +uint32_t pcg32_random(void); +uint32_t pcg32_random_r(pcg32_random_t* rng); + +// pcg32_boundedrand(bound): +// pcg32_boundedrand_r(rng, bound): +// Generate a uniformly distributed number, r, where 0 <= r < bound + +uint32_t pcg32_boundedrand(uint32_t bound); +uint32_t pcg32_boundedrand_r(pcg32_random_t* rng, uint32_t bound); + +#if __cplusplus +} +#endif + +#endif // PCG_BASIC_H_INCLUDED diff --git a/bin/defocus.c b/bin/defocus.c index a7ada1c..e96f4b4 100644 --- a/bin/defocus.c +++ b/bin/defocus.c @@ -5,6 +5,33 @@ #define STB_IMAGE_WRITE_IMPLEMENTATION #include +df_plane get_image_filling_plane(int width, int height, df_image_handle image, float focal_length) +{ + float aspect = (float)width / (float)height; + + df_plane plane; + plane.base_x = 0.f; + plane.base_y = 0.f; + plane.base_z = -2 * focal_length; + plane.normal_x = 0.f; + plane.normal_y = 0.f; + plane.normal_z = 1.f; + plane.img_p0_x = -aspect / focal_length; + plane.img_p0_y = -1.f / focal_length; + plane.img_p0_z = plane.base_z; + plane.img_w = 2.f * aspect / focal_length; + plane.img_h = 2.f / focal_length; + 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 = image; + + return plane; +} + int main(int argc, char **argv) { df_sphere spheres[2]; spheres[0].center_x = 0.f; @@ -16,39 +43,24 @@ int main(int argc, char **argv) { spheres[1].center_z = -1.5f; spheres[1].radius = .75f; - int image_width = 1024; - int image_height = 512; + int image_width = 0; + int image_height = 0; + float focal_length = 1.f; df_image_handle input_image = df_load_image("../test_image.png", &image_width, &image_height); float aspect = (float)image_width / (float)image_height; - - 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_w = 2.f * aspect; - plane.img_h = 2.f; - 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; + df_plane plane = get_image_filling_plane(image_width, image_height, input_image, focal_length); uint8_t *image = NULL; + df_thin_lense_camera_data camera_data = df_create_thin_lense_camera_data(image_width, image_height, 42.f, focal_length); df_trace_rays((df_trace_rays_settings){ - .focal_length = 1.f, .image_width = image_width, .image_height = image_height, + .samples_per_pixel = 64, + .camera = df_thin_lense_camera, + .camera_data = &camera_data, }, spheres, 0, &plane, 1, diff --git a/include/defocus/defocus.h b/include/defocus/defocus.h index 397935c..5c1cb47 100644 --- a/include/defocus/defocus.h +++ b/include/defocus/defocus.h @@ -21,25 +21,87 @@ #endif /* __cplusplus */ +/* Useful constants */ #define DF_EPSF32 0.00001f +#define DF_PIF32 3.1415926f +#define DF_PI_OVER_4F32 (DF_PIF32 / 4.f) +#define DF_PI_OVER_2F32 (DF_PIF32 / 2.f) #define DF_ARRAY_COUNT(A) (sizeof((A)) / sizeof((A)[0])) +#define DF_UNUSED(x) ((void)sizeof((x))) + +/* Math stuff */ + +typedef struct +{ + float x; + float y; + float z; +} df_v3; + +typedef struct +{ + float x; + float y; +} df_v2; + +typedef struct +{ + int x; + int y; +} df_v2i; /* Images */ typedef unsigned int df_image_handle; DF_API df_image_handle df_load_image(const char *path, int *w, int *h); +/* cameras */ + +typedef struct +{ + df_v3 origin; + df_v3 dir; +} df_ray; + +/* A camera function takes the uv coordinates of the output pixel and generates a ray. + * The sample index can be used to introduce some randomness for multi-sampling. */ +typedef df_ray (*df_camera_fn)(float u, float v, unsigned int sample_idx, void *camera_data); + +typedef struct +{ + df_v2 viewport_size; + df_v3 lower_left; +} df_pinhole_camera_data; + +DF_API df_pinhole_camera_data df_create_pinhole_camera_data(int image_width, int image_height, float focal_length); + +DF_API df_ray df_pinhole_camera(float u, float v, unsigned int sample_idx, void *camera_data); + +typedef struct +{ + df_v2 viewport_size; + df_v3 lower_left; + float lens_radius; +} df_thin_lense_camera_data; + +DF_API df_thin_lense_camera_data df_create_thin_lense_camera_data(int image_width, + int image_height, + float aperture, + float focal_distance); + +DF_API df_ray df_thin_lense_camera(float u, float v, unsigned int sample_idx, void *camera_data); + /* Settings for the basic ray tracing function. */ typedef struct { /* Width needs to be divisible by 4! */ int image_width; int image_height; + int samples_per_pixel; - /* Distance between lense and image plane */ - float focal_length; - float lens_radius; + df_camera_fn camera; + void *camera_data; } df_trace_rays_settings; /* Sphere shape */ @@ -67,7 +129,7 @@ typedef struct float img_w; float img_h; - /* TODO(Kevin): These could be calculated from p0 and p1... */ + /* TODO(Kevin): These could be calculated from p0 and p1 (p0+wh)... */ float img_ax0_x; float img_ax0_y; float img_ax0_z; @@ -77,7 +139,6 @@ typedef struct df_image_handle image; } df_plane; - DF_API float df_max_f32(const float *list, unsigned int count); /* Core function, implements raytracing. diff --git a/lib/raytracer.c b/lib/raytracer.c index 9795d94..3e67ff3 100644 --- a/lib/raytracer.c +++ b/lib/raytracer.c @@ -7,19 +7,13 @@ #define STB_IMAGE_IMPLEMENTATION #include "stb_image.h" -/* *********** * - * * - * Vec 3 * - * * - * *********** */ +#include "pcg/pcg_basic.h" -/* We can later use this to store 4 vecs for simd */ -typedef struct -{ - float x; - float y; - float z; -} df_v3; +/* ************* * + * * + * Vectors * + * * + * ************* */ static df_v3 normalize(df_v3 v) { @@ -29,67 +23,8 @@ static df_v3 normalize(df_v3 v) 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 * - * * - * ****************** */ +static float v2_len(df_v2 v) { return sqrtf(v.x * v.x + v.y * v.y); } -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; -} /* ********************* * * * @@ -131,10 +66,115 @@ DF_API df_image_handle df_load_image(const char *path, int *w, int *h) return handle; } +/* *************************** * + * * + * Builtin camera models * + * * + * *************************** */ + +DF_API df_pinhole_camera_data df_create_pinhole_camera_data(int image_width, int image_height, float focal_length) +{ + float aspect_ratio = (float)image_width / (float)image_height; + float viewport_height = 2.f; + float viewport_width = aspect_ratio * viewport_height; + + /* 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; + + return (df_pinhole_camera_data){ + .viewport_size = {.x = viewport_width, .y = viewport_height}, + .lower_left = { .x = lower_left_x, .y = lower_left_y, .z = lower_left_z }, + }; +} + +DF_API df_ray df_pinhole_camera(float u, float v, unsigned int sample_idx, void *camera_data) +{ + DF_UNUSED(sample_idx); + df_pinhole_camera_data *data = camera_data; + + df_v3 origin = {0.f, 0.f, 0.f}; + df_v3 target = { + data->lower_left.x + u * data->viewport_size.x, + data->lower_left.y + v * data->viewport_size.y, + data->lower_left.z + }; + + return (df_ray){ + .origin = origin, + .dir = target, + }; +} + +static df_v2 random_point_on_disk(unsigned int i) +{ + for (unsigned int j = 0; j < i; ++j) + pcg32_random(); + + while (1) { + df_v2 p = {.x = (float)pcg32_boundedrand(1024) / 1023.f, .y = (float)pcg32_boundedrand(1024) / 1023.f}; + if (v2_len(p) <= 1.f) + return p; + } +} + +DF_API df_thin_lense_camera_data df_create_thin_lense_camera_data(int image_width, + int image_height, + float aperture, + float focal_distance) +{ + float aspect_ratio = (float)image_width / (float)image_height; + float viewport_height = 2.f; + float viewport_width = aspect_ratio * viewport_height; + + /* 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_distance; + + float lens_radius = aperture / 2.f; + + return (df_thin_lense_camera_data){ + .lens_radius = lens_radius, + .lower_left = { + .x = lower_left_x, + .y = lower_left_y, + .z = lower_left_z, + }, + .viewport_size = { + .x = viewport_width, + .y = viewport_height + } + }; +} + +DF_API df_ray df_thin_lense_camera(float u, float v, unsigned int sample_idx, void *camera_data) +{ + df_thin_lense_camera_data *data = camera_data; + + df_v2 lensp = random_point_on_disk(sample_idx); + lensp.x *= data->lens_radius; + lensp.y *= data->lens_radius; + + df_v3 offset = {u * lensp.x, v * lensp.y, .z = 0.f}; + + df_v3 origin = offset; + df_v3 target = {data->lower_left.x + u * data->viewport_size.x + offset.x, + data->lower_left.y + v * data->viewport_size.y + offset.y, + data->lower_left.z + offset.z}; + return (df_ray){.origin = origin, .dir = target}; +} + + /* ********************************* * * * * Intersection test functions * - * + * * * ********************************* */ typedef struct @@ -191,12 +231,8 @@ static float sphere_test(float ray_origin_x, 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, +static df_hit plane_test(df_v3 ray_origin, + df_v3 ray_d, const df_plane *planes, unsigned int plane_count) { @@ -204,11 +240,11 @@ static df_hit plane_test(float ray_origin_x, .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); + float dot = (ray_d.x * planes[i].normal_x + ray_d.y * planes[i].normal_y + ray_d.z * 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 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; @@ -216,9 +252,9 @@ static df_hit plane_test(float ray_origin_x, 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 px = ray_origin.x + t * ray_d.x; + float py = ray_origin.y + t * ray_d.y; + float pz = ray_origin.z + t * ray_d.z; float img_p3_x = px - planes[i].img_p0_x; float img_p3_y = py - planes[i].img_p0_y; @@ -241,6 +277,17 @@ static df_hit plane_test(float ray_origin_x, return result; } +/* Ray packets to support more than one ray (=sample) per pixel */ +typedef struct +{ + unsigned int samples_per_pixel; + unsigned int ray_count; + + df_v3 *ray_origins; /* (0, 0, 0) is the center of the lens */ + df_v3 *ray_deltas; /* not necessarily normalized! */ + df_v2i *ray_pixels; /* pixel coordinates in the output image */ +} df_ray_packet; + DF_API int df_trace_rays(df_trace_rays_settings settings, const df_sphere *spheres, unsigned int sphere_count, @@ -250,47 +297,77 @@ DF_API int df_trace_rays(df_trace_rays_settings settings, { 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; + unsigned int samples_per_pixel = (unsigned int)settings.samples_per_pixel; + if (samples_per_pixel == 0) + samples_per_pixel = 1; + df_camera_fn camera = settings.camera; + if (!camera) + return 0; + void *camera_data = settings.camera_data; uint8_t *pixels = malloc(image_width * image_height * 3); if (!pixels) return 0; + memset(pixels, 0, image_width * image_height * 3); - float max_img_u = 0.f; - float max_img_v = 0.f; + df_ray_packet packet; + packet.samples_per_pixel = samples_per_pixel; + packet.ray_count = image_width * image_height * samples_per_pixel; + void *packet_mem = malloc(packet.ray_count * (2 * sizeof(df_v3) + sizeof(df_v2i))); + if (!packet_mem) { + free(pixels); + return 0; + } + packet.ray_origins = (df_v3 *)packet_mem; + packet.ray_deltas = packet.ray_origins + packet.ray_count; + packet.ray_pixels = (df_v2i *)(packet.ray_deltas + packet.ray_count); + /* FIXME(Kevin): Replace with dynamic (time() ?) initializer */ + pcg32_srandom(1523789, 901842398023); + + /* Generate the rays */ + unsigned int ray_idx = 0; 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); + for (unsigned int samplei = 0; samplei < samples_per_pixel; ++samplei) { + packet.ray_pixels[ray_idx] = (df_v2i){x, y}; + + df_ray ray = camera(u, v, samplei, camera_data); + + packet.ray_origins[ray_idx] = ray.origin; + packet.ray_deltas[ray_idx] = ray.dir; + ++ray_idx; + } +#if 0 /* 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; + df_v3 origin = {0.f, 0.f, 0.f}; + df_v3 target = {lower_left_x + u * viewport_width, lower_left_y + v * viewport_height, 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); + /* Adjust for lens effect */ + if (settings.lens_radius > 0.f) { + df_v2 lensp = concentric_sample_disk((df_v2){u, v}); + lensp.x *= settings.lens_radius; + lensp.y *= settings.lens_radius; + float ft = -focal_length / target.z; + df_v3 focusp = {target.x * ft, target.y * ft, target.z * ft}; - #if 0 + origin.x = lensp.x; + origin.y = lensp.y; + origin.z = 0.f; + + target = (df_v3){focusp.x - lensp.x, focusp.y - lensp.y, focusp.z}; + } + + /* Raycast against all planes */ + df_hit plane_hit = plane_test(origin, target, 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)); @@ -303,7 +380,7 @@ DF_API int df_trace_rays(df_trace_rays_settings settings, 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 +#endif if (plane_hit.ray_t >= 0) { float img_u = plane_hit.img_u; float img_v = plane_hit.img_v; @@ -326,8 +403,8 @@ DF_API int df_trace_rays(df_trace_rays_settings settings, } 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 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; @@ -336,9 +413,39 @@ DF_API int df_trace_rays(df_trace_rays_settings settings, row[x * 3 + 1] = (uint8_t)(g * 255); row[x * 3 + 2] = (uint8_t)(b * 255); } + +#endif + } + } + + /* Trace the rays */ + for (unsigned int i = 0; i < packet.ray_count; ++i) { + df_v2i pixelp = packet.ray_pixels[i]; + + uint8_t *row = pixels + pixelp.y * image_width * 3; + + df_hit plane_hit = plane_test(packet.ray_origins[i], packet.ray_deltas[i], planes, plane_count); + + if (plane_hit.ray_t >= 0) { + float img_u = plane_hit.img_u; + float img_v = plane_hit.img_v; + + 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[pixelp.x * 3 + 0] += pixel[0] / samples_per_pixel; + row[pixelp.x * 3 + 1] += pixel[1] / samples_per_pixel; + row[pixelp.x * 3 + 2] += pixel[2] / samples_per_pixel; + } } } + free(packet_mem); *result = pixels; return 1; } diff --git a/meson.build b/meson.build index 4502ba9..3e909bf 100644 --- a/meson.build +++ b/meson.build @@ -15,6 +15,8 @@ lib = library('df', 'lib/utils.c', 'include/defocus/defocus.h', '3p/stb_image.h', + '3p/pcg/pcg_basic.c', + '3p/pcg/pcg_basic.h', include_directories: incdir, dependencies: m_dep)