Path: news.mathworks.com!not-for-mail
From: "Sadik " <sadik.hava@gmail.com>
Newsgroups: comp.soft-sys.matlab
Subject: Re: How can i plot filled round
Date: Sat, 23 May 2009 20:03:02 +0000 (UTC)
Organization: The MathWorks, Inc.
Lines: 76
Message-ID: <gv9kpm$t51$1@fred.mathworks.com>
References: <gv9ajd$q6$1@fred.mathworks.com> <15576475.11032.1243101515834.JavaMail.jakarta@nitrogen.mathforum.org>
Reply-To: "Sadik " <sadik.hava@gmail.com>
NNTP-Posting-Host: webapp-03-blr.mathworks.com
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1243108982 29857 172.30.248.38 (23 May 2009 20:03:02 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sat, 23 May 2009 20:03:02 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 1666517
Xref: news.mathworks.com comp.soft-sys.matlab:542066


I had done something similar in the past and fortunately I found the code. At a first glance, it may look complex but in deed it is not. What it does is, it first assumes the major axis is along the x-axis and the minor is along the y-axis [although you can still input values such that major < minor. This is just the naming convention of the code]. It generates points inside this horizontal ellipse [satisfying a smaller-than-or-equal-to inequality] and the number of points is determined by the user. The other inputs are the center (x0,y0) of the ellipse, and the angle it makes with the x-axis [angleOfRotation], as well as the major axis and minor axis lengths. The final input is just to control drawing. If it is 1, the points are plotted with the actual elliptic border surrounding and the two axes. If it is 0, no drawing is performed. In both cases, the output of the function is the set 
of generated points. 

An example usage:

pointsGenerated = generatePointsInsideEllipse(3,2,5,3,pi/4,1000,1);

A set of points in an ellipse at an angle of pi/4 with the x-axis, centered at (3,2) with axis lengths major = 5, minor = 3. The total number of points to be generated is 1000, but not all of them may fall into the ellipse. Finally, we request a draw.

The code is as follows. Please copy and paste the entire code to a single blank m-file. [If you have any questions, please feel free to ask me at sadik.hava@gmail.com]


% CODE STARTS HERE
function pointsGenerated = generatePointsInsideEllipse(centerX0,centerY0,majorAxis,minorAxis,angleOfRotation,totalNumberOfPoints,drawAlso)

giveEllipse(centerX0,centerY0,majorAxis,minorAxis,angleOfRotation,drawAlso);

rotationMatrix = [cos(angleOfRotation) -sin(angleOfRotation);sin(angleOfRotation) cos(angleOfRotation)];

generatedX = unifrnd(-majorAxis,majorAxis,1,totalNumberOfPoints);
generatedY = unifrnd(-minorAxis,minorAxis,1,totalNumberOfPoints);

pointsGenerated = zeros(2,totalNumberOfPoints);

numberInside = 0;

for n = 1:totalNumberOfPoints
    
    if zeroMeanEllipseEquation(generatedX(n),generatedY(n),majorAxis,minorAxis) <= 0
        
        numberInside = numberInside + 1;
        pointsGenerated(1,numberInside) = generatedX(n);
        pointsGenerated(2,numberInside) = generatedY(n);
        
    end
    
end

pointsGenerated(:,numberInside+1:end) = [];

pointsGenerated = rotationMatrix*pointsGenerated+[centerX0*ones(1,length(pointsGenerated));centerY0*ones(1,length(pointsGenerated))];


if drawAlso
    
    scatter(pointsGenerated(1,:),pointsGenerated(2,:),'.')
    [u,s,v] = svd(cov(pointsGenerated'));
    line([centerX0,centerX0+2*sqrt(s(2,2))*u(1,2)],[centerY0,centerY0+2*sqrt(s(2,2))*u(2,2)],'Color','Red','Linewidth',3)
    line([centerX0,centerX0+2*sqrt(s(1,1))*u(1,1)],[centerY0,centerY0+2*sqrt(s(1,1))*u(2,1)],'Color','Red','Linewidth',3)
    
end

function [xPrimeUp,xPrimeDown,yPrimeUp,yPrimeDown] = giveEllipse(centerX0,centerY0,majorAxis,minorAxis,angleOfRotation,drawAlso)

x = (-majorAxis:0.01:majorAxis)+centerX0;
yUp = minorAxis*sqrt(1 - ((x-centerX0)/majorAxis).^2)+centerY0;
yDown = -minorAxis*sqrt(1 - ((x-centerX0)/majorAxis).^2)+centerY0;
yPrimeUp = (x-centerX0)*sin(angleOfRotation)+(yUp-centerY0)*cos(angleOfRotation)+centerY0;
yPrimeDown = (x-centerX0)*sin(angleOfRotation)+(yDown-centerY0)*cos(angleOfRotation)+centerY0;
xPrimeUp = (x-centerX0)*cos(angleOfRotation)-(yUp-centerY0)*sin(angleOfRotation)+centerX0;
xPrimeDown = (x-centerX0)*cos(angleOfRotation)-(yDown-centerY0)*sin(angleOfRotation)+centerX0;

if drawAlso
        
    figure;plot(xPrimeUp,yPrimeUp)
    hold
    plot(xPrimeDown,yPrimeDown)
    
end


function result = zeroMeanEllipseEquation(givenX,givenY,majorAxis,minorAxis)

result = (givenX/majorAxis)^2 + (givenY/minorAxis)^2 - 1;

% CODE ENDS HERE