Finish 3D level set

include/taichi/math/levelset_2d.cpp

@@ -56,7 +56,7 @@ Vector2 LevelSet2D::get_normalized_gradient(const Vector2 &pos) const
 
 real LevelSet2D::get(const Vector2 &pos) const
 {
-    assert_info(inside(pos), "LevelSet Gradient Query out of Bound! (" + std::to_string(pos.x) + ", " + std::to_string(pos.y) + ")");
+    assert_info(inside(pos), "LevelSet Query out of Bound! (" + std::to_string(pos.x) + ", " + std::to_string(pos.y) + ")");
     real x = pos.x, y = pos.y;
     x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps);
     y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps);

include/taichi/math/levelset_3d.cpp

@@ -0,0 +1,113 @@
+/*******************************************************************************
+    Taichi - Physically based Computer Graphics Library
+
+    Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
+
+    All rights reserved. Use of this source code is governed by
+    the MIT license as written in the LICENSE file.
+*******************************************************************************/
+
+#include "levelset_3d.h"
+
+TC_NAMESPACE_BEGIN
+
+
+void LevelSet3D::add_sphere(Vector3 center, real radius, bool inside_out) {
+    for (auto &ind : get_region()) {
+        Vector3 sample = ind.get_pos();
+        real dist = (inside_out ? -1 : 1) * (length(center - sample) - radius);
+        set(ind, std::min(Array3D::get(ind), dist));
+    }
+}
+
+Vector3 LevelSet3D::get_gradient(const Vector3 &pos) const
+{
+    assert_info(inside(pos), "LevelSet Gradient Query out of Bound! ("
+                             + std::to_string(pos.x) + ", "
+                             + std::to_string(pos.y) + ", "
+                             + std::to_string(pos.z) + ")");
+    real x = pos.x, y = pos.y, z = pos.z;
+    x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps);
+    y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps);
+    z = clamp(z - storage_offset.z, 0.f, depth - 1.f - eps);
+    const int x_i = clamp(int(x), 0, width - 2);
+    const int y_i = clamp(int(y), 0, height - 2);
+    const int z_i = clamp(int(z), 0, depth - 2);
+    const real x_r = x - x_i;
+    const real y_r = y - y_i;
+    const real z_r = z - z_i;
+    const real gx = lerp(y_r,
+                         lerp(z_r, Array3D::get(x_i + 1, y_i, z_i) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i + 1, y_i, z_i + 1) - Array3D::get(x_i, y_i, z_i + 1)),
+                         lerp(z_r, Array3D::get(x_i + 1, y_i + 1, z_i) - Array3D::get(x_i, y_i + 1, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i + 1, z_i + 1)));
+    const real gy = lerp(z_r,
+                         lerp(x_r, Array3D::get(x_i, y_i + 1, z_i) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i + 1, y_i + 1, z_i) - Array3D::get(x_i + 1, y_i, z_i)),
+                         lerp(x_r, Array3D::get(x_i, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i, z_i + 1), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i + 1, y_i, z_i + 1)));
+    const real gz = lerp(x_r,
+                         lerp(y_r, Array3D::get(x_i, y_i, z_i + 1) - Array3D::get(x_i, y_i, z_i), Array3D::get(x_i, y_i + 1, z_i + 1) - Array3D::get(x_i, y_i + 1, z_i)),
+                         lerp(y_r, Array3D::get(x_i + 1, y_i, z_i + 1) - Array3D::get(x_i + 1, y_i, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1) - Array3D::get(x_i + 1, y_i + 1, z_i)));
+    return Vector3(gx, gy, gz);
+}
+
+Vector3 LevelSet3D::get_normalized_gradient(const Vector3 &pos) const
+{
+    Vector3 gradient = get_gradient(pos);
+    if (length(gradient) < 1e-10f)
+        return Vector3(1, 0, 0);
+    else
+        return normalize(gradient);
+}
+
+real LevelSet3D::get(const Vector3 &pos) const
+{
+    assert_info(inside(pos), "LevelSet Query out of Bound! ("
+                             + std::to_string(pos.x) + ", "
+                             + std::to_string(pos.y) + ", "
+                             + std::to_string(pos.z) + ")");
+    real x = pos.x, y = pos.y, z = pos.z;
+    x = clamp(x - storage_offset.x, 0.f, width - 1.f - eps);
+    y = clamp(y - storage_offset.y, 0.f, height - 1.f - eps);
+    z = clamp(z - storage_offset.z, 0.f, depth - 1.f - eps);
+    const int x_i = clamp(int(x), 0, width - 2);
+    const int y_i = clamp(int(y), 0, height - 2);
+    const int z_i = clamp(int(z), 0, depth - 2);
+    const real x_r = x - x_i;
+    const real y_r = y - y_i;
+    const real z_r = z - z_i;
+    return lerp(x_r,
+                lerp(y_r,
+                     lerp(z_r, Array3D::get(x_i, y_i, z_i), Array3D::get(x_i, y_i, z_i + 1)),
+                     lerp(z_r, Array3D::get(x_i, y_i + 1, z_i), Array3D::get(x_i, y_i + 1, z_i + 1))),
+                lerp(y_r,
+                     lerp(z_r, Array3D::get(x_i + 1, y_i, z_i), Array3D::get(x_i + 1, y_i, z_i + 1)),
+                     lerp(z_r, Array3D::get(x_i + 1, y_i + 1, z_i), Array3D::get(x_i + 1, y_i + 1, z_i + 1))));
+}
+
+
+Array3D<real> LevelSet3D::rasterize(int width, int height, int depth) {
+    for (auto &p : (*this)) {
+        if (std::isnan(p)) {
+            printf("Warning: nan in levelset.");
+        }
+    }
+    Array3D<real> out(width, height, depth);
+    Vector3 actual_size;
+    if (storage_offset == Vector3(0.0f, 0.0f, 0.0f)) {
+        actual_size = Vector3(this->width - 1, this->height - 1, this->depth - 1);
+    }
+    else {
+        actual_size = Vector3(this->width, this->height, this->depth);
+    }
+
+    Vector3 scale_factor = actual_size / Vector3(width, height, depth);
+
+    for (auto &ind : Region3D(0, width, 0, height, 0, depth, Vector3(0.5f, 0.5f, 0.5f))) {
+        Vector3 p = scale_factor * ind.get_pos();
+        out[ind] = sample(p);
+        if (std::isnan(out[ind])) {
+            out[ind] = std::numeric_limits<real>::infinity();
+        }
+    }
+    return out;
+}
+
+TC_NAMESPACE_END

include/taichi/math/levelset_3d.h

@@ -0,0 +1,71 @@
+/*******************************************************************************
+    Taichi - Physically based Computer Graphics Library
+
+    Copyright (c) 2017 Yu Fang <squarefk@gmail.com>
+
+    All rights reserved. Use of this source code is governed by
+    the MIT license as written in the LICENSE file.
+*******************************************************************************/
+
+#pragma once
+
+#include <limits>
+#include <taichi/math/array_3d.h>
+
+
+TC_NAMESPACE_BEGIN
+
+class LevelSet3D : public Array3D<real> {
+public:
+    real friction = 1.0f;
+
+    LevelSet3D() : LevelSet3D(0, 0, 0) {}
+
+    LevelSet3D(int width, int height, int depth, Vector3 offset = Vector3(0.5f, 0.5f, 0.5f)) {
+        initialize(width, height, depth, offset);
+    }
+
+    void initialize(int width, int height, int depth, Vector3 offset) {
+        Array3D<real>::initialize(width, height, depth, std::numeric_limits<real>::infinity(), offset);
+    }
+
+    void initialize(const Vector3i &res, Vector3 offset) {
+        Array3D<real>::initialize(res[0], res[1], res[2], std::numeric_limits<real>::infinity(), offset);
+    }
+
+    void initialize(int width, int height, int depth, Vector3 offset, real value) {
+        Array3D<real>::initialize(width, height, depth, value, offset);
+    }
+
+    void initialize(const Vector3i &res, Vector3 offset, real value) {
+        Array3D<real>::initialize(res[0], res[1], res[2], value, offset);
+    }
+
+    void add_sphere(Vector3 center, real radius, bool inside_out = false);
+
+    Vector3 get_gradient(const Vector3 &pos) const; // Note this is not normalized!
+
+    Vector3 get_normalized_gradient(const Vector3 &pos) const;
+
+    real get(const Vector3 &pos) const;
+
+    static real fraction_outside(real phi_a, real phi_b) {
+        return 1.0f - fraction_inside(phi_a, phi_b);
+    }
+
+    static real fraction_inside(real phi_a, real phi_b) {
+        if (phi_a < 0 && phi_b < 0)
+            return 1;
+        if (phi_a < 0 && phi_b >= 0)
+            return phi_a / (phi_a - phi_b);
+        if (phi_a >= 0 && phi_b < 0)
+            return phi_b / (phi_b - phi_a);
+        else
+            return 0;
+    }
+
+    Array3D<real> rasterize(int width, int height, int depth);
+};
+
+TC_NAMESPACE_END
+