How to sum (row,col) until the first red pixel

Hi

Any ways to use the sum() command in an image, in a way to scan the whole image until it finds the first red pixel where it should stop and proceed to the nect column starting the scanning again.

Any ideas? Thank you.

[Rik]

Your definitions/data structure are not entirely clear. What do you mean with red pixel? Do you have an indexed image or is it an RGB?

Is every row guaranteed to have a red pixel? If not, what should happen in that case?

[Stelios Fanourakis]

@Rik

I attach the Image.

The user needs to define the start and end point on the image (row, columns). The sum() command will start scanning vertically, until it finds the first red pixel on the contour you see in the image. Then it will stop and proceed to the next column.

Answer 1

[Image Analyst]

I believe I answered the image in your duplicate question

Here it is again:

I think this does what you want.

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;

%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'A2.jpg';
folder = pwd;
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

%=======================================================================================
% Read in demo image.
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display image.
subplot(2, 3, 1);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\n%s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0.05 1 0.95]);
% 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.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;

%===================================================================================
% FIND THE RED MASK
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
% get a mask for the red line.
mask = redChannel > greenChannel;
% Display the image.
subplot(2, 3, 2);
imshow(mask);
caption = sprintf('Color Segmentation Mask Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
drawnow;
% fill holes.
mask = imfill(mask, 'holes');
% Take the largest blob only.
mask = bwareafilt(mask, 1);
% Shrink two layers of pixels to come inside the red lines.
mask = imerode(mask, true(5));
% Display the image.
subplot(2, 3, 3);
imshow(mask);
title('Final Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
% Get a masked image.
% Mask the image using bsxfun() function to multiply the mask by each channel individually.
maskedRgbImage = bsxfun(@times, rgbImage, cast(mask, 'like', rgbImage));
% Display the image.
subplot(2, 3, 4);
imshow(maskedRgbImage);
title('Masked RGB Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');

% A = contour(mask,'b');
set(gca,'ydir','reverse')
conversion = 0.9183; % in mm/pixel
addMM = @(x) sprintf('%.3fmm', x * conversion);
xticklabels(cellfun(addMM, num2cell(xticks'), 'UniformOutput', false));
yticklabels(cellfun(addMM, num2cell(yticks'), 'UniformOutput', false));
    
%===================================================================================
% FIND THE HEIGHT OF THE MASK, AND SUM OF GRAY LEVELS FROM ROW 1 TO RED ROW.
topRow = 1; % Whatever you want.
grayLevelSum = zeros(1, columns); % Initialize sums for every column.
heights = zeros(1, columns); % Initialize heights for every column.
% Find the left and right column of the mask so we know where to sum over.
[rows, columns] = find(mask);
col1 = min(columns);
col2 = max(columns);
for col = col1 : col2
    row2 = find(mask(:, col), 1, 'last'); % Get bottom row of mask
    if ~isempty(row2)
        % Get the sum in the GREEN channel, from row1 to row2.
        grayLevelSum(col)  = sum(greenChannel(topRow:row2, col));
		% Compute the height of the mask at every column.
		heights(col) = sum(mask(topRow:row2, col));
    end
end
    
% Plot integrated gray level.
subplot(2, 3, 5);
plot(grayLevelSum, 'b-', 'LineWidth', 2);
xlabel('Column', 'FontSize', 20);
ylabel('Sum of gray levels in column', 'FontSize', 20);
title('Sum of gray levels in column', 'FontSize', 20);
grid on;
% Plot heights
subplot(2, 3, 6);
plot(heights, 'b-', 'LineWidth', 2);
xlabel('Column', 'FontSize', 20);
ylabel('Height of mask in column', 'FontSize', 20);
title('Mask Height vs. Column', 'FontSize', 20);
grid on;

% Tell user the answer.
message = sprintf('Done!');
uiwait(helpdlg(message));

[Stelios Fanourakis]

That's what I call life saving. No words to describe my grace.

Although, I believe the height of the mask is a bit more than it should. The real heigh should be around 10 mm maximum and not exceeding them at all.

In my case, it shows it around 13 mm, but I believe it is attributed to the row I chose, since, it is out of the mask boundaries where the sum begins, thus the overall graph should be a bit bigger.

Thanks a ton sir Analyst

[Image Analyst]

I counted from line 1. Make sure that is what should be the line to start counting downwards from.

I also don't know what your spatial calibration iis in mm per pixel. Maybe it's not right.

[Stelios Fanourakis]

Actually, I managed to get the difference using that code.

open('A2.fig');
a = get(gca,'Children');
xdata = get(a, 'XData');
ydata = get(a, 'YData');
zdata = get(a, 'ZData');

diff = heights - ydata

I attach both picture plots and its plot of difference. I would assume the difference plot (example.jpg) to be around 2-3 mm as it is obvious a difference among the other two plots (second,second2.jpg). But I struggle see the difference of 2-3mm depicted.

[Stelios Fanourakis]

Image Analyst. I used your perfect code to do my job and it worked for the first image ('A1.jpg') I attach.

Consequently, I used a small extra script to export the columns and row pixel values in xls file format. For first picture it worked.

When I try to do the same thing, to the second image ('A1b.jpg'), it is not working. The exporting pixel values in an xls column is only 0. Also the graphs are not correct.

Can you please tell me what is wrong?

[Stelios Fanourakis]

Also, see attached image (question.jpg), can you please tell me what is actually this?

I use this script to export the heights to an xls column

data = [heights(:)];
sheet = 1;
range = 'B';
xlswrite('YourFile.xls', data, sheet, range)

Although, this time I get some values, they are not correct. They should be more or less the same as the previous image but they are indeed awkward. For instance for column 136 I should take 83 pixels and Matlab gives me only 4 which is not correct.

I believe it has to do something with the image I attach and only the upper part with those jaggy spikes is calculated for heights and not the whole image area where it should.

[Image Analyst]

Looks like the red outline is not in the same place? Why not? Can you make sure it has the red in the same location as the other image? Who or what is burning the red line into the image.

Do not use JPG for image analysis since it blurs the red line.

[Stelios Fanourakis]

Actually, it is done manually. An operator segments with red vectors both images. At A1.jpg he closes the whole area making a closed contour. At the A1b.jpg image he just draws the bottom line. He doesn't make the whole area. But both bottom red lines are placed at the same position in both images. Maybe there are 1-2 pixels difference at some points. But 98-99% they are identical, at least for the bottom line.

What format should I export images?

[Image Analyst]

The images should be saved as PNG (or BMP or TIFF), but never JPG.

Why does the user only draw the top of the region in some of the images, but not all of them? Sees like a user training issue to me. Just tell them to draw the top of the region. Do they not want to for some reason? It can be done to do it both ways but you have to first check to see where the top red line is. If it's on line 1 and flat, do ti one way. If it's on some other line and flat, then do it a different way. And if the top line is not flat at all, I guess do it a third way. Do you really want to spend the time to make your program so flexible when it's probably no big deal to your users just to tell them where to draw the top? They'll probably happily do it.

[Stelios Fanourakis]

Actually, there is a reason for users NOT to draw the top line. It will be more flexible to have just the bottom line and that's because this algorithm is supposed to implemented in different ultrasound units, which means, that the penetration depth of the beam will not be the same. We care only for the bottom line because it is the interface where skin and a standoff material meet. If the user draws the top line we limit the algorithm to trace for a certain depth. But that will not be the case for all ultrasound units.

So, NO it is not a training issue. We want for the user, to define by choice the row and column for the sum scanning to start no matter the height, cause the depth obviously will be different. But we care to stop where the bottom red line is, since it is the interface we care, and provide the sum of pixels.

Can you modify the algorithm please? Obviously I don't get the right values. There should be slighlty differences among the two images, but they are way too much irrelevant.

[Image Analyst]

I'm heading out for most of the day. While I'm out, tell me how the users area drawing the line. Are they using imfreehand? Some other method? Can we at least assume that the boundary will be the lower of the coordinates? In other words, the flat, non-boundary part of the red (that closes off the region) won't be along the bottom of the image or below the interface.

[Stelios Fanourakis]

The users use any vector software to draw the red line. There won't be any straight upper part of it. Only the lower part which is in 'A1b.jpg' I attach

[Stelios Fanourakis]

The lower red line at the middle part of the image. Discard any red straight lines that might appear up or side of the image. We only care about the imaging window.

Any help?

[Image Analyst]

Try this:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 16;
%===============================================================================
% Get the name of the image the user wants to use.
baseFileName = 'A1b.jpg';
folder = pwd;
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
%=======================================================================================
% Read in demo image.
rgbImage = imread(fullFileName);
% Get the dimensions of the image.
[rows, columns, numberOfColorChannels] = size(rgbImage)
% Display image.
subplot(2, 2, 1);
imshow(rgbImage, []);
impixelinfo;
axis on;
caption = sprintf('Original Color Image\n%s', baseFileName);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
hp = impixelinfo(); % Set up status line to see values when you mouse over the image.
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0.05 1 0.95]);
% 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.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
drawnow;
%===================================================================================
% Mask out everything except the signal.
mask = rgbImage(:, :, 2) > 0;
% Get the largest region
mask = bwareafilt(mask, 1);
mask = imfill(mask, 'holes'); % Fill any holes.
% Mask the image using bsxfun() function to multiply the mask by each channel individually.
rgbImage = bsxfun(@times, rgbImage, cast(mask, 'like', rgbImage));
% imshow(rgbImage, []);
% Find the topmost row
[maskRows, maskColumns] = find(mask);
topRow = min(maskRows);
% Draw a yellow line there
line([1, columns], [topRow, topRow], 'Color', 'y', 'LineWidth', 2);

%===================================================================================
% FIND THE RED MASK
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
% get a mask for the red line.
mask = redChannel > greenChannel;
% Display the image.
subplot(2, 2, 2);
imshow(mask);
% Draw a yellow line there
line([1, columns], [topRow, topRow], 'Color', 'y', 'LineWidth', 2);
caption = sprintf('Mask Image');
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
axis('on', 'image');
drawnow;
    
%===================================================================================
% FIND THE HEIGHT OF THE MASK, AND SUM OF GRAY LEVELS FROM TOP ROW TO RED ROW.
grayLevelSum = zeros(1, columns); % Initialize sums for every column.
heights = zeros(1, columns); % Initialize heights for every column.
% Find the left and right column of the mask so we know where to sum over.
[rows, columns] = find(mask);
col1 = min(columns);
col2 = max(columns);
for col = col1 : col2
    row2 = find(mask(:, col), 1, 'last'); % Get bottom row of mask
    if ~isempty(row2)
        % Get the sum in the GREEN channel, from row1 to row2.
        grayLevelSum(col)  = sum(greenChannel(topRow:row2, col));
		% Compute the height of the mask at every column.
		heights(col) = row2 - topRow;
    end
end
    
% Plot integrated gray level.
subplot(2, 2, 3);
plot(grayLevelSum, 'b-', 'LineWidth', 2);
xlabel('Column', 'FontSize', 20);
ylabel('Sum of gray levels in column', 'FontSize', 20);
title('Sum of gray levels in column', 'FontSize', 20);
grid on;
% Plot heights
subplot(2, 2, 4);
plot(heights, 'b-', 'LineWidth', 2);
xlabel('Column', 'FontSize', 20);
ylabel('Height of mask in column', 'FontSize', 20);
title('Mask Height vs. Column', 'FontSize', 20);
grid on;
% Tell user the answer.
message = sprintf('Done!');
uiwait(helpdlg(message));

[Stelios Fanourakis]

Under a first evaluation, it seems to work correctly and thank you very much. I am just trying to understand the changes that have been done.

The "heights" subtraction gives the difference between the mask and the topRow which is manually can be defined by the user. Right? And I get the pixel values between this area. Am I right?

I think it's brilliantly genious.

[Image Analyst]

The topRow is computed as the top row of all that signal stuff, after the annotations have been masked out. If you don't like that row and want a different one, feel free to redefine it to be some other row. The sum of the pixel values is computed, and the number of pixels are computed -- both are per row, so they're arrays where the index is the column number. I guess you could say the number of pixels in a column is the area of that column. If you sum them all, you'd get the area of the signal region.

[Stelios Fanourakis]

Yes and thanks a lot. One more question is how do I make it exclude certain excel rows when storing data.

For example, I want to store the data from B2 (second row) as I want the first row as a headline. How do I do this programmatically?

I use that code. I want to adjust the range

data = [heights(:)];
sheet = 1;
range = 'B';
xlswrite('YourFile.xls', data, sheet, range)

[Image Analyst]

Make a cell array for your header, then call xlswrite once for the header and once for the numerical data

header = cell(1, 2); % Row vector
header{1} = 'Column 1';
header{2} = 'Column 2';
sheet = 1;
xlswrite('YourFile.xls', header, sheet, 'A1')
xlswrite('YourFile.xls', data, sheet, 'A2') % Write one row below the header.