dump state

does not take plane rotation into account!

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.

3p/pcg/pcg_basic.c

@@ -0,0 +1,116 @@
+/*
+ * PCG Random Number Generation for C.
+ *
+ * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
+ *
+ * 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);
+}
+

3p/pcg/pcg_basic.h

@@ -0,0 +1,78 @@
+/*
+ * PCG Random Number Generation for C.
+ *
+ * Copyright 2014 Melissa O'Neill <oneill@pcg-random.org>
+ *
+ * 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 <inttypes.h>
+
+#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

bin/defocus.c

@@ -1,242 +1,100 @@
-#include <stdio.h>
-#include <string.h>
-#include <stdlib.h>
+/* Main command line executable. */
+
 #include <defocus/defocus.h>
+#include <defocus/df_math.h>
 
-#include <defocus/camera.h>
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+#include <stb_image_write.h>
 
-int pinhole_fn(int argc, char **argv);
-int thin_lense_fn(int argc, char **argv);
-int test_fn(int argc, char **argv);
-
-static const char *_model_names[] = {"pinhole", "thin_lense", "test"};
-typedef int (*model_fn)(int argc, char **argv);
-
-static model_fn _model_fns[] = {
-    pinhole_fn,
-    thin_lense_fn,
-    test_fn,
-};
-
-void usage(const char *pname)
+df_plane get_image_filling_plane(int width, int height, df_image_handle image, float focal_length)
 {
-    printf("Usage: %s <model-name> [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]);
-    }
+    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)
-{
-    if (argc < 2) {
-        fprintf(stderr, "Missing model name!\n");
-        usage(argv[0]);
-        return 1;
-    }
+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;
 
-    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;
-    }
+    int image_width = 0;
+    int image_height = 0;
+    float focal_length = 1.f;
 
-    return fn(argc - 2, &argv[2]);
-}
+    df_image_handle input_image = df_load_image("../test_image.png", &image_width, &image_height);
+    float aspect = (float)image_width / (float)image_height;
     
-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];
-        }
-    }
+    df_plane plane;
+    plane.base_x = 0.f;
+    plane.base_y = 0.f;
+    plane.base_z = -2;
+    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;
 
-    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 <path>\t\tPath to the input image.\n");
-        printf("  -o  | --output <path>\t\tPath to the output image (Default: out.png).\n");
-        printf("  -iz | --input-z <number>\t\tDistance to the input image plane.\n");
-        printf("  -oz | --output-z <number>\t\tDistance to the output image plane.\n");
-        printf("  -f  | --focal-length <number>\t\tThe focal length of the camera lens.\n");
-        return 1;
-    }
+    uint8_t *image = NULL;
     
-    float focal, in_z, out_z;
-    focal = atof(focal_str);
-    in_z = atof(in_z_str);
-    out_z = atof(out_z_str);
+    df_v3 look_from = {1, 1.5, 5};
+    df_v3 look_at = {0, 0, -2};
+    df_v3 up = {0, 1, 0};
 
-    df_image *in_image, *out_image;
+    float dist_to_focus = df_len_v3(df_sub_v3(look_at, look_from));
 
-    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;
-    }
+    df_thin_lense_camera_data camera_data = df_create_thin_lense_camera_data(look_from, look_at, up, 45.f, aspect, 2.0, dist_to_focus);
 
-    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_trace_rays((df_trace_rays_settings){
+            .image_width = image_width,
+            .image_height = image_height,
+            .samples_per_pixel = 64,
+            .camera = df_pinhole_camera,
+            .camera_data = &camera_data,
+        },
+        spheres, 0,
+        &plane, 1,
+        &image);
 
-    df_pinhole((df_pinhole_params){.orig_z = in_z, .new_z = out_z, .focal_length = focal}, in_image, out_image);
+    stbi_write_png("raytracer.png", image_width, image_height, 3, image, image_width * 3);
 
-    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 <path>\t\tPath to the input image.\n");
-        printf("  -o  | --output <path>\t\tPath to the output image (Default: out.png).\n");
-        printf("  -iz | --input-z <number>\t\tDistance to the input image plane.\n");
-        printf("  -oz | --output-z <number>\t\tDistance to the output image plane.\n");
-        printf("  -f  | --focal-length <number>\t\tThe focal length of the camera lens.\n");
-        printf("  -a  | --aperture <number>\t\tThe size of the camera lens.\n");
-        printf("  -s  | --samples <number>\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);
-    df_evaluate_ray_packet(&packet);
-
-    df_release_ray_packet(&packet);
-    
-    cam.release(cam.o);
     return 0;
 }

include/defocus/base.h

@@ -1,278 +0,0 @@
-#ifndef DEFOCUS_BASE_H
-#define DEFOCUS_BASE_H
-
-#include <inttypes.h>
-#include <string.h>
-
-#include <immintrin.h>
-
-/** @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

include/defocus/camera.h

@@ -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

include/defocus/defocus.h

@@ -1,11 +1,172 @@
-#ifndef DEFOCUS_DEFOCUS_H
-#define DEFOCUS_DEFOCUS_H
+#ifndef DEFOCUS_H
+#define DEFOCUS_H
 
-/** @file defocus.h
- * @brief meta header
+#include <stdint.h>
+
+#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 */
+
+/* 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_PI_OVER_180F32 (DF_PIF32 / 180.f)
+
+#define DF_DEGREES_TO_RADIANS(deg) ((deg) * DF_PI_OVER_180F32)
+
+#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 */
+
+/* The camera data structs (df_pinhole_camera_data etc.) are defined here
+ * to enable you to store these on the stack.
+ * However - unless you know what you are doing (-: - you should treat these
+ * structs as opaque.
  */
-#include "base.h"
-#include "image.h"
-#include "models.h"
+
+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_v3 horizontal;
+    df_v3 vertical;
+    df_v3 lower_left;
+    df_v3 origin;
+} df_pinhole_camera_data;
+
+/* Creates data for a pinhole camera.
+ * The vertical fov (vfov) is given in degrees.
+ */
+DF_API df_pinhole_camera_data df_create_pinhole_camera_data(df_v3 look_from, df_v3 look_at, df_v3 up,
+    float vfov, float aspect);
+
+DF_API df_ray df_pinhole_camera(float u, float v, unsigned int sample_idx, void *camera_data);
+
+typedef struct
+{
+    df_v3 horizontal;
+    df_v3 vertical;
+    df_v3 lower_left;
+    df_v3 origin;
+    float lens_radius;
+} df_thin_lense_camera_data;
+
+DF_API df_thin_lense_camera_data df_create_thin_lense_camera_data(df_v3 look_from, df_v3 look_at, df_v3 up,
+    float vfov, float aspect, 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;
+
+    df_v3 look_from;
+    df_v3 look_at;
+
+    df_camera_fn camera;
+    void *camera_data;
+} 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_w;
+    float img_h;
+
+    /* TODO(Kevin): These could be calculated from p0 and p1 (p0+wh)... */
+    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.
+ */
+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

include/defocus/image.h

@@ -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

include/defocus/intrinsic_helper.h

@@ -1,38 +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 <immintrin.h>
-#include <pmmintrin.h>
-
-/* 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)
-
-#endif

include/defocus/models.h

@@ -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

include/defocus/raytracing.h

@@ -1,29 +0,0 @@
-#ifndef DF_RAYTRACING_H
-#define DF_RAYTRACING_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);
-
-#endif

include/defocus/scene.h

@@ -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

lib/camera.c

@@ -1,152 +0,0 @@
-#include <defocus/camera.h>
-
-#include <string.h>
-#include <stdlib.h>
-#include <assert.h>
-
-#ifdef _MSC_VER
-#include <malloc.h>
-#endif
-
-/* ********************************************************
- *
- * Ray Packet Functions
- *
- * ********************************************************/
-
-DF_API void df_release_ray_packet(df_ray_packet *rays)
-{
-    free(rays->simd_mem);
-    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;
-}

lib/color.c

@@ -1,14 +0,0 @@
-#include <defocus/base.h>
-
-#include <math.h>
-
-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;
-}

lib/image.c

@@ -1,110 +0,0 @@
-#include <defocus/image.h>
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-
-#define STB_IMAGE_IMPLEMENTATION
-#include <stb_image.h>
-
-#define STB_IMAGE_WRITE_IMPLEMENTATION
-#include <stb_image_write.h>
-
-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);
-    }
-}

lib/log.c

@@ -1,15 +0,0 @@
-#include <defocus/base.h>
-#include <stdio.h>
-#include <stdarg.h>
-
-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);
-}

lib/math.c

@@ -1,326 +0,0 @@
-#include <defocus/base.h>
-#include <defocus/intrinsic_helper.h>
-
-#include <math.h>
-#include <immintrin.h>
-#include <pmmintrin.h>
-
-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;
-}

lib/pinhole.c

@@ -1,76 +0,0 @@
-#include <defocus/models.h>
-#include <defocus/base.h>
-#include <defocus/image.h>
-
-#include <math.h>
-
-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;
-    }
-}

lib/raytracer.c

@@ -0,0 +1,471 @@
+#include <defocus/defocus.h>
+#include <defocus/df_math.h>
+
+#include <math.h>
+#include <stdlib.h>
+#include <string.h>
+
+#define STB_IMAGE_IMPLEMENTATION
+#include "stb_image.h"
+
+#include "pcg/pcg_basic.h"
+
+/* ********************* *
+ *                       *
+ *    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;
+}
+
+/* *************************** *
+ *                             *
+ *    Builtin camera models    *
+ *                             *
+ * *************************** */
+
+DF_API df_pinhole_camera_data
+df_create_pinhole_camera_data(df_v3 look_from, df_v3 look_at, df_v3 up, float vfov, float aspect)
+{
+    float theta = DF_DEGREES_TO_RADIANS(vfov);
+    float h = tanf(theta * .5f);
+    float viewport_height = 2.f * h;
+    float viewport_width = aspect * viewport_height;
+
+    df_v3 view_vec = df_sub_v3(look_from, look_at);
+    df_v3 w = df_normalize_v3(view_vec);
+    df_v3 u = df_normalize_v3(df_cross(up, w));
+    df_v3 v = df_cross(w, u);
+
+    df_v3 horizontal = df_mul_v3(viewport_width, u);
+    df_v3 vertical = df_mul_v3(viewport_height, v);
+    df_v3 lower_left = {
+        .x = look_from.x - horizontal.x / 2.f - vertical.x / 2.f - w.x,
+        .y = look_from.y - horizontal.y / 2.f - vertical.y / 2.f - w.y,
+        .z = look_from.z - horizontal.z / 2.f - vertical.z / 2.f - w.z,
+    };
+
+    return (df_pinhole_camera_data){
+        .origin = look_from,
+        .horizontal = horizontal,
+        .vertical = vertical,
+        .lower_left = lower_left,
+    };
+}
+
+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 = data->origin;
+    df_v3 horizontal = data->horizontal;
+    df_v3 vertical = data->vertical;
+
+    df_v3 target = {
+        data->lower_left.x + u * horizontal.x + v * vertical.x - origin.x, 
+        data->lower_left.y + u * horizontal.y + v * vertical.y - origin.y,
+        data->lower_left.z + u * horizontal.z + v * vertical.z - origin.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 (df_len_v2(p) <= 1.f)
+            return p;
+    }
+}
+
+DF_API df_thin_lense_camera_data df_create_thin_lense_camera_data(
+    df_v3 look_from, df_v3 look_at, df_v3 up, float vfov, float aspect, float aperture, float focal_distance)
+{
+    float theta = DF_DEGREES_TO_RADIANS(vfov);
+    float h = tanf(theta * .5f);
+    float viewport_height = 2.f * h;
+    float viewport_width = aspect * viewport_height;
+
+    df_v3 view_vec = df_sub_v3(look_from, look_at);
+    df_v3 w = df_normalize_v3(view_vec);
+    df_v3 u = df_normalize_v3(df_cross(up, w));
+    df_v3 v = df_cross(w, u);
+
+    df_v3 horizontal = df_mul_v3(viewport_width, u);
+    df_v3 vertical = df_mul_v3(viewport_height, v);
+    df_v3 lower_left = {
+        .x = look_from.x - horizontal.x / 2.f - vertical.x / 2.f - focal_distance * w.x,
+        .y = look_from.y - horizontal.y / 2.f - vertical.y / 2.f - focal_distance * w.y,
+        .z = look_from.z - horizontal.z / 2.f - vertical.z / 2.f - focal_distance * w.z,
+    };
+
+    return (df_thin_lense_camera_data){
+        .origin = look_from,
+        .horizontal = horizontal,
+        .vertical = vertical,
+        .lower_left = lower_left,
+        .lens_radius = aperture / 2.f,
+    };
+}
+
+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 = df_add_v3(data->origin, offset);
+    df_v3 horizontal = data->horizontal;
+    df_v3 vertical = data->vertical;
+
+    df_v3 target = {
+        data->lower_left.x + u * horizontal.x + v * vertical.x - origin.x,
+        data->lower_left.y + u * horizontal.y + v * vertical.y - origin.y,
+        data->lower_left.z + u * horizontal.z + v * vertical.z - origin.z,
+    };
+
+    return (df_ray){
+        .origin = origin,
+        .dir = target,
+    };
+}
+
+
+/* ********************************* *
+ *                                   *
+ *    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(df_v3 ray_origin,
+                        df_v3 ray_d,
+                        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_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 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_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;
+                float img_p3_z = pz - planes[i].img_p0_z;
+
+                float w = planes[i].img_w;
+                float h = planes[i].img_h;
+
+                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;
+}
+
+/* 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,
+                         const df_plane *planes,
+                         unsigned int plane_count,
+                         uint8_t **result)
+{
+    int image_width = settings.image_width;
+    int image_height = settings.image_height;
+    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;
+
+    float *fpixels = malloc(sizeof(float) * image_width * image_height * 3);
+    if (!fpixels)
+        return 0;
+    memset(fpixels, 0, sizeof(float) * image_width * image_height * 3);
+
+    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(fpixels);
+        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 */
+        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) */
+            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};
+
+            /* 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};
+
+                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));
+
+            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);
+            }
+
+#endif
+        }
+    }
+    
+    float fsamples = (float)samples_per_pixel;
+    /* Trace the rays */
+    for (unsigned int i = 0; i < packet.ray_count; ++i) {
+        df_v2i pixelp = packet.ray_pixels[i];
+
+        float *row = fpixels + 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] += (float)pixel[0];
+                row[pixelp.x * 3 + 1] += (float)pixel[1];
+                row[pixelp.x * 3 + 2] += (float)pixel[2];
+            }
+        }
+    }
+
+    /* Map down to 8bit rgb */
+    uint8_t *pixels = malloc(image_width * image_height * 3);
+    if (!pixels) {
+        free(packet_mem);
+        free(fpixels);
+        return 0;
+    }
+    for (unsigned int i = 0; i < image_width * image_height; ++i) {
+        float *fp = &fpixels[i * 3];
+        uint8_t *p = &pixels[i * 3];
+
+        p[0] = (uint8_t)(fp[0] / fsamples);
+        p[1] = (uint8_t)(fp[1] / fsamples);
+        p[2] = (uint8_t)(fp[2] / fsamples);
+    }
+
+    free(fpixels);
+    free(packet_mem);
+    *result = pixels;
+    return 1;
+}

