Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I have a stack of images (usually ca 100) of the same sample. The images have intrinsic variation of the sample, which is my signal, and a lot of statistical noise. I did a principal components analysis (PCA) on the whole stack and found that components 2-5 are just random noise, whereas the rest is fine. How can I produce a new stack of images where the noise components are filtered out?

EDIT:

I am sorry I was not as active as you yesterday. I must admit I am bit overwhelm by the depth and yet simplicity of your answers. It is hard for me to choose one, since all of them work great and give what I actually wanted.

I feel that I need to elaborate a bit more the problem I am working on. Unfortunately, my supervisor does not allow me to upload nay data before we have published the final results, so I have to work in abstract terms. We have an atomic cloud cooled to a temperature of 10 µK. Due to inter-atomic and laser interaction, the atomic cloud (all of the atoms a whole) is excited and starts to oscillates in different vibrational modes. This dynamic behavior is of great interest to us, since it provides and insight to the inter-atomic physics.

The problem is that most of the relevant variations are obscured by noise due to the imaging process. The noise usually is greatly suppressed if you take two images one with Noise+Signal and Noise only and then subtract them. However, this does not work if the noise in the two images is not correlated, which sadly is our case. Therefore, we decided to use PCA, because there you can clearly see the oscillation modes and filter everything that is crap. If you are interested in using PCA to visualize dynamics, you can have a look at this paper by different group:

http://iopscience.iop.org/article/10.1088/1367-2630/16/12/122001

I deeply thank everybody who contributed.

share|improve this question
3  
PCA is basically an SVD process. You take the SVD and keep only the singular vectors that correspond to the desired singular values. So in Mathematica, you should look at SingularValueDecomposition[k,m] – bill s May 6 at 12:39
2  
Also this post should have a concrete example as input, along with whatever code you have tried thus far. – Daniel Lichtblau May 6 at 22:20
2  
I'm curious as to why components 6-9 would be useful, from an SVD won't they be explaining less of the variance than components 2-5? – blochwave May 6 at 23:10
2  
I have some sympathy for the situation you are in, constrained by prepublication requirements, research supervisor, etc. All the same, it seems quite difficult to get agreement by very experienced users as to what exactly comprises a viable example or methodology. I'll vote to reopen if maybe that much can be sorted out. – Daniel Lichtblau May 8 at 14:51
2  
I mean a data set that you are allowed to use (obviously not the one your supervisor is disallowing), one that has similar noise characteristics. That way people can try the PCA approch on a set that is known to share the salient features of your problem. Without that there has been a lot of guesswork and dispute as to what constitutes a representative example or code to handle such. – Daniel Lichtblau May 8 at 21:20
up vote 5 down vote accepted

This answer compares two dimension reduction techniques SVD and Non-Negative Matrix Factorization (NNMF) over a set of images with two different classes of signals (two digits below) produced by different generators and overlaid with different types of noise.

Note that question states that the images have one class of signals:

I have a stack of images (usually ca 100) of the same sample.

PCA/SVD produces somewhat good results, but NNMF often provides great results. The factors of NNMF allow interpretation of the basis vectors of the dimension reduction, PCA in general does not. This can be seen in the example below.

Data

The data set-up is explained in more detail in my previous answer. Here is the code used in this one:

MNISTdigits = ExampleData[{"MachineLearning", "MNIST"}, "TestData"];
{testImages, testImageLabels} = 
  Transpose[List @@@ RandomSample[Cases[MNISTdigits, HoldPattern[(im_ -> 0 | 4)]], 100]];
testImages

enter image description here

See the breakdown of signal classes:

Tally[testImageLabels]
(* {{4, 48}, {0, 52}} *)

Verify the images have the same sizes:

Tally[ImageDimensions /@ testImages]
dims = %[[1, 1]]

Add different kinds of noise to the images:

noisyTestImages6 = 
  Table[ImageEffect[
    testImages[[i]], 
    {RandomChoice[{"GaussianNoise", "PoissonNoise", "SaltPepperNoise"}], RandomReal[1]}], {i, Length[testImages]}];
RandomSample[Thread[{testImages, noisyTestImages6}], 15]

enter image description here

Since the important values of the signals are 0 or close to 0 we negate the noisy images:

negNoisyTestImages6 = ImageAdjust@*ColorNegate /@ noisyTestImages6

enter image description here

Linear vector space representation

We unfold the images into vectors and stack them into a matrix:

noisyTestImagesMat = (Flatten@*ImageData) /@ negNoisyTestImages6;
Dimensions[noisyTestImagesMat]

(* {100, 784} *)

Here is centralized version of the matrix to be used with PCA/SVD:

cNoisyTestImagesMat = 
  Map[# - Mean[noisyTestImagesMat] &, noisyTestImagesMat];

(With NNMF we want to use the non-centralized one.)

Here confirm the values in those matrices:

Grid[{Histogram[Flatten[#], 40, PlotRange -> All, 
     ImageSize -> Medium] & /@ {noisyTestImagesMat, 
    cNoisyTestImagesMat}}]

enter image description here

SVD dimension reduction

For more details see the previous answers.

{U, S, V} = SingularValueDecomposition[cNoisyTestImagesMat, 100];
ListPlot[Diagonal[S], PlotRange -> All, PlotTheme -> "Detailed"]
dS = S;
Do[dS[[i, i]] = 0, {i, Range[10, Length[S], 1]}]
newMat = U.dS.Transpose[V];
denoisedImages = 
  Map[Image[Partition[# + Mean[noisyTestImagesMat], dims[[2]]]] &, newMat];

enter image description here

Here are how the new basis vectors look like:

Take[#, 50] &@
 MapThread[{#1, Norm[#2], 
    ImageAdjust@Image[Partition[Rescale[#3], dims[[1]]]]} &, {Range[
    Dimensions[V][[2]]], Diagonal[S], Transpose[V]}]

enter image description here

Basically, we cannot tell much from these SVD basis vectors images.

Load packages

Here we load the packages for the what is computed next:

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/NonNegativeMatrixFactorization.m"]

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/OutlierIdentifiers.m"]

NNMF

This command factorizes the image matrix into the product $W H$ :

{W, H} = GDCLS[noisyTestImagesMat, 20, "MaxSteps" -> 200];
{W, H} = LeftNormalizeMatrixProduct[W, H];

Dimensions[W]
Dimensions[H]

(* {100, 20} *)
(* {20, 784} *)

The rows of $H$ are interpreted as new basis vectors and the rows of $W$ are the coordinates of the images in that new basis. Some appropriate normalization was also done for that interpretation. Note that we are using the non-normalized image matrix.

Let us see the norms of $H$ and mark the top outliers:

norms = Norm /@ H;
ListPlot[norms, PlotRange -> All, PlotLabel -> "Norms of H rows", 
  PlotTheme -> "Detailed"] // 
 ColorPlotOutliers[TopOutliers@*HampelIdentifierParameters]
OutlierPosition[norms, TopOutliers@*HampelIdentifierParameters]
OutlierPosition[norms, TopOutliers@*SPLUSQuartileIdentifierParameters]

enter image description here

Here is the interpretation of the new basis vectors (the outliers are marked in red):

MapIndexed[{#2[[1]], Norm[#], Image[Partition[#, dims[[1]]]]} &, H] /. (# -> Style[#, Red] & /@ 
   OutlierPosition[norms, TopOutliers@*HampelIdentifierParameters])

enter image description here

Using only the outliers of $H$ let us reconstruct the image matrix and the de-noised images:

pos = {1, 6, 10}
dHN = Total[norms]/Total[norms[[pos]]]*
   DiagonalMatrix[
    ReplacePart[ConstantArray[0, Length[norms]], 
     Map[List, pos] -> 1]];
newMatNNMF = W.dHN.H;
denoisedImagesNNMF = 
  Map[Image[Partition[#, dims[[2]]]] &, newMatNNMF];

Comparison

At this point we can plot all images together for comparison:

imgRows = 
  Transpose[{testImages, noisyTestImages6, 
    ImageAdjust@*ColorNegate /@ denoisedImages, 
    ImageAdjust@*ColorNegate /@ denoisedImagesNNMF}];
With[{ncol = 5}, 
 Grid[Prepend[Partition[imgRows, ncol], 
   Style[#, Blue, FontFamily -> "Times"] & /@ 
    Table[{"original", "noised", "SVD", "NNMF"}, ncol]]]]

enter image description here

We can see that NNMF produces cleaner images. This can be also observed/confirmed using threshold binarization -- the NNMF images are much cleaner.

imgRows =
  With[{th = 0.5},
   MapThread[{#1, #2, Binarize[#3, th], 
      Binarize[#4, th]} &, {testImageLabels, noisyTestImages6, 
     ImageAdjust@*ColorNegate /@ denoisedImages, 
     ImageAdjust@*ColorNegate /@ denoisedImagesNNMF}]];
With[{ncol = 5}, 
 Grid[Prepend[Partition[imgRows, ncol], 
   Style[#, Blue, FontFamily -> "Times"] & /@ 
    Table[{"label", "noised", "SVD", "NNMF"}, ncol]]]]

enter image description here

Usually with NNMF in order to get good results we have to do more that one iteration of the whole factorization and reconstruction process. And of course NNMF is much slower. Nevertheless, we can see clear advantages of NNMF's interpretability and leveraging it.

Further experiments

Gallery of experiments other digit pairs

See the gallery with screenshots from other experiments in

Using Classify

Further comparisons can be done using Classify -- for more details see the last sections of the posts linked above.

Here is an image of such comparison:

enter image description here

share|improve this answer

Start data

First let us get some images. I am going to use the MNIST dataset for clarity. (And because I experimented with similar data some time ago.)

MNISTdigits = ExampleData[{"MachineLearning", "MNIST"}, "TestData"];
testImages = RandomSample[Cases[MNISTdigits, (im_ -> 0) :> im], 100]

enter image description here

Let us convince ourselves that all images have the same dimension:

Tally[ImageDimensions /@ testImages]
dims = %[[1, 1]];

(* {{{28, 28}, 100}} *)

Let us create noisy images. Here are the options:

enter image description here

Below I will show results with 2 and 5.

The listing/selection of these options comes from asking ourselves what we know about the data. What does "of the same sample" in the question mean? In terms of digit writing:

  1. is it the same instance of digit writing but obtained by different channels that add noise, or

  2. is it different instances of digit writing obtained by a channel that adds noise?

Noisy images

noisyTestImages = 
  Table[ImageEffect[
    testImages[[12]], {"GaussianNoise", RandomReal[1]}], {100}];

RandomSample[Thread[{testImages[[12]], noisyTestImages}], 15]

enter image description here

noisyTestImages5 = 
  Table[ImageEffect[
    testImages[[i]], {"GaussianNoise", RandomReal[1]}], {i, 100}];
RandomSample[Thread[{testImages, noisyTestImages5}], 15]

enter image description here

Linear vector space representation

Unfold the images into vectors and stack together into a matrix:

noisyTestImagesMat = (Flatten@*ImageData) /@ noisyTestImages;
Dimensions[noisyTestImagesMat]

(* {100, 784} *)

No centralizing and standardizing

The vector representations of the images I use have 0's as the important signal values. It is better to use the negated images as in my second answer (not done in this one).

Initially I thought that centralizing and normalizing is detrimental to the described process. Thanks to arguments in the comments and Rahul's answer now I think that centralizing is just useless for the data in my answer -- see this comparison of centralized and non-centralized results using the code in Rahul's answer.

My other reason for not using centralization is that I would like to have interpretable basis when using matrix factorizations of data. (See the NNMF basis in my second answer.) We will have at least one such vector with SVD if we do not centralize. The question does not specify which of the results of PCA are going to be used: the de-noised images or the obtained components.

De-noising

Apply PCA (SVD):

{U, S, V} = SingularValueDecomposition[noisyTestImagesMat, 12];

Let us look at the singular values:

ListPlot[Diagonal[S], PlotRange -> All, PlotTheme -> "Detailed"]

enter image description here

The main step of the denoising process -- zeroing the basis vectors we think correspond to noise:

dS = S;
Do[dS[[i, i]] = 0, {i, Range[2, 12, 1]}]

Reconstruct the stack of images matrix and plot it:

newMat = U.dS.Transpose[V];

MatrixPlot[newMat]

enter image description here

Convert to denoised images:

denoisedImages = Map[Image[Partition[#, dims[[2]]]] &, newMat];

Compare noisy and denoised images:

Transpose[{noisyTestImages, denoisedImages}]

enter image description here

The harder set

Repeating the above steps for option 5 (all images same noise type with different parameters) we get this:

enter image description here

Not as de-noised as with the more regular set but still good.

Post processing

As it can be seen with the last result the procedure finds the "mean image" ("mean digit" in this case.) This can be exploited with some sort of image normalization:

ImageMultiply[ImageAdd[denoisedImages], 1/Length[denoisedImages]]
Binarize[%, 0.7]

enter image description here

share|improve this answer
1  
Aren't you supposed to subtract the mean vector before applying SVD? – nikie May 6 at 17:10
    
Well, it depends what are your goals. This procedure kind of finds the "mean digit", that being the main vector and the rest noise vectors. Although, some row-wise normalization probably should be done or verified. I was going to add comments and code illustrating that... – Anton Antonov May 6 at 17:28
3  
I don't get it. If I take an uncentered point set, and I look for the vector that will give me the lowest average reconstruction error for that point set (that's what the SVD does), the vector will always point more or less towards the center of the point set, right? If I want that, why not just take the mean? – nikie May 6 at 19:19
    
@nikie 1. Read the question and comments by OP -- his boss wants PCA applied (not mean, or Fourier, or whatever). 2. The question is not clear about this: are we looking for a component or are we trying to de-noise each image separately? The described procedure conveniently de-noises each image and finds good approximation of the signal. 3. The vector representations of the images I use have 0's as the important signal parts. If we normalize, things get messier -- we have to use the negated images if we want to apply normalization. This is not something I wanted to deal with... – Anton Antonov May 6 at 21:29
3  
No, what nikie is suggesting is the standard approach for PCA/SVD-based data processing. I have just tried it on this data and it works fine. Maybe you should try it too. – Rahul May 6 at 23:07

As requested by Anton. I halved the amount of noise because otherwise some images have barely any signal left. As you can see below, we are still putting in a significant amount of noise.

(To conserve space I'm only visualizing the first ten images in this answer, but the denoising is happening over all 100 test images.)

SeedRandom[2016] (* for reproducibility *)
MNISTdigits = ExampleData[{"MachineLearning", "MNIST"}, "TestData"];
images = RandomSample[Cases[MNISTdigits, (im_ -> 0) :> im], 100];
n = Length@images;
noisyImages = (ImageEffect[#, {"GaussianNoise", RandomReal[0.5]}] &) /@ images;
Take[noisyImages, 10]

enter image description here

Convert the data to vectors and find their mean:

toVector = Flatten@*ImageData;
fromVector = Image[Partition[#, 28]] &;
data = toVector /@ noisyImages;
mean = Mean[data];
fromVector[mean]

enter image description here

Center the data by subtracting the mean from all the data points. Now, what we have is how much each image differs from the mean, i.e. the variations in how the letters are drawn, plus the noise in each image.

centeredData = (# - mean &) /@ data;
(ImageAdjust@*fromVector) /@ Take[centeredData, 10]

enter image description here

The singular value decomposition extracts the patterns in these variations, from most to least common.

{u, s, v} = SingularValueDecomposition[centeredData, n];

Apparently you already know which components are important, but if you didn't, you could look at the distribution of singular values:

ListPlot[Diagonal@s, PlotRange -> All]

enter image description here

It looks like the first five stand out from the rest. Let's look at what variations they encode:

(ImageAdjust@*fromVector) /@ Take[Transpose[v], 5]

enter image description here

Zero out all the rest of the singular values and reconstruct the modified data matrix, then the images themselves:

snew = DiagonalMatrix@Table[If[i <= 5, s[[i, i]], 0], {i, n}]
denoisedData = (# + mean &) /@ (u.snew.Transpose[v])
denoisedImages = fromVector /@ denoisedData;
Take[denoisedImages, 10]

enter image description here

Let's review, because the actual denoising code is quite short and simple:

data = toVector /@ noisyImages;
mean = Mean[data];
centeredData = (# - mean &) /@ data;
{u, s, v} = SingularValueDecomposition[centeredData, n];
snew = DiagonalMatrix@Table[If[i <= 5, s[[i, i]], 0], {i, n}]
denoisedData = (# + mean &) /@ (u.snew.Transpose[v])
denoisedImages = fromVector /@ denoisedData;

enter image description here enter image description here

share|improve this answer
4  
Now this is the low-rank approximation I'm accustomed to… :) – J. M. May 7 at 6:59
    
As I mentioned in another comment, I run this answer computations with and without centralizing -- see this image. I do not see any benefit of the centralization. – Anton Antonov May 7 at 22:07
1  
@Anton: The benefit of centralization is that it saves you exactly one singular vector -- the one needed to represent the mean. What would be more illuminating is if you compared centralizing vs. not centralizing for just 1 nonzero singular value, then 2, then 3, and so on. I think your input images have so much noise that by the 5th singular value the components are almost entirely swamped by noise. – Rahul May 8 at 1:09
    
@Rahul I agree with what you are say in principle. Please see the section "No centralizing and standardizing" of my first answer where my point of view in this case is explained. – Anton Antonov May 8 at 13:46
    
@Anton: I've seen that already, it doesn't address what I said in my comment. Anyway I don't think we are going to get anywhere with this discussion. I will refer you instead to the paper the OP has added to the question. Have a nice day. – Rahul May 8 at 14:15

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.