Восстановление размытого изображения
artpgchЛинейный метод наименьших квадратов с ограничениями на примере восстановления размытого изображения
Введение
В данном примере рассмотрим, как восстановить размытое изображение, решая задачу оптимизации методом линейных наименьших квадратов с ограничениями.
Присоединим необходимые библиотеки.
#import Pkg
#Pkg.add(["FileIO", "Images", "ImageShow", "LinearAlgebra", "SparseArrays", "JuMP", "Ipopt", "Printf"])
using FileIO, Images, ImageShow, LinearAlgebra, SparseArrays, JuMP, Ipopt, Printf
Постановка задачи
Имеется фотография, на которой изображён стенд Engee.
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="Исходное изображение $m x $n ($mn пикселей)")
display(p1)
Задача состоит в том, чтобы взять размытую версию этой фотографии и попытаться восстановить её чёткость. Исходное изображение является чёрно-белым, что означает, что его можно представить в виде матрицы размером , состоящей из значений пикселей в диапазоне от до .
Размытие изображения
Смоделируем эффект вертикального размытия в движении, усреднив каждый пиксель с 15 пикселями сверху и снизу. Построим разреженную матрицу D для выполнения размытия умножением матриц.
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);
Отобразим первые 30×30 элементов матрицы D.
println("\nОператор размытия (первые 30x30 элементов):")
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)
Умножив матрицу D на исходное изображение, преобразованное в вектор-столбец, получим размытое изображение.
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="Размытое изображение")
display(p2)
Изображение стало нечётким, надписи нечитаемыми.
Восстановление изображения
Чтобы восстановить чёткость, предположим, что известен оператор размытия . Насколько хорошо можно устранить размытие и восстановить исходное изображение?
Самый простой подход — решить задачу метода наименьших квадратов для :
при условии .
В этой задаче матрица размытия задана, и требуется найти такой вектор , при котором будет максимально близко к . Чтобы решение представляло собой осмысленные значения пикселей, его ограничивают диапазоном от до .
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)
Решим задачу оптимизации.
optimize!(model)
x_val = value.(x)
xpic = reshape(x_val, m, n)
Отобразим восстановленное изображение.
println("\nВосстановленное изображение")
p3 = heatmap(reverse(xpic, dims=1), aspect_ratio=:equal, color=:grays, axis=false, colorbar=false, title="Восстановленное изображение")
display(p3)
Применив решение задачи оптимизации, нам удалось практически идеально восстановить исходное изображение до такой степени, что надписи на нём снова стали различимы и читаемы.
Заключение
В данной работе была рассмотрена практическая задача восстановления изображений на примере устранения вертикального размытия. Экспериментальные результаты убедительно демонстрируют эффективность метода линейных наименьших квадратов с ограничениями для решения обратных задач в обработке изображений, наглядно демонстрируя, как абстрактные математические концепции находят конкретное применение в решении практических проблем цифровой обработки визуальной информации.