# The Gradient Operator of a Vectorized Image in Matrix Form

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.

• I'm confused how you can write your in- and output images as vectors $x$ and $y$ and still differentiate them. I'm missing something here, or you're doing a transformation that doesn't work. What exactly is the set from which $x$ and $y$ come? (also, bad choice re-using letters that you already used as coordinates in the first equation) Jul 6, 2018 at 21:48
• i edited my question, also here is the reference dsp.stackexchange.com/questions/50329/… Jul 6, 2018 at 23:57
• Why did you remove thee better half of your question? That was your own attempt! Without that, your question really is just "solve this problem for me"! Instead of deleting, you should have explained your approach! Jul 7, 2018 at 0:06
• i add the resolution details Jul 7, 2018 at 0:52
• That's just now copies of the answer you've got there. What is your question, precisely? Jul 7, 2018 at 0:53

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.

• what’s the role of lambda ? Jul 7, 2018 at 12:34
• @lafiraed, That belongs to the other question you posted not here.
– Royi
Jul 7, 2018 at 13:18