mitsuba-af602c6fd98a/src/shapes/hair.cpp
/**
	 * \brief Intersect an infinite cylinder with an
	 * AABB face and bound the resulting clipped ellipse
	 */
\brief Intersect an infinite cylinder with an
AABB face and bound the resulting clipped ellipse
\ brief Intersect an infinite cylinder with an AABB face and bound the resulting clipped ellipse
AABB intersectCylFace(int axis,
			const Point &min, const Point &max,
			const Point &cylPt, const Vector &cylD) const {
		int axis1 = (axis + 1) % 3;
		int axis2 = (axis + 2) % 3;
		Normal planeNrml(0.0f);
		planeNrml[axis] = 1;
		Point ellipseCenter;
		Vector ellipseAxes[2];
		Float ellipseLengths[2];
		AABB aabb;
		if (!intersectCylPlane(min, planeNrml, cylPt, cylD, m_radius * (1 + Epsilon),
			ellipseCenter, ellipseAxes, ellipseLengths)) {
			return aabb;
		}
		for (int i=0; i<4; ++i) {
			Point p1, p2;
			p1[axis] = p2[axis] = min[axis];
			p1[axis1] = ((i+1) & 2) ? min[axis1] : max[axis1];
			p1[axis2] = ((i+0) & 2) ? min[axis2] : max[axis2];
			p2[axis1] = ((i+2) & 2) ? min[axis1] : max[axis1];
			p2[axis2] = ((i+1) & 2) ? min[axis2] : max[axis2];
			Point2 p1l(
				dot(p1 - ellipseCenter, ellipseAxes[0]) / ellipseLengths[0],
				dot(p1 - ellipseCenter, ellipseAxes[1]) / ellipseLengths[1]);
			Point2 p2l(
				dot(p2 - ellipseCenter, ellipseAxes[0]) / ellipseLengths[0],
				dot(p2 - ellipseCenter, ellipseAxes[1]) / ellipseLengths[1]);
			Vector2 rel = p2l-p1l;
			Float A = dot(rel, rel);
			Float B = 2*dot(Vector2(p1l), rel);
			Float C = dot(Vector2(p1l), Vector2(p1l))-1;
			Float x0, x1;
			if (solveQuadratic(A, B, C, x0, x1)) {
				if (x0 >= 0 && x0 <= 1)
					aabb.expandBy(p1+(p2-p1)*x0);
				if (x1 >= 0 && x1 <= 1)
					aabb.expandBy(p1+(p2-p1)*x1);
			}
		}
		ellipseAxes[0] *= ellipseLengths[0];
		ellipseAxes[1] *= ellipseLengths[1];
		AABB faceBounds(min, max);
		for (int i=0; i<2; ++i) {
			int j = (i==0) ? axis1 : axis2;
			Float alpha = ellipseAxes[0][j];
			Float beta = ellipseAxes[1][j];
			Float tmp = 1 / std::sqrt(alpha*alpha + beta*beta);
			Float cosTheta = alpha * tmp, sinTheta = beta*tmp;
			Point p1 = ellipseCenter + cosTheta*ellipseAxes[0] + sinTheta*ellipseAxes[1];
			Point p2 = ellipseCenter - cosTheta*ellipseAxes[0] - sinTheta*ellipseAxes[1];
			if (faceBounds.contains(p1))
				aabb.expandBy(p1);
			if (faceBounds.contains(p2))
				aabb.expandBy(p2);
		}
		return aabb;
	}
AABB intersectCylFace ( int axis , const Point & min , const Point & max , const Point & cylPt , const Vector & cylD ) const { int axis1 = ( axis + 1 ) % 3 ; int axis2 = ( axis + 2 ) % 3 ; Normal planeNrml ( 0.0f ) ; planeNrml [ axis ] = 1 ; Point ellipseCenter ; Vector ellipseAxes [ 2 ] ; Float ellipseLengths [ 2 ] ; AABB aabb ; if ( ! intersectCylPlane ( min , planeNrml , cylPt , cylD , m_radius * ( 1 + Epsilon ) , ellipseCenter , ellipseAxes , ellipseLengths ) ) { return aabb ; } for ( int i = 0 ; i < 4 ; ++ i ) { Point p1 , p2 ; p1 [ axis ] = p2 [ axis ] = min [ axis ] ; p1 [ axis1 ] = ( ( i + 1 ) & 2 ) ? min [ axis1 ] : max [ axis1 ] ; p1 [ axis2 ] = ( ( i + 0 ) & 2 ) ? min [ axis2 ] : max [ axis2 ] ; p2 [ axis1 ] = ( ( i + 2 ) & 2 ) ? min [ axis1 ] : max [ axis1 ] ; p2 [ axis2 ] = ( ( i + 1 ) & 2 ) ? min [ axis2 ] : max [ axis2 ] ; Point2 p1l ( dot ( p1 - ellipseCenter , ellipseAxes [ 0 ] ) / ellipseLengths [ 0 ] , dot ( p1 - ellipseCenter , ellipseAxes [ 1 ] ) / ellipseLengths [ 1 ] ) ; Point2 p2l ( dot ( p2 - ellipseCenter , ellipseAxes [ 0 ] ) / ellipseLengths [ 0 ] , dot ( p2 - ellipseCenter , ellipseAxes [ 1 ] ) / ellipseLengths [ 1 ] ) ; Vector2 rel = p2l - p1l ; Float A = dot ( rel , rel ) ; Float B = 2 * dot ( Vector2 ( p1l ) , rel ) ; Float C = dot ( Vector2 ( p1l ) , Vector2 ( p1l ) ) - 1 ; Float x0 , x1 ; if ( solveQuadratic ( A , B , C , x0 , x1 ) ) { if ( x0 >= 0 && x0 <= 1 ) aabb . expandBy ( p1 + ( p2 - p1 ) * x0 ) ; if ( x1 >= 0 && x1 <= 1 ) aabb . expandBy ( p1 + ( p2 - p1 ) * x1 ) ; } } ellipseAxes [ 0 ] *= ellipseLengths [ 0 ] ; ellipseAxes [ 1 ] *= ellipseLengths [ 1 ] ; AABB faceBounds ( min , max ) ; for ( int i = 0 ; i < 2 ; ++ i ) { int j = ( i == 0 ) ? axis1 : axis2 ; Float alpha = ellipseAxes [ 0 ] [ j ] ; Float beta = ellipseAxes [ 1 ] [ j ] ; Float tmp = 1 / std :: sqrt ( alpha * alpha + beta * beta ) ; Float cosTheta = alpha * tmp , sinTheta = beta * tmp ; Point p1 = ellipseCenter + cosTheta * ellipseAxes [ 0 ] + sinTheta * ellipseAxes [ 1 ] ; Point p2 = ellipseCenter - cosTheta * ellipseAxes [ 0 ] - sinTheta * ellipseAxes [ 1 ] ; if ( faceBounds . contains ( p1 ) ) aabb . expandBy ( p1 ) ; if ( faceBounds . contains ( p2 ) ) aabb . expandBy ( p2 ) ; } return aabb ; }
\brief Intersect an infinite cylinder with an
 AABB face and bound the resulting clipped ellipse
\ brief Intersect an infinite cylinder with an AABB face and bound the resulting clipped ellipse
/* Degenerate case -- return an invalid AABB. This is
			   not a problem, since one of the other faces will provide
			   enough information to arrive at a correct clipped AABB */ /* Intersect the ellipse against the sides of the AABB face */ /* Find the componentwise maxima of the ellipse */