function [qm,qsd] = symPC(imageRef,imageDis)
% =============================================================
% Copyright(c) 2023 Congmin Chen, Xuanqin Mou
% All Rights Reserved.
% 
% symPC- Measure the image quality of distorted image with the reference image using circular symmetric Phase Congruency.
% 
% inputs:
% 
% imageRef - the reference image
% imageDis - the distorted image
% The dynamic range of input images should be 0-255.
% 
% outputs:
% 
% qm - image score with average pooling
% qsd - image score with standard deviation pooling
% 
% This is an implementation of the following algorithm:
% Congmin Chen, Xuanqin Mou,
% "Phase Congruency Based on Derivatives of Circular Symmetric Gaussian Function:
% An Efficient Feature Map for Image Quality Assessment"
% =============================================================

imageRef = double(imageRef);
imageDis = double(imageDis);
[rows, cols] = size(imageRef(:,:,1));

%%%%%%%%% settings %%%%%%%%%%%
minDimension = min(rows,cols);
F = max(1,round(minDimension / 256));
aveKernel = fspecial('average',F); % downsampling filter
SIG = [0.3,0.6]; % scale selection
ep = 25;
c0=120; % constant in divisive normalization
c1 = 3e-05; % constant in quality map computing
%%%%%%%%%%%%%%%%%%%%%%%%

if ndims(imageRef) == 3 % color images
    Y1 = 0.299 * double(imageRef(:,:,1)) + 0.587 * double(imageRef(:,:,2)) + 0.114 * double(imageRef(:,:,3));
    Y2 = 0.299 * double(imageDis(:,:,1)) + 0.587 * double(imageDis(:,:,2)) + 0.114 * double(imageDis(:,:,3));
    I1 = 0.596 * double(imageRef(:,:,1)) - 0.274 * double(imageRef(:,:,2)) - 0.322 * double(imageRef(:,:,3));
    I2 = 0.596 * double(imageDis(:,:,1)) - 0.274 * double(imageDis(:,:,2)) - 0.322 * double(imageDis(:,:,3));
    Q1 = 0.211 * double(imageRef(:,:,1)) - 0.523 * double(imageRef(:,:,2)) + 0.312 * double(imageRef(:,:,3));
    Q2 = 0.211 * double(imageDis(:,:,1)) - 0.523 * double(imageDis(:,:,2)) + 0.312 * double(imageDis(:,:,3));
    
    % downsample on color components
    aveI1 = conv2(I1, aveKernel,'same');
    aveI2 = conv2(I2, aveKernel,'same');
    I1 = aveI1(1:F:rows,1:F:cols);
    I2 = aveI2(1:F:rows,1:F:cols);
    
    aveQ1 = conv2(Q1, aveKernel,'same');
    aveQ2 = conv2(Q2, aveKernel,'same');
    Q1 = aveQ1(1:F:rows,1:F:cols);
    Q2 = aveQ2(1:F:rows,1:F:cols);    

else % grayscale images
    Y1 = imageRef;
    Y2 = imageDis;
end

% downsample on Luma channel
aveY1 = conv2(Y1, aveKernel, 'same');
aveY2 = conv2(Y2, aveKernel, 'same');
Y1 = aveY1(1:F:end,1:F:end);
Y2 = aveY2(1:F:end,1:F:end);

[rows, cols] = size(Y1);
GM1 = ones(rows,cols,length(SIG));   GM2 = GM1;
L1 = ones(rows,cols,length(SIG));        L2 = L1;
A1 = ones(rows,cols,length(SIG));       A2 = A1;

% compute PC using GM and LoG
for lev = 1:length(SIG)
    
    sigma = SIG(lev);
    win = 2*ceil(3*sigma)+1;
    
    % GM
    winG = fspecial('gaussian', win, sigma);
    imG1 = filter2(winG, Y1, 'same');
    imG2 = filter2(winG, Y2, 'same');
    [Gx1, Gy1] = gradient(imG1);
    [Gx2, Gy2] = gradient(imG2);
    GM1(:,:,lev) = sqrt(Gx1.^2 + Gy1.^2);
    GM2(:,:,lev) = sqrt(Gx2.^2 + Gy2.^2);
    
    % LoG
    winL = fspecial('log', win, sigma); %LoG  filter
    winL = winL/sum(abs(winL(:)));
    L1(:,:,lev) = filter2(winL, Y1, 'same');
    L2(:,:,lev) = filter2(winL, Y2, 'same');
    
    % divisive normalization
    win2 = fspecial('gaussian', 2*win+3, 2*sigma);
    win2 = win2/sum(abs(win2(:)));
    GM1(:,:,lev) = GM1(:,:,lev)./(c0+sqrt(filter2(win2, GM1(:,:,lev).^2, 'same')));
    GM2(:,:,lev) = GM2(:,:,lev)./(c0+sqrt(filter2(win2, GM2(:,:,lev).^2, 'same')));
    L1(:,:,lev) = L1(:,:,lev)./(c0+sqrt(filter2(win2, L1(:,:,lev).^2, 'same')));
    L2(:,:,lev) = L2(:,:,lev)./(c0+sqrt(filter2(win2, L2(:,:,lev).^2, 'same')));    
    
    % local amplitude
    A1(:,:,lev) = sqrt(GM1(:,:,lev).^2 +L1(:,:,lev).^2);
    A2(:,:,lev) = sqrt(GM2(:,:,lev).^2 + L2(:,:,lev).^2);
end

% local energy 
H1 = sum(GM1,3); H2 = sum(GM2,3);
F1 = sum(abs(L1),3);  F2 = sum(abs(L2),3);
E1 = sqrt( F1.^2 + H1.^2 );
E2 = sqrt( F2.^2 + H2.^2 );

% Phase Congruency
PC1 = E1./(ep+sum(A1,3));
PC2 = E2./(ep+sum(A2,3));

% quality map
quality_map = (2 * PC1.*PC2 + c1)./(PC1.^2 + PC2.^2 + c1);

if ndims(imageRef) == 3 %color images
    c2 = 200; %fixed
    c3 = 200; %fixed
    ISimMatrix = (2 * I1 .* I2 + c2) ./ (I1.^2 + I2.^2 + c2);
    QSimMatrix = (2 * Q1 .* Q2 + c3) ./ (Q1.^2 + Q2.^2 + c3);
    quality_map = quality_map .* real((ISimMatrix .* QSimMatrix).^0.03);
end

qm =  (1-mean2(quality_map)).^(1/3);
qsd = (std2(quality_map)).^(1/3);

end




