Dilatation and Erosion as custom function for signal processing

5 views (last 30 days)
I want to define the two function for accomplish the morphological operation of erosion and dilatation on a 1-D signal (in this case an ecg signal). I know that I can use functions like imerode and imdilate, but I want to define them by myself, obtaining a result comparable to the functions mentioned before.
To design the two function I have to followed the example propose below:
Based on this knowledge I have defined two functions: dilatation and erosion (below). However when I apply the algorithm (show below in the code) the result is something that I do not understand. Based on the methods presented above I can't understand what I is missing or what is wrong in the two functions.
My Code:
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 15;
% Signal info.
G = 200; % Gain
Fs = 360; % [Hz]
L = 3600; % lenght of ECG signals
T = linspace(0,L/Fs,L); % time axis
F = linspace(-Fs/2, Fs/2, L); % Frequency axis
% Load Signal
load("100m.mat");
ecg = val/G;
% Plot the ecg signal
figure(1); plot(T, ecg); grid on;
title("ECG Signal","FontSize",fontSize);
xlabel("Time (sec)", "FontSize", fontSize);
ylabel("voltage [mV]", "FontSize", fontSize);
% Apply MMF algorithm
% Declaration of the two Structuring Elements
B1 = [0 1 5 1 0]; B2 = [1 1 1 1 1];
% modified morphological filtering algorithm
d1 = dilatation(ecg, B1); % d1 = imdilate(ecg, B1);
e1 = erosion(d1,B2); % e1 = imerode(ecg, B2);
e2 = erosion(ecg,B1); % e2 = imerode(ecg, B1);
d2 = dilatation(e2,B2); % d2 = imdilate(e2, B2);
mmf = (e1+d2)/2; % average
% Plot
figure(2); plot(T,mmf,"b-"); grid on;
title("MMF Denoised Ecg signal", "FontSize", fontSize);
xlabel("Time (sec)", "FontSize", fontSize);
ylabel("Voltage (Hz)", "FontSize", fontSize);
Dilatation Function:
function [output] = dilatation(inputSignal,structuringElement)
% Dilatation operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
output = zeros(length(inputSignal));
m = length(structuringElement);
n = length(inputSignal);
for i = abs((m-1)/2):n-abs((m+1)/2)-1
for j=1:m
output(i) = max(output(i-1), inputSignal(i-abs((m-1)/2)+j)+structuringElement(j));
end
end
end
Erosion Function:
function [output] = erosion(inputSignal,structuringElement)
% Erosion operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
output = zeros(length(inputSignal));
m = length(structuringElement);
n = length(inputSignal);
for i = abs((m-1)/2):n-abs((m+1)/2)-1
for j=1:m
output(i) = min(output(i-1), inputSignal(i-abs((m-1)/2)+j)-structuringElement(j));
end
end
end

Accepted Answer

DGM
DGM on 10 Jul 2022
Edited: DGM on 10 Jul 2022
The loop part can be cleaned up, but a big part of why the output looks messed up starts here:
output = zeros(length(inputSignal));
This will create a square matrix, not a vector. So then this assignment inside the loop:
output(i) = min(...);
... will populate the first column. Then if you feed the result to plot(), it plots each row as a series.
So just make sure to preallocate a vector.
I'm going to assume that this is closer to what you're looking for.
insignal = [0 0 0 0 0 1 1 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0];
se = [1 1 1];
outsignal = erosion(insignal,se);
subplot(2,1,1)
plot(insignal)
subplot(2,1,2)
plot(outsignal)
function [output] = erosion(inputSignal,structuringElement)
% Erosion operation on a 1-D signal
% insputSignal = signal in input
% StructuringElement = structuring element in input
% get sizes
M = length(structuringElement);
N = length(inputSignal);
% this relies on the strel being odd width
if mod(M,2)==0
error('strel length must be odd')
end
% perform dilation
hw = floor(M/2); % half-width of strel
output = zeros(1,N);
for n = ceil(M/2):N-ceil(M/2)+1
output(n) = min(inputSignal(n-hw:n+hw) - structuringElement);
end
end
Bear in mind that the behavior at the ends is due to the way the operation is constrained there. The alternative to avoiding the ends of the signal vector is to add enough padding to them that you can filter the whole vector.
  9 Comments

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!