Any ideas on debugging this segmentation code?

I have (with input from various people) developed this code for segmentation of angiographic images (already filtered, windowed and thresholded). However it doesn't work correctly. It works as follows:

1. The images are preloaded Jnorm (please ask me to send you the sample matrix containing an idealised angiogram if needed).
2. When the code is run it first displays a movie of all the slices so that the feature of interest to seed (in theory an aneurysm) can be found.
3. The movie viewer is closed and a window asking for the slice of interest is displayed. The number of the slice containing the feature is entered.
4. An image viewer containing said slice is displayed and the mouse contains crosshairs. The feature is clicked on to "seed" it. Press enter after clicking the feature to continue.
5. The segmentation runs...
6. The waitbar (should) close and a movie (should) display with the now segmented image.

The problem is that the segmentation progresses downwards (i.e. to lower number slices) from where the seed originates and never moves upwards. Also it creates a segmentation "path" before it actually starts segmenting. It seems to need to hit a "wall" before it starts segmenting all the connected features (i.e. grow) in that slice.

dcmmov = implay(Jnorm); %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(Jnorm(:,:,z)); %shows the slice of interst
[y,x] = ginput; %permits the user to seed a point within the aneurysm to start the growing from.

x = round(x); %makes sure the co-ordinates are integers
y = round(y);
error = .15; %current error threshold to be used when compairing neighbouring pixels.

Jlinear = reshape(Jnorm,1,[]);
V = zeros(size(Jnorm));

queueMask=false(size(Jnorm)); %binary masks of voxels left to process
queueMask(x,y,z)=true; %seed
queuelength=nnz(queueMask);  
[i,j,k]=meshgrid(1:3,1:3,1:3);
jumps=sub2ind(size(Jnorm),i,j,k)-sub2ind(size(Jnorm),2,2,2);
remaining = size(find(Jnorm>0),1);
%changed it due to 'Undefined function 'lenght' for input arguments of type 'double'.' error originally read: 'remaining=length(find(Jnorm>0))'.
voxelqueue = [];
bar = waitbar(0,'Segmenting image. Please wait...'); %starts a progress bar.

while queuelength>0
    voxel = find(queueMask,1);
    neighbours = voxel + jumps;
    if sum(ismember(voxelqueue,voxel)) == 0 && size(find(neighbours > 0),1) == 27 %~ismember(neighbours,voxels)
        T=double(Jlinear(neighbours)); %remember to do this (make sure T is positive).
        S=T(14);
        V(neighbours) = (abs(T - S)/S < error);
        queueMask(neighbours) = (abs(T - S)/S < error);
        %disp(['Next interation: voxel ', int2str(voxel), ' ; length of queue ', int2str(queuelength)]);
    end
    queueMask(voxel)=false;
    queuelength = nnz(queueMask);
    voxelqueue(end+1,:) = voxel;
    waitbar(size((voxelqueue),1)/remaining);
end

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

Answer 1

[Image Analyst]

Identify which slice is the "wall" slice. Then examine the values of that slice. It's probably anomalous for some reason.