Sentinel-1影像数据陆地掩模与dB值转换(IDL Programming)

;+
; :Author: Baikal
;-
; Request: .zip format source file should be stored in the same directory
; with the processed folder using the procSentinelRTC_recipe.py script.
;
; Purpose:
; Mask land areas and convert pixel values to dB.
; -Version: 1.0
; Time: Sep.24,2017
;
; First Modification:
; Purpose:
; To sovle abnormal display of pixel values in elevation file and HH&VV output img.
; Time : Sept.25,2017
;
; Second Modification:
; Purpose:
; To clip the image based on Mask_Sea file rather than elevation image.
; -Version: 2.0
; Time : Oct.20,2017
; Third Modification:
; Purpose:
; To add modification of header file.
; -Version: 2.0
; Time : Oct.21,2017


PRO Mask_Land_To_dB_V2

  COMPILE_OPT idl2

  ;给定文件路径
  directory = 'C:\Users\Baikal\Desktop\test\'
  ;查找zip文件
  file = file_search(directory,count = num,'*.zip')
  for id = 0,num-1 DO BEGIN
    ;获取文件名
    file_name = strmid(file[id],70,67,/Reverse_Offset)

    ;获取HH极化与HV极化影像路径
    Image_HH = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\" + "Beta0_HH.img"
    Image_HV = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\" + "Beta0_HV.img"

    ;获取HH极化与HV极化影像头文件路径
    HH_hdr = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\" + "Beta0_HH.hdr"
    HV_hdr = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\" + "Beta0_HV.hdr"

    ;获取HH极化影像的掩膜文件路径（掩膜文件中仅含有陆地区域）
    Image_HH_Mask = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC(Mask_Sea).data\" + "Beta0_HH.img"

    ;查找影像中是否有已经存在的转为dB值的文件，若有，则删除
    ;    HH_file = file_search(directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\",'Aim_Beta0_HH.hdr')
    ;    IF HH_file THEN  FILE_DELETE,HH_file
    ;    HV_file = file_search(directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\",'Aim_Beta0_HV.hdr')
    ;    IF HV_file THEN  FILE_DELETE,HV_file
    ;    HH_HV_file = file_search(directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\",'Aim_Beta0_HH_HV.hdr')
    ;    IF HH_HV_file THEN  FILE_DELETE,HH_HV_file

    ;对HH、HV、HH/HV极化影像的三个波段进行波段合成
    header_path = directory + file_name + '\' + file_name + "_TNR_OB_CAL_TC.data\"
    FILE_COPY, HH_hdr, header_path + 'Aim_Beta0_HH.hdr'
    FILE_COPY, HV_hdr, header_path + 'Aim_Beta0_HV.hdr'
    FILE_COPY, HV_hdr, header_path + 'Aim_Beta0_HH_HV.hdr'
    Aim_HH_hdr    = header_path + 'Aim_Beta0_HH.hdr'
    Aim_HV_hdr    = header_path + 'Aim_Beta0_HV.hdr'
    Aim_HH_HV_hdr = header_path + 'Aim_Beta0_HH_HV.hdr'
    ;复制生成转为dB值文件的头文件
    Output_Directory = directory + file_name + '\Composition\'
    file_mkdir, Output_Directory
    FILE_COPY, HV_hdr, Output_Directory + 'Layer Stacking.hdr'
    LS_hdr = Output_Directory + 'Layer Stacking.hdr'

    ;将头文件中的内容进行修改
    OPENU,hdr_lun,Aim_HH_hdr,/get_lun
    line    = 'band names = { Beta0_HH }'
    replace = 'band names={Aim_Beta0_HH}'
    Skip_Lun,hdr_lun,10,/lines
    Point_lun,-hdr_lun,line
    printf,hdr_lun,replace
    CLOSE,hdr_lun & FREE_LUN, hdr_lun

    OPENU,hdr_lun,Aim_HV_hdr,/get_lun
    line    = 'band names = { Beta0_HV }'
    replace = 'band names={Aim_Beta0_HV}'
    Skip_Lun,hdr_lun,10,/lines 
    Point_lun,-hdr_lun,line
    printf,hdr_lun,replace
    CLOSE,hdr_lun & FREE_LUN, hdr_lun

    OPENU,hdr_lun,Aim_HH_HV_hdr,/get_lun
    line    = 'band names = { Beta0_HV }'
    replace = 'band names={Beta0_HH_HV}'
    Skip_Lun,hdr_lun,10,/lines
    Point_lun,-hdr_lun,line
    printf,hdr_lun,replace
    CLOSE,hdr_lun & FREE_LUN, hdr_lun

    OPENU,hdr_lun,LS_hdr,/get_lun
    line_Band  = 'bands = 1'
    replace_line_Band = 'bands = 3'
    Skip_Lun,hdr_lun,4,/lines
    Point_lun,-hdr_lun,line_Band
    printf,hdr_lun,replace_line_Band
    CLOSE,hdr_lun & FREE_LUN, hdr_lun
   
    OPENU,hdr_lun,LS_hdr,/get_lun
    line_Band_Name    = 'band names = { Beta0_HV }'
    replace_Band_Name = 'band names={HH,HV,HH_HV}'
    Skip_Lun,hdr_lun,10,/lines
    Point_lun,-hdr_lun,line_Band_Name
    printf,hdr_lun,replace_Band_Name
    CLOSE,hdr_lun & FREE_LUN, hdr_lun

    ;打开HH极化的头文件
    OPENR,hdr_lun,HH_hdr,/get_lun
    ;设立一个存储头文件行数值的变量
    line = ''
    WHILE (~EOF(hdr_lun)) DO BEGIN
      ;读取一行数据
      READF,hdr_lun,line
      tmp = strsplit(line[0],'=',/extract)
      ;利用空格( )获得影像的头文件参数信息
      header_keyword = Strsplit(tmp[0], ' ', /extract)
      ;按照头文件中变量名字获得相应的参数
      IF header_keyword[0] EQ 'samples' THEN xs = LONG(tmp[1])
      IF header_keyword[0] EQ 'lines' THEN ys = LONG(tmp[1])
      IF header_keyword[0] EQ 'data type' THEN type = LONG(tmp[1])
    ENDWHILE
    ;关闭并释放头文件指针
    CLOSE,hdr_lun & FREE_LUN, hdr_lun

    ;读取HH极化影像数据
    OPENR,HH_unit,Image_HH,/get_lun
    HH = MAKE_ARRAY(xs,ys,type = type)
    READU,HH_unit,HH
    ;!!!修改影像读取的存储方式
    BYTEORDER,HH,/FTOXDR

    ;读取HV极化影像数据
    OPENR,HV_unit,Image_HV,/get_lun
    HV = MAKE_ARRAY(xs,ys,type = type)
    READU,HV_unit,HV
    BYTEORDER,HV,/FTOXDR

    ;打开掩膜文件的头文件
    OPENR,HH_Mask_unit,Image_HH_Mask,/get_lun
    HH_Mask = MAKE_ARRAY(xs,ys,type = type)
    READU,HH_Mask_unit,HH_Mask
    BYTEORDER,HH_Mask,/FTOXDR

    ;生成HH/HV波段数据
    HH_HV = MAKE_ARRAY(xs,ys,type = type)
    ；修改影像的存储方式
    BYTEORDER,HH_HV,/FTOXDR

    for i = 0,xs-1 DO BEGIN
      for j = 0,ys-1 DO BEGIN
        ;如果掩膜文件为非零值或者HH极化数据为零值(背景值)，则将输出的文件设置为255
        IF HH_Mask[i,j] NE 0 OR HH[i,j] EQ 0 THEN BEGIN
          HH[i,j]    =  255
          HV[i,j]    =  255
          HH_HV[i,j] =  255
        ENDIF ELSE BEGIN
          HH[i,j]    = 10*ALOG10(HH[i,j])
          HV[i,j]    = 10*ALOG10(HV[i,j])
          HH_HV[i,j] = HH[i,j]/HV[i,j]
        ENDELSE
      endfor
    endfor
    FREE_LUN, HH_unit
    FREE_LUN, HV_unit
    FREE_LUN, HH_Mask_unit

   
    ;!!!!!将文件的存储方式进行修改
    BYTEORDER,HH,/XDRTOF
    BYTEORDER,HV,/XDRTOF
    BYTEORDER,HH_HV,/XDRTOF

    ;检查程序是否存在数学错误
    IF check_math() THEN Print, check_math()
    Beta0_HH_Img = header_path + 'Aim_Beta0_HH.img'
    OPENW,Beta0_HH_lun,Beta0_HH_Img,/get_lun
    ;生成HH极化转为dB值文件的的头文件信息
    WriteU,Beta0_HH_lun,HH
    Close,Beta0_HH_lun
    Free_LUN,Beta0_HH_lun

    Beta0_HV_Img = header_path + 'Aim_Beta0_HV.img'
    OPENW,Beta0_HV_lun,Beta0_HV_Img,/get_lun
    WriteU,Beta0_HV_lun,HV
    Close,Beta0_HV_lun
    Free_LUN,Beta0_HV_lun

    Beta0_HH_HV_Img = header_path + 'Aim_Beta0_HH_HV.img'
    OPENW,Beta0_HH_HV_lun,Beta0_HH_HV_Img,/get_lun
    WriteU,Beta0_HH_HV_lun,HH_HV
    Close,Beta0_HH_HV_lun
    Free_LUN,Beta0_HH_HV_lun

    ;HH、HV、HH_HV影像3个波段进行波段合成
    Aim_Beta0_HH_HV = [[HH],[HV],[HH_HV]]
    Aim_Beta0_HH_HV = Reform(Aim_Beta0_HH_HV,[3,xs,ys])
    Output_File = Output_Directory + 'Layer Stacking.img'
    OPENW,lun,Output_File,/get_lun
    WriteU,lun,Aim_Beta0_HH_HV
    Close,lun
    Free_Lun,lun

    ;输出显示处理状态
    Print,'还剩',strtrim(num-id-1,1),'景影像待处理，已处理',strtrim(id+1,1),'景影像'
  ENDfor
  print,'处理完毕'
END

move cursor in DOS/Windows command line（DOS/Windows窗口快捷鍵）

参考链接：http://www.contrib.andrew.cmu.edu/~minlix/wordpress/?p=238

word公式编辑

1.       打开视图—标尺；

2.       “开始”选项卡下，“样式”栏选择“创建样式”，点击修改；

3.       选择格式，点击制表位；

4.       输入39.5，选择右对齐，前导符无，点击设置；

5.       再输入17.95，居中，前导符无，设置，确定；

6.       编辑一个公式，后面编号，选中我们的编号公式格式；

7.       在公式前键入一个Tab，编号前键入一个Tab，biu～自动就好了

实例：

                                           F=ma                                                        (1)

Textural Features of Gray Level Co-occurrence Matrix (IDL)

使用方法参考链接：
http://www.harrisgeospatial.com/docs/TEXTURE_COOCCUR_DOIT.html

GLCM with Sliding Window （带有滑动窗口的GLCM程序,Matlab)

具有滑动窗口的灰度共生矩阵Matlab代码

来源:Github

代码下载链接
https://github.com/palmerc/GLCM

其中，Sliding GLCM.m为主程序，其他程序（GLCM.m，PadMatrix.m，QuadrantEnergy.m以及Quadrants.m）均为子程序。

Calculation of texture features of GLCM(MATLab)

function [out] = GLCM_Features(glcmin,pairs)
% 
% GLCM_Features1 helps to calculate the features from the different GLCMs
% that are input to the function. The GLCMs are stored in a i x j x n
% matrix, where n is the number of GLCMs calculated usually due to the
% different orientation and displacements used in the algorithm. Usually
% the values i and j are equal to 'NumLevels' parameter of the GLCM
% computing function graycomatrix(). Note that matlab quantization values
% belong to the set {1,..., NumLevels} and not from {0,...,(NumLevels-1)}
% as provided in some references
% http://www.mathworks.com/access/helpdesk/help/toolbox/images/graycomatrix
% .html
% 
% Although there is a function graycoprops() in Matlab Image Processing
% Toolbox that computes four parameters Contrast, Correlation, Energy,
% and Homogeneity. The paper by Haralick suggests a few more parameters
% that are also computed here. The code is not fully vectorized and hence
% is not an efficient implementation but it is easy to add new features
% based on the GLCM using this code. Takes care of 3 dimensional glcms
% (multiple glcms in a single 3D array)
% 
% If you find that the values obtained are different from what you expect 
% or if you think there is a different formula that needs to be used 
% from the ones used in this code please let me know. 
% A few questions which I have are listed in the link 
% http://www.mathworks.com/matlabcentral/newsreader/view_thread/239608
%
% I plan to submit a vectorized version of the code later and provide 
% updates based on replies to the above link and this initial code. 
%
% Features computed 
% Autocorrelation: [2]                      (out.autoc)
% Contrast: matlab/[1,2]                    (out.contr)
% Correlation: matlab                       (out.corrm)
% Correlation: [1,2]                        (out.corrp)
% Cluster Prominence: [2]                   (out.cprom)
% Cluster Shade: [2]                        (out.cshad)
% Dissimilarity: [2]                        (out.dissi)
% Energy: matlab / [1,2]                    (out.energ)
% Entropy: [2]                              (out.entro)
% Homogeneity: matlab                       (out.homom)
% Homogeneity: [2]                          (out.homop)
% Maximum probability: [2]                  (out.maxpr)
% Sum of sqaures: Variance [1]              (out.sosvh)
% Sum average [1]                           (out.savgh)
% Sum variance [1]                          (out.svarh)
% Sum entropy [1]                           (out.senth)
% Difference variance [1]                   (out.dvarh)
% Difference entropy [1]                    (out.denth)
% Information measure of correlation1 [1]   (out.inf1h)
% Informaiton measure of correlation2 [1]   (out.inf2h)
% Inverse difference (INV) is homom [3]     (out.homom)
% Inverse difference normalized (INN) [3]   (out.indnc) 
% Inverse difference moment normalized [3]  (out.idmnc)
%
% The maximal correlation coefficient was not calculated due to
% computational instability 
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
% Formulae from MATLAB site (some look different from
% the paper by Haralick but are equivalent and give same results)
% Example formulae: 
% Contrast = sum_i(sum_j(  (i-j)^2 * p(i,j) ) ) (same in matlab/paper)
% Correlation = sum_i( sum_j( (i - u_i)(j - u_j)p(i,j)/(s_i.s_j) ) ) (m)
% Correlation = sum_i( sum_j( ((ij)p(i,j) - u_x.u_y) / (s_x.s_y) ) ) (p[2])
% Energy = sum_i( sum_j( p(i,j)^2 ) )           (same in matlab/paper)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + |i-j|) ) ) (as in matlab)
% Homogeneity = sum_i( sum_j( p(i,j) / (1 + (i-j)^2) ) ) (as in paper)
% 
% Where:
% u_i = u_x = sum_i( sum_j( i.p(i,j) ) ) (in paper [2])
% u_j = u_y = sum_i( sum_j( j.p(i,j) ) ) (in paper [2])
% s_i = s_x = sum_i( sum_j( (i - u_x)^2.p(i,j) ) ) (in paper [2])
% s_j = s_y = sum_i( sum_j( (j - u_y)^2.p(i,j) ) ) (in paper [2])
%
% 
% Normalize the glcm:
% Compute the sum of all the values in each glcm in the array and divide 
% each element by it sum
%
% Haralick uses 'Symmetric' = true in computing the glcm
% There is no Symmetric flag in the Matlab version I use hence
% I add the diagonally opposite pairs to obtain the Haralick glcm
% Here it is assumed that the diagonally opposite orientations are paired
% one after the other in the matrix.
% If the above assumption is true with respect to the input glcm then
% setting the flag 'pairs' to 1 will compute the final glcms that would result 
% by setting 'Symmetric' to true. If your glcm is computed using the
% Matlab version with 'Symmetric' flag you can set the flag 'pairs' to 0
%
% References:
% 1. R. M. Haralick, K. Shanmugam, and I. Dinstein, Textural Features of
% Image Classification, IEEE Transactions on Systems, Man and Cybernetics,
% vol. SMC-3, no. 6, Nov. 1973
% 2. L. Soh and C. Tsatsoulis, Texture Analysis of SAR Sea Ice Imagery
% Using Gray Level Co-Occurrence Matrices, IEEE Transactions on Geoscience
% and Remote Sensing, vol. 37, no. 2, March 1999.
% 3. D A. Clausi, An analysis of co-occurrence texture statistics as a
% function of grey level quantization, Can. J. Remote Sensing, vol. 28, no.
% 1, pp. 45-62, 2002
% 4. http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
%
%
% Example:
%
% Usage is similar to graycoprops() but needs extra parameter 'pairs' apart
% from the GLCM as input
% I = imread('circuit.tif');
% GLCM2 = graycomatrix(I,'Offset',[2 0;0 2]);
% stats = GLCM_features1(GLCM2,0)
% The output is a structure containing all the parameters for the different
% GLCMs
%
% [Avinash Uppuluri: avinash_uv@yahoo.com: Last modified: 11/20/08]

