从层析成像中提取正弦波

Question

我叫洛伦佐，是意大利的博士后研究员。我的工作是利用同步辐射进行断层成像。这是我的一个新领域，我开始用Matlab编写一些代码。

我是全新的断层成像和Matlab给我带来了新的挑战。我的实际问题是从一堆平行的图像投影中创建汉字图。对于不在现场的情况下，正弦波图是一幅地图，它以X射线束和样品之间的角度作为函数，在探测器上探测到一个特征在样品中的位置。

我从一个实验中得到的是一系列不同角度的二维射线摄影，你可以看到它是一个矩形的“体积”，其中尺寸分别是单个投影中的行数和列数，而体积是从角度的数目给出的。正弦波只是这个体积的横切。这意味着不是从侧面而是从顶部读取卷，所以我将创建一个新的图像数组，其尺寸是列数和投影数，数组的长度是投影中的行数。为了给出这个实验中的数字，我有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秒。我已经找到了几个解决办法来改善这一点，我将这样做。

第二步是对一个平面场(没有样本的图像)进行分割。我的代码获取了平面，并使用了对平面图像数组中的每个图像进行划分的方法：

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

你能帮我加速这次行动吗？我希望我的问题是清楚的，否则请问，我会尽量解释得更好。

Answer 1

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

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

改进第一个循环：

Pre-allocation是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秒。然而，你的案子可能不会像这样。您甚至可能会发现，由于第一次启动共享池和在以后运行时分配员工所造成的高开销成本而造成的损失。这可能值得一试取决于你的工作。我只是保存我的CT-数据集在垫-文件后，第一次运行，它是更快的工作，这些。

平坦场的划分：

由于现在使用JIT编译器的MATLAB，Vectorization可能是费力的，有时收益很小。但是，在某些情况下，它仍然可以提供帮助，在这里，它只是一个函数调用：

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]);