图像质量评估系统的Matlab GUIDE实现(四)——相关算法

系统的大体流程和介绍在（一）~（三）中大致介绍完了。

该系统采用的方法主要有两种，一种是全参考（RF），另一种是无参考（NF）

贴出的算法代码是用于系统调用的返回值版：

一、全参考（RF）

RF中主要包括了MSE\PSNR\JND\PDM\SSIM\SIExt\QILV\VIF几个算法,在这里就不逐个介绍了。

(1)MSE(Mean Square Error)均方误差

该算法直接比较得到两幅图像 I、J 每一个像素位置对应的绝对差值，最后总平均差值即为 MSE。

(2)PSNR(Peak Signal to Noise Ratio)峰值信噪比

可由MSE值计算得到。

function [PSNR, MSE] = mseandpsnr(X, Y)
D = X-Y;
[m n d]=size(D);
MSE = sum(D(:).^2)/(m*n*d);
PSNR = 10*log10(258^2/MSE);

(3)JND(Just Noticeable Difference)最小可觉差/差别感知阈限

这个在网上找到了一个版本的代码，但是无法搬到系统中，搜了半天看到有说这个算法基于mallat，于是就用mallat替代了。

%离散快速正交小波变换——Mallat算法
function [psnr]=mallat(I)
    order = 3; %变换级数为3
    I=rgb2gray(I);
    I=imresize(I,[256,256]);
    wname = 'db2';
    ordata=double(I);
    data=ordata;
    [row,col]=size(data);
    imdata=zeros(size(data,1),size(data,2));
    %分解
    [lo_d,hi_d]=wfilters(wname,'d');
    for i=1:order
        data=decomp(lo_d,hi_d,i,data);
        MIN=min(min(data));
        MAX=max(max(data));
        for j=1:size(data,1)
            for k=1:size(data,2)
                imdata(j,k)=((data(j,k)-MIN)/(MAX-MIN))*255;
            end
        end
        %figure;
        %imshow(uint8(imdata));    
    end

    %计算最后一级分解后能量分布
    energy=zeros(1,4);
    energy=cal_energy(data,order);

    %阈值化处理
    [thr,sorh]=ddencmp('den','wv',ordata);  %获得阈值
    ddata=wthresh(data,sorh,thr);           %阈值化处理
    MIN=min(min(ddata));
    MAX=max(max(ddata));
        for j=1:size(ddata,1)
            for k=1:size(ddata,2)
                imdata(j,k)=((ddata(j,k)-MIN)/(MAX-MIN))*255;
            end
        end
    %figure;    
    %imshow(uint8(imdata));  

    %统计0的个数
    sum=0;
    for i=1:row
        for j=1:col
            if (ddata(i,j)==0)
                sum=sum+1;
            end
        end
    end
    Zero=sum/(row*col);

    %重建
    [lo_r,hi_r]=wfilters(wname,'r');
    redata=zeros(size(data,1),size(data,2));
    data=ddata;
    for i=order:-1:1
        data=reconstrn(lo_r,hi_r,i,data);
        MIN=min(min(data));
        MAX=max(max(data));
        for j=1:size(data,1)
            for k=1:size(data,2)
                redata(j,k)=((data(j,k)-MIN)/(MAX-MIN))*255;
            end
        end
        %figure;
        %imshow(uint8(redata));
    end

    %计算重构后的均值和方差
    [mean,var]=cal_mv(data);

    %峰值信噪比PSNR
    [row,col]=size(data);
    mse=0;
    for i=1:row
        for j=1:col
            mse=mse+(data(i,j)-ordata(i,j))^2;
        end
    end
    mse=mse/(row*col);
    psnr=10*log10(255^2/mse);
end

调用到的其他函数：

cal_mv.m

%求均值、方差
function [mean var]=cal_mv(data)
[row col]=size(data);
sum=0;
for i=1:row
    for j=1:col
        sum=sum+data(i,j);
    end
end
mean = sum/(row*col);
sum=0;
for i=1:row
    for j=1:col
        sum=sum+(data(i,j)-mean)^2;
    end
end
var = sum/(row*col);
        

        

decomp.m 

%分解变换
function data_d=decomp(lo_d,hi_d,order,data)
[row,col]=size(data);
%一维行变换
firowlen=size(lo_d,2);
l=fix(firowlen/2);
%excoldata=wextend('ac','sym',data(1:row/2^(order-1),1:col/2^(order-1)),l);
excoldata=wextend('ac','ppd',data(1:row/2^(order-1),1:col/2^(order-1)),firowlen-1,'r');%往右周期延拓
for i=1:row/2^(order-1);
    for j=1:row/2^order;
        data(i,j)=lo_d*excoldata(i,(2*j-1):(2*j+firowlen-2))';
        data(i,j+row/2^order)=hi_d*excoldata(i,(2*j-1):(2*j+firowlen-2))';
    end
end

%一维列变换
ficollen=size(hi_d,2);
l=fix(ficollen/2);
%exrowdata=wextend('ar','sym',data(1:row/2^(order-1),1:col/2^(order-1)),l);
exrowdata=wextend('ar','ppd',data(1:row/2^(order-1),1:col/2^(order-1)),ficollen-1,'d');%往下周期延拓
for i=1:col/2^(order-1);
    for j=1:col/2^order;
        data(j,i)=lo_d*exrowdata((2*j-1):(2*j+ficollen-2),i);
        data(j+col/2^order,i)=hi_d*exrowdata((2*j-1):(2*j+ficollen-2),i);
    end
end
data_d= data;

reconstrn.m 

%重建
function data_r=reconstrn(lo_r,hi_r,order,data)
[row,col]=size(data);
lo_ro=lo_r(1:2:end);%提取lo_r奇数位上的值
lo_re=lo_r(2:2:end);%提取lo_r偶数位上的值
hi_ro=hi_r(1:2:end);%提取hi_r奇数位上的值
hi_re=hi_r(2:2:end);%提取hi_r偶数位上的值
leno=size(lo_ro,2);
lene=size(lo_re,2);
%列重构
exloodata=wextend('ar','ppd',data(1:row/2^order,1:col/2^(order-1)),(leno-1),'u');%向上周期延拓
exloedata=wextend('ar','ppd',data(1:row/2^order,1:col/2^(order-1)),(lene-1),'u');%向上周期延拓
exhiodata=wextend('ar','ppd',data((row/2^order+1):row/2^(order-1),1:col/2^(order-1)),(leno-1),'u');%向上周期延拓
exhiedata=wextend('ar','ppd',data((row/2^order+1):row/2^(order-1),1:col/2^(order-1)),(lene-1),'u');%向上周期延拓
for i=1:row/2^(order-1);
    for j=1:col/2^(order-1);
        k=fix((i+1)/2);
        if(mod(i,2)==0)
            lodata=lo_ro*exloodata(k:k+leno-1,j);
            hidata=hi_ro*exhiodata(k:k+leno-1,j);
            data(i,j)=lodata+hidata;
        else
            lodata=lo_re*exloedata(k:k+lene-1,j);
            hidata=hi_re*exhiedata(k:k+lene-1,j);
            data(i,j)=lodata+hidata;
        end
    end
end

%行重构
exloodata=wextend('ac','ppd',data(1:row/2^(order-1),1:col/2^order),(leno-1),'l');%向左周期延拓
exloedata=wextend('ac','ppd',data(1:row/2^(order-1),1:col/2^order),(lene-1),'l');%向左周期延拓
exhiodata=wextend('ac','ppd',data(1:row/2^(order-1),(col/2^order+