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