Is MATLAB's PIQE function wrong?

I'm trying to train a Deep Learning model for image super resolution, and I wanted to implement the PIQE score as a loss function. Since I will be training the model with pytorch, I was trying to make an own Python implementation of the algorithm to compute the PIQE score.

As a first step, I looked at the MATLAB implementation of piqe (the link takes you to the main page but I am looking at the source code) to see how it's done and then adapt it to Python. There is one thing that bothers me, however.

The PIQE score starts off by calculating the Mean-Substracted Contrast-Normalized coefficients with the following formula:

But the matlab code at that step looks like this:

mu = imgaussfilt(ipImage,7/6,'FilterSize',7,'Padding','replicate');
sigma = sqrt(abs(imgaussfilt(ipImage.*ipImage,7/6,'FilterSize',7,'Padding','replicate') - mu.*mu));
imnorm = (ipImage-mu)./(sigma+1);

I'm puzzled about the calculation of the variance, sigma. In the algorithm of the paper, at each pixel, the mean of the 7x7 neighborhood is calculated and then subtracted from each value of said 7x7 neighborhood. Then, the differences are squared and multiplied by its corresponding Gaussian weight w(k,l)

Instead, the MATLAB algorithm multiplies the Gaussian weighting (by using imgausssfilt) with the squared pixel values, and then subtracts the squared means from that matrix, taking the absolute values of that operation. Correct me if I'm wrong, but isn't this a case of mistakenly using (a-b)² = a² - b² ?

Basically my question is if you could kindly confirm whether what I said before is true, or I misinterpreted the MATLAB code. Thanks in advance!

I understand why you are confused, but both are right. It uses the classical identity

Var(X) = E [(X-E(X))^2]= E(X^2) - E(X)^2

Just multiply your (I-mu)^2 out and compare the result with the definition of mu, then you will see that they cancel.