radiant/brush/winding.cpp

/*
 Copyright (C) 1999-2006 Id Software, Inc. and contributors.
 For a list of contributors, see the accompanying CONTRIBUTORS file.

 This file is part of GtkRadiant.

 GtkRadiant is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 2 of the License, or
 (at your option) any later version.

 GtkRadiant 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 GtkRadiant; if not, write to the Free Software
 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

#include "winding.h"

#include <algorithm>

#include "math/line.h"

inline double plane3_distance_to_point (const Plane3& plane, const DoubleVector3& point)
{
	return point.dot(plane.normal()) - plane.dist();
}

inline double plane3_distance_to_point (const Plane3& plane, const Vector3& point)
{
	return point.dot(plane.normal()) - plane.dist();
}

/// \brief Returns the point at which \p line intersects \p plane, or an undefined value if there is no intersection.
inline DoubleVector3 line_intersect_plane (const DoubleLine& line, const Plane3& plane)
{
	return line.origin + line.direction * (-plane3_distance_to_point(plane, line.origin) / line.direction.dot(
			plane.normal()));
}

inline bool float_is_largest_absolute (double axis, double other)
{
	return fabs(axis) > fabs(other);
}

/// \brief Returns the index of the component of \p v that has the largest absolute value.
inline int vector3_largest_absolute_component_index (const DoubleVector3& v)
{
	return (float_is_largest_absolute(v[1], v[0])) ? (float_is_largest_absolute(v[1], v[2])) ? 1 : 2
			: (float_is_largest_absolute(v[0], v[2])) ? 0 : 2;
}

/// \brief Returns the infinite line that is the intersection of \p plane and \p other.
inline DoubleLine plane3_intersect_plane3 (const Plane3& plane, const Plane3& other)
{
	DoubleLine line;
	line.direction = plane.normal().crossProduct(other.normal());
	switch (vector3_largest_absolute_component_index(line.direction)) {
	case 0:
		line.origin.x() = 0;
		line.origin.y() = (-other.dist() * plane.normal().z() - -plane.dist() * other.normal().z())
				/ line.direction.x();
		line.origin.z() = (-plane.dist() * other.normal().y() - -other.dist() * plane.normal().y())
				/ line.direction.x();
		break;
	case 1:
		line.origin.x() = (-plane.dist() * other.normal().z() - -other.dist() * plane.normal().z())
				/ line.direction.y();
		line.origin.y() = 0;
		line.origin.z() = (-other.dist() * plane.normal().x() - -plane.dist() * other.normal().x())
				/ line.direction.y();
		break;
	case 2:
		line.origin.x() = (-other.dist() * plane.normal().y() - -plane.dist() * other.normal().y())
				/ line.direction.z();
		line.origin.y() = (-plane.dist() * other.normal().x() - -other.dist() * plane.normal().x())
				/ line.direction.z();
		line.origin.z() = 0;
		break;
	default:
		break;
	}

	return line;
}

/// \brief Keep the value of \p infinity as small as possible to improve precision in Winding_Clip.
void Winding_createInfinite (FixedWinding& winding, const Plane3& plane, double infinity)
{
	double max = -infinity;
	int x = -1;
	for (int i = 0; i < 3; i++) {
		const double d = fabs(plane.normal()[i]);
		if (d > max) {
			x = i;
			max = d;
		}
	}
	if (x == -1) {
		globalErrorStream() << "invalid plane\n";
		return;
	}

	DoubleVector3 vup = g_vector3_identity;
	switch (x) {
	case 0:
	case 1:
		vup[2] = 1;
		break;
	case 2:
		vup[0] = 1;
		break;
	}

	vup += plane.normal() * (-vup.dot(plane.normal()));
	vector3_normalise(vup);

	DoubleVector3 org = plane.normal() * plane.dist();

	DoubleVector3 vright = vup.crossProduct(plane.normal());

	vup *= infinity;
	vright *= infinity;

	// project a really big  axis aligned box onto the plane

	DoubleLine r1, r2, r3, r4;
	r1.origin = (org - vright) + vup;
	r1.direction = vright.getNormalised();
	winding.push_back(FixedWindingVertex(r1.origin, r1, c_brush_maxFaces));
	r2.origin = org + vright + vup;
	r2.direction = (-vup).getNormalised();
	winding.push_back(FixedWindingVertex(r2.origin, r2, c_brush_maxFaces));
	r3.origin = (org + vright) - vup;
	r3.direction = (-vright).getNormalised();
	winding.push_back(FixedWindingVertex(r3.origin, r3, c_brush_maxFaces));
	r4.origin = (org - vright) - vup;
	r4.direction = vup.getNormalised();
	winding.push_back(FixedWindingVertex(r4.origin, r4, c_brush_maxFaces));
}

inline PlaneClassification Winding_ClassifyDistance (const double distance, const double epsilon)
{
	if (distance > epsilon) {
		return ePlaneFront;
	}
	if (distance < -epsilon) {
		return ePlaneBack;
	}
	return ePlaneOn;
}

/// \brief Returns true if
/// !flipped && winding is completely BACK or ON
/// or flipped && winding is completely FRONT or ON
bool Winding_TestPlane (const Winding& winding, const Plane3& plane, bool flipped)
{
	const int test = (flipped) ? ePlaneBack : ePlaneFront;
	for (Winding::const_iterator i = winding.begin(); i != winding.end(); ++i) {
		if (test == Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON)) {
			return false;
		}
	}
	return true;
}

/// \brief Returns true if any point in \p w1 is in front of plane2, or any point in \p w2 is in front of plane1
bool Winding_PlanesConcave (const Winding& w1, const Winding& w2, const Plane3& plane1, const Plane3& plane2)
{
	return !Winding_TestPlane(w1, plane2, false) || !Winding_TestPlane(w2, plane1, false);
}

BrushSplitType Winding_ClassifyPlane (const Winding& winding, const Plane3& plane)
{
	BrushSplitType split;
	for (Winding::const_iterator i = winding.begin(); i != winding.end(); ++i) {
		++split.counts[Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON)];
	}
	return split;
}

const double DEBUG_EPSILON_SQUARED = ON_EPSILON * ON_EPSILON;

#define WINDING_DEBUG 0

/// \brief Clip \p winding which lies on \p plane by \p clipPlane, resulting in \p clipped.
/// If \p winding is completely in front of the plane, \p clipped will be identical to \p winding.
/// If \p winding is completely in back of the plane, \p clipped will be empty.
/// If \p winding intersects the plane, the edge of \p clipped which lies on \p clipPlane will store the value of \p adjacent.
void Winding_Clip (const FixedWinding& winding, const Plane3& plane, const Plane3& clipPlane, std::size_t adjacent,
		FixedWinding& clipped)
{
	PlaneClassification classification = Winding_ClassifyDistance(plane3_distance_to_point(clipPlane,
			winding.back().vertex), ON_EPSILON);
	PlaneClassification nextClassification;
	// for each edge
	for (std::size_t next = 0, i = winding.size() - 1; next != winding.size(); i = next, ++next, classification
			= nextClassification) {
		nextClassification = Winding_ClassifyDistance(plane3_distance_to_point(clipPlane, winding[next].vertex),
				ON_EPSILON);
		const FixedWindingVertex& vertex = winding[i];

		// if first vertex of edge is ON
		if (classification == ePlaneOn) {
			// append first vertex to output winding
			if (nextClassification == ePlaneBack) {
				// this edge lies on the clip plane
				clipped.push_back(
						FixedWindingVertex(vertex.vertex, plane3_intersect_plane3(plane, clipPlane), adjacent));
			} else {
				clipped.push_back(vertex);
			}
			continue;
		}

		// if first vertex of edge is FRONT
		if (classification == ePlaneFront) {
			// add first vertex to output winding
			clipped.push_back(vertex);
		}
		// if second vertex of edge is ON
		if (nextClassification == ePlaneOn) {
			continue;
		}
		// else if second vertex of edge is same as first
		else if (nextClassification == classification) {
			continue;
		}
		// else if first vertex of edge is FRONT and there are only two edges
		else if (classification == ePlaneFront && winding.size() == 2) {
			continue;
		}
		// else first vertex is FRONT and second is BACK or vice versa
		else {
			// append intersection point of line and plane to output winding
			DoubleVector3 mid(line_intersect_plane(vertex.edge, clipPlane));

			if (classification == ePlaneFront) {
				// this edge lies on the clip plane
				clipped.push_back(FixedWindingVertex(mid, plane3_intersect_plane3(plane, clipPlane), adjacent));
			} else {
				clipped.push_back(FixedWindingVertex(mid, vertex.edge, vertex.adjacent));
			}
		}
	}
}

std::size_t Winding_FindAdjacent (const Winding& winding, std::size_t face)
{
	for (std::size_t i = 0; i < winding.size(); ++i) {
		ASSERT_MESSAGE(winding[i].adjacent != c_brush_maxFaces, "edge connectivity data is invalid");
		if (winding[i].adjacent == face) {
			return i;
		}
	}
	return c_brush_maxFaces;
}

std::size_t Winding_Opposite (const Winding& winding, const std::size_t& index, const std::size_t& other)
{
	ASSERT_MESSAGE(index < winding.size() && other < winding.size(), "Winding_Opposite: index out of range");

	double dist_best = 0;
	std::size_t index_best = c_brush_maxFaces;

	Ray edge(ray_for_points(winding[index].vertex, winding[other].vertex));

	for (std::size_t i = 0; i < winding.size(); ++i) {
		if (i == index || i == other) {
			continue;
		}

		double dist_squared = ray_squared_distance_to_point(edge, winding[i].vertex);

		if (dist_squared > dist_best) {
			dist_best = dist_squared;
			index_best = i;
		}
	}
	return index_best;
}

std::size_t Winding_Opposite (const Winding& winding, const std::size_t index)
{
	return Winding_Opposite(winding, index, winding.next(index));
}

/// \brief Calculate the \p centroid of the polygon defined by \p winding which lies on plane \p plane.
void Winding_Centroid (const Winding& winding, const Plane3& plane, Vector3& centroid)
{
	double area2 = 0, x_sum = 0, y_sum = 0;
	const ProjectionAxis axis = projectionaxis_for_normal(plane.normal());
	const indexremap_t remap = indexremap_for_projectionaxis(axis);
	for (std::size_t i = winding.size() - 1, j = 0; j < winding.size(); i = j, ++j) {
		const double ai = winding[i].vertex[remap.x] * winding[j].vertex[remap.y] - winding[j].vertex[remap.x]
				* winding[i].vertex[remap.y];
		area2 += ai;
		x_sum += (winding[j].vertex[remap.x] + winding[i].vertex[remap.x]) * ai;
		y_sum += (winding[j].vertex[remap.y] + winding[i].vertex[remap.y]) * ai;
	}

	centroid[remap.x] = static_cast<float> (x_sum / (3 * area2));
	centroid[remap.y] = static_cast<float> (y_sum / (3 * area2));
	{
		Ray ray(Vector3(0, 0, 0), Vector3(0, 0, 0));
		ray.origin[remap.x] = centroid[remap.x];
		ray.origin[remap.y] = centroid[remap.y];
		ray.direction[remap.z] = 1;
		centroid[remap.z] = static_cast<float> (ray_distance_to_plane(ray, plane));
	}
}