1022 lines
53 KiB
Matlab
1022 lines
53 KiB
Matlab
function f = SRM(IMAGE)
|
|
% -------------------------------------------------------------------------
|
|
% Copyright (c) 2011 DDE Lab, Binghamton University, NY.
|
|
% All Rights Reserved.
|
|
% -------------------------------------------------------------------------
|
|
% Permission to use, copy, modify, and distribute this software for
|
|
% educational, research and non-profit purposes, without fee, and without a
|
|
% written agreement is hereby granted, provided that this copyright notice
|
|
% appears in all copies. The program is supplied "as is," without any
|
|
% accompanying services from DDE Lab. DDE Lab does not warrant the
|
|
% operation of the program will be uninterrupted or error-free. The
|
|
% end-user understands that the program was developed for research purposes
|
|
% and is advised not to rely exclusively on the program for any reason. In
|
|
% no event shall Binghamton University or DDE Lab be liable to any party
|
|
% for direct, indirect, special, incidental, or consequential damages,
|
|
% including lost profits, arising out of the use of this software. DDE Lab
|
|
% disclaims any warranties, and has no obligations to provide maintenance,
|
|
% support, updates, enhancements or modifications.
|
|
% -------------------------------------------------------------------------
|
|
% Contact: jan@kodovsky.com | fridrich@binghamton.edu | October 2011
|
|
% http://dde.binghamton.edu/download/feature_extractors
|
|
% -------------------------------------------------------------------------
|
|
% Extracts all 106 submodels presented in [1] as part of a rich model for
|
|
% steganalysis of digital images. All features are calculated in the
|
|
% spatial domain and are stored in a structured variable 'f'. For more
|
|
% deatils about the individual submodels, please see the publication [1].
|
|
% Total dimensionality of all 106 submodels is 34671.
|
|
% -------------------------------------------------------------------------
|
|
% Input: IMAGE ... path to the image (can be JPEG)
|
|
% Output: f ....... extracted SRM features in a structured format
|
|
% -------------------------------------------------------------------------
|
|
% [1] Rich Models for Steganalysis of Digital Images, J. Fridrich and J.
|
|
% Kodovsky, IEEE Transactions on Information Forensics and Security, 2011.
|
|
% Under review.
|
|
% -------------------------------------------------------------------------
|
|
|
|
X = double(imread(IMAGE));
|
|
|
|
f = post_processing(all1st(X,1),'f1',1); % 1st order, q=1
|
|
f = post_processing(all1st(X,2),'f1',2,f); % 1st order, q=2
|
|
for q=[1 1.5 2], f = post_processing(all2nd(X,q*2),'f2',q,f); end % 2nd order
|
|
for q=[1 1.5 2], f = post_processing(all3rd(X,q*3),'f3',q,f); end % 3rd order
|
|
for q=[1 1.5 2], f = post_processing(all3x3(X,q*4),'f3x3',q,f); end % 3x3
|
|
for q=[1 1.5 2], f = post_processing(all5x5(X,q*12),'f5x5',q,f); end % 5x5
|
|
|
|
function RESULT = post_processing(DATA,f,q,RESULT)
|
|
|
|
Ss = fieldnames(DATA);
|
|
for Sid = 1:length(Ss)
|
|
VARNAME = [f '_' Ss{Sid} '_q' strrep(num2str(q),'.','')];
|
|
eval(['RESULT.' VARNAME ' = reshape(single(DATA.' Ss{Sid} '),1,[]);' ])
|
|
end
|
|
|
|
% symmetrize
|
|
L = fieldnames(RESULT);
|
|
for i=1:length(L)
|
|
name = L{i}; % feature name
|
|
if name(1)=='s', continue; end
|
|
[T,N,Q] = parse_feaname(name);
|
|
if strcmp(T,''), continue; end
|
|
% symmetrization
|
|
if strcmp(N(1:3),'min') || strcmp(N(1:3),'max')
|
|
% minmax symmetrization
|
|
OUT = ['s' T(2:end) '_minmax' N(4:end) '_' Q];
|
|
if isfield(RESULT,OUT), continue; end
|
|
eval(['Fmin = RESULT.' strrep(name,'max','min') ';']);
|
|
eval(['Fmax = RESULT.' strrep(name,'min','max') ';']);
|
|
F = symfea([Fmin Fmax]',2,4,'mnmx')'; %#ok<*NASGU>
|
|
eval(['RESULT.' OUT ' = single(F);' ]);
|
|
elseif strcmp(N(1:4),'spam')
|
|
% spam symmetrization
|
|
OUT = ['s' T(2:end) '_' N '_' Q];
|
|
if isfield(RESULT,OUT), continue; end
|
|
eval(['Fold = RESULT.' name ';']);
|
|
F = symm1(Fold',2,4)';
|
|
eval(['RESULT.' OUT ' = single(F);' ]);
|
|
end
|
|
end
|
|
% delete RESULT.f*
|
|
L = fieldnames(RESULT);
|
|
for i=1:length(L)
|
|
name = L{i}; % feature name
|
|
if name(1)=='f'
|
|
RESULT = rmfield(RESULT,name);
|
|
end
|
|
end
|
|
% merge spam features
|
|
L = fieldnames(RESULT);
|
|
for i=1:length(L)
|
|
name = L{i}; % feature name
|
|
[T,N,Q] = parse_feaname(name);
|
|
if ~strcmp(N(1:4),'spam'), continue; end
|
|
if strcmp(T,''), continue; end
|
|
if strcmp(N(end),'v')||(strcmp(N,'spam11')&&strcmp(T,'s5x5'))
|
|
elseif strcmp(N(end),'h')
|
|
% h+v union
|
|
OUT = [T '_' N 'v_' Q ];
|
|
if isfield(RESULT,OUT), continue; end
|
|
name2 = strrep(name,'h_','v_');
|
|
eval(['Fh = RESULT.' name ';']);
|
|
eval(['Fv = RESULT.' name2 ';']);
|
|
eval(['RESULT.' OUT ' = [Fh Fv];']);
|
|
RESULT = rmfield(RESULT,name);
|
|
RESULT = rmfield(RESULT,name2);
|
|
elseif strcmp(N,'spam11')
|
|
% KBKV creation
|
|
OUT = ['s35_' N '_' Q];
|
|
if isfield(RESULT,OUT), continue; end
|
|
name1 = strrep(name,'5x5','3x3');
|
|
name2 = strrep(name,'3x3','5x5');
|
|
if ~isfield(RESULT,name1), continue; end
|
|
if ~isfield(RESULT,name2), continue; end
|
|
eval(['F_KB = RESULT.' name1 ';']);
|
|
eval(['F_KV = RESULT.' name2 ';']);
|
|
eval(['RESULT.' OUT ' = [F_KB F_KV];']);
|
|
RESULT = rmfield(RESULT,name1);
|
|
RESULT = rmfield(RESULT,name2);
|
|
end
|
|
end
|
|
function [T,N,Q] = parse_feaname(name)
|
|
[T,N,Q] = deal('');
|
|
P = strfind(name,'_'); if length(P)~=2, return; end
|
|
T = name(1:P(1)-1);
|
|
N = name(P(1)+1:P(2)-1);
|
|
Q = name(P(2)+1:end);
|
|
function g = all1st(X,q)
|
|
%
|
|
% X must be a matrix of doubles or singles (the image) and q is the
|
|
% quantization step (any positive number).
|
|
%
|
|
% Recommended values of q are c, 1.5c, 2c, where c is the central
|
|
% coefficient in the differential (at X(I,J)).
|
|
%
|
|
% This function outputs co-occurrences of ALL 1st-order residuals
|
|
% listed in Figure 1 in our journal HUGO paper (version from June 14),
|
|
% including the naming convention.
|
|
%
|
|
% List of outputted features:
|
|
%
|
|
% 1a) spam14h
|
|
% 1b) spam14v (orthogonal-spam)
|
|
% 1c) minmax22v
|
|
% 1d) minmax24
|
|
% 1e) minmax34v
|
|
% 1f) minmax41
|
|
% 1g) minmax34
|
|
% 1h) minmax48h
|
|
% 1i) minmax54
|
|
%
|
|
% Naming convention:
|
|
%
|
|
% name = {type}{f}{sigma}{scan}
|
|
% type \in {spam, minmax}
|
|
% f \in {1,2,3,4,5} number of filters that are "minmaxed"
|
|
% sigma \in {1,2,3,4,8} symmetry index
|
|
% scan \in {h,v,\emptyset} scan of the cooc matrix (empty = sum of both
|
|
% h and v scans).
|
|
%
|
|
% All odd residuals are implemented the same way simply by
|
|
% narrowing the range for I and J and replacing the residuals --
|
|
% -- they should "stick out" (trcet) in the same direction as
|
|
% the 1st order ones. For example, for the 3rd order:
|
|
%
|
|
% RU = -X(I-2,J+2)+3*X(I-1,J+1)-3*X(I,J)+X(I+1,J-1); ... etc.
|
|
%
|
|
% Note1: The term X(I,J) should always have the "-" sign.
|
|
% Note2: This script does not include s, so, cout, cin versions (weak).
|
|
% This function calls Cooc.m and Quant.m
|
|
|
|
[M N] = size(X); [I,J,T,order] = deal(2:M-1,2:N-1,2,4);
|
|
% Variable names are self-explanatory (R = right, U = up, L = left, D = down)
|
|
[R,L,U,D] = deal(X(I,J+1)-X(I,J),X(I,J-1)-X(I,J),X(I-1,J)-X(I,J),X(I+1,J)-X(I,J));
|
|
[Rq,Lq,Uq,Dq] = deal(Quant(R,q,T),Quant(L,q,T),Quant(U,q,T),Quant(D,q,T));
|
|
[RU,LU,RD,LD] = deal(X(I-1,J+1)-X(I,J),X(I-1,J-1)-X(I,J),X(I+1,J+1)-X(I,J),X(I+1,J-1)-X(I,J));
|
|
[RUq,RDq,LUq,LDq] = deal(Quant(RU,q,T),Quant(RD,q,T),Quant(LU,q,T),Quant(LD,q,T));
|
|
% minmax22h -- to be symmetrized as mnmx, directional, hv-nonsymmetrical.
|
|
[RLq_min,UDq_min,RLq_max,UDq_max] = deal(min(Rq,Lq),min(Uq,Dq),max(Rq,Lq),max(Uq,Dq));
|
|
g.min22h = reshape(Cooc(RLq_min,order,'hor',T) + Cooc(UDq_min,order,'ver',T),[],1);
|
|
g.max22h = reshape(Cooc(RLq_max,order,'hor',T) + Cooc(UDq_max,order,'ver',T),[],1);
|
|
% minmax34h -- to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
[Uq_min,Rq_min,Dq_min,Lq_min] = deal(min(min(Lq,Uq),Rq),min(min(Uq,Rq),Dq),min(min(Rq,Dq),Lq),min(min(Dq,Lq),Uq));
|
|
[Uq_max,Rq_max,Dq_max,Lq_max] = deal(max(max(Lq,Uq),Rq),max(max(Uq,Rq),Dq),max(max(Rq,Dq),Lq),max(max(Dq,Lq),Uq));
|
|
g.min34h = reshape(Cooc([Uq_min;Dq_min],order,'hor',T) + Cooc([Lq_min Rq_min],order,'ver',T),[],1);
|
|
g.max34h = reshape(Cooc([Uq_max;Dq_max],order,'hor',T) + Cooc([Rq_max Lq_max],order,'ver',T),[],1);
|
|
% spam14h/v -- to be symmetrized as spam, directional, hv-nonsymmetrical
|
|
g.spam14h = reshape(Cooc(Rq,order,'hor',T) + Cooc(Uq,order,'ver',T),[],1);
|
|
g.spam14v = reshape(Cooc(Rq,order,'ver',T) + Cooc(Uq,order,'hor',T),[],1);
|
|
% minmax22v -- to be symmetrized as mnmx, directional, hv-nonsymmetrical. Good with higher-order residuals! Note: 22h is bad (too much neighborhood overlap).
|
|
g.min22v = reshape(Cooc(RLq_min,order,'ver',T) + Cooc(UDq_min,order,'hor',T),[],1);
|
|
g.max22v = reshape(Cooc(RLq_max,order,'ver',T) + Cooc(UDq_max,order,'hor',T),[],1);
|
|
% minmax24 -- to be symmetrized as mnmx, directional, hv-symmetrical. Darn good, too.
|
|
[RUq_min,RDq_min,LUq_min,LDq_min] = deal(min(Rq,Uq),min(Rq,Dq),min(Lq,Uq),min(Lq,Dq));
|
|
[RUq_max,RDq_max,LUq_max,LDq_max] = deal(max(Rq,Uq),max(Rq,Dq),max(Lq,Uq),max(Lq,Dq));
|
|
g.min24 = reshape(Cooc([RUq_min;RDq_min;LUq_min;LDq_min],order,'hor',T) + Cooc([RUq_min RDq_min LUq_min LDq_min],order,'ver',T),[],1);
|
|
g.max24 = reshape(Cooc([RUq_max;RDq_max;LUq_max;LDq_max],order,'hor',T) + Cooc([RUq_max RDq_max LUq_max LDq_max],order,'ver',T),[],1);
|
|
% minmax34v -- v works well, h does not, to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
g.min34v = reshape(Cooc([Uq_min Dq_min],order,'ver',T) + Cooc([Rq_min;Lq_min],order,'hor',T),[],1);
|
|
g.max34v = reshape(Cooc([Uq_max Dq_max],order,'ver',T) + Cooc([Rq_max;Lq_max],order,'hor',T),[],1);
|
|
% minmax41 -- to be symmetrized as mnmx, non-directional, hv-symmetrical
|
|
[R_min,R_max] = deal(min(RLq_min,UDq_min),max(RLq_max,UDq_max));
|
|
g.min41 = reshape(Cooc(R_min,order,'hor',T) + Cooc(R_min,order,'ver',T),[],1);
|
|
g.max41 = reshape(Cooc(R_max,order,'hor',T) + Cooc(R_max,order,'ver',T),[],1);
|
|
% minmax34 -- good, to be symmetrized as mnmx, directional, hv-symmetrical
|
|
[RUq_min,RDq_min,LUq_min,LDq_min] = deal(min(RUq_min,RUq),min(RDq_min,RDq),min(LUq_min,LUq),min(LDq_min,LDq));
|
|
[RUq_max,RDq_max,LUq_max,LDq_max] = deal(max(RUq_max,RUq),max(RDq_max,RDq),max(LUq_max,LUq),max(LDq_max,LDq));
|
|
g.min34 = reshape(Cooc([RUq_min;RDq_min;LUq_min;LDq_min],order,'hor',T) + Cooc([RUq_min RDq_min LUq_min LDq_min],order,'ver',T),[],1);
|
|
g.max34 = reshape(Cooc([RUq_max;RDq_max;LUq_max;LDq_max],order,'hor',T) + Cooc([RUq_max RDq_max LUq_max LDq_max],order,'ver',T),[],1);
|
|
% minmax48h -- h better than v, to be symmetrized as mnmx, directional, hv-nonsymmetrical. 48v is almost as good as 48h; for 3rd-order but weaker for 1st-order. Here, I am outputting both but Figure 1 in our paper lists only 48h.
|
|
[RUq_min2,RDq_min2,LDq_min2,LUq_min2] = deal(min(RUq_min,LUq),min(RDq_min,RUq),min(LDq_min,RDq),min(LUq_min,LDq));
|
|
[RUq_min3,RDq_min3,LDq_min3,LUq_min3] = deal(min(RUq_min,RDq),min(RDq_min,LDq),min(LDq_min,LUq),min(LUq_min,RUq));
|
|
g.min48h = reshape(Cooc([RUq_min2;LDq_min2;RDq_min3;LUq_min3],order,'hor',T) + Cooc([RDq_min2 LUq_min2 RUq_min3 LDq_min3],order,'ver',T),[],1);
|
|
g.min48v = reshape(Cooc([RDq_min2;LUq_min2;RUq_min3;LDq_min3],order,'hor',T) + Cooc([RUq_min2 LDq_min2 RDq_min3 LUq_min3],order,'ver',T),[],1);
|
|
[RUq_max2,RDq_max2,LDq_max2,LUq_max2] = deal(max(RUq_max,LUq),max(RDq_max,RUq),max(LDq_max,RDq),max(LUq_max,LDq));
|
|
[RUq_max3,RDq_max3,LDq_max3,LUq_max3] = deal(max(RUq_max,RDq),max(RDq_max,LDq),max(LDq_max,LUq),max(LUq_max,RUq));
|
|
g.max48h = reshape(Cooc([RUq_max2;LDq_max2;RDq_max3;LUq_max3],order,'hor',T) + Cooc([RDq_max2 LUq_max2 RUq_max3 LDq_max3],order,'ver',T),[],1);
|
|
g.max48v = reshape(Cooc([RDq_max2;LUq_max2;RUq_max3;LDq_max3],order,'hor',T) + Cooc([RUq_max2 LDq_max2 RDq_max3 LUq_max3],order,'ver',T),[],1);
|
|
% minmax54 -- to be symmetrized as mnmx, directional, hv-symmetrical
|
|
[RUq_min4,RDq_min4,LDq_min4,LUq_min4] = deal(min(RUq_min2,RDq),min(RDq_min2,LDq),min(LDq_min2,LUq),min(LUq_min2,RUq));
|
|
[RUq_min5,RDq_min5,LDq_min5,LUq_min5] = deal(min(RUq_min3,LUq),min(RDq_min3,RUq),min(LDq_min3,RDq),min(LUq_min3,LDq));
|
|
g.min54 = reshape(Cooc([RUq_min4;LDq_min4;RDq_min5;LUq_min5],order,'hor',T) + Cooc([RDq_min4 LUq_min4 RUq_min5 LDq_min5],order,'ver',T),[],1);
|
|
[RUq_max4,RDq_max4,LDq_max4,LUq_max4] = deal(max(RUq_max2,RDq),max(RDq_max2,LDq),max(LDq_max2,LUq),max(LUq_max2,RUq));
|
|
[RUq_max5,RDq_max5,LDq_max5,LUq_max5] = deal(max(RUq_max3,LUq),max(RDq_max3,RUq),max(LDq_max3,RDq),max(LUq_max3,LDq));
|
|
g.max54 = reshape(Cooc([RUq_max4;LDq_max4;RDq_max5;LUq_max5],order,'hor',T) + Cooc([RDq_max4 LUq_max4 RUq_max5 LDq_max5],order,'ver',T),[],1);
|
|
function g = all2nd(X,q)
|
|
%
|
|
% X must be a matrix of doubles or singles (the image) and q is the
|
|
% quantization step (any positive number).
|
|
%
|
|
% Recommended values of q are c, 1.5c, 2c, where c is the central
|
|
% coefficient in the differential (at X(I,J)).
|
|
%
|
|
% This function outputs co-occurrences of ALL 2nd-order residuals
|
|
% listed in Figure 1 in our journal HUGO paper (version from June 14),
|
|
% including the naming convention.
|
|
%
|
|
% List of outputted features:
|
|
%
|
|
% 1a) spam12h
|
|
% 1b) spam12v (orthogonal-spam)
|
|
% 1c) minmax21
|
|
% 1d) minmax41
|
|
% 1e) minmax24h (24v is also outputted but not listed in Figure 1)
|
|
% 1f) minmax32
|
|
%
|
|
% Naming convention:
|
|
%
|
|
% name = {type}{f}{sigma}{scan}
|
|
% type \in {spam, minmax}
|
|
% f \in {1,2,3,4,5} number of filters that are "minmaxed"
|
|
% sigma \in {1,2,3,4,8} symmetry index
|
|
% scan \in {h,v,\emptyset} scan of the cooc matrix (empty = sum of both
|
|
% h and v scans).
|
|
%
|
|
% All even residuals are implemented the same way simply by
|
|
% narrowing the range for I and J and replacing the residuals.
|
|
%
|
|
% Note1: The term X(I,J) should always have the "-" sign.
|
|
% Note2: This script does not include s, so, cout, cin versions (weak).
|
|
%
|
|
% This function calls Residual.m, Cooc.m, and Quant.m
|
|
|
|
[T,order] = deal(2,4);
|
|
% 2nd-order residuals are implemented using Residual.m
|
|
[Dh,Dv,Dd,Dm] = deal(Residual(X,2,'hor'),Residual(X,2,'ver'),Residual(X,2,'diag'),Residual(X,2,'mdiag'));
|
|
[Yh,Yv,Yd,Ym] = deal(Quant(Dh,q,T),Quant(Dv,q,T),Quant(Dd,q,T),Quant(Dm,q,T));
|
|
% spam12h/v
|
|
g.spam12h = reshape(Cooc(Yh,order,'hor',T) + Cooc(Yv,order,'ver',T),[],1);
|
|
g.spam12v = reshape(Cooc(Yh,order,'ver',T) + Cooc(Yv,order,'hor',T),[],1);
|
|
% minmax21
|
|
[Dmin,Dmax] = deal(min(Yh,Yv),max(Yh,Yv));
|
|
g.min21 = reshape(Cooc(Dmin,order,'hor',T) + Cooc(Dmin,order,'ver',T),[],1);
|
|
g.max21 = reshape(Cooc(Dmax,order,'hor',T) + Cooc(Dmax,order,'ver',T),[],1);
|
|
% minmax41
|
|
[Dmin2,Dmax2] = deal(min(Dmin,min(Yd,Ym)),max(Dmax,max(Yd,Ym)));
|
|
g.min41 = reshape(Cooc(Dmin2,order,'hor',T) + Cooc(Dmin2,order,'ver',T),[],1);
|
|
g.max41 = reshape(Cooc(Dmax2,order,'hor',T) + Cooc(Dmax2,order,'ver',T),[],1);
|
|
% minmax32 -- good, directional, hv-symmetrical, to be symmetrized as mnmx
|
|
[RUq_min,RDq_min] = deal(min(Dmin,Ym),min(Dmin,Yd));
|
|
[RUq_max,RDq_max] = deal(max(Dmax,Ym),max(Dmax,Yd));
|
|
g.min32 = reshape(Cooc([RUq_min;RDq_min],order,'hor',T) + Cooc([RUq_min RDq_min],order,'ver',T),[],1);
|
|
g.max32 = reshape(Cooc([RUq_max;RDq_max],order,'hor',T) + Cooc([RUq_max RDq_max],order,'ver',T),[],1);
|
|
% minmax24h,v -- both "not bad," h slightly better, directional, hv-nonsymmetrical, to be symmetrized as mnmx
|
|
[RUq_min2,RDq_min2,RUq_min3,LUq_min3] = deal(min(Ym,Yh),min(Yd,Yh),min(Ym,Yv),min(Yd,Yv));
|
|
g.min24h = reshape(Cooc([RUq_min2;RDq_min2],order,'hor',T)+Cooc([RUq_min3 LUq_min3],order,'ver',T),[],1);
|
|
g.min24v = reshape(Cooc([RUq_min2 RDq_min2],order,'ver',T)+Cooc([RUq_min3;LUq_min3],order,'hor',T),[],1);
|
|
[RUq_max2,RDq_max2,RUq_max3,LUq_max3] = deal(max(Ym,Yh),max(Yd,Yh),max(Ym,Yv),max(Yd,Yv));
|
|
g.max24h = reshape(Cooc([RUq_max2;RDq_max2],order,'hor',T)+Cooc([RUq_max3 LUq_max3],order,'ver',T),[],1);
|
|
g.max24v = reshape(Cooc([RUq_max2 RDq_max2],order,'ver',T)+Cooc([RUq_max3;LUq_max3],order,'hor',T),[],1);
|
|
function g = all3rd(X,q)
|
|
%
|
|
% X must be a matrix of doubles or singles (the image) and q is the
|
|
% quantization step (any positive number).
|
|
%
|
|
% Recommended values of q are c, 1.5c, 2c, where c is the central
|
|
% coefficient in the differential (at X(I,J)).
|
|
%
|
|
% This function outputs co-occurrences of ALL 3rd-order residuals
|
|
% listed in Figure 1 in our journal HUGO paper (version from June 14),
|
|
% including the naming convention.
|
|
%
|
|
% List of outputted features:
|
|
%
|
|
% 1a) spam14h
|
|
% 1b) spam14v (orthogonal-spam)
|
|
% 1c) minmax22v
|
|
% 1d) minmax24
|
|
% 1e) minmax34v
|
|
% 1f) minmax41
|
|
% 1g) minmax34
|
|
% 1h) minmax48h
|
|
% 1i) minmax54
|
|
%
|
|
% Naming convention:
|
|
%
|
|
% name = {type}{f}{sigma}{scan}
|
|
% type \in {spam, minmax}
|
|
% f \in {1,2,3,4,5} number of filters that are "minmaxed"
|
|
% sigma \in {1,2,3,4,8} symmetry index
|
|
% scan \in {h,v,\emptyset} scan of the cooc matrix (empty = sum of both
|
|
% h and v scans).
|
|
%
|
|
% All odd residuals are implemented the same way simply by
|
|
% narrowing the range for I and J and replacing the residuals --
|
|
% -- they should "stick out" (trcet) in the same direction as
|
|
% the 1st order ones. For example, for the 3rd order:
|
|
%
|
|
% RU = -X(I-2,J+2)+3*X(I-1,J+1)-3*X(I,J)+X(I+1,J-1); ... etc.
|
|
%
|
|
% Note1: The term X(I,J) should always have the "-" sign.
|
|
% Note2: This script does not include s, so, cout, cin versions (weak).
|
|
|
|
[M N] = size(X); [I,J,T,order] = deal(3:M-2,3:N-2,2,4);
|
|
[R,L,U,D] = deal(-X(I,J+2)+3*X(I,J+1)-3*X(I,J)+X(I,J-1),-X(I,J-2)+3*X(I,J-1)-3*X(I,J)+X(I,J+1),-X(I-2,J)+3*X(I-1,J)-3*X(I,J)+X(I+1,J),-X(I+2,J)+3*X(I+1,J)-3*X(I,J)+X(I-1,J));
|
|
[Rq,Lq,Uq,Dq] = deal(Quant(R,q,T),Quant(L,q,T),Quant(U,q,T),Quant(D,q,T));
|
|
[RU,LU,RD,LD] = deal(-X(I-2,J+2)+3*X(I-1,J+1)-3*X(I,J)+X(I+1,J-1),-X(I-2,J-2)+3*X(I-1,J-1)-3*X(I,J)+X(I+1,J+1),-X(I+2,J+2)+3*X(I+1,J+1)-3*X(I,J)+X(I-1,J-1),-X(I+2,J-2)+3*X(I+1,J-1)-3*X(I,J)+X(I-1,J+1));
|
|
[RUq,RDq,LUq,LDq] = deal(Quant(RU,q,T),Quant(RD,q,T),Quant(LU,q,T),Quant(LD,q,T));
|
|
% minmax22h -- to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
[RLq_min,UDq_min] = deal(min(Rq,Lq),min(Uq,Dq));
|
|
[RLq_max,UDq_max] = deal(max(Rq,Lq),max(Uq,Dq));
|
|
g.min22h = reshape(Cooc(RLq_min,order,'hor',T) + Cooc(UDq_min,order,'ver',T),[],1);
|
|
g.max22h = reshape(Cooc(RLq_max,order,'hor',T) + Cooc(UDq_max,order,'ver',T),[],1);
|
|
% minmax34h -- to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
[Uq_min,Rq_min,Dq_min,Lq_min] = deal(min(RLq_min,Uq),min(UDq_min,Rq),min(RLq_min,Dq),min(UDq_min,Lq));
|
|
[Uq_max,Rq_max,Dq_max,Lq_max] = deal(max(RLq_max,Uq),max(UDq_max,Rq),max(RLq_max,Dq),max(UDq_max,Lq));
|
|
g.min34h = reshape(Cooc([Uq_min;Dq_min],order,'hor',T)+Cooc([Rq_min Lq_min],order,'ver',T),[],1);
|
|
g.max34h = reshape(Cooc([Uq_max;Dq_max],order,'hor',T)+Cooc([Rq_max Lq_max],order,'ver',T),[],1);
|
|
% spam14h,v -- to be symmetrized as spam, directional, hv-nonsymmetrical
|
|
g.spam14h = reshape(Cooc(Rq,order,'hor',T) + Cooc(Uq,order,'ver',T),[],1);
|
|
g.spam14v = reshape(Cooc(Rq,order,'ver',T) + Cooc(Uq,order,'hor',T),[],1);
|
|
% minmax22v -- to be symmetrized as mnmx, directional, hv-nonsymmetrical. Good with higher-order residuals! Note: 22h is bad (too much neighborhood overlap).
|
|
g.min22v = reshape(Cooc(RLq_min,order,'ver',T) + Cooc(UDq_min,order,'hor',T),[],1);
|
|
g.max22v = reshape(Cooc(RLq_max,order,'ver',T) + Cooc(UDq_max,order,'hor',T),[],1);
|
|
% minmax24 -- to be symmetrized as mnmx, directional, hv-symmetrical Note: Darn good, too.
|
|
[RUq_min,RDq_min,LUq_min,LDq_min] = deal(min(Rq,Uq),min(Rq,Dq),min(Lq,Uq),min(Lq,Dq));
|
|
[RUq_max,RDq_max,LUq_max,LDq_max] = deal(max(Rq,Uq),max(Rq,Dq),max(Lq,Uq),max(Lq,Dq));
|
|
g.min24 = reshape(Cooc([RUq_min;RDq_min;LUq_min;LDq_min],order,'hor',T) + Cooc([RUq_min RDq_min LUq_min LDq_min],order,'ver',T),[],1);
|
|
g.max24 = reshape(Cooc([RUq_max;RDq_max;LUq_max;LDq_max],order,'hor',T) + Cooc([RUq_max RDq_max LUq_max LDq_max],order,'ver',T),[],1);
|
|
% minmax34v -- v works well, h does not, to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
g.min34v = reshape(Cooc([Uq_min Dq_min],order,'ver',T) + Cooc([Rq_min;Lq_min],order,'hor',T),[],1);
|
|
g.max34v = reshape(Cooc([Uq_max Dq_max],order,'ver',T) + Cooc([Rq_max;Lq_max],order,'hor',T),[],1);
|
|
% minmax41 -- unknown performance as of 6/14/11, to be symmetrized as mnmx, non-directional, hv-symmetrical
|
|
[R_min,R_max] = deal(min(RUq_min,LDq_min),max(RUq_max,LDq_max));
|
|
g.min41 = reshape(Cooc(R_min,order,'hor',T) + Cooc(R_min,order,'ver',T),[],1);
|
|
g.max41 = reshape(Cooc(R_max,order,'hor',T) + Cooc(R_max,order,'ver',T),[],1);
|
|
% minmax34 -- good, to be symmetrized as mnmx, directional, hv-symmetrical
|
|
[RUq_min2,RDq_min2,LUq_min2,LDq_min2] = deal(min(RUq_min,RUq),min(RDq_min,RDq),min(LUq_min,LUq),min(LDq_min,LDq));
|
|
[RUq_max2,RDq_max2,LUq_max2,LDq_max2] = deal(max(RUq_max,RUq),max(RDq_max,RDq),max(LUq_max,LUq),max(LDq_max,LDq));
|
|
g.min34 = reshape(Cooc([RUq_min2;RDq_min2;LUq_min2;LDq_min2],order,'hor',T) + Cooc([RUq_min2 RDq_min2 LUq_min2 LDq_min2],order,'ver',T),[],1);
|
|
g.max34 = reshape(Cooc([RUq_max2;RDq_max2;LUq_max2;LDq_max2],order,'hor',T) + Cooc([RUq_max2 RDq_max2 LUq_max2 LDq_max2],order,'ver',T),[],1);
|
|
% minmax48h -- h better than v, to be symmetrized as mnmx, directional, hv-nonsymmetrical. 48v is almost as good as 48h for 3rd-order but weaker for 1st-order. Here, I am outputting both but Figure 1 in our paper lists only 48h.
|
|
[RUq_min3,RDq_min3,LDq_min3,LUq_min3] = deal(min(RUq_min2,LUq),min(RDq_min2,RUq),min(LDq_min2,RDq),min(LUq_min2,LDq));
|
|
[RUq_min4,RDq_min4,LDq_min4,LUq_min4] = deal(min(RUq_min2,RDq),min(RDq_min2,LDq),min(LDq_min2,LUq),min(LUq_min2,RUq));
|
|
g.min48h = reshape(Cooc([RUq_min3;LDq_min3;RDq_min4;LUq_min4],order,'hor',T)+Cooc([RDq_min3 LUq_min3 RUq_min4 LDq_min4],order,'ver',T),[],1);
|
|
g.min48v = reshape(Cooc([RUq_min3 LDq_min3 RDq_min4 LUq_min4],order,'ver',T)+Cooc([RDq_min3;LUq_min3;RUq_min4;LDq_min4],order,'hor',T),[],1);
|
|
[RUq_max3,RDq_max3,LDq_max3,LUq_max3] = deal(max(RUq_max2,LUq),max(RDq_max2,RUq),max(LDq_max2,RDq),max(LUq_max2,LDq));
|
|
[RUq_max4,RDq_max4,LDq_max4,LUq_max4] = deal(max(RUq_max2,RDq),max(RDq_max2,LDq),max(LDq_max2,LUq),max(LUq_max2,RUq));
|
|
g.max48h = reshape(Cooc([RUq_max3;LDq_max3;RDq_max4;LUq_max4],order,'hor',T)+Cooc([RDq_max3 LUq_max3 RUq_max4 LDq_max4],order,'ver',T),[],1);
|
|
g.max48v = reshape(Cooc([RUq_max3 LDq_max3 RDq_max4 LUq_max4],order,'ver',T)+Cooc([RDq_max3;LUq_max3;RUq_max4;LDq_max4],order,'hor',T),[],1);
|
|
% minmax54 -- to be symmetrized as mnmx, directional, hv-symmetrical
|
|
[RUq_min5,RDq_min5,LDq_min5,LUq_min5] = deal(min(RUq_min3,RDq),min(RDq_min3,LDq),min(LDq_min3,LUq),min(LUq_min3,RUq));
|
|
[RUq_max5,RDq_max5,LDq_max5,LUq_max5] = deal(max(RUq_max3,RDq),max(RDq_max3,LDq),max(LDq_max3,LUq),max(LUq_max3,RUq));
|
|
g.min54 = reshape(Cooc([RUq_min5;LDq_min5;RDq_min5;LUq_min5],order,'hor',T) + Cooc([RDq_min5 LUq_min5 RUq_min5 LDq_min5],order,'ver',T),[],1);
|
|
g.max54 = reshape(Cooc([RUq_max5;LDq_max5;RDq_max5;LUq_max5],order,'hor',T) + Cooc([RDq_max5 LUq_max5 RUq_max5 LDq_max5],order,'ver',T),[],1);
|
|
function g = all3x3(X,q)
|
|
% This function outputs co-occurrences of ALL residuals based on the
|
|
% KB kernel and its "halves" (EDGE residuals) as listed in Figure 1
|
|
% in our journal HUGO paper (version from June 14), including the naming
|
|
% convention.
|
|
[T,order] = deal(2,4);
|
|
% spam11 (old name KB residual), good, non-directional, hv-symmetrical, to be symmetrized as spam
|
|
D = Residual(X,2,'KB'); Y = Quant(D,q,T);
|
|
g.spam11 = reshape(Cooc(Y,order,'hor',T) + Cooc(Y,order,'ver',T),[],1);
|
|
% EDGE residuals
|
|
D = Residual(X,2,'edge-h');Du = D(:,1:size(D,2)/2);Db = D(:,size(D,2)/2+1:end);
|
|
D = Residual(X,2,'edge-v');Dl = D(:,1:size(D,2)/2);Dr = D(:,size(D,2)/2+1:end);
|
|
[Yu,Yb,Yl,Yr] = deal(Quant(Du,q,T),Quant(Db,q,T),Quant(Dl,q,T),Quant(Dr,q,T));
|
|
% spam14h,v not bad, directional, hv-nonsym, to be symmetrized as spam
|
|
g.spam14v = reshape(Cooc([Yu Yb],order,'ver',T) + Cooc([Yl;Yr],order,'hor',T),[],1);
|
|
g.spam14h = reshape(Cooc([Yu;Yb],order,'hor',T) + Cooc([Yl Yr],order,'ver',T),[],1);
|
|
% minmax24 -- EXCELLENT, directional, hv-sym, to be symmetrized as mnmx
|
|
[Dmin1,Dmin2,Dmin3,Dmin4] = deal(min(Yu,Yl),min(Yb,Yr),min(Yu,Yr),min(Yb,Yl));
|
|
g.min24 = reshape(Cooc([Dmin1 Dmin2 Dmin3 Dmin4],order,'ver',T) + Cooc([Dmin1;Dmin2;Dmin3;Dmin4],order,'hor',T),[],1);
|
|
[Dmax1,Dmax2,Dmax3,Dmax4] = deal(max(Yu,Yl),max(Yb,Yr),max(Yu,Yr),max(Yb,Yl));
|
|
g.max24 = reshape(Cooc([Dmax1 Dmax2 Dmax3 Dmax4],order,'ver',T) + Cooc([Dmax1;Dmax2;Dmax3;Dmax4],order,'hor',T),[],1);
|
|
% minmax22 - hv-nonsymmetrical
|
|
% min22h -- good, to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
% min22v -- EXCELLENT - to be symmetrized as mnmx, directional,
|
|
[UEq_min,REq_min] = deal(min(Yu,Yb),min(Yr,Yl));
|
|
g.min22h = reshape(Cooc(UEq_min,order,'hor',T) + Cooc(REq_min,order,'ver',T),[],1);
|
|
g.min22v = reshape(Cooc(UEq_min,order,'ver',T) + Cooc(REq_min,order,'hor',T),[],1);
|
|
[UEq_max,REq_max] = deal(max(Yu,Yb),max(Yr,Yl));
|
|
g.max22h = reshape(Cooc(UEq_max,order,'hor',T) + Cooc(REq_max,order,'ver',T),[],1);
|
|
g.max22v = reshape(Cooc(UEq_max,order,'ver',T) + Cooc(REq_max,order,'hor',T),[],1);
|
|
% minmax41 -- good, non-directional, hv-sym, to be symmetrized as mnmx
|
|
[Dmin5,Dmax5] = deal(min(Dmin1,Dmin2),max(Dmax1,Dmax2));
|
|
g.min41 = reshape(Cooc(Dmin5,order,'ver',T) + Cooc(Dmin5,order,'hor',T),[],1);
|
|
g.max41 = reshape(Cooc(Dmax5,order,'ver',T) + Cooc(Dmax5,order,'hor',T),[],1);
|
|
function g = all5x5(X,q)
|
|
% This function outputs co-occurrences of ALL residuals based on the
|
|
% KV kernel and its "halves" (EDGE residuals) as listed in Figure 1
|
|
% in our journal HUGO paper (version from June 14), including the naming
|
|
% convention.
|
|
[M N] = size(X); [I,J,T,order] = deal(3:M-2,3:N-2,2,4);
|
|
% spam11 (old name KV residual), good, non-directional, hv-symmetrical, to be symmetrized as spam
|
|
D = Residual(X,3,'KV'); Y = Quant(D,q,T);
|
|
g.spam11 = reshape(Cooc(Y,order,'hor',T) + Cooc(Y,order,'ver',T),[],1);
|
|
% EDGE residuals
|
|
Du = 8*X(I,J-1)+8*X(I-1,J)+8*X(I,J+1)-6*X(I-1,J-1)-6*X(I-1,J+1)-2*X(I,J-2)-2*X(I,J+2)-2*X(I-2,J)+2*X(I-1,J-2)+2*X(I-2,J-1)+2*X(I-2,J+1)+2*X(I-1,J+2)-X(I-2,J-2)-X(I-2,J+2)-12*X(I,J);
|
|
Dr = 8*X(I-1,J)+8*X(I,J+1)+8*X(I+1,J)-6*X(I-1,J+1)-6*X(I+1,J+1)-2*X(I-2,J)-2*X(I+2,J)-2*X(I,J+2)+2*X(I-2,J+1)+2*X(I-1,J+2)+2*X(I+1,J+2)+2*X(I+2,J+1)-X(I-2,J+2)-X(I+2,J+2)-12*X(I,J);
|
|
Db = 8*X(I,J+1)+8*X(I+1,J)+8*X(I,J-1)-6*X(I+1,J+1)-6*X(I+1,J-1)-2*X(I,J-2)-2*X(I,J+2)-2*X(I+2,J)+2*X(I+1,J+2)+2*X(I+2,J+1)+2*X(I+2,J-1)+2*X(I+1,J-2)-X(I+2,J+2)-X(I+2,J-2)-12*X(I,J);
|
|
Dl = 8*X(I+1,J)+8*X(I,J-1)+8*X(I-1,J)-6*X(I+1,J-1)-6*X(I-1,J-1)-2*X(I-2,J)-2*X(I+2,J)-2*X(I,J-2)+2*X(I+2,J-1)+2*X(I+1,J-2)+2*X(I-1,J-2)+2*X(I-2,J-1)-X(I+2,J-2)-X(I-2,J-2)-12*X(I,J);
|
|
[Yu,Yb,Yl,Yr] = deal(Quant(Du,q,T),Quant(Db,q,T),Quant(Dl,q,T),Quant(Dr,q,T));
|
|
% spam14v not bad, directional, hv-nonsym, to be symmetrized as spam
|
|
g.spam14v = reshape(Cooc([Yu Yb],order,'ver',T) + Cooc([Yl;Yr],order,'hor',T),[],1);
|
|
g.spam14h = reshape(Cooc([Yu;Yb],order,'hor',T) + Cooc([Yl Yr],order,'ver',T),[],1);
|
|
% minmax24 -- EXCELLENT, directional, hv-sym, to be symmetrized as mnmx
|
|
[Dmin1,Dmin2,Dmin3,Dmin4] = deal(min(Yu,Yl),min(Yb,Yr),min(Yu,Yr),min(Yb,Yl));
|
|
g.min24 = reshape(Cooc([Dmin1 Dmin2 Dmin3 Dmin4],order,'ver',T) + Cooc([Dmin1;Dmin2;Dmin3;Dmin4],order,'hor',T),[],1);
|
|
[Dmax1,Dmax2,Dmax3,Dmax4] = deal(max(Yu,Yl),max(Yb,Yr),max(Yu,Yr),max(Yb,Yl));
|
|
g.max24 = reshape(Cooc([Dmax1 Dmax2 Dmax3 Dmax4],order,'ver',T) + Cooc([Dmax1;Dmax2;Dmax3;Dmax4],order,'hor',T),[],1);
|
|
% minmax22 - hv-nonsymmetrical
|
|
% min22h -- good, to be symmetrized as mnmx, directional, hv-nonsymmetrical
|
|
% min22v -- EXCELLENT - to be symmetrized as mnmx, directional,
|
|
[UEq_min,REq_min] = deal(min(Yu,Yb),min(Yr,Yl));
|
|
g.min22h = reshape(Cooc(UEq_min,order,'hor',T) + Cooc(REq_min,order,'ver',T),[],1);
|
|
g.min22v = reshape(Cooc(UEq_min,order,'ver',T) + Cooc(REq_min,order,'hor',T),[],1);
|
|
[UEq_max,REq_max] = deal(max(Yu,Yb),max(Yr,Yl));
|
|
g.max22h = reshape(Cooc(UEq_max,order,'hor',T) + Cooc(REq_max,order,'ver',T),[],1);
|
|
g.max22v = reshape(Cooc(UEq_max,order,'ver',T) + Cooc(REq_max,order,'hor',T),[],1);
|
|
% minmax41 -- good, non-directional, hv-sym, to be symmetrized as mnmx
|
|
[Dmin5,Dmax5] = deal(min(Dmin1,Dmin2),max(Dmax1,Dmax2));
|
|
g.min41 = reshape(Cooc(Dmin5,order,'ver',T) + Cooc(Dmin5,order,'hor',T),[],1);
|
|
g.max41 = reshape(Cooc(Dmax5,order,'ver',T) + Cooc(Dmax5,order,'hor',T),[],1);
|
|
function f = Cooc(D,order,type,T)
|
|
% Co-occurrence operator to be appied to a 2D array of residuals D \in [-T,T]
|
|
% T ... threshold
|
|
% order ... cooc order \in {1,2,3,4,5}
|
|
% type ... cooc type \in {hor,ver,diag,mdiag,square,square-ori,hvdm}
|
|
% f ... an array of size (2T+1)^order
|
|
|
|
B = 2*T+1;
|
|
if max(abs(D(:))) > T, fprintf('*** ERROR in Cooc.m: Residual out of range ***\n'), end
|
|
|
|
switch order
|
|
case 1
|
|
f = hist(D(:),-T:T);
|
|
case 2
|
|
f = zeros(B,B);
|
|
if strcmp(type,'hor'), L = D(:,1:end-1); R = D(:,2:end);end
|
|
if strcmp(type,'ver'), L = D(1:end-1,:); R = D(2:end,:);end
|
|
if strcmp(type,'diag'), L = D(1:end-1,1:end-1); R = D(2:end,2:end);end
|
|
if strcmp(type,'mdiag'), L = D(1:end-1,2:end); R = D(2:end,1:end-1);end
|
|
for i = -T : T
|
|
R2 = R(L(:)==i);
|
|
for j = -T : T
|
|
f(i+T+1,j+T+1) = sum(R2(:)==j);
|
|
end
|
|
end
|
|
case 3
|
|
f = zeros(B,B,B);
|
|
if strcmp(type,'hor'), L = D(:,1:end-2); C = D(:,2:end-1); R = D(:,3:end);end
|
|
if strcmp(type,'ver'), L = D(1:end-2,:); C = D(2:end-1,:); R = D(3:end,:);end
|
|
if strcmp(type,'diag'), L = D(1:end-2,1:end-2); C = D(2:end-1,2:end-1); R = D(3:end,3:end);end
|
|
if strcmp(type,'mdiag'), L = D(1:end-2,3:end); C = D(2:end-1,2:end-1); R = D(3:end,1:end-2);end
|
|
for i = -T : T
|
|
C2 = C(L(:)==i);
|
|
R2 = R(L(:)==i);
|
|
for j = -T : T
|
|
R3 = R2(C2(:)==j);
|
|
for k = -T : T
|
|
f(i+T+1,j+T+1,k+T+1) = sum(R3(:)==k);
|
|
end
|
|
end
|
|
end
|
|
case 4
|
|
f = zeros(B,B,B,B);
|
|
if strcmp(type,'hor'), L = D(:,1:end-3); C = D(:,2:end-2); E = D(:,3:end-1); R = D(:,4:end);end
|
|
if strcmp(type,'ver'), L = D(1:end-3,:); C = D(2:end-2,:); E = D(3:end-1,:); R = D(4:end,:);end
|
|
if strcmp(type,'diag'), L = D(1:end-3,1:end-3); C = D(2:end-2,2:end-2); E = D(3:end-1,3:end-1); R = D(4:end,4:end);end
|
|
if strcmp(type,'mdiag'), L = D(4:end,1:end-3); C = D(3:end-1,2:end-2); E = D(2:end-2,3:end-1); R = D(1:end-3,4:end);end
|
|
if strcmp(type,'square'), L = D(2:end,1:end-1); C = D(2:end,2:end); E = D(1:end-1,2:end); R = D(1:end-1,1:end-1);end
|
|
if strcmp(type,'square-ori'), [M, N] = size(D); Dh = D(:,1:M); Dv = D(:,M+1:2*M);
|
|
L = Dh(2:end,1:end-1); C = Dv(2:end,2:end); E = Dh(1:end-1,2:end); R = Dv(1:end-1,1:end-1);end
|
|
if strcmp(type,'hvdm'), [M, N] = size(D); L = D(:,1:M); C = D(:,M+1:2*M); E = D(:,2*M+1:3*M); R = D(:,3*M+1:4*M);end
|
|
if strcmp(type,'s-in'), [M, N] = size(D); Dh = D(:,1:M); Dv = D(:,M+1:2*M); Dh1 = D(:,2*M+1:3*M); Dv1 = D(:,3*M+1:4*M);
|
|
L = Dh(2:end,1:end-1); C = Dh1(2:end,2:end); E = Dh1(1:end-1,2:end); R = Dh(1:end-1,1:end-1);end
|
|
if strcmp(type,'s-out'), [M, N] = size(D); Dh = D(:,1:M); Dv = D(:,M+1:2*M); Dh1 = D(:,2*M+1:3*M); Dv1 = D(:,3*M+1:4*M);
|
|
L = Dh1(2:end,1:end-1); C = Dh(2:end,2:end); E = Dh(1:end-1,2:end); R = Dh1(1:end-1,1:end-1);end
|
|
if strcmp(type,'ori-in'), [M, N] = size(D); Dh = D(:,1:M); Dv = D(:,M+1:2*M); Dh1 = D(:,2*M+1:3*M); Dv1 = D(:,3*M+1:4*M);
|
|
L = Dh(2:end,1:end-1); C = Dv1(2:end,2:end); E = Dh1(1:end-1,2:end); R = Dv(1:end-1,1:end-1);end
|
|
if strcmp(type,'ori-out'),[M, N] = size(D); Dh = D(:,1:M); Dv = D(:,M+1:2*M); Dh1 = D(:,2*M+1:3*M); Dv1 = D(:,3*M+1:4*M);
|
|
L = Dh1(2:end,1:end-1); C = Dv(2:end,2:end); E = Dh(1:end-1,2:end); R = Dv1(1:end-1,1:end-1);end
|
|
for i = -T : T
|
|
ind = (L(:)==i); C2 = C(ind); E2 = E(ind); R2 = R(ind);
|
|
for j = -T : T
|
|
ind = (C2(:)==j); E3 = E2(ind); R3 = R2(ind);
|
|
for k = -T : T
|
|
R4 = R3(E3(:)==k);
|
|
for l = -T : T
|
|
f(i+T+1,j+T+1,k+T+1,l+T+1) = sum(R4(:)==l);
|
|
end
|
|
end
|
|
end
|
|
end
|
|
case 5
|
|
f = zeros(B,B,B,B,B);
|
|
if strcmp(type,'hor'),L = D(:,1:end-4); C = D(:,2:end-3); E = D(:,3:end-2); F = D(:,4:end-1); R = D(:,5:end);end
|
|
if strcmp(type,'ver'),L = D(1:end-4,:); C = D(2:end-3,:); E = D(3:end-2,:); F = D(4:end-1,:); R = D(5:end,:);end
|
|
|
|
for i = -T : T
|
|
ind = (L(:)==i); C2 = C(ind); E2 = E(ind); F2 = F(ind); R2 = R(ind);
|
|
for j = -T : T
|
|
ind = (C2(:)==j); E3 = E2(ind); F3 = F2(ind); R3 = R2(ind);
|
|
for k = -T : T
|
|
ind = (E3(:)==k); F4 = F3(ind); R4 = R3(ind);
|
|
for l = -T : T
|
|
R5 = R4(F4(:)==l);
|
|
for m = -T : T
|
|
f(i+T+1,j+T+1,k+T+1,l+T+1,m+T+1) = sum(R5(:)==m);
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
|
|
% normalization
|
|
f = double(f);
|
|
fsum = sum(f(:));
|
|
if fsum>0, f = f/fsum; end
|
|
|
|
function Y = Quant(X,q,T)
|
|
% Quantization routine
|
|
% X ... variable to be quantized/truncated
|
|
% T ... threshold
|
|
% q ... quantization step (with type = 'scalar') or a vector of increasing
|
|
% non-negative integers outlining the quantization process.
|
|
% Y ... quantized/truncated variable
|
|
% Example 0: when q is a positive scalar, Y = trunc(round(X/q),T)
|
|
% Example 1: q = [0 1 2 3] quantizes 0 to 0, 1 to 1, 2 to 2, [3,Inf) to 3,
|
|
% (-Inf,-3] to -3, -2 to -2, -1 to -1. It is equivalent to Quant(.,3,1).
|
|
% Example 2: q = [0 2 4 5] quantizes 0 to 0, {1,2} to 1, {3,4} to 2,
|
|
% [5,Inf) to 3, and similarly with the negatives.
|
|
% Example 3: q = [1 2] quantizes {-1,0,1} to 0, [2,Inf) to 1, (-Inf,-2] to -1.
|
|
% Example 4: q = [1 3 7 15 16] quantizes {-1,0,1} to 0, {2,3} to 1, {4,5,6,7}
|
|
% to 2, {8,9,10,11,12,13,14,15} to 3, [16,Inf) to 4, and similarly the
|
|
% negatives.
|
|
|
|
if numel(q) == 1
|
|
if q > 0, Y = trunc(round(X/q),T);
|
|
else fprintf('*** ERROR: Attempt to quantize with non-positive step. ***\n'),end
|
|
else
|
|
q = round(q); % Making sure the vector q is made of integers
|
|
if min(q(2:end)-q(1:end-1)) <= 0
|
|
fprintf('*** ERROR: quantization vector not strictly increasing. ***\n')
|
|
end
|
|
if min(q) < 0, fprintf('*** ERROR: Attempt to quantize with negative step. ***\n'),end
|
|
|
|
T = q(end); % The last value determines the truncation threshold
|
|
v = zeros(1,2*T+1); % value-substitution vector
|
|
Y = trunc(X,T)+T+1; % Truncated X and shifted to positive values
|
|
if q(1) == 0
|
|
v(T+1) = 0; z = 1; ind = T+2;
|
|
for i = 2 : numel(q)
|
|
v(ind:ind+q(i)-q(i-1)-1) = z;
|
|
ind = ind+q(i)-q(i-1);
|
|
z = z+1;
|
|
end
|
|
v(1:T) = -v(end:-1:T+2);
|
|
else
|
|
v(T+1-q(1):T+1+q(1)) = 0; z = 1; ind = T+2+q(1);
|
|
for i = 2 : numel(q)
|
|
v(ind:ind+q(i)-q(i-1)-1) = z;
|
|
ind = ind+q(i)-q(i-1);
|
|
z = z+1;
|
|
end
|
|
v(1:T-q(1)) = -v(end:-1:T+2+q(1));
|
|
end
|
|
Y = v(Y); % The actual quantization :)
|
|
end
|
|
function Z = trunc(X,T)
|
|
% Truncation to [-T,T]
|
|
Z = X;
|
|
Z(Z > T) = T;
|
|
Z(Z < -T) = -T;
|
|
function fsym = symfea(f,T,order,type)
|
|
% Marginalization by sign and directional symmetry for a feature vector
|
|
% stored as one of our 2*(2T+1)^order-dimensional feature vectors. This
|
|
% routine should be used for features possiessing sign and directional
|
|
% symmetry, such as spam-like features or 3x3 features. It should NOT be
|
|
% used for features from MINMAX residuals. Use the alternative
|
|
% symfea_minmax for this purpose.
|
|
% The feature f is assumed to be a 2dim x database_size matrix of features
|
|
% stored as columns (e.g., hor+ver, diag+minor_diag), with dim =
|
|
% 2(2T+1)^order.
|
|
|
|
[dim,N] = size(f);
|
|
B = 2*T+1;
|
|
c = B^order;
|
|
ERR = 1;
|
|
|
|
if strcmp(type,'spam')
|
|
if dim == 2*c
|
|
switch order % Reduced dimensionality for a B^order dimensional feature vector
|
|
case 1, red = T + 1;
|
|
case 2, red = (T + 1)^2;
|
|
case 3, red = 1 + 3*T + 4*T^2 + 2*T^3;
|
|
case 4, red = B^2 + 4*T^2*(T + 1)^2;
|
|
case 5, red = 1/4*(B^2 + 1)*(B^3 + 1);
|
|
end
|
|
fsym = zeros(2*red,N);
|
|
|
|
for i = 1 : N
|
|
switch order
|
|
case 1, cube = f(1:c,i);
|
|
case 2, cube = reshape(f(1:c,i),[B B]);
|
|
case 3, cube = reshape(f(1:c,i),[B B B]);
|
|
case 4, cube = reshape(f(1:c,i),[B B B B]);
|
|
case 5, cube = reshape(f(1:c,i),[B B B B B]);
|
|
end
|
|
% [size(symm_dir(cube,T,order)) red]
|
|
fsym(1:red,i) = symm(cube,T,order);
|
|
switch order
|
|
case 1, cube = f(c+1:2*c,i);
|
|
case 2, cube = reshape(f(c+1:2*c,i),[B B]);
|
|
case 3, cube = reshape(f(c+1:2*c,i),[B B B]);
|
|
case 4, cube = reshape(f(c+1:2*c,i),[B B B B]);
|
|
case 5, cube = reshape(f(c+1:2*c,i),[B B B B B]);
|
|
end
|
|
fsym(red+1:2*red,i) = symm(cube,T,order);
|
|
end
|
|
else
|
|
fsym = [];
|
|
fprintf('*** ERROR: feature dimension is not 2x(2T+1)^order. ***\n')
|
|
end
|
|
ERR = 0;
|
|
end
|
|
|
|
if strcmp(type,'mnmx')
|
|
if dim == 2*c
|
|
switch order
|
|
case 3, red = B^3 - T*B^2; % Dim of the marginalized set is (2T+1)^3-T*(2T+1)^2
|
|
case 4, red = B^4 - 2*T*(T+1)*B^2; % Dim of the marginalized set is (2T+1)^4-2T*(T+1)*(2T+1)^2
|
|
end
|
|
fsym = zeros(red, N);
|
|
for i = 1 : N
|
|
switch order
|
|
case 1, cube_min = f(1:c,i); cube_max = f(c+1:2*c,i);
|
|
case 2, cube_min = reshape(f(1:c,i),[B B]); cube_max = reshape(f(c+1:2*c,i),[B B]); f_signsym = cube_min + cube_max(end:-1:1,end:-1:1);
|
|
case 3, cube_min = reshape(f(1:c,i),[B B B]); cube_max = reshape(f(c+1:2*c,i),[B B B]); f_signsym = cube_min + cube_max(end:-1:1,end:-1:1,end:-1:1);
|
|
case 4, cube_min = reshape(f(1:c,i),[B B B B]); cube_max = reshape(f(c+1:2*c,i),[B B B B]); f_signsym = cube_min + cube_max(end:-1:1,end:-1:1,end:-1:1,end:-1:1);
|
|
case 5, cube_min = reshape(f(1:c,i),[B B B B B]); cube_max = reshape(f(c+1:2*c,i),[B B B B B]); f_signsym = cube_min + cube_max(end:-1:1,end:-1:1,end:-1:1,end:-1:1,end:-1:1);
|
|
end
|
|
% f_signsym = cube_min + cube_max(end:-1:1,end:-1:1,end:-1:1);
|
|
fsym(:,i) = symm_dir(f_signsym,T,order);
|
|
end
|
|
end
|
|
ERR = 0;
|
|
end
|
|
|
|
if ERR == 1, fprintf('*** ERROR: Feature dimension and T, order incompatible. ***\n'), end
|
|
function As = symm_dir(A,T,order)
|
|
% Symmetry marginalization routine. The purpose is to reduce the feature
|
|
% dimensionality and make the features more populated.
|
|
% A is an array of features of size (2*T+1)^order, otherwise error is outputted.
|
|
%
|
|
% Directional marginalization pertains to the fact that the
|
|
% differences d1, d2, d3, ... in a natural (both cover and stego) image
|
|
% are as likely to occur as ..., d3, d2, d1.
|
|
|
|
% Basically, we merge all pairs of bins (i,j,k, ...) and (..., k,j,i)
|
|
% as long as they are two different bins. Thus, instead of dim =
|
|
% (2T+1)^order, we decrease the dim by 1/2*{# of non-symmetric bins}.
|
|
% For order = 3, the reduced dim is (2T+1)^order - T(2T+1)^(order-1),
|
|
% for order = 4, it is (2T+1)^4 - 2T(T+1)(2T+1)^2.
|
|
|
|
B = 2*T+1;
|
|
done = zeros(size(A));
|
|
switch order
|
|
case 3, red = B^3 - T*B^2; % Dim of the marginalized set is (2T+1)^3-T*(2T+1)^2
|
|
case 4, red = B^4 - 2*T*(T+1)*B^2; % Dim of the marginalized set is (2T+1)^4-2T*(T+1)*(2T+1)^2
|
|
case 5, red = B^5 - 2*T*(T+1)*B^3;
|
|
end
|
|
As = zeros(red, 1);
|
|
m = 1;
|
|
|
|
switch order
|
|
case 3
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
if k ~= i % Asymmetric bin
|
|
if done(i+T+1,j+T+1,k+T+1) == 0
|
|
As(m) = A(i+T+1,j+T+1,k+T+1) + A(k+T+1,j+T+1,i+T+1); % Two mirror-bins are merged
|
|
done(i+T+1,j+T+1,k+T+1) = 1;
|
|
done(k+T+1,j+T+1,i+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
else % Symmetric bin is just copied
|
|
As(m) = A(i+T+1,j+T+1,k+T+1);
|
|
done(i+T+1,j+T+1,k+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
case 4
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
for n = -T : T
|
|
if (i ~= n) || (j ~= k) % Asymmetric bin
|
|
if done(i+T+1,j+T+1,k+T+1,n+T+1) == 0
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,n+T+1) + A(n+T+1,k+T+1,j+T+1,i+T+1); % Two mirror-bins are merged
|
|
done(i+T+1,j+T+1,k+T+1,n+T+1) = 1;
|
|
done(n+T+1,k+T+1,j+T+1,i+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
else % Symmetric bin is just copied
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,n+T+1);
|
|
done(i+T+1,j+T+1,k+T+1,n+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
case 5
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
for l = -T : T
|
|
for n = -T : T
|
|
if (i ~= n) || (j ~= l) % Asymmetric bin
|
|
if done(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) == 0
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) + A(n+T+1,l+T+1,k+T+1,j+T+1,i+T+1); % Two mirror-bins are merged
|
|
done(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) = 1;
|
|
done(n+T+1,l+T+1,k+T+1,j+T+1,i+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
else % Symmetric bin is just copied
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1);
|
|
done(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) = 1;
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
otherwise
|
|
fprintf('*** ERROR: Order not equal to 3 or 4 or 5! ***\n')
|
|
end
|
|
function fsym = symm1(f,T,order)
|
|
% Marginalization by sign and directional symmetry for a feature vector
|
|
% stored as a (2T+1)^order-dimensional array. The input feature f is
|
|
% assumed to be a dim x database_size matrix of features stored as columns.
|
|
|
|
[dim,N] = size(f);
|
|
B = 2*T+1;
|
|
c = B^order;
|
|
ERR = 1;
|
|
|
|
if dim == c
|
|
ERR = 0;
|
|
switch order % Reduced dimensionality for a c-dimensional feature vector
|
|
case 1, red = T + 1;
|
|
case 2, red = (T + 1)^2;
|
|
case 3, red = 1 + 3*T + 4*T^2 + 2*T^3;
|
|
case 4, red = B^2 + 4*T^2*(T + 1)^2;
|
|
case 5, red = 1/4*(B^2 + 1)*(B^3 + 1);
|
|
end
|
|
fsym = zeros(red,N);
|
|
|
|
for i = 1 : N
|
|
switch order
|
|
case 1, cube = f(:,i);
|
|
case 2, cube = reshape(f(:,i),[B B]);
|
|
case 3, cube = reshape(f(:,i),[B B B]);
|
|
case 4, cube = reshape(f(:,i),[B B B B]);
|
|
case 5, cube = reshape(f(:,i),[B B B B B]);
|
|
end
|
|
% [size(symm_dir(cube,T,order)) red]
|
|
fsym(:,i) = symm(cube,T,order);
|
|
end
|
|
end
|
|
|
|
if ERR == 1, fprintf('*** ERROR in symm1: Feature dimension and T, order incompatible. ***\n'), end
|
|
function As = symm(A,T,order)
|
|
% Symmetry marginalization routine. The purpose is to reduce the feature
|
|
% dimensionality and make the features more populated. It can be applied to
|
|
% 1D -- 5D co-occurrence matrices (order \in {1,2,3,4,5}) with sign and
|
|
% directional symmetries (explained below).
|
|
% A must be an array of (2*T+1)^order, otherwise error is outputted.
|
|
%
|
|
% Marginalization by symmetry pertains to the fact that, fundamentally,
|
|
% the differences between consecutive pixels in a natural image (both cover
|
|
% and stego) d1, d2, d3, ..., have the same probability of occurrence as the
|
|
% triple -d1, -d2, -d3, ...
|
|
%
|
|
% Directional marginalization pertains to the fact that the
|
|
% differences d1, d2, d3, ... in a natural (cover and stego) image are as
|
|
% likely to occur as ..., d3, d2, d1.
|
|
|
|
ERR = 1; % Flag denoting when size of A is incompatible with the input parameters T and order
|
|
m = 2;
|
|
B = 2*T + 1;
|
|
|
|
switch order
|
|
case 1 % First-order coocs are only symmetrized
|
|
if numel(A) == 2*T+1
|
|
As(1) = A(T+1); % The only non-marginalized bin is the origin 0
|
|
As(2:T+1) = A(1:T) + A(T+2:end);
|
|
As = As(:);
|
|
ERR = 0;
|
|
end
|
|
case 2
|
|
if numel(A) == (2*T+1)^2
|
|
As = zeros((T+1)^2, 1);
|
|
As(1) = A(T+1,T+1); % The only non-marginalized bin is the origin (0,0)
|
|
for i = -T : T
|
|
for j = -T : T
|
|
if (done(i+T+1,j+T+1) == 0) && (abs(i)+abs(j) ~= 0)
|
|
As(m) = A(i+T+1,j+T+1) + A(T+1-i,T+1-j);
|
|
done(i+T+1,j+T+1) = 1;
|
|
done(T+1-i,T+1-j) = 1;
|
|
if (i ~= j) && (done(j+T+1,i+T+1) == 0)
|
|
As(m) = As(m) + A(j+T+1,i+T+1) + A(T+1-j,T+1-i);
|
|
done(j+T+1,i+T+1) = 1;
|
|
done(T+1-j,T+1-i) = 1;
|
|
end
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
ERR = 0;
|
|
end
|
|
case 3
|
|
if numel(A) == B^3
|
|
done = zeros(size(A));
|
|
As = zeros(1+3*T+4*T^2+2*T^3, 1);
|
|
As(1) = A(T+1,T+1,T+1); % The only non-marginalized bin is the origin (0,0,0)
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
if (done(i+T+1,j+T+1,k+T+1) == 0) && (abs(i)+abs(j)+abs(k) ~= 0)
|
|
As(m) = A(i+T+1,j+T+1,k+T+1) + A(T+1-i,T+1-j,T+1-k);
|
|
done(i+T+1,j+T+1,k+T+1) = 1;
|
|
done(T+1-i,T+1-j,T+1-k) = 1;
|
|
if (i ~= k) && (done(k+T+1,j+T+1,i+T+1) == 0)
|
|
As(m) = As(m) + A(k+T+1,j+T+1,i+T+1) + A(T+1-k,T+1-j,T+1-i);
|
|
done(k+T+1,j+T+1,i+T+1) = 1;
|
|
done(T+1-k,T+1-j,T+1-i) = 1;
|
|
end
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
ERR = 0;
|
|
end
|
|
case 4
|
|
if numel(A) == (2*T+1)^4
|
|
done = zeros(size(A));
|
|
As = zeros(B^2 + 4*T^2*(T+1)^2, 1);
|
|
As(1) = A(T+1,T+1,T+1,T+1); % The only non-marginalized bin is the origin (0,0,0,0)
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
for n = -T : T
|
|
if (done(i+T+1,j+T+1,k+T+1,n+T+1) == 0) && (abs(i)+abs(j)+abs(k)+abs(n)~=0)
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,n+T+1) + A(T+1-i,T+1-j,T+1-k,T+1-n);
|
|
done(i+T+1,j+T+1,k+T+1,n+T+1) = 1;
|
|
done(T+1-i,T+1-j,T+1-k,T+1-n) = 1;
|
|
if ((i ~= n) || (j ~= k)) && (done(n+T+1,k+T+1,j+T+1,i+T+1) == 0)
|
|
As(m) = As(m) + A(n+T+1,k+T+1,j+T+1,i+T+1) + A(T+1-n,T+1-k,T+1-j,T+1-i);
|
|
done(n+T+1,k+T+1,j+T+1,i+T+1) = 1;
|
|
done(T+1-n,T+1-k,T+1-j,T+1-i) = 1;
|
|
end
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
ERR = 0;
|
|
end
|
|
case 5
|
|
if numel(A) == (2*T+1)^5
|
|
done = zeros(size(A));
|
|
As = zeros(1/4*(B^2 + 1)*(B^3 + 1), 1);
|
|
As(1) = A(T+1,T+1,T+1,T+1,T+1); % The only non-marginalized bin is the origin (0,0,0,0,0)
|
|
for i = -T : T
|
|
for j = -T : T
|
|
for k = -T : T
|
|
for l = -T : T
|
|
for n = -T : T
|
|
if (done(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) == 0) && (abs(i)+abs(j)+abs(k)+abs(l)+abs(n)~=0)
|
|
As(m) = A(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) + A(T+1-i,T+1-j,T+1-k,T+1-l,T+1-n);
|
|
done(i+T+1,j+T+1,k+T+1,l+T+1,n+T+1) = 1;
|
|
done(T+1-i,T+1-j,T+1-k,T+1-l,T+1-n) = 1;
|
|
if ((i ~= n) || (j ~= l)) && (done(n+T+1,l+T+1,k+T+1,j+T+1,i+T+1) == 0)
|
|
As(m) = As(m) + A(n+T+1,l+T+1,k+T+1,j+T+1,i+T+1) + A(T+1-n,T+1-l,T+1-k,T+1-j,T+1-i);
|
|
done(n+T+1,l+T+1,k+T+1,j+T+1,i+T+1) = 1;
|
|
done(T+1-n,T+1-l,T+1-k,T+1-j,T+1-i) = 1;
|
|
end
|
|
m = m + 1;
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
end
|
|
ERR = 0;
|
|
end
|
|
otherwise
|
|
As = [];
|
|
fprintf(' Order of cooc is not in {1,2,3,4,5}.\n')
|
|
end
|
|
|
|
if ERR == 1
|
|
As = [];
|
|
fprintf('*** ERROR in symm: The number of elements in the array is not (2T+1)^order. ***\n')
|
|
end
|
|
|
|
As = As(:);
|
|
function D = Residual(X,order,type)
|
|
% Computes the noise residual of a given type and order from MxN image X.
|
|
% residual order \in {1,2,3,4,5,6}
|
|
% type \in {hor,ver,diag,mdiag,KB,edge-h,edge-v,edge-d,edge-m}
|
|
% The resulting residual is an (M-b)x(N-b) array of the specified order,
|
|
% where b = ceil(order/2). This cropping is little more than it needs to
|
|
% be to make sure all the residuals are easily "synchronized".
|
|
% !!!!!!!!!!!!! Use order = 2 with KB and all edge residuals !!!!!!!!!!!!!
|
|
|
|
[M N] = size(X);
|
|
I = 1+ceil(order/2) : M-ceil(order/2);
|
|
J = 1+ceil(order/2) : N-ceil(order/2);
|
|
|
|
switch type
|
|
case 'hor'
|
|
switch order
|
|
case 1, D = - X(I,J) + X(I,J+1);
|
|
case 2, D = X(I,J-1) - 2*X(I,J) + X(I,J+1);
|
|
case 3, D = X(I,J-1) - 3*X(I,J) + 3*X(I,J+1) - X(I,J+2);
|
|
case 4, D = -X(I,J-2) + 4*X(I,J-1) - 6*X(I,J) + 4*X(I,J+1) - X(I,J+2);
|
|
case 5, D = -X(I,J-2) + 5*X(I,J-1) - 10*X(I,J) + 10*X(I,J+1) - 5*X(I,J+2) + X(I,J+3);
|
|
case 6, D = X(I,J-3) - 6*X(I,J-2) + 15*X(I,J-1) - 20*X(I,J) + 15*X(I,J+1) - 6*X(I,J+2) + X(I,J+3);
|
|
end
|
|
case 'ver'
|
|
switch order
|
|
case 1, D = - X(I,J) + X(I+1,J);
|
|
case 2, D = X(I-1,J) - 2*X(I,J) + X(I+1,J);
|
|
case 3, D = X(I-1,J) - 3*X(I,J) + 3*X(I+1,J) - X(I+2,J);
|
|
case 4, D = -X(I-2,J) + 4*X(I-1,J) - 6*X(I,J) + 4*X(I+1,J) - X(I+2,J);
|
|
case 5, D = -X(I-2,J) + 5*X(I-1,J) - 10*X(I,J) + 10*X(I+1,J) - 5*X(I+2,J) + X(I+3,J);
|
|
case 6, D = X(I-3,J) - 6*X(I-2,J) + 15*X(I-1,J) - 20*X(I,J) + 15*X(I+1,J) - 6*X(I+2,J) + X(I+3,J);
|
|
end
|
|
case 'diag'
|
|
switch order
|
|
case 1, D = - X(I,J) + X(I+1,J+1);
|
|
case 2, D = X(I-1,J-1) - 2*X(I,J) + X(I+1,J+1);
|
|
case 3, D = X(I-1,J-1) - 3*X(I,J) + 3*X(I+1,J+1) - X(I+2,J+2);
|
|
case 4, D = -X(I-2,J-2) + 4*X(I-1,J-1) - 6*X(I,J) + 4*X(I+1,J+1) - X(I+2,J+2);
|
|
case 5, D = -X(I-2,J-2) + 5*X(I-1,J-1) - 10*X(I,J) + 10*X(I+1,J+1) - 5*X(I+2,J+2) + X(I+3,J+3);
|
|
case 6, D = X(I-3,J-3) - 6*X(I-2,J-2) + 15*X(I-1,J-1) - 20*X(I,J) + 15*X(I+1,J+1) - 6*X(I+2,J+2) + X(I+3,J+3);
|
|
end
|
|
case 'mdiag'
|
|
switch order
|
|
case 1, D = - X(I,J) + X(I-1,J+1);
|
|
case 2, D = X(I-1,J+1) - 2*X(I,J) + X(I+1,J-1);
|
|
case 3, D = X(I-1,J+1) - 3*X(I,J) + 3*X(I+1,J-1) - X(I+2,J-2);
|
|
case 4, D = -X(I-2,J+2) + 4*X(I-1,J+1) - 6*X(I,J) + 4*X(I+1,J-1) - X(I+2,J-2);
|
|
case 5, D = -X(I-2,J+2) + 5*X(I-1,J+1) - 10*X(I,J) + 10*X(I+1,J-1) - 5*X(I+2,J-2) + X(I+3,J-3);
|
|
case 6, D = X(I-3,J+3) - 6*X(I-2,J+2) + 15*X(I-1,J+1) - 20*X(I,J) + 15*X(I+1,J-1) - 6*X(I+2,J-2) + X(I+3,J-3);
|
|
end
|
|
case 'KB'
|
|
D = -X(I-1,J-1) + 2*X(I-1,J) - X(I-1,J+1) + 2*X(I,J-1) - 4*X(I,J) + 2*X(I,J+1) - X(I+1,J-1) + 2*X(I+1,J) - X(I+1,J+1);
|
|
case 'edge-h'
|
|
Du = 2*X(I-1,J) + 2*X(I,J-1) + 2*X(I,J+1) - X(I-1,J-1) - X(I-1,J+1) - 4*X(I,J); % -1 2 -1
|
|
Db = 2*X(I+1,J) + 2*X(I,J-1) + 2*X(I,J+1) - X(I+1,J-1) - X(I+1,J+1) - 4*X(I,J); % 2 C 2 + flipped vertically
|
|
D = [Du,Db];
|
|
case 'edge-v'
|
|
Dl = 2*X(I,J-1) + 2*X(I-1,J) + 2*X(I+1,J) - X(I-1,J-1) - X(I+1,J-1) - 4*X(I,J); % -1 2
|
|
Dr = 2*X(I,J+1) + 2*X(I-1,J) + 2*X(I+1,J) - X(I-1,J+1) - X(I+1,J+1) - 4*X(I,J); % 2 C + flipped horizontally
|
|
D = [Dl,Dr]; % -1 2
|
|
case 'edge-m'
|
|
Dlu = 2*X(I,J-1) + 2*X(I-1,J) - X(I-1,J-1) - X(I+1,J-1) - X(I-1,J+1) - X(I,J); % -1 2 -1
|
|
Drb = 2*X(I,J+1) + 2*X(I+1,J) - X(I+1,J+1) - X(I+1,J-1) - X(I-1,J+1) - X(I,J); % 2 C + flipped mdiag
|
|
D = [Dlu,Drb]; % -1
|
|
case 'edge-d'
|
|
Dru = 2*X(I-1,J) + 2*X(I,J+1) - X(I-1,J+1) - X(I-1,J-1) - X(I+1,J+1) - X(I,J); % -1 2 -1
|
|
Dlb = 2*X(I,J-1) + 2*X(I+1,J) - X(I+1,J-1) - X(I+1,J+1) - X(I-1,J-1) - X(I,J); % C 2 + flipped diag
|
|
D = [Dru,Dlb]; % -1
|
|
case 'KV'
|
|
D = 8*X(I-1,J) + 8*X(I+1,J) + 8*X(I,J-1) + 8*X(I,J+1);
|
|
D = D - 6*X(I-1,J+1) - 6*X(I-1,J-1) - 6*X(I+1,J-1) - 6*X(I+1,J+1);
|
|
D = D - 2*X(I-2,J) - 2*X(I+2,J) - 2*X(I,J+2) - 2*X(I,J-2);
|
|
D = D + 2*X(I-1,J-2) + 2*X(I-2,J-1) + 2*X(I-2,J+1) + 2*X(I-1,J+2) + 2*X(I+1,J+2) + 2*X(I+2,J+1) + 2*X(I+2,J-1) + 2*X(I+1,J-2);
|
|
D = D - X(I-2,J-2) - X(I-2,J+2) - X(I+2,J-2) - X(I+2,J+2) - 12*X(I,J);
|
|
end
|