# Retrieval of Similar Evolution Patterns from Satellite Image Time Series +
## 2 Materials and Methods
### Determining the Optimal Threshold Using Expectation-Maximization
_Maximization step (M-step)_. Update the model parameters such that the log-likelihood approaches its maximum, i.e., the convergence of the EM algorithm guarantees that the log-likelihood value is increased with each iteration [19]. It is usually convenient to introduce a mapping between the new model parameters \(\theta^{(t+1)}\) and the previous ones \(\theta^{(t)}\). According to [18] and following the notations in [19], the means, squared standard deviations, and mixture probabilities for each class \(k\in\{s,n\}\) can be re-estimated using the following set of relations: \[\pi_{k}^{(t+1)} = \frac{\sum\limits_{x}\zeta_{k}^{(t)}(x)}{N}\] (9) \[\mu_{k}^{(t+1)} = \frac{\sum\limits_{x}\zeta_{k}^{(t)}(x)\cdot x}{\sum\limits_{x} \zeta_{k}^{(t)}(x)}\] (10) \[\sigma_{k}^{(t+1)} = \left(\frac{\sum\limits_{x}\zeta_{k}^{(t)}(x)\cdot\left(x-\mu_{k }^{(t)}\right)^{2}}{\sum\limits_{x}\zeta_{k}^{(t)}(x)}\right)^{\frac{1}{2}}\] (11) where \(N=W\times H\) is the total number of pixels in the distance image (i.e., the number of spatio-temporal evolutions) and \(\zeta_{k}^{(t)}(x)\) values are derived as \[\zeta_{k}^{(t)}(x)=\frac{\pi_{k}^{(t)}\mathcal{N}(x|\mu_{k}^{(t)},\sigma_{k}^ {(t)})}{\sum\
olimits_{x}\pi_{k^{\prime}}^{(t)}\mathcal{N}(x|\mu_{k^{\prime}} ^{(t)},\sigma_{k^{\prime}}^{(t)})}.\] (12)  
The initialization of the EM algorithm is done by separating the similarity scores into two subsets, \(\mathcal{S}_{n}^{(0)}\) and \(\mathcal{S}_{c}^{(0)}\). Any unsupervised clustering method can be applied, but we chose the \(K\)-means (\(K=2\)) clustering method due to its fast response time when the number of clusters is small. However, being an unsupervised clustering method, a constraint must still be fulfilled. Namely, the cluster of similar evolutions must be characterized by DTW similarity scores that are smaller than the DTW similarity scores obtained for the temporal sequences pertaining to the cluster of non-similar evolutions with respect to the query. The initial values for the prior probabilities, means, and squared standard deviations can be easily computed as statistics over the two subsets, \(\mathcal{S}_{n}^{(0)}\) and \(\mathcal{S}_{c}^{(0)}\). Compared to the initialization proposed in [1], the K-means initialization speeds up the convergence of the EM algorithm [19], and thus, the speed of the whole algorithm. After the estimation of the \(\theta\) model parameters, the optimal threshold \(T_{o}\) that separates the two classes (i.e., _similar_ and _non-similar_) is determined from the equality  
\[\pi_{s}\mathcal{N}(T_{o}|\mu_{s},\sigma_{s})=\pi_{n}\mathcal{N}(T_{o}|\mu_{n},\sigma_{n}) \tag{13}\]  
that naturally follows from the Maximum Likelihood rule  
\[p_{\theta}(x)=\pi_{s}\mathcal{N}(x|\mu_{s},\sigma_{s})+\pi_{n}\mathcal{N}(x| \mu_{n},\sigma_{n}) \tag{6}\]0Taking the logarithm of both sides of Equation (13) yields a quadratic equation in \(T_{o}\)  
\[p_{\theta}(x)=\pi_{s}\mathcal{N}(x|\mu_{s},\sigma_{s})+\pi_{n}\mathcal{N}(x| \mu_{n},\sigma_{n}) \tag{6}\]1  
that has two possible solutions. However, only one of these solutions lies between \(\mu_{s}\) and \(\mu_{n}\) and this solution represents the optimal threshold value. After determining the optimal threshold value \(T_{o}\), the DTW distance image **DI** is transformed into a binary image **S**, delimiting the locations of evolutions that are similar to the given query:  
\[\textbf{S}(w,h)= \begin{cases}1,&\textbf{DI}(w,h)\leq T_{o}\\ 0,&\text{otherwise}\end{cases}  \tag{16}\]  
where \(1\leq w\leq W\) and \(1\leq h\leq H\), \(W\) and \(H\) being the width and height of the images in SITS.