Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Sparse octree. Work in progress. #1

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Sparse octree. Work in progress. #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
a9a43ef
More cleanup and improvements
chevtche Mar 15, 2017
12fba2f
Added CGAL decimation
chevtche Mar 17, 2017
087230c
CGAL self intersection test. Automatic parameters computation.
Mar 19, 2017
5664b9d
Fixing self intersection. Work in progress.
chevtche Mar 21, 2017
1042dfd
Fixing self intersection. Work in progress 2.
chevtche Mar 22, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
5 changes: 4 additions & 1 deletion apps/Mesher/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,10 @@
# along with this library; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.

find_package(CGAL)

set(MESHER_SOURCES main.cpp)
set(MESHER_LINK_LIBRARIES PUBLIC SIMDVoxelizerCommon)
set(MESHER_LINK_LIBRARIES PUBLIC SIMDVoxelizerCommon CGAL CGAL_Core CGAL_ImageIO
boost_thread boost_program_options gmp mpfr)

common_application(Mesher)
260 changes: 75 additions & 185 deletions apps/Mesher/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,202 +19,92 @@
* 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*/

#include <simdvoxelizer/SIMDSparseVoxelizer_ispc.h>
#include <simdvoxelizer/SIMDVoxelizer_ispc.h>
#include <simdvoxelizer/Octree.h>
#include <string.h>
#include <cmath>
#include <fstream>
#include <vector>
#include <map>
#include <limits>
#include <iostream>
#include <sstream>
#include <ctime>
#include <cstdlib>
#include <boost/program_options.hpp>

#if 0 // WORK IN PROGRESS
#include <CGAL/Surface_mesh_default_triangulation_3.h>
#include <CGAL/Surface_mesh_default_criteria_3.h>
#include <CGAL/Complex_2_in_triangulation_3.h>
#include <CGAL/IO/Complex_2_in_triangulation_3_file_writer.h>
#include <fstream>
#include <CGAL/make_surface_mesh.h>
#include <CGAL/Gray_level_image_3.h>
#include <CGAL/Implicit_surface_3.h>

// default triangulation for Surface_mesher
typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
// c2t3
typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
typedef Tr::Geom_traits GT;
typedef CGAL::Gray_level_image_3<GT::FT, GT::Point_3> Gray_level_image;
typedef CGAL::Implicit_surface_3<GT, Gray_level_image> Surface_3;

using namespace ispc;

int span = 32;
float voxelSize = 2.f;
float threshold = 15.f;
std::string inputFile;
std::string outputFile;
float triangleSize = 5;
int generateVolume = 1;
#include <simdvoxelizer/Mesher.h>

