How to draw circle in a 3D space?

Hi,

I would like to draw 2 circles. One in plane [P1 P6 P5], the other in plane [P1 P8 P5].

RGB values of points P1, P5, P6 and P8 are known. E.g., the RGB values of P1 is [R1 G1 B1]; the RGB values of P5 is [R5 G5 B5].

Requirements:

1. Using line P1P5 as the diameter
2. Using the mid point of line P1P5 as the center
3. P6 and P8 are points only used for deciding the location of the planes, P6 and P8 does not necessarily need to be contained in the circles. They could be outside of the circles.

Any thoughts please?

Answer 1

[Fabio Freschi]

The solution requires the definition of a rotation matrix. There are several ways to use it. Matlab has the rotx, roty and rotz functions, but they only work with one rotation at time. My implementation (see attachment) works by defining the new coordinate system identified by the position of the new center P0, any point along the new z axis Pz, and any point along the new x axis Px. So the call is R = rotationMatrix(P0,Pz,Px); I understand that these few words can be complex to understand but if you need details, I can explain further. Now the code

clear variables, close all

% points
P1 = [0 0 0]; P8 = [0 0 1]; P5 = [1 .5 2]; P6 = [2 2 3];

% plot
h = figure; hold on
plot3(P1(1),P1(2),P1(3),'o')
plot3(P8(1),P8(2),P8(3),'o')
plot3(P5(1),P5(2),P5(3),'o')
plot3(P6(1),P6(2),P6(3),'o')
axis equal, view([1 1 1])

% center
C0 = (P1+P5)/2;  % center
figure(h), plot3(C0(1),C0(2),C0(3),'*')

% radius
R = sqrt(sum((P5-P1).^2,2))/2;

% circumpherence in local coordinates
nPt = 100;
alpha = linspace(0,2*pi,nPt).';
P = [R.*cos(alpha) R.*sin(alpha) zeros(nPt,1)];

% direction orthogonal to circumpherence 1
v1 = cross(P8-P1,P5-P1);
figure(h), quiver3(C0(1),C0(2),C0(3),v1(1),v1(2),v1(3))
% rotation matrix
R1 = rotationMatrix(C0,C0+v1,P5);
% circumpherence 1 rotate points
Pc1 = P*R1'+C0;
plot3(Pc1(:,1),Pc1(:,2),Pc1(:,3),'-');


% direction orthogonal to circumpherence 2
v2 = cross(P6-P1,P5-P1);
figure(h), quiver3(C0(1),C0(2),C0(3),v2(1),v2(2),v2(3))
% rotation matrix
R2 = rotationMatrix(C0,C0+v2,P5);
% circumpherence 2 rotate points
Pc2 = P*R2'+C0;
plot3(Pc2(:,1),Pc2(:,2),Pc2(:,3),'-');

Please check if I have provided everything you need, because I have many personal utility functions in matalb path

[rough93]

rotationmatrix is undefined here, are you using a toolbox or predefined personal function?

[Fabio Freschi]

the file is attached

[Salad Box]

Thanks Fabio. I modified RGB values of P1,5,6 and 8 based on my own situation. Netherless it wouldn't change your method.

Also to make the case simpler, I used P1, P5, P6 as one example to draw a circle.

I also replaced plot3 with scatter3 due to my own preference.

clear variables
% points
P1 = [1,0,0];  P5 = [239,185,138]; P6 = [1,185,138];

% plot
figure;
 
hp1 = scatter3(P1(1),P1(2),P1(3),'o');
tp1 = text(P1(1), P1(2), P1(3), 'P1');
xlim([0 255]); ylim([0 255]); zlim([0 255]);

hold on
hp5 = scatter3(P5(1),P5(2),P5(3),'o')
tp5 = text(P5(1), P5(2), P5(3), 'P5');

hold on
hp6 = scatter3(P6(1),P6(2),P6(3),'o')
tp6 = text(P6(1), P6(2), P6(3), 'P6');

% connect p1, p5 and p6 with lines
hd15 = line([P1(1),P5(1)],[P1(2),P5(2)],[P1(3),P5(3)],'linestyle','--');
hd16 = line([P1(1),P6(1)],[P1(2),P6(2)],[P1(3),P6(3)],'linestyle','--');
hd56 = line([P6(1),P5(1)],[P6(2),P5(2)],[P6(3),P5(3)],'linestyle','--');

% center
C0 = (P1+P5)/2;

hold on
scatter3(C0(1), C0(2), C0(3), '*');

% radius
R = sqrt(sum((P5-P1).^2))/2;

% circumpherence in local coordinates
nPt = 100;
alpha = linspace(0,2*pi,nPt)';
P = [R.*cos(alpha) R.*sin(alpha) zeros(nPt,1)];
hold on
hc_local = scatter3(P(:,1), P(:,2), P(:,3),'o');

