16

I would like to show a Torus with three paths on its surface:

enter image description here

  • Both points ($a$ and $b$) are on the surface of the torus.
  • All three paths (orange, red, green) are on the surface of the torus
  • All three paths start in $a$ and end in $b$
  • The orange and the red path are pretty "straight" (one goes the left side, one the right side, but the green path is more interesting. It makes a "curve" (I really don't know how to describe this better.

The torus was created with sketch:

def torus {
    def n_segs 40
    sweep [draw=black, fill=lightgray, fill opacity=0.75] {n_segs, rotate(360/n_segs, (0,0,0), [0,1,0])}
        sweep {n_segs, rotate(360/n_segs, (1.5,0,0), [0,0,1])}
        (2,0,0)
}

put { view((10,4,2)) } {{torus}}

global { language tikz }

How can I print the paths on the torus?

  • 1
    Check out the related )to the right). Expccially tex.stackexchange.com/questions/70090/… – John Kormylo Dec 14 '13 at 19:46
  • Dou you have a parametrization of your curves? I would like to try to draw that directly in tikz tomorrow. – Ronny Dec 23 '13 at 22:43
  • @Ronny: No, I don't have a parametrization of my curves. – Martin Thoma Dec 24 '13 at 16:03
19
+50

Here is an answer that uses Asymptote to produce a vector graphic result: enter image description here

The actual pdf file may be found, for now, at this location; but I don't think you will have any trouble compiling it (although it may take a while--about 77 seconds on my computer). I've omitted the LaTeX wrapper and divided the code into two code blocks for readability, but you can just copy and paste them one after the other to make a coherent .asy file.

The first code block is really an Asymptote module I am writing that is in its very early stages:

settings.outformat="pdf";

import graph3;
import contour;

// A bunch of auxiliary functions.

real fuzz = .001;

real umin(surface s) { return 0; }
real vmin(surface s) { return 0; }
pair uvmin(surface s) { return (umin(s), vmin(s)); }
real umax(surface s, real fuzz=fuzz) {
  if (s.ucyclic()) return s.index.length;
  else return s.index.length - fuzz;
}
real vmax(surface s, real fuzz=fuzz) {
  if (s.vcyclic()) return s.index[0].length;
  return s.index[0].length - fuzz;
}
pair uvmax(surface s, real fuzz=fuzz) { return (umax(s,fuzz), vmax(s,fuzz)); }

typedef real function(real, real);

function normalDot(surface s, triple eyedir(triple)) {
  real toreturn(real u, real v) {
    return dot(s.normal(u, v), eyedir(s.point(u,v)));
  }
  return toreturn;
}

struct patchWithCoords {
  patch p;
  real u;
  real v;
  void operator init(patch p, real u, real v) {
    this.p = p;
    this.u = u;
    this.v = v;
  }
  void operator init(surface s, real u, real v) {
    int U=floor(u);
    int V=floor(v);
    int index = (s.index.length == 0 ? U+V : s.index[U][V]);

    this.p = s.s[index];
    this.u = u-U;
    this.v = v-V;
  }
  triple partialu() {
    return p.partialu(u,v);
  }
  triple partialv() {
    return p.partialv(u,v);
  }
}

typedef triple paramsurface(pair);

paramsurface tangentplane(surface s, pair pt) {
  patchWithCoords thepatch = patchWithCoords(s, pt.x, pt.y);
  triple partialu = thepatch.partialu();
  triple partialv = thepatch.partialv();
  return new triple(pair tangentvector) {
    return s.point(pt.x, pt.y) + (tangentvector.x * partialu) + (tangentvector.y * partialv);
  };
}

guide[] normalpathuv(surface s, projection P = currentprojection, int n = ngraph) {
  triple eyedir(triple a);
  if (P.infinity) eyedir = new triple(triple) { return P.camera; };
  else eyedir = new triple(triple pt) { return P.camera - pt; };
  return contour(normalDot(s, eyedir), uvmin(s), uvmax(s), new real[] {0}, nx=n)[0];
}

