Skip to main content

Image Quality Assessment Index VSI: A Visual Saliency-Induced Index MATLAB Code (VSI.m)

For Full Article Refer http://ieeexplore.ieee.org/abstract/document/6873260/


Perceptual image quality assessment (IQA) aims to use computational models to measure the image quality in consistent with subjective evaluations. Visual saliency (VS) has been widely studied by psychologists, neurobiologists, and computer scientists during the last decade to investigate, which areas of an image will attract the most attention of the human visual system. Intuitively, VS is closely related to IQA in that suprathreshold distortions can largely affect VS maps of images. With this consideration, we propose a simple but very effective full reference IQA method using VS. In our proposed IQA model, the role of VS is twofold. First, VS is used as a feature when computing the local quality map of the distorted image. Second, when pooling the quality score, VS is employed as a weighting function to reflect the importance of a local region. The proposed IQA index is called visual saliency-based index (VSI). Several prominent computational VS models have been investigated in the context of IQA and the best one is chosen for VSI. Extensive experiments performed on four largescale benchmark databases demonstrate that the proposed IQA index VSI works better in terms of the prediction accuracy than all state-of-the-art IQA indices we can find while maintaining a moderate computational complexity.

At times, i have faced problem in obtaining the code from official site. Here is the MATLAB Code: 

function sim = VSI(image1, image2)
% ========================================================================
% Visual Saliency based Index
% Copyright(c) 2013 Lin ZHANG
% All Rights Reserved.
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is hereQ
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
% For more information, please refer to
% Lin Zhang et al., "VSI: A Visual Saliency Induced Index for Perceptual
% Image Quality Assessment", submitted to TIP
%----------------------------------------------------------------------
%
%Input : (1) image1: the first image being compared, which is a RGB image
%        (2) image2: the second image being compared, which is a RGB image
%
%Output: sim: the similarity score between two images, a real number
%        
%-----------------------------------------------------------------------

constForVS = 1.27;%fixed
constForGM = 386;%fixed
constForChrom = 130;%fixed
alpha = 0.40;%fixed
lambda = 0.020;%fixed
sigmaF = 1.34;%fixed donot change
omega0 = 0.0210;%fixed
sigmaD = 145;%fixed
sigmaC = 0.001;%fixed

%compute the visual saliency map using SDSP
saliencyMap1 = SDSP(image1,sigmaF,omega0,sigmaD,sigmaC);
saliencyMap2 = SDSP(image2,sigmaF,omega0,sigmaD,sigmaC);

[rows, cols, junk] = size(image1);
%transform into an opponent color space
L1 = 0.06 * double(image1(:,:,1)) + 0.63 * double(image1(:,:,2)) + 0.27 * double(image1(:,:,3));
L2 = 0.06 * double(image2(:,:,1)) + 0.63 * double(image2(:,:,2)) + 0.27 * double(image2(:,:,3));
M1 = 0.30 * double(image1(:,:,1)) + 0.04 * double(image1(:,:,2)) - 0.35 * double(image1(:,:,3));
M2 = 0.30 * double(image2(:,:,1)) + 0.04 * double(image2(:,:,2)) - 0.35 * double(image2(:,:,3));
N1 = 0.34 * double(image1(:,:,1)) - 0.60 * double(image1(:,:,2)) + 0.17 * double(image1(:,:,3));
N2 = 0.34 * double(image2(:,:,1)) - 0.60 * double(image2(:,:,2)) + 0.17 * double(image2(:,:,3));
%%%%%%%%%%%%%%%%%%%%%%%%%
% Downsample the image
%%%%%%%%%%%%%%%%%%%%%%%%%
minDimension = min(rows,cols);
F = max(1,round(minDimension / 256));
aveKernel = fspecial('average',F);

aveM1 = conv2(M1, aveKernel,'same');
aveM2 = conv2(M2, aveKernel,'same');
M1 = aveM1(1:F:rows,1:F:cols);
M2 = aveM2(1:F:rows,1:F:cols);

aveN1 = conv2(N1, aveKernel,'same');
aveN2 = conv2(N2, aveKernel,'same');
N1 = aveN1(1:F:rows,1:F:cols);
N2 = aveN2(1:F:rows,1:F:cols);

aveL1 = conv2(L1, aveKernel,'same');
aveL2 = conv2(L2, aveKernel,'same');
L1 = aveL1(1:F:rows,1:F:cols);
L2 = aveL2(1:F:rows,1:F:cols);

aveSM1 = conv2(saliencyMap1, aveKernel,'same');
aveSM2 = conv2(saliencyMap2, aveKernel,'same');
saliencyMap1 = aveSM1(1:F:rows,1:F:cols);
saliencyMap2 = aveSM2(1:F:rows,1:F:cols);
%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the gradient map
%%%%%%%%%%%%%%%%%%%%%%%%%
dx = [3 0 -3; 10 0 -10;  3  0 -3]/16;
dy = [3 10 3; 0  0   0; -3 -10 -3]/16;

IxL1 = conv2(L1, dx, 'same');     
IyL1 = conv2(L1, dy, 'same');    
gradientMap1 = sqrt(IxL1.^2 + IyL1.^2);

IxL2 = conv2(L2, dx, 'same');     
IyL2 = conv2(L2, dy, 'same');    
gradientMap2 = sqrt(IxL2.^2 + IyL2.^2);

%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate the VSI
%%%%%%%%%%%%%%%%%%%%%%%%%
VSSimMatrix = (2 * saliencyMap1 .* saliencyMap2 + constForVS) ./ (saliencyMap1.^2 + saliencyMap2.^2 + constForVS);
gradientSimMatrix = (2*gradientMap1.*gradientMap2 + constForGM) ./(gradientMap1.^2 + gradientMap2.^2 + constForGM);

weight = max(saliencyMap1, saliencyMap2);

ISimMatrix = (2 * M1 .* M2 + constForChrom) ./ (M1.^2 + M2.^2 + constForChrom);
QSimMatrix = (2 * N1 .* N2 + constForChrom) ./ (N1.^2 + N2.^2 + constForChrom);

SimMatrixC = (gradientSimMatrix .^ alpha) .* VSSimMatrix .* real((ISimMatrix .* QSimMatrix) .^ lambda) .* weight;
sim = sum(sum(SimMatrixC)) / sum(weight(:));

return;
%===================================
function VSMap = SDSP(image,sigmaF,omega0,sigmaD,sigmaC)
% ========================================================================
% SDSP algorithm for salient region detection from a given image.
% Copyright(c) 2013 Lin ZHANG, School of Software Engineering, Tongji
% University
% All Rights Reserved.
% ----------------------------------------------------------------------
% Permission to use, copy, or modify this software and its documentation
% for educational and research purposes only and without fee is here
% granted, provided that this copyright notice and the original authors'
% names appear on all copies and supporting documentation. This program
% shall not be used, rewritten, or adapted as the basis of a commercial
% software or hardware product without first obtaining permission of the
% authors. The authors make no representations about the suitability of
% this software for any purpose. It is provided "as is" without express
% or implied warranty.
%----------------------------------------------------------------------
%
% This is an implementation of the algorithm for calculating the
% SDSP (Saliency Detection by combining Simple Priors).
%
% Please refer to the following paper
%
% Lin Zhang, Zhongyi Gu, and Hongyu Li,"SDSP: a novel saliency detection 
% method by combining simple priors", ICIP, 2013.
%----------------------------------------------------------------------
%
%Input : image: an uint8 RGB image with dynamic range [0, 255] for each
%channel
%        
%Output: VSMap: the visual saliency map extracted by the SDSP algorithm.
%Data range for VSMap is [0, 255]. So, it can be regarded as a common
%gray-scale image.
%        
%-----------------------------------------------------------------------
%convert the image into LAB color space
[oriRows, oriCols, junk] = size(image);
image = double(image);
dsImage(:,:,1) = imresize(image(:,:,1), [256, 256],'bilinear');
dsImage(:,:,2) = imresize(image(:,:,2), [256, 256],'bilinear');
dsImage(:,:,3) = imresize(image(:,:,3), [256, 256],'bilinear');
lab = RGB2Lab(dsImage); 

LChannel = lab(:,:,1);
AChannel = lab(:,:,2);
BChannel = lab(:,:,3);

LFFT = fft2(double(LChannel));
AFFT = fft2(double(AChannel));
BFFT = fft2(double(BChannel));

