Quality as an image-specific characteristic perceived by an average human observer

end

bm_v = 1/sum(prof_vert(1:n-1))*sqrt((1/(bl-1))*sum_FPH_vert);

blockness = sqrt((bm_v^2)*0.5 + (bm_h^2)*0.5);

end

BlurOne

% BLUR METRICS

%

% The Blur Effect: Perception and Estimation with a New No-Reference

% Perceptual Blur Metric

%

% metric NO 001

function blur = blurOne(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

% lets create blurred image B:

hv = [1 1 1 1 1 1 1 1 1] * 1/9;

hh = hv';

Bvertical = imfilter(imageF, hv);

Bhorizontal = imfilter(imageF, hh);

% compute absolute difference images of B and F:

diff_Fvertical = abs(imageF(2:m, :) - imageF(1:m-1, :));

diff_Fhorizontal = abs(imageF(:, 2:n) - imageF(:, 1:n-1));

diff_Bvertical = abs(Bvertical(2:m, :) - Bvertical(1:m-1, :));

diff_Bhorizontal = abs(Bhorizontal(:, 2:n) - Bhorizontal(:, 1:n-1));

% compute variation

var_vertical = max(0, diff_Fvertical - diff_Bvertical);

var_horizontal = max(0, diff_Fhorizontal - diff_Bhorizontal);

% compute sum of coefficients of differences (sum of all matrix elements):

sum_Fvertical = sum(sum(diff_Fvertical(2:m-1, 2:n-1)));

sum_Fhorizontal = sum(sum(diff_Fhorizontal(2:m-1, 2:n-1)));

sum_var_vertical = sum(sum(var_vertical(2:m-1, 2:n-1)));

sum_var_horizontal = sum(sum(var_horizontal(2:m-1, 2:n-1)));

% normalize the result:

norm_Fvertical = (sum_Fvertical - sum_var_vertical)/sum_Fvertical;

norm_Fhorizontal = (sum_Fhorizontal - sum_var_horizontal)/sum_Fhorizontal;

blur = max(norm_Fvertical, norm_Fhorizontal);

end

BlurTwo

% BLUR METRICS

%

% A No-Reference Blur Image Quality measure

% Based on Wavelet Transform

%

% Min Goo Choi method

% metric NO 002

function blur = blurTwo(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

% treshold = 0.1

% compute absolute horizontal difference of image F (HADV):

diff_horizontal = abs(imageF(:, 1:n-2) - imageF(:, 3:n));

sum_horizontal = sum(sum(diff_horizontal));

diff_hor_mean = sum_horizontal /(m*(n-2));

% find edge candidates:

Ch = zeros(m,n-2);

for i=1:m

for j=1:n-2

if diff_horizontal(i,j) > diff_hor_mean

Ch(i,j) = diff_horizontal(i,j);

else

Ch(i,j) = 0;

end

end

end

Eh = zeros(m, n-2);

for i=1:m

for j=2:n-3

if (Ch(i,j) > Ch(i, j-1)) & (Ch(i,j) > Ch(i, j+1))

Eh(i,j) = 1;

else

Eh(i,j) = 0;

end

end

end

% compute for vertical:

diff_vertical = abs(imageF(1:m-2, :) - imageF(3:m, :));

sum_vertical = sum(sum(diff_vertical));

diff_vert_mean = sum_vertical /(m*(n-2));

Cv = zeros(m-2,n);

for i=1:m-2

for j=1:n

if diff_vertical(i,j) > diff_vert_mean

Cv(i,j) = diff_vertical(i,j);

else

Cv(i,j) = 0;

end

end

end

Ev = zeros(m-2, n);

for i=2:m-3

for j=1:n

if (Cv(i,j) > Cv(i-1, j)) & (Cv(i,j) > Cv(i+1, j))

Ev(i,j) = 1;

else

Ev(i,j) = 0;

end

end

end

% does detected edge pixel correspond to a blurred edge?

A_vertical = imageF;

for j=1:n

for i=2:m-1

A_vertical(i,j) = (abs(imageF(i+1, j) + imageF(i-1, j)))/2.0;

end

for i = 1

A_vertical(i,j) = A_vertical(i,j);

end

for i = m

A_vertical(i,j) = A_vertical(i,j);

end

end

BR_vertical = (abs(imageF(:,:) - A_vertical(:,:)))./max(A_vertical(:,:),zeros(m,n) + 1e-10);

A_horizontal = imageF;

for i=1:m

for j=2:n-1

A_horizontal(i,j) = (abs(imageF(i, j-1) + imageF(i, j+1)))/2.0;

end

for j = 1

A_horizontal(i,j) = A_horizontal(i,j);

end

for j = n

A_horizontal(i,j) = A_horizontal(i,j);

end

end

BR_horizontal = (abs(imageF(:,:) - A_horizontal(:,:)))./max(A_horizontal(:,:),zeros(m,n) + 1e-10) ;

B = zeros(m, n);

for i=1:m

for j=1:n

if max(BR_vertical(i,j), BR_horizontal(i,j)) < .1;

B(i,j) = 1;

else

B(i,j) = 0;

end

end

end

% how many edge and blurred pixels?

blur_cnt_h = nnz(B(:, 2:n-1) .* Eh);

blur_cnt_v = nnz(B(2:m-1, :) .* Ev);

edge_cnt_h = nnz(Eh);

edge_cnt_v = nnz(Ev);

ratio_h = blur_cnt_h/edge_cnt_h;

ratio_v = blur_cnt_v/edge_cnt_v;

blur = max(ratio_v, ratio_h);

end

Contrast

function contr = contrast(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

kernel = [-1, -1, -1, -1, 8, -1, -1, -1]/8;

diffImage = convn(imageF, kernel, 'same');

cpp = mean2(diffImage);

contr = cpp;

Contrast 2

function contr2 = contrast2(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

image = imageF(:); % make it vector

m = mean(image);

image = image - m;

contr2=norm(image)/sqrt(numel(image));

EI

function EI = edgeIntens(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

Gx = zeros(m,n);

Gy = Gx;

F = Gy;

for i=2:m-1

for j=2:n-1

a1 = imageF(i+1,j-1);

a2 = 2*imageF(i+1,j);

a3 = imageF(i+1, j+1);

sum1 = abs(a1+a2+a3);

a4 = imageF(i-1,j);

a5 = 2*imageF(i-1,j);

a6 =imageF(i-1,j+1);

sum2=abs(a4+a5+a6);

Gx(i,j) = sum1-sum2;

b1 = imageF(i-1,j+1);

b2 = 2*imageF(i,j+1);

b3 = imageF(i+1, j+1);

bum1 = abs(a1+a2+a3);

b4 = imageF(i-1,j-1);

b5 = 2*imageF(i,j-1);

b6 =imageF(i+1,j-1);

bum2=abs(a4+a5+a6);

Gy(i,j) = bum1-bum2;

F(i,j) = abs(Gx(i,j)^2 + Gy(i,j)^2);

end

end

EI = sum(F(:))/(m*n);

Fractal

% Fractal dimension

function fractal = fractal(imageF)

imageF = im2double(imageF);

[m, n] = size(imageF);

% const

e= 0.04;

%edge pixels

edge_x = abs(imageF(1:m-2, :) - imageF(3:m, :));

max_edge_x = max(edge_x);

edge_y = abs(imageF(:, 1:n-2) - imageF(:, 3:n));

max_edge_y = max(edge_y);

edge_pix_x= zeros(m,n);

for i=2:m-1

for j=1:n

if edge_x(i-1,j) > (max_edge_x*e)

edge_pix_x(i,j) = 1;

else

edge_pix_x(i,j) = 0;

end

end

end

edge_pix_y= zeros(m,n);

for j=2:n-1

for i=1:m

if edge_y(i,j-1) > (max_edge_y*e)

edge_pix_y(i,j) = 1;

else

edge_pix_y(i,j) = 0;

end

end

end

[n,r] = boxcount(edge_pix_x);

s = -gradient(log(n))./gradient(log(r));

fractal = mean(s(2:5));

end

Separability

% Segmentation

% object separability

%segmentation by edge

function separability = separability(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

medianF = mean(imageF(:));

% foreground

for i=1:m

for j=1:n

if imageF(i,j) <= medianF

imageF(i,j) = 0;

else imageF(i,j) = imageF(i,j);

end

end

end

imageB = im2double(image);

% background

for i=1:m

for j=1:n

if imageB(i,j) >= medianF

imageB(i,j) = 0;

else imageB(i,j) = imageB(i,j);

end

end

end

%size of segment

[ff,e] = size(nonzeros(imageF));

[bb,e] = size(nonzeros(imageB));

[m, n] = size(imageF);

%avg diff for central pixel in segment

Uk = zeros(m,n);

for i=2:m-1

for j=2:n-1

a1 = imageF(i,j);

a2 = imageF(i-1,j-1);

a3 = imageF(i-1, j);

a4 = imageF(i-1,j+1);

a5 = imageF(i,j+1);

a6 =imageF(i,j-1);

a7 = imageF(i+1,j-1);

a8 = imageF(i+1, j);

a9 =imageF(i+1,j+1);

if a1 > 0

vec = [a2,a3,a4,a5,a6,a7,a8,a9];

ind = find(vec>0);

[m,n]=size(ind);

pix = zeros(1,n);

for k=1:n

pix(1,k) = (vec(ind(k))-a1)^2;

end

Uk(i,j) = (sum(pix))/8;

else

continue;

end

end

end

%avg diff for central pixel in segment

Uk2 = zeros(m,n);

for i=2:m-1

for j=2:n-1

a1 = imageB(i,j);

a2 = imageB(i-1,j-1);

a3 = imageB(i-1, j);

a4 = imageB(i-1,j+1);

a5 = imageB(i,j+1);

a6 =imageB(i,j-1);

a7 = imageB(i+1,j-1);

a8 = imageB(i+1, j);

a9 =imageB(i+1,j+1);

if a1 > 0

vec = [a2,a3,a4,a5,a6,a7,a8,a9];

ind = find(vec>0);

[m,n]=size(ind);

pix = zeros(1,n);

for k=1:n

pix(1,k) = (vec(ind(k))-a1)^2;

end

Uk2(i,j) = (sum(pix))/8;

else

continue;

end

end

end

PkF = sum(Uk(:))/ff;

PkB = sum(

Uk2(:))/bb;

W = (PkF+PkB)/2;

%average pixel intensity in segment

IF = (sum(nonzeros(imageF)))/ff;

IB = (sum(nonzeros(imageB)))/bb;

G = (IF-IB)^2;

B = 1/G;

IM = 1000*W+B;

separability = IM;

end

Sharpness

% Sharpness

% Sharpness Estimation for Document and Scene Images

% metric NO 005

function sharpness1 = sharpness1(image)

% original image F:

% width of block

w = 10;

% sharpness threshold

p = 2;

% edge threshold

e = 0.09;

%imageF = im2double(image);

imageF = rgb2gray(image);

[m, n] = size(imageF);

% median filter for image

image_med = medfilt2(double(imageF), [3,3]);

% DoM:diff of diff for median filtered image

diff_med_vertical = abs(abs(image_med(3:m-2, :) - image_med(1:m-4, :)) - abs(image_med(3:m-2,:) - image_med(5:m,:)));

diff_med_horizontal = abs(abs(image_med(:, 3:n-2) - image_med(:,1:n-4)) - abs(image_med(:, 3:n-2) - image_med(:, 5:n)));

% diff for original image

diff_vertical = abs(image_med(1:m-1, :) - image_med(2:m, :));

diff_horizontal = abs(image_med(:, 1:n-1) - image_med(:, 2:n));

s_x = zeros(m,n);

for i = 4+w:m-3-w

for j=1+w:n-w

podmatrix = diff_med_vertical(i-w-2:i+w-2, j);

up = sum(podmatrix(:));

podmatrix_original = diff_vertical(i-w-1:i+w-1, j);

down = sum(podmatrix_original(:));

if down == 0

s_x(i,j) = 1;

else

s_x(i,j) = up/(down+1e-10);

end

end

end

s_y = zeros(m,n);

for j = 4+w:n-3-w

for i=1+w:m-w

podmatrix = diff_med_horizontal(i, j-w-2:j+w-2);

up = sum(podmatrix(:));

podmatrix_original = diff_horizontal(i, j-w-1:j+w-1);

down = sum(podmatrix_original(:));

if down == 0

s_y(i,j) = 1;

else

s_y(i,j) = up/(down+1e-10);

end

end

end

%edge pixels

edge_x = abs(imageF(1:m-2, :) - imageF(3:m, :));

max_edge_x = max(edge_x);

edge_y = abs(imageF(:, 1:n-2) - imageF(:, 3:n));

max_edge_y = max(edge_y);

edge_pix_x= zeros(m,n);

for i=2:m-1

for j=1:n

if edge_x(i-1,j) > (max_edge_x*e)

edge_pix_x(i,j) = 1;

else

edge_pix_x(i,j) = 0;

end

end

end

edge_pix_y= zeros(m,n);

for j=2:n-1

for i=1:m

if edge_y(i,j-1) > (max_edge_y*e)

edge_pix_y(i,j) = 1;

else

edge_pix_y(i,j) = 0;

end

end

end

matrix_edge_sharp_x = zeros(m,n);

for i = 4+w:m-3-w

for j=1+w:n-w

if edge_pix_x(i,j) == 1 && s_x(i,j) > p

matrix_edge_sharp_x(i,j) = 1;

end

end

end

matrix_edge_sharp_y = zeros(m,n);

for j = 4+w:n-3-w

for i=1+w:m-w

if edge_pix_y(i,j) == 1 && s_y(i,j) > p

matrix_edge_sharp_y(i,j) = 1;

end

end

end

R_x = sum(matrix_edge_sharp_x(:))/sum(edge_pix_x(:));

R_y = sum(matrix_edge_sharp_y(:))/sum(edge_pix_y(:));

sharpness1 = sqrt(R_x^2 + R_y^2);

%sharpness1 = imshow(matrix_edge_sharp_x);

%sharpness1= imshow(s_x);

end

Spectral

% Spectral flatness

% No-Reference Image Quality Metrics for Structural MRI

% metric NO 004

function spec = spectral(image)

% original image F:

imageF = im2double(image);

[m, n] = size(imageF);

medianF = mean(imageF(:));

% background

for i=1:m

for j=1:n

if imageF(i,j) <= medianF

imageF(i,j) = 0;

else imageF(i,j) = imageF(i,j);

end

end

end

%image as 2D signal

% 2D descrete Furier transform

F= fft2(imageF);

% reshape

F_v1 = reshape(F, m*n, 1);

function s = a(x)

s = abs(x)^2;

end

F_v = arrayfun(@a, F_v1(:));

geom_mean = geomean(F_v);

arith_mean = mean2(F_v);

%spectral flatness quality measure

FF = geom_mean / arith_mean;

%spatial flatness quality measure

I_v1 = reshape(imageF, m*n, 1);

I_v = arrayfun(@a, I_v1(:));

geom_mean = geomean(I_v);

arith_mean = mean2(I_v);

II = geom_mean / arith_mean;

%final measure

var_I = var(imageF(:));

FE = FF * var_I;

spec = FE;

end

C. Finding best linear regression models

function [cBest, fitBest, CombsF] = OM_BestRegression(Y, X, nVars, nSets, nRegressionType)

% Find the best subset of predictor variables X

% containing exactly nVars variables

% to predict the dependent variable Y

% nSets - the number of best regression parameter combinations

% nRegressionType - defines different regression types

% to return in CombsF

Y_trees = Y(1:50,:);

X_trees = X(1:50,:);

Y_med = Y(50:104,:);

X_med = X(50:104,:);

% Verify input

if nargin < 4

nSets = 20;

end

if nargin < 5

nRegressionType = 2; % least squares

end

% Define residual function for L1 regression

NormL1=@(r)mean(abs(r));

NormL2=@(r)std(r);

NormLM=@(r)median(abs(r));

opts = optimset('MaxIter', 10000000, 'MaxFunEvals', 1000000);

% Find the sizes of independent variables X

[~, nx] = size(X);

% Initialize all possible combinations of nVars variables

Combs = FindBestCombs(Y_trees,X_trees,nVars,50*nVars);

Combs2 = FindBestCombs(Y_med,X_med,nVars,50*nVars);

%Combs = [Combs1;Combs2]

nAllCombs = length(Combs);

% If we have too many combinations, branch and bound

% to use 20 best

% nUseBranchAndBound = 1000000;

% if (nRegressionType~=2)

% nUseBranchAndBound = nUseBranchAndBound/10;

% end;

% if(nAllCombs>nUseBranchAndBound)

% fprintf('Using branch-and-bound for speedup\n');

% [~,Combs,~,~] = bbdireg(Y,X,nVars,10*nSets);

% nAllCombs = length(Combs);

% end;

% Find full model regression, using L2 as a baseline

r0_trees = Y_trees-X_trees*(X_trees\Y_trees);

stdr0_trees = std(r0_trees);

r0_med = Y_med-X_med*(X_med\Y_med);

stdr0_med = std(r0_med);

% Set combinatorial parameters

CombsF = [zeros(nAllCombs,4) Combs];

nSteps = max(1, round(nAllCombs/20));

% Try all possible combinations, find the best

for nComb = 1:nAllCombs

% Output current progress

if mod(nComb, nSteps)==0 && nSteps>50

fprintf(' %d', round(100*nComb/nAllCombs));

end;

% Build predictor for the current subset

Xtmp_trees = X_trees(:, Combs(nComb,:));

Xtmp_med = X_med(:, Combs(nComb,:));

% Compute least squares regression L2 asa baseline

[pL2_trees, ~, r_trees, ~, stats_trees] = regress(Y_trees, Xtmp_trees);

[pL2_med, ~, r_med, ~, stats_med] = regress(Y_med, Xtmp_med);

% Compute the required regression type

if nRegressionType==0 % median regression

p=fminsearch(@(p)NormLM(Y-Xtmp*p),pL2,opts);

r = Y-Xtmp*p;

Err0 = median(abs(r0));

Err = median(abs(r));

R2 = 1-Err/median(abs(Y));

elseif nRegressionType==1 % L1 regression

p_trees=fminsearch(@(p)NormL1(Y_trees-Xtmp_trees*p),pL2_trees,opts);

p_med=fminsearch(@(p)NormL1(Y_med-Xtmp_med*p),pL2_med,opts);

r_trees = Y_trees-Xtmp_trees*p_trees;

Err0_trees = mean(abs(r0_trees));

Err_trees = mean(abs(r_trees));

R2_trees = 1-(Err_trees/mean(abs(Y_trees)));

r_med = Y_med-Xtmp_med*p_med;

Err0_med = mean(abs(r0_med));

Err_med = mean(abs(r_med));

R2_med = 1-(Err_med/mean(abs(Y_med)));

%weights

Err_complex = 1*Err_trees + 0*Err_med;

R2_avg = (R2_trees+R2_med)/2;

elseif nRegressionType==5 % L1 regression

Err0_trees = mean(abs(r0_trees));

Err_trees = mean(abs(r_trees));

R2_trees = 1-(Err_trees/mean(abs(Y_trees)));

Err0_med = mean(abs(r0_med));

Err_med = mean(abs(r_med));

R2_med = 1-(Err_med/mean(abs(Y_med)));

%weights

Err_complex = 1*Err_trees + 0*Err_med;

R2_avg = (R2_trees+R2_med)/2;

end;

% Record regression statistics

CombsF(nComb,1) = Err_complex;

CombsF(nComb,2) = R2_avg;

CombsF(nComb,3) = R2_med; % R2, similar to stats(1)

CombsF(nComb,4) = R2_trees+1; % R2, similar to stats(1)

% CombsF(nComb,5) = Err_trees; % Error for this regression type %log10(stats(3)); % p-val, log10 of

% CombsF(nComb,6) = Err_med+1; % Error for this regression type %log10(stats(3)); % p-val, log10 of

end;

fprintf('\n');

% Record best predictor sets

CombsF = sortrows(CombsF,1);

fitBest = CombsF(1,1);

cBest = CombsF(1, 4:nVars+3);

len = min(size(CombsF,1), 50);

%len = size(CombsF,1);

CombsF = CombsF(1:len, :);

end

% Helper function to check ALL combinations and return the

% most promising NCombs. Use it when the number of combinations

% is too large to fit into memory

function Combs = FindBestCombs(Y,X,nVars,NCombs)

% Initialize

Combs = [];

% Find the sizes of independent variables X

[~, nX] = size(X);

if(nX<=nVars)

return;

end;

nCombsLog = sum(log10(1:nX)) - sum(log10(1:nVars)) - sum(log10(1:(nX-nVars)));

% If we have realtively small combination count, enumerate all of them

if(nCombsLog<1)

Combs = combnk(1:nX, nVars);

return;

end;

% We have lots of combinaitons. Evaluate each one

NCombsMax = 10*NCombs;

Combs = 1.0e40*ones(NCombsMax, nVars+1);

nc = 0;

H = nextchoose(nX, nVars);

nAll = prod((nX-nVars+1):nX)/prod(1:nVars);

nSteps = max(1, round(nAll/20));

bSorted = 0;

for n=1:nAll-1

% Output current progress

if mod(n, nSteps)==0 && nSteps>50

fprintf(' %d', round(100*n/nAll));

end;

% Find next combination

com = H();

xt = X(:, com);

e = std(Y-xt*(xt\Y));

if(nc>=NCombsMax) % cleanup

Combs = sortrows(Combs, 1);

nc = NCombs+1;

bSorted = 1;

elseif (bSorted>0 && nc>NCombs && e>Combs(NCombs,1))

continue; % no improvement

end;

% OK, good candidate, add

nc = nc+1;

Combs(nc, :) = [e com];

end;

% Cleanup

Combs = sortrows(Combs, 1);

Combs = Combs(1:NCombs, 2:nVars+1);

end

Размещено на Allbest.ru