path3 onSurface(surface s, path p) {
  triple f(int t) {
    pair point = point(p,t);
    return s.point(point.x, point.y);
  }

  guide3 toreturn = f(0);
  paramsurface thetangentplane = tangentplane(s, point(p,0));
  triple oldcontrol, newcontrol;
  int size = length(p);
  for (int i = 1; i <= size; ++i) {
    oldcontrol = thetangentplane(postcontrol(p,i-1) - point(p,i-1));
    thetangentplane = tangentplane(s, point(p,i));
    newcontrol = thetangentplane(precontrol(p, i) - point(p,i));
    toreturn = toreturn .. controls oldcontrol and newcontrol .. f(i);
  }

  if (cyclic(p)) toreturn = toreturn & cycle;

  return toreturn;

}

/*
 * This method returns an array of paths that trace out all the
 * points on s at which s is parallel to eyedir.
 */

path[] paramSilhouetteNoEdges(surface s, projection P = currentprojection, int n = ngraph) {
   guide[] uvpaths = normalpathuv(s, P, n);
  //Reduce the number of segments to conserve memory
  for (int i = 0; i < uvpaths.length; ++i) {
    real len = length(uvpaths[i]);
    uvpaths[i] = graph(new pair(real t) {return point(uvpaths[i],t);}, 0, len, n=n);
  }
  return uvpaths;
}   

private typedef real function2(real, real);
private typedef real function3(triple);

triple[] normalVectors(triple dir, triple surfacen) {
  dir = unit(dir);
  surfacen = unit(surfacen);
  triple v1, v2;
  int i = 0;
  do {
    v1 = unit(cross(dir, (unitrand(), unitrand(), unitrand())));
    v2 = unit(cross(dir, (unitrand(), unitrand(), unitrand())));
    ++i;
  } while ((abs(dot(v1,v2)) > Cos(10) || abs(dot(v1,surfacen)) > Cos(5) || abs(dot(v2,surfacen)) > Cos(5)) && i < 1000);
  if (i >= 1000) {
    write("problem: Unable to comply.");
    write(" dir = " + (string)dir);
    write(" surface normal = " + (string)surfacen);
  }
  return new triple[] {v1, v2};
}

function3 planeEqn(triple pt, triple normal) {
  return new real(triple r) {
    return dot(normal, r - pt);
  };
}

function2 pullback(function3 eqn, surface s) {
  return new real(real u, real v) {
    return eqn(s.point(u,v));
  };
}

/*
 * returns the distinct points in which the surface intersects
 * the line through the point pt in the direction dir
 */

triple[] intersectionPoints(surface s, pair parampt, triple dir) {
  triple pt = s.point(parampt.x, parampt.y);
  triple[] lineNormals = normalVectors(dir, s.normal(parampt.x, parampt.y));
  path[][] curves;
  for (triple n : lineNormals) {
    function3 planeEn = planeEqn(pt, n);
    function2 pullback = pullback(planeEn, s);
    guide[] contour = contour(pullback, uvmin(s), uvmax(s), new real[]{0})[0];

    curves.push(contour);
  }
  pair[] intersectionPoints;
  for (path c1 : curves[0])
    for (path c2 : curves[1])
      intersectionPoints.append(intersectionpoints(c1, c2));
  triple[] toreturn;
  for (pair P : intersectionPoints)
    toreturn.push(s.point(P.x, P.y));
  return toreturn;
}



/*
 * Returns those intersection points for which the vector from pt forms an
 * acute angle with dir.
 */
 int numPointsInDirection(surface s, pair parampt, triple dir, real fuzz=.05) {
  triple pt = s.point(parampt.x, parampt.y);
  dir = unit(dir);
  triple[] intersections = intersectionPoints(s, parampt, dir);
  int num = 0;
  for (triple isection: intersections)
    if (dot(isection - pt, dir) > fuzz) ++num;
  return num;
}

bool3 increasing(real t0, real t1) {
  if (t0 < t1) return true;
  if (t0 > t1) return false;
  return default;
}

