MATLAB Examples

%------------------------------------------------------------------------------------------------
% 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