高光谱图像分类
任务目标要求
阅读论文《HybridSN: Exploring 3-D–2-DCNN Feature Hierarchy for Hyperspectral Image Classification》，进行相关代码实验，完成对高光谱图像的分类，进而思考以下问题
思考
1、3D卷积和2D卷积的区别？
2、对于训练好的同一个模型，为什么多次测试的结果不同？
3、如果想要进一步提升高光谱图像的分类性能，可以如何使用注意力机制？
代码实现部分
取得数据，并引入基本函数库。
! 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
! pip install spectral
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
定义 HybridSN 类
模型网络结构图如下
三维卷积部分：
conv1：（1, 30, 25, 25）， 8个 7x3x3 的卷积核 ==>（8, 24, 23, 23）
conv2：（8, 24, 23, 23）， 16个 5x3x3 的卷积核 ==>（16, 20, 21, 21）
conv3：（16, 20, 21, 21），32个 3x3x3 的卷积核 ==>（32, 18, 19, 19）
接下来要进行二维卷积，因此把前面的 32*18 reshape 一下，得到 （576, 19, 19）
二维卷积：（576, 19, 19） 64个 3x3 的卷积核，得到 （64, 17, 17）
接下来是一个 flatten 操作，变为 18496 维的向量，
接下来依次为256，128节点的全连接层，都使用比例为0.4的 Dropout，
最后输出为 16 个节点，是最终的分类类别数。
下面是 HybridSN 类的代码：
class HybridSN(nn.Module):
def __init__(self):
super(HybridSN, self).__init__()
self.conv3d_1 = nn.Sequential(
nn.Conv3d(1, 8, kernel_size=(7, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(8),
nn.ReLU(inplace = True),
)
self.conv3d_2 = nn.Sequential(
nn.Conv3d(8, 16, kernel_size=(5, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(16),
nn.ReLU(inplace = True),
)
self.conv3d_3 = nn.Sequential(
nn.Conv3d(16, 32, kernel_size=(3, 3, 3), stride=1, padding=0),
nn.BatchNorm3d(32),
nn.ReLU(inplace = True)
)
self.conv2d_4 = nn.Sequential(
nn.Conv2d(576, 64, kernel_size=(3, 3), stride=1, padding=0),
nn.BatchNorm2d(64),
nn.ReLU(inplace = True),
)
self.fc1 = nn.Linear(18496,256)
self.fc2 = nn.Linear(256,128)
self.fc3 = nn.Linear(128,16)
self.dropout = nn.Dropout(p = 0.4)
def forward(self,x):
out = self.conv3d_1(x)
out = self.conv3d_2(out)
out = self.conv3d_3(out)
out = self.conv2d_4(out.reshape(out.shape[0],-1,19,19))
out = out.reshape(out.shape[0],-1)
out = F.relu(self.dropout(self.fc1(out)))
out = F.relu(self.dropout(self.fc2(out)))
out = self.fc3(out)
return out
创建数据集
# 对高光谱数据 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
# 对单个像素周围提取 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
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
# 地物类别
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
print('Hyperspectral data shape: ', X.shape)
print('Label shape: ', y.shape)
print('\n... ... PCA tranformation ... ...')
X_pca = applyPCA(X, numComponents=pca_components)
print('Data shape after PCA: ', X_pca.shape)
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)
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)
# 改变 Xtrain, Ytrain 的形状，以符合 keras 的要求
Xtrain = Xtrain.reshape(-1, patch_size, patch_size, pca_components, 1)
Xtest = Xtest.reshape(-1, patch_size, patch_size, pca_components, 1)
print('before transpose: Xtrain shape: ', Xtrain.shape)
print('before transpose: Xtest shape: ', Xtest.shape)
# 为了适应 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
# 创建 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)
开始训练
# 使用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')
测试
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)
显示分类效果
# load the original image
X = sio.loadmat('Indian_pines_corrected.mat')['indian_pines_corrected']
y = sio.loadmat('Indian_pines_gt.mat')['indian_pines_gt']
height = y.shape[0]
width = y.shape[1]
X = applyPCA(X, numComponents= pca_components)
X = padWithZeros(X, patch_size//2)
# 逐像素预测类别
outputs = np.zeros((height,width))
for i in range(height):
for j in range(width):
if int(y[i,j]) == 0:
continue
else :
image_patch = X[i:i+patch_size, j:j+patch_size, :]
image_patch = image_patch.reshape(1,image_patch.shape[0],image_patch.shape[1], image_patch.shape[2], 1)
X_test_image = torch.FloatTensor(image_patch.transpose(0, 4, 3, 1, 2)).to(device)
prediction = net(X_test_image)
prediction = np.argmax(prediction.detach().cpu().numpy(), axis=1)
outputs[i][j] = prediction+1
if i % 20 == 0:
print('... ... row ', i, ' handling ... ...')
predict_image = spectral.imshow(classes = outputs.astype(int),figsize =(5,5))
模型测试准确度很高，但是却得不到理想的分类结果，下文尝试使用注意力机制解决问题
回答思考问题
1、3D卷积和2D卷积的区别？
2D卷积的卷积核有两个自由度，在本文中是对图像每个通道的数据面进行操作，2个维度分别是水平与垂直方向的滑动，代表空间维度。
3D卷积的卷积核有三个自由度，对于不同的问题空间，3D卷积代表不同的含义，本文中的3D卷积是在空间维度的基础上加入了一个光谱波段维度，引入光谱波段的维度来表示同一物质（如山、水、房屋）在不同光谱中的表现，从而完成分类工作。3D卷积对数据块进行操作，得到的仍是一个数据块。
正如论文作者所提到的，2D卷积和3D卷积在高光谱图像分类工作中各有优劣：
仅2-D-CNN不能从光谱维度中提取良好的区分特征图。 相似地，较深的3-D-CNN在计算上更加复杂，这也是坐着提出混合CNN模型的动机。 3-D-CNN和2-D-CNN层组装在一起，使得它们充分利用光谱图和空间特征图，以实现最大可能的精度，得到理想的分类结果。
2、对于训练好的同一个模型，为什么多次测试的结果不同？
网络的全连接层中使用了 nn.Dropout(p=0.4)，网络层的节点会随机失活，所以每次测试结果不同
3、如果想要进一步提升高光谱图像的分类性能，可以如何使用注意力机制？
这里加入空间注意力机制测试高光谱图像的分类结果
# 空间注意力机制
class SpatialAttention(nn.Module):
def __init__(self, kernel_size=7):
super(SpatialAttention, self).__init__()
assert kernel_size in (3, 7), 'kernel size must be 3 or 7'
padding = 3 if kernel_size == 7 else 1
self.conv1 = nn.Conv2d(2, 1, kernel_size, padding=padding, bias=False)
self.sigmoid = nn.Sigmoid()
def forward(self, x):
avg_out = torch.mean(x, dim=1, keepdim=True)
max_out, _ = torch.max(x, dim=1, keepdim=True)
x = torch.cat([avg_out, max_out], dim=1)
x = self.conv1(x)
return self.sigmoid(x)
# 网络骨架
class HybridSN(nn.Module):
def __init__(self, num_classes, self_attention=False):
super(HybridSN, self).__init__()
# out = (width - kernel_size + 2*padding)/stride + 1
# => padding = ( stride * (out-1) + kernel_size - width)
# 这里因为 stride == 1 所有卷积计算得到的padding都为 0
#默认不使用注意力机制
self.self_attention = self_attention
# 3D卷积块
self.block_1_3D = nn.Sequential(
nn.Conv3d(
in_channels=1,
out_channels=8,
kernel_size=(7, 3, 3),
stride=1,
padding=0
),
nn.ReLU(inplace=True),
nn.Conv3d(
in_channels=8,
out_channels=16,
kernel_size=(5, 3, 3),
stride=1,
padding=0
),
nn.ReLU(inplace=True),
nn.Conv3d(
in_channels=16,
out_channels=32,
kernel_size=(3, 3, 3),
stride=1,
padding=0
),
nn.ReLU(inplace=True)
)
if self_attention:
self.spatial_attention_1 = SpatialAttention(kernel_size=7)
# 2D卷积块
self.block_2_2D = nn.Sequential(
nn.Conv2d(
in_channels=576,
out_channels=64,
kernel_size=(3, 3)
),
nn.ReLU(inplace=True)
)
if self_attention:
self.spatial_attention_2 = SpatialAttention(kernel_size=7)
# 全连接层
self.classifier = nn.Sequential(
nn.Linear(
in_features=18496,
out_features=256
),
nn.Dropout(p=0.4),
nn.Linear(
in_features=256,
out_features=128
),
nn.Dropout(p=0.4),
nn.Linear(
in_features=128,
out_features=num_classes
)
# pytorch交叉熵损失函数是混合了softmax的。不需要再使用softmax
)
def forward(self, x):
y = self.block_1_3D(x)
y = y.view(-1, y.shape[1] * y.shape[2], y.shape[3], y.shape[4])
if self.self_attention:
# y = self.channel_attention_1(y) * y
y = self.spatial_attention_1(y) * y
y = self.block_2_2D(y)
if self.self_attention:
# y = self.channel_attention_2(y) * y
y = self.spatial_attention_2(y) * y
y = y.view(y.size(0), -1)
y = self.classifier(y)
return y
仅使用空间注意力机制效果有了很大的提升，推测进一步还可以在通道、光谱波段使用注意力机制来提高网络的分类能力。
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
