bni normal integration optimized for 0 discontinuity, few small bug fixed.

Author: ponchio

relight-cli/rtibuilder.cpp

@@ -1068,6 +1068,8 @@ size_t RtiBuilder::savePTM(const std::string &output) {
 	//second reading.
 	imageset.restart();
 
+	nworkers = QThread::idealThreadCount();
+
 	vector<Worker *> workers(height, nullptr);
 	for(size_t i = 0; i < nworkers; i++) {
 		workers[i] = new Worker(*this);

relightlab/CMakeLists.txt

@@ -15,13 +15,14 @@ set(CMAKE_AUTOUIC ON)
 set(CMAKE_AUTORCC ON)
 
 set(OPENCV_RELIGHT_LIBS opencv_core  opencv_imgcodecs opencv_imgproc opencv_video)
-	
+add_compile_definitions(WITH_OPENCV)
+
 if (APPLE)
 	set(CMAKE_OSX_DEPLOYMENT_TARGET "11.0" CACHE STRING "Minimum OS X deployment version" FORCE)
 	SET(CMAKE_INSTALL_RPATH $ORIGIN/../Frameworks)
 	set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
 #else()
-#	add_compile_definitions(WITH_OPENCV)
+#	
 endif()
 
 

relightlab/main.cpp

@@ -1,16 +1,16 @@
+
+#include "relightapp.h"
+#include "mainwindow.h"
+
 #include <QApplication>
 
 #include <QVBoxLayout>
 #include <QPushButton>
 #include <QLabel>
 #include <QImageReader>
 
-#include "relightapp.h"
-#include "mainwindow.h"
-
-#include <locale.h>
-
 #include <iostream>
+#include <locale.h>
 using namespace std;
 
 #define RELIGHT_STRINGIFY0(v) #v

relightlab/normalsplan.cpp

@@ -164,7 +164,7 @@ NormalsSurfaceRow::NormalsSurfaceRow(NormalsParameters &_parameters, QFrame *par
 			QHBoxLayout *bni_parameter = new QHBoxLayout;
 			bni_parameter->addWidget(new QLabel("Discontinuity propensity."));
 			bni_parameter->addWidget(bni_k = new QDoubleSpinBox);
-			bni_k->setRange(0.01, 50);
+			bni_k->setRange(0.00, 50);
 			bni_k->setValue(parameters.bni_k);
 
 			bni_layout->addLayout(bni_parameter);

relightlab/processqueue.h

@@ -19,7 +19,6 @@ class ProcessQueue: public QThread {
 	Q_OBJECT
 public:
 	~ProcessQueue();
-	int nthreads = 4;
 
 	Task *task = nullptr; //task is removed from the queue when executing.
 	bool stopped = false;

relightlab/relightlab.pro

@@ -6,7 +6,7 @@ CONFIG += c++17
 #target_link_libraries( ${APP_NAME} Qt5::Svg )
 
 DEFINES += QT_DEPRECATED_WARNINGS
-DEFINES += _USE_MATH_DEFINES #WITH_OPENCV
+DEFINES += _USE_MATH_DEFINES WITH_OPENCV
 DEFINES += NOMINMAX
 
 INCLUDEPATH += ../external/
@@ -16,17 +16,17 @@ win32:INCLUDEPATH += ../external/libjpeg-turbo-2.0.6/include \
     ../src/
 win32:LIBS += ../external/libjpeg-turbo-2.0.6/lib/jpeg-static.lib
 
+unix::QMAKE_CXXFLAGS = -fopenmp
 unix:INCLUDEPATH += ../external/eigen-3.3.9/ /usr/include/opencv4
-unix:LIBS += -ljpeg -ltiff -lgomp
-#unix:LIBS += -lopencv_core -lopencv_imgcodecs -lopencv_imgproc -lopencv_video
-
+unix:LIBS += -ljpeg -ltiff
+unix:LIBS += -lopencv_core -lopencv_imgcodecs -lopencv_imgproc -lopencv_video
+unix::LIBS += -fopenmp #-lgomp
 
 mac:INCLUDEPATH += /usr/local/Cellar/jpeg-turbo/3.1.0/include \
     /usr/local/include \
     /usr/local/include/eigen3
 mac:LIBS += -L/usr/local/Cellar/jpeg-turbo/3.1.0/lib/ -ljpeg
 mac:LIBS += -framework Accelerate
-#mac:QMAKE_CXXFLAGS += -fopenmp
 mac:QMAKE_CXXFLAGS += -Xpreprocessor -I/usr/local/include #fopenmp -lomp
 # mac:QMAKE_LFLAGS += -lomp
 mac:LIBS += -L /usr/local/lib /usr/local/lib/libomp.dylib

src/bni_normal_integration.cpp

@@ -283,6 +283,92 @@ std::vector<float> bni_pyramid(std::function<bool(QString s, int n)> progressed,
 	h = result.h;
 	return result.heights;
 }
+#include <Eigen/SparseCholesky>
+
+
+Eigen::VectorXd bni_integrate_iterative(std::function<bool(QString s, int n)> progressed, Eigen::SparseMatrix<double> &A, Eigen::VectorXd &b, Eigen::SparseMatrix<double> W,
+							 double k,
+							 double tolerance, double solver_tolerance,
+							 int max_iterations, int max_solver_iterations) {
+
+	int n = b.size()/4;
+
+	Eigen::VectorXd z = Eigen::VectorXd::Zero(n);
+
+	Eigen::MatrixXd tmp = A*z - b;
+	double energy = (tmp.transpose() * W * tmp)(0);
+	double start_energy = energy;
+	if(isnan(energy)) {
+		throw "Computational problems.";
+	}
+
+	cout << "Energy : " << energy << endl;
+
+	vector<Triple> triples;
+	for(int i = 0; i < max_iterations; i++) {
+
+		Eigen::SparseMatrix<double> A_mat = A.transpose()*W*A;
+		Eigen::VectorXd b_vec = A.transpose()*W*b;
+
+		Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
+
+		solver.compute(A_mat);
+		solver.setTolerance(solver_tolerance);
+		solver.setMaxIterations(max_solver_iterations);
+		z = solver.solveWithGuess(b_vec, z);
+		if (solver.info() != Eigen::Success) {
+			double finalResidual = solver.error();
+			cerr << "Max iter reached with error: " << finalResidual << endl;
+		}
+
+		// Get the number of iterations
+		int numIterations = solver.iterations();
+		std::cout << "Number of iterations: " << numIterations << " error: " << solver.error() << " tolerance: " << solver_tolerance << std::endl;
+
+
+		if(k == 0)
+			break;
+
+		Eigen::VectorXd wu = ((A.block(n, 0, n, n)*z).array().pow(2) -
+							  (A.block(0, 0, n, n)*z).array().pow(2)).unaryExpr([k](double x) { return sigmoid(x, k); });
+		Eigen::VectorXd wv = ((A.block(n*3, 0, n, n)*z).array().pow(2) -
+							  (A.block(n*2, 0, n, n)*z).array().pow(2)).unaryExpr([k](double x) { return sigmoid(x, k); });
+
+		triples.clear();
+		for(int i = 0; i < n; i++) {
+			triples.push_back(Triple(i, i, wu(i)));
+			triples.push_back(Triple(i+n, i+n, 1.0 - wu(i)));
+			triples.push_back(Triple(i+n*2, i+n*2, wv(i)));
+			triples.push_back(Triple(i+n*3, i+n*3, 1 - wv(i)));
+		}
+		W.setFromTriplets(triples.begin(), triples.end());
+
+		double energy_old = energy;
+		tmp = A*z - b;
+		energy = (tmp.transpose() * W * tmp)(0);
+		cout << "Energy: " << energy << endl;
+
+		double relative_energy = fabs(energy - energy_old) / energy_old;
+		double total_progress = fabs(energy - start_energy) / start_energy;
+		if(progressed) {
+			bool proceed = progressed("Integrating normals...", 100*(1 - (log(relative_energy) - log(tolerance))/(log(total_progress) - log(tolerance))));
+			if(!proceed) break;
+		}
+		if(relative_energy < tolerance)
+			break;
+	}
+	return z;
+}
+
+Eigen::VectorXd bni_integrate_direct(Eigen::SparseMatrix<double> &A, Eigen::VectorXd &b, Eigen::SparseMatrix<double> &W) {
+	Eigen::SparseMatrix<double> A_mat = A.transpose()*W*A;
+	Eigen::VectorXd b_vec = A.transpose()*W*b;
+
+	Eigen::SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+	solver.compute(A_mat);
+	return solver.solve(b_vec);
+}
+
 
 void bni_integrate(std::function<bool(QString s, int n)> progressed, int w, int h, std::vector<float> &normalmap, std::vector<float> &heights,
 				   double k,
@@ -361,77 +447,14 @@ void bni_integrate(std::function<bool(QString s, int n)> progressed, int w, int
 		triples.push_back(Triple(i, i, 0.5));
 	W.setFromTriplets(triples.begin(), triples.end());
 
-	Eigen::VectorXd z = Eigen::VectorXd::Zero(n);
-	//for(int i = 0; i < n; i++)
-	//	z(i) = heights[i];
-
-	Eigen::MatrixXd tmp = A*z - b;
-	double energy = (tmp.transpose() * W * tmp)(0);
-	double start_energy = energy;
-	if(isnan(energy)) {
-		throw "Accidentaccio!";
-	}
-
-	cout << "Energy : " << energy << endl;
-
-	for(int i = 0; i < max_iterations; i++) {
-
-		Eigen::SparseMatrix<double> A_mat = A.transpose()*W*A;
-		Eigen::VectorXd b_vec = A.transpose()*W*b;
-
-		Eigen::ConjugateGradient<Eigen::SparseMatrix<double>> solver;
-
-		//Eigen::LeastSquaresConjugateGradient<Eigen::SparseMatrix<double>> solver;
-		//Eigen::BiCGSTAB<Eigen::SparseMatrix<double> > solver;
-		solver.compute(A_mat);
-		solver.setTolerance(solver_tolerance);
-		solver.setMaxIterations(max_solver_iterations);
-
-		z = solver.solveWithGuess(b_vec, z);
-		if (solver.info() != Eigen::Success) {
-			double finalResidual = solver.error();
-			cerr << "Max iter reached with error: " << finalResidual << endl;
-			//			return std::vector<double>();
-		}
-
-		// Get the number of iterations
-		int numIterations = solver.iterations();
-		std::cout << "Number of iterations: " << numIterations << " error: " << solver.error() << " tolerance: " << solver_tolerance << std::endl;
-
+	Eigen::VectorXd z;
 
-		if(k == 0)
-			break;
+	if(k == 0.0)
+		z = bni_integrate_direct(A, b, W);
+	else
+		z = bni_integrate_iterative(progressed, A, b, W, k, tolerance, solver_tolerance, max_iterations, max_solver_iterations);
 
-		Eigen::VectorXd wu = ((A.block(n, 0, n, n)*z).array().pow(2) -
-							  (A.block(0, 0, n, n)*z).array().pow(2)).unaryExpr([k](double x) { return sigmoid(x, k); });
-		Eigen::VectorXd wv = ((A.block(n*3, 0, n, n)*z).array().pow(2) -
-							  (A.block(n*2, 0, n, n)*z).array().pow(2)).unaryExpr([k](double x) { return sigmoid(x, k); });
-
-		triples.clear();
-		for(int i = 0; i < n; i++) {
-			triples.push_back(Triple(i, i, wu(i)));
-			triples.push_back(Triple(i+n, i+n, 1.0 - wu(i)));
-			triples.push_back(Triple(i+n*2, i+n*2, wv(i)));
-			triples.push_back(Triple(i+n*3, i+n*3, 1 - wv(i)));
-		}
-		W.setFromTriplets(triples.begin(), triples.end());
-
-		double energy_old = energy;
-		tmp = A*z - b;
-		energy = (tmp.transpose() * W * tmp)(0);
-		cout << "Energy: " << energy << endl;
-
-		double relative_energy = fabs(energy - energy_old) / energy_old;
-		double total_progress = fabs(energy - start_energy) / start_energy;
-		if(progressed) {
-			bool proceed = progressed("Integrating normals...", 100*(1 - (log(relative_energy) - log(tolerance))/(log(total_progress) - log(tolerance))));
-			if(!proceed) break;
-		}
-		if(relative_energy < tolerance)
-			break;
-	}
 	heights.resize(w*h);
 	for(int i = 0; i < w*h; i++)
 		heights[i] = z(i);
-	//memcpy(&depthmap[0], z.data(), n*8);
 }

src/bni_normal_integration.h

@@ -15,15 +15,15 @@ void bni_integrate(std::function<bool(QString s, int n)> progressed,
 								  double tolerance = 1e-5,
 								  double solver_tolerance = 1e-5,
 								  int max_iterations = 10,
-								  int max_solver_iterations = 5000);
+								  int max_solver_iterations = 500);
 
 std::vector<float> bni_pyramid(std::function<bool(QString s, int n)> progressed,
 								  int &w, int &h, std::vector<float> &normalmap,
 								  double k = 2.0,
 								  double tolerance = 1e-5,
 								  double solver_tolerance = 1e-5,
-								  int max_iterations = 150,
-								  int max_solver_iterations = 5000,
+								  int max_iterations = 15,
+								  int max_solver_iterations = 500,
 								  int scale = 0);
 
 bool savePly(const QString &filename, size_t w, size_t h, std::vector<float> &z);

src/dome.cpp

@@ -69,6 +69,7 @@ void Dome::fromSpheres(std::vector<Image> &images, std::vector<Sphere *> &sphere
 	case Dome::LIGHTS3D:
 		computeParallaxPositions(images, spheres, lens, positions3d);
 		directions = positions3d;
+		positionsSphere.resize(directions.size());
 		for(size_t i = 0; i < directions.size(); i++) {
 			float len = directions[i].norm();
 			directions[i] /= len;

src/project.cpp

@@ -280,6 +280,12 @@ void Project::load(QString filename) {
 	if(!setDir(folder))
 		throw QString("The folder " + obj["folder"].toString() + " does not exists.");
 
+	//ensure resources folder has been created (for old project)
+	QDir resources = dir.filePath("resources");
+	if (!resources.exists() && !dir.mkpath("resources")) {
+		throw QString("Could not create the resources folder in %1").arg(dir.absolutePath());
+	}
+
 	lens.fromJson(obj["lens"].toObject());
 	lens.width = imgsize.width();
 	lens.height = imgsize.height();