2

我叫 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

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

4

1 回答 1

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 的相当旧的双核笔记本电脑上的 100 张图像的基准测试（假设 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]);
于 2016-08-14T14:43:36.657 回答