Review May 2023

An efficient iterative thresholding method for image segmentation

The paper: An efficient iterative thresholding method for image segmentation was referenced by one of the reviewers to indicate that it is possible to succesfully include a data term in a threshold dynamic framework.

Indeed, this work uses threshold dynamics to optimize the length component of the Chan-Vese model. Nonetheless, the data term is optimization of the data term is straightforward and does not make use of threshold dynamics. In fact, the energy is splitted in two optimization problems that are evaluated separately (taken as input the solution of the other problem).

Paper Summary

The paper proposes a model of image segmentation. It is based on the Chan-Vese model. The formulation differs in how the length component is computed. The multi-phase Chan-Vese model can be written as

\[ \displaystyle \begin{align*} (u^{\star}, C^{\star}) &= \text{argmin}_{u \in S, C \in \mathbb{R}^n}{ \sum_{i=1}^{n}\left[ \int_{\Omega} u_i(x)g_i(x)d\Omega + \lambda |\partial \Omega_i| \right]. } \end{align*} \]

, where we have

\[ \begin{align*} S &= \left\{ u = (u_1,u_2,\cdots,u_n) \in BV(\Omega) \; : \; u_i(x) = 0,1 \text{ and } \sum_{i=1}^{n} = 1 \right\}, \\[1em] g_i &= || C_i - f ||^2_2. \end{align*} \]

As pointed out by Threshold Dynamics for Networks with Arbitrary Surface Tensions, one can approximate the intersection between two sub-regions of the image by evaluating the result of a difussion process whenever \(\delta t \ll 1\). That is,

\[ \begin{align*} |\partial \Omega_i \cap \partial \Omega_j| & \approx \sqrt{\frac{\pi}{\delta t}}\int_{\Omega}u_iG_{\delta t} * u_jd\Omega. \end{align*} \]

where \(*\) represents convolution and

\[ G_{\delta t}(x) = \frac{1}{4\pi \delta t}exp(-\frac{|x|^2}{4 \delta t}) \]

is the heat kernel. Citing the paper:

The above integral measures the amount of heat that escapes from \(\Omega_j\) to \(\Omega_i\).

Therefore

\[ |\Omega_i| \approx \sum_{j=1, j \neq i}^{n}{\sqrt{\frac{\pi}{\delta t}}\int_{\Omega}{u_iG_{\delta t}*u_jd\Omega}}. \]

Finally, the energy to be minimized is defined as

\[ \mathcal{E}(u,C) = \sum_{i=1}^{n}\int_{\Omega}{\left( u_ig_i + \lambda \sum_{j=1,j \neq i}^{n} \sqrt{\frac{\pi}{\delta t}}u_iG_{\delta t}*u_j\right)d\Omega}. \]

The energy is optimized by splitting the energy in two optimization problems. Given some initial solution \(u^0\) and \(C^0\):

\[ \begin{align*} u^{k+1} &= \text{argmin}_{u \in S} \mathcal{E}^{\delta t} (u,C^k), \\[1em] C^{k+1} &= \text{argmin}_{C \in \mathbb{R}^n} \mathcal{E}^{\delta t}(u^{k+1}, C). \end{align*} \]

The second optimization problem above has solution

\[ C^{k+1}_{i} = \frac{\int_{\Omega}{u_i^{k+1}fd\Omega}}{\int_{\Omega}u_i^{k+1}d\Omega} \]

The first problem is a non-convex energy functional. Surprisingly, its solution is equal to the solution of its convex relaxation counterpart. It is there that enters in scene threshold dynamics.

The idea is to obtain a linear approximation \(\mathcal{L}(u,u^k)\) of the energy and then minimize this linear approximation to obtain \(u^{k+1}\). Following the computations in the paper, we obtain:

\[ \mathcal{L}(u,u^k) = \sum_{i=1}^{n}\int_{\Omega}u_i\phi_i^kd\Omega = \int_{\Omega}\sum_{i=1}^{n}u_i\phi_i^kd\Omega. \]

Given the problem constraints (each pixel can only belong to a single region), it is easy to see that:

\[ u_i^{k+1}(x) = \begin{cases} 1 & \text{if } \phi_i^k(x) = \min_j \phi_j^k(x),\\ 0 & \text{otherwise}. \end{cases} \]

This is equivalent to

\[ \Omega_i^{k+1} = \left\{ x \; : \; \phi_i^k(x) < \min_{j \neq i}\phi_j^k(x)\right\}. \]