Message Boards Message Boards

Repair damaged image pixel using Laplace Interpolation

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).

enter image description here

The Laplace Interpolation

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

enter image description here

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: enter image description 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];
  repl = Thread[unknowns -> sol];
  Partition[yList /. repl, xdim]
  ]

Let's add some convenience functions:

 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];
  black = ConstantArray[bad, {xdim, ydim}];
  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]]

enter image description here

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]

enter image description here

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]

enter image description here

Attachments:
POSTED BY: Markus Roellig
4 Replies
Sort By:
Alex Alex
Alex Alex
Posted 7 years ago
POSTED BY: Alex Alex
POSTED BY: Moderation Team

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.

enter image description here
POSTED BY: Matthias Odisio

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.
POSTED BY: Markus Roellig
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract