Крупномасштабный линейный метод наименьших квадратов с ограничениями, основанный на проблеме

Этот пример показов, как восстановить размытое изображение путем решения крупномасштабной задачи линейного метода наименьших квадратов оптимизации с ограничениями. В примере используется основанный на проблеме подход. Для основанного на решателе подхода смотрите Крупномасштабный Ограниченный Линейный метод наименьших квадратов, Основанный на решателе.

Проблема

Вот фото людей, сидящих в машине с интересным номерным знаком.

load optdeblur
[m,n] = size(P);
mn = m*n;
figure
imshow(P);
colormap(gray);
axis off image;
title([int2str(m) ' x ' int2str(n) ' (' int2str(mn) ') pixels'])

Figure contains an axes. The axes with title 149 x 311 (46339) pixels contains an object of type image.

Проблема в том, чтобы взять размытую версию этой фотографии и попытаться развенчать ее. Начальное изображение является черно-белым, что означает, что оно состоит из пиксельных значений от 0 до 1 в матрице m x n P.

Добавить движение

Симулируйте эффект размытия вертикального движения путем усреднения каждого пикселя с 5 пикселями выше и ниже. Создайте разреженную матрицу D чтобы размыть с одной матрицей умножить.

blur = 5;
mindex = 1:mn;
nindex = 1:mn;
for i = 1:blur
  mindex=[mindex i+1:mn 1:mn-i];
  nindex=[nindex 1:mn-i i+1:mn];
end
D = sparse(mindex,nindex,1/(2*blur+1));

Нарисуйте картину Д.

cla
axis off ij
xs = 31;
ys = 15;
xlim([0,xs+1]);
ylim([0,ys+1]);
[ix,iy] = meshgrid(1:(xs-1),1:(ys-1));
l = abs(ix-iy) <= blur;
text(ix(l),iy(l),'x')
text(ix(~l),iy(~l),'0')
text(xs*ones(ys,1),1:ys,'...');
text(1:xs,ys*ones(xs,1),'...');
title('Blurring Operator D (x = 1/11)')

Умножьте изображение P на матрицу D, чтобы создать размытое изображение G.

G = D*(P(:));
figure
imshow(reshape(G,m,n));

Figure contains an axes. The axes contains an object of type image.

Изображение намного менее отчётливо; вы больше не можете считать номерной знак.

Отлаженное изображение

Чтобы деблур, предположим, что вы знаете размытого оператора D. Как хорошо можно удалить размытие и восстановить оригинальное изображение P?

Самый простой подход состоит в том, чтобы решить задачу методом наименьших квадратов для x:

min(Dx-G2) при условии, что 0x1.

Эта задача принимает размытую матрицу D, как задано, и пытается найти x, который делает Dx ближе всего к G = DP. В порядок для решения, чтобы представлять разумные пиксельные значения, ограничьте решение от 0 до 1.

x = optimvar('x',mn,'LowerBound',0,'UpperBound',1);
expr = D*x-G;
objec = expr'*expr;
blurprob = optimproblem('Objective',objec);
sol = solve(blurprob);
Solving problem using quadprog.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
xpic = reshape(sol.x,m,n);
figure
imshow(xpic)
title('Deblurred Image')

Figure contains an axes. The axes with title Deblurred Image contains an object of type image.

Окрашенное изображение намного яснее, чем размытое изображение. Можно еще раз прочитать номерной знак. Однако отложенное изображение имеет некоторые программные продукты, такие как горизонтальные полосы в нижней правой области дорожного покрытия. Возможно, эти программные продукты могут быть удалены регуляризацией.

Регуляризация

Регуляризация является способом сглаживания решения. Существует много методов регуляризации. Для простого подхода добавьте термин к целевой функции следующим образом:

min((D+εI)x-G2) при условии, что 0x1.

ТерминεI делает полученную квадратичную задачу более стабильной. Бери ε=0.02 и решить проблему снова.

addI = speye(mn);
expr2 = (D + 0.02*addI)*x - G;
objec2 = expr2'*expr2;
blurprob2 = optimproblem('Objective',objec2);
sol2 = solve(blurprob2);
Solving problem using quadprog.

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.
xpic2 = reshape(sol2.x,m,n);
figure
imshow(xpic2)
title('Deblurred Regularized Image')

Figure contains an axes. The axes with title Deblurred Regularized Image contains an object of type image.

Судя по всему, эта простая регуляризация не удаляет программные продукты.

Похожие темы