How to indentify a broken pill package?
1 次查看(过去 30 天)
显示 更早的评论
Hello, I am writing a program, which should read in an image of a pill package and count how many pills are left over in the package. The code uses a combination of image processing functions including the canny transform into a gradient image and the hough transform for circles. Here are some images to understand what I am talking about (The numbers in the green circles show the order of the circles' metric values, from high to low)



Now I am trying to write some code, which can identify the difference between pills, which are already used and those, which are still left over in the package. Unfortunately the metric value does not give a good indication if the pill pocket is broken or not. I also calculated a ratio where I divided the number of gradient pixels by the total area of the circle. Unfortunately this value is also very unreliable for this problem.
Do you have any other ideas how to identify which of these pill pockets are broken? If it helps, the full pill package is always analyzed first and then, some informations are saved and can be loaded into the code again as a possible comparision. However, the images are never 100% the same, since they are all done by hand within an app.
I'm thankful for any tips, ideas and help!
采纳的回答
Image Analyst
2020-10-9
I'd probably start with just having a mask/template with ROIs over each pill pocket and take the standard deviation of the image intensities. If the standard deviation is high, then it's been used. If it's low, then the intensity if pretty uniform and the pill has not been pushed out.
Or you could try a deep learning app.
17 个评论
Jannis Holtkoetter
2020-10-9
Thank you for this quick reply. I already have a mask, which includes the information about the ROIs and other regionprops information. It looks like this:

I tried the calculation of the standard deviation in different image channels such as the value, saturation channel and even in a grayscale image of the original, however there is unfortunately no peak for the circles of the broken pills.
Do you have any other ideas? Or any advice how to implement a deep learning app into this. That would be a completely new topic for me.
Image Analyst
2020-10-9
I did not say to do a histogram on the Canny edge image. Of course that would be not very helpful since that image is binary. You need to do it on the original gray scale image. Perhaps something like (untested)
[labeledImage, numPills] = bwlabel(mask);
props = regionprops(mask, 'Centroid');
% Get all centroids.
xy = vertcat(props.Centroid)
stdDevs = zeros(numPills, 1);
for k = 1 : numPills
% Get a mask of just this one pill.
thisPill = ismember(labeledImage, k);
imshow(thisPill);
drawnow;
stdDevs(k) = std(grayImage(thisPill))
end
% If it's standard deviation is higher than some number, say 10, then put a red X over it.
imshow(grayImage, []);
hold on;
for k = 1 : numPills
if stdDevs(k) > 10
% Plot red X over it.
plot(xy(k, 1), xy(k, 2), 'rx', 'MarkerSize', 30, 'LineWidth', 2);
else
% Plot green spot over it.
plot(xy(k, 1), xy(k, 2), 'g.', 'MarkerSize', 30, 'LineWidth', 2);
end
end
Jannis Holtkoetter
2020-10-9
编辑:Jannis Holtkoetter
2020-10-10
Thanks again for your help! I tried that with the image I already attached previously.
I have a questions for the usage of the grayscale image. I was using my RGB image and used the function "grayImage = rgb2gray(RGBImage)" in order to transform it form RGB into grayscale. Then the std function gave me the error, that it can only work with a single or double data type. So, I changed the grayscale image into double with the im2double function. is there any flaw in that?
Otherwise I show you the result of your offered code. I changed the std limit to 0.2, since the std Values where somewhere between 0.17 and 0.23. Unfortunately the highest value are not the broken ones.


So, are there any other ideas?
Image Analyst
2020-10-9
编辑:Image Analyst
2020-10-9
I would have thought the SD's would be more different than that. But the broken ones are above 0.2 while the others are less. You just don't understand how the masks are ordered. It goes top down then left to right. In other words, whichever one is leftmost is #1, regardless if it's in the top row or not. For example we know that blobs 5, 6, and 7 are in the right column, but I see that 5 and 6 have 0.2 so that means they're in row 1 and 4 and are very slightly positioned to the left of the blob in row 2, column 2.
I suggest you add the region number to the overlay to see how they're ordered:
for k = 1 : numPills
if stdDevs(k) > 10
% Plot red X over it.
plot(xy(k, 1), xy(k, 2), 'rx', 'MarkerSize', 30, 'LineWidth', 2);
else
% Plot green spot over it.
plot(xy(k, 1), xy(k, 2), 'g.', 'MarkerSize', 30, 'LineWidth', 2);
end
txt = sprintf(' %d', k);
text(xy(k, 1), xy(k, 2), txt, 'Color', 'b', 'FontSize', 20, 'FontWeight', 'bold');
end
If you want to you could sort them based on rows and columns. Use kmeans() to categorize the x centroids into 2 classes, and the y into 4 classes. Then use sort to resort them if you want. However that really had no bearing on the count of broken or unbroken.
To explain labeling, run this code:
% Illustrates how regions labels (i.e. numbers) are assigned to blobs randomly placed in an image.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 22;
%-------------------------------------------------------------------------------------------------------------------------------------
% First create image with a bunch of regions in it.
% Define the number of regions
numRegions = 20;
% Define a max region size for the circles.
maxRadius = 30;
minRadius = 6;
% Define the image size
rows = 600;
columns = 800;
% Create a binary image
binaryImage = false(rows, columns);
% Get a list of all coordinates in the image.
% First create the image.
imageSizeX = 640;
imageSizeY = 480;
[columnsInImage, rowsInImage] = meshgrid(1:columns, 1:rows);
% Start adding circles of random size.
for k = 1 : numRegions
% Create a logical image of a circle with specified
% diameter, center, and image size.
% Next create the circle in the image.
radius = minRadius + (maxRadius - minRadius) * rand(1);
centerX = maxRadius + (columns - 2 * maxRadius) * rand(1);
centerY = maxRadius + (rows - 2 * maxRadius) * rand(1);
fprintf('Pasting circle #%d of radius %.1f centered at (%.1f, %.1f)\n', k, radius, centerX, centerY);
circlePixels = (rowsInImage - centerY).^2 ...
+ (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Paste this circle onto binaryImage
binaryImage = binaryImage | circlePixels;
% imshow(binaryImage);
% drawnow;
end
%-------------------------------------------------------------------------------------------------------------------------------------
% Next let's do connected components labeling of the regions.
[labeledImage, numRegions] = bwlabel(binaryImage); % bwlabel() is where it actually does the connected components labeling.
% Let's assign each blob a different color to visually show the user the distinct blobs,
% just in case some overlapped which would cause them to merge into one single blob instead of two.
coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels
% coloredLabels is an RGB image. We could have applied a colormap instead (but only with R2014b and later)
imshow(coloredLabels);
axis('on', 'image');
caption = sprintf('Binary image of %d distinct, separate regions', numRegions);
title(caption, 'FontSize', 15);
%-------------------------------------------------------------------------------------------------------------------------------------
% Now find out their pixel locations so we can put up a line there along the left edge, which is how it determines the label number.
props = regionprops(binaryImage, 'PixelList', 'Centroid');
for k = 1 : numRegions
thesePixels = props(k).PixelList;
[leftColumn(k), index] = min(thesePixels(:, 1)); % Min x value.
% Get the y value for the left edge of the blob.
yt = thesePixels(index, 2);
% Put up a line there
xline(leftColumn(k), 'Color', 'w');
% Put a numerical label at the centroid.
xt = leftColumn(k); %props(k).Centroid(1);
txt = sprintf(' %d', k);
text(xt, yt, txt, 'Color', 'k', 'FontWeight', 'bold', 'FontSize', 12, 'VerticalAlignment', 'middle');
end
% Show all the left columns in the command window so you'll see they are sorted from left to right.
leftColumn
% Maximize the window.
g = gcf;
g.WindowState = 'maximized';

Now scroll over to the right and look at the special case where blobs 17 and 18 both started at column 683. Because the blue blob (17) is higher than the pink blob (18), it was assigned the number 17 and the lower pink blob was assigned number 18. But in general you can see that blobs are numbered from left to right. And your mask is not a perfect rectangular array so it's quite possible some blobs have numbers different than what you might think at first.
Jannis Holtkoetter
2020-10-10
Hello Image Analyst, yes I was aware of the ordering of the blobs in my mask image. I should have made it clearer in my post. And I applied your code to my image and it confirms the ordering from left to right and then from top to bottom:

However, this doesn't change my problem, since there is no special behaviour of the standard deviation in the broken pills (no. 4 & 6).
I was also looking into texture filters based on the local entropy, local stadard deviation and local range. I found the information about it here: Texture Segmentation Using Texture Filters (https://www.mathworks.com/help/images/texture-segmentation-using-texture-filters.html).
The application of these filters to my image gives me this result:

Do you have any ideas how I could apply this maybe in combination with the information I already have about the location of the pills in order to get a better solution?
Image Analyst
2020-10-10
You could use those images as the gray scale image in regionprops and compute the mean intensity
stdImage = stdfilt(grayImage, ones(9));
props = regionprops(mask, stdImage, 'MeanIntensity');
Since you're not sharing the mask or code how you got it, with us, you'll have to do it yourself.
Jannis Holtkoetter
2020-10-10
Thanks again!
I got the mask by using the hough transform for circles and then using different information about the pills (all same radius, rough X-axis-location, etc. ) in order to filter out all possible circles which were detected in the Hough Transform but are not pills. This works for the main part of my images, however I can't share any code, since it is quite long and there are different bits of code involved in the complete code.
However, I tried the calculation of the mean intensity, but it seems to give me unclear results for some images. I mean you can clearly see, that the mean intensity is generally higher in the broken pills, however the values are very close to each other:

After doing several tests, the limit which differentiates a non broken pill from a broken pill seems to be somewhere between 0.13 and 0.14. However as you can see in this example, sometimes there can be a little dent or some light reflection which bumps up the mean intensity.
Do you have any idea how to improve the differentiation between the values, so that it is a little bit more clear which ones are the real broken ones?
Image Analyst
2020-10-11
You can combine several metrics, like one based on texture and one based on histogram and one based on something else, like perhaps how radially symmetric the region is. Then make a decision based on the whole collection of measurements you made, not just one.
Jannis Holtkoetter
2020-10-11
Yes, I actually had that in mind to combine several metrics. So far, I just couldn't figure out which ones.
Could you maybe how you would proceed based on the histogram and based on the radially symmetry of the region? I haven't used those approaches so far.
Image Analyst
2020-10-11
Attached are demos where I compute the average radial profile. Adapt as needed.
Jannis Holtkoetter
2020-10-11
Thank you very much for that! I don't quite understand, what is shown in the Euclidean distance transform and what I can do with the average radial profile then. Should I compute the profile for each of my ROIs? And if so, how do I do that?
Image Analyst
2020-10-11
Yes, for each blob you can compute the average radial profile from the center using the Euclidean distance transform.
Jannis Holtkoetter
2020-10-11
Okay, and how do I do that? Do I have to include the code from 'average_radial_profile.m' into a loop where I use an overlay of the grayscale image with one blob of the mask for each iteration? That seems like a pretty cumbersome process.
Do you know a better plan how to compute the profiles only for the ROIs?
Image Analyst
2020-10-11
In your function that analyzes a single image, AnalyzeSingleImage(), you could have a loop over all blob regions. And in that loop you could call a function ComputeAverageRadialProfile() where it computes and returns a 1-D list of the intensity values from the center out to the edge of the blob. Not that cumbersome. Most of my programs have between 3000 and 7000 lines of code, and this is nowhere near that complicated/cumbersome.
Jannis Holtkoetter
2020-10-12
Yes, I am sure there are many way more complicated things out there. However, I still can't really connect, how I can implement the idea. Could you maybe give me some advice?
Should I use an image like this, where I applied the mask onto the local standard deviation filter image?

And then how do I calculate the radial profile for each pill individually with the center of the pill as the starting point for the calculation? In bwdist function I couldn't figure out how to change the starting (center) point.
Or is it maybe smarter to use an image, where I cropped the image which each BoundingBox of the pills, such as following (this would be for the upper left pill from the image):

And then I guess I also have to find a proper way, how to automatically binarize this image, so that the threshold always suits the purpose. How would you do that?
I hope you can help me with some of my questions, since I really want to implement the idea, but I don't know how.
Image Analyst
2020-10-12
Yes, I'd just start with the cropped bounding box image. And of course you'll need the cropped circular mask image also so you can get the EDT of that so you know what radius the pixel's intensity is measured at.
Jannis Holtkoetter
2020-10-12
I applied your recommendation of using the crop and feeding these images into a loop with the averageRadialProfile function. I have to say, that I used the grayscale image based on the local standard deviation filter (values between 0 and 1). Is that a problem and should I have used just a grayscale version of the original jpg image (values between 0 and 255? Does it change anything?
The caluclation is being done with the code, which you provided in the file "average_radial_profile.m", as shown in the following. I was looking into the documentation of the different functions, but had problems to fully understand what the averageRadialProfile actually describes...
edtImage = bwdist(binaryImage);
maxDistance = ceil(max(edtImage(:)));
% Allocate an array for the profile
profileSums = zeros(1, maxDistance);
profileCounts = zeros(1, maxDistance);
% Scan the original image getting gray level, and scan edtImage getting distance.
% Then add those values to the profile.
[rows, columns] = size(edtImage);
for column = 1 : columns
for row = 1 : rows
thisDistance = round(edtImage(row, column));
if thisDistance <= 0
continue;
end
profileSums(thisDistance) = profileSums(thisDistance) + double(grayImage(row, column));
profileCounts(thisDistance) = profileCounts(thisDistance) + 1;
end
end
% Divide the sums by the counts at each distance to get the average profile
averageRadialProfile = profileSums ./ profileCounts;
In the loop I put each averageRadialProfile array into one cell array which now contains the averageRadialProfile for each pill. If I plot it for the example above, it gives me this result (the plots with the red circle are broken pills). How would can I analyze these results in a quantified way, so that I can atleast draw some conclusion between the average radial profile and a broken pill?

Again, I would be very thankful for any help!
更多回答(1 个)
Jannis Holtkoetter
2020-10-22
I just accepted the answer with regard to the fact, that I will be able to use the mean intensity of the local standard deviation image for each pill in combination with the mean intensity of a canny edge image for each pill to define if a pill can be marked as used or not. It is definetly not working perfect but it works for a lot of specific looking blisters.
I still do not really know, how to interpret the radial profile of each pills and therefore I have no use for it, however the provided code worked for me and I could compute each radial profile in my example. I'm still open for any explanation or suggestion for implementation.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 3-D Volumetric Image Processing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom(English)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)