How to extract color parameters of each rice grain and then how should I compare them?

I have an image named 50, I have segmented the iamge and obtained label of each object (rice grain). This image contains rice of different varieties and I want to separate them based on their colour properties because they are similar in shape and size. Can anyone explain me which parmateters will help me to distinguish them.

[KSSV]

You can go by area, perimeter. ???

[Ankit Rathore]

I would have, but the grains are of same size and another problem is that both varities contains some grains of same size?

Answer 1

[DGM]

Consider what information is available. The exposure doesn't seem too bad. There are some grains which seem maybe a bit more colorful than others, but it's hard to tell what's meaningful.

First, strip out the background and adjust the levels so that the grain color is more pronounced.

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/890025/50.jpg');
A = rot90(A); % i'm just rotating this so it fits better in the web-view

A = imadjust(A,[0.7 1]);
mask = bwareaopen(rgb2gray(A)>128,1000);
A(~mask(:,:,[1 1 1])) = 0;

imshow(A)

Is that enough information to identify which are which? I'm not sure.

Select only the pixels representing the rice grains. Out of this entire population of pixels, is there any salient group information?

mask = imerode(mask,strel('disk',5)); % erode the mask to ignore edge shading
Am = reshape(A(mask(:,:,[1 1 1])),[],1,3);

Amlab = rgb2lab(Am);
for c = 1:3
    subplot(3,1,c)
    histogram(Amlab(:,:,c))
end

Maybe? It's hard to tell. Of course, that's the entire population of pixels. What if we consider the average color of each grain?

Alab = rgb2lab(A);

[L N] = bwlabel(mask);
graincolor = zeros(N,3);
for b = 1:N
    grainmk = L==b;
    Agrain = Alab(grainmk(:,:,[1 1 1]));
    graincolor(b,:) = mean(reshape(Agrain,[],3),1);
end
graincolor % average color [L a b] for each grain

graincolor = 50×3

69.0426 -9.3356 27.3126 80.1833 -7.5283 30.9029 76.8120 -8.9078 25.2107 78.9900 -7.3287 23.0446 75.2678 -6.4827 22.0959 77.4466 -8.9433 26.8382 82.7059 -8.4911 24.3281 88.5259 -6.4093 18.3634 85.4462 -7.5029 22.2742 87.3065 -9.6987 29.2395

figure
for c = 1:3
    subplot(3,1,c)
    histogram(graincolor(:,c))
end

Now it looks even less like there's a discernable pattern. What if we just look at where those color points are and try to find a pattern by eye?

figure
plot3(graincolor(:,2),graincolor(:,3),graincolor(:,1),'b.')
grid on

Well, I'm not sure if that really helps either. It's hard to say if I'm getting anywhere since I'm not sure which are which to begin with. Let's just say that there's a group defined by L>90. Let's just select that.

figure
grainmask = graincolor(:,1)>90; 

S = regionprops(mask,'centroid');
C = vertcat(S.Centroid);
imshow(A); hold on
for b = 1:N
    if grainmask(b)
        plot(C(b,1),C(b,2),'bo','markersize',25);
    else
        plot(C(b,1),C(b,2),'ro','markersize',25);
    end
end

Are those the right ones?

[Image Analyst]

My guess is that he forgot to take a photo with both types of rice on it, or forgot to attach the photo that had only the other type of rice on it.

Not sure what the real world context is. I'm guessing there is none. No one in the real world is going to sort rice by mixing different types together and then tossing them on a black paper to be sorted. Seems like something you'd only do in a class assignment. By the way, for another way to do it, you could use regionprops to measure the color in each blob.

Alab = rgb2lab(A);
% Get L mean intensities
propsL = regionprops(mask, Alab(:,:,1), 'MeanIntensity')
% Get A mean intensities
propsA = regionprops(mask, Alab(:,:,2), 'MeanIntensity')
% Get B mean intensities
propsB = regionprops(mask, Alab(:,:,3), 'MeanIntensity')
% Get all the colors in an N by 3 array.  Average color [L a b] for each grain
graincolor = [[propsL.MeanIntensity]; [propsA.MeanIntensity]; [propsB.MeanIntensity]]';

[DGM]

I didn't even think about that.

I swear. I need to keep a cached copy of the regionprops webdocs. The help synopsis is too terse for its features, and loading webdocs is slow enough that it becomes a hassle I avoid.

[yanqi liu]

may be use KMeans to make color cluster

[Image Analyst]

kmeans is not useful or robust in this situation. The algorithm should be able to detect the percentage of each type of rice regardless if it's 100% of type 1 and 0% of type 2, or 0% of type 1 and 100% of type 2, or any amount in between. If you say k=2 but it's like >95% of one type, kmeans will force it to find two classes and it will probably say it's like 50% of type 1 and 50% of type 2, when that is not what you have. It's better to use fixed thresholds on some feature, or else some kind of model that predicts each grain of rice by itself and is not influenced by how many of each type there are (for example measure 20 things about each rice grain and train a model using the Regression Learner app).