xlim([-255 255]); ylim([-255 255]); zlim([-255 255]);

% direction orthogonal to circumpherence 1 (P1 P6 P5)
    v1 = cross(P6 - P1, P5 - P1);
    hold on
    quiver3(C0(1),C0(2), C0(3), v1(1), v1(2), v1(3))

    % rotation matrix
    R1 = rotationMatrix(C0, C0 + v1, P5);
    % circumpherence 1 rotate points
    Pc1 = P*R1' + C0;
    hold on
    scatter3(Pc1(:,1), Pc1(:,2), Pc1(:,3), 'ro')

    xlim([-300 300]); ylim([-300 300]); zlim([-300 300]);

My questions are:

1. I guess the sizes of red circle and green circle won't be the same (different radius values)?
2. Why do we need the green circle in order to draw the red circle?
3. Why do we need v1 = cross(P6 - P1, P5 - P1); and why do we need to move the center using quiver3(C0(1),C0(2), C0(3), v1(1), v1(2), v1(3))?
4. About the rotationMatrix function you wrote, R = rotationMatrix(P0,Pz,Px) INPUT: 'P0': system origin; 'Pz': point on z-axis; 'Px': point on x-axis; R1 = rotationMatrix(C0, C0 + v1, P5); why Pz = C0 + v1 and why Px = P5?
5. What is the equation of the red circle? Pc1 = P*R1' + C0? What is the equation of the red circle in terms of variables R G B? Because the purpose of this whole thing is trying to select points that are evenly distributed on the circumpherence of this red circle.

[Fabio Freschi]

My answers

1) the radius of the green circle equals the radius of the red one: in the original post I understood that the diameter is P1-P5, so I calculated the radius as

R = sqrt(sum((P5-P1).^2))/2;

2) The key of the method is to create the green circle in a local reference frame (simple in local coordinates, e.g. x-y plane, centered) then rotate according to R1 or R2 and finally translate it to the midpoint of P1-P5

3) v1 is the vector orthogonal to the plane wher you want to lay your circle. P1-P5-P6 define a plane and the cross product of two vectors lying on a plane gives a vector orthogonal to that plane. Quiver does not move the the center. I wanted to visualize the plane with its orthogonal vector (note here that any direction is good, so both v1 and -v1 are fine).

4) We have to define a coordinate system x'y'z' lying on the plane of the circle, identified by P0, Pz and Px. P0 is the origin of x'y'z', so in our case it is C0. Pz is any point along z'. Because the z' axis is orthogonal to the plane, we can optain Pz as C0+v1 (actually any C0+k*v1, with real k is fine). Px is any point on the x' axis. We have a circle, so the x' axis can be any point in the plane of the circle. I have chosen P5 because it is ready (but also P1 is good). See the picture below

5) I think your equation is correct. In any case, the point on the red circumherence are evenly distributed along the azimuthal coordinate

alpha = linspace(0,2*pi,nPt).';

Hope it helps. If the answer fits your request, please accept the answer

Fabio

[Salad Box]

Thanks for your detailed response Fabio. I have acceped your answer.

I followed all your answers apart from 4):

"4) We have to define a coordinate system x'y'z' lying on the plane of the circle"...

To me x'y'z' doesn't seem to be lying on the same plane as the circle. z' seems more orthogonal to the circle. So I wonder how x'y'z' is accurately defined. How to decide axis x', and axis y', if axis z' is the axis where line p1 to p5 sits on?

[Fabio Freschi]

if you look at the green circle, it lays in the xy plane and the orthogonal vector is along the z axis. We want a similar situation for the red circle in the x'y'z' coordinate system: it must lay in the x'y' plane and the z' axis must be orthogonal to that plane. so, for x'y'z':

• the origin O' is in C0 (P0 = C0)
• the z' axis is along the v1 vector. Accordingly, a point along z' is Pz = C0+v1
• the x' and y' axes are on the plane orthogonal to v1. In this case the position of the x' axis is not important; for simplicity I decided to align x' with the line passing from P5-C0, so a point along the x' axis is Px = P5;

[Salad Box]

Thanks Fabio for your precise response.

I'm hoping that I understand it correctly.

• x' and y' are in the same plane as the red circle while axis x' is perpendicular to axis y'.
• z' is orthogonal to plane x'y'.

Thanks a million for your help!

[Fabio Freschi]

exact!

[Salad Box]

One more question:

Why Pc1 = P*R1' + C0?

[Fabio Freschi]

P are the points along the green circle R1 is the rotation matrix, C0 the center of the roto-translated circle. So: the green circumpherence is rotated (P*R1') and then translated to the new center (P*R1'+C0)

Answer 2

[darova]

Another version

Rotation matrix from HERE