Problems with the FastICA filter
Description
There is a FastICAFilter
in the DimensionalityReduction
module. The filter is not tested (there is actually a test that is compiled but not executed). The DimensionalityReduction
application uses this filter, but it also is untested.
In MR !336 (merged) we tried to add this test to CTest
, but it produces different results depending on the number of threads used, which makes Dashboard testing hard. One solution would be to force the number of threads to 1 in the test, however, I don't think the result of FastICA should depend on the number of thread used and that this filter does not do what it is supposed to do.
Note that the result of the test also changes if we use float instead of double as pixel type.
The FastICA
algorithm used is the one described in [1]. It is implemented as an ImageToImageFilter
and uses the same framework as the other Dimensionality Reduction filters (pca, mnf ...) i.e. :

They can work both in "Forward" and "Backward" mode, forward meaning applying the reduction transformation and backward applying its inverse. This is a template parameter.

The
GenerateData()
method of the filter only apply the transformation, i.e. the result of the preprocessing PCA (whitening) followed by the result of the ICA. this transformation can be supplied to the filter via Setters, but if it is not, it will be computed in theUpdateOutputInformation()
method using the FastICA algorithm (the filter also computes the PCA whitening transform if needed). This means that most of the processing of this filter is done inUpdateOutputInformation()
. 
The
FastICA
filter uses (in theUpdateOutputInformation()
method!) aFastICAInternalOptimizerVectorImageFilter
which is a "kind of" Persistent filter (meaning it is not persistent but maybe it should be), this filter takes as input the output of the PCAFilter, and computes output pixels but also accumulators and is plugged to a StreamingStatisticsVectorImageFilter. This pipeline is executed at each iteration of the FastICA algorithm.
I found many problems in this filter, but nothing that really explains why the result depends on the number of thread :

The "persistent computation" of the
FastICAInternalOptimizerVectorImageFilter
i.e. computing quantities (m_Den
andm_Beta
) based on the accumulators is done in theAfterThreadedGenerateData()
method.AfterThreadedGenerateData()
is called every timeGenerateData()
is called, i.e. once for each piece of stream. So instead of computingm_Den
andm_Beta
using the whole image, it is computed using only the last stream of the input image. This means that the filter does not support streaming ... Note that it supports threading, and this does not explain the bug. 
I think the formula used for the transformation matrix update step is wrong. Currently it is :
W(band, bd) = m_Mu * ( estimator>GetMean()[bd]
 optimizer>GetBeta() * W(band, bd) / optimizer>GetDen() );
but according to [1] it should be :
W(band, bd) = m_Mu * ( estimator>GetMean()[bd]
 optimizer>GetBeta() * W(band, bd) ) / optimizer>GetDen() ;
(note the parenthesis)

At the first step of the FastICA algorithm, we use the input Image as input of the filter instead of the PCA output. This is a problem because FastICA expects whitened input (uncorrelated, zero mean and unit variance).

The
PCA
is computed with the optionuseNormalization
to true, The input image is normalized in mean and variance before the PCA. While removing the mean is required, I'm not sure about the variance normalization before the PCA for ICA in general and it doesn't work in OTB case. I'm not sure what thePCAFilter
computes when the variance normalization is on, but it is not a whitening transformation (i.e. D^1/2 * E, where D is the eigenvalue matrix and E the (transposed) eigenvector matrix). If we setUseVarianceForNormalization
to false the data is correctly whitened. 
The non linear function
f
which is stored as a function pointer can be set in the filter by the user. By default it isstd::tanh
. The problem is that we also need its derivative, which is hardcoded as1f^2
, so the filter will only work for non linear function that verifiesf'=1f^2
, i.etanh
. 
The dimensionality reduction is done at the end of the algorithm, by selecting the five first independent components. But these components are in no specific order, so this is equivalent to taking five components at random (ICA is not really a dimensionality reduction algorithm ...). Instead the reduction should be done in the PCA step. This would also reduce computation time!
[1] A. Hyvarinen. Fast and robust fixedpoint algorithms for independent component analysis. Neural Networks, IEEE Transactions on, 10(3):626–634, 1999