distance measurement in image based on intensity difference

Hi,

I have a serious of images similar to the one attached here. I intend to measure the distance of particular portion of the image from the bottom as shown as yellow arrow lines in the attached figure.. If somebody came across this kind of problem, please shed some light on this

Answer 1

[Image Analyst]

Try this. If it works you can "Accept this answer" and you might want to consider going over your prior posts and "Accepting" some of them to let people know they're solved and to give the people who helped you "reputation points.":

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 = 22;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'B00015.jpg';
% Get the full filename, with path prepended.
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
grayImage = imread(fullFileName);


% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Display histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of original gray image', 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION OF IMAGE
threshold = 60;
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
mask = bwareafilt(mask, 1);
% Display the image.
subplot(2, 2, 3);
imshow(mask, []);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Measure Bounding Box.
props = regionprops(mask, 'BoundingBox');
% Display the original image again so we can overlay graphics on it.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Overlays', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Now show the bounding box
hold on;
rectangle('Position', props.BoundingBox, 'LineWidth', 2, 'EdgeColor', 'r');
% Find the last row
lastRow = floor(props.BoundingBox(2) + props.BoundingBox(4))
% Compute the distance between the last row and the bottom of the image.
distance = rows - lastRow	% Add 1 if you want.
% Find out all the columns that are in the blob.
[blobRows, blobColumns] = find(mask);
lastRow = max(blobRows) % Should be the same as above.
% Find out columns where the blob is in the last row
lastColumns = find(mask(lastRow, :))
% Get the mean of these to get the mean x (column) value.
meanX = mean(lastColumns);
% Draw a yellow line up to there.
line([meanX, meanX], [columns, lastRow], 'Color', 'y', 'LineWidth', 2);
% Put up text next to the line half way up.
textLabel = sprintf('  Last Row at %d.\n  Distance = %d', lastRow, distance);
yt = (rows + lastRow)/2
text(meanX, yt, textLabel, 'FontSize', 13, 'FontWeight', 'bold', ...
	'Color', 'y', 'VerticalAlignment', 'middle');

[Turbulence Analysis]

Hi Image Analyst,

Thank you so much for your detailed response. It's worked well for me..

Actually, I have to implement in this in set of files, rather than one file. In that case, obvisiouly i have to choose optimal thresold value for all the images.. Besides, to deal with multiple images present inside the folder, I tried as follows, but got the error msg

folder = pwd;

baseFileName = dir('*.jpg');

% Get the full filename, with path prepended.

fullFileName = fullfile(folder, baseFileName);

Error using fullfile (line 67)

All inputs must be strings, character vectors, or cell arrays of character vectors.

[Image Analyst]

The optimal threshold should be chosen somehow but then it should be fixed for all images. As long as all your images are taken under the same conditions, a fixed threshold is the only way to make sure your mask will detect everything from 0% coverage to 100% coverage and be consistent from one snapshot to the next. You certainly don't want an automatic threshold. The belief (usually by novices) that an automatic threshold is the best approach is a misinformed myth. Just think about it some and you'll realize I'm right. OK, let's take an example. Let's say you have a light spray or cloud (whatever it is) and you pick a low threshold automatically so it gets all the light spray. Now let's say the next image has a heavy white cloud and you pick an automatic threshold, and it picks a higher threshold. Now the lighter spready that was surrounding the heavy white spray will not get detected, whereas before it was detected. Now let's say you have an image where it's mostly white (almost completely covered) or all black (mostly uncovered). You'd be forcing it to choose a threshold that would split the histogram in two (foregound and background) and that threshold would like in the middle of the narrow histogram. So it might say the covereage was around 50% whereas in reality it was 0% or 100%. So an automatic threshold will not be consistently detecting the same levels of spray density across multiple images.

To process a sequence of files, see the FAQ: https://matlab.fandom.com/wiki/FAQ#How_can_I_process_a_sequence_of_files.3F

This is not correct:

baseFileName = dir('*.jpg')

look up the help for dir() and you can see it does not produce a base file name. It returns a structure of files - a file listing. To get the k'th base file name and full name (in a loop over k) you should do

fileList = dir('*.jpg');
for k = 1 : length(fileList)
    baseFileName = fileList(k).name;
    fullFileName = fullfile(fileList(k).folder, baseFileName);
    grayImage = imread(fullFileName);
    % etc.
