Using Finite Element Method for Partial Differential Equation Methods in Image Inpainting

Question

How does one traverse through a corrupted grayscale image fixing all the corrupted pixels in it without checking the pixels that are not corrupted?

Any help would be greatly appreciated.

P.S. I have an algorithm that fixes corrupted pixels by interpolating with its surrounding pixels but I cannot use this algorithm on the entire image as it will take too long.

Here is the algorithm (in pseudocode):

Ek[m, n] = Ik[m − 1, n] + Ik[m + 1, n]+ Ik[m, n − 1] + Ik [m, n + 1] − 4I[m, n]
Ik+1[m, n] = Ik[m, n] + α Ek[m, n]

where error is Ek, α is a coeficient and Ik is the image at the kth iteration.

Answer 1

As far as I understand your question, you're searching for a way of iterating though each pixel of a 2D array and performing some action whenever you find a "corrupted" pixel. Not knowing what if your definition of "corrupted", I'll assume that by "corrupted" you mean isnan - Not A Number. Your algorithm makes no sense then (because you're averaging NaN into your corrected pixel value which will give you a NaN). I will therefore assume that the algorithm is simply replacing the pixel completely with the average of it's nearest non-diagonal neighbors:

Ak+1[m, n] = 1/4 * ( Ak[m − 1, n] + Ak[m + 1, n]+ Ak[m, n − 1] + Ak [m, n + 1] )

and If I've in fact misunderstood you, you'll hopefully be able to work form my answer anyway.

First of all:

% Make a sample 'corrupted' array
A = rand(100,100);
A(A>0.9) = NaN; % ~10% corrupted

Now, the problem is complicated by the fact that this algorithm fails if there are two corrupted pixels next to one another, as well as at the edges of the array. We will worry about that later.

Simple (but not working) aproach

You can find all the corrupted pixels with

nA = isnan(A); % Boolean array the size of A where each corrupted pixel value is equal to 1

and, if you want to, get their indices with

I = find(nA); % Finds the indices of the above using linear indexing
[m,n] = find(nA); % Finds the indices using subscript indexing

If you're new to Matlab, an interesting reading here might be: "logical arrays", "find" and "linear indexing". We will also need the size of A so lets put that into workspace

sA = size(A);

You then iterate though all the indices and apply the algorithm.

for j = 1:length(I)
    i = I(j); %j'th index
    [m,n] = ind2sub(sA,i); % Change linear indices to subscripts
    A(m,n) = 1/4 * ( A(m−1,n) + A(m+1,n)+ A(m,n−1) + A(m,n+1) );
end

More complicated (and actually working) method

As mentioned, the algorithm has problems when there are edges around, and when there are two or more NaN's around (which you can think of as a tiny little edge). In order to deal with the edge problem we will wrap the border around and make your image into a torus (although a 1 pixel padding would work too). We will then iterate though A replacing the pixels according to the nearest neighbor method. I had to allow for 3 and 2 neighbor replacement to make it work on heavily corrupted images. The code can easily be modified to stop that. Anyway, here is the working example:

A = peaks(100);
A(rand(size(A))>0.3) = NaN; % 70% corrupted image
sA = size(A);

nrep = [1 1 1]; % Initialise as nonzero
imagesc(A); drawnow; % Show what is happening
while any(isnan(A(:)))
    getA = @(m,n) A(mod(m-1,sA(1)-1)+1,mod(n-1,sA(2)-1)+1); % Wrap the image around the edges
    numA = @(m,n) sum(isnan([getA(m+1,n),getA(m-1,n),getA(m,n+1),getA(m,n-1)])); % Number of neighbouring NaNs
    for j = 1:numel(A);
        if isnan(A(j));
            [m,n] = ind2sub(sA,j);
            switch numA(m,n)
                case 0
                    cA = 1/4 * ( getA(m-1,n) + getA(m+1,n) + getA(m,n-1) + getA(m,n+1) );
                    nrep(1) = nrep(1) + 1;
                case 1
                    cA1 = 1/3 * ( getA(m+1,n) + getA(m,n-1) + getA(m,n+1) );
                    cA2 = 1/3 * ( getA(m-1,n) + getA(m,n-1) + getA(m,n+1) );
                    cA3 = 1/3 * ( getA(m-1,n) + getA(m+1,n) + getA(m,n+1) );
                    cA4 = 1/3 * ( getA(m-1,n) + getA(m+1,n) + getA(m,n-1) );
                    cAs = [cA1 cA2 cA3 cA4];
                    ok = ~isnan(cAs); ok = find(ok,1); % find first cA# which is not a NaN
                    cA = cAs(ok);
                    nrep(2) = nrep(2) + 1;
                case 2
                        cA1 = 1/2 * ( getA(m+1,n) + getA(m,n+1) );
                        cA2 = 1/2 * ( getA(m+1,n) + getA(m,n-1) );
                        cA3 = 1/2 * ( getA(m-1,n) + getA(m,n+1) );
                        cA4 = 1/2 * ( getA(m-1,n) + getA(m,n-1) );
                        cA5 = 1/2 * ( getA(m+1,n) + getA(m-1,n) );
                        cA6 = 1/2 * ( getA(m,n+1) + getA(m,n-1) );
                        cAs = [cA1 cA2 cA3 cA4 cA5 cA6];
                        ok = ~isnan(cAs); ok = find(ok,1); % find first cA# which is not a NaN
                        cA = cAs(ok);
                        nrep(3) = nrep(3) + 1;
                case 3
                    continue
                case 4
                    continue
            end
        A(j) = cA; % Replace element
        end
    end
    imagesc(A); drawnow; pause(0.01); % Show what is happening
end

Comment out the lines with the imagesc to suppress plotting the figures (cool as it may look ;) ). The code is somewhat ramshackled and could be improved a lot, but it works pretty well as is:

Original

Reconstruction form 50% corruption

Reconstruction form 70% corruption

Hive-mind method (my go-to method)

Know what to google for and hence arrive at this piece of code. Adapt as needed.