mitsuba-af602c6fd98a/src/utils/cylclip.cpp
/**
	 * Compute the ellipse created by the intersection of an infinite
	 * cylinder and a plane. Returns false in the degenerate case.
	 * Based on:
	 * www.geometrictools.com/Documentation/IntersectionCylinderPlane.pdf
	 */
Compute the ellipse created by the intersection of an infinite
cylinder and a plane. Returns false in the degenerate case.
Compute the ellipse created by the intersection of an infinite cylinder and a plane . Returns false in the degenerate case .
bool intersectCylPlane(Point planePt, Normal planeNrml,
			Point cylPt, Vector cylD, Float radius, Point &center,
			Vector *axes, Float *lengths) const {
		if (absDot(planeNrml, cylD) < Epsilon)
			return false;
		Vector B, A = cylD - dot(cylD, planeNrml)*planeNrml;
		Float length = A.length();
		if (length != 0) {
			A /= length;
			B = cross(planeNrml, A);
		} else {
			coordinateSystem(planeNrml, A, B);
		}
		Vector delta = planePt - cylPt,
			   deltaProj = delta - cylD*dot(delta, cylD);
		Float aDotD = dot(A, cylD);
		Float bDotD = dot(B, cylD);
		Float c0 = 1-aDotD*aDotD;
		Float c1 = 1-bDotD*bDotD;
		Float c2 = 2*dot(A, deltaProj);
		Float c3 = 2*dot(B, deltaProj);
		Float c4 = dot(delta, deltaProj) - radius*radius;
		Float lambda = (c2*c2/(4*c0) + c3*c3/(4*c1) - c4)/(c0*c1);
		Float alpha0 = -c2/(2*c0),
			  beta0 = -c3/(2*c1);
		lengths[0] = std::sqrt(c1*lambda),
		lengths[1] = std::sqrt(c0*lambda);
		center = planePt + alpha0 * A + beta0 * B;
		axes[0] = A;
		axes[1] = B;
		return true;
	}
bool intersectCylPlane ( Point planePt , Normal planeNrml , Point cylPt , Vector cylD , Float radius , Point & center , Vector * axes , Float * lengths ) const { if ( absDot ( planeNrml , cylD ) < Epsilon ) return false ; Vector B , A = cylD - dot ( cylD , planeNrml ) * planeNrml ; Float length = A . length ( ) ; if ( length != 0 ) { A /= length ; B = cross ( planeNrml , A ) ; } else { coordinateSystem ( planeNrml , A , B ) ; } Vector delta = planePt - cylPt , deltaProj = delta - cylD * dot ( delta , cylD ) ; Float aDotD = dot ( A , cylD ) ; Float bDotD = dot ( B , cylD ) ; Float c0 = 1 - aDotD * aDotD ; Float c1 = 1 - bDotD * bDotD ; Float c2 = 2 * dot ( A , deltaProj ) ; Float c3 = 2 * dot ( B , deltaProj ) ; Float c4 = dot ( delta , deltaProj ) - radius * radius ; Float lambda = ( c2 * c2 / ( 4 * c0 ) + c3 * c3 / ( 4 * c1 ) - c4 ) / ( c0 * c1 ) ; Float alpha0 = - c2 / ( 2 * c0 ) , beta0 = - c3 / ( 2 * c1 ) ; lengths [ 0 ] = std :: sqrt ( c1 * lambda ) , lengths [ 1 ] = std :: sqrt ( c0 * lambda ) ; center = planePt + alpha0 * A + beta0 * B ; axes [ 0 ] = A ; axes [ 1 ] = B ; return true ; }
Compute the ellipse created by the intersection of an infinite
 cylinder and a plane.
Compute the ellipse created by the intersection of an infinite cylinder and a plane .