diff --git a/.gitignore b/.gitignore
index f59c916d00..dc37045e73 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,5 +1,6 @@
 # temporary files
 *~
+.nfs*
 .*.swp # vim
 *.flc # xemacs
 .DS_Store # MacOS
diff --git a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
index 3acba733f5..91af70f754 100644
--- a/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
+++ b/src/aliceVision/fuseCut/DelaunayGraphCut.cpp
@@ -3271,7 +3271,7 @@ void DelaunayGraphCut::voteFullEmptyScore(const StaticVector<int>& cams, const s
         if(false)
         {
             std::unique_ptr<mesh::Mesh> meshf(createTetrahedralMesh(false, 0.9f, [](const fuseCut::GC_cellInfo& c) { return c.emptinessScore; }));
-            meshf->saveToObj(folderName + "tetrahedralMesh_beforeForceTEdge_emptiness.obj");
+            meshf->save(folderName + "tetrahedralMesh_beforeForceTEdge_emptiness");
         }
 
         if(forceTEdge)
@@ -3789,47 +3789,47 @@ void DelaunayGraphCut::exportDebugMesh(const std::string& filename, const Point3
     }
 
     const std::string tempDirPath = boost::filesystem::temp_directory_path().generic_string();
-    mesh->saveToObj(tempDirPath + "/" + filename + ".obj");
-    meshf->saveToObj(tempDirPath + "/" + filename + "Filtered.obj");
+    mesh->save(tempDirPath + "/" + filename);
+    meshf->save(tempDirPath + "/" + filename);
 }
 
 void DelaunayGraphCut::exportFullScoreMeshs(const std::string& outputFolder, const std::string& name) const
 {
-    const std::string nameExt = (name.empty() ? "" : "_" + name) + ".obj";
+    const std::string nameExt = (name.empty() ? "" : "_" + name);
     {
         std::unique_ptr<mesh::Mesh> meshEmptiness(
             createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.emptinessScore; }));
-        meshEmptiness->saveToObj(outputFolder + "/mesh_emptiness" + nameExt);
+        meshEmptiness->save(outputFolder + "/mesh_emptiness" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> meshFullness(
             createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.fullnessScore; }));
-        meshFullness->saveToObj(outputFolder + "/mesh_fullness" + nameExt);
+        meshFullness->save(outputFolder + "/mesh_fullness" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> meshSWeight(
             createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellSWeight; }));
-        meshSWeight->saveToObj(outputFolder + "/mesh_sWeight" + nameExt);
+        meshSWeight->save(outputFolder + "/mesh_sWeight" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> meshTWeight(
             createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.cellTWeight; }));
-        meshTWeight->saveToObj(outputFolder + "/mesh_tWeight" + nameExt);
+        meshTWeight->save(outputFolder + "/mesh_tWeight" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> meshOn(
             createTetrahedralMesh(false, 0.999f, [](const GC_cellInfo& c) { return c.on; }));
-        meshOn->saveToObj(outputFolder + "/mesh_on" + nameExt);
+        meshOn->save(outputFolder + "/mesh_on" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> mesh(createTetrahedralMesh(
             false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.fullnessScore - c.emptinessScore; }));
-        mesh->saveToObj(outputFolder + "/mesh_fullness-emptiness" + nameExt);
+        mesh->save(outputFolder + "/mesh_fullness-emptiness" + nameExt);
     }
     {
         std::unique_ptr<mesh::Mesh> mesh(createTetrahedralMesh(
             false, 0.99f, [](const fuseCut::GC_cellInfo& c) { return c.cellSWeight - c.cellTWeight; }));
-        mesh->saveToObj(outputFolder + "/mesh_s-t" + nameExt);
+        mesh->save(outputFolder + "/mesh_s-t" + nameExt);
     }
 }
 
diff --git a/src/aliceVision/mesh/Mesh.cpp b/src/aliceVision/mesh/Mesh.cpp
index 8c918735c6..24f158fd18 100644
--- a/src/aliceVision/mesh/Mesh.cpp
+++ b/src/aliceVision/mesh/Mesh.cpp
@@ -14,6 +14,7 @@
 #include <geogram/points/kd_tree.h>
 
 #include <boost/filesystem.hpp>
+#include <boost/algorithm/string/case_conv.hpp> 
 
 #include <assimp/Importer.hpp>
 #include <assimp/Exporter.hpp>
@@ -36,9 +37,56 @@ Mesh::~Mesh()
 {
 }
 
-void Mesh::saveToObj(const std::string& filename)
+
+std::string EFileType_enumToString(const EFileType meshFileType)
 {
-    ALICEVISION_LOG_INFO("Writing obj and mtl file.");
+    switch(meshFileType)
+    {
+        case EFileType::OBJ:
+            return "obj";
+        case EFileType::FBX:
+            return "fbx";
+        case EFileType::STL:
+            return "stl";
+        case EFileType::GLTF:
+            return "gltf";
+    }
+    throw std::out_of_range("Unrecognized EMeshFileType");
+}
+
+EFileType EFileType_stringToEnum(const std::string& meshFileType)
+{
+    std::string m = meshFileType;
+    boost::to_lower(m);
+
+    if(m == "obj")
+        return EFileType::OBJ;
+    if(m == "fbx")
+        return EFileType::FBX;
+    if(m == "stl")
+        return EFileType::STL;
+    if(m == "gltf")
+        return EFileType::GLTF;
+    throw std::out_of_range("Invalid mesh file type " + meshFileType);
+}
+
+std::ostream& operator<<(std::ostream& os, EFileType meshFileType)
+{
+    return os << EFileType_enumToString(meshFileType);
+}
+std::istream& operator>>(std::istream& in, EFileType& meshFileType)
+{
+    std::string token;
+    in >> token;
+    meshFileType = EFileType_stringToEnum(token);
+    return in;
+}
+
+void Mesh::save(const std::string& filepath, EFileType fileType)
+{
+    const std::string fileTypeStr = EFileType_enumToString(fileType);
+
+    ALICEVISION_LOG_INFO("Save " << fileTypeStr << " mesh file");
 
     aiScene scene;
 
@@ -85,15 +133,36 @@ void Mesh::saveToObj(const std::string& filename)
         }
     }
 
+    std::string formatId = fileTypeStr;
+    unsigned int pPreprocessing = 0u;
+    // If gltf, use gltf 2.0
+    if (fileType == EFileType::GLTF)
+    {
+        formatId = "gltf2";
+        // gen normals in order to have correct shading in Qt 3D Scene
+        // but cause problems with assimp importer
+        pPreprocessing |= aiProcess_GenNormals;
+    }
+    // If obj, do not use material
+    else if (fileType == EFileType::OBJ)
+    {
+        formatId = "objnomtl";
+    }
+
     Assimp::Exporter exporter;
-    exporter.Export(&scene, "objnomtl", filename);
+    exporter.Export(&scene, formatId, filepath, pPreprocessing);
+
+    ALICEVISION_LOG_INFO("Save mesh to " << fileTypeStr << " done.");
 
-    ALICEVISION_LOG_INFO("Save mesh to obj done.");
+    ALICEVISION_LOG_DEBUG("Vertices: " << pts.size());
+    ALICEVISION_LOG_DEBUG("Triangles: " << tris.size());
+    ALICEVISION_LOG_DEBUG("UVs: " << uvCoords.size());  
+    ALICEVISION_LOG_DEBUG("Normals: " << normals.size());
 }
 
-bool Mesh::loadFromBin(const std::string& binFileName)
+bool Mesh::loadFromBin(const std::string& binFilepath)
 {
-    FILE* f = fopen(binFileName.c_str(), "rb");
+    FILE* f = fopen(binFilepath.c_str(), "rb");
 
     if(f == nullptr)
         return false;
@@ -114,12 +183,12 @@ bool Mesh::loadFromBin(const std::string& binFileName)
     return true;
 }
 
-void Mesh::saveToBin(const std::string& binFileName)
+void Mesh::saveToBin(const std::string& binFilepath)
 {
     long t = std::clock();
     ALICEVISION_LOG_DEBUG("Save mesh to bin.");
     // printf("open\n");
-    FILE* f = fopen(binFileName.c_str(), "wb");
+    FILE* f = fopen(binFilepath.c_str(), "wb");
 
     int npts = pts.size();
     // printf("write npts %i\n",npts);
@@ -885,26 +954,26 @@ void Mesh::getDepthMap(StaticVector<float>& depthMap, StaticVector<StaticVector<
     }     // for pix.x
 }
 
-void Mesh::getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& depthMapFileName, const std::string& trisMapFileName,
+void Mesh::getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& depthMapFilepath, const std::string& trisMapFilepath,
                                                        const mvsUtils::MultiViewParams& mp, int rc, int w, int h)
 {
     StaticVector<float> depthMap;
-    loadArrayFromFile<float>(depthMap, depthMapFileName);
+    loadArrayFromFile<float>(depthMap, depthMapFilepath);
     StaticVector<StaticVector<int>> trisMap;
-    loadArrayOfArraysFromFile<int>(trisMap, trisMapFileName);
+    loadArrayOfArraysFromFile<int>(trisMap, trisMapFilepath);
 
     getVisibleTrianglesIndexes(out_visTri, trisMap, depthMap, mp, rc, w, h);
 }
 
 void Mesh::getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& tmpDir, const mvsUtils::MultiViewParams& mp, int rc, int w, int h)
 {
-    std::string depthMapFileName = tmpDir + "depthMap" + std::to_string(mp.getViewId(rc)) + ".bin";
-    std::string trisMapFileName = tmpDir + "trisMap" + std::to_string(mp.getViewId(rc)) + ".bin";
+    std::string depthMapFilepath = tmpDir + "depthMap" + std::to_string(mp.getViewId(rc)) + ".bin";
+    std::string trisMapFilepath = tmpDir + "trisMap" + std::to_string(mp.getViewId(rc)) + ".bin";
 
     StaticVector<float> depthMap;
-    loadArrayFromFile<float>(depthMap, depthMapFileName);
+    loadArrayFromFile<float>(depthMap, depthMapFilepath);
     StaticVector<StaticVector<int>> trisMap;
-    loadArrayOfArraysFromFile<int>(trisMap, trisMapFileName);
+    loadArrayOfArraysFromFile<int>(trisMap, trisMapFilepath);
 
     getVisibleTrianglesIndexes(out_visTri, trisMap, depthMap, mp, rc, w, h);
 }
