Plotting data from a loop

 for k = 1:length(s)
        xbar = s(k).Centroid(1);
        ybar = s(k).Centroid(2);

        a = s(k).MajorAxisLength/2.6;
        b = s(k).MinorAxisLength/2.7;

        theta = pi*s(k).Orientation/180;
        R = [ cos(theta)   sin(theta)
            -sin(theta)   cos(theta)];

        xy = [a*cosphi; b*sinphi];
        xy = R*xy;

        x = xy(1,:) + xbar;
        y = xy(2,:) + ybar;

        plot(x,y,'r','LineWidth',2);

    end

My code shown above is part of a greater program which plots an ellipse onto a set of 4 images. The ellipse axes (stored in a and b) are growing throughout my loop, and I want to store them in an array to plot them when the loop is finished. Whenever I've been trying to plot, I've just been able to plot one data point.

How do I plot a and b after the loop has finished?

Answer 1

[Azzi Abdelmalek]

Use a(k) and b(k)

Answer 2

[Image Analyst]

You just need to use hold on, and to keep track of the rotation angle and the angle of points around the ellipse (you need to keep them separate). Try it like this:

imshow(originalImage); % Display image.
hold on;
axis on;
theta = 0 : 0.01 : 2*pi;
% Plot ellipses over image.
for k = 1:length(s)
  xCenter = s(k).Centroid(1);
  yCenter = s(k).Centroid(2);
  xRadius = s(k).MajorAxisLength/2;
  yRadius = s(k).MinorAxisLength/2;
  x = xRadius * cos(theta);
  y = yRadius * sin(theta);
  rotationAngle = pi*s(k).Orientation/180;
  R = [cos(rotationAngle)   sin(rotationAngle)
    -sin(rotationAngle)   cos(rotationAngle)];
  
  xy = [x;y];
  xy = R*xy;
  
  x = xy(1,:) + xCenter;
  y = xy(2,:) + yCenter;
  
  plot(x,y,'r','LineWidth',2);
end

If you want the full-blown demo, where I used the FAQ to create an ellipse and my BlobsDemo, here it is:

%------------------------------------------------------------------------------------------------
% Demo to illustrate simple blob detection, measurement, and filtering.
% Requires the Image Processing Toolbox (IPT) because it demonstates some functions
% supplied by that toolbox, plus it uses the "coins" demo image supplied with that toolbox.
% If you have the IPT (you can check by typing ver on the command line), you should be able to
% run this demo code simply by copying and pasting this code into a new editor window,
% and then clicking the green "run" triangle on the toolbar.
% Running time = 7.5 seconds the first run and 2.5 seconds on subsequent runs.
% A similar Mathworks demo:
% http://www.mathworks.com/products/image/demos.html?file=/products/demos/shipping/images/ipexprops.html
% Code written and posted by ImageAnalyst, July 2009.
%------------------------------------------------------------------------------------------------
% function BlobsDemo()
% echo on;
% Startup code.
tic; % Start timer.
clc; % Clear command window.
clearvars; % Get rid of variables from prior run of this m-file.
disp('Running BlobsDemo.m...'); % Message sent to command window.
workspace; % Make sure the workspace panel with all the variables is showing.
imtool close all;  % Close all imtool figures.
format longg;
format compact;
fontSize = 20;

% Check that user has the Image Processing Toolbox installed.
hasIPT = license('test', 'image_toolbox');
if ~hasIPT
  % User does not have the toolbox installed.
  message = sprintf('Sorry, but you do not seem to have the Image Processing Toolbox.\nDo you want to try to continue anyway?');
  reply = questdlg(message, 'Toolbox missing', 'Yes', 'No', 'Yes');
  if strcmpi(reply, 'No')
    % User said No, so exit.
    return;
  end
end

% Read in standard MATLAB demo image
baseFileName = 'coins.png';
folder = pwd;  % Look in the current folder.
fullFileName = fullfile(folder, baseFileName);
if ~exist(fullFileName, 'file')
  % It doesn't exist in the current folder.
  % Look on the search path.
  if ~exist(baseFileName, 'file')
    % It doesn't exist on the search path either.
    % Alert user that we can't find the image.
    warningMessage = sprintf('Error: the input image file\n%s\nwas not found.\nClick OK to exit the demo.', fullFileName);
    uiwait(warndlg(warningMessage));
    fprintf(1, 'Finished running BlobsDemo.m.\n');
    return;
  end
  % Found it on the search path.  Construct the file name.
  fullFileName = baseFileName; % Note: don't prepend the folder.
