Procedure is as follows:
Void CMyDlg: : OnHuiDuGongSheng ()
{
//TODO: Add your command handler code here

EngEvalString (ep, "Gray=imread (' 121. BMP);" );
EngEvalString (ep, "[M, N, O]=size (Gray);" );
EngEvalString (ep, M=128; "" );
EngEvalString (ep, N=128; "" );

//in order to reduce the amount of calculation, the original image is compressed grayscale, quantification into Gray level 16
EngEvalString (ep, "for I=1: M; \
For j=1: N. \
For n=25 6/16; \
If (n - 1) * 16 & lt;=Gray (I, j) & amp; & Gray (I, j) & lt;=(n - 1) * 16 + 15; \
Gray (I, j)=n - 1; \
end; \
end; \
end; \
End." );

Four co-occurrence matrix P//calculation, take a distance of 1, Angle respectively 0,45,90,135
EngEvalString (ep, "P=zeros (16,16,4);" );
EngEvalString (ep, "the for m=1:16; \
For n=1:16; \
For I=1: M; \
For j=1: N. \
If jP (m, n, 1)=P (m, n, 1) + 1; \
P (n, m, 1)=P (m, n, 1); \
end; \
If i> 1& & jP (m, n, 2)=P (m, n, 2) + 1; \
P (n, m, 2)=P (m, n, 2); \
end; \
If iP (m, n, 3)=P (m, n, 3) + 1; \
P (n, m, 3)=P (m, n, 3); \
end; \
If iP (m, n, 4)=P (m, n, 4) + 1; \
P (n, m, 4)=P (m, n, 4); \
end; \
end; \
end; \
If m==n; \
P (m, n, :)=P (m, n, :) * 2; \
end; \
end; \
End." );

//the symbiotic matrix normalized
EngEvalString (ep, "for n=1:4; \
P (:, :, n)=P (:, :, n)/sum (the sum (P (:, :, n))); \
End." );
//the symbiotic matrix to calculate energy, entropy, moment of inertia, four texture parameters related to
EngEvalString (ep, "H=zeros (1, 4);" );
EngEvalString (ep, I=H; "" );
EngEvalString (ep, Ux=H; "" );
EngEvalString (ep, "Uy=H;" );
EngEvalString (ep, "deltaX=H;" );
EngEvalString (ep, deltaY=H; "" );
EngEvalString (ep, C=H; "" );
EngEvalString (ep, "for n=1:4; \
E (n)=sum (the sum (P (:, :, n). ^ 2)); \
For I=1:16. \
For j=1:16; \
If P (I, j, n) ~=0; \
H (n)=P (I, j, n) * log (P (I, j, n)) + H (n); \
end; \
I (n)=(I - j) ^ 2 * P + I (I, j, n) (n); \
Ux (n)=I * P (I, j, n) + Ux (n); \
Uy (n)=j * P (I, j, n) + Uy (n); \
end; \
end; \
End." );
EngEvalString (ep, "for n=1:4; \
For I=1:16. \
For j=1:16; \
DeltaX (n)=(I - Ux (n)) ^ 2 * P (I, j, n) + deltaX (n); \
DeltaY (n)=(j - Uy (n)) ^ 2 * P (I, j, n) + deltaY (n); \
C (n)=I * j * P (I, j, n) + C (n); \
end; \
end; \
C (n)=(C (n) - Ux Uy (n) * (n))/deltaX (n)/deltaY (n); \
End." );

//energy, entropy, mean and standard deviation of the moment of inertia, relevant as the final 8 d texture feature
EngEvalString (ep, "T (1)=mean (E);" );
EngEvalString (ep, "T (2)=mean (H);" );
EngEvalString (ep, "T (3)=mean (I);" );
EngEvalString (ep, "T (4)=mean (C);" );
//engEvalString (ep, "fprintf (" % f, % f, % f, % f", T (1), T (2), T (3), T (4));" );
EngEvalString (ep, "xlswrite (' testdata. XLS, [T (1) (2) T T T (4) and (3)], 'B2: H2);" );