int[] extremes(real[] f, bool cyclic = f.cyclic) {
  bool3 lastIncreasing;
  bool3 nextIncreasing;
  int max;
  if (cyclic) {
    lastIncreasing = increasing(f[-1], f[0]);
    max = f.length - 1;
  } else {
    max = f.length - 2;
    if (increasing(f[0], f[1])) lastIncreasing = false;
    else lastIncreasing = true;
  }
  int[] toreturn;
  for (int i = 0; i <= max; ++i) {
    nextIncreasing = increasing(f[i], f[i+1]);
    if (lastIncreasing != nextIncreasing) {
      toreturn.push(i);
    }
    lastIncreasing = nextIncreasing;
  }
  if (!cyclic) toreturn.push(f.length - 1);
  toreturn.cyclic = cyclic;
  return toreturn;
}

int[] extremes(path path, real f(pair) = new real(pair P) {return P.x;})
{
  real[] fvalues = new real[size(path)];
  for (int i = 0; i < fvalues.length; ++i) {
    fvalues[i] = f(point(path, i));
  }
  fvalues.cyclic = cyclic(path);
  int[] toreturn = extremes(fvalues);
  fvalues.delete();
  return toreturn;
}

path[] splitAtExtremes(path path, real f(pair) = new real(pair P) {return P.x;})
{
  int[] splittingTimes = extremes(path, f);
  path[] toreturn;
  if (cyclic(path)) toreturn.push(subpath(path, splittingTimes[-1], splittingTimes[0]));
  for (int i = 0; i+1 < splittingTimes.length; ++i) {
    toreturn.push(subpath(path, splittingTimes[i], splittingTimes[i+1]));
  }
  return toreturn;
}

path[] splitAtExtremes(path[] paths, real f(pair P) = new real(pair P) {return P.x;})
{
  path[] toreturn;
  for (path path : paths) {
    toreturn.append(splitAtExtremes(path, f));
  }
  return toreturn;
}

path3 toCamera(triple p, projection P=currentprojection, real fuzz = .01, real upperLimit = 100) {
  if (!P.infinity) {
    triple directionToCamera = unit(P.camera - p);
    triple startingPoint = p + fuzz*directionToCamera;
    return startingPoint -- P.camera;
  }
  else {
    triple directionToCamera = unit(P.camera);
    triple startingPoint = p + fuzz*directionToCamera;

    return startingPoint -- (p + upperLimit*directionToCamera);
  }
}

int numSheetsHiding(surface s, pair parampt, projection P = currentprojection) {
  triple p = s.point(parampt.x, parampt.y);
  path3 tocamera = toCamera(p, P);
  triple pt = beginpoint(tocamera);
  triple dir = endpoint(tocamera) - pt;
  return numPointsInDirection(s, parampt, dir);
}

struct coloredPath {
  path path;
  pen pen;
  void operator init(path path, pen p=currentpen) {
    this.path = path;
    this.pen = p;
  }
  /* draws the path with the pen having the specified weight (using colors)*/
  void draw(real weight) {
    draw(path, p=weight*pen + (1-weight)*white);
  }
}
coloredPath[][] layeredPaths;
// onTop indicates whether the path should be added at the top or bottom of the specified layer
void addPath(path path, pen p=currentpen, int layer, bool onTop=true) {
  coloredPath toAdd = coloredPath(path, p);
  if (layer >= layeredPaths.length) {
    layeredPaths[layer] = new coloredPath[] {toAdd};
  } else if (onTop) {
    layeredPaths[layer].push(toAdd);
  } else layeredPaths[layer].insert(0, toAdd);
}

void drawLayeredPaths() {
  for (int layer = layeredPaths.length - 1; layer >= 0; --layer) {
    real layerfactor = (1/3)^layer;
    for (coloredPath toDraw : layeredPaths[layer]) {
      toDraw.draw(layerfactor);
    }
  }
}

real[] cutTimes(path tocut, path[] knives) {
  real[] intersectionTimes = new real[] {0, length(tocut)};
  for (path knife : knives) {
    real[][] complexIntersections = intersections(tocut, knife);
    for (real[] times : complexIntersections) {
      intersectionTimes.push(times[0]);
    }
  }
  return sort(intersectionTimes);
}

path[] cut(path tocut, path[] knives) {
  real[] cutTimes = cutTimes(tocut, knives);
  path[] toreturn;
  for (int i = 0; i + 1 < cutTimes.length; ++i) {
    toreturn.push(subpath(tocut,cutTimes[i], cutTimes[i+1]));
  }
  return toreturn;
}

real[] condense(real[] values, real fuzz=.001) {
  values = sort(values);
  values.push(infinity);
  real previous = -infinity;
  real lastMin;
  real[] toReturn;
  for (real t : values) {
    if (t - fuzz > previous) {
      if (previous > -infinity) toReturn.push((lastMin + previous) / 2);
      lastMin = t;
    }
    previous = t;
  }
  return toReturn;
}

/*
 * A smooth surface parametrized by the domain [0,1] x [0,1]
 */
struct SmoothSurface {
  surface s;
  private real sumax;
  private real svmax;
  path[] paramSilhouette;
  path[] projectedSilhouette;
  projection theProjection;

  path3 onSurface(path paramPath) {
    return onSurface(s, scale(sumax,svmax)*paramPath);
  }

  triple point(real u, real v) { return s.point(sumax*u, svmax*v); }
  triple point(pair uv) { return point(uv.x, uv.y); }
  triple normal(real u, real v) { return s.normal(sumax*u, svmax*v); }
  triple normal(pair uv) { return normal(uv.x, uv.y); }

  void operator init(surface s, projection P=currentprojection) {
    this.s = s;
    this.sumax = umax(s);
    this.svmax = vmax(s);
    this.theProjection = P;
    this.paramSilhouette = scale(1/sumax, 1/svmax) * paramSilhouetteNoEdges(s,P);
    this.projectedSilhouette = sequence(new path(int i) {
    path3 truePath = onSurface(paramSilhouette[i]);
    path projectedPath = project(truePath, theProjection, ninterpolate=1);
    return projectedPath;
      }, paramSilhouette.length);
  }

  int numSheetsHiding(pair parampt) {
    return numSheetsHiding(s, scale(sumax,svmax)*parampt);
  }

  void drawSilhouette(pen p=currentpen, bool includePathsBehind=false, bool onTop = true) {
    int[][] extremes;
    for (path path : projectedSilhouette) {
      extremes.push(extremes(path));
    }

    path[] splitSilhouette;
    path[] paramSplitSilhouette;

    /*
     * First, split at extremes to ensure that there are no
     * self-intersections of any one subpath in the projected silhouette.
     */

    for (int j = 0; j < paramSilhouette.length; ++j) {
      path current = projectedSilhouette[j];

      path currentParam = paramSilhouette[j];

      int[] dividers = extremes[j];
      for (int i = 0; i + 1 < dividers.length; ++i) {
    int start = dividers[i];
    int end = dividers[i+1];
    splitSilhouette.push(subpath(current,start,end));
    paramSplitSilhouette.push(subpath(currentParam, start, end));
      }
    }

    /*
     * Now, split at intersections of distinct subpaths.
     */

    for (int j = 0; j < splitSilhouette.length; ++j) {
      path current = splitSilhouette[j];
      path currentParam = paramSplitSilhouette[j];

      real[] splittingTimes = new real[] {0,length(current)};
      for (int k = 0; k < splitSilhouette.length; ++k) {
    if (j == k) continue;
    real[][] times = intersections(current, splitSilhouette[k]);
    for (real[] time : times) {
      real relevantTime = time[0];
      if (.01 < relevantTime && relevantTime < length(current) - .01) splittingTimes.push(relevantTime);
    }
      }
      splittingTimes = sort(splittingTimes);
      for (int i = 0; i + 1 < splittingTimes.length; ++i) {
    real start = splittingTimes[i];
    real end = splittingTimes[i+1];
    real mid = start + ((end-start) / (2+0.1*unitrand()));
    pair theparampoint = point(currentParam, mid);
    int sheets = numSheetsHiding(theparampoint);

    if (sheets == 0 || includePathsBehind) {
      path currentSubpath = subpath(current, start, end);
      addPath(currentSubpath, p=p, onTop=onTop, layer=sheets);
    }

      }
    }
  }

  /*
    Splits a parametrized path along the parametrized silhouette,
    taking [0,1] x [0,1] as the
    fundamental domain.  Could be implemented more efficiently.
  */
  private real[] splitTimes(path thepath) {
    pair min = min(thepath);
    pair max = max(thepath);
    path[] baseknives = paramSilhouette;
    path[] knives;
    for (int u = floor(min.x); u < max.x + .001; ++u) {
      for (int v = floor(min.y); v < max.y + .001; ++v) {
    knives.append(shift(u,v)*baseknives);
      }
    }
    return cutTimes(thepath, knives);
  }

  /*
    Returns the times at which the projection of the given path3 intersects
    the projection of the surface silhouette. This may miss unstable
    intersections that can be detected by the previous method.
  */
  private real[] silhouetteCrossingTimes(path3 thepath, real fuzz = .01) {
    path projectedpath = project(thepath, theProjection, ninterpolate=1);
    real[] crossingTimes = cutTimes(projectedpath, projectedSilhouette);
    if (crossingTimes.length == 0) return crossingTimes;
    real current = 0;
    real[] toReturn = new real[] {0};
    for (real prospective : crossingTimes) {
      if (prospective > current + fuzz
      && prospective < length(thepath) - fuzz) {
    toReturn.push(prospective);
    current = prospective;
      }
    }
    toReturn.push(length(thepath));
    return toReturn;
  }

  void drawSurfacePath(path parampath, pen p=currentpen, bool onTop=true) {
    path[] toDraw;
    real[] crossingTimes = splitTimes(parampath);
    crossingTimes.append(silhouetteCrossingTimes(onSurface(parampath)));
    crossingTimes = condense(crossingTimes);
    for (int i = 0; i+1 < crossingTimes.length; ++i) {
      toDraw.push(subpath(parampath, crossingTimes[i], crossingTimes[i+1]));
    }
    for (path thepath : toDraw) {
      pair midpoint = point(thepath, length(thepath) / (2+.1*unitrand()));
      int sheets = numSheetsHiding(midpoint);
      path path3d = project(onSurface(thepath), theProjection, ninterpolate = 1);
      addPath(path3d, p=p, onTop=onTop, layer=sheets);
    }
  }
}

The second code block is the code that uses the utilities defined above to actually draw a torus. It bears a certain resemblance to the code from my previous (rasterized-only) answer.

real unit = 4cm;
unitsize(unit);
triple eye = (10,1,4);
//currentprojection=perspective(2*eye);
currentprojection=orthographic(eye);

surface torus = surface(Circle(c=2Y, r=0.6, normal=X, n=32), c=O, axis=Z, n=32);
torus.ucyclic(true);
torus.vcyclic(true);

SmoothSurface Torus = SmoothSurface(torus);

Torus.drawSilhouette(p=black, includePathsBehind=true);

pair a = (22/40, 3/40);
pair b = (5/40, .5);

path abpathparam(int ucycles, int vcycles) {
  pair bshift = (ucycles, vcycles);
  pair f(real t) {
    return (1-t)*a + t*(b+bshift);
  }
  return graph(f, 0, 1, n=10);
}

real linewidth = 0.8pt;

Torus.drawSurfacePath(abpathparam(0,0), p=linewidth + orange);
Torus.drawSurfacePath(abpathparam(1,0), p=linewidth + red);
Torus.drawSurfacePath(abpathparam(1,-1), p=linewidth + (darkgreen + 0.2blue));

pen meshpen = gray(0.6);
for (real u = 0; u < 1; u += 1/40) {
  Torus.drawSurfacePath(graph(new pair(real v) {return (u,v);}, 0,1,n=5), p=meshpen, onTop=false);
}
for (real v = 0; v < 1; v += 1/20) {
  Torus.drawSurfacePath(graph(new pair(real u) {return (u,v);}, 0, 1, n=5), p=meshpen, onTop=false);
}

drawLayeredPaths();

