# Robust Registration of Medical Images in the Presence of Spatially-Varying Noise
## 2 Materials and Methods  
### Empirical Mode Decomposition  
Decomposing nonlinear and non-stationary signals into their intrinsic modes have been a challenging task in signal processing. Most of the previously proposed decomposition is based on time-frequency transformations such as Short-Time Fourier Transform [34] and wavelet transform [35] that expands the signal into a set of basis functions. Unlike these transformations, algorithms based on empirical mode decomposition (EMD) [29; 36] are able to decompose the linear or nonlinear signal into a set of functions defined by the signal itself [37]. EMD consists of a set of spatial and temporal processes that decomposes a signal into Intrinsic Mode Functions (IMFs). Each IMF is an oscillating mode of the signal. In contrast with harmonic modes, oscillating modes have time-variant amplitude and frequency. EMD has been widely used in applications of signal decomposition such as prediction [38], classification [39] and denoising [40].  
In particular, bidimensional empirical mode decomposition [41] is our focus in this paper since it deals with two-dimensional signals such as images. Each IMF in the bidimensional empirical mode decomposition could be extracted from the image using a sifting process. Let \(\mathbb{I}\) denote the image, \(\mathbb{IMF}_{i}\) and \(\mathbb{RES}_{i}\) denote the \(i\)th IMF and residual, respectively. Here \(i\) represents the level or scale of the decomposition. For simplicity, let us assume that the residual at level 0 is the original image. That is:  
\[\mathbb{RES}_{0}=\mathbb{I}.\]  
To find the IMF in level \(i\), we iterate through the following steps: First, identify all the local minimum and maximum of the \(\mathbb{RES}_{i-1}\). Then, generate the bounding maximum and minimum envelopes, \(\mathbb{E}^{max}\) and \(\mathbb{E}^{min}\), using an interpolation or a surface fitting technique. The IMF in this iteration is equal to \(\mathbb{RES}_{i-1}\) subtracted by the mean of maximum and minimum envelope:  
\[\mathbb{IMF}_{i}=\mathbb{RES}_{i-1}-\frac{\mathbb{E}^{max}+\mathbb{E}^{min}}{2}\]  
We repeat this process until a convergence condition is met. At this point, the residual is defined as:  
\[\mathbb{RES}_{i}=\mathbb{RES}_{i-1}-\mathbb{IMF}_{i}\]  
The first IMF has the finest and the most detailed patterns while the rest of the IMFs have patterns with coarser content. The extracted IMFs have a number of crucial characteristics. First, the number of zero-crossing points is equal to or one less than the number of extremes in each IMF. Second, each IMF has one minimum with zero value. Finally, the mean value of the bounding maximum and minimum envelopes is equal to zero for each IMF.  
It is known that the EMD has limited capacity in processing narrow-band signals [42]. These signals have concentrated energy in a narrow frequency band, which causes the decomposition to fail in extracting the intrinsic modes. Thus, the accuracy and efficiency of the EMD could be limited in this situation. However, this is not an issue in our study for MRI and retina images. These images do not have a narrow-band frequency component and therefore, are suited to be analyzed through EMD.