[Photometric Stereo] MultiView fusion in Texturing

src/CMakeLists.txt

@@ -51,7 +51,7 @@ option(ALICEVISION_BUILD_SFM "Build AliceVision SfM part" ON)
 option(ALICEVISION_BUILD_MVS "Build AliceVision MVS part" ON)
 option(ALICEVISION_BUILD_HDR "Build AliceVision HDR part" ON)
 option(ALICEVISION_BUILD_SEGMENTATION "Build AliceVision Segmentation part" ON)
-option(ALICEVISION_BUILD_STEREOPHOTOMETRY "Build AliceVision StereoPhotometry part" ON)
+option(ALICEVISION_BUILD_PHOTOMETRICSTEREO "Build AliceVision Photometric Stereo part" ON)
 option(ALICEVISION_BUILD_PANORAMA "Build AliceVision panorama part" ON)
 option(ALICEVISION_BUILD_SOFTWARE "Build AliceVision command line tools." ON)
 option(ALICEVISION_BUILD_COVERAGE "Enable code coverage generation (gcc only)" OFF)
@@ -94,7 +94,7 @@ if(NOT ALICEVISION_BUILD_SFM)
     SET(ALICEVISION_BUILD_MVS OFF)
     SET(ALICEVISION_BUILD_HDR OFF)
     SET(ALICEVISION_BUILD_SEGMENTATION OFF)
-    SET(ALICEVISION_BUILD_STEREOPHOTOMETRY OFF)
+    SET(ALICEVISION_BUILD_PHOTOMETRICSTEREO OFF)
     SET(ALICEVISION_BUILD_PANORAMA OFF)
     SET(ALICEVISION_BUILD_LIDAR OFF)
 endif()

src/aliceVision/CMakeLists.txt

@@ -51,7 +51,6 @@ endif()
 # MVS modules
 
 if(ALICEVISION_BUILD_MVS)
-  add_subdirectory(lightingEstimation)
   add_subdirectory(mesh)
   add_subdirectory(mvsData)
   add_subdirectory(mvsUtils)
@@ -72,12 +71,13 @@ if(ALICEVISION_BUILD_SFM AND ALICEVISION_BUILD_MVS)
   add_subdirectory(sfmMvsUtils)
 endif()
 
-if (ALICEVISION_BUILD_STEREOPHOTOMETRY)
+if (ALICEVISION_BUILD_PHOTOMETRICSTEREO)
     if(ALICEVISION_HAVE_OPENCV)
         add_subdirectory(photometricStereo)
         if(ALICEVISION_HAVE_ONNX)
             add_subdirectory(sphereDetection)
         endif()
+        add_subdirectory(lightingEstimation)
     endif()
 endif()
 

src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp

@@ -67,7 +67,7 @@ BOOST_AUTO_TEST_CASE(fuseCut_delaunayGraphCut)
     const NViewDatasetConfigurator config(1000, 1000, 500, 500, 1, 0);
     SfMData sfmData = generateSfm(config, 6);
 
-    mvsUtils::MultiViewParams mp(sfmData, "", "", "", false);
+    mvsUtils::MultiViewParams mp(sfmData, "", "", "");
 
     mp.userParams.put("LargeScale.universePercentile", 0.999);
     mp.userParams.put("delaunaycut.seed", 1);

src/aliceVision/lightingEstimation/CMakeLists.txt

@@ -3,20 +3,24 @@ set(lightingEstimation_files_headers
     augmentedNormals.hpp
     lightingEstimation.hpp
     lightingCalibration.hpp
+    ellipseGeometry.hpp
 )
 
 # Sources
 set(lightingEstimation_files_sources
     augmentedNormals.cpp
     lightingEstimation.cpp
     lightingCalibration.cpp
+    ellipseGeometry.cpp
 )
 
 alicevision_add_library(aliceVision_lightingEstimation
   SOURCES ${lightingEstimation_files_headers} ${lightingEstimation_files_sources}
   PUBLIC_LINKS
+    ${OpenCV_LIBS}
     aliceVision_image
     aliceVision_system
+    aliceVision_photometricStereo
 )
 
 

src/aliceVision/lightingEstimation/augmentedNormals.hpp

@@ -67,8 +67,8 @@ class AugmentedNormal : public Eigen::Matrix<float, 9, 1, 0, 9, 1>
     inline const T& nz() const { return (*this)(2); }
     inline T& nz() { return (*this)(2); }
 
-    inline const T& nambiant() const { return (*this)(3); }
-    inline T& nambiant() { return (*this)(3); }
+    inline const T& nambient() const { return (*this)(3); }
+    inline T& nambient() { return (*this)(3); }
 
     inline const T& nx_ny() const { return (*this)(4); }
     inline T& nx_ny() { return (*this)(4); }

src/aliceVision/lightingEstimation/ellipseGeometry.cpp

@@ -0,0 +1,232 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2023 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#include "ellipseGeometry.hpp"
+#include <aliceVision/photometricStereo/photometricDataIO.hpp>
+#include <aliceVision/image/io.hpp>
+
+#include <opencv2/imgcodecs.hpp>
+#include <opencv2/imgproc.hpp>
+
+#include <algorithm>
+#include <limits>
+
+namespace aliceVision {
+namespace lightingEstimation {
+
+#define PI 3.14159265
+
+void quadraticFromEllipseParameters(const std::array<float, 5>& ellipseParameters,
+                                    Eigen::Matrix3f& Q)
+{
+    float phi = ellipseParameters[0]*PI/180.0;
+    float c_x = ellipseParameters[1];
+    float c_y = ellipseParameters[2];
+    float b = ellipseParameters[3];
+    float a = ellipseParameters[4];
+
+    float A = a*a*sin(phi)*sin(phi) + b*b*cos(phi)*cos(phi);
+    float B = 2*(b*b-a*a)*sin(phi)*cos(phi);
+    float C = a*a*cos(phi)*cos(phi) + b*b*sin(phi)*sin(phi);
+    float D = -2*A*c_x - B*c_y;
+    float E = -B*c_x - 2*C*c_y;
+    float F = A*c_x*c_x + C*c_y*c_y - a*a*b*b + B*c_x*c_y;
+
+    Q << A, B/2, D/2,
+         B/2, C, E/2,
+         D/2, E/2, F;
+
+    Q /= Q.norm();
+}
+
+int findUniqueIndex(const std::vector<int>& vec) {
+
+    std::unordered_map<int, int> countMap;
+
+    // Compter le nombre d'occurrences de chaque élément dans le vecteur
+    for (int num : vec) {
+        countMap[num]++;
+    }
+
+    int toReturn = -1;
+    // Trouver l'élément avec une seule occurrence
+    for (size_t i = 0; i < vec.size(); ++i) {
+        if (countMap[vec[i]] == 1) {
+            toReturn = i;
+        }
+    }
+
+    return toReturn;
+}
+
+void estimateSphereCenter(const std::array<float, 5>& ellipseParameters,
+                          const float sphereRadius,
+                          const Eigen::Matrix3f& K,
+                          std::array<float, 3>& sphereCenter)
+{
+    Eigen::Matrix3f Q;
+    quadraticFromEllipseParameters(ellipseParameters, Q);
+
+    Eigen::Matrix3f M = K.transpose()*Q*K;
+    Eigen::EigenSolver<Eigen::MatrixXf> es(M);
+
+    Eigen::Vector3f eigval = es.eigenvalues().real();
+
+    std::vector<int> eigvalSign;
+
+    for (int i = 0; i < 3; ++i) {
+        eigvalSign.push_back((eigval[i] > 0) ? 1 : -1);
+    }
+
+    int index = findUniqueIndex(eigvalSign);
+    float uniqueEigval = eigval[index];
+
+    float dist = sqrt(1 - ((eigval[0] + eigval[1] + eigval[2] - uniqueEigval)/2) / uniqueEigval);
+
+    Eigen::Vector3f eigvect = es.eigenvectors().col(index).real();
+    float sign = eigvect[2] > 0 ? 1 : -1;
+
+    float norm = eigvect.norm();
+    float C_factor = sphereRadius*dist*sign/norm;
+    Eigen::Vector3f C = C_factor*eigvect;
+    sphereCenter[0]=C[0];
+    sphereCenter[1]=C[1];
+    sphereCenter[2]=C[2];
+}
+
+void sphereRayIntersection(const Eigen::Vector3f& direction,
+                           const std::array<float, 3>& sphereCenter,
+                           const float sphereRadius,
+                           float& delta,
+                           Eigen::Vector3f& normal)
+{
+    float a = direction.dot(direction);
+
+    Eigen::Vector3f spCenter;
+    spCenter << sphereCenter[0], sphereCenter[1], sphereCenter[2];
+
+    float b = -2*direction.dot(spCenter);
+    float c = spCenter.dot(spCenter) - sphereRadius*sphereRadius;
+
+    delta = b*b - 4*a*c;
+
+    float factor;
+
+    if (delta >= 0) {
+        factor = (-b - sqrt(delta))/(2*a);
+        normal = (direction*factor - spCenter)/sphereRadius;
+    }
+    else {
+        delta = -1;
+    }
+}
+
+void estimateSphereNormals(const std::array<float, 3>& sphereCenter,
+                           const float sphereRadius,
+                           const Eigen::Matrix3f& K,
+                           image::Image<image::RGBfColor>& normals,
+                           image::Image<float>& newMask) 
+{
+    Eigen::Matrix3f invK = K.inverse();
+
+    for (int i = 0; i < normals.rows(); ++i) {
+        for (int j = 0; j < normals.cols(); ++j) {
+            // Get homogeneous coordinates of the pixel :
+            Eigen::Vector3f coordinates = Eigen::Vector3f(j, i, 1.0f);
+
+            // Get the direction of the ray :
+            Eigen::Vector3f direction = invK*coordinates;
+
+            // Estimate the interception of the ray with the sphere :
+            float delta = 0.0;
+            Eigen::Vector3f normal = Eigen::Vector3f(0.0f, 0.0f, 0.0f);
+
+            sphereRayIntersection(direction, sphereCenter, sphereRadius, delta, normal);
+
+            // If the ray intercepts the sphere we add this pixel to the new mask :
+            if(delta > 0)
+            {
+                newMask(i, j) = 1.0;
+                normals(i, j) = image::RGBfColor(normal[0], normal[1], normal[2]);
+            }
+            else
+            {
+                newMask(i, j) = 0.0;
+                normals(i, j) = image::RGBfColor(0.0, 0.0, 0.0);
+            }
+        }
+    }
+}
+
+void getRealNormalOnSphere(const cv::Mat& maskCV,
+                       const Eigen::Matrix3f& K,
+                       const float sphereRadius,
+                       image::Image<image::RGBfColor>& normals,
+                       image::Image<float>& newMask)
+{
+
+    // Apply a threshold to the image
+    int thresholdValue = 150;
+    cv::Mat thresh;
+    cv::threshold(maskCV, thresh, thresholdValue, 255, 0);
+
+    // Find contours
+    std::vector<std::vector<cv::Point>> contours;
+    std::vector<cv::Vec4i> hierarchy;
+    cv::findContours(thresh, contours, hierarchy, cv::RETR_TREE, cv::CHAIN_APPROX_SIMPLE);
+
+    // Fit the ellipse to the first contour
+    std::vector<cv::Point> cnt = contours[0];
+    cv::RotatedRect ellipse = cv::fitEllipse(cnt);
+    float angle = ellipse.angle;
+    cv::Size2f size = ellipse.size;
+    cv::Point2f center = ellipse.center;
+
+    // Ellipse is converted as five-parameter array
+    std::array<float, 5> ellipseParameters;
+
+    ellipseParameters[0] = angle;
+    ellipseParameters[1] = center.x;
+    ellipseParameters[2] = center.y;
+    ellipseParameters[3] = size.width/2;
+    ellipseParameters[4] = size.height/2;
+
+    std::array<float, 3> sphereCenter;
+    estimateSphereCenter(ellipseParameters, sphereRadius, K, sphereCenter);
+    estimateSphereNormals(sphereCenter, sphereRadius, K, normals, newMask);
+}
+
+void getEllipseMaskFromSphereParameters(const std::array<float, 3>& sphereParam, const Eigen::Matrix3f& K, std::array<float, 5>& ellipseParameters, cv::Mat maskCV)
+{
+
+    // Distance between center of image and center of disk :
+    float imX = K(0, 2);
+    float imY = K(1, 2);
+    float f = K(0, 0);
+    float delta = sqrt((sphereParam[0])*(sphereParam[0]) + (sphereParam[1])*(sphereParam[1]));
+    // sphere params = x,y, radius
+    // ellipse params  = angle, x, y, semi minor axe, semi major axe
+
+    // Ellipse from sphere parameters
+    // semi minor axe = radius;
+    // semi major axe = formula from paper
+    // main direction = disc center to picture center
+    // ellipse center = disc center
+
+    float radians =  atan((sphereParam[0]) / (sphereParam[1]));
+    ellipseParameters[0] = radians * (180.0/3.141592653589793238463);
+
+    ellipseParameters[1] = sphereParam[0];
+    ellipseParameters[2] = sphereParam[1];
+    ellipseParameters[3] = sphereParam[2];
+
+
+    // a² = b² ((distance between image center and disc center)² + f² + b²)/(f² + b²)
+    ellipseParameters[4] = sqrt((ellipseParameters[3]*ellipseParameters[3]) * (delta * delta + f*f + ellipseParameters[3]*ellipseParameters[3])/(f*f + ellipseParameters[3]*ellipseParameters[3]));
+}
+
+} // namespace lightingEstimation
+} // namespace aliceVision

