how to draw a 3d torus

Update:

This new solution is a bit longer, but overcomes all 3 difficulties I listing for the previous solution (below). In particular, this is a vector graphic.

Since you lot expressed an interest in imitating a 3D helix with hidden lines, I also produced an instance of that (in grayscale vector graphics with no shading):

To compile either of these, first salvage the following code in a file called surfacepaths.asy:

          import graph3; import contour;  // A bunch of auxiliary functions.  real fuzz = .001;  real umin(surface s) { return 0; } real vmin(surface s) { render 0; } pair uvmin(surface s) { return (umin(s), vmin(due south)); } real umax(surface due south, existent fuzz=fuzz) {   if (south.ucyclic()) render due south.alphabetize.length;   else render s.index.length - fuzz; } existent vmax(surface southward, real fuzz=fuzz) {   if (south.vcyclic()) return s.index[0].length;   return s.index[0].length - fuzz; } pair uvmax(surface s, real fuzz=fuzz) { render (umax(s,fuzz), vmax(s,fuzz)); }  typedef existent part(real, real);  function normalDot(surface southward, triple eyedir(triple)) {   real toreturn(existent u, real five) {     return dot(s.normal(u, five), eyedir(southward.point(u,5)));   }   return toreturn; }  struct patchWithCoords {   patch p;   existent u;   real 5;   void operator init(patch p, real u, real v) {     this.p = p;     this.u = u;     this.5 = v;   }   void operator init(surface due south, real u, existent v) {     int U=flooring(u);     int V=floor(v);     int alphabetize = (southward.index.length == 0 ? U+5 : s.alphabetize[U][V]);      this.p = s.s[index];     this.u = u-U;     this.5 = v-V;   }   triple partialu() {     render p.partialu(u,five);   }   triple partialv() {     return p.partialv(u,v);   } }  triple[] derivative(surface s, pair pt) {   patchWithCoords thepatch = patchWithCoords(due south, pt.x, pt.y);   return new triple[] {thepatch.partialu(), thepatch.partialv()}; }  typedef triple paramsurface(pair);  paramsurface tangentplane(surface s, pair pt) {   patchWithCoords thepatch = patchWithCoords(southward, pt.x, pt.y);   triple partialu = thepatch.partialu();   triple partialv = thepatch.partialv();   return new triple(pair tangentvector) {     return s.signal(pt.ten, pt.y) + (tangentvector.x * partialu) + (tangentvector.y * partialv);   }; }  guide[] normalpathuv(surface due south, project 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(due south, eyedir), uvmin(south), uvmax(southward), new real[] {0}, nx=n)[0]; }  path3 onSurface(surface due south, path p) {   triple f(int t) {     pair point = point(p,t);     return south.bespeak(point.x, bespeak.y);   }    guide3 toreturn = f(0);   paramsurface thetangentplane = tangentplane(s, indicate(p,0));   triple oldcontrol, newcontrol;   int size = length(p);   for (int i = 1; i <= size; ++i) {     oldcontrol = thetangentplane(postcontrol(p,i-ane) - bespeak(p,i-i));     thetangentplane = tangentplane(s, signal(p,i));     newcontrol = thetangentplane(precontrol(p, i) - point(p,i));     toreturn = toreturn .. controls oldcontrol and newcontrol .. f(i);   }    if (cyclic(p)) toreturn = toreturn & cycle;    render toreturn;  }  /*  * This method returns an assortment of paths that trace out all the  * points on s at which south is parallel to eyedir.  */  path[] paramSilhouetteNoEdges(surface s, projection P = currentprojection, int n = ngraph) {    guide[] uvpaths = normalpathuv(south, 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) {render point(uvpaths[i],t);}, 0, len, n=northward);   }   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(cantankerous(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 < g);   if (i >= m) {     write("problem: Unable to comply.");     write(" dir = " + (string)dir);     write(" surface normal = " + (cord)surfacen);   }   return new triple[] {v1, v2}; }  function3 planeEqn(triple pt, triple normal) {   render new real(triple r) {     render dot(normal, r - pt);   }; }  function2 pullback(function3 eqn, surface s) {   render new real(real u, real v) {     render eqn(s.point(u,five));   }; }  /*  * 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 = south.point(parampt.x, parampt.y);   triple[] lineNormals = normalVectors(dir, s.normal(parampt.10, parampt.y));   path[][] curves;   for (triple northward : lineNormals) {     function3 planeEn = planeEqn(pt, north);     function2 pullback = pullback(planeEn, s);     guide[] contour = contour(pullback, uvmin(south), uvmax(s), new real[]{0})[0];      curves.push(profile);   }   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 button(s.point(P.x, P.y));   return toreturn; }    /*  * Returns those intersection points for which the vector from pt forms an  * astute angle with dir.  */  int numPointsInDirection(surface s, pair parampt, triple dir, real fuzz=.05) {   triple pt = s.point(parampt.x, parampt.y);   dir = unit of measurement(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(existent t0, real t1) {   if (t0 < t1) return true;   if (t0 > t1) return fake;   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.button(f.length - 1);   toreturn.circadian = cyclic;   return toreturn; }  int[] extremes(path path, existent f(pair) = new existent(pair P) {return P.10;}) {   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();   render 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[-i], 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, existent f(pair P) = new real(pair P) {render P.10;}) {   path[] toreturn;   for (path path : paths) {     toreturn.append(splitAtExtremes(path, f));   }   render toreturn; }  path3 toCamera(triple p, projection P=currentprojection, real fuzz = .01, existent upperLimit = 100) {   if (!P.infinity) {     triple directionToCamera = unit of measurement(P.camera - p);     triple startingPoint = p + fuzz*directionToCamera;     return startingPoint -- P.camera;   }   else {     triple directionToCamera = unit(P.camera);     triple startingPoint = p + fuzz*directionToCamera;      render 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(due south, 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) {     describe(path, p=weight*pen + (1-weight)*white);   } } coloredPath[][] layeredPaths; // onTop indicates whether the path should be added at the elevation 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 - i; layer >= 0; --layer) {     real layerfactor = (ane/3)^layer;     for (coloredPath toDraw : layeredPaths[layer]) {       toDraw.draw(layerfactor);     }   }   layeredPaths.delete(); }  existent[] cutTimes(path tocut, path[] knives) {   real[] intersectionTimes = new real[] {0, length(tocut)};   for (path knife : knives) {     real[][] complexIntersections = intersections(tocut, knife);     for (existent[] times : complexIntersections) {       intersectionTimes.push button(times[0]);     }   }   return sort(intersectionTimes); }  path[] cut(path tocut, path[] knives) {   existent[] cutTimes = cutTimes(tocut, knives);   path[] toreturn;   for (int i = 0; i + one < cutTimes.length; ++i) {     toreturn.push button(subpath(tocut,cutTimes[i], cutTimes[i+ane]));   }   return toreturn; }  real[] condense(existent[] values, real fuzz=.001) {   values = sort(values);   values.push button(infinity);   real previous = -infinity;   real lastMin;   real[] toReturn;   for (real t : values) {     if (t - fuzz > previous) {       if (previous > -infinity) toReturn.push button((lastMin + previous) / 2);       lastMin = t;     }     previous = t;   }   return toReturn; }  /*  * A smooth surface parametrized past the domain [0,ane] 10 [0,ane]  */ struct SmoothSurface {   surface due south;   private existent sumax;   private real svmax;   path[] paramSilhouette;   path[] projectedSilhouette;   projection theProjection;    triple photographic camera;    path3 onSurface(path paramPath) {     return onSurface(s, scale(sumax,svmax)*paramPath);   }    triple point(existent u, real v) { return s.point(sumax*u, svmax*5); }   triple point(pair uv) { return signal(uv.x, uv.y); }   triple normal(real u, existent 5) { return s.normal(sumax*u, svmax*v); }   triple normal(pair uv) { render normal(uv.x, uv.y); }    triple[] derivative(real u, real v) { return derivative(s, (u/sumax, v/svmax)); }   triple[] derivative(pair uv) { return derivative(uv.x, uv.y); }    pen colour(real u, real 5, material chiliad, low-cal lite) {     return colour(normal=normal(u,v), m=m, low-cal=light);   }   pen color(pair uv, material thou, light low-cal) {     return this.colour(uv.10, uv.y, m, calorie-free);   }   pair projection(triple pt) { return project(pt, theProjection); }    void operator init(surface south, projection P=currentprojection) {     this.southward = s;     this.sumax = umax(s);     this.svmax = vmax(s);     this.theProjection = P;     this.paramSilhouette = scale(1/sumax, 1/svmax) * paramSilhouetteNoEdges(southward,P);     this.projectedSilhouette = sequence(new path(int i) {     path3 truePath = onSurface(paramSilhouette[i]);     path projectedPath = projection(truePath, theProjection, ninterpolate=1);     return projectedPath;       }, paramSilhouette.length);     camera=theProjection.photographic camera;     if(theProjection.infinity) {       triple m=min(s);       triple M=max(due south);       camera = theProjection.target     + camerafactor*(abs(M-m)+abs(yard-theProjection.target))*unit(theProjection.vector());     }    }    int numSheetsHiding(pair parampt) {     return numSheetsHiding(s, scale(sumax,svmax)*parampt);   }   int numSheetsHiding(real u, real five) {     render numSheetsHiding((u,five));   }    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 in that location are no      * self-intersections of whatever 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,finish));     paramSplitSilhouette.push(subpath(currentParam, outset, end));       }     }      /*      * At present, split at intersections of distinct subpaths.      */      for (int j = 0; j < splitSilhouette.length; ++j) {       path current = splitSilhouette[j];       path currentParam = paramSplitSilhouette[j];        existent[] splittingTimes = new existent[] {0,length(current)};       for (int k = 0; k < splitSilhouette.length; ++one thousand) {     if (j == k) keep;     existent[][] times = intersections(current, splitSilhouette[chiliad]);     for (real[] fourth dimension : times) {       real relevantTime = time[0];       if (.01 < relevantTime && relevantTime < length(current) - .01) splittingTimes.button(relevantTime);     }       }       splittingTimes = sort(splittingTimes);       for (int i = 0; i + 1 < splittingTimes.length; ++i) {     existent first = splittingTimes[i];     real end = splittingTimes[i+1];     real mid = starting time + ((stop-offset) / (2+0.1*unitrand()));     pair theparampoint = point(currentParam, mid);     int sheets = numSheetsHiding(theparampoint);      if (sheets == 0 || includePathsBehind) {       path currentSubpath = subpath(current, starting time, end);       addPath(currentSubpath, p=p, onTop=onTop, layer=sheets);     }        }     }   }    /*     Splits a parametrized path along the parametrized silhouette,     taking [0,1] x [0,1] every bit 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.ten); u < max.x + .001; ++u) {       for (int 5 = floor(min.y); five < max.y + .001; ++five) {     knives.append(shift(u,5)*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;     existent 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;   }    path3 drawSurfacePath(path parampath, pen p=currentpen, bool onTop=true) {     path[] toDraw;     existent[] crossingTimes = splitTimes(parampath);     crossingTimes.append(silhouetteCrossingTimes(onSurface(parampath)));     crossingTimes = condense(crossingTimes);     for (int i = 0; i+1 < crossingTimes.length; ++i) {       toDraw.button(subpath(parampath, crossingTimes[i], crossingTimes[i+1]));     }     for (path thepath : toDraw) {       pair midpoint = point(thepath, length(thepath) / (2+.i*unitrand()));       int sheets = numSheetsHiding(midpoint);       path path3d = project(onSurface(thepath), theProjection, ninterpolate = 1);       addPath(path3d, p=p, onTop=onTop, layer=sheets);     }     render onSurface(parampath);   }  }  SmoothSurface operator *(transform3 t, SmoothSurface southward) {   return SmoothSurface(t*south.s); }                  

