ICALAB for Image Processing 
:%s 
The ICALAB Package:
has been designed and developed by:
Andrzej Cichocki, Shunichi Amari, Krzysztof Siwek,
Sergio Cruces, Toshihisa Tanaka, Pando Georgiev, Zbigniew Leonowicz, Tomasz Rutkowski, Seungjin Choi, Adel Belouchrani, Allan Barros, Ruck Thawonmas, Tetsuya Hoya, Wakako Hashimoto, Yasushi Terazono and Tomomi Watanabe in cooperation with other members of the Laboratory for Advanced Brain Signal Processing.
The graphic design, data visualization, user interface, extensive testing and integration of most of the existing algorithms have been implemented in MATLAB^{®} by: Krzysztof Siwek, Andrzej Cichocki, Toshihisa Tanaka and Tomasz Rutkowski.
The current version is 2.0, as of March 20, 2004.
MATLAB^{®} is a registered trademark of The MathWorks, Inc.
A similar package has been developed for Signal Processing.
The comprehensive reference for both packages is the following book:
The important and unique features of our ICALAB toolboxes are preprocessing and postprocessing tools (Fig. 0).
Actual optional PREPROCESSING tools include: Principal Component Analysis (PCA), prewhitening, filtering: High Pass Filtering (HPF), Low Pass Filtering (LPF), Subband filters (Butterworth, Chebyshev, Elliptic) with adjustable order of filters, frequency subbands and the number of subbands).
POSTPROCESSING tools actually includes: Deflation and Reconstruction
("cleaning") of original raw data by removing undesirable components, noise or
artifacts.
Moreover, the ICALAB Toolboxes have flexible and extendable structure with the
possibility to extend the toolbox by the users by adding their own algorithms.
The algorithms can perform not only ICA ;but also Second Order Statistics Blind
Source Separation (BSS) Sparse Component Analysis (SCA), Nonnegative Matrix
Factorization (NMF), Smooth Component Analysis (SmoCA), Factor Analysis (FA)
and any other possible matrix factorization of the form X=HS+N or
Y=WX where H=W^{+} is a mixing matrix or a matrix of basis
vectors.
The ICA/BSS algorithms are pure mathematical formulas, powerful, but rather mechanical procedures:
There is not very much left for the user to do after the
machinery has been optimally implemented. The successful and efficient use of
the ICALAB strongly depends on a priori knowledge, common sense and
appropriate use of the preprocessing and postprocessing tools. In other words,
it is preprocessing of data and postprocessing of models where expertise is
truly needed (see the
book).
On the other hand, the assumed linear mixing models must be valid at least
approximately and original sources signals should have specified statistical
properties.
ICALAB can be useful in the following tasks:
The package contains a collection of algorithms for whitening, robust orthogonalization, ICA, BSS and BSE. The user can easily compare various algorithms for Blind Source Separation (BSS) employing the second order statistics (SOS) and ICA using the higher order statistics (HOS). This package is hence quite versatile and extendable for a user algorithm.
Several benchmarks are included to illustrate the performance of the various algorithms for a selection of synthetic and real world images (see Benchmarks).
Limitation of version 2.0:
The version 2.0 of the package is limited to a maximum of 16 images. In the future, we plan to extend ICALAB for Image Processing for processing up to 256 images, which might be useful for applications such as computer tomography and functional neuroimaging.
NEITHER THE AUTHORS NOR THEIR EMPLOYERS ACCEPT ANY RESPONSIBILITY OR LIABILITY FOR LOSS OR DAMAGE OCCASIONED TO ANY PERSON OR PROPERTY THROUGH USING SOFTWARE, MATERIALS, INSTRUCTIONS, METHODS OR IDEAS CONTAINED HEREIN, OR ACTING OR REFRAINING FROM ACTING AS A RESULT OF SUCH USE. THE AUTHORS EXPRESSLY DISCLAIM ALL IMPLIED WARRANTIES, INCLUDING MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE. THERE WILL BE NO DUTY ON THE AUTHORS TO CORRECT ANY ERRORS OR DEFECTS IN THE SOFTWARE. THIS SOFTWARE AND THE DOCUMENTATIONS ARE THE PROPERTY OF THE AUTHORS AND SHOULD BE ONLY USED FOR SCIENTIFIC AND EDUCATIONAL PURPOSES. ALL SOFTWARE IS PROVIDED FREE AND IT IS NOT SUPPORTED. THE AUTHORS ARE, HOWEVER, HAPPY TO RECEIVE COMMENTS, CRITICISM AND SUGGESTIONS ADDRESSED TO
To start ICALAB for Image Processing type:
icalab
in the MATLAB command window (Note: this package runs on MATLAB 6.0 or higher).
ICALAB for Image Processing was developed under MATLAB version 6.0 and tested under MATLAB versions: 6.0 , 6.1 and 6.5 (Note: Previous versions (i.e. 5.3) may not work properly due to some unsupported functions.)
To load new images for further processing
jpeg
,
png
, bmp
, pcx
, tiff
.
Moreover, you can save all preloaded images in one single
file in MATLAB
format *.mat
(Press on image to enlarge).
The processing speed of many algorithms depends on the size of images and the amount of memory and processor speed of your computer. In order to process large images, you might need to increase the system memory.
mat
file, which contains all images in a specific benchmark.
After loading all images, click on Cancel.
MATLAB
file name.mat
for quick loading and processing of images in the future.The window Select images to load will appear. You can now select the images you want to load. You can discard (delete) a specific image by clicking on it (the color scale of the image will become darker). After selecting images, you can save the set of images into a single MAT file for more convenient loading in future sessions.
The button Save as MAT allows you to choose the destination and file name. Later, you can read such exported image sets directly from Select images to load window.
Before you proceed to the final step, you may choose the loading parameters such as:
Direction of image scanning
Processing area
Finally, you can load so encoded images by clicking on the LOAD button.
You can mix images synthetically in case they are not originally mixed. Leave the option with identity (unit matrix) for real world data.
The option to mix the source signals is applied only for testing and comparing the performance of various algorithms. In ICALAB, eight alternative ways to generate such mixing matrices are preset:
The last option is limited to 15x15 mixing matrix. For this option, click on the EDIT button and the following window will appears. Every element of the mixing matrix may now be edited. After typing in the entries, you will see that both the determinant and condition numbers of the mixing matrix H are automatically updated.
Optionally, the mixing matrix H can be rectangular, you can select a size of the matrix H by clicking on subwindow No. of sensor and choosing the desired number of observations.
You can also add noise to the images before performing ICA or BSS with the following options:
This option can be used e.g., to investigate the robustness of a specific algorithm with respect to additive noise (please try, for example, SOBIRO or SONS algorithm).
Once you have loaded data and optionally mixed them, you can select one of the available algorithms. There is a long list of algorithms which you can apply. You can find detailed descriptions of each algorithm either in the Book or through the online help here.
You can also add your own algorithm(s) to test and compare their performance with other algorithms. Please refer to the example mfiles: user_algk.m to see how ICALAB calls the userdefined algorithms. The user algorithm should return only demixing (separating) matrix W.
user_algk.m
file.Most of the algorithms within the package are supplied with default parameters. Thus, you can start testing the algorithms without the need for adjusting or preselecting any parameter. The default parameters are already tuned to close to approximate optimum values for typical data. Otherwise, you can tune the parameters for most of the algorithms by clicking on the advanced option button ADV. OPTIONS. It is recommended that you use this option and tune the parameters if you are already familiar with the algorithm (see the Book for derivation, properties, and description of the algorithms).
After selecting an algorithm and optionally adjusting its parameters, you can click on the button RUN ALGORITHM. The learning procedure will begin and the algorithm specific messages will appear in the main MATLAB command window. During computation, an additional window displaying the algorithm name will appear. You can stop the learning process by clicking on the INTERRUPT button.
By definition, standard ICA algorithms are not able to estimate statistically dependent sources, that is, when the assumption of independence does not hold. In many cases, however, we may be able to reconstruct the original sources using simple preprocessing techniques and to estimate mixing and separating matrices, even if the sources are not independent.
The ICALAB Toolbox enables blind separation of sources for a wide class of signals that do not satisfy the independence assumption. This can be achieved by applying second order statistics (SOS), exploiting spatiotemporal decorrelation (STD) of sources, or applying linear predictability (LP) and smoothness criteria (see the book) and for some preprocessing, such as: differentiation, high lowpass filtering, sparsification or subband decomposition.
Moreover, each unknown source can be modeled or represented as a sum of narrowband subsignals (components). Provided that for some of the subbands (at least one) subcomponents are mutually independent or temporally decorrelated, suitably designed subband filters can be used in the preprocessing stage to extract mixture of them assuming that these subbands can be identified by a priori knowledge. The standard ICA or BSS algorithms for such transformed (filtered) mixed signals can then be applied. In one of the simplest case, the source signals can be modeled or decomposed into their low and high frequency components. In practice, the high frequency components are often found to be mutually independent. In this case, we can use a High Pass Filter (HPF) to extract high frequency subcomponents and then apply any standard ICA algorithm to such preprocessed sensor (observed) signals.
In the new, updated version ICALAB 2.0 this optional preprocessing has been implemented. To use this option, click on the Preprocessing button at the main ICALAB window.This option is particularly useful for blind separation of dependent or correlated source signals or images, such as faces or natural images, where you will notice significant improvements in the performance of the algorithms. In the preprocessing stage, more sophisticated methods, such as band pass filters or wavelet transforms, can also be applied. Optimal choice of a transformation depends on a specific application and optimal parameters are problem dependent. Experiments are necessary to choose the optimal parameters.
Click on the Preprocessing button in order to perform preprocessing of sensor data. In the first step two window appear (Fig. A), when you can select different preprocessing techniques.
You can choose one from the following options:
Every preprocessing procedure is performed before the ICA algorithm and after pressing of the button OK. Preprocessing options remain active as long as you do not change the mixing matrix, the noise level or one of the preprocessing options.
This option allows to compute differences between values of subsequent pixels and enhance edges in the image.
This option allows comprehensive design of IIR and FIR filters with visualization of parameters. The Fig. C below shows the options, which include:
This option provides a powerful preprocessing method for the ICA/BSS. The subband transform decomposes the input signal into several subbands by applying the corresponding bandpass filters. The figure shows the subband decomposition structure and the frequency responses for the filters. The number of subbands and the specific filter (Butterworth, Chebyshev I/II, Elliptic) can be selected by the user.
Let be m signal mixtures x_{i}(k); (i= 1, ..., n). Let L be the number of subbands. Then, every mixture x_{i}(k) is decomposed into L subsignals x_{i}^{(l)}(k) ; (l = 1, ..., L). It is expected that if we select one or preferably several subsignal(s) from them (including the original mixture x_{i}(k) (denoted as x_{i}^{(0)}(k) in the figure B), based on an appropriate criterion, we can achieve better separation.
You can set parameters listed as follows:
Number of subbands: This parameter corresponds to L in the above figure.
Filter name: To construct a bank of filters, you can choose a filter from Butterworth, Chebyshev I/II and Elliptic.
Order of the filter
Number L of subbands to be selected
Subband selection criterion: It can be chosen the following cost functions : l_{}_{1}norm, l_{}_{p}norm or kurtosis
It is possible to view the spectrum of the data (FFT) by pressing the View FFT button in the preprocessing window.
If you check Display channel signals and select nodes manually you can display a subband signal for each mixture channel and check the value of the cost function. The detail of this option will be explained below.
If this option is selected (see Fig. F), after you click the APPLY button, the window appears where it is possible to choose channels and subbands for further processing:
Subband 3

Subband 5 
The lowerleft figure (in Fig. G) shows the value of the chosen criterion for each subband.
If the mixture i.e. in original, observed data (ROOT in the figure) does not contain higher frequency components, those components are automatically discarded, because the cost function is normalized by the variance of filtered signals, which implies that even if the amplitudes are negligibly small, the value of the cost function can be large. The cost values for those spurious or undesirable subbands are indicated by yellow bars. The selected subband(s) is/ are indicated by red bar(s).
Although the ICALAB can select the subband(s) automatically, the user can preselect any subband manually with possibility to display the filtered signals (Fig. G). The upper plot in figure G presents the subband decomposition tree. The subband images x_{i}^{(l)}(k) are plotted in the right part of Fig. G. To change the target mixture, please choose the mixture number i by selecting the corresponding checkbox below the subband number.
After selecting the suitable subband filtering of sensor signals, press OK in the main preprocessing window and ICA/BSS algorithm starts to run. If you want to reset some options, press CANCEL.
The main references are:
When the program finishes computation, we have three sets of parameters:
After convergence, you can visualize the results by clicking on the PLOT button. The decomposed images are shown in the Independent components window. In addition, you can see the loaded images in the Sources window, which displays the original images. In case you mixed the images by the randomly generated mixing matrix H a I, we recommend that you look at the Sensors (mixtures) window. Essentially, there is no need for this procedure if your original images were already mixed or they represent real (measured or observed) superimposed sequence of images prepared for decomposition into independent components and further processing. You can look at the performance of the selected algorithm in the Global Matrix (Performance Index) window for the artificially mixed sources.
(a) 
(b) 
(c) 
(d) 
If you choose the 3D option, the global matrix G = W H (or G = U Q H) is also plotted in the three dimensional plot (where Q is a prewhitening matrix and U is a rotation matrix in the twostage procedure with orthogonalization or prewhitening. Please see the Book for more details).
Note: Rendering 3D plots of large number of signals might take a long processing time.
The window for Estimated mixing matrix H = inv(W) shows distribution of entries of estimated mixing matrix.
Optionally, you can specify which plots to show. At the bottom of the main ICALAB window, just check/uncheck the plots you want to display. Instead of closing the windows separately, just click on the CLOSE button. Plot windows are closed automatically if you rerun any algorithm.
You can save the separation results for further processing by clicking on the SAVE button. Note that all of the signals will be saved into a single MAT file (with the specified name). In the file, the variables S, X, Y, H, W, G are saved (with the same names as the variables), i.e., S  sources, X  sensors (mixture), Y  independent components or estimated sources, H  mixing matrix, W  demixing matrix, and G  global (mixingdemixing) matrix.
The ICALAB package generally assumes the following generative mixing model
where v(k) is a vector of additive noise,
or in batch mode
where:
V is a matrix representing the additive noise,
and demixing model:
or
or
where
Variable  Description  Dimension of matrix 
S  Sources or independent components  n x N 
X  Observations (sensor signals)  m x N (m >= n) 
Y  Estimated sources or independent components  n x N 
H  Mixing matrix  m x m or n x m 
W  Demixing matrix  m x n or n x n 
G  Global (mixingdemixing) matrix  m x n or n x n 
Remark: In this version some variables are noncapitalized (s, x, y).
Remark: Some algorithms automatically detect the number of sources and reduce the number of outputs to n. The other algorithms assume that the demixing matrix W is an m x m square matrix (see the Book for explanation).
After extracting the independent components or performing blind separation of signals (from the mixture), you can discard some of the components and then project the remaining components back to the sensor domain. This procedure is called deflation and allows you to remove the unnecessary (or undesirable) components that are hidden in the mixture (or overlapped data). In other words, the deflation procedure allows you to extract and remove one or more independent components (or uncorrelated sources) from the mixture x(k). To perform this operation, open the Deflation procedure window by clicking on the DEFLATION button in the main ICALAB window. The DEFLATION button becomes activated only after estimation of the demixing matrix W is completed using any of the builtin or userdefined algorithms.
The deflation procedure is carried out in two steps.
In batch format, the reconstruction procedure can be written as
where
X_{r} = [ x_{r}(1), x_{r}(2), ..., x_{r}(N) ]  reconstructed sensor images,
Y_{r} = [ y_{r}(1), y_{r}(2), ..., y_{r}(N) ]  reduced (selected or filtered) independent components.
In the deflation procedure, you can eliminate the undesirable components by clicking on them in the left plane. In general, you may want to discard more than one component by choosing specific components, representing for example noise, artifacts or undesirable interference. After selecting the appropriate components you can perform deflation by pressing the button Do deflation. An example of removing a component (image no 9) is shown in the Fig. 16.
Note that the deflation procedure may be slow for large data sets.
After completing the deflation procedure, the reconstructed sensor signals can be saved for further processing by clicking on the Save results button. This allows you to process reconstructed sensor signals by applying the same or different ICA/BSS/BSE algorithms. For example, in the first stage you can apply a secondorder statistics BSS algorithm to recover sources with temporal structures. In the second stage, you may use any higherorder statistics ICA algorithm to recover independent components.
In the same way, you can save the reconstructed images obtained by the deflation procedure.
The reconstructed data can be saved in name_XR.mat
file.
It is recommended that one exit from the ICALAB program by clicking on the EXIT button in the main program window below:
Algorithms SANG, NGFICA, NGOnLine, ERICA belong to the family of natural gradient (NG) algorithms and are described in detail in Chapters 6  8.
The algorithms have been originally presented in the following papers:
The wide class of NG algorithms can be represented in a general form as
where is a symmetric positive definite matrix, typically, = _{0} I,
F(y) is a suitably chosen n by n matrix of nonlinearly transformed output signals, typically with entries
or
and
is a diagonal positive definite matrix, typically with entries
or
SANG is a batch algorithm which can be represented in the following form:
where the entries of the diagonal matrix D satisfy the nonholonomic constraints, i.e.
and the learning rate matrix is a diagonal matrix with entries
Main references:
The NG  FlexICA (Natural Gradient Flexible ICA) algorithm has been developed by S. Choi, A. Cichocki and S. Amari (with MATLAB implementation by Seungjin Choi).
The main references are:
The algorithm is described in detail in Section 6.4.
TICA Algorithm
Thin algorithm for Independent Component Analysis (TICA) has been developed by Sergio Cruces and Andrzej Cichocki [1].
The TICA algorithm is able to extract simultaneously arbitrary number of components specified by the user. The algorithm is based on criteria that jointly perform maximization of several cumulants of the outputs and/or second order time delay covariance matrices. This employed contrast function combines the robustness of the joint approximate diagonalization techniques with the flexibility of the methods for blind signal extraction. Its maximization leads to hierarchical and simultaneous ICA extraction algorithms which are respectively based on the thin QR and thin SVD factorizations. The TICA algorithm can be regarded as hierarchical/simultaneous extensions of the fast fixed point algorithms [2] and has also close links to [3] and [4,5].
The main references are:
ERICA Algorithm  Equivariant Robust Independent Component Analysis algorithm (which is asymptotically equivariant in the presence of Gaussian noise) has been developed by Sergio Cruces, Luis Castedo and Andrzej Cichocki.
This algorithm separates the signals from m mixture of n sources (with non zero kurtosis) in the presence of Gaussian noise.
This algorithm is a quasiNewton iteration that will converge to a saddle point with locally isotropic convergence, regardless of the distributions of sources. The use of prewhitening is not necessary for this algorithm to converge.
The main references are:
SIMBEC algorithm (Simultaneous Blind signal Extraction using Cumulants) was implemented by Sergio Cruces, Andrzej Cichocki, and Shunichi Amari.
It performs simultaneous blind signal extraction of an arbitrary group of sources from a rather large number of observations. Amari proposed in [1] a gradient algorithm that optimizes the ML criteria on the Stiefel manifold and solves the problem when the approximate (or hypothetical) densities of the desired signals are a priori known. This algorithm [24] extends this result to other contrast functions that do not require explicit knowledge of the source densities. The necessary and sufficient local stability conditions of the algorithm are exploited to obtain fast convergence.
The main references are:
The algorithm can extract an arbitrary group of sources specified by the user.
The JADE algorithm has been originally developed and implemented by JeanFrancois Cardoso and Antoine Souloumiac.
The main references are:
The related references are:
In the ICALAB package, we have used a modified version of this algorithm with reduced number of eigen matrices as described in Chapter 4. We have also optimized some numerical procedures in MATLAB to speed up the algorithm. This algorithm is free of any adjustable parameters (there is no parameter tuning).
The JADEopt has been modified and implemented by Y. Terazono.
The Fixed Point or Fast ICA algorithm has been originally developed by Aapo Hyvärinen and Erkki Oja:
In the ICALAB package, we have implemented the FPICA sequential blind extraction algorithm which extracts independent nonGaussian distributed sources one by one as described in Chapter 5, Section 5.2.4.
Pearson system is an HOS ICA algorithm originally developed by Juha Karvanen, Jan Eriksson and Visa Koivunen:
The Pearson system is employed for modeling a wide class of both symmetric and asymmetric source distributions.
In fact, the algorithm is similar to the natural gradient (NG) and its optimized parallel version Fast ICA algorithm for maximum likelihood estimation. The main difference is the selection of an activation function f(y) (score function) on the basis of estimated online second, third and fourth order moments (see chapter 6, section 6.4).
The BSS SVD algorithm is similar to the AMUSE algorithm (see the Book). However, it uses a different prewhitening procedure and singular value decomposition (SVD) instead of the symmetric EVD for a single timedelayed covariance matrix. The algorithm is very fast but unfortunately, sensitive to additive noise.
See Chapter 4 and the related references:
The SOBIRO (Robust Second Order Blind Identification with Robust Orthogonalization) algorithm is described in detail in Chapter 4 (Section 4.4).
The SOBIRO algorithm employs the robust orthogonalization preprocessing as described in Section 4.3.1.
The related references are:
SONS (Second Order Nonstationary Source Separation) algorithm is described in detail in Chapter 4 (Section 4.4.1).
The algorithm allows you to perform both ICA (for arbitrary distributed nonstationary sources with temporal structures) and BSS (for colored sources with different spectra), depending on the time delays, the number of time windows, and the size of each window.
In the Advanced option, you can select these parameters to optimize the algorithm for a specific problem.
If you want to obtain smooth signals the subwindow should have 500 samples or more. In order to achieve sparse independent signals choose a short subwindow in the range 1050 samples.
The SONS algorithm has been developed by S. Choi and A. Cichocki and presented in the following:
Related references are:
You can integrate your own algorithms or any ICA/BSS, BSE algorithm available
in MATLAB within the ICALAB package, by simply inserting the code(s) in each of:
user_alg1.m
 user_alg10.m
files.
Note: the algorithm must be coded in such a way that it returns only one single variable: demixing matrix W.
If the algorithm estimates the mixing matrix H, you may use W = pinv(H).