MATLAB Answers

i am currently trying to detect sickle cell using this code which i got from the MATLAB answer. The code run well, but it didn't detect the exact sickle cell. maybe because of the value parameter. where should i adjust the parameter value?

11 views (last 30 days)
haniz azwa
haniz azwa on 2 Oct 2017
Commented: John D'Errico on 2 Oct 2017
% Startup code.
tic; % Start timer.
clc; % Clear command window.
clearvars; % Get rid of variables from prior run of this m-file.
fprintf('Running SickleDemo.m...\n'); % Message sent to command window.
workspace; % Make sure the workspace panel with all the variables is showing.
imtool close all; % Close all imtool figures.
format long g;
format compact;
captionFontSize = 14;
[fn, pn] = uigetfile('sickle cell.bmp','select bitmap file');
fullFileName = strcat(pn,fn);
if ~exist(fullFileName, 'file')
% It doesn't exist in the current folder.
% Look on the search path.
if ~exist(baseFileName, 'file')
% It doesn't exist on the search path either.
% Alert user that we can't find the image.
warningMessage = sprintf('Error: the input image file\n%s\nwas not found.\nClick OK to exit the demo.', fullFileName);
fprintf(1, 'Finished running SickleDemo.m.\n');
% Found it on the search path. Construct the file name.
fullFileName = baseFileName; % Note: don't prepend the folder.
% If we get here, we should have found the image file.
originalImage = imread(fullFileName);
% Check to make sure that it is grayscale, just in case the user substituted their own image.
[rows, columns, numberOfColorChannels] = size(originalImage);
if numberOfColorChannels > 1
promptMessage = sprintf('Your image file has %d color channels.\nThis program was designed for grayscale images.\nDo you want me to convert it to grayscale for you so you can continue?', numberOfColorChannels);
button = questdlg(promptMessage, 'Continue', 'Convert and Continue', 'Cancel', 'Convert and Continue');
if strcmp(button, 'Cancel')
fprintf(1, 'Finished running SickleDemo.m.\n');
% Do the conversion using standard book formula
hsv = rgb2hsv(im2double(originalImage));
% Display the grayscale image.
mask = hsv(:,:,2) > 0.2;
mask_remove = bwareaopen(mask, 100);
mask_fill = imfill(mask_remove, 'holes');
se = strel('square', 7);
mask_final = imdilate(mask_fill, se);
mask_final = repmat(mask_final, [1 1 3]);
originalImage= originalImage .* uint8(~mask_final);
originalImage(mask_final) = 255;
originalImage = rgb2gray(originalImage);
% Maximize the figure window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Force it to display RIGHT NOW (otherwise it might not display until it's all done, unless you've stopped at a breakpoint.)
caption = sprintf('Grey image');
title(caption, 'FontSize', captionFontSize);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
% its histogram and display it.
[pixelCount, grayLevels] = imhist(originalImage);
title('Histogram of original image', 'FontSize', captionFontSize);
xlim([0 grayLevels(end)]); % Scale x axis manually.
grid on;
% Threshold the image to get a binary image (only 0's and 1's) of class "logical."
% using a logical operation.
thresholdValue = 200;
%binaryImage = originalImage > thresholdValue; % Bright objects will be chosen if you use >.
% ========== IMPORTANT OPTION ============================================================
% Use < if you want to find dark objects instead of bright objects.
binaryImage = originalImage < thresholdValue; % Dark objects will be chosen if you use <.
% Do a "hole fill" to get rid of any background pixels or "holes" inside the blobs.
binaryImage = imfill(binaryImage, 'holes');
binaryImage = imclearborder (binaryImage);
binaryImage = bwareaopen(binaryImage, 100);
% Show the threshold as a vertical red bar on the histogram.
hold on;
maxYValue = ylim;
line([thresholdValue, thresholdValue], maxYValue, 'Color', 'r');
% 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]);
% Display the binary image.
title('Binary Image, after thresholding', 'FontSize', captionFontSize);
% Identify individual blobs by seeing which pixels are connected to each other.
% Each group of connected pixels will be given a label, a number, to identify it and distinguish it from the other blobs.
% Do connected components labeling with either bwlabel() or bwconncomp().
labeledImage = bwlabel(binaryImage, 8); % Label each blob so we can make measurements of it
% labeledImage is an integer-valued image where all pixels in the blobs have values of 1, or 2, or 3, or ... etc.
imshow(labeledImage, []); % Show the gray scale image.
title('Labeled Image, from bwlabel()', 'FontSize', captionFontSize);
% Get all the blob properties.
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 coins on the original grayscale image using the coordinates returned by bwboundaries.
title('Outlines, from bwboundaries()', 'FontSize', captionFontSize);
axis image; % Make sure image is not artificially stretched because of screen's aspect ratio.
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
thisBoundary = boundaries{k};
plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
hold off;
textFontSize = 10; % 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.
% directly into regionprops.
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 one at a time
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', textFontSize, 'FontWeight', 'Bold');
%another way to get centroids.
% We can get the centroids of ALL the blobs into 2 arrays,
% one for the centroid x values and one for the centroid y values.
allBlobCentroids = [blobMeasurements.Centroid];
centroidsX = allBlobCentroids(1:2:end-1);
centroidsY = allBlobCentroids(2:2:end);
% Put the labels on the rgb labeled image also.
%subplot(2, 2, 3);
for k = 1 : numberOfBlobs % Loop through all blobs.
text(centroidsX(k) + labelShiftX, centroidsY(k), num2str(k), 'FontSize', textFontSize, 'FontWeight', 'Bold');
% Now I'll demonstrate how to select certain blobs based using the ismember() function.
% Let's say that we wanted to find only those blobs
% with an intensity between 150 and 220 and an area less than 2000 pixels.
% This would give us the three brightest dimes (the smaller coin type).
allBlobIntensities = [blobMeasurements.MeanIntensity];
allBlobAreas = [blobMeasurements.Area];
% Get a list of the blobs that meet our criteria and we need to keep.
% These will be logical indices - lists of true or false depending on whether the feature meets the criteria or not.
% for example [1, 0, 0, 1, 1, 0, 1, .....]. Elements 1, 4, 5, 7, ... are true, others are false.
allowableIntensityIndexes = (allBlobIntensities > 150) & (allBlobIntensities < 220);
allowableAreaIndexes = allBlobAreas < 2000; % Take the small objects.
% Now let's get actual indexes, rather than logical indexes, of the features that meet the criteria.
% for example [1, 4, 5, 7, .....] to continue using the example from above.
keeperIndexes = find(allowableIntensityIndexes & allowableAreaIndexes);
% 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. Result will be an image - the same as labeledImage but with only the blobs listed in keeperIndexes in it.
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.
% Plot the centroids in the original image in the upper left.
% Dimes will have a red cross, nickels will have a blue X.
message = sprintf('Now I will plot the centroids over the original image.');
reply = questdlg(message, 'Plot Centroids?', 'OK', 'Cancel', 'Cancel');
% Note: reply will = '' for Upper right X, 'OK' for OK, and 'Cancel' for Cancel.
if strcmpi(reply, 'Cancel')
%subplot(2, 2, 1);
RBC = 0;
Sickle = 0;
hold on; % Don't blow away image.
for k = 1 : numberOfBlobs % Loop through all keeper blobs.
% Identify if blob #k is a sickle or nickel.
itsADime = allBlobAreas(k) < 240; % sickle are small.
if itsADime
% Plot dimes with a green +.
plot(centroidsX(k), centroidsY(k), 'g+', 'MarkerSize', 10, 'LineWidth', 2);
Sickle = Sickle + 1;
% Plot dimes with a yellow x.
plot(centroidsX(k), centroidsY(k), 'yx', 'MarkerSize', 10, 'LineWidth', 2);
RBC = RBC + 1;
disp(['RBC = ' num2str(RBC)]) ;
disp(['Sickle = ' num2str(Sickle)]) ;
h = msgbox(cat(1, {'Cells counting succeeded.'}, {'No. of Red Blood cells :'}, num2str(RBC), {'No. of Sickle Cell :'}, num2str(Sickle)), 'Success');
% Now use the keeper blobs 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.
% Now let's get the nickels (the larger coin type).
keeperIndexes = find(allBlobAreas > 2000); % Take the larger objects.
% Note how we use ismember to select the blobs that meet our criteria.
nickelBinaryImage = ismember(labeledImage, keeperIndexes);
% Let's get the nickels from the original grayscale image, with the other non-nickel pixels blackened.
% In other words, we will create a "masked" image.
maskedImageNickel = originalImage; % Simply a copy at first.
maskedImageNickel(~nickelBinaryImage) = 0; % Set all non-nickel pixels to zero.
elapsedTime = toc;
% 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; % Create a new figure window.
% Maximize the figure window.
set(gcf, 'Units','Normalized','OuterPosition',[0 0 1 1]);
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);
% Determine if it's a dime (small) or a nickel (large coin).
if blobMeasurements(k).Area > 230
coinType = 'Red Blood Cell';
coinType = 'Sickle Cell';
% Display the image with informative caption.
subplot(3, 4, k);
caption = sprintf('Coin #%d is a %s.\nDiameter = %.1f pixels\nArea = %d pixels', ...
k, coinType, blobECD(k), blobMeasurements(k).Area);
title(caption, 'FontSize', textFontSize);


Rik on 2 Oct 2017
This is the worst of both worlds: we don't have all relevant files to run this code, but it is a wall of text. Can you condense this to a MWE (minimal working example)? Also, you mention that you got this code from somewhere else. Where did you get this? That could give an important clue.
John D'Errico
John D'Errico on 2 Oct 2017
Contact the person you got the code from. It is virtually impossible to answer this question for someone else who did not write the code.

Sign in to comment.

Answers (0)

Translated by