@@ -1855,9 +1924,9 @@ void Mesh::computeTrisCams(StaticVector<StaticVector<int>>& trisCams, const mvsU
     long t1 = mvsUtils::initEstimate();
     for(int rc = 0; rc < mp.ncams; ++rc)
     {
-        std::string visTrisFileName = tmpDir + "visTris" + std::to_string(mp.getViewId(rc)) + ".bin";
+        std::string visTrisFilepath = tmpDir + "visTris" + std::to_string(mp.getViewId(rc)) + ".bin";
         StaticVector<int> visTris;
-        loadArrayFromFile<int>(visTris, visTrisFileName);
+        loadArrayFromFile<int>(visTris, visTrisFilepath);
         if(!visTris.empty())
         {
             for(int i = 0; i < visTris.size(); ++i)
@@ -1883,9 +1952,9 @@ void Mesh::computeTrisCams(StaticVector<StaticVector<int>>& trisCams, const mvsU
     t1 = mvsUtils::initEstimate();
     for(int rc = 0; rc < mp.ncams; ++rc)
     {
-        std::string visTrisFileName = tmpDir + "visTris" + std::to_string(mp.getViewId(rc)) + ".bin";
+        std::string visTrisFilepath = tmpDir + "visTris" + std::to_string(mp.getViewId(rc)) + ".bin";
         StaticVector<int> visTris;
-        loadArrayFromFile<int>(visTris, visTrisFileName);
+        loadArrayFromFile<int>(visTris, visTrisFilepath);
         if(!visTris.empty())
         {
             for(int i = 0; i < visTris.size(); ++i)
@@ -2279,7 +2348,7 @@ void Mesh::getLargestConnectedComponentTrisIds(StaticVector<int>& out) const
     }
 }
 
-void Mesh::loadFromObjAscii(const std::string& objAsciiFileName)
+void Mesh::load(const std::string& filepath)
 {
     Assimp::Importer importer;
 
@@ -2294,16 +2363,58 @@ void Mesh::loadFromObjAscii(const std::string& objAsciiFileName)
     normals.clear();
     pointsVisibilities.clear();
 
-    if(!boost::filesystem::exists(objAsciiFileName))
-    {
-        ALICEVISION_THROW_ERROR("Mesh::loadFromObjAscii: no such file: " << objAsciiFileName);
-    }
-
-    unsigned int pFlags = aiProcessPreset_TargetRealtime_MaxQuality & (~aiProcess_SplitLargeMeshes);
-    const aiScene * scene = importer.ReadFile(objAsciiFileName, pFlags);
+    if(!boost::filesystem::exists(filepath))
+    {
+        ALICEVISION_THROW_ERROR("Mesh::load: no such file: " << filepath);
+    }
+
+    // see https://github.com/assimp/assimp/blob/master/include/assimp/postprocess.h#L85
+    const unsigned int pFlags =
+        // If this flag is not specified, no vertices are referenced by more than one face
+        aiProcess_JoinIdenticalVertices |
+        // if a face contain more than 3 vertices, split it in triangles
+        aiProcess_Triangulate |
+        // Removes the node graph and pre-transforms all vertices with the local transformation matrices of their nodes.
+        //aiProcess_PreTransformVertices |
+        // This is intended to get rid of some common exporter errors
+        //// aiProcess_FindInvalidData |
+        // Face normals are shared between all points of a single face,
+        // so a single point can have multiple normals, which forces the library to duplicate vertices in some cases.
+        aiProcess_DropNormals |
+        // This makes sure that all indices are valid
+        //aiProcess_ValidateDataStructure |
+        aiProcess_RemoveComponent |
+        // This step searches all meshes for degenerate primitives and converts them to proper lines or points.
+        // A face is 'degenerate' if one or more of its points are identical.
+        aiProcess_FindDegenerates |
+        // aiProcess_SortByPType needed for aiProcess_FindDegenerates.
+        // This step splits meshes with more than one primitive type in homogeneous sub-meshes
+        // (point and line in different meshes, so we can remove them with AI_CONFIG_PP_SBP_REMOVE).
+        aiProcess_SortByPType |
+        0;
+    importer.SetPropertyInteger(
+        AI_CONFIG_PP_RVC_FLAGS, 
+        aiComponent_NORMALS | 
+        aiComponent_TANGENTS_AND_BITANGENTS
+        // We do not remove texture coords as we need it for texturing,
+        // but it causes vertices duplicates with assimp. This is problematic for mesh post-processing
+        // but we should face this problem only for mesh coming from retopology,
+        // in which case we will only do texturing.
+        //aiComponent_TEXCOORDS
+        );
+
+    // aiProcess_FindDegenerates will convert degenerate triangles.
+    // As we don't want lines and points, we set the AI_CONFIG_PP_SBP_REMOVE.
+    importer.SetPropertyInteger(AI_CONFIG_PP_SBP_REMOVE, aiPrimitiveType_POINT | aiPrimitiveType_LINE);
+
+    // aiProcess_FindDegenerates will also remove very small triangles with a surface area smaller than 10^-6.
+    // As we don't want this extra-behavior, we set the property AI_CONFIG_PP_FD_CHECKAREA to false.
+    importer.SetPropertyBool(AI_CONFIG_PP_FD_CHECKAREA, false);
+
+    const aiScene* scene = importer.ReadFile(filepath, pFlags);
     if (!scene)
     {
-        ALICEVISION_THROW_ERROR("Failed loading mesh from file: " << objAsciiFileName);
+        ALICEVISION_THROW_ERROR("Failed loading mesh from file: " << filepath);
     }
 
     std::list<aiNode *> nodes;
@@ -2416,6 +2527,9 @@ void Mesh::loadFromObjAscii(const std::string& objAsciiFileName)
             nodes.push_back(node->mChildren[idChild]);
         }
     }
+    ALICEVISION_LOG_DEBUG("Vertices: " << pts.size());
+    ALICEVISION_LOG_DEBUG("Triangles: " << tris.size());
+    ALICEVISION_LOG_DEBUG("UVs: " << uvCoords.size());
 }
 
 bool Mesh::getEdgeNeighTrisInterval(Pixel& itr, Pixel& edge, StaticVector<Voxel>& edgesXStat,
diff --git a/src/aliceVision/mesh/Mesh.hpp b/src/aliceVision/mesh/Mesh.hpp
index 2cc52d20c2..063538b2fa 100644
--- a/src/aliceVision/mesh/Mesh.hpp
+++ b/src/aliceVision/mesh/Mesh.hpp
@@ -41,6 +41,22 @@ ALICEVISION_BITMASK(EVisibilityRemappingMethod);
 EVisibilityRemappingMethod EVisibilityRemappingMethod_stringToEnum(const std::string& method);
 std::string EVisibilityRemappingMethod_enumToString(EVisibilityRemappingMethod method);
 
+/**
+ * @brief File type available for exporting mesh
+ */
+enum class EFileType
+{
+    OBJ = 0,
+    FBX,
+    GLTF,
+    STL
+};
+
+EFileType EFileType_stringToEnum(const std::string& filetype);
+std::string EFileType_enumToString(const EFileType filetype);
+std::istream& operator>>(std::istream& in, EFileType& meshFileType);
+std::ostream& operator<<(std::ostream& os, EFileType meshFileType);
+
 
 class Mesh
 {
@@ -148,11 +164,11 @@ class Mesh
     Mesh();
     ~Mesh();
 
-    void saveToObj(const std::string& filename);
+    void save(const std::string& filepath, EFileType filetype = EFileType::OBJ);
 
-    bool loadFromBin(const std::string& binFileName);
-    void saveToBin(const std::string& binFileName);
-    void loadFromObjAscii(const std::string& objAsciiFileName);
+    bool loadFromBin(const std::string& binFilepath);
+    void saveToBin(const std::string& binFilepath);
+    void load(const std::string& filepath);
 
     void addMesh(const Mesh& mesh);
 
@@ -177,7 +193,7 @@ class Mesh
     void getPtsNeighPtsOrdered(StaticVector<StaticVector<int>>& out_ptsNeighTris) const;
 
     void getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& tmpDir, const mvsUtils::MultiViewParams& mp, int rc, int w, int h);
-    void getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& depthMapFileName, const std::string& trisMapFileName,
+    void getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, const std::string& depthMapFilepath, const std::string& trisMapFilepath,
                                                   const mvsUtils::MultiViewParams& mp, int rc, int w, int h);
     void getVisibleTrianglesIndexes(StaticVector<int>& out_visTri, StaticVector<StaticVector<int>>& trisMap,
                                                   StaticVector<float>& depthMap, const mvsUtils::MultiViewParams& mp, int rc, int w,
diff --git a/src/aliceVision/mesh/MeshEnergyOpt.cpp b/src/aliceVision/mesh/MeshEnergyOpt.cpp
index d7f1da5e3c..0acb35a02a 100644
--- a/src/aliceVision/mesh/MeshEnergyOpt.cpp
+++ b/src/aliceVision/mesh/MeshEnergyOpt.cpp
@@ -105,7 +105,7 @@ bool MeshEnergyOpt::optimizeSmooth(float lambda, int niter, StaticVectorBool& pt
         ALICEVISION_LOG_INFO("Optimizing mesh smooth: iteration " << i);
         updateGradientParallel(lambda, LU, RD, ptsCanMove);
         //if(saveDebug)
-        //    saveToObj(folder + "mesh_smoothed_" + std::to_string(i) + ".obj");
+        //    save(folder + "mesh_smoothed_" + std::to_string(i));
     }
 
     return true;
diff --git a/src/aliceVision/mesh/Texturing.cpp b/src/aliceVision/mesh/Texturing.cpp
index 91c9dbb49f..7891b146c7 100644
--- a/src/aliceVision/mesh/Texturing.cpp
+++ b/src/aliceVision/mesh/Texturing.cpp
@@ -24,10 +24,19 @@
 #include <geogram/mesh/mesh_io.h>
 #include <geogram/parameterization/mesh_atlas_maker.h>
 
+#include <geogram/basic/permutation.h>
+#include <geogram/basic/attributes.h>
+#include <geogram/basic/geometry_nd.h>
+#include <geogram/points/kd_tree.h>
+#include <geogram/mesh/mesh_AABB.h>
+#include <geogram/mesh/mesh_reorder.h>
+#include <geogram/mesh/mesh_geometry.h>
+
 #include <boost/algorithm/string/case_conv.hpp> 
 
 #include <assimp/Exporter.hpp>
 #include <assimp/scene.h>
+#include <assimp/postprocess.h>
 
 #include <map>
 #include <set>
@@ -96,6 +105,37 @@ std::string EVisibilityRemappingMethod_enumToString(EVisibilityRemappingMethod m
     throw std::out_of_range("Unrecognized EVisibilityRemappingMethod");
 }
 
+EBumpMappingType EBumpMappingType_stringToEnum(const std::string& type) 
+{
+    if(type == "Height")
+        return EBumpMappingType::Height;
+    if(type == "Normal")
+        return EBumpMappingType::Normal;
+    throw std::out_of_range("Invalid bump mapping type " + type);
+}
+std::string EBumpMappingType_enumToString(EBumpMappingType type) 
+{
+    switch(type)
+    {
+        case EBumpMappingType::Height:
+            return "Height";
+        case EBumpMappingType::Normal:
+            return "Normal";
+    }
+    throw std::out_of_range("Invalid bump mapping type enum");
+}
+std::ostream& operator<<(std::ostream& os, EBumpMappingType bumpMappingType)
+{
+    return os << EBumpMappingType_enumToString(bumpMappingType);
+}
+std::istream& operator>>(std::istream& in, EBumpMappingType& bumpMappingType)
+{
+    std::string token;
+    in >> token;
+    bumpMappingType = EBumpMappingType_stringToEnum(token);
+    return in;
+}
+
 /**
  * @brief Return whether a pixel is contained in or intersected by a 2D triangle.
  * @param[in] triangle the triangle as an array of 3 point2Ds
@@ -698,6 +738,26 @@ void Texturing::generateTexturesSubSet(const mvsUtils::MultiViewParams& mp,
     }
 }
 
+void Texturing::generateNormalAndHeightMaps(const mvsUtils::MultiViewParams& mp, const Mesh& denseMesh,
+                                            const bfs::path& outPath, const mesh::BumpMappingParams& bumpMappingParams)
+{
+    GEO::Mesh geoDenseMesh;
+    toGeoMesh(denseMesh, geoDenseMesh);
+    GEO::compute_normals(geoDenseMesh);
+    GEO::MeshFacetsAABB denseMeshAABB(geoDenseMesh); // warning: mesh_reorder called inside
+
+    GEO::Mesh geoSparseMesh;
+    toGeoMesh(*mesh, geoSparseMesh);
+    GEO::compute_normals(geoSparseMesh);
+
+    mvsUtils::ImagesCache imageCache(&mp, imageIO::EImageColorSpace::NO_CONVERSION);
+    
+    for(size_t atlasID = 0; atlasID < _atlases.size(); ++atlasID)
+        _generateNormalAndHeightMaps(mp, denseMeshAABB, geoSparseMesh, atlasID, imageCache, outPath, bumpMappingParams);
+}
+
+
+
 void Texturing::writeTexture(AccuImage& atlasTexture, const std::size_t atlasID, const boost::filesystem::path &outPath,
                              imageIO::EImageFileType textureFileType, const int level)
 {
@@ -845,13 +905,13 @@ void Texturing::clear()
     mesh = nullptr;
 }
 
-void Texturing::loadOBJWithAtlas(const std::string& filename, bool flipNormals)
+void Texturing::loadWithAtlas(const std::string& filepath, bool flipNormals)
 {
     // Clear internal data
     clear();
     mesh = new Mesh();
     // Load .obj
-    mesh->loadFromObjAscii(filename);
+    mesh->load(filepath);
 
     // Handle normals flipping
     if(flipNormals)
@@ -905,7 +965,7 @@ void Texturing::replaceMesh(const std::string& otherMeshPath, bool flipNormals)
     mesh = nullptr;
     
     // load input obj file
-    loadOBJWithAtlas(otherMeshPath, flipNormals);
+    loadWithAtlas(otherMeshPath, flipNormals);
     // allocate pointsVisibilities for new internal mesh
     mesh->pointsVisibilities = PointsVisibility();
     // remap visibilities from reconstruction onto input mesh
@@ -959,11 +1019,15 @@ void Texturing::unwrap(mvsUtils::MultiViewParams& mp, EUnwrapMethod method)
     }
 }
 
-void Texturing::saveAsOBJ(const bfs::path& dir, const std::string& basename, imageIO::EImageFileType textureFileType)
+void Texturing::saveAs(const bfs::path& dir, const std::string& basename, 
+    EFileType meshFileType, 
+    imageIO::EImageFileType textureFileType,
+    const BumpMappingParams& bumpMappingParams)
 {
-    ALICEVISION_LOG_INFO("Writing obj and mtl file.");
+    const std::string meshFileTypeStr = EFileType_enumToString(meshFileType);
+    const std::string filepath = (dir / (basename + "." + meshFileTypeStr)).string();
 
-    const std::string objFilename = (dir / (basename + ".obj")).string();
+    ALICEVISION_LOG_INFO("Save " << filepath << " mesh file");
 
     if (_atlases.empty())
     {
@@ -1002,7 +1066,33 @@ void Texturing::saveAsOBJ(const bfs::path& dir, const std::string& basename, ima
         scene.mMaterials[atlasId]->AddProperty(&shininess, 1, AI_MATKEY_SHININESS);
         scene.mMaterials[atlasId]->AddProperty(&texFile, AI_MATKEY_TEXTURE_DIFFUSE(0));
         scene.mMaterials[atlasId]->AddProperty(&texName, AI_MATKEY_NAME);
+
+        // Color Mapping
+        if(textureFileType != imageIO::EImageFileType::NONE)
+        {
+            const aiString texFile(texturePath);
+            scene.mMaterials[atlasId]->AddProperty(&texFile, AI_MATKEY_TEXTURE_DIFFUSE(0));
+        }
+
+        // Displacement Mapping
+        if(bumpMappingParams.displacementFileType != imageIO::EImageFileType::NONE)
+        {
+            const aiString texFileHeightMap("Displacement_" + std::to_string(textureId) + "." +EImageFileType_enumToString(bumpMappingParams.bumpMappingFileType));
+            scene.mMaterials[atlasId]->AddProperty(&texFileHeightMap, AI_MATKEY_TEXTURE_DISPLACEMENT(0));
+        }
         
+        // Bump Mapping
+        if(bumpMappingParams.bumpType == EBumpMappingType::Normal && bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::NONE)
+        {
+            const aiString texFileNormalMap("Normal_" + std::to_string(textureId) + "." + EImageFileType_enumToString(bumpMappingParams.bumpMappingFileType));
+            scene.mMaterials[atlasId]->AddProperty(&texFileNormalMap, AI_MATKEY_TEXTURE_NORMALS(0));
+        }
+        else if(bumpMappingParams.bumpType == EBumpMappingType::Height && bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::NONE)
+        {
+            const aiString texFileHeightMap("Bump_" + std::to_string(textureId) + "." + EImageFileType_enumToString(bumpMappingParams.displacementFileType));
+            scene.mMaterials[atlasId]->AddProperty(&texFileHeightMap, AI_MATKEY_TEXTURE_HEIGHT(0));
+        }
+
         scene.mRootNode->mMeshes[atlasId] = atlasId;
         scene.mMeshes[atlasId] = new aiMesh;
         aiMesh * aimesh = scene.mMeshes[atlasId];
@@ -1068,8 +1158,339 @@ void Texturing::saveAsOBJ(const bfs::path& dir, const std::string& basename, ima
         }
     }
 
+    std::string formatId = meshFileTypeStr;
+    unsigned int pPreprocessing = 0u;
+    // If gltf, use gltf 2.0
+    if(meshFileType == EFileType::GLTF)
+    {
+        formatId = "gltf2";
+        // Flip UVs when exporting (issue with UV origin for gltf2)
+        // https://github.com/around-media/ue4-custom-prompto/commit/044dbad90fc2172f4c5a8b67c779b80ceace5e1e
+        pPreprocessing |= aiPostProcessSteps::aiProcess_FlipUVs | aiProcess_GenNormals;
+    }
+
     Assimp::Exporter exporter;
-    exporter.Export(&scene, "obj", objFilename);
+    exporter.Export(&scene, formatId, filepath, pPreprocessing);
+
+    ALICEVISION_LOG_INFO("Save mesh to " << meshFileTypeStr << " done.");
+}
+
+
+inline GEO::vec3 mesh_facet_interpolate_normal_at_point(const GEO::Mesh& mesh, GEO::index_t f, const GEO::vec3& p)
+{
+    const GEO::index_t v0 = mesh.facets.vertex(f, 0);
+    const GEO::index_t v1 = mesh.facets.vertex(f, 1);
+    const GEO::index_t v2 = mesh.facets.vertex(f, 2);
+
+    const GEO::vec3 p0 = mesh.vertices.point(v0);
+    const GEO::vec3 p1 = mesh.vertices.point(v1);
+    const GEO::vec3 p2 = mesh.vertices.point(v2);
+
+    const GEO::vec3 n0 = GEO::normalize(GEO::Geom::mesh_vertex_normal(mesh, v0));
+    const GEO::vec3 n1 = GEO::normalize(GEO::Geom::mesh_vertex_normal(mesh, v1));
+    const GEO::vec3 n2 = GEO::normalize(GEO::Geom::mesh_vertex_normal(mesh, v2));
+
+    GEO::vec3 barycCoords;
+    GEO::vec3 closestPoint;
+    GEO::Geom::point_triangle_squared_distance<GEO::vec3>(p, p0, p1, p2, closestPoint, barycCoords.x, barycCoords.y,
+                                                          barycCoords.z);
+
+    const GEO::vec3 n = barycCoords.x * n0 + barycCoords.y * n1 + barycCoords.z * n2;
+
+    return GEO::normalize(n);
+}
+
+inline GEO::vec3 mesh_facet_interpolate_normal_at_point(const StaticVector<Point3d>& ptsNormals, const Mesh& mesh,
+                                                        GEO::index_t f, const GEO::vec3& p)
+{
+    const GEO::index_t v0 = (mesh.tris)[f].v[0];
+    const GEO::index_t v1 = (mesh.tris)[f].v[1];
+    const GEO::index_t v2 = (mesh.tris)[f].v[2];
+
+    const GEO::vec3 p0 ((mesh.pts)[v0].x, (mesh.pts)[v0].y, (mesh.pts)[v0].z);
+    const GEO::vec3 p1 ((mesh.pts)[v1].x, (mesh.pts)[v1].y, (mesh.pts)[v1].z);
+    const GEO::vec3 p2 ((mesh.pts)[v2].x, (mesh.pts)[v2].y, (mesh.pts)[v2].z);
+
+    const GEO::vec3 n0 (ptsNormals[v0].x, ptsNormals[v0].y, ptsNormals[v0].z);
+    const GEO::vec3 n1 (ptsNormals[v1].x, ptsNormals[v1].y, ptsNormals[v1].z);
+    const GEO::vec3 n2 (ptsNormals[v2].x, ptsNormals[v2].y, ptsNormals[v2].z);
+
+    GEO::vec3 barycCoords;
+    GEO::vec3 closestPoint;
+    GEO::Geom::point_triangle_squared_distance<GEO::vec3>(p, p0, p1, p2, closestPoint, barycCoords.x, barycCoords.y,
+                                                          barycCoords.z);
+
+    const GEO::vec3 n = barycCoords.x * n0 + barycCoords.y * n1 + barycCoords.z * n2;
+
+    return GEO::normalize(n);
+}
+
+template <class T, GEO::index_t DIM>
+inline Eigen::Matrix<T, DIM, 1> toEigen(const GEO::vecng<DIM, T>& v)
+{
+    return Eigen::Matrix<T, DIM, 1>(v.data());
+}
+
+/**
+ * @brief Compute a transformation matrix to convert coordinates in world space coordinates into the triangle space.
+ * The triangle space is define by the Z-axis as the normal of the triangle,
+ * the X-axis aligned with the horizontal line in the texture file (using texture/UV coordinates).
+ *
+ * @param[in] mesh: input mesh
+ * @param[in] f: facet/triangle index
+ * @param[in] triPts: UV Coordinates
+ * @return Rotation matrix to convert from world space coordinates in the triangle space
+ */
+inline Eigen::Matrix3d computeTriangleTransform(const Mesh& mesh, int f, const Point2d* triPts)
+{
+    const Eigen::Vector3d p0 = toEigen((mesh.pts)[(mesh.tris)[f].v[0]]);
+    const Eigen::Vector3d p1 = toEigen((mesh.pts)[(mesh.tris)[f].v[1]]);
+    const Eigen::Vector3d p2 = toEigen((mesh.pts)[(mesh.tris)[f].v[2]]);
+
+    const Eigen::Vector3d tX = (p1 - p0).normalized();                       // edge0 => local triangle X-axis
+    const Eigen::Vector3d N = tX.cross((p2 - p0).normalized()).normalized(); // cross(edge0, edge1) => Z-axis
+
+    // Correct triangle X-axis to be align with X-axis in the texture
+    const GEO::vec2 t0 = GEO::vec2(triPts[0].m);
+    const GEO::vec2 t1 = GEO::vec2(triPts[1].m);
+    const GEO::vec2 tV = GEO::normalize(t1 - t0);
+    const GEO::vec2 origNormal(1.0, 0.0); // X-axis in the texture
+    const double tAngle = GEO::Geom::angle(tV, origNormal);
+    Eigen::Matrix3d transform(Eigen::AngleAxisd(tAngle, N).toRotationMatrix());
+    // Rotate triangle v0v1 axis around Z-axis, to get a X axis aligned with the 2d texture
+    Eigen::Vector3d X = (transform * tX).normalized();
+
+    const Eigen::Vector3d Y = N.cross(X).normalized(); // Y-axis
+
+    Eigen::Matrix3d m;
+    m.col(0) = X;
+    m.col(1) = Y;
+    m.col(2) = N;
+    // const Eigen::Matrix3d mInv = m.inverse();
+    const Eigen::Matrix3d mT = m.transpose();
+
+    return mT;
+}
+
+inline void computeNormalHeight(const GEO::Mesh& mesh, double orientation, double t, GEO::index_t f,
+                                const Eigen::Matrix3d& m, const GEO::vec3& q, const GEO::vec3& qA, const GEO::vec3& qB,
+                                float& out_height, Color& out_normal)
+{
+    GEO::vec3 intersectionPoint = t * qB + (1.0 - t) * qA;
+    out_height = q.distance(intersectionPoint) * orientation;
+
+    // Use facet normal
+    // GEO::vec3 denseMeshNormal_f = normalize(GEO::Geom::mesh_facet_normal(mesh, f));
+    // Use per pixel normal using weighted interpolation of the facet vertex normals
+    const GEO::vec3 denseMeshNormal = mesh_facet_interpolate_normal_at_point(mesh, f, intersectionPoint);
+
+    Eigen::Vector3d dNormal = m * toEigen(denseMeshNormal);
+    dNormal.normalize();
+    out_normal = Color(dNormal(0), dNormal(1), dNormal(2));
+}
+
+void Texturing::_generateNormalAndHeightMaps(const mvsUtils::MultiViewParams& mp,
+                                             const GEO::MeshFacetsAABB& denseMeshAABB, const GEO::Mesh& sparseMesh,
+                                             size_t atlasID, mvsUtils::ImagesCache& imageCache,
+                                             const bfs::path& outPath, const mesh::BumpMappingParams& bumpMappingParams)
+{
+    ALICEVISION_LOG_INFO("Generating Height and Normal Maps for atlas " << atlasID + 1 << "/" << _atlases.size() << " ("
+                                                                        << _atlases[atlasID].size() << " triangles).");
+
+    std::vector<Color> normalMap(texParams.textureSide * texParams.textureSide);
+    std::vector<float> heightMap(texParams.textureSide * texParams.textureSide);
+    const auto& triangles = _atlases[atlasID];
+
+    // iterate over atlas' triangles
+#pragma omp parallel for
+    for(int ti = 0; ti < triangles.size(); ++ti)
+    {
+        const unsigned int triangleId = triangles[ti];
+        //  const Point3d __triangleNormal_ = me->computeTriangleNormal(triangleId).normalize();
+        //  const GEO::vec3 __triangleNormal(__triangleNormal_.x, __triangleNormal_.y, __triangleNormal_.z);
+        const double minEdgeLength = mesh->computeTriangleMinEdgeLength(triangleId);
+        // const GEO::vec3 scaledTriangleNormal = triangleNormal * minEdgeLength;
+
+        // retrieve triangle 3D and UV coordinates
+        Point2d triPixs[3];
+        Point3d triPts[3];
+        auto& triangleUvIds = mesh->trisUvIds[triangleId];
+
+        // compute the Bottom-Left minima of the current UDIM for [0,1] range remapping
+        Point2d udimBL;
+        StaticVector<Point2d>& uvCoords = mesh->uvCoords;
+        udimBL.x = std::floor(std::min(std::min(uvCoords[triangleUvIds[0]].x, uvCoords[triangleUvIds[1]].x),uvCoords[triangleUvIds[2]].x));
+        udimBL.y = std::floor(std::min(std::min(uvCoords[triangleUvIds[0]].y, uvCoords[triangleUvIds[1]].y),uvCoords[triangleUvIds[2]].y));
+
+        for(int k = 0; k < 3; k++)
+        {
+            const int pointIndex = (mesh->tris)[triangleId].v[k];
+            triPts[k] = (mesh->pts)[pointIndex]; // 3D coordinates
+            const int uvPointIndex = triangleUvIds.m[k];
+
+            Point2d uv = uvCoords[uvPointIndex];
+            // UDIM: remap coordinates between [0,1]
+            uv = uv - udimBL;
+
+            triPixs[k] = uv * texParams.textureSide; // UV coordinates
+        }
+
+        // compute triangle bounding box in pixel indexes
+        // min values: floor(value)
+        // max values: ceil(value)
+        Pixel LU, RD;
+        LU.x = static_cast<int>(std::floor(std::min(std::min(triPixs[0].x, triPixs[1].x), triPixs[2].x)));
+        LU.y = static_cast<int>(std::floor(std::min(std::min(triPixs[0].y, triPixs[1].y), triPixs[2].y)));
+        RD.x = static_cast<int>(std::ceil(std::max(std::max(triPixs[0].x, triPixs[1].x), triPixs[2].x)));
+        RD.y = static_cast<int>(std::ceil(std::max(std::max(triPixs[0].y, triPixs[1].y), triPixs[2].y)));
+
+        // sanity check: clamp values to [0; textureSide]
+        int texSide = static_cast<int>(texParams.textureSide);
+        LU.x = clamp(LU.x, 0, texSide);
+        LU.y = clamp(LU.y, 0, texSide);
+        RD.x = clamp(RD.x, 0, texSide);
+        RD.y = clamp(RD.y, 0, texSide);
+
+        const Eigen::Matrix3d worldToTriangleMatrix = computeTriangleTransform(*mesh, triangleId, triPixs);
+        // const Point3d triangleNormal = me->computeTriangleNormal(triangleId);
+
+        // iterate over bounding box's pixels
+        for(int y = LU.y; y < RD.y; ++y)
+        {
+            for(int x = LU.x; x < RD.x; ++x)
+            {
+                Pixel pix(x, y); // top-left corner of the pixel
+                Point2d barycCoords;
+
+                // test if the pixel is inside triangle
+                // and retrieve its barycentric coordinates
+                if(!isPixelInTriangle(triPixs, pix, barycCoords))
+                {
+                    continue;
+                }
+
+                // remap 'y' to image coordinates system (inverted Y axis)
+                const unsigned int y_ = (texParams.textureSide - 1) - y;
+                // 1D pixel index
+                unsigned int xyoffset = y_ * texParams.textureSide + x;
+                // get 3D coordinates
+                // Point3d pt3d = barycentricToCartesian(triPts, Point2d(barycCoords.z, barycCoords.y));
+                Point3d pt3d = barycentricToCartesian(triPts, barycCoords);
+                GEO::vec3 q(pt3d.x, pt3d.y, pt3d.z);
+
+                // Texel normal (weighted normal from the 3 vertices normals), instead of face normal for better
+                // transitions (reduce seams)
+                const GEO::vec3 triangleNormal_p = mesh_facet_interpolate_normal_at_point(sparseMesh, triangleId, q);
+                // const GEO::vec3 triangleNormal_p = GEO::vec3(triangleNormal.m); // to use the triangle normal instead
+                const GEO::vec3 scaledTriangleNormal = triangleNormal_p * minEdgeLength * 10;
+
+                const double epsilon = 0.00001;
+                GEO::vec3 qA1 = q - (scaledTriangleNormal * epsilon);
+                GEO::vec3 qB1 = q + scaledTriangleNormal;
+                double t = 0.0;
+                GEO::index_t f = 0.0;
+                bool intersection = denseMeshAABB.segment_nearest_intersection(qA1, qB1, t, f);
+                if(intersection)
+                {
+                    computeNormalHeight(*denseMeshAABB.mesh(), 1.0, t, f, worldToTriangleMatrix, q, qA1, qB1,
+                                        heightMap[xyoffset], normalMap[xyoffset]);
+                }
+                else
+                {
+                    GEO::vec3 qA2 = q + (scaledTriangleNormal * epsilon);
+                    GEO::vec3 qB2 = q - scaledTriangleNormal;
+                    bool intersection = denseMeshAABB.segment_nearest_intersection(qA2, qB2, t, f);
+                    if(intersection)
+                    {
+                        computeNormalHeight(*denseMeshAABB.mesh(), -1.0, t, f, worldToTriangleMatrix, q, qA2, qB2,
+                                            heightMap[xyoffset], normalMap[xyoffset]);
+                    }
+                    else
+                    {
+                        heightMap[xyoffset] = 0.0f;
+                        normalMap[xyoffset] = Color(0.0f, 0.0f, 0.0f);
+                    }
+                }
+            }
+        }
+    }
+
+    // Save Normal Map
+    if(bumpMappingParams.bumpType == EBumpMappingType::Normal && bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::NONE)
+    {
+        unsigned int outTextureSide = texParams.textureSide;
+        // downscale texture if required
+        if(texParams.downscale > 1)
+        {
+            ALICEVISION_LOG_INFO("Downscaling normal map (" << texParams.downscale << "x).");
+            std::vector<Color> resizedBuffer;
+            outTextureSide = texParams.textureSide / texParams.downscale;
+            // use nearest-neighbor interpolation to avoid meaningless interpolation of normals on edges.
+            const std::string interpolation = "box";
+            imageAlgo::resizeImage(texParams.textureSide, texParams.textureSide, texParams.downscale, normalMap,
+                                   resizedBuffer, interpolation);
+
+            std::swap(resizedBuffer, normalMap);
+        }
+
+        // X: -1 to +1 : Red : 0 to 255
+        // Y: -1 to +1 : Green : 0 to 255
+        // Z: 0 to -1 : Blue : 128 to 255 OR 0 to 255 (like Blender)
+        for(unsigned int i = 0; i < normalMap.size(); ++i)
+            // normalMap[i] = Color(normalMap[i].r * 0.5 + 0.5, normalMap[i].g * 0.5 + 0.5, normalMap[i].b); // B:
+            // 0:+1 => 0-255
+            normalMap[i] = Color(normalMap[i].r * 0.5 + 0.5, normalMap[i].g * 0.5 + 0.5,
+                                    normalMap[i].b * 0.5 + 0.5); // B: -1:+1 => 0-255 which means 0:+1 => 128-255
+
+        const std::string name = "Normal_" + std::to_string(1001 + atlasID) + "." + EImageFileType_enumToString(bumpMappingParams.bumpMappingFileType);
+        bfs::path normalMapPath = outPath / name;
+        ALICEVISION_LOG_INFO("Writing normal map: " << normalMapPath.string());
+
+        imageIO::OutputFileColorSpace outputColorSpace(imageIO::EImageColorSpace::NO_CONVERSION,imageIO::EImageColorSpace::NO_CONVERSION);
+        imageIO::writeImage(normalMapPath.string(), outTextureSide, outTextureSide, normalMap, imageIO::EImageQuality::OPTIMIZED, outputColorSpace);
+    }
+
+    // Save Height Maps
+    if(bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::NONE || bumpMappingParams.displacementFileType != imageIO::EImageFileType::NONE)
+    {
+        unsigned int outTextureSide = texParams.textureSide;
+        if(texParams.downscale > 1)
+        {
+            ALICEVISION_LOG_INFO("Downscaling height map (" << texParams.downscale << "x).");
+            std::vector<float> resizedBuffer;
+            outTextureSide = texParams.textureSide / texParams.downscale;
+            imageAlgo::resizeImage(texParams.textureSide, texParams.textureSide, texParams.downscale, heightMap,
+                                   resizedBuffer);
+            std::swap(resizedBuffer, heightMap);
+        }
+
+        // Height maps are only in .EXR at this time, so this will never be executed.
+        // 
+        //if(bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::EXR)
+        //{
+        //    // Y: [-1, 0, +1] => [0, 128, 255]
+        //    for(unsigned int i = 0; i < heightMap.size(); ++i)
+        //        heightMap[i] = heightMap[i] * 0.5 + 0.5;
+        //}
+
+        // Save Bump Map
+        imageIO::OutputFileColorSpace outputColorSpace(imageIO::EImageColorSpace::AUTO);
+        if(bumpMappingParams.bumpType == EBumpMappingType::Height)
+        {
+            const std::string bumpName = "Bump_" + std::to_string(1001 + atlasID) + "." + EImageFileType_enumToString(bumpMappingParams.bumpMappingFileType);
+            bfs::path bumpMapPath = outPath / bumpName;
+            ALICEVISION_LOG_INFO("Writing bump map: " << bumpMapPath);
+            imageIO::writeImage(bumpMapPath.string(), outTextureSide, outTextureSide, heightMap, imageIO::EImageQuality::OPTIMIZED, outputColorSpace);
+        }
+        // Save Displacement Map
+        if(bumpMappingParams.displacementFileType != imageIO::EImageFileType::NONE)
+        {
+            const std::string dispName = "Displacement_" + std::to_string(1001 + atlasID) + "." + EImageFileType_enumToString(bumpMappingParams.displacementFileType);
+            bfs::path dispMapPath = outPath / dispName;
+            ALICEVISION_LOG_INFO("Writing displacement map: " << dispMapPath);
+            imageIO::writeImage(dispMapPath.string(), outTextureSide, outTextureSide, heightMap, imageIO::EImageQuality::OPTIMIZED, outputColorSpace);
+        }
+    }
 }
 
 } // namespace mesh
diff --git a/src/aliceVision/mesh/Texturing.hpp b/src/aliceVision/mesh/Texturing.hpp
index bfe9c9d372..33ce081651 100644
--- a/src/aliceVision/mesh/Texturing.hpp
+++ b/src/aliceVision/mesh/Texturing.hpp
@@ -20,6 +20,11 @@
 
 namespace bfs = boost::filesystem;
 
+namespace GEO {
+    class MeshFacetsAABB;
+    class Mesh;
+}
+
 namespace aliceVision {
 namespace mesh {
 
@@ -46,6 +51,24 @@ EUnwrapMethod EUnwrapMethod_stringToEnum(const std::string& method);
  */
 std::string EUnwrapMethod_enumToString(EUnwrapMethod method);
 
+enum class EBumpMappingType
+{
+    Height = 0,
+    Normal
+};
+EBumpMappingType EBumpMappingType_stringToEnum(const std::string& type);
+std::string EBumpMappingType_enumToString(EBumpMappingType type);
+std::istream& operator>>(std::istream& in, EBumpMappingType& meshFileType);
+std::ostream& operator<<(std::ostream& os, EBumpMappingType meshFileType);
+
+
+struct BumpMappingParams
+{
+    imageIO::EImageFileType bumpMappingFileType = imageIO::EImageFileType::NONE;
+    imageIO::EImageFileType displacementFileType = imageIO::EImageFileType::NONE;
+
+    EBumpMappingType bumpType = EBumpMappingType::Normal;
+};
 
 struct TexturingParams
 {
@@ -64,6 +87,7 @@ struct TexturingParams
     double bestScoreThreshold = 0.1; //< 0.0 to disable filtering based on threshold to relative best score
     double angleHardThreshold = 90.0; //< 0.0 to disable angle hard threshold filtering
 
+    imageIO::EImageFileType textureFileType = imageIO::EImageFileType::NONE;
     imageIO::EImageColorSpace processColorspace = imageIO::EImageColorSpace::SRGB; // colorspace for the texturing internal computation
     mvsUtils::ImagesCache::ECorrectEV correctEV{mvsUtils::ImagesCache::ECorrectEV::NO_CORRECTION};
 
@@ -92,7 +116,7 @@ struct Texturing
     void clear();
 
     /// Load a mesh from a .obj file and initialize internal structures
-    void loadOBJWithAtlas(const std::string& filename, bool flipNormals=false);
+    void loadWithAtlas(const std::string& filepath, bool flipNormals=false);
 
     /**
      * @brief Remap visibilities
@@ -173,12 +197,22 @@ struct Texturing
                          const std::vector<size_t>& atlasIDs, mvsUtils::ImagesCache& imageCache,
                          const bfs::path &outPath, imageIO::EImageFileType textureFileType = imageIO::EImageFileType::PNG);
 
+    void generateNormalAndHeightMaps(const mvsUtils::MultiViewParams& mp, const Mesh& denseMesh,
+                                     const bfs::path& outPath, const mesh::BumpMappingParams& bumpMappingParams);
+
+    void _generateNormalAndHeightMaps(const mvsUtils::MultiViewParams& mp, const GEO::MeshFacetsAABB& denseMeshAABB,
+                                      const GEO::Mesh& sparseMesh, size_t atlasID, mvsUtils::ImagesCache& imageCache,
+                                      const bfs::path& outPath, const mesh::BumpMappingParams& bumpMappingParams);
+
     ///Fill holes and write texture files for the given texture atlas
     void writeTexture(AccuImage& atlasTexture, const std::size_t atlasID, const bfs::path& outPath,
                       imageIO::EImageFileType textureFileType, const int level);
 
     /// Save textured mesh as an OBJ + MTL file
-    void saveAsOBJ(const bfs::path& dir, const std::string& basename, imageIO::EImageFileType textureFileType = imageIO::EImageFileType::PNG);
+    void saveAs(const bfs::path& dir, const std::string& basename,
+                aliceVision::mesh::EFileType meshFileType = aliceVision::mesh::EFileType::OBJ,
+                imageIO::EImageFileType textureFileType = imageIO::EImageFileType::EXR,
+                const BumpMappingParams& bumpMappingParams = BumpMappingParams());
 };
 
 } // namespace mesh
diff --git a/src/aliceVision/mesh/meshPostProcessing.cpp b/src/aliceVision/mesh/meshPostProcessing.cpp
index 31f158f50a..6ec1a6e07a 100644
--- a/src/aliceVision/mesh/meshPostProcessing.cpp
+++ b/src/aliceVision/mesh/meshPostProcessing.cpp
@@ -28,7 +28,7 @@ void meshPostProcessing(Mesh*& inout_mesh, StaticVector<StaticVector<int>>& inou
     bool exportDebug = (float)mp.userParams.get<bool>("delaunaycut.exportDebugGC", false);
 
     if(exportDebug)
-        inout_mesh->saveToObj(debugFolderName + "rawGraphCut.obj");
+        inout_mesh->save(debugFolderName + "rawGraphCut");
 
     // copy ptsCams
     {
@@ -74,7 +74,7 @@ void meshPostProcessing(Mesh*& inout_mesh, StaticVector<StaticVector<int>>& inou
         meOpt.cleanMesh(10);
 
         if(exportDebug)
-            meOpt.saveToObj(debugFolderName + "MeshClean.obj");
+            meOpt.save(debugFolderName + "MeshClean");
 
         /////////////////////////////
         {
@@ -128,7 +128,7 @@ void meshPostProcessing(Mesh*& inout_mesh, StaticVector<StaticVector<int>>& inou
             meOpt.optimizeSmooth(lambda, smoothNIter, ptsCanMove);
 
             if(exportDebug)
-                meOpt.saveToObj(debugFolderName + "mesh_smoothed.obj");
+                meOpt.save(debugFolderName + "mesh_smoothed");
         }
 
         meOpt.deallocateCleaningAttributes();
diff --git a/src/aliceVision/mvsData/Point3d.hpp b/src/aliceVision/mvsData/Point3d.hpp
index 52273e4716..a3204da05a 100644
--- a/src/aliceVision/mvsData/Point3d.hpp
+++ b/src/aliceVision/mvsData/Point3d.hpp
@@ -10,6 +10,14 @@
 #include <cmath>
 #include <ostream>
 
+#include <Eigen/Core>
+
+// #include <geogram/basic/vecg.h>
+namespace GEO {
+    template <unsigned int DIM, class T>
+    class vecng;
+};
+
 namespace aliceVision {
 
 class Point3d
@@ -128,6 +136,12 @@ class Point3d
         return x * x + y * y + z * z;
     }
 
+    template <class T>
+    operator GEO::vecng<3, T>() const
+    {
+      return GEO::vecng<3, T>(x, y, z);
+    }
+
     friend double dist(const Point3d& p1, const Point3d& p2);
 
     friend double dot(const Point3d& p1, const Point3d& p2);
@@ -186,4 +200,9 @@ inline std::ostream& operator<<(std::ostream& stream, const Point3d& p)
     return stream;
 }
 
+inline Eigen::Matrix<double, 3, 1> toEigen(const Point3d& v)
+{
+  return Eigen::Matrix<double, 3, 1>(v.m);
+}
+
 } // namespace aliceVision
diff --git a/src/aliceVision/mvsData/imageIO.cpp b/src/aliceVision/mvsData/imageIO.cpp
index e970617b28..d5975c5fef 100644
--- a/src/aliceVision/mvsData/imageIO.cpp
+++ b/src/aliceVision/mvsData/imageIO.cpp
@@ -132,6 +132,7 @@ std::string EImageFileType_enumToString(const EImageFileType imageFileType)
     case EImageFileType::PNG:   return "png";
     case EImageFileType::TIFF:  return "tif";
     case EImageFileType::EXR:   return "exr";
+    case EImageFileType::NONE:  return "none";
   }
   throw std::out_of_range("Invalid EImageType enum");
 }
diff --git a/src/aliceVision/mvsData/imageIO.hpp b/src/aliceVision/mvsData/imageIO.hpp
index 032df29211..e255507d1e 100644
--- a/src/aliceVision/mvsData/imageIO.hpp
+++ b/src/aliceVision/mvsData/imageIO.hpp
@@ -28,7 +28,8 @@ enum class EImageFileType
   JPEG,
   PNG,
   TIFF,
-  EXR
+  EXR,
+  NONE
 };
 
 /**
diff --git a/src/software/convert/CMakeLists.txt b/src/software/convert/CMakeLists.txt
index 28b7f082ad..142d8a7deb 100644
--- a/src/software/convert/CMakeLists.txt
+++ b/src/software/convert/CMakeLists.txt
@@ -62,7 +62,7 @@ alicevision_add_software(aliceVision_convertMesh
   FOLDER ${FOLDER_SOFTWARE_CONVERT}
   LINKS aliceVision_system
         aliceVision_numeric
-	Geogram::geogram
+	aliceVision_mesh
         ${Boost_LIBRARIES}
 )
 
diff --git a/src/software/convert/main_convertMesh.cpp b/src/software/convert/main_convertMesh.cpp
index 9928100774..2303c1f6c1 100644
--- a/src/software/convert/main_convertMesh.cpp
+++ b/src/software/convert/main_convertMesh.cpp
@@ -8,16 +8,12 @@
 #include <aliceVision/system/cmdline.hpp>
 #include <aliceVision/system/main.hpp>
 #include <aliceVision/system/Timer.hpp>
+#include <aliceVision/mesh/Texturing.hpp>
+#include <aliceVision/mesh/Mesh.hpp>
 
 #include <boost/program_options.hpp>
 #include <boost/filesystem.hpp>
 
-#include <geogram/mesh/mesh.h>
-#include <geogram/mesh/mesh_io.h>
-
-#include <geogram/basic/command_line.h>
-#include <geogram/basic/command_line_args.h>
-
 #include <iostream>
 #include <fstream>
 #include <ostream>
@@ -58,14 +54,12 @@ int aliceVision_main(int argc, char** argv)
       ("output,o", po::value<std::string>(&outputFilePath)->default_value(outputFilePath),
         "Output file path for the new mesh file (*.obj, *.mesh, *.meshb, *.ply, *.off, *.stl)");
 
-    po::options_description optionalParams("Optional parameters");
-
     po::options_description logParams("Log parameters");
     logParams.add_options()
       ("verboseLevel,v", po::value<std::string>(&verboseLevel)->default_value(verboseLevel),
         "verbosity level (fatal, error, warning, info, debug, trace).");
 
-    allParams.add(requiredParams).add(optionalParams).add(logParams);
+    allParams.add(requiredParams).add(logParams);
 
     po::variables_map vm;
     try
@@ -126,31 +120,31 @@ int aliceVision_main(int argc, char** argv)
         }
     }
 
-    GEO::initialize();
-    GEO::CmdLine::import_arg_group("standard");
-    GEO::CmdLine::import_arg_group("algo");
-
-    ALICEVISION_LOG_INFO("Geogram initialized.");
-
-    GEO::Mesh inputMesh;
-
     // load input mesh
-    if(!GEO::mesh_load(inputMeshPath, inputMesh))
+    mesh::Texturing texturing;
+    texturing.loadWithAtlas(inputMeshPath);
+    mesh::Mesh* inputMesh = texturing.mesh;
+
+    if(!inputMesh)
     {
-        ALICEVISION_LOG_ERROR("Failed to load mesh file: \"" << inputMeshPath << "\".");
+        ALICEVISION_LOG_ERROR("Unable to read input mesh from the file: " << inputMeshPath);
         return EXIT_FAILURE;
     }
 
-    // save output mesh
-    ALICEVISION_LOG_INFO("Convert mesh.");
-    if(!GEO::mesh_save(inputMesh, outputFilePath))
+    if(inputMesh->pts.empty() || inputMesh->tris.empty())
     {
-        ALICEVISION_LOG_ERROR("Failed to save mesh file: \"" << outputFilePath << "\".");
+        ALICEVISION_LOG_ERROR("Error: empty mesh from the file " << inputMeshPath);
+        ALICEVISION_LOG_ERROR("Input mesh: " << inputMesh->pts.size() << " vertices and " << inputMesh->tris.size()
+                                             << " facets.");
         return EXIT_FAILURE;
     }
 
-    ALICEVISION_LOG_INFO("Mesh file: \"" << outputFilePath << "\" saved.");
+    // save output mesh
+    ALICEVISION_LOG_INFO("Convert mesh.");
+    mesh::EFileType outputMeshFileType = mesh::EFileType_stringToEnum(fs::path(outputFilePath).extension().string().substr(1));
+    inputMesh->save(outputFilePath, outputMeshFileType);
+
     ALICEVISION_LOG_INFO("Task done in (s): " + std::to_string(timer.elapsed()));
 
     return EXIT_SUCCESS;
-}
\ No newline at end of file
+}
diff --git a/src/software/pipeline/main_meshFiltering.cpp b/src/software/pipeline/main_meshFiltering.cpp
index e4bad38bd3..debaa8183b 100644
--- a/src/software/pipeline/main_meshFiltering.cpp
+++ b/src/software/pipeline/main_meshFiltering.cpp
@@ -101,6 +101,7 @@ int aliceVision_main(int argc, char* argv[])
     std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel());
     std::string inputMeshPath;
     std::string outputMeshPath;
+    aliceVision::mesh::EFileType outputMeshFileType;
 
     bool keepLargestMeshOnly = false;
 
@@ -123,12 +124,14 @@ int aliceVision_main(int argc, char* argv[])
     po::options_description requiredParams("Required parameters");
     requiredParams.add_options()
         ("inputMesh,i", po::value<std::string>(&inputMeshPath)->required(),
-            "Input Mesh (OBJ file format).")
+            "Input Mesh")
         ("outputMesh,o", po::value<std::string>(&outputMeshPath)->required(),
-            "Output mesh (OBJ file format).");
+            "Output mesh");
 
     po::options_description optionalParams("Optional parameters");
     optionalParams.add_options()
+        ("outputMeshFileType", po::value<aliceVision::mesh::EFileType>(&outputMeshFileType)->default_value(aliceVision::mesh::EFileType::GLTF),
+            "output mesh file type")
         ("keepLargestMeshOnly", po::value<bool>(&keepLargestMeshOnly)->default_value(keepLargestMeshOnly),
             "Keep only the largest connected triangles group.")
         ("smoothingSubset",po::value<std::string>(&smoothingSubsetTypeName)->default_value(smoothingSubsetTypeName),
@@ -199,7 +202,7 @@ int aliceVision_main(int argc, char* argv[])
         bfs::create_directory(outDirectory);
 
     mesh::Texturing texturing;
-    texturing.loadOBJWithAtlas(inputMeshPath);
+    texturing.loadWithAtlas(inputMeshPath);
     mesh::Mesh* mesh = texturing.mesh;
 
     if(!mesh)
@@ -310,7 +313,7 @@ int aliceVision_main(int argc, char* argv[])
     ALICEVISION_LOG_INFO("Save mesh.");
 
     // Save output mesh
-    outMesh.saveToObj(outputMeshPath);
+    outMesh.save(outputMeshPath, outputMeshFileType);
 
     ALICEVISION_LOG_INFO("Mesh file: \"" << outputMeshPath << "\" saved.");
 
diff --git a/src/software/pipeline/main_meshMasking.cpp b/src/software/pipeline/main_meshMasking.cpp
index 261edea021..83d838fee1 100644
--- a/src/software/pipeline/main_meshMasking.cpp
+++ b/src/software/pipeline/main_meshMasking.cpp
@@ -378,7 +378,8 @@ void meshMasking(
     const bool invert,
     const bool smoothBoundary,
     const bool undistortMasks,
-    const bool usePointsVisibilities
+    const bool usePointsVisibilities,
+    aliceVision::mesh::EFileType outputMeshFileType
     )
 {
     MaskCache maskCache(mp, masksFolders, undistortMasks);
@@ -517,7 +518,7 @@ void meshMasking(
     }
 
     // Save output mesh
-    filteredMesh.saveToObj(outputMeshPath);
+    filteredMesh.save(outputMeshPath, outputMeshFileType);
 
     ALICEVISION_LOG_INFO("Mesh file: \"" << outputMeshPath << "\" saved.");
 }
@@ -533,6 +534,7 @@ int main(int argc, char **argv)
     std::string inputMeshPath;
     std::vector<std::string> masksFolders;
     std::string outputMeshPath;
+    aliceVision::mesh::EFileType outputMeshFileType;
     std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel());
 
     int threshold = 1;
@@ -548,18 +550,20 @@ int main(int argc, char **argv)
         ("input,i", po::value<std::string>(&sfmFilePath)->default_value(sfmFilePath)->required(),
             "A SfMData file (*.sfm).")
         ("inputMesh,i", po::value<std::string>(&inputMeshPath)->required(),
-            "Input Mesh (OBJ file format).")
+            "Input Mesh")
         ("masksFolders", po::value<std::vector<std::string>>(&masksFolders)->multitoken(),
             "Use masks from specific folder(s).\n"
             "Filename should be the same or the image uid.")
         ("outputMesh,o", po::value<std::string>(&outputMeshPath)->required(),
-            "Output mesh (OBJ file format).")
+            "Output mesh")
         ("threshold", po::value<int>(&threshold)->default_value(threshold)->notifier(optInRange(1, INT_MAX, "threshold"))->required(),
             "The minimum number of visibility to keep a vertex.")
         ;
 
     po::options_description optionalParams("Optional parameters");
     optionalParams.add_options()
+        ("outputMeshFileType", po::value<aliceVision::mesh::EFileType>(&outputMeshFileType)->default_value(aliceVision::mesh::EFileType::GLTF),
+            "output mesh file type")
         ("invert", po::value<bool>(&invert)->default_value(invert),
             "Invert the mask.")
         ("smoothBoundary", po::value<bool>(&smoothBoundary)->default_value(smoothBoundary),
@@ -618,7 +622,7 @@ int main(int argc, char **argv)
     // check input mesh
     ALICEVISION_LOG_INFO("Load input mesh.");
     mesh::Mesh inputMesh;
-    inputMesh.loadFromObjAscii(inputMeshPath);
+    inputMesh.load(inputMeshPath);
 
     // check sfm file
     if(!sfmFilePath.empty() && !fs::exists(sfmFilePath) && !fs::is_regular_file(sfmFilePath))
@@ -672,7 +676,7 @@ int main(int argc, char **argv)
     }
 
     ALICEVISION_LOG_INFO("Mask mesh");
-    meshMasking(mp, inputMesh, masksFolders, outputMeshPath, threshold, invert, smoothBoundary, undistortMasks, usePointsVisibilities);
+    meshMasking(mp, inputMesh, masksFolders, outputMeshPath, threshold, invert, smoothBoundary, undistortMasks, usePointsVisibilities, outputMeshFileType);
     ALICEVISION_LOG_INFO("Task done in (s): " + std::to_string(timer.elapsed()));
     return EXIT_SUCCESS;
 }
diff --git a/src/software/pipeline/main_meshing.cpp b/src/software/pipeline/main_meshing.cpp
index 23443fa360..029f8a1243 100644
--- a/src/software/pipeline/main_meshing.cpp
+++ b/src/software/pipeline/main_meshing.cpp
@@ -254,6 +254,7 @@ int aliceVision_main(int argc, char* argv[])
     std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel());
     std::string sfmDataFilename;
     std::string outputMesh;
+    aliceVision::mesh::EFileType outputMeshFileType;
     std::string outputDensePointCloud;
     std::string depthMapsFolder;
     EPartitioningMode partitioningMode = ePartitioningSingleBlock;
@@ -294,10 +295,12 @@ int aliceVision_main(int argc, char* argv[])
         ("output,o", po::value<std::string>(&outputDensePointCloud)->required(),
           "Output Dense SfMData file.")
         ("outputMesh,o", po::value<std::string>(&outputMesh)->required(),
-          "Output mesh (OBJ file format).");
+          "Output mesh");
 
     po::options_description optionalParams("Optional parameters");
     optionalParams.add_options()
+        ("outputMeshFileType", po::value<aliceVision::mesh::EFileType>(&outputMeshFileType)->default_value(aliceVision::mesh::EFileType::OBJ),
+            "output mesh file type")
         ("depthMapsFolder", po::value<std::string>(&depthMapsFolder),
             "Input filtered depth maps folder.")
         ("boundingBox", po::value<BoundingBox>(&boundingBox),
@@ -602,7 +605,8 @@ int aliceVision_main(int argc, char* argv[])
     sfmDataIO::Save(densePointCloud, outputDensePointCloud, sfmDataIO::ESfMData::ALL_DENSE);
 
     ALICEVISION_LOG_INFO("Save obj mesh file.");
-    mesh->saveToObj(outputMesh);
+    ALICEVISION_LOG_INFO("OUTPUT MESH " << outputMesh);
+    mesh->save(outputMesh, outputMeshFileType);
     delete mesh;
 
 
diff --git a/src/software/pipeline/main_texturing.cpp b/src/software/pipeline/main_texturing.cpp
index c231627aae..3f0ea2583c 100644
--- a/src/software/pipeline/main_texturing.cpp
+++ b/src/software/pipeline/main_texturing.cpp
@@ -20,6 +20,8 @@
 #include <aliceVision/system/main.hpp>
 #include <aliceVision/system/Timer.hpp>
 
+#include <geogram/basic/common.h>
+
 #include <boost/program_options.hpp>
 #include <boost/filesystem.hpp>
 
@@ -44,10 +46,13 @@ int aliceVision_main(int argc, char* argv[])
 
     std::string verboseLevel = system::EVerboseLevel_enumToString(system::Logger::getDefaultVerboseLevel());
     std::string sfmDataFilename;
-    std::string inputMeshFilepath;
+
+    std::string inputMeshFilepath;    // Model to texture (HighPoly for diffuse, LowPoly for Diffuse+Normal)
+    std::string inputRefMeshFilepath; // HighPoly for NormalMap
+    aliceVision::mesh::EFileType outputMeshFileType;
+
     std::string outputFolder;
     std::string imagesFolder;
-    std::string outTextureFileTypeName = imageIO::EImageFileType_enumToString(imageIO::EImageFileType::PNG);
     std::string processColorspaceName = imageIO::EImageColorSpace_enumToString(imageIO::EImageColorSpace::SRGB);
     bool flipNormals = false;
     bool correctEV = false;
@@ -56,6 +61,10 @@ int aliceVision_main(int argc, char* argv[])
     std::string unwrapMethod = mesh::EUnwrapMethod_enumToString(mesh::EUnwrapMethod::Basic);
     std::string visibilityRemappingMethod = mesh::EVisibilityRemappingMethod_enumToString(texParams.visibilityRemappingMethod);
 
+    mesh::BumpMappingParams bumpMappingParams;
+    imageIO::EImageFileType normalFileType;
+    imageIO::EImageFileType heightFileType;
+
     po::options_description allParams("AliceVision texturing");
 
     po::options_description requiredParams("Required parameters");
@@ -65,10 +74,12 @@ int aliceVision_main(int argc, char* argv[])
         ("inputMesh", po::value<std::string>(&inputMeshFilepath)->required(),
             "Input mesh to texture.")
         ("output,o", po::value<std::string>(&outputFolder)->required(),
-            "Folder for output mesh: OBJ, material and texture files.");
+            "Folder for output mesh");
 
     po::options_description optionalParams("Optional parameters");
     optionalParams.add_options()
+        ("inputRefMesh", po::value<std::string>(&inputRefMeshFilepath),
+            "Optional input mesh to compute height maps and normal maps. If not provided, no additional maps with geometric information will be generated.")
         ("imagesFolder", po::value<std::string>(&imagesFolder),
           "Use images from a specific folder instead of those specify in the SfMData file.\n"
           "Filename should be the image uid.")
@@ -76,8 +87,18 @@ int aliceVision_main(int argc, char* argv[])
             "Output texture size")
         ("downscale", po::value<unsigned int>(&texParams.downscale)->default_value(texParams.downscale),
             "Texture downscale factor")
-        ("outputTextureFileType", po::value<std::string>(&outTextureFileTypeName)->default_value(outTextureFileTypeName),
-          imageIO::EImageFileType_informations().c_str())
+        ("outputMeshFileType", po::value<aliceVision::mesh::EFileType>(&outputMeshFileType)->default_value(aliceVision::mesh::EFileType::OBJ),
+            "output mesh file type")
+        ("colorMappingFileType", po::value<imageIO::EImageFileType>(&texParams.textureFileType)->default_value(texParams.textureFileType),
+          imageIO::EImageFileType_informations().c_str())(
+        "heightFileType", po::value<imageIO::EImageFileType>(&heightFileType)->default_value(imageIO::EImageFileType::NONE),
+            imageIO::EImageFileType_informations().c_str())
+        ("normalFileType", po::value<imageIO::EImageFileType>(&normalFileType)->default_value(imageIO::EImageFileType::NONE),
+            imageIO::EImageFileType_informations().c_str())
+        ("displacementMappingFileType", po::value<imageIO::EImageFileType>(&bumpMappingParams.displacementFileType)->default_value(bumpMappingParams.displacementFileType),
+            imageIO::EImageFileType_informations().c_str())
+        ("bumpType", po::value<mesh::EBumpMappingType>(&bumpMappingParams.bumpType)->default_value(bumpMappingParams.bumpType),
+            "Use HeightMap for displacement or bump mapping")
         ("unwrapMethod", po::value<std::string>(&unwrapMethod)->default_value(unwrapMethod),
             "Method to unwrap input mesh if it does not have UV coordinates.\n"
             " * Basic (> 600k faces) fast and simple. Can generate multiple atlases.\n"
@@ -152,24 +173,30 @@ int aliceVision_main(int argc, char* argv[])
     ALICEVISION_COUT("Program called with the following parameters:");
     ALICEVISION_COUT(vm);
 
+    // set bump mapping file type
+    bumpMappingParams.bumpMappingFileType = (bumpMappingParams.bumpType == mesh::EBumpMappingType::Normal) ? normalFileType : heightFileType;
+
     // set verbose level
     system::Logger::get()->setLogLevel(verboseLevel);
 
+    GEO::initialize();
+
     texParams.visibilityRemappingMethod = mesh::EVisibilityRemappingMethod_stringToEnum(visibilityRemappingMethod);
     texParams.processColorspace = imageIO::EImageColorSpace_stringToEnum(processColorspaceName);
-    // set output texture file type
-    const imageIO::EImageFileType outputTextureFileType = imageIO::EImageFileType_stringToEnum(outTextureFileTypeName);
-
+    
     texParams.correctEV = mvsUtils::ImagesCache::ECorrectEV::NO_CORRECTION;
     if(correctEV) { texParams.correctEV = mvsUtils::ImagesCache::ECorrectEV::APPLY_CORRECTION; }
 
     // read the input SfM scene
-    ALICEVISION_LOG_INFO("Load dense point cloud.");
     sfmData::SfMData sfmData;
-    if(!sfmDataIO::Load(sfmData, sfmDataFilename, sfmDataIO::ESfMData::ALL_DENSE))
+    if(!sfmDataFilename.empty())
     {
-      ALICEVISION_LOG_ERROR("The input SfMData file '" << sfmDataFilename << "' cannot be read.");
-      return EXIT_FAILURE;
+        ALICEVISION_LOG_INFO("Load dense point cloud.");
+        if(!sfmDataIO::Load(sfmData, sfmDataFilename, sfmDataIO::ESfMData::ALL_DENSE))
+        {
+          ALICEVISION_LOG_ERROR("The input SfMData file '" << sfmDataFilename << "' cannot be read.");
+          return EXIT_FAILURE;
+        }
     }
 
     // initialization
@@ -183,7 +210,7 @@ int aliceVision_main(int argc, char* argv[])
     // load input mesh (to texture) obj file
     ALICEVISION_LOG_INFO("Load input mesh.");
     mesh.clear();
-    mesh.loadOBJWithAtlas(inputMeshFilepath, flipNormals);
+    mesh.loadWithAtlas(inputMeshFilepath, flipNormals);
 
     // load reference dense point cloud with visibilities
     ALICEVISION_LOG_INFO("Convert dense point cloud into ref mesh");
@@ -195,13 +222,16 @@ int aliceVision_main(int argc, char* argv[])
     {
         // Need visibilities to compute unwrap
         mesh.remapVisibilities(texParams.visibilityRemappingMethod, mp, refMesh);
-        ALICEVISION_LOG_INFO("Input mesh has no UV coordinates, start unwrapping (" + unwrapMethod +")");
+        ALICEVISION_LOG_INFO("Input mesh has no UV coordinates, start unwrapping (" + unwrapMethod + ")");
         mesh.unwrap(mp, mesh::EUnwrapMethod_stringToEnum(unwrapMethod));
         ALICEVISION_LOG_INFO("Unwrapping done.");
     }
 
     // save final obj file
-    mesh.saveAsOBJ(outputFolder, "texturedMesh", outputTextureFileType);
+    if(!inputMeshFilepath.empty())
+    {
+        mesh.saveAs(outputFolder, "texturedMesh", outputMeshFileType, texParams.textureFileType, bumpMappingParams);
+    }
 
     if(texParams.subdivisionTargetRatio > 0)
     {
@@ -223,9 +253,25 @@ int aliceVision_main(int argc, char* argv[])
         // mesh.saveAsOBJ(outputFolder, "subdividedMesh", outputTextureFileType);
     }
 
-    // generate textures
-    ALICEVISION_LOG_INFO("Generate textures.");
-    mesh.generateTextures(mp, outputFolder, outputTextureFileType);
+    // generate diffuse textures
+    if(!inputMeshFilepath.empty() && !sfmDataFilename.empty() && texParams.textureFileType != imageIO::EImageFileType::NONE)
+    {
+        ALICEVISION_LOG_INFO("Generate textures.");
+        mesh.generateTextures(mp, outputFolder, texParams.textureFileType);
+    }
+
+
+    if(!inputRefMeshFilepath.empty() && !inputMeshFilepath.empty() &&
+       (bumpMappingParams.bumpMappingFileType != imageIO::EImageFileType::NONE ||
+        bumpMappingParams.displacementFileType != imageIO::EImageFileType::NONE))
+    {
+        ALICEVISION_LOG_INFO("Generate height and normal maps.");
+
+        mesh::Mesh denseMesh;
+        denseMesh.load(inputRefMeshFilepath);
+
+        mesh.generateNormalAndHeightMaps(mp, denseMesh, outputFolder, bumpMappingParams);
+    }
 
     ALICEVISION_LOG_INFO("Task done in (s): " + std::to_string(timer.elapsed()));
     return EXIT_SUCCESS;