Lagrange
compute_tangent_bitangent.h
1/*
2 * Copyright 2020 Adobe. All rights reserved.
3 * This file is licensed to you under the Apache License, Version 2.0 (the "License");
4 * you may not use this file except in compliance with the License. You may obtain a copy
5 * of the License at http://www.apache.org/licenses/LICENSE-2.0
6 *
7 * Unless required by applicable law or agreed to in writing, software distributed under
8 * the License is distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR REPRESENTATIONS
9 * OF ANY KIND, either express or implied. See the License for the specific language
10 * governing permissions and limitations under the License.
11 */
12#pragma once
13
14#include <lagrange/Mesh.h>
15#include <lagrange/MeshTrait.h>
16#include <lagrange/chain_corners_around_edges.h>
17#include <lagrange/chain_corners_around_vertices.h>
18#include <lagrange/common.h>
19#include <lagrange/corner_to_edge_mapping.h>
20#include <lagrange/legacy/inline.h>
21#include <lagrange/utils/DisjointSets.h>
22#include <lagrange/utils/geometry3d.h>
23
24// clang-format off
25#include <lagrange/utils/warnoff.h>
26#include <tbb/parallel_for.h>
27#include <lagrange/utils/warnon.h>
28// clang-format on
29
30#include <cmath>
31
32namespace lagrange {
33LAGRANGE_LEGACY_INLINE
34namespace legacy {
35
36namespace tangent_bitangent_detail {
37
49template <typename Scalar>
50Eigen::Matrix<Scalar, 2, 3> triangle_tangent_bitangent(
51 const Eigen::Matrix<Scalar, 3, 3>& vertices,
52 const Eigen::Matrix<Scalar, 3, 2>& uvs,
53 bool& sign)
54{
55 const auto edge0 = (vertices.row(1) - vertices.row(0)).eval();
56 const auto edge1 = (vertices.row(2) - vertices.row(0)).eval();
57
58 const auto edge_uv0 = (uvs.row(1) - uvs.row(0)).eval();
59 const auto edge_uv1 = (uvs.row(2) - uvs.row(0)).eval();
60
61 const Scalar uv_signed_area = (edge_uv0.x() * edge_uv1.y() - edge_uv0.y() * edge_uv1.x());
62 const Scalar r = (uv_signed_area == Scalar(0) ? Scalar(1) : Scalar(1) / uv_signed_area);
63
64 sign = (uv_signed_area > 0);
65
66 Eigen::Matrix<Scalar, 2, 3> T_BT;
67 T_BT.row(0) = ((edge0 * edge_uv1.y() - edge1 * edge_uv0.y()) * r).stableNormalized();
68 T_BT.row(1) = ((edge1 * edge_uv0.x() - edge0 * edge_uv1.x()) * r).stableNormalized();
69 return T_BT;
70}
71
81template <typename Scalar>
82bool triangle_uv_orientation(const Eigen::Matrix<Scalar, 3, 2>& uvs)
83{
84 const auto edge_uv0 = (uvs.row(1) - uvs.row(0)).eval();
85 const auto edge_uv1 = (uvs.row(2) - uvs.row(0)).eval();
86
87 const Scalar uv_signed_area = (edge_uv0.x() * edge_uv1.y() - edge_uv0.y() * edge_uv1.x());
88
89 return uv_signed_area > 0;
90}
91
104template <typename Scalar>
105Eigen::Matrix<Scalar, 2, 3> quad_tangent_bitangent(
106 const Eigen::Matrix<Scalar, 4, 3>& vertices,
107 const Eigen::Matrix<Scalar, 4, 2>& uvs,
108 bool& sign)
109{
110 // Check whether there are colocated vertices in UV space
111 Eigen::Vector3i indices = Eigen::Vector3i(0, 1, 2);
112 if ((uvs.row(0) - uvs.row(1)).isZero()) {
113 indices = Eigen::Vector3i(1, 2, 3);
114 } else if ((uvs.row(1) - uvs.row(2)).isZero()) {
115 indices = Eigen::Vector3i(2, 3, 0);
116 } else if ((uvs.row(2) - uvs.row(3)).isZero()) {
117 indices = Eigen::Vector3i(3, 0, 1);
118 }
119
120 Eigen::Matrix<Scalar, 3, 3> triangle_vertices;
121 triangle_vertices.row(0) = vertices.row(indices(0));
122 triangle_vertices.row(1) = vertices.row(indices(1));
123 triangle_vertices.row(2) = vertices.row(indices(2));
124
125 Eigen::Matrix<Scalar, 3, 2> triangle_uvs;
126 triangle_uvs.row(0) = uvs.row(indices(0));
127 triangle_uvs.row(1) = uvs.row(indices(1));
128 triangle_uvs.row(2) = uvs.row(indices(2));
129
130 // Calculate from triangle
131 return triangle_tangent_bitangent(triangle_vertices, triangle_uvs, sign);
132}
133
148template <typename MeshType, typename DerivedT, typename DerivedB>
149void accumulate_tangent_bitangent(
150 MeshType& mesh,
151 const Eigen::MatrixBase<DerivedT>& tangents,
152 const Eigen::MatrixBase<DerivedB>& bitangents)
153{
154 static_assert(MeshTrait<MeshType>::is_mesh(), "Input type is not Mesh");
155 using Scalar = typename MeshType::Scalar;
156 using Index = typename MeshType::Index;
157 using FacetArray = typename MeshType::FacetArray;
158 using AttributeArray = typename MeshType::AttributeArray;
159 using IndexArray = Eigen::Matrix<Index, Eigen::Dynamic, 1>;
160 using RowVector3r = Eigen::Matrix<Scalar, 1, 3>;
161
162 la_runtime_assert(
163 mesh.get_vertex_per_facet() == 3,
164 "Only triangle meshes are supported for this.");
165 la_runtime_assert(safe_cast<Index>(tangents.rows()) == mesh.get_num_facets() * 3);
166 la_runtime_assert(safe_cast<Index>(bitangents.rows()) == mesh.get_num_facets() * 3);
167 la_runtime_assert(tangents.cols() == bitangents.cols());
168 la_runtime_assert(tangents.cols() == 3 || tangents.cols() == 4);
169
170 // Compute edge information
171 logger().trace("Corner to edge mapping");
172 IndexArray c2e;
173 corner_to_edge_mapping(mesh.get_facets(), c2e);
174 logger().trace("Chain corners around edges");
175 IndexArray e2c;
176 IndexArray next_corner_around_edge;
177 chain_corners_around_edges(mesh.get_facets(), c2e, e2c, next_corner_around_edge);
178 logger().trace("Chain corners around vertices");
179 IndexArray v2c;
180 IndexArray next_corner_around_vertex;
181 chain_corners_around_vertices(
182 mesh.get_num_vertices(),
183 mesh.get_facets(),
184 v2c,
185 next_corner_around_vertex);
186
187 const auto nvpf = mesh.get_vertex_per_facet();
188 const auto& vertices = mesh.get_vertices();
189 const auto& facets = mesh.get_facets();
190 const auto& uv_values = std::get<0>(mesh.get_indexed_attribute("uv"));
191 const auto& uv_indices = std::get<1>(mesh.get_indexed_attribute("uv"));
192 const auto& normal_values = std::get<0>(mesh.get_indexed_attribute("normal"));
193 const auto& normal_indices = std::get<1>(mesh.get_indexed_attribute("normal"));
194
195 // Check the uv triangle orientation
196 logger().trace("Is orient preserving?");
197 std::vector<unsigned char> is_orient_preserving(mesh.get_num_facets());
198 tbb::parallel_for(Index(0), mesh.get_num_facets(), [&](Index f) {
199 Eigen::Matrix<Scalar, 3, 2> uvs;
200 uvs.row(0) = uv_values.row(uv_indices(f, 0));
201 uvs.row(1) = uv_values.row(uv_indices(f, 1));
202 uvs.row(2) = uv_values.row(uv_indices(f, 2));
203 is_orient_preserving[f] = triangle_uv_orientation(uvs);
204 });
205
206 // Iterate over facets incident to eij until fi and fj are found. The edge is considered smooth
207 // if both corners at (eij, fi) and (eij, fj) share the same uv and normal indices.
208 auto is_edge_smooth = [&](Index eij, Index fi, Index fj) {
209 std::array<Index, 2> v, n, uv;
210 std::fill(v.begin(), v.end(), invalid<Index>());
211 for (Index c = e2c[eij]; c != invalid<Index>(); c = next_corner_around_edge[c]) {
212 Index f = c / 3;
213 Index lv = c % 3;
214 if (f == fi || f == fj) {
215 Index v0 = facets(f, lv);
216 Index v1 = facets(f, (lv + 1) % 3);
217 Index n0 = normal_indices(f, lv);
218 Index n1 = normal_indices(f, (lv + 1) % 3);
219 Index uv0 = uv_indices(f, lv);
220 Index uv1 = uv_indices(f, (lv + 1) % 3);
221 if (v[0] == invalid<Index>()) {
222 v[0] = v0;
223 v[1] = v1;
224 n[0] = n0;
225 n[1] = n1;
226 uv[0] = uv0;
227 uv[1] = uv1;
228 continue;
229 }
230 if (v0 != v[0]) {
231 std::swap(v0, v1);
232 std::swap(n0, n1);
233 std::swap(uv0, uv1);
234 }
235 assert(v0 == v[0]);
236 return (n0 == n[0] && n1 == n[1] && uv0 == uv[0] && uv1 == uv[1]);
237 }
238 }
239 assert(false);
240 return false;
241 };
242
243 // Check if two face vertices are collapsed
244 auto is_face_degenerate = [&](Index f) {
245 for (Index lv = 0; lv < nvpf; ++lv) {
246 if (facets(f, lv) == facets(f, (lv + 1) % nvpf)) {
247 return true;
248 }
249 }
250 return false;
251 };
252
253 // STEP 1: For each vertex v, iterate over each pair of incident facet (fi, fj) that share a
254 // common edge eij.
255 logger().trace("Loop to unify corner indices");
256 DisjointSets<Index> unified_indices((Index)tangents.rows());
257 tbb::parallel_for(Index(0), Index(mesh.get_num_vertices()), [&](Index v) {
258 for (Index ci = v2c[v]; ci != invalid<Index>(); ci = next_corner_around_vertex[ci]) {
259 Index eij = c2e[ci];
260 Index fi = ci / 3;
261 Index lvi = ci % 3;
262 Index vi = facets(fi, lvi);
263 if (is_face_degenerate(fi)) continue;
264
265 auto si = is_orient_preserving[fi];
266 for (Index cj = e2c[eij]; cj != invalid<Index>(); cj = next_corner_around_edge[cj]) {
267 Index fj = cj / 3;
268 Index lvj = cj % 3;
269 auto sj = is_orient_preserving[fj];
270 if (fi == fj) continue;
271 if (si != sj) continue;
272 Index vj = facets(fj, lvj);
273 if (vi != vj) {
274 lvj = (lvj + 1) % 3;
275 vj = facets(fj, lvj);
276 assert(vi == vj);
277 }
278
279 if (is_edge_smooth(eij, fi, fj)) {
280 unified_indices.merge(ci, fj * 3 + lvj);
281 }
282 }
283 }
284 });
285
286 // STEP 2: Perform averaging and reindex attribute
287 logger().trace("Compute new indices");
288 const Index num_corners = (Index)tangents.rows();
289 std::vector<Index> repr(num_corners, invalid<Index>());
290 Index num_indices = 0;
291 for (Index n = 0; n < num_corners; ++n) {
292 Index r = unified_indices.find(n);
293 if (repr[r] == invalid<Index>()) {
294 repr[r] = num_indices++;
295 }
296 repr[n] = repr[r];
297 }
298 logger().trace("Compute offsets");
299 std::vector<Index> indices(num_corners);
300 std::vector<Index> offsets(num_indices + 1, 0);
301 for (Index c = 0; c < num_corners; ++c) {
302 offsets[repr[c] + 1]++;
303 }
304 for (Index r = 1; r <= num_indices; ++r) {
305 offsets[r] += offsets[r - 1];
306 }
307 {
308 // Bucket sort for corner indices
309 std::vector<Index> counter = offsets;
310 for (Index c = 0; c < num_corners; ++c) {
311 indices[counter[repr[c]]++] = c;
312 }
313 }
314
315 logger().trace("Project and average tangent/bitangent");
316 FacetArray tangent_indices(mesh.get_num_facets(), 3);
317 AttributeArray tangent_values(num_indices, tangents.cols());
318 AttributeArray bitangent_values(num_indices, bitangents.cols());
319 tangent_values.setZero();
320 bitangent_values.setZero();
321 tbb::parallel_for(Index(0), Index(offsets.size() - 1), [&](Index r) {
322 for (Index i = offsets[r]; i < offsets[r + 1]; ++i) {
323 Index c = indices[i];
324 Index f = c / 3;
325 Index lv = c % 3;
326 assert(repr[c] == r);
327
328 // Project tangent frame onto the tangent plane
329 RowVector3r nrm = normal_values.row(normal_indices(f, lv));
330 RowVector3r t = tangents.row(c).template head<3>();
331 RowVector3r bt = bitangents.row(c).template head<3>();
332 t = project_on_plane(t, nrm).stableNormalized();
333 bt = project_on_plane(bt, nrm).stableNormalized();
334
335 // Compute weight as the angle between the (projected) edges
336 RowVector3r v0 = vertices.row(facets(f, lv));
337 RowVector3r v1 = vertices.row(facets(f, (lv + 1) % 3));
338 RowVector3r v2 = vertices.row(facets(f, (lv + 2) % 3));
339 RowVector3r e01 = v1 - v0;
340 RowVector3r e02 = v2 - v0;
341 Scalar w = projected_angle_between(e01, e02, nrm);
342
343 tangent_values.row(r).template head<3>() += t * w;
344 bitangent_values.row(r).template head<3>() += bt * w;
345 if (tangents.cols() == 4) {
346 if (tangent_values(r, 3) == 0) {
347 tangent_values(r, 3) = tangents(c, 3);
348 bitangent_values(r, 3) = bitangents(c, 3);
349 } else {
350 assert(tangent_values(r, 3) == tangents(c, 3));
351 assert(bitangent_values(r, 3) == bitangents(c, 3));
352 }
353 }
354
355 tangent_indices(f, lv) = r;
356 }
357 });
358 logger().trace("Normalize vectors");
359 tbb::parallel_for(Index(0), Index(tangent_values.rows()), [&](Index c) {
360 tangent_values.row(c).template head<3>().stableNormalize();
361 bitangent_values.row(c).template head<3>().stableNormalize();
362 });
363
364 logger().trace("Write attributes");
365 FacetArray bitangent_indices = tangent_indices;
366 mesh.add_indexed_attribute("tangent");
367 mesh.add_indexed_attribute("bitangent");
368 mesh.import_indexed_attribute("tangent", tangent_values, tangent_indices);
369 mesh.import_indexed_attribute("bitangent", bitangent_values, bitangent_indices);
370}
371
386template <typename MeshType, typename DerivedT, typename DerivedB>
387void corner_tangent_bitangent_raw(
388 const MeshType& mesh,
389 Eigen::PlainObjectBase<DerivedT>& T,
390 Eigen::PlainObjectBase<DerivedB>& BT,
391 bool pad_with_sign)
392{
393 static_assert(MeshTrait<MeshType>::is_mesh(), "Input type is not Mesh");
394
395 la_runtime_assert(mesh.get_vertices().cols() == 3, "Mesh must be 3D");
396 la_runtime_assert(mesh.is_uv_initialized(), "UVs must be initialized");
397
398 using Scalar = typename MeshType::Scalar;
399 using Index = typename MeshType::Index;
400
401 const auto& F = mesh.get_facets();
402 const auto& V = mesh.get_vertices();
403 const auto& UV_val = mesh.get_uv();
404 const auto& UV_indices = mesh.get_uv_indices();
405
406 const auto n_faces = mesh.get_num_facets();
407 const auto per_facet = mesh.get_vertex_per_facet();
408
409 T.resize(n_faces * per_facet, pad_with_sign ? 4 : 3);
410 BT.resize(n_faces * per_facet, pad_with_sign ? 4 : 3);
411
412 // Triangular mesh
413 if (per_facet == 3) {
414 tbb::parallel_for(Index(0), Index(F.rows()), [&](Index i) {
415 const auto& indices = F.row(i);
416 Eigen::Matrix<Scalar, 3, 3> vertices;
417 vertices.row(0) = V.row(indices(0));
418 vertices.row(1) = V.row(indices(1));
419 vertices.row(2) = V.row(indices(2));
420
421 Eigen::Matrix<Scalar, 3, 2> uvs;
422 uvs.row(0) = UV_val.row(UV_indices(i, 0));
423 uvs.row(1) = UV_val.row(UV_indices(i, 1));
424 uvs.row(2) = UV_val.row(UV_indices(i, 2));
425
426 bool sign;
427 auto T_BT = triangle_tangent_bitangent(vertices, uvs, sign);
428
429 for (Index k = 0; k < per_facet; k++) {
430 T.row(i * per_facet + k).template head<3>() = T_BT.row(0);
431 BT.row(i * per_facet + k).template head<3>() = T_BT.row(1);
432 if (pad_with_sign) {
433 T.row(i * per_facet + k).w() = Scalar(sign ? 1 : -1);
434 BT.row(i * per_facet + k).w() = Scalar(sign ? 1 : -1);
435 }
436 }
437 });
438 }
439 // Quad or mixed mesh
440 else if (per_facet == 4) {
441 tbb::parallel_for(Index(0), Index(F.rows()), [&](Index i) {
442 const auto& indices = F.row(i);
443
444 bool sign;
445 Eigen::Matrix<Scalar, 2, 3> T_BT;
446
447 // Triangle in mixed mesh
448 if (indices(3) == invalid<typename MeshType::Index>()) {
449 Eigen::Matrix<Scalar, 3, 3> vertices;
450 vertices.row(0) = V.row(indices(0));
451 vertices.row(1) = V.row(indices(1));
452 vertices.row(2) = V.row(indices(2));
453
454 Eigen::Matrix<Scalar, 3, 2> uvs;
455 uvs.row(0) = UV_val.row(UV_indices(i, 0));
456 uvs.row(1) = UV_val.row(UV_indices(i, 1));
457 uvs.row(2) = UV_val.row(UV_indices(i, 2));
458
459 T_BT = triangle_tangent_bitangent(vertices, uvs, sign);
460 }
461 // Quad
462 else {
463 Eigen::Matrix<Scalar, 4, 3> vertices;
464 vertices.row(0) = V.row(indices(0));
465 vertices.row(1) = V.row(indices(1));
466 vertices.row(2) = V.row(indices(2));
467 vertices.row(3) = V.row(indices(3));
468
469 Eigen::Matrix<Scalar, 4, 2> uvs;
470 uvs.row(0) = UV_val.row(UV_indices(i, 0));
471 uvs.row(1) = UV_val.row(UV_indices(i, 1));
472 uvs.row(2) = UV_val.row(UV_indices(i, 2));
473 uvs.row(3) = UV_val.row(UV_indices(i, 3));
474
475 T_BT = quad_tangent_bitangent(vertices, uvs, sign);
476 }
477
478 for (Index k = 0; k < per_facet; k++) {
479 if (indices(k) == invalid<typename MeshType::Index>()) break;
480 T.row(i * per_facet + k).template head<3>() = T_BT.row(0);
481 BT.row(i * per_facet + k).template head<3>() = T_BT.row(1);
482 if (pad_with_sign) {
483 T.row(i * per_facet + k).w() = Scalar(sign ? 1 : -1);
484 BT.row(i * per_facet + k).w() = Scalar(sign ? 1 : -1);
485 }
486 }
487 });
488 }
489}
490
491} // namespace tangent_bitangent_detail
492
494
506template <typename MeshType>
507void compute_corner_tangent_bitangent(MeshType& mesh, bool pad_with_sign = false)
508{
509 static_assert(MeshTrait<MeshType>::is_mesh(), "Input type is not Mesh");
510
511 la_runtime_assert(mesh.get_vertices().cols() == 3, "Mesh must be 3D");
512 la_runtime_assert(mesh.is_uv_initialized(), "UVs must be initialized");
513
514 typename MeshType::AttributeArray T;
515 typename MeshType::AttributeArray BT;
516 tangent_bitangent_detail::corner_tangent_bitangent_raw(mesh, T, BT, pad_with_sign);
517
518 mesh.add_corner_attribute("tangent");
519 mesh.add_corner_attribute("bitangent");
520 mesh.import_corner_attribute("tangent", T);
521 mesh.import_corner_attribute("bitangent", BT);
522}
523
535template <typename MeshType>
536void compute_indexed_tangent_bitangent(MeshType& mesh, bool pad_with_sign = false)
537{
538 static_assert(MeshTrait<MeshType>::is_mesh(), "Input type is not Mesh");
539
540 la_runtime_assert(mesh.get_vertices().cols() == 3, "Mesh must be 3D");
541 la_runtime_assert(mesh.is_uv_initialized(), "UVs must be initialized");
542 la_runtime_assert(mesh.get_vertex_per_facet() == 3, "Input must be a triangle mesh");
543 la_runtime_assert(mesh.has_indexed_attribute("normal"), "Mesh must have indexed normals");
544
545 typename MeshType::AttributeArray T;
546 typename MeshType::AttributeArray BT;
547 logger().debug("Compute corner tangent info");
548 tangent_bitangent_detail::corner_tangent_bitangent_raw(mesh, T, BT, pad_with_sign);
549 logger().debug("Accumulate tangent info");
550 tangent_bitangent_detail::accumulate_tangent_bitangent(mesh, T, BT);
551}
552
553} // namespace legacy
554} // namespace lagrange
lagrange::Mesh
Definition: Mesh.h:48
lagrange::logger
LA_CORE_API spdlog::logger & logger()
Retrieves the current logger.
Definition: Logger.cpp:40
la_runtime_assert
#define la_runtime_assert(...)
Runtime assertion check.
Definition: assert.h:169
lagrange::swap
constexpr void swap(function_ref< R(Args...)> &lhs, function_ref< R(Args...)> &rhs) noexcept
Swaps the referred callables of lhs and rhs.
Definition: function_ref.h:114
lagrange
Main namespace for Lagrange.
Definition: AABBIGL.h:30
lagrange::sign
int sign(T val)
Get the sign of the value Returns either -1, 0, or 1.
Definition: utils.h:43
lagrange::corner_to_edge_mapping
Eigen::Index corner_to_edge_mapping(const Eigen::MatrixBase< DerivedF > &F, Eigen::PlainObjectBase< DerivedC > &C2E)
Computes a mapping from mesh corners (k*f+i) to unique edge ids.
Definition: corner_to_edge_mapping.impl.h:36
lagrange::chain_corners_around_edges
void chain_corners_around_edges(const Eigen::MatrixBase< DerivedF > &facets, const Eigen::MatrixBase< DerivedC > &corner_to_edge, Eigen::PlainObjectBase< DerivedE > &edge_to_corner, Eigen::PlainObjectBase< DerivedN > &next_corner_around_edge)
Chains facet corners around edges of a mesh.
Definition: chain_corners_around_edges.h:39
lagrange::chain_corners_around_vertices
void chain_corners_around_vertices(typename DerivedF::Scalar num_vertices, const Eigen::MatrixBase< DerivedF > &facets, Eigen::PlainObjectBase< DerivedE > &vertex_to_corner, Eigen::PlainObjectBase< DerivedN > &next_corner_around_vertex)
Chains facet corners around vertices of a mesh.
Definition: chain_corners_around_vertices.h:41