Dilatation and Erosion as custom function for signal processing

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

Answer 1

[DGM]

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.

[Sabaudian]

Thanks you for the clarification.

Just one more thing. What kind of ending behaviour are you referring to?

[DGM]

The output array is initialized with zeros, but the calculated values do not extend all the way to the end of the output vector due to the width of the strel. As a result, the output vector has (because the strel was width 3) one element at each end that still has the original zero in it.

Assume that the strel was width 5. There would be two elements at each end which don't get processed:

The alternative is to pad the array with floor(M/2) extra values at each end. That way you can process all the values in the signal vector. The padding can be done either by simple concatenation, or by using padarray(). By adjusting the loop indices accordingly, the output array size remains the same.

There are other alternatives that can be imagined, and what's appropriate depends on the intent. Just pay attention to how the strel overlaps the signal vector.

[Sabaudian]

I am sorry to bother you again but,

I have tried simply to do something like this (add the padding to the signal before callling the function to erode and dilatate the signal), but does not seems to work

new_signal = [zeros(round(M/2)) my_signal zeros(round(M/2))] 

Even with the padding at the start and at the end of the signal there is this kind of "distortion"

I can't understand what I am supposed to do with the two function (referring to this: "By adjusting the loop indices accordingly, the output array size remains the same.") because if I change the size before use them the indices are just right as they are.

% 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;
% Distort the Baseline aadding drift
X = linspace(0,20*pi,L);
A = 2; % ampl. of upward or downward drift
N = 30; % severity of baseline roll
sn = A*cos(X./N);
minDriftOffset = 1;
maxDriftOffset = 5;
slantedLine = linspace(minDriftOffset,maxDriftOffset,L);
drift = sn + slantedLine;
% Apply drift to all signals
dft_ecg = ecg + drift;
%% Baseline correction
close all;
% Defining the two structuring element
Bo = ones(1,73);
Bc = ones(1,145);
M = length(Bc);

% Padding before conditioning the signal
new_signal = [zeros(round(M/2)) dft_ecg zeros(round(M/2))];

% Detection and Correction of the Wandering Baseline 
peaksSuppression = dilatation(erosion(new_signal,Bo),Bo);
pitsSuppression = erosion(dilatation(peaksSuppression,Bc),Bc);
detectedDrift = pitsSuppression;
% Correction subtracting the drift from signals
Correction = new_signal - detectedDrift;
% Final Baseline result
finalBaseline = erosion(dilatation(dilatation(erosion(Correction,Bo),Bo),Bc),Bc);
% Plot Corrected Signal and Baseline
figure(1); subplot(2,1,1); hold on;
plot(detectedDrift,"g-","LineWidth",3);
plot(new_signal, "b-","LineWidth",0.5);
title("Detected Baseline Drift", "FontSize", fontSize); 
xlabel("Time (sec)", "FontSize", fontSize);  
ylabel("Voltage (Hz)", "FontSize", fontSize);
grid on; hold off;
subplot(2,1,2); hold on;
plot(finalBaseline,"g-","LineWidth",3);
plot(Correction, "b-","LineWidth",0.5);
title("Drift Correction", "FontSize", fontSize); 
xlabel("Time (sec)", "FontSize", fontSize);  
ylabel("Voltage (Hz)", "FontSize", fontSize);
grid on; hold off;

[DGM]

Something like this.

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_padded(insignal,se);
subplot(2,1,1)
plot(insignal)
subplot(2,1,2)
plot(outsignal)

function [output] = erosion_padded(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
    
    % pad vector
    hw = floor(M/2); % half-width of strel
    inputSignal = padarray(inputSignal,[0 hw],0,'both');
    
    % perform dilation
    output = zeros(1,N); % same length as input
    for n = hw+1:hw+N % n = (1:N)+hw (indexes over the same length as the input, but offset by +hw)
        % the indexing on the LHS removes that +hw offset
        output(n-hw) = min(inputSignal(n-hw:n+hw) - structuringElement);
    end
end

Also bear in mind that everything so far assumes that the input signal and strel are row vectors.

[Sabaudian]

Ok, I have tried your solution and it seems to give a better result but not optimal lets say, there is still some part of the signal that is not properly included into the operations (erosion and dilatation), causing a "distortion" in the output at the start and at the end, as it is possible to see in the image below

So to resolve this, what I am supposed to to? Enlarge the Padding? or something else?

[DGM]

If you're still padding prior to calling the erosion/dilation functions, don't. The padding is handled within the functions themselves.

[Sabaudian]

Yes I have already delete that part when I applied your solution. the only padding is in the functions.

The code is the same as before (minus the padding part) and as you can see the problem persist.

% 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;
% Distort the Baseline aadding drift
X = linspace(0,20*pi,L);
A = 2; % ampl. of upward or downward drift
N = 30; % severity of baseline roll
sn = A*cos(X./N);
minDriftOffset = 1;
maxDriftOffset = 5;
slantedLine = linspace(minDriftOffset,maxDriftOffset,L);
drift = sn + slantedLine;
% Apply drift to all signals
dft_ecg = ecg + drift;
%% Baseline correction
close all;
% Defining the two structuring element
Bo = ones(1,73);
Bc = ones(1,145);
% Detection and Correction of the Wandering Baseline 
peaksSuppression = dilatation(erosion(dft_ecg,Bo),Bo);
pitsSuppression = erosion(dilatation(peaksSuppression,Bc),Bc);
detectedDrift = pitsSuppression;
% Correction subtracting the drift from signals
Correction = dft_ecg - detectedDrift;
% Final Baseline result
finalBaseline = erosion(dilatation(dilatation(erosion(Correction,Bo),Bo),Bc),Bc);
% Plot Corrected Signal and Baseline
figure(1); subplot(2,1,1); hold on;
plot(T,detectedDrift,"g-","LineWidth",3);
plot(T,dft_ecg, "b-","LineWidth",0.5);
title("Detected Baseline Drift");
% legend("baseline","signal");
xlabel("Time (sec)");  
ylabel("Voltage (Hz)");
grid on; hold off;
subplot(2,1,2); hold on;
plot(T,finalBaseline,"g-","LineWidth",3);
plot(T,Correction, "b-","LineWidth",0.5);
title("Drift Correction");
% legend("baseline","signal");
xlabel("Time (sec)");  
ylabel("Voltage (Hz)");
grid on; hold off;

[DGM]

Augh. It's picking up the padding anyway. Try doing this instead of padding with zeros:

input = padarray(input,[0 hw],'replicate','both');

I am kind of confused by the behavior of this method of dilation/erosion. I'm used to doing this with images, where the strel is multiplied instead of added/subtracted. It's just not a convention I'm used to thinking about. At least with the above change, it seems to work okay. See what you think.

[Sabaudian]

It works! Thank you so much for the help and the patience