Blind Source Separation for MRI Images

The MRI images at the focus of this project are of the brain's axial slices. The brain mostly contains two types of substances: water and fat.

The MRI images at the focus of this project are of the brain’s axial slices. The brain mostly contains two types of substances: water and fat. An examination of the properties of MRI images reveals that the images of the brain are a linear combination of the images of the brain’s components: water and fat tissues ( [1], [2]). Thus, we regard the MRI images as the observed mixtures, and the fat and water tissues pictures as the sources.

The project’s goal is to recover the images of the sources from a given set of MRI images. The sources would give us an indication of any abnormality in the tissues’ structure in the brain, and could help in early diagnosis of cancer tumors. Since the MRI images are a linear combination of the sources, the problem of separating the observed mixtures into their sources is a well-known problem from the signal processing field, called Blind Source Separation (BSS).

The linear equation that describes the process is:


Where X1, X2, S1, S2 are in row stack image representation.
The noise component is added to account for other types of subdominant tissues such as skull tissue. However, the noise component is typically small in MRI images, and it will be ignored to simplify the analysis. Further references on noise impact and practical solutions can be found in the project report.


The key to solving the equation is to estimate 2and then calculate 3. In order to test the accuracy of the BSS method and examine our accuracy in finding A, we have created several test-cases of linear combinations of images and implemented our method on the mixtures.


BSS method ([3])
We will use BSS in order to estimate the mixing matrix A, given X. This method requires sparse sources, so we used the Wavelet decomposition in order to transform the sources into their sparse representations. An image is composed of sparse sources if most pixels in the image contain one dominant source. In MRI pictures this translates into only one dominant tissue substance in each pixel of the axial cut picture.
As mentioned, the Wavelet transform ([4], [5]) turns S into sparse sources. Due to the linearity of the transform, performing the transform on X is equal to performing it on S. In this fashion, we get the following expressions after the transform:


The wavelet transform creates sparse sources due to the fact that it decomposes the image both in the time and frequency domains. Therefore, the transform extracts local variations in the source image, such as contours. The contours vary in different source images because they are the main difference between two images. Therefore, representing an image through its contours creates a sparse representation in the mixtures when we linearly combine the sources.
If the sources are approximately sparse, a graph of X1_w vs. X2_ will produce two lines with differing slopes:






Thus, by finding the slopes of the two lines, we can find A – though the columns of the matrix may be swapped, and the columns may be multiplied by a constant.


Finding the slopes from the graph
The task of devising an algorithm to find the slopes of the two lines is a central requirement of the project, and relates to many problems in the field of image processing – such as the peak picking problem. This process is a simple task for the human observer, since for the human visual system the dominant lines are usually obvious. However, for a computer, the task is much more complicated and challenging, yet it offers the advantages of accuracy, time efficiency, and protection from human errors.
In order to find the slopes, we have created a polar histogram of the angle from the origin of all the points in the graph above.



Addressing this problem, we chose to take two different steps:
1. Noise reduction: around the origin we have many points which function as a large noise source, because a small error at a point near the origin usually translates into a meaningful error in the points’ slope. We have created a function that finds a radius of a large density of points around the origin, and extracts all the points in that radius from the histogram. This tool has significantly improved the accuracy of our solutions.
2. Peak picking: we found the slopes from the histogram with a combination of three filters (using [6]):
– Smoothing the graph with a low pass Gaussian filter to avoid a decrease of the multiple peaks phenomena due to noise.
– Selecting all the suspected peaks, where a peak is defined as a value that is larger than its neighboring values.
– Checking if a peak is a local maximum: by checking that the average of the values that surround the peak is lower than the value of the peak. We used Taurus distance to make a circular action in the averaging near (-90) and 90 (since they are in fact the same slope). Next, we chose the two maximum values of the remaining peaks.
An example of this process can be seen in the following graphs. Notice that the graphs use a linear scale, and that the peaks in the bottom two graphs are emphasized (for visualization).


The result of the two steps in the example above is:



We can observe that the expected slope and the automatically detected slopes are the same!


Working with three observed mixtures

Up until now, we had only two mixtures provided. Now let us assume that we are given three mixtures. In this case, we are provided more information, and we can reconstruct A more efficiently.


In order to find S1 and S2, we must estimate A. From A, we can find the best solution for S in MS manner by performing a least square optimization ([7]):


For each pixel in the space x1,x2,x3:


Therefore, in this space, the mixture’s pixels are a linear combination of two constant vectors, so they would all be on a plane in space containing the origin. After we perform a Wavelet transform on x1, x2, x3 we will get, like before, sparse sources, and the sparse graph will have two dominant lines, such as in the following example:


As in the previous case, by finding the directions of the two dominant lines, we will find A, and then S.


In our project we have demonstrated a couple of important points:

– In order to use our method of finding slopes with two mixtures, we have projected all the points to the common plane, and then used our 2D algorithm. In this stage, we had to define a new 2D axis system.
– From the 2D slopes we found the direction in space using the new axis system.
– We have proved that we can perform the Wavelet transform AFTER the projection on the 2D points, and obtain similar results. In this way we can perform only two transforms, and refrain from having to perform a complex, expensive and redundant transformation.
– In case of noise, we will perform a plane fitting, project all the points to that plane, and determine a solution using the same algorithm as the linear case.


The following is an example of the method on two MRI mixtures:


The source images are:


The sources are the images with the white background. We have also shown the negatives from the reconstructed images, because we found the images were multiplied by a constant, which could be negative. In order to reduce noise, we focused on a rectangular area that did not contain any skull tissue (subdominant source).

An example of image mixtures is show below:


Their associated reconstructed images are:


As one can notice, the reconstruction is accurate.

Further 2D and 3D Examples, mathematical explanations and proofs, implementation details, and a theoretical background can be found in the project report.
We would like to thank our project supervisor Ehud Orian for the guidance, patience, inspiration and efforts. We are also grateful to the Ollendorf Minerva Center Fund for supporting this project.


[2] Basics of Magnetic Resonance Imaging (Chapter 1)
[3] Multiscale Framework for Blind Separation of Linearly Mixed Signals by Pavel Kisilev,

Michael Zibulevsky, Yehoshua Y. Zeevi
[4] “”A Really Friendly Guide for Wavelets”” by C.Valence
[5] Wavelets lab (from visl lab)
[6] Complex Domain Onset Detection For Musical Signals by Chris Duxbury, Juan Pablo Bello,

Mike Davies and Mark Sandler