A consistent analytical formulation for volume-estimation of geometries enclosed by implicitly defined surfaces

TL;DR

This paper introduces a general analytical method for accurately estimating the volume of geometries enclosed by implicit surfaces, ensuring consistency across mesh refinements with demonstrated second-order accuracy.

Contribution

It provides a unified analytical formulation for 2D and 3D volume estimation that addresses mesh refinement inconsistencies.

Findings

Achieves 2nd-order accuracy in volume estimation

Applicable to all 2D cases and elementary 3D cases

Ensures consistency across mesh refinements

Abstract

We have derived an analytical formulation for estimating the volume of geometries enclosed by implicitly defined surfaces. The novelty of this work is due to two aspects. First we provide a general analytical formulation for all two-dimensional cases, and for elementary three three-dimensional cases by which the volume of general three-dimensional cases can be computed. Second, our method addresses the inconsistency issue due to mesh refinement. It is demonstrated by several two-dimensional and three-dimensional cases that this analytical formulation exhibits 2nd-order accuracy.

Equations

$$\Gamma =\left\{\left(x,y\right)\in \Omega :\varphi (x,y)={\beta }_0+{\beta }_{1}x+{\beta }_{2}y+{\beta }_{3}xy=0\right\},$$

$${\beta }_0={\varphi }_{00},{\beta }_1={\varphi }_{10}-{\varphi }_{00},{\beta }_2={\varphi }_{01}-{\varphi }_{00},{\beta }_3={\varphi }_{00}+{\varphi }_{11}-{\varphi }_{01}-{\varphi }_{10},$$

$$\zeta =\frac{-a\eta -b}{c\eta +d}$$

$${\mathbf{x}}_0=\underset{\mathbf{x}\in \mathbf{X}}{\text{arg max }}\#\left({\mathbf{x}}^{\prime }\in \mathbf{X}|\varphi (\mathbf{x})\varphi ({\mathbf{x}}^{\prime })<0\wedge \Vert {\mathbf{x}}^{\prime }-\mathbf{x}{\Vert }_2=1\right)$$

$$l_x=-\frac{{\beta }_2{i}_y+{\beta }_0}{{\beta }_3{i}_y+{\beta }_1}\text{mbox}if{\varphi }_{i_x{i}_y}{\varphi }_{|1-i_x|i_y}<0\text{mbox}else1.0,l_y=-\frac{{\beta }_1{i}_x+{\beta }_0}{{\beta }_3{i}_x+{\beta }_2}\text{mbox}if{\varphi }_{i_x{i}_y}{\varphi }_{i_x|1-i_y|}<0\text{mbox}else1.0,$$

$$\alpha ={(1-2i_{\varsigma })\left[\frac{(ad-bc)\log |c\eta +d|}{c^2}-\frac{a\eta }{c}\right]|}_{{\eta }_0}^{{\eta }_1}+i_{\varsigma }({\eta }_1-{\eta }_0).$$

$$\Gamma =\{\left(x,y,z\right)\in \Omega :\varphi (x,y,z)={\beta }_0+{\beta }_{1}x+{\beta }_{2}y+{\beta }_{3}z+{\beta }_{4}xy+{\beta }_{5}yz+{\beta }_{6}xz+{\beta }_{7}xyz=0\},$$

$${\beta }_0={\varphi }_{000},{\beta }_1={\varphi }_{100}-{\varphi }_{000},{\beta }_2={\varphi }_{010}-{\varphi }_{000},{\beta }_3={\varphi }_{001}-{\varphi }_{000},$$

$${\beta }_4={\varphi }_{110}-{\varphi }_{100}-{\varphi }_{010}+{\varphi }_{000},{\beta }_5={\varphi }_{011}-{\varphi }_{010}-{\varphi }_{001}+{\varphi }_{000},$$

$${\beta }_6={\varphi }_{101}-{\varphi }_{100}-{\varphi }_{001}+{\varphi }_{000},{\beta }_7={\varphi }_{111}-{\varphi }_{110}-{\varphi }_{101}+{\varphi }_{100}-{\varphi }_{011}+{\varphi }_{010}+{\varphi }_{001}-{\varphi }_{000}.$$

$$l_x$$

$$l_y$$

$$l_z$$

$$\zeta (\xi ,\eta )=\frac{-\xi (a\eta +b)-c\eta -d}{\xi (e\eta +f)+g\eta +h},$$