int main( int argc, char* argv[] )
{
if( argc != 8 )
{
std::cerr << "usage: SIMDVoxelizer <voxel_size> <span> <input_file> <output_file> <threshold> <triangle_size> <generate_volume>" << std::endl;
exit(1);
}

voxelSize = atof(argv[1]);
span = atoi(argv[2]);
inputFile = argv[3];
outputFile = argv[4];
threshold = atof(argv[5]);
triangleSize = atof(argv[6]);
generateVolume = atoi(argv[7]);

glm::uvec3 volumeDim( 534, 829, 253);
if( generateVolume == 1 )
float alpha = 3.0f;
float angularBound = -1.0f;
float radialBound = -1.0f;
float distanceBound = -1.0f;
float leafSize = -1.0f;
float minRadius = 0.1f;
float decimationRatio = 1.0f;
float epsilon = 0.0000033f;
uint32_t targetFaceNumber = 0u;

std::string outputFile;
std::string inputFile;

namespace po = boost::program_options;
po::options_description desc("");

desc.add_options()
("help,h", "Print this help messages.")
("input,i", po::value< std::string >( &inputFile )->required(),
"Input file containing the morphology.")
("output,o", po::value< std::string >( &outputFile )->required(),
"Output file that will contain the mesh.\n")
("simplification-ratio,s", po::value< float >( &decimationRatio )->default_value( 1.0f ),
"The decimation ratio by which the mesh must be decimated.\n")
("target-faces,t", po::value< uint32_t >( &targetFaceNumber ),
"The target number of faces for the decimated mesh. This parameter overwrite the"
" --simplification-ratio parameter.\n")
("min-radius", po::value< float >( &minRadius )->default_value( 0.1f ),
"All radii in the morphology smaller than the specified value will be "
"increased to that value.")
("angular-bound,a", po::value< float >( &angularBound ),
"A lower bound in degrees for the angles of mesh facets. If not specified "
"this parameter is automatically computed.")
("radial-bound,r", po::value< float >( &radialBound ),
"An upper bound on the radii of surface Delaunay balls. A surface "
"Delaunay ball is a ball circumscribing a mesh facet and centered on the "
"surface. If not specified this parameter is automatically computed.")
("distance-bound,d", po::value< float >( &distanceBound ),
"An upper bound for the distance between the circumcenter of a mesh facet "
"and the center of a surface Delaunay ball of this facet. If not specified "
"this parameter is automatically computed.\n")
("epsilon,e", po::value< float >( &epsilon )->default_value( 0.0000033f ),
"Epsilon multiplyed by the radius of the bounding sphere is the precision at"
"which the dichotomy will stop.\n")
("leaf-size,l", po::value< float >( &leafSize ),
"The size of the octree's leaves in morphology space. If not specified "
"this parameter is automatically computed.")
("alpha", po::value< float >( &alpha )->default_value( 3.0f ),
"Octree's nodes further than closest radius * alpha are discarded");

po::variables_map vm;

try
{
std::vector<float> events;
std::ifstream file(inputFile.c_str(), std::ios::in | std::ios::binary);
if( file.is_open( ))
{
while( !file.eof( ))
{
float v;
file.read( (char *)&v, sizeof(float));
events.push_back(v);
}
file.close();
}

std::cout << "--------------------------------------------" << std::endl;
std::cout << "SIMDVoxelizer" << std::endl;
std::cout << "--------------------------------------------" << std::endl;
std::cout << "Voxel size : " << voxelSize << std::endl;
std::cout << "Span : " << span << std::endl;
std::cout << "Input file : " << inputFile << std::endl;
std::cout << "Output file : " << outputFile << std::endl;
std::cout << "--------------------------------------------" << std::endl;

Octree morphoOctree( events, voxelSize );
uint64_t volumeSize = morphoOctree.getVolumeSize();
volumeDim = morphoOctree.getVolumeDim();

float* volume = new float[ volumeSize ];

std::cout<< "Volume dims: " << volumeDim.x << " " << volumeDim.y << " " << volumeDim.z << " " << volumeSize << std::endl;

uint32_t zLenght = 32;
uint32_t nPasses = std::ceil( volumeDim.z / zLenght );

for( uint32_t zOffset = 0; zOffset < volumeDim.z; zOffset += zLenght )
{
std::cout << "z: " << zOffset << std::endl;
SIMDSparseVoxelizer_ispc( zOffset, zOffset + zLenght, span, voxelSize,
morphoOctree.getOctreeSize(), volumeDim.x,
volumeDim.y, volumeDim.z,
morphoOctree.getFlatIndexes(),
morphoOctree.getFlatData(), volume );
}

float minValue = std::numeric_limits<float>::max();
float maxValue = -std::numeric_limits<float>::max();
for( int i = 0; i < volumeSize; ++i )
{
if( volume[i] != 0.f )
{
minValue = std::min( minValue, volume[i] );
maxValue = std::max( maxValue, volume[i] );
}
}
po::store( po::parse_command_line( argc, argv, desc ), vm );

//float a = 255.f / ( maxValue - minValue );
float a = 9.f / ( maxValue - minValue );
std::cout << "Normalization [" << minValue << " - " << maxValue << "] " << a << std::endl;

char* volumeAsChar = new char[ volumeSize ];
for( int i = 0; i < volumeSize; ++i )
if( vm.count( "help" ))
{
//float normalizedValue = (volume[i] - minValue) * a;
float normalizedValue = 255.0f * std::log((volume[i] - minValue) * a + 1.0);
volumeAsChar[i] = (uint8_t)normalizedValue;
std::cout << desc << std::endl;
return 0;
}
std::cout << std::endl;
std::cout << "--------------------------------------------" << std::endl;

std::ofstream volumeFile(
outputFile.c_str(), std::ios::out | std::ios::binary);
std::stringstream header;
header << "#INRIMAGE-4#{\n";
header << "XDIM=" << volumeDim.x << "\n";
header << "YDIM=" << volumeDim.y << "\n";
header << "ZDIM=" << volumeDim.z << "\n";
header << "VDIM=1\n";
header << "TYPE=unsigned fixed\n";
header << "SCALE=1**0\n";
header << "PIXSIZE=8 bits\n";
header << "CPU=pc\n";
header << "VX=1\n";
header << "VY=1\n";
header << "VZ=1\n";
const size_t len = 256 - 4 - header.str().length();
for( size_t i =0 ; i < len; ++i )
header << "\n";
header << "##}\n";

volumeFile.write( (char*)header.str().c_str(), sizeof(char) * header.str().length());
volumeFile.write( (char*)&volumeAsChar[0], sizeof(char) * volumeSize );
volumeFile.close();
delete [] volume;
delete [] volumeAsChar;
po::notify( vm );
}
catch( boost::program_options::required_option& e )
{
std::cerr << "ERROR: " << e.what() << std::endl << std::endl;
std::cout << desc << std::endl;
return 1;
}

std::cout << "Meshing..." << std::endl;
Tr tr; // 3D-Delaunay triangulation
C2t3 c2t3 (tr); // 2D-complex in 3D-Delaunay triangulation
// the 'function' is a 3D gray level image
Gray_level_image image( outputFile.c_str(), threshold );

// Carefully choosen bounding sphere: the center must be inside the
// surface defined by 'image' and the radius must be high enough so that
// the sphere actually bounds the whole image.
GT::Point_3 bounding_sphere_center(
0.498567f * volumeDim.x,
0.199998f * volumeDim.y,
0.609048f * volumeDim.z );
const float radius =
std::max( volumeDim.x, std::max( volumeDim.y, volumeDim.z ));
std::cout << "Radius: " << radius << std::endl;
GT::FT bounding_sphere_squared_radius = radius * radius;
GT::Sphere_3 bounding_sphere(bounding_sphere_center,
bounding_sphere_squared_radius);

// definition of the surface, with 10^-5 as relative precision
Surface_3 surface( image, bounding_sphere, 1e-7 );
if( minRadius <= 0 )
{
std::cerr << " ERROR: min-radius is: " << minRadius
<<". The algorithm will not converge." << std::endl;
return 1;
}

// defining meshing criteria
std::cout << "Triangle Size: " << triangleSize << std::endl;
CGAL::Surface_mesh_default_criteria_3<Tr> criteria( 30, triangleSize, triangleSize );
Mesher mesher( inputFile, alpha, angularBound, radialBound, distanceBound,
leafSize, minRadius, decimationRatio, epsilon, targetFaceNumber );
mesher.mesh( outputFile );

// meshing surface, with the "manifold without boundary" algorithm
//CGAL::make_surface_mesh( c2t3, surface, criteria, CGAL::Manifold_tag( ));
CGAL::make_surface_mesh( c2t3, surface, criteria, CGAL::Non_manifold_tag());
std::ofstream out( "out.off" );
CGAL::output_surface_facets_to_off( out, c2t3 );
std::cout << "Final number of points: " << tr.number_of_vertices() << "\n";
return 0;
}
#else
int main( int , char* argv[] )
{
std::cout << argv << std::endl;
return 1;
}

