Skip to content

Commit

Permalink
More cleanup and improvements
Browse files Browse the repository at this point in the history
  • Loading branch information
@chevtche
chevtche committed Mar 15, 2017
1 parent 5b6080a commit a9a43ef
Show file tree
Hide file tree
Showing 10 changed files with 541 additions and 361 deletions.
5 changes: 4 additions & 1 deletion apps/Mesher/CMakeLists.txt
View file
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)
245 changes: 58 additions & 187 deletions apps/Mesher/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,202 +19,73 @@
* 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 = -1.0f;

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

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

desc.add_options()
("help,h", "Print 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")
("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")
("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.")
("min-radius", po::value< float >( &minRadius ),
"All radii in the morphology smaller than the specified value will be "
"increased to that value.")
("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 0;
}
Mesher mesher( inputFile, alpha, angularBound, radialBound, distanceBound,
leafSize, minRadius );
mesher.mesh( outputFile );

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 );

// defining meshing criteria
std::cout << "Triangle Size: " << triangleSize << std::endl;
CGAL::Surface_mesh_default_criteria_3<Tr> criteria( 30, triangleSize, triangleSize );

// 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
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
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

0 comments on commit a9a43ef

Please sign in to comment.