14

In the code below, I define a cylinder (as a surface of revolution) and a plane (as a cyclic path). From both, I can get a surface that renders in 3D. Now, I would like to hide (or better remove) the part of the cylinder, which is below the surface (in this example at z < 0, rendered in a darker shade of red), while the plane should remain semi-transparent.

How can this be done elegantly?

import three;
import solids;

currentprojection = obliqueY();

path3 xyplane = path3(scale(10) * box((-1,-1),(1,1)));
revolution c =  rotate(-45,Y) * shift((0,0,-5)) *cylinder(O,1,15);

draw(surface(xyplane),black+opacity(.5));
draw(xyplane,black+linewidth(.1));

draw(surface(c),red);
draw(c,red);

Rendered image

3
  • Why not drawing a white opaque surface before to draw the semi-transparent one? Commented Dec 14, 2012 at 13:03
  • @LionelMANSUY This will be part of a complex scene. I need to draw another cylinder below the xy-plane.
    – Jan
    Commented Dec 14, 2012 at 14:01
  • You could perhaps arrange the order of all the drawings? Commented Dec 14, 2012 at 14:56

2 Answers 2

12

Here's a way to cut a surface by another surface that is defined by an equation.

First, save the following code in a file called crop3D.asy:

import three;

/**********************************************/
/* Code for splitting surfaces: */

struct possibleInt {
  int value;
  bool holds;
}

// Get versions of hsplit and vsplit with no extra optional
// argument.
triple[][][] old_hsplit(triple[][] P) { return hsplit(P); }
triple[][][] old_vsplit(triple[][] P) { return vsplit(P); }

int operator cast(possibleInt i) { return i.value; }

restricted int maxdepth = 20;
restricted void maxdepth(int n) { maxdepth = n; }

surface[] divide(surface s, possibleInt region(patch), int numregions,
         bool keepregion(int) = null) {

  if (keepregion == null) keepregion = new bool(int region) {
      return (0 <= region && region < numregions);
    };

  surface[] toreturn = new surface[numregions];
  for (int i = 0; i < numregions; ++i)
    toreturn[i] = new surface;

  void addPatch(patch P, int region) {
    if (keepregion(region)) toreturn[region].push(P);
  }

  void divide(patch P, int depth) {
    if (depth == 0) {
      addPatch(P, region(P));
      return;
    }

    possibleInt region = region(P);
    if (region.holds) {
      addPatch(P, region);
      return;
    }

    // Choose the splitting function based on the parity of the recursion depth.
    triple[][][] Split(triple[][] P) =
        (depth % 2 == 0 ? old_hsplit : old_vsplit);

    patch[] Split(patch P) {
      triple[][][] patches = Split(P.P);
      return sequence(new patch(int i) {return patch(patches[i]);}, patches.length);
    }

    patch[] patches = Split(P);
    for (patch PP : patches)
      divide(PP, depth-1);
  }

  for (patch P : s.s)
    divide(P, maxdepth);

  return toreturn;
}


surface[] divide(surface s, int region(triple), int numregions,
         bool keepregion(int) = null) {
  possibleInt patchregion(patch P) {
    triple[][] controlpoints = P.P;
    possibleInt theRegion;
    theRegion.value = region(controlpoints[0][0]);
    theRegion.holds = true;
    for (triple[] ta : controlpoints) {
      for (triple t : ta) {
    if (region(t) != theRegion.value) {
      theRegion.holds = false;
      break;
    }
      }
      if (!theRegion.holds) break;
    }
    return theRegion;
  }

  return divide(s, patchregion, numregions, keepregion);
}


/**************************************************/
/* Code for cropping surfaces */

// Return 0 iff the point lies in box(a,b).
int cropregion(triple pt, triple a=O, triple b=(1,1,1)) {
  real x=pt.x, y=pt.y, z=pt.z;
  int toreturn=0;
  real xmin=a.x, xmax=b.x, ymin = a.y, ymax=b.y, zmin=a.z, zmax=b.z;
  if (xmin > xmax) { xmin = b.x; xmax = a.x; }
  if (ymin > ymax) { ymin = b.y; ymax = a.y; }
  if (zmin > zmax) { zmin = b.z; zmax = a.z; }
  if (x < xmin) --toreturn;
  else if (x > xmax) ++toreturn;
  toreturn *= 2;
  if (y < ymin) --toreturn;
  else if (y > ymax) ++toreturn;
  toreturn *= 2;
  if (z < zmin) --toreturn;
  else if (z > zmax) ++toreturn;
  return toreturn;
}

// Crop the surface to box(a,b).
surface crop(surface s, triple a, triple b) {
  int region(triple pt) {
    return cropregion(pt, a, b);
  }
  return divide(s, region=region, numregions=1)[0];
}

// Crop the surface to things contained in a region described by a bool(triple) function
surface crop(surface s, bool allow(triple)) {
  int region(triple pt) {
    if (allow(pt)) return 0;
    else return -1;
  }
  return divide(s, region=region, numregions=1)[0];
}