To get the first image (with the two circles), in the same directory, run the following *.asy code:

          settings.outformat="pdf"; import surfacepaths;  unitsize(1.5cm); triple eye = 4*(0,1,3.5); currentprojection=perspective(eye, up=Y);  surface torus = surface(O, shift(2X)*path3(calibration(0.4)*unitcircle), axis=Y); SmoothSurface Torus = SmoothSurface(torus);  Torus.drawSilhouette(); real unequal = .07; real midpoint = .75; path circle1 = project(Torus.drawSurfacePath((midpoint-unequal, 0) -- (midpoint-diff, 1)), ninterpolate=i); path circle2 = projection(Torus.drawSurfacePath((midpoint+unequal,0) -- (midpoint+diff,1)),ninterpolate=one);  drawLayeredPaths();    // Depict the paths (with hidden path removal) that were computed in the previous code block.  real circle1mintime = mintimes(circle1)[ane];    // the time at which circle1 reaches its minimum vertical extent pair p1 = betoken(circle1, circle1mintime); characterization("$\mathcal C$", position=p1, align=SW);  existent circle2mintime = mintimes(circle2)[1]; pair p2 = point(circle2, circle2mintime); label("$\mathcal C'$", position=p2, align=SE);  unequal += .01; path3 arc = Torus.onSurface((midpoint-diff,.75) -- (midpoint+diff,.75)); draw(shift(0,-one)*project(arc), arrow=ArcArrow());                  