% If 'pairs' not entered: set pairs to 0 
if ((nargin > 2) || (nargin == 0))
   error('Too many or too few input arguments. Enter GLCM and pairs.');
elseif ( (nargin == 2) ) 
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
        error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end    
elseif (nargin == 1) % only GLCM is entered
    pairs = 0; % default is numbers and input 1 for percentage
    if ((size(glcmin,1) <= 1) || (size(glcmin,2) <= 1))
       error('The GLCM should be a 2-D or 3-D matrix.');
    elseif ( size(glcmin,1) ~= size(glcmin,2) )
       error('Each GLCM should be square with NumLevels rows and NumLevels cols');
    end    
end


format long e
if (pairs == 1)
    newn = 1;
    for nglcm = 1:2:size(glcmin,3)
        glcm(:,:,newn)  = glcmin(:,:,nglcm) + glcmin(:,:,nglcm+1);
        newn = newn + 1;
    end
elseif (pairs == 0)
    glcm = glcmin;
end

size_glcm_1 = size(glcm,1);
size_glcm_2 = size(glcm,2);
size_glcm_3 = size(glcm,3);

% checked 
out.autoc = zeros(1,size_glcm_3); % Autocorrelation: [2] 
out.contr = zeros(1,size_glcm_3); % Contrast: matlab/[1,2]
out.corrm = zeros(1,size_glcm_3); % Correlation: matlab
out.corrp = zeros(1,size_glcm_3); % Correlation: [1,2]
out.cprom = zeros(1,size_glcm_3); % Cluster Prominence: [2]
out.cshad = zeros(1,size_glcm_3); % Cluster Shade: [2]
out.dissi = zeros(1,size_glcm_3); % Dissimilarity: [2]
out.energ = zeros(1,size_glcm_3); % Energy: matlab / [1,2]
out.entro = zeros(1,size_glcm_3); % Entropy: [2]
out.homom = zeros(1,size_glcm_3); % Homogeneity: matlab
out.homop = zeros(1,size_glcm_3); % Homogeneity: [2]
out.maxpr = zeros(1,size_glcm_3); % Maximum probability: [2]

