22
33#include " embedding.hh"
44
5- namespace CurvatureMetric
6- {
5+ namespace CurvatureMetric {
76
8- Scalar
9- squared_area (Scalar li, Scalar lj, Scalar lk)
10- {
7+ Scalar squared_area (Scalar li, Scalar lj, Scalar lk)
8+ {
119 // Sort the lengths for numerical stability
1210 Scalar a = li;
1311 Scalar b = lj;
1412 Scalar c = lk;
15- if (a < b)
16- swap (a, b);
17- if (a < c)
18- swap (a, c);
19- if (b < c)
20- swap (b, c);
13+ if (a < b) swap (a, b);
14+ if (a < c) swap (a, c);
15+ if (b < c) swap (b, c);
2116
2217 // Compute the area
2318 Scalar A = a + (b + c);
@@ -26,33 +21,29 @@ namespace CurvatureMetric
2621 Scalar D = a + (b - c);
2722
2823 return A * B * C * D / 16.0 ;
29- }
24+ }
3025
31- Scalar
32- squared_area_length_derivative (Scalar variable_length, Scalar lj, Scalar lk)
33- {
26+ Scalar squared_area_length_derivative (Scalar variable_length, Scalar lj, Scalar lk)
27+ {
3428 // Sort the lengths and keep track of the derivative edge.
3529 Scalar a = variable_length;
3630 Scalar b = lj;
3731 Scalar c = lk;
3832 bool deriv_a = true ;
3933 bool deriv_b = false ;
40- if (a < b)
41- {
42- swap (a, b);
43- deriv_a = false ;
44- deriv_b = true ;
34+ if (a < b) {
35+ swap (a, b);
36+ deriv_a = false ;
37+ deriv_b = true ;
4538 }
46- if (a < c)
47- {
48- swap (a, c);
49- deriv_a = false ;
39+ if (a < c) {
40+ swap (a, c);
41+ deriv_a = false ;
5042 }
51- if (b < c)
52- {
53- swap (b, c);
54- // The derivative is b iff the derivate was c before the swap
55- deriv_b = !(deriv_a || deriv_b);
43+ if (b < c) {
44+ swap (b, c);
45+ // The derivative is b iff the derivate was c before the swap
46+ deriv_b = !(deriv_a || deriv_b);
5647 }
5748
5849 // Compute stable factors for Heron's formula
@@ -68,51 +59,45 @@ namespace CurvatureMetric
6859 Scalar TmC = (S * A * B) / 16.0 ;
6960
7061 // Compute the derivative for li
71- if (deriv_a)
72- {
73- return TmS - TmA + TmB + TmC;
62+ if (deriv_a) {
63+ return TmS - TmA + TmB + TmC;
64+ } else if (deriv_b) {
65+ return TmS + TmA - TmB + TmC;
66+ } else {
67+ return TmS + TmA + TmB - TmC;
7468 }
75- else if (deriv_b)
76- {
77- return TmS + TmA - TmB + TmC;
78- }
79- else
80- {
81- return TmS + TmA + TmB - TmC;
82- }
83- }
69+ }
8470
85- VectorX squared_areas (const DifferentiableConeMetric & cone_metric)
86- {
71+ VectorX squared_areas (const DifferentiableConeMetric& cone_metric)
72+ {
8773 int num_halfedges = cone_metric.n_halfedges ();
8874 VectorX he2areasq (num_halfedges);
8975
9076 int num_faces = cone_metric.h .size ();
9177 // #pragma omp parallel for
92- for (int f = 0 ; f < num_faces; f++)
93- {
94- // Get halfedges of face f
95- int hi = cone_metric.h [f];
96- int hj = cone_metric.n [hi];
97- int hk = cone_metric.n [hj];
98-
99- // Get lengths of the halfedges
100- Scalar li = cone_metric.l [hi];
101- Scalar lj = cone_metric.l [hj];
102- Scalar lk = cone_metric.l [hk];
103-
104- // Compute the area of the face adjacent to the halfedges
105- Scalar areasq = squared_area (li, lj, lk);
106- he2areasq[hi] = areasq;
107- he2areasq[hj] = areasq;
108- he2areasq[hk] = areasq;
78+ for (int f = 0 ; f < num_faces; f++) {
79+ // Get halfedges of face f
80+ int hi = cone_metric.h [f];
81+ int hj = cone_metric.n [hi];
82+ int hk = cone_metric.n [hj];
83+
84+ // Get lengths of the halfedges
85+ Scalar li = cone_metric.l [hi];
86+ Scalar lj = cone_metric.l [hj];
87+ Scalar lk = cone_metric.l [hk];
88+
89+ // Compute the area of the face adjacent to the halfedges
90+ Scalar areasq = squared_area (li, lj, lk);
91+ he2areasq[hi] = areasq;
92+ he2areasq[hj] = areasq;
93+ he2areasq[hk] = areasq;
10994 }
11095
11196 return he2areasq;
112- }
97+ }
11398
114- VectorX areas (const DifferentiableConeMetric & cone_metric)
115- {
99+ VectorX areas (const DifferentiableConeMetric& cone_metric)
100+ {
116101 int num_halfedges = cone_metric.n_halfedges ();
117102
118103 // Compute squared areas
@@ -121,46 +106,44 @@ namespace CurvatureMetric
121106
122107 // Take square roots
123108 VectorX he2area (num_halfedges);
124- for (int h = 0 ; h < num_halfedges; ++h)
125- {
126- he2area[h] = sqrt (std::max<Scalar>(he2areasq[h], 0.0 ));
109+ for (int h = 0 ; h < num_halfedges; ++h) {
110+ he2area[h] = sqrt (std::max<Scalar>(he2areasq[h], 0.0 ));
127111 }
128112
129113 return he2area;
130- }
114+ }
131115
132- VectorX squared_area_length_derivatives (const DifferentiableConeMetric & cone_metric)
133- {
116+ VectorX squared_area_length_derivatives (const DifferentiableConeMetric& cone_metric)
117+ {
134118 int num_halfedges = cone_metric.n_halfedges ();
135119 VectorX he2areasqderiv (num_halfedges);
136120
137121 // Compute maps from halfedges to derivatives of area with respect to the edge
138122 // length
139123 int num_faces = cone_metric.h .size ();
140124 // #pragma omp parallel for
141- for (int f = 0 ; f < num_faces; f++)
142- {
143- // Get halfedges of face f
144- int hi = cone_metric.h [f];
145- int hj = cone_metric.n [hi];
146- int hk = cone_metric.n [hj];
147-
148- // Get lengths of the halfedges
149- Scalar li = cone_metric.l [hi];
150- Scalar lj = cone_metric.l [hj];
151- Scalar lk = cone_metric.l [hk];
152-
153- // Compute the derivative of the area of f with respect to each halfedge
154- he2areasqderiv[hi] = squared_area_length_derivative (li, lj, lk);
155- he2areasqderiv[hj] = squared_area_length_derivative (lj, lk, li);
156- he2areasqderiv[hk] = squared_area_length_derivative (lk, li, lj);
125+ for (int f = 0 ; f < num_faces; f++) {
126+ // Get halfedges of face f
127+ int hi = cone_metric.h [f];
128+ int hj = cone_metric.n [hi];
129+ int hk = cone_metric.n [hj];
130+
131+ // Get lengths of the halfedges
132+ Scalar li = cone_metric.l [hi];
133+ Scalar lj = cone_metric.l [hj];
134+ Scalar lk = cone_metric.l [hk];
135+
136+ // Compute the derivative of the area of f with respect to each halfedge
137+ he2areasqderiv[hi] = squared_area_length_derivative (li, lj, lk);
138+ he2areasqderiv[hj] = squared_area_length_derivative (lj, lk, li);
139+ he2areasqderiv[hk] = squared_area_length_derivative (lk, li, lj);
157140 }
158141
159142 return he2areasqderiv;
160- }
143+ }
161144
162- VectorX squared_area_log_length_derivatives (const DifferentiableConeMetric & cone_metric)
163- {
145+ VectorX squared_area_log_length_derivatives (const DifferentiableConeMetric& cone_metric)
146+ {
164147 int num_halfedges = cone_metric.n_halfedges ();
165148
166149 // Compute squared areas length derivatives
@@ -169,12 +152,11 @@ namespace CurvatureMetric
169152
170153 // Apply chain rule to A(l) = A(e^(lambda/2))
171154 VectorX he2areasq_log_deriv (num_halfedges);
172- for (int h = 0 ; h < num_halfedges; ++h)
173- {
174- he2areasq_log_deriv[h] = he2areasq_deriv[h] * cone_metric.l [h] / 2.0 ;
155+ for (int h = 0 ; h < num_halfedges; ++h) {
156+ he2areasq_log_deriv[h] = he2areasq_deriv[h] * cone_metric.l [h] / 2.0 ;
175157 }
176158
177159 return he2areasq_log_deriv;
178- }
179-
180160}
161+
162+ } // namespace CurvatureMetric
0 commit comments