end

Finally, don't use JPG format for imaging if you can help it. It introduces JPEG compression artifacts. Use uncompressed or lossless compression formats like PNG, TIFF, or BMP.

[Turbulence Analysis]

Thanks for an excellent insights on selection of optimal thresold..

for getting the distance for the multiple images , I have added a below line inside for loop

D1(k) = distance;

and outside the For loop I used

D1 = zeros (1, length(fileList));

Its working well, however, beyond some number of images, it throws below error..

Error using horzcat
Out of memory. Type HELP MEMORY for your options.
Error in regionprops>ComputePixelList (line 913)
            stats(k).PixelList = [In{:}];
Error in regionprops>ComputeBoundingBox (line 648)
        ComputePixelList(imageSize,stats,statsAlreadyComputed);
Error in regionprops (line 257)
                ComputeBoundingBox(imageSize,stats,statsAlreadyComputed); 

[Image Analyst]

I'd need more info than that. Presumably you called regionprops() but you're only showing the errors in regionprops(). You chopped off the later part of the error message traceback that shows where (line number) and how you called regionprops() from your program. Attach your m-file so I can fix it.

[Turbulence Analysis]

Hi,

It's working normal, I guess one of the image inside the folder has got some issue.. That's why while running in a loop it throws some error...

Attached you the code I am using with loop..

Besides, I am having some images (show in the attachements), here I need to find the distance on both left and right sides.. So, for this do i need to crop images into two halves and execute the code to get the distance on both the sides ???

Thanks in advance for your time and kind support..

[Image Analyst]

Turb: You just need to make a few modifications to be able to handle multiple blobs. I'm attaching that 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 = 22;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'original img.jpg';
% Get the full filename, with path prepended.
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
grayImage = imread(fullFileName);


% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Display histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of original gray image', 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION OF IMAGE
threshold = 60;
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
% Extract only blobs with an area of 100 or more.
mask = bwareafilt(mask, [100, inf], 4); % Connectivity of 4
% Apply a variety of pseudo-colors to the regions.
[labeledImage, numberOfBlobs] = bwlabel(mask, 4);
coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
% Display the pseudo-colored image.
subplot(2, 2, 3);
imshow(coloredLabelsImage);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Measure Bounding Boxes.
props = regionprops(labeledImage, 'Area', 'BoundingBox', 'PixelList');
caption = sprintf('Mask with %d blobs', length(props));
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
allAreas = [props.Area]
% Display the original image again so we can overlay graphics on it.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Overlays', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

for k = 1 : length(props)	
	% Now show the bounding box
	hold on;
	rectangle('Position', props(k).BoundingBox, 'LineWidth', 2, 'EdgeColor', 'r');
	% Find the last row
	lastRow(k) = floor(props(k).BoundingBox(2) + props(k).BoundingBox(4));
	% Compute the distance between the last row and the bottom of the image.
	distances(k) = rows - lastRow(k);	% Add 1 if you want.
	% Find out all the columns that are in the blob.
	thisBlobsMaskXYCoordinates = props(k).PixelList;
	lastRow(k) = max(thisBlobsMaskXYCoordinates(:, 2));
	% Find out columns where the blob is in the last row.
	% First get the indexes (rows in the thisBlobsMaskXYCoordinates matrix) where y (column 2) equals the last row.
	yIndexes = find(thisBlobsMaskXYCoordinates(:, 2) == lastRow(k));
	% Get the mean of these to get the mean x (column) value.
	meanX(k) = mean(thisBlobsMaskXYCoordinates(yIndexes, 1));
	% Draw a yellow line up to there.
	line([meanX(k), meanX(k)], [columns, lastRow(k)], 'Color', 'y', 'LineWidth', 2);
	% Put up text next to the line half way up.
	textLabel = sprintf('Blob #%d : Last Row at %d.\n          Distance from bottom = %d', k, lastRow(k), distances(k));
	fprintf('%s\n', textLabel); % Print to command window.
	yt = (rows + lastRow(k))/2;
	text(meanX(k), yt, textLabel, 'FontSize', 8, 'FontWeight', 'bold', ...
		'Color', 'y', 'VerticalAlignment', 'middle');
end

