## 3 Self-intersection Penalty Term
### Gradients Calculation  
The self-intersection penalty term \(E_{SPT}(V)\) is defined as:  
\[E_{SPT}(V)=V_{intersect} \tag{1}\]  
where \(V\) denotes the set of coordinates of all vertices, \(V_{intersection}\) is the volume of the self-intersection region. For the sake of further discussion, a vertex \(p_{t}(t=1,\dots,N)\) is randomly chosen from the triangular mesh with \(N\) vertices. Coordinates of all vertices are frozen, except for \(p_{t}\). The coordinate of \(p_{t}\) is denoted as \((x_{t},y_{t},z_{t})\). Under the preceding assumptions, the penalty term \(E_{SPT}(V)\) can be regarded as a function of \(x_{t}\),\(y_{t}\) and \(z_{t}\). \[E_{SPT}(V)=f(x_{t},y_{t},z_{t}) \tag{2}\]  
To compute \(\frac{\partial E_{SPT}(V)}{\partial x_{t}}\), \(\frac{\partial E_{SPT}(V)}{\partial y_{t}}\) and \(\frac{\partial E_{SPT}(V)}{\partial z_{t}}\), \(p_{t}\) is displaced from \((x_{t},y_{t},z_{t})\) to \((x_{t}+\Delta x,y_{t}+\Delta y,\)\(z_{t}+\Delta z)\). Then the change in \(E_{STP}(V)\) due to the displacement could be represented as:  
\[\Delta E_{SPT}(V)=f(x_{t}+\Delta x,y_{t}+\Delta y,z_{t}+\Delta z)-f(x_{t},y_{t },z_{t}) \tag{3}\]  
Further computation is hard to continue without imposing constraints on \(p_{t}\). First, for the simplest case, \(p_{t}\) is assumed to belong to \(V_{0}\), that is to say \(p_{t}\in V_{0}\). Since tiny displacement of \(p_{t}\) brings no effect to the volume of self-intersection, it is obvious that:  
\[\Delta E_{SPT}(V)=\Delta f=0,p_{t}\in V_{0} \tag{4}\]  
The gradient of \(E_{SPT}(V)\) could be represented as:  
\[\left( \begin{array}{c}\frac{\partial E_{SPT}(V)}{\partial x_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial y_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial z_{t}}\end{array} \right)=\left( \begin{array} []{c}\lim_{\Delta x\to 0}\frac{\Delta f}{\Delta x}\\ \lim_{\Delta y\to 0}\frac{\Delta f}{\Delta y}\\ \lim_{\Delta z\to 0}\frac{\Delta f}{\Delta z}\end{array} \right)=\left( \begin{array}{c}0\\ 0\\ 0\end{array} \right),p_{t}\in V_{0} \tag{5}\]  
For the second case \(p_{j}\in V_{out}\), more assumptions need to be made for a detailed discussion. We assume that there are \(n\) neighboring triangles sharing \(p_{t}\) as a common vertex. One of the neighboring triangles is denoted as \(T_{l}(l=1,\ldots,n)\), the area of \(T_{l}\) is represented as \(S_{l}\) and the unit normal vector of \(T_{l}\) is denoted as \(\mathbf{n}_{l}\). An intuitive representation of this situation is shown in Figure 5.  
The change in the volume of self-intersection due to tiny displacement of vertex \(p_{t}\) can be represented as:  
\[\Delta f=\Delta V_{intersection}=\frac{1}{3}\sum_{l=1}^{n}S_{l}\mathbf{n}_{l}( \Delta x\mathbf{i}+\Delta y\mathbf{j}+\Delta z\mathbf{k}),p_{t}\in V_{out} \tag{6}\]  
where \(\mathbf{i}\),\(\mathbf{j}\) and \(\mathbf{k}\) are unit vectors in the same directions as the positive directions of \(x\), \(y\) and \(z\) axes. The gradient can be obtained as:  
\[\left( \begin{array}{c}\frac{\partial E_{SPT}(V)}{\partial x_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial y_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial z_{t}}\end{array} \right)=\left( \begin{array} []{c}\lim_{\Delta x\to 0}\frac{\Delta f}{\Delta x}\\ \lim_{\Delta y\to 0}\frac{\Delta f}{\Delta y}\\ \lim_{\Delta z\to 0}\frac{\Delta f}{\Delta z}\end{array} \right)=\left( \begin{array}{c}\frac{1}{3}\sum_{l=1}^{n}S_{l}\mathbf{n}_{l}\mathbf{i}\\ \frac{1}{3}\sum_{l=1}^{n}S_{l}\mathbf{n}_{l}\mathbf{k}\end{array} \right),p_{t}\in V_{out} \tag{7}\]  
Figure 5: A tiny displacement on the vertex \(p_{t}\). The change in the volume of self-intersection is equal to the sum of volume of several neighboring tetrahedrons. The equation above can be simplified as:  
\[\left( \begin{array}{c}\frac{\partial E_{SPT}(V)}{\partial x_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial y_{t}}\\ \frac{\partial E_{SPT}(V)}{\partial z_{t}}\end{array}