$$\alpha =(1-2i_{\varsigma }){\int }_{{\xi }_0}^{{\xi }_1}\left[{\int }_{{\eta }_0}^{{\eta }_1}\zeta (\xi ,\eta )d\eta \right]d\xi +i_{\varsigma }{\int }_{{\xi }_0}^{{\xi }_1}q(\xi )d\xi =(1-2i_{\varsigma })\left[\mathbf{F}\left({\xi }_1\right)-\mathbf{F}\left({\xi }_0\right)\right]+i_{\varsigma }\left[\mathbf{G}\left({\xi }_1\right)-\mathbf{G}\left({\xi }_0\right)\right]$$

$$q(\xi )=\frac{-(f{i}_{\varsigma }+b)\xi -h{i}_{\varsigma }-d}{(e{i}_{\varsigma }+a)\xi +g{i}_{\varsigma }+c}.$$

$${\eta }_0=-\frac{\left(f{i}_{\zeta }+b\right)\xi +h{i}_{\zeta }+d}{\left(e{i}_{\zeta }+a\right)\xi +g{i}_{\zeta }+c}\text{mbox}if{i}_{\eta }=1\text{mbox}else0.0,{\eta }_1=-\frac{\left(f{i}_{\zeta }+b\right)\xi +h{i}_{\zeta }+d}{\left(e{i}_{\zeta }+a\right)\xi +g{i}_{\zeta }+c}\text{mbox}if{i}_{\eta }=0\text{mbox}else1.0,$$

$$\mathbf{F}(\xi )=\int \left[-\frac{\left(t_2+\xi t_1+{\xi }^2{t}_0\right)}{{(e\xi +g)}^2}\left[\log \left(f\xi +h\right)+\log \left(a\xi +c\right)-\log \left(t_2+\xi t_1+{\xi }^2{t}_0\right)\right]+\frac{b\xi +d}{e\xi +g}\right]d\xi$$

$$\mathbf{G}(\xi )=-\frac{1}{(i_{\zeta }e+a{)}^2}\left[\left(eh-fg\right)i_{\zeta }^2+\left(ah-bg-cf+de\right)i_{\zeta }+ad-bc\right]\log \left(i_{\zeta }{\xi }_0+{\xi }_4\right)-\frac{\left(i_{\zeta }f+b\right)\xi }{i_{\zeta }e+a}$$

$${\mathbf{F}}_1(\xi )$$

$${\mathbf{F}}_2(\xi )$$

$${\mathbf{F}}_1(\xi )=\underset{\mathbf{A}}{\underbrace{\log ({\Xi }_1)\left[\frac{{\Xi }_0\log ({\Xi }_0)t_6-t_3}{e^3{\Xi }_0}+\frac{\xi t_0}{e^2}\right]}}+\underset{{\mathbf{B}}_1}{\underbrace{\frac{t_6}{e^3}\int \frac{(t_1+2\xi t_0)\log ({\Xi }_0)}{{\Xi }_1}d\xi }}+\underset{{\mathbf{B}}_2}{\underbrace{\int \frac{(t_1+2\xi t_0)\left[-t_3+e\xi {\Xi }_0{t}_0\right]}{e^3{\Xi }_0{\Xi }_1}d\xi }},$$

$${\mathbf{B}}_1=\frac{t_6}{e^3}\text{bigintsss}\log ({\Xi }_0)\left(\frac{2t_0}{2t_0\xi +\sqrt{t_1^2-4t_0{t}_2}+t_1}+\frac{2t_0}{2t_0\xi -\sqrt{t_1^2-4t_0{t}_2}-t_1}\right)d\xi ,$$

$${\mathbf{B}}_1=\frac{t_6}{e^3}\text{bigintsss}\log ({\Xi }_0)\left(\frac{2t_0}{2t_0\xi +i\sqrt{-t_1^2+4t_0{t}_2}+t_1}+\frac{2t_0}{2t_0\xi -i\sqrt{-t_1^2+4t_0-t_1}2t_0}\right)d\xi .$$

$${\mathbf{B}}_1$$

$$\mathbf{C}$$

$$\mathbf{D}$$

$$\mathbf{E}=\int \frac{b\xi +d}{e\xi +g}d\xi =\frac{b\xi }{e}-\frac{(bg-de)\log ({\Xi }_0)}{e^2},$$

$$\alpha$$

$${\mathbf{C}}^{\ast }$$

$${\mathbf{D}}^{\ast }$$

