Merge pull request #39 from squarefk/ampm

MeshTexture

include/taichi/geometry/primitives.h

@@ -169,6 +169,12 @@ struct Triangle {
     }
 };
 
+class BoundingBox {
+public:
+    Vector3 lower_boundary;
+    Vector3 upper_boundary;
+};
+
 class Instance {
     Matrix4 transform;
     int id;

include/taichi/visual/scene.cpp

@@ -123,7 +123,7 @@ void Scene::add_mesh(std::shared_ptr<Mesh> mesh) {
     meshes.push_back(*mesh);
 }
 
-void Scene::finalize() {
+void Scene::finalize_geometry() {
     int triangle_count = 0;
     for (auto &mesh : meshes) {
         triangle_id_start[&mesh] = triangle_count;
@@ -142,6 +142,9 @@ void Scene::finalize() {
     }
     num_triangles = triangle_count;
     printf("Scene loaded. Triangle count: %d\n", triangle_count);
+};
+
+void Scene::finalize_lighting() {
     if (!emissive_triangles.empty()) {
         update_emission_cdf();
         update_light_emission_cdf();
@@ -152,4 +155,9 @@ void Scene::finalize() {
     }
 }
 
+void Scene::finalize() {
+    finalize_geometry();
+    finalize_lighting();
+}
+
 TC_NAMESPACE_END

include/taichi/visual/scene.h

@@ -41,6 +41,26 @@ class Mesh {
     void set_untransformed_triangles(const std::vector<Triangle> &triangles) {
         untransformed_triangles = triangles;
     }
+    // bounding box
+    BoundingBox get_bounding_box() {
+        BoundingBox bb;
+        bool is_first_triangle = true;
+        for (auto t : untransformed_triangles) {
+            if (is_first_triangle) {
+                bb.lower_boundary = t.v[0];
+                bb.upper_boundary = t.v[0];
+                is_first_triangle = false;
+            }
+            // loop for xyz
+            for (int i = 0; i < 3; ++i)
+                // loop for each vertex
+                for (int j = 0; j < 3; ++j){
+                    bb.lower_boundary[i] = std::min(bb.lower_boundary[i], t.v[j][i]);
+                    bb.upper_boundary[i] = std::max(bb.upper_boundary[i], t.v[j][i]);
+                }
+        }
+        return bb;
+    }
     std::vector<Triangle> get_triangles() {
         std::vector<Triangle> triangles;
         for (auto t : untransformed_triangles) {
@@ -141,6 +161,10 @@ class Scene {
 
     void add_mesh(std::shared_ptr<Mesh> mesh);
 
+    void finalize_geometry();
+
+    void finalize_lighting();
+
     void finalize();
 
     void update_light_emission_cdf() {

python/examples/simulation/3d/mpm_snow_ball_hit_3d.py

@@ -0,0 +1,106 @@
+import math
+
+from taichi.dynamics.mpm import MPM3
+from taichi.core import tc_core
+from taichi.misc.util import Vector
+from taichi.visual import *
+from taichi.visual.post_process import *
+from taichi.visual.texture import Texture
+import taichi as tc
+
+gi_render = True
+step_number = 400
+# step_number = 1
+# total_frames = 1
+grid_downsample = 2
+output_downsample = 1
+render_epoch = 20
+
+
+def create_mpm_sand_block(fn):
+    particles = tc_core.RenderParticles()
+    assert particles.read(fn)
+    downsample = grid_downsample
+    tex = Texture.from_render_particles((255 / downsample, 255 / downsample, 255 / downsample), particles) * 5
+    # mesh_transform = tc_core.Matrix4(1.0).scale_s(0.5).translate(Vector(0.5, 0.5, 0.5))
+    # transform = tc_core.Matrix4(1.0).scale_s(2).scale(Vector(2.0, 0.5, 1.0)).translate(Vector(-2, -0.99, -1))
+    mesh_transform = tc_core.Matrix4(1.0).translate(Vector(0, 0.01, 0))
+    transform = tc_core.Matrix4(1.0).scale_s(2).translate(Vector(-1, -1, -1))
+    vol = VolumeMaterial('sdf_voxel', scattering=5, absorption=0, tex=tex,
+                         resolution=(255 / downsample, 255 / downsample, 255 / downsample),
+                         transform_ptr=transform.get_ptr_string())
+    material = SurfaceMaterial('plain_interface')
+    material.set_internal_material(vol)
+    return Mesh('cube', material=material, transform=transform * mesh_transform)
+
+
+def create_sand_scene(frame, d, t):
+    downsample = output_downsample
+    width, height = 540 / downsample, 540 / downsample
+    camera = Camera('thinlens', width=width, height=height, fov=75,
+                    origin=(0, 1, 4), look_at=(0.0, -0.9, 0.0), up=(0, 1, 0), aperture=0.01)
+
+    scene = Scene()
+    with scene:
+        scene.set_camera(camera)
+        rep = Texture.create_taichi_wallpaper(10, rotation=0, scale=0.95) * Texture('const', value=(0.7, 0.5, 0.5))
+        material = SurfaceMaterial('pbr', diffuse_map=rep)
+        scene.add_mesh(Mesh('holder', material=material, translate=(0, -1, -6), scale=2))
+
+        mesh = Mesh('plane', SurfaceMaterial('emissive', color=(1, 1, 1)),
+                    translate=(1.0, 1.0, -1), scale=(0.1, 0.1, 0.1), rotation=(180, 0, 0))
+        scene.add_mesh(mesh)
+
+        material = tc.SurfaceMaterial('microfacet', color=(1, 1, 0.5), roughness=(0.1, 0, 0, 0), f0=1)
+        sphere = tc.Mesh('sphere', material,
+                         translate=((t+0.05) * 0.5 - 0.35, -0.61, 0), scale=0.1, rotation=(0, 0, 0))
+        scene.add_mesh(sphere)
+
+        # Change this line to your particle output path pls.
+        # fn = r'../sand-sim/particles%05d.bin' % frame
+        fn = d + r'/particles%05d.bin' % frame
+        mesh = create_mpm_sand_block(fn)
+        scene.add_mesh(mesh)
+
+    return scene
+
+
+def render_sand_frame(frame, d, t):
+    renderer = Renderer(output_dir='volumetric', overwrite=True, frame=frame)
+    renderer.initialize(preset='pt', scene=create_sand_scene(frame, d, t), sampler='prand')
+    renderer.set_post_processor(LDRDisplay(exposure=0.6, bloom_radius=0.0, bloom_threshold=1.0))
+    renderer.render(render_epoch)
+
+
+if __name__ == '__main__':
+    downsample = grid_downsample
+    resolution = (255 / downsample, 255 / downsample, 255 / downsample)
+    # tex = Texture('ring', outer=0.15) * 4
+    # tex = Texture('bound', tex=tex, axis=2, bounds=(0.0, 0.4), outside_val=(0, 0, 0))
+    # tex = Texture('rotate', tex=tex, rotate_axis=0, rotate_times=1)
+    tex = Texture('mesh', resolution=resolution, filename=tc.get_asset_path('meshes/suzanne.obj')) * 8
+    tex = tex.zoom((0.5, 0.5, 0.5), (0.5, 0.5, 0.5), False)
+    # tex = Texture('rotate', tex=tex, rotate_axis=1, rotate_times=1)
+    mpm = MPM3(resolution=resolution, gravity=(0, -10, 0), initial_velocity=(0, 0, 0), delta_t=0.0005, num_threads=8,
+               density_tex=tex.id)
+
+
+    # Dynamic Levelset
+    def levelset_generator(t):
+        levelset = mpm.create_levelset()
+        levelset.add_sphere(Vector(0.325 + 0.25 * (t+0.05), 0.2, 0.5), 0, False)
+        # levelset.add_sphere(Vector(0.5, 0.2, 0.5), t, False)
+        return levelset
+    mpm.set_levelset(levelset_generator, True)
+
+    t = 0
+    for i in range(step_number):
+        print 'process(%d/%d)' % (i, step_number)
+        mpm.step(0.01)
+        t += 0.01
+        if gi_render:
+            d = mpm.get_directory()
+            if i % 10 == 0:
+                render_sand_frame(i, d, t)
+                pass
+    mpm.make_video()

python/taichi/dynamics/mpm.py

@@ -52,14 +52,16 @@ def step(self, step_t):
         print 'Simulation time:', t
         T = time.time()
         self.update_levelset(t, t + step_t)
+        print 'Update Leveset Time:', time.time() - T
+        T = time.time()
         self.c.step(step_t)
-        print 'Time:', time.time() - T
+        print 'Step Time:', time.time() - T
         image_buffer = tc_core.Array2DVector3(self.video_manager.width, self.video_manager.height, Vector(0, 0, 0.0))
         particles = self.c.get_render_particles()
         particles.write(self.directory + '/particles%05d.bin' % self.frame)
         res = map(float, self.resolution)
         camera = Camera('pinhole', origin=(0, res[1] * 0.4, res[2] * 1.4),
-                        look_at=(0, -res[1] * 0.5, 0), up=(0, 1, 0), fov=70,
+                        look_at=(0, -res[1] * 0.5, 0), up=(0, 1, 0), fov=90,
                         width=10, height=10)
         self.particle_renderer.set_camera(camera)
         self.particle_renderer.render(image_buffer, particles)

src/texture/sources.cpp

@@ -8,6 +8,7 @@
 *******************************************************************************/
 
 #include <taichi/visual/texture.h>
+#include <taichi/visual/scene_geometry.h>
 #include <taichi/visualization/image_buffer.h>
 #include <taichi/math/array_3d.h>
 #include <taichi/common/asset_manager.h>
@@ -284,5 +285,77 @@ class SlicedTexture : public Texture {
 
 TC_IMPLEMENTATION(Texture, SlicedTexture, "sliced");
 
+class MeshTexture : public Texture {
+protected:
+    Array3D<char> arr;
+    Vector3i resolution;
+public:
+    // parameter name for mesh path is 'filename'
+    void initialize(const Config &config) override {
+        Texture::initialize(config);
+        Mesh mesh;
+        mesh.initialize(config);
+        Scene scene;
+        scene.add_mesh(std::make_shared<Mesh>(mesh));
+        scene.finalize_geometry();
+        auto ray_intersection = create_instance<RayIntersection>("embree");
+        std::shared_ptr<SceneGeometry> scene_geometry(new SceneGeometry(std::make_shared<Scene>(scene), ray_intersection));
+        BoundingBox bb = mesh.get_bounding_box();
+        // calc offset and delta_x to keep the mesh in the center of texture
+        resolution = config.get_vec3i("resolution");
+        arr = Array3D<char>(resolution.x, resolution.y, resolution.z);
+        real delta_x = 0;
+        Vector3i offset;
+        for (int i = 0; i < 3; ++i)
+            delta_x = std::max(delta_x, (bb.upper_boundary[i] - bb.lower_boundary[i]) / resolution[i]);
+        for (int i = 0; i < 3; ++i)
+            offset[i] = int(resolution[i] - (bb.upper_boundary[i] - bb.lower_boundary[i])/ delta_x) / 2;
+        // cast ray for each i, j
+        for (int i = 0; i < resolution.x; ++i)
+            for (int j = 0; j < resolution.y; ++j) {
+                real x = (i - offset.x) * delta_x + bb.lower_boundary.x;
+                real y = (j - offset.y) * delta_x + bb.lower_boundary.y;
+                real z = (0 - offset.z) * delta_x + bb.lower_boundary.z;
+                real base_z = z - eps;
+                int k = 0;
+                bool inside = false;
+                while (k < resolution.z) {
+                    int next_k;
+                    Ray ray(Vector3(x, y, z), Vector3(0, 0, 1));
+                    scene_geometry->query(ray);
+                    if (ray.dist == Ray::DIST_INFINITE) {
+                        next_k = resolution.z;
+                        inside = false;
+                    }
+                    else {
+                        z += ray.dist;
+                        next_k = std::min(int((z - base_z) / delta_x), resolution.z);
+                    }
+                    while (k < next_k) {
+                        arr.set(i, j, k, char(inside));
+                        ++k;
+                    }
+                    inside = !inside;
+                }
+            }
+    }
+
+    bool inside(const Vector3 &coord) const {
+        return
+            0 < coord.x && coord.x < 1.f &&
+            0 < coord.y && coord.y < 1.f &&
+            0 < coord.z && coord.z < 1.f;
+    }
+
+    virtual Vector4 sample(const Vector3 &coord) const override {
+        if (inside(coord))
+            return Vector4(real(arr.sample_relative_coord(coord)));
+        else
+            return Vector4(0);
+    }
+};
+
+TC_IMPLEMENTATION(Texture, MeshTexture, "mesh")
+
 TC_NAMESPACE_END