原文来自Perra C , Massidda F , Giusto D D . Image blockiness evaluation based on Sobel operator[C]// IEEE International Conference on Image Processing. IEEE, 2005.
摘要：该文章提出了一种新颖的无参考图像块效应评价方法，该方法利用Sobel算子，块效应分数是图像块边缘像素及中心像素的亮度变化的组合。 详细方法：
对图像I和Sobel算子的蒙版Mx和My分别卷积： Mx=[-1,0,1;-2,0,2;-1,0,1] My=[-1,-2,-1;0,0,0;1,2,1] 得到边缘图像Dx和Dy: Dx=IMx,Dy=IMy 其中"*"代表卷积操作符
将Dx和Dy分解为若干个不重叠的N*N大小的块，当块效应为Jpeg压缩产生时，N为8。
对于每个N*N的块(Dx和Dy分解后的块)，做如下分割： 从左到右的区域（图中灰色部分）依次记作Ω1V、Ω1H和Ω2.即每个图像块中位于第一列和最后一列的所有像素属于Ω1V，第一行和最后一行的所有像素属于Ω1H，满足位于第二行到倒数第二行、第二列到倒数第二列的所有像素属于Ω2.
s1和s2定义如下： 其中N1和N2分别为右侧和式的长度，也就是外圈的像素数和内部的像素数。
块效应指标S： β取2.
Matlab实现：
首先将图像I转为灰度图像，若I已经是灰度图像则不操作： if numel(size(I))>2 %如果不为灰度图像 I=rgb2gray(I); end 通过conv2函数或者edge获得水平、垂直边缘图像dx、dy以及总边缘图像d。为方便之后的操作，将三个边缘图像转为dobule类型。 dx=double(edge(I,'Sobel','horizontal')); dy=double(edge(I,'Sobel','vertical')); d=double(edge(I,'Sobel')); 分块：设lenx和leny为块的高和宽，默认设置为8，num表示为已经分块的个数， 初始值为0。水平方向以leny个单位为间隔、垂直以lenx个单位为间隔进行两层遍历，遍历终点分别为图像的宽(size(I,2))和高(size(I,1))。对于垂直和水平遍历当前值i和j，该块的左上角像素在原始图像中的位置就是(i,j)，右下角像素为(xend,yend)，xend为高和i+lenx-1的最小值，即设置为该块的高最大为lenx，最小的话取决于到底部的距离(因为NN分块，最底部和最右部的块不一定满足NN大小)，yend同理，把对dx、dy、d分块后的结果存在结构体Block内。两层遍历结束时，Block的元素number设置为变量num,表示为块数。 lenx=8; leny=8; num=0; for i=1:lenx:size(I,1) for j=1:leny:size(I,2) xend=min(i+lenx-1,size(I,1)); yend=min(j+leny-1,size(I,2)); num=num+1; Block.images_dx{num}=dx(i:xend,j:yend); Block.images_dy{num}=dy(i:xend,j:yend); Block.images_d{num}=d(i:xend,j:yend); end end Block.number=num;4.获得水平、垂直边缘、内部的像素：以水平边缘像素为例，编写函数get_h()返回图像块的水平边缘像素在该块内的坐标。首先获得图像块的大小m*n,如果m大于等于2，则返回的是两行(第一行和最后一行)，设置horizontal矩阵为2n（两行像素个数）2，每一行为一个二元坐标。带如果为1说明只有一行，则返回这一行，horizontal设置为n2。 首先i遍历1到n，把坐标(1,i)存入horizontal的第i行，如果m大于等于2，则继续存入：再次i遍历1到n，把坐标 （n,i）存入horizontal的第i+n行。
function [horizontal]=get_h(dy) [m,n]=size(dy); if m>=2 horizontal=zeros(2*n,2); else horizontal=zeros(n,2); end for i=1:n horizontal(i,1)=1; horizontal(i,2)=i; end if m>=2 for i=1:n horizontal(i+n,1)=n; horizontal(i+n,2)=i; end end end垂直边缘像素同理编写函数get_v()。对于内部像素，获得图像块的大小mn,由于最上下左右侧的像素不计入，所以内部像素个数为(m-2)(n-2)。垂直方向和水平方向分别i，j从2遍历到m-1（或n-1），则结果的第k行为(i,j)，k为目前的坐标数，初始化为0。
function [inner]=get_inner(d) [m,n]=size(d); inner=zeros((m-2)*(n-2),2); k=0; for i=2:m-1 for j=2:n-1 k=k+1; inner(k,1)=i; inner(k,2)=j; end end end i从1遍历到Block.number，分别取当前第i个dx、dy、d图像块，即Block.images_dx{i}，其他两个同理，调入之前编写的函数，获得存储垂直、水平和内部像素坐标的矩阵，记作v、h、in。则N1就是v的高和h的高之和（v和h两个矩阵所存储的点数值和），N2为in的高(in矩阵存储的点数)。对于当前块而言，初始化s1、s2为0，取max_x、max_y、max_in为当前块Block.images_dx{i}、Block.images_dy{i}、Block.images_d{i}的最大值。j遍历v的所有行，根据公式,s1累加abs(Block.images_dx{i}(v(j,1),v(j,2)))/max_x，同理遍历h、in。最后将s1的值除以N1保存在Block的s1矩阵中，s2的值除以N2保存在s2矩阵中。 for i=1:Block.number v=get_v(Block.images_dx{i}); h=get_h(Block.images_dy{i}); in=get_inner(Block.images_d{i}); N1=size(v,1)+size(h,1); N2=size(in,1); s1=0; s2=0; max_x=max(max(Block.images_dx{i})); max_y=max(max(Block.images_dy{i})); max_in=max(max(Block.images_d{i})); for j=1:size(v,1) s1=s1+abs(Block.images_dx{i}(v(j,1),v(j,2)))/max_x; end for j=1:size(h,1) s1=s1+abs(Block.images_dy{i}(h(j,1),h(j,2)))/max_y; end for j=1:N2 s2=s2+abs(Block.images_d{i}(in(j,1),in(j,2)))/max_in; end Block.s1(i)=s1/N1; Block.s2(i)=s2/N2; end 计算S：首先，由于部分块的N1或N2为0，Block的s1、s2矩阵中存在NaN值，首先对其中NaN的值替换为0：Block.s1(isnan(Block.s1))=0;然后求s1、s2的均值S1、S2，最后根据公式求S。 Block.s1(isnan(Block.s1))=0; Block.s2(isnan(Block.s2))=0; S1=mean(Block.s1); S2=mean(Block.s2); lambda=2; S=abs((S1^lambda-S2^lambda)/(S1^lambda+S2^lambda));完整代码：
I=imread('1.jpg'); if numel(size(I))>2 I=rgb2gray(I); end dx=double(edge(I,'Sobel','horizontal')); dy=double(edge(I,'Sobel','vertical')); d=double(edge(I,'Sobel')); lenx=8; leny=8; num=0; for i=1:lenx:size(I,1) for j=1:leny:size(I,2) xend=min(i+lenx-1,size(I,1)); yend=min(j+leny-1,size(I,2)); num=num+1; Block.images_dx{num}=dx(i:xend,j:yend); Block.images_dy{num}=dy(i:xend,j:yend); Block.images_d{num}=d(i:xend,j:yend); end end Block.number=num; for i=1:Block.number v=get_v(Block.images_dx{i}); h=get_h(Block.images_dy{i}); in=get_inner(Block.images_d{i}); N1=size(v,1)+size(h,1); N2=size(in,1); s1=0; s2=0; max_x=max(max(Block.images_dx{i})); max_y=max(max(Block.images_dy{i})); max_in=max(max(Block.images_d{i})); for j=1:size(v,1) s1=s1+abs(Block.images_dx{i}(v(j,1),v(j,2)))/max_x; end for j=1:size(h,1) s1=s1+abs(Block.images_dy{i}(h(j,1),h(j,2)))/max_y; end for j=1:N2 s2=s2+abs(Block.images_d{i}(in(j,1),in(j,2)))/max_in; end Block.s1(i)=s1/N1; Block.s2(i)=s2/N2; end Block.s1(isnan(Block.s1))=0; Block.s2(isnan(Block.s2))=0; S1=mean(Block.s1); S2=mean(Block.s2); lambda=2; S=abs((S1^lambda-S2^lambda)/(S1^lambda+S2^lambda)); function [vertical]=get_v(dx) [m,n]=size(dx); if n>=2 vertical=zeros(2*m,2); else vertical=zeros(m,2); end for i=1:m vertical(i,1)=i; vertical(i,2)=1; end if n>=2 for i=1:m vertical(i+m,1)=i; vertical(i+m,2)=n; end end end function [horizontal]=get_h(dy) [m,n]=size(dy); if m>=2 horizontal=zeros(2*n,2); else horizontal=zeros(n,2); end for i=1:n horizontal(i,1)=1; horizontal(i,2)=i; end if m>=2 for i=1:n horizontal(i+n,1)=n; horizontal(i+n,2)=i; end end end function [inner]=get_inner(d) [m,n]=size(d); inner=zeros((m-2)*(n-2),2); k=0; for i=2:m-1 for j=2:n-1 k=k+1; inner(k,1)=i; inner(k,2)=j; end end end