Remove noise from screw image

The screw image shown below has some dust/fur on the screw thread. What are the ways to remove those noise and produce a new image with clear outline of the screw?

I tried it myself and came out result as shown below.

Some of the noises are still not removed. I have also attached my codes as shown below.

clc
clear
close all

%% Read
A = imread('20.bmp');

%% Turn into grayscale
I0 = rgb2gray(A);

%% Crop
I1 = imcrop(I0,[1,78,710-1,410-78]);
%% Histogram
[pixelCount, grayLevels] = imhist(I0);

%% Adjustment
I2 = imadjust(I1);
I3 = medfilt2(I2);
ThresholdValue = graythresh(I3);
I4 = imbinarize(I3, ThresholdValue);
I5 = edge(I4,'Sobel');

Answer 1

[Image Analyst]

The screw and the threads are the black thing and the dust is the lighter things on it. I think you are just using the wrong threshold. You can adjust it in imbinarize(), or you can do it interactively with my interactive thresholding app in my File Exchange.

[Ang Xian Jia]

Hi Image Analyst. How to correctly determine the desired threshold value through your interactive file because it needs to be done manually. Does this mean that I should manually determine threshold value for different image with same type of screw?

[Image Analyst]

If you have good control over your exposure and illumination and camera settings then you can do it interactively ONCE to find out what the threshold should be, then just use that value as a fixed threshold from that point on.

[Ang Xian Jia]

Hi Image Analyst,

Is the image below considered ok? I cant get a very good fixed threshold using the interactive tool. Would you mind to help me take a look?

Answer 2

[Stephan]

Replace the medfilt2 by a gaussian blur using imgaussfilt:

%% read image and find points that should lie on a straight line
A = imread('20.bmp');
% Turn into grayscale
I0 = rgb2gray(A);
% Crop
I1 = imcrop(I0,[1,78,710-1,410-78]);
% Histogram
[pixelCount, grayLevels] = imhist(I0);
% Adjustment
I2 = imadjust(I1);
I3 = imgaussfilt(I2,4);
ThresholdValue = graythresh(I3);
I4 = imbinarize(I3, ThresholdValue);
I5 = edge(I4,'Sobel');
imshowpair(I4,I5,'montage')
gaussblur.PNG

For adjusting i suggest to play around with the standard deviation of imgaussfilt. Here i used a value of 4.

BTW:

Did you notice that you can accept and/or vote for useful answers? If they dont help you - leave a comment to discuss what you expect to get.

Answer 3

[Image Analyst]

Here is a full demo.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
imtool close all;  % Close all imtool figures if you have the Image Processing Toolbox.
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;
% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end
%===============================================================================
% Read in a standard MATLAB color demo image.
folder = pwd;
baseFileName = 'screw.bmp';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % Didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, 'file')
    % Still didn't find it.  Alert user.
    errorMessage = sprintf('Error: %s does not exist.', fullFileName);
    uiwait(warndlg(errorMessage));
    return;
  end
end
grayImage = imread(fullFileName);
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows, columns, numberOfColorBands] = size(grayImage);
if numberOfColorBands > 1
	grayImage = rgb2gray(grayImage);
end
% Display the original color image.
subplot(2, 2, 1);
imshow(grayImage);
axis on;
title('Original Image', 'FontSize', fontSize);
impixelinfo;
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'Outerposition', [0, 0.1, 1, 0.9]);

subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of Gray Scale Image', 'FontSize', fontSize);

% Get the binaryImage
binaryImage = grayImage < 200;
% There is a white frame around it.  Get rid of that.
% binaryImage = imclearborder(binaryImage);
% Take the largest blob only
binaryImage = bwareafilt(binaryImage, 1);
% Display the image.
subplot(2, 2, 3);
imshow(binaryImage);
axis('on', 'image');
title('Binary Image', 'FontSize', fontSize);

% Get the top row
topRows = nan(1, columns);
bottomRows = nan(1, columns);
for col = 1 : columns
	thisColumn = binaryImage(:, col);
	t = find(thisColumn);
	if ~isempty(t) % If there is some white in this column
		topRows(col) = t(1);
		bottomRows(col) = t(end);
	end
end

% Smooth them a bit by smoothing with a Savitzky-Golay filter.
% This will scan with a moving quadratic fit window.
topRows = sgolayfilt(topRows, 2, 11);

% Plot boundary over gray scale image.
subplot(2, 2, 4);
imshow(binaryImage);
axis('on', 'image');
hold on;
plot(1:columns, topRows, 'r-', 'LineWidth', 2);
plot(1:columns, bottomRows, 'r-', 'LineWidth', 2);

caption = sprintf('With top and bottom rows overlaid.');
title(caption, 'FontSize', fontSize);

% Now find peaks of the top rows.
x = 1 : columns;
badColumns = isnan(topRows);
x(badColumns) = [];
topRows(badColumns) = [];
% Find peaks.
[topPeaksY, indexes] = findpeaks(topRows);
topPeaksX = x(indexes)
% Plot a circle over them.
plot(topPeaksX, topPeaksY, 'bo', 'LineWidth', 2, 'MarkerSize', 18);
% Find valleys.
[topValleysY, indexes] = findpeaks(-topRows);
topValleysX = x(indexes)
% Plot a circle over them.
plot(topValleysX, -topValleysY, 'co', 'LineWidth', 2, 'MarkerSize', 18);

% Now find peaks of the top rows.
x = 1 : columns;
badColumns = isnan(bottomRows);
x(badColumns) = [];
bottomRows(badColumns) = [];
% Find peaks.
[topPeaksY, indexes] = findpeaks(bottomRows);
topPeaksX = x(indexes)
% Plot a circle over them.
plot(topPeaksX, topPeaksY, 'bo', 'LineWidth', 2, 'MarkerSize', 18);
% Find valleys.
[topValleysY, indexes] = findpeaks(-bottomRows);
topValleysX = x(indexes)
% Plot a circle over them.
plot(topValleysX, -topValleysY, 'co', 'LineWidth', 2, 'MarkerSize', 18);

I used a threshold of 200.

[Ang Xian Jia]

Hi Image Analyst,

May I know the effect of this sentence?

binaryImage = bwareafilt(binaryImage, 1);

Does it mean retaining the largest object, which is the screw? But when i put together the image before and after the filter, I cannot see any differences. Is it to remove the white spot on the top left corner?

[Image Analyst]

It means take the largest object. I do this just in case there might be some small noise in the image, maybe even some that you can't see if they are too small so that they were subsampled away when it went to display the image. If there are no other noise blobs, then it will just return the same image. However, if there WERE other tiny blobs in the image, then that would complicate the subsequent processing so I just make sure there is only one because I KNOW there is supposed to be only one.

[Ang Xian Jia]

Greetings Image Analyst,

Thank you for the explanation. The code identify some of the peaks and valleys incorrectly. This is due to the blurry image(mechanical limitation) and cannott be corrected through software right?

[Ang Xian Jia]

The code correctly identify the peaks and valleys if I use 'MinimumPeakDistance'. Is this considered a good way to solve the problem?

[topPeaksY, indexes] = findpeaks(topRows,'MinPeakDistance',50);

[topValleysY, indexes] = findpeaks(-topRows,'MinPeakDistance',50);

However, the way the code identify the peaks and valleys is on the middle of every pixel.

So, I'm considering using subpixel edge pixel detection, extract the coordinates of edge, and similarly find peaks and valleys after that. Is this a good way to start to improve accuracy of identify peaks and valleys? Btw, thank you so much for the code! I learn a lot.

[Image Analyst]

As you can see, the quantized/digitized boundary is not very accurate, at least I don't believe it is. Sure you can get the actual location on the boundary by not doing the smoothing filter, but who's to say which location is more accurate? You have no ground truth. So why do you think the quantized one is more accurate? However, if you want that, rather than the denoised version, so ahead and use it.

If you want subpixel accuracy you'll have to use the gray scale image and use smoothing (like I did) or some other kind of modeling. With thresholding the best you can get is pixel accuracy. You can only identify the location on an exact pixel, not in between pixels like my code did.

[Ang Xian Jia]

Hi Image Analyst, since you have also identified the top & bottom peak and valley, is your code also suitable to be used to calculate major and minor diameter of the screw? Thank you once again.