MATLAB: How to speed this up (it takes days to run currently)

for loopImage Processing Toolboximage segmentationmedical imagesregion growingspeedvolume growing

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

Best Answer

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
Related Question