src/aliceVision/lightingEstimation/ellipseGeometry.hpp

@@ -0,0 +1,75 @@
+// This file is part of the AliceVision project.
+// Copyright (c) 2024 AliceVision contributors.
+// This Source Code Form is subject to the terms of the Mozilla Public License,
+// v. 2.0. If a copy of the MPL was not distributed with this file,
+// You can obtain one at https://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include <aliceVision/image/Image.hpp>
+
+#include <opencv2/imgcodecs.hpp>
+#include <opencv2/imgproc.hpp>
+
+#include <string>
+#include <vector>
+#include <array>
+
+namespace aliceVision {
+namespace lightingEstimation {
+
+/**
+ * @brief Estimate the center of a sphere from the ellipse parameters
+ * @param[in] ellipseParameters An array of 5 floating-point: the parameters of the ellipse
+ * @param[in] sphereRadius The radius of the sphere
+ * @param[in] K Intrinsic parameters of the camera
+ * @param[out] sphereCenter An array of 3 floating-point: the coordinates of the sphere center in the picture frame
+ */
+void estimateSphereCenter(const std::array<float, 5>& ellipseParameters,
+                          const float sphereRadius,
+                          const Eigen::Matrix3f& K,
+                          std::array<float, 3>& sphereCenter);
+
+void sphereRayIntersection(const Eigen::Vector3f& direction,
+                           const std::array<float, 3>& sphereCenter,
+                           const float sphereRadius,
+                           float& delta,
+                           Eigen::Vector3f& normal);
+
+void quadraticFromEllipseParameters(const std::array<float, 5>& ellipseParameters,
+                                       Eigen::Matrix3f& Q);
+
+int findUniqueIndex(const std::vector<int>& vec);
+
+/**
+ * @brief Estimate the normals of a sphere from a mask
+ * @param[in] sphereCenter An array of 3 floating-point: the coordinates of the sphere center in the picture frame
+ * @param[in] ellipseMask The binary mask of the sphere in the picture
+ * @param[in] K Intrinsic parameters of the camera
+ * @param[out] normals Normals on the sphere
+ * @param[out] newMask The mask of the sphere after ray tracing
+ */
+void estimateSphereNormals(const std::array<float, 3>& sphereCenter,
+                           const float sphereRadius,
+                           const Eigen::Matrix3f& K,
+                           image::Image<image::RGBfColor>& normals,
+                           image::Image<float>& newMask);
+
+/**
+ * @brief Estimate the normals of a sphere from a mask
+ * @param[in] maskCV The openCV image of the binary mask of the ellipse in the picture
+ * @param[in] K Intrinsic parameters of the camera
+ * @param[out] normals Normals on the sphere in camera frame
+ * @param[out] newMask The mask of the sphere after ray tracing
+*/
+void getRealNormalOnSphere(const cv::Mat& maskCV,
+                       const Eigen::Matrix3f& K,
+                       const float sphereRadius,
+                       image::Image<image::RGBfColor>& normals,
+                       image::Image<float>& newMask);
+
+
+void getEllipseMaskFromSphereParameters(const std::array<float, 3>& sphereParam, const Eigen::Matrix3f& K, std::array<float, 5>& ellipseParameters, cv::Mat maskCV);
+
+} // namespace lightingEstimation
+} // namespace aliceVision