From 831bfbd3459d268f89d7445f7aafaaa42176314a Mon Sep 17 00:00:00 2001
From: David Coeurjolly <david.coeurjolly@liris.cnrs.fr>
Date: Tue, 31 Oct 2017 15:55:21 +0100
Subject: [PATCH] better scaling

---
 converters/mesh2vol.cpp | 79 ++++++++++++++++++++++++-----------------
 1 file changed, 47 insertions(+), 32 deletions(-)

diff --git a/converters/mesh2vol.cpp b/converters/mesh2vol.cpp
index fb63b4d6..4891e180 100644
--- a/converters/mesh2vol.cpp
+++ b/converters/mesh2vol.cpp
@@ -28,7 +28,7 @@
 #include <iostream>
 #include <fstream>
 #include <chrono>
-
+#include "DGtal/helpers/StdDefs.h"
 #include "DGtal/shapes/MeshVoxelizer.h"
 #include "DGtal/kernel/sets/CDigitalSet.h"
 #include "DGtal/kernel/domains/HyperRectDomain.h"
@@ -61,14 +61,13 @@ namespace po = boost::program_options;
   -o [ --output ] arg             filename of ouput voxelization file (.off)
                                   (auto-generated by argument values if empty)
   -s [ --separation ] arg (=6)    voxelization 6-separated or 26-separated.
-  -d [ --domainSize ] arg (=64)   digitization domain size.
-  -f [ --scaleFactor ] arg (=1.0) scale factor to apply to the mesh
+  -r [ --resolution ]             digitization domain size (e.g. 128). The mesh will be scaled such that its bounding box maps to [0,resolution)^3.
 
 @endcode
 
 @b Example:
 @code
-  $ mesh2vol -i ${DGtal}/examples/samples/tref.off --separation 26 --resolution 256
+  $ mesh2vol -i ${DGtal}/examples/samples/tref.off --separation 26 --resolution 256 -o output.vol
 @endcode
 
 @see mesh2vol.cpp
@@ -78,41 +77,50 @@ namespace po = boost::program_options;
 template< unsigned int SEP >
 void voxelizeAndExport(const std::string inputFilename,
                        const std::string outputFilename,
-                       const unsigned int resolution,
-                       const double scaleFactor)
+                       const unsigned int resolution)
 {
+  using Domain   = Z3i::Domain;
+  using PointR3  = Z3i::RealPoint;
+  using PointZ3  = Z3i::Point;
 
-  using Space3Dint = SpaceND<3>;
-  using Domain   = HyperRectDomain<Space3Dint>;
-  using PointR3  = PointVector<3, double>;
-  using PointZ3  = PointVector<3, int>;
-
+  trace.beginBlock("Preparing the mesh");
   trace.info() << "Reading input file: " << inputFilename;
-
   Mesh<PointR3> inputMesh;
-
-  unsigned int base = resolution/2; // move center of digital set to (0;0;0)
-  Domain aDomain( PointZ3(-base, -base, -base),
-                  PointZ3( base,  base,  base) );
-
   MeshReader<PointR3>::importOFFFile(inputFilename.c_str(), inputMesh);
   trace.info() << " [done]" << std::endl;
+  const std::pair<PointR3, PointR3> bbox = inputMesh.getBoundingBox();
+  trace.info()<< "Mesh bounding box: "<<bbox.first <<" "<<bbox.second<<std::endl;
+
+  const double smax = (bbox.second - bbox.first).max();
+  const double factor = resolution / smax;
+  const PointR3 translate = -bbox.first;
+  trace.info() << "Scale = "<<factor<<" translate = "<<translate<<std::endl;
+  for(auto it = inputMesh.vertexBegin(), itend = inputMesh.vertexEnd();
+      it != itend; ++it)
+  {
+    //scale + translation
+    *it += translate;
+    *it *= factor;
+  }
+  trace.endBlock();
+  
+  trace.beginBlock("Voxelization");
   trace.info() << "Voxelization " << SEP << "-separated ; " << resolution << "^3 ";
-
+  Domain aDomain(PointZ3().diagonal(0), PointZ3().diagonal(resolution));
+  
   //Digitization step
-  DigitalSetBySTLSet<Domain> mySet(aDomain);
-  MeshVoxelizer<DigitalSetBySTLSet<Domain>, SEP> aVoxelizer;
-  double factor = scaleFactor;
-  if(factor == 1.0)
-    factor = resolution / inputMesh.getBoundingBox().second[0];
+  Z3i::DigitalSet mySet(aDomain);
+  MeshVoxelizer<Z3i::DigitalSet, SEP> aVoxelizer;
   
   auto start = std::chrono::high_resolution_clock::now();
-  aVoxelizer.voxelize(mySet, inputMesh, factor);
+  aVoxelizer.voxelize(mySet, inputMesh, 1.0);
   auto end = std::chrono::high_resolution_clock::now();
   auto diff = end - start;
 
   trace.info() << " [done] (" << std::chrono::duration <double, std::milli> (diff).count() << " ms)" << std::endl;
-
+  trace.endBlock();
+  
+  trace.beginBlock("Exporting");
   // Export the digital set to a vol file
   trace.info()<<aDomain<<std::endl;
   ImageContainerBySTLVector<Domain, unsigned char> image(aDomain);
@@ -122,6 +130,7 @@ void voxelizeAndExport(const std::string inputFilename,
     image.setValue(p, 128);
   }
   image >> outputFilename.c_str();
+  trace.endBlock();
 }
 
 int main( int argc, char** argv )