Full text

††journal: Elsevier

1 Introduction

Accurate and efficient volume estimation is important for investigating complex geometries encountered in a variety of scientific and engineering problems. In astronomy, the estimation of an asteroid volume is helpful to determine the bulk density for assessing the material composition. Except for some large spheroidal asteroids whose volumes can be estimated by their mean diameters, most asteroids have irregular shapes which require more complicated volume estimation methods. Usually, three-dimensional (3D) reconstructed shape models are employed to approximate such irregular geometries [7, 1]. In computational fluid dynamics, a volume-fraction estimation is required for the transport of multiple fluids by the volume-of-fluid (VOF) method or the level-set method [12]. The volume fraction of a specific fluid in a uniform mesh cell can be computed by the Heaviside [23, 25] or color function [9], or combination of either one with the interface normal direction [11]. The interface reconstruction step of the VOF method may also require a volume-estimation method to determine the interface location [22, 5, 4]. Another example is medical engineering, where in order to detect the development of lung hypoplasia and investigate the correlation between lung growth and fetal presentation, accurate measurement of fetal lung volume is a pressing need [21]. Tumor volume estimation is also widely used in many cancer treatments, such as prostate cancer [3], brain tumors [13] and pelvic neoplasms [18]. Most of such volume estimations are based on magnetic-resonance (MR) imaging. A three-dimensional reconstruction method is considered to be much more reliable and accurate than other more traditional approximation methods, such as ellipsoid and disc-summation approximations, and automated image segmentation [6, 8, 2, 14]. For example, volume estimation of the articular cartilage is performed by the Marching Cube algorithm [17, 19] applied on the three-dimensional reconstruction [20]. The disadvantage of this method is that one first needs to reconstruct the iso-surface before volume estimation. However, there are too many different cases ($16$ for 2D and $256$ for 3D) to deal with. The volume estimation method is consistent during mesh refinement in the sense that the volume fractions are conservative with respect to the coarsest level. It is possible to obtain a unique volume approximation during mesh refinement by a consistent volume estimation method. Any numerical method introduces inconsistencies during mesh refinement, irrespectively of its order of accuracy, as the truncation error depends on mesh resolution. Although not critical for many applications, consistency is desirable, and it is important for numerical simulations of multiphase flows [24].

The objective of the present paper is to develop a consistent volume estimation method for calculating the volumes of geometries enclosed by implicitly defined surfaces. To achieve this property, we derive analytical formulations for generic 2D and 3D geometries, which avoid truncation errors, unlike numerical methods. The formulations are then subsumed to derive a general formulation that can be applied to more complex cases. This method is made consistent during mesh refinement by employing bilinear/trilinear interpolation. The paper is organized as follows. Sec. 2 describes the derivation of 2D and 3D analytical volume formulation for arbitrary geometries. We also discuss how to preserve consistency during mesh refinement. Sec. 3 is dedicated to assessing the capability of the present method to calculate the volumes for various cases, followed by conclusions in Sec. 4.

2 Analytical Volume Fraction

2.1 2D formulation

Before presenting the 3D formulation, we first consider the 2D case for its simplicity and give a comprehensive description of the underlying concept. Suppose that we have volumetric data, say discrete data of a level-set field $\phi$, which implicitly defines the interface in the domain $\Omega$, the interface in any cell can be expressed by the zero value of a reconstructed field $\phi(x,y)$ which is obtained by bilinear interpolation. The interface is hereby represented by piecewise parabolas for two dimensions (2D), i.e., the interface inside the cell $[i,i+1]\times[j,j+1]$ is defined as

[TABLE]

where the coefficients $\beta_{0}$, $\beta_{1}$, $\beta_{2}$ and $\beta_{3}$ are determined by the interpolation conditions,

[TABLE]

with $\phi_{00}$, $\phi_{10}$, $\phi_{01}$ and $\phi_{11}$ defined in Fig. 1(a). The volume of this cell can be easily calculated by integrating the reconstructed field in $x$ direction or $y$ direction. In order to simplify the integration, we reformulate Eq. (1) with a coordinate mapping,

[TABLE]

whose symbols need careful definition to simplify the calculations. The definition is based on the intersection points at the edges of a mesh cell. Specifically in 2D, we first detect the grid point

[TABLE]

