From 3bb1f87d539cccbb96c79e50f3d7803b974ebd13 Mon Sep 17 00:00:00 2001 From: Grant Karapetyan Date: Mon, 29 Aug 2022 18:40:03 +0300 Subject: [PATCH 1/3] Smooth voxel mask meshing --- source/MRMesh/MRVolumeSegment.cpp | 70 +++++++++++++++++++++++++++++++ source/MRMesh/MRVolumeSegment.h | 7 ++++ 2 files changed, 77 insertions(+) diff --git a/source/MRMesh/MRVolumeSegment.cpp b/source/MRMesh/MRVolumeSegment.cpp index 517a26b2f8fe..c348967b2562 100644 --- a/source/MRMesh/MRVolumeSegment.cpp +++ b/source/MRMesh/MRVolumeSegment.cpp @@ -11,6 +11,76 @@ namespace MR { +tl::expected meshFromVoxelsMask( const ObjectVoxels& volume, const VoxelBitSet& mask ) +{ + if ( !volume.grid() ) + return tl::make_unexpected( "Cannot create mesh from empty volume." ); + if ( mask.none() ) + return tl::make_unexpected( "Cannot create mesh from empty mask." ); + + const auto& indexer = volume.getVolumeIndexer(); + auto expandedMask = mask; + expandVoxelsMask( expandedMask, indexer, 25 ); // 25 should be in params + + auto accessor = volume.grid()->getConstAccessor(); + double insideAvg = 0.0; + double outsideAvg = 0.0; + Box3i partBox; + for ( auto voxelId : expandedMask ) + { + auto pos = indexer.toPos( voxelId ); + partBox.include( pos ); + if ( mask.test( voxelId ) ) + insideAvg += double( accessor.getValue( { pos.x,pos.y,pos.z } ) ); + else + outsideAvg += double( accessor.getValue( { pos.x,pos.y,pos.z } ) ); + } + insideAvg /= double( mask.count() ); + outsideAvg /= double( expandedMask.count() - mask.count() ); + auto range = float( insideAvg - outsideAvg ); + + SimpleVolume voulmePart; + voulmePart.dims = partBox.size() + Vector3i::diagonal( 1 ); + + auto smallExpMask = mask; + auto smallShrMask = mask; + expandVoxelsMask( smallExpMask, indexer, 3 ); + shrinkVoxelsMask( smallShrMask, indexer, 3 ); + + + const int newDimX = voulmePart.dims.x; + const size_t newDimXY = size_t( newDimX ) * voulmePart.dims.y; + voulmePart.data.resize( newDimXY * voulmePart.dims.z ); + for ( int z = partBox.min.z; z <= partBox.max.z; ++z ) + for ( int y = partBox.min.y; y <= partBox.max.y; ++y ) + for ( int x = partBox.min.x; x <= partBox.max.x; ++x ) + { + auto voxId = indexer.toVoxelId( { x,y,z } ); + if ( smallShrMask.test( voxId ) ) + voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = 1.0f; + else if ( smallExpMask.test( voxId ) ) + voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = + std::clamp( ( accessor.getValue( { x,y,z } ) - float( outsideAvg ) ) / range, 0.0f, 1.0f ); + else + voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = 0.0f; + } + + auto voxelSize = volume.voxelSize(); + auto grid = simpleVolumeToDenseGrid( voulmePart ); + auto mesh = gridToMesh( grid, voxelSize, 0.5f ).value(); // no callback so cannot b stopped + + for ( auto& p : mesh.points ) + { + p = p + mult( Vector3f( partBox.min ), voxelSize ); + } + + if ( mesh.topology.numValidFaces() == 0 ) + return tl::make_unexpected( "Failed to create segmented mesh" ); + + return mesh; + +} + tl::expected segmentVolume( const ObjectVoxels& volume, const std::vector>& pairs, const VolumeSegmentationParameters& params ) { diff --git a/source/MRMesh/MRVolumeSegment.h b/source/MRMesh/MRVolumeSegment.h index d51725103b52..6cefd08ac4fe 100644 --- a/source/MRMesh/MRVolumeSegment.h +++ b/source/MRMesh/MRVolumeSegment.h @@ -18,6 +18,13 @@ class ObjectVoxels; * \{ */ +/** + * \brief Creates mesh from voxels mask + * \param mask in space of whole volume + * density inside mask is expected to be higher then outside + */ +MRMESH_API tl::expected meshFromVoxelsMask( const ObjectVoxels& volume, const VoxelBitSet& mask ); + /** * \brief Parameters for volume segmentation * From ad7bea12cd19c793fd11fbfb4eef309b33fdef95 Mon Sep 17 00:00:00 2001 From: Grant Karapetyan Date: Tue, 30 Aug 2022 14:31:11 +0300 Subject: [PATCH 2/3] Remove duplicated code --- source/MRMesh/MRVolumeSegment.cpp | 173 ++++++++++++++++-------------- 1 file changed, 95 insertions(+), 78 deletions(-) diff --git a/source/MRMesh/MRVolumeSegment.cpp b/source/MRMesh/MRVolumeSegment.cpp index c348967b2562..4ff18eeeefbb 100644 --- a/source/MRMesh/MRVolumeSegment.cpp +++ b/source/MRMesh/MRVolumeSegment.cpp @@ -11,74 +11,119 @@ namespace MR { -tl::expected meshFromVoxelsMask( const ObjectVoxels& volume, const VoxelBitSet& mask ) +// creates mesh from simple volume as 0.5 iso-surface +tl::expected meshFromSimpleVolume( const SimpleVolume& volumePart, const Vector3i& shift ) { - if ( !volume.grid() ) - return tl::make_unexpected( "Cannot create mesh from empty volume." ); - if ( mask.none() ) - return tl::make_unexpected( "Cannot create mesh from empty mask." ); + auto grid = simpleVolumeToDenseGrid( volumePart ); + auto mesh = gridToMesh( grid, volumePart.voxelSize, 0.5f ).value(); // no callback so cannot b stopped - const auto& indexer = volume.getVolumeIndexer(); + auto minCorner = mult( Vector3f( shift ), volumePart.voxelSize ); + for ( auto& p : mesh.points ) + { + p = p + minCorner; + } + + if ( mesh.topology.numValidFaces() == 0 ) + return tl::make_unexpected( "Failed to create mesh from mask" ); + + return mesh; +} + +// returns block of volume and new mask in space of new block, and shift of block +std::tuple simpleVolumeFromVoxelsMask( const ObjectVoxels& volume, const VoxelBitSet& mask, int expansion ) +{ + assert( volume.grid() ); + assert( mask.any() ); + + SimpleVolume volumePart; + volumePart.voxelSize = volume.voxelSize(); auto expandedMask = mask; - expandVoxelsMask( expandedMask, indexer, 25 ); // 25 should be in params + const auto& indexer = volume.getVolumeIndexer(); + expandVoxelsMask( expandedMask, indexer, expansion ); - auto accessor = volume.grid()->getConstAccessor(); - double insideAvg = 0.0; - double outsideAvg = 0.0; Box3i partBox; for ( auto voxelId : expandedMask ) { auto pos = indexer.toPos( voxelId ); partBox.include( pos ); - if ( mask.test( voxelId ) ) - insideAvg += double( accessor.getValue( { pos.x,pos.y,pos.z } ) ); + } + + volumePart.dims = partBox.size() + Vector3i::diagonal( 1 ); + const int newDimX = volumePart.dims.x; + const int newDimXY = newDimX * volumePart.dims.y; + volumePart.data.resize( newDimXY * volumePart.dims.z ); + + auto partIndexer = VolumeIndexer( volumePart.dims ); + VoxelBitSet volumePartMask( volumePart.data.size() ); + auto accessor = volume.grid()->getConstAccessor(); + for ( VoxelId i = VoxelId( 0 ); i < volumePart.data.size(); ++i ) + { + auto pos = partIndexer.toPos( i ) + partBox.min; + if ( mask.test( indexer.toVoxelId( pos ) ) ) + volumePartMask.set( i ); + volumePart.data[i] = accessor.getValue( { pos.x,pos.y,pos.z } ); + } + return { volumePart,volumePartMask,partBox.min }; +} + +enum class VolumeMaskMeshingMode +{ + Simple, + Smooth +}; + +void prepareVolumePart( SimpleVolume& volumePart, const VoxelBitSet& mask, VolumeMaskMeshingMode mode ) +{ + if ( mode == VolumeMaskMeshingMode::Simple ) + { + for ( VoxelId i = VoxelId( 0 ); i < volumePart.data.size(); ++i ) + { + volumePart.data[i] = mask.test( i ) ? 1.0f : 0.0f; + } + return; + } + // else: mode == VolumeMaskMeshingMode::Smooth + double insideAvg = 0.0; + double outsideAvg = 0.0; + for ( VoxelId i = VoxelId( 0 ); i < volumePart.data.size(); ++i ) + { + if ( mask.test( i ) ) + insideAvg += volumePart.data[i]; else - outsideAvg += double( accessor.getValue( { pos.x,pos.y,pos.z } ) ); + outsideAvg += volumePart.data[i]; } insideAvg /= double( mask.count() ); - outsideAvg /= double( expandedMask.count() - mask.count() ); + outsideAvg /= double( volumePart.data.size() - mask.count() ); auto range = float( insideAvg - outsideAvg ); - SimpleVolume voulmePart; - voulmePart.dims = partBox.size() + Vector3i::diagonal( 1 ); - + auto partIndexer = VolumeIndexer( volumePart.dims ); auto smallExpMask = mask; auto smallShrMask = mask; - expandVoxelsMask( smallExpMask, indexer, 3 ); - shrinkVoxelsMask( smallShrMask, indexer, 3 ); - - - const int newDimX = voulmePart.dims.x; - const size_t newDimXY = size_t( newDimX ) * voulmePart.dims.y; - voulmePart.data.resize( newDimXY * voulmePart.dims.z ); - for ( int z = partBox.min.z; z <= partBox.max.z; ++z ) - for ( int y = partBox.min.y; y <= partBox.max.y; ++y ) - for ( int x = partBox.min.x; x <= partBox.max.x; ++x ) - { - auto voxId = indexer.toVoxelId( { x,y,z } ); - if ( smallShrMask.test( voxId ) ) - voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = 1.0f; - else if ( smallExpMask.test( voxId ) ) - voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = - std::clamp( ( accessor.getValue( { x,y,z } ) - float( outsideAvg ) ) / range, 0.0f, 1.0f ); - else - voulmePart.data[x - partBox.min.x + ( y - partBox.min.y ) * newDimX + ( z - partBox.min.z ) * newDimXY] = 0.0f; - } - - auto voxelSize = volume.voxelSize(); - auto grid = simpleVolumeToDenseGrid( voulmePart ); - auto mesh = gridToMesh( grid, voxelSize, 0.5f ).value(); // no callback so cannot b stopped - - for ( auto& p : mesh.points ) + expandVoxelsMask( smallExpMask, partIndexer, 3 ); + shrinkVoxelsMask( smallShrMask, partIndexer, 3 ); + for ( VoxelId i = VoxelId( 0 ); i < volumePart.data.size(); ++i ) { - p = p + mult( Vector3f( partBox.min ), voxelSize ); + if ( smallShrMask.test( i ) ) + volumePart.data[i] = 1.0f; + else if ( smallExpMask.test( i ) ) + volumePart.data[i] = std::clamp( float( volumePart.data[i] - outsideAvg ) / range, 0.0f, 1.0f ); + else + volumePart.data[i] = 0.0f; } +} - if ( mesh.topology.numValidFaces() == 0 ) - return tl::make_unexpected( "Failed to create segmented mesh" ); +tl::expected meshFromVoxelsMask( const ObjectVoxels& volume, const VoxelBitSet& mask ) +{ + if ( !volume.grid() ) + return tl::make_unexpected( "Cannot create mesh from empty volume." ); + if ( mask.none() ) + return tl::make_unexpected( "Cannot create mesh from empty mask." ); - return mesh; + auto [volumePart, partMask, minVoxel] = simpleVolumeFromVoxelsMask( volume, mask, 25 ); + + prepareVolumePart( volumePart, partMask, VolumeMaskMeshingMode::Smooth ); + return meshFromSimpleVolume( volumePart, minVoxel ); } tl::expected segmentVolume( const ObjectVoxels& volume, const std::vector>& pairs, @@ -164,38 +209,10 @@ tl::expected VolumeSegmenter::segmentVolume( float seg tl::expected VolumeSegmenter::createMeshFromSegmentation( const VoxelBitSet& segmentation ) const { - auto newDims = volumePart_.dims; - auto newDimXY = newDims.x * newDims.y; - auto CoordToNewVoxelId = [&] ( const Vector3i& coord )->VoxelId - { - return VoxelId( coord.x + coord.y * newDims.x + coord.z * newDimXY ); - }; - auto segmentBlockCopy = volumePart_; - // Prepare block for meshing - for ( int z = 0; z < newDims.z; ++z ) - for ( int y = 0; y < newDims.y; ++y ) - for ( int x = 0; x < newDims.x; ++x ) - { - auto index = CoordToNewVoxelId( Vector3i( x, y, z ) ); - segmentBlockCopy.data[index] = segmentation.test( index ) ? 1.0f : 0.0f; - } - - - auto voxelSize = volume_.voxelSize(); - auto grid = simpleVolumeToDenseGrid( segmentBlockCopy ); - auto mesh = gridToMesh( grid, voxelSize, 0.5f ).value(); // no callback so cannot b stopped - - - for ( auto& p : mesh.points ) - { - p = p + mult( Vector3f( minVoxel_ ), voxelSize ); - } - - if ( mesh.topology.numValidFaces() == 0 ) - return tl::make_unexpected( "Failed to create segmented mesh" ); - - return mesh; + segmentBlockCopy.voxelSize = volume_.voxelSize(); + prepareVolumePart( segmentBlockCopy, segmentation, VolumeMaskMeshingMode::Simple ); + return meshFromSimpleVolume( segmentBlockCopy, minVoxel_ ); } const MR::Vector3i& VolumeSegmenter::getVolumePartDimensions() const From a7b2df681519924412340581580529ee9cd8b27a Mon Sep 17 00:00:00 2001 From: Grant Karapetyan Date: Tue, 30 Aug 2022 14:35:02 +0300 Subject: [PATCH 3/3] Add some comments --- source/MRMesh/MRVolumeSegment.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/source/MRMesh/MRVolumeSegment.cpp b/source/MRMesh/MRVolumeSegment.cpp index 4ff18eeeefbb..9d0a42dd3973 100644 --- a/source/MRMesh/MRVolumeSegment.cpp +++ b/source/MRMesh/MRVolumeSegment.cpp @@ -15,7 +15,7 @@ namespace MR tl::expected meshFromSimpleVolume( const SimpleVolume& volumePart, const Vector3i& shift ) { auto grid = simpleVolumeToDenseGrid( volumePart ); - auto mesh = gridToMesh( grid, volumePart.voxelSize, 0.5f ).value(); // no callback so cannot b stopped + auto mesh = gridToMesh( grid, volumePart.voxelSize, 0.5f ).value(); // no callback so cannot be stopped auto minCorner = mult( Vector3f( shift ), volumePart.voxelSize ); for ( auto& p : mesh.points ) @@ -66,12 +66,14 @@ std::tuple simpleVolumeFromVoxelsMask( cons return { volumePart,volumePartMask,partBox.min }; } +// Mode of meshing enum class VolumeMaskMeshingMode { - Simple, - Smooth + Simple, // 1 inside, 0 outside, 0.5 iso + Smooth // 1 deep inside, (density - outsideAvgDensity)/(insideAvgDensity - outsideAvgDensity) on the edge, 0 - far outside, 0.5 iso }; +// changes volume part due to meshing mode void prepareVolumePart( SimpleVolume& volumePart, const VoxelBitSet& mask, VolumeMaskMeshingMode mode ) { if ( mode == VolumeMaskMeshingMode::Simple )