Geometry Processing on Voxel Objects

David Coeurjolly
CNRS, Université de Lyon

           

\( %%%% If you see this, MathJax wasn't loaded ... \definecolor{myred}{RGB}{231,76,60} \definecolor{mydarkred}{RGB}{166,51,38} \definecolor{mygreen}{RGB}{38,166,115} \definecolor{myblue}{RGB}{23,110,179} \definecolor{mylightblue}{RGB}{104,171,216} \definecolor{myyellow}{RGB}{101,103,10} \definecolor{mygrey}{RGB}{101,123,131} \definecolor{mynormal}{RGB}{41,128,185} \)

Digital geometry model in one slide

Topology and geometry processing of objects defined on regular lattices

  • Digital Objects = subsets of $\mathbb{Z}^d$
  • Digital Topology = adjacency relationship induced by the lattice
  • Digital surfaces = cells of a Cartesian cubical complex

  • Geometrical predicates = integer only computations
  • $\Rightarrow$ strong arithmetical results when doing geometry on grids

Why voxel objects?




  • Widely used in geometric modeling and rendering to represent complex and interactive scenes
  • Nice mathematical modeling framework

$64K^3$ voxel grid!
Villanueva, Alberto Jaspe, Fabio Marton, and Enrico Gobbetti. "SSVDAGs: symmetry-aware sparse voxel DAGs." In Proceedings of the 20th ACM SIGGRAPH Symposium on Interactive 3D Graphics and Games, pp. 7-14. ACM, 2016.

Material Sciences Applications

Non-invasive snow micro-tomographic images
X-ray Micro-tomograph
3SR Lab and CEN/CNRM - GAME URA 1357/Météo-France - CNRS
Up to $2048^3$

Outline




  • Curvature tensor estimation

  • Laplace-Beltrami operator on digital surfaces

  • Variational geometry processing: Ambrosio-Tortorelli functional

Collaborators: Jacques-Olivier Lachaud (Chambéry), Tristan Roussillon (Lyon), Nicolas Bonneel (Lyon), Jérémy Levallois (Lyon), Thomas Caissard (Lyon), Marion Foare (Lyon)


C., Jacques-Olivier Lachaud, Jérémy Levallois. "Multigrid Convergent Principal Curvature Estimators in Digital Geometry". Computer Vision and Image Understanding, 129(1):27-41, June 2014.
Thomas Caissard, C., Jacques-Olivier Lachaud, Tristan Roussillon. "Laplace–Beltrami Operator on Digital Surfaces". Journal of Mathematical Imaging and Vision, 2018.
C., Marion Foare, Pierre Gueth, Jacques-Olivier Lachaud. "Piecewise smooth reconstruction of normal vector field on digital data". Computer Graphics Forum (Proceedings of Pacific Graphics), 35(7), September 2016.
Nicolas Bonneel, C., Pierre Gueth, Jacques-Olivier Lachaud. "Mumford-Shah Mesh Processing using the Ambrosio-Tortorelli Functional". Computer Graphics Forum (Proceedings of Pacific Graphics), 37(7), October 2018.

Curvature Tensor Estimation

Digitization model

Given $\Shape \subset \R^d$, its digitization at gridstep $h$ is $$\DigF{\Shape}{h} = \left(\frac{1}{h} \cdot \Shape\right) \cap \Z^d$$

$\Shape \in \Shapes$
$\DigF{\Shape}{1}$
$\DigF{\Shape}{0.5}$
$\DigF{\Shape}{0.25}$



Example: $h^2|\DigF{\Shape}{h}|$ converges to the measure of $\Shape$ as $h\rightarrow 0$

[Gauss, Dirichlet, Huxley...]

Multigrid convergence of a local estimator


Given a digitization process $\AnyDig$, a local discrete geometric estimator $\hat{E}$ of some geometric quantity $E$ is multigrid convergent for the family of shapes $\Shapes$ if and only if, for any $\Shape \in \Shapes$, there exists a grid step $h_\Shape>0$ such that the estimate $\hat{E}(\AnyDig_{\Shape}(h),\hat{\vx},h)$ is defined for all $\hat{\vx} \in \Bd{\Body{\AnyDig_{\Shape}(h)}{h}}$ with $0< h< h_\Shape$, and for any $\vx \in \dS$,

$$ \forall \hat{\vx} \in \Bd{\Body{\AnyDig_{\Shape}(h)}{h}} \text{ with } \| \hat{\vx} -\vx\|_\infty \le h, \quad\quad \color{myred}{| \hat{E}(\AnyDig_{\Shape}(h),\hat{\vx},h) - E(\Shape,\vx) | \le \tau_{\Shape,\vx}(h)}, $$
where $\tau_{\Shape,\vx}: \R^{+} \setminus\{0\} \rightarrow \R^+$ has null limit at $0$. The convergence is uniform for $\Shape$ when every $\tau_{\Shape,\vx}$ is bounded from above by a function $\tau_\Shape$ independent of $\vx \in \Bd{\Shape}$ with null limit at $0$.

Digital/Continuous mapping

    Let $\Shape$ be a compact domain of $\R^d$ such that $\partial X$ has positive reach greater that $R$. Let $\partial_h\Shape\EqDef\partial\Body{\Dig_h(\Shape))}{h}$. Then for any $0 < h < 2R/\sqrt{2}$,

  • The Hausdorff between $\partial \Shape$ and $\partial_h\Shape$ is bounded by $\sqrt{d}h/2$
  • For $d=2$, there exists an homeomorphism between $\partial X$ and $\partial_h X$
  • For $d\geq 3$, no homeomorphism ☹, but
    • Projection operator $\xi:\,\partial_h\Shape\rightarrow\partial \Shape$ is surjective
    • Area of non-injective parts of $\xi$ the tends to zero

Digitization as an Hausdorff sampling of the continuous object

Can we estimate the curvature tensor on digital surfaces with multigrid convergence properties ?

Huge literature on differential quantity estimators




→ accuracy depends on the mesh/point cloud quality
→ incompatible constraints on convergence theorem w.r.t. digital surface

Main contributions: Integral Invariant approach


If $\Shape$ is compact and $\partial \Shape$ has positive reach $\rho$ and $C^3$-continuity then

  • Mean curvature and principal curvatures estimations converge in $O(h^{\frac{1}{3}})$
  • Normal vector field estimation converges in $O(h^{\frac{2}{3}})$
  • Principal Curvature directions converge in $\frac{1}{|\PrincCurv{1}(\Shape,\vx)-\PrincCurv{2}(\Shape,\vx)|} O(h^{\frac{1}{3}})$

Overall proof scheme

$\vx$ $\cdot$

SVG-Viewer needed.

${\kappa}(\Shape,\vx) \EqDef \underbrace{\frac{3\pi}{2{{R}}} - \frac{3\cdot { A_R(\Shape,\vx)}}{{{R}}^3}}_{ \tilde{\kappa}^{R}(\Shape,\vx) } + O(R)$ [Pottmann et al. 2007]

SVG-Viewer needed.

SVG-Viewer needed.

$A_R(\Shape,\vx) \rightarrow \AreaC(\Ball{R/h}{\vx/h} \cap \DigF{\Shape}{h})$

SVG-Viewer needed.

+ [Pottmann et al. 2007] $\quad\CurvH{R}(\DSh,\vx,h)$

SVG-Viewer needed.

SVG-Viewer needed.

$\CurvH{R}(\DSh,\hat{\vx},h) \rightarrow \Curv(\Shape,\vx)$

Area estimator "à la" Gauss

\[ {\color{myblue}\AreaC(\DSh)} \EqDef h^2 \Card (\DSh) \]
$\quad \AreaC(Y) = \Area(\Shape) + {\color{myblue}O(h^\beta)}\quad\quad\text{with}\quad \beta \ge 1$

SVG-Viewer needed.

  • $\beta = 1$ for $\Shapes^C$ [Gauss, Dirichlet]
  • $\beta = \frac{15}{11}-\epsilon$ for $\Shapes^{3-SC}$ (non null curvature) [Huxley 1990]

Area estimator "à la" Gauss

\[ {\color{myblue}\AreaC(Y)} \EqDef h^2 \Card (Y) \quad\quad\text{with}\quad Y = \Ball{R/h}{\vx/h} \cap \DigF{\Shape}{h} \]
$\quad \AreaC(Y) = A_R(\Shape,\vx) + {\color{myblue}O(h^\beta R^{2-\beta})} \quad\quad\text{with}\quad \beta \ge 1$
$\vx$
$\vx$
\( \AreaC(Y) = 12 \times 1^2 = {\color{myblue}12} \)
\( \AreaC(Y) = 50 \times (0.5)^2 = {\color{myblue}12.5} \)
\( \color{myred} A_R(\Shape,\vx) \approx 12.9 \)

