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

【问题标题】：Extracting sinograms from tomography projections从断层扫描投影中提取正弦图
【发布时间】：2016-01-18 19:00:49
【问题描述】：

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

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

我从实验中得到的是一系列不同角度的 2D 射线照相，您可以将其视为一个矩形“体积”，其中尺寸分别是单个投影中的行数和列数，体积来自角的数量。正弦图只是这个体积的横向切割。这意味着不是从侧面而是从顶部读取体积，因此我将创建一个新的图像数组，其尺寸是列数和投影，数组长度是投影中的行数。为了在这个实验中给出数字，我有 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 次。你有什么建议吗？有没有更好的方法来执行它？

现在我最大的问题是如何以最快的方式在投影体积中进行水平切割？我的基本想法是从“顶部”读取音量。报告如下，它完成了它的工作，但是需要很长时间：

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

        end
end

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

【问题讨论】：

我建议在 Computational Science 交流中提问：scicomp.stackexchange.com 这是一个相当复杂的数学问题，你可能会在那儿获得更好的运气。
我还是很乐观，我已经在这里找到了很多我的问题的答案:) 我会听从你的建议，我也会在那里问...
也只是在 CT-area 开始一个项目。有没有办法从点云对象 STL 数据创建图像投影？

标签： performance matlab tomography-reconstruction

【解决方案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 似乎通常受到 I/O 限制，但我已经看到了使用它的好处。

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

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 并在以后的运行中分配工作人员，您甚至可能会看到由于高昂的间接成本而造成的损失。根据您的工作，这可能值得一试。我只是在第一次运行后将我的 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]);

【讨论】：