how to separate the wave from ecg image

Attachment is original color image of ECG in png version.

The task is to separate the wave from the rest (removing red dots, black column and row bars, the little black rectangle over the peak of wave).

Ang way to do that? Following is a portion of the ECG image

Answer 1

[Image Analyst]

Try this:

% Demo by Image Analyst to extract an ECG signal from an image or a paper strip chart. 
% Feb. 15, 2022
clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clearvars;
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
fprintf('Beginning to run %s.m ...\n', mfilename);
%-----------------------------------------------------------------------------------------------------------------------------------
% Read in image.
folder = pwd;
baseFileName = 'ecg.png';
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
    % The file doesn't exist -- didn't find it there in that folder.
    % Check the entire search path (other folders) for the file by stripping off the folder.
    fullFileNameOnSearchPath = baseFileName; % No path this time.
    if ~exist(fullFileNameOnSearchPath, 'file')
        % Still didn't find it.  Alert user.
        errorMessage = sprintf('Error: %s does not exist in the search path folders.', fullFileName);
        uiwait(warndlg(errorMessage));
        return;
    end
end
rgbImage = imread(fullFileName);
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display the image.
subplot(3, 2, 1);
imshow(rgbImage, []);
axis('on', 'image');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
caption = sprintf('Original RGB Image : "%s"\n%d rows by %d columns', baseFileName, rows, columns);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
drawnow;
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
hFig1 = gcf;
hFig1.Units = 'Normalized';
hFig1.WindowState = 'maximized';
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
hFig1.Name = 'Demo by Image Analyst';
%-----------------------------------------------------------------------------------------------------------------------------------
% Get the red channel and threshold it at 107.
mask = rgbImage(:, :, 1) < 107;
% Display the image.
subplot(3, 2, 2);
imshow(mask, []);
axis('on', 'image');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
caption = sprintf('Binary Mask Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Find the vertical axis by summing vertically
horizontalProfile = sum(mask, 1);
% Display the image.
subplot(3, 2, 3);
plot(horizontalProfile, 'b-', 'LineWidth', 2)
grid on;
caption = sprintf('Mask Summed Vertically');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Find out which columns have more than 100 pixels in them and erase those from the mask.
badColumns = horizontalProfile > 100;
mask(:, badColumns) = false;
% Now take the largest blob.
mask = bwareafilt(mask, 1);
% Now skeletonize down to a single pixel wide line.
mask = bwskel(mask);
% Display the image.
subplot(3, 2, 4);
imshow(mask, []);
axis('on', 'image');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
caption = sprintf('Final Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
% Get the signal
signal = nan(1, columns);
for col = 1 : columns
    t = find(mask(:, col), 1, "first");
    if ~isempty(t)
        signal(col) = rows - t;
    end
end
% Get rid of any nans
badIndexes = isnan(signal);
x = 1 : columns; % Initialize to entire width.
x(badIndexes) = []; % Remove nans
signal(badIndexes) = [];
% Now we can rescale the signal from units of rows to whatever we want, like to the 0-1 range.
signal = rescale(signal, 0, 1);
% Now plot the final signal.
subplot(3, 2, 5:6);
plot(x, signal, 'b-', 'LineWidth', 3);
grid on;
title('Final Signal', 'FontSize', fontSize);
xlabel('Column', 'FontSize', fontSize);
ylabel('Row', 'FontSize', fontSize);

[Image Analyst]

By the way, can you tell me what the purpose of the red dots is? Why do they put them on the paper anyway? Is it like to provide some sort of grid, and they don't want to use lines?

[lt c]

Great, appreciate your hard work and timely response.

You are right about red dots being grid, which provides unit length for myocardial potential and time interval. I will apply your code in the project. Hope you do not mind that there might be follow-up questions.

The parameters:

'''

mask = rgbImage(:, :, 1) < 107;

badColumns = horizontalProfile > 100;

'''

How could I pick 107 and 100 for thresholds? A few extra steps to come up with these figures?

Could you shed some light on the segmentation of the original ECG. The attached original image has 10 rows of ECG tracing, which is a continuous recording of one patient's heartbeats.

Thank you.

Cheers, LT

[Image Analyst]

I just picked a threshold that worked for that particular image. If your exposure varies all over the place, then you'd have to have a different threshold selection method that takes that into account.

The thresholding looks at the red channel because in the red channel the red dots will be the same brightness as the paper. If we were to look at the blue or green channel, then there would be cark dots where the red dots are and that would just complicate things because we now could not say that all dark stuff is the trace.

The threshold for the y axis was taken by looking at the horizontal profile. Where there is a trace, the vertical sum of dark stuff is a few pixels, while at the y axis, there are lots of dark pixels. So when there are lots of dark pixels in a column, I know it must be the y axis and I can exclude that part.

[lt c]

Nice move, I got it now.

Your approach is way more effective. Back then, I was struggled with an idea of building a clean background as mask (such as using a n*n square of red dots without curve, then multiply by m to fit every area between y-axis), and subtract by the original image, hoping backgrounds cancel each other out and leave the curve.

Looking forward to your doing of segmentation, thank you.

Cheers

[Image Analyst]

Yeah, just taking the red channel is a much easier way to exclude the red dots. So if we're done here, could you please click the "Accept this Answer" link? Thanks in advance. 🙂

[lt c]

Yes, certainly. Before that, I do need your help with the segmentation. Maybe I did not make myself clear.

Your code works excellently on this figure which is only a portion of the original one. As you may notice the original image (following) has multiple lines of signals that we previously worked on.

Since it is a continuous recording, the logical solution being to segment every line of tracing and put them back together in nose to tail, as the red arrow indicated. Does this make sense to you?

The original image was posted as the attachment(s.png). I was trying to change the threshold to 400, and nest a row loop in the "getting signal" section in order to fit the whole image. But I failed. Could you help me with this problem? And how to form these pieces of tracing into one line as it should be, in natural condition.

Thank you.

Cheers, LT

[Image Analyst]

Can you just ask the user how many strips there are and then extract them into subimages? Like if N = 10

N = 10;
[rows, columns, numberOfColorChannels] = size(rgbImage)
startRows = linspace(1, rows+1, N+1);
endRows = startRows(2:end) - 1;
startRows = startRows(1:end-1);
for k = 1 : length(startRows)
    row1 = startRows(k);
    row2 = endRows(k);
    thisImage = rgbImage(row1:row2, :, :);
    % Now process this one image.
end

If N varies, you can determine it by taking the red channel and summing it horizontally:

verticalProfile = mean(rgbImage(:, :, 1), 2);

[lt c]

I kind of get your idea. Basically, I use imcrop to cut the upper bound (index of patient id) out. Then, use linspace to equally divide the image into several pieces (dependent on how many strips). Finally, use startrow and endrow to extract every single strip for further preprocessing, right?

Unfortunately, I had a hard time with the 2nd solution.

verticalProfile = mean(rgbImage(:, :, 1), 2); This code returns a column vector, containing the mean of every row, right? Say there are 11 strips, how does the mean of row help with the N issue? Collecting rows with mean < 255, and see the distribution of row index?

Also, don't mind me asking about the assembling problem. I took two adjacent strips: thisImage(2) and thisImage(3), assembled.

testImage = [thisImage(2) thisImage(3)];

imshow(testImage,[])

These strips do not line up evenly (the curves that blue and red arrows point were supposed to meet).

Also, somehow there is a space between two strips.

Any way to solve these problems? Thanks

Cheers

[Image Analyst]

By looking at the vertical profile you can find out where the signal is lowest, which is where those solid black lines go all the way across the image. If you can't do it that way, just divide it evenly using startRows and endRows.

[lt c]

Quite right. And the number of strips equals to the number of solid black lines minus 1.

Thank you. It is great to have you in this community. I have accepted your answer.

Hope I could learn more from you.