MATLAB: Curve fitting and scatter plots

scatter plots polynomials splines

Hi all,

This may not be possible but is there a function or a way you can fit a curve to a select number of points within a scatterplot? Looking specifically to fit a curve around the 'lower boundary' of a scatterplot. For example say I had a scatterplot of hundreds of points that made up a circle. Is there a way to fit a curve around the bottom perimeter?

If you need any more info please ask.

Best

Best Answer

Craig, try this to see if it's what you want. I create a bunch of noisy points around a circle. Then I select the bottom half of those. Then I fit a circle through the points, using only those points in the bottom half of the circle to determine the fitted equation. I think that's what you've been saying you want. Save the following in test.m and run it.

function test
clc;    % Clear the command window.
workspace;  % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
xCenter = 12;
yCenter = 10;
theta = 0 : 0.01 : 2*pi;
radius = 5 + 2*rand(1, length(theta));
x = radius .* cos(theta) + xCenter;
y = radius .* sin(theta) + yCenter;
plot(x, y, 'b*');
axis square;
xlim([0 20]);
ylim([0 20]);
grid on;
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Get the bottom half of them
selector = y < yCenter;
bottomHalfY = y(selector);
bottomHalfX = x(selector);
hold on;
plot(bottomHalfX, bottomHalfY, 'ro');
% Fit a circle to the bottom half noisy data:
[xc,yc,R,a] = circfit(x,y)
xFit = R .* cos(theta) + xc;
yFit = R .* sin(theta) + yc;
plot(xFit, yFit, 'r-', 'LineWidth', 3);
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
% http://matlab.wikia.com/wiki/FAQ#How_can_I_fit_a_circle_to_a_set_of_XY_data.3F
% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
% by Bucher izhak 25/oct/1991
n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));

Related Solutions

MATLAB: How to create a circles (smallest and biggest circle) based on points in an image given, then find the center and radius

To get the minimal bounding outer circle, try this:

John D'Errico's File Exchange

I asked John for the largest interior circle, and he said that's a much more difficult problem and he doesn't have code for that.

MATLAB: Drawing a circle around optic disk

I would not do an edge detection. That's a typical "new bee" thing to do and rarely works. Since your optic disk is a bright blob on the left part of the image and known/assumed to be a perfect circle, you can just use that knowledge to your advantage and threshold the image to find the disc's region (or at least most of it), then fit the blob's "good" right side to a circle using code from the FAQ. Here is the code:

% Routine to find the optic disc in one particular image.
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 = 20;
%===============================================================================
% Read in gray scale demo image.
folder = pwd; % Determine where demo folder is (works with all versions).
baseFileName = 'i.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
rgbImage = 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(rgbImage)
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(rgbImage);
  % 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 = rgbImage(:, :, 2); % Take green channel.
else
  grayImage = rgbImage; % It's already gray scale.
end
% Now it's gray scale with range of 0 to 255.
% Display the image.


subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();
%------------------------------------------------------------------------------
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% 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;
% Optic disk is known to be on the left, so erase the image on the right.
clippedImage = grayImage; % Initialize.
clippedImage(:, 1200:end) = grayImage(1,1); % Set equal to surrounding template.
% Show the histogram
subplot(2, 2, 2);
[counts, grayLevels] = histcounts(clippedImage(:));
counts(counts == max(counts)) = 0; % Suppress spike
bar(grayLevels(1:end-1), counts);
grid on;
title('Histogram of Original Image', 'FontSize', fontSize, 'Interpreter', 'None');
% Binarize the image
binaryImage = clippedImage >= 245;
% Take largest blob
binaryImage = bwareafilt(binaryImage, 1);
% Take convex hull
binaryImage = bwconvhull(binaryImage);
% Display the image.
subplot(2, 2, 3);
imshow(binaryImage, []);
title('Binary Image', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();
% It is known/assumed that the optic disk's "good part" is on the right,
% so find pixels on the right side.
% First find boundaries
boundaries = bwboundaries(binaryImage);
xBoundary = boundaries{1}(:, 2);
yBoundary = boundaries{1}(:, 1);
% Scan down finding right edge
xRight = zeros(1, columns);
yRight = zeros(1, columns);
for row = 1 : rows
  x = find(binaryImage(row, :), 1, 'last');
  if ~isempty(x)
    xRight(row) = x;
    yRight(row) = row;
  end
end
% Get rid of rows with no object in them.
noObjectRows = xRight == 0;
xRight(noObjectRows) = [];
yRight(noObjectRows) = [];
% Fit xRight and yRight to a circle
[xCenter, yCenter, R, a] = circfit(xRight,yRight);
% Plot over image.
% Display the image.
subplot(2, 2, 4);
imshow(grayImage, []);
title('Original Image with Disc Circled', 'FontSize', fontSize, 'Interpreter', 'None');
axis('on', 'image');
hp = impixelinfo();
viscircles([xCenter, yCenter], R);
% Circle fitting routine from the FAQ.
function [xc,yc,R,a] = circfit(x,y)
%CIRCFIT  Fits a circle in x,y plane
%

% [XC, YC, R, A] = CIRCFIT(X,Y)
% Result is center point (yc,xc) and radius R.  A is an optional
% output describing the circle's equation:
%
%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
% by Bucher izhak 25/oct/1991
n=length(x);  xx=x.*x; yy=y.*y; xy=x.*y;
A=[sum(x) sum(y) n;sum(xy) sum(yy) sum(y);sum(xx) sum(xy) sum(x)];
B=[-sum(xx+yy) ; -sum(xx.*y+yy.*y) ; -sum(xx.*x+xy.*y)];
a=A\B;
xc = -.5*a(1);
yc = -.5*a(2);
R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));
end

Best Answer

Related Solutions

MATLAB: How to create a circles (smallest and biggest circle) based on points in an image given, then find the center and radius

MATLAB: Drawing a circle around optic disk

Related Question