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 yList
we 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]]
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: