Extracting sinograms from tomography projections

Question

my name is Lorenzo and I am a postdoc researcher in Italy. My job is related to tomography imaging using synchrotron radiation. This is a new field for me and I start to face with Matlab to write some code.

I am completely new with tomography imaging and Matlab shows new challenges for me. My actual problem is to create sinograms from a stacking of parallel image projections. For whose are not in the field, a sinogram is a map detecting the projection on the detector of the position of a feature in the sample, as function of angle between the X-Ray beam and the sample.

What I get from an experiment is a series of 2D radiography at different angles, you can see it as a rectangular "volume" where dimensions are the number of rows and columns in the single projection, respectively, and volume is given from the number of angles. A sinogram is just a transversal cut of this volume. This means reading the volume not from lateral side but from the top, so I will create a new array of images, whose dimensions are the number of columns and projection, and array lenght the number of rows in the projections. To give numbers in this experiment I have 4000 projections with size 2048x1370 pixels, so for my computer skills this is a huge computational problem.

I need your help to perform some operations faster. My code in a first part allocates an array to contain all images, this is a 34 Gb array but I have 130 Gb RAM, so no problem. The code to perform this operation uses a loop of imread:

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

This is not the fastest way off course, now it takes 332 seconds to create this array. I have found several solutions to improve this and I will do it.

Second step is to divide for a flat field (image taken without the sample). My code takes the flat and using imdivide divides each image in the array for the flat image:

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

this step seems to be fast however it is called 4000 times. Do you have suggestion? Is there a better way to perform it?

Now my greatest problem is how to take in the fastest way the horizontal cut in the projection volume? My basic idea to read the volume from the "top". is reported below, it does its job, however it takes long time:

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

        end
end

Can you help me to speed up also this operation? I hope my question is clear, otherwise please ask and I will try to explain it better.

Answer 1

Did you eventually find the answers you are looking for? I'll go through some here for future reference:

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

Improving the first loop:

Pre-allocation is critical in MATLAB. If you allocate an image dynamically, MATLAB will have to create a new imArray for each time you add an image to it. I'm assuming that you do it already as the loop starts from 2 and is fairly fast. If not, you should look into it.

parfor can be used to accelerate imread to a degree. In my experience imread seems to usually be I/O limited, but I have seen gains using it.

For my benchmark of a 100 images on a fairly old dual-core laptop with a fast ssd (assumedly CPU-limited):

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

took 47.506717 seconds.

vs. a serial form:

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

Which took 55.928427 seconds. However, your case will probably not look like this. You might even see losses due to the high overhead costs from starting up the parpool the first time and allocating the workers on later runs. It might be worth a try depending on your work. I just save my CT-datasets in mat-files after the first run, it's much faster to work with those.

dividing for flat field:

Vectorization can be laborous and sometimes with little gains, due to MATLAB using a JIT-compiler nowadays. However, in some cases it can still help and here it is only a single function call away:

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

took 3.809438 seconds

While:

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

took 1.268413 seconds

Or even better:

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

took 0.970662 seconds

bsxfun seems to be popular, for a reason. It's multithreaded by default.

The sinogram:

Are you simply switching the dimensions here? If so, try:

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