[rows, cols, junk] = size(dsImage);
LG = logGabor(rows,cols,omega0,sigmaF);
FinalLResult = real(ifft2(LFFT.*LG));
FinalAResult = real(ifft2(AFFT.*LG));
FinalBResult = real(ifft2(BFFT.*LG));

SFMap = sqrt(FinalLResult.^2 + FinalAResult.^2 + FinalBResult.^2);

%the central areas will have a bias towards attention
coordinateMtx = zeros(rows, cols, 2);
coordinateMtx(:,:,1) = repmat((1:1:rows)', 1, cols);
coordinateMtx(:,:,2) = repmat(1:1:cols, rows, 1);

centerY = rows / 2;
centerX = cols / 2;
centerMtx(:,:,1) = ones(rows, cols) * centerY;
centerMtx(:,:,2) = ones(rows, cols) * centerX;
SDMap = exp(-sum((coordinateMtx - centerMtx).^2,3) / sigmaD^2);

%warm colors have a bias towards attention
maxA = max(AChannel(:));
minA = min(AChannel(:));
normalizedA = (AChannel - minA) / (maxA - minA);

maxB = max(BChannel(:));
minB = min(BChannel(:));
normalizedB = (BChannel - minB) / (maxB - minB);

labDistSquare = normalizedA.^2 + normalizedB.^2;
SCMap = 1 - exp(-labDistSquare / (sigmaC^2));
% VSMap = SFMap .* SDMap;
VSMap = SFMap .* SDMap .* SCMap;

VSMap =  imresize(VSMap, [oriRows, oriCols],'bilinear');
VSMap = mat2gray(VSMap);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function labImage = RGB2Lab(image)

image = double(image);
normalizedR = image(:,:,1) / 255;
normalizedG = image(:,:,2) / 255;
normalizedB = image(:,:,3) / 255;

RSmallerOrEqualto4045 = normalizedR <= 0.04045;
RGreaterThan4045 = 1 - RSmallerOrEqualto4045;
tmpR = (normalizedR / 12.92) .* RSmallerOrEqualto4045;
tmpR = tmpR + power((normalizedR + 0.055)/1.055,2.4) .* RGreaterThan4045;

GSmallerOrEqualto4045 = normalizedG <= 0.04045;
GGreaterThan4045 = 1 - GSmallerOrEqualto4045;
tmpG = (normalizedG / 12.92) .* GSmallerOrEqualto4045;
tmpG = tmpG + power((normalizedG + 0.055)/1.055,2.4) .* GGreaterThan4045;

BSmallerOrEqualto4045 = normalizedB <= 0.04045;
BGreaterThan4045 = 1 - BSmallerOrEqualto4045;
tmpB = (normalizedB / 12.92) .* BSmallerOrEqualto4045;
tmpB = tmpB + power((normalizedB + 0.055)/1.055,2.4) .* BGreaterThan4045;

X = tmpR*0.4124564 + tmpG*0.3575761 + tmpB*0.1804375;
Y = tmpR*0.2126729 + tmpG*0.7151522 + tmpB*0.0721750;
Z = tmpR*0.0193339 + tmpG*0.1191920 + tmpB*0.9503041;

epsilon = 0.008856; %actual CIE standard
kappa   = 903.3; %actual CIE standard
Xr = 0.9642; %reference white D50
Yr = 1.0; %reference white
Zr = 0.8251; %reference white

xr = X/Xr;
yr = Y/Yr;
zr = Z/Zr;

xrGreaterThanEpsilon = xr > epsilon;
xrSmallerOrEqualtoEpsilon = 1 - xrGreaterThanEpsilon;
fx = power(xr, 1.0/3.0) .* xrGreaterThanEpsilon;
fx = fx + (kappa*xr + 16.0)/116.0 .* xrSmallerOrEqualtoEpsilon;

yrGreaterThanEpsilon = yr > epsilon;
yrSmallerOrEqualtoEpsilon = 1 - yrGreaterThanEpsilon;
fy = power(yr, 1.0/3.0) .* yrGreaterThanEpsilon;
fy = fy + (kappa*yr + 16.0)/116.0 .* yrSmallerOrEqualtoEpsilon;

zrGreaterThanEpsilon = zr > epsilon;
zrSmallerOrEqualtoEpsilon = 1 - zrGreaterThanEpsilon;
fz = power(zr, 1.0/3.0) .* zrGreaterThanEpsilon;
fz = fz + (kappa*zr + 16.0)/116.0 .* zrSmallerOrEqualtoEpsilon;

[rows,cols,junk] = size(image);
labImage = zeros(rows,cols,3);
labImage(:,:,1) = 116.0 * fy - 16.0;
labImage(:,:,2) = 500.0 * (fx - fy);
labImage(:,:,3) = 200.0 * (fy - fz);
return;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function LG = logGabor(rows,cols,omega0,sigmaF)
     [u1, u2] = meshgrid(([1:cols]-(fix(cols/2)+1))/(cols-mod(cols,2)), ...
            ([1:rows]-(fix(rows/2)+1))/(rows-mod(rows,2)));
     mask = ones(rows, cols);
     for rowIndex = 1:rows
         for colIndex = 1:cols
             if u1(rowIndex, colIndex)^2 + u2(rowIndex, colIndex)^2 > 0.25
                 mask(rowIndex, colIndex) = 0;
             end
         end
     end
     u1 = u1 .* mask;
     u2 = u2 .* mask;
     
     u1 = ifftshift(u1);  
     u2 = ifftshift(u2);
     
     radius = sqrt(u1.^2 + u2.^2);    
     radius(1,1) = 1;
            
     LG = exp((-(log(radius/omega0)).^2) / (2 * (sigmaF^2)));  
     LG(1,1) = 0; 
return;

Comments

Popular posts from this blog

Image RGB to CMYK : Image Conversion - MATLAB Code (RGB to CMYK, CMYK to RGB)

MATAB use makecform and applycform functions to perform this color model conversion  Ref: makecform:  https://www.mathworks.com/help/images/ref/makecform.html Ref: applycform:  https://www.mathworks.com/help/images/ref/applycform.html RGB to CMYK: I_rgb = imread('lena_color.jpg'); %Read the color image  C = makecform('srgb2cmyk');     %srgb (standard rgb) to cmyk  I_CMYK= applycform(I_rgb,C); TO SAVE IN CMYK Format imwrite(I_CMYK, 'test.tiff'); USE IMFINFO to confirm the conversion imfinfo('test.tiff') the command gives the output ... Width: 512 Height: 512 BitDepth: 32 ColorType: 'CMYK' ... CMYK to RGB: I_CMYK = imread('lena_color.jpg'); %Read the color image  C = makecform('cmyk2srgb');     %srgb (standard rgb) to cmyk  I_rgb= applycform(I_CMYK, C);

Digital Half Toning - Ordered Dithering - MATLAB Code Bayer/ Ulichney

Digital Half-toning is a technique to convert the gray scale / color image (in the range 0-255) to binary images (in the range 0-1)  that is useful for printing (especially black and white printers). Though at first glance dots are clearly visible, as Human Vision System is a low pass filter, when dots are printed closely and viewed from the specific distance it can be perceived as a continuous tone image (i.e. gray scale). In general, half toning techniques is widely used in newspapers, journals and magazines printing.  Halftoning algorithms aims for better #Image Quality #Low Complexity #Optimal ink usage Coming to the topic, Ordered dithering is one of the standard technique to produce halftone image from continuous (gray scale) images. Let us assume, the first block containing 'B' is an image matrix and the  'H' is the threshold matrix. And the process is simply a element wise (pixel wise) comparison; if element of B > H (g...

Direct Binary Search Halftoning (DBS) - MATLAB Code (Efficient Direct Binary Search, Halftoning)

Direct binary approach is a heuristic optimization method and is proven to be very powerful to  obtain optimized binary patterns. The approach achieves the least minimum square error between the perceived halftone and original image, though swap and toggle operation. The swap operation consists of switching the current pixel with eight of its neighborhood pixels. Toggle operation is about switching the values between 0 to 1 or vice versa (as shown in Fig. below)      With each iteration, the perceived error starts to reduce and in finally a superior halftone quality is achieved. The method is computationally very expensive and difficult to implement in hardware. The code for Efficient DBS is provided in the link: (MATLAB Code) https://github.com/SankarSrin/EDBS-Halftoning