Add OpenSfM support

CMakeLists.txt

@@ -124,7 +124,7 @@ endif()
 add_library(gsplat_cpu vendor/gsplat-cpu/gsplat_cpu.cpp)
 target_include_directories(gsplat_cpu PRIVATE ${TORCH_INCLUDE_DIRS})
 
-add_executable(opensplat opensplat.cpp point_io.cpp nerfstudio.cpp model.cpp kdtree_tensor.cpp spherical_harmonics.cpp cv_utils.cpp utils.cpp project_gaussians.cpp rasterize_gaussians.cpp ssim.cpp optim_scheduler.cpp colmap.cpp input_data.cpp tensor_math.cpp)
+add_executable(opensplat opensplat.cpp point_io.cpp nerfstudio.cpp model.cpp kdtree_tensor.cpp spherical_harmonics.cpp cv_utils.cpp utils.cpp project_gaussians.cpp rasterize_gaussians.cpp ssim.cpp optim_scheduler.cpp colmap.cpp opensfm.cpp input_data.cpp tensor_math.cpp)
 set_property(TARGET opensplat PROPERTY CXX_STANDARD 17)
 target_include_directories(opensplat PRIVATE ${PROJECT_SOURCE_DIR}/vendor/glm ${GPU_INCLUDE_DIRS})
 target_link_libraries(opensplat PUBLIC ${STDPPFS_LIBRARY} ${GPU_LIBRARIES} ${GSPLAT_LIBS} ${TORCH_LIBRARIES} ${OpenCV_LIBS})

README.md

