# 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) – Marcus Müller Jul 6 '18 at 21:48
• i edited my question, also here is the reference dsp.stackexchange.com/questions/50329/… – lafi raed Jul 6 '18 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! – Marcus Müller Jul 7 '18 at 0:06
• i add the resolution details – lafi raed Jul 7 '18 at 0:52
• That's just now copies of the answer you've got there. What is your question, precisely? – Marcus Müller Jul 7 '18 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 Question 50329 - GitHub Repository.

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 ? – lafi raed Jul 7 '18 at 12:34
• @lafiraed, That belongs to the other question you posted not here. – Royi Jul 7 '18 at 13:18