and set is as the origin of a local Cartesian coordinate system. $\mathbf{X}$ is the set of all $4$ vertices ($8$ in 3D) of the current cell. Then the coordinate transformation is determined by $\eta=\underset{\gamma\in\{x,y\}}{\text{arg max }}l_{\gamma}$ and $\zeta=\left\{x,y\right\}\backslash\left\{\eta\right\}$, where $l_{\gamma}$ is the distance from the vertex $(i,j)$ to the intersection point on a neighbour edge in $x$-direction or $y$-direction,

[TABLE]

where $i_{x}$ and $i_{y}$ are the coordinates of the origin in the $x$ and $y$ directions, respectively, see Fig. 1(a). Once the origin and the local transformed coordinates ($\eta,\zeta$) are defined, the coefficients of Eq. (3), $a$, $b$, $c$ and $d$, are determined according to the relation between the $(x,y)$ and $(\eta,\zeta)$ coordinate systems. Then the cut volume of cell $[i,i+1]\times[j,j+1]$ can be obtained by integration of $\zeta(\eta)$ as

[TABLE]

The integration range is determined by $\eta_{0}=\mbox{min}(l_{\eta},i_{\eta})$ and $\eta_{1}=\mbox{max}(l_{\eta},i_{\eta})$, where ${l}_{\eta}$ is the length from the origin of the coordinate system to the intersection point in $\eta$ direction.

2.2 3D formulation

The 2D formulation can be easily extended to three dimensions. Again, the interface inside a cell, say $[i,i+1]\times[j,j+1]\times[k,k+1]$ in Fig. 1(b), can be represented by trilinear interpolation,

[TABLE]

where the coefficients are uniquely determined by the interpolation conditions, see Fig. 1(b),

[TABLE]

Upon identification of the origin ($O$) by Eq. (4), the local coordinates are defined as $\xi=\underset{\gamma\in\{x,y,z\}}{\text{arg max }}l_{\gamma}$, $\eta=\underset{\gamma\in\{x,y,z\}\setminus\{\xi\}}{\text{arg max }}l_{\gamma}$ and $\zeta=\{x,y,z\}\setminus\{\xi,\eta\}$, where the distance from the grid $(i,j,k)$ to the intersection points, $l_{\gamma}$, are calculated by

[TABLE]

The coordinates of the origin in the $x$, $y$ and $z$ directions are $i_{x}$, $i_{y}$ and $i_{z}$, respectively, as shown in Fig. 1(b). Eq. (7) is rewritten as

[TABLE]

with the coefficients, $a,b,\dots,f,g$, determined by the relation between the $(x,y,z)$ and $(\xi,\eta,\zeta)$ coordinate systems. Thus the volume of each cell $\alpha$ can be calculated by the double integral of $\zeta(\xi,\eta)$,

[TABLE]

with

[TABLE]

The integration ranges are $\xi_{0}=\mbox{min}(l_{\xi},i_{\xi})$ and $\xi_{1}=\mbox{max}(l_{\xi},i_{\xi})$ in $\xi$-direction, and

[TABLE]

in $\eta$-direction.

For all cells containing a resolved interface, Eq. (11) for the three cases sketched in Fig. 2 can be expressed by elementary integrals. First we subsume those elementary cases by a single formulation. Other cases can be expressed as combination of these elementary cases, and their volume estimation will be discussed later.

First we consider $i_{\zeta}=0$ such that $\mathbf{F}(\xi)$ in Eq.(11) becomes

[TABLE]

with $t_{0}=af-be$, $t_{1}=ah-bg+cf-de$, and $t_{2}=ch-dg$, while $\mathbf{G}(\xi)$ is

[TABLE]

We define $\mathbf{F}(\xi)=\mathbf{F}_{1}(\xi)+\mathbf{F}_{2}(\xi)$, where $\mathbf{F}_{1}$ and $\mathbf{F}_{2}$ are

[TABLE]

and

[TABLE]

respectively. Substituting the variables, ${{\Xi}_{0}}=e\xi+g$, ${{\Xi}_{1}}={{t}_{2}}+\xi{{t}_{1}}+{{\xi}^{2}}{{t}_{0}}$, ${{t}_{3}}={{e}^{2}}{{t}_{2}}-eg{{t}_{1}}+{{g}^{2}}{{t}_{0}}$ and ${{t}_{6}}=e{{t}_{1}}-2g{{t}_{0}}$, into the above formula, we obtain

[TABLE]

where if $t_{1}^{2}-4{{t}_{0}}{{t}_{2}}\geq 0$,

