跳到主要内容

火山下的文明——维苏威古卷

在上一讲中，我们探讨了自定义二值化的各种策略，这些策略有助于我们在默认二值化效果不理想的情况下找到解决方案。

今天，我们将以此为基础，挑战一个真实且持续进行的技术难题：《火山下的文明——维苏威古卷》

我希望通过这个项目，让大家了解到机器学习的完整流程，并对其中的关键概念给出代码并做出解释。

项目概述​

公元79年10月24日，意大利的维苏威火山爆发，一天之内就毁灭了两万多人的庞贝古城。

火山灰掩盖了整座城市，其中有一栋房子存放了各种书籍。直到18世纪，这栋房子才重新被发现，下面是考古学家的建筑复原图。

1711416683059

房子里面的1800多卷纸莎草古书，都已经完全碳化。

1711416773698

1711416781220

考古学家试图展开这些烧焦的古书，但是纸张一碰就碎。

1711416794510

每一页的碎片就像下面这样。

1711416808202

没有人能从这些烧焦的古书上，读出一个字。

2019年，一位美国科学家提出了一种新方法，使用短红外高光谱成像技术，对这些古卷进行 3D 扫描，然后虚拟建模，在不接触古卷的情况下，就可以展开或压平每一页，从而复现上面的笔迹。

1711416835510

1711416847661

下面是扫描后的页面结果。

1711416871566

问题是上面的笔迹还是很模糊，无法确定是哪一个字母，更不要说读出句子了。

2023年3月，在一些风险投资者的资助下，古卷保管单位举办了一次比赛，邀请全世界的程序员参加，要求使用机器学习，认出上面是什么字母，破解古卷的内容，奖金总额高达100万美元。

提示

该项目首要赞助商为马斯克基金会。

发掘地的图书馆已被证实属于当时罗马最有权势的人之一：凯撒的岳父。

红外成像部分使用的是蔡司的镜头（光刻机用的也是蔡司的镜头）

主办单位原以为，一年之内成功的可能性不到30%。但是，2023年10月，21岁的内布拉斯加大学的学生卢克·法里托（Luke Farritor）就读出了第一个单词 ΠΟΡΦΥΡΑϹ（紫色）。

1711418146719

后来，他与正在柏林读博士的纳德（Youssef Nader）、瑞士苏黎世理工学院的机器人专业的学生席里格（Julian Schilliger），组成了一个团队，致力于建立一个完整的 AI 模型识别这些古书。

他们最终在比赛的截止日期前（2024年1月1日）识别出了2000个字符。

下面就是采用他们的模型，结合红外成像识别出来的碎片，内容已经清晰还原出来了。

1711416997765

2月5日，主办单位宣布，他们获得了本次比赛的第一名。

1711417019492

由于比赛结果令人鼓舞，主办单位扩展了自己的目标，接下来将对90%的维苏威古卷，进行扫描和识别，彻底破解两千年前的罗马人在书里写了什么。

他们所用的模型，已经全部开源，就放在 GitHub。任何人都可以安装和运行，尝试改进他们的模型，获得更好的结果。

现代科技的发展，真像神话一样，烧成灰的纸都能辨认出文字。

项目官网： https://scrollprize.org/

该项目在B站有纪录片。

项目分析​

在2024年，文本识别已经是一项非常成熟的技术。然而，每个项目都有其独特性，阅读碳化卷轴的过程分为三个步骤：

扫描​

使用 X 射线断层扫描创建卷轴或片段的 3D 扫描

3D扫描的原理是：墨水区域与纸张区域对粒子的穿透率不同，因此在 X 射线扫描中墨水会呈现白色，纸张会呈现黑色。

提示

这个步骤通常由专业的机构完成，我们只需要获取扫描结果即可。 不同的粒子可以获得不同的建模，但是目前来说这一步并不会影响结果。

分割和展平​

在 3D 扫描中找到卷起的纸莎草纸的各层，然后将它们展开成展平的“表面体积”

卷轴。 从技术上讲，我们可以为整个卷轴制作一个巨大的“片段”，但实际上卷轴包装可能很难区分，因此我们将其分成更易于管理的部分。

分割可能具有挑战性，因为纸莎草的不同层可能会被损坏、扭曲或磨损。碳化的纸莎草会起泡，不同的层甚至可以相互融合。

1711419035325

在碎片上，这个过程要容易一些，因为它们已经相当平坦，并且有一个暴露的表面，我们实际上可以在上面看到墨水。尽管如此，这些碎片通常并不完全平坦，并且可能在可见层下方附着有纸莎草的“隐藏层”。

1711416794510

提示

新手可以从碎片入手，因为它们更容易处理。

墨水检测​

使用机器学习模型识别展平表面体积中的墨水区域

演示视频：https://scrollprize.org/img/tutorials/ink-detection-anim2-dark.webm

数据集：https://www.kaggle.com/competitions/vesuvius-challenge-ink-detection/data

test数据集包含 tif （一种切片图片格式）每个片段包含 65 个切片。将这些图像叠加在一起，每个片段的体素数量为宽＊高＊65。

提示

也就是说：00.tif -- 64.tif 是一个片段,表示一个炭化碎片。通过这个数据集，我们可以训练一个模型，识别出墨水区域。

想象一块豆腐被片成65片

  • mask.png — 二进制蒙版，表示哪些像素包含数据。

其中，train 数据集额外包含

  • ir.png - 人工拍摄的红外成像。 ::: 用于参考的图片 :::

  • inklabels.png — 一个区分墨水与无墨水标签的二值掩膜（二值化的图像，用于表示图像中的某个区域是否有墨水）。

提示

我们需要的结果图像

  • inklabels_rle.csv — 一个包含 inklabels 的 run-length 编码的 CSV 文件， 返回的数据是成对的： 连续的1开始的位置 ，长度。我们需要提交这样的结果给到主办方。这里主办方要求从1开始，不是从0开始。
提示

[0,1,1,0,0,1,1,1] 返回 2 2 表示从2个数开始，连续2个1 6 3 表示从6个数开始，连续3个1

# 导入 NumPy 库，用于进行数值计算
import numpy as np

# 导入 Pillow 库，用于处理图像
from PIL import Image

# 定义一个函数，名为 rle，该函数用于实现图像的行程长度编码
def rle (img):
    # 将输入的图像（假设为二维数组）展平为一维数组
    flat_img = img.flatten()
    
    # 对展平后的图像进行二值化处理，所有大于0.5的元素转换为1，其他元素转换为0
    # 数组的比较，与我们上节课内容一致
    flat_img = np.where(flat_img > 0.5, 1, 0).astype(np.uint8)

    # 找出所有从0变为1的位置，这些位置表示一段新的连续1的开始
    # 数组的逻辑运算
    starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))
'''
    [0,1,1,0,0,1,1,1,0] -> [T,F,F,T,T,F,F,F,*]
    [0,1,1,0,0,1,1,1,0] -> [*,T,T,F,F,T,T,T,F]

    [T,F,F,T,T,F,F,F] 
    & 
    [T,T,F,F,T,T,T,F]

print(True & True) # True
print(True & False) # False
print(False & True) # False
print(False & False) # False

    [T,F,F,F,T,F,F,F]
    
'''
    # 找出所有从1变为0的位置，这些位置表示一段连续1的结束
    ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))
    
    # 找出所有开始和结束的位置的索引
    starts_ix = np.where(starts)[0] + 2
    ends_ix = np.where(ends)[0] + 2
'''
starts = np.array([1,0,0,0,1,0,0,0])

print(np.where(starts))
# (array([0, 4], dtype=int64),)

print(np.where(starts)[0])
# [0 4]

print(np.where(starts)[0] + 2)
# [2 6]
'''
    # 计算每一段连续1的长度
    # 数组的减法
    lengths = ends_ix - starts_ix

    # 返回每一段连续1的开始位置和长度
    return starts_ix, lengths

# 打开名为 'inklabels.png' 的图像，并将其转换为 NumPy 数组
inklabels = np.array(Image.open('inklabels.png'), dtype=np.uint8)

# 对图像进行行程长度编码
starts_ix, lengths = rle(inklabels)

# 将开始位置和长度配对，转换为字符串，然后用空格连接
inklabels_rle = " ".join(map(str, sum(zip(starts_ix, lengths), ())))

# 将结果写入 'inklabels_rle.csv' 文件中，文件的第一行是 "Id,Predicted"，第二行开始是行程长度编码的结果
print("Id,Predicted\n1," + inklabels_rle, file=open('inklabels_rle.csv', 'w'))

可以用下面的代码来测试

import numpy as np

# 导入 Pillow 库，用于处理图像
from PIL import Image

# 定义一个函数，名为 rle，该函数用于实现图像的行程长度编码
def rle ():
    # 将输入的图像（假设为二维数组）展平为一维数组
    flat_img = np.array([0,1,1,0,0,1,1,1,0])
    
    # 对展平后的图像进行二值化处理，所有大于0.5的元素转换为1，其他元素转换为0
    # 数组的比较，与我们上节课内容一致
    flat_img = np.where(flat_img > 0.5, 1, 0).astype(np.uint8)

    # 找出所有从0变为1的位置，这些位置表示一段新的连续1的开始
    # 数组的逻辑运算
    starts = np.array((flat_img[:-1] == 0) & (flat_img[1:] == 1))

    # 找出所有从1变为0的位置，这些位置表示一段连续1的结束
    ends = np.array((flat_img[:-1] == 1) & (flat_img[1:] == 0))
    print(ends)
    
    # 找出所有开始和结束的位置的索引
    starts_ix = np.where(starts)[0] + 2
    ends_ix = np.where(ends)[0] + 2

    # 计算每一段连续1的长度
    # 数组的减法
    lengths = ends_ix - starts_ix

    # 返回每一段连续1的开始位置和长度
    return starts_ix, lengths

# 对图像进行行程长度编码
starts_ix, lengths = rle()

print(starts_ix, lengths)
# [2 6] [2 3]

任务聚焦：墨水检测（高阶学生需要掌握）​

到这里，我们需要暂停一下。思考下怎么编写代码。或者说，代码应该包含哪些部分？

import torch
import torch.nn as nn
import torch.optim as optim
import numpy as np
import glob
import PIL.Image as Image
import torch.utils.data as data
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from tqdm import tqdm
from ipywidgets import interact, fixed

PREFIX = '/kaggle/input/vesuvius-challenge/train/1/'
BUFFER = 30  # Buffer size in x and y direction
Z_START = 27 # First slice in the z direction to use
Z_DIM = 10   # Number of slices in the z direction
TRAINING_STEPS = 30000
LEARNING_RATE = 0.03
BATCH_SIZE = 32
DEVICE = torch.device("cuda" if torch.cuda.is_available() else "cpu")

plt.imshow(Image.open(PREFIX+"ir.png"), cmap="gray")

mask = np.array(Image.open(PREFIX+"mask.png").convert('1'))
label = torch.from_numpy(np.array(Image.open(PREFIX+"inklabels.png"))).gt(0).float().to(DEVICE)
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_title("mask.png")
ax1.imshow(mask, cmap='gray')
ax2.set_title("inklabels.png")
ax2.imshow(label.cpu(), cmap='gray')
plt.show()


# Load the 3d x-ray scan, one slice at a time
images = [np.array(Image.open(filename), dtype=np.float32)/65535.0 for filename in tqdm(sorted(glob.glob(PREFIX+"surface_volume/*.tif"))[Z_START:Z_START+Z_DIM])]
image_stack = torch.stack([torch.from_numpy(image) for image in images], dim=0).to(DEVICE)