out.sosvh = zeros(1,size_glcm_3); % Sum of sqaures: Variance [1]
out.savgh = zeros(1,size_glcm_3); % Sum average [1]
out.svarh = zeros(1,size_glcm_3); % Sum variance [1]
out.senth = zeros(1,size_glcm_3); % Sum entropy [1]
out.dvarh = zeros(1,size_glcm_3); % Difference variance [4]
%out.dvarh2 = zeros(1,size_glcm_3); % Difference variance [1]
out.denth = zeros(1,size_glcm_3); % Difference entropy [1]
out.inf1h = zeros(1,size_glcm_3); % Information measure of correlation1 [1]
out.inf2h = zeros(1,size_glcm_3); % Informaiton measure of correlation2 [1]
%out.mxcch = zeros(1,size_glcm_3);% maximal correlation coefficient [1]
%out.invdc = zeros(1,size_glcm_3);% Inverse difference (INV) is homom [3]
out.indnc = zeros(1,size_glcm_3); % Inverse difference normalized (INN) [3]
out.idmnc = zeros(1,size_glcm_3); % Inverse difference moment normalized [3]

% correlation with alternate definition of u and s
%out.corrm2 = zeros(1,size_glcm_3); % Correlation: matlab
%out.corrp2 = zeros(1,size_glcm_3); % Correlation: [1,2]