dot(project(Torus.point(a.x,a.y)), L="$a$", align=W);
dot(project(Torus.point(b.x,b.y)), L="$b$", align=N);
  • You are amazing! Thanks for putting so much work in this image (it's now part of my free lecture notes for geometry and topology). In the evening, I will try to compile it. If it works, I'll give you a bounty. (By the way, the image at your university site looks great.) – Martin Thoma Dec 23 '13 at 10:25
  • 1
    When I try to compile this, I get error: out of memory. – Martin Thoma Dec 23 '13 at 17:52
  • @moose: I've put a fair amount of effort into getting rid of those errors, but I guess it makes sense that different installations will behave differently. If you use /*...*/ to comment out all the lines strictly between Torus.drawSilhouette and drawLayeredPaths();, what happens? – Charles Staats Dec 23 '13 at 19:03
  • It compiles (to a quite crappy image) and gives: next.asy: 539.26: no matching variable 'a.x' next.asy: 539.30: no matching variable 'a.y' – Martin Thoma Dec 23 '13 at 21:26
  • 1
    @moose: The example may be found here. The code is here; it requires this auxiliary code to be saved in a file in the same directory. – Charles Staats Feb 20 '14 at 15:29
17

How's this?

The code using asymptote:

\documentclass[margin=1cm]{standalone}
\usepackage{asymptote}
\begin{document}
\begin{asy}
settings.render = 8;
settings.prc = false;

import graph3;
import contour;
size3(8cm);

currentprojection = orthographic(10,1,4);
defaultrender = render(merge = true);

// create torus as surface of rotation
int umax = 40;
int vmax = 40;
surface torus = surface(Circle(c=2Y, r=0.6, normal=X, n=vmax), c=O, axis=Z, n=umax);
torus.ucyclic(true);
torus.vcyclic(true);

pen meshpen = 0.3pt + gray;

draw(torus, surfacepen=material(diffusepen=white+opacity(0.6), emissivepen=white));
for (int u = 0; u < umax; ++u)
  draw(torus.uequals(u), p=meshpen);
for (int v = 0; v < vmax; ++v)
  draw(graph(new triple(real u) {return torus.point(u,v); }, 0, umax, operator ..),
       p=meshpen);

pair a = (floor(umax/2) + 2, 3);
dot(torus.point(a.x, a.y), L="$a$", align=W);
pair b = (5, floor(vmax/2));
dot(torus.point(b.x, b.y), L="$b$", align=2Z + X);

path3 abpath(int ucycles, int vcycles) {
  pair bshift = (ucycles*umax, vcycles*vmax);
  triple f(real t) {
    pair uv = (1-t)*a + t*(b+bshift);
    return torus.point(uv.x, uv.y);
  }
  return graph(f, 0, 1, operator ..);
}

real linewidth = 0.8pt;

draw(abpath(0,0), p=linewidth + orange);
draw(abpath(1,0), p=linewidth + red);
draw(abpath(1,-1), p=linewidth + darkgreen);
\end{asy}
\end{document}
  • A great answer! But when I try to compile it, I get an error: /usr/local/texlive/2013/texmf-dist/asymptote/three.asy: 2906.13: runtime: to support onscreen rendering, please install glut library, run ./configure, and recompile. Currently, I update Texlive and hope it fixes the problem, because I didn't find glut within the repository (I use Linux Mint). – Martin Thoma Dec 15 '13 at 18:03
  • Updating texlive didn't help. freeglut3 is installed. What can I do to fix the problem? – Martin Thoma Dec 15 '13 at 18:16
  • @moose: I wish I knew enough to help, in part because if I did, then I might have some chance of manually updating Asymptote on my own machine. All I know is that for me, after installing MacTeX, Asymptote magically works. – Charles Staats Dec 15 '13 at 20:39
  • @moose: I do have one suggestion, which may or may not be helpful. Make sure that whatever workflow you are using is calling asy without the -V option; the error you are getting may be an error in the preview mechanism rather than the mechanism for producing a pdf file. – Charles Staats Dec 16 '13 at 1:27
  • 1
    @moose: When rendering 3d graphics as vector graphics, Asymptote has much more limited capabilities for hidden surface removal, at least for the moment. Try the same source code with settings.render = 0; instead of 8; you will probably be able to compile it successfully, and the result will be a vector graphic, but the green path will look like it's on top of the torus all the way. – Charles Staats Dec 17 '13 at 0:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.