[Turbulence Analysis]

Thanks a lot image analyst... Your contribution to the image prcessing community is amazing...

[Turbulence Analysis]

Dear Image analyst,

In refrence to our earlier conversation, I have noticed some issues while processing the images.. For instant here I have attached three raw and processed images. As you see, the distance is correctly detected in image 10, however, it detects in wrong place for image 8, 9.. Actually teh idea is to locate the first white portion from the bottom of the images.. I played with the thresold values, but not successful.

Could you please shed some light on this....

[Image Analyst]

You did something to change it, but you did not tell me what it was. You forgot to attach your code. My code detects all 13 blobs and finds the lowest point in each one. If you want just the lowest one instead of all of them, let me know. That's fairly easy to put in but we first need to figure out why your code is detecting only one, instead of all of them, and why that one is not the lowest one.

[Turbulence Analysis]

Dear Image analyst,

Thanks for your reply. I have attached the code below. Actually this is the first code you suggested (i.e. detection of only one position). Yes, I am intrested in the lowest one..

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 = 22;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'B00008.bmp';
% Get the full filename, with path prepended.
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
grayImage = imread(fullFileName);


% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Display histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of original gray image', 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION OF IMAGE
threshold = 75;
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
mask = bwareafilt(mask, 1);
% Display the image.
subplot(2, 2, 3);
imshow(mask, []);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Measure Bounding Box.
props = regionprops(mask, 'BoundingBox');
% Display the original image again so we can overlay graphics on it.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Overlays', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Now show the bounding box
hold on;
rectangle('Position', props.BoundingBox, 'LineWidth', 2, 'EdgeColor', 'r');
% Find the last row
lastRow = floor(props.BoundingBox(2) + props.BoundingBox(4))
% Compute the distance between the last row and the bottom of the image.
distance = rows - lastRow	% Add 1 if you want.
% Find out all the columns that are in the blob.
[blobRows, blobColumns] = find(mask);
lastRow = max(blobRows) % Should be the same as above.
% Find out columns where the blob is in the last row
lastColumns = find(mask(lastRow, :))
% Get the mean of these to get the mean x (column) value.
meanX = mean(lastColumns);
% Draw a yellow line up to there.
line([meanX, meanX], [columns, lastRow], 'Color', 'y', 'LineWidth', 2);
% Put up text next to the line half way up.
textLabel = sprintf('  Last Row at %d.\n  Distance = %d', lastRow, distance);
yt = (rows + lastRow)/2
text(meanX, yt, textLabel, 'FontSize', 13, 'FontWeight', 'bold', ...
	'Color', 'y', 'VerticalAlignment', 'middle');

[Image Analyst]

I'm sure you got it figured out by now (in the last 5 hours) since it was only a trivial 5 lines of code, but for what it's worth, here is what I did.

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 = 22;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
folder = pwd;
baseFileName = 'b00008.bmp';
% Get the full filename, with path prepended.
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
grayImage = imread(fullFileName);


% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Display histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of original gray image', 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION OF IMAGE
threshold = 60;
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
% Extract only blobs with an area of 100 or more.
mask = bwareafilt(mask, [100, inf], 4); % Connectivity of 4
% Apply a variety of pseudo-colors to the regions.
[labeledImage, numberOfBlobs] = bwlabel(mask, 4);
coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
% Display the pseudo-colored image.
subplot(2, 2, 3);
imshow(coloredLabelsImage);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Measure Bounding Boxes.
props = regionprops(labeledImage, 'Area', 'BoundingBox', 'PixelList');
caption = sprintf('Mask with %d blobs', length(props));
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
allAreas = [props.Area]
% Display the original image again so we can overlay graphics on it.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Overlays', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Find the last row of any of the blobs:
[r, c] = find(mask);
finalRow = max(r);

for k = 1 : length(props)	
	% Now show the bounding box
	hold on;
	rectangle('Position', props(k).BoundingBox, 'LineWidth', 2, 'EdgeColor', 'r');
	% Find the last row
	lastRow(k) = floor(props(k).BoundingBox(2) + props(k).BoundingBox(4));
	% If this one's last row is not the final row, skip it.
	if lastRow(k) < finalRow
		continue; % Skip this blob.
	end
	% Compute the distance between the last row and the bottom of the image.
	distances(k) = rows - lastRow(k);	% Add 1 if you want.
	% Find out all the columns that are in the blob.
	thisBlobsMaskXYCoordinates = props(k).PixelList;
	lastRow(k) = max(thisBlobsMaskXYCoordinates(:, 2));
	% Find out columns where the blob is in the last row.
	% First get the indexes (rows in the thisBlobsMaskXYCoordinates matrix) where y (column 2) equals the last row.
	yIndexes = find(thisBlobsMaskXYCoordinates(:, 2) == lastRow(k));
	% Get the mean of these to get the mean x (column) value.
	meanX(k) = mean(thisBlobsMaskXYCoordinates(yIndexes, 1));
	% Draw a yellow line up to there.
	line([meanX(k), meanX(k)], [columns, lastRow(k)], 'Color', 'y', 'LineWidth', 2);
	% Put up text next to the line half way up.
	textLabel = sprintf('Blob #%d : Last Row at %d.\n          Distance from bottom = %d', k, lastRow(k), distances(k));
	fprintf('%s\n', textLabel); % Print to command window.
	yt = (rows + lastRow(k))/2;
	text(meanX(k), yt, textLabel, 'FontSize', 8, 'FontWeight', 'bold', ...
		'Color', 'y', 'VerticalAlignment', 'middle');
end

[Turbulence Analysis]

It's amazing. I have tried with series of images in the directory, with the below code. This directory has 5 images inside, but in ''distances'' matrix I got only 1 x 8 values. I guess for each image, the size of the ''distances'' matrix depends on the no. of blobs. For instance, if the image has five blobs, it generates 1 x5 array with last column contains the distance of last row. Hence,here, just to keep only non zero entry i.e. last blobs distance I used distances( :, all(~distances,1) ) = [ ], but I am not succesful..

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 = 22;

%--------------------------------------------------------------------------------------------------------
%    READ IN IMAGE
fileList = dir('*.bmp');
% Get the full filename, with path prepended.
% D1 = zeros (1, length(fileList));
% D2=D1';
for kk = 1 : length(fileList)
    baseFileName = fileList(kk).name;
    fullFileName = fullfile(fileList(kk).folder, baseFileName);
    grayImage = imread(fullFileName);
    % etc.
% 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

grayImage = imread(fullFileName);


% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
% Display the image.
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

% Display histogram
subplot(2, 2, 2);
imhist(grayImage);
grid on;
title('Histogram of original gray image', 'FontSize', fontSize);

