Category: UV Unwrapping symmetric dirichlet energy with ceres Integrationuv_unwrapping_slim_and_ceres

adding mesh_unwrapper_ceres

6 files changed, 313 insertions, 0 deletions

diff --git a/intern/SLIM/src/ceres_mesh_unwrapper.cpp b/intern/SLIM/src/ceres_mesh_unwrapper.cpp
new file mode 100644
index 00000000000..b0ef67628b5
--- /dev/null
+++ b/intern/SLIM/src/ceres_mesh_unwrapper.cpp
@@ -0,0 +1,302 @@
+//
+//  ceres_mesh_unwrapper.cpp
+//  Blender
+//
+//  Created by Aurel Gruber on 24.02.17.
+//
+//
+
+#include "ceres_mesh_unwrapper.h"
+#include "ceres/ceres.h"
+#include <iostream>
+#include <math.h>       /* cos */
+
+
+using namespace std;
+
+template <typename T>
+T cross_2D(const Eigen::Matrix<T, 1, 2> &a, const Eigen::Matrix<T, 1, 2> &b){
+	return -T(a(0) * b(1) - a(1) * b(0));
+}
+
+/**
+ *find SVD according to http://scicomp.stackexchange.com/questions/8899/robust-algorithm-for-2x2-svd
+ *	-> Alex Eftimiades' answer
+ */
+template <typename T>
+void computeSingularValues(T j11, T j12, T j21, T j22, T& s1, T& s2) {
+
+	T y1 = j12 + j21;
+	T x1 = j11 - j22;
+	T y2 = j12 - j21;
+	T x2 = j11 + j22;
+
+	T h1 = T(sqrt(y1*y1 + x1*x1));
+	T h2 = T(sqrt(y2*y2 + x2*x2));
+
+	s1 = (h1 + h2) / T(2.0);
+	s2 = (h1 - h2) / T(2.0);
+
+	//cout << "s1: " << s1 << endl;
+	//cout << "s2: " << s2 << endl;
+
+}
+
+template <typename T>
+void findJacobian(const double* orig_v2D_A, const double* orig_v2D_B, const double* orig_v2D_C,
+				  T* new_v2D_A, T* new_v2D_B, T* new_v2D_C,
+				  T& j11, T& j12, T& j21, T& j22) {
+	/*
+
+	 To find J we can use the points B and C:
+
+	 J = j11	j12
+	 j21	j22
+
+	 |j11  j12|   |b1|   |b'1|
+	 JB = B'	 ->  |		  | * |  | = |   |
+	 |j21  j22|   |b2|   |b'2|
+
+	 -->
+
+	 |j11  j12|   |c1|   |c'1|
+	 JC = C'	 ->  |        | * |  | = |   |
+	 |j21  j22|   |c2|   |c'2|
+
+
+	 We know, that c2 is zero, due to the nature of the embedding of the original triangle in 2d
+
+	 j11 * b1 + j12 * b2 = b'1
+	 j21 * b1 + j22 * b2 = b'2
+	 j11 * c1 = c'1
+	 j21 * c1 = c'2
+
+	 Thereofore:
+
+	 j11 = c'1 / c1
+	 j21 = c'2 / c1
+
+	 j12 = ( b'1 - j11*b1 ) / b2
+	 j22 = ( b'2 - j21*b1 ) / b2
+
+	 */
+
+	j11 = new_v2D_C[0] / T(orig_v2D_C[0]);
+	j21 = new_v2D_C[1] / T(orig_v2D_C[0]);
+
+	j12 = ( new_v2D_B[0] - j11 * T(orig_v2D_B[0]) ) / T(orig_v2D_B[1]);
+	j22 = ( new_v2D_B[1] - j21 * T(orig_v2D_B[0]) ) / T(orig_v2D_B[1]);
+}
+
+void map_3D_triangles_to_2D_undistorded(
+										const Eigen::Vector3d &v3DA,
+										const Eigen::Vector3d &v3DB,
+										const Eigen::Vector3d &v3DC,
+										Eigen::Vector2d &v2DAembedded,
+										Eigen::Vector2d &v2DBembedded,
+										Eigen::Vector2d &v2DCembedded
+										) {
+	/*
+	     b
+	 A ______ B
+	   \    /
+	  a \  / c
+	     \/
+	     C
+
+	 We place A on 0,0
+	 A = 0,0
+
+	 We place C on the x axis
+	 C = ||C-A||, 0
+
+	 We compute B via angle
+
+	 b = B - A
+	 a = C - A
+
+	 phi =  arccos(a*b)/|a|*|b|
+
+	 B = cos(phi)*b, sin(phi)*b
+	 */
+
+	Eigen::Vector3d xUnitVector = Eigen::Vector3d::UnitX();
+	Eigen::Vector3d yUnitVector = Eigen::Vector3d::UnitY();
+
+	//translate to origin
+	Eigen::Vector3d v3DAembedded;
+	v3DAembedded<< 0,0,0;
+
+	Eigen::Vector3d v3DBembedded = v3DB - v3DA;
+	Eigen::Vector3d v3DCembedded = v3DC - v3DA;
+
+	// rotate triangle, s.t. C is on x axis
+	double angleC_X = std::acos( (v3DCembedded.dot(xUnitVector)) / v3DCembedded.norm() );
+	Eigen::Vector3d axis = v3DCembedded.cross(xUnitVector);
+
+	if (axis.norm() > 0) {
+		axis.normalize();
+		Eigen::AngleAxisd rotationY(angleC_X, axis);
+
+		v3DCembedded = rotationY * v3DCembedded;
+		v3DBembedded = rotationY * v3DBembedded;
+	}
+
+
+	// rotate triangle, such that B is in x-y plane
+	Eigen::Vector3d v3DBembeddedProjectedZYPlane;
+	v3DBembeddedProjectedZYPlane << 0, v3DBembedded(1), v3DBembedded(2);
+	double angleB_XY = std::acos( (v3DBembeddedProjectedZYPlane.dot(yUnitVector)) / v3DBembeddedProjectedZYPlane.norm() );
+	axis = v3DBembeddedProjectedZYPlane.cross(yUnitVector);
+
+	if (axis.norm() > 0) {
+
+		axis.normalize();
+		Eigen::AngleAxisd rotationX(angleB_XY, axis);
+
+		v3DBembedded = rotationX * v3DBembedded;
+	}
+
+	v2DAembedded << v3DAembedded(0), v3DAembedded(1);
+	v2DBembedded << v3DBembedded(0), v3DBembedded(1);
+	v2DCembedded << v3DCembedded(0), v3DCembedded(1);
+}
+
+struct DistortionResidual {
+	DistortionResidual(const double* v2D_A, const double* v2D_B, const double* v2D_C)
+	: orig_v2D_A(v2D_A), orig_v2D_B(v2D_B), orig_v2D_C(v2D_C) {}
+	template <typename T>
+	bool operator()(
+					const T* const new_v2D_A_offset,
+					const T* const new_v2D_B_offset,
+					const T* const new_v2D_C_offset,
+					T* residuals) const {
+
+
+		/*
+		 The Jacobian J represents the mapping from the embedded, undistorted triangles to the distorted triangles of the current iterate:
+
+            A ______ B         A'__________ B'
+              \    /	  J		|	    .
+               \  /		-->     |.  ^
+                \/				C'
+                 C
+
+		 Given an ortogonal frame per triangle, we can express any point in this local coordinatesystem. We first place the triangle of the current iterate at the origin:
+		 */
+
+		T new_v2D_A[2];
+		T new_v2D_B[2];
+		T new_v2D_C[2];
+
+		new_v2D_A[0] = T(0);
+		new_v2D_A[1] = T(0);
+		new_v2D_B[0] = new_v2D_B_offset[0] - new_v2D_A_offset[0];
+		new_v2D_B[1] = new_v2D_B_offset[1] - new_v2D_A_offset[1];
+		new_v2D_C[0] = new_v2D_C_offset[0] - new_v2D_A_offset[0];
+		new_v2D_C[1] = new_v2D_C_offset[1] - new_v2D_A_offset[1];
+
+		/*
+		 The unwrapping can be characterised as a linear function w.r.t. the undistorted triangle. Given a point p1 "on" the undistorted triangle, it must hold that
+
+		 J*p1 = q1
+
+		 where q1 is the position of the point after the mapping. Ideally, no distortion is induced. This means J == R, that is some rotation matrix.
+
+		 That means, the singular values of J should all be 1.
+		 */
+
+		T j11, j12, j21, j22;
+		findJacobian(orig_v2D_A, orig_v2D_B, orig_v2D_C,
+					 new_v2D_A, new_v2D_B, new_v2D_C,
+					 j11, j12, j21, j22);
+
+		/*
+			Given J, we can find the singular values s1, s2:
+		 */
+
+		T s1, s2;
+		computeSingularValues(j11, j12, j21, j22, s1, s2);
+
+
+		const Eigen::Map<const Eigen::Matrix<T, 1, 2> > ev1(new_v2D_A_offset);
+		const Eigen::Map<const Eigen::Matrix<T, 1, 2> > ev2(new_v2D_B_offset);
+		const Eigen::Map<const Eigen::Matrix<T, 1, 2> > ev3(new_v2D_C_offset);
+
+		Eigen::Matrix<T, 1, 2> e21 = ev1 - ev2;
+		Eigen::Matrix<T, 1, 2> e23 = ev3 - ev2;
+
+		T area_2 = cross_2D(e21, e23);
+
+		if (area_2 < T(0.0)){
+			cout << "area negative! " << endl;
+			return false;
+		}
+
+		residuals[0] = s1;
+		residuals[1] = s2;
+		residuals[2] = T(1.0) / s1;
+		residuals[3] = T(1.0) / s2;
+
+		return true;
+	}
+	const double* orig_v2D_A;
+	const double* orig_v2D_B;
+	const double* orig_v2D_C;
+};
+
+bool solve_map_with_ceres(const Eigen::MatrixXd &Vertex3D, const Eigen::MatrixXi &FaceIndices, Eigen::MatrixXd &UV, int nIterations){
+	int num_faces = FaceIndices.rows();
+
+	Eigen::Matrix<double, Eigen::Dynamic, 2, Eigen::RowMajor> vertices2d = UV;
+	Eigen::Matrix<double, Eigen::Dynamic, 3, Eigen::RowMajor> vertices3d = Vertex3D;
+	Eigen::Matrix<double, Eigen::Dynamic, 2, Eigen::RowMajor> vertices2DEmbedded(num_faces*3, 2);
+
+	ceres::Problem problem;
+	for (int i = 0; i < num_faces; ++i) {
+
+		// computed undistorted embedding
+		Eigen::Vector2d vertex2DAEmbedded;
+		Eigen::Vector2d vertex2DBEmbedded;
+		Eigen::Vector2d vertex2DCEmbedded;
+
+		map_3D_triangles_to_2D_undistorded(vertices3d.row(FaceIndices(i, 0)),
+										   vertices3d.row(FaceIndices(i, 1)),
+										   vertices3d.row(FaceIndices(i, 2)),
+										   vertex2DAEmbedded,
+										   vertex2DBEmbedded,
+										   vertex2DCEmbedded);
+
+		vertices2DEmbedded.row(3*i) = vertex2DAEmbedded;
+		vertices2DEmbedded.row(3*i+1) = vertex2DBEmbedded;
+		vertices2DEmbedded.row(3*i+2) = vertex2DCEmbedded;
+
+		problem.AddResidualBlock(
+			new ceres::AutoDiffCostFunction<DistortionResidual, 4, 2, 2, 2>(
+				new DistortionResidual(
+						&vertices2DEmbedded(3*i, 0),
+						&vertices2DEmbedded(3*i+1, 0),
+						&vertices2DEmbedded(3*i+2, 0)
+				)
+			),
+			NULL,
+			&vertices2d(FaceIndices(i, 0), 0),
+			&vertices2d(FaceIndices(i, 1), 0),
+			&vertices2d(FaceIndices(i, 2), 0));
+	}
+
+	// Make Ceres automatically detect the bundle structure. Note that the
+	// standard solver, SPARSE_NORMAL_CHOLESKY, also works fine but it is slower
+	// for standard bundle adjustment problems.
+	ceres::Solver::Options options;
+	options.max_num_iterations = nIterations;
+	options.linear_solver_type = ceres::CGNR;
+	options.sparse_linear_algebra_library_type = ceres::EIGEN_SPARSE;
+	options.minimizer_progress_to_stdout = true;
+	ceres::Solver::Summary summary;
+	ceres::Solve(options, &problem, &summary);
+	std::cout << summary.FullReport() << "\n";
+	
+	UV = vertices2d;
+	return true;
+}
diff --git a/intern/SLIM/src/ceres_mesh_unwrapper.h b/intern/SLIM/src/ceres_mesh_unwrapper.h
new file mode 100644
index 00000000000..59a0ca4b5bb
--- /dev/null
+++ b/intern/SLIM/src/ceres_mesh_unwrapper.h
@@ -0,0 +1,11 @@
+//
+//  ceres_mesh_unwrapper.hpp
+//  Blender
+//
+//  Created by Aurel Gruber on 24.02.17.
+//
+//
+
+#pragma once
+#include <Eigen/Dense>
+bool solve_map_with_ceres(const Eigen::MatrixXd &Vertex3D, const Eigen::MatrixXi &FaceIndices, Eigen::MatrixXd &UV, int nIterations);
diff --git a/release/datafiles/locale b/release/datafiles/locale
-Subproject dc166057192ea882b5cc70484d4c8bacd7cb41b
+Subproject 507eacde9608ce0190c6087fb338dd63f7a1649
diff --git a/release/scripts/addons b/release/scripts/addons
-Subproject 06dad53c80801e0e0919f086040e3d9c31bbd0a
+Subproject ff8ef5cb121e590ef1a5f3e5f57b9bf99c46559
diff --git a/release/scripts/addons_contrib b/release/scripts/addons_contrib
-Subproject 04af69be141a5757fc60b44cc1a5b72db524af3
+Subproject e0e0fe92738e345746ef566b6ecfb79d28a28d6
diff --git a/source/tools b/source/tools
-Subproject 896c5f78952adb2d091d28c65086d46992dabda
+Subproject b11375e89061303401376f7aeae42ac2fd64692