# How to cut a surface at the intersection with another surface in Asymptote?

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);
``````

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

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) {
return;
}

possibleInt region = region(P);
if (region.holds) {
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);

}

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;
}

// 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));
}
}

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

/***************************************/
/* 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;
}
}

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

(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

• 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 Mar 8 '17 at 15:43
• @GMueller: Thanks for pointing this out. I think I've fixed it. Try it now. – Charles Staats Mar 9 '17 at 6:05

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);
``````