@@ -5,7 +5,7 @@ A free and open source implementation of 3D [gaussian splatting](https://www.you
 <img src="https://github.com/pierotofy/OpenSplat/assets/1951843/c9327c7c-31ad-402d-a5a5-04f7602ca5f5" width="49%" />
 <img src="https://github.com/pierotofy/OpenSplat/assets/1951843/eba4ae75-2c88-4c9e-a66b-608b574d085f" width="49%" />
 
-OpenSplat takes camera poses + sparse points in [COLMAP](https://colmap.github.io/) or [nerfstudio](https://docs.nerf.studio/quickstart/custom_dataset.html) project format and computes a [scene file](https://drive.google.com/file/d/1w-CBxyWNXF3omA8B_IeOsRmSJel3iwyr/view?usp=sharing) (.ply) that can be later imported for [viewing](https://github.com/shg8/3DGS.cpp), editing and rendering in other [software](https://github.com/MrNeRF/awesome-3D-gaussian-splatting?tab=readme-ov-file#open-source-implementations).
+OpenSplat takes camera poses + sparse points in [COLMAP](https://colmap.github.io/), [OpenSfM](https://github.com/mapillary/OpenSfM), [ODM](https://github.com/OpenDroneMap/ODM) or [nerfstudio](https://docs.nerf.studio/quickstart/custom_dataset.html) project format and computes a [scene file](https://drive.google.com/file/d/1w-CBxyWNXF3omA8B_IeOsRmSJel3iwyr/view?usp=sharing) (.ply) that can be later imported for [viewing](https://github.com/shg8/3DGS.cpp), editing and rendering in other [software](https://github.com/MrNeRF/awesome-3D-gaussian-splatting?tab=readme-ov-file#open-source-implementations).
 
 Graphics card recommended, but not required! OpenSplat runs the fastest on NVIDIA and AMD GPUs, but can also run entirely on the CPU (~100x slower).
 

constants.hpp

@@ -2,5 +2,6 @@
 #define CONSTANTS_H
 
 #define PI 3.14159265358979323846
+#define FLOAT_EPS 1e-9f
 
 #endif

input_data.cpp

@@ -7,6 +7,7 @@ using namespace torch::indexing;
 
 namespace ns{ InputData inputDataFromNerfStudio(const std::string &projectRoot); }
 namespace cm{ InputData inputDataFromColmap(const std::string &projectRoot); }
+namespace osfm { InputData inputDataFromOpenSfM(const std::string &projectRoot); }
 
 InputData inputDataFromX(const std::string &projectRoot){
     fs::path root(projectRoot);
@@ -15,6 +16,10 @@ InputData inputDataFromX(const std::string &projectRoot){
         return ns::inputDataFromNerfStudio(projectRoot);
     }else if (fs::exists(root / "sparse") || fs::exists(root / "cameras.bin")){
         return cm::inputDataFromColmap(projectRoot);
+    }else if (fs::exists(root / "reconstruction.json")){
+        return osfm::inputDataFromOpenSfM(projectRoot);
+    }else if (fs::exists(root / "opensfm" / "reconstruction.json")){
+        return osfm::inputDataFromOpenSfM((root / "opensfm").string());
     }else{
         throw std::runtime_error("Invalid project folder (must be either a colmap or nerfstudio project folder)");
     }

model.cpp

@@ -498,10 +498,10 @@ void Model::savePlySplat(const std::string &filename){
 
     float zeros[] = { 0.0f, 0.0f, 0.0f };
 
-    torch::Tensor meansCpu = (means.cpu() / scale) + translation;
+    torch::Tensor meansCpu = keepCrs ? (means.cpu() / scale) + translation : means.cpu();
     torch::Tensor featuresDcCpu = featuresDc.cpu();
     torch::Tensor opacitiesCpu = opacities.cpu();
-    torch::Tensor scalesCpu = torch::log((torch::exp(scales.cpu()) / scale));
+    torch::Tensor scalesCpu = keepCrs ? torch::log((torch::exp(scales.cpu()) / scale)) : scales.cpu();
     torch::Tensor quatsCpu = quats.cpu();
 
     for (size_t i = 0; i < numPoints; i++) {
@@ -535,7 +535,7 @@ void Model::saveDebugPly(const std::string &filename){
     o << "property uchar blue" << std::endl;
     o << "end_header" << std::endl;
 
-    torch::Tensor meansCpu = (means.cpu() / scale) + translation;
+    torch::Tensor meansCpu = keepCrs ? (means.cpu() / scale) + translation : means.cpu();
     torch::Tensor rgbsCpu = (sh2rgb(featuresDc.cpu()) * 255.0f).toType(torch::kUInt8);
 
     for (size_t i = 0; i < numPoints; i++) {

model.hpp

@@ -23,12 +23,12 @@ struct Model{
   Model(const InputData &inputData, int numCameras,
         int numDownscales, int resolutionSchedule, int shDegree, int shDegreeInterval, 
         int refineEvery, int warmupLength, int resetAlphaEvery, float densifyGradThresh, float densifySizeThresh, int stopScreenSizeAt, float splitScreenSize,
-        int maxSteps,
+        int maxSteps, bool keepCrs,
         const torch::Device &device) :
     numCameras(numCameras),
     numDownscales(numDownscales), resolutionSchedule(resolutionSchedule), shDegree(shDegree), shDegreeInterval(shDegreeInterval), 
     refineEvery(refineEvery), warmupLength(warmupLength), resetAlphaEvery(resetAlphaEvery), stopSplitAt(maxSteps / 2), densifyGradThresh(densifyGradThresh), densifySizeThresh(densifySizeThresh), stopScreenSizeAt(stopScreenSizeAt), splitScreenSize(splitScreenSize),
-    maxSteps(maxSteps),
+    maxSteps(maxSteps), keepCrs(keepCrs),
     device(device), ssim(11, 3){
 
     long long numPoints = inputData.points.xyz.size(0);
@@ -131,6 +131,7 @@ struct Model{
   int stopScreenSizeAt;
   float splitScreenSize;
   int maxSteps;
+  bool keepCrs;
 
   float scale;
   torch::Tensor translation;

nerfstudio.cpp

@@ -109,6 +109,7 @@ void from_json(const json& j, Transforms &t){
 Transforms readTransforms(const std::string &filename){
     std::ifstream f(filename);
     json data = json::parse(f);
+    f.close();
     return data.template get<Transforms>();
 }
 

opensfm.cpp

@@ -0,0 +1,153 @@
+#include <filesystem>
+#include <cstdlib>
+#include "vendor/json/json.hpp"
+#include "opensfm.hpp"
+#include "point_io.hpp"
+#include "cv_utils.hpp"
+#include "tensor_math.hpp"
+
+namespace fs = std::filesystem;
+
+using json = nlohmann::json;
+using namespace torch::indexing;
+
+namespace osfm{
+
+void from_json(const json& j, Cam &c){
+    j.at("projection_type").get_to(c.projectionType);
+    if (j.contains("width")) j.at("width").get_to(c.width);
+    if (j.contains("height")) j.at("height").get_to(c.height);
+    if (j.contains("focal_x")) j.at("focal_x").get_to(c.fx);
+    if (j.contains("focal_y")) j.at("focal_y").get_to(c.fy);
+    if (j.contains("focal")){
+        j.at("focal").get_to(c.fx);
+        j.at("focal").get_to(c.fy);
+    }
+    if (j.contains("c_x")) j.at("c_x").get_to(c.cx);
+    if (j.contains("c_y")) j.at("c_y").get_to(c.cy);
+    if (j.contains("k1")) j.at("k1").get_to(c.k1);
+    if (j.contains("k2")) j.at("k2").get_to(c.k2);
+    if (j.contains("p1")) j.at("p1").get_to(c.p1);
+    if (j.contains("p2")) j.at("p2").get_to(c.p2);
+    if (j.contains("k3")) j.at("k3").get_to(c.k3);
+}
+
+void from_json(const json& j, Shot &s){
+    j.at("rotation").get_to(s.rotation);
+    j.at("translation").get_to(s.translation);
+    j.at("camera").get_to(s.camera);
+}    
+
+void from_json(const json& j, Point &p){
+    j.at("coordinates").get_to(p.coordinates);
+    j.at("color").get_to(p.color);
+}
+
+void from_json(const json& j, Reconstruction &r){
+    j.at("cameras").get_to(r.cameras);
+    j.at("shots").get_to(r.shots);
+    j.at("points").get_to(r.points);
+
+}
+
+InputData inputDataFromOpenSfM(const std::string &projectRoot){
+    InputData ret;
+    fs::path nsRoot(projectRoot);
+    fs::path reconstructionPath = nsRoot / "reconstruction.json";
+    fs::path imageListPath = nsRoot / "image_list.txt";
+
+    if (!fs::exists(reconstructionPath)) throw std::runtime_error(reconstructionPath.string() + " does not exist");
+    if (!fs::exists(imageListPath)) throw std::runtime_error(imageListPath.string() + " does not exist");
+
+    std::ifstream f(reconstructionPath.string());
+    json data = json::parse(f);
+    f.close();
+
+    std::unordered_map<std::string, std::string> images;
+    f.open(imageListPath.string());
+    std::string line;
+    while(std::getline(f, line)){
+        fs::path p(line);
+        if (p.is_absolute()) images[p.filename().string()] = line;
+        else images[p.filename().string()] = fs::absolute(nsRoot / p).string();
+    }
+    f.close();
+    auto reconstructions = data.template get<std::vector<Reconstruction>>();
+
+    if (reconstructions.size() == 0) throw std::runtime_error("No reconstructions found");
+    if (reconstructions.size() > 1) std::cout << "Warning: multiple OpenSfM reconstructions found, choosing the first" << std::endl;
+
+    auto reconstruction = reconstructions[0];
+    auto shots = reconstruction.shots;
+    auto cameras = reconstruction.cameras;
+    auto points = reconstruction.points;
+
+    torch::Tensor unorientedPoses = torch::zeros({static_cast<long int>(shots.size()), 4, 4}, torch::kFloat32);
+    size_t i = 0;
+    for (const auto &s : shots){
+        Shot shot = s.second;
+
+        torch::Tensor rotation = rodriguesToRotation(torch::from_blob(shot.rotation.data(), {static_cast<long>(shot.rotation.size())}, torch::kFloat32));
+        torch::Tensor translation = torch::from_blob(shot.translation.data(), {static_cast<long>(shot.translation.size())}, torch::kFloat32);
+        torch::Tensor w2c = torch::eye(4, torch::kFloat32);
+        w2c.index_put_({Slice(None, 3), Slice(None, 3)}, rotation);
+        w2c.index_put_({Slice(None, 3), Slice(3,4)}, translation.reshape({3, 1}));
+
+        unorientedPoses[i] = torch::linalg::inv(w2c);
+
+        // Convert OpenSfM's camera CRS (OpenCV) to OpenGL
+        unorientedPoses[i].index_put_({Slice(0, 3), Slice(1,3)}, unorientedPoses[i].index({Slice(0, 3), Slice(1,3)}) * -1.0f);
+        i++;
+    }
+
+    auto r = autoScaleAndCenterPoses(unorientedPoses);
+    torch::Tensor poses = std::get<0>(r);
+    ret.translation = std::get<1>(r);
+    ret.scale = std::get<2>(r);
+
+    i = 0;
+    for (const auto &s : shots){
+        std::string filename = s.first;
+        Shot shot = s.second;
+
+        Cam &c = cameras[shot.camera];
+        if (c.projectionType != "perspective" && c.projectionType != "brown"){
+            throw std::runtime_error("Camera projection type " + c.projectionType + " is not supported");
+        }
+
+        float normalizer = static_cast<float>((std::max)(c.width, c.height));
+        ret.cameras.emplace_back(Camera(c.width, c.height, 
+                            static_cast<float>(c.fx * normalizer), static_cast<float>(c.fy * normalizer), 
+                            static_cast<float>(static_cast<float>(c.width) / 2.0f + normalizer * c.cx), static_cast<float>(static_cast<float>(c.height) / 2.0f + normalizer * c.cy), 
+                            static_cast<float>(c.k1), static_cast<float>(c.k2), static_cast<float>(c.k3), 
+                            static_cast<float>(c.p1), static_cast<float>(c.p2),  
+
+                            poses[i++], images[filename]));
+    }
+
+    size_t numPoints = points.size();
+    torch::Tensor xyz = torch::zeros({static_cast<long>(numPoints), 3}, torch::kFloat32);
+    torch::Tensor rgb = torch::zeros({static_cast<long>(numPoints), 3}, torch::kUInt8);
+
+    i = 0;
+    for (const auto &pt: points){
+        Point p = pt.second;
+
+        xyz[i][0] = p.coordinates[0];
+        xyz[i][1] = p.coordinates[1];
+        xyz[i][2] = p.coordinates[2];
+
+        rgb[i][0] = static_cast<uint8_t>(p.color[0]);
+        rgb[i][1] = static_cast<uint8_t>(p.color[1]);
+        rgb[i][2] = static_cast<uint8_t>(p.color[2]);
+
+        i++;
+    }
+
+    ret.points.xyz = (xyz - ret.translation) * ret.scale;
+    ret.points.rgb = rgb;
+
+    return ret;
+}
+
+}

opensfm.hpp

@@ -0,0 +1,54 @@
+#ifndef OPENSFM_H
+#define OPENSFM_H
+
+#include <iostream>
+#include <string>
+#include <fstream>
+#include <unordered_map>
+#include <torch/torch.h>
+#include "input_data.hpp"
+#include "vendor/json/json_fwd.hpp"
+
+using json = nlohmann::json;
+
+namespace osfm{
+    struct Cam{
+        std::string projectionType = "";
+        int width = 0;
+        int height = 0;
+        double fx = 0;
+        double fy = 0;
+        double cx = 0;
+        double cy = 0;
+        double k1 = 0;
+        double k2 = 0;
+        double p1 = 0;
+        double p2 = 0;
+        double k3 = 0;
+    };
+    void from_json(const json& j, Cam &c);
+
+    struct Shot{
+        std::vector<float> rotation = {0.0f, 0.0f, 0.0f};
+        std::vector<float> translation = {0.0f, 0.0f, 0.0f};
+        std::string camera = "";
+    };
+    void from_json(const json& j, Shot &s);
+
+    struct Point{
+        std::vector<float> color;
+        std::vector<float> coordinates;
+    };
+    void from_json(const json& j, Point &p);
+
+    struct Reconstruction{
+        std::unordered_map<std::string, Cam> cameras;
+        std::unordered_map<std::string, Shot> shots;
+        std::unordered_map<std::string, Point> points;
+    };
+    void from_json(const json& j, Reconstruction &r);
+
+    InputData inputDataFromOpenSfM(const std::string &projectRoot);
+}   
+
+#endif

opensplat.cpp

@@ -18,6 +18,7 @@ int main(int argc, char *argv[]){
         ("val", "Withhold a camera shot for validating the scene loss")
         ("val-image", "Filename of the image to withhold for validating scene loss", cxxopts::value<std::string>()->default_value("random"))
         ("val-render", "Path of the directory where to render validation images", cxxopts::value<std::string>()->default_value(""))
+        ("keep-crs", "Retain the project input's coordinate reference system")
         ("cpu", "Force CPU execution")
 
         ("n,num-iters", "Number of iterations to run", cxxopts::value<int>()->default_value("30000"))
@@ -38,7 +39,7 @@ int main(int argc, char *argv[]){
         ("h,help", "Print usage")
         ;
     options.parse_positional({ "input" });
-    options.positional_help("[colmap or nerfstudio project path]");
+    options.positional_help("[colmap/nerfstudio/opensfm/odm project path]");
     cxxopts::ParseResult result;
     try {
         result = options.parse(argc, argv);
@@ -61,7 +62,7 @@ int main(int argc, char *argv[]){
     const std::string valImage = result["val-image"].as<std::string>();
     const std::string valRender = result["val-render"].as<std::string>();
     if (!valRender.empty() && !fs::exists(valRender)) fs::create_directories(valRender);
-
+    const bool keepCrs = result.count("keep-crs") > 0;
     const float downScaleFactor = (std::max)(result["downscale-factor"].as<float>(), 1.0f);
     const int numIters = result["num-iters"].as<int>();
     const int numDownscales = result["num-downscales"].as<int>();
@@ -103,7 +104,7 @@ int main(int argc, char *argv[]){
                     cams.size(),
                     numDownscales, resolutionSchedule, shDegree, shDegreeInterval, 
                     refineEvery, warmupLength, resetAlphaEvery, densifyGradThresh, densifySizeThresh, stopScreenSizeAt, splitScreenSize,
-                    numIters,
+                    numIters, keepCrs,
                     device);
 
         std::vector< size_t > camIndices( cams.size() );

tensor_math.cpp

@@ -25,7 +25,6 @@ torch::Tensor quatToRotMat(const torch::Tensor &quat){
             1.0 - 2.0 * (x.pow(2) + y.pow(2))
         }, -1)
     }, -2);
-
 }
 
 std::tuple<torch::Tensor, torch::Tensor, float> autoScaleAndCenterPoses(const torch::Tensor &poses){
@@ -66,4 +65,29 @@ torch::Tensor rotationMatrix(const torch::Tensor &a, const torch::Tensor &b){
     skew[2][1] = v[0];
 
     return torch::eye(3) + skew + torch::matmul(skew, skew * ((1 - c) / (s.pow(2) + EPS)));
+}
+
+torch::Tensor rodriguesToRotation(const torch::Tensor &rodrigues){
+    float theta = torch::linalg::vector_norm(rodrigues, 2, { -1 }, true, torch::kFloat32).item<float>();
+    if (theta < FLOAT_EPS){
+        return torch::eye(3, torch::kFloat32);
+    }
+    torch::Tensor r = rodrigues / theta;
+    torch::Tensor ident = torch::eye(3, torch::kFloat32);
+    float a = r[0].item<float>();
+    float b = r[1].item<float>();
+    float c = r[2].item<float>();
+    torch::Tensor rrT = torch::tensor({
+        {a * a, a * b, a * c},
+        {b * a, b * b, b * c},
+        {c * a, c * b, c * c}
+    }, torch::kFloat32);
+    torch::Tensor rCross = torch::tensor({
+        {0.0f, -c, b},
+        {c, 0.0f, -a},
+        {-b, a, 0.0f}
+    }, torch::kFloat32);
+    float cosTheta = std::cos(theta);
+
+    return cosTheta * ident + (1 - cosTheta) * rrT + std::sin(theta) * rCross;
 }

tensor_math.hpp

@@ -3,10 +3,11 @@
 
 #include <torch/torch.h>
 #include <tuple>
+#include "constants.hpp"
 
 torch::Tensor quatToRotMat(const torch::Tensor &quat);
 std::tuple<torch::Tensor, torch::Tensor, float> autoScaleAndCenterPoses(const torch::Tensor &poses);
 torch::Tensor rotationMatrix(const torch::Tensor &a, const torch::Tensor &b);
-
+torch::Tensor rodriguesToRotation(const torch::Tensor &rodrigues);
 
 #endif