glcm_sum  = zeros(size_glcm_3,1);
glcm_mean = zeros(size_glcm_3,1);
glcm_var  = zeros(size_glcm_3,1);

% http://www.fp.ucalgary.ca/mhallbey/glcm_mean.htm confuses the range of 
% i and j used in calculating the means and standard deviations.
% As of now I am not sure if the range of i and j should be [1:Ng] or
% [0:Ng-1]. I am working on obtaining the values of mean and std that get
% the values of correlation that are provided by matlab.
u_x = zeros(size_glcm_3,1);
u_y = zeros(size_glcm_3,1);
s_x = zeros(size_glcm_3,1);
s_y = zeros(size_glcm_3,1);

% % alternate values of u and s
% u_x2 = zeros(size_glcm_3,1);
% u_y2 = zeros(size_glcm_3,1);
% s_x2 = zeros(size_glcm_3,1);
% s_y2 = zeros(size_glcm_3,1);

% checked p_x p_y p_xplusy p_xminusy
p_x = zeros(size_glcm_1,size_glcm_3); % Ng x #glcms[1]  
p_y = zeros(size_glcm_2,size_glcm_3); % Ng x #glcms[1]
p_xplusy = zeros((size_glcm_1*2 - 1),size_glcm_3); %[1]
p_xminusy = zeros((size_glcm_1),size_glcm_3); %[1]
% checked hxy hxy1 hxy2 hx hy
hxy  = zeros(size_glcm_3,1);
hxy1 = zeros(size_glcm_3,1);
hx   = zeros(size_glcm_3,1);
hy   = zeros(size_glcm_3,1);
hxy2 = zeros(size_glcm_3,1);

