# Repair damaged image pixel using Laplace Interpolation

GROUPS:

I came across Laplace Interpolation as a specialized interpolation method for restoring missing data on a grid and was so delighted by the results that I thought it might be worth sharing. (see e.g. here: http://numerical.recipes/CS395T/lectures2010/201019LaplaceInterpolation.pdf)

So the basic idea of Laplace interpolation is to set $y(x_i)=y_i$ at every known data point and solve $\nabla^2 y=0$ at every unknown point.

Before I show the code let's take a look at the results. The following shows, the original image (left), the damaged image with 40% healthy pixels (middle) and the repaired image (right).

## The Laplace Interpolation

Each pixel gives one equation. Consider an arbitrary pixel $y_0$. Its value is determined by the neighborhood:

To compute $y_0$

If the pixel is ok:
$$y_0=y_{0,measured}$$
If the pixel is damaged;
$$y_0-\frac{1}{4}y_u-\frac{1}{4}y_d-\frac{1}{4}y_l-\frac{1}{4}y_r=0$$

and the following special cases (boundary):

(left and right boundary)
$$y_0-\frac{1}{2}y_u-\frac{1}{2}y_d=0$$
(top and bottom boundary)
$$y_0-\frac{1}{2}y_l-\frac{1}{2}y_r=0$$
(top left corner)
$$y_0-\frac{1}{2}y_r-\frac{1}{2}y_d=0$$
(top right corner)
$$y_0-\frac{1}{2}y_l-\frac{1}{2}y_d=0$$ (bottom left corner)
$$y_0-\frac{1}{2}y_r--\frac{1}{2}y_u=0$$
(bottom right corner)
$$y_0-\frac{1}{2}y_l-\frac{1}{2}y_u=0$$ There is exactly one equation for each grid point, so we can solve this as a giant (sparse!) linear system,

## The Algorithm

I decided to flatten the image into a 1-D list enumerating as shown here:

Taking this scheme I set up the test functions to decide whether a given index belongs to the image edges or corners;

topLeftCornerQ[index_Integer,{xdim_Integer,ydim_Integer}]:=index==1
topRightCornerQ[index_Integer,{xdim_Integer,ydim_Integer}]:=index==xdim
bottomRightCornerQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(index==(xdim*ydim))
bottomLeftCornerQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(index==xdim*(ydim-1)+1)
leftEdgeQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(Mod[index,xdim]-1==0)&&Not[bottomLeftCornerQ[index,{xdim,ydim}]]&&Not[topLeftCornerQ[index,{xdim,ydim}]]
rightEdgeQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(Mod[index,xdim]==0)&&Not[bottomRightCornerQ[index,{xdim,ydim}]]&&Not[topRightCornerQ[index,{xdim,ydim}]]
topEdgeQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(2<= index<=xdim-1)
bottomEdgeQ[index_Integer,{xdim_Integer,ydim_Integer}]:=(xdim*(ydim-1)+2<= index<=xdim*ydim-1)


We then need some functions to select the boundary values:

yLeft[index_,{xdim_,ydim_}]:=Module[{pos},
Which[
leftEdgeQ[index,{xdim ,ydim }],{},
topLeftCornerQ[index,{xdim ,ydim }],{},
bottomLeftCornerQ[index,{xdim ,ydim }],{},
True,index-1
]
]
yRight[index_,{xdim_,ydim_}]:=Module[{pos},
Which[
rightEdgeQ[index,{xdim ,ydim }],{},
topRightCornerQ[index,{xdim ,ydim }],{},
bottomRightCornerQ[index,{xdim ,ydim }],{},
True,index+1
]
]
yDown[index_,{xdim_,ydim_}]:=Module[{pos},
Which[
bottomEdgeQ[index,{xdim ,ydim }],{},
bottomLeftCornerQ[index,{xdim ,ydim }],{},
bottomRightCornerQ[index,{xdim ,ydim }],{},
True,index+xdim
]
]
yUp[index_,{xdim_,ydim_}]:=Module[{pos},
Which[
topEdgeQ[index,{xdim ,ydim }],{},
topLeftCornerQ[index,{xdim ,ydim }],{},
topRightCornerQ[index,{xdim ,ydim }],{},
True,index-xdim
]
]


(The sanity checks inside can also be removed).

The input to the interpolation is an array where the damaged pixel are indicated by their value $10^{99}$ Every instance of $10^{99}$ is replaced by a enumerated variable $y_i$, where $i$ is the pixel index.

  {xdim , ydim} = Dimensions[array];
arrayList = Flatten[array, 1];
symbolList = ParallelTable[ToExpression["y" <> ToString[i]], {i, 1, xdim*ydim}];
(*replace unknown pixel with y-variable*)
yList = MapIndexed[If[#1 == 10^99, symbolList[[First[#2]]], #1] &,arrayList];


With yListwe can now construct the list of equations to solve:

createEquations[yList_, {xdim_, ydim_}] := Module[{value},
DeleteCases[ParallelTable[
value = yList[[i]];
Which[
NumberQ[value], {},
topLeftCornerQ[i, {xdim, ydim}], value - 0.5 yList[[yRight[i, {xdim, ydim}]]] - 0.5 yList[[yDown[i, {xdim, ydim}]]] == 0.,
topRightCornerQ[i, {xdim, ydim}], value - 0.5 yList[[yLeft[i, {xdim, ydim}]]] - 0.5 yList[[yDown[i, {xdim, ydim}]]] == 0.,
bottomLeftCornerQ[i, {xdim, ydim}], value - 0.5 yList[[yRight[i, {xdim, ydim}]]] - 0.5 yList[[yUp[i, {xdim, ydim}]]] == 0.,
bottomRightCornerQ[i, {xdim, ydim}], value - 0.5 yList[[yLeft[i, {xdim, ydim}]]] - 0.5 yList[[yUp[i, {xdim, ydim}]]] == 0.,
leftEdgeQ[i, {xdim, ydim}], value - 0.5 yList[[yDown[i, {xdim, ydim}]]] - 0.5 yList[[yUp[i, {xdim, ydim}]]] == 0.,
rightEdgeQ[i, {xdim, ydim}], value - 0.5 yList[[yDown[i, {xdim, ydim}]]] - 0.5 yList[[yUp[i, {xdim, ydim}]]] == 0.,
topEdgeQ[i, {xdim, ydim}], value - 0.5 yList[[yLeft[i, {xdim, ydim}]]] - 0.5 yList[[yRight[i, {xdim, ydim}]]] == 0.,
bottomEdgeQ[i, {xdim, ydim}], value - 0.5 yList[[yLeft[i, {xdim, ydim}]]] - 0.5 yList[[yRight[i, {xdim, ydim}]]] == 0.,
True, value - 0.25 yList[[yLeft[i, {xdim, ydim}]]] - 0.25 yList[[yRight[i, {xdim, ydim}]]] - 0.25 yList[[yUp[i, {xdim, ydim}]]] - 0.25 yList[[yDown[i, {xdim, ydim}]]]], {i, 1,
Length[yList]}], {}]]


Putting it all together and adding the LinearSolver as well as reformate the output array:

laplaceInterpolate[array_] :=
Module[{xdim , ydim , list, yList, symbolList, arrayList, unknowns,
equations, coeffArrays, rhs, m, sol, repl},
{xdim , ydim} = Dimensions[array];
arrayList = Flatten[array, 1];
symbolList =
ParallelTable[ToExpression["y" <> ToString[i]], {i, 1, xdim*ydim}];
(*replace unknown pixel with y-variable*)
yList =
MapIndexed[If[#1 == 10^99, symbolList[[First[#2]]], #1] &,
arrayList];
unknowns = DeleteCases[yList, _?NumberQ];
equations = createEquations[yList, {xdim , ydim}];
(*Print[Row[{Length[unknowns],": damaged pixel,   ",Length[
equations],": equations"}]];*)
{rhs, m} = CoefficientArrays[equations, unknowns];
sol = Abs@LinearSolve[m, rhs];
Partition[yList /. repl, xdim]
]


 damageImage[im_Image, healthyPixel_Integer] :=
Module[{imDat, xdim, ydim, black, points, bad},
imDat = ImageData[im];
bad = If[ImageChannels[im] == 3, {10^99, 10^99, 10^99}, 10^99];
{xdim, ydim} = ImageDimensions[im];
points =
RandomSample[Flatten[Table[{i, j}, {i, 1, xdim}, {j, 1, ydim}], 1],
healthyPixel];
Table[black[[points[[i, 1]], points[[i, 2]]]] =
imDat[[points[[i, 1]], points[[i, 2]]]], {i, 1, Length[points]}];
black]
repairImage[im_Image] := Module[{repaired, channels},
If[ImageChannels[im] == 1,
repaired = Image@laplaceInterpolate[ImageData[im]],
channels = ImageData /@ ColorSeparate[im];
repaired =
ColorCombine[Image /@ (laplaceInterpolate[#] & /@ channels),
"RGB"]];
repaired]


## Action

We can now test the procedure. The function damageImage takes an input image and keeps healthyPixel random pixel positions and destroys all other image pixel. I resize the image to 128x128 to keep the computing times low. Any performance improvements are appreciated.

 im = ImageResize[ExampleData[{"TestImage", "Man"}], {128}];
dam = Image@damageImage[im, Round[0.9 128*128]];
rep = repairImage[dam];
GraphicsRow[{
Labeled[Image[im, ImageSize -> 250], "Original", ImageMargins -> 10],
Labeled[Image[dam, ImageSize -> 250], "Damaged (10%)", ImageMargins -> 10],
Labeled[Image[rep, ImageSize -> 250], "Repaired", ImageMargins -> 10]}, ImageSize -> 800, Spacings -> 0]]


It is amazing how much information can be restored:

 im = ImageResize[ExampleData[{"TestImage", "Man"}], {128}];
dam = Image@damageImage[im, Round[0.1 128*128]];
rep = repairImage[dam];
GraphicsRow[{
Labeled[Image[im, ImageSize -> 250], "Original", ImageMargins -> 10],
Labeled[Image[dam, ImageSize -> 250], "Damaged (90%)",ImageMargins -> 10],
Labeled[Image[rep, ImageSize -> 250], "Repaired", ImageMargins -> 10]}, ImageSize -> 800, Spacings -> 0]


In case of an RGB image each channel has to be interpolated seperately.

 im = ImageResize[ExampleData[{"TestImage", "Lena"}], {128}];
dam = Image@damageImage[im, Round[0.25 128*128]];
rep = repairImage[dam];
GraphicsRow[{
Labeled[Image[im, ImageSize -> 250], "Original", ImageMargins -> 10],
Labeled[Image[dam, ImageSize -> 250], "Damaged (75%)", ImageMargins -> 10],
Labeled[Image[rep, ImageSize -> 250], "Repaired", ImageMargins -> 10]}, ImageSize -> 800, Spacings -> 0]


Attachments:
2 years ago
4 Replies
 Matthias Odisio 3 Votes That looks like good restoration, Markus. Thanks for sharing it!For comparison, see below the results of the function Inpaint with two different methods yielding similar results. You'll also be pleased with the speed.
 Markus Roellig 1 Vote Matthias, thanks for pointing out Inpaint, of course there is already a built-in function ;) I still would be curious how much one could speed up my code - probably not up to compete with the internal functions, but who knows.