end
% If we get here, we should have found the image file.
originalImage = imread(fullFileName);
% Check to make sure that it is grayscale.
[rows, columns, numberOfColorBands] = size(originalImage);
if numberOfColorBands > 1
  promptMessage = sprintf('Your image file has %d color bands.\nThis demo was designed for grayscale images.\nDo you want me to convert it to grayscale for you so you can continue?', numberOfColorBands);
  button = questdlg(promptMessage, 'Continue', 'Convert and Continue', 'Cancel', 'Convert and Continue');
  if strcmp(button, 'Cancel')
    fprintf(1, 'Finished running BlobsDemo.m.\n');
    return;
  end
  % Do the conversion using standard book formula
  originalImage = rgb2gray(originalImage);
end

% Display the grayscale image.
subplot(3, 3, 1);
imshow(originalImage);
% Maximize the figure window.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Force it to display RIGHT NOW (otherwise it might not display until it's all done, unless you've stopped at a breakpoint.)
drawnow;
caption = sprintf('Original "coins" image showing\n6 nickels (the larger coins) and 4 dimes (the smaller coins).');
title(caption);
axis square; % Make sure image is not artificially stretched because of screen's aspect ratio.

% Just for fun, let's get its histogram.
[pixelCount, grayLevels] = imhist(originalImage);
subplot(3, 3, 2);
bar(pixelCount); title('Histogram of original image');
xlim([0 grayLevels(end)]); % Scale x axis manually.
grid on;

% Threshold the image to get a binary image (only 0's and 1's) of class "logical."
% Method #1: using im2bw()
%   normalizedThresholdValue = 0.4; % In range 0 to 1.
%   thresholdValue = normalizedThresholdValue * max(max(originalImage)); % Gray Levels.
%   binaryImage = im2bw(originalImage, normalizedThresholdValue);       % One way to threshold to binary
% Method #2: using a logical operation.
thresholdValue = 100;
binaryImage = originalImage > thresholdValue; % Bright objects will be the chosen if you use >.
%   binaryImage = originalImage < thresholdValue; % Dark objects will be the chosen if you use <.

% Do a "hole fill" to get rid of any background pixels inside the blobs.
binaryImage = imfill(binaryImage, 'holes');

% Show the threshold as a vertical red bar on the histogram.
hold on;
maxYValue = ylim;
hStemLines = stem(thresholdValue, maxYValue(2), 'r');
children = get(hStemLines, 'children');
set(children(2),'visible', 'off');
% Place a text label on the bar chart showing the threshold.
annotationText = sprintf('Thresholded at %d gray levels', thresholdValue);
% For text(), the x and y need to be of the data class "double" so let's cast both to double.
text(double(thresholdValue + 5), double(0.5 * maxYValue(2)), annotationText, 'FontSize', 10, 'Color', [0 .5 0]);
text(double(thresholdValue - 70), double(0.94 * maxYValue(2)), 'Background', 'FontSize', 10, 'Color', [0 0 .5]);
text(double(thresholdValue + 50), double(0.94 * maxYValue(2)), 'Foreground', 'FontSize', 10, 'Color', [0 0 .5]);

% Display the binary image.
subplot(3, 3, 3); imagesc(binaryImage); colormap(gray(256)); title('Binary Image, obtained by thresholding'); axis square;

labeledImage = bwlabel(binaryImage, 8);     % Label each blob so we can make measurements of it
coloredLabels = label2rgb (labeledImage, 'hsv', 'k', 'shuffle'); % pseudo random color labels

subplot(3, 3, 4);
imshow(labeledImage, []);
title('Labeled Image, from bwlabel()');
axis square;
subplot(3, 3, 5);
imagesc(coloredLabels);
axis square;
caption = sprintf('Pseudo colored labels, from label2rgb().\nBlobs are numbered from top to bottom, then from left to right.');
title(caption);

% Get all the blob properties.  Can only pass in originalImage in version R2008a and later.
blobMeasurements = regionprops(labeledImage, originalImage, 'all');
numberOfBlobs = size(blobMeasurements, 1);

% bwboundaries() returns a cell array, where each cell contains the row/column coordinates for an object in the image.
% Plot the borders of all the coins on the original grayscale image using the coordinates returned by bwboundaries.
subplot(3, 3, 6);
imagesc(originalImage);
title('Outlines, from bwboundaries()'); axis square;
hold on;
boundaries = bwboundaries(binaryImage);
numberOfBoundaries = size(boundaries, 1);
for k = 1 : numberOfBoundaries
  thisBoundary = boundaries{k};
  plot(thisBoundary(:,2), thisBoundary(:,1), 'g', 'LineWidth', 2);
end
hold off;

fontSize = 14;  % Used to control size of "blob number" labels put atop the image.
labelShiftX = -7;  % Used to align the labels in the centers of the coins.
blobECD = zeros(1, numberOfBlobs);
% Print header line in the command window.
fprintf(1,'Blob #      Mean Intensity  Area   Perimeter    Centroid       Diameter\n');
% Loop over all blobs printing their measurements to the command window.
for k = 1 : numberOfBlobs           % Loop through all blobs.
  % Find the mean of each blob.  (R2008a has a better way where you can pass the original image
  % directly into regionprops.  The way below works for all versions including earlier versions.)
  thisBlobsPixels = blobMeasurements(k).PixelIdxList;  % Get list of pixels in current blob.
  meanGL = mean(originalImage(thisBlobsPixels)); % Find mean intensity (in original image!)
  meanGL2008a = blobMeasurements(k).MeanIntensity; % Mean again, but only for version >= R2008a
  
  blobArea = blobMeasurements(k).Area;    % Get area.
  blobPerimeter = blobMeasurements(k).Perimeter;    % Get perimeter.
  blobCentroid = blobMeasurements(k).Centroid;    % Get centroid one at a time
  blobECD(k) = sqrt(4 * blobArea / pi);          % Compute ECD - Equivalent Circular Diameter.
  fprintf(1,'#%2d %17.1f %11.1f %8.1f %8.1f %8.1f % 8.1f\n', k, meanGL, blobArea, blobPerimeter, blobCentroid, blobECD(k));
  % Put the "blob number" labels on the "boundaries" grayscale image.
  text(blobCentroid(1) + labelShiftX, blobCentroid(2), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');
end

% Now, I'll show you another way to get centroids.
% We can get the centroids of ALL the blobs into 2 arrays,
% one for the centroid x values and one for the centroid y values.
allBlobCentroids = [blobMeasurements.Centroid];
centroidsX = allBlobCentroids(1:2:end-1);
centroidsY = allBlobCentroids(2:2:end);
% Put the labels on the rgb labeled image also.
subplot(3, 3, 5);
for k = 1 : numberOfBlobs           % Loop through all blobs.
  text(centroidsX(k) + labelShiftX, centroidsY(k), num2str(k), 'FontSize', fontSize, 'FontWeight', 'Bold');
end

elapsedTime = toc;
subplot(3, 3, 8);
s = blobMeasurements;
imshow(originalImage); % Display image.
hold on;
axis on;
theta = 0 : 0.01 : 2*pi;
% Plot ellipses over image.
for k = 1:length(s)
  xCenter = s(k).Centroid(1);
  yCenter = s(k).Centroid(2);
  xRadius = s(k).MajorAxisLength/2;
  yRadius = s(k).MinorAxisLength/2;
  x = xRadius * cos(theta);
  y = yRadius * sin(theta);
  rotationAngle = pi*s(k).Orientation/180;
  R = [cos(rotationAngle)   sin(rotationAngle)
    -sin(rotationAngle)   cos(rotationAngle)];
  
  xy = [x;y];
  xy = R*xy;
  
  x = xy(1,:) + xCenter;
  y = xy(2,:) + yCenter;
  
  plot(x,y,'r','LineWidth',2);
end

[N.B.]

Thanks, Image Analyst.

Now, in order to plot these Major and Minor axes with respect to another variable, should I write a plot(xRadius, k) after the final end?

[Image Analyst]

I don't know what "plot these Major and Minor axes with respect to another variable" means. What is the other variable, and what does "with respect to" mean in that case? I plotted the ellipses that fit each blob, so what else do you need?

[N.B.]

Sorry for the confusion. I am looking to plot the length of the Major and Minor axes of each ellipse on an xy-scatter graph in a new figure.

[Image Analyst]

You can use the line() function. You know the centroid, and the lengths, and the angle, so it's simple algebra/trigonometry to determine the line endpoints.

[N.B.]

Could I do something like:

plot(xRadius, k)

after the loop has ended?

[Image Analyst]

I doubt it - what would that do?

[N.B.]

I'm just throwing it out here. I'm not trying to create a line object. All I want to do is plot the Major or Minor axes on a graph and see the pattern over the course of the 4 images that I am applying the ellipse to.

[Image Analyst]

You could certain plot any variable, such as xRadius, as a function of iteration if that's what you want.

[N.B.]

Yep, that's what I'm looking to do. I'm just not sure the best way to go about it.