My custom implementation of MATLAB's imdilate

Question

I have been assigned to implement the morphological dilation operation without using MATLAB's imdilate, and I came up with the following solution:

% -------------------------------------------------------------------------
% -> Morphological Operations
% -> Image Dilation
% -------------------------------------------------------------------------
% -------------------------------------------------------------------------
% - [Pre-Operation Code: Getting Everything Ready] -
% -------------------------------------------------------------------------
% Starting With A Clean (Workspace) And (Command Window)
clear;
clc;
% Creating A Matrix Consisting Of 0s And 1s Representing A Binary Image.
binaryImage = [ ...
 0 0 0 0 0 0
 0 0 1 1 0 0
 0 0 1 1 0 0
 0 0 0 0 0 0];
% Creating A Matrix Representing Our Structuring Element
structuringElement = [ ...
 0 1 0
 1 1 1
 0 1 0];
% Getting The Number Of Rows (Height) And The Number Of Columns (Width) Of The Binary Image
[imageRows, imageColumns] = size(binaryImage);
% Getting The Number Of Rows (Height) And The Number Of Columns (Width) Of The Structuring Element
[structuringRows, structuringColumns] = size(structuringElement);
% Creating An Empty Matrix That Will Be Used To Store The Final Processed Image
dilatedImage = zeros(imageRows, imageColumns);
ref = imdilate(binaryImage, structuringElement);
% -------------------------------------------------------------------------
% - [Dilation Operation] -
% -------------------------------------------------------------------------
% Going Over Each Row In The Binary Image
for i = 1:imageRows
 % Going Over Each Column In The Binary Image
 for j = 1:imageColumns
 % If The Current Pixel Is A Foreground Pixel (1)
 if (binaryImage(i, j) == 1)
 % Going Over Each Row In The Structuring Element
 for k = 1:structuringRows
 % Going Over Each Column In The Structuring Element
 for l = 1:structuringColumns
 % If The Current Pixel In The Structuring Element Is A Foreground Pixel (1)
 if (structuringElement(k, l) == 1)
 dilatedImage(i + k - 2, j + l - 2) = 1;
 end
 end
 end
 end
 end
end
subplot(1, 3, 1), imshow(binaryImage), title('Original Image');
subplot(1, 3, 2), imshow(dilatedImage), title('Dilated Image');
subplot(1, 3, 3), imshow(ref), title('MATLAB imdilate');

I am open to any suggestions! Thank you in advance.

Answer 1

I have a few different comments.

Firstly, for 2D arrays like in your example, you can simply do conv2(binaryImage,structuringElement,'same'), it's equivalent to imdilate (and about 4 times faster). That might not be acceptable for your case if you are supposed to write the code yourself though.

Next, and importantly, your code has a logic error. If your binary image has a white pixel along any edge, the results are incorrect or you get an error (if it's in the top row). This is due to your indexing into dilatedImage where i+k-2 can be less than 1, or can be larger than the number of rows in dilatedImage.

You have a lot of comments, which is better than none, but I think it's excessive here. It is very clear what this line of code does [imageRows, imageColumns] = size(binaryImage);, for example, even for people who don't know Matlab. Some of your comments make the lines very long as well.

A minor issue, I'd move the imdilate command somewhere else, that part of the code is where you are re-creating imdilate, I think it's nicer to separate the verification part.

Another minor issue is in using i and j as loop indices. Since i and j can be used for complex numbers in Matlab, by convention they are avoided as loop indices (although I do not usually worry about this).

I would put the code that does the work into its own function that takes binaryImage and structuringElement as inputs and returns dilatedImage.

The big thing here is you 6 nested for/if statements, which is hard to read and hard to code. Historically, such nested loops would have been very slow in Matlab, but nowadays it's not such a concern. One thing you should do is swap the order of the first two for loops. Matlab stores arrays in a column major format, so loop over the columns first, then the rows. This gives roughly a 20-25% speedup in my tests. Additionally, checking to remove the out-of-bounds error is also a speedup, since we only set pixels to one if they are inside the image, irregardless of the value of the structuring array. Check out this version of your loop:

for j = 1:imageColumns
 for i = 1:imageRows
 if binaryImage(i, j) == 1
 % Loop through the entries of the structuring element
 for l = 1:structuringColumns
 for k = 1:structuringRows
 % make sure the new pixel is inside the image
 if i+k-2>0 && j+l-2>0 && i+k-2<=imageRows && j+l-2<=imageColumns
 % make a white pixel if necessary
 if structuringElement(k, l)==1
 dilatedImage(i + k - 2, j + l - 2) = 1;
 end
 end
 end
 end
 end
 end
end

Of course it would be nice to avoid writing all these loops, and that can be done. This next method is only slightly faster, but to me is cleaner. The idea is to first find the white pixels, and just loop over them, for each one working out which surrounding pixels should be changed.

[rowIdx,colIdx] = find(binaryImage); % the indices of white pixels in the image
[maskRowIdx,maskColIdx] = find(structuringElement); % indices of structuring elements
% now get indices of structuring elements relative to the centre of the array
maskRowIdx = maskRowIdx - floor(structuringRows/2) - 1;
maskColIdx = maskColIdx - floor(structuringColumns/2) - 1;
dilatedImage = zeros(imageRows, imageColumns);
% loop over each white pixel in the image
for i = 1:numel(rowIdx)
 % these are just the indices of the pixel
 rI = rowIdx(i);
 cI = colIdx(i);
 % now loop over each non-zero element of the structuring array
 for j = 1:numel(maskRowIdx)
 % the position of the pixel to change is (r,c)
 r = rI+maskRowIdx(j);
 c = cI+maskColIdx(j);
 % if the pixel is inside the image, we make it a 1
 if r>0 && c>0 && r<=imageRows && c<=imageColumns
 dilatedImage(r,c) = 1;
 end
 end
end

• saad.sawash: Thank you a lot for the tips! By the way, what about the erosion operation implementation? Is it possible to implement it following the same code you provided just with some modifications? And what will I modify?
• David: The simplest thing to do is just use the same code but operate on ~binaryImage instead of binaryImage, that is, swap the black and white pixels, do a dilation, swap the pixels back again.