Linear least squares method with constraints on the example of blurred image restoration
Introduction
In this example, let's look at how to restore a blurred image by solving an optimization problem using the linear least squares method with constraints.
We will attach the necessary libraries.
# import Pkg
# Pkg.add(["FileIO", "Images", "ImageShow", "LinearAlgebra", "SparseArrays", "JuMP", "Ipopt", "Printf"])
using FileIO, Images, ImageShow, LinearAlgebra, SparseArrays, JuMP, Ipopt, Printf
Setting the task
There is a photo showing the Engee booth.
P = load("image.jpg")
P = Gray.(P)
P = Float64.(channelview(P))
m, n = size(P)
mn = m * n
p1 = heatmap(reverse(P, dims=1), aspect_ratio=:equal, color=:grays, axis=false, colorbar=false, title="The original image is $m x $n ($mn pixels)")
display(p1)
The task is to take a blurred version of this photo and try to restore its clarity. The original image is black and white, which means that it can be represented as a matrix. size , consisting of pixel values ranging from before .
Blurring the image
Let's simulate the effect of vertical motion blur by averaging each pixel with 15 pixels at the top and bottom. Let's construct a sparse matrix D to perform blurring by matrix multiplication.
blur = 15
weight = 1/(2*blur + 1)
mindex = Int[]
nindex = Int[]
values = Float64[]
for i = 1:mn
push!(mindex, i)
push!(nindex, i)
push!(values, weight)
for k = 1:blur
if i + k <= mn
push!(mindex, i + k)
push!(nindex, i)
push!(values, weight)
end
if i - k >= 1
push!(mindex, i - k)
push!(nindex, i)
push!(values, weight)
end
end
end
D = sparse(mindex, nindex, values, mn, mn);
Let's display the first 30×30 elements of the matrix D.
println("\Blur iterator (first 30x30 elements):")
for i in 1:30
for j in 1:30
if abs(i-j) <= blur
print("x ")
else
print("0 ")
end
end
println("...")
end
println("..."^10)
Multiplying the matrix D by the original image transformed into a column vector, we get a blurred image.
G_vec = D * vec(P)
G = reshape(G_vec, m, n)
p2 = heatmap(reverse(G, dims=1), aspect_ratio=:equal, color=:grays, axis=false, colorbar=false, title="Blurred image")
display(p2)
The image has become blurry, and the labels are unreadable.
Image Recovery
To restore clarity, assume that the blur operator is known. . How well can the blur be eliminated and the original image restored?
The simplest approach is to solve the least squares problem for :
on condition .
In this problem, the blur matrix is set, and it is required to find such a vector , in which it will be as close as possible to . In order for the solution to represent meaningful pixel values, it is limited to a range from before .
model = Model(Ipopt.Optimizer)
set_silent(model)
@variable(model, 0 <= x[1:mn] <= 1)
@expression(model, expr, D * x - G_vec)
@objective(model, Min, expr' * expr)
Let's solve the optimization problem.
optimize!(model)
x_val = value.(x)
xpic = reshape(x_val, m, n)
We will display the restored image.
println("\The restored image")
p3 = heatmap(reverse(xpic, dims=1), aspect_ratio=:equal, color=:grays, axis=false, colorbar=false, title="The restored image")
display(p3)
By applying the optimization problem solution, we were able to almost perfectly restore the original image to such an extent that the inscriptions on it became distinguishable and readable again.
Conclusion
In this paper, the practical task of image restoration was considered using the example of eliminating vertical blurring. The experimental results convincingly demonstrate the effectiveness of the linear least squares method with constraints for solving inverse problems in image processing, clearly demonstrating how abstract mathematical concepts find concrete application in solving practical problems of digital visual information processing.