Port point-to-point-cloud distance algorithm to CUDA (#4230)

* WIP: Port point-to-point cloud algorithms to CUDA

* Fix

* Fix

* Add heap bytes functions

* Add comment

* Add option to skip point with same index

* WIP: Add generalized interface

* WIP: Change API

* WIP: Add progress callback

* Add comments

* Fix build

* WIP: Add unit tests for MRCuda

* Fix error

* Fix tests

* Increase verbosity

* Refactor

* Fix GCC build

* Fix GCC build

* Fix VS project

* Refactor

* Fix bindings build

* Fix MSVC build

* Fix build

* Fix build

* Fix MSVC build

* Fix VS project

Author: oitel

CMakeLists.txt

@@ -314,6 +314,9 @@ ENDIF()
 IF(BUILD_TESTING)
   enable_testing()
   add_subdirectory(${PROJECT_SOURCE_DIR}/MRTest ./MRTest)
+  IF(MESHLIB_BUILD_MRCUDA)
+    add_subdirectory(${PROJECT_SOURCE_DIR}/MRTestCuda ./MRTestCuda)
+  ENDIF()
   IF(MESHLIB_EXPERIMENTAL_BUILD_C_BINDING)
     add_subdirectory(${PROJECT_SOURCE_DIR}/MRTestC ./MRTestC)
   ENDIF()

scripts/check_vcxproj.py

@@ -36,10 +36,18 @@ def find_missing_entries(vcxproj_path):
             VCXPROJ_NAMESPACES,
         )
     }
+    cuda_compiles = {
+        item.attrib['Include'].lower()
+        for item in project.iterfind(
+            'msbuild:ItemGroup/msbuild:CudaCompile',
+            VCXPROJ_NAMESPACES,
+        )
+    }
 
     result = {
         'ClInclude': [],
         'ClCompile': [],
+        'CudaCompile': [],
     }
     for path in vcxproj_dir.iterdir():
         name, suffix = path.name.lower(), path.suffix.lower()
@@ -48,22 +56,23 @@ def find_missing_entries(vcxproj_path):
         if suffix in {".cpp"}:
             if name not in compiles:
                 result['ClCompile'].append(path)
-        elif suffix in {".h", ".hpp"}:
+        elif suffix in {".h", ".hpp", ".cuh"}:
             if name not in includes:
                 result['ClInclude'].append(path)
+        elif suffix in {".cu"}:
+            if name not in cuda_compiles:
+                result['CudaCompile'].append(path)
 
     return result
 
 
 def process_file(vcxproj_path):
     result = find_missing_entries(vcxproj_path)
     ok = True
-    for path in result['ClInclude']:
-        print(f"{vcxproj_path}: missing ClInclude item: {path.name}", file=sys.stderr)
-        ok = False
-    for path in result['ClCompile']:
-        print(f"{vcxproj_path}: missing ClCompile entry: {path.name}", file=sys.stderr)
-        ok = False
+    for group_name, group in result.items():
+        for path in group:
+            print(f"{vcxproj_path}: missing {group_name} item: {path.name}", file=sys.stderr)
+            ok = False
     return ok
 
 

source/MRCuda/MRCuda.vcxproj

@@ -92,6 +92,7 @@
     <CudaCompile Include="MRCudaSolarRadiation.cu" />
     <CudaCompile Include="MRCudaTest.cu" />
     <CudaCompile Include="MRCudaContoursDistanceMap.cu" />
+    <CudaCompile Include="MRCudaPointsProject.cu" />
   </ItemGroup>
   <ItemGroup>
     <ClInclude Include="exports.h" />
@@ -114,6 +115,11 @@
     <ClInclude Include="MRCudaTest.cuh" />
     <ClInclude Include="MRCudaTest.h" />
     <ClInclude Include="MRCudaContoursDistanceMap.h" />
+    <ClInclude Include="MRCudaMath.h" />
+    <ClInclude Include="MRCudaPointCloud.h" />
+    <ClInclude Include="MRCudaPointCloud.cuh" />
+    <ClInclude Include="MRCudaPointsProject.h" />
+    <ClInclude Include="MRCudaPointsProject.cuh" />
   </ItemGroup>
   <ItemGroup>
     <ClCompile Include="MRCudaBasic.cpp" />
@@ -124,6 +130,9 @@
     <ClCompile Include="MRCudaPointsToMeshProjector.cpp" />
     <ClCompile Include="MRCudaSolarRadiation.cpp" />
     <ClCompile Include="MRCudaTest.cpp" />
+    <ClCompile Include="MRCudaMath.cpp" />
+    <ClCompile Include="MRCudaPointCloud.cpp" />
+    <ClCompile Include="MRCudaPointsProject.cpp" />
   </ItemGroup>
   <ItemGroup>
     <ProjectReference Include="..\MRMesh\MRMesh.vcxproj">

source/MRCuda/MRCudaFloat.cuh

@@ -1,6 +1,10 @@
 #pragma once
+
 #include "cuda_runtime.h"
 
+#include <cassert>
+#include <cmath>
+
 namespace MR
 {
 

source/MRCuda/MRCudaMath.cpp

@@ -0,0 +1,32 @@
+#include "MRCudaMath.h"
+#include "MRCudaMath.cuh"
+
+#include "MRMesh/MRAffineXf3.h"
+
+namespace MR::Cuda
+{
+
+float3 fromVec( const Vector3f& v )
+{
+    return {
+        .x = v.x,
+        .y = v.y,
+        .z = v.z,
+    };
+}
+
+Matrix4 fromXf( const MR::AffineXf3f& xf )
+{
+    if ( xf == AffineXf3f{} )
+        return Matrix4 { .isIdentity = true };
+
+    return {
+        .x = fromVec( xf.A.x ),
+        .y = fromVec( xf.A.y ),
+        .z = fromVec( xf.A.z ),
+        .b = fromVec( xf.b ),
+        .isIdentity = false,
+    };
+}
+
+} // namespace MR::Cuda

source/MRCuda/MRCudaMath.h

@@ -0,0 +1,21 @@
+#pragma once
+#ifndef MR_PARSING_FOR_PB11_BINDINGS
+
+#include "exports.h"
+
+#include "MRMesh/MRMeshFwd.h"
+
+struct float3;
+
+namespace MR::Cuda
+{
+
+// structs from MRCudaMath.cuh
+struct Matrix4;
+
+// copy from CPU to GPU structs
+MRCUDA_API float3 fromVec( const Vector3f& v );
+MRCUDA_API Matrix4 fromXf( const MR::AffineXf3f& xf );
+
+} // namespace MR::Cuda
+#endif

source/MRCuda/MRCudaPointCloud.cpp

@@ -0,0 +1,40 @@
+#include "MRCudaPointCloud.h"
+#include "MRCudaPointCloud.cuh"
+
+#include "MRMesh/MRAABBTreePoints.h"
+#include "MRMesh/MRPointCloud.h"
+
+namespace MR::Cuda
+{
+
+Expected<std::unique_ptr<PointCloudDataHolder>> copyDataFrom( const PointCloud& pc, bool copyNormals,
+    const std::vector<Vector3f>* normals )
+{
+    const auto& tree = pc.getAABBTree();
+    const auto& nodes = tree.nodes();
+    const auto& points = tree.orderedPoints();
+
+    auto result = std::make_unique<PointCloudDataHolder>();
+
+    CUDA_LOGE_RETURN_UNEXPECTED( result->nodes.fromVector( nodes.vec_ ) );
+    CUDA_LOGE_RETURN_UNEXPECTED( result->points.fromVector( points ) );
+    if ( copyNormals )
+        CUDA_LOGE_RETURN_UNEXPECTED( result->normals.fromVector( normals ? *normals : pc.normals.vec_ ) );
+
+    return result;
+}
+
+size_t pointCloudHeapBytes( const PointCloud& pc, bool copyNormals, const std::vector<Vector3f>* normals )
+{
+    const auto& tree = pc.getAABBTree();
+    const auto& nodes = tree.nodes();
+    const auto& points = tree.orderedPoints();
+
+    return
+          nodes.size() * sizeof( Node3 )
+        + points.size() * sizeof( OrderedPoint )
+        + ( copyNormals ? ( normals ? normals->size() : pc.normals.size() ) * sizeof( float3 ) : 0 )
+    ;
+}
+
+} // namespace MR::Cuda

source/MRCuda/MRCudaPointCloud.cuh

@@ -0,0 +1,41 @@
+#pragma once
+
+#include "MRCudaBasic.cuh"
+#include "MRCudaMath.cuh"
+
+namespace MR::Cuda
+{
+
+// struct similar to MR::Point
+struct OrderedPoint
+{
+    float3 coord;
+    int id;
+};
+
+// point cloud data required for algorithms
+struct PointCloudData
+{
+    const Node3* __restrict__ nodes;
+    const OrderedPoint* __restrict__ points;
+    const float3* __restrict__ normals;
+};
+
+// GPU memory holder for point cloud data
+struct PointCloudDataHolder
+{
+    DynamicArray<Node3> nodes;
+    DynamicArray<OrderedPoint> points;
+    DynamicArray<float3> normals;
+
+    [[nodiscard]] PointCloudData data() const
+    {
+        return {
+            nodes.data(),
+            points.data(),
+            normals.data(),
+        };
+    }
+};
+
+} // namespace MR::Cuda

source/MRCuda/MRCudaPointCloud.h

@@ -0,0 +1,23 @@
+#pragma once
+#ifndef MR_PARSING_FOR_PB11_BINDINGS
+
+#include "exports.h"
+
+#include "MRMesh/MRMeshFwd.h"
+#include "MRMesh/MRExpected.h"
+
+namespace MR::Cuda
+{
+
+struct PointCloudDataHolder;
+
+/// copy point cloud-related data to the GPU memory
+MRCUDA_API Expected<std::unique_ptr<PointCloudDataHolder>> copyDataFrom( const PointCloud& pc, bool copyNormals = false,
+    const std::vector<Vector3f>* normals = nullptr );
+
+/// return the amount of GPU memory required for \ref MR::Cuda::PointCloudDataHolder
+MRCUDA_API size_t pointCloudHeapBytes( const PointCloud& pc, bool copyNormals = false,
+    const std::vector<Vector3f>* normals = nullptr );
+
+} // namespace MR::Cuda
+#endif

source/MRCuda/MRCudaPointsProject.cpp

@@ -0,0 +1,104 @@
+#include "MRCudaPointsProject.cuh"
+#include "MRCudaPointsProject.h"
+
+#include "MRCudaBasic.cuh"
+#include "MRCudaBasic.h"
+#include "MRCudaMath.h"
+#include "MRCudaPointCloud.h"
+
+#include "MRMesh/MRAffineXf3.h"
+#include "MRMesh/MRBitSet.h"
+#include "MRMesh/MRChunkIterator.h"
+#include "MRMesh/MRProgressCallback.h"
+
+static_assert( sizeof( MR::Cuda::PointsProjectionResult ) == sizeof( MR::PointsProjectionResult ) );
+
+namespace MR::Cuda
+{
+
+Expected<std::vector<MR::PointsProjectionResult>> findProjectionOnPoints( const PointCloud& pointCloud,
+    const std::vector<Vector3f>& points, const FindProjectionOnPointsSettings& settings )
+{
+    std::vector<MR::PointsProjectionResult> results;
+    PointsProjector projector;
+    return projector.setPointCloud( pointCloud )
+        .and_then( [&] { return projector.findProjections( results, points, settings ); } )
+        .transform( [&] { return results; } );
+}
+
+Expected<void> PointsProjector::setPointCloud( const PointCloud& pointCloud )
+{
+    if ( auto res = copyDataFrom( pointCloud ) )
+    {
+        data_ = std::move( *res );
+        return {};
+    }
+    else
+    {
+        return unexpected( std::move( res.error() ) );
+    }
+}
+
+Expected<void> PointsProjector::findProjections( std::vector<MR::PointsProjectionResult>& results,
+    const std::vector<Vector3f>& points, const FindProjectionOnPointsSettings& settings ) const
+{
+    if ( !data_ )
+        return unexpected( "No reference point cloud is set" );
+
+    const auto totalSize = points.size();
+    const auto bufferSize = maxBufferSize( getCudaSafeMemoryLimit(), totalSize, sizeof( float3 ) + sizeof( PointsProjectionResult ) );
+
+    DynamicArray<float3> cudaPoints;
+    CUDA_LOGE_RETURN_UNEXPECTED( cudaPoints.resize( bufferSize ) );
+
+    DynamicArray<PointsProjectionResult> cudaResult;
+    CUDA_LOGE_RETURN_UNEXPECTED( cudaResult.resize( bufferSize ) );
+
+    results.resize( totalSize );
+
+    DynamicArray<uint64_t> cudaValid;
+    if ( settings.valid )
+    {
+        assert( points.size() <= settings.valid->size() );
+        std::vector<uint64_t> validVec;
+        boost::to_block_range( *settings.valid, std::back_inserter( validVec ) );
+        CUDA_LOGE_RETURN_UNEXPECTED( cudaValid.fromVector( validVec ) );
+    }
+
+    const auto cudaXf = settings.xf ? fromXf( *settings.xf ) : Matrix4{};
+
+    if ( !reportProgress( settings.cb, 0.60f ) )
+        return unexpectedOperationCanceled();
+
+    const auto cb1 = subprogress( settings.cb, 0.60f, 1.00f );
+    const auto iterCount = chunkCount( totalSize, bufferSize );
+    size_t iterIndex = 0;
+
+    for ( const auto [offset, size] : splitByChunks( totalSize, bufferSize ) )
+    {
+        const auto cb2 = subprogress( cb1, iterIndex++, iterCount );
+
+        CUDA_LOGE_RETURN_UNEXPECTED( cudaPoints.copyFrom( points.data() + offset, size ) );
+
+        findProjectionOnPointsKernel( cudaResult.data(), data_->data(), cudaPoints.data(), settings.valid ? cudaValid.data() : nullptr, cudaXf, settings.upDistLimitSq, settings.loDistLimitSq, settings.skipSameIndex, size, offset );
+        CUDA_LOGE_RETURN_UNEXPECTED( cudaGetLastError() );
+        if ( !reportProgress( cb2, 0.33f ) )
+            return unexpectedOperationCanceled();
+
+        CUDA_LOGE_RETURN_UNEXPECTED( cudaResult.copyTo( results.data() + offset, size ) );
+        if ( !reportProgress( cb2, 1.00f ) )
+            return unexpectedOperationCanceled();
+    }
+
+    return {};
+}
+
+size_t findProjectionOnPointsHeapBytes( const PointCloud& pointCloud, size_t pointsCount )
+{
+    constexpr size_t cMinCudaBufferSize = 1 << 24; // 16 MiB
+    return
+          pointCloudHeapBytes( pointCloud )
+        + std::min( ( sizeof( float3 ) + sizeof( PointsProjectionResult ) ) * pointsCount, cMinCudaBufferSize );
+}
+
+} // namespace MR::Cuda