To become the torus wrapped in a helix, run the following Asymptote lawmaking:

          settings.outformat="pdf"; import surfacepaths;  size(400,0); currentprojection=perspective(camera=(0,4,14), up=Y);  //function to return a line, but split into northward parts (for better approximation when the line is transformed to lie on the surface) path splitline(pair a, pair b, int n) {   pair f(real t) { return (1-t)*a + t*b; }   return graph(f, 0, 1, n=n); }  surface torus = surface(O, Circle(c=2X,r=0.4,normal=Z,n=32), centrality=Y, north=32); /* Let Asymptote know (in example information technology did not already) that the u and 5  * coordinates of this surface are both cyclic.  */ torus.ucyclic(true); torus.vcyclic(truthful); SmoothSurface Torus = SmoothSurface(torus);  Torus.drawSilhouette(includePathsBehind=truthful); Torus.drawSurfacePath(splitline((0,0), (ane,20), n=100), p=greyness);  drawLayeredPaths();                  

---

Onetime (rasterized-simply) solution:

Hither'southward an Asymptote solution that sort of works.

Note that as currently formulated, it has the following difficulties:

1. The torus must be opaque.
2. Anything that should be partially obscured by the torus has to be given in rasterized output (although I've got a scaling factor included to requite high resolution). The rest of the image (i.e., the labels and the arrow) is vector output.
3. Information technology only works with an orthographic projection.

Anyway, here's the code:

          import graph3; import contour;  // A bunch of auxiliary functions for drawing silhouettes.  real fuzz = .001;  real umin(surface south) { render 0; } real vmin(surface south) { return 0; } pair uvmin(surface due south) { return (umin(s), vmin(due south)); } existent umax(surface south, real fuzz=fuzz) {   if (south.ucyclic()) render s.index.length;   else return s.index.length - fuzz; } existent vmax(surface southward, existent fuzz=fuzz) {   if (due south.vcyclic()) return s.alphabetize[0].length;   render due south.index[0].length - fuzz; } pair uvmax(surface s, existent fuzz=fuzz) { return (umax(southward,fuzz), vmax(south,fuzz)); }  typedef real function(real, real);  function normalDot(surface southward, triple eyedir) {   existent toreturn(real u, real v) {     render dot(southward.normal(u, v), eyedir);   }   return toreturn; }  guide[] normalpathuv(surface southward, triple eyedir) {   return profile(normalDot(southward, eyedir), uvmin(s), uvmax(s), new existent[] {0})[0]; }  path3 onSurface(surface s, path p) {   triple f(real t) {     pair point = betoken(p,t);     return due south.point(point.ten, indicate.y);   }   if (circadian(p)) {     guide3 toreturn = f(0);     for (int i = ane; i < size(p); ++i)       toreturn = toreturn -- f(i);     toreturn = toreturn -- cycle;     return toreturn;   }   return graph(f, 0, length(p)); }  /*  * This method returns an assortment of paths that trace out all the  * points on s at which southward is normal to eyedir.  */ path3[] silhouetteNormal(surface southward, triple eyedir) {   guide[] uvpaths = normalpathuv(due south, eyedir);   path3[] toreturn = new path3[uvpaths.length];   for (int i = 0; i < uvpaths.length; ++i) {     toreturn[i] = onSurface(due south, uvpaths[i]);   }   return toreturn; }    /*  * At present, let's actually draw the picture.  */  settings.outformat="pdf"; int resolutionfactor = 4; settings.return=2*resolutionfactor; settings.china=false;  unitsize(1.5cm); triple center = (0,1,iii); currentprojection=orthographic(heart, up=Y);  surface torus = surface(O, shift(2X)*path3(scale(0.five)*unitcircle), axis=Y);  draw(silhouetteNormal(torus, eye)); depict(torus, surfacepen=emissive(white));  existent umin = umin(torus); real umax = umax(torus); real diff = 0.07; path3 circle1 = torus.uequals((.75-diff)*umax); draw(circle1); pair p1 = project(bespeak(circle1, 3)); label("$\mathcal C$", position=p1, align=2SW);  path3 circle2 = torus.uequals((.75+diff)*umax); draw(circle2); pair p2 = project(point(circle2,3)); label("$\mathcal C'$", position=p2, align=2SE);  pair midpoint = (p1 + p2)/two + (0,-.1); path arc = p1 .. midpoint .. p2;  draw(shift(0,-one)*arc, arrow=ArcArrow());  shipout(scale(resolutionfactor)*currentpicture.fit());                  

lunatheable.blogspot.com

Source: https://tex.stackexchange.com/questions/119386/drawing-the-contour-of-a-3d-torus

Share this post