Message Boards Message Boards

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

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

POSTED BY: Markus Roellig
Answer
1 year ago

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
Answer
1 year ago

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
Answer
1 year ago

enter image description here - you earned "Featured Contributor" badge, congratulations !

This is a great post and it has been selected for the curated Staff Picks group. Your profile is now distinguished by a "Featured Contributor" badge and displayed on the "Featured Contributor" board.

POSTED BY: Moderation Team
Answer
1 year ago

Finite difference method is used for image compression ( http://jre.cplire.ru/iso/nov12/1/text.pdf ) and ( http://elibrary.ru/item.asp?id=27147482 ) and ( http://www.freepatent.ru/images/patents/497/2500067/patent-2500067.pdf ). First formed boundary conditions (this pattern), then pattern compressed through Huffman or arithmetic coding. The method of differential compression works well on gradient images or images containing large fields of the same color or brightness. The essential question is the choice of color space YCrCb or RGB or. Also proposed ( http://jre.cplire.ru/iso/nov12/1/text.pdf ) methods of improving the quality of image reconstruction using prefiltration, with a pattern formed on the original image. enter image description here In the image: a -- the original image (a photo of the Earth remote sensing); б -- boundary conditions (pattern) in red denotes excluded elements; в -- is the restored image. enter image description here On the image: a - the original image (photo from the International Space Station); B - boundary conditions (figure) in red indicate the excluded elements; B is the restored image.

POSTED BY: Alex Alex
Answer
8 months ago

Group Abstract Group Abstract