第四周 HybridSN 高光谱分类

一、论文阅读

参考论文：HybridSN: Exploring 3D-2D CNN FeatureHierarchy for Hyperspectral Image Classification
论文pdf链接：https://arxiv.org/pdf/1902.06701.pdf

二、HybridSN网络结构的搭建

1、数据集下载

有两种方法：
第一种方法使用wget方法
! wget http://www.ehu.eus/ccwintco/uploads/6/67/Indian_pines_corrected.mat
! wget http://www.ehu.eus/ccwintco/uploads/c/c4/Indian_pines_gt.mat
第二种方法下载到本地并上传到google drive进行colab与drive的链接
from google.colab import drive
drive.mount('/content/drive')

class_num=16
X=sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']
y=sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']
#用于测试样本的比例
test_ratio=0.90
#每个像素周围提取patch的尺寸
patch_size=25
#使用PCA降维，得到主成分的数量
pca_components=30

2、导入所需要的库

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
from sklearn.decomposition import PCA
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix, accuracy_score, classification_report, cohen_kappa_score
import spectral
import torch
import torchvision
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim

3、搭建HybridSN网络结构

分为三维卷积、二维卷积以及FC部分。
1、三维卷积：四个参数为（Channel,D,M,N）
（1）conv1:(1,30,25,25)  Kernel_size:(7,3,3)  output:(8,24,23,23)
（2）conv2:(8,24,23,23)  Kernel_size:(5,3,3)  output:(16,20,21,21)
（3）conv3:(16,20,21,21) Kernel_size:(3,3,3)  output:(32,18,19,19)
2、二维卷积：进行Conv2d之前需要先把数据reshape为（576，19，19）
	conv4：（576, 19, 19） Kernel_size(3,3) output：(64,17,17)
3、FC：需要变成18496维向量。经过三层FC层，输出16个种类。

class HybridSN(nn.Module):
  def __init__(self):
    super(HybridSN,self).__init__()
    self.conv1=nn.Conv3d(in_channels=1,out_channels=8,kernel_size=(7,3,3))
    self.bn1=nn.BatchNorm3d(8)
    self.relu1=nn.ReLU(inplace=True)
    

    self.conv2=nn.Conv3d(in_channels=8,out_channels=16,kernel_size=(5,3,3))
    self.bn2=nn.BatchNorm3d(16)
    self.relu2=nn.ReLU(inplace=True)

    self.conv3=nn.Conv3d(in_channels=16,out_channels=32,kernel_size=(3,3,3))
    self.bn3=nn.BatchNorm3d(32)
    self.relu3=nn.ReLU(inplace=True)

    self.conv4=nn.Conv2d(in_channels=576,out_channels=64,kernel_size=(3,3),stride=1)
    self.bn4=nn.BatchNorm2d(64)
    self.relu4=nn.ReLU(inplace=True)

    self.fc1=nn.Linear(in_features=18496,out_features=256)
    self.fc2=nn.Linear(in_features=256,out_features=128)
    self.fc3=nn.Linear(in_features=128,out_features=class_num)
    self.dropout1=nn.Dropout(p=0.4)
    self.dropout2=nn.Dropout(p=0.4)

  def forward(self,x):
    out=self.conv1(x)
    out=self.relu1(out)
    out=self.bn1(out)

    out=self.conv2(out)
    out=self.bn2(out)
    out=self.relu2(out)

    out=self.conv3(out)
    out=self.bn3(out)
    out=self.relu3(out)

    out=out.reshape((out.shape[0],-1,19,19))
    out=self.conv4(out)
    out=self.bn4(out)
    out=self.relu4(out)
    
    out = out.reshape(out.shape[0],-1)
    #print(out.shape)
    out=self.fc1(out)
    out=self.dropout1(out)
    out=self.fc2(out)
    out=self.dropout2(out)
    out=self.fc3(out)

    return out

随机输入，测试网络结构是否连通

x = torch.randn(1, 1, 30, 25, 25)
net = HybridSN()
y = net(x)
print(y.shape)

4、定义需要用到的函数

1、PCA函数
去除光谱冗余，M×N×D ===> M×N×B

# 对高光谱数据 X 应用 PCA 变换
def applyPCA(X, numComponents):
    newX = np.reshape(X, (-1, X.shape[2]))
    pca = PCA(n_components=numComponents, whiten=True)
    newX = pca.fit_transform(newX)
    newX = np.reshape(newX, (X.shape[0], X.shape[1], numComponents))
    return newX

print('\n... ... PCA tranformation ... ...')
X_pca = applyPCA(X, numComponents=pca_components)
print('Data shape after PCA: ', X_pca.shape)

2、提取patch

# 对单个像素周围提取 patch 时，边缘像素就无法取了，因此，给这部分像素进行 padding 操作
def padWithZeros(X, margin=2):
    newX = np.zeros((X.shape[0] + 2 * margin, X.shape[1] + 2* margin, X.shape[2]))
    x_offset = margin
    y_offset = margin
    newX[x_offset:X.shape[0] + x_offset, y_offset:X.shape[1] + y_offset, :] = X
    return newX

# 在每个像素周围提取 patch ，然后创建成符合 keras 处理的格式
def createImageCubes(X, y, windowSize=5, removeZeroLabels = True):
    # 给 X 做 padding
    margin = int((windowSize - 1) / 2)
    zeroPaddedX = padWithZeros(X, margin=margin)
    # split patches
    patchesData = np.zeros((X.shape[0] * X.shape[1], windowSize, windowSize, X.shape[2]))
    patchesLabels = np.zeros((X.shape[0] * X.shape[1]))
    patchIndex = 0
    for r in range(margin, zeroPaddedX.shape[0] - margin):
        for c in range(margin, zeroPaddedX.shape[1] - margin):
            patch = zeroPaddedX[r - margin:r + margin + 1, c - margin:c + margin + 1]   
            patchesData[patchIndex, :, :, :] = patch
            patchesLabels[patchIndex] = y[r-margin, c-margin]
            patchIndex = patchIndex + 1
    if removeZeroLabels:
        patchesData = patchesData[patchesLabels>0,:,:,:]
        patchesLabels = patchesLabels[patchesLabels>0]
        patchesLabels -= 1
    return patchesData, patchesLabels

print('\n... ... create data cubes ... ...')
X_pca, y = createImageCubes(X_pca, y, windowSize=patch_size)
print('Data cube X shape: ', X_pca.shape)
print('Data cube y shape: ', y.shape)

3、划分训练集以及测试集

def splitTrainTestSet(X, y, testRatio, randomState=345):
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=testRatio, random_state=randomState, stratify=y)
    return X_train, X_test, y_train, y_test

print('\n... ... create train & test data ... ...')
Xtrain, Xtest, ytrain, ytest = splitTrainTestSet(X_pca, y, test_ratio)
print('Xtrain shape: ', Xtrain.shape)
print('Xtest  shape: ', Xtest.shape)

4、数据处理
维度转置，变成（channel，D，W，H）格式；Tensor处理

# 为了适应 pytorch 结构，数据要做 transpose
Xtrain = Xtrain.transpose(0, 4, 3, 1, 2)
Xtest  = Xtest.transpose(0, 4, 3, 1, 2)
print('after transpose: Xtrain shape: ', Xtrain.shape) 
print('after transpose: Xtest  shape: ', Xtest.shape) 

""" Training dataset"""
class TrainDS(torch.utils.data.Dataset): 
    def __init__(self):
        self.len = Xtrain.shape[0]
        self.x_data = torch.FloatTensor(Xtrain)
        self.y_data = torch.LongTensor(ytrain)        
    def __getitem__(self, index):
        # 根据索引返回数据和对应的标签
        return self.x_data[index], self.y_data[index]
    def __len__(self): 
        # 返回文件数据的数目
        return self.len

""" Testing dataset"""
class TestDS(torch.utils.data.Dataset): 
    def __init__(self):
        self.len = Xtest.shape[0]
        self.x_data = torch.FloatTensor(Xtest)
        self.y_data = torch.LongTensor(ytest)
    def __getitem__(self, index):
        # 根据索引返回数据和对应的标签
        return self.x_data[index], self.y_data[index]
    def __len__(self): 
        # 返回文件数据的数目
        return self.len

5、加载数据	

# 创建 trainloader 和 testloader
trainset = TrainDS()
testset  = TestDS()
train_loader = torch.utils.data.DataLoader(dataset=trainset, batch_size=128, shuffle=True, num_workers=2)
test_loader  = torch.utils.data.DataLoader(dataset=testset,  batch_size=128, shuffle=False, num_workers=2)

5、训练模型

# 使用GPU训练，可以在菜单 "代码执行工具" -> "更改运行时类型" 里进行设置
device = torch.device("cuda:0" if torch.cuda.is_available() else "cpu")

# 网络放到GPU上
net = HybridSN().to(device)
criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(net.parameters(), lr=0.001)

# 开始训练
total_loss = 0
for epoch in range(100):
    for i, (inputs, labels) in enumerate(train_loader):
        inputs = inputs.to(device)
        labels = labels.to(device)
        # 优化器梯度归零
        optimizer.zero_grad()
        # 正向传播 +　反向传播 + 优化 
        outputs = net(inputs)
        loss = criterion(outputs, labels)
        loss.backward()
        optimizer.step()
        total_loss += loss.item()
    print('[Epoch: %d]   [loss avg: %.4f]   [current loss: %.4f]' %(epoch + 1, total_loss/(epoch+1), loss.item()))

print('Finished Training')

6、测试模型

count = 0
# 模型测试
for inputs, _ in test_loader:
    inputs = inputs.to(device)
    outputs = net(inputs)
    outputs = np.argmax(outputs.detach().cpu().numpy(), axis=1)
    if count == 0:
        y_pred_test =  outputs
        count = 1
    else:
        y_pred_test = np.concatenate( (y_pred_test, outputs) )

# 生成分类报告
classification = classification_report(ytest, y_pred_test, digits=4)
print(classification)

准确率为97.18%

三、问题总结

1、训练HybridSN，然后多测试几次，会发现每次分类的结果都不一样，请思考为什么？
可能由于不同物体的相似性，相同物体的互异性。
2、如果想要进一步提升高光谱图像的分类性能，可以如何改进？
引入注意力机制；数据增强；
3、depth-wise conv 和 分组卷积有什么区别与联系？
区别：DW卷积，每个卷积核的Channel=1；分组卷积（加入分为3组），每个卷积核的Channel变为原来的$\frac{1}{3}$。
联系：DW卷积可以看成一组只有一个Channel，是特殊的分组卷积。
4、SENet 的注意力是不是可以加在空间位置上？
可以加在空间位置上。
5、在 ShuffleNet 中，通道的 shuffle 如何用代码实现？

def channel_shuffle(x:Tensor,groups:int):
	batch_size,num_channels,height,width=x.size()
	channels_per_group=num_channels//groups
	x=x.view(batch_size,groups,channels_per_group,height,width)
	#把groups维度与channels_per_group维度进行调换。可以理解为按照channel进行了分类。
	x=troch.transpose(x,1,2).contiguous()
	x=x.view(batch_size,-1,height,width)

来源：https://www.bilibili.com/video/BV1dh411r76X/?spm_id_from=333.788&vd_source=dfc6d3291de08a14fa590df61c75fb19