mitsuba/src/libcore/triangle.cpp

/*
    This file is part of Mitsuba, a physically based rendering system.

    Copyright (c) 2007-2011 by Wenzel Jakob and others.

    Mitsuba is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License Version 3
    as published by the Free Software Foundation.

    Mitsuba is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#include <mitsuba/core/triangle.h>

MTS_NAMESPACE_BEGIN

Point Triangle::sample(const Point *positions, const Normal *normals,
		Normal &normal, const Point2 &sample) const {
	const Point &p0 = positions[idx[0]];
	const Point &p1 = positions[idx[1]];
	const Point &p2 = positions[idx[2]];

	Point2 bary = squareToTriangle(sample);
	Vector sideA = p1 - p0, sideB = p2 - p0;
	Point p = p0 + (sideA * bary.x) + (sideB * bary.y);	
	
	if (normals) {
		const Normal &n0 = normals[idx[0]];
		const Normal &n1 = normals[idx[1]];
		const Normal &n2 = normals[idx[2]];

		normal = Normal(normalize(
			n0 * (1.0f - bary.x - bary.y) +
			n1 * bary.x + n2 * bary.y
		));
	} else {
		normal = Normal(normalize(cross(sideA, sideB)));
	}

	return p;
}

Float Triangle::surfaceArea(const Point *positions) const {
	const Point &p0 = positions[idx[0]];
	const Point &p1 = positions[idx[1]];
	const Point &p2 = positions[idx[2]];
	Vector sideA = p1 - p0, sideB = p2 - p0;
	return 0.5f * cross(sideA, sideB).length();
}

#define MAX_VERTS 10

static int sutherlandHodgman(Point3d *input, int inCount, Point3d *output, int axis, 
		double splitPos, bool isMinimum) {
	if (inCount < 3)
		return 0;

	Point3d cur       = input[0];
	double sign       = isMinimum ? 1.0f : -1.0f;
	double distance   = sign * (cur[axis] - splitPos);
	bool  curIsInside = (distance >= 0);
	int   outCount    = 0;

	for (int i=0; i<inCount; ++i) {
		int nextIdx = i+1;
		if (nextIdx == inCount)
			nextIdx = 0;
		Point3d next = input[nextIdx];
		distance = sign * (next[axis] - splitPos);
		bool nextIsInside = (distance >= 0);

		if (curIsInside && nextIsInside) {
			/* Both this and the next vertex are inside, add to the list */
			SAssertEx(outCount + 1 < MAX_VERTS, "Overflow in sutherlandHodgman()!");
			output[outCount++] = next;
		} else if (curIsInside && !nextIsInside) {
			/* Going outside -- add the intersection */
			double t = (splitPos - cur[axis]) / (next[axis] - cur[axis]);
			SAssertEx(outCount + 1 < MAX_VERTS, "Overflow in sutherlandHodgman()!");
			Point3d p = cur + (next - cur) * t;
			p[axis] = splitPos; // Avoid roundoff errors
			output[outCount++] = p;
		} else if (!curIsInside && nextIsInside) {
			/* Coming back inside -- add the intersection + next vertex */
			double t = (splitPos - cur[axis]) / (next[axis] - cur[axis]);
			SAssertEx(outCount + 2 < MAX_VERTS, "Overflow in sutherlandHodgman()!");
			Point3d p = cur + (next - cur) * t;
			p[axis] = splitPos; // Avoid roundoff errors
			output[outCount++] = p;
			output[outCount++] = next;
		} else {
			/* Entirely outside - do not add anything */
		}
		cur = next;
		curIsInside = nextIsInside;
	}
	return outCount;
}

AABB Triangle::getClippedAABB(const Point *positions, const AABB &aabb) const {
	/* Reserve room for some additional vertices */
	Point3d vertices1[MAX_VERTS], vertices2[MAX_VERTS];
	int nVertices = 3;

	/* The kd-tree code will frequently call this function with
	   almost-collapsed AABBs. It's extremely important not to introduce
	   errors in such cases, otherwise the resulting tree will incorrectly
	   remove triangles from the associated nodes. Hence, do the
	   following computation in double precision! */
	for (int i=0; i<3; ++i) 
		vertices1[i] = Point3d(positions[idx[i]]);

	for (int axis=0; axis<3; ++axis) {
		nVertices = sutherlandHodgman(vertices1, nVertices, vertices2, axis, aabb.min[axis], true);
		nVertices = sutherlandHodgman(vertices2, nVertices, vertices1, axis, aabb.max[axis], false);
	}

	AABB result;
	for (int i=0; i<nVertices; ++i) {
		for (int j=0; j<3; ++j) {
			/* Now this is really paranoid! */
			double pos_d = vertices1[i][j];
			float pos_f = (float) pos_d;
			Float pos_roundedDown, pos_roundedUp;

			if (pos_f < pos_d) {
				/* Float value is too small */
				pos_roundedDown = pos_f;
				pos_roundedUp = nextafterf(pos_f, 
					std::numeric_limits<float>::infinity());
			} else if (pos_f > pos_d) {
				/* Float value is too large */
				pos_roundedUp = pos_f;
				pos_roundedDown = nextafterf(pos_f, 
					-std::numeric_limits<float>::infinity());
			} else {
				/* Double value is exactly representable */
				pos_roundedDown = pos_roundedUp = pos_f;
			}

			result.min[j] = std::min(result.min[j], pos_roundedDown);
			result.max[j] = std::max(result.max[j], pos_roundedUp);
		}
	}
	result.clip(aabb);

	return result;
}

MTS_NAMESPACE_END