Overall proof scheme - 1/3

SVG-Viewer needed.

SVG-Viewer needed.

$A_R(\Shape,\vx) \rightarrow \AreaC(\Ball{R/h}{\vx/h} \cap \DigF{\Shape}{h})$

SVG-Viewer needed.

+ [Pottmann et al. 2007] $\quad\CurvH{R}(\DSh,\vx,h)$

SVG-Viewer needed.

SVG-Viewer needed.

$\CurvH{R}(\DSh,\hat{\vx},h) \rightarrow \Curv(\Shape,\vx)$

Curvature estimation by digital integration

\[ \hat{\kappa}^{R}(\DigF{\Shape}{h},\vx,h) \EqDef \frac{3\pi}{2{{R}}} - \frac{3 {\color{myblue}\AreaC(Y)}}{{{R}}^3} {\color{mydarkred}+O(R)} \quad\text{with}\quad Y = \Ball{R/h}{\vx/h} \cap \DigF{\Shape}{h} \]
$\vx$
$\vx$
\( \hat{\kappa}^{3}(\DigF{\Shape}{1},\vx,1) = \frac{3\pi}{2\times 3} - \frac{3 \times {\color{myblue}12}}{3^3} \approx 0.237 \)
\( \hat{\kappa}^{3}(\DigF{\Shape}{0.5},\vx,0.5) = \frac{3\pi}{2 \times 3} - \frac{3 \times {\color{myblue}12.5}}{3^3} \approx 0.182 \)
\( \color{myred} {\kappa}(\Shape,\vx) = \frac{1}{7.5} \approx 0.133 \)

Overall proof scheme - 2/3

SVG-Viewer needed.

SVG-Viewer needed.

$A_R(\Shape,\vx) \rightarrow \AreaC(\Ball{R/h}{\vx/h} \cap \DigF{\Shape}{h})$

SVG-Viewer needed.

+ [Pottmann et al. 2007] $\quad\CurvH{R}(\DSh,\vx,h)$

SVG-Viewer needed.

SVG-Viewer needed.

$\CurvH{R}(\DSh,\hat{\vx},h) \rightarrow \Curv(\Shape,\vx)$

Positioning Error


$\pi^\Shape_h(\hat{\vx})$
For any 2D shape with positive reach, given $0 < h \le \reach{\Shape}$, let $\mathbf{n}(\Shape, \vx, l)$ be the segment centred at $\vx$, aligned with the normal vector at $\vx$ of length $2l$. Then : \[\begin{align} \pi^\Shape_h : \partial\Body{\DSh}{h} &\rightarrow \dS \,,\\ \hat{\vx} &\mapsto \vx = \pi^\Shape_h(\hat{\vx}) \,, \end{align} \] where $\vx$ is the unique point such that $\hat{\vx} \in \mathbf{n}\left(\Shape,\vx, {\color{myred}\frac{\sqrt{2}}{2}h}\right)$.
$\vx$
$\hat{\vx}$
${\color{mygreen}|| \vx - \hat{\vx} ||_\infty = O(h^{\alpha'})} \quad\quad\text{with}\quad\alpha' \ge 1$

Multigrid convergence of the digital curvature estimator





Let $\Shape$ be a convex shape in $\R^2$ with $C^3$ bounded positive curvature boundary. \[ \begin{split} \forall \vx \in \partial\Shape, \forall \hat{\vx} \in \partial\Body{\DSh}{h}, \| \hat{x} -x\|_\infty \le h \Implies\\ | \CurvH{R}(\DSh,\hat{\vx},h) - \Curv(\Shape,\vx) | =\quad & {\color{mydarkred} O(R)}\\ +\quad & {\color{myblue} O\left(\frac{h^\beta}{R^{1+\beta}}\right)}\\ +\quad & {\color{mygreen} O\left(\frac{h^{\alpha'}}{R^2}\right) + O\left(h^{\alpha'}\right) + O\left(\frac{h^{2\alpha'}}{R^2}\right)} \end{split} \]

→ Setting $R=kh^\alpha$, we select $\alpha$ to minimize all errors.

\[ \color{myred} \left| \CurvH{R}(\DSh,\hat{\vx},h) - \Curv(\Shape,\vx) \right| \le O\left(h^\frac{1}{3}\right) \quad\text{setting}\quad R=kh^\frac{1}{3} \]

Curvature tensor on digital surfaces

\[ \begin{align} &\PrincCurvH{1}{R}(\DSh,\vx,h) \EqDef \frac{6}{\pi R^6}\left({\color{myblue}\hat{\lambda}_2} - 3{\color{myblue}\hat{\lambda}_1}\right) + \frac{8}{5R}\\ &\PrincCurvH{2}{R}(\DSh,\vx,h) \EqDef \frac{6}{\pi R^6}\left({\color{myblue}\hat{\lambda}_1} - 3{\color{myblue}\hat{\lambda}_2}\right) + \frac{8}{5R}\\ &\PrincDirH{1}{R}(\DSh,\vx,h) \EqDef {\color{myblue}\hat{\nu}_1}\\ &\PrincDirH{2}{R}(\DSh,\vx,h) \EqDef {\color{myblue}\hat{\nu}_2}\\ &\NormalDirH{R}(\DSh,\vx,h) \EqDef {\color{myblue}\hat{\nu}_3} \end{align} \]

$\{{\color{myblue}\hat\nu_i,\hat\lambda_i}\}$ are the eigenvalues/eigenvectors of the covariance matrix of $B_r(\vx)\cap \Shape$

\[ \color{myred} \left| \PrincCurvH{i}{R}(\DSh,\hat{\vx},h) - \PrincCurv{i}(\Shape,\vx) \right| \le O\left(h^\frac{1}{3}\right) \quad\text{setting}\quad R=kh^\frac{1}{3} \]

\[ \begin{align} &\exists h_\Shape \in \R^+,\,\, \forall h \in \R,\,\, 0 < h < h_\Shape,\,\, \,\,\nonumber\quad \quad\quad \quad\\ &\forall \vx \in \Bd{\Shape}, \quad \forall \hat{\vx} \in \Bd{\Body{\DSh}{h}} \text{ avec } \| \hat{\vx} -\vx\|_\infty \le h, \nonumber\quad \quad\quad \quad \\ &\quad \quad\quad \quad\| \PrincDirH{1}{R}( \DSh, \hat{\vx}, h ) - \PrincDir{1}(\Shape, \vx) \| \le \frac{1}{|\PrincCurv{1}(\Shape,\vx)-\PrincCurv{2}(\Shape,\vx)|} O(h^{\frac{1}{3}})\,, \\ &\quad \quad\quad \quad\| \PrincDirH{2}{R}( \DSh, \hat{\vx}, h ) - \PrincDir{2}(\Shape, \vx) \| \le \frac{1}{|\PrincCurv{1}(\Shape,\vx)-\PrincCurv{2}(\Shape,\vx)|} O(h^{\frac{1}{3}})\,, \\ &\quad \quad\quad \quad\| \NormalDirH{R}( \DSh, \hat{\vx}, h ) - \NormalDir(\Shape, \vx) \| \le O(h^{\frac{2}{3}})\,. \end{align} \]

In summary





  • Robust, Efficient implementation (convolutions)
  • Parametrized by a integration radius $R$ or a grid step $h$
  • Proof relies on digital/continuous relationships and geometrical moment estimation

Laplace-Beltrami Operator

Motivations

\[ \Delta u = \nabla\cdot\nabla u \]

Key operator for many geometry processing tasks


Discretization of the Laplace-Beltrami operator

Many discretization scheme for triangular meshes, polygonal meshes, point clouds...
SYM LOC LIN POS PSD $C^2$-CON
Mean Value x $\checkmark$ $\checkmark$ $\checkmark$ x x
Intrinsic Del $\checkmark$ x $\checkmark$ $\checkmark$ $\checkmark$ x
Combinatorial $\checkmark$ $\checkmark$ x $\checkmark$ $\checkmark$ x
Cotan x $\checkmark$ $\checkmark$ x $\checkmark$ x
Polygonal Lap. x $\checkmark$ $\checkmark$ x $\checkmark$ x
Convolutional x x ? $\checkmark$ ? $\checkmark$
$r-$local $\checkmark$ x ? $\checkmark$ ? $\checkmark$

(update of "Discrete Laplace operators: No free lunch" [Wardetzky et al., 2007])

Strong consistency of the operator: $\quad\text{lim}_{\varepsilon \rightarrow 0} || \Delta_\varepsilon v - \Delta v ||_{L^\infty} = \text{lim}_{\varepsilon \rightarrow 0} \,\, \text{sup}_{x \in \partial M} | (\Delta_\varepsilon v) (x) - (\Delta v) (x) | = 0,\quad \forall v\in C^2(\partial M).$

What about Digital Surfaces ?



  • No Laplace-Beltrami operator with strong consistency property exists on digital surfaces
  • Anisotropic nature of digital surfaces may lead to geometrical inconsistencies

Heat equation based Laplace-Beltrami operator on meshes

$\Delta g(x,t) = \frac{\partial }{\partial t} g(x,t), u=g(•,0)\quad\rightarrow \quad $ $\Delta g(x,t) = \text{lim}_{t\rightarrow 0} \frac{1}{t}\int_{\partial M} p(t,x,y)(u(y)-u(x)) dy$

$$(\mathcal{L}_t\, u)(x) \EqDef \frac{1}{t (4 \pi t)^{\frac{d}{2}}} \int_{y \in \partial M} e^{- \frac{ ||y - x||^2 }{4 t}} (u(y) - u(x))dy,$$

As $t\rightarrow 0$, $(\mathcal{L}_t\, u)$ converges to $\Delta\,u$.

$$(\mathcal{L}_{MESH}\, u)(p) \EqDef \frac{1}{4 \pi t^2} \sum_{f\in F} \frac{A_f}{3} \sum_{q \in V(f)} e^{- \frac{||p-q||^2}{4t}}\left(u(q) - u(p)\right)$$

If the mesh is a nice triangulation of a smooth manifold, $(\mathcal{L}_{MESH}\, u)$ converges to $(\mathcal{L}_t\, u)$
($t\approx$ 1/density).

(multigrid) Digital Surfaces are not nice triangulations!

Convolution based Laplace-Beltrami operator on Digital Surfaces

$$(L_h\, \tilde u)({\mathbf{s}}) \EqDef \frac{1}{t_h(4 \pi t_h)^{\frac{d}{2}}} {\color{mygreen}\sum_{\mathbf{r}\in S}} e^{- \frac{||{\mathbf{r}} - {\mathbf{s}}||^2}{4t_h}} [\tilde u({{\mathbf{r}}}) - \tilde u({{\mathbf{s}}})]{\color{myred} \mu(\mathbf{r})}$$

As $h\rightarrow 0$, $(L_h\, \tilde{u})$ converges to $\Delta\,u$.


($\tilde{u}$ is the extension of $u$ from $M$ to $\mathbb{R}^3$ along the normal vectors)

Measure of a surface element $\quad\mu(s)\EqDef \vn_s \cdot \vn_s^e$



Convolution based Laplace-Beltrami operator on Digital Surfaces

$$(L_h\, \tilde u)({\mathbf{s}}) \EqDef \frac{1}{t_h(4 \pi t_h)^{\frac{d}{2}}} {\color{mygreen}\sum_{\mathbf{r}\in S}} e^{- \frac{||{\mathbf{r}} - {\mathbf{s}}||^2}{4t_h}} [\tilde u({{\mathbf{r}}}) - \tilde u({{\mathbf{s}}})]{\color{myred} \mu(\mathbf{r})}$$

As $h\rightarrow 0$, $(L_h\, \tilde{u})$ converges to $\Delta\,u$.


($\tilde{u}$ is the extension of $u$ from $M$ to $\mathbb{R}^3$ along the normal vectors)

Sketch of the proof



$$ \Big| \, (\Delta u) (\xi({s})) - (L_h \tilde{u})({s})\,\Big| \leq {\Big| \, (\Delta u)(\xi({s})) - (\mathcal{L}_t u)(\xi({s})) \Big|} + {\Big|\, (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L} \tilde{u})({s}) \Big|} + {\Big|\, (\mathcal{L}_t \tilde{u})({s}) - (L_h \tilde{u})({s}) \Big| } $$

Sketch of the proof



$$ \Big| \, (\Delta u) (\xi({s})) - (L_h \tilde{u})({s})\,\Big| \leq \underbrace{\Big| \, (\Delta u)(\xi({s})) - (\mathcal{L}_t u)(\xi({s})) \Big|}_{\text{[Belkin et al]}} + \underbrace{\Big|\, (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L} \tilde{u})({s}) \Big|}_{\text{Projection error}} + { \underbrace{\Big|\, (\mathcal{L}_t \tilde{u})({s}) - (L_h \tilde{u})({s}) \Big| }_{\text{Digital integration error}}} $$

Sketch of the proof



$$ \Big| \, (\Delta u) (\xi({s})) - (L_h \tilde{u})({s})\,\Big| \leq \underbrace{\Big| \, (\Delta u)(\xi({s})) - (\mathcal{L}_t u)(\xi({s})) \Big|}_{\text{[Belkin et al]}} + {\color{mygreen} \underbrace{\Big|\, (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L}_t \tilde{u})({s}) \Big|}_{\text{Projection error}} } + { \underbrace{\Big|\, (\mathcal{L}_t \tilde{u})({s}) - (L_h \tilde{u})({s}) \Big| }_{\text{Digital integration error}}} $$

Let $\mathbf{s} \in \partial_h M$, a function $u\in C^2(\partial M)$ and its extension $\tilde{u}$. For $t_h = h^{\alpha}$, $0 < \alpha \leq \frac{2}{2 + d}$ and $h \leq h_{max}$ with $h_{max}$ the minimum between $\text{Diam}(\partial M)$, $K_3(d, \alpha, \text{Diam}(\partial M))$ and $R / \sqrt{d+1}$, we have \begin{aligned}\\ | (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L}_t \tilde{u})({s}) | \leq \text{Area}(\partial M) \,||\nabla u||_\infty \left[ K_1(d) h^{1 - \alpha(1 + \frac{d}{2})} + K_2(d) h^{2 - \alpha \frac{3 + d}{2}} \right] \end{aligned} with \[ K_1(d) \EqDef \frac{\sqrt{d+1}}{ 2^{d - 1} e \pi^{\frac{d}{2}}} \text{ and } K_2(d) \EqDef \frac{3 (d+1)}{ 2^{d + \frac{5}{2}} \sqrt{e} \pi^{\frac{d}{2}}} . \]

Technical proof using the regularity of $u$ and Hausdorff distance between $\partial M$ and $\partial_h M$.

Sketch of the proof



$$ \Big| \, (\Delta u) (\xi({s})) - (L_h \tilde{u})({s})\,\Big| \leq \underbrace{\Big| \, (\Delta u)(\xi({s})) - (\mathcal{L}_t u)(\xi({s})) \Big|}_{\text{[Belkin et al]}} + \underbrace{\Big|\, (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L} \tilde{u})({s}) \Big|}_{\text{Projection error}} + {\color{myred} \underbrace{\Big|\, (\mathcal{L}_t \tilde{u})({s}) - (L_h \tilde{u})({s}) \Big| }_{\text{Digital integration error}}} $$

Sketch of the proof



$$ \Big| \, (\Delta u) (\xi({s})) - (L_h \tilde{u})({s})\,\Big| \leq \underbrace{\Big| \, (\Delta u)(\xi({s})) - (\mathcal{L}_t u)(\xi({s})) \Big|}_{\text{[Belkin et al]}} + \underbrace{\Big|\, (\mathcal{L}_t u)(\xi({s})) - (\mathcal{L} \tilde{u})({s}) \Big|}_{\text{Projection error}} + {\color{myred} \underbrace{\Big|\, (\mathcal{L}_t \tilde{u})({s}) - (L_h \tilde{u})({s}) \Big| }_{\text{Digital integration error}}} $$

Let $ M $ be a compact domain whose boundary has positive reach $R$. For $h \leq \frac{R}{\sqrt{d+1}}$, the digital integral is multigrid convergent toward the integral over $\partial M $. More precisely, for any measurable function $f : \mathbb{R}^{d+1} \rightarrow \mathbb{R}$, one gets
$$ \,\\ \Bigg| \int_{\partial M } f(x)dx - \text{DI}_h(f,{ M }_{h},\hat{\mathbf{n}}) \Bigg| \leq 2^{d + 3} (d+1)^{\frac{3}{2}} \ Area (\partial M ) \Big( \text{Lip}(f) \sqrt{d + 1} \ h \\ \phantom{dfqfijeiofjqzdjiqjziieqioj} + ||f||_\infty \cdot {\color{mygreen}||\hat{\mathbf{n}} - \mathbf{n}||_{est}}\Big), $$ ($\text{DI}_h(f,{ M }_{h},\hat{\mathbf{n}})\approx$ summation of $f$ evaluated at each surfel $s$ and weighted by $\mu(s)$)

$\Rightarrow$ We need a multigrid convergent normal vector estimation

Remaining steps: We set $f(x) \EqDef \frac{1}{t_h (4 \pi t_h)^{\frac{d}{2}}} e^{- \frac{||{x} - {s}||^2}{4t_h}}(\tilde u({x}) - \tilde u({s}))$ and we derive bounds for $\text{Lip}(f)$ and $||f||_\infty$.

Main result



Let $\mathbf{s}$ be a surfel in $\partial_h M$, a function $u \in C^2(\partial M)$ and its extension $\tilde u$. Let $t_h = h^\alpha$ and let the convergence speed of the normal estimator be in $O(h^\beta)$. Let $h_0$ be the minimum between $\text{Diam}(\partial M)$, $R / \sqrt{d+1}$ and $K_3(d, \alpha, \text{Diam}(\partial M))$. For $0 < h \leq h_0$ we have
$$\\ \,\\ \lim\limits_{h \rightarrow 0} |{(\Delta u)(\xi({s})) - (L \tilde{u})({s})}| = 0 $$
if $0 < \alpha < \min\left( \frac{2}{d + 2}, \frac{2 \beta}{d + 1} \right)$.

In summary





  • Strong consistency thanks to a multigrid convergent normal vector field

  • Discrete operator is not as sparse as the cotangent one
  • Efficient implementation (convolutions on compact support)

Piecewise smooth reconstruction

Problem statement


$\Rightarrow$ Piecewise smooth reconstruction of the normal vector field

Ambrosio-Tortorelli functional





$$ \mathcal{AT}_\eps(u,v) = {\alpha} \underbrace{\dint_{{M}} |u - g|^{{2}} \dx}_{\text{attachment term}} + \underbrace{\dint_{{M}} |{v} \grad u|^{{2}} \dx}_{\text{smoothness term}} + \underbrace{{\lambda} \dint_{{M}} {\eps} |\grad {v}|^{{2}} + \frac{1}{4 {\eps}} |1 - {v}|^{{2}} \dx}_{\text{discontinuities length}}$$

  • Two functions to optimize: $u$ (scalar map) and $v$ (feature scalar map, $\mathcal{M}\rightarrow [0,1]$)
  • $v\approx 1$ on smooth parts, $v\approx 0$ near features
  • Quadratic terms
  • the AT functional $\Gamma-$converges to Mumford-Shah's functional $\mathcal{AT}_\eps \gammaconv{\eps \to 0} \mathcal{MS}$
  • (integration domain does not change, no Hausdorff measure)

Ambrosio-Tortorelli functional





$$ \mathcal{AT}_\eps(u,v) = {\alpha} \underbrace{\dint_{{M}} |u - g|^{{2}} \dx}_{\text{attachment term}} + \underbrace{\dint_{{M}} |{v} \grad u|^{{2}} \dx}_{\text{smoothness term}} + \underbrace{{\lambda} \dint_{{M}} {\eps} |\grad {v}|^{{2}} + \frac{1}{4 {\eps}} |1 - {v}|^{{2}} \dx}_{\text{discontinuities length}}$$

  • $\eps$ -- thickness of the feature set ($[distance]$)
  • ${\alpha}$ -- attachment coefficient to control the smoothing strength ($[area^{-1}]$)
  • $\lambda$ -- proportional to the length of discontinuities ($[distance^{-1}]$)

Discretization

"à la" Discrete Exterior Calculus [Hirani, Desbrun, Grady...]:

$\mathcal{M}$ is a cellular complex ($\mathcal{M}'$ its dual), $\sigma^k$ are $k-$cells of $\mathcal{M}$ (resp. $\sigma'^k$ of $\mathcal{M}'$)

      

  • $k-$forms are vectors of $|\{\sigma^k\}|$ scalars
  • Linear operators are matrices
  • e.g.
    • $d_k$ (exterior derivative) maps primal $k-$forms to primal $(k+1)-$forms
    • wedge product $\alpha\wedge \beta$ maps $k-$forms and $l-$forms to $(k+l)-$forms
    • Hodge-star $\star_k$ operator to maps primal $k-$forms to dual $k-$forms
    • ...

Discretization (bis)



$$ \mathcal{AT}_\eps(u,v) = {\alpha} \underbrace{\dint_{{M}} |u - g|^{{2}} \dx}_{\text{attachment term}} + \underbrace{\dint_{{M}} |{v} \grad u|^{{2}} \dx}_{\text{smoothness term}} + \underbrace{{\lambda} \dint_{{M}} {\eps} |\grad {v}|^{{2}} + \frac{1}{4 {\eps}} |1 - {v}|^{{2}} \dx}_{\text{discontinuities length}}$$

$v$ is a primal $0$-form, $u$ is a triple of dual $0$-forms $(u_1,u_2,u_3)$ and we discretize $\mathcal{AT}_\epsilon$

Discretization (ter)

$$\mathcal{AT}^d_\epsilon(u,v) \EqDef \alpha \sum_{i=1}^{3} \langle{u_i - g_i},{u_i - g_i}\rangle_\bar{0} + \sum_{i=1}^{3} \langle{v \wedge d_\bar{0} u_i},{v \wedge d_\bar{0} u_i}_\bar{1}\rangle_\bar{1} \\ \quad\quad\quad\quad\quad\quad\quad\quad\quad+ \lambda \epsilon \langle{ d_{0} v },{ d_{0} v }\rangle_1 + \frac{\lambda}{4\epsilon} \langle{1 - v},{1 - v}\rangle_0$$

  • if $\gamma$ is a primal $0-$form and $\beta$ a dual $1-$form, then $\gamma\wedge\beta=diag(\beta)M\gamma$ with $M\EqDef\frac{1}{2}|d_0|$
  • $A$ is the matrix form of $d_0$
  • $B$ is the matrix form of $d_\bar{0}$
  • $\vu_i,\,\vv,\,\vg$ are column vectors containing associated $k-$form scalars
  • $S_i$ is a diagonal matrix encoding the Hodge star $\star_i$

$$\mathcal{AT}^d_\epsilon(\vu,\vv)= \alpha (u-g)^T S_{\bar{0}} (u-g) + u^T B^T \operatorname{Diag}(M v) S_{\bar{1}} \operatorname{Diag}(M v) B u \\ + \lambda\,\epsilon\, v^T A^T S_1 A v + \frac{\lambda}{4\epsilon} (1-v)^T S_0 (1-v)$$

Optimization

$\min_{u,v} \mathcal{AT}_\eps \Leftrightarrow \left( \nabla_u \mathcal{AT}_\eps = 0 \right) \wedge \left(\nabla_v \mathcal{AT}_\eps = 0 \right)$

For a given $\epsilon$:


$$ \nabla_u AT_\epsilon[u, v] = 0 \Leftrightarrow \left[ \alpha\,S_{\bar{0}} - B^T \operatorname{Diag}(M v) S_{\bar{1}} \operatorname{Diag}(M v) B \right] u = \alpha\,S_{0}\,g $$

$$\nabla_v AT_\epsilon[u, v] = 0 \Leftrightarrow \left[ \frac{\lambda}{4\epsilon} S_0 + \lambda\,\epsilon\,A^T S_1 A + M^T \operatorname{Diag}(B u) S_{\bar{1}} \operatorname{Diag}(B u) M \right] v = \frac{\lambda}{4\,\epsilon} S_0$$

$\Rightarrow$ only linear system solves on sparse matrices !

$\lambda$ parameter





$\epsilon$ parameter

Noise level w.r.t. $\alpha$ parameter

In summary





  • Anisotropic normal vector field regularization with feature selection
  • Sharp features
  • Parameters make sense :)
  • Variational problem discretization using a combinatorial representation of the digital surface

Chicken/egg problem: measure of quads $\mu(s)$ used in the DEC operators

Ambrosio-Tortorelli on meshes





Nicolas Bonneel, C., Pierre Gueth, Jacques-Olivier Lachaud. Mumford-Shah Mesh Processing using the Ambrosio-Tortorelli Functional. Computer Graphics Forum (Proceedings of Pacific Graphics), 37(7), October 2018.

DGtal / dgtal.org

github.com/DGtal-team

@libdgtal

Conclusion




  • Nice geometrical model with many interactions (arithmetic's, theory of words, computational geometry, discrete mathematics...)
  • Very specific discrete/continuous properties
  • Related to various areas (image processing, material sciences, geometrical modeling, rendering...) data