radiant/winding.cpp

   1 /*
   2    Copyright (C) 1999-2006 Id Software, Inc. and contributors.
   3    For a list of contributors, see the accompanying CONTRIBUTORS file.
   4
   5    This file is part of GtkRadiant.
   6
   7    GtkRadiant is free software; you can redistribute it and/or modify
   8    it under the terms of the GNU General Public License as published by
   9    the Free Software Foundation; either version 2 of the License, or
  10    (at your option) any later version.
  11
  12    GtkRadiant is distributed in the hope that it will be useful,
  13    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15    GNU General Public License for more details.
  16
  17    You should have received a copy of the GNU General Public License
  18    along with GtkRadiant; if not, write to the Free Software
  19    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20  */
  21
  22 #include "winding.h"
  23
  24 #include <algorithm>
  25
  26 #include "math/line.h"
  27
  28
  29 inline double plane3_distance_to_point(const Plane3 &plane, const DoubleVector3 &point)
  30 {
  31     return vector3_dot(point, plane.normal()) - plane.dist();
  32 }
  33
  34 inline double plane3_distance_to_point(const Plane3 &plane, const Vector3 &point)
  35 {
  36     return vector3_dot(point, plane.normal()) - plane.dist();
  37 }
  38
  39 /// \brief Returns the point at which \p line intersects \p plane, or an undefined value if there is no intersection.
  40 inline DoubleVector3 line_intersect_plane(const DoubleLine &line, const Plane3 &plane)
  41 {
  42     return line.origin + vector3_scaled(
  43             line.direction,
  44             -plane3_distance_to_point(plane, line.origin)
  45             / vector3_dot(line.direction, plane.normal())
  46     );
  47 }
  48
  49 inline bool float_is_largest_absolute(double axis, double other)
  50 {
  51     return fabs(axis) > fabs(other);
  52 }
  53
  54 /// \brief Returns the index of the component of \p v that has the largest absolute value.
  55 inline int vector3_largest_absolute_component_index(const DoubleVector3 &v)
  56 {
  57     return (float_is_largest_absolute(v[1], v[0]))
  58            ? (float_is_largest_absolute(v[1], v[2]))
  59              ? 1
  60              : 2
  61            : (float_is_largest_absolute(v[0], v[2]))
  62              ? 0
  63              : 2;
  64 }
  65
  66 /// \brief Returns the infinite line that is the intersection of \p plane and \p other.
  67 inline DoubleLine plane3_intersect_plane3(const Plane3 &plane, const Plane3 &other)
  68 {
  69     DoubleLine line;
  70     line.direction = vector3_cross(plane.normal(), other.normal());
  71     switch (vector3_largest_absolute_component_index(line.direction)) {
  72         case 0:
  73             line.origin.x() = 0;
  74             line.origin.y() =
  75                     (-other.dist() * plane.normal().z() - -plane.dist() * other.normal().z()) / line.direction.x();
  76             line.origin.z() =
  77                     (-plane.dist() * other.normal().y() - -other.dist() * plane.normal().y()) / line.direction.x();
  78             break;
  79         case 1:
  80             line.origin.x() =
  81                     (-plane.dist() * other.normal().z() - -other.dist() * plane.normal().z()) / line.direction.y();
  82             line.origin.y() = 0;
  83             line.origin.z() =
  84                     (-other.dist() * plane.normal().x() - -plane.dist() * other.normal().x()) / line.direction.y();
  85             break;
  86         case 2:
  87             line.origin.x() =
  88                     (-other.dist() * plane.normal().y() - -plane.dist() * other.normal().y()) / line.direction.z();
  89             line.origin.y() =
  90                     (-plane.dist() * other.normal().x() - -other.dist() * plane.normal().x()) / line.direction.z();
  91             line.origin.z() = 0;
  92             break;
  93         default:
  94             break;
  95     }
  96
  97     return line;
  98 }
  99
 100
 101 /// \brief Keep the value of \p infinity as small as possible to improve precision in Winding_Clip.
 102 void Winding_createInfinite(FixedWinding &winding, const Plane3 &plane, double infinity)
 103 {
 104     double max = -infinity;
 105     int x = -1;
 106     for (int i = 0; i < 3; i++) {
 107         double d = fabs(plane.normal()[i]);
 108         if (d > max) {
 109             x = i;
 110             max = d;
 111         }
 112     }
 113     if (x == -1) {
 114         globalErrorStream() << "invalid plane\n";
 115         return;
 116     }
 117
 118     DoubleVector3 vup = g_vector3_identity;
 119     switch (x) {
 120         case 0:
 121         case 1:
 122             vup[2] = 1;
 123             break;
 124         case 2:
 125             vup[0] = 1;
 126             break;
 127     }
 128
 129
 130     vector3_add(vup, vector3_scaled(plane.normal(), -vector3_dot(vup, plane.normal())));
 131     vector3_normalise(vup);
 132
 133     DoubleVector3 org = vector3_scaled(plane.normal(), plane.dist());
 134
 135     DoubleVector3 vright = vector3_cross(vup, plane.normal());
 136
 137     vector3_scale(vup, infinity);
 138     vector3_scale(vright, infinity);
 139
 140     // project a really big  axis aligned box onto the plane
 141
 142     DoubleLine r1, r2, r3, r4;
 143     r1.origin = vector3_added(vector3_subtracted(org, vright), vup);
 144     r1.direction = vector3_normalised(vright);
 145     winding.push_back(FixedWindingVertex(r1.origin, r1, c_brush_maxFaces));
 146     r2.origin = vector3_added(vector3_added(org, vright), vup);
 147     r2.direction = vector3_normalised(vector3_negated(vup));
 148     winding.push_back(FixedWindingVertex(r2.origin, r2, c_brush_maxFaces));
 149     r3.origin = vector3_subtracted(vector3_added(org, vright), vup);
 150     r3.direction = vector3_normalised(vector3_negated(vright));
 151     winding.push_back(FixedWindingVertex(r3.origin, r3, c_brush_maxFaces));
 152     r4.origin = vector3_subtracted(vector3_subtracted(org, vright), vup);
 153     r4.direction = vector3_normalised(vup);
 154     winding.push_back(FixedWindingVertex(r4.origin, r4, c_brush_maxFaces));
 155 }
 156
 157
 158 inline PlaneClassification Winding_ClassifyDistance(const double distance, const double epsilon)
 159 {
 160     if (distance > epsilon) {
 161         return ePlaneFront;
 162     }
 163     if (distance < -epsilon) {
 164         return ePlaneBack;
 165     }
 166     return ePlaneOn;
 167 }
 168
 169 /// \brief Returns true if
 170 /// !flipped && winding is completely BACK or ON
 171 /// or flipped && winding is completely FRONT or ON
 172 bool Winding_TestPlane(const Winding &winding, const Plane3 &plane, bool flipped)
 173 {
 174     const int test = (flipped) ? ePlaneBack : ePlaneFront;
 175     for (Winding::const_iterator i = winding.begin(); i != winding.end(); ++i) {
 176         if (test == Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON)) {
 177             return false;
 178         }
 179     }
 180     return true;
 181 }
 182
 183 /// \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
 184 bool Winding_PlanesConcave(const Winding &w1, const Winding &w2, const Plane3 &plane1, const Plane3 &plane2)
 185 {
 186     return !Winding_TestPlane(w1, plane2, false) || !Winding_TestPlane(w2, plane1, false);
 187 }
 188
 189 brushsplit_t Winding_ClassifyPlane(const Winding &winding, const Plane3 &plane)
 190 {
 191     brushsplit_t split;
 192     for (Winding::const_iterator i = winding.begin(); i != winding.end(); ++i) {
 193         ++split.counts[Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON)];
 194     }
 195     return split;
 196 }
 197
 198
 199 #define DEBUG_EPSILON ON_EPSILON
 200 const double DEBUG_EPSILON_SQUARED = DEBUG_EPSILON * DEBUG_EPSILON;
 201
 202 #define WINDING_DEBUG 0
 203
 204 /// \brief Clip \p winding which lies on \p plane by \p clipPlane, resulting in \p clipped.
 205 /// If \p winding is completely in front of the plane, \p clipped will be identical to \p winding.
 206 /// If \p winding is completely in back of the plane, \p clipped will be empty.
 207 /// If \p winding intersects the plane, the edge of \p clipped which lies on \p clipPlane will store the value of \p adjacent.
 208 void Winding_Clip(const FixedWinding &winding, const Plane3 &plane, const Plane3 &clipPlane, std::size_t adjacent,
 209                   FixedWinding &clipped)
 210 {
 211     PlaneClassification classification = Winding_ClassifyDistance(
 212             plane3_distance_to_point(clipPlane, winding.back().vertex), ON_EPSILON);
 213     PlaneClassification nextClassification;
 214     // for each edge
 215     for (std::size_t next = 0, i = winding.size() - 1;
 216          next != winding.size(); i = next, ++next, classification = nextClassification) {
 217         nextClassification = Winding_ClassifyDistance(plane3_distance_to_point(clipPlane, winding[next].vertex),
 218                                                       ON_EPSILON);
 219         const FixedWindingVertex &vertex = winding[i];
 220
 221         // if first vertex of edge is ON
 222         if (classification == ePlaneOn) {
 223             // append first vertex to output winding
 224             if (nextClassification == ePlaneBack) {
 225                 // this edge lies on the clip plane
 226                 clipped.push_back(
 227                         FixedWindingVertex(vertex.vertex, plane3_intersect_plane3(plane, clipPlane), adjacent));
 228             } else {
 229                 clipped.push_back(vertex);
 230             }
 231             continue;
 232         }
 233
 234         // if first vertex of edge is FRONT
 235         if (classification == ePlaneFront) {
 236             // add first vertex to output winding
 237             clipped.push_back(vertex);
 238         }
 239         // if second vertex of edge is ON
 240         if (nextClassification == ePlaneOn) {
 241             continue;
 242         }
 243             // else if second vertex of edge is same as first
 244         else if (nextClassification == classification) {
 245             continue;
 246         }
 247             // else if first vertex of edge is FRONT and there are only two edges
 248         else if (classification == ePlaneFront && winding.size() == 2) {
 249             continue;
 250         }
 251             // else first vertex is FRONT and second is BACK or vice versa
 252         else {
 253             // append intersection point of line and plane to output winding
 254             DoubleVector3 mid(line_intersect_plane(vertex.edge, clipPlane));
 255
 256             if (classification == ePlaneFront) {
 257                 // this edge lies on the clip plane
 258                 clipped.push_back(FixedWindingVertex(mid, plane3_intersect_plane3(plane, clipPlane), adjacent));
 259             } else {
 260                 clipped.push_back(FixedWindingVertex(mid, vertex.edge, vertex.adjacent));
 261             }
 262         }
 263     }
 264 }
 265
 266 std::size_t Winding_FindAdjacent(const Winding &winding, std::size_t face)
 267 {
 268     for (std::size_t i = 0; i < winding.numpoints; ++i) {
 269         ASSERT_MESSAGE(winding[i].adjacent != c_brush_maxFaces, "edge connectivity data is invalid");
 270         if (winding[i].adjacent == face) {
 271             return i;
 272         }
 273     }
 274     return c_brush_maxFaces;
 275 }
 276
 277 std::size_t Winding_Opposite(const Winding &winding, const std::size_t index, const std::size_t other)
 278 {
 279     ASSERT_MESSAGE(index < winding.numpoints && other < winding.numpoints, "Winding_Opposite: index out of range");
 280
 281     double dist_best = 0;
 282     std::size_t index_best = c_brush_maxFaces;
 283
 284     Ray edge(ray_for_points(winding[index].vertex, winding[other].vertex));
 285
 286     for (std::size_t i = 0; i < winding.numpoints; ++i) {
 287         if (i == index || i == other) {
 288             continue;
 289         }
 290
 291         double dist_squared = ray_squared_distance_to_point(edge, winding[i].vertex);
 292
 293         if (dist_squared > dist_best) {
 294             dist_best = dist_squared;
 295             index_best = i;
 296         }
 297     }
 298     return index_best;
 299 }
 300
 301 std::size_t Winding_Opposite(const Winding &winding, const std::size_t index)
 302 {
 303     return Winding_Opposite(winding, index, Winding_next(winding, index));
 304 }
 305
 306 /// \brief Calculate the \p centroid of the polygon defined by \p winding which lies on plane \p plane.
 307 void Winding_Centroid(const Winding &winding, const Plane3 &plane, Vector3 &centroid)
 308 {
 309     double area2 = 0, x_sum = 0, y_sum = 0;
 310     const ProjectionAxis axis = projectionaxis_for_normal(plane.normal());
 311     const indexremap_t remap = indexremap_for_projectionaxis(axis);
 312     for (std::size_t i = winding.numpoints - 1, j = 0; j < winding.numpoints; i = j, ++j) {
 313         const double ai = winding[i].vertex[remap.x] * winding[j].vertex[remap.y] -
 314                           winding[j].vertex[remap.x] * winding[i].vertex[remap.y];
 315         area2 += ai;
 316         x_sum += (winding[j].vertex[remap.x] + winding[i].vertex[remap.x]) * ai;
 317         y_sum += (winding[j].vertex[remap.y] + winding[i].vertex[remap.y]) * ai;
 318     }
 319
 320     centroid[remap.x] = static_cast<float>( x_sum / (3 * area2));
 321     centroid[remap.y] = static_cast<float>( y_sum / (3 * area2));
 322     {
 323         Ray ray(Vector3(0, 0, 0), Vector3(0, 0, 0));
 324         ray.origin[remap.x] = centroid[remap.x];
 325         ray.origin[remap.y] = centroid[remap.y];
 326         ray.direction[remap.z] = 1;
 327         centroid[remap.z] = static_cast<float>( ray_distance_to_plane(ray, plane));
 328     }
 329 }