[Ankit Rathore]

I have not forgot to take photo of both the rice . And coming to the use case of this . Its actually a big problem to separate the grains while sampling and studying on it further. I would have explained usecase but this is not the plaform for that. For Better understaning I am attaching some new image where I have also marked the rice grain of different type. Image "Main" contains original mixture of rice grains while Image "Result" contains the grains marked of other type. I want the to separate out grains which I have marked by circle. Thanks in Advance

[Image Analyst]

There is no significant color difference between the ones you circled and the ones you didn't. They all range from 240 to 255 with many of the colors nearly pegged aat 255 so we don't know what color they were really mean to be because they were clipped.

rgbImage = imread('result.jpg');
rgbImage = rescale(double(rgbImage), 0, 1);
rgbImage2 = imadjust(rgbImage, [240, 255]/255, [0, 255]/255);
imshow(rgbImage2);
impixelinfo;

Above is an image scaled so that the range is changed from 240-255 to 0 to 255. It includes both ones you circled and didn't, so you can see that you can't really tell which is which.

[DGM]

I'm not really sure I can identify the difference either. Like I said, some of the grains look lighter or darker than others, but I'm not sure if that's charactersistic of the variety, or something else (moisture content, external abrasion, orientation WRT illuminator). Some of the grains seem subtly shaped different, but I'm not sure what that means either.

Regarding lightness and shape, I don't see a pattern in what's been annotated -- at least not offhand.

Lemme see if I can see anything after processing.

[DGM]

This is what I got:

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/893335/Main.jpg');
A = imresize(A,0.25); % i'm just going to work with a smaller copy for testing
A = imadjust(A,[0.7 1]);

mask = bwareaopen(rgb2gray(A)>128,1000);
A(~mask(:,:,[1 1 1])) = 0;

% get mean color of each grain
Alab = rgb2lab(A);

[L N] = bwlabel(mask);
graincolor = zeros(N,3);
for b = 1:N
    grainmk = L==b;
    Agrain = Alab(grainmk(:,:,[1 1 1]));
    graincolor(b,:) = mean(reshape(Agrain,[],3),1);
end

% plot grain color points
r51marked = [5 6 24 30]; % indices for marked grains in rice51.jpg
mk = false(N,1);
mk(r51marked) = true;

plot3(graincolor(mk,2),graincolor(mk,3),graincolor(mk,1),'r.'); hold on
plot3(graincolor(~mk,2),graincolor(~mk,3),graincolor(~mk,1),'b.')
grid on

figure

% make a selection and plot
grainmask = graincolor(:,1)>93 ...
          & graincolor(:,2)<-6 ...
          & graincolor(:,3)>19 ...
          & sum(graincolor(:,2:3),2)>12; 

S = regionprops(mask,'centroid');
C = vertcat(S.Centroid);
imshow(A); hold on
for b = 1:N
    if grainmask(b)
        plot(C(b,1),C(b,2),'ro','markersize',15);
    else
        plot(C(b,1),C(b,2),'bo','markersize',15);
    end
end
hold off

Note that while I was able to find those four grains, the mask needed to be fairly elaborate to exclude nearby color points. With a larger sample size, it may prove impractical to find a geometry which can reliably make this exclusion.

Note also that the mask thresholds depend on the exposure and the input levels given to imadjust().

[Ankit Rathore]

Thank You @DGM you made it. It helped me a lot. Should I integrate it with HSV YCrCb to increase the efficiency of it?

[DGM]

You can try to plot the color points in HSV or anything you want and see if they're more easily separable, but in my experiments, I didn't see much difference. If you're after computational efficiency, using YCbCr might speed up that one line by quite a bit, but it only contribues to a small amount of the total time anyway.

Separation might be a bit easier in LCH simply due to the axis alignment with the apparent group geometry, but the benefit is uncertain at this point and I didn't want to overcomplicate the example.

Further optimizing the mask really warrants a larger sample size so that the shape and variability of the ROI is more clearly known.

[Image Analyst]

Ycbcr takes longer:

rgbImage = imread('peppers.png');
tic
for k = 1 : 1000
    hsvImage = rgb2hsv(rgbImage);
end
hsvTime = toc

hsvTime = 0.6810

tic
for k = 1 : 1000
    hsvImage = rgb2ycbcr(rgbImage);
end
ycbcrTime = toc

ycbcrTime = 3.2302

Here is the gamut in HSV space:

and here it is in YCbCr space:

You can see that the gamut shapes are extremely similar. Also, we don't know where in that gamut your circled grains lie. Like @DGM and I have been saying, your circled grains are somewhere in the midst of all that. It's not like they are all at the top, middle, or outside edge where they can easily be thresholded.

@Ankit Rathore I think it might benefit you if you learned something about color science first. I don't mean to offend you. It's a very very complicated subject that takes years to learn and it's very confusing even for experts (I know, I've asked some of the top people in the world in person) because it mixes physics with physiology and psychology. All this stuff about different color spaces and delta E, etc. is totally baffling for beginners (even me when I started decades ago). This might be one place to start:

https://www.hunterlab.com/en/color-measurement-learning/

[DGM]

I meant that rgb2ycbcr() would be quicker than the demonstrated method, which was rgb2lab(), which it is about 10x as fast for me (though that depends). I didn't think that HSV really offered much margin for error in a rectangular mask region. To be fair, the demonstrated method didn't either, hence the diagonal face.

One can of course extend the color point viewing routine to provide more information about the margin around a simple box selection in various models. As before, this plots the target grains in red, and the unselected grains in blue. Any false positives that are unintentionally selected with the target grains are in green. For a scale factor of 1, this creates the smallest axis-aligned box selection that can include the target points. By simply scaling the box geometry, one can get an idea how much wiggle room there is. I've included some numbers for a handful of tests.

EDIT: fixed box scaling factor, updated

A = imread('https://www.mathworks.com/matlabcentral/answers/uploaded_files/893335/Main.jpg');
A = imresize(A,0.25);
A = imadjust(A,[0.7 1]);

mask = bwareaopen(rgb2gray(A)>128,1000);
A(~mask(:,:,[1 1 1])) = 0;

% get mean color of each grain
Alab = rgb2hsv(A); % yeah yeah i didn't change the variable name

[L N] = bwlabel(mask);
graincolor = zeros(N,3);
for b = 1:N
    grainmk = L==b;
    Agrain = Alab(grainmk(:,:,[1 1 1]));
    graincolor(b,:) = mean(reshape(Agrain,[],3),1);
end

% plot grain color points
% find roi grains
r51marked = [5 6 24 30]; % indices for marked grains in rice51
mk1 = false(N,1);
mk1(r51marked) = true;

% maximum allowable box size
% NATIVE TOOLS
% +8% for lab (rgb2lab(A))
% +13% for hsv (rgb2hsv(A))
% +44% for ycbcr (rgb2ycbcr(A))
% +47% for yiq (rgb2ntsc(A))
% MIMT TOOLS
% +56% for hsl (rgb2hsl(A))
% +16% for hsi (rgb2hsi(A))
% +45% for lchab (rgb2lch(A,'lab'))
% +32% for lchuv (rgb2lch(A,'luv'))
% +45% for lchsr (rgb2lch(A,'srlab')) (srlab2)
% +46% for lchok (rgb2lch(A,'oklab')) (oklab)
% +49% for ych (rgb2lch(A,'ypbpr')) (polar ypbpr)
% -26% for hsy (rgb2hsy(A)) (chroma-normalized ypbpr)

% find false matches, isolate masks
boxscale = 1.15; % scale selection box
roibox = graincolor(mk1,:);
roibox = [min(roibox,[],1); max(roibox,[],1)];
bwdelta = abs(diff(roibox,1,1))*(boxscale-1);
roibox = roibox + [-1; 1].*bwdelta/2;
mk3 = all(graincolor>=roibox(1,:) & graincolor<=roibox(2,:),2);
mk2 = ~mk1 & ~mk3;
mk3 = mk3 & ~mk1;
nnz(mk3) % number of false positives

ans = 1

% plot color points
ms = 20;
scatter3(graincolor(mk1,2),graincolor(mk1,3),graincolor(mk1,1),ms,'r','filled'); hold on
scatter3(graincolor(mk2,2),graincolor(mk2,3),graincolor(mk2,1),ms,'b','filled')
ph = scatter3(graincolor(mk3,2),graincolor(mk3,3),graincolor(mk3,1),ms,'g','filled');
ph.MarkerEdgeColor = 'k';

% plot selection box
fa = 0.1;
patch(roibox([1 2 2 1],2),roibox([1 1 2 2],3),roibox([1 1 1 1],1),'k','facealpha',fa)
patch(roibox([1 2 2 1],2),roibox([1 1 2 2],3),roibox([2 2 2 2],1),'k','facealpha',fa)
patch(roibox([1 1 1 1],2),roibox([1 2 2 1],3),roibox([1 1 2 2],1),'k','facealpha',fa)
patch(roibox([2 2 2 2],2),roibox([1 2 2 1],3),roibox([1 1 2 2],1),'k','facealpha',fa)
patch(roibox([1 1 2 2],2),roibox([1 1 1 1],3),roibox([1 2 2 1],1),'k','facealpha',fa)
patch(roibox([1 1 2 2],2),roibox([2 2 2 2],3),roibox([1 2 2 1],1),'k','facealpha',fa)

grid on

The first group of conversions can be done with regular MATLAB tools. The rest are done using tools from MIMT, though LCHab and LCHuv can be done with existing tools. It's just not as convenient. Generally, simple transformations like HSV/HSL/HSI or luma-chroma models like YCbCr/YPbPr/YIQ will be much faster than LAB/LUV/etc. I don't think anything fancy is needed here. I'm just including a bunch of things to demonstrate that something as simple as YCbCr will probably work fine, and that there aren't really any esoteric miracle solutions that do remarkably better.

The particular properties of the various model don't matter that much right now, so long as the color points are distributed in a way that makes them (more) easily separated from their neighbors using a simple 3D geometry (a rectangular prism in this example). RIght now, the strategy is still squarely in the realm of "trying to make the best of a bad situation".

We're still trying to make a low-tolerance estimation of an unknown 3D boundary with only a very sparse test set. The apparent margins provided by the various models should be considered very much subject to change as the number of samples increases. What seems best now may have little margins with some other image, or like HSY in this case, it may simply not be able to make the selection with a simple box.

[Image Analyst]

I think the main takeaway from all this discussion is that the class 1 rice (red spots) and class 2 rice (blue spots) are not distinguishable since the red spots lie in the middle of a huge cluster of colors from class 2. They're totally intermingled and we can't get them isolated with reasonable segmentation methods. It's highly unlikely that we'd need to have a mask of some bizarre, irregular shape to isolate the red points. I don't thing that would represent the real world situation at all. So, to me, the bottom line is that the two rice classes are not different colors despite what the original poster wants to believe. At least not in the photo. Maybe with the naked eye they are, but if they are not distinguishable in the photo then the photo may have a problem like it's over exposed or not enough resolution.

[DGM]

I think that the need for an elaborate mask and the strong need to exploit a particular color model are both good signs that automated segmentation is not going to be robust. I think I already mentioned that. Maybe it was a different thread. I can't remember anymore.

[Image Analyst]

Yes you did mention it. I just wanted to concur and drive the point home. Hopefully @Ankit Rathore is still monitoring his question (but I doubt it).

[Ankit Rathore]

@Image Analyst I am still on this problem. As of now computational efficiency is not of major concern for me , currently my focus is to identify grains with utmost efficiency( ability to separate the grains).

[Image Analyst]

OK, good luck with that. I've done all I can do with this image. And I don't think there's anything we can do short of putting in ridiculously restrictive and specific criteria that will work to find the circled grains in this image only, but likely won't work for any other images.

[Ankit Rathore]

What if make a setup in which the lightining conditions will be same for every image i take. In that case I can apply those specific criteria to other images also.. @Image Analyst I am open to any other suggestion also . If you were have been at my place what would you have do?? I am restricted to image processing only. I have hyperspectral imaging here but I have to find a solution where a simple laptop and matlab would help.

[Image Analyst]

Hyper spectral could help. Post those images. Take a picture with all one type of rice on the left and all the other type on the right.

I've been designing imaging systems for about 40 years so I know a little bit about it. 🙂 Other than that (using your hyperspectral images) you should have a light booth that excludes ambient light. You should also check the illumination level with a lux meter to make sure it's the same each time. Aim for somewhere between 1000 and 6000 lux. Take a "blank shot" of just the background so you know how the exposure varies with location in the image due to lens shading and light non-uniformities. Then divide your test images by that background image. See attached demo.

But the big problem is your sample presentation. Your rice you said is all mixed together and you don't know which is which. So, using image analysis, let's say you figure out which type each grain of rice is. Now you need to collect those into bags or buckets (each of one type) to separate them. I guess you'd have to do that by hand, which is subject to error. Plus it would be slow and tedious - definitely not a high throughput method where you could sort through hundreds of kg per day. So maybe you could have the rice on a moving belt and analyze just one at a time. Then, depending on what type of rice it is, a blower tube blasts each grain of rice to the right or left (off the belt) into the proper bucket. Then you'd need to have a bucket fill indicator so you can either automatically switch out the bucket or shutdown the system until the user can put in a new bucket.

[Ankit Rathore]

@Image Analyst The attached demo is missing

Answer 2

[Benjamin Thompson]

You will have to decide how to separate them based on your intended purpose. See the documentation on the regionprops function in the Image Processing Toolbox for a variety of shape measurements that can be calculated. Your image also is grayscale or very nearly so, which means your ability to separate by color will be difficult. Try running the colorThresholder tool on it to create some masks.

Answer 3

[Image Analyst]

How many colors are you seeing there? Because most of us see one : white. In fact it looks like some of the grains might be saturated (overexposed) so they are clipped at 255. If there are two types of rice there you might try hyperspectral imaging. Maybe as some narrow band of light there will be some brightness differentiation.