Image quality not reducing when I try to implement the given Mean structural similarity best-case/worst-case validation model in a paper

Question

Im currently trying to implement the mathematical models in the paper " Image Quality Assessment: From Error Visibility to Structural Similarity" by Zhou Wang.

At section IV.A the author introduced a mathematical model to validate the MSSIM model. The model modify the distroted image Y by increasing/decreasing MSSIM by descending/ascending the gradient of MSSIM w.r.t Y.

When I try to implement image quality descent, the output image quality is not what I expected as the image's quality go better instead of worse.

For the implementation, I used skimage's structural_similariy function imported as ssim. I defined load_gray_img function which loads an image and convert it to grayscale. Function P correspond to the Projection operation in the paper P(X,Y) and unit correspond to Unit vector operation E(X,Y).

Here is my implementation, I want to ask where the problem is and how can it be resolved? I have tried multiple different step_size and MSE constraints it all end up with the similar problem.

import numpy as np
import matplotlib.pyplot as plt
from skimage.metrics import structural_similarity as ssim
from skimage.metrics import mean_squared_error as mse


def load_gray_img(image_path):
    # Load the image
    image = cv.imread(image_path)
    if image is None:
        raise ValueError("Image not found or unable to load.")
    
    # Convert to grayscale
    gray_image = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
    return gray_image

X = load_gray_img("mos/reference_images/I01.BMP")
Y = load_gray_img("mos/distorted_images/I01_14_4.bmp")

plt.imshow(Y, cmap='gray')

final_mssim = ssim(X, Y, win_size=11, gaussian_weights=True, data_range=255)
final_mse = mse(X, Y)
print(f"Initial MSSIM: {final_mssim}")
print(f"Initial MSE: {final_mse}")

def P(X, Y):
    E = unit(X, Y)
    E = np.array(E, dtype=np.float32)
    E_T = E.T
    I = np.eye(E.shape[0])
    return I - np.dot(E, E_T)

def unit(X, Y):
    return (Y - X) / np.linalg.norm(Y - X, ord='fro')

lmbda = 0.5  # step size
sigma = 50  # MSE constraint

def quality_descent(X, Y, iterations=100):
    mse_list = []
    mssim_list = []
    for _ in range(iterations):
        current_mssim, grad = ssim(X, Y, gradient=True, win_size=11, gaussian_weights=True, data_range=255)
        Y = Y + lmbda * np.dot(P(X, Y), grad)
        Y = X + sigma * unit(X, Y)
        current_mse = mse(X, Y)
        mse_list.append(current_mse)
        mssim_list.append(current_mssim)
    return Y, mse_list, mssim_list

# Perform quality descent
Y, mse_list, mssim_list = quality_descent(X, Y)

# Plotting MSE and MSSIM over iterations
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(mse_list, label="MSE")
plt.xlabel("Iterations")
plt.ylabel("MSE")
plt.title("MSE over Iterations")
plt.legend()

plt.subplot(1, 2, 2)
plt.plot(mssim_list, label="MSSIM", color="orange")
plt.xlabel("Iterations")
plt.ylabel("MSSIM")
plt.title("MSSIM over Iterations")
plt.legend()
plt.tight_layout()
plt.show()

# Display the final image
plt.imshow(Y, cmap='gray')
plt.title("Final Image after Quality change")
plt.show()

# Print final MSSIM and MSE values
final_mssim = ssim(X, Y, win_size=11, gaussian_weights=True, data_range=255)
final_mse = mse(X, Y)
print(f"Final MSSIM: {final_mssim}")
print(f"Final MSE: {final_mse}")

Output:
Initial MSSIM: 0.8459229680970184
Initial MSE: 295.9006754557292

Final MSSIM: 0.9999909113990414
Final MSE: 0.01271565755208334