mingzailao
/
mitsuba
mirror of https://github.com/Quartenia/mitsuba.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

added clipping for cylinders

metadata
Wenzel Jakob 2010-10-22 19:13:45 +02:00
parent f2387696d1
commit 58eb27730f
2 changed files with 163 additions and 50 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

150
src/shapes/cylinder.cpp
View File

@ -174,20 +174,15 @@ public:
return result;
}
inline Float sqr(Float f) const {
return f*f;
}
/**
* 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
*/
FINLINE bool intersectCylPlane(Point planePt, Normal planeNrml,
bool intersectCylPlane(Point planePt, Normal planeNrml,
Point cylPt, Vector cylD, Float radius, Point &center,
Vector &axis1, Vector &axis2) const {
Vector *axes, Float *lengths) const {
if (absDot(planeNrml, cylD) < Epsilon)
return false;
@ -198,17 +193,16 @@ public:
A /= length;
B = cross(planeNrml, A);
} else {
Log(EInfo, "interesting: %s %s", cylD.toString().c_str(),
planeNrml.toString().c_str());
coordinateSystem(planeNrml, A, B);
}
Vector delta = planePt - cylPt,
deltaProj = delta - cylD*dot(delta, cylD);
Float c0 = 1-sqr(dot(A, cylD));
Float c1 = 1-sqr(dot(B, 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;
@ -216,35 +210,135 @@ public:
Float lambda = (c2*c2/(4*c0) + c3*c3/(4*c1) - c4)/(c0*c1);
Float alpha0 = -c2/(2*c0),
beta0 = -c3/(2*c1),
L_alpha = std::sqrt(c1*lambda),
L_beta = std::sqrt(c0*lambda);
beta0 = -c3/(2*c1);
lengths[0] = std::sqrt(c1*lambda),
lengths[1] = std::sqrt(c0*lambda);
center = planePt + alpha0 * A + beta0 * B;
axis1 = L_alpha * A;
axis2 = L_beta * B;
axes[0] = A;
axes[1] = B;
return true;
}
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,
ellipseCenter, ellipseAxes, ellipseLengths)) {
/* 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 */
return aabb;
}
/* Intersect the ellipse against the sides of the AABB face */
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);
/* Find the componentwise maxima of the ellipse */
for (int i=0; i<2; ++i) {
int j = (i==0) ? axis1 : axis2;
Float alpha = ellipseAxes[0][j];
Float beta = ellipseAxes[1][j];
Float ratio = beta/alpha, tmp = std::sqrt(1+ratio*ratio);
Float cosTheta = 1/tmp, sinTheta = ratio/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 getClippedAABB(const AABB &box) const {
/* Compute a base bounding box */
AABB base(getAABB());
base.clip(box);
Point p0 = m_objectToWorld(Point(0, 0, 0));
Point p1 = m_objectToWorld(Point(0, 0, m_length));
Vector cylD(normalize(p1-p0));
Point center;
Vector axis1, axis2;
intersectCylPlane(box.min, Normal(1, 0, 0), p0, cylD, m_radius, center, axis1, axis2);
// cout << center.toString() << endl;
Point cylPt = m_objectToWorld(Point(0, 0, 0));
Vector cylD(m_objectToWorld(Vector(0, 0, 1)));
/* Now forget about the cylinder ends and
intersect an infinite cylinder with each AABB face */
return base;
AABB clippedAABB;
clippedAABB.expandBy(intersectCylFace(0,
Point(base.min.x, base.min.y, base.min.z),
Point(base.min.x, base.max.y, base.max.z),
cylPt, cylD));
clippedAABB.expandBy(intersectCylFace(0,
Point(base.max.x, base.min.y, base.min.z),
Point(base.max.x, base.max.y, base.max.z),
cylPt, cylD));
clippedAABB.expandBy(intersectCylFace(1,
Point(base.min.x, base.min.y, base.min.z),
Point(base.max.x, base.min.y, base.max.z),
cylPt, cylD));
clippedAABB.expandBy(intersectCylFace(1,
Point(base.min.x, base.max.y, base.min.z),
Point(base.max.x, base.max.y, base.max.z),
cylPt, cylD));
clippedAABB.expandBy(intersectCylFace(2,
Point(base.min.x, base.min.y, base.min.z),
Point(base.max.x, base.max.y, base.min.z),
cylPt, cylD));
clippedAABB.expandBy(intersectCylFace(2,
Point(base.min.x, base.min.y, base.max.z),
Point(base.max.x, base.max.y, base.max.z),
cylPt, cylD));
clippedAABB.clip(box);
return clippedAABB;
}
#if 0

63
src/utils/cylclip.cpp
View File

@ -27,31 +27,35 @@ public:
m_viewTransform = Transform::lookAt(Point(10*std::sin(m_angle), 0, std::cos(m_angle)*10),
Point(0, 0, 0), Vector(0, 1, 0));
m_lineParams = Point2(M_PI/2, 0.28f);
m_cylPos = Point(0.0f);
m_device->setFSAA(4);
m_red[0] = 1.0f;
m_blue[2] = 1.0f;
m_showEllipses = false;
m_showRectangles = false;
m_showClippedAABB = false;
m_radius = .2f;
}
void mouseDragged(const DeviceEvent &event) {
if (event.getMouseButton() == 1) {
if (event.getMouseButton() == Device::ELeftButton) {
m_angle += event.getMouseRelative().x / 100.0f;
m_viewTransform = Transform::lookAt(Point(10*std::sin(m_angle), 0, std::cos(m_angle)*10),
Point(0, 0, 0), Vector(0, 1, 0));
} else {
} else if (event.getMouseButton() == Device::ERightButton) {
m_lineParams += Vector2(
event.getMouseRelative().x / 500.0f,
event.getMouseRelative().y / 500.0f
);
} else if (event.getMouseButton() == Device::EMiddleButton) {
m_cylPos += Vector3(
event.getMouseRelative().x / 300.0f,
0,
event.getMouseRelative().y / 300.0f
);
}
}
inline Float sqr(Float f) const {
return f*f;
}
/**
* Compute the ellipse created by the intersection of an infinite
* cylinder and a plane. Returns false in the degenerate case.
@ -77,8 +81,10 @@ public:
Vector delta = planePt - cylPt,
deltaProj = delta - cylD*dot(delta, cylD);
Float c0 = 1-sqr(dot(A, cylD));
Float c1 = 1-sqr(dot(B, 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;
@ -111,9 +117,13 @@ public:
Float ellipseLengths[2];
AABB aabb;
if (!intersectCylPlane(min, planeNrml, cylPt, cylD, .2,
ellipseCenter, ellipseAxes, ellipseLengths))
return aabb; // return an invalid AABB
if (!intersectCylPlane(min, planeNrml, cylPt, cylD, m_radius,
ellipseCenter, ellipseAxes, ellipseLengths)) {
/* 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 */
return aabb;
}
if (m_showEllipses) {
m_renderer->setColor(m_gray);
@ -200,6 +210,12 @@ public:
case 'a':
m_showClippedAABB = !m_showClippedAABB;
break;
case '[':
m_radius *= 1.1;
break;
case ']':
m_radius /= 1.1;
break;
}
}
@ -218,44 +234,46 @@ public:
m_renderer->setColor(m_gray);
Vector cylD(sphericalDirection(m_lineParams.x, m_lineParams.y));
Point cylP(0, 0, 0);
m_renderer->drawLine(cylP-cylD*1e4, cylP+cylD*1e4);
m_renderer->drawLine(m_cylPos-cylD*1e4, m_cylPos+cylD*1e4);
AABB clippedAABB;
clippedAABB.expandBy(intersectCylFace(0,
Point(aabb.min.x, aabb.min.y, aabb.min.z),
Point(aabb.min.x, aabb.max.y, aabb.max.z),
cylP, cylD));
m_cylPos, cylD));
clippedAABB.expandBy(intersectCylFace(0,
Point(aabb.max.x, aabb.min.y, aabb.min.z),
Point(aabb.max.x, aabb.max.y, aabb.max.z),
cylP, cylD));
m_cylPos, cylD));
clippedAABB.expandBy(intersectCylFace(1,
Point(aabb.min.x, aabb.min.y, aabb.min.z),
Point(aabb.max.x, aabb.min.y, aabb.max.z),
cylP, cylD));
m_cylPos, cylD));
clippedAABB.expandBy(intersectCylFace(1,
Point(aabb.min.x, aabb.max.y, aabb.min.z),
Point(aabb.max.x, aabb.max.y, aabb.max.z),
cylP, cylD));
m_cylPos, cylD));
clippedAABB.expandBy(intersectCylFace(2,
Point(aabb.min.x, aabb.min.y, aabb.min.z),
Point(aabb.max.x, aabb.max.y, aabb.min.z),
cylP, cylD));
m_cylPos, cylD));
clippedAABB.expandBy(intersectCylFace(2,
Point(aabb.min.x, aabb.min.y, aabb.max.z),
Point(aabb.max.x, aabb.max.y, aabb.max.z),
cylP, cylD));
m_cylPos, cylD));
m_renderer->setColor(m_gray);
if (m_showClippedAABB)
m_renderer->drawAABB(clippedAABB);
if (m_showClippedAABB) {
if (clippedAABB.isValid())
m_renderer->drawAABB(clippedAABB);
}
m_renderer->setDepthTest(false);
drawHUD(formatString("Cylinder clipping test. LMB-dragging moves the camera, RMB-dragging rotates the cylinder\n"
@ -273,7 +291,8 @@ private:
Transform m_projTransform, m_viewTransform;
Spectrum m_red, m_blue, m_gray;
Point2 m_lineParams;
Float m_angle;
Point m_cylPos;
Float m_angle, m_radius;
bool m_showEllipses;
bool m_showRectangles;
bool m_showClippedAABB;