{"id":9918,"date":"2018-07-06T09:00:53","date_gmt":"2018-07-06T13:00:53","guid":{"rendered":"https:\/\/blogs.mathworks.com\/pick\/?p=9918"},"modified":"2018-07-05T16:20:18","modified_gmt":"2018-07-05T20:20:18","slug":"a-couple-of-my-favorite-new-image-processing-toolbox-functions","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/pick\/2018\/07\/06\/a-couple-of-my-favorite-new-image-processing-toolbox-functions\/","title":{"rendered":"A couple of my favorite new Image Processing Toolbox functions"},"content":{"rendered":"
<\/p>\n

In July, 2016, I selected Philip Kollmansberger<\/a>’s useful File Exchange contribution Skeleton3D<\/a> as the Pick of the Week. I’d like to revisit that particular blog post today, if I may. For those of you who don’t yet have access to R2018a, Phillip’s contribution is still perhaps your best option for skeletonizing volumetric imagery. But a new function added to the Image Processing Toolbox<\/a> in the current release (that is, R2018a) gives you provides an alternative. In the Release Notes documenting modifications to the Toolbox, we learn that there is a new function called bwskel<\/tt><\/a> that supports skeletonization of 2D and 3D images. This is big!<\/p>\n

In light of this new function, it is interesting to take a fresh look at the 3D skeletonization of my earlier Pick.<\/p>\n

Consider:<\/p>\n

%% Load and display sample 3D data<\/span>\r\nload spiralVol.mat<\/span>;\r\nstepSize = 25;\r\n% Detect all spiral points:<\/span>\r\ninds = find(spiralVol);\r\n% Note the order reversal for x,y (col\/row):<\/span>\r\n[y,x,z] = ind2sub(size(spiralVol),inds);\r\nsubplot(2,2,1)\r\nmyTitle = sprintf('Full Spiral (Downsampled):\\n %i Points'<\/span>,...<\/span>\r\n    numel(inds));\r\nrefreshSpiralVisualization(x,y,z,stepSize,myTitle);\r\n<\/pre>\n
%% First, skeletonization using Skeleton3D (grayscale):<\/span>\r\ntic\r\nskel = Skeleton3D(spiralVol);\r\nt = toc;\r\ninds = find(skel);\r\nmyTitle = sprintf('Using Skeleton3D:\\n%i \"Skeleton\" Points; t = %0.2f sec '<\/span>,...<\/span>\r\n    numel(inds),t);\r\ndisp('Using Skeleton3D, grayscale'<\/span>)\r\nsubplot(2,2,2)\r\nrefreshSpiralVisualization(x,y,z,stepSize,myTitle);\r\nhold on<\/span>\r\n[yskel,xskel,zskel] = ind2sub(size(skel),inds);\r\nscatter3(xskel,yskel,zskel,5,...<\/span>\r\n    'Marker'<\/span>,'o'<\/span>,...<\/span>\r\n    'MarkerEdgeColor'<\/span>,'none'<\/span>,...<\/span>\r\n    'MarkerFaceColor'<\/span>,'r'<\/span>);\r\n<\/pre>\n
%% Next, skeletonization using Skeleton3D (binary):<\/span>\r\ntic\r\nspiralVolBinary = imbinarize(spiralVol);\r\nskel = Skeleton3D(spiralVolBinary);\r\nt = toc;\r\ninds = find(skel);\r\nmyTitle = sprintf('Using Skeleton3D (Binary):\\n%i \"Skeleton\" Points; t = %0.2f sec '<\/span>,...<\/span>\r\n    numel(inds),t);\r\n[yskel,xskel,zskel] = ind2sub(size(skel),inds);\r\nsubplot(2,2,3)\r\nrefreshSpiralVisualization(x,y,z,stepSize,myTitle);\r\nhold on<\/span>\r\nscatter3(xskel,yskel,zskel,5,...<\/span>\r\n    'Marker'<\/span>,'o'<\/span>,...<\/span>\r\n    'MarkerEdgeColor'<\/span>,'none'<\/span>,...<\/span>\r\n    'MarkerFaceColor'<\/span>,'r'<\/span>);\r\n<\/pre>\n
%% Finally, using |bwskel()| (Binary Only)<\/span>\r\ntic\r\nspiralVolLogical = imbinarize(spiralVol);\r\nskel = bwskel(spiralVolLogical);\r\nt = toc;\r\ninds = find(skel);\r\nmyTitle = sprintf('Using bwskel() (Binary):\\n%i \"Skeleton\" Points; t = %0.2f sec '<\/span>,...<\/span>\r\n    numel(inds),t);\r\n[yskel,xskel,zskel] = ind2sub(size(skel),inds);\r\nsubplot(2,2,4)\r\nrefreshSpiralVisualization(x,y,z,stepSize,myTitle);\r\nhold on<\/span>\r\nscatter3(xskel,yskel,zskel,5,...<\/span>\r\n    'Marker'<\/span>,'o'<\/span>,...<\/span>\r\n    'MarkerEdgeColor'<\/span>,'none'<\/span>,...<\/span>\r\n    'MarkerFaceColor'<\/span>,'r'<\/span>);\r\n<\/pre>\n
%%<\/span>\r\nfunction<\/span> refreshSpiralVisualization(x,y,z,stepVal,myTitle)\r\nscatter3(x(1:stepVal:end),y(1:stepVal:end),z(1:stepVal:end),0.1,...<\/span>\r\n    'Marker'<\/span>,'o'<\/span>,...<\/span>\r\n    'MarkerEdgeColor'<\/span>,'c'<\/span>,...<\/span>\r\n   'MarkerFaceAlpha'<\/span>,0.3,...<\/span>\r\n    'MarkerFaceColor'<\/span>,'c'<\/span>);\r\nxlabel('x'<\/span>);\r\nylabel('y'<\/span>);\r\nzlabel('z'<\/span>);\r\nset(gca,'view'<\/span>,[-26.94  23.6])\r\ntitle(myTitle)\r\nend<\/span>\r\n<\/pre>\n
\"\" <\/p>\n
A few things come become apparent looking at these results. First, Skeleton3D()<\/tt> operates without error or warning on grayscale images (upper right axes), but appears to label every nonzero input pixel as a “skeleton” point. More usefully, as suggested by Phillip’s description of the algorithm, the function returns more reasonable skeleton points when operating on a binarization of the input data (lower left axes). However, that calculation took over 25 seconds (including binarization) on my computer. In contrast, bwskel()<\/tt> (lower right) errors (usefully, in my opinion) if you try to operate directly on the grayscale data, and takes less than a second (including binarization) operating on the segmented data. (The outputs of the two functions are not identical, but they both seem reasonable.)<\/p>\n
In the 2016 post, I also mentioned that the calculation of blood vessel tortuosity<\/a> in 3D was compounded by the challenge of determining endpoints of 3D skeletons. That difficulty, too, has fallen by the wayside. The new function bwmorph3<\/tt><\/a> (also in the Image Processing Toolbox) now facilitates the calculation of a half-dozen operations on 3D binary images; fortunately, that list of operations includes endpoint-detection. The net result: if we had a volumetric representation of the vasculature we wanted to evaluate, we could now do so relatively easily.<\/p>\n
One final comment: note in the code above that the script contains a local function, or subfunction. The capacity to include subfunctions in scripts was added in R2016b.<\/p>\n
As always, I welcome your thoughts and comments<\/a>.<\/p>\n