/******************************************/
/* Code for cropping paths */

// A rectangular solid with opposite vertices a, b:
surface surfacebox(triple a, triple b) {
  return shift(a)*scale((b-a).x,(b-a).y,(b-a).z)*unitcube;
}

bool containedInBox(triple pt, triple a, triple b) {
  return cropregion(pt, a, b) == 0;
}

// Crop a path3 to box(a,b).
path3[] crop(path3 g, triple a, triple b) {
  surface thebox = surfacebox(a,b);
  path3[] toreturn;
  real[] times = new real[] {0};
  real[][] alltimes = intersections(g, thebox);
  for (real[] threetimes : alltimes)
    times.push(threetimes[0]);
  times.push(length(g));
  for (int i = 1; i < times.length; ++i) {
    real mintime = times[i-1];
    real maxtime = times[i];
    triple midpoint = point(g, (mintime+maxtime)/2);
    if (containedInBox(midpoint, a, b))
      toreturn.push(subpath(g, mintime, maxtime));
  }
  return toreturn;
}

path3[] crop(path3[] g, triple a, triple b) {
  path3[] toreturn;
  for (path3 gi : g)
    toreturn.append(crop(gi, a, b));
  return toreturn;
}

/***************************************/
/* Code to return only the portion of the surface facing the camera */

bool facingCamera(triple vec, triple pt=O, projection P = currentprojection, bool towardsCamera = true) {
  triple normal = P.camera;
  if (!P.infinity) {
    normal = P.camera - pt;
  }
  if (towardsCamera) return (dot(vec, normal) >= 0);
  else return (dot(vec, normal) <= 0);
}

surface facingCamera(surface s, bool towardsCamera = true, int maxdepth = 10) {

  int oldmaxdepth = maxdepth;
  maxdepth(maxdepth);

  possibleInt facingregion(patch P) {
    int n = 2;
    possibleInt toreturn;
    unravel toreturn;
    bool facingcamera = facingCamera(P.normal(1/2, 1/2), pt=P.point(1/2,1/2), towardsCamera);
    value = facingcamera ? 0 : 1;
    holds = true;
    for (int i = 0; i <= n; ++i) {
      real u = i/n;
      for (int j = 0; j <= n; ++j) {
    real v = j/n;
    if (facingCamera(P.normal(u,v), P.point(u,v), towardsCamera) != facingcamera) {
      holds = false;
      break;
    }
      }
      if (!holds) break;
    }
    return toreturn;
  }

  surface toreturn = divide(s, facingregion, numregions=1)[0];
  maxdepth(oldmaxdepth);
  return toreturn;
}

(This is essentially defining a new module; only a portion of the code is actually required for this example.) Then, run the Asymptote code

settings.outformat="png";
settings.render=16;

import three;
import solids;
import crop3D;

currentprojection = obliqueY();

path3 xyplane = path3(scale(10) * box((-1,-1),(1,1)));
surface c =  surface( rotate(-45,Y) * shift((0,0,-5)) * cylinder(O,1,15) );

bool zpositive(triple pt) { return pt.z > 0; }
c = crop(c, zpositive);

draw(surface(xyplane),black+opacity(.5));
draw(xyplane,black+linewidth(.1));

draw(c,red);

to produce the image

2
  • Currently your example throws me errors: crop3D.asy: 44.67: cannot cast 'triple[][][](triple[][] P, real u=<default>)' to 'triple[][][](triple[][] P)' crop3D.asy: 44.58: cannot cast 'triple[][][](triple[][] P, real v=<default>)' to 'triple[][][](triple[][] P)' Unfortunately I am not expert enough to fix this - any help would be much appreciated! This is with asymptote v2.38
    – GPMueller
    Commented Mar 8, 2017 at 15:43
  • 1
    @GMueller: Thanks for pointing this out. I think I've fixed it. Try it now. Commented Mar 9, 2017 at 6:05
7

Here's the workaround I came up with for this particular example. If the surface, you want to cut is a surface of revolution, you can define it as a parametric surface excluding the part below the plane. However, this is not a generic answer to the question.

import three;
import solids;

currentprojection = obliqueY();

path3 xyplane = path3(scale(10) * box((-1,-1),(1,1)));

surface cylinderSurfaceTiltedPlane(real R, real z0, real planeAlpha, real planePhi) {
    triple parametricCylinder(pair p) {
        real Phi = p.x;
        real Z = p.y;
        real x = R*sin(Phi);
        real y = R*cos(Phi);
        real z = Z * (tan(planeAlpha*pi/180)*sin(Phi - planePhi*pi/180)*R + z0)/z0;
        return (x,y,z);
    }
    return surface(parametricCylinder, (0,0), (2pi,z0), Spline);
}


draw(surface(xyplane),black+opacity(.5));
draw(xyplane,black+linewidth(.1));

surface cF = rotate(180-45,Y)  * shift((0,0,-10)) *     cylinderSurfaceTiltedPlane(1, 10, -45, 0);
draw(cF,red);

rendered image

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .