从断层扫描投影中提取正弦图

Question

我叫Lorenzo，是意大利的一名博士后研究员。我的工作与使用同步辐射的断层扫描成像有关。这对我来说是一个新的领域，我开始面对Matlab写一些代码。

我对层析成像完全陌生，Matlab 向我展示了新的挑战。我的实际问题是从一堆平行图像投影中创建正弦图。对于不在现场的人，正弦图是检测样本中特征位置在检测器上的投影的图，作为 X 射线束和样本之间角度的函数。

我从实验中得到的是一系列不同角度的 2D 射线照相，您可以将其视为一个矩形 "volume"，其中维度分别是单个投影中的行数和列数，以及体积由角度数给出。正弦图只是这个体积的横向切割。这意味着不是从侧面而是从顶部读取体积，所以我将创建一个新的图像数组，其维度是列数和投影，数组长度是投影中的行数。为了在这个实验中给出数字，我有 4000 个大小为 2048x1370 像素的投影，所以对于我的计算机技能来说，这是一个巨大的计算问题。

我需要你的帮助才能更快地执行一些操作。我在第一部分的代码分配了一个数组来包含所有图像，这是一个 34 Gb 的数组，但我有 130 Gb RAM，所以没问题。执行这个操作的代码使用了一个循环 imread:

for i=2:num_proj
    filename=strcat(path_im,list_proj(i).name);
    image=imread(filename);
    imArray(:,:,i)=image;
end

这当然不是最快的方法，现在创建这个数组需要 332 秒。我已经找到了几种解决方案来改善这一点，我会这样做的。

第二步是划分平场（没有样本拍摄的图像）。我的代码采用平面并使用 imdivide 将数组中的每个图像划分为平面图像：

for i=1:size(imArray, 3);
    imArray(:,:,i)=imdivide(imArray(:,:,i), flat);
end

这一步貌似很快，但是调用了4000次。你有什么建议吗？有没有更好的实现方式？

现在我最大的问题是如何在投影体积中以最快的方式进行水平切割？我的基本思路是从"top"读卷。报告如下，它完成了它的工作，但是需要很长时间：

for i=1:size(sinogram, 3)
        for k=1:size(sinogram, 1)
            sinogram(k,:,i)=imArray(i,:,k);

        end
end

你能帮我加快这个操作吗？我希望我的问题很清楚，否则请提问，我会尽力解释得更好。

Answer 1

你最终找到你想要的答案了吗？我将在这里介绍一些以供将来参考：

for i=2:num_proj
     filename=strcat(path_im,list_proj(i).name);
     image=imread(filename);
     imArray(:,:,i)=image;
end

改进第一个循环：

预分配 在 MATLAB 中很重要。如果您动态分配图像，每次向其中添加图像时，MATLAB 都必须创建一个新的 imArray。我假设您已经这样做了，因为循环从 2 开始并且速度相当快。如果没有，你应该调查一下。

parfor可以在一定程度上加速imread。根据我的经验，imread 似乎通常 I/O 受到限制，但我已经看到使用它的好处。

对于我在具有快速 ssd（假设 CPU-有限）的相当旧的双核笔记本电脑上的 100 张图像的基准测试：

imArray = cell([100,1])
parfor i=1:100
    imArray{i} = imread(filenames{i})
end
imArray = cat(3, imArray{:});

用了 47.506717 秒。

对比连续形式：

firstI = imread(files{1});
imArray = zeros(size(firstI,1),size(firstI,2),100, 'uint8');
imArray(:,:,1) = firstI;
for i=2:100
    imArray(:,:,i)=imread(filenames{i})
end

耗时 55.928427 秒。但是，您的情况可能不会像这样。您甚至可能会看到由于第一次启动 parpool 并在稍后运行s 分配工作人员的高开销成本而造成的损失。根据您的工作，这可能值得一试。我只是在第一个运行之后将我的 CT 数据集保存在 mat 文件中，使用这些数据集要快得多。

平场分割：

矢量化 可能很费力，有时收效甚微，因为现在 MATLAB 使用 JIT 编译器。然而，在某些情况下它仍然可以提供帮助，这里它只是一个函数调用：

for i=1:size(imArray, 3);
    imArray(:,:,i)=imdivide(imArray(:,:,i), flat);
end

用了 3.809438 秒

同时：

imArray = imArray ./ repmat(flat,1,1,size(imArray,3));

用了 1.268413 秒

甚至更好：

imArray = bsxfun(@rdivide, imArray, flat);

用了 0.970662 秒

bsxfun 似乎很受欢迎，这是有原因的。它默认是多线程的。

正弦图：

你只是在这里切换尺寸吗？如果是这样，请尝试：

sinogram = permute(imArray, [3,2,1]);

从断层扫描投影中提取正弦图

Extracting sinograms from tomography projections

performance

matlab

tomography-reconstruction