Drawing circles on a torus in matlab

Question

Using the answer found here, I was able to draw a nice torus. I would like to draw two "circles" on this torus, the first of which is standard. The second should be somewhat random, perhaps like a ladder, ideally on the opposite side from the first. It might be kind of jagged, but should remain on the edges of the lattice and should connect back to itself. Here is what I have so far

am = 1.;
rm = 3.;
t = linspace(-pi,pi,16);
p = linspace(0.,2.*pi,16);
[t,p] = meshgrid(p,t);
x = (rm + am.*cos(p)).*cos(t);
y = (rm + am.*cos(p)).*sin(t);
z = am.*sin(p);
hsurf = surf(x,y,z);
axis equal;
set(hsurf, 'FaceColor','interp','FaceAlpha',0.5,'EdgeAlpha',0.25);
hold on
plot3((rm+am.*cos(p)).*cos(t(8)),(rm+am.*cos(p)).*sin(t(8)),am.*sin(p), 'b', 'LineWidth',2);

I have the standard loop, but can't get the jagged one. I have been trying to specify the other loop by specifying the vertices, but it doesn't seem to be working. Can I specify a curve on this mesh by specifying its vertices? If not, are there any suggestions? Thanks for any help you can provide.

EDIT: Sorry for being vague earlier. The picture below shows what I would like. Note that the right circle is given by the above code, whereas the left hand 'circle' is drawn by hand, and is what I would like to automate.

Answer 1

If you want the vertices of the (jagged or not) circle to match the points of the surface, you do not need to recalculate any points. You can simply re-use some of the points of your surface, choosing the path which define the shape you want.

I'll use your code to define the initial tore, with just a slight modification: I'll use the upper-case letter for variable which define a matrix (and lower case for variables which are simple 1D vector). So your initial tore becomes:

%% // define initial tore
am = 1.; rm = 3.;
t = linspace(-pi,pi,16);
p = linspace(0,2.*pi,16);
[T,P] = meshgrid(p,t);
X = (rm + am.*cos(P)).*cos(T);
Y = (rm + am.*cos(P)).*sin(T);
Z = am.*sin(P);
hsurf = surf(X,Y,Z,'FaceColor','interp','FaceAlpha',0.5,'EdgeAlpha',0.25);
axis equal ; hold on

Not a big deal but this allow to have only 1D vector coordinates in the case where you define your circle by equation as you were doing. It keeps your variables and workspace cleaner and reduce the risk of mistakes and dimension errors.

For the simple circle, you could redefine the coordinate by equation as you did (developed here for clarity)

%% // your simple circle - calculated
colID = 8 ;
xs = (rm+am.*cos(p)).*cos(t(colID)) ;
ys = (rm+am.*cos(p)).*sin(t(colID)) ;
zs = am.*sin(p) ;
hp = plot3(xs,ys,zs, 'b', 'LineWidth',2);
%// This was generating 2D arrays (matrices) for xs, ys and zs when you
%// were using the meshed version of T and P. Using only the vector version of
%// 't' and 'p' generate only vector coordinates. => cleaner and safer.

But I'd like to bring your attention to another way, just use the points of the surface already defined. This (i) eliminates the need for new calculations, and (ii) insure that the vertices of the circle will be matching exactly some vertices of the tore.

%% // another simple circle - extracted from the tore vertices
colID = 8 ;
xs = X(:,colID) ;
ys = Y(:,colID) ;
zs = Z(:,colID) ;
hpc = plot3(xs,ys,zs, 'b', 'LineWidth',2);

Again, there are extra lines for clarity but you can compact the code once you understand what is going on. These 2 different ways of producing a circle result in the following figures:

For illustration purpose, I highlighted the vertices of the tore (small black dots), so you can see if the circles vertices match or not.

---

Now for the jagged circle, the last method presented will have even more interest. As you observed, to get a circle coordinate we simply extracted one column of each of the tore matrix coordinates. We used a fixed value for the colID to retrieve. To "jagg" your circle, we just need to introduce a small perturbation in the column number, so we will occasionally retrieve the adjacent vertex coordinates.

The following code generate a column ID varying around the staring column ID. You can define the maximum total spread as well as the maximum single increment:

%% // generate index of columns to be retrieved
idxcol = zeros(size(X,1),1) ;
maxDeviation = 5 ;      %// total maximum deviation to the initial center line
maxPerturbation = 2 ;   %// max deviation for 1 increment

deviation = cumsum(randi([-maxPerturbation maxPerturbation],size(X,1),1)) ;
%// bound deviations to maximum values
deviation(deviation>maxDeviation)  = maxDeviation ;
deviation(deviation<-maxDeviation) = -maxDeviation ;

startColID = 8 ;
colID = startColID + deviation ;
%// close the profile by repeating first point in the end if it's not closed already
if colID(end) ~= colID(1) ; colID(end) = colID(1) ; end

If we now use coordinate from theses columns, we get:

npt = size(colID,1) ;
xcj = zeros(npt) ; ycj = xcj ; zcj = xcj ;
for k=1:npt
    xcj(k) = X( k , colID(k) ) ;
    ycj(k) = Y( k , colID(k) ) ;
    zcj(k) = Z( k , colID(k) ) ;
end
hpc = plot3(xcj,ycj,zcj, 'b', 'LineWidth',2);

It is a closed profile turning around the tore but the lines do not match the lines of the surface as in your drawn example. This is because the transition from one column ID to another happen at each change of x index.

---

To rectify that, we can use the stairs function:

%% // create a stepped profile
[xd,yd] = stairs(colID) ;

%% // display rectified jagged circle
npt = size(yd,1) ;
xcj = zeros(npt) ; ycj = xcj ; zcj = xcj ;
for k=1:npt
    xcj(k) = X( xd(k) , yd(k) ) ;
    ycj(k) = Y( xd(k) , yd(k) ) ;
    zcj(k) = Z( xd(k) , yd(k) ) ;
end
hpc = plot3(xcj,ycj,zcj, 'b', 'LineWidth',2);

This will now generate the following profile:

Much better ... However, when the column ID moves by more than one index for the same x, we still miss anchor points to keep the jagged circle profile matching the tore profile. If your basic column deviation increment is more than one, you will even have many more of these artefact.

---

To completely rectify that, we need to add the missing anchor points. I couldn't find a built-in function to do that so I went the good old loop way. Feel free to optimize if you find some way.

%% // recreate index vectors with all the necessary anchor points
colX(1,1) = xd(1) ;
colY(1,1) = yd(1) ;
k = 2 ;
for t=2:size(yd,1)
    inc = sign(yd(t)-yd(t-1)) ;         %// define increment of 1 with correct sign
    ids = yd(t-1):inc:yd(t)-inc ;       %// range of index to be covered
    if isempty(ids) ; ids = yd(t) ; end %// catch corner cases

    %// now add these points to the list
    n = length(ids) ;
    colX(k:k+n-1,1) = xd(t-1) ;
    colY(k:k+n-1,1) = ids ;
    k=k+n;
end
%// close profile (add last point in the cases where it is necessary)
if inc ~= 0
    colX(end+1,1) = xd(end) ;
    colY(end+1,1) = yd(end) ;
end

%% // display fully rectified jagged circle
npt = size(colX,1) ;
xcj = zeros(npt) ; ycj = xcj ; zcj = xcj ;
for k=1:npt
    xcj(k) = X( colX(k) , colY(k) ) ;
    ycj(k) = Y( colX(k) , colY(k) ) ;
    zcj(k) = Z( colX(k) , colY(k) ) ;
end

And now we don't have gaps between two column ID any more, all the points of the jagged circle match a point of the surface profile: