From ed4e61a8a4b0700e98ee66cf2877e63043fa55b0 Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Mon, 8 Apr 2024 14:05:59 -0700 Subject: [PATCH 1/8] Initial draft seg faults --- cpp/open3d/core/TensorKey.h | 4 +- cpp/open3d/t/geometry/RaycastingScene.cpp | 4 +- cpp/open3d/t/geometry/RaycastingScene.h | 2 +- cpp/open3d/t/geometry/TriangleMesh.cpp | 241 ++++++++++++++++++- cpp/open3d/t/geometry/TriangleMesh.h | 21 +- cpp/open3d/t/geometry/Utility.h | 2 +- cpp/open3d/t/geometry/VoxelBlockGrid.cpp | 4 +- cpp/open3d/t/geometry/kernel/IPPImage.cpp | 78 +++++- cpp/open3d/t/geometry/kernel/IPPImage.h | 34 +-- cpp/open3d/visualization/gui/Application.cpp | 2 + cpp/pybind/t/geometry/trianglemesh.cpp | 19 +- cpp/tests/t/geometry/TriangleMesh.cpp | 37 +++ 12 files changed, 409 insertions(+), 39 deletions(-) diff --git a/cpp/open3d/core/TensorKey.h b/cpp/open3d/core/TensorKey.h index 3c90b10cf59..eb203baea73 100644 --- a/cpp/open3d/core/TensorKey.h +++ b/cpp/open3d/core/TensorKey.h @@ -82,11 +82,11 @@ class TensorKey { /// For TensorKeyMode::Slice only. int64_t GetStart() const; - /// Get stop index. Throws exception if start is None. + /// Get stop index. Throws exception if stop is None. /// For TensorKeyMode::Slice only. int64_t GetStop() const; - /// Get step index. Throws exception if start is None. + /// Get step index. Throws exception if step is None. /// For TensorKeyMode::Slice only. int64_t GetStep() const; diff --git a/cpp/open3d/t/geometry/RaycastingScene.cpp b/cpp/open3d/t/geometry/RaycastingScene.cpp index 14f9962c26c..8040d327978 100644 --- a/cpp/open3d/t/geometry/RaycastingScene.cpp +++ b/cpp/open3d/t/geometry/RaycastingScene.cpp @@ -777,7 +777,7 @@ uint32_t RaycastingScene::AddTriangles(const TriangleMesh& mesh) { } std::unordered_map RaycastingScene::CastRays( - const core::Tensor& rays, const int nthreads) { + const core::Tensor& rays, const int nthreads) const { AssertTensorDtypeLastDimDeviceMinNDim(rays, "rays", 6, impl_->tensor_device_); auto shape = rays.GetShape(); @@ -1173,4 +1173,4 @@ uint32_t RaycastingScene::INVALID_ID() { return RTC_INVALID_GEOMETRY_ID; } } // namespace geometry } // namespace t -} // namespace open3d \ No newline at end of file +} // namespace open3d diff --git a/cpp/open3d/t/geometry/RaycastingScene.h b/cpp/open3d/t/geometry/RaycastingScene.h index f25d994b0b5..9e3e83f0c1e 100644 --- a/cpp/open3d/t/geometry/RaycastingScene.h +++ b/cpp/open3d/t/geometry/RaycastingScene.h @@ -70,7 +70,7 @@ class RaycastingScene { /// - \b primitive_normals A tensor with the normals of the hit /// triangles. The shape is {.., 3}. std::unordered_map CastRays( - const core::Tensor &rays, const int nthreads = 0); + const core::Tensor &rays, const int nthreads = 0) const; /// \brief Checks if the rays have any intersection with the scene. /// \param rays A tensor with >=2 dims, shape {.., 6}, and Dtype Float32 diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 2bf329838cc..b7bc2e641b4 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -16,23 +16,36 @@ #include #include +#include +#include #include +#include #include +#include +#include #include "open3d/core/CUDAUtils.h" +#include "open3d/core/Device.h" +#include "open3d/core/Dtype.h" #include "open3d/core/EigenConverter.h" #include "open3d/core/ShapeUtil.h" #include "open3d/core/Tensor.h" #include "open3d/core/TensorCheck.h" +#include "open3d/core/TensorKey.h" +#include "open3d/core/linalg/AddMM.h" +#include "open3d/core/linalg/Matmul.h" #include "open3d/t/geometry/LineSet.h" #include "open3d/t/geometry/PointCloud.h" #include "open3d/t/geometry/RaycastingScene.h" #include "open3d/t/geometry/VtkUtils.h" +#include "open3d/t/geometry/kernel/IPPImage.h" #include "open3d/t/geometry/kernel/PCAPartition.h" #include "open3d/t/geometry/kernel/PointCloud.h" #include "open3d/t/geometry/kernel/Transform.h" #include "open3d/t/geometry/kernel/TriangleMesh.h" #include "open3d/t/geometry/kernel/UVUnwrapping.h" +#include "open3d/t/io/ImageIO.h" +#include "open3d/utility/Optional.h" #include "open3d/utility/ParallelScan.h" namespace open3d { @@ -972,7 +985,7 @@ TriangleMesh::BakeTriangleAttrTextures( } core::Tensor tensor = triangle_attr_.at(attr).To(core::Device()).Contiguous(); - DISPATCH_DTYPE_TO_TEMPLATE(tensor.GetDtype(), [&]() { + DISPATCH_DTYPE_TO_TEMPLATE_WITH_BOOL(tensor.GetDtype(), [&]() { core::Tensor tex; if (GetTriangleIndices().GetDtype() == core::Int32) { tex = BakeAttribute( @@ -1149,7 +1162,8 @@ static bool IsNegative(T val) { return false; } -TriangleMesh TriangleMesh::SelectByIndex(const core::Tensor &indices) const { +TriangleMesh TriangleMesh::SelectByIndex(const core::Tensor &indices, + bool copy_attributes) const { core::AssertTensorShape(indices, {indices.GetLength()}); if (indices.NumElements() == 0) { return {}; @@ -1249,7 +1263,8 @@ TriangleMesh TriangleMesh::SelectByIndex(const core::Tensor &indices) const { if (tris_cpu.NumElements() > 0) { // To() needs non-empty tensor result.SetTriangleIndices(tris_cpu.To(GetDevice())); } - CopyAttributesByMasks(result, *this, vertex_mask, tri_mask); + if (copy_attributes) + CopyAttributesByMasks(result, *this, vertex_mask, tri_mask); return result; } @@ -1304,11 +1319,229 @@ TriangleMesh TriangleMesh::RemoveUnreferencedVertices() { utility::LogDebug( "[RemoveUnreferencedVertices] {:d} vertices have been removed.", - (int)(num_verts_old - GetVertexPositions().GetLength())); + static_cast(num_verts_old - GetVertexPositions().GetLength())); return *this; } +namespace { + +core::Tensor Project( + const core::Tensor &t_xyz, // contiguous {...,3} + const core::Tensor &t_intrinsic_matrix, // contiguous {3,3} + const core::Tensor &t_extrinsic_matrix) { // contiguous {4,4} + auto xy_shape = t_xyz.GetShape(); + xy_shape[xy_shape.GetLength() - 1] = 2; + core::Tensor t_xy(xy_shape, t_xyz.GetDtype()); + DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(t_xyz.GetDtype(), [&]() { + Eigen::Map> xy(t_xy.GetDataPtr(), 2, + t_xy.NumElements() / 2); + Eigen::Map> xyz( + t_xyz.GetDataPtr(), 3, t_xyz.NumElements() / 3); + Eigen::Map> intrinsic_matrix( + t_intrinsic_matrix.GetDataPtr()); + Eigen::Map> extrinsic_matrix( + t_extrinsic_matrix.GetDataPtr()); + + Eigen::ArrayX pxyz = + intrinsic_matrix * + (extrinsic_matrix.topLeftCorner<3, 3>() * xyz + + extrinsic_matrix.topRightCorner<3, 1>()); // {...,3} + pxyz.leftCols<2>().colwise() /= pxyz.rightCols<1>(); + xy = pxyz.leftCols<2>(); + }); + return t_xy; // contiguous {...,2} +} +} // namespace + +std::pair TriangleMesh::SelectVisible( + const core::Tensor &intrinsic_matrix, + const core::Tensor &extrinsic_matrix, + int width_px, + int height_px, + const RaycastingScene &rcs) { + using tk = core::TensorKey; + using core::None; + constexpr double EPS = 1e-6; + const core::Tensor vertices = GetVertexPositions().To(core::Device()); + core::Tensor imxy = Project(vertices, intrinsic_matrix, extrinsic_matrix); + core::Tensor vis_vert = + core::Tensor::Empty({vertices.GetLength()}, core::Bool); + bool *vis_vert_ptr = vis_vert.GetDataPtr(); + int n_vis = 0; + DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(vertices.GetDtype(), [&]() { + scalar_t *imxy_ptr = imxy.GetDataPtr(); + for (int i = 0; i < vertices.GetLength(); ++i, imxy_ptr += 3) { + vis_vert_ptr[i] = (imxy_ptr[0] > 0 && imxy_ptr[0] < width_px && + imxy_ptr[1] > 0 && imxy_ptr[1] < height_px); + n_vis += vis_vert_ptr[i]; + } + }); + core::Tensor cam_loc = + -extrinsic_matrix + .GetItem({tk::Slice(0, 3, None), tk::Slice(0, 3, None)}) + .T() + .Matmul(extrinsic_matrix.GetItem( + {tk::Slice(0, 3, None), tk::Index(3)})); + utility::LogInfo("cam_loc: {}", cam_loc.ToString()); + core::Tensor rays = core::Tensor::Empty({n_vis, 6}, vertices.GetDtype()); + rays.SetItem({tk::Slice(None, None, None), tk::Slice(0, 3, 1)}, + cam_loc.T()); + core::Tensor vis_vert_idx(vis_vert.NonZero().Reshape({-1})); + rays.SetItem({tk::Slice(None, None, None), tk::Slice(3, None, 1)}, + vertices.GetItem({tk::IndexTensor(vis_vert_idx)})); + + auto result = rcs.CastRays(rays); + float *t_hit_ptr = result["t_hit"].GetDataPtr(); + int64_t *vis_vert_idx_ptr = vis_vert_idx.GetDataPtr(); + for (int i = 0; i < n_vis; ++i) { + vis_vert_ptr[vis_vert_idx_ptr[i]] = t_hit_ptr[i] > 1 - EPS; + } + core::Tensor vis_tris = core::Tensor::Zeros( + {GetTriangleIndices().GetShape(0)}, core::Dtype::Bool); + DISPATCH_INT_DTYPE_PREFIX_TO_TEMPLATE( + GetTriangleIndices().GetDtype(), tri, [&]() { + auto tris_ptr = GetTriangleIndices() + .To(core::Device()) + .Contiguous() + .GetDataPtr(); + auto vis_tris_ptr = vis_tris.GetDataPtr(); + for (int i = 0; i < GetTriangleIndices().GetShape(0); ++i) { + vis_tris_ptr[i] = (vis_vert_ptr[tris_ptr[3 * i + 0]] && + vis_vert_ptr[tris_ptr[3 * i + 1]] && + vis_vert_ptr[tris_ptr[3 * i + 2]]); + } + }); + return std::make_pair(vis_vert, vis_tris); +} + +Image TriangleMesh::ProjectImagesToAlbedo( + const std::vector &images, + const std::vector &intrinsic_matrices, + const std::vector &extrinsic_matrices, + int tex_size /*=1024*/, + bool update_material /*=true*/) { + using core::None; + using tk = core::TensorKey; + constexpr double EPS = 1e-6; + if (!triangle_attr_.Contains("texture_uvs")) { + utility::LogError("Cannot find triangle attribute 'texture_uvs'"); + } + core::Tensor texture_uvs = + triangle_attr_.at("texture_uvs").To(core::Device()).Contiguous(); + core::AssertTensorShape(texture_uvs, {core::None, 3, 2}); + core::AssertTensorDtype(texture_uvs, {core::Float32}); + + if (images.size() != extrinsic_matrices.size() || + images.size() != intrinsic_matrices.size()) { + utility::LogError( + "Received {} images, but {} extrinsic matrices and {} " + "intrinsic matrices.", + images.size(), extrinsic_matrices.size(), + intrinsic_matrices.size()); + } + // (u,v) -> (x,y,z) : {tex_size, tex_size, 3} + core::Tensor position_map = BakeVertexAttrTextures( + tex_size, {"positions"}, 1, 0, false)["positions"]; + core::Tensor albedo = + core::Tensor::Zeros({tex_size, tex_size, 4}, core::Float32); + albedo.Slice(2, 3, 4, 1).Fill(EPS); // regularize + core::Tensor this_albedo = + core::Tensor::Zeros({tex_size, tex_size, 4}, core::Float32); + core::Tensor weighted_image; + + /* For each image, project vertices onto 2D image and do occlusion + * culling. This creates (xim, yim) -> (u,v) mapping. Invert this + * mapping with ipp::remap to get (u, v) -> (xim, yim) -> rgb to create + * albedo texture rgb(u,v) with bilinear interpolation. Merge albedos + * from all images OR (u,v) -> (x, y, z) -> project to (xim, yim) -> rgb + * Discard if occlusion test fails / hit test for ray fails + */ + RaycastingScene rcs; + rcs.AddTriangles(*this); + + for (size_t i = 0; i < images.size(); ++i) { + auto width = images[i].GetCols(), height = images[i].GetRows(); + core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); + core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); + + // A. Get image space weight matrix + auto rays = RaycastingScene::CreateRaysPinhole( + intrinsic_matrices[i], extrinsic_matrices[i], + images[i].GetCols(), images[i].GetRows()); + auto result = rcs.CastRays(rays); + // Eigen is column-major order + Eigen::Map normals_e( + result["primitive_normals"].GetDataPtr(), 3, + width * height); + Eigen::Map rays_e(rays.GetDataPtr(), 6, + width * height); + Eigen::Map t_hit(result["t_hit"].GetDataPtr(), 3, + width * height); + // Weight map based on dot product between ray normals and triangle + // normals. + auto depth = t_hit * (rays_e.bottomRows<3>() - rays_e.topRows<3>()) + .colwise() + .norm() + .array(); + auto footprint = depth * normals_e + .cwiseProduct((rays_e.bottomRows<3>() - + rays_e.topRows<3>()) + .colwise() + .normalized()) + .sum(); + if (!weighted_image.GetShape().IsCompatible({height, width, 4})) { + weighted_image = core::Tensor({height, width, 4}, core::Float32); + } + weighted_image.Slice(2, 0, 3, 1) = images[i].AsTensor(); + Eigen::Map weighted_image_e( + weighted_image.GetDataPtr(), 4, width * height); + // footprint is inf if depth or t_hit is inf. Never zero. + weighted_image_e.bottomRows<3>() = 1. / footprint; + /* weighted_image.Slice(2, 3, 4, 1) = 1. / footprint; */ + + // B. Get texture space (u,v) -> (x,y) map and valid domain in uv space. + core::Tensor uv2xy = Project(position_map, intrinsic_matrices[i], + extrinsic_matrices[i]); // {ts, ts, 2} + core::Tensor visible_triangles = + SelectVisible(intrinsic_matrices[i], extrinsic_matrices[i], + images[i].GetCols(), images[i].GetRows(), rcs) + .second; + SetTriangleAttr("visible", visible_triangles); + core::Tensor visible_map = BakeTriangleAttrTextures( + tex_size, {"visible"}, 2, 0, + false)["visible"]; // uv space {ts, ts, 1} + /* RemoveTriangleAttr("visible"); */ + // do not remap occluded / out of view faces + io::WriteImage("visible_map.png", Image(visible_map.To(core::UInt8))); + uv2xy.Reshape({-1, 2}).GetItem( + tk::IndexTensor(visible_map.Reshape({-1}).NonZero())) = -1.f; + // albedo[u,v] = image[ i[u,v], j[u,v] ] + // ij[u,v] ? = identity[ uv[i,j] ] + + // C. Interpolate weighted image to weighted texture + ipp::Remap(weighted_image /*{height, width, 4} f32*/, + uv2xy /* {texsz, texsz, 2} f32*/, + this_albedo /*{texsz, texsz, 4} f32*/, + t::geometry::Image::InterpType::Cubic); + // weighted sum + albedo.Slice(2, 0, 3, 1) += + this_albedo.Slice(2, 0, 3, 1) * this_albedo.Slice(2, 3, 4, 1); + albedo.Slice(2, 3, 4, 1) += this_albedo.Slice(2, 3, 4, 1); + } + albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); + + Image albedo_texture(albedo.Slice(2, 0, 3, 1)); + if (update_material) { + if (!HasMaterial()) { + SetMaterial(visualization::rendering::Material()); + GetMaterial().SetDefaultProperties(); // defaultUnlit + } + GetMaterial().SetTextureMap("albedo", albedo_texture); + } + return albedo_texture; +} + } // namespace geometry } // namespace t } // namespace open3d diff --git a/cpp/open3d/t/geometry/TriangleMesh.h b/cpp/open3d/t/geometry/TriangleMesh.h index 4a1f9138fd6..a3a2616e931 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.h +++ b/cpp/open3d/t/geometry/TriangleMesh.h @@ -22,6 +22,7 @@ namespace t { namespace geometry { class LineSet; +class RaycastingScene; /// \class TriangleMesh /// \brief A triangle mesh contains vertices and triangles. @@ -940,15 +941,33 @@ class TriangleMesh : public Geometry, public DrawableGeometry { /// the mesh or has a negative value, it is ignored. /// \param indices An integer list of indices. Duplicates are /// allowed, but ignored. Signed and unsigned integral types are allowed. + /// \param copy_attributes Indicates if vertex attributes (other than + /// positions) and triangle attributes (other than indices) should be copied + /// to the returned mesh. /// \return A new mesh with the selected vertices and faces built /// from the selected vertices. If the original mesh is empty, return /// an empty mesh. - TriangleMesh SelectByIndex(const core::Tensor &indices) const; + TriangleMesh SelectByIndex(const core::Tensor &indices, + bool copy_attributes = true) const; /// Removes unreferenced vertices from the mesh. /// \return The reference to itself. TriangleMesh RemoveUnreferencedVertices(); + std::pair SelectVisible( + const core::Tensor &intrinsic_matrix, + const core::Tensor &extrinsic_matrix, + int width_px, + int height_px, + const RaycastingScene &rcs); + + Image ProjectImagesToAlbedo( + const std::vector &images, + const std::vector &intrinsic_matrices, + const std::vector &extrinsic_matrices, + int tex_size = 1024, + bool update_material = true); + protected: core::Device device_ = core::Device("CPU:0"); TensorMap vertex_attr_; diff --git a/cpp/open3d/t/geometry/Utility.h b/cpp/open3d/t/geometry/Utility.h index a8a73ae2895..343d5d77aaf 100644 --- a/cpp/open3d/t/geometry/Utility.h +++ b/cpp/open3d/t/geometry/Utility.h @@ -66,7 +66,7 @@ inline void CheckExtrinsicTensor(const core::Tensor& extrinsic) { } } -inline void CheckBlockCoorinates(const core::Tensor& block_coords) { +inline void CheckBlockCoordinates(const core::Tensor& block_coords) { if (block_coords.GetDtype() != core::Dtype::Int32) { utility::LogError("Unsupported block coordinate dtype {}", block_coords.GetDtype().ToString()); diff --git a/cpp/open3d/t/geometry/VoxelBlockGrid.cpp b/cpp/open3d/t/geometry/VoxelBlockGrid.cpp index b98c09b95f2..535e9423fd9 100644 --- a/cpp/open3d/t/geometry/VoxelBlockGrid.cpp +++ b/cpp/open3d/t/geometry/VoxelBlockGrid.cpp @@ -301,7 +301,7 @@ void VoxelBlockGrid::Integrate(const core::Tensor &block_coords, AssertInitialized(); bool integrate_color = color.AsTensor().NumElements() > 0; - CheckBlockCoorinates(block_coords); + CheckBlockCoordinates(block_coords); CheckDepthTensor(depth.AsTensor()); if (integrate_color) { CheckColorTensor(color.AsTensor()); @@ -338,7 +338,7 @@ TensorMap VoxelBlockGrid::RayCast(const core::Tensor &block_coords, float trunc_voxel_multiplier, int range_map_down_factor) { AssertInitialized(); - CheckBlockCoorinates(block_coords); + CheckBlockCoordinates(block_coords); CheckIntrinsicTensor(intrinsic); CheckExtrinsicTensor(extrinsic); diff --git a/cpp/open3d/t/geometry/kernel/IPPImage.cpp b/cpp/open3d/t/geometry/kernel/IPPImage.cpp index 275d15b560b..e0ba8167401 100644 --- a/cpp/open3d/t/geometry/kernel/IPPImage.cpp +++ b/cpp/open3d/t/geometry/kernel/IPPImage.cpp @@ -7,10 +7,14 @@ #include "open3d/t/geometry/kernel/IPPImage.h" +#include +#include + #include #include #include #include +#include #include "open3d/core/Dtype.h" #include "open3d/core/ParallelFor.h" @@ -76,8 +80,8 @@ void RGBToGray(const core::Tensor &src_im, core::Tensor &dst_im) { } } -void Resize(const open3d::core::Tensor &src_im, - open3d::core::Tensor &dst_im, +void Resize(const core::Tensor &src_im, + core::Tensor &dst_im, t::geometry::Image::InterpType interp_type) { auto dtype = src_im.GetDtype(); // Create IPP wrappers for all Open3D tensors @@ -153,9 +157,9 @@ void Dilate(const core::Tensor &src_im, core::Tensor &dst_im, int kernel_size) { } } -void Filter(const open3d::core::Tensor &src_im, - open3d::core::Tensor &dst_im, - const open3d::core::Tensor &kernel) { +void Filter(const core::Tensor &src_im, + core::Tensor &dst_im, + const core::Tensor &kernel) { // Supported device and datatype checking happens in calling code and will // result in an exception if there are errors. auto dtype = src_im.GetDtype(); @@ -291,7 +295,71 @@ void FilterSobel(const core::Tensor &src_im, utility::LogError("IPP-IW error {}: {}", e.m_status, e.m_string); } } + +void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, + const core::Tensor &dst2src_map /*{Wd, Hd, 2}, float*/, + core::Tensor &dst_im /*{Wd, Hd, 2}*/, + t::geometry::Image::InterpType interp_type) { + auto dtype = src_im.GetDtype(); + if (dtype != dst_im.GetDtype()) { + utility::LogError( + "Source ({}) and destination ({}) image dtypes are different!", + dtype.ToString(), dst_im.GetDtype().ToString()); + } + if (dst2src_map.GetDtype() != core::Float32) { + utility::LogError("dst2src_map dtype ({}) must be Float32.", + dst2src_map.GetDtype().ToString()); + } + + static const std::unordered_map + interp_dict = { + {t::geometry::Image::InterpType::Nearest, IPPI_INTER_NN}, + {t::geometry::Image::InterpType::Linear, IPPI_INTER_LINEAR}, + {t::geometry::Image::InterpType::Cubic, IPPI_INTER_CUBIC}, + {t::geometry::Image::InterpType::Lanczos, + IPPI_INTER_LANCZOS}, + /* {t::geometry::Image::InterpType::Cubic2p_CatmullRom, */ + /* IPPI_INTER_CUBIC2P_CATMULLROM}, */ + }; + + auto interp_it = interp_dict.find(interp_type); + if (interp_it == interp_dict.end()) { + utility::LogError("Unsupported interp type {}", + static_cast(interp_type)); + } + + if (src_im.GetDtype() == core::Float32 && src_im.GetShape(2) == 4) { + /* IPPAPI(IppStatus, ippiRemap_32f_C4R, (const Ipp32f* pSrc, IppiSize + * srcSize, */ + /* int srcStep, IppiRect srcROI, const Ipp32f* pxMap, int xMapStep, + */ + /* const Ipp32f* pyMap, int yMapStep, Ipp32f* pDst, int dstStep, */ + /* IppiSize dstRoiSize, int interpolation)) */ + IppiSize src_size{static_cast(src_im.GetShape(1)), + static_cast(src_im.GetShape(0))}, + dst_roi_size{static_cast(dst_im.GetShape(1)), + static_cast(dst_im.GetShape(0))}; + IppiRect src_roi{0, 0, static_cast(src_im.GetShape(1)), + static_cast(src_im.GetShape(0))}; + + IppStatus sts = ippiRemap_32f_C4R( + src_im.GetDataPtr(), src_size, 1, src_roi, + dst2src_map.GetDataPtr(), 2, + dst2src_map.GetDataPtr() + 1, 2, + dst_im.GetDataPtr(), 1, dst_roi_size, interp_it->second); + if (sts != ippStsNoErr) { + // See comments in icv/include/ippicv_types.h for meaning + utility::LogError("IPP remap error {}", sts); + } + + } else { + utility::LogError( + "Remap not implemented for dtype ({}) and channels ({}).", + src_im.GetDtype().ToString(), src_im.GetShape(2)); + } +} } // namespace ipp + } // namespace geometry } // namespace t } // namespace open3d diff --git a/cpp/open3d/t/geometry/kernel/IPPImage.h b/cpp/open3d/t/geometry/kernel/IPPImage.h index b861dd4d920..b5c7fddab6c 100644 --- a/cpp/open3d/t/geometry/kernel/IPPImage.h +++ b/cpp/open3d/t/geometry/kernel/IPPImage.h @@ -50,33 +50,37 @@ void To(const core::Tensor &src_im, void RGBToGray(const core::Tensor &src_im, core::Tensor &dst_im); -void Dilate(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim, - int kernel_size); +void Dilate(const core::Tensor &srcim, core::Tensor &dstim, int kernel_size); -void Resize(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim, +void Resize(const core::Tensor &srcim, + core::Tensor &dstim, t::geometry::Image::InterpType interp_type); -void Filter(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim, - const open3d::core::Tensor &kernel); +void Filter(const core::Tensor &srcim, + core::Tensor &dstim, + const core::Tensor &kernel); -void FilterBilateral(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim, +void FilterBilateral(const core::Tensor &srcim, + core::Tensor &dstim, int kernel_size, float value_sigma, float distance_sigma); -void FilterGaussian(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim, +void FilterGaussian(const core::Tensor &srcim, + core::Tensor &dstim, int kernel_size, float sigma); -void FilterSobel(const open3d::core::Tensor &srcim, - open3d::core::Tensor &dstim_dx, - open3d::core::Tensor &dstim_dy, +void FilterSobel(const core::Tensor &srcim, + core::Tensor &dstim_dx, + core::Tensor &dstim_dy, int kernel_size); + +void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, + const core::Tensor &dst2src_map /*{Wd, Hd, 2}, float*/, + core::Tensor &dst_im /*{Wd, Hd, 2}*/, + Image::InterpType interp_type); + } // namespace ipp } // namespace geometry } // namespace t diff --git a/cpp/open3d/visualization/gui/Application.cpp b/cpp/open3d/visualization/gui/Application.cpp index 89cfa9746b2..288b0fd2caa 100644 --- a/cpp/open3d/visualization/gui/Application.cpp +++ b/cpp/open3d/visualization/gui/Application.cpp @@ -89,6 +89,8 @@ std::string FindResourcePath(int argc, const char *argv[]) { for (auto &subpath : {"/resources", "/../resources" /*building with Xcode */, "/share/resources" /* GNU */, "/share/Open3D/resources" /* GNU */}) { + open3d::utility::LogInfo("Checking for resources in {}", + path + subpath); if (o3dfs::DirectoryExists(path + subpath)) { return path + subpath; } diff --git a/cpp/pybind/t/geometry/trianglemesh.cpp b/cpp/pybind/t/geometry/trianglemesh.cpp index cf245426d58..5698753faf4 100644 --- a/cpp/pybind/t/geometry/trianglemesh.cpp +++ b/cpp/pybind/t/geometry/trianglemesh.cpp @@ -706,7 +706,7 @@ This function always uses the CPU device. Returns: This function creates a face attribute "texture_uvs" and returns a tuple - with (max stretch, num_charts, num_partitions) storing the + with (max stretch, num_charts, num_partitions) storing the actual amount of stretch, the number of created charts, and the number of parallel partitions created. @@ -836,6 +836,13 @@ This function always uses the CPU device. plt.imshow(texture_tensors['albedo'].numpy()) )"); + triangle_mesh.def( + "project_images_to_albedo", &TriangleMesh::ProjectImagesToAlbedo, + "images"_a, "intrinsic_matrices"_a, "extrinsic_matrices"_a, + "tex_size"_a = 1024, "update_material"_a = true, + py::call_guard(), + R"(Create an albedo texture from images of an object taken with a calibrated camera.)"); + triangle_mesh.def("extrude_rotation", &TriangleMesh::ExtrudeRotation, "angle"_a, "axis"_a, "resolution"_a = 16, "translation"_a = 0.0, "capping"_a = true, @@ -883,7 +890,7 @@ This function always uses the CPU device. "max_faces"_a, R"(Partition the mesh by recursively doing PCA. -This function creates a new face attribute with the name "partition_ids" storing +This function creates a new face attribute with the name "partition_ids" storing the partition id for each face. Args: @@ -892,7 +899,7 @@ the partition id for each face. Example: - This code partitions a mesh such that each partition contains at most 20k + This code partitions a mesh such that each partition contains at most 20k faces:: import open3d as o3d @@ -911,15 +918,15 @@ the partition id for each face. R"(Returns a new mesh with the faces selected by a boolean mask. Args: - mask (open3d.core.Tensor): A boolean mask with the shape (N) with N as the + mask (open3d.core.Tensor): A boolean mask with the shape (N) with N as the number of faces in the mesh. - + Returns: A new mesh with the selected faces. If the original mesh is empty, return an empty mesh. Example: - This code partitions the mesh using PCA and then visualized the individual + This code partitions the mesh using PCA and then visualized the individual parts:: import open3d as o3d diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index b891d6b6986..6198179ede4 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -11,7 +11,12 @@ #include "core/CoreTest.h" #include "open3d/core/EigenConverter.h" +#include "open3d/core/Tensor.h" #include "open3d/core/TensorCheck.h" +#include "open3d/geometry/LineSet.h" +#include "open3d/t/io/ImageIO.h" +#include "open3d/t/io/TriangleMeshIO.h" +#include "open3d/visualization/utility/Draw.h" #include "tests/Tests.h" namespace open3d { @@ -1342,5 +1347,37 @@ TEST_P(TriangleMeshPermuteDevices, RemoveUnreferencedVertices) { EXPECT_TRUE(torus.GetTriangleNormals().AllClose(expected_tri_normals)); } +TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { + using namespace t::geometry; + core::Device device = GetParam(); + TriangleMesh sphere = + TriangleMesh::FromLegacy(*geometry::TriangleMesh::CreateSphere( + 1.0, 10, /*create_uv_map=*/true)); + core::Tensor view = core::Tensor::Zeros({192, 256, 3}, core::Float32); + view.Slice(2, 0, 1, 1).Fill(1.0); // red + core::Tensor intrinsic_matrix = core::Tensor::Init( + {{256, 0, 96}, {0, 256, 128}, {0, 0, 1}}, device); + core::Tensor extrinsic_matrix = core::Tensor::Init( + {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 1}, {0, 0, 0, 1}}, device); + + Eigen::Map e_intrinsic( + intrinsic_matrix.GetDataPtr()); + Eigen::Map e_extrinsic( + extrinsic_matrix.GetDataPtr()); + auto p_camera = open3d::geometry::LineSet::CreateCameraVisualization( + 256, 192, e_intrinsic.cast(), e_extrinsic.cast()); + + Image texture = sphere.ProjectImagesToAlbedo( + {Image(view)}, {intrinsic_matrix}, {extrinsic_matrix}, 256, true); + t::io::WriteImage("texture.png", texture); + t::io::WriteTriangleMesh("sphere-projected.obj", sphere); + + visualization::Draw( + {visualization::DrawObject("camera", p_camera, true), + visualization::DrawObject{ + "mesh", std::shared_ptr(&sphere), true}}, + "ProjectImagesToAlbedo", 1024, 768); +} + } // namespace tests } // namespace open3d From 9ac925c21fc7680926f2f752981be6fda8299f83 Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Mon, 8 Apr 2024 23:42:00 -0700 Subject: [PATCH 2/8] Bug fixes --- cpp/open3d/t/geometry/TriangleMesh.cpp | 28 ++++++++++++++------------ 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index b7bc2e641b4..48201716ead 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -1331,24 +1331,26 @@ core::Tensor Project( const core::Tensor &t_intrinsic_matrix, // contiguous {3,3} const core::Tensor &t_extrinsic_matrix) { // contiguous {4,4} auto xy_shape = t_xyz.GetShape(); - xy_shape[xy_shape.GetLength() - 1] = 2; + xy_shape[t_xyz.NumDims() - 1] = 2; core::Tensor t_xy(xy_shape, t_xyz.GetDtype()); DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(t_xyz.GetDtype(), [&]() { + // Eigen is column major Eigen::Map> xy(t_xy.GetDataPtr(), 2, t_xy.NumElements() / 2); Eigen::Map> xyz( t_xyz.GetDataPtr(), 3, t_xyz.NumElements() / 3); - Eigen::Map> intrinsic_matrix( + Eigen::Map> KT( t_intrinsic_matrix.GetDataPtr()); - Eigen::Map> extrinsic_matrix( + Eigen::Map> TT( t_extrinsic_matrix.GetDataPtr()); - Eigen::ArrayX pxyz = - intrinsic_matrix * - (extrinsic_matrix.topLeftCorner<3, 3>() * xyz + - extrinsic_matrix.topRightCorner<3, 1>()); // {...,3} - pxyz.leftCols<2>().colwise() /= pxyz.rightCols<1>(); - xy = pxyz.leftCols<2>(); + auto K = KT.transpose(); + auto T = TT.transpose(); + auto R = T.topLeftCorner<3, 3>(); + auto t = T.topRightCorner<3, 1>(); + auto pxyz = (K * ((R * xyz).colwise()+t)).array(); + pxyz.topRows<2>().rowwise() /= pxyz.bottomRows<1>(); + xy = pxyz.topRows<2>(); }); return t_xy; // contiguous {...,2} } @@ -1476,8 +1478,8 @@ Image TriangleMesh::ProjectImagesToAlbedo( width * height); Eigen::Map rays_e(rays.GetDataPtr(), 6, width * height); - Eigen::Map t_hit(result["t_hit"].GetDataPtr(), 3, - width * height); + Eigen::Map t_hit(result["t_hit"].GetDataPtr(), + 1, width * height); // Weight map based on dot product between ray normals and triangle // normals. auto depth = t_hit * (rays_e.bottomRows<3>() - rays_e.topRows<3>()) @@ -1497,7 +1499,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( Eigen::Map weighted_image_e( weighted_image.GetDataPtr(), 4, width * height); // footprint is inf if depth or t_hit is inf. Never zero. - weighted_image_e.bottomRows<3>() = 1. / footprint; + weighted_image_e.bottomRows<1>() = 1. / footprint; /* weighted_image.Slice(2, 3, 4, 1) = 1. / footprint; */ // B. Get texture space (u,v) -> (x,y) map and valid domain in uv space. @@ -1531,7 +1533,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( } albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); - Image albedo_texture(albedo.Slice(2, 0, 3, 1)); + Image albedo_texture(albedo.Slice(2, 0, 3, 1).Contiguous()); if (update_material) { if (!HasMaterial()) { SetMaterial(visualization::rendering::Material()); From a7125a7f0db61b2ee019b6493d64fa9cc193e86b Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Wed, 17 Apr 2024 17:15:15 -0700 Subject: [PATCH 3/8] Debug wip Raycasting fails at vertices. (t_hit is inf) Switch from vertices to triangle centers - Raycasting still fails (t_hit is inf) --- cpp/open3d/t/geometry/TriangleMesh.cpp | 151 +++++++++++++--------- cpp/open3d/t/geometry/kernel/IPPImage.cpp | 96 ++++++++------ cpp/open3d/t/geometry/kernel/IPPImage.h | 7 +- 3 files changed, 151 insertions(+), 103 deletions(-) diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 48201716ead..bcdaca926b8 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -1327,12 +1327,15 @@ TriangleMesh TriangleMesh::RemoveUnreferencedVertices() { namespace { core::Tensor Project( - const core::Tensor &t_xyz, // contiguous {...,3} - const core::Tensor &t_intrinsic_matrix, // contiguous {3,3} - const core::Tensor &t_extrinsic_matrix) { // contiguous {4,4} + const core::Tensor &t_xyz, // contiguous {...,3} + const core::Tensor &t_intrinsic_matrix, // contiguous {3,3} + const core::Tensor &t_extrinsic_matrix, // contiguous {4,4} + core::Tensor &t_xy) { // contiguous {...,2} auto xy_shape = t_xyz.GetShape(); xy_shape[t_xyz.NumDims() - 1] = 2; - core::Tensor t_xy(xy_shape, t_xyz.GetDtype()); + if (t_xy.GetDtype() != t_xyz.GetDtype() || t_xy.GetShape() != xy_shape) { + t_xy = core::Tensor(xy_shape, t_xyz.GetDtype()); + } DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(t_xyz.GetDtype(), [&]() { // Eigen is column major Eigen::Map> xy(t_xy.GetDataPtr(), 2, @@ -1348,9 +1351,8 @@ core::Tensor Project( auto T = TT.transpose(); auto R = T.topLeftCorner<3, 3>(); auto t = T.topRightCorner<3, 1>(); - auto pxyz = (K * ((R * xyz).colwise()+t)).array(); - pxyz.topRows<2>().rowwise() /= pxyz.bottomRows<1>(); - xy = pxyz.topRows<2>(); + auto pxyz = (K * ((R * xyz).colwise() + t)).array(); + xy = pxyz.topRows<2>().rowwise() / pxyz.bottomRows<1>(); }); return t_xy; // contiguous {...,2} } @@ -1366,19 +1368,29 @@ std::pair TriangleMesh::SelectVisible( using core::None; constexpr double EPS = 1e-6; const core::Tensor vertices = GetVertexPositions().To(core::Device()); - core::Tensor imxy = Project(vertices, intrinsic_matrix, extrinsic_matrix); - core::Tensor vis_vert = - core::Tensor::Empty({vertices.GetLength()}, core::Bool); - bool *vis_vert_ptr = vis_vert.GetDataPtr(); + const core::Tensor tris = GetTriangleIndices().To(core::Device()), + trisT = tris.T(); + core::Tensor tri_centers = + (vertices.IndexGet({trisT[0]}) + vertices.IndexGet({trisT[1]}) + + vertices.IndexGet({trisT[2]})) / + 3.f + + 0.01; + core::Tensor imxy; + Project(tri_centers, intrinsic_matrix, extrinsic_matrix, imxy); + core::Tensor vis_tris = + core::Tensor::Empty({tri_centers.GetLength()}, core::Bool); + bool *vis_vert_ptr = vis_tris.GetDataPtr(); int n_vis = 0; - DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(vertices.GetDtype(), [&]() { + DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(tri_centers.GetDtype(), [&]() { scalar_t *imxy_ptr = imxy.GetDataPtr(); - for (int i = 0; i < vertices.GetLength(); ++i, imxy_ptr += 3) { + for (int i = 0; i < tri_centers.GetLength(); ++i, imxy_ptr += 2) { vis_vert_ptr[i] = (imxy_ptr[0] > 0 && imxy_ptr[0] < width_px && imxy_ptr[1] > 0 && imxy_ptr[1] < height_px); n_vis += vis_vert_ptr[i]; } }); + utility::LogInfo("tri_centers: {}\nprojections: {}", tri_centers.ToString(), + imxy.ToString()); core::Tensor cam_loc = -extrinsic_matrix .GetItem({tk::Slice(0, 3, None), tk::Slice(0, 3, None)}) @@ -1386,34 +1398,25 @@ std::pair TriangleMesh::SelectVisible( .Matmul(extrinsic_matrix.GetItem( {tk::Slice(0, 3, None), tk::Index(3)})); utility::LogInfo("cam_loc: {}", cam_loc.ToString()); - core::Tensor rays = core::Tensor::Empty({n_vis, 6}, vertices.GetDtype()); - rays.SetItem({tk::Slice(None, None, None), tk::Slice(0, 3, 1)}, - cam_loc.T()); - core::Tensor vis_vert_idx(vis_vert.NonZero().Reshape({-1})); - rays.SetItem({tk::Slice(None, None, None), tk::Slice(3, None, 1)}, - vertices.GetItem({tk::IndexTensor(vis_vert_idx)})); + core::Tensor rays = core::Tensor::Empty({n_vis, 6}, tri_centers.GetDtype()); + rays.Slice(1, 0, 3, 1) = cam_loc.T(); + core::Tensor vis_vert_idx(vis_tris.NonZero().View({-1})); + rays.Slice(1, 3, 6, 1) = + tri_centers.GetItem({tk::IndexTensor(vis_vert_idx)}); auto result = rcs.CastRays(rays); float *t_hit_ptr = result["t_hit"].GetDataPtr(); int64_t *vis_vert_idx_ptr = vis_vert_idx.GetDataPtr(); + // t_hit is in units of ray segments, so <1 means ray hit something else. for (int i = 0; i < n_vis; ++i) { - vis_vert_ptr[vis_vert_idx_ptr[i]] = t_hit_ptr[i] > 1 - EPS; - } - core::Tensor vis_tris = core::Tensor::Zeros( - {GetTriangleIndices().GetShape(0)}, core::Dtype::Bool); - DISPATCH_INT_DTYPE_PREFIX_TO_TEMPLATE( - GetTriangleIndices().GetDtype(), tri, [&]() { - auto tris_ptr = GetTriangleIndices() - .To(core::Device()) - .Contiguous() - .GetDataPtr(); - auto vis_tris_ptr = vis_tris.GetDataPtr(); - for (int i = 0; i < GetTriangleIndices().GetShape(0); ++i) { - vis_tris_ptr[i] = (vis_vert_ptr[tris_ptr[3 * i + 0]] && - vis_vert_ptr[tris_ptr[3 * i + 1]] && - vis_vert_ptr[tris_ptr[3 * i + 2]]); - } - }); + vis_vert_ptr[vis_vert_idx_ptr[i]] = (t_hit_ptr[i] > 1 - EPS); + } + utility::LogInfo("t_hit = {}", result["t_hit"].ToString()); + utility::LogInfo("vis_tris = {}", vis_tris.ToString()); + core::Tensor vis_vert = core::Tensor::Zeros( + {GetVertexPositions().GetShape(0)}, core::Dtype::Bool); + vis_vert.IndexGet({tris.IndexGet({vis_tris}).View({-1})}) = true; + utility::LogInfo("vis_vert = {}", vis_vert.ToString()); return std::make_pair(vis_vert, vis_tris); } @@ -1424,7 +1427,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( int tex_size /*=1024*/, bool update_material /*=true*/) { using core::None; - using tk = core::TensorKey; + /* using tk = core::TensorKey; */ constexpr double EPS = 1e-6; if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); @@ -1445,12 +1448,14 @@ Image TriangleMesh::ProjectImagesToAlbedo( // (u,v) -> (x,y,z) : {tex_size, tex_size, 3} core::Tensor position_map = BakeVertexAttrTextures( tex_size, {"positions"}, 1, 0, false)["positions"]; + io::WriteImage("position_map.png", + Image(((position_map + 1) * 127.5).To(core::UInt8))); core::Tensor albedo = core::Tensor::Zeros({tex_size, tex_size, 4}, core::Float32); albedo.Slice(2, 3, 4, 1).Fill(EPS); // regularize core::Tensor this_albedo = - core::Tensor::Zeros({tex_size, tex_size, 4}, core::Float32); - core::Tensor weighted_image; + core::Tensor::Empty({tex_size, tex_size, 4}, core::Float32); + core::Tensor weighted_image({}, core::Float32), uv2xy({}, core::Float32); /* For each image, project vertices onto 2D image and do occlusion * culling. This creates (xim, yim) -> (u,v) mapping. Invert this @@ -1467,7 +1472,8 @@ Image TriangleMesh::ProjectImagesToAlbedo( core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); - // A. Get image space weight matrix + // A. Get image space weight matrix, as inverse of pixel footprint on + // the mesh. auto rays = RaycastingScene::CreateRaysPinhole( intrinsic_matrices[i], extrinsic_matrices[i], images[i].GetCols(), images[i].GetRows()); @@ -1480,18 +1486,17 @@ Image TriangleMesh::ProjectImagesToAlbedo( width * height); Eigen::Map t_hit(result["t_hit"].GetDataPtr(), 1, width * height); - // Weight map based on dot product between ray normals and triangle - // normals. auto depth = t_hit * (rays_e.bottomRows<3>() - rays_e.topRows<3>()) .colwise() .norm() .array(); - auto footprint = depth * normals_e - .cwiseProduct((rays_e.bottomRows<3>() - - rays_e.topRows<3>()) - .colwise() - .normalized()) - .sum(); + auto unit_area_proj = (normals_e.cwiseProduct((rays_e.bottomRows<3>() - + rays_e.topRows<3>()) + .colwise() + .normalized())) + .sum(); + auto footprint = + (depth * unit_area_proj).abs(); // ignore face orientation if (!weighted_image.GetShape().IsCompatible({height, width, 4})) { weighted_image = core::Tensor({height, width, 4}, core::Float32); } @@ -1500,11 +1505,10 @@ Image TriangleMesh::ProjectImagesToAlbedo( weighted_image.GetDataPtr(), 4, width * height); // footprint is inf if depth or t_hit is inf. Never zero. weighted_image_e.bottomRows<1>() = 1. / footprint; - /* weighted_image.Slice(2, 3, 4, 1) = 1. / footprint; */ // B. Get texture space (u,v) -> (x,y) map and valid domain in uv space. - core::Tensor uv2xy = Project(position_map, intrinsic_matrices[i], - extrinsic_matrices[i]); // {ts, ts, 2} + Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], + uv2xy); // {ts, ts, 2} core::Tensor visible_triangles = SelectVisible(intrinsic_matrices[i], extrinsic_matrices[i], images[i].GetCols(), images[i].GetRows(), rcs) @@ -1515,17 +1519,48 @@ Image TriangleMesh::ProjectImagesToAlbedo( false)["visible"]; // uv space {ts, ts, 1} /* RemoveTriangleAttr("visible"); */ // do not remap occluded / out of view faces - io::WriteImage("visible_map.png", Image(visible_map.To(core::UInt8))); - uv2xy.Reshape({-1, 2}).GetItem( - tk::IndexTensor(visible_map.Reshape({-1}).NonZero())) = -1.f; + io::WriteImage("visible_map.png", + Image(visible_map.To(core::UInt8) * 255)); + bool *p_visible_map = visible_map.GetDataPtr(); + for (float *p_uv2xy = uv2xy.GetDataPtr(); + p_uv2xy < uv2xy.GetDataPtr() + uv2xy.NumElements(); + p_uv2xy += 2, ++p_visible_map) { + if (!*p_visible_map) *p_uv2xy = *(p_uv2xy + 1) = -1.f; + } + uv2xy = uv2xy.Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} + io::WriteImage("uv2x.png", Image((uv2xy[0].To(core::UInt8)))); + io::WriteImage("uv2y.png", Image((uv2xy[1].To(core::UInt8)))); // albedo[u,v] = image[ i[u,v], j[u,v] ] // ij[u,v] ? = identity[ uv[i,j] ] // C. Interpolate weighted image to weighted texture - ipp::Remap(weighted_image /*{height, width, 4} f32*/, - uv2xy /* {texsz, texsz, 2} f32*/, - this_albedo /*{texsz, texsz, 4} f32*/, + this_albedo.Fill(0.f); + ipp::Remap(weighted_image, /*{height, width, 4} f32*/ + uv2xy[0], /* {texsz, texsz} f32*/ + uv2xy[1], /* {texsz, texsz} f32*/ + this_albedo, /*{texsz, texsz, 4} f32*/ t::geometry::Image::InterpType::Cubic); + io::WriteImage( + "this_albedo.png", + Image((this_albedo.Slice(2, 0, 3, 1) * 255).To(core::UInt8))); + float wtmin = this_albedo.Slice(2, 3, 4, 1) + .Min({0, 1, 2}) + .Item(), + wtmax = this_albedo.Slice(2, 3, 4, 1) + .Max({0, 1, 2}) + .Item(); + io::WriteImage( + "image_weights.png", + Image(weighted_image.Slice(2, 3, 4, 1).Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); + io::WriteImage( + "this_albedo_weight.png", + Image(this_albedo.Slice(2, 3, 4, 1).Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); // weighted sum albedo.Slice(2, 0, 3, 1) += this_albedo.Slice(2, 0, 3, 1) * this_albedo.Slice(2, 3, 4, 1); @@ -1533,7 +1568,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( } albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); - Image albedo_texture(albedo.Slice(2, 0, 3, 1).Contiguous()); + Image albedo_texture((albedo.Slice(2, 0, 3, 1) * 255).To(core::UInt8)); if (update_material) { if (!HasMaterial()) { SetMaterial(visualization::rendering::Material()); diff --git a/cpp/open3d/t/geometry/kernel/IPPImage.cpp b/cpp/open3d/t/geometry/kernel/IPPImage.cpp index e0ba8167401..b4d9350b616 100644 --- a/cpp/open3d/t/geometry/kernel/IPPImage.cpp +++ b/cpp/open3d/t/geometry/kernel/IPPImage.cpp @@ -16,6 +16,7 @@ #include #include +#include "open3d/core/Dispatch.h" #include "open3d/core/Dtype.h" #include "open3d/core/ParallelFor.h" #include "open3d/core/ShapeUtil.h" @@ -82,7 +83,7 @@ void RGBToGray(const core::Tensor &src_im, core::Tensor &dst_im) { void Resize(const core::Tensor &src_im, core::Tensor &dst_im, - t::geometry::Image::InterpType interp_type) { + Image::InterpType interp_type) { auto dtype = src_im.GetDtype(); // Create IPP wrappers for all Open3D tensors const ::ipp::IwiImage ipp_src_im( @@ -96,14 +97,13 @@ void Resize(const core::Tensor &src_im, 0 /* border buffer size */, dst_im.GetDataPtr(), dst_im.GetStride(0) * dtype.ByteSize()); - static const std::unordered_map + static const std::unordered_map type_dict = { - {t::geometry::Image::InterpType::Nearest, ippNearest}, - {t::geometry::Image::InterpType::Linear, ippLinear}, - {t::geometry::Image::InterpType::Cubic, ippCubic}, - {t::geometry::Image::InterpType::Lanczos, ippLanczos}, - {t::geometry::Image::InterpType::Super, ippSuper}, + {Image::InterpType::Nearest, ippNearest}, + {Image::InterpType::Linear, ippLinear}, + {Image::InterpType::Cubic, ippCubic}, + {Image::InterpType::Lanczos, ippLanczos}, + {Image::InterpType::Super, ippSuper}, }; auto it = type_dict.find(interp_type); @@ -296,31 +296,34 @@ void FilterSobel(const core::Tensor &src_im, } } -void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, - const core::Tensor &dst2src_map /*{Wd, Hd, 2}, float*/, - core::Tensor &dst_im /*{Wd, Hd, 2}*/, - t::geometry::Image::InterpType interp_type) { +void Remap(const core::Tensor &src_im, /*{Ws, Hs, C}*/ + const core::Tensor &dst2src_xmap, /*{Wd, Hd}, float*/ + const core::Tensor &dst2src_ymap, /*{Wd, Hd}, float*/ + core::Tensor &dst_im, /*{Wd, Hd, C}*/ + Image::InterpType interp_type) { auto dtype = src_im.GetDtype(); if (dtype != dst_im.GetDtype()) { utility::LogError( "Source ({}) and destination ({}) image dtypes are different!", dtype.ToString(), dst_im.GetDtype().ToString()); } - if (dst2src_map.GetDtype() != core::Float32) { - utility::LogError("dst2src_map dtype ({}) must be Float32.", - dst2src_map.GetDtype().ToString()); + if (dst2src_xmap.GetDtype() != core::Float32) { + utility::LogError("dst2src_xmap dtype ({}) must be Float32.", + dst2src_xmap.GetDtype().ToString()); + } + if (dst2src_ymap.GetDtype() != core::Float32) { + utility::LogError("dst2src_ymap dtype ({}) must be Float32.", + dst2src_ymap.GetDtype().ToString()); } - static const std::unordered_map - interp_dict = { - {t::geometry::Image::InterpType::Nearest, IPPI_INTER_NN}, - {t::geometry::Image::InterpType::Linear, IPPI_INTER_LINEAR}, - {t::geometry::Image::InterpType::Cubic, IPPI_INTER_CUBIC}, - {t::geometry::Image::InterpType::Lanczos, - IPPI_INTER_LANCZOS}, - /* {t::geometry::Image::InterpType::Cubic2p_CatmullRom, */ - /* IPPI_INTER_CUBIC2P_CATMULLROM}, */ - }; + static const std::unordered_map interp_dict = { + {Image::InterpType::Nearest, IPPI_INTER_NN}, + {Image::InterpType::Linear, IPPI_INTER_LINEAR}, + {Image::InterpType::Cubic, IPPI_INTER_CUBIC}, + {Image::InterpType::Lanczos, IPPI_INTER_LANCZOS}, + /* {Image::InterpType::Cubic2p_CatmullRom, */ + /* IPPI_INTER_CUBIC2P_CATMULLROM}, */ + }; auto interp_it = interp_dict.find(interp_type); if (interp_it == interp_dict.end()) { @@ -328,6 +331,20 @@ void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, static_cast(interp_type)); } + IppiSize src_size{static_cast(src_im.GetShape(1)), + static_cast(src_im.GetShape(0))}, + dst_roi_size{static_cast(dst_im.GetShape(1)), + static_cast(dst_im.GetShape(0))}; + IppiRect src_roi{0, 0, static_cast(src_im.GetShape(1)), + static_cast(src_im.GetShape(0))}; + IppStatus sts = ippStsNoErr; + + int src_step = src_im.GetDtype().ByteSize() * src_im.GetStride(0); + int dst_step = dst_im.GetDtype().ByteSize() * dst_im.GetStride(0); + int xmap_step = + dst2src_xmap.GetDtype().ByteSize() * dst2src_xmap.GetStride(0); + int ymap_step = + dst2src_ymap.GetDtype().ByteSize() * dst2src_ymap.GetStride(0); if (src_im.GetDtype() == core::Float32 && src_im.GetShape(2) == 4) { /* IPPAPI(IppStatus, ippiRemap_32f_C4R, (const Ipp32f* pSrc, IppiSize * srcSize, */ @@ -335,28 +352,23 @@ void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, */ /* const Ipp32f* pyMap, int yMapStep, Ipp32f* pDst, int dstStep, */ /* IppiSize dstRoiSize, int interpolation)) */ - IppiSize src_size{static_cast(src_im.GetShape(1)), - static_cast(src_im.GetShape(0))}, - dst_roi_size{static_cast(dst_im.GetShape(1)), - static_cast(dst_im.GetShape(0))}; - IppiRect src_roi{0, 0, static_cast(src_im.GetShape(1)), - static_cast(src_im.GetShape(0))}; - - IppStatus sts = ippiRemap_32f_C4R( - src_im.GetDataPtr(), src_size, 1, src_roi, - dst2src_map.GetDataPtr(), 2, - dst2src_map.GetDataPtr() + 1, 2, - dst_im.GetDataPtr(), 1, dst_roi_size, interp_it->second); - if (sts != ippStsNoErr) { - // See comments in icv/include/ippicv_types.h for meaning - utility::LogError("IPP remap error {}", sts); - } - + const auto p_src_im = src_im.GetDataPtr(); + auto p_dst_im = dst_im.GetDataPtr(); + const auto p_dst2src_xmap = dst2src_xmap.GetDataPtr(); + const auto p_dst2src_ymap = dst2src_ymap.GetDataPtr(); + sts = ippiRemap_32f_C4R(p_src_im, src_size, src_step, src_roi, + p_dst2src_xmap, xmap_step, p_dst2src_ymap, + ymap_step, p_dst_im, dst_step, dst_roi_size, + interp_it->second); } else { utility::LogError( "Remap not implemented for dtype ({}) and channels ({}).", src_im.GetDtype().ToString(), src_im.GetShape(2)); } + if (sts != ippStsNoErr) { + // See comments in icv/include/ippicv_types.h for meaning + utility::LogError("IPP remap error {}", sts); + } } } // namespace ipp diff --git a/cpp/open3d/t/geometry/kernel/IPPImage.h b/cpp/open3d/t/geometry/kernel/IPPImage.h index b5c7fddab6c..774a1953416 100644 --- a/cpp/open3d/t/geometry/kernel/IPPImage.h +++ b/cpp/open3d/t/geometry/kernel/IPPImage.h @@ -76,9 +76,10 @@ void FilterSobel(const core::Tensor &srcim, core::Tensor &dstim_dy, int kernel_size); -void Remap(const core::Tensor &src_im /*{Ws, Hs, C}*/, - const core::Tensor &dst2src_map /*{Wd, Hd, 2}, float*/, - core::Tensor &dst_im /*{Wd, Hd, 2}*/, +void Remap(const core::Tensor &src_im, /*{Ws, Hs, C}*/ + const core::Tensor &dst2src_xmap, /*{Wd, Hd}, float*/ + const core::Tensor &dst2src_ymap, /*{Wd, Hd, 2}, float*/ + core::Tensor &dst_im, /*{Wd, Hd, 2}*/ Image::InterpType interp_type); } // namespace ipp From 5a8f7d2d2aafed08f21a673d04f8d081ec2bbbc4 Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Thu, 2 May 2024 23:25:47 -0700 Subject: [PATCH 4/8] WIP changes Remove SelectVisible (implementation was incorrect). Instead use visibility map for each texel directly. Add pixel_foreshortening Python test example by rendering images and using them to project albedo. --- cpp/open3d/io/ImageIO.cpp | 3 +- cpp/open3d/t/geometry/TriangleMesh.cpp | 175 +++++++----------- cpp/open3d/t/geometry/TriangleMesh.h | 7 - cpp/open3d/t/io/ImageIO.cpp | 3 +- cpp/pybind/t/geometry/trianglemesh.cpp | 9 + cpp/pybind/visualization/o3dvisualizer.cpp | 4 +- cpp/tests/t/geometry/TriangleMesh.cpp | 72 +++++-- .../triangle_mesh_project_to_albedo.py | 88 +++++++++ python/open3d/visualization/draw.py | 4 +- 9 files changed, 225 insertions(+), 140 deletions(-) create mode 100644 examples/python/geometry/triangle_mesh_project_to_albedo.py diff --git a/cpp/open3d/io/ImageIO.cpp b/cpp/open3d/io/ImageIO.cpp index a4126df14e6..68e1af09de0 100644 --- a/cpp/open3d/io/ImageIO.cpp +++ b/cpp/open3d/io/ImageIO.cpp @@ -77,7 +77,8 @@ bool WriteImage(const std::string &filename, auto map_itr = file_extension_to_image_write_function.find(filename_ext); if (map_itr == file_extension_to_image_write_function.end()) { utility::LogWarning( - "Write geometry::Image failed: unknown file extension."); + "Write geometry::Image failed: file extension {} unknown.", + filename_ext); return false; } return map_itr->second(filename, image, quality); diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 944cca59f18..55b848b69b6 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -17,6 +17,7 @@ #include #include +#include #include #include #include @@ -1180,7 +1181,7 @@ static bool IsNegative(T val) { } TriangleMesh TriangleMesh::SelectByIndex(const core::Tensor &indices, - bool copy_attributes) const { + bool copy_attributes /*=true*/) const { core::AssertTensorShape(indices, {indices.GetLength()}); if (indices.NumElements() == 0) { return {}; @@ -1343,15 +1344,15 @@ TriangleMesh TriangleMesh::RemoveUnreferencedVertices() { namespace { -core::Tensor Project( - const core::Tensor &t_xyz, // contiguous {...,3} - const core::Tensor &t_intrinsic_matrix, // contiguous {3,3} - const core::Tensor &t_extrinsic_matrix, // contiguous {4,4} - core::Tensor &t_xy) { // contiguous {...,2} +core::Tensor Project(const core::Tensor &t_xyz, // contiguous {...,3} + const core::Tensor &t_intrinsic_matrix, // {3,3} + const core::Tensor &t_extrinsic_matrix, // {4,4} + core::Tensor &t_xy) { // contiguous {...,2} auto xy_shape = t_xyz.GetShape(); + auto dt = t_xyz.GetDtype(); xy_shape[t_xyz.NumDims() - 1] = 2; - if (t_xy.GetDtype() != t_xyz.GetDtype() || t_xy.GetShape() != xy_shape) { - t_xy = core::Tensor(xy_shape, t_xyz.GetDtype()); + if (t_xy.GetDtype() != dt || t_xy.GetShape() != xy_shape) { + t_xy = core::Tensor(xy_shape, dt); } DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(t_xyz.GetDtype(), [&]() { // Eigen is column major @@ -1360,9 +1361,9 @@ core::Tensor Project( Eigen::Map> xyz( t_xyz.GetDataPtr(), 3, t_xyz.NumElements() / 3); Eigen::Map> KT( - t_intrinsic_matrix.GetDataPtr()); + t_intrinsic_matrix.To(dt).Contiguous().GetDataPtr()); Eigen::Map> TT( - t_extrinsic_matrix.GetDataPtr()); + t_extrinsic_matrix.To(dt).Contiguous().GetDataPtr()); auto K = KT.transpose(); auto T = TT.transpose(); @@ -1375,68 +1376,6 @@ core::Tensor Project( } } // namespace -std::pair TriangleMesh::SelectVisible( - const core::Tensor &intrinsic_matrix, - const core::Tensor &extrinsic_matrix, - int width_px, - int height_px, - const RaycastingScene &rcs) { - using tk = core::TensorKey; - using core::None; - constexpr double EPS = 1e-6; - const core::Tensor vertices = GetVertexPositions().To(core::Device()); - const core::Tensor tris = GetTriangleIndices().To(core::Device()), - trisT = tris.T(); - core::Tensor tri_centers = - (vertices.IndexGet({trisT[0]}) + vertices.IndexGet({trisT[1]}) + - vertices.IndexGet({trisT[2]})) / - 3.f + - 0.01; - core::Tensor imxy; - Project(tri_centers, intrinsic_matrix, extrinsic_matrix, imxy); - core::Tensor vis_tris = - core::Tensor::Empty({tri_centers.GetLength()}, core::Bool); - bool *vis_vert_ptr = vis_tris.GetDataPtr(); - int n_vis = 0; - DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(tri_centers.GetDtype(), [&]() { - scalar_t *imxy_ptr = imxy.GetDataPtr(); - for (int i = 0; i < tri_centers.GetLength(); ++i, imxy_ptr += 2) { - vis_vert_ptr[i] = (imxy_ptr[0] > 0 && imxy_ptr[0] < width_px && - imxy_ptr[1] > 0 && imxy_ptr[1] < height_px); - n_vis += vis_vert_ptr[i]; - } - }); - utility::LogInfo("tri_centers: {}\nprojections: {}", tri_centers.ToString(), - imxy.ToString()); - core::Tensor cam_loc = - -extrinsic_matrix - .GetItem({tk::Slice(0, 3, None), tk::Slice(0, 3, None)}) - .T() - .Matmul(extrinsic_matrix.GetItem( - {tk::Slice(0, 3, None), tk::Index(3)})); - utility::LogInfo("cam_loc: {}", cam_loc.ToString()); - core::Tensor rays = core::Tensor::Empty({n_vis, 6}, tri_centers.GetDtype()); - rays.Slice(1, 0, 3, 1) = cam_loc.T(); - core::Tensor vis_vert_idx(vis_tris.NonZero().View({-1})); - rays.Slice(1, 3, 6, 1) = - tri_centers.GetItem({tk::IndexTensor(vis_vert_idx)}); - - auto result = rcs.CastRays(rays); - float *t_hit_ptr = result["t_hit"].GetDataPtr(); - int64_t *vis_vert_idx_ptr = vis_vert_idx.GetDataPtr(); - // t_hit is in units of ray segments, so <1 means ray hit something else. - for (int i = 0; i < n_vis; ++i) { - vis_vert_ptr[vis_vert_idx_ptr[i]] = (t_hit_ptr[i] > 1 - EPS); - } - utility::LogInfo("t_hit = {}", result["t_hit"].ToString()); - utility::LogInfo("vis_tris = {}", vis_tris.ToString()); - core::Tensor vis_vert = core::Tensor::Zeros( - {GetVertexPositions().GetShape(0)}, core::Dtype::Bool); - vis_vert.IndexGet({tris.IndexGet({vis_tris}).View({-1})}) = true; - utility::LogInfo("vis_vert = {}", vis_vert.ToString()); - return std::make_pair(vis_vert, vis_tris); -} - Image TriangleMesh::ProjectImagesToAlbedo( const std::vector &images, const std::vector &intrinsic_matrices, @@ -1444,8 +1383,9 @@ Image TriangleMesh::ProjectImagesToAlbedo( int tex_size /*=1024*/, bool update_material /*=true*/) { using core::None; - /* using tk = core::TensorKey; */ - constexpr double EPS = 1e-6; + using tk = core::TensorKey; + float pixel_foreshortening_threshold = 0.0f; + constexpr float EPS = 1e-6; if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); } @@ -1473,6 +1413,8 @@ Image TriangleMesh::ProjectImagesToAlbedo( core::Tensor this_albedo = core::Tensor::Empty({tex_size, tex_size, 4}, core::Float32); core::Tensor weighted_image({}, core::Float32), uv2xy({}, core::Float32); + core::Tensor uvrays = + core::Tensor::Empty({tex_size, tex_size, 6}, core::Float32); /* For each image, project vertices onto 2D image and do occlusion * culling. This creates (xim, yim) -> (u,v) mapping. Invert this @@ -1485,6 +1427,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( rcs.AddTriangles(*this); for (size_t i = 0; i < images.size(); ++i) { + auto i_str = std::to_string(i); auto width = images[i].GetCols(), height = images[i].GetRows(); core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); @@ -1492,61 +1435,69 @@ Image TriangleMesh::ProjectImagesToAlbedo( // A. Get image space weight matrix, as inverse of pixel footprint on // the mesh. auto rays = RaycastingScene::CreateRaysPinhole( - intrinsic_matrices[i], extrinsic_matrices[i], - images[i].GetCols(), images[i].GetRows()); + intrinsic_matrices[i], extrinsic_matrices[i], width, height); + core::Tensor cam_loc = + rays.GetItem({tk::Index(0), tk::Index(0), tk::Slice(0, 3, 1)}); auto result = rcs.CastRays(rays); // Eigen is column-major order - Eigen::Map normals_e( + Eigen::Map normals_e( result["primitive_normals"].GetDataPtr(), 3, width * height); - Eigen::Map rays_e(rays.GetDataPtr(), 6, + Eigen::Map rays_e(rays.GetDataPtr(), 6, width * height); Eigen::Map t_hit(result["t_hit"].GetDataPtr(), 1, width * height); - auto depth = t_hit * (rays_e.bottomRows<3>() - rays_e.topRows<3>()) - .colwise() - .norm() - .array(); - auto unit_area_proj = (normals_e.cwiseProduct((rays_e.bottomRows<3>() - - rays_e.topRows<3>()) - .colwise() - .normalized())) - .sum(); + auto depth = t_hit * rays_e.bottomRows<3>().colwise().norm().array(); + // removing this eval() increase runtime a lot (?) + auto rays_dir = rays_e.bottomRows<3>().colwise().normalized().eval(); + auto pixel_foreshortening = (normals_e * rays_dir) + .colwise() + .sum() + .abs(); // ignore face orientation auto footprint = - (depth * unit_area_proj).abs(); // ignore face orientation + (depth * (pixel_foreshortening < pixel_foreshortening_threshold) + .select(INFINITY, pixel_foreshortening)) + .eval(); + // fix for bad normals + footprint = footprint.isNaN().select(INFINITY, footprint).eval(); + + utility::LogInfo("Image {}, footprint range: {}-{}", i, + footprint.minCoeff(), footprint.maxCoeff()); if (!weighted_image.GetShape().IsCompatible({height, width, 4})) { weighted_image = core::Tensor({height, width, 4}, core::Float32); } - weighted_image.Slice(2, 0, 3, 1) = images[i].AsTensor(); + weighted_image.Slice(2, 0, 3, 1) = + images[i].To(core::Float32).AsTensor(); // range: [0,1] Eigen::Map weighted_image_e( weighted_image.GetDataPtr(), 4, width * height); // footprint is inf if depth or t_hit is inf. Never zero. weighted_image_e.bottomRows<1>() = 1. / footprint; // B. Get texture space (u,v) -> (x,y) map and valid domain in uv space. + utility::LogInfo("cam_loc={}", cam_loc.ToString()); + uvrays.GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(0, 3, 1)}) = cam_loc; + uvrays.GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(3, 6, 1)}) = position_map - cam_loc; + auto uvresult = rcs.CastRays(uvrays); + auto &t_hit_uv = uvresult["t_hit"]; + io::WriteImage("t_hit_uv_" + i_str + ".png", + Image((t_hit_uv * 255).To(core::UInt8))); + Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], uv2xy); // {ts, ts, 2} - core::Tensor visible_triangles = - SelectVisible(intrinsic_matrices[i], extrinsic_matrices[i], - images[i].GetCols(), images[i].GetRows(), rcs) - .second; - SetTriangleAttr("visible", visible_triangles); - core::Tensor visible_map = BakeTriangleAttrTextures( - tex_size, {"visible"}, 2, 0, - false)["visible"]; // uv space {ts, ts, 1} - /* RemoveTriangleAttr("visible"); */ - // do not remap occluded / out of view faces - io::WriteImage("visible_map.png", - Image(visible_map.To(core::UInt8) * 255)); - bool *p_visible_map = visible_map.GetDataPtr(); - for (float *p_uv2xy = uv2xy.GetDataPtr(); + // Disable self-occluded points + for (float *p_uv2xy = uv2xy.GetDataPtr(), + *p_t_hit = t_hit_uv.GetDataPtr(); p_uv2xy < uv2xy.GetDataPtr() + uv2xy.NumElements(); - p_uv2xy += 2, ++p_visible_map) { - if (!*p_visible_map) *p_uv2xy = *(p_uv2xy + 1) = -1.f; + p_uv2xy += 2, ++p_t_hit) { + if (*p_t_hit < 1 - EPS) *p_uv2xy = *(p_uv2xy + 1) = -1.f; } uv2xy = uv2xy.Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} - io::WriteImage("uv2x.png", Image((uv2xy[0].To(core::UInt8)))); - io::WriteImage("uv2y.png", Image((uv2xy[1].To(core::UInt8)))); + io::WriteImage("uv2x_" + i_str + ".png", + Image((uv2xy[0].To(core::UInt8)))); + io::WriteImage("uv2y_" + i_str + ".png", + Image((uv2xy[1].To(core::UInt8)))); // albedo[u,v] = image[ i[u,v], j[u,v] ] // ij[u,v] ? = identity[ uv[i,j] ] @@ -1556,9 +1507,10 @@ Image TriangleMesh::ProjectImagesToAlbedo( uv2xy[0], /* {texsz, texsz} f32*/ uv2xy[1], /* {texsz, texsz} f32*/ this_albedo, /*{texsz, texsz, 4} f32*/ - t::geometry::Image::InterpType::Cubic); + t::geometry::Image::InterpType::Linear); + // Weights can become negative with higher order interpolation io::WriteImage( - "this_albedo.png", + "this_albedo_" + i_str + ".png", Image((this_albedo.Slice(2, 0, 3, 1) * 255).To(core::UInt8))); float wtmin = this_albedo.Slice(2, 3, 4, 1) .Min({0, 1, 2}) @@ -1566,14 +1518,15 @@ Image TriangleMesh::ProjectImagesToAlbedo( wtmax = this_albedo.Slice(2, 3, 4, 1) .Max({0, 1, 2}) .Item(); + utility::LogInfo("Image {}, weight range: {}-{}", i, wtmin, wtmax); io::WriteImage( - "image_weights.png", + "image_weights_" + i_str + ".png", Image(weighted_image.Slice(2, 3, 4, 1).Contiguous()) .To(core::UInt8, /*copy=*/false, /*scale=*/255.f / (wtmax - wtmin), /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); io::WriteImage( - "this_albedo_weight.png", + "this_albedo_weight_" + i_str + ".png", Image(this_albedo.Slice(2, 3, 4, 1).Contiguous()) .To(core::UInt8, /*copy=*/false, /*scale=*/255.f / (wtmax - wtmin), diff --git a/cpp/open3d/t/geometry/TriangleMesh.h b/cpp/open3d/t/geometry/TriangleMesh.h index 04953574beb..2dbb71ec55d 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.h +++ b/cpp/open3d/t/geometry/TriangleMesh.h @@ -979,13 +979,6 @@ class TriangleMesh : public Geometry, public DrawableGeometry { /// \return The reference to itself. TriangleMesh RemoveUnreferencedVertices(); - std::pair SelectVisible( - const core::Tensor &intrinsic_matrix, - const core::Tensor &extrinsic_matrix, - int width_px, - int height_px, - const RaycastingScene &rcs); - Image ProjectImagesToAlbedo( const std::vector &images, const std::vector &intrinsic_matrices, diff --git a/cpp/open3d/t/io/ImageIO.cpp b/cpp/open3d/t/io/ImageIO.cpp index df1ca3625ce..52d0305b6c3 100644 --- a/cpp/open3d/t/io/ImageIO.cpp +++ b/cpp/open3d/t/io/ImageIO.cpp @@ -79,7 +79,8 @@ bool WriteImage(const std::string &filename, auto map_itr = file_extension_to_image_write_function.find(filename_ext); if (map_itr == file_extension_to_image_write_function.end()) { utility::LogWarning( - "Write geometry::Image failed: unknown file extension."); + "Write geometry::Image failed: file extension {} unknown.", + filename_ext); return false; } return map_itr->second(filename, image.To(core::Device("CPU:0")), quality); diff --git a/cpp/pybind/t/geometry/trianglemesh.cpp b/cpp/pybind/t/geometry/trianglemesh.cpp index 213b4e11e09..fabc79dc38e 100644 --- a/cpp/pybind/t/geometry/trianglemesh.cpp +++ b/cpp/pybind/t/geometry/trianglemesh.cpp @@ -970,6 +970,7 @@ the partition id for each face. triangle_mesh.def( "select_by_index", &TriangleMesh::SelectByIndex, "indices"_a, + "copy_attributes"_a = true, R"(Returns a new mesh with the vertices selected according to the indices list. If an item from the indices list exceeds the max vertex number of the mesh or has a negative value, it is ignored. @@ -977,6 +978,9 @@ or has a negative value, it is ignored. Args: indices (open3d.core.Tensor): An integer list of indices. Duplicates are allowed, but ignored. Signed and unsigned integral types are accepted. + copy_attributes (bool): Indicates if vertex attributes (other than + positions) and triangle attributes (other than indices) should be copied to + the returned mesh. Returns: A new mesh with the selected vertices and faces built from these vertices. @@ -995,6 +999,11 @@ or has a negative value, it is ignored. triangle_mesh.def("remove_unreferenced_vertices", &TriangleMesh::RemoveUnreferencedVertices, "Removes unreferenced vertices from the mesh in-place."); + + triangle_mesh.def("project_images_to_albedo", + &TriangleMesh::ProjectImagesToAlbedo, "images"_a, + "intrinsic_matrices"_a, "extrinsic_matrices"_a, + "tex_size"_a = 1024, "update_material"_a = true); } } // namespace geometry diff --git a/cpp/pybind/visualization/o3dvisualizer.cpp b/cpp/pybind/visualization/o3dvisualizer.cpp index a97d108f918..9f8e22a1217 100644 --- a/cpp/pybind/visualization/o3dvisualizer.cpp +++ b/cpp/pybind/visualization/o3dvisualizer.cpp @@ -315,8 +315,8 @@ void pybind_o3dvisualizer(py::module& m) { "redraw is required.", "callback"_a) .def("export_current_image", &O3DVisualizer::ExportCurrentImage, - "Exports a PNG image of what is currently displayed to the " - "given path.", + "Exports a PNG or JPEG image of what is currently displayed " + "to the given path.", "path"_a) .def("start_rpc_interface", &O3DVisualizer::StartRPCInterface, "address"_a, "timeout"_a, diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index 79ca6af4e55..37d43134b0e 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -1350,31 +1350,71 @@ TEST_P(TriangleMeshPermuteDevices, RemoveUnreferencedVertices) { TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { using namespace t::geometry; + using ls = open3d::geometry::LineSet; core::Device device = GetParam(); TriangleMesh sphere = TriangleMesh::FromLegacy(*geometry::TriangleMesh::CreateSphere( - 1.0, 10, /*create_uv_map=*/true)); - core::Tensor view = core::Tensor::Zeros({192, 256, 3}, core::Float32); - view.Slice(2, 0, 1, 1).Fill(1.0); // red + 1.0, 20, /*create_uv_map=*/true)); + core::Tensor view[3] = {core::Tensor::Zeros({192, 256, 3}, core::Float32), + core::Tensor::Zeros({192, 256, 3}, core::Float32), + core::Tensor::Zeros({192, 256, 3}, core::Float32)}; + view[0].Slice(2, 0, 1, 1).Fill(1.0); // red + view[1].Slice(2, 1, 2, 1).Fill(1.0); // green + view[2].Slice(2, 2, 3, 1).Fill(1.0); // blue core::Tensor intrinsic_matrix = core::Tensor::Init( - {{256, 0, 96}, {0, 256, 128}, {0, 0, 1}}, device); - core::Tensor extrinsic_matrix = core::Tensor::Init( - {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 1}, {0, 0, 0, 1}}, device); + {{256, 0, 128}, {0, 256, 96}, {0, 0, 1}}, device); + core::Tensor extrinsic_matrix[3] = { + core::Tensor::Init( + {{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, 1, 3}, {0, 0, 0, 1}}, + device), + core::Tensor::Init({{-0.5, 0, -0.8660254, 0}, + {0, 1, 0, 0}, + {0.8660254, 0, -0.5, 3}, + {0, 0, 0, 1}}, + device), + core::Tensor::Init({{-0.5, 0, 0.8660254, 0}, + {0, 1, 0, 0}, + {-0.8660254, 0, -0.5, 3}, + {0, 0, 0, 1}}, + device), + }; Eigen::Map e_intrinsic( intrinsic_matrix.GetDataPtr()); - Eigen::Map e_extrinsic( - extrinsic_matrix.GetDataPtr()); - auto p_camera = open3d::geometry::LineSet::CreateCameraVisualization( - 256, 192, e_intrinsic.cast(), e_extrinsic.cast()); - - Image texture = sphere.ProjectImagesToAlbedo( - {Image(view)}, {intrinsic_matrix}, {extrinsic_matrix}, 256, true); - t::io::WriteImage("texture.png", texture); - t::io::WriteTriangleMesh("sphere-projected.obj", sphere); + Eigen::Map e_extrinsic[3] = { + Eigen::Map( + extrinsic_matrix[0].GetDataPtr()), + Eigen::Map( + extrinsic_matrix[1].GetDataPtr()), + Eigen::Map( + extrinsic_matrix[2].GetDataPtr())}; + std::shared_ptr p_camera[3] = { + ls::CreateCameraVisualization( + 256, 192, e_intrinsic.transpose().cast(), + e_extrinsic[0].transpose().cast()), + ls::CreateCameraVisualization( + 256, 192, e_intrinsic.transpose().cast(), + e_extrinsic[1].transpose().cast()), + ls::CreateCameraVisualization( + 256, 192, e_intrinsic.transpose().cast(), + e_extrinsic[2].transpose().cast())}; + + Image albedo = sphere.ProjectImagesToAlbedo( + {Image(view[0]), Image(view[1]), Image(view[2])}, + {intrinsic_matrix, intrinsic_matrix, intrinsic_matrix}, + {extrinsic_matrix[0], extrinsic_matrix[1], extrinsic_matrix[2]}, + 256, true); + utility::LogInfo("Mesh: {}", sphere.ToString()); + utility::LogInfo("Texture: {}", albedo.ToString()); + t::io::WriteImage("albedo.png", albedo); + /* t::io::WriteTriangleMesh("sphere-projected.obj", sphere); */ + t::io::WriteTriangleMesh("sphere-projected.glb", sphere); + t::io::WriteTriangleMesh("sphere-projected.npz", sphere); visualization::Draw( - {visualization::DrawObject("camera", p_camera, true), + {visualization::DrawObject("camera_0", p_camera[0], true), + visualization::DrawObject("camera_1", p_camera[1], true), + visualization::DrawObject("camera_2", p_camera[2], true), visualization::DrawObject{ "mesh", std::shared_ptr(&sphere), true}}, "ProjectImagesToAlbedo", 1024, 768); diff --git a/examples/python/geometry/triangle_mesh_project_to_albedo.py b/examples/python/geometry/triangle_mesh_project_to_albedo.py new file mode 100644 index 00000000000..4981db7769c --- /dev/null +++ b/examples/python/geometry/triangle_mesh_project_to_albedo.py @@ -0,0 +1,88 @@ +import argparse +from pathlib import Path +import numpy as np +import open3d as o3d +from open3d.visualization import gui, rendering, O3DVisualizer +from open3d.core import Tensor + + +def create_dataset(meshfile): + width, height = 1024, 768 + K = np.array([[512, 0, 512], [0, -512, 384], [0, 0, 1]]) + model = o3d.io.read_triangle_model(meshfile) + t = np.array([0, 0, 0.5]) # origin / object in camera ref frame + + def move_camera_and_shoot(o3dvis): + Rts = [] + images = [] + o3dvis.scene.scene.enable_sun_light(False) + for n in range(12): + theta = n * (2 * np.pi) / 12 + Rt = np.eye(4) + Rt[:3, 3] = t + Rt[:3, : + 3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_axis_angle( + [0, theta, 0]) + print(f"{Rt=}") + Rts.append(Rt) + o3dvis.setup_camera(K, Rt, width, height) + o3dvis.post_redraw() + o3dvis.export_current_image(f"render-{n}.jpg") + images.append(f"render-{n}.jpg") + np.savez("cameras.npz", + width=width, + height=height, + K=K, + Rts=Rts, + images=images) + + o3d.visualization.draw([model], + show_ui=False, + width=width, + height=height, + actions=[("Save Images", move_camera_and_shoot)]) + + # Linux only :-( + # render = rendering.OffscreenRenderer(width, height) + # render.scene.add_geometry(model) + # img = render.render_to_image() + # o3d.io.write_image("render-image.jpg", img) + + +def albedo_from_images(meshfile, calib_data_file): + + model = o3d.io.read_triangle_model(meshfile) + tmeshes = o3d.t.geometry.TriangleMesh.from_triangle_mesh_model(model) + tmeshes = list(tmeshes.values()) + calib = np.load(calib_data_file) + Ks = list(Tensor(calib["K"]) for _ in range(len(calib["Rts"])))[:2] + Rts = list(Tensor(Rt) for Rt in calib["Rts"])[:2] + images = list(o3d.t.io.read_image(imfile) for imfile in calib["images"])[:2] + calib.close() + tmeshes[0].project_images_to_albedo(images, Ks, Rts, 256) + cam_vis = list({ + "name": + f"camera-{i}", + "geometry": + o3d.geometry.LineSet.create_camera_visualization( + images[0].columns, images[0].rows, K.numpy(), Rt.numpy(), 0.1) + } for i, (K, Rt) in enumerate(zip(Ks, Rts))) + o3d.visualization.draw(cam_vis + [{ + "name": meshfile.name, + "geometry": tmeshes[0] + }], + show_ui=True) + + +if __name__ == "__main__": + + parser = argparse.ArgumentParser() + parser.add_argument("action", + choices=('create_dataset', 'albedo_from_images')) + parser.add_argument("meshfile", type=Path) + args = parser.parse_args() + + if args.action == "create_dataset": + create_dataset(args.meshfile) + else: + albedo_from_images(args.meshfile, "cameras.npz") diff --git a/python/open3d/visualization/draw.py b/python/open3d/visualization/draw.py index f7764b41400..96396a0ffd0 100644 --- a/python/open3d/visualization/draw.py +++ b/python/open3d/visualization/draw.py @@ -92,12 +92,12 @@ def draw(geometry=None, on_animation_frame (Callable): Callback for each animation frame update with signature:: - Callback(O3DVisualizer, double time) -> None + Callback(O3DVisualizer o3dvis, double time) -> None on_animation_tick (Callable): Callback for each animation time step with signature:: - Callback(O3DVisualizer, double tick_duration, double time) -> TickResult + Callback(O3DVisualizer o3dvis, double tick_duration, double time) -> TickResult If the callback returns ``TickResult.REDRAW``, the scene is redrawn. It should return ``TickResult.NOCHANGE`` if redraw is not required. From 08a9e69970a51611deda7f80c5e192b45b1a0e02 Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Sat, 18 May 2024 08:22:20 -0700 Subject: [PATCH 5/8] Fix memory error --- cpp/open3d/t/geometry/TriangleMesh.cpp | 25 +++--- cpp/tests/t/geometry/TriangleMesh.cpp | 3 +- .../triangle_mesh_project_to_albedo.py | 79 +++++++++++++------ python/open3d/visualization/draw.py | 2 +- 4 files changed, 74 insertions(+), 35 deletions(-) diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 55b848b69b6..dfa4659b749 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -46,6 +46,7 @@ #include "open3d/t/geometry/kernel/TriangleMesh.h" #include "open3d/t/geometry/kernel/UVUnwrapping.h" #include "open3d/t/io/ImageIO.h" +#include "open3d/t/io/NumpyIO.h" #include "open3d/utility/Optional.h" #include "open3d/utility/ParallelScan.h" @@ -1348,22 +1349,26 @@ core::Tensor Project(const core::Tensor &t_xyz, // contiguous {...,3} const core::Tensor &t_intrinsic_matrix, // {3,3} const core::Tensor &t_extrinsic_matrix, // {4,4} core::Tensor &t_xy) { // contiguous {...,2} + utility::LogInfo("K: {}\ncTw: {}", t_intrinsic_matrix.ToString(), + t_extrinsic_matrix.ToString()); auto xy_shape = t_xyz.GetShape(); auto dt = t_xyz.GetDtype(); + auto t_K = t_intrinsic_matrix.To(dt).Contiguous(), + t_T = t_extrinsic_matrix.To(dt).Contiguous(); xy_shape[t_xyz.NumDims() - 1] = 2; if (t_xy.GetDtype() != dt || t_xy.GetShape() != xy_shape) { t_xy = core::Tensor(xy_shape, dt); } - DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(t_xyz.GetDtype(), [&]() { + DISPATCH_FLOAT_DTYPE_TO_TEMPLATE(dt, [&]() { // Eigen is column major Eigen::Map> xy(t_xy.GetDataPtr(), 2, t_xy.NumElements() / 2); Eigen::Map> xyz( t_xyz.GetDataPtr(), 3, t_xyz.NumElements() / 3); Eigen::Map> KT( - t_intrinsic_matrix.To(dt).Contiguous().GetDataPtr()); + t_K.GetDataPtr()); Eigen::Map> TT( - t_extrinsic_matrix.To(dt).Contiguous().GetDataPtr()); + t_T.GetDataPtr()); auto K = KT.transpose(); auto T = TT.transpose(); @@ -1384,7 +1389,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( bool update_material /*=true*/) { using core::None; using tk = core::TensorKey; - float pixel_foreshortening_threshold = 0.0f; + float pixel_foreshortening_threshold = 0.25f; constexpr float EPS = 1e-6; if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); @@ -1454,12 +1459,12 @@ Image TriangleMesh::ProjectImagesToAlbedo( .colwise() .sum() .abs(); // ignore face orientation - auto footprint = + auto footprint_err = (depth * (pixel_foreshortening < pixel_foreshortening_threshold) - .select(INFINITY, pixel_foreshortening)) - .eval(); + .select(INFINITY, pixel_foreshortening)); // fix for bad normals - footprint = footprint.isNaN().select(INFINITY, footprint).eval(); + auto footprint = + footprint_err.isNaN().select(INFINITY, footprint_err).eval(); utility::LogInfo("Image {}, footprint range: {}-{}", i, footprint.minCoeff(), footprint.maxCoeff()); @@ -1495,9 +1500,9 @@ Image TriangleMesh::ProjectImagesToAlbedo( } uv2xy = uv2xy.Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} io::WriteImage("uv2x_" + i_str + ".png", - Image((uv2xy[0].To(core::UInt8)))); + Image((uv2xy[0].To(core::UInt16)))); io::WriteImage("uv2y_" + i_str + ".png", - Image((uv2xy[1].To(core::UInt8)))); + Image((uv2xy[1].To(core::UInt16)))); // albedo[u,v] = image[ i[u,v], j[u,v] ] // ij[u,v] ? = identity[ uv[i,j] ] diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index 37d43134b0e..a431804da9e 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -1416,7 +1416,8 @@ TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { visualization::DrawObject("camera_1", p_camera[1], true), visualization::DrawObject("camera_2", p_camera[2], true), visualization::DrawObject{ - "mesh", std::shared_ptr(&sphere), true}}, + "mesh", std::make_shared(std::move(sphere)), + true}}, "ProjectImagesToAlbedo", 1024, 768); } diff --git a/examples/python/geometry/triangle_mesh_project_to_albedo.py b/examples/python/geometry/triangle_mesh_project_to_albedo.py index 4981db7769c..a8ed06e552c 100644 --- a/examples/python/geometry/triangle_mesh_project_to_albedo.py +++ b/examples/python/geometry/triangle_mesh_project_to_albedo.py @@ -1,46 +1,66 @@ import argparse from pathlib import Path +import subprocess as sp import numpy as np import open3d as o3d from open3d.visualization import gui, rendering, O3DVisualizer from open3d.core import Tensor -def create_dataset(meshfile): - width, height = 1024, 768 - K = np.array([[512, 0, 512], [0, -512, 384], [0, 0, 1]]) +def create_dataset(meshfile, n_images=10, movie=False): + SCALING = 2 + width, height = 1024, 1024 + focal_length = 512 + K = np.array([[focal_length, 0, width / 2], [0, focal_length, height / 2], + [0, 0, 1]]) model = o3d.io.read_triangle_model(meshfile) - t = np.array([0, 0, 0.5]) # origin / object in camera ref frame + unlit = rendering.MaterialRecord() + unlit.shader = "unlit" + t = np.array([0, 0, 0.3]) # origin / object in camera ref frame - def move_camera_and_shoot(o3dvis): + def rotate_camera_and_shoot(o3dvis): Rts = [] images = [] o3dvis.scene.scene.enable_sun_light(False) - for n in range(12): - theta = n * (2 * np.pi) / 12 + print("Rendering images: ", end='', flush=True) + for n in range(n_images): + theta = n * (2 * np.pi) / n_images Rt = np.eye(4) Rt[:3, 3] = t - Rt[:3, : - 3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_axis_angle( - [0, theta, 0]) - print(f"{Rt=}") + Rt[:3, :3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_zyx( + [np.pi, theta, 0]) Rts.append(Rt) o3dvis.setup_camera(K, Rt, width, height) o3dvis.post_redraw() - o3dvis.export_current_image(f"render-{n}.jpg") - images.append(f"render-{n}.jpg") + o3dvis.export_current_image(f"render-{n:02}.jpg") + images.append(f"render-{n:02}.jpg") + print('.', end='', flush=True) np.savez("cameras.npz", width=width, height=height, K=K, Rts=Rts, images=images) + # Now create a movie from the saved images by calling ffmpeg with + # subprocess + if movie: + print("\nCreating movie...", end='', flush=True) + sp.run([ + "ffmpeg", "-framerate", f"{n_images/6}", "-pattern_type", + "glob", "-i", "render-*.jpg", "-y", meshfile.stem + ".mp4" + ], + check=True) + print("\nDone.") - o3d.visualization.draw([model], + o3d.visualization.draw([{ + 'geometry': model, + 'name': meshfile.name, + 'material': unlit + }], show_ui=False, - width=width, - height=height, - actions=[("Save Images", move_camera_and_shoot)]) + width=int(width / SCALING), + height=int(height / SCALING), + actions=[("Save Images", rotate_camera_and_shoot)]) # Linux only :-( # render = rendering.OffscreenRenderer(width, height) @@ -55,14 +75,16 @@ def albedo_from_images(meshfile, calib_data_file): tmeshes = o3d.t.geometry.TriangleMesh.from_triangle_mesh_model(model) tmeshes = list(tmeshes.values()) calib = np.load(calib_data_file) - Ks = list(Tensor(calib["K"]) for _ in range(len(calib["Rts"])))[:2] - Rts = list(Tensor(Rt) for Rt in calib["Rts"])[:2] - images = list(o3d.t.io.read_image(imfile) for imfile in calib["images"])[:2] + Ks = list(Tensor(calib["K"]) for _ in range(len(calib["Rts"]))) + Rts = list(Tensor(Rt) for Rt in calib["Rts"]) + images = list(o3d.t.io.read_image(imfile) for imfile in calib["images"]) calib.close() - tmeshes[0].project_images_to_albedo(images, Ks, Rts, 256) + # breakpoint() + albedo = tmeshes[0].project_images_to_albedo(images, Ks, Rts, 1024) + o3d.t.io.write_image("albedo.png", albedo) cam_vis = list({ "name": - f"camera-{i}", + f"camera-{i:02}", "geometry": o3d.geometry.LineSet.create_camera_visualization( images[0].columns, images[0].rows, K.numpy(), Rt.numpy(), 0.1) @@ -72,6 +94,7 @@ def albedo_from_images(meshfile, calib_data_file): "geometry": tmeshes[0] }], show_ui=True) + o3d.t.io.write_triangle_mesh(meshfile.stem + "_albedo.glb", tmeshes[0]) if __name__ == "__main__": @@ -80,9 +103,19 @@ def albedo_from_images(meshfile, calib_data_file): parser.add_argument("action", choices=('create_dataset', 'albedo_from_images')) parser.add_argument("meshfile", type=Path) + parser.add_argument("--n-images", + type=int, + default=10, + help="Number of images to render.") + parser.add_argument( + "--movie", + action="store_true", + help= + "Create movie from rendered images with ffmpeg. ffmpeg must be installed and in path." + ) args = parser.parse_args() if args.action == "create_dataset": - create_dataset(args.meshfile) + create_dataset(args.meshfile, n_images=args.n_images, movie=args.movie) else: albedo_from_images(args.meshfile, "cameras.npz") diff --git a/python/open3d/visualization/draw.py b/python/open3d/visualization/draw.py index 96396a0ffd0..16027786de9 100644 --- a/python/open3d/visualization/draw.py +++ b/python/open3d/visualization/draw.py @@ -152,7 +152,7 @@ def toggle_result(o3dvis): ("Toggle truth/result", toggle_result)]) """ gui.Application.instance.initialize() - w = O3DVisualizer(title, width, height) + w = O3DVisualizer(title, int(width), int(height)) w.set_background(bg_color, bg_image) if actions is not None: From 91a242cf103f4e0f35b0a6988a25be84fe81e8ec Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Sat, 25 May 2024 23:42:46 -0700 Subject: [PATCH 6/8] MAX and AVERAGE BlendingMethods Documentation Cleanup Add download for Baluster vase model in python example. --- cpp/open3d/t/geometry/TriangleMesh.cpp | 330 +++++++++++------- cpp/open3d/t/geometry/TriangleMesh.h | 31 +- cpp/pybind/t/geometry/trianglemesh.cpp | 43 ++- cpp/tests/t/geometry/TriangleMesh.cpp | 49 +-- docs/tutorial/data/index.rst | 10 +- .../triangle_mesh_project_to_albedo.py | 73 +++- 6 files changed, 345 insertions(+), 191 deletions(-) diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index dfa4659b749..709cb71cb84 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -7,6 +7,8 @@ #include "open3d/t/geometry/TriangleMesh.h" +#include +#include #include #include #include @@ -18,7 +20,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -29,6 +33,7 @@ #include "open3d/core/Device.h" #include "open3d/core/Dtype.h" #include "open3d/core/EigenConverter.h" +#include "open3d/core/ParallelFor.h" #include "open3d/core/ShapeUtil.h" #include "open3d/core/Tensor.h" #include "open3d/core/TensorCheck.h" @@ -1349,8 +1354,6 @@ core::Tensor Project(const core::Tensor &t_xyz, // contiguous {...,3} const core::Tensor &t_intrinsic_matrix, // {3,3} const core::Tensor &t_extrinsic_matrix, // {4,4} core::Tensor &t_xy) { // contiguous {...,2} - utility::LogInfo("K: {}\ncTw: {}", t_intrinsic_matrix.ToString(), - t_extrinsic_matrix.ToString()); auto xy_shape = t_xyz.GetShape(); auto dt = t_xyz.GetDtype(); auto t_K = t_intrinsic_matrix.To(dt).Contiguous(), @@ -1386,10 +1389,11 @@ Image TriangleMesh::ProjectImagesToAlbedo( const std::vector &intrinsic_matrices, const std::vector &extrinsic_matrices, int tex_size /*=1024*/, - bool update_material /*=true*/) { + bool update_material /*=true*/, + BlendingMethod blending_method /*=MAX*/) { + const bool DEBUG = false; using core::None; using tk = core::TensorKey; - float pixel_foreshortening_threshold = 0.25f; constexpr float EPS = 1e-6; if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); @@ -1410,138 +1414,206 @@ Image TriangleMesh::ProjectImagesToAlbedo( // (u,v) -> (x,y,z) : {tex_size, tex_size, 3} core::Tensor position_map = BakeVertexAttrTextures( tex_size, {"positions"}, 1, 0, false)["positions"]; - io::WriteImage("position_map.png", - Image(((position_map + 1) * 127.5).To(core::UInt8))); + if (DEBUG) { + io::WriteImage("position_map.png", + Image(((position_map + 1) * 127.5).To(core::UInt8))); + } core::Tensor albedo = core::Tensor::Zeros({tex_size, tex_size, 4}, core::Float32); + std::mutex albedo_mutex; albedo.Slice(2, 3, 4, 1).Fill(EPS); // regularize - core::Tensor this_albedo = - core::Tensor::Empty({tex_size, tex_size, 4}, core::Float32); - core::Tensor weighted_image({}, core::Float32), uv2xy({}, core::Float32); - core::Tensor uvrays = - core::Tensor::Empty({tex_size, tex_size, 6}, core::Float32); - - /* For each image, project vertices onto 2D image and do occlusion - * culling. This creates (xim, yim) -> (u,v) mapping. Invert this - * mapping with ipp::remap to get (u, v) -> (xim, yim) -> rgb to create - * albedo texture rgb(u,v) with bilinear interpolation. Merge albedos - * from all images OR (u,v) -> (x, y, z) -> project to (xim, yim) -> rgb - * Discard if occlusion test fails / hit test for ray fails - */ + RaycastingScene rcs; rcs.AddTriangles(*this); - for (size_t i = 0; i < images.size(); ++i) { - auto i_str = std::to_string(i); - auto width = images[i].GetCols(), height = images[i].GetRows(); - core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); - core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); - - // A. Get image space weight matrix, as inverse of pixel footprint on - // the mesh. - auto rays = RaycastingScene::CreateRaysPinhole( - intrinsic_matrices[i], extrinsic_matrices[i], width, height); - core::Tensor cam_loc = - rays.GetItem({tk::Index(0), tk::Index(0), tk::Slice(0, 3, 1)}); - auto result = rcs.CastRays(rays); - // Eigen is column-major order - Eigen::Map normals_e( - result["primitive_normals"].GetDataPtr(), 3, - width * height); - Eigen::Map rays_e(rays.GetDataPtr(), 6, - width * height); - Eigen::Map t_hit(result["t_hit"].GetDataPtr(), - 1, width * height); - auto depth = t_hit * rays_e.bottomRows<3>().colwise().norm().array(); - // removing this eval() increase runtime a lot (?) - auto rays_dir = rays_e.bottomRows<3>().colwise().normalized().eval(); - auto pixel_foreshortening = (normals_e * rays_dir) - .colwise() - .sum() - .abs(); // ignore face orientation - auto footprint_err = - (depth * (pixel_foreshortening < pixel_foreshortening_threshold) - .select(INFINITY, pixel_foreshortening)); - // fix for bad normals - auto footprint = - footprint_err.isNaN().select(INFINITY, footprint_err).eval(); - - utility::LogInfo("Image {}, footprint range: {}-{}", i, - footprint.minCoeff(), footprint.maxCoeff()); - if (!weighted_image.GetShape().IsCompatible({height, width, 4})) { - weighted_image = core::Tensor({height, width, 4}, core::Float32); + // Simulate thread_local Tensors with a vector of Tensors. thread_local + // variables are only destructed when the TBB thread pool is finalized which + // can cause memory leaks in Python code / long running processes. + // Also TBB can schedule multiple tasks on one thread at the same time + size_t max_workers = tbb::this_task_arena::max_concurrency(); + // Tensor copy ctor does shallow copies - OK for empty tensors. + std::vector this_albedo(max_workers, + core::Tensor({}, core::Float32)), + weighted_image(max_workers, core::Tensor({}, core::Float32)), + uv2xy(max_workers, core::Tensor({}, core::Float32)), + uvrays(max_workers, core::Tensor({}, core::Float32)); + + auto project_one_image = [&](const tbb::blocked_range &range) { + size_t widx = tbb::task_arena::current_thread_index(); + // initialize thread variables + if (!this_albedo[widx].GetShape().IsCompatible( + {tex_size, tex_size, 4})) { + this_albedo[widx] = + core::Tensor::Empty({tex_size, tex_size, 4}, core::Float32); + uvrays[widx] = + core::Tensor::Empty({tex_size, tex_size, 6}, core::Float32); } - weighted_image.Slice(2, 0, 3, 1) = - images[i].To(core::Float32).AsTensor(); // range: [0,1] - Eigen::Map weighted_image_e( - weighted_image.GetDataPtr(), 4, width * height); - // footprint is inf if depth or t_hit is inf. Never zero. - weighted_image_e.bottomRows<1>() = 1. / footprint; - - // B. Get texture space (u,v) -> (x,y) map and valid domain in uv space. - utility::LogInfo("cam_loc={}", cam_loc.ToString()); - uvrays.GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), - tk::Slice(0, 3, 1)}) = cam_loc; - uvrays.GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), - tk::Slice(3, 6, 1)}) = position_map - cam_loc; - auto uvresult = rcs.CastRays(uvrays); - auto &t_hit_uv = uvresult["t_hit"]; - io::WriteImage("t_hit_uv_" + i_str + ".png", - Image((t_hit_uv * 255).To(core::UInt8))); - - Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], - uv2xy); // {ts, ts, 2} - // Disable self-occluded points - for (float *p_uv2xy = uv2xy.GetDataPtr(), - *p_t_hit = t_hit_uv.GetDataPtr(); - p_uv2xy < uv2xy.GetDataPtr() + uv2xy.NumElements(); - p_uv2xy += 2, ++p_t_hit) { - if (*p_t_hit < 1 - EPS) *p_uv2xy = *(p_uv2xy + 1) = -1.f; + for (size_t i = range.begin(); i < range.end(); ++i) { + auto i_str = std::to_string(i); + auto width = images[i].GetCols(), height = images[i].GetRows(); + if (!weighted_image[widx].GetShape().IsCompatible( + {height, width, 4})) { + weighted_image[widx] = + core::Tensor({height, width, 4}, core::Float32); + } + core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); + core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); + + // A. Get image space weight matrix, as inverse of pixel + // footprint on the mesh. + auto rays = RaycastingScene::CreateRaysPinhole( + intrinsic_matrices[i], extrinsic_matrices[i], width, + height); + core::Tensor cam_loc = rays.GetItem( + {tk::Index(0), tk::Index(0), tk::Slice(0, 3, 1)}); + + // A nested parallel_for's threads must be isolated from the threads + // running this paralel_for, else we get BAD ACCESS errors. + auto result = tbb::this_task_arena::isolate( + [&rays, &rcs]() { return rcs.CastRays(rays); }); + // Eigen is column-major order + Eigen::Map normals_e( + result["primitive_normals"].GetDataPtr(), 3, + width * height); + Eigen::Map rays_e(rays.GetDataPtr(), 6, + width * height); + Eigen::Map t_hit( + result["t_hit"].GetDataPtr(), 1, width * height); + auto depth = + t_hit * rays_e.bottomRows<3>().colwise().norm().array(); + // removing this eval() increase runtime a lot (?) + auto rays_dir = + rays_e.bottomRows<3>().colwise().normalized().eval(); + auto pixel_foreshortening = + (normals_e * rays_dir) + .colwise() + .sum() + .abs(); // ignore face orientation + // fix for bad normals + auto inv_footprint = pixel_foreshortening.isNaN().select( + 0, pixel_foreshortening) / + (depth * depth); + + utility::LogDebug( + "[ProjectImagesToAlbedo] Image {}, weight (inv_footprint) " + "range: {}-{}", + i, inv_footprint.minCoeff(), inv_footprint.maxCoeff()); + weighted_image[widx].Slice(2, 0, 3, 1) = + images[i].To(core::Float32).AsTensor(); // range: [0,1] + Eigen::Map weighted_image_e( + weighted_image[widx].GetDataPtr(), 4, + width * height); + weighted_image_e.bottomRows<1>() = inv_footprint; + + // B. Get texture space (u,v) -> (x,y) map and valid domain in + // uv space. + uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(0, 3, 1)}) = cam_loc; + uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(3, 6, 1)}) = position_map - cam_loc; + // A nested parallel_for's threads must be isolated from the threads + // running this paralel_for, else we get BAD ACCESS errors. + result = tbb::this_task_arena::isolate([&rcs, &uvrays, widx]() { + return rcs.CastRays(uvrays[widx]); + }); + auto &t_hit_uv = result["t_hit"]; + if (DEBUG) { + io::WriteImage("t_hit_uv_" + i_str + ".png", + Image((t_hit_uv * 255).To(core::UInt8))); + } + + Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], + uv2xy[widx]); // {ts, ts, 2} + // Disable self-occluded points + for (float *p_uv2xy = uv2xy[widx].GetDataPtr(), + *p_t_hit = t_hit_uv.GetDataPtr(); + p_uv2xy < + uv2xy[widx].GetDataPtr() + uv2xy[widx].NumElements(); + p_uv2xy += 2, ++p_t_hit) { + if (*p_t_hit < 1 - EPS) *p_uv2xy = *(p_uv2xy + 1) = -1.f; + } + core::Tensor uv2xy2 = + uv2xy[widx].Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} + if (DEBUG) { + io::WriteImage("uv2x_" + i_str + ".png", + Image((uv2xy2[0].To(core::UInt16)))); + io::WriteImage("uv2y_" + i_str + ".png", + Image((uv2xy2[1].To(core::UInt16)))); + } + // albedo[u,v] = image[ i[u,v], j[u,v] ] + // ij[u,v] ? = identity[ uv[i,j] ] + + // C. Interpolate weighted image to weighted texture + this_albedo[widx].Fill(0.f); + ipp::Remap(weighted_image[widx], /*{height, width, 4} f32*/ + uv2xy2[0], /* {texsz, texsz} f32*/ + uv2xy2[1], /* {texsz, texsz} f32*/ + this_albedo[widx], /*{texsz, texsz, 4} f32*/ + t::geometry::Image::InterpType::Linear); + // Weights can become negative with higher order interpolation + if (DEBUG) { + io::WriteImage("this_albedo_" + i_str + ".png", + Image((this_albedo[widx].Slice(2, 0, 3, 1) * 255) + .To(core::UInt8))); + float wtmin = this_albedo[widx] + .Slice(2, 3, 4, 1) + .Min({0, 1, 2}) + .Item(), + wtmax = this_albedo[widx] + .Slice(2, 3, 4, 1) + .Max({0, 1, 2}) + .Item(); + io::WriteImage("image_weights_" + i_str + ".png", + Image(weighted_image[widx] + .Slice(2, 3, 4, 1) + .Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / + (wtmax - wtmin))); + io::WriteImage( + "this_albedo_weight_" + i_str + ".png", + Image(this_albedo[widx].Slice(2, 3, 4, 1).Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / + (wtmax - wtmin))); + } + { + std::lock_guard albedo_lock{albedo_mutex}; + if (blending_method == BlendingMethod::MAX) { + // Select albedo value with max weight + for (auto p_albedo = albedo.GetDataPtr(), + p_this_albedo = + this_albedo[widx].GetDataPtr(); + p_albedo < + albedo.GetDataPtr() + albedo.NumElements(); + p_albedo += 4, p_this_albedo += 4) { + if (p_albedo[3] < p_this_albedo[3]) { + for (auto k = 0; k < 4; ++k) + p_albedo[k] = p_this_albedo[k]; + } + } + } else if (blending_method == BlendingMethod::AVERAGE) { + for (auto p_albedo = albedo.GetDataPtr(), + p_this_albedo = + this_albedo[widx].GetDataPtr(); + p_albedo < + albedo.GetDataPtr() + albedo.NumElements(); + p_albedo += 4, p_this_albedo += 4) { + for (auto k = 0; k < 3; ++k) + p_albedo[k] += p_this_albedo[k] * p_this_albedo[3]; + p_albedo[3] += p_this_albedo[3]; + } + } + } } - uv2xy = uv2xy.Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} - io::WriteImage("uv2x_" + i_str + ".png", - Image((uv2xy[0].To(core::UInt16)))); - io::WriteImage("uv2y_" + i_str + ".png", - Image((uv2xy[1].To(core::UInt16)))); - // albedo[u,v] = image[ i[u,v], j[u,v] ] - // ij[u,v] ? = identity[ uv[i,j] ] - - // C. Interpolate weighted image to weighted texture - this_albedo.Fill(0.f); - ipp::Remap(weighted_image, /*{height, width, 4} f32*/ - uv2xy[0], /* {texsz, texsz} f32*/ - uv2xy[1], /* {texsz, texsz} f32*/ - this_albedo, /*{texsz, texsz, 4} f32*/ - t::geometry::Image::InterpType::Linear); - // Weights can become negative with higher order interpolation - io::WriteImage( - "this_albedo_" + i_str + ".png", - Image((this_albedo.Slice(2, 0, 3, 1) * 255).To(core::UInt8))); - float wtmin = this_albedo.Slice(2, 3, 4, 1) - .Min({0, 1, 2}) - .Item(), - wtmax = this_albedo.Slice(2, 3, 4, 1) - .Max({0, 1, 2}) - .Item(); - utility::LogInfo("Image {}, weight range: {}-{}", i, wtmin, wtmax); - io::WriteImage( - "image_weights_" + i_str + ".png", - Image(weighted_image.Slice(2, 3, 4, 1).Contiguous()) - .To(core::UInt8, /*copy=*/false, - /*scale=*/255.f / (wtmax - wtmin), - /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); - io::WriteImage( - "this_albedo_weight_" + i_str + ".png", - Image(this_albedo.Slice(2, 3, 4, 1).Contiguous()) - .To(core::UInt8, /*copy=*/false, - /*scale=*/255.f / (wtmax - wtmin), - /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); - // weighted sum - albedo.Slice(2, 0, 3, 1) += - this_albedo.Slice(2, 0, 3, 1) * this_albedo.Slice(2, 3, 4, 1); - albedo.Slice(2, 3, 4, 1) += this_albedo.Slice(2, 3, 4, 1); - } - albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); + }; + tbb::parallel_for(tbb::blocked_range(0, images.size()), + project_one_image); + if (blending_method == BlendingMethod::AVERAGE) { + albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); + } Image albedo_texture((albedo.Slice(2, 0, 3, 1) * 255).To(core::UInt8)); if (update_material) { @@ -1549,7 +1621,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( SetMaterial(visualization::rendering::Material()); GetMaterial().SetDefaultProperties(); // defaultUnlit } - GetMaterial().SetTextureMap("albedo", albedo_texture); + GetMaterial().SetAlbedoMap(albedo_texture); } return albedo_texture; } diff --git a/cpp/open3d/t/geometry/TriangleMesh.h b/cpp/open3d/t/geometry/TriangleMesh.h index 2dbb71ec55d..a8402bbd4b1 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.h +++ b/cpp/open3d/t/geometry/TriangleMesh.h @@ -26,6 +26,9 @@ namespace geometry { class LineSet; class RaycastingScene; +/// Texture Blending method for ProjectImagesToAlbedo() from overlapping images. +enum class BlendingMethod { MAX, AVERAGE }; + /// \class TriangleMesh /// \brief A triangle mesh contains vertices and triangles. /// @@ -979,12 +982,38 @@ class TriangleMesh : public Geometry, public DrawableGeometry { /// \return The reference to itself. TriangleMesh RemoveUnreferencedVertices(); + /// Create an albedo for the triangle mesh using calibrated images. The + /// triangle mesh must have texture coordinates ("texture_uvs" triangle + /// attribute). This works by back projecting the images onto the texture + /// surface. Overlapping images are blended together in the resulting + /// albedo. For best results, use images captured with exposure and white + /// balance lock to reduce the chance of seams in the output texture. + /// + /// \param images vector of images. + /// \param intrinsic_matrices vector of {3,3} intrinsic matrices describing + /// the pinhole camera. + /// \param extrinsic_matrices vector of {4,4} extrinsic matrices describing + /// the position and orientation of the camera. + /// \param tex_size Output albedo texture size. This is a square image, so + /// only one side is needed. + /// \param update_material Whether to update the material of the triangle + /// mesh, possibly overwriting an existing albedo texture. + /// \param blending_method BlendingMethod enum specifying the blending + /// method for overlapping images: + /// - MAX: For each texel, pick the input pixel with the max weight from + /// all overlapping images. This creates sharp textures but may have + /// visible seams. + /// - AVERAGE: The output texel value is the weighted sum of input + /// pixels. This creates smooth blending without seams, but the results + /// may be blurry. + /// \return Image with albedo texture Image ProjectImagesToAlbedo( const std::vector &images, const std::vector &intrinsic_matrices, const std::vector &extrinsic_matrices, int tex_size = 1024, - bool update_material = true); + bool update_material = true, + BlendingMethod blending_method = BlendingMethod::MAX); protected: core::Device device_ = core::Device("CPU:0"); diff --git a/cpp/pybind/t/geometry/trianglemesh.cpp b/cpp/pybind/t/geometry/trianglemesh.cpp index fabc79dc38e..1c741e4f025 100644 --- a/cpp/pybind/t/geometry/trianglemesh.cpp +++ b/cpp/pybind/t/geometry/trianglemesh.cpp @@ -858,13 +858,6 @@ This function always uses the CPU device. plt.imshow(texture_tensors['albedo'].numpy()) )"); - triangle_mesh.def( - "project_images_to_albedo", &TriangleMesh::ProjectImagesToAlbedo, - "images"_a, "intrinsic_matrices"_a, "extrinsic_matrices"_a, - "tex_size"_a = 1024, "update_material"_a = true, - py::call_guard(), - R"(Create an albedo texture from images of an object taken with a calibrated camera.)"); - triangle_mesh.def("extrude_rotation", &TriangleMesh::ExtrudeRotation, "angle"_a, "axis"_a, "resolution"_a = 16, "translation"_a = 0.0, "capping"_a = true, @@ -1000,10 +993,44 @@ or has a negative value, it is ignored. &TriangleMesh::RemoveUnreferencedVertices, "Removes unreferenced vertices from the mesh in-place."); + py::enum_(m, "BlendingMethod") + .value("MAX", BlendingMethod::MAX) + .value("AVERAGE", BlendingMethod::AVERAGE); triangle_mesh.def("project_images_to_albedo", &TriangleMesh::ProjectImagesToAlbedo, "images"_a, "intrinsic_matrices"_a, "extrinsic_matrices"_a, - "tex_size"_a = 1024, "update_material"_a = true); + "tex_size"_a = 1024, "update_material"_a = true, + "blending_method"_a = BlendingMethod::MAX, + py::call_guard(), R"( +Create an albedo for the triangle mesh using calibrated images. The triangle +mesh must have texture coordinates ("texture_uvs" triangle attribute). This works +by back projecting the images onto the texture surface. Overlapping images are +blended together in the resulting albedo. For best results, use images captured +with exposure and white balance lock to reduce the chance of seams in the output +texture. + +Args: + images (List[open3d.t.geometry.Image]): List of images. + intrinsic_matrices (List[open3d.core.Tensor]): List of (3,3) intrinsic matrices describing + the pinhole camera. + extrinsic_matrices (List[open3d.core.Tensor]): List of (4,4) extrinsic matrices describing + the position and orientation of the camera. + tex_size (int): Output albedo texture size. This is a square image, so + only one side is needed. + update_material (bool): Whether to update the material of the triangle + mesh, possibly overwriting an existing albedo texture. + blending_method (BlendingMethod) enum specifying the blending + method for overlapping images:: + + - `MAX`: For each texel, pick the input pixel with the max weight from + all overlapping images. This creates sharp textures but may have + visible seams. + - `AVERAGE`: The output texel value is the weighted sum of input + pixels. This creates smooth blending without seams, but the results + may be blurry. + +Returns: + Image with albedo texture.)"); } } // namespace geometry diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index a431804da9e..f2640d67cdf 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -8,10 +8,12 @@ #include "open3d/t/geometry/TriangleMesh.h" #include +#include #include "core/CoreTest.h" #include "open3d/core/Dtype.h" #include "open3d/core/EigenConverter.h" +#include "open3d/core/SizeVector.h" #include "open3d/core/Tensor.h" #include "open3d/core/TensorCheck.h" #include "open3d/geometry/LineSet.h" @@ -1350,7 +1352,6 @@ TEST_P(TriangleMeshPermuteDevices, RemoveUnreferencedVertices) { TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { using namespace t::geometry; - using ls = open3d::geometry::LineSet; core::Device device = GetParam(); TriangleMesh sphere = TriangleMesh::FromLegacy(*geometry::TriangleMesh::CreateSphere( @@ -1379,46 +1380,22 @@ TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { device), }; - Eigen::Map e_intrinsic( - intrinsic_matrix.GetDataPtr()); - Eigen::Map e_extrinsic[3] = { - Eigen::Map( - extrinsic_matrix[0].GetDataPtr()), - Eigen::Map( - extrinsic_matrix[1].GetDataPtr()), - Eigen::Map( - extrinsic_matrix[2].GetDataPtr())}; - std::shared_ptr p_camera[3] = { - ls::CreateCameraVisualization( - 256, 192, e_intrinsic.transpose().cast(), - e_extrinsic[0].transpose().cast()), - ls::CreateCameraVisualization( - 256, 192, e_intrinsic.transpose().cast(), - e_extrinsic[1].transpose().cast()), - ls::CreateCameraVisualization( - 256, 192, e_intrinsic.transpose().cast(), - e_extrinsic[2].transpose().cast())}; - Image albedo = sphere.ProjectImagesToAlbedo( {Image(view[0]), Image(view[1]), Image(view[2])}, {intrinsic_matrix, intrinsic_matrix, intrinsic_matrix}, {extrinsic_matrix[0], extrinsic_matrix[1], extrinsic_matrix[2]}, 256, true); - utility::LogInfo("Mesh: {}", sphere.ToString()); - utility::LogInfo("Texture: {}", albedo.ToString()); - t::io::WriteImage("albedo.png", albedo); - /* t::io::WriteTriangleMesh("sphere-projected.obj", sphere); */ - t::io::WriteTriangleMesh("sphere-projected.glb", sphere); - t::io::WriteTriangleMesh("sphere-projected.npz", sphere); - - visualization::Draw( - {visualization::DrawObject("camera_0", p_camera[0], true), - visualization::DrawObject("camera_1", p_camera[1], true), - visualization::DrawObject("camera_2", p_camera[2], true), - visualization::DrawObject{ - "mesh", std::make_shared(std::move(sphere)), - true}}, - "ProjectImagesToAlbedo", 1024, 768); + + EXPECT_TRUE(sphere.HasMaterial()); + EXPECT_TRUE(sphere.GetMaterial().HasAlbedoMap()); + EXPECT_TRUE(albedo.AsTensor().GetShape().IsCompatible({256, 256, 3})); + EXPECT_TRUE(albedo.GetDtype() == core::UInt8); + core::Tensor mean_color_ref = + core::Tensor::Init({92.465515, 71.62926, 67.55928}); + EXPECT_TRUE(albedo.AsTensor() + .To(core::Float32) + .Mean({0, 1}) + .AllClose(mean_color_ref)); } } // namespace tests diff --git a/docs/tutorial/data/index.rst b/docs/tutorial/data/index.rst index e17c8d4619e..f1fcf2dd21c 100644 --- a/docs/tutorial/data/index.rst +++ b/docs/tutorial/data/index.rst @@ -175,13 +175,13 @@ A 3D Mobius knot mesh in PLY format. data::KnotMesh dataset; auto mesh = io::CreateMeshFromFile(dataset.GetPath()); -TriangleModel with PRB texture +TriangleModel with PBR texture ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ MonkeyModel ----------- -The monkey model with PRB texture. +The monkey model with PBR texture. .. code-block:: python @@ -197,7 +197,7 @@ The monkey model with PRB texture. SwordModel ---------- -The sword model with PRB texture. +The sword model with PBR texture. .. code-block:: python @@ -213,7 +213,7 @@ The sword model with PRB texture. CrateModel ---------- -The crate model with PRB texture. +The crate model with PBR texture. .. code-block:: python @@ -229,7 +229,7 @@ The crate model with PRB texture. FlightHelmetModel ----------------- -The flight helmet gltf model with PRB texture. +The flight helmet gltf model with PBR texture. .. code-block:: python diff --git a/examples/python/geometry/triangle_mesh_project_to_albedo.py b/examples/python/geometry/triangle_mesh_project_to_albedo.py index a8ed06e552c..0f057da3c8c 100644 --- a/examples/python/geometry/triangle_mesh_project_to_albedo.py +++ b/examples/python/geometry/triangle_mesh_project_to_albedo.py @@ -1,22 +1,62 @@ +# ---------------------------------------------------------------------------- +# - Open3D: www.open3d.org - +# ---------------------------------------------------------------------------- +# Copyright (c) 2018-2023 www.open3d.org +# SPDX-License-Identifier: MIT +# ---------------------------------------------------------------------------- +"""This example demonstrates project_image_to_albedo. Use create_dataset mode to +render images of a 3D mesh or model from different viewpoints. +albedo_from_dataset mode then uses the calibrated images to re-create the albedo +texture for the mesh. +""" import argparse from pathlib import Path import subprocess as sp +import time import numpy as np import open3d as o3d from open3d.visualization import gui, rendering, O3DVisualizer from open3d.core import Tensor +def download_smithsonian_baluster_vase(): + """Download the Smithsonian Baluster Vase 3D model.""" + vase_url = 'https://3d-api.si.edu/content/document/3d_package:d8c62634-4ebc-11ea-b77f-2e728ce88125/resources/F1980.190%E2%80%93194_baluster_vase-150k-4096.glb' + import urllib.request + + def show_progress(block_num, block_size, total_size): + total_size = total_size >> 20 if total_size > 0 else "??" # Convert to MB if known + print( + "Downloading F1980_baluster_vase.glb... " + f"{(block_num * block_size) >>20}MB / {total_size}MB", + end="\r") + + urllib.request.urlretrieve(vase_url, + filename="F1980_baluster_vase.glb", + reporthook=show_progress) + print("\nDownload complete.") + + def create_dataset(meshfile, n_images=10, movie=False): + """Render images of a 3D mesh from different viewpoints. These form a + synthetic dataset to test the project_images_to_albedo function. + """ + # Adjust these parameters to properly frame your model. + # Window system pixel scaling (e.g. 1 for normal, 2 for HiDPI / retina display) SCALING = 2 - width, height = 1024, 1024 + width, height = 1024, 1024 # image width, height focal_length = 512 + d_camera_obj = 0.3 # distance from camera to object K = np.array([[focal_length, 0, width / 2], [0, focal_length, height / 2], [0, 0, 1]]) + t = np.array([0, 0, d_camera_obj]) # origin / object in camera ref frame + model = o3d.io.read_triangle_model(meshfile) + # DefaultLit shader will produce non-uniform images with specular + # highlights, etc. These should be avoided to accurately capture the diffuse + # albedo unlit = rendering.MaterialRecord() unlit.shader = "unlit" - t = np.array([0, 0, 0.3]) # origin / object in camera ref frame def rotate_camera_and_shoot(o3dvis): Rts = [] @@ -52,6 +92,8 @@ def rotate_camera_and_shoot(o3dvis): check=True) print("\nDone.") + print("If the object is properly framed in the GUI window, click on the " + "'Save Images' action in the menu.") o3d.visualization.draw([{ 'geometry': model, 'name': meshfile.name, @@ -62,12 +104,6 @@ def rotate_camera_and_shoot(o3dvis): height=int(height / SCALING), actions=[("Save Images", rotate_camera_and_shoot)]) - # Linux only :-( - # render = rendering.OffscreenRenderer(width, height) - # render.scene.add_geometry(model) - # img = render.render_to_image() - # o3d.io.write_image("render-image.jpg", img) - def albedo_from_images(meshfile, calib_data_file): @@ -79,9 +115,11 @@ def albedo_from_images(meshfile, calib_data_file): Rts = list(Tensor(Rt) for Rt in calib["Rts"]) images = list(o3d.t.io.read_image(imfile) for imfile in calib["images"]) calib.close() - # breakpoint() + start = time.time() albedo = tmeshes[0].project_images_to_albedo(images, Ks, Rts, 1024) + print(f"project_images_to_albedo ran in {time.time()-start:.2f}s") o3d.t.io.write_image("albedo.png", albedo) + o3d.t.io.write_triangle_mesh(meshfile.stem + "_albedo.glb", tmeshes[0]) cam_vis = list({ "name": f"camera-{i:02}", @@ -94,19 +132,24 @@ def albedo_from_images(meshfile, calib_data_file): "geometry": tmeshes[0] }], show_ui=True) - o3d.t.io.write_triangle_mesh(meshfile.stem + "_albedo.glb", tmeshes[0]) if __name__ == "__main__": - parser = argparse.ArgumentParser() + parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("action", choices=('create_dataset', 'albedo_from_images')) - parser.add_argument("meshfile", type=Path) + parser.add_argument("--meshfile", + type=Path, + default=".", + help="Path to mesh file.") parser.add_argument("--n-images", type=int, default=10, help="Number of images to render.") + parser.add_argument("--download_sample_model", + help="Download a sample 3D model for this example.", + action="store_true") parser.add_argument( "--movie", action="store_true", @@ -116,6 +159,12 @@ def albedo_from_images(meshfile, calib_data_file): args = parser.parse_args() if args.action == "create_dataset": + if args.download_sample_model: + download_smithsonian_baluster_vase() + args.meshfile = "F1980_baluster_vase.glb" + if args.meshfile == Path("."): + parser.error("Please provide a path to a mesh file, or use " + "--download_sample_model.") create_dataset(args.meshfile, n_images=args.n_images, movie=args.movie) else: albedo_from_images(args.meshfile, "cameras.npz") From 86a0bf4b0bacb29109685626bc339f9c3d2d42e8 Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Fri, 31 May 2024 14:44:40 -0700 Subject: [PATCH 7/8] Color Correction (WIP) Needs tbb::feeder from tbb::parallel_for_each --- cpp/open3d/t/geometry/TriangleMesh.cpp | 439 +++++++++++------- cpp/open3d/t/geometry/TriangleMesh.h | 24 +- cpp/pybind/t/geometry/trianglemesh.cpp | 23 +- cpp/tests/t/geometry/TriangleMesh.cpp | 10 +- .../triangle_mesh_project_to_albedo.py | 34 +- .../python/visualization/remove_geometry.py | 3 +- 6 files changed, 359 insertions(+), 174 deletions(-) diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 709cb71cb84..0b583c85861 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -9,6 +9,8 @@ #include #include +#include +#include #include #include #include @@ -19,6 +21,7 @@ #include #include +#include #include #include #include @@ -1382,6 +1385,98 @@ core::Tensor Project(const core::Tensor &t_xyz, // contiguous {...,3} }); return t_xy; // contiguous {...,2} } + +/// Estimate contrast and brightness in each color channel for this_albedo to +/// match albedo, based on the overlapping area in the texture image. +std::tuple, + std::array, + core::Tensor, + core::Tensor> +get_color_correction(const core::Tensor &albedo, + const core::Tensor &this_albedo, + float SPECULAR_THRESHOLD) { + constexpr double EPS = 1e-6; + std::array this_contrast{}, this_brightness{}, sum{}, this_sum{}, + absdiffsum{EPS, EPS, EPS}, this_absdiffsum{EPS, EPS, EPS}; + auto q_albedo = + Image(albedo).Resize(0.25, Image::InterpType::Linear).AsTensor(); + auto q_this_albedo = Image(this_albedo) + .Resize(0.25, Image::InterpType::Linear) + .AsTensor(); + int64_t rowstride = albedo.GetStride(0), row = 0, col = 0; + for (float *pq_albedo = q_albedo.GetDataPtr(), + *pq_this_albedo = q_this_albedo.GetDataPtr(); + pq_albedo < q_albedo.GetDataPtr() + q_albedo.NumElements(); + pq_albedo += 4, pq_this_albedo += 4, ++col) { + if (col == albedo.GetShape(1)) { + col = 0; + ++row; + } + if (pq_albedo[3] <= EPS || pq_this_albedo[3] <= EPS) continue; + bool update_absdiffsum = row > 0 && row < albedo.GetShape(0) - 1 && + col > 0 && col < albedo.GetShape(1) - 1; + for (int c = 0; c < 3; ++c) { + sum[c] += pq_albedo[c]; + this_sum[c] += pq_this_albedo[c]; + // Absolute central difference + if (update_absdiffsum) { + absdiffsum[c] += std::fabs( + pq_albedo[c] - 0.25 * (pq_albedo[c - 4] + + pq_albedo[c - 4 - rowstride] + + pq_albedo[c + 4] + + pq_albedo[c + rowstride + 4])); + this_absdiffsum[c] += + std::fabs(pq_this_albedo[c] - + 0.25 * (pq_this_albedo[c - 4] + + pq_this_albedo[c - 4 - rowstride] + + pq_this_albedo[c + 4] + + pq_this_albedo[c + rowstride + 4])); + } + } + } + for (int c = 0; c < 3; ++c) { + this_contrast[c] = absdiffsum[c] / this_absdiffsum[c]; + this_brightness[c] = sum[c] - this_contrast[c] * this_sum[c]; + } + + // Detect specular highlights from outliers. + core::Tensor highlights_mask = core::Tensor::Empty( + {albedo.GetShape(0) / 4, albedo.GetShape(1) / 4, 1}, + core::Bool), + this_highlights_mask = core::Tensor::Empty( + {albedo.GetShape(0) / 4, albedo.GetShape(1) / 4, 1}, + core::Bool); + std::array err{}; + bool *p_highlights_mask = highlights_mask.GetDataPtr(), + *p_this_highlights_mask = this_highlights_mask.GetDataPtr(); + for (float *pq_albedo = q_albedo.GetDataPtr(), + *pq_this_albedo = q_this_albedo.GetDataPtr(); + pq_albedo < q_albedo.GetDataPtr() + q_albedo.NumElements(); + pq_albedo += 4, pq_this_albedo += 4, ++p_highlights_mask, + ++p_this_highlights_mask) { + if (pq_albedo[3] <= EPS || pq_this_albedo[3] <= EPS) continue; + for (int c = 0; c < 3; ++c) { + err[c] = 1.f - (pq_this_albedo[c] * this_contrast[c] + + this_brightness[c]) / + pq_albedo[c]; + } + *p_highlights_mask = + std::min({err[0], err[1], err[2]}) > SPECULAR_THRESHOLD; + *p_this_highlights_mask = + std::max({err[0], err[1], err[2]}) < -SPECULAR_THRESHOLD; + } + highlights_mask = Image(highlights_mask) + .Dilate(/*kernel_size=*/3) + .Resize(4.f) + .AsTensor(); + this_highlights_mask = Image(this_highlights_mask) + .Dilate(/*kernel_size=*/3) + .Resize(4.f) + .AsTensor(); + return std::make_tuple(this_contrast, this_brightness, highlights_mask, + this_highlights_mask); +} + } // namespace Image TriangleMesh::ProjectImagesToAlbedo( @@ -1395,6 +1490,7 @@ Image TriangleMesh::ProjectImagesToAlbedo( using core::None; using tk = core::TensorKey; constexpr float EPS = 1e-6; + constexpr float SPECULAR_THRESHOLD = 0.1f; if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); } @@ -1438,8 +1534,13 @@ Image TriangleMesh::ProjectImagesToAlbedo( uv2xy(max_workers, core::Tensor({}, core::Float32)), uvrays(max_workers, core::Tensor({}, core::Float32)); - auto project_one_image = [&](const tbb::blocked_range &range) { - size_t widx = tbb::task_arena::current_thread_index(); + // Used to control order of blending projected images into the texture. This + // ensures correct computation of color matching (only neded for + // COLOR_CORRECTION). + std::condition_variable cv_next_blend_image; + std::atomic_uint next_blend_image{0}; + auto project_one_image = [&](size_t i, tbb::feeder &feeder) { + size_t widx = tbb::this_task_arena::current_thread_index(); // initialize thread variables if (!this_albedo[widx].GetShape().IsCompatible( {tex_size, tex_size, 4})) { @@ -1448,170 +1549,192 @@ Image TriangleMesh::ProjectImagesToAlbedo( uvrays[widx] = core::Tensor::Empty({tex_size, tex_size, 6}, core::Float32); } - for (size_t i = range.begin(); i < range.end(); ++i) { - auto i_str = std::to_string(i); - auto width = images[i].GetCols(), height = images[i].GetRows(); - if (!weighted_image[widx].GetShape().IsCompatible( - {height, width, 4})) { - weighted_image[widx] = - core::Tensor({height, width, 4}, core::Float32); - } - core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); - core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); - - // A. Get image space weight matrix, as inverse of pixel - // footprint on the mesh. - auto rays = RaycastingScene::CreateRaysPinhole( - intrinsic_matrices[i], extrinsic_matrices[i], width, - height); - core::Tensor cam_loc = rays.GetItem( - {tk::Index(0), tk::Index(0), tk::Slice(0, 3, 1)}); - - // A nested parallel_for's threads must be isolated from the threads - // running this paralel_for, else we get BAD ACCESS errors. - auto result = tbb::this_task_arena::isolate( - [&rays, &rcs]() { return rcs.CastRays(rays); }); - // Eigen is column-major order - Eigen::Map normals_e( - result["primitive_normals"].GetDataPtr(), 3, - width * height); - Eigen::Map rays_e(rays.GetDataPtr(), 6, - width * height); - Eigen::Map t_hit( - result["t_hit"].GetDataPtr(), 1, width * height); - auto depth = - t_hit * rays_e.bottomRows<3>().colwise().norm().array(); - // removing this eval() increase runtime a lot (?) - auto rays_dir = - rays_e.bottomRows<3>().colwise().normalized().eval(); - auto pixel_foreshortening = - (normals_e * rays_dir) - .colwise() - .sum() - .abs(); // ignore face orientation - // fix for bad normals - auto inv_footprint = pixel_foreshortening.isNaN().select( - 0, pixel_foreshortening) / - (depth * depth); + auto i_str = std::to_string(i); + auto width = images[i].GetCols(), height = images[i].GetRows(); + if (!weighted_image[widx].GetShape().IsCompatible({height, width, 4})) { + weighted_image[widx] = + core::Tensor({height, width, 4}, core::Float32); + } + core::AssertTensorShape(intrinsic_matrices[i], {3, 3}); + core::AssertTensorShape(extrinsic_matrices[i], {4, 4}); + + // A. Get image space weight matrix, as inverse of pixel + // footprint on the mesh. + auto rays = RaycastingScene::CreateRaysPinhole( + intrinsic_matrices[i], extrinsic_matrices[i], width, height); + core::Tensor cam_loc = + rays.GetItem({tk::Index(0), tk::Index(0), tk::Slice(0, 3, 1)}); + + // A nested parallel_for's threads must be isolated from the threads + // running this paralel_for, else we get BAD ACCESS errors. + auto result = tbb::this_task_arena::isolate( + [&rays, &rcs]() { return rcs.CastRays(rays); }); + // Eigen is column-major order + Eigen::Map normals_e( + result["primitive_normals"].GetDataPtr(), 3, + width * height); + Eigen::Map rays_e(rays.GetDataPtr(), 6, + width * height); + Eigen::Map t_hit(result["t_hit"].GetDataPtr(), + 1, width * height); + auto depth = t_hit * rays_e.bottomRows<3>().colwise().norm().array(); + // removing this eval() increase runtime a lot (?) + auto rays_dir = rays_e.bottomRows<3>().colwise().normalized().eval(); + auto pixel_foreshortening = (normals_e * rays_dir) + .colwise() + .sum() + .abs(); // ignore face orientation + // fix for bad normals + auto inv_footprint = + pixel_foreshortening.isNaN().select(0, pixel_foreshortening) / + (depth * depth); + + utility::LogDebug( + "[ProjectImagesToAlbedo] Image {}, weight (inv_footprint) " + "range: {}-{}", + i, inv_footprint.minCoeff(), inv_footprint.maxCoeff()); + weighted_image[widx].Slice(2, 0, 3, 1) = + images[i].To(core::Float32).AsTensor(); // range: [0,1] + Eigen::Map weighted_image_e( + weighted_image[widx].GetDataPtr(), 4, width * height); + weighted_image_e.bottomRows<1>() = inv_footprint; + + // B. Get texture space (u,v) -> (x,y) map and valid domain in + // uv space. + uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(0, 3, 1)}) = cam_loc; + uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), + tk::Slice(3, 6, 1)}) = position_map - cam_loc; + // A nested parallel_for's threads must be isolated from the threads + // running this paralel_for, else we get BAD ACCESS errors. + result = tbb::this_task_arena::isolate( + [&rcs, &uvrays, widx]() { return rcs.CastRays(uvrays[widx]); }); + auto &t_hit_uv = result["t_hit"]; + if (DEBUG) { + io::WriteImage("t_hit_uv_" + i_str + ".png", + Image((t_hit_uv * 255).To(core::UInt8))); + } + Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], + uv2xy[widx]); // {ts, ts, 2} + // Disable self-occluded points + for (float *p_uv2xy = uv2xy[widx].GetDataPtr(), + *p_t_hit = t_hit_uv.GetDataPtr(); + p_uv2xy < + uv2xy[widx].GetDataPtr() + uv2xy[widx].NumElements(); + p_uv2xy += 2, ++p_t_hit) { + if (*p_t_hit < 1 - EPS) *p_uv2xy = *(p_uv2xy + 1) = -1.f; + } + core::Tensor uv2xy2 = + uv2xy[widx].Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} + if (DEBUG) { + io::WriteImage("uv2x_" + i_str + ".png", + Image((uv2xy2[0].To(core::UInt16)))); + io::WriteImage("uv2y_" + i_str + ".png", + Image((uv2xy2[1].To(core::UInt16)))); + } + // albedo[u,v] = image[ i[u,v], j[u,v] ] + // ij[u,v] ? = identity[ uv[i,j] ] + + // C. Interpolate weighted image to weighted texture + this_albedo[widx].Fill(0.f); + ipp::Remap(weighted_image[widx], /*{height, width, 4} f32*/ + uv2xy2[0], /* {texsz, texsz} f32*/ + uv2xy2[1], /* {texsz, texsz} f32*/ + this_albedo[widx], /*{texsz, texsz, 4} f32*/ + t::geometry::Image::InterpType::Linear); + // Weights can become negative with higher order interpolation + if (DEBUG) { + io::WriteImage("this_albedo_" + i_str + ".png", + Image((this_albedo[widx].Slice(2, 0, 3, 1) * 255) + .To(core::UInt8))); + float wtmin = this_albedo[widx] + .Slice(2, 3, 4, 1) + .Min({0, 1, 2}) + .Item(), + wtmax = this_albedo[widx] + .Slice(2, 3, 4, 1) + .Max({0, 1, 2}) + .Item(); + io::WriteImage( + "image_weights_" + i_str + ".png", + Image(weighted_image[widx].Slice(2, 3, 4, 1).Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); + io::WriteImage( + "this_albedo_weight_" + i_str + ".png", + Image(this_albedo[widx].Slice(2, 3, 4, 1).Contiguous()) + .To(core::UInt8, /*copy=*/false, + /*scale=*/255.f / (wtmax - wtmin), + /*offset=*/-wtmin * 255.f / (wtmax - wtmin))); + } + std::array this_contrast{1.f, 1.f, 1.f}, + this_brightness{0.f, 0.f, 0.f}; + core::Tensor highlights_mask, this_highlights_mask; + if (TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION)) { + std::tie(this_contrast, this_brightness, highlights_mask, + this_highlights_mask) = + get_color_correction(albedo, this_albedo[widx], + SPECULAR_THRESHOLD); utility::LogDebug( - "[ProjectImagesToAlbedo] Image {}, weight (inv_footprint) " - "range: {}-{}", - i, inv_footprint.minCoeff(), inv_footprint.maxCoeff()); - weighted_image[widx].Slice(2, 0, 3, 1) = - images[i].To(core::Float32).AsTensor(); // range: [0,1] - Eigen::Map weighted_image_e( - weighted_image[widx].GetDataPtr(), 4, - width * height); - weighted_image_e.bottomRows<1>() = inv_footprint; - - // B. Get texture space (u,v) -> (x,y) map and valid domain in - // uv space. - uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), - tk::Slice(0, 3, 1)}) = cam_loc; - uvrays[widx].GetItem({tk::Slice(0, None, 1), tk::Slice(0, None, 1), - tk::Slice(3, 6, 1)}) = position_map - cam_loc; - // A nested parallel_for's threads must be isolated from the threads - // running this paralel_for, else we get BAD ACCESS errors. - result = tbb::this_task_arena::isolate([&rcs, &uvrays, widx]() { - return rcs.CastRays(uvrays[widx]); - }); - auto &t_hit_uv = result["t_hit"]; - if (DEBUG) { - io::WriteImage("t_hit_uv_" + i_str + ".png", - Image((t_hit_uv * 255).To(core::UInt8))); - } + "[ProjectImagesToAlbedo] Image {}, contrast: {}, " + "brightness {}", + i_str, this_contrast, this_brightness); + } - Project(position_map, intrinsic_matrices[i], extrinsic_matrices[i], - uv2xy[widx]); // {ts, ts, 2} - // Disable self-occluded points - for (float *p_uv2xy = uv2xy[widx].GetDataPtr(), - *p_t_hit = t_hit_uv.GetDataPtr(); - p_uv2xy < - uv2xy[widx].GetDataPtr() + uv2xy[widx].NumElements(); - p_uv2xy += 2, ++p_t_hit) { - if (*p_t_hit < 1 - EPS) *p_uv2xy = *(p_uv2xy + 1) = -1.f; - } - core::Tensor uv2xy2 = - uv2xy[widx].Permute({2, 0, 1}).Contiguous(); // {2, ts, ts} - if (DEBUG) { - io::WriteImage("uv2x_" + i_str + ".png", - Image((uv2xy2[0].To(core::UInt16)))); - io::WriteImage("uv2y_" + i_str + ".png", - Image((uv2xy2[1].To(core::UInt16)))); + std::unique_lock albedo_lock{albedo_mutex}; + if (TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION)) { + // Ensure images are blended in order to correctly calculate + // color correction + cv_next_blend_image.wait(albedo_lock, [&i, &next_blend_image]() { + return next_blend_image == i; + }); + } + if (TEST_ENUM_FLAG(BlendingMethod, blending_method, MAX)) { + utility::LogInfo("Blending image {} with method MAX", i_str); + // Select albedo value with max weight + for (auto p_albedo = albedo.GetDataPtr(), + p_this_albedo = this_albedo[widx].GetDataPtr(); + p_albedo < albedo.GetDataPtr() + albedo.NumElements(); + p_albedo += 4, p_this_albedo += 4) { + if (p_albedo[3] < p_this_albedo[3]) { + for (auto k = 0; k < 3; ++k) + p_albedo[k] = this_contrast[k] * p_this_albedo[k] + + this_brightness[k]; + p_albedo[3] = p_this_albedo[3]; + } } - // albedo[u,v] = image[ i[u,v], j[u,v] ] - // ij[u,v] ? = identity[ uv[i,j] ] - - // C. Interpolate weighted image to weighted texture - this_albedo[widx].Fill(0.f); - ipp::Remap(weighted_image[widx], /*{height, width, 4} f32*/ - uv2xy2[0], /* {texsz, texsz} f32*/ - uv2xy2[1], /* {texsz, texsz} f32*/ - this_albedo[widx], /*{texsz, texsz, 4} f32*/ - t::geometry::Image::InterpType::Linear); - // Weights can become negative with higher order interpolation - if (DEBUG) { - io::WriteImage("this_albedo_" + i_str + ".png", - Image((this_albedo[widx].Slice(2, 0, 3, 1) * 255) - .To(core::UInt8))); - float wtmin = this_albedo[widx] - .Slice(2, 3, 4, 1) - .Min({0, 1, 2}) - .Item(), - wtmax = this_albedo[widx] - .Slice(2, 3, 4, 1) - .Max({0, 1, 2}) - .Item(); - io::WriteImage("image_weights_" + i_str + ".png", - Image(weighted_image[widx] - .Slice(2, 3, 4, 1) - .Contiguous()) - .To(core::UInt8, /*copy=*/false, - /*scale=*/255.f / (wtmax - wtmin), - /*offset=*/-wtmin * 255.f / - (wtmax - wtmin))); - io::WriteImage( - "this_albedo_weight_" + i_str + ".png", - Image(this_albedo[widx].Slice(2, 3, 4, 1).Contiguous()) - .To(core::UInt8, /*copy=*/false, - /*scale=*/255.f / (wtmax - wtmin), - /*offset=*/-wtmin * 255.f / - (wtmax - wtmin))); + } else if (TEST_ENUM_FLAG(BlendingMethod, blending_method, AVERAGE)) { + utility::LogInfo("Blending image {} with method AVERAGE", i_str); + for (auto p_albedo = albedo.GetDataPtr(), + p_this_albedo = this_albedo[widx].GetDataPtr(); + p_albedo < albedo.GetDataPtr() + albedo.NumElements(); + p_albedo += 4, p_this_albedo += 4) { + for (auto k = 0; k < 3; ++k) + p_albedo[k] += (this_contrast[k] * p_this_albedo[k] + + this_brightness[k]) * + p_this_albedo[3]; + p_albedo[3] += p_this_albedo[3]; } - { - std::lock_guard albedo_lock{albedo_mutex}; - if (blending_method == BlendingMethod::MAX) { - // Select albedo value with max weight - for (auto p_albedo = albedo.GetDataPtr(), - p_this_albedo = - this_albedo[widx].GetDataPtr(); - p_albedo < - albedo.GetDataPtr() + albedo.NumElements(); - p_albedo += 4, p_this_albedo += 4) { - if (p_albedo[3] < p_this_albedo[3]) { - for (auto k = 0; k < 4; ++k) - p_albedo[k] = p_this_albedo[k]; - } - } - } else if (blending_method == BlendingMethod::AVERAGE) { - for (auto p_albedo = albedo.GetDataPtr(), - p_this_albedo = - this_albedo[widx].GetDataPtr(); - p_albedo < - albedo.GetDataPtr() + albedo.NumElements(); - p_albedo += 4, p_this_albedo += 4) { - for (auto k = 0; k < 3; ++k) - p_albedo[k] += p_this_albedo[k] * p_this_albedo[3]; - p_albedo[3] += p_this_albedo[3]; - } - } + } + if (TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION)) { + cv_next_blend_image.notify_all(); + if (next_blend_image + max_workers < images.size()) { + feeder.add(next_blend_image + max_workers); } + ++next_blend_image; } }; - tbb::parallel_for(tbb::blocked_range(0, images.size()), - project_one_image); - if (blending_method == BlendingMethod::AVERAGE) { + + size_t n_init_images = + TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION) + ? std::min(max_workers, images.size()) + : images.size(); + std::vector range{0, n_init_images}; + std::iota(range.begin(), range.end(), 0); + tbb::parallel_for_each(range, project_one_image); + if (TEST_ENUM_FLAG(BlendingMethod, blending_method, AVERAGE)) { albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); } diff --git a/cpp/open3d/t/geometry/TriangleMesh.h b/cpp/open3d/t/geometry/TriangleMesh.h index a8402bbd4b1..735608e96c6 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.h +++ b/cpp/open3d/t/geometry/TriangleMesh.h @@ -8,6 +8,7 @@ #pragma once #include +#include #include #include "open3d/core/Tensor.h" @@ -27,7 +28,23 @@ class LineSet; class RaycastingScene; /// Texture Blending method for ProjectImagesToAlbedo() from overlapping images. -enum class BlendingMethod { MAX, AVERAGE }; +enum class BlendingMethod : uint8_t { + /// For each texel, pick the input pixel with the max weight from all + /// overlapping images. This creates sharp textures but may have visible + /// seams. + MAX = 1 << 0, + /// The output texel value is the weighted sum of input pixels. This + /// creates smooth blending without seams, but the results may be blurry. + AVERAGE = 1 << 1, + /// Try to match white balance and exposure settings for the images with a + /// linear 3 channel contrast + brightness correction for each image. Also + /// detects and masks specular highlights and saturated regions. This can + /// be combined with either the MAX or AVERAGE blending methods. + COLOR_CORRECTION = 1 << 2 +}; +#define TEST_ENUM_FLAG(ENUM, VALUE, FLAG) \ + (static_cast>(VALUE) & \ + static_cast>(ENUM::FLAG)) /// \class TriangleMesh /// \brief A triangle mesh contains vertices and triangles. @@ -1006,6 +1023,11 @@ class TriangleMesh : public Geometry, public DrawableGeometry { /// - AVERAGE: The output texel value is the weighted sum of input /// pixels. This creates smooth blending without seams, but the results /// may be blurry. + /// - COLOR_CORRECTION: Try to match white balance and exposure + /// settings for the images with a linear 3 channel contrast + + /// brightness correction for each image. Also detects and masks + /// specular highlights and saturated regions. This can be combined with + /// either the MAX or AVERAGE blending methods. /// \return Image with albedo texture Image ProjectImagesToAlbedo( const std::vector &images, diff --git a/cpp/pybind/t/geometry/trianglemesh.cpp b/cpp/pybind/t/geometry/trianglemesh.cpp index 1c741e4f025..8a69563864a 100644 --- a/cpp/pybind/t/geometry/trianglemesh.cpp +++ b/cpp/pybind/t/geometry/trianglemesh.cpp @@ -993,9 +993,21 @@ or has a negative value, it is ignored. &TriangleMesh::RemoveUnreferencedVertices, "Removes unreferenced vertices from the mesh in-place."); - py::enum_(m, "BlendingMethod") - .value("MAX", BlendingMethod::MAX) - .value("AVERAGE", BlendingMethod::AVERAGE); + py::enum_(m, "BlendingMethod", py::arithmetic()) + .value("MAX", BlendingMethod::MAX, + "For each texel, pick the input pixel with the max weight " + "from all overlapping images. This creates sharp textures " + "but may have visible seams.") + .value("AVERAGE", BlendingMethod::AVERAGE, + "The output texel value is the weighted sum of input " + "pixels. This creates smooth blending without seams, but " + "the results may be blurry.") + .value("COLOR_CORRECTION", BlendingMethod::COLOR_CORRECTION, + "Try to match white balance and exposure settings for the " + "images with a linear 3 channel contrast + brightness " + "correction for each image. Also detects and masks specular " + "highlights and saturated regions. This can be combined " + "with either the `MAX` or `AVERAGE` blending methods."); triangle_mesh.def("project_images_to_albedo", &TriangleMesh::ProjectImagesToAlbedo, "images"_a, "intrinsic_matrices"_a, "extrinsic_matrices"_a, @@ -1028,6 +1040,11 @@ texture. - `AVERAGE`: The output texel value is the weighted sum of input pixels. This creates smooth blending without seams, but the results may be blurry. + - `COLOR_CORRECTION`: Try to match white balance and exposure + settings for the images with a linear 3 channel contrast + + brightness correction for each image. Also detects and masks + specular highlights and saturated regions. This can be combined with + either the MAX or AVERAGE blending methods. Returns: Image with albedo texture.)"); diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index f2640d67cdf..4c39dea8da8 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -1352,6 +1352,8 @@ TEST_P(TriangleMeshPermuteDevices, RemoveUnreferencedVertices) { TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { using namespace t::geometry; + using ::testing::ElementsAre; + using ::testing::FloatEq; core::Device device = GetParam(); TriangleMesh sphere = TriangleMesh::FromLegacy(*geometry::TriangleMesh::CreateSphere( @@ -1390,12 +1392,12 @@ TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { EXPECT_TRUE(sphere.GetMaterial().HasAlbedoMap()); EXPECT_TRUE(albedo.AsTensor().GetShape().IsCompatible({256, 256, 3})); EXPECT_TRUE(albedo.GetDtype() == core::UInt8); - core::Tensor mean_color_ref = - core::Tensor::Init({92.465515, 71.62926, 67.55928}); - EXPECT_TRUE(albedo.AsTensor() + EXPECT_THAT(albedo.AsTensor() .To(core::Float32) .Mean({0, 1}) - .AllClose(mean_color_ref)); + .ToFlatVector(), + ElementsAre(FloatEq(92.465515), FloatEq(71.62926), + FloatEq(67.55928))); } } // namespace tests diff --git a/examples/python/geometry/triangle_mesh_project_to_albedo.py b/examples/python/geometry/triangle_mesh_project_to_albedo.py index 0f057da3c8c..92b40651fde 100644 --- a/examples/python/geometry/triangle_mesh_project_to_albedo.py +++ b/examples/python/geometry/triangle_mesh_project_to_albedo.py @@ -38,8 +38,9 @@ def show_progress(block_num, block_size, total_size): def create_dataset(meshfile, n_images=10, movie=False): - """Render images of a 3D mesh from different viewpoints. These form a - synthetic dataset to test the project_images_to_albedo function. + """Render images of a 3D mesh from different viewpoints, covering the + northern hemisphere. These form a synthetic dataset to test the + project_images_to_albedo function. """ # Adjust these parameters to properly frame your model. # Window system pixel scaling (e.g. 1 for normal, 2 for HiDPI / retina display) @@ -63,12 +64,25 @@ def rotate_camera_and_shoot(o3dvis): images = [] o3dvis.scene.scene.enable_sun_light(False) print("Rendering images: ", end='', flush=True) + n_0 = 2 * n_images // 3 + n_1 = n_images - n_0 - 1 for n in range(n_images): - theta = n * (2 * np.pi) / n_images Rt = np.eye(4) Rt[:3, 3] = t - Rt[:3, :3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_zyx( - [np.pi, theta, 0]) + if n < n_0: + theta = n * (2 * np.pi) / n_0 + Rt[:3, : + 3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_zyx( + [np.pi, theta, 0]) + elif n < n_images - 1: + theta = (n - n_0) * (2 * np.pi) / n_1 + Rt[:3, : + 3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_xyz( + [np.pi / 4, theta, np.pi]) + else: # one image from the top + Rt[:3, : + 3] = o3d.geometry.Geometry3D.get_rotation_matrix_from_zyx( + [np.pi, 0, -np.pi / 2]) Rts.append(Rt) o3dvis.setup_camera(K, Rt, width, height) o3dvis.post_redraw() @@ -90,6 +104,7 @@ def rotate_camera_and_shoot(o3dvis): "glob", "-i", "render-*.jpg", "-y", meshfile.stem + ".mp4" ], check=True) + o3dvis.close() print("\nDone.") print("If the object is properly framed in the GUI window, click on the " @@ -105,7 +120,7 @@ def rotate_camera_and_shoot(o3dvis): actions=[("Save Images", rotate_camera_and_shoot)]) -def albedo_from_images(meshfile, calib_data_file): +def albedo_from_images(meshfile, calib_data_file, albedo_contrast=1.25): model = o3d.io.read_triangle_model(meshfile) tmeshes = o3d.t.geometry.TriangleMesh.from_triangle_mesh_model(model) @@ -116,7 +131,10 @@ def albedo_from_images(meshfile, calib_data_file): images = list(o3d.t.io.read_image(imfile) for imfile in calib["images"]) calib.close() start = time.time() - albedo = tmeshes[0].project_images_to_albedo(images, Ks, Rts, 1024) + albedo = tmeshes[0].project_images_to_albedo( + images, Ks, Rts, 1024, True, o3d.t.geometry.BlendingMethod.MAX) + albedo = albedo.linear_transform(scale=albedo_contrast) # brighten albedo + tmeshes[0].material.texture_maps["albedo"] = albedo print(f"project_images_to_albedo ran in {time.time()-start:.2f}s") o3d.t.io.write_image("albedo.png", albedo) o3d.t.io.write_triangle_mesh(meshfile.stem + "_albedo.glb", tmeshes[0]) @@ -165,6 +183,8 @@ def albedo_from_images(meshfile, calib_data_file): if args.meshfile == Path("."): parser.error("Please provide a path to a mesh file, or use " "--download_sample_model.") + if args.n_images < 10: + parser.error("Atleast 10 images should be used!") create_dataset(args.meshfile, n_images=args.n_images, movie=args.movie) else: albedo_from_images(args.meshfile, "cameras.npz") diff --git a/examples/python/visualization/remove_geometry.py b/examples/python/visualization/remove_geometry.py index 2353462174f..81a3245ffbb 100644 --- a/examples/python/visualization/remove_geometry.py +++ b/examples/python/visualization/remove_geometry.py @@ -38,7 +38,7 @@ def visualize_non_blocking(vis, pcds): curr_sec = int(time.time() - start_time) prev_sec = curr_sec - 1 -while True: +while curr_sec < 10: curr_sec = int(time.time() - start_time) if curr_sec - prev_sec == 1: prev_sec = curr_sec @@ -54,3 +54,4 @@ def visualize_non_blocking(vis, pcds): print("Removing %d" % i) visualize_non_blocking(vis, pcds) + time.sleep(0.025) # yield CPU to others while maintaining responsiveness From 54088bbf27329ebda1d621b479996e0fe8c9d60a Mon Sep 17 00:00:00 2001 From: Sameer Sheorey Date: Tue, 4 Jun 2024 08:13:37 -0700 Subject: [PATCH 8/8] COLOR CORRECTION WIP2 Mostly works - slight issues Image::Resize for Bool has bus error / seg fault. highlights detection doesn't work. --- cpp/open3d/t/geometry/Image.cpp | 6 + cpp/open3d/t/geometry/TriangleMesh.cpp | 142 +++++++++++------- cpp/open3d/t/geometry/TriangleMesh.h | 5 + cpp/pybind/t/geometry/trianglemesh.cpp | 6 +- cpp/tests/t/geometry/TriangleMesh.cpp | 2 +- .../triangle_mesh_project_to_albedo.py | 19 ++- 6 files changed, 118 insertions(+), 62 deletions(-) diff --git a/cpp/open3d/t/geometry/Image.cpp b/cpp/open3d/t/geometry/Image.cpp index 3af4fc20b42..84cc0ec6f98 100644 --- a/cpp/open3d/t/geometry/Image.cpp +++ b/cpp/open3d/t/geometry/Image.cpp @@ -156,6 +156,12 @@ Image Image::Resize(float sampling_rate, InterpType interp_type) const { if (sampling_rate == 1.0f) { return *this; } + if (GetDtype() == core::Bool) { // Resize via UInt8 + return Image(Image(data_.ReinterpretCast(core::Bool)) + .Resize(sampling_rate, interp_type) + .AsTensor() + .ReinterpretCast(core::UInt8)); + } static const dtype_channels_pairs npp_supported{ {core::UInt8, 1}, {core::UInt16, 1}, {core::Float32, 1}, diff --git a/cpp/open3d/t/geometry/TriangleMesh.cpp b/cpp/open3d/t/geometry/TriangleMesh.cpp index 0b583c85861..4d71f9a54ff 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.cpp +++ b/cpp/open3d/t/geometry/TriangleMesh.cpp @@ -1387,7 +1387,10 @@ core::Tensor Project(const core::Tensor &t_xyz, // contiguous {...,3} } /// Estimate contrast and brightness in each color channel for this_albedo to -/// match albedo, based on the overlapping area in the texture image. +/// match albedo, based on the overlapping area in the texture image. Contrast +/// and brightness are estimated by matching the second and first moments +/// (variance and mean) of the pixel colors, respectively. +/// albedo and this_albedo are float images with range [0,1] std::tuple, std::array, core::Tensor, @@ -1395,48 +1398,59 @@ std::tuple, get_color_correction(const core::Tensor &albedo, const core::Tensor &this_albedo, float SPECULAR_THRESHOLD) { - constexpr double EPS = 1e-6; - std::array this_contrast{}, this_brightness{}, sum{}, this_sum{}, - absdiffsum{EPS, EPS, EPS}, this_absdiffsum{EPS, EPS, EPS}; - auto q_albedo = - Image(albedo).Resize(0.25, Image::InterpType::Linear).AsTensor(); + const double EPS = 1e-6; + const float shift = 0.5f; // compute shifted sum / sumsqr for stability. + const unsigned MIN_OVERLAP_PIXELS = 32 * 32; + // Perform the color correction with albedo downsampled to size 256x256. + float resize_down = 256.f / albedo.GetShape(0) + /*,resize_up = albedo.GetShape(0) / 256.f*/; + std::array this_contrast{1.f, 1.f, 1.f}, + this_brightness{0.f, 0.f, 0.f}, sum{}, this_sum{}, sumsqr{}, + this_sumsqr{}; + auto q_albedo = Image(albedo) + .Resize(resize_down, Image::InterpType::Linear) + .AsTensor(); auto q_this_albedo = Image(this_albedo) - .Resize(0.25, Image::InterpType::Linear) + .Resize(resize_down, Image::InterpType::Linear) .AsTensor(); - int64_t rowstride = albedo.GetStride(0), row = 0, col = 0; + unsigned count = 0; for (float *pq_albedo = q_albedo.GetDataPtr(), *pq_this_albedo = q_this_albedo.GetDataPtr(); pq_albedo < q_albedo.GetDataPtr() + q_albedo.NumElements(); - pq_albedo += 4, pq_this_albedo += 4, ++col) { - if (col == albedo.GetShape(1)) { - col = 0; - ++row; - } + pq_albedo += 4, pq_this_albedo += 4) { if (pq_albedo[3] <= EPS || pq_this_albedo[3] <= EPS) continue; - bool update_absdiffsum = row > 0 && row < albedo.GetShape(0) - 1 && - col > 0 && col < albedo.GetShape(1) - 1; + ++count; for (int c = 0; c < 3; ++c) { - sum[c] += pq_albedo[c]; - this_sum[c] += pq_this_albedo[c]; - // Absolute central difference - if (update_absdiffsum) { - absdiffsum[c] += std::fabs( - pq_albedo[c] - 0.25 * (pq_albedo[c - 4] + - pq_albedo[c - 4 - rowstride] + - pq_albedo[c + 4] + - pq_albedo[c + rowstride + 4])); - this_absdiffsum[c] += - std::fabs(pq_this_albedo[c] - - 0.25 * (pq_this_albedo[c - 4] + - pq_this_albedo[c - 4 - rowstride] + - pq_this_albedo[c + 4] + - pq_this_albedo[c + rowstride + 4])); - } + sum[c] += pq_albedo[c] - shift; + sumsqr[c] += (pq_albedo[c] - shift) * (pq_albedo[c] - shift); + this_sum[c] += pq_this_albedo[c] - shift; + this_sumsqr[c] += + (pq_this_albedo[c] - shift) * (pq_this_albedo[c] - shift); } } + if (count <= MIN_OVERLAP_PIXELS) { + utility::LogWarning( + "[ProjectImagesToAlbedo] Too few overlapping pixels ({}/{}) " + "found for color correction.", + count, MIN_OVERLAP_PIXELS); + return std::make_tuple( + this_contrast, this_brightness, + core::Tensor::Zeros({albedo.GetShape(0), albedo.GetShape(1), 1}, + core::Bool), + core::Tensor::Zeros({albedo.GetShape(0), albedo.GetShape(1), 1}, + core::Bool)); + } for (int c = 0; c < 3; ++c) { - this_contrast[c] = absdiffsum[c] / this_absdiffsum[c]; - this_brightness[c] = sum[c] - this_contrast[c] * this_sum[c]; + float variance = (sumsqr[c] - sum[c] * sum[c] / count) / count, + this_variance = + (this_sumsqr[c] - this_sum[c] * this_sum[c] / count) / + count; + utility::LogWarning("count: {}, variance: {}, this_variance: {}", count, + variance, this_variance); + this_contrast[c] = sqrt((variance + EPS) / (this_variance + EPS)); + sum[c] += count * shift; // get un-shifted sum for brightness. + this_sum[c] += count * shift; + this_brightness[c] = (sum[c] - this_contrast[c] * this_sum[c]) / count; } // Detect specular highlights from outliers. @@ -1465,14 +1479,24 @@ get_color_correction(const core::Tensor &albedo, *p_this_highlights_mask = std::max({err[0], err[1], err[2]}) < -SPECULAR_THRESHOLD; } - highlights_mask = Image(highlights_mask) - .Dilate(/*kernel_size=*/3) - .Resize(4.f) - .AsTensor(); - this_highlights_mask = Image(this_highlights_mask) - .Dilate(/*kernel_size=*/3) - .Resize(4.f) - .AsTensor(); + utility::LogInfo( + "n_highlights: {}, n_this_highlights: {}", + highlights_mask.To(core::UInt32).Sum({0, 1, 2}).Item(), + this_highlights_mask.To(core::UInt32) + .Sum({0, 1, 2}) + .Item()); + // highlights_mask = Image(highlights_mask) + // .Dilate(/*kernel_size=*/3) + // .Resize(resize_up) + // .AsTensor(); + // this_highlights_mask = Image(this_highlights_mask) + // .Dilate(/*kernel_size=*/3) + // .Resize(resize_up) + // .AsTensor(); + highlights_mask = core::Tensor::Zeros( + {albedo.GetShape(0), albedo.GetShape(1), 1}, core::Bool), + this_highlights_mask = core::Tensor::Zeros( + {albedo.GetShape(0), albedo.GetShape(1), 1}, core::Bool); return std::make_tuple(this_contrast, this_brightness, highlights_mask, this_highlights_mask); } @@ -1494,6 +1518,10 @@ Image TriangleMesh::ProjectImagesToAlbedo( if (!triangle_attr_.Contains("texture_uvs")) { utility::LogError("Cannot find triangle attribute 'texture_uvs'"); } + if (!TEST_ENUM_FLAG(BlendingMethod, blending_method, MAX) && + !TEST_ENUM_FLAG(BlendingMethod, blending_method, AVERAGE)) { + utility::LogError("Select one of MAX and AVERAGE BlendingMethod s."); + } core::Tensor texture_uvs = triangle_attr_.at("texture_uvs").To(core::Device()).Contiguous(); core::AssertTensorShape(texture_uvs, {core::None, 3, 2}); @@ -1538,8 +1566,10 @@ Image TriangleMesh::ProjectImagesToAlbedo( // ensures correct computation of color matching (only neded for // COLOR_CORRECTION). std::condition_variable cv_next_blend_image; - std::atomic_uint next_blend_image{0}; - auto project_one_image = [&](size_t i, tbb::feeder &feeder) { + size_t next_blend_image{0}; + /* auto project_one_image = [&](size_t i, tbb::feeder &feeder) { */ + auto project_one_image = [&](size_t i, + tbb::parallel_do_feeder &feeder) { size_t widx = tbb::this_task_arena::current_thread_index(); // initialize thread variables if (!this_albedo[widx].GetShape().IsCompatible( @@ -1633,10 +1663,9 @@ Image TriangleMesh::ProjectImagesToAlbedo( io::WriteImage("uv2y_" + i_str + ".png", Image((uv2xy2[1].To(core::UInt16)))); } - // albedo[u,v] = image[ i[u,v], j[u,v] ] - // ij[u,v] ? = identity[ uv[i,j] ] // C. Interpolate weighted image to weighted texture + // albedo[u,v] = image[ i[u,v], j[u,v] ] this_albedo[widx].Fill(0.f); ipp::Remap(weighted_image[widx], /*{height, width, 4} f32*/ uv2xy2[0], /* {texsz, texsz} f32*/ @@ -1672,16 +1701,6 @@ Image TriangleMesh::ProjectImagesToAlbedo( std::array this_contrast{1.f, 1.f, 1.f}, this_brightness{0.f, 0.f, 0.f}; core::Tensor highlights_mask, this_highlights_mask; - if (TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION)) { - std::tie(this_contrast, this_brightness, highlights_mask, - this_highlights_mask) = - get_color_correction(albedo, this_albedo[widx], - SPECULAR_THRESHOLD); - utility::LogDebug( - "[ProjectImagesToAlbedo] Image {}, contrast: {}, " - "brightness {}", - i_str, this_contrast, this_brightness); - } std::unique_lock albedo_lock{albedo_mutex}; if (TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION)) { @@ -1690,6 +1709,14 @@ Image TriangleMesh::ProjectImagesToAlbedo( cv_next_blend_image.wait(albedo_lock, [&i, &next_blend_image]() { return next_blend_image == i; }); + std::tie(this_contrast, this_brightness, highlights_mask, + this_highlights_mask) = + get_color_correction(albedo, this_albedo[widx], + SPECULAR_THRESHOLD); + utility::LogInfo( + "[ProjectImagesToAlbedo] Image {}, contrast: {}, " + "brightness {}", + i_str, this_contrast, this_brightness); } if (TEST_ENUM_FLAG(BlendingMethod, blending_method, MAX)) { utility::LogInfo("Blending image {} with method MAX", i_str); @@ -1731,9 +1758,10 @@ Image TriangleMesh::ProjectImagesToAlbedo( TEST_ENUM_FLAG(BlendingMethod, blending_method, COLOR_CORRECTION) ? std::min(max_workers, images.size()) : images.size(); - std::vector range{0, n_init_images}; + std::vector range(n_init_images, 0); std::iota(range.begin(), range.end(), 0); - tbb::parallel_for_each(range, project_one_image); + /* tbb::parallel_for_each(range, project_one_image); */ + tbb::parallel_do(range.begin(), range.end(), project_one_image); if (TEST_ENUM_FLAG(BlendingMethod, blending_method, AVERAGE)) { albedo.Slice(2, 0, 3, 1) /= albedo.Slice(2, 3, 4, 1); } diff --git a/cpp/open3d/t/geometry/TriangleMesh.h b/cpp/open3d/t/geometry/TriangleMesh.h index 735608e96c6..eeabc3c524f 100644 --- a/cpp/open3d/t/geometry/TriangleMesh.h +++ b/cpp/open3d/t/geometry/TriangleMesh.h @@ -42,6 +42,11 @@ enum class BlendingMethod : uint8_t { /// be combined with either the MAX or AVERAGE blending methods. COLOR_CORRECTION = 1 << 2 }; +inline BlendingMethod operator|(BlendingMethod a, BlendingMethod b) { + return static_cast( + static_cast>(a) | + static_cast>(b)); +} #define TEST_ENUM_FLAG(ENUM, VALUE, FLAG) \ (static_cast>(VALUE) & \ static_cast>(ENUM::FLAG)) diff --git a/cpp/pybind/t/geometry/trianglemesh.cpp b/cpp/pybind/t/geometry/trianglemesh.cpp index 8a69563864a..575e5c43287 100644 --- a/cpp/pybind/t/geometry/trianglemesh.cpp +++ b/cpp/pybind/t/geometry/trianglemesh.cpp @@ -7,6 +7,8 @@ #include "open3d/t/geometry/TriangleMesh.h" +#include + #include #include @@ -1007,7 +1009,9 @@ or has a negative value, it is ignored. "images with a linear 3 channel contrast + brightness " "correction for each image. Also detects and masks specular " "highlights and saturated regions. This can be combined " - "with either the `MAX` or `AVERAGE` blending methods."); + "with either the `MAX` or `AVERAGE` blending methods.") + .def("__or__", + py::overload_cast(&operator|)); triangle_mesh.def("project_images_to_albedo", &TriangleMesh::ProjectImagesToAlbedo, "images"_a, "intrinsic_matrices"_a, "extrinsic_matrices"_a, diff --git a/cpp/tests/t/geometry/TriangleMesh.cpp b/cpp/tests/t/geometry/TriangleMesh.cpp index 4c39dea8da8..72ecb6f9631 100644 --- a/cpp/tests/t/geometry/TriangleMesh.cpp +++ b/cpp/tests/t/geometry/TriangleMesh.cpp @@ -1386,7 +1386,7 @@ TEST_P(TriangleMeshPermuteDevices, ProjectImagesToAlbedo) { {Image(view[0]), Image(view[1]), Image(view[2])}, {intrinsic_matrix, intrinsic_matrix, intrinsic_matrix}, {extrinsic_matrix[0], extrinsic_matrix[1], extrinsic_matrix[2]}, - 256, true); + 256, true, BlendingMethod::MAX | BlendingMethod::COLOR_CORRECTION); EXPECT_TRUE(sphere.HasMaterial()); EXPECT_TRUE(sphere.GetMaterial().HasAlbedoMap()); diff --git a/examples/python/geometry/triangle_mesh_project_to_albedo.py b/examples/python/geometry/triangle_mesh_project_to_albedo.py index 92b40651fde..390c086c55a 100644 --- a/examples/python/geometry/triangle_mesh_project_to_albedo.py +++ b/examples/python/geometry/triangle_mesh_project_to_albedo.py @@ -37,7 +37,7 @@ def show_progress(block_num, block_size, total_size): print("\nDownload complete.") -def create_dataset(meshfile, n_images=10, movie=False): +def create_dataset(meshfile, n_images=10, movie=False, vary_exposure=False): """Render images of a 3D mesh from different viewpoints, covering the northern hemisphere. These form a synthetic dataset to test the project_images_to_albedo function. @@ -59,6 +59,10 @@ def create_dataset(meshfile, n_images=10, movie=False): unlit = rendering.MaterialRecord() unlit.shader = "unlit" + def triangle_wave(n, period=1): + """Triangle wave function between [0,1] with given period.""" + return abs(n % period - period / 2) / (period / 2) + def rotate_camera_and_shoot(o3dvis): Rts = [] images = [] @@ -85,6 +89,11 @@ def rotate_camera_and_shoot(o3dvis): [np.pi, 0, -np.pi / 2]) Rts.append(Rt) o3dvis.setup_camera(K, Rt, width, height) + # Vary IBL intensity as a poxy for exposure value. IBL ranges from + # [0,150000]. We vary it between 20000 and 100000. + if vary_exposure: + o3dvis.set_ibl_intensity(20000 + + 80000 * triangle_wave(n, n_images / 4)) o3dvis.post_redraw() o3dvis.export_current_image(f"render-{n:02}.jpg") images.append(f"render-{n:02}.jpg") @@ -132,7 +141,8 @@ def albedo_from_images(meshfile, calib_data_file, albedo_contrast=1.25): calib.close() start = time.time() albedo = tmeshes[0].project_images_to_albedo( - images, Ks, Rts, 1024, True, o3d.t.geometry.BlendingMethod.MAX) + images, Ks, Rts, 1024, True, o3d.t.geometry.BlendingMethod.MAX | + o3d.t.geometry.BlendingMethod.COLOR_CORRECTION) albedo = albedo.linear_transform(scale=albedo_contrast) # brighten albedo tmeshes[0].material.texture_maps["albedo"] = albedo print(f"project_images_to_albedo ran in {time.time()-start:.2f}s") @@ -185,6 +195,9 @@ def albedo_from_images(meshfile, calib_data_file, albedo_contrast=1.25): "--download_sample_model.") if args.n_images < 10: parser.error("Atleast 10 images should be used!") - create_dataset(args.meshfile, n_images=args.n_images, movie=args.movie) + create_dataset(args.meshfile, + n_images=args.n_images, + movie=args.movie, + vary_exposure=True) else: albedo_from_images(args.meshfile, "cameras.npz")