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   {
  73   case 0:
  74     line.origin.x() = 0;
  75     line.origin.y() = (-other.dist() * plane.normal().z() - -plane.dist() * other.normal().z()) / line.direction.x();
  76     line.origin.z() = (-plane.dist() * other.normal().y() - -other.dist() * plane.normal().y()) / line.direction.x();
  77     break;
  78   case 1:
  79     line.origin.x() = (-plane.dist() * other.normal().z() - -other.dist() * plane.normal().z()) / line.direction.y();
  80     line.origin.y() = 0;
  81     line.origin.z() = (-other.dist() * plane.normal().x() - -plane.dist() * other.normal().x()) / line.direction.y();
  82     break;
  83   case 2:
  84     line.origin.x() = (-other.dist() * plane.normal().y() - -plane.dist() * other.normal().y()) / line.direction.z();
  85     line.origin.y() = (-plane.dist() * other.normal().x() - -other.dist() * plane.normal().x()) / line.direction.z();
  86     line.origin.z() = 0;
  87     break;
  88   default:
  89     break;
  90   }
  91
  92   return line;
  93 }
  94
  95
  96 /// \brief Keep the value of \p infinity as small as possible to improve precision in Winding_Clip.
  97 void Winding_createInfinite(FixedWinding& winding, const Plane3& plane, double infinity)
  98 {
  99   double max = -infinity;
 100   int x = -1;
 101   for (int i=0 ; i<3; i++)
 102   {
 103     double d = fabs(plane.normal()[i]);
 104     if (d > max)
 105     {
 106       x = i;
 107       max = d;
 108     }
 109   }
 110   if(x == -1)
 111   {
 112     globalErrorStream() << "invalid plane\n";
 113     return;
 114   }
 115
 116   DoubleVector3 vup = g_vector3_identity;
 117   switch (x)
 118   {
 119   case 0:
 120   case 1:
 121     vup[2] = 1;
 122     break;
 123   case 2:
 124     vup[0] = 1;
 125     break;
 126   }
 127
 128
 129   vector3_add(vup, vector3_scaled(plane.normal(), -vector3_dot(vup, plane.normal())));
 130   vector3_normalise(vup);
 131
 132   DoubleVector3 org = vector3_scaled(plane.normal(), plane.dist());
 133
 134   DoubleVector3 vright = vector3_cross(vup, plane.normal());
 135
 136   vector3_scale(vup, infinity);
 137   vector3_scale(vright, infinity);
 138
 139   // project a really big  axis aligned box onto the plane
 140
 141   DoubleLine r1, r2, r3, r4;
 142   r1.origin = vector3_added(vector3_subtracted(org, vright), vup);
 143   r1.direction = vector3_normalised(vright);
 144   winding.push_back(FixedWindingVertex(r1.origin, r1, c_brush_maxFaces));
 145   r2.origin = vector3_added(vector3_added(org, vright), vup);
 146   r2.direction = vector3_normalised(vector3_negated(vup));
 147   winding.push_back(FixedWindingVertex(r2.origin, r2, c_brush_maxFaces));
 148   r3.origin = vector3_subtracted(vector3_added(org, vright), vup);
 149   r3.direction = vector3_normalised(vector3_negated(vright));
 150   winding.push_back(FixedWindingVertex(r3.origin, r3, c_brush_maxFaces));
 151   r4.origin = vector3_subtracted(vector3_subtracted(org, vright), vup);
 152   r4.direction = vector3_normalised(vup);
 153   winding.push_back(FixedWindingVertex(r4.origin, r4, c_brush_maxFaces));
 154 }
 155
 156
 157 inline PlaneClassification Winding_ClassifyDistance(const double distance, const double epsilon)
 158 {
 159   if(distance > epsilon)
 160   {
 161     return ePlaneFront;
 162   }
 163   if(distance < -epsilon)
 164   {
 165     return ePlaneBack;
 166   }
 167   return ePlaneOn;
 168 }
 169
 170 /// \brief Returns true if
 171 /// !flipped && winding is completely BACK or ON
 172 /// or flipped && winding is completely FRONT or ON
 173 bool Winding_TestPlane(const Winding& winding, const Plane3& plane, bool flipped)
 174 {
 175   const int test = (flipped) ? ePlaneBack : ePlaneFront;
 176   for(Winding::const_iterator i = winding.begin(); i != winding.end(); ++i)
 177   {
 178     if(test == Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON))
 179     {
 180       return false;
 181     }
 182   }
 183   return true;
 184 }
 185
 186 /// \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
 187 bool Winding_PlanesConcave(const Winding& w1, const Winding& w2, const Plane3& plane1, const Plane3& plane2)
 188 {
 189   return !Winding_TestPlane(w1, plane2, false) || !Winding_TestPlane(w2, plane1, false);
 190 }
 191
 192 brushsplit_t Winding_ClassifyPlane(const Winding& winding, const Plane3& plane)
 193 {
 194   brushsplit_t split;
 195   for(Winding::const_iterator i = winding.begin(); i != winding.end(); ++i)
 196   {
 197     ++split.counts[Winding_ClassifyDistance(plane3_distance_to_point(plane, (*i).vertex), ON_EPSILON)];
 198   }
 199   return split;
 200 }
 201
 202
 203 #define DEBUG_EPSILON ON_EPSILON
 204 const double DEBUG_EPSILON_SQUARED = DEBUG_EPSILON * DEBUG_EPSILON;
 205
 206 #define WINDING_DEBUG 0
 207
 208 /// \brief Clip \p winding which lies on \p plane by \p clipPlane, resulting in \p clipped.
 209 /// If \p winding is completely in front of the plane, \p clipped will be identical to \p winding.
 210 /// If \p winding is completely in back of the plane, \p clipped will be empty.
 211 /// If \p winding intersects the plane, the edge of \p clipped which lies on \p clipPlane will store the value of \p adjacent.
 212 void Winding_Clip(const FixedWinding& winding, const Plane3& plane, const Plane3& clipPlane, std::size_t adjacent, FixedWinding& clipped)
 213 {
 214   PlaneClassification classification = Winding_ClassifyDistance(plane3_distance_to_point(clipPlane, winding.back().vertex), ON_EPSILON);
 215   PlaneClassification nextClassification;
 216   // for each edge
 217   for(std::size_t next = 0, i = winding.size()-1; next != winding.size(); i = next, ++next, classification = nextClassification)
 218   {
 219     nextClassification = Winding_ClassifyDistance(plane3_distance_to_point(clipPlane, winding[next].vertex), ON_EPSILON);
 220     const FixedWindingVertex& vertex = winding[i];
 221
 222     // if first vertex of edge is ON
 223     if(classification == ePlaneOn)
 224     {
 225       // append first vertex to output winding
 226       if(nextClassification == ePlaneBack)
 227       {
 228         // this edge lies on the clip plane
 229         clipped.push_back(FixedWindingVertex(vertex.vertex, plane3_intersect_plane3(plane, clipPlane), adjacent));
 230       }
 231       else
 232       {
 233         clipped.push_back(vertex);
 234       }
 235       continue;
 236     }
 237
 238     // if first vertex of edge is FRONT
 239     if(classification == ePlaneFront)
 240     {
 241       // add first vertex to output winding
 242       clipped.push_back(vertex);
 243     }
 244     // if second vertex of edge is ON
 245     if(nextClassification == ePlaneOn)
 246     {
 247       continue;
 248     }
 249     // else if second vertex of edge is same as first
 250     else if(nextClassification == classification)
 251     {
 252       continue;
 253     }
 254     // else if first vertex of edge is FRONT and there are only two edges
 255     else if(classification == ePlaneFront && winding.size() == 2)
 256     {
 257       continue;
 258     }
 259     // else first vertex is FRONT and second is BACK or vice versa
 260     else
 261     {
 262       // append intersection point of line and plane to output winding
 263       DoubleVector3 mid(line_intersect_plane(vertex.edge, clipPlane));
 264
 265       if(classification == ePlaneFront)
 266       {
 267         // this edge lies on the clip plane
 268         clipped.push_back(FixedWindingVertex(mid, plane3_intersect_plane3(plane, clipPlane), adjacent));
 269       }
 270       else
 271       {
 272         clipped.push_back(FixedWindingVertex(mid, vertex.edge, vertex.adjacent));
 273       }
 274     }
 275   }
 276 }
 277
 278 std::size_t Winding_FindAdjacent(const Winding& winding, std::size_t face)
 279 {
 280   for(std::size_t i=0; i<winding.numpoints; ++i)
 281   {
 282     ASSERT_MESSAGE(winding[i].adjacent != c_brush_maxFaces, "edge connectivity data is invalid");
 283     if(winding[i].adjacent == face)
 284     {
 285       return i;
 286     }
 287   }
 288   return c_brush_maxFaces;
 289 }
 290
 291 std::size_t Winding_Opposite(const Winding& winding, const std::size_t index, const std::size_t other)
 292 {
 293   ASSERT_MESSAGE(index < winding.numpoints && other < winding.numpoints, "Winding_Opposite: index out of range");
 294
 295   double dist_best = 0;
 296   std::size_t index_best = c_brush_maxFaces;
 297
 298   Ray edge(ray_for_points(winding[index].vertex, winding[other].vertex));
 299
 300   for(std::size_t i=0; i<winding.numpoints; ++i)
 301   {
 302     if(i == index || i == other)
 303     {
 304       continue;
 305     }
 306
 307     double dist_squared = ray_squared_distance_to_point(edge, winding[i].vertex);
 308
 309     if(dist_squared > dist_best)
 310     {
 311       dist_best = dist_squared;
 312       index_best = i;
 313     }
 314   }
 315   return index_best;
 316 }
 317
 318 std::size_t Winding_Opposite(const Winding& winding, const std::size_t index)
 319 {
 320   return Winding_Opposite(winding, index, Winding_next(winding, index));
 321 }
 322
 323 /// \brief Calculate the \p centroid of the polygon defined by \p winding which lies on plane \p plane.
 324 void Winding_Centroid(const Winding& winding, const Plane3& plane, Vector3& centroid)
 325 {
 326   double area2 = 0, x_sum = 0, y_sum = 0;
 327   const ProjectionAxis axis = projectionaxis_for_normal(plane.normal());
 328   const indexremap_t remap = indexremap_for_projectionaxis(axis);
 329   for(std::size_t i = winding.numpoints-1, j = 0; j < winding.numpoints; i = j, ++j)
 330   {
 331     const double ai = winding[i].vertex[remap.x] * winding[j].vertex[remap.y] - winding[j].vertex[remap.x] * winding[i].vertex[remap.y];
 332     area2 += ai;
 333     x_sum += (winding[j].vertex[remap.x] + winding[i].vertex[remap.x]) * ai;
 334     y_sum += (winding[j].vertex[remap.y] + winding[i].vertex[remap.y]) * ai;
 335   }
 336
 337   centroid[remap.x] = static_cast<float>(x_sum / (3 * area2));
 338   centroid[remap.y] = static_cast<float>(y_sum / (3 * area2));
 339   {
 340     Ray ray(Vector3(0, 0, 0), Vector3(0, 0, 0));
 341     ray.origin[remap.x] = centroid[remap.x];
 342     ray.origin[remap.y] = centroid[remap.y];
 343     ray.direction[remap.z] = 1;
 344     centroid[remap.z] = static_cast<float>(ray_distance_to_plane(ray, plane));
 345   }
 346 }