lib/raytracing.c

@@ -1,51 +0,0 @@
-#include <defocus/camera.h>
-
-#include <immintrin.h>
-
-DF_API void df_evaluate_ray_packet(const df_ray_packet *rays) {
-    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) {
-        __m128 rays_base_x = base_x[i];
-        __m128 rays_base_y = base_y[i];
-        __m128 rays_base_z = base_z[i];
-        __m128 rays_dir_x = dir_x[i];
-        __m128 rays_dir_y = dir_y[i];
-        __m128 rays_dir_z = dir_z[i];
-
-        /* 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);
-
-    }
-
-    /* Handle remaining (< 4) rays */
-    if ((ray_count % 4) != 0) {
-    
-    }
-}

lib/thin_lense.c

@@ -1,129 +0,0 @@
-#include <defocus/base.h>
-#include <defocus/models.h>
-#include <defocus/image.h>
-
-#include <math.h>
-#include <stdlib.h>
-
-/* 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
-}

lib/utils.c

@@ -0,0 +1,10 @@
+#include <defocus/defocus.h>
+#include <math.h>
+
+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;
+}

meson.build

@@ -9,28 +9,24 @@ if cc.get_id() == 'gcc' or cc.get_id() == 'clang'
   add_project_arguments([ '-msse3', '-msse4.1', '-Wno-missing-braces' ], language: 'c')
 endif
 
+# Main library
 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',
+              'lib/raytracer.c',
+              'lib/utils.c',
+              'lib/df_math.c',
 	          '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',
+              'include/defocus/df_math.h',
+              '3p/stb_image.h',
+              '3p/pcg/pcg_basic.c',
+              '3p/pcg/pcg_basic.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'])

subprojects/.gitignore

@@ -0,0 +1 @@
+munit