/* SPDX-License-Identifier: Apache-2.0 * Copyright 2020-2022 Blender Foundation */ #include "scene/volume.h" #include "scene/attribute.h" #include "scene/image_vdb.h" #include "scene/scene.h" #ifdef WITH_OPENVDB # include # include # include # include #endif #include "util/hash.h" #include "util/log.h" #include "util/openvdb.h" #include "util/progress.h" #include "util/types.h" CCL_NAMESPACE_BEGIN NODE_DEFINE(Volume) { NodeType *type = NodeType::add("volume", create, NodeType::NONE, Mesh::get_node_type()); SOCKET_FLOAT(clipping, "Clipping", 0.001f); SOCKET_FLOAT(step_size, "Step Size", 0.0f); SOCKET_BOOLEAN(object_space, "Object Space", false); SOCKET_FLOAT(velocity_scale, "Velocity Scale", 1.0f); return type; } Volume::Volume() : Mesh(get_node_type(), Geometry::VOLUME) { clipping = 0.001f; step_size = 0.0f; object_space = false; } void Volume::clear(bool preserve_shaders) { Mesh::clear(preserve_shaders, true); } struct QuadData { int v0, v1, v2, v3; float3 normal; }; enum { QUAD_X_MIN = 0, QUAD_X_MAX = 1, QUAD_Y_MIN = 2, QUAD_Y_MAX = 3, QUAD_Z_MIN = 4, QUAD_Z_MAX = 5, }; #ifdef WITH_OPENVDB const int quads_indices[6][4] = { /* QUAD_X_MIN */ {4, 0, 3, 7}, /* QUAD_X_MAX */ {1, 5, 6, 2}, /* QUAD_Y_MIN */ {4, 5, 1, 0}, /* QUAD_Y_MAX */ {3, 2, 6, 7}, /* QUAD_Z_MIN */ {0, 1, 2, 3}, /* QUAD_Z_MAX */ {5, 4, 7, 6}, }; const float3 quads_normals[6] = { /* QUAD_X_MIN */ make_float3(-1.0f, 0.0f, 0.0f), /* QUAD_X_MAX */ make_float3(1.0f, 0.0f, 0.0f), /* QUAD_Y_MIN */ make_float3(0.0f, -1.0f, 0.0f), /* QUAD_Y_MAX */ make_float3(0.0f, 1.0f, 0.0f), /* QUAD_Z_MIN */ make_float3(0.0f, 0.0f, -1.0f), /* QUAD_Z_MAX */ make_float3(0.0f, 0.0f, 1.0f), }; static int add_vertex(int3 v, vector &vertices, int3 res, unordered_map &used_verts) { size_t vert_key = v.x + v.y * (res.x + 1) + v.z * (res.x + 1) * (res.y + 1); unordered_map::iterator it = used_verts.find(vert_key); if (it != used_verts.end()) { return it->second; } int vertex_offset = vertices.size(); used_verts[vert_key] = vertex_offset; vertices.push_back(v); return vertex_offset; } static void create_quad(int3 corners[8], vector &vertices, vector &quads, int3 res, unordered_map &used_verts, int face_index) { QuadData quad; quad.v0 = add_vertex(corners[quads_indices[face_index][0]], vertices, res, used_verts); quad.v1 = add_vertex(corners[quads_indices[face_index][1]], vertices, res, used_verts); quad.v2 = add_vertex(corners[quads_indices[face_index][2]], vertices, res, used_verts); quad.v3 = add_vertex(corners[quads_indices[face_index][3]], vertices, res, used_verts); quad.normal = quads_normals[face_index]; quads.push_back(quad); } #endif /* Create a mesh from a volume. * * The way the algorithm works is as follows: * * - The topologies of input OpenVDB grids are merged into a temporary grid. * - Voxels of the temporary grid are dilated to account for the padding necessary for volume * sampling. * - Quads are created on the boundary between active and inactive leaf nodes of the temporary * grid. */ class VolumeMeshBuilder { public: #ifdef WITH_OPENVDB /* use a MaskGrid to store the topology to save memory */ openvdb::MaskGrid::Ptr topology_grid; openvdb::CoordBBox bbox; #endif bool first_grid; VolumeMeshBuilder(); #ifdef WITH_OPENVDB void add_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, float volume_clipping); #endif void add_padding(int pad_size); void create_mesh(vector &vertices, vector &indices, vector &face_normals, const float face_overlap_avoidance); void generate_vertices_and_quads(vector &vertices_is, vector &quads); void convert_object_space(const vector &vertices, vector &out_vertices, const float face_overlap_avoidance); void convert_quads_to_tris(const vector &quads, vector &tris, vector &face_normals); bool empty_grid() const; #ifdef WITH_OPENVDB template void merge_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, float volume_clipping) { typename GridType::ConstPtr typed_grid = openvdb::gridConstPtrCast(grid); if (do_clipping) { using ValueType = typename GridType::ValueType; typename GridType::Ptr copy = typed_grid->deepCopy(); typename GridType::ValueOnIter iter = copy->beginValueOn(); for (; iter; ++iter) { if (openvdb::math::Abs(iter.getValue()) < ValueType(volume_clipping)) { iter.setValueOff(); } } typed_grid = copy; } topology_grid->topologyUnion(*typed_grid); } #endif }; VolumeMeshBuilder::VolumeMeshBuilder() { first_grid = true; } #ifdef WITH_OPENVDB void VolumeMeshBuilder::add_grid(openvdb::GridBase::ConstPtr grid, bool do_clipping, float volume_clipping) { /* set the transform of our grid from the first one */ if (first_grid) { topology_grid = openvdb::MaskGrid::create(); topology_grid->setTransform(grid->transform().copy()); first_grid = false; } /* if the transforms do not match, we need to resample one of the grids so that * its index space registers with that of the other, here we resample our mask * grid so memory usage is kept low */ else if (topology_grid->transform() != grid->transform()) { openvdb::MaskGrid::Ptr temp_grid = topology_grid->copyWithNewTree(); temp_grid->setTransform(grid->transform().copy()); openvdb::tools::resampleToMatch(*topology_grid, *temp_grid); topology_grid = temp_grid; topology_grid->setTransform(grid->transform().copy()); } if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { merge_grid(grid, do_clipping, volume_clipping); } else if (grid->isType()) { topology_grid->topologyUnion(*openvdb::gridConstPtrCast(grid)); } } #endif void VolumeMeshBuilder::add_padding(int pad_size) { #ifdef WITH_OPENVDB openvdb::tools::dilateActiveValues( topology_grid->tree(), pad_size, openvdb::tools::NN_FACE, openvdb::tools::IGNORE_TILES); #else (void)pad_size; #endif } void VolumeMeshBuilder::create_mesh(vector &vertices, vector &indices, vector &face_normals, const float face_overlap_avoidance) { #ifdef WITH_OPENVDB /* We create vertices in index space (is), and only convert them to object * space when done. */ vector vertices_is; vector quads; /* make sure we only have leaf nodes in the tree, as tiles are not handled by * this algorithm */ topology_grid->tree().voxelizeActiveTiles(); generate_vertices_and_quads(vertices_is, quads); convert_object_space(vertices_is, vertices, face_overlap_avoidance); convert_quads_to_tris(quads, indices, face_normals); #else (void)vertices; (void)indices; (void)face_normals; (void)face_overlap_avoidance; #endif } #ifdef WITH_OPENVDB static bool is_non_empty_leaf(const openvdb::MaskGrid::TreeType &tree, const openvdb::Coord coord) { auto *leaf_node = tree.probeLeaf(coord); return (leaf_node && !leaf_node->isEmpty()); } #endif void VolumeMeshBuilder::generate_vertices_and_quads(vector &vertices_is, vector &quads) { #ifdef WITH_OPENVDB const openvdb::MaskGrid::TreeType &tree = topology_grid->tree(); tree.evalLeafBoundingBox(bbox); const int3 resolution = make_int3(bbox.dim().x(), bbox.dim().y(), bbox.dim().z()); unordered_map used_verts; for (auto iter = tree.cbeginLeaf(); iter; ++iter) { if (iter->isEmpty()) { continue; } openvdb::CoordBBox leaf_bbox = iter->getNodeBoundingBox(); /* +1 to convert from exclusive to include bounds. */ leaf_bbox.max() = leaf_bbox.max().offsetBy(1); int3 min = make_int3(leaf_bbox.min().x(), leaf_bbox.min().y(), leaf_bbox.min().z()); int3 max = make_int3(leaf_bbox.max().x(), leaf_bbox.max().y(), leaf_bbox.max().z()); int3 corners[8] = { make_int3(min[0], min[1], min[2]), make_int3(max[0], min[1], min[2]), make_int3(max[0], max[1], min[2]), make_int3(min[0], max[1], min[2]), make_int3(min[0], min[1], max[2]), make_int3(max[0], min[1], max[2]), make_int3(max[0], max[1], max[2]), make_int3(min[0], max[1], max[2]), }; /* Only create a quad if on the border between an active and an inactive leaf. * * We verify that a leaf exists by probing a coordinate that is at its center, * to do so we compute the center of the current leaf and offset this coordinate * by the size of a leaf in each direction. */ static const int LEAF_DIM = openvdb::MaskGrid::TreeType::LeafNodeType::DIM; auto center = leaf_bbox.min() + openvdb::Coord(LEAF_DIM / 2); if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() - LEAF_DIM, center.y(), center.z()))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MIN); } if (!is_non_empty_leaf(tree, openvdb::Coord(center.x() + LEAF_DIM, center.y(), center.z()))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_X_MAX); } if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() - LEAF_DIM, center.z()))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MIN); } if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y() + LEAF_DIM, center.z()))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Y_MAX); } if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() - LEAF_DIM))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MIN); } if (!is_non_empty_leaf(tree, openvdb::Coord(center.x(), center.y(), center.z() + LEAF_DIM))) { create_quad(corners, vertices_is, quads, resolution, used_verts, QUAD_Z_MAX); } } #else (void)vertices_is; (void)quads; #endif } void VolumeMeshBuilder::convert_object_space(const vector &vertices, vector &out_vertices, const float face_overlap_avoidance) { #ifdef WITH_OPENVDB /* compute the offset for the face overlap avoidance */ bbox = topology_grid->evalActiveVoxelBoundingBox(); openvdb::Coord dim = bbox.dim(); float3 cell_size = make_float3(1.0f / dim.x(), 1.0f / dim.y(), 1.0f / dim.z()); float3 point_offset = cell_size * face_overlap_avoidance; out_vertices.reserve(vertices.size()); for (size_t i = 0; i < vertices.size(); ++i) { openvdb::math::Vec3d p = topology_grid->indexToWorld( openvdb::math::Vec3d(vertices[i].x, vertices[i].y, vertices[i].z)); float3 vertex = make_float3((float)p.x(), (float)p.y(), (float)p.z()); out_vertices.push_back(vertex + point_offset); } #else (void)vertices; (void)out_vertices; (void)face_overlap_avoidance; #endif } void VolumeMeshBuilder::convert_quads_to_tris(const vector &quads, vector &tris, vector &face_normals) { int index_offset = 0; tris.resize(quads.size() * 6); face_normals.reserve(quads.size() * 2); for (size_t i = 0; i < quads.size(); ++i) { tris[index_offset++] = quads[i].v0; tris[index_offset++] = quads[i].v2; tris[index_offset++] = quads[i].v1; face_normals.push_back(quads[i].normal); tris[index_offset++] = quads[i].v0; tris[index_offset++] = quads[i].v3; tris[index_offset++] = quads[i].v2; face_normals.push_back(quads[i].normal); } } bool VolumeMeshBuilder::empty_grid() const { #ifdef WITH_OPENVDB return !topology_grid || topology_grid->tree().leafCount() == 0; #else return true; #endif } #ifdef WITH_OPENVDB template static openvdb::GridBase::ConstPtr openvdb_grid_from_device_texture(device_texture *image_memory, float volume_clipping, Transform transform_3d) { using ValueType = typename GridType::ValueType; openvdb::CoordBBox dense_bbox(0, 0, 0, image_memory->data_width - 1, image_memory->data_height - 1, image_memory->data_depth - 1); typename GridType::Ptr sparse = GridType::create(ValueType(0.0f)); if (dense_bbox.empty()) { return sparse; } openvdb::tools::Dense dense( dense_bbox, static_cast(image_memory->host_pointer)); openvdb::tools::copyFromDense(dense, *sparse, ValueType(volume_clipping)); /* #copyFromDense will remove any leaf node that contains constant data and replace it with a * tile, however, we need to preserve the leaves in order to generate the mesh, so re-voxelize * the leaves that were pruned. This should not affect areas that were skipped due to the * volume_clipping parameter. */ sparse->tree().voxelizeActiveTiles(); /* Compute index to world matrix. */ float3 voxel_size = make_float3(1.0f / image_memory->data_width, 1.0f / image_memory->data_height, 1.0f / image_memory->data_depth); transform_3d = transform_inverse(transform_3d); openvdb::Mat4R index_to_world_mat((double)(voxel_size.x * transform_3d[0][0]), 0.0, 0.0, 0.0, 0.0, (double)(voxel_size.y * transform_3d[1][1]), 0.0, 0.0, 0.0, 0.0, (double)(voxel_size.z * transform_3d[2][2]), 0.0, (double)transform_3d[0][3], (double)transform_3d[1][3], (double)transform_3d[2][3], 1.0); openvdb::math::Transform::Ptr index_to_world_tfm = openvdb::math::Transform::createLinearTransform(index_to_world_mat); sparse->setTransform(index_to_world_tfm); return sparse; } static int estimate_required_velocity_padding(openvdb::GridBase::ConstPtr grid, float velocity_scale) { /* TODO: we may need to also find outliers and clamp them to avoid adding too much padding. */ openvdb::math::Extrema extrema; openvdb::Vec3d voxel_size; /* External .vdb files have a vec3 type for velocity, but the Blender exporter creates a vec4. */ if (grid->isType()) { openvdb::Vec3fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast(grid); extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn()); voxel_size = vel_grid->voxelSize(); } else if (grid->isType()) { openvdb::Vec4fGrid::ConstPtr vel_grid = openvdb::gridConstPtrCast(grid); extrema = openvdb::tools::extrema(vel_grid->cbeginValueOn()); voxel_size = vel_grid->voxelSize(); } else { assert(0); return 0; } /* We should only have uniform grids, so x = y = z, but we never know. */ const double max_voxel_size = openvdb::math::Max(voxel_size.x(), voxel_size.y(), voxel_size.z()); if (max_voxel_size == 0.0) { return 0; } const double estimated_padding = extrema.max() * static_cast(velocity_scale) / max_voxel_size; return static_cast(std::ceil(estimated_padding)); } static openvdb::FloatGrid::ConstPtr get_vdb_for_attribute(Volume *volume, AttributeStandard std) { Attribute *attr = volume->attributes.find(std); if (!attr) { return nullptr; } ImageHandle &handle = attr->data_voxel(); VDBImageLoader *vdb_loader = handle.vdb_loader(); if (!vdb_loader) { return nullptr; } openvdb::GridBase::ConstPtr grid = vdb_loader->get_grid(); if (!grid) { return nullptr; } if (!grid->isType()) { return nullptr; } return openvdb::gridConstPtrCast(grid); } class MergeScalarGrids { typedef openvdb::FloatTree ScalarTree; openvdb::tree::ValueAccessor m_acc_x, m_acc_y, m_acc_z; public: MergeScalarGrids(const ScalarTree *x_tree, const ScalarTree *y_tree, const ScalarTree *z_tree) : m_acc_x(*x_tree), m_acc_y(*y_tree), m_acc_z(*z_tree) { } MergeScalarGrids(const MergeScalarGrids &other) : m_acc_x(other.m_acc_x), m_acc_y(other.m_acc_y), m_acc_z(other.m_acc_z) { } void operator()(const openvdb::Vec3STree::ValueOnIter &it) const { using namespace openvdb; const math::Coord xyz = it.getCoord(); float x = m_acc_x.getValue(xyz); float y = m_acc_y.getValue(xyz); float z = m_acc_z.getValue(xyz); it.setValue(math::Vec3s(x, y, z)); } }; static void merge_scalar_grids_for_velocity(const Scene *scene, Volume *volume) { if (volume->attributes.find(ATTR_STD_VOLUME_VELOCITY)) { /* A vector grid for velocity is already available. */ return; } openvdb::FloatGrid::ConstPtr vel_x_grid = get_vdb_for_attribute(volume, ATTR_STD_VOLUME_VELOCITY_X); openvdb::FloatGrid::ConstPtr vel_y_grid = get_vdb_for_attribute(volume, ATTR_STD_VOLUME_VELOCITY_Y); openvdb::FloatGrid::ConstPtr vel_z_grid = get_vdb_for_attribute(volume, ATTR_STD_VOLUME_VELOCITY_Z); if (!(vel_x_grid && vel_y_grid && vel_z_grid)) { return; } openvdb::Vec3fGrid::Ptr vecgrid = openvdb::Vec3SGrid::create(openvdb::Vec3s(0.0f)); /* Activate voxels in the vector grid based on the scalar grids to ensure thread safety during * the merge. */ vecgrid->tree().topologyUnion(vel_x_grid->tree()); vecgrid->tree().topologyUnion(vel_y_grid->tree()); vecgrid->tree().topologyUnion(vel_z_grid->tree()); MergeScalarGrids op(&vel_x_grid->tree(), &vel_y_grid->tree(), &vel_z_grid->tree()); openvdb::tools::foreach (vecgrid->beginValueOn(), op, true, false); /* Assume all grids have the same transformation. */ openvdb::math::Transform::Ptr transform = openvdb::ConstPtrCast( vel_x_grid->transformPtr()); vecgrid->setTransform(transform); /* Make an attribute for it. */ Attribute *attr = volume->attributes.add(ATTR_STD_VOLUME_VELOCITY); ImageLoader *loader = new VDBImageLoader(vecgrid, "merged_velocity"); ImageParams params; attr->data_voxel() = scene->image_manager->add_image(loader, params); } #endif /* ************************************************************************** */ void GeometryManager::create_volume_mesh(const Scene *scene, Volume *volume, Progress &progress) { string msg = string_printf("Computing Volume Mesh %s", volume->name.c_str()); progress.set_status("Updating Mesh", msg); /* Find shader and compute padding based on volume shader interpolation settings. */ Shader *volume_shader = NULL; int pad_size = 0; for (Node *node : volume->get_used_shaders()) { Shader *shader = static_cast(node); if (!shader->has_volume) { continue; } volume_shader = shader; if (shader->get_volume_interpolation_method() == VOLUME_INTERPOLATION_LINEAR) { pad_size = max(1, pad_size); } else if (shader->get_volume_interpolation_method() == VOLUME_INTERPOLATION_CUBIC) { pad_size = max(2, pad_size); } break; } /* Clear existing volume mesh, done here in case we early out due to * empty grid or missing volume shader. * Also keep the shaders to avoid infinite loops when synchronizing, as this will tag the shaders * as having changed. */ volume->clear(true); volume->need_update_rebuild = true; if (!volume_shader) { return; } /* Create volume mesh builder. */ VolumeMeshBuilder builder; #ifdef WITH_OPENVDB merge_scalar_grids_for_velocity(scene, volume); for (Attribute &attr : volume->attributes.attributes) { if (attr.element != ATTR_ELEMENT_VOXEL) { continue; } bool do_clipping = false; ImageHandle &handle = attr.data_voxel(); /* Try building from OpenVDB grid directly. */ VDBImageLoader *vdb_loader = handle.vdb_loader(); openvdb::GridBase::ConstPtr grid; if (vdb_loader) { grid = vdb_loader->get_grid(); /* If building from an OpenVDB grid, we need to manually clip the values. */ do_clipping = true; } /* Else fall back to creating an OpenVDB grid from the dense volume data. */ if (!grid) { device_texture *image_memory = handle.image_memory(); if (image_memory->data_elements == 1) { grid = openvdb_grid_from_device_texture( image_memory, volume->get_clipping(), handle.metadata().transform_3d); } else if (image_memory->data_elements == 3) { grid = openvdb_grid_from_device_texture( image_memory, volume->get_clipping(), handle.metadata().transform_3d); } else if (image_memory->data_elements == 4) { grid = openvdb_grid_from_device_texture( image_memory, volume->get_clipping(), handle.metadata().transform_3d); } } if (grid) { /* Add padding based on the maximum velocity vector. */ if (attr.std == ATTR_STD_VOLUME_VELOCITY && scene->need_motion() != Scene::MOTION_NONE) { pad_size = max(pad_size, estimate_required_velocity_padding(grid, volume->get_velocity_scale())); } builder.add_grid(grid, do_clipping, volume->get_clipping()); } } #else (void)scene; #endif /* If nothing to build, early out. */ if (builder.empty_grid()) { return; } builder.add_padding(pad_size); /* Slightly offset vertex coordinates to avoid overlapping faces with other * volumes or meshes. The proper solution would be to improve intersection in * the kernel to support robust handling of multiple overlapping faces or use * an all-hit intersection similar to shadows. */ const float face_overlap_avoidance = 0.1f * hash_uint_to_float(hash_string(volume->name.c_str())); /* Create mesh. */ vector vertices; vector indices; vector face_normals; builder.create_mesh(vertices, indices, face_normals, face_overlap_avoidance); volume->reserve_mesh(vertices.size(), indices.size() / 3); volume->used_shaders.clear(); volume->used_shaders.push_back_slow(volume_shader); for (size_t i = 0; i < vertices.size(); ++i) { volume->add_vertex(vertices[i]); } for (size_t i = 0; i < indices.size(); i += 3) { volume->add_triangle(indices[i], indices[i + 1], indices[i + 2], 0, false); } Attribute *attr_fN = volume->attributes.add(ATTR_STD_FACE_NORMAL); float3 *fN = attr_fN->data_float3(); for (size_t i = 0; i < face_normals.size(); ++i) { fN[i] = face_normals[i]; } /* Print stats. */ VLOG_WORK << "Memory usage volume mesh: " << ((vertices.size() + face_normals.size()) * sizeof(float3) + indices.size() * sizeof(int)) / (1024.0 * 1024.0) << "Mb."; } CCL_NAMESPACE_END