[TABLE]

and otherwise

[TABLE]

The integration of Eq. (19) results in

[TABLE]

with $\Xi_{2}=t_{1}+2\xi t_{0}$, $s=\sqrt{\left|t_{1}^{2}-4{{t}_{0}}{{t}_{2}}\right|}$, $s_{1}=t_{6}+es$ and $s_{2}=t_{6}-es$. The dilogarithm (Spence’s function) is defined as $\text{Li}_{2}(z)=-\int_{0}^{z}{\log(1-t)dt}$. Eq. (20) becomes

[TABLE]

Similarly, in Eq. (18) we have $\mathbf{B}_{2}=\frac{1}{{{e}^{3}}}\left[-\log({{\Xi}_{0}}){{t}_{6}}-g\log({{\Xi}_{1}}){{t}_{0}}+2e\xi{{t}_{0}}+es\,\log\left(\frac{{{\Xi}_{2}}-s}{{{\Xi}_{2}}+s}\right)\right]$ if $t_{1}^{2}-4{{t}_{0}}{{t}_{2}}\geq 0$, otherwise $\mathbf{B}_{2}=\frac{1}{{{e}^{3}}}\left[-\log({{\Xi}_{0}}){{t}_{6}}-g\log({{\Xi}_{1}}){{t}_{0}}+2e\xi{{t}_{0}}-2es\,\text{atan}(\frac{\Xi_{2}}{s})\right]$.

To calculate $\mathbf{F}_{2}(\xi)$, we define $\mathbf{F}_{2}(\xi)=\mathbf{C}+\mathbf{D}+\mathbf{E}$, with $\mathbf{C}$, $\mathbf{D}$ and $\mathbf{E}$ being

[TABLE]

[TABLE]

and

[TABLE]

respectively, with ${{t}_{4}}=eh-fg$, ${{t}_{5}}=ce-ag$, ${{\Xi}_{3}}=f\xi+h$ and ${{\Xi}_{4}}=a\xi+c$. We can generalize Eq. (11) for cases with different $i_{\zeta}$,

[TABLE]

where $\mathbf{F}^{*}\left(\xi\right)=\left(1-2{{i}_{\eta}}\right)\left(\mathbf{A^{*}}-{{\mathbf{B}}_{\mathbf{1}}}^{*}-{{\mathbf{B}}_{\mathbf{2}}}^{*}\right)+\left(1-2{{i}_{\eta}}\right)\left(\mathbf{C}^{*}+\mathbf{D}^{*}+\mathbf{E}^{*}\right)+\mathbf{G}^{*}$ and $\mathbf{G}\left(\xi\right)^{*}=\left(1-2{{i}_{\eta}}\right)\mathbf{G}\left(\xi\right)$, where ${{\mathbf{A}}}^{*}$, ${{\mathbf{B}}_{\mathbf{1}}}^{*}$ and ${{\mathbf{B}}_{\mathbf{2}}}^{*}$ are the same as ${{\mathbf{A}}}$, ${{\mathbf{B}}_{\mathbf{1}}}$ and ${{\mathbf{B}}_{\mathbf{2}}}$, respectively. ${\mathbf{C}^{*}}$, ${\mathbf{D}^{*}}$, ${\mathbf{E}^{*}}$ and ${\mathbf{G}^{*}}$ are defined as

[TABLE]

[TABLE]

[TABLE]

and ${\mathbf{G}^{*}}=\left(1-2{{i}_{\eta}}\right){\mathbf{G}}+\xi{{i}_{\eta}}$, respectively. Now, we can use this general formulation to calculate 3D volume fraction of all cases sketched in Fig. 2. These elementary cases generate more complex cases as indicated in Fig. 3, and are defined as type I, II, and III, respectively. These complex cases can be split them into elementary cases in Fig. 2, as illustrated in each subgraph. For instance, the case in Fig. 3(a) is composed of type I and II, whereas cases in Fig. 3(b) and (c) are combinations of all three elementary types. In order to compute the volumes of those cases, the first step is to decompose the cell into elementary types, which is easy to accomplish. For example, in Fig. 3(b), along the $\xi$-direction we can find two intersection points $\xi_{1}$ and $\xi_{2}$ which occur at two different edges parallel to $\xi$-direction. By these intersection points, the whole cell is cut into three parts in $\xi$ directions, as shown in Fig. 3(b).

2.3 Accuracy test

We calculate the volumes of a circle and a sphere to assess the accuracy of our analytic formulations in 2D and 3D. The results of our method are compared to those of a linear approximation method [15], as plotted in Fig. 4.

For 2D, as shown in Fig. 4(a), the cumulative volumes along $x$-direction indicate that both linear and analytical results are in good agreement with the exact result while locally the analytical formulation is more accurate. In Fig. 4(b) 2nd-order convergence rate is demonstrated for our formulation because it is based on bilinear/triliner interpolation while the linear method is 1st-order as expected. The magnitude of the error is much smaller in our method. The result of the 3D case in Fig. 4(c) exhibits similar error distribution with Fig. 4b. As a consequence, for 2D and 3D, our analytical formulation provides more accurate volume-estimation than the linear method due to a more accurate interface representation.

2.4 Consistency test

Definition. The consistency condition of the volume-estimation during mesh refinement is $\int_{\Omega}\alpha^{0}dV=\int_{\Omega}\alpha^{\ell}dV$, i.e., volumes calculated from the initial volumetric field and the refined field are exactly the same.

The significance of maintaining consistency of volume estimation is obvious. For example, for an adaptive multiphase solver [24] conservation errors can only be prevented if the volumes have an identical value at each resolution level. This can not be achieved by a linear representation. All such numerical volume estimations are not consistent as the errors introduced by numerical methods depend on grid resolution.

Theorem. The analytical volume-estimation formulation is consistent if the level-set field is refined by bilinear/trilinear interpolation.

Proof 1

The necessary and sufficient condition for consistency is that the parameters of Eq. (1) or Eq. (7) remain invariant upon mesh refinement as no numerical error is introduced by the analytical formulation and the coordinate transformation doest not change the computed volume. This can be easily proved. For simplicity, we consider the 2D cases for description. Assume the grid point level-set values $\phi^{\ell}$ of a cell $\left[i,i+1\right]\times\left[j,j+1\right]$ on level $\ell$ are given. During refinement, a cell is divided into $4$ subcells. The level-set values $\phi^{\ell+1}$ at level $\ell+1$ are computed by bilinear interpolation. Thus the parameters of the interface equation for the subcell $\left[2i+1,2i+2\right]\times\left[2j+1,2j+2\right]$ are

[TABLE]

which can be expressed by $\phi^{\ell}$ with bilinear interpolation,

[TABLE]

This implies that the interface equation at $\ell+1$,

[TABLE]

is the same as Eq. (1), as $x^{\prime}=2\,x-1$ and $y^{\prime}=2\,y-1$. For the other subcell we can obtain the same conclusion, which completes the proof. $\square$

For example, a single cell test with specific level-set values clearly shows the consistency feature of our method. After prescribing the grid values, $\phi^{\ell}_{00}$ = 0.1, $\phi^{\ell}_{10}$ = 0.6, $\phi^{\ell}_{01}$ = -0.3, and $\phi^{\ell}_{11}$ = -0.1, we obtain the coefficients $a=0.5$, $b=0.1$, $c=-0.3$, and $d=-0.4$. The origin is placed at $(0,0)$ and the integration variable $\eta$ is $x$, hence the interface profile becomes $\zeta=(5\eta+1)/(3\eta+4)$ with an integration range of $[0,1]$. After substituting the parameters into the Eq. (6), the volume of this piece on this level is $\alpha=\frac{1}{9}(17\,\mathrm{log}(4)-17\,\mathrm{log}(7)+15)$. The bilinear interpolation leads to $\phi^{\ell+1}_{00}$ = 0.1, $\phi^{\ell+1}_{10}$ = 0.35, $\phi^{\ell+1}_{01}$ = -0.1, and $\phi^{\ell+1}_{11}$ = 0.075 for subcell cell $\left[2i,2i+1\right]\times\left[2j,2j+1\right]$. Substitute these values into Eq. (6), we can get a volume of $\alpha_{00}=\frac{1}{63}(120-476\,\mathrm{log}\left(\frac{68}{7}\right))+\frac{68}{9}\,\mathrm{log}\left(8\right)+\frac{3}{7}$. Analogously, the volumes of the other three subcells are $\alpha_{10}=1.0$, $\alpha_{01}=\frac{68}{9}(\,\mathrm{log}\left(\frac{68}{7}\right)-12-68\,\mathrm{log}\left(11\right)+21)$ and $\alpha_{11}=\frac{68}{9}(\,\mathrm{log}\left(11\right)-\,\mathrm{log}\left(14\right)+21)$, respectively. Obviously, the consistency is satisfied as $\alpha=(\alpha_{11}+\alpha_{10}+\alpha_{01}+\alpha_{00})/4$. Analogously, we can proof the consistency of our volume-estimation method during coarsening.

In addition, simple cases used in the accuracy test are employed to test the consistency of the present method. As shown in Fig. 5, the initial resolution is $h=0.05$ on $\ell=0$, then the level-set fields are refined from $\ell_{max}=1$ to $5$. The dash-dotted, solid and dashed lines stand for error norms,

[TABLE]

respectively, with the superscripts ‘[math]’ and ‘$\ell$’ being the volume fraction of the cell $[i,j]\times[i+1,j+1]$ on $\ell=0$ and the summation of volume fractions on $\ell>0$ within the same cell. In 2D, with $\ell_{\text{max}}$ increasing, the discrepancy between $\ell=0$ and $\ell_{\text{max}}$ for linear approximation method is large, ranging from $5\times 10^{-2}$ to $1\times 10^{-1}$ for $L_{\infty}$. The linear reconstruction inside each cell is responsible for this distinct errors as it does not have a consistent representation of interface between different resolutions. With respect to the analytical formulation the magnitudes are much lower than those of the linear method. So the bilinear/trilinear representation and interpolation are beneficial to preserving the consistency of the volume estimation. The 2D and 3D errors produced in our method can be neglected. Therefore it is concluded that for both 2D and 3D our method automatically achieves consistency which is attributed to the analytical formulation and the bilinear/trilinear interpolation.

2.5 Additional remarks

Note that only resolved cases have an analytical formulation. If there exists ambiguity such as multiple intersections in one cell edge, pure analytical method may not provide the correct results as the formulation is not able to resolve this cell. One can easily address this issue by adaptive refinement in the ambiguity regions or performing an interface scale separation [10]. 2. 2.

Our method achieves 2nd-order accuracy which is generally acceptable for most of the applications. However the goal of our work is not to achieve high-order accuracy, but to maintain the consistency which is the unique feature of our method.

3 Applications

In this section, we apply the analytical formulation to more complicated cases. Firstly we give two cases which have exact volumes to test both accuracy and consistency. The first one is refered to as “double-circle” case where two circles with an identical radius $r=0.25$ and different centers ($\textbf{x}_{0}=(0.3,0.5)$ and $\textbf{x}_{1}=(0.7,0.5)$) are merged. While the Zalesak disk with a radius of $0.4$, a notch width of $0.2$ and a notch height of $0.6$ is placed in at $(0.5,0.5)$. As shown in Figs. 6a and 6c, the results converge with a 2nd-order rate in our method and a 1st-order rate in the liner method. Both these cases show distinguished consistency in the analytical formulation and large inconsistency in the linear method. A 3D double-sphere case with the same parameters in the double-circle case confirms the conclusion, as shown in Fig. 7. In Fig. 8, A case with $15$ randomly generated circles provides similar convergence and consistence results with the above cases.

Further more, a more practical case is tested to manifest the features of the present method. As shown in Fig. 10, this case is the image segmentation of a human brain wihch has more complicated boundaries than the above cases. The level-set fields are generated by the code in Ref. [16] with the same MR image. The red line represents the air-tissue interface while the blue one stands for the interface of two different tissues. The consistency results are presented in Fig. 11. Apparently, the errors in both areas are significantly reduced by the analytical formulation, indicating the consistency of our method in this brain image segmentation.

4 Concluding remarks

In this paper we have proposed a volume-estimation method for implicitly defined geometries. We have derived a general analytical formulation for 2D cases and three elementary 3D cases which can be employed for more complicated 3D cases. The analytical formulation achieves 2nd-order accuracy for arbitrary geometries. Based on this analytical formulation, the consistency issue during mesh refinement is imposed by the bilinear/trilinear interpolation. The accuracy and consistency are demonstrated by 2D and 3D examples. This consistent analytical volume-estimation method has the potential to be applied in many scientific and engineering fields, such as astronomy, computational fluid dynamics and medical imaging. Based on the idea developed in this paper, a consistent analytical volume-estimation formulation with high-order accuracy and its extension to arbitrary convex polyhedra is subject to further work.

Acknowledgment

This work is supported by China Scholarship Council under No. 201306290030.