#endif
38 changes: 32 additions & 6 deletions apps/Voxelizer/main.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -62,15 +62,15 @@ void flattenChildren( const OctreeNode* node, uint32_t* offsetPerLevel,
flatOctreeData[ offsetPerLevel[ level ] * 4u ] = node->getCenter().x;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 1 ] = node->getCenter().y;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 2 ] = node->getCenter().z;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 3 ] = node->getValue();
flatOctreeData[ offsetPerLevel[ level ] * 4u + 3 ] = node->getMaxValue();

offsetPerLevel[ level ] += 1u;
return;
}
flatOctreeData[ offsetPerLevel[ level ] * 4u ] = node->getCenter().x;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 1 ] = node->getCenter().y;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 2 ] = node->getCenter().z;
flatOctreeData[ offsetPerLevel[ level ] * 4u + 3 ] = node->getValue();
flatOctreeData[ offsetPerLevel[ level ] * 4u + 3 ] = node->getMaxValue();

flatOctreeIndex[ offsetPerLevel[ level ] * 2u ] =
offsetPerLevel[ level - 1 ];
Expand All @@ -83,9 +83,9 @@ void flattenChildren( const OctreeNode* node, uint32_t* offsetPerLevel,
level - 1u );
}

int main( int argc, char* argv[] )
int main( /*int argc, char* argv[]*/ )
{
if( argc != 5 )
/*if( argc != 5 )
{
std::cerr << "usage: SIMDVoxelizer <voxel_size> <span> "
<< "<input_file> <output_file>" << std::endl;
Expand All @@ -109,6 +109,32 @@ int main( int argc, char* argv[] )
file.close();
}

// Bounding box
glm::vec3 min;
min.x = std::numeric_limits<float>::max();
min.y = std::numeric_limits<float>::max();
min.z = std::numeric_limits<float>::max();

glm::vec3 max;
max.x = -std::numeric_limits<float>::max();
max.y = -std::numeric_limits<float>::max();
max.z = -std::numeric_limits<float>::max();

std::cout << "Nb of events :" << events.size() / 5 << std::endl;
for( size_t i = 0; i < events.size(); i += 5 )
{
min.x = std::min( min.x, events[i] );
min.y = std::min( min.y, events[i+1] );
min.z = std::min( min.z, events[i+2] );

max.x = std::max( max.x, events[i] );
max.y = std::max( max.y, events[i+1] );
max.z = std::max( max.z, events[i+2] );
}

//erase the soma position ( padded with zeros )-> posX posY posZ 0.0f 0.0f
events.erase( events.begin(), events.begin() + 5 );

std::cout << "--------------------------------------------" << std::endl;
std::cout << "SIMDVoxelizer" << std::endl;
std::cout << "--------------------------------------------" << std::endl;
Expand All @@ -118,7 +144,7 @@ int main( int argc, char* argv[] )
std::cout << "Output file : " << outputFile << std::endl;
std::cout << "--------------------------------------------" << std::endl;

Octree morphoOctree( events, voxelSize );
Octree morphoOctree( events, voxelSize, min, max );

uint64_t volumeSize = morphoOctree.getVolumeSize();
glm::uvec3 volumeDim = morphoOctree.getVolumeDim();
Expand Down Expand Up @@ -168,7 +194,7 @@ int main( int argc, char* argv[] )
std::ofstream volumeFile(
outputFile.c_str(), std::ios::out | std::ios::binary);
volumeFile.write( (char*)volumeAsChar.data(), sizeof(char) * volumeSize );
volumeFile.close();
volumeFile.close();*/

return 0;
}
5 changes: 4 additions & 1 deletion simdvoxelizer/CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,11 +86,13 @@ endmacro()

set(SIMDVOXELIZERCOMMON_SOURCES
ispc/tasksys.cpp
Mesher.cpp
OctreeNode.cpp
Octree.cpp
)

set(SIMDVOXELIZERCOMMON_PUBLIC_HEADERS
Mesher.h
OctreeNode.h
Octree.h
)
Expand All @@ -106,6 +108,7 @@ INCLUDE_DIRECTORIES_ISPC(${SIMDVOXELIZER_INCLUDE_DIRS})
ISPC_COMPILE(${SIMDVOXELIZERCOMMON_ISPC_SOURCES})
list(APPEND SIMDVOXELIZERCOMMON_SOURCES ${ISPC_OBJECTS})

set(SIMDVOXELIZERCOMMON_LINK_LIBRARIES PUBLIC glm pthread stdc++)
set(SIMDVOXELIZERCOMMON_LINK_LIBRARIES PUBLIC glm pthread stdc++ CGAL
CGAL_Core CGAL_ImageIO boost_thread gmp mpfr)

common_library(SIMDVoxelizerCommon)
Loading