The Gradient Operator of a Vectorized Image in Matrix Form

Question

I have this optimization problem:

$$ \arg \min_{ X \left( i, j \right) } \sum_{i, j} \left\| X \left( i, j \right) - 255 \right\|_{2}^{2} + \lambda \sum_{i, j} \left\| \nabla X \left( i, j \right) - \nabla Y \left( i, j \right) \right\|_{2}^{2} $$

Where $ X $ is the output image and $ Y $ is the input image.

Let's say the input image is $ y $ and the output image is $ x $ (Transform image to vector) then the problem can be rewritten:

$$ \hat{x} = \arg \min_{x} \frac{1}{2} {\left\| x - 255 \cdot \boldsymbol{1} \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{h} \left( x - y \right) \right\|}_{2}^{2} + \frac{\lambda}{2} {\left\| {D}_{v} \left( x - y \right) \right\|}_{2}^{2} $$

Where $D_h$ is the horizontal Derivative Operator, $D_v$ is the vertical Derivative Operator and $1$ is vector of ones.

Then the solution is given by:

$$ \hat{x} = { \left( I + \lambda {D}_{h}^{T} {D}_{h} + \lambda {D}_{v}^{T} {D}_{v} \right) }^{-1} \left( \lambda {D}_{h}^{T} {D}_{h} y + \lambda {D}_{v}^{T} {D}_{v} y + 255 \cdot \boldsymbol{1} \right) $$

My question is given the input $ y $ how to apply $ {D}_{h} $ and $ {D}_{v} $ to this specific equation.

Thanks for your reply.

Answer 1

It is pretty simple to create those Matrices.
The real issue with them is their size which is enormous for real world images.

For small kernels they are sparse which saves the day.
Indeed for the Derivative Operator, which has only 2 elements, they are highly sparse.

I built them in MATLAB using:

mI = im2double(imread(imageFileName));
mI = mI(11:410, 201:600, 1);

% mI = mI(11:20, 201:210, 1);

numRows     = size(mI, 1);
numCols     = size(mI, 2);
numPixels   = numRows * numCols;

mO = zeros([numRows, numCols, length(vParamLambda)]); %<! Output

vDx = [1, -1]; %<! Matrix is doing Correlation, this for Convolution
vDy = [1; -1];

% Sanity Check
mIxRef = conv2(mI, vDx, 'valid');
mIyRef = conv2(mI, vDy, 'valid');

% mIx = reshape(mDx * mI(:), numRows, numCols - 1);
% mIy = reshape(mDy * mI(:), numRows - 1, numCols);

mDh = sparse(numPixels - numCols, numPixels);
mDv = sparse(numPixels - numRows, numPixels);

% mDx = zeros(numPixels - numCols, numPixels);
% mDy = zeros(numPixels - numRows, numPixels);

tic();
colShift = 0;
for ii = 1:(numPixels - numRows)
    if(mod(ii + colShift, numRows) == 0)
        colShift = colShift + 1;
    end
    mDv(ii, ii + colShift) = -1;
    mDv(ii, ii + colShift + 1) = 1;
end

for ii = 1:(numPixels - numCols)
    mDh(ii, ii) = -1;
    mDh(ii, ii + numCols) = 1;
end
toc();


mIx = reshape(mDh * mI(:), numRows, numCols - 1);
mIy = reshape(mDv * mI(:), numRows - 1, numCols);

mE = mIy - mIyRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Vertical Derivative) - ', num2str(max(abs(mE(:))))]);

mE = mIx - mIxRef;
disp(['Maximum Absolute Error Between Matrix Form to Convolution (Horizontal Derivative) - ', num2str(max(abs(mE(:))))]);

If you load image into mI you see the result is identical to the convolution operator.

The full MATLAB code which implements the solution is at my Signal Processing StackExchange Q50329 - GitHub Repository (Look at the SignalProcessing\Q50329 folder).

Remark
In practice usually those kind of equation are solved using Preconditioned Conjugate Gradient Method.
The nice trick is that method requires only the result of the operators applied on a vectorized image. This can be done using the convolution operator itself instead of using the matrix form.