%Q    = zeros(size(glcm));

for k = 1:size_glcm_3 % number glcms

    glcm_sum(k) = sum(sum(glcm(:,:,k)));
    glcm(:,:,k) = glcm(:,:,k)./glcm_sum(k); % Normalize each glcm
    glcm_mean(k) = mean2(glcm(:,:,k)); % compute mean after norm
    glcm_var(k)  = (std2(glcm(:,:,k)))^2;
    
    for i = 1:size_glcm_1

        for j = 1:size_glcm_2

            out.contr(k) = out.contr(k) + (abs(i - j))^2.*glcm(i,j,k);
            out.dissi(k) = out.dissi(k) + (abs(i - j)*glcm(i,j,k));
            out.energ(k) = out.energ(k) + (glcm(i,j,k).^2);
            out.entro(k) = out.entro(k) - (glcm(i,j,k)*log(glcm(i,j,k) + eps));
            out.homom(k) = out.homom(k) + (glcm(i,j,k)/( 1 + abs(i-j) ));
            out.homop(k) = out.homop(k) + (glcm(i,j,k)/( 1 + (i - j)^2));
            % [1] explains sum of squares variance with a mean value;
            % the exact definition for mean has not been provided in 
            % the reference: I use the mean of the entire normalized glcm 
            out.sosvh(k) = out.sosvh(k) + glcm(i,j,k)*((i - glcm_mean(k))^2);
            
            %out.invdc(k) = out.homom(k);
            out.indnc(k) = out.indnc(k) + (glcm(i,j,k)/( 1 + (abs(i-j)/size_glcm_1) ));
            out.idmnc(k) = out.idmnc(k) + (glcm(i,j,k)/( 1 + ((i - j)/size_glcm_1)^2));
            u_x(k)          = u_x(k) + (i)*glcm(i,j,k); % changed 10/26/08
            u_y(k)          = u_y(k) + (j)*glcm(i,j,k); % changed 10/26/08
            % code requires that Nx = Ny 
            % the values of the grey levels range from 1 to (Ng) 
        end
        
    end
    out.maxpr(k) = max(max(glcm(:,:,k)));
