OIF值的计算

from PIL import Image
import matplotlib.pyplot as plt
import scipy.stats as sts
import numpy as np
import cv2 as cv
from numpy import cov,corrcoef
from scipy.special import comb, perm
from itertools import combinations, permutations
img1 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B1.TIF",1))
img2 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B2.TIF",1))
img3 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B3.TIF",1))
img4 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B4.TIF",1))
img5 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B5.TIF",1))
img6 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B6.TIF",1))
img7 = np.array(cv.imread("LC08_L1TP_121034_20170626_20170713_01_T1_B7.TIF",1))
IMGN2P=[img1,img2,img3,img4,img5,img6,img7]   ### 二维矩阵数组

img1np=np.array(img1).flatten()   
img2np=np.array(img2).flatten() 
img3np=np.array(img3).flatten() 
img4np=np.array(img4).flatten() 
img5np=np.array(img5).flatten() 
img6np=np.array(img6).flatten() 
img7np=np.array(img7).flatten() 
IMGNP=[img1np,img2np,img3np,img4np,img5np,img6np,img7np]   ###  一维矩阵数组
### 创建 数组
mean=[0.00,1.00,2.00,3.00,4.00,5.00,6.00]   # 均值数组
ptp=range(7) #   极差
var=range(7) #  方差
std=range(7)# 标准差
mean=np.array(mean)          ###  均值矩阵
ptp=np.array(ptp,dtype=float)   ###   极差矩阵
var=np.array(var,dtype=float)   ###   方差矩阵
std=np.array(std,dtype=float)   ###   标准差
###cov corrcof   ###### 协方差矩阵    
data1 =  np.array([img1np,img2np]) ##波段1.2
cove = cov(data1,bias=1)
print("波段1 2  协方差:",cove)
print("波段1 2 相关系数: ",corrcoef(data1))

###1 2 3 波段OIF
Sk=std[0]+std[1]+std[2]
data2= np.array([img1np,img3np])  ## 波段 1  3
data3=np.array([img2np,img3np])  ##波段 2 3
R1=abs(corrcoef(data1).min())
R2=abs(corrcoef(data2).min())
R3=abs(corrcoef(data3).min())
R=R1+R2+R3
OIF=Sk/R
print("OIF(1 2 3 ) : \n",OIF)

# ncb73=comb(7,3)
nstd=list(range(100))
# ncb72=comb(7,2)
ncorrcoef=range(21)
k=0
npcrroef=np.zeros((7,7),dtype=float)
for x in range(6):   ###相关系数
    mxk=x
    for y in range(1,6-x):
        mxk=mxk+1
        ndata=np.array([IMGNP[x],IMGNP[mxk]])   ### ju zhen  qiu xiang guan xi shu 相关系数矩阵  12  13 ....
        npcrroef[x][mxk]=abs(corrcoef(ndata).min())
        k=k+1
m2=0
mk=0
for x in  range(5):
    m2=m2+1
    men=m2
    for y in range(5-m2):
        men=men+1
        nstd[mk]=std[x]+std[m2]+std[men]
for x in range(35):
    if x < 5:
        nco = npcrroef[0][x+2]+npcrroef[0][1]+npcrroef[1][x+2]
        OIF = nstd[x]/nco
        print("OIF( 1 2 ",x+3," )",OIF)
    elif x < 9:
        nco = npcrroef[0][2]+npcrroef[0][x-2]+npcrroef[2][x-2]
        OIF = nstd[x]/nco
        print("OIF(1 3 ",x-1," )",OIF)
    elif x < 12:
        nco = npcrroef[0][3]+npcrroef[0][x-5]+npcrroef[3][x-5]
        OIF = nstd[x]/nco
        print("OIF ( 1 4 ",x-4," )",OIF)
    elif x < 14:
        nco = npcrroef[0][4]+npcrroef[0][x-7]+npcrroef[4][x-7]
        OIF= nstd[x]/nco
        print("OIF ( 1 5 ",x-7," )",OIF)
    elif x <15:
        nco = npcrroef[0][5]+npcrroef[0][6]+npcrroef[5][6]
        OIF= nstd[x]/nco
        print("OIF ( 1 6 7)",OIF)
    elif x < 19:
        nco = npcrroef[1][2]+npcrroef[1][x-12]+npcrroef[2][x-12]
        OIF=nstd[x]/nco
        print("OIF ( 2 3 ",x-11," )",OIF)
    elif x < 22:
        nco=npcrroef[1][3]+npcrroef[1][x-15]+npcrroef[3][x-15]
        OIF=nstd[x]/nco
        print("OIF ( 2 4 ",x-14," )",OIF)
    elif x< 24:
        nco=npcrroef[1][4]+npcrroef[1][x-17]+npcrroef[4][x-17]
        OIF=nstd[x]/nco
        print("OIF (2 5 ",x-16," )",OIF)
    elif x <25:
        nco=npcrroef[1][5]+npcrroef[1][x-18]+npcrroef[5][x-18]
        OIF=nstd[x]/nco
        print("OIF (2 6 7)",OIF)
    elif x < 28:
        nco = npcrroef[2][3]+npcrroef[2][x-21]+npcrroef[3][x-21]
        OIF=nstd[x]/nco
        print("OIF (3 4 ",x-20,OIF)
    elif x < 30:
         nco=npcrroef[2][4]+npcrroef[2][x-23]+npcrroef[4][x-23]
         OIF=nstd[x]/nco
         print("OIF (3 5 ",x-22," )",OIF)
    elif x < 31:
        nco=npcrroef[2][5]+npcrroef[2][6]+npcrroef[5][6]
        OIF=nstd[x]/nco
        print("OIF （3 6 7）",OIF)
    elif x < 33:
        nco=npcrroef[3][4]+npcrroef[3][x-26]+npcrroef[4][x-26]
        OIF=nstd[x]/nco
        print("OIF (4 5 ",x-25," )",OIF)
    elif x<34:
        nco=npcrroef[3][5]+npcrroef[3][6]+npcrroef[5][6]
        OIF=nstd[x]/nco
        print("OIF (4 6 7)",OIF)
    elif x <35:
        nco=npcrroef[4][5]+npcrroef[4][6]+npcrroef[5][6]
        OIF=nstd[x]/nco
        print("OIF( 5 6 7)",OIF)