A toy example for using pyrpca

Author: Roland Kwitt, Kitware Inc., 2013 (roland.kwitt@kitware.com, rkwitt@gmx.at)

Our problem setting is as follows: We want to load an image (which we know is low-rank due to it's structure), then corrupt the image with random outliers and eventually recover the low-rank part (i.e., the original image).

Let's load SimpleITK and numpy first:



In [78]:

    
import SimpleITK as sitk
import numpy as np

We'll tackle the task using ialm, i.e., the Inexact Lagrangian Multipier (IALM) approach:



In [79]:

    
import core.ialm as ialm

We next load a checkerboard image (our low-rank signal):



In [80]:

    
I = sitk.ReadImage("examples/checkerboard.png")



In [81]:

    
Inp = sitk.GetArrayFromImage(I)



In [82]:

    
imshow(Inp, cm.gray)









    Out[82]:





<matplotlib.image.AxesImage at 0x112214190>

As we can see, the image has a very distinctive structure. In the following, we'll corrupt 30% of the pixel values by pertubing the signal with values in the range [190, 210]. The outliers will be uniformly distributed.



In [83]:

    
p = 0.3
how_many = np.round(N*p)
outliers = np.round(np.random.uniform(-10, 10, how_many))



In [84]:

    
N = np.prod(Inp.shape)
where = np.random.random_integers(0, N-1, how_many)

We create a copy of the original image and set the outliers:



In [85]:

    
Inp_corrupted = np.copy(Inp)



In [97]:

    
Inp_corrupted.ravel()[where] = np.array(200 + outliers, dtype=np.uint8)



In [98]:

    
imshow(Inp_corrupted, cm.gray)









    Out[98]:





<matplotlib.image.AxesImage at 0x1124c1910>

We'll run IALM (with it's default settings) on the image next which gives us the low-rank and the sparse (hopefully the outliers) signal parts:



In [99]:

    
low_rank, sparse, _ = ialm.recover(Inp_corrupted)









    



[iter: 0000]: rank(P) = 0002, |C|_0 = 0000, crit=12.7699549168
[iter: 0010]: rank(P) = 0059, |C|_0 = 61785, crit=0.4638434471
[iter: 0020]: rank(P) = 0002, |C|_0 = 67923, crit=0.0001199796
[iter: 0030]: rank(P) = 0002, |C|_0 = 67923, crit=0.0000000527



In [100]:

    
imshow(low_rank, cm.gray)









    Out[100]:





<matplotlib.image.AxesImage at 0x1124ee6d0>



In [93]:

    
imshow(sparse, cm.gray)









    Out[93]:





<matplotlib.image.AxesImage at 0x111d3f2d0>



In [102]:

    
# compute Frobenius norm of low-rank part minus the original signal
error = np.linalg.norm(low_rank-Inp,ord='fro')
print error









    



0.000109736534477



In [114]:

    
# compute Frobenius norm
error = np.linalg.norm(np.abs(np.round(low_rank))-Inp,'fro')
print error

0.0

So, in the end, we could perfectly separate the low-rank and sparse part of this signal.