end
% glcms have been normalized:
% The contrast has been computed for each glcm in the 3D matrix
% (tested) gives similar results to the matlab function

for k = 1:size_glcm_3
    
    for i = 1:size_glcm_1
        
        for j = 1:size_glcm_2
            p_x(i,k) = p_x(i,k) + glcm(i,j,k); 
            p_y(i,k) = p_y(i,k) + glcm(j,i,k); % taking i for j and j for i
            if (ismember((i + j),[2:2*size_glcm_1])) 
                p_xplusy((i+j)-1,k) = p_xplusy((i+j)-1,k) + glcm(i,j,k);
            end
            if (ismember(abs(i-j),[0:(size_glcm_1-1)])) 
                p_xminusy((abs(i-j))+1,k) = p_xminusy((abs(i-j))+1,k) +...
                    glcm(i,j,k);
            end
        end
    end
    
%     % consider u_x and u_y and s_x and s_y as means and standard deviations
%     % of p_x and p_y
%     u_x2(k) = mean(p_x(:,k));
%     u_y2(k) = mean(p_y(:,k));
%     s_x2(k) = std(p_x(:,k));
%     s_y2(k) = std(p_y(:,k));
    
end

% marginal probabilities are now available [1]
% p_xminusy has +1 in index for matlab (no 0 index)
% computing sum average, sum variance and sum entropy:
for k = 1:(size_glcm_3)
    
    for i = 1:(2*(size_glcm_1)-1)
        out.savgh(k) = out.savgh(k) + (i+1)*p_xplusy(i,k);
        % the summation for savgh is for i from 2 to 2*Ng hence (i+1)
        out.senth(k) = out.senth(k) - (p_xplusy(i,k)*log(p_xplusy(i,k) + eps));
    end

end
% compute sum variance with the help of sum entropy
for k = 1:(size_glcm_3)
    
    for i = 1:(2*(size_glcm_1)-1)
        out.svarh(k) = out.svarh(k) + (((i+1) - out.senth(k))^2)*p_xplusy(i,k);
        % the summation for savgh is for i from 2 to 2*Ng hence (i+1)
    end

