How can I speed this up (it takes days to run currently)?

I have a piece of code that I wrote that takes a 3D matrix containing imaging data of brain vessels and uses region growing starting at a user-defined seed. However it takes to long. I'm sure there's a way to speed this up as some pixels are reprocessed (found in the queue) multiple times. Also I think it's taking to long to process the pixels (my image contains approx. 1.6E7 pixels) however only a small percentage of them are "of interest", i.e. those that are at the minimum should be excluded from the computation as the probability that they belong to the volume of interest is 0%. Following is the code:

clc; clear all; close all;

load('Jmatrix.mat'); %this file contains a 3D matrix with the grayscale vessel slices (from DICOM).

dcmmov = implay(J); %allows the user to scan the images in order to select the slice containing the aneurysm.

waitfor(dcmmov); %waits for the user to close the image scanner (movie) above.

prompt = {'Enter slice number:'};
dlg_title = 'Slice Selection';
num_lines = 1;
slice = inputdlg(prompt,dlg_title,num_lines); %user inputs slice of interest for seeding.

z = str2double(slice{1});

imshow(J(:,:,z)); %shows the slice of interst
[y,x] = ginput; %permits the user to seed a point within the aneurysm to start the growing from.

% seed = [x,y];

% seed = [165,85];

x = round(x); %makes sure the co-ordinates are integers
y = round(y);

[h, w, d] = size(J); %gets the dimensions of the original image matrix.

V = zeros(w,h,d); %creates a prealocated matrix of the same size as the original (will be binary).
V(x,y,z)= 1; %assigns a 1 (i.e. true) to the seed pixel
error = .10; %current error threshold to be used when compairing neighbouring pixels.

queue = [x, y, z]; %creates the que for the pixels to be analysed.
bar = waitbar(0,'Segmenting image. Please wait...'); %starts a progress bar.

while size(queue, 1)

    voxelx = queue(1,1);  %indicate co-ordinates of the next pixel (voxel) in the queue.
    voxely = queue(1,2);
    voxelz = queue(1,3);

    S = double(J(voxelx, voxely, voxelz)); %this is the pixel value to be compared with neighbours (originally the seed).

    queue(1,:) = []; %deletes the queue

    for i = -1:1
      for j = -1:1
          for k = -1:1
              if  J(voxelx+i,voxely+j,voxelz+k) == 1;      %skips to the next text pixel if the value of this pixel is already 1 (removes redundancy).
                  V(voxelx+i,voxely+j,voxelz+k)=1;
              elseif J(voxelx+i,voxely+j,voxelz+k) == 0;   %skips to the next text pixel if the value of this pixel is already 1 (removes redundancy).
                  V(voxelx+i,voxely+j,voxelz+k)=1;
              else
                  T = double(J(voxelx+i, voxely+j,voxelz+k)); %defines the neighbouring pixel to be tested.
                  if abs(T - S)/S < error   %this is the formula that determines if the tested pixel is sufficiently similar (within an error threshold) to the comparison pixel S. 
                   V(voxelx+i,voxely+j,voxelz+k) = 1; %if sufficiently similar assigns a 1 (i.e. true).
                  end

              end
              queue(end+1,:) = [voxelx+i, voxely+j,voxelz+k]; %adds the next pixel to the queue.
              sizeq = size(queue); %these lines allow monitoring of the process
              waitbar(sizeq(1) / (w*h*d));
              pixel = [voxelx+i,voxely+j,voxelz+k];
              value = V(voxelx+i,voxely+j,voxelz+k);
              [sizeq(1) pixel value]
          end
      end
    end
  end
  close(bar);          

implay(V); %permits viewing of the final segmented volume.
  end

Thanks, Daniel

[Matt J]

It doesn't look like there's any way that your queue can empty. In one line you delete a voxel

    queue(1,:) = []; %deletes the queue

but then later, in your triple for-loop, you re-include the same voxel plus 26 more in

queue(end+1,:) = [voxelx+i, voxely+j,voxelz+k];

Aside from that, when you say that there are some voxels that are not of interest, you should expect that only you (not us) are going to know which voxels those are.

[Daniel]

Hi Matt, thanks for your answer. I am not quite sure what you mean by saying that only I know which voxels are of interest. What I mean by voxels of interest are those which are part of the geometry (i.e. blood vessels with radiological contrast inside). The images that I am using are generated by digital subtraction angiography which are similar to CT images, hence tomographic (in axial slices). The vessels "light up" and everything surrounding them is black (or nearly 0 pixel value) since it has been subtracted from a non-contrast run, therefore the only feature visible should be the contrast. Nevertheless, these images contain a significant amount of noise, hence my need to segment the vessels with a volume growing algorithm. The black (or nearly black) pixels are nothing or noise and are therefore not of interest, those with high(er) pixel values are likely part of the feature of interest (i.e. vessel).

[Matt J]

I am not quite sure what you mean by saying that only I know which voxels are of interest.

I meant, "how would we know better than you"? Where in your code do you specify the thresholding criteria that determines whether a particular voxel belongs to a vessel?

It sort of looks like this is the criterion

 abs(T - S)/S < error

but then why do you also set V=1 when J=1 and also V=1 when J=0?

And why do you always add the currently examined voxel to the queue regardless of what's inside it?

queue(end+1,:) = [voxelx+i, voxely+j,voxelz+k]; %adds the next pixel to the queue.

Answer 1

[Matt J]

Below is an idea for getting rid of your triple for-loop. You'll notice that I don't use (i,j,k) voxel coordinates, but rather linear indices instead.

Also, your queue is not pre-allocated and both shrinks and grows throughout the while loop. This is a notoriously bad thing. My suggestion is to instead use a binary queueMask which is true in voxels that need processing and false elsewhere.

Finally, you should be aware that displaying things on the screen is time consuming. Removing statements like

[sizeq(1) pixel value]

will make things go a lot faster.

queueMask=false(size(J)); %binary masks of voxels left to process
  queueMask(x,y,z)=true; %seed

queuelength=nnz(queueMask);  

[i,j,k]=ndgrid(1:3,1:3,1:3);
jumps=sub2ind(size(J),i,j,k)-sub2ind(size(J),2,2,2);

while queuelength>0

      voxel = find(queueMask,1);
      neighbours = voxel + jumps;

      T=double(J(neighbours));
      S=T(14);
      V(neighbours) = J==1|J==0|(abs(T - S)/S < error);

      %I don't know if this is how the queue should be updated, but it's an
      %example
      queueMask(neighbours)=true;
      queueMask(voxel)=false;
      queuelength=queuelength+...

end

[Daniel]

Thanks again Matt for your valuable input. I have been analysing what you did and get the general idea. If I am not mistaken you are moving a test 3x3x3 mask around and comparing that with the central "voxel". Nevertheless, I am unable to make this work. I see that V(neighbours) is a logical statement but I cannot make sense of the relationship or what this variable is supposed to represent. The error that I get when running this is related to:

size(J==1) = size(J==0) ≠ size(abs(T-S)/S)

Finally, I myself can't think of a way to update the queue

queuelength=queueLength+...

is incomplete I assume).

Thanks again!

[Matt J]

Sorry, V(neighbours) should be updated as

 V(neighbours) = T==1|T==0|(abs(T - S)/S < error);

or at least this looks equivalent to what you had in your code.

Answer 2

[Image Analyst]

The things that strikes me right ott the bat is that your indexes are in the reverse order. MATLAB is in "column major order". so that it goes down rows in a column before it moves over to the next column. So in imageArray(row, column) the left-most index, row, varies fastest and the rightmost one varies slowest. You have your loops in the opposite order so that is going to make it the most memory inefficient that it can be. Try putting your "i" loop inner-most and your k loop outer-most. Time it and see if it's faster. I bet it will be.

[Daniel]

Thanks for your answer. I tested changing the loops as you recommended and clocked a set number of iterations. Unfortunately the difference in computation time is insignificant. However, I'll keep this suggestion into consideration when nesting for loops.