@@ -133,8 +142,7 @@ int main( int argc, char** argv )
     ("input,i", po::value<std::string>(), "mesh file (.off) " )
     ("output,o", po::value<std::string>(), "filename of ouput volumetric file (Vol,PGM3D,..).")
     ("separation,s", po::value<unsigned int>()->default_value(6), "voxelization 6-separated or 26-separated." )
-    ("scaleFactor,f", po::value<double>()->default_value(1.0), "scale factor to apply to the mesh." )
-    ("domainSize,d", po::value<unsigned int>()->default_value(64), "digitization domain size." );
+    ("resolution,r", po::value<unsigned int>(), "digitization domain size (e.g. 128). The mesh will be scaled such that its bounding box maps to [0,resolution)^3." );
 
   bool parseOK=true;
   po::variables_map vm;
@@ -153,7 +161,7 @@ int main( int argc, char** argv )
               << "Convert a mesh file into a 26-separated or 6-separated volumetric voxelization in a given resolution grid."
               << general_opt << "\n";
     std::cout << "Example:\n"
-              << "mesh2vol -i ${DGtal}/examples/samples/tref.off --separation 26 --resolution 256 \n";
+              << "mesh2vol -i ${DGtal}/examples/samples/tref.off -o output.vol --separation 26 --resolution 256 \n";
     return -1;
   }
   
@@ -171,10 +179,17 @@ int main( int argc, char** argv )
   }
 
   unsigned int separation = vm["separation"].as<unsigned int>();
-  unsigned int resolution = vm["domainSize"].as<unsigned int>();
+  unsigned int resolution;
+  if (vm.count("resolution"))
+    resolution = vm["resolution"].as<unsigned int>();
+  else
+  {
+    trace.error()<<"Missing output parameter '--resolution'"<<std::endl;
+    return -1;
+  }
+  
   string inputFilename = vm["input"].as<std::string>();
   string outputFilename;
-  double scaleFactor = vm["scaleFactor"].as<double>();
   
   if( vm.count("output") )
     outputFilename = vm["output"].as<std::string>();
@@ -186,9 +201,9 @@ int main( int argc, char** argv )
   
   
   if(separation == 6)
-    voxelizeAndExport<6>(inputFilename, outputFilename, resolution, scaleFactor);
+    voxelizeAndExport<6>(inputFilename, outputFilename, resolution);
   else if(separation == 26)
-    voxelizeAndExport<26>(inputFilename, outputFilename, resolution, scaleFactor);
+    voxelizeAndExport<26>(inputFilename, outputFilename, resolution);
 
   return 0;
 }