end
% compute difference variance, difference entropy, 
for k = 1:size_glcm_3
% out.dvarh2(k) = var(p_xminusy(:,k));
% but using the formula in 
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
% we have for dvarh
    for i = 0:(size_glcm_1-1)
        out.denth(k) = out.denth(k) - (p_xminusy(i+1,k)*log(p_xminusy(i+1,k) + eps));
        out.dvarh(k) = out.dvarh(k) + (i^2)*p_xminusy(i+1,k);
    end
end

% compute information measure of correlation(1,2) [1]
for k = 1:size_glcm_3
    hxy(k) = out.entro(k);
    for i = 1:size_glcm_1
        
        for j = 1:size_glcm_2
            hxy1(k) = hxy1(k) - (glcm(i,j,k)*log(p_x(i,k)*p_y(j,k) + eps));
            hxy2(k) = hxy2(k) - (p_x(i,k)*p_y(j,k)*log(p_x(i,k)*p_y(j,k) + eps));
%             for Qind = 1:(size_glcm_1)
%                 Q(i,j,k) = Q(i,j,k) +...
%                     ( glcm(i,Qind,k)*glcm(j,Qind,k) / (p_x(i,k)*p_y(Qind,k)) ); 
%             end
        end
        hx(k) = hx(k) - (p_x(i,k)*log(p_x(i,k) + eps));
        hy(k) = hy(k) - (p_y(i,k)*log(p_y(i,k) + eps));
    end
    out.inf1h(k) = ( hxy(k) - hxy1(k) ) / ( max([hx(k),hy(k)]) );
    out.inf2h(k) = ( 1 - exp( -2*( hxy2(k) - hxy(k) ) ) )^0.5;
%     eig_Q(k,:)   = eig(Q(:,:,k));
%     sort_eig(k,:)= sort(eig_Q(k,:),'descend');
%     out.mxcch(k) = sort_eig(k,2)^0.5;
% The maximal correlation coefficient was not calculated due to
% computational instability 
% http://murphylab.web.cmu.edu/publications/boland/boland_node26.html
end

corm = zeros(size_glcm_3,1);
corp = zeros(size_glcm_3,1);
% using http://www.fp.ucalgary.ca/mhallbey/glcm_variance.htm for s_x s_y
for k = 1:size_glcm_3
    for i = 1:size_glcm_1
        for j = 1:size_glcm_2
            s_x(k)  = s_x(k)  + (((i) - u_x(k))^2)*glcm(i,j,k);
            s_y(k)  = s_y(k)  + (((j) - u_y(k))^2)*glcm(i,j,k);
            corp(k) = corp(k) + ((i)*(j)*glcm(i,j,k));
            corm(k) = corm(k) + (((i) - u_x(k))*((j) - u_y(k))*glcm(i,j,k));
            out.cprom(k) = out.cprom(k) + (((i + j - u_x(k) - u_y(k))^4)*...
                glcm(i,j,k));
            out.cshad(k) = out.cshad(k) + (((i + j - u_x(k) - u_y(k))^3)*...
                glcm(i,j,k));
        end
    end
    % using http://www.fp.ucalgary.ca/mhallbey/glcm_variance.htm for s_x
    % s_y : This solves the difference in value of correlation and might be
    % the right value of standard deviations required 
    % According to this website there is a typo in [2] which provides
    % values of variance instead of the standard deviation hence a square
    % root is required as done below:
    s_x(k) = s_x(k) ^ 0.5;
    s_y(k) = s_y(k) ^ 0.5;
    out.autoc(k) = corp(k);
    out.corrp(k) = (corp(k) - u_x(k)*u_y(k))/(s_x(k)*s_y(k));
    out.corrm(k) = corm(k) / (s_x(k)*s_y(k));
%     % alternate values of u and s
%     out.corrp2(k) = (corp(k) - u_x2(k)*u_y2(k))/(s_x2(k)*s_y2(k));
%     out.corrm2(k) = corm(k) / (s_x2(k)*s_y2(k));
end
% Here the formula in the paper out.corrp and the formula in matlab
% out.corrm are equivalent as confirmed by the similar results obtained

% % The papers have a slightly different formular for Contrast
% % I have tested here to find this formula in the papers provides the 
% % same results as the formula provided by the matlab function for 
% % Contrast (Hence this part has been commented)
% out.contrp = zeros(size_glcm_3,1);
% contp = 0;
% Ng = size_glcm_1;
% for k = 1:size_glcm_3
%     for n = 0:(Ng-1)
%         for i = 1:Ng
%             for j = 1:Ng
%                 if (abs(i-j) == n)
%                     contp = contp + glcm(i,j,k);
%                 end
%             end
%         end
%         out.contrp(k) = out.contrp(k) + n^2*contp;
%         contp = 0;
%     end
%     
% end

