|
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
|