I have a 3D array with curved lines (skeletonized data) and need to find the shortest connecting path between two points on those lines.
How do I find the shortest path between two points along a 3D skeleton array in Matlab?
There are several people asking a similar question both here and on the mathworks forum. But most answers point to this blog post on finding shortest paths, which does not search along a skeleton. If I try to force it along the skeleton it also returns paths that are not between the two points but farther away from both. Also some matlab functions that people suggest can only handle 2D data (like bwdistgeodesic, which could be used to evade non-skeleton voxels).
I use two functions from the file exchange. Here's what I do:
%get data
load mri
D = double(squeeze(D)); D = D./max(D(:)); D = real(ifft(ifft(ifft(ifftshift(ifftshift(ifftshift( padarray( fftshift(fftshift(fftshift(fft(fft(fft(D,[],1),[],2),[],3),1),2),3), size(D),0,'both' ) ,1),2),3),[],1),[],2),[],3));
figure; ax(1) = subplot(1,3,1); imshow(mean(D,3),[]); title('mean image')
%find vasculature
spacing = [1 1 1];sigmas = [0.5:0.5:2];tau = 1;whiteondark = false;Dvess = vesselness3D(D, sigmas, spacing, tau, whiteondark);
ax(2) = subplot(1,3,2); imshow(imfuse(squeeze(max(Dvess,[],3)), squeeze(max(Dvess>0.5,[],3))),[]); title('masked "vessels"')
%skeletonize
Dskel = Skeleton3D(Dvess>0.5);
ax(3) = subplot(1,3,3); imshow(imfuse(squeeze(max(Dskel,[],3)), squeeze(max(Dvess>0.5,[],3))),[]); title('skeletonized mask'); linkaxes(ax); zoom(2);
%select seeds
corline1 = imline(); idx1 = find(repmat(corline1.createMask, [1 1 size(Dskel,3)]).*Dskel); [sub11, sub21, sub31] = ind2sub(size(Dskel), idx1);
corline2 = imline(); idx2 = find(repmat(corline2.createMask, [1 1 size(Dskel,3)]).*Dskel); [sub12, sub22, sub32] = ind2sub(size(Dskel), idx2);
So far so good with these results:
The white lines in the right image show the skeleton and the blue lines show the seeds. But then if I try to apply distance transforms I cant find the shortest path:
%distance map
dist1 = zeros(size(Dskel)); dist1(sub11, sub21, sub31) = 1; dist1 = bwdist(dist1);
dist2 = zeros(size(Dskel)); dist2(sub12, sub22, sub32) = 1; dist2 = bwdist(dist2);
dist12 = dist1+dist2;
%path
path1 = imregionalmin(Dskel.*dist12);
path2 = imregionalmin(Dskel.*dist12+10000.*(~Dskel));
%path3 = imregionalmin(Dskel.*dist12+Inf.*Dskel); %function cant handle inf
figure; ax(1) = subplot(1,3,1); imshow(squeeze(sum(double(path1),3)),[]); title('')
ax(2) = subplot(1,3,2); imshow(squeeze(sum(double(path2),3)),[]); title('')
%ax(2) = subplot(1,3,2); imshow(squeeze(sum(double(path3),3)),[]); title('')
linkaxes(ax); zoom(3);
The first finds finds pretty much the entire skeleton again and the second looks like bogus. The regional minimum definitely is not appropriate to find the path, but not sure what else would be. Also my distance calculation does not account for the skeleton, which I tried to circumvent in the second path but that also did not work nicely.
What should I do different?
Just found something called 'fast marching', maybe that would work. Looking into it now.