fig, axes = plt.subplots(1, len(images), figsize=(15, 3))
for image, ax in zip(images, axes):
  ax.imshow(np.array(Image.fromarray(image).resize((image.shape[1]//20, image.shape[0]//20)), dtype=np.float32), cmap='gray')
  ax.set_xticks([]); ax.set_yticks([])
fig.tight_layout()
plt.show()

rect = (1100, 3500, 700, 950)
fig, ax = plt.subplots()
ax.imshow(label.cpu())
patch = patches.Rectangle((rect[0], rect[1]), rect[2], rect[3], linewidth=2, edgecolor='r', facecolor='none')
ax.add_patch(patch)
plt.show()


class SubvolumeDataset(data.Dataset):
    def __init__(self, image_stack, label, pixels):
        self.image_stack = image_stack
        self.label = label
        self.pixels = pixels
    def __len__(self):
        return len(self.pixels)
    def __getitem__(self, index):
        y, x = self.pixels[index]
        subvolume = self.image_stack[:, y-BUFFER:y+BUFFER+1, x-BUFFER:x+BUFFER+1].view(1, Z_DIM, BUFFER*2+1, BUFFER*2+1)
        inklabel = self.label[y, x].view(1)
        return subvolume, inklabel

model = nn.Sequential(
    nn.Conv3d(1, 16, 3, 1, 1), nn.MaxPool3d(2, 2),
    nn.Conv3d(16, 32, 3, 1, 1), nn.MaxPool3d(2, 2),
    nn.Conv3d(32, 64, 3, 1, 1), nn.MaxPool3d(2, 2),
    nn.Flatten(start_dim=1),
    nn.LazyLinear(128), nn.ReLU(),
    nn.LazyLinear(1), nn.Sigmoid()
).to(DEVICE)

print("Generating pixel lists...")
# Split our dataset into train and val. The pixels inside the rect are the 
# val set, and the pixels outside the rect are the train set.
# Adapted from https://www.kaggle.com/code/jamesdavey/100x-faster-pixel-coordinate-generator-1s-runtime
# Create a Boolean array of the same shape as the bitmask, initially all True
not_border = np.zeros(mask.shape, dtype=bool)
not_border[BUFFER:mask.shape[0]-BUFFER, BUFFER:mask.shape[1]-BUFFER] = True
arr_mask = np.array(mask) * not_border
inside_rect = np.zeros(mask.shape, dtype=bool) * arr_mask
# Sets all indexes with inside_rect array to True
inside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = True
# Set the pixels within the inside_rect to False
outside_rect = np.ones(mask.shape, dtype=bool) * arr_mask
outside_rect[rect[1]:rect[1]+rect[3]+1, rect[0]:rect[0]+rect[2]+1] = False
pixels_inside_rect = np.argwhere(inside_rect)
pixels_outside_rect = np.argwhere(outside_rect)

print("Training...")
train_dataset = SubvolumeDataset(image_stack, label, pixels_outside_rect)
train_loader = data.DataLoader(train_dataset, batch_size=BATCH_SIZE, shuffle=True)
criterion = nn.BCELoss()
optimizer = optim.SGD(model.parameters(), lr=LEARNING_RATE)
scheduler = torch.optim.lr_scheduler.OneCycleLR(optimizer, max_lr=LEARNING_RATE, total_steps=TRAINING_STEPS)
model.train()
# running_loss = 0.0
for i, (subvolumes, inklabels) in tqdm(enumerate(train_loader), total=TRAINING_STEPS):
    if i >= TRAINING_STEPS:
        break
    optimizer.zero_grad()
    outputs = model(subvolumes.to(DEVICE))
    loss = criterion(outputs, inklabels.to(DEVICE))
    loss.backward()
    optimizer.step()
    scheduler.step()
#     running_loss += loss.item()
#     if i % 3000 == 3000-1:
#         print("Loss:", running_loss / 3000)
#         running_loss = 0.0

eval_dataset = SubvolumeDataset(image_stack, label, pixels_inside_rect)
eval_loader = data.DataLoader(eval_dataset, batch_size=BATCH_SIZE, shuffle=False)
output = torch.zeros_like(label).float()
model.eval()
with torch.no_grad():
    for i, (subvolumes, _) in enumerate(tqdm(eval_loader)):
        for j, value in enumerate(model(subvolumes.to(DEVICE))):
            output[tuple(pixels_inside_rect[i*BATCH_SIZE+j])] = value

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(output.cpu(), cmap='gray')
ax2.imshow(label.cpu(), cmap='gray')
plt.show()

THRESHOLD = 0.4
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(output.gt(THRESHOLD).cpu(), cmap='gray')
ax2.imshow(label.cpu(), cmap='gray')
plt.show()