Shade/color matching of the input image with reference shade chart

Dear Friends,

I have to segment different shades of colour from an image. There is name for each reference color. I have utlisied many suggested methods here like colorlab/edge detection etc. But they did not work as my color of interest is all white based and closer (in colour) to each other.

I am looking for a method to match the image with reference shade as in the figure below.

Please someone help me

Thanks in advance.

Answer 1

[Image Analyst]

You need to calculate the color difference. See lengthy discussions I had with this poster who wanted to match grass to standard colored tiles. http://www.mathworks.com/matlabcentral/answers/contributors/3595115-elvin/questions

[Kumar]

OMG Thanks a lot. Infact I used your method and got some positve results. Actually my task is recurring with many input images with the same standard shade chart. This method necessitates me having both chart and image as an input and it enables me to find only single colour. Is it possible to mark multiple colour and map into label?

Thank you so much for your kind help

[Image Analyst]

Yes it is. If you have simple RGB graphics, you can probably just calculate the Euclidean distance in RGB color space and whichever standard is closest, you assign it to that. Are your standard colors in an image? Do you have the RGB values for them? If not, and you have an image, you can measure them from the image. Then for each standard color, calculate the Euclidean rgb distance:

deltaR = rgbImage(:,:,1) - standardR;
deltaG = rgbImage(:,:,2) - standardG;
deltaB = rgbImage(:,:,3) - standardB;
colorDiffImage = sqrt(deltaR .^ 2 + deltaG .^ 2 + deltaB .^ 2);

Do that for every standard then find out which one each pixel is closest to. Let me see your code for that because there are efficient (min function) and inefficient (for loop) ways of doing that.

Did you see my File Exchange with color segmentation demos: http://www.mathworks.com/matlabcentral/fileexchange/?term=authorid%3A31862

[Kumar]

Thank you so much for valuable input. Yes I am seeing your demo files and learning. It's a great resource.

Unfortunately my knowledge in this domain is very limited. :-( wish I learn and improve.

Interestingly I have L*a*b values for my reference shades as follow:

Shade = (L,a,b)

A1= (82.56,1.37,12.33)

B4=(78.43,3.82, 21.35)

C4= (74.44,4.02,20.23)

D4= (73.68,5.41,23.09)

Is it possible to get compared all these colours simultaneously from an image?

Thanks again Have a nice day

[Image Analyst]

Yes, if you put all the L into one array, all the b into another and all the b into another.

standardLs = [A1(1), B4(1), C4(1), D4(1)];
standardAs = [A1(2), B4(2), C4(2), D4(2)];
standardBs = [A1(3), B4(3), C4(3), D4(3)];
deltaLs = standardLs - testL;
deltaAs = standardAs - testA;
deltaBs = standardBs - testB;
deltaEs = sqrt(deltaLs .^ 2 + deltaAs  .^ 2 + deltaBs  .^ 2);
[closestDeltaE, index] = min(deltaEs);

[Kumar]

Thanks again. Now I am getting dimensional error. I tried to implement your methedology of deltaE.It turns dimensional error.

Please forgive me if I am foolish.

I read image using imread and then tried to convert into Lab colourspace using your deltaE example. but getting error:

>> fullImageFileName = 'teeth.jpg';
>> % Read in image into an array.
  [rgbImage storedColorMap] = imread(fullImageFileName); 
  [rows columns numberOfColorBands] = size(rgbImage); 
  % If it's monochrome (indexed), convert it to color. 
  % Check to see if it's an 8-bit image needed later for scaling).
  if strcmpi(class(rgbImage), 'uint8')
    % Flag for 256 gray levels.
    eightBit = true;
  else
    eightBit = false;
  end
  if numberOfColorBands == 1
    if isempty(storedColorMap)
      % Just a simple gray level image, not indexed with a stored color map.
      % Create a 3D true color image where we copy the monochrome image into all 3 (R, G, & B) color planes.
      rgbImage = cat(3, rgbImage, rgbImage, rgbImage);
    else
      % It's an indexed image.
      rgbImage = ind2rgb(rgbImage, storedColorMap);
      % ind2rgb() will convert it to double and normalize it to the range 0-1.
      % Convert back to uint8 in the range 0-255, if needed.
      if eightBit
        rgbImage = uint8(255 * rgbImage);
      end
    end
  end
  
  % Display the original image.
  h1 = subplot(3, 4, 1);
  imshow(rgbImage);
  drawnow; % Make it display immediately. 
  if numberOfColorBands > 1 
    title('Original Color Image', 'FontSize', fontSize); 
  else 
    caption = sprintf('Original Indexed Image\n(converted to true color with its stored colormap)');
    title(caption, 'FontSize', fontSize);
  end
Undefined function or variable 'fontSize'.

>> cform = makecform('srgb2lab');
  lab_Image = applycform(im2double(rgbImage),cform);
  
  % Extract out the color bands from the original image
  % into 3 separate 2D arrays, one for each color component.
  LChannel = lab_Image(:, :, 1); 
  aChannel = lab_Image(:, :, 2); 
  bChannel = lab_Image(:, :, 3); 
>> % Display the lab images.
  subplot(3, 4, 2);
  imshow(LChannel, []);
  title('L Channel', 'FontSize', fontSize);
  subplot(3, 4, 3);
  imshow(aChannel, []);
  title('a Channel', 'FontSize', fontSize);
  subplot(3, 4, 4);
  imshow(bChannel, []);
  title('b Channel', 'FontSize', fontSize);

  % Get the average lab color value.
  [LMean, aMean, bMean] = GetMeanLABValues(LChannel, aChannel, bChannel, mask);
Undefined function or variable 'fontSize'.

>>  LStandard = [82.56,78.43,74.44,73.68];
    aStandard = [1.37,3.82,4.02,5.41];
    bStandard = [12.33,21.35,20.23,23.09];
  
  
  % Create the delta images: delta L, delta A, and delta B.
  deltaL = LChannel - LStandard;
  deltaa = aChannel - aStandard;
  deltab = bChannel - bStandard;
  
  % Create the Delta E image.
  % This is an image that represents the color difference.
  % Delta E is the square root of the sum of the squares of the delta images.
  deltaE = sqrt(deltaL .^ 2 + deltaa .^ 2 + deltab .^ 2);
  
  % Mask it to get the Delta E in the mask region only.
  maskedDeltaE = deltaE .* mask;
  % Get the mean delta E in the mask region
  % Note: deltaE(mask) is a 1D vector of ONLY the pixel values within the masked area.
   meanMaskedDeltaE = mean(deltaE(mask));
  % Get the standard deviation of the delta E in the mask region
  stDevMaskedDeltaE = std(deltaE(mask));
  message = sprintf('The mean LAB = (%.2f, %.2f, %.2f).\nThe mean Delta E in the masked region is %.2f +/- %.2f',...
    LMean, aMean, bMean, meanMaskedDeltaE, stDevMaskedDeltaE);  
  
  % Display the masked Delta E image - the delta E within the masked region only.
  subplot(3, 4, 6);
  imshow(maskedDeltaE, []);
  caption = sprintf('Delta E between image within masked region\nand mean color within masked region.\n(With amplified intensity)');
  title(caption, 'FontSize', fontSize);

  % Display the Delta E image - the delta E over the entire image.
  subplot(3, 4, 7);
  imshow(deltaE, []);
  caption = sprintf('Delta E Image\n(Darker = Better Match)');
  title(caption, 'FontSize', fontSize);

  % Plot the histograms of the Delta E color difference image,
  % both within the masked region, and for the entire image.
  PlotHistogram(deltaE(mask), deltaE, [3 4 8], 'Histograms of the 2 Delta E Images');
  
  message = sprintf('%s\n\nRegions close in color to the color you picked\nwill be dark in the Delta E image.\n', message);
  msgboxw(message);

  % Find out how close the user wants to match the colors.
  prompt = {sprintf('First, examine the histogram.\nThen find pixels within this Delta E (from the average color in the region you drew):')};
  dialogTitle = 'Enter Delta E Tolerance';
  numberOfLines = 1;
  % Set the default tolerance to be the mean delta E in the masked region plus two standard deviations.
  strTolerance = sprintf('%.1f', meanMaskedDeltaE + 3 * stDevMaskedDeltaE);
  defaultAnswer = {strTolerance};  % Suggest this number to the user.
  response = inputdlg(prompt, dialogTitle, numberOfLines, defaultAnswer);
  % Update tolerance with user's response.
  tolerance = str2double(cell2mat(response));

  % Let them interactively select the threshold with the threshold() m-file.
  % (Note: This is a separate function in a separate file in my File Exchange.)
%   threshold(deltaE);
  
  % Place a vertical bar at the threshold location.
  handleToSubPlot8 = subplot(3, 4, 8);  % Get the handle to the plot.
  PlaceVerticalBarOnPlot(handleToSubPlot8, tolerance, [0 .5 0]);  % Put a vertical red line there.
  
  % Find pixels within that delta E.
  binaryImage = deltaE <= tolerance;
  subplot(3, 4, 9);
  imshow(binaryImage, []);
  title('Matching Colors Mask', 'FontSize', fontSize);
  
  % Mask the image with the matching colors and extract those pixels.
  matchingColors = bsxfun(@times, rgbImage, cast(binaryImage, class(rgbImage)));
  subplot(3, 4, 10);
  imshow(matchingColors);
  caption = sprintf('Matching Colors (Delta E <= %.1f)', tolerance);
  title(caption, 'FontSize', fontSize);

  % Mask the image with the NON-matching colors and extract those pixels.
  nonMatchingColors = bsxfun(@times, rgbImage, cast(~binaryImage, class(rgbImage)));
  subplot(3, 4, 11);
  imshow(nonMatchingColors);
  caption = sprintf('Non-Matching Colors (Delta E > %.1f)', tolerance);
  title(caption, 'FontSize', fontSize);

  % Display credits: the MATLAB logo and my name.
  ShowCredits(); % Display logo in plot position #12.
  
  % Alert user that the demo has finished.
  message = sprintf('Done!\n\nThe demo has finished.\nRegions close in color to the color you picked\nwill be dark in the Delta E image.\n');
  msgbox(message);
  
catch ME
  errorMessage = sprintf('Error running this m-file:\n%s\n\nThe error message is:\n%s', ...
    mfilename('fullpath'), ME.message);
  errordlg(errorMessage);
end
return; % from SimpleColorDetection()
% ---------- End of main function ---------------------------------

%----------------------------------------------------------------------------
% Display the MATLAB logo.
function ShowCredits()
try
%   xpklein;
%   surf(peaks(30));
  logoFig = subplot(3, 4, 12);
  caption = sprintf('A MATLAB Demo\nby ImageAnalyst');
  text(0.5,1.15, caption, 'Color','r', 'FontSize', 18, 'FontWeight','b', 'HorizontalAlignment', 'Center') ;
  positionOfLowerRightPlot = get(logoFig, 'position');
  L = 40*membrane(1,25);
  logoax = axes('CameraPosition', [-193.4013 -265.1546  220.4819],...
    'CameraTarget',[26 26 10], ...
    'CameraUpVector',[0 0 1], ...
    'CameraViewAngle',9.5, ...
    'DataAspectRatio', [1 1 .9],...
    'Position', positionOfLowerRightPlot, ...
    'Visible','off', ...
    'XLim',[1 51], ...
    'YLim',[1 51], ...
    'ZLim',[-13 40], ...
    'parent',gcf);
  s = surface(L, ...
    'EdgeColor','none', ...
    'FaceColor',[0.9 0.2 0.2], ...
    'FaceLighting','phong', ...
    'AmbientStrength',0.3, ...
    'DiffuseStrength',0.6, ... 
    'Clipping','off',...
    'BackFaceLighting','lit', ...
    'SpecularStrength',1.1, ...
    'SpecularColorReflectance',1, ...
    'SpecularExponent',7, ...
    'Tag','TheMathWorksLogo', ...
    'parent',logoax);
  l1 = light('Position',[40 100 20], ...
    'Style','local', ...
    'Color',[0 0.8 0.8], ...
    'parent',logoax);
  l2 = light('Position',[.5 -1 .4], ...
    'Color',[0.8 0.8 0], ...
    'parent',logoax);

catch ME
  errorMessage = sprintf('Error running ShowCredits().\n\nThe error message is:\n%s', ...
    ME.message);
  errordlg(errorMessage);
end
return; % from ShowCredits()

%-----------------------------------------------------------------------------
function [xCoords, yCoords, roiPosition] = DrawBoxRegion(handleToImage)
  try
  % Open a temporary full-screen figure if requested.
  enlargeForDrawing = true;
  axes(handleToImage);
  if enlargeForDrawing
    hImage = findobj(gca,'Type','image');
    numberOfImagesInside = length(hImage);
    if numberOfImagesInside > 1
      imageInside = get(hImage(1), 'CData');
    else
      imageInside = get(hImage, 'CData');
    end
    hTemp = figure;
    hImage2 = imshow(imageInside, []);
    [rows columns NumberOfColorBands] = size(imageInside);
    set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
  end
  
  txtInfo = sprintf('Draw a box over the unstained fabric by clicking and dragging over the image.\nDouble click inside the box to finish drawing.');
  text(10, 40, txtInfo, 'color', 'r', 'FontSize', 24);

      % Prompt user to draw a region on the image.
    msgboxw(txtInfo);
    
    % Erase all previous lines.
    if ~enlargeForDrawing
      axes(handleToImage);
  %     ClearLinesFromAxes(handles);
    end
    
    hBox = imrect;
    roiPosition = wait(hBox);
    roiPosition
    % Erase all previous lines.
    if ~enlargeForDrawing
      axes(handleToImage);
  %     ClearLinesFromAxes(handles);
    end

  xCoords = [roiPosition(1), roiPosition(1)+roiPosition(3), roiPosition(1)+roiPosition(3), roiPosition(1), roiPosition(1)];
  yCoords = [roiPosition(2), roiPosition(2), roiPosition(2)+roiPosition(4), roiPosition(2)+roiPosition(4), roiPosition(2)];

  % Plot the mask as an outline over the image.
  hold on;
  plot(xCoords, yCoords, 'linewidth', 2);
  close(hTemp);
  catch ME
    errorMessage = sprintf('Error running DrawRegion:\n\n\nThe error message is:\n%s', ...
      ME.message);
    WarnUser(errorMessage);
  end
  return; % from DrawRegion
  
%-----------------------------------------------------------------------------
function [mask] = DrawFreehandRegion(handleToImage, rgbImage)
try
  fontSize = 14;
  % Open a temporary full-screen figure if requested.
  enlargeForDrawing = true;
  axes(handleToImage);
  if enlargeForDrawing
    hImage = findobj(gca,'Type','image');
    numberOfImagesInside = length(hImage);
    if numberOfImagesInside > 1
      imageInside = get(hImage(1), 'CData');
    else
      imageInside = get(hImage, 'CData');
    end
    hTemp = figure;
    hImage2 = imshow(imageInside, []);
    [rows columns NumberOfColorBands] = size(imageInside);
    set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
  end
  
  message = sprintf('Left click and hold to begin drawing.\nSimply lift the mouse button to finish');
  text(10, 40, message, 'color', 'r', 'FontSize', fontSize);

      % Prompt user to draw a region on the image.
    uiwait(msgbox(message));
    
    % Now, finally, have the user freehand draw the mask in the image.
    hFH = imfreehand();

  % Once we get here, the user has finished drawing the region.
  % Create a binary image ("mask") from the ROI object.
  mask = hFH.createMask();
  
  % Close the maximized figure because we're done with it.
  close(hTemp);

  % Display the freehand mask.
  subplot(3, 4, 5);
  imshow(mask);
  title('Binary mask of the region', 'FontSize', fontSize);
  
  % Mask the image.
  maskedRgbImage = bsxfun(@times, rgbImage, cast(mask,class(rgbImage)));
  % Display it.
  subplot(3, 4, 6);
  imshow(maskedRgbImage);
catch ME
  errorMessage = sprintf('Error running DrawFreehandRegion:\n\n\nThe error message is:\n%s', ...
    ME.message);
  WarnUser(errorMessage);
end
return; % from DrawFreehandRegion

%-----------------------------------------------------------------------------
% Get the average lab within the mask region.
function [LMean, aMean, bMean] = GetMeanLABValues(LChannel, aChannel, bChannel, mask)
try
  LVector = LChannel(mask); % 1D vector of only the pixels within the masked area.
  LMean = mean(LVector);
  aVector = aChannel(mask); % 1D vector of only the pixels within the masked area.
  aMean = mean(aVector);
  bVector = bChannel(mask); % 1D vector of only the pixels within the masked area.
  bMean = mean(bVector);
catch ME
  errorMessage = sprintf('Error running GetMeanLABValues:\n\n\nThe error message is:\n%s', ...
    ME.message);
  WarnUser(errorMessage);
end
return; % from GetMeanLABValues

%==========================================================================================================================
function WarnUser(warningMessage)
  uiwait(warndlg(warningMessage));
  return; % from WarnUser()
  
%==========================================================================================================================
function msgboxw(message)
  uiwait(msgbox(message));
  return; % from msgboxw()
  
%==========================================================================================================================
% Plots the histograms of the pixels in both the masked region and the entire image.
function PlotHistogram(maskedRegion, doubleImage, plotNumber, caption)
try
  fontSize = 14;
  subplot(plotNumber(1), plotNumber(2), plotNumber(3)); 

  % Find out where the edges of the histogram bins should be.
  maxValue1 = max(maskedRegion(:));
  maxValue2 = max(doubleImage(:));
  maxOverallValue = max([maxValue1 maxValue2]);
  edges = linspace(0, maxOverallValue, 100);

  % Get the histogram of the masked region into 100 bins.
  pixelCount1 = histc(maskedRegion(:), edges);

  % Get the histogram of the entire image into 100 bins.
  pixelCount2 = histc(doubleImage(:), edges);

  % Plot the  histogram of the entire image.
  plot(edges, pixelCount2, 'b-');
  
  % Now plot the histogram of the masked region.
  % However there will likely be so few pixels that this plot will be so low and flat compared to the histogram of the entire
  % image that you probably won't be able to see it.  To get around this, let's scale it to make it higher so we can see it.
  gainFactor = 1.0;
  maxValue3 = max(max(pixelCount2));
  pixelCount3 = gainFactor * maxValue3 * pixelCount1 / max(pixelCount1);
  hold on;
  plot(edges, pixelCount3, 'r-');
  title(caption, 'FontSize', fontSize);
  
  % Scale x axis manually.
  xlim([0 edges(end)]);
  legend('Entire', 'Masked');
  
catch ME
  errorMessage = sprintf('Error running PlotHistogram:\n\n\nThe error message is:\n%s', ...
    ME.message);
  WarnUser(errorMessage);
end
return; % from PlotHistogram

%=====================================================================
% Shows vertical lines going up from the X axis to the curve on the plot.
function lineHandle = PlaceVerticalBarOnPlot(handleToPlot, x, lineColor)
  try
    % If the plot is visible, plot the line.
    if get(handleToPlot, 'visible')
      axes(handleToPlot);  % makes existing axes handles.axesPlot the current axes.
      % Make sure x location is in the valid range along the horizontal X axis.
      XRange = get(handleToPlot, 'XLim');
      maxXValue = XRange(2);
      if x > maxXValue
        x = maxXValue;
      end
      % Erase the old line.
      %hOldBar=findobj('type', 'hggroup');
      %delete(hOldBar);
      % Draw a vertical line at the X location.
      hold on;
      yLimits = ylim;
      lineHandle = line([x x], [yLimits(1) yLimits(2)], 'Color', lineColor, 'LineWidth', 3);
      hold off;
    end
  catch ME
    errorMessage = sprintf('Error running PlaceVerticalBarOnPlot:\n\n\nThe error message is:\n%s', ...
      ME.message);
    WarnUser(errorMessage);
end
  return;  % End of PlaceVerticalBarOnPlot
Error using  - 
Matrix dimensions must agree.

>>

[Image Analyst]

Assign fontSize somewhere in the beginning

fontSize = 13; % or whatever...

[Kumar]

Thanks still it returns error. I simplified further through to understand the basic. It gives me too many error:

imread('teeth.jpg')
%load rgb image
src = 'teeth.jpg';
rgbI = imread(src);

%convert to lab
labTransformation = makecform('srgb2lab');
labI = applycform(rgbI,labTransformation);

%seperate l,a,b
l = labI(:,:,1);
a = labI(:,:,2);
b = labI(:,:,3);

figure, imshow(l) , title('l');
figure, imshow(a) , title('a');
figure, imshow(b) , title('b');
 LStandard = [82.56,78.43,74.44,73.68];
 aStandard = [1.37,3.82,4.02,5.41];
 bStandard = [12.33,21.35,20.23,23.09];
  
  
  % Create the delta images: delta L, delta A, and delta B.
  deltaL = l - LStandard;
  deltaa = a - aStandard;
  deltab = b - bStandard;
  
  % Create the Delta E image.
  % This is an image that represents the color difference.
  % Delta E is the square root of the sum of the squares of the delta images.
  deltaE = sqrt(deltaL .^ 2 + deltaa .^ 2 + deltab .^ 2);
    figure, imshow(deltaE)

This gives me error like:

Error using - Integers can only be combined with integers of the same class, or scalar doubles.

Error in Untitled2 (line 24) deltaL = l - LStandard;

[Image Analyst]

Make sure everything is double before you use it.

l = double(labI(:,:,1));

same for a and b.

[Kumar]

You are my saviour. I have done that. now again: Error using - Matrix dimensions must agree. Error using - Integers can only be combined with integers of the same class, or scalar doubles.

Error in final (line 23) deltaL = l - LStandard;

[Image Analyst]

Since l is an entire image, you're going to have to subtract one LStandard value at a time:

deltaL = l - LStandard(1); % deltaL is also an image.

repeat for LStandard(2), LStandard(3), and LStandard(4).

[Rifayet Ashraf]

Would you please tell me , how to run this full code? I cannot run because of error. Is there any specific initialization or any default settings? the m.file is here what I am trying to run.

[Image Analyst]

You forgot to attach 'C:\Users\user\Desktop\FYP pic\Single teeth.JPG'

[Walter Roberson]

Addressed in https://www.mathworks.com/matlabcentral/answers/381436-code-error-of-shade-matching-image-processing