%       GLCM Features (Soh, 1999; Haralick, 1973; Clausi 2002)
%           f1. Uniformity / Energy / Angular Second Moment (done)
%           f2. Entropy (done)
%           f3. Dissimilarity (done)
%           f4. Contrast / Inertia (done)
%           f5. Inverse difference    
%           f6. correlation
%           f7. Homogeneity / Inverse difference moment
%           f8. Autocorrelation
%           f9. Cluster Shade
%          f10. Cluster Prominence
%          f11. Maximum probability
%          f12. Sum of Squares
%          f13. Sum Average
%          f14. Sum Variance
%          f15. Sum Entropy
%          f16. Difference variance
%          f17. Difference entropy
%          f18. Information measures of correlation (1)
%          f19. Information measures of correlation (2)
%          f20. Maximal correlation coefficient
%          f21. Inverse difference normalized (INN)
%          f22. Inverse difference moment normalized (IDN)

SVM Classifier

https://github.com/Abusamra/LIBSVM-multi-classification

Matlab Plot Functions(Matlab 画图函数）

All kinds of plot functions can be easily found in the reference link.

Good luck!

Reference：
https://cn.mathworks.com/help/matlab/creating_plots/types-of-matlab-plots.html

关于灰度共生矩阵的理解

首先来介绍灰度共生矩阵的公式，这是最重要也是最基本的。

Fig.1 GLCM 计算公式

Fig. 2 GLCM计算参数实验

其中，W为滑动窗口的大小，Moving step为滑动窗口移动的距离，而d为共生距离，也就是灰度共生矩阵中的d，最后的K为量化后的灰度级别。

Fig.3 滑动窗口的理解

Note: Fig.3 所示为水平方向（0度）获得共生矩阵的演示图。其他3个方向（45度、90度和135度）情况类似。

每一个窗口均获得一个GLCM，然后依据这些GLCM可以获得特征参数，然后即可获得如Fig.4 的散点关系图。

Fig.4 不同纹理特征的散点图

注意

需要注意的一点是，如果选定4个方向为GLCM的计算方向，那么在根据GLCM计算纹理特征时，需要首先根据4个方向的GLCM计算其纹理特征，然而将这4个方向的纹理特征取平均，获得最终的纹理特征计算结果[1]。

Reference:

[1] https://cn.mathworks.com/help/images/specify-offset-used-in-glcm-calculation.html

灰度共生矩阵

（先来补充）ENVI中关于灰度共生矩阵（又称为二阶矩阵，Second-Order Metrics）的介绍（另包含一阶矩阵的介绍）：https://harrisgeospatial.com/docs/backgroundtexturemetrics.html

灰度共生矩阵（Gray Level Co-occurrence Matrix）是影像像元间角度以及距离二者的函数。该矩阵中的数值反应的是影像中某个像元与某一个方向与距离处指定像元像元值出现的次数。该算法可以很好的进行对影像进行纹理分析，来源是Haralick, Shanmugan, and Dinstein (1973)的文章。

灰度共生矩阵的运算举例说明：

图1 数组实例

图1数组的像元数值范围为0-6，ENVI在计算的过程中，假定定义一个3×3的目标窗口，在上面的这步移动过程(X 方向右移1个像素，Y方向不移动)中，仅仅对于这个目标区域以及这一次的移动过程，可以获得以下共生矩阵：

图2 灰度共生矩阵

以图2中（6，4）处的数值为2进行解释，即（x，y）相对于（x+1，y）像元为（6，4）出现的次数为2次，所以是2。同理，可以理解图2中其他数值的含义。
如果该3×3窗口继续移动，遍历整个图像，则灰度共生矩阵不单单是4×4的数组，应该是5×5的数组，这是因为原示例数组的数值范围为2-6，同时，灰度共生矩阵中的数值在窗口移动的过程中不断被改写，由于需要遍历整幅影像，因而该过程比较耗时。

参考网页链接：

https://harrisgeospatial.com/docs/backgroundtexturemetrics.html

Reference(参考文献)：
Haralick, R., Shanmugan, K., and Dinstein, I. "Textural Features for Image Classification." IEEE Transactions on Systems, Man, and Cybernetics 3, no. 6 (1973): 610-621.

Import Snappy Module in Python（在Python中导入Snappy模块）

Reference：

https://www.math.uic.edu/t3m/SnapPy/installing.html