Ensure correctness of `polycube_distortion`

[sebmestrallet]

Context

The ultimate goal is to make good hex-mehses. One of the metrics for good hex-meshes is the per-cell distortion like the Scaled Jacobian. However, to generate a good hex-mesh of a given 3D model with an iterative algorithm, it is too costly to go as far as the hex-mesh in the feedback loop.

Polycube-based hex-meshing proposes to iterate over polycubes (orthogonal polyhedra), a good polycube being a proxy for a good hex-mesh. Some of the metrics for good polycubes are the angle and area distortions of the mapping.

The pipeline can be further simplified by iterating on polycube labelings (association of surface triangles to ±X, ±Y or ±Z) instead of handling a (volume) polycube, a good labeling being a proxy for a good polycube. A labeling must be valid (i.e. being the representation of a polycube boundary), and a good labeling is compact, has high geometric fidelity (but not always) and all-monotone boundaries.

Currently in the code base we can evaluate labelings and hex-meshes (with the Scaled Jacobian). Contrary to some papers on polycube-based hex-meshing, we cannot provide values on polycube distortion. Evocube¹ leverages distortion metrics, on one hand for results analysis (Tables 1 & 3 + supplemental material), on the other hand as score components (section 2.1, workability).

Starting point

Like Evocube¹, we have two input triangle meshes, with the same number of vertices & triangles, with the same vertices connectivity (= the same triangles) but of different vertices coordinates.

Evocube section 2.1:

The per-triangle distortion is measured using the singular values $\sigma_1$ and $\sigma_2$ of the Jacobian of the mapping from the initial triangle to its equivalent in [the polycube domain].

On this point, the authors refer to Tarini et al. 2004 ².

On the score side (section 2.1), the per-triangle workability $e_w$ is:

$$e_w = \sigma_1 + \sigma_2 + \frac{1}{\sigma_1 \sigma_2} + \frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} - 4~~~~(eq.1)$$

Questions

Where does this -4 come from? $e_w$ seems in 3 parts:

• $\sigma_1 + \sigma_2 + \frac{1}{\sigma_1 \sigma_2}$ looks like the area distortion ($\sigma_1 \sigma_2 + \frac{1}{\sigma_1 \sigma_2}$) we will see below. Is there a typo?
• $\frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1}$ is the angle distortion we will see below.
• lastly $-4$. The ideal value for each distortion being $1$, why not subtract 2 so that the cost is null when the distortion is null?

As for the results analysis (section 4), Evocube uses the area distortion $D_a$ (Table 1), referring to Tarini et al. 2004 ². In Table 3, the angle distortion also appears. Their calculation is not explained.

In the source code of Evocube, we can find out how the computation of these metrics has been implemented.

On the score side, the function include/evaluator.h > evaluate() as indeed 4 components:

• the (in)valididy score invalid_score
• the compactness score compact_score
• the (geometric) fidelity score fidelity_score
• the workability fast_poly_score, obtained via src/quick_label_ev.cpp > QuickLabelEv::evaluate(), which is the average distortion + a cost for inverted triangles

Let's take the average distortion calculation step-by-step:

1. A polycube is obtained with src/quick_label_ev.cpp QuickLabelEv::LDLTDeformation(), which looks like a transcoding of fastbndpolycube and/or polycube_withHexEx.
2. The normals N_def of the polycube triangles are computed
3. Per-triangle distortions are computed via src/distortion.cpp > computeDisto(), which:
   1. Computes the Jacobian jacobians[f_id]
   2. Computes singular values s1 and s2 ($\sigma_1$ et $\sigma_2$) using Eigen::JacobiSVD
   3. Computes the area distortion s1 * s2 + 1.0 / (s1 * s2) - 2.0 (this goes along with a typo in the article, see above)
   4. Computes the angle distortion s1 / s2 + s2 / s1 - 2.0
   5. Writes the sum of the two distortions as the per-triangle distorsion
4. Squaring of each per-triangle distorsion
5. Computation of the distorsion integral over the surface via src/distortion.cpp > integrateDistortion(), which :
   1. Caps distorsions with a threshold at $10^8$, giving the distop vector
   2. Computes (A.array() * distop.array()).sum() / A.sum(). A being the triangles area vector of the input mesh, we have the weighted sum of distortions, normalized by the total area:

$$\frac{ \displaystyle\sum_{t \in \delta T_\Omega}{A_t \left( \sigma_1 \sigma_2 + \frac{1}{\sigma_1 \sigma_2} + \frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} -4 \right)} }{ \displaystyle\sum_{t \in \delta T_\Omega}{A_t} }~~~~(eq.2)$$

6. (For the final value, taking into consideration of the number of inverted triangles, i.e. polycube triangles for which the dot product with the label's direction is negative)

Question

Why 2 is subtracted to the area and angle distortions? It's consistent with the -4 seen above, but not consistent with the equations we'll see below.

As for the results analysis side, Evocube has a measurement app (app/measurement.cpp) which reads a triangle mesh and a polycube, to write in the logs (JSON files) the distortion metrics value.

The first operations are:

1. The computation of per-triangle areas for the input mesh (A1) and the polycube (A2)
2. The computation of the total areas for the input mesh (A_m) and the polycube (A_d)
3. The polycube is scaled to have the same total surface as the input mesh: V2 = V2 * std::sqrt(A_m / A_d). I suppose there is a square root because we're going from a scalar to dimension 2 (area).
4. Areas are recomputed
5. Jacobians are computed
6. Singular values are computed using Eigen::JacobiSVD
7. The stretch (another distortion metric) is computed via src/distortion.cpp > computeStretch()
8. The area distortion is computed via src/distortion.cpp > computeAreaDisto()
9. The angle distorsion is computed via src/distortion.cpp > computeAngleDisto()
10. The isometric distortion (another distortion metric) is computed via src/distortion.cpp > computeIsometricDisto()

Let's take these distortion computations one by one.

// Spherical Parametrization and Remeshing
// Praun & Hoppe
double computeStretch(const Eigen::VectorXd& A, double A_m, double A_d,
                      std::vector<std::pair<double, double>> per_tri_singular_values){
    
    double sum = 0;
    for (int f_id = 0; f_id < per_tri_singular_values.size(); f_id++){
        double s1 = per_tri_singular_values[f_id].first;
        double s2 = per_tri_singular_values[f_id].second;
        sum += A(f_id) * (s1 * s1 + s2 * s2) / 2.0;
    }
    sum /= A.sum();

    return (A_m / A_d) * (1.0 / std::pow(sum, 2)); 
}

So we have:

$$stretch = \frac{ \displaystyle\sum_{t \in \delta T_\Omega}{A_t} }{ \displaystyle\sum_{t \in \delta T_\Omega}{{A'}_t} } \frac{1}{\displaystyle\sum_{t \in \delta T_\Omega}{\left(A_t \frac{\sigma_1^2+\sigma_2^2}{2}\right)^2}}~~~~(eq.3)$$

// PolyCube-Maps
// Tarini & Hormann & Cignoni & Montani
double computeAreaDisto(const Eigen::VectorXd& A, std::vector<std::pair<double, double>> per_tri_singular_values){
    double sum = 0;
    for (int f_id = 0; f_id < per_tri_singular_values.size(); f_id++){
        double s1 = per_tri_singular_values[f_id].first;
        double s2 = per_tri_singular_values[f_id].second;
        sum += A(f_id) * 0.5 * (s1 * s2 + 1.0 / (s1 * s2));
    }

    return sum / A.sum();
}

This gives:

$$area\ distortion = \frac{ \displaystyle\sum_{t \in \delta T_\Omega}{0.5 \mathcal{A}_t \left( \sigma_1 \sigma_2 + \frac{1}{\sigma_1 \sigma_2} \right) } }{ \displaystyle\sum_{t \in \delta T_\Omega}{\mathcal{A}_t} }~~~~(eq.4)$$

double computeAngleDisto(const Eigen::VectorXd& A, std::vector<std::pair<double, double>> per_tri_singular_values){
    double sum = 0;
    for (int f_id = 0; f_id < per_tri_singular_values.size(); f_id++){
        double s1 = per_tri_singular_values[f_id].first;
        double s2 = per_tri_singular_values[f_id].second;
        sum += A(f_id) * 0.5 * (s1 / s2 + s2 / s1);
    }

    return sum / A.sum();
}

This gives:

$$angle\ distortion = \frac{ \displaystyle\sum_{t \in \delta T_\Omega}{0.5 \mathcal{A}_t \left( \frac{\sigma_1}{\sigma_2} + \frac{\sigma_2}{\sigma_1} \right) } }{ \displaystyle\sum_{t \in \delta T_\Omega}{\mathcal{A}_t} }~~~~(eq.5)$$

Lastly:

// Computing Surface PolyCube-Maps by Constrained Voxelization
// Yang, Fu, Liu
double computeIsometricDisto(const Eigen::VectorXd& A, std::vector<std::pair<double, double>> per_tri_singular_values){
    double sum = 0;
    for (int f_id = 0; f_id < per_tri_singular_values.size(); f_id++){
        double s1 = per_tri_singular_values[f_id].first;
        double s2 = per_tri_singular_values[f_id].second;
        sum += std::max(s1, 1.0 / s2) * A(f_id);
    }
    return sum / A.sum();
}

This results in:

$$isometric\ distortion = \frac{ \displaystyle\sum_{t \in \delta T_\Omega}{ \max \left(\sigma_1,\frac{1}{\sigma_2} \right) \mathcal{A}_t } }{ \displaystyle\sum_{t \in \delta T_\Omega}{\mathcal{A}_t} }~~~~(eq.6)$$

Transcoding

In validity-first-polycube-labeling, I simply copied the Evocube code, cf src/geometry_distortion.cpp, but using the STL & Geogram, instead of libigl & Eigen.

Result: our executable polycube_distortion gives the same values as measurement from Evocube.

However, I don't find the same values as the Evocube output data. When Evocube was executed on the whole input dataset, measurement, via init_from_folder, created a logs.json file for each 3D model. Then supplemental_generator browsed all sub-folders and generated a PDF report. By executing our polycube_distortion on the Evocube output mesh, I don't obtain the same values as the ones in logs.json, the latter being different from the one in the PDF (!). Was the distortions calculation modified in the meantime? The code history seems to say no.

Dive into referred articles

Let's start with the stretch, In the code, Evocube refer to Praun and Hoppe 2003 ³. Indeed, section 3.3 "Review of planar-domain stretch metrics", we find a local definition, at a point $(s,t)$, of the root-mean-square stretch $L^2(s,t)$:

$$L^2(s,t) = \sqrt{\frac{1}{2} \left( \sigma_1^2 + \sigma_2^2 \right)}~~~~(eq.7)$$

The $L^2$ norm of the stretch is the $L^2(s,t)$ integral over the surface $M$, to obtain the stretch metric $L^2(M)$:

$$L^2(M) = \sqrt{ \frac{1}{\mathcal{A}_M} \iint_{(s,t) \in D}{ \left( L^2(s,t) \right)^2~~~d\mathcal{A}_M(s,t) } }\\ = \sqrt{ \frac{1}{\mathcal{A}_M} \iint_{(s,t) \in D}{ \frac{1}{2} \left( \sigma_1^2 + \sigma_2^2 \right)~~~d\mathcal{A}_M(s,t) } }~~~~(eq.8)$$

With $D$ a planar domain and $\mathcal{A}_M$ the total area (of the input mesh? of the polycube?)

In case of a triangle mesh, the parametrization $\phi$ is piecewise linear and its Jacobian $J_\phi$ is constant over each triangle. This way, the integral of the metric can be rewritten has a finite sum.

Question

How do you go from the integral to the sum?

$$\sqrt{ \frac{1}{\mathcal{A}_M} \displaystyle\sum_{t \in \delta T_\Omega}{ \frac{1}{2} \left( \sigma_1^2 + \sigma_2^2 \right) } }~~~~~~?$$

Question

Why Evocube (cf eq.3 of this issue) uses an inverse function? And why the sum is squared, whereas here there is an encompassing square root?

Praun and Hoppe 2003 ³ refers to Sander et al. 2001 ⁴ about the stretch. This article talks about it in section 3 "Texture stretch metric". We find the eq.7 of this issue (calculation of $L^2(T)$ ). However, the integral over a discrete surface $M={T_i}$ is made clear:

$$L^2(M) = \sqrt{\frac{ \displaystyle\sum_{T_i \in M}{ \left( L^2(T) \right)^2 \mathcal{A}'(T_i) } }{ \displaystyle\sum_{T_i \in M}{\mathcal{A}'(T_i)} }}\\\ = \sqrt{\frac{ \displaystyle\sum_{T_i \in M}{ \frac{\sigma_1^2 + \sigma_2^2}{2} \mathcal{A}'(T_i) } }{ \displaystyle\sum_{T_i \in M}{\mathcal{A}'(T_i)} }}~~~~(eq.9)$$

With $\mathcal{A}'(T_i)$ the area of triangle $T_i$ in 3D.

Question

That doesn't look like what Evocube does, does it?

The normalization, so that 1.0 is the minimum, is doable in two ways:

• by the scaling, upstream, of the texture domain, so that its area is the same as the area of the 3D surface
• by the multiplication of the stretch by the factor $\sqrt{\frac{ \displaystyle\sum_{T_i \in M}{A(T_i)} }{ \displaystyle\sum_{T_i \in M}{A'(T_i)} }}$

Question

Evocube does both, doesn't it?

After the stretch, let's see area and angle distortions. In the source code, Evocube referred to Tarini et al. 2004 ². In the Table 1 legend, we can read they are "measured by integrating and normalizing the values $\sigma_1 \sigma_2 + 1/(\sigma_1 \sigma_2)$ and $\sigma_1/\sigma_2 + \sigma_2/\sigma_1$ , respectively". For more details, the reader is referred to Degener et al. 2003 ⁵ and Floater & Hormann 2004 ⁶.

In Degener et al. 2003 ⁵, we have in section 2.5:

$$E_\text{area}(\omega) := \sqrt{\det I_\phi(\omega)} + \frac{1}{\sqrt{\det I_\phi(\omega)}}~~~~(eq.10)$$ $$E_\text{angle}(\omega) := \sqrt{\frac{\lambda_\text{max}}{\lambda_\text{min}}} + \sqrt{\frac{\lambda_\text{min}}{\lambda_\text{max}}} ~~~~(eq.11)$$

I suppose eq.10 is equivalent to $\sigma_1 \sigma_2 + \frac{1}{\sigma_1 \sigma_2}$?

The next section (2.6) talks about discretization and what are becoming $E_\text{area}(\omega)$ and $E_\text{angle}(\omega)$ in case of a triangle mesh, but I didn't understand the reasoning sequence.

In Floater & Hormann 2004 ⁶, we find (section 8):

$$E_A(G) = \det J + \frac{1}{\det J} = \sigma_1 \sigma_2 + \frac{1}{\sigma_1 \sigma_2} ~~~~(eq.12)$$

I didn't find a similar formulation of the angle distortion.

Lastly, the isometric distortion, for which the Evocube source code is referring Yang et al. 2019 ⁷. This article defines this distortion on a triangle $t$ as:

$$d_{t} = \{\sigma_{max}, 1/ \sigma_{min}\} ~~~~(eq.13)$$

with $\sigma_{max}$ and $\sigma_{min}$ the singular values of the Jacobian on $t$. When $\sigma_{max} = \sigma_{min} = 1$, $d_{t}$ reaches its minimum of $1$, indicating there is no distortion.

Question

What do the braces mean?

Yang et al. 2019 ⁷ refers (in section 4) to Fu et al. 2015 ⁸, which present the isometric distortion as penalizing both the conformal distortion (angle?) and the area distortion:

$$(\sigma_1 \sigma_2^{-1} + \sigma_2 \sigma_1^{-1})(\sigma_1 \sigma_2 + \sigma_1^{-1} \sigma_2^{-1})^\theta ~~~~(eq.14)$$

Question

What is $\theta$ here?

They refer to Degener et al. 2003 ⁵ as well, mentioned above for another distortion. In section 2.6 we can find:

$$E_T := E_{angle}(T) \cdot E_{area}(T)^\theta ~~~~(eq.15)$$

Which looks like eq.14?

But neither eq.13 nor eq.14 nor eq.15 look like the Evocube code.

Conclusion

• the workability of Evocube is, although not clearly transcribed in the article, the sum of the angle and area distortions, from which we subtract 4 (not explained).
• the stretch calculation in the Evocube code seems very different from referred articles
• the isometric distortion in the Evocube code does not look like the equations of referred articles
• the metrics normalization does not seem homogeneous
• I can't find the same values that Evocube recorded

Footnotes

1. Dumery, Protais, Mestrallet, Bourcier, Ledoux, "Evocube: a Genetic Labeling Framework for Polycube-Maps", Computer Graphics Forum, 2022, https://onlinelibrary.wiley.com/doi/10.1111/cgf.14649 ↩ ↩²

2. Tarini, Hormann, Cignoni, Montani, "PolyCube-Maps", Proceeding of SIGGRAPH, 2004, https://vcg.isti.cnr.it/polycubemaps/resources/sigg04.pdf ↩ ↩² ↩³

3. Praun, Hoppe, "Spherical Parametrization and Remeshing", ACM Transactions on Graphics, 2003, https://hhoppe.com/sphereparam.pdf ↩ ↩²

4. Sander, Snyder, Gortier, Hoppe, "Texture Mapping Progressive Meshes", SIGGRAPH, 2001, https://hhoppe.com/tmpm.pdf ↩

5. Degener, Meseth, Klein, "An adaptable surface parametrization method", International Meshing Roundtable, 2003, https://cg.cs.uni-bonn.de/backend/v1/files/publications/degener-2003-adaptable.pdf ↩ ↩² ↩³

6. Floater, Hormann, "Surface parametrization: a tutorial and survey", Advances in Multiresolution for Geometric Modelling. Mathematics and Visualization, 2004, https://graphics.stanford.edu/courses/cs468-05-fall/Papers/param-survey.pdf ↩ ↩²

7. Yang, Fu, Liu, "Computing Surface PolyCube-Maps by Constrained Voxelization", Computer Graphics Forum, 2019, https://onlinelibrary.wiley.com/doi/full/10.1111/cgf.13838 ↩ ↩²

8. Fu, Liu, Guo, "Computing locally injective mappings by advanced MIPS", ACM Transactions on Graphics, Volume 34, Issue 4, 2015, https://dl.acm.org/doi/10.1145/2766938 ↩