MATLAB: Can I more accurately measure this pixel streak

image processingImage Processing Toolbox

For reference, here is the image I am looking at: http://i630.photobucket.com/albums/uu30/Coulton/Image1.jpg

I am trying to find the orientation and center of a streak in an image. The image I am looking at can be found at the link about. For those inquiring, it is a aluminum strip falling in water.

I want to find the orientation and center as mentioned above. I have been able to measure theta fairly accurately. However, when trying to find the center of the streak I obtain a slight difference in the length calculation from the boundingbox cropping image and the overall binary image. Can anyone explain the miscalculation? I am sure it is user error, but I only have a year or so of Matlab experience so I am still learning. Also, is there a way I can measure the center without using for loops and if structures? Let me know if you need more info. The code is below:

clc;    % Clear the command window.
close all;  % Close all figures (except those of imtool.)
clear;  % Erase all existing variables.
workspace;  % Make sure the workspace panel is showing.
format longg;
format compact;
fontSize = 20;
Length = input('Enter length of strip (in inches): ');
% Read in a standard MATLAB gray scale demo image.
folder = 'K:\Matlab_Research';
baseFileName = 'image1.jpg';
% Get the full filename, with path prepended.
fullFileName = fullfile(folder, baseFileName);
% Check if file exists.
if ~exist(fullFileName, 'file')
  % File doesn't exist -- didn't find it there.  Check the search path for it.
  fullFileName = baseFileName; % No path this time.
  if ~exist(fullFileName, '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);
% Crop away tick marks, colorbar, etc..
grayImage = rgb2gray(imcrop(grayImage, [1 1 725 725]));
% Get the binary image
% Use an algorithm to find the threshold.
binaryImage = grayImage > 100;
% Smooth it out some and clean it up.
binaryImage1 = bwareaopen(binaryImage, 80); 
labeledImage = bwlabel(binaryImage1);
measurements = regionprops(labeledImage, 'all');
% Crop it to the bounding box.
bb = measurements.BoundingBox;
% Crop the image.
binaryImage = imcrop(binaryImage1, bb);
% Skeletonize the image.
skeleton = bwmorph(binaryImage, 'skel', 'inf');
% Find the endpoints of the skeleton.
skeleton = bwmorph(skeleton, 'endpoints');
% Get the coordinates of all endpoints
[rows cols] = find(skeleton);
% Find the two that are farthest apart
maxDistance = -1;
for k1 = 1 : length(rows)
  row1 = rows(k1);
  col1 = cols(k1);
  for k2 = 1 : length(rows)
    row2 = rows(k2);
    col2 = cols(k2);
    distance = sqrt((row1-row2)^2 + (col1-col2)^2);
    if distance > maxDistance
      bestIndex1 = k1;
      bestIndex2 = k2;
      maxDistance = distance;
              % Finding x and y component distances
              y = abs(row2-row1);
              x = abs(col2-col1);
      end
    end
  end
[rows cols] = find(binaryImage1);
% Find the two that are farthest apart from uncropped image
maxDistance1 = -1;
for k11 = 1 : length(rows)
  row11 = rows(k11);
  col11 = cols(k11);
  for k22 = 1 : length(rows)
    row22 = rows(k22);
    col22 = cols(k22);
    distance = sqrt((row11-row22)^2 + (col11-col22)^2);
    if distance > maxDistance1
      maxDistance1 = distance;
              % Store max point locations
              y1 = row11;
              y2 = row22;
              x1 = col11;
              x2 = col22;
              % Find length of strip from uncropped image
              maxDistanceloc = sqrt((y1-y2)^2 + (x1-x2)^2);
              % Determine horizontal location
              if row11 > row22
                  ybar = row11 - abs(row11-row22)/2;
              else
                  ybar = row22 - abs(row11-row22)/2;
              end
              % Deteremine vertical location
              if col11 > col22
                  xbar = col11 - abs(col11-col22)/2;
              else
                  xbar = row22 - abs(row11-row22)/2;
              end
      end
    end
  end
% Angle from horizontal
Theta = atan(y/x)*180/pi;
% Determine pixels per inch
PPI = maxDistance/Length;
% Determine center location from top left of image
xloc = xbar/PPI;
yloc = ybar/PPI;
% Difference in strip length from cropped to uncropped image
SLD = (maxDistance1 - maxDistance)/PPI;
fprintf('\nStrip center location from the top left of the image: ') 
fprintf('y = %2.4f inches and x = %2.4f inches\n\n',yloc,xloc')
fprintf('Theta = %2.2f degrees\n\n',Theta)
fprintf('Difference in strip lengths: %.4f inches\n',SLD)

Best Answer

You can get the centroid location of the binary image before you skeletonize it.

measurements = regionprops(binaryImage, 'Centroid');
centroid = measurements(1).Centriod;

Related Solutions

MATLAB: Finding tortuosity of parasite

Use bwskel() to get the longest spine. Then use bwmorph() to get the endpoints of the spine. Then use find() and sqrt() to get the straight line distance between the endpoints, then divide the two.

% This demo computes the tortuosity of a snake-like shape by dividing the actual longest spine distance by the straight line distance between the endpoints.
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;
%--------------------------------------------------------------------------------------------------------

%    READ IN IMAGE
folder = pwd;
baseFileName = 'worm.png';
% 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
% Read in the image from disk.  If storedColorMap is not empty, it's an indexed image with a stored colormap.
[grayImage, storedColorMap] = imread(fullFileName);
if ~isempty(storedColorMap)
	grayImage = ind2rgb(grayImage, storedColorMap);
end
% 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.
hFig = figure;
subplot(1, 2, 1);
imshow(grayImage, []);
title('Original Grayscale Image', 'FontSize', fontSize, 'Interpreter', 'None');
impixelinfo;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;
%--------------------------------------------------------------------------------------------------------
%    SEGMENTATION
% Binarize the image
binaryImage = imbinarize(grayImage);
% Take the largest blob only.
binaryImage = bwareafilt(binaryImage, 1);
% Fill holes, if any.
binaryImage = imfill(binaryImage, 'holes');
% Display the binary image.
subplot(1, 2, 1);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
% Now skeletonize
skelImage = bwskel(binaryImage); % First guess.  Will have short spurs.
% First get the area of the skeleton and take half of it to be the longest spine length.
% This will get rid of all spurs.
MinBranchLength = round(sum(skelImage(:))/2)
% Now take the skeleton again using that min branch length we just computed.
skelImage = bwskel(binaryImage, 'MinBranchLength', MinBranchLength);
% Display the skeleton image.
subplot(1, 2, 2);
imshow(skelImage, []);
title('Skeleton Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
% Find the endpoints
endpointImage = bwmorph(skelImage, 'endpoints');
[rows, columns] = find(endpointImage)
% Draw a red line between them.
hold on;
plot(columns, rows, 'r-', 'LineWidth', 2);
% Find the spine length.
% (Would be more accurate to find the spine length by traversing it.)
spineLength = sum(skelImage(:))
% Compute the straight line distance.
straightLineDistance = sqrt((columns(2) - columns(1))^2 + (rows(2) - rows(1))^2)
% Compute the tortuosity
tortuosity = spineLength / straightLineDistance
caption = sprintf('Skeleton Image.  Tortuosity = %f', tortuosity);
title(caption, 'FontSize', fontSize, 'Interpreter', 'None');

MATLAB: Error when reading data

If more than one region is found in the bwlabel, then the regionprops is going to return a structure array, and then measurements.BoundingBox would be attempting to expand a structure array, leading to the error.

Best Answer

Related Solutions

MATLAB: Finding tortuosity of parasite

MATLAB: Error when reading data

Related Question