## 2 Methods
### \"K-Means + Gaussian Mixture Model\" Clustering
The GMM can fit data through arbitrary ellipsoidal clusters, which are more suitable for dealing with the regional division of complex biological tissues. However, the GMM is very sensitive to the setting of initial parameters, and the number of clusters needs to be set in advance. The clustering result of K-means can provide it with the initial values of these parameters, avoiding the influence of the randomness of the initial parameters on the GMM. The GMM is a parametric probabilistic model. This model is given as follows:  
\[f(p_{i})=\sum_{k=1}^{K}a_{k}N_{k}(p_{i},\mu_{k},\Sigma_{k}). \tag{9}\]  
where \(N_{k}(p_{i},\mu_{k},\Sigma_{k})\) stands for the probability density function of the kth Gaussian distribution. \(a_{k}\) is the weight coefficient, which represents the proportion of each type of tissue. \(\mu_{k}\) is the mean value, which corresponds to the central position of the distribution of each tissue. \(\Sigma_{k}\) is the variance, which describes the distribution of data around the mean value. When the GMM is applied to fit the photoacoustic pressure data, the data need to be allocated to \(K\) Gaussian distributions, and each Gaussian distribution corresponds to the distribution information of a class of tissues. The initial parameters are set based on the clustering result of K-means. The cluster center point of the \(k\)th cluster can be set directly to the initial mean value of the \(k\)th Gaussian distribution. The initial weight coefficient of the \(k\)th Gaussian distribution is obtained according to the proportion of the number of data \(M\) in the \(k\)th cluster to the total data. The initial mean and the number of data \(M\) are  
Figure 2: (**a**) Four sets of randomly generated data. (**b**) Cluster deviation map for different K-values. brought into Equation (10). The initial variance parameter of the \(k\)th Gaussian distribution can be obtained. \[\Sigma_{0,k}=\frac{\sum_{i}^{M}\left(p_{i}-\mu_{0,k}\right)\left(p_{i}-\mu_{0,k} \right)^{T}}{M}. \tag{10}\]  
After the initial parameters of the GMM are obtained, the expectation-maximization (EM) algorithm is used to iteratively optimize these parameters to obtain a set of optimal model parameters. The posteriori probability is estimated by the current model parameters in E-step of the EM algorithm, namely, the probability that the data \(p_{i}\) belong to the kth Gaussian distribution as follows [20]:  
\[\gamma_{i,k}^{new}=\frac{a_{k}N_{k}(p_{i},\mu_{k},\Sigma_{k})}{\sum_{k=1}^{K}a _{k}N_{k}(p_{i},\mu_{k},\Sigma_{k})}. \tag{11}\]  
The updated model parameters are calculated through the M-step of the EM algorithm as follows:  
\[\mu_{k}^{new}=\frac{\sum_{i}^{N}\left(\gamma_{i,k}^{new}p_{i}\right)}{\sum_{i }^{N}\gamma_{i,k}^{new}}, \tag{12}\]  
\[\Sigma_{k}^{new}=\frac{\sum_{i}^{N}\gamma_{i,k}^{new}(p_{i}-\mu_{k}^{new})(p_{ i}-\mu_{k}^{new})^{T}}{\sum_{i}^{N}\gamma_{i,k}^{new}}, \tag{13}\]  
\[a_{k}^{new}=\frac{\sum_{i}^{N}\gamma_{i,k}^{new}}{N}. \tag{14}\]  
In Equation (14), \(N\) represents the total number of photoacoustic pressure data. From Equations (12)-(14), a set of more accurate model parameters can be obtained. The GMM described by the set of model parameters is closer to the real distribution of biological tissues. These parameters are used to calculate the log-likelihood function of the GMM as follows:  
\[\ln f(p_{i})=\sum_{i=1}^{N}\ln\Bigl{\{}\sum_{k=1}^{K}a_{k}^{new}N_{k}^{new}(p_ {i},\mu_{k}^{new},\Sigma_{k}^{new})\Bigr{\}}. \tag{15}\]  
If the log-likelihood function does not converge, the model parameters need to return to E-step again for iterative updates. When the log-likelihood function of the model reaches the convergence state, the fitting effect of the GMM to the data also reaches its best. The current model parameter is the optimal model parameter. Each photoacoustic pressure datum is assigned to the Gaussian distribution with the largest probability according to the posteriori probability, and the clustering result can be obtained.