%------------------------------------------------------------------------------------------------ % Demo to illustrate simple blob detection, measurement, and filtering. % Requires the Image Processing Toolbox (IPT) because it demonstates some functions % supplied by that toolbox, plus it uses the "coins" demo image supplied with that toolbox. % If you have the IPT (you can check by typing ver on the command line), you should be able to % run this demo code simply by copying and pasting this code into a new editor window, % and then clicking the green "run" triangle on the toolbar. % Running time = 7.5 seconds the first run and 2.5 seconds on subsequent runs. % A similar Mathworks demo: % http://www.mathworks.com/products/image/demos.html?file=/products/demos/shipping/images/ipexprops.html % Code written and posted by ImageAnalyst, July 2009. % Modified for MRI scans by C Rahul, February 2014. %------------------------------------------------------------------------------------------------ % function BlobsDemo() % echo on; % Startup code. tic; % Start timer. clc; % Clear command window. clear all; % Get rid of variables from prior run of this m-file. disp('Running BlobsDemo.m...'); % Message sent to command window. workspace; % Show panel with all the variables. % Change the current folder to the folder of this m-file. % (The line of code below is from Brett Shoelson of The Mathworks.) if(~isdeployed) cd(fileparts(which(mfilename))); end % Read in standard MATLAB demo image originalImage = rgb2gray(imread('4.jpg')); subplot(3, 3, 1); imshow(originalImage); % Maximize the figure window. set(gcf, 'Position', get(0, 'ScreenSize')); % Force it to display RIGHT NOW (otherwise it might not display until it's all done, unless you've stopped at a breakpoint.) drawnow; caption = sprintf('Original MRI Image'); title(caption); axis square; % Make sure image is not artificially stretched because of screen's aspect ratio. % Just for fun, let's get its histogram. [pixelCount grayLevels] = imhist(originalImage); subplot(3, 3, 2); bar(pixelCount); title('Histogram of original image'); xlim([0 grayLevels(end)]); % Scale x axis manually. % Filtering starts % Setting all pixels more than 180 and less than 60 to 0 to create a clear distinction % between spine fluid and spine for i=1:512 for j=1:512 if originalImage(i,j)> 100 || originalImage(i,j)< 63 originalImage(i,j) = 0; end end end % Threshold the image to get a binary image (only 0's and 1's) of class "logical." thresholdValue = 63; binaryImage = originalImage > thresholdValue; % Bright objects will be the chosen if you use >. % binaryImage = originalImage < thresholdValue; % Dark objects will be the chosen if you use <. % Do a "hole fill" to get rid of any background pixels inside the blobs. binaryImage = imfill(binaryImage, 'holes'); % Show the threshold as a vertical red bar on the histogram. hold on; maxYValue = ylim; hStemLines = stem(thresholdValue, maxYValue(2), 'r'); children = get(hStemLines, 'children'); set(children(2),'visible', 'off'); % Place a text label on the bar chart showing the threshold. annotationText = sprintf('Thresholded at %d gray levels', thresholdValue); % For text(), the x and y need to be of the data class "double" so let's cast both to double. text(double(thresholdValue + 5), double(0.5 * maxYValue(2)), annotationText, 'FontSize', 10, 'Color', [0 .5 0]); text(double(thresholdValue - 70), double(0.94 * maxYValue(2)), 'Background', 'FontSize', 10, 'Color', [0 0 .5]); text(double(thresholdValue + 50), double(0.94 * maxYValue(2)), 'Foreground', 'FontSize', 10, 'Color', [0 0 .5]); subplot(3, 3, 3); imagesc(originalImage); title('Binary Image, obtained after initial spine fluid segmentation'); axis square; % Display the binary image. subplot(3, 3, 4); imagesc(binaryImage); colormap(gray(256)); title('Binary Image, obtained by thresholding'); axis square; labeledImage = bwlabel(binaryImage, 8); % Label each blob so we can make measurements of it coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels subplot(3, 3, 5); imshow(labeledImage, []); title('Labeled Image, from bwlabel()'); axis square; subplot(3, 3, 6); imagesc(coloredLabels); axis square; caption = sprintf('Pseudo colored labels, from label2rgb().\nBlobs are numbered from top to bottom, then from left to right.'); title(caption); % Get all the blob properties. Can only pass in originalImage in version R2008a and later. blobMeasurements = regionprops(labeledImage, originalImage, 'all'); numberOfBlobs = size(blobMeasurements, 1); % bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image. % Plot the borders of all the blobs on the original grayscale image using the coordinates returned by bwboundaries. subplot(3, 3, 7); imagesc(originalImage); title('Outlines, from bwboundaries()'); axis square; hold on; [boundaries ,L,N,A] = bwboundaries(binaryImage,8,'noholes'); colors=['b' 'g' 'r' 'c' 'm' 'y']; numberOfBoundaries = size(boundaries); poly_eq = zeros(numberOfBoundaries(1),3); % Setting the equations matrix for k = 1 : numberOfBoundaries boundary = boundaries{k}; cidx = mod(k,length(colors))+1; plot(boundary(:,2), boundary(:,1),... colors(cidx),'LineWidth',2); %randomize text position for better visibility rndRow = ceil(length(boundary)/(mod(rand*k,7)+1)); col = boundary(rndRow,2); row = boundary(rndRow,1); h = text(col+1, row-1, num2str(L(row,col))); set(h,'Color',colors(cidx),... 'FontSize',6,'FontWeight','bold'); % getting the polynomial equation for each blob num_of_points_in_blob = size(boundary); if(num_of_points_in_blob(1) >=4) eq = polyfit(boundary(:, 1),boundary(:, 2),2); poly_eq(k,:) = [eq(1,1) eq(1,2) eq(1,3)]; end end spy(A); hold on; fontSize = 6; % Used to control size of "blob number" labels put atop the image. labelShiftX = -7; % Used to align the labels in the centers of the coins. blobECD = zeros(1, numberOfBlobs); % Print header line in the command window. fprintf(1,'Blob # Mean Intensity Area Perimeter Centroid Diameter\n'); % Loop over all blobs printing their measurements to the command window. for k = 1 : numberOfBlobs % Loop through all blobs. % Find the mean of each blob. (R2008a has a better way where you can pass the original image % directly into regionprops. The way below works for all versions including earlier versions.) thisBlobsPixels = blobMeasurements(k).PixelIdxList; % Get list of pixels in current blob. meanGL = mean(originalImage(thisBlobsPixels)); % Find mean intensity (in original image!) meanGL2008a = blobMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a blobArea = blobMeasurements(k).Area; % Get area. blobPerimeter = blobMeasurements(k).Perimeter; % Get perimeter. blobCentroid = blobMeasurements(k).Centroid; % Get centroid. blobECD(k) = sqrt(4 * blobArea / pi); % Compute ECD - Equivalent Circular Diameter. fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, blobArea, blobPerimeter, blobCentroid, blobECD(k)); % Put the "blob number" labels on the "boundaries" grayscale image. text(blobCentroid(1) + labelShiftX, blobCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold'); end % Put the labels on the rgb labeled image also. subplot(3, 3, 5); for k = 1 : numberOfBlobs % Loop through all blobs. blobCentroid = blobMeasurements(k).Centroid; % Get centroid. text(blobCentroid(1) + labelShiftX, blobCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold'); end % FILTER % Find only those blobs % with a diameter between 150 and 220 and an area more than 110. allBlobDiameters = [blobMeasurements.EquivDiameter]; allBlobAreas = [blobMeasurements.Area]; % Get a list of the blobs that meet our criteria and we need to keep. allowableDiameterIndexes = (allBlobDiameters > 20) & (allBlobDiameters < 120); allowableAreaIndexes = (allBlobAreas > 110) & (allBlobAreas < 600) ; % Take the small objects. keeperIndexes = find(allowableAreaIndexes & allowableDiameterIndexes); % Extract only those blobs that meet our criteria, and % eliminate those blobs that don't meet our criteria. % Note how we use ismember() to do this. keeperBlobsImage = ismember(labeledImage, keeperIndexes); % Re-label with only the keeper blobs kept. labeledDimeImage = bwlabel(keeperBlobsImage, 8); % Label each blob so we can make measurements of it % Now we're done. We have a labeled image of blobs that meet our specified criteria. subplot(3, 3, 8); imshow(labeledDimeImage, []); axis square; title('Blobs after area-filtering'); % Now use the keeper blob as a mask on the original image. % This will let us display the original image in the regions of the keeper blobs. maskedImageDime = originalImage; % Simply a copy at first. maskedImageDime(~keeperBlobsImage) = 0; % Set all non-keeper pixels to zero. subplot(3, 3, 9); imshow(maskedImageDime); axis square; title('Only the 3 brightest dimes from the original image'); % Now let's get the nickels (the larger coin type) keeperIndexes = find(allBlobAreas > 300); % Take the larger objects. % Note how we use ismember to select the blobs that meet our criteria. nickelBinaryImage = ismember(labeledImage, keeperIndexes); maskedImageNickel = originalImage; % Simply a copy at first. maskedImageNickel(~nickelBinaryImage) = 0; % Set all non-nickel pixels to zero. %subplot(3, 3, 9); %imshow(maskedImageNickel, []); %axis square; %title('Only the nickels from the original image'); elapsedTime = toc; % Alert user that the demo is done and give them the option to save an image. % message = sprintf('Finished running BlobsDemo.m.\n\nElapsed time = %.2f seconds.', elapsedTime); % message = sprintf('%s\n\nCheck out the figure window for the images.\nCheck out the command window for the numerical results.', message); % message = sprintf('%s\n\nDo you want to save the pseudo-colored image?', message); % reply = questdlg(message, 'Save image?', 'Yes', 'No', 'No'); % % Note: reply will = '' for Upper right X, 'Yes' for Yes, and 'No' for No. % if strcmpi(reply, 'Yes') % % Ask user for a filename. % FilterSpec = {'*.tif', 'TIFF images (*.tif)'; '*.*', 'All Files (*.*)'}; % DialogTitle = 'Save image file name'; % % Get the default filename. Make sure it's in the folder where this m-file lives. % % (If they run this file but the cd is another folder then pwd will show that folder, not this one. % thisFile = mfilename('fullpath'); % [thisFolder, baseFileName, ext, version] = fileparts(thisFile); % DefaultName = sprintf('%s/%s.tif', thisFolder, baseFileName); % [fileName, specifiedFolder] = uiputfile(FilterSpec, DialogTitle, DefaultName); % % Parse what they actually specified. % [folder, baseFileName, ext, version] = fileparts(fileName); % % Create the full filename, making sure it has a tif filename. % fullImageFileName = fullfile(specifiedFolder, [baseFileName '.tif']); % % Save the labeled image as a tif image. % imwrite(uint8(coloredLabels), fullImageFileName); % % Just for fun, read image back into the imtool utility to demonstrate that tool. % tifimage = imread(fullImageFileName); % imtool(tifimage, []); % end % % message = sprintf('Would you like to crop out each coin to individual images?'); % reply = questdlg(message, 'Extract Individual Images?', 'Yes', 'No', 'Yes'); % % Note: reply will = '' for Upper right X, 'Yes' for Yes, and 'No' for No. % if strcmpi(reply, 'Yes') % figure; % % Maximize the figure window. % set(gcf, 'Position', get(0, 'ScreenSize')); % for k = 1 : numberOfBlobs % Loop through all blobs. % % Find the bounding box of each blob. % thisBlobsBoundingBox = blobMeasurements(k).BoundingBox; % Get list of pixels in current blob. % % Extract out this coin into it's own image. % subImage = imcrop(originalImage, thisBlobsBoundingBox); % % Display the image with informative caption. % subplot(3, 4, k); % imshow(subImage); % caption = sprintf('Coin #%d\nDiameter = %.1f pixels\nArea = %d pixels', k, blobECD(k), blobMeasurements(k).Area); % title(caption, 'FontSize', 14); % end %end

Running BlobsDemo.m... Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Warning: Polynomial is badly conditioned. Add points with distinct X values, reduce the degree of the polynomial, or try centering and scaling as described in HELP POLYFIT. Blob # Mean Intensity Area Perimeter Centroid Diameter # 1 76.1 191.0 60.6 41.2 230.2 15.6 # 2 50.2 3648.0 953.5 136.9 375.0 68.2 # 3 69.6 11.0 10.7 54.6 242.9 3.7 # 4 65.1 11.0 9.6 68.5 122.0 3.7 # 5 67.3 20.0 16.8 68.2 191.3 5.0 # 6 66.1 90.0 53.5 76.1 126.5 10.7 # 7 95.6 36.0 28.4 103.7 252.2 6.8 # 8 68.4 21.0 17.0 122.2 225.7 5.2 # 9 79.6 91.0 35.0 131.6 227.9 10.8 #10 69.1 97.0 40.9 151.3 434.9 11.1 #11 74.6 70.0 31.3 164.6 427.2 9.4 #12 65.0 4.0 5.9 172.0 135.5 2.3 #13 66.3 7.0 6.9 172.6 188.3 3.0 #14 67.3 8.0 10.8 179.1 155.4 3.2 #15 65.4 5.0 8.3 178.6 169.0 2.5 #16 68.0 14.0 15.0 186.0 120.1 4.2 #17 48.3 297.0 85.9 193.9 240.6 19.4 #18 74.3 34.0 25.5 189.0 255.1 6.6 #19 78.6 24.0 27.7 190.0 164.9 5.5 #20 65.0 1.0 0.0 189.0 187.0 1.1 #21 72.4 458.0 118.8 202.7 176.1 24.1 #22 55.0 109.0 37.1 197.0 119.5 11.8 #23 74.8 137.0 97.8 197.8 215.6 13.2 #24 65.0 4.0 4.4 193.5 95.0 2.3 #25 71.7 45.0 37.8 195.6 294.4 7.6 #26 83.6 521.0 88.7 215.4 209.0 25.8 #27 64.0 2.0 2.0 203.5 157.0 1.6 #28 73.6 311.0 135.9 206.3 297.0 19.9 #29 74.3 23.0 19.3 207.5 245.5 5.4 #30 69.2 13.0 13.6 211.9 191.5 4.1 #31 65.8 5.0 6.4 209.4 268.6 2.5 #32 80.5 1804.0 573.4 238.5 154.7 47.9 #33 70.4 523.0 93.7 225.3 242.3 25.8 #34 75.6 19.0 14.4 213.1 275.8 4.9 #35 67.0 4.0 3.6 214.5 307.5 2.3 #36 53.2 508.0 86.5 231.7 275.1 25.4 #37 66.4 583.0 193.7 234.2 309.7 27.2 #38 82.7 533.0 92.5 233.9 341.8 26.1 #39 65.8 5.0 7.0 220.8 374.2 2.5 #40 69.9 10.0 9.5 224.3 291.9 3.6 #41 73.3 32.0 22.6 229.5 324.2 6.4 #42 70.0 280.0 99.5 235.1 374.6 18.9 #43 79.0 40.0 25.9 234.8 358.3 7.1 #44 82.2 6.0 10.3 230.5 195.5 2.8 #45 81.5 2.0 2.0 232.0 200.5 1.6 #46 82.5 8.0 14.7 233.9 206.5 3.2 #47 66.5 12.0 11.3 236.3 403.9 3.9 #48 82.0 1.0 0.0 236.0 212.0 1.1 #49 76.2 9.0 17.1 238.7 218.0 3.4 #50 87.3 8.0 13.7 241.0 227.5 3.2 #51 84.0 5.0 7.8 242.0 236.0 2.5 #52 83.4 11.0 21.1 244.4 245.0 3.7 #53 64.6 8.0 10.3 246.8 384.6 3.2 #54 72.5 80.0 33.4 249.8 374.8 10.1 #55 83.9 13.0 24.0 247.8 258.0 4.1 #56 74.3 13.0 22.0 247.5 272.5 4.1 #57 66.3 2817.0 1152.0 289.8 233.5 59.9 #58 80.8 5.0 7.8 249.0 282.0 2.5 #59 82.0 4.0 5.9 250.0 287.5 2.3 #60 65.6 30.0 21.9 254.9 409.6 6.2 #61 93.4 7.0 11.8 251.0 294.0 3.0 #62 79.5 24.0 18.0 253.3 149.5 5.5 #63 79.5 1172.0 469.7 267.5 361.4 38.6 #64 69.9 8.0 7.7 253.9 162.3 3.2 #65 100.0 2.0 2.0 253.0 272.5 1.6 #66 75.2 176.0 52.7 262.3 171.5 15.0 #67 97.1 10.0 12.2 254.5 285.0 3.6 #68 95.9 20.0 24.5 255.2 347.7 5.0 #69 96.9 14.0 14.8 255.6 296.9 4.2 #70 99.7 3.0 3.9 257.0 359.0 2.0 #71 49.3 765.0 145.2 277.4 246.4 31.2 #72 100.0 1.0 0.0 267.0 258.0 1.1 #73 59.9 167.0 95.0 274.2 273.7 14.6 #74 83.8 89.0 117.3 271.3 309.4 10.6 #75 65.1 7.0 6.9 272.7 149.6 3.0 #76 69.6 41.0 24.2 278.6 144.2 7.2 #77 74.7 239.0 67.0 284.9 339.9 17.4 #78 67.8 59.0 49.0 284.1 160.1 8.7 #79 65.3 4.0 4.4 280.5 136.0 2.3 #80 77.3 236.0 83.1 301.3 359.1 17.3 #81 72.4 41.0 29.7 290.9 157.9 7.2 #82 65.4 8.0 10.0 289.5 377.0 3.2 #83 66.2 14.0 11.0 292.0 372.5 4.2 #84 65.6 8.0 7.7 297.8 392.9 3.2 #85 64.3 4.0 5.9 303.0 149.5 2.3 #86 67.3 4.0 5.9 312.0 281.5 2.3 #87 86.0 5.0 7.8 312.0 290.0 2.5 #88 80.2 280.0 164.2 332.8 320.3 18.9 #89 100.0 1.0 0.0 314.0 191.0 1.1 #90 64.9 7.0 8.6 317.3 396.7 3.0 #91 99.0 2.0 2.0 318.0 300.5 1.6 #92 98.7 7.0 11.4 319.7 305.9 3.0 #93 85.3 4.0 5.9 320.0 204.5 2.3 #94 81.8 4.0 5.9 321.0 210.5 2.3 #95 81.0 5.0 7.8 322.0 216.0 2.5 #96 64.8 9.0 11.0 324.2 363.0 3.4 #97 79.3 3.0 3.9 323.0 222.0 2.0 #98 97.5 14.0 13.4 325.1 314.8 4.2 #99 83.0 3.0 3.9 324.0 227.0 2.0 #100 83.3 4.0 5.9 325.0 232.5 2.3 #101 80.0 5.0 7.8 326.0 238.0 2.5 #102 81.7 29.0 55.9 328.3 256.0 6.1 #103 64.0 4.0 4.4 328.5 363.0 2.3 #104 76.1 20.0 37.2 330.0 288.5 5.0 #105 75.5 2.0 2.0 331.0 300.5 1.6 #106 67.1 10.0 9.9 340.1 343.0 3.6