%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION OF IMAGE
threshold = 60;
mask = grayImage > threshold;
mask = imfill(mask, 'holes');
% Extract only blobs with an area of 100 or more.
mask = bwareafilt(mask, [100, inf], 4); % Connectivity of 4
% Apply a variety of pseudo-colors to the regions.
[labeledImage, numberOfBlobs] = bwlabel(mask, 4);
coloredLabelsImage = label2rgb (labeledImage, 'hsv', 'k', 'shuffle');
% Display the pseudo-colored image.
subplot(2, 2, 3);
imshow(coloredLabelsImage);
title('Mask', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Measure Bounding Boxes.
props = regionprops(labeledImage, 'Area', 'BoundingBox', 'PixelList');
caption = sprintf('Mask with %d blobs', length(props));
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');
allAreas = [props.Area]
% Display the original image again so we can overlay graphics on it.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Overlays', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;

% Find the last row of any of the blobs:
[r, c] = find(mask);
finalRow = max(r);
for k = 1 : length(props)	
	% Now show the bounding box
	hold on;
	rectangle('Position', props(k).BoundingBox, 'LineWidth', 2, 'EdgeColor', 'r');
	% Find the last row
	lastRow(k) = floor(props(k).BoundingBox(2) + props(k).BoundingBox(4));
	% If this one's last row is not the final row, skip it.
	if lastRow(k) < finalRow
		continue; % Skip this blob.
	end
	% Compute the distance between the last row and the bottom of the image.
	distances(k) = rows - lastRow(k);	% Add 1 if you want.
%     distances( :, all(~distances,1) ) = [];
%     D1(kk) = distances;
	% Find out all the columns that are in the blob.
	thisBlobsMaskXYCoordinates = props(k).PixelList;
	lastRow(k) = max(thisBlobsMaskXYCoordinates(:, 2));
	% Find out columns where the blob is in the last row.
	% First get the indexes (rows in the thisBlobsMaskXYCoordinates matrix) where y (column 2) equals the last row.
	yIndexes = find(thisBlobsMaskXYCoordinates(:, 2) == lastRow(k));
	% Get the mean of these to get the mean x (column) value.
	meanX(k) = mean(thisBlobsMaskXYCoordinates(yIndexes, 1));
	% Draw a yellow line up to there.
	line([meanX(k), meanX(k)], [columns, lastRow(k)], 'Color', 'y', 'LineWidth', 2);
	% Put up text next to the line half way up.
	textLabel = sprintf('Blob #%d : Last Row at %d.\n          Distance from bottom = %d', k, lastRow(k), distances(k));
	fprintf('%s\n', textLabel); % Print to command window.
	yt = (rows + lastRow(k))/2;
	text(meanX(k), yt, textLabel, 'FontSize', 8, 'FontWeight', 'bold', ...
		'Color', 'y', 'VerticalAlignment', 'middle');
end
end

[Image Analyst]

You need to index yt, like

yt(kk) = lastRow(k);

or whatever you want to store.

[Turbulence Analysis]

Thanks, I have used distances(kk) = rows - lastRow(k), now for 1000 images I am getting exactly 1 x 1000 array.. Thanks again for your timely support..

[Turbulence Analysis]

Hi,

Is there a way in which I can calculate surface area of the white structure shown in the attached image.

Thanks in advnace for time and support.

[Image Analyst]

Why is the surface area not the Area returned by regionprops()???

props = regionprops(mask, 'BoundingBox', 'Area');

Do you have some other definition?

[Turbulence Analysis]

Hi,

Thanks.. I thought of this. As this surface posessess irregular shape, hence, the bounding box may not yields the correct value. I will look into this later on..

I intend to calculate the distance of discontinuity as shown as yellow arrow arrow marks in the attached figures 1,2. Please help me with this..

[Image Analyst]

The area is the area of the irregularly-shaped region, not of the bounding box.

It looks like the arrow goes to the top of the region. The last (fourth) element of the bounding box is the height of the box and region from bottom to top.

[Turbulence Analysis]

Yes, youe are correct it returns the surface area.

Regarding the distance, as you pointed out, the last values in the bounding box has the height information. And this is saved in the props struct file, may i know how to read this value. I have attached props struct file for your refrence.

Furthermore, I tried some edge detection on the single image

Imageinput = imread('B00002.bmp');

ED = edge(Imageinput ,'sobel');

It pushes below error message

Error using edge

Expected input number 1, I, to be two-dimensional.

Error in edge>parse_inputs (line 482)

validateattributes(I,{'numeric','logical'},{'real','nonsparse','2d'},mfilename,'I',1);

Error in edge (line 213)

[a,method,thresh,sigma,thinning,H,kx,ky] = parse_inputs(args{:});

[Image Analyst]

To get the height of the second blob in the image

height = props(2).BoundingBox(4);

Your image is not gray scale - it's 3-D RGB color, even if it looks gray scale to you. You should do this to check the number of color channels, like I showed you before:

grayImage = imread('B00002.bmp');
% Get the dimensions of the image.
% numberOfColorChannels should be = 1 for a gray scale image, and 3 for an RGB color image.
[rows, columns, numberOfColorChannels] = size(grayImage);
if numberOfColorChannels > 1
	% It's not really gray scale like we expected - it's color.
	% Use weighted sum of ALL channels to create a gray scale image.
	grayImage = rgb2gray(grayImage);
	% ALTERNATE METHOD: Convert it to gray scale by taking only the green channel,
	% which in a typical snapshot will be the least noisy channel.
	% grayImage = grayImage(:, :, 2); % Take green channel.
end
ED = edge(grayImage, 'sobel');

[Turbulence Analysis]

Thanks a lot.

For edge detection, I have tried imerode () function as well, I got the attaced image. I am wondering is there any any way to keep only outer boundary....

[Image Analyst]

I don't even know why you're trying edge detection in the first place. It does not look like it's necessary or appropriate. Many/most beginners think that because you can see edges in the image you need to do edge detection. That is a false assumption in most cases. In your case I would not do it, unless you can expain to me why you think it's necessary. I think just using a threshold, probably a constant/fixed threshold, would be best.

[Turbulence Analysis]

Sorry, actually I thought of extracting outer profile, this i need for some comparision with other process. Anyway, I will play with the thresolding as you said. I am having the images from two cameras, I have to first, superimpose the images and make the videos then. Could you please help me with this.. Below is my code for superimposing just two images....

%first camera image

I=loadvec('Image2.im7');
% Non dimenionalize axis
x=(I.x);
y=(I.y);
F =imagesc(x,y,(I.w'));
ax=gca;
set(gca, 'YDir','normal')
colormap(jet);
F.AlphaData = 0.5;
c1=colorbar;
axis equal
hold on;


% second camera image
I1=loadvec('strain1.im7');
% Non dimenionalize axis
x=(I1.x);
y=(I1.y);
G=imagesc(x,y,(I1.w'));
ax=gca;
set(gca, 'YDir','normal')
colormap(jet);
G.AlphaData = 0.5;
axis equal

[Image Analyst]

You can just average them if you want

meanImage = uint8((double(I) + double(I1)) / 2);

[Turbulence Analysis]

Thanks. I have fixed the issue..

[Turbulence Analysis]

Hi,

I have superimposed two images, could you please let me know, how to have different contour color maps for two images ??. For instance, here I defined 'bone' and 'jet' for first , second images respcctively, but still I end up with 'jet' map for both the images.

% First image
I=loadvec('Image1.im7');
x=(I.x);
y=(I.y);
h = figure (2);
G= imagesc(x,y,(I.w'));
ax=gca;
G.AlphaData = 0.8;
colormap(bone);
set(gca, 'YDir','normal')
axis equal
hold on;

% Second image
V=loadvec('VectorField1.vc7');
x=(V.x);
y=(V.y);
curl = vec2scal(V,'rot');
% S = showf(curl);
vort=curl.w';
HHH = imagesc(x,y,vort);
HHH.AlphaData = 0.5;
colormap(jet);
set(gca, 'YDir','normal')
axis equal

[Image Analyst]

What I'd do it to convert both indexed images to RGB images with ind2rgb() with the colormap you want. Then average the two RGB images together.

[Turbulence Analysis]

Thanks.. Actually my image here is displayed by imagesc(), May I know how to implement the ind2grb () funcction over the imagesc....

[Image Analyst]

Something like this:

colorMap1 = jet(256);
colorMap2 = bone(256);
rgb1 = rgb2ind(grayImage1, colorMap1);
rgb2 = rgb2ind(grayImage2, colorMap2);
rgbAverage = uint8((double(rgb1) + double(rgb2)) / 2;)

Obviously replace the variable names with whatever you have or want.

[Turbulence Analysis]

Thanks, Actually I need to retain the actual intensity values in the original image, hence, I am using .im7 fies instead of image files. I tried your above suggestion with the images instead of .im7 files. But I couldn't get the actual intensity values (capture 2) as shown in image accquired from .im7 files (capture 1)

%%%% code for running with images..
I1 = imread ('B00003.bmp');
%imshow (I1);
I2 = imread ('B00045.bmp');
% I3 = ind2rgb(I1);
% I4 = ind2rgb(I2);
colorMap1 = jet(256);
colorMap2 = bone(256);
rgb1 = rgb2ind(I1, colorMap1);
rgb2 = rgb2ind(I2, colorMap2);
hhhh = imshow (rgb1);
hhhh.AlphaData = 0.5;
hold on
gggg=imshow (rgb2);
gggg.AlphaData = 0.5;
rgbAverage = uint8((double(rgb1) + double(rgb2)) / 2;);

% Code for running with .im7 files 

I=loadvec('B00003.im7');
x=(I.x);
y=(I.y);
h = figure (2);
colormap(jet);
G = imagesc(x,y,(I.w'));
ax1=gca;
G.AlphaData = 0.8;
set(gca, 'YDir','normal')
axis equal
hold on;

I=loadvec('B00045.im7');
x=(I.x);
y=(I.y);
colormap(hot);
GH = imagesc(x,y,(I.w'));
ax1=gca;
GH.AlphaData = 0.5;
set(gca, 'YDir','normal')
axis equal
hold on;

[Image Analyst]

Please reply and attach loadvec.m and the .im7 files again.

[Turbulence Analysis]

Hi, Here it is ....

Please let me know if you face any issues while reading .im7 files.

https://www.dropbox.com/sh/lqqy4twpnya1o7l/AADIMQW0qFuxYDxnwViFNZw6a?dl=0

[Turbulence Analysis]

Hi,

Hope you have recived the files in proper format.

[Turbulence Analysis]

Hi,

This in continuation with our previous discussions on having the two color bars during the image superimposition.

Now, I able to display the two images with different colormap by converting the image file read from .im7 to rgb. However, I am not succesful in showing two colorbars.

For instance, in the attached file, one image is diplayed with 'hot' map and other with 'jet'. But, at the end, both the colorbar switches to 'hot' pattern. Could you please hlep me with this ??

%%% First image
I=readimx('B00003.im7');
h = (I.Frames{1}.Components{1}.Planes {1,1})';
h1 =fix(h/10);
h1(h1==0) = 1 ;
rgb1 = ind2rgb(h1,jet(256));
G = imagesc(flipud(rgb1));
ax1=gca;
colormap(ax1,'jet');
c = colorbar(ax1,'Location', 'east','FontSize',20,'TickLabelInterpreter', 'latex');
c.Title.Interpreter = 'latex';
c.Title.FontSize = 16;
G.AlphaData = 0.5;
set(gca, 'YDir','normal')
axis equal
hold on;
linkaxes([ax1,ax2]);
%%%% second image
I1=readimx('B00045.im7');
h2 = (I1.Frames{2,1}.Components{1,1}.Planes {1,1})';
h3 =fix(h2/10);
h3(h3==0) = 1 ;
rgb2 = ind2rgb(h3,hot(256));
GH = imagesc(flipud(rgb2));
ax2=gca;
colormap(ax2,'hot')
c1 = colorbar(ax2,'Location', 'west','FontSize',20,'TickLabelInterpreter', 'latex');
GH.AlphaData = 0.5;
set(gca, 'YDir','normal')
axis equal

Answer 2

[jonas]

I would first binarize the image and then go columnwise to find the last white pixel.

RGB = imread('B00001.jpg');
GRAY = rgb2gray(RGB);
BW = imbinarize(GRAY);
imshow(BW);hold on
y = nan(size(BW,2),1);
for i = 1:size(BW,2)
    lastcol = find(BW(:,i) == true,1,'last');
    if ~isempty(lastcol)
        y(i) = lastcol;
    end
end
x = 1:size(BW,2);
plot(x,y,'r','linewidth',2);

[Turbulence Analysis]

Hi Jonas,

Thanks for your kind support.. I tried with your code, got the attached plot.. Here, y reprsents the intensity.. I couldn't able to see the binarized image, while executing, imshow(BW), I am getting the below error..

Error using images.internal.imageDisplayValidateParams>validateCData (line 119)
If input is logical (binary), it must be two-dimensional.

Error in images.internal.imageDisplayValidateParams (line 27)
common_args.CData = validateCData(common_args.CData,image_type);

Error in images.internal.imageDisplayParseInputs (line 78)
common_args = images.internal.imageDisplayValidateParams(common_args);

Error in imshow (line 245)
    images.internal.imageDisplayParseInputs({'Parent','Border','Reduce'},preparsed_varargin{:});

[jonas]

If you provide the image (not graph with axes) and the exact code then I can see what is wrong.

Or do you mean that you got the error with my cameraman.tif example?

[jonas]

Thought your original image was grayscale.. I've edited the answer for future reference, but honestly IA's answer is better. Cheers.

[Image Analyst]

I never use imagesc(). Use imshow() and pass in the axes handle and the colormap you want for that axes. For example:

GH = imshow(flipud(rgb2), 'ColorMap', hot(256), 'Parent', ax2);

By the way, RGB color images do not have colormaps. It simply does not apply. They already have an RGB color for every pixel, unlike a grayscale (indexed) image that has a single value that we can apply a color to.

[Turbulence Analysis]

Thanks for your response.

Actually, I need to display image with x,y axes values i.e. imagesc(x,y,flipud(rgb2), thats why I prefer imagesc. If I use x,y, in imshow command it throws error...

[Image Analyst]

What error? You can have the axes in pixels, going either up or down. Or you can have it in arbitrary units if you specify the 'XData' and 'YData' options. What tick labels do you want along the axes of your image, and how did you call imshow()? And what was the error message?

[Turbulence Analysis]

Hi,

I used as shown below,

I am sure that, I am not defining the x, y values inside imshow () in proper way..

Attached here is the rgb1, x, y, matrix

ax1 = gca;
GH = imshow(x,y,flipud(rgb1), 'ColorMap', hot(256), 'Parent', ax1);

%% error meassage

Error using images.internal.imageDisplayParsePVPairs (line 86)
Too many input arguments were passed to imageDisplayParsePVPairs.
Error in images.internal.imageDisplayParseInputs (line 69)
[common_args,specific_args] = images.internal.imageDisplayParsePVPairs(varargin{:});
Error in imshow (line 245)
    images.internal.imageDisplayParseInputs({'Parent','Border','Reduce'},preparsed_varargin{:}); 

[Image Analyst]

Like I said, you need to use the 'XData' and 'YData' options:

GH = imshow(flipud(rgb1), 'ColorMap', hot(256), 'Parent', ax1, 'XData', x, 'YData', y);

[Turbulence Analysis]

Thanks, I am getting axis values now. I just added 'axis on'' just below this..

However, Still I am not getting two different color bars..

Below is the script I used..Attached data matrix for your refrence...

% first image
ax1=gca;
G = imshow(flipud(rgb1), 'ColorMap', jet(256), 'Parent', ax1, 'XData', x, 'YData', y);
axis on
colormap(ax1,'jet');
%caxis ([0 3000]);
%c = colorbar('FontSize',20,'TickLabelInterpreter', 'latex');
c = colorbar(ax1,'Location', 'east','FontSize',20,'TickLabelInterpreter', 'latex');
G.AlphaData = 0.5;
set(gca, 'YDir','normal')
axis equal
hold on;
% linkaxes([ax1,ax2])
% Second image
ax2=gca;      
GH = imshow(flipud(rgb2), 'ColorMap', hot(256), 'Parent', ax2, 'XData', x, 'YData', y);
axis on
colormap(ax2,'hot')
% caxis ([0 3000]);
%c1 = colorbar('FontSize',20,'TickLabelInterpreter', 'latex');
c1 = colorbar(ax2,'Location', 'west','FontSize',20,'TickLabelInterpreter', 'latex');
GH.AlphaData = 0.5;
set(gca, 'YDir','normal')
axis equal

[Image Analyst]

Like I said, colormaps and colorbars don't apply to RGB images.

[Turbulence Analysis]

Okay, thanks. Atleast, with RGB I can able to display images with two different color patterns..

Furthermore, I intend to read images placed inside folder and subfolder within a for loop.. Is there a way to achieve this within a single for loop..

For e.g, the first part of code needs to read image from folder named 'A' and second part needs to read from subfolder 'A1'

%%%%%%% Needs to be read from folder A

for f = 1:1:1 if (f>=1) && (f<=9) fname_strt = 'B0000' ; elseif (f>=10) && (f<=99) fname_strt='B000'; elseif (f>=100) && (f<=999) fname_strt='B00'; else fname_strt='B0'; end fname_end = num2str(f); fname = strcat(fname_strt,fname_end,'.jpg'); I=imread(fname); G = imagesc(x,y,flipud(I)); G.AlphaData = 0.5; hold on;

%%%%%%% Needs to be read from subfolder A1 which is inside A

    if (f>=1) && (f<=9)
        fname_strt = 'B0000' ;
    elseif (f>=10) && (f<=99)
        fname_strt='B000';
    elseif (f>=100) && (f<=999)
        fname_strt='B00';
    else
        fname_strt='B0';
    end
    fname_end = num2str(f);
    fname = strcat(fname_strt,fname_end,'.jpg');
    h=readimx(fname);
H1 = imagesc(flipud(h));
H1.AlphaData = 0.5;
axis equal
hold on;
end

[Image Analyst]

Yes, see the FAQ:

https://matlab.fandom.com/wiki/FAQ#How_can_I_process_a_sequence_of_files.3F

[Turbulence Analysis]

Thanks for sharing the link..

I will work out this..