Skip to content

Commit

Permalink
Merge pull request #1590 from alicevision/dev/ps_sh
Browse files Browse the repository at this point in the history
[Photometric Stereo] MultiView fusion in Texturing
  • Loading branch information
fabiencastan authored Sep 6, 2024
2 parents 084b4c6 + 977f4cc commit 54cdaca
Show file tree
Hide file tree
Showing 31 changed files with 1,334 additions and 291 deletions.
4 changes: 2 additions & 2 deletions src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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()
Expand Down
4 changes: 2 additions & 2 deletions src/aliceVision/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -51,7 +51,6 @@ endif()
# MVS modules

if(ALICEVISION_BUILD_MVS)
add_subdirectory(lightingEstimation)
add_subdirectory(mesh)
add_subdirectory(mvsData)
add_subdirectory(mvsUtils)
Expand All @@ -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()

Expand Down
2 changes: 1 addition & 1 deletion src/aliceVision/fuseCut/DelaunayGraphCut_test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand Down
4 changes: 4 additions & 0 deletions src/aliceVision/lightingEstimation/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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
)


Expand Down
4 changes: 2 additions & 2 deletions src/aliceVision/lightingEstimation/augmentedNormals.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -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); }
Expand Down
232 changes: 232 additions & 0 deletions src/aliceVision/lightingEstimation/ellipseGeometry.cpp
Original file line number Diff line number Diff line change
@@ -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
75 changes: 75 additions & 0 deletions src/aliceVision/lightingEstimation/ellipseGeometry.hpp
Original file line number Diff line number Diff line change
@@ -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
Loading

0 comments on commit 54cdaca

Please sign in to comment.