19 March, 2016

# A step-by-step guide of an image segmentation task

Some weeks ago I received an e-mail from a reader pointing to an interesting task in the field of image processing. He kindly let me to publish the problem and the images we were working on, so I decided to explain the problem and write a small guide showing how such a task can be solved.

Have a look at the input image below. The goal was to roughly identify the scratched region found in the lower part of the image. The solution also had to work for similar inputs.

As the first step it is important consider the possible different approaches and to find such properties of the image or a-priori knowledge which may simplify the problem. For example here are some of the pieces of information I was interested in:

• Do all input images have the same size?
• Do the scratch marks have the same or similar width as on the example image?
• Is the run-time important?
• Is the scratch marked region always contiguous?

It turned out, that all images will have the same size and the scratch marks also have almost the same width on each image. The scratchy region is always contiguous. I was also told that the run-time is not so important. (Do not believe it to anyone: it is always important... :)

By knowing all the required information the skeleton of a possible solution is the following:

• As color information is not important, convert the input image to gray-scale.
• Apply some preprocessing if required.
• Use a matched filter which is sensitive to horizontal stripes.
• Identify the high-response regions.
• Select the largest contiguous region as the result.

One important thing to note: the size of the input image was 1294x930 pixels. Do we really need this high resolution? Horizontal stripes remain horizontal stripes even in smaller images. Because the accuracy of the identified regions is not so important, it is worth to scale the input by a factor of 1/4. This way the information to process is reduced.

In this post I used Octave, which is very similar to MATLAB, but the source code may require some tuning...

The preprocessing part is the following:

% Read the image
imageFile = 'image002.png';

% Select the green channel
% Note: There is no need to convert the image to gray as it is almost gray already
image = image(:, :, 2);

% Convert to double: required if you use Octave
image = double(image);

% Resize the image to the 1/4th size
image = imresize(image, 0.25);

% Some noise smoothing
image = imsmooth(image, 'Gaussian', 1.00);


After these steps we have a smaller and noise filtered image:

Now we have to apply a filter which is sensitive to horizontal stripes. In this post we will use a rather simple stripe filter. The stripeness can be measured with the following formula which works on a single column or row of an image:

$$R_{c,m} = 2 \cdot I_c - I_{c-m} - I_{c+m}$$

R is the response of the processed pixel in the current row or column: we subtract the I intensity of the two side pixels at distance m from the central, cth pixel. This approach works because of the following reasons:

• If the central pixel is bright while the ones on its side are dark, the response is a large positive number.
• If the central pixel is dark while the ones on its side are bright, the response is a large negative number.
• If the central pixel and the ones on its side have similar intensity, we get a result around zero.

This filter can be easily applied by using convolution. In addition because the kernel is one-dimensional, it will be very fast. We can improve the results by taking the absolute value of the response: between two bright stripes there is a dark stripe too, so all extreme values are important for us despite of their sign.

Because the widths of the stripes vary in a range we have to apply filters created for different widths. This is another advantage of reducing the image size: the possible widths of stripes has also decreased, which results in lower computational effort. Let us apply the filters for different widths and store the maximum response for each pixel:

$$R_{c,final} = max(R_{c,m}), m \in [1, 8]$$

Because we know that the scratch marks look very similar on each input image, their width – measured on the example image – must be always in the range of [1, 8] pixels. The kernel family we use for convolution are demonstrated below. Please note the transpose operator as we have to use vertical kernels:

% Kernel family to enhance horizontal stripes
k1 = [-1 2 -1]';
k2 = [-1 0 2 0 -1]';
k3 = [-1 0 0 2 0 0 -1]';
k4 = [-1 0 0 0 2 0 0 0 -1]';
k5 = [-1 0 0 0 0 2 0 0 0 0 -1]';


The piece of code for filtering is the following, the kernels are generated on the fly:

% Create an accumulator to store maximum filter response
accum = zeros(size(image));

% Create kernels of different widths
for w = 0 : 4
% Create the current kernel
kernel = double([-1; zeros(w, 1); 2; zeros(w, 1); -1]);

% Calculate absolute response
response = abs(conv2(image, kernel,  'same'));

% Store the maximum filter response in the accumulator
accum = max(accum, response);
end


On the result it can be seen, that the filter highlights the scratchy regions, but also the borders of the objects:

In the next step the regions having high response shall be segmented. First use a median filter to smooth the results of the filters, run a threshold on it, then apply a median filter again to smooth the binary regions. The code example is:

% Smoothen the result by median filtering
accum = medfilt2(accum, [17, 17]);

% Thresholding can be improved later, now use a fixed one
binary = accum > 22;

% Median filter again for noise and shape smoothing
binary = medfilt2(binary, [17, 17]);


The next three images show the progress:

As we know the size of the input image, we can use a median filter of a fixed size. But the applied threshold value needs to be tuned later or be calculated by a much clever method. For demonstration purposes it worked well.

As it can be seen on the final picture, border regions are also highlighted. We can remove these thin regions by a morphological eroding and dilating:

% Apply morphological operations to remove thin regions
binary = imerode(binary, ones(17));
binary = imdilate(binary, ones(13));


The result – highlighted on the original image – is the following:

The scratchy region is well identified. Although there are some false regions too, they can be easily removed for example by finding the largest blob. This time I let this step to the reader. :)