Problems with the FastICA filter
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.
FastICA algorithm used is the one described in . 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.
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 the
UpdateOutputInformation()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 in
FastICAfilter uses (in the
FastICAInternalOptimizerVectorImageFilterwhich 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
FastICAInternalOptimizerVectorImageFilteri.e. computing quantities (
m_Beta) based on the accumulators is done in the
AfterThreadedGenerateData()is called every time
GenerateData()is called, i.e. once for each piece of stream. So instead of computing
m_Betausing 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  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).
PCAis computed with the option
useNormalizationto 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 the
PCAFiltercomputes 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 set
UseVarianceForNormalizationto false the data is correctly whitened.
The non linear function
fwhich is stored as a function pointer can be set in the filter by the user. By default it is
std::tanh. The problem is that we also need its derivative, which is hard-coded as
1-f^2, so the filter will only work for non linear function that verifies
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!
 A. Hyvarinen. Fast and robust fixed-point algorithms for independent component analysis. Neural Networks, IEEE Transactions on, 10(3):626–634, 1999