QGIS API Documentation  3.8.0-Zanzibar (11aff65)
qgsraycastingutils_p.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsraycastingutils_h.cpp
3  --------------------------------------
4  Date : June 2018
5  Copyright : (C) 2018 by Martin Dobias
6  Email : wonder dot sk at gmail dot com
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsraycastingutils_p.h"
17 
18 #include "qgis.h"
19 #include "qgsaabb.h"
20 
21 #include <Qt3DRender/QCamera>
22 
24 
25 
26 namespace QgsRayCastingUtils
27 {
28 
29 // copied from qt3d/src/render/raycasting/qray3d_p.h + qray3d.cpp
30 // by KDAB, licensed under the terms of LGPL
31 
32 
33  Ray3D::Ray3D()
34  : m_direction( 0.0f, 0.0f, 1.0f )
35  {
36  }
37 
38  Ray3D::Ray3D( QVector3D origin, QVector3D direction, float distance )
39  : m_origin( origin )
40  , m_direction( direction )
41  , m_distance( distance )
42  {}
43 
44  QVector3D Ray3D::origin() const
45  {
46  return m_origin;
47  }
48 
49  void Ray3D::setOrigin( QVector3D value )
50  {
51  m_origin = value;
52  }
53 
54  QVector3D Ray3D::direction() const
55  {
56  return m_direction;
57  }
58 
59  void Ray3D::setDirection( QVector3D value )
60  {
61  if ( value.isNull() )
62  return;
63 
64  m_direction = value;
65  }
66 
67  float Ray3D::distance() const
68  {
69  return m_distance;
70  }
71 
72  void Ray3D::setDistance( float distance )
73  {
74  m_distance = distance;
75  }
76 
77  QVector3D Ray3D::point( float t ) const
78  {
79  return m_origin + t * m_direction;
80  }
81 
82  Ray3D &Ray3D::transform( const QMatrix4x4 &matrix )
83  {
84  m_origin = matrix * m_origin;
85  m_direction = matrix.mapVector( m_direction );
86 
87  return *this;
88  }
89 
90  Ray3D Ray3D::transformed( const QMatrix4x4 &matrix ) const
91  {
92  return Ray3D( matrix * m_origin, matrix.mapVector( m_direction ) );
93  }
94 
95  bool Ray3D::operator==( const Ray3D &other ) const
96  {
97  return m_origin == other.origin() && m_direction == other.direction();
98  }
99 
100  bool Ray3D::operator!=( const Ray3D &other ) const
101  {
102  return !( *this == other );
103  }
104 
105  bool Ray3D::contains( QVector3D point ) const
106  {
107  QVector3D ppVec( point - m_origin );
108  if ( ppVec.isNull() ) // point coincides with origin
109  return true;
110  const float dot = QVector3D::dotProduct( ppVec, m_direction );
111  if ( qFuzzyIsNull( dot ) )
112  return false;
113  return qFuzzyCompare( dot * dot, ppVec.lengthSquared() * m_direction.lengthSquared() );
114  }
115 
116  bool Ray3D::contains( const Ray3D &ray ) const
117  {
118  const float dot = QVector3D::dotProduct( m_direction, ray.direction() );
119  if ( !qFuzzyCompare( dot * dot, m_direction.lengthSquared() * ray.direction().lengthSquared() ) )
120  return false;
121  return contains( ray.origin() );
122  }
123 
124  float Ray3D::projectedDistance( QVector3D point ) const
125  {
126  Q_ASSERT( !m_direction.isNull() );
127 
128  return QVector3D::dotProduct( point - m_origin, m_direction ) /
129  m_direction.lengthSquared();
130  }
131 
132  QVector3D Ray3D::project( QVector3D vector ) const
133  {
134  QVector3D norm = m_direction.normalized();
135  return QVector3D::dotProduct( vector, norm ) * norm;
136  }
137 
138 
139  float Ray3D::distance( QVector3D point ) const
140  {
141  float t = projectedDistance( point );
142  return ( point - ( m_origin + t * m_direction ) ).length();
143  }
144 
145  QDebug operator<<( QDebug dbg, const Ray3D &ray )
146  {
147  QDebugStateSaver saver( dbg );
148  dbg.nospace() << "QRay3D(origin("
149  << ray.origin().x() << ", " << ray.origin().y() << ", "
150  << ray.origin().z() << ") - direction("
151  << ray.direction().x() << ", " << ray.direction().y() << ", "
152  << ray.direction().z() << "))";
153  return dbg;
154  }
155 
156 }
157 
158 
160 
161 
162 struct box
163 {
164  box( const QgsAABB &b )
165  {
166  min[0] = b.xMin; min[1] = b.yMin; min[2] = b.zMin;
167  max[0] = b.xMax; max[1] = b.yMax; max[2] = b.zMax;
168  }
169  double min[3];
170  double max[3];
171 };
172 
173 struct ray
174 {
175  ray( double xO, double yO, double zO, double xD, double yD, double zD )
176  {
177  origin[0] = xO; origin[1] = yO; origin[2] = zO;
178  dir[0] = xD; dir[1] = yD; dir[2] = zD;
179  dir_inv[0] = 1 / dir[0]; dir_inv[1] = 1 / dir[1]; dir_inv[2] = 1 / dir[2];
180  }
181  ray( const QgsRayCastingUtils::Ray3D &r )
182  {
183  // ignoring length...
184  origin[0] = r.origin().x(); origin[1] = r.origin().y(); origin[2] = r.origin().z();
185  dir[0] = r.direction().x(); dir[1] = r.direction().y(); dir[2] = r.direction().z();
186  dir_inv[0] = 1 / dir[0]; dir_inv[1] = 1 / dir[1]; dir_inv[2] = 1 / dir[2];
187  }
188  double origin[3];
189  double dir[3];
190  double dir_inv[3];
191 };
192 
193 // https://tavianator.com/fast-branchless-raybounding-box-intersections/
194 // https://tavianator.com/fast-branchless-raybounding-box-intersections-part-2-nans/
195 
196 bool intersection( const box &b, const ray &r )
197 {
198  double t1 = ( b.min[0] - r.origin[0] ) * r.dir_inv[0];
199  double t2 = ( b.max[0] - r.origin[0] ) * r.dir_inv[0];
200 
201  double tmin = std::min( t1, t2 );
202  double tmax = std::max( t1, t2 );
203 
204  for ( int i = 1; i < 3; ++i )
205  {
206  t1 = ( b.min[i] - r.origin[i] ) * r.dir_inv[i];
207  t2 = ( b.max[i] - r.origin[i] ) * r.dir_inv[i];
208 
209  tmin = std::max( tmin, std::min( std::min( t1, t2 ), tmax ) );
210  tmax = std::min( tmax, std::max( std::max( t1, t2 ), tmin ) );
211  }
212 
213  return tmax > std::max( tmin, 0.0 );
214 }
215 
216 
217 namespace QgsRayCastingUtils
218 {
219 
220  bool rayBoxIntersection( const Ray3D &r, const QgsAABB &aabb )
221  {
222  box b( aabb );
223 
224  // intersection() does not like yMin==yMax (excludes borders)
225  if ( b.min[0] == b.max[0] ) b.max[0] += 0.1;
226  if ( b.min[1] == b.max[1] ) b.max[1] += 0.1;
227  if ( b.min[2] == b.max[2] ) b.max[2] += 0.1;
228 
229  return intersection( b, ray( r ) );
230  }
231 
232  // copied from https://stackoverflow.com/questions/23975555/how-to-do-ray-plane-intersection
233  bool rayPlaneIntersection( const Ray3D &r, const Plane3D &plane, QVector3D &pt )
234  {
235  float denom = QVector3D::dotProduct( plane.normal, r.direction() );
236  if ( std::abs( denom ) > 0.0001f ) // your favorite epsilon
237  {
238  float t = QVector3D::dotProduct( plane.center - r.origin(), plane.normal ) / denom;
239  if ( t >= 0 )
240  {
241  pt = r.point( t );
242  return true; // you might want to allow an epsilon here too
243  }
244  }
245  return false;
246  }
247 
248 // copied from intersectsSegmentTriangle() from qt3d/src/render/backend/triangleboundingvolume.cpp
249 // by KDAB, licensed under the terms of LGPL
250  bool rayTriangleIntersection( const Ray3D &ray,
251  QVector3D a,
252  QVector3D b,
253  QVector3D c,
254  QVector3D &uvw,
255  float &t )
256  {
257  // Note: a, b, c in clockwise order
258  // RealTime Collision Detection page 192
259 
260  const QVector3D ab = b - a;
261  const QVector3D ac = c - a;
262  const QVector3D qp = ( ray.origin() - ray.point( ray.distance() ) );
263 
264  const QVector3D n = QVector3D::crossProduct( ab, ac );
265  const float d = QVector3D::dotProduct( qp, n );
266 
267  if ( d <= 0.0f )
268  return false;
269 
270  const QVector3D ap = ray.origin() - a;
271  t = QVector3D::dotProduct( ap, n );
272 
273  if ( t < 0.0f || t > d )
274  return false;
275 
276  const QVector3D e = QVector3D::crossProduct( qp, ap );
277  uvw.setY( QVector3D::dotProduct( ac, e ) );
278 
279  if ( uvw.y() < 0.0f || uvw.y() > d )
280  return false;
281 
282  uvw.setZ( -QVector3D::dotProduct( ab, e ) );
283 
284  if ( uvw.z() < 0.0f || uvw.y() + uvw.z() > d )
285  return false;
286 
287  const float ood = 1.0f / d;
288  t *= ood;
289  uvw.setY( uvw.y() * ood );
290  uvw.setZ( uvw.z() * ood );
291  uvw.setX( 1.0f - uvw.y() - uvw.z() );
292 
293  return true;
294  }
295 
296 }
297 
298 
300 
301 
302 
303 static QRect windowViewport( QSize area, const QRectF &relativeViewport )
304 {
305  if ( area.isValid() )
306  {
307  const int areaWidth = area.width();
308  const int areaHeight = area.height();
309  return QRect( relativeViewport.x() * areaWidth,
310  ( 1.0 - relativeViewport.y() - relativeViewport.height() ) * areaHeight,
311  relativeViewport.width() * areaWidth,
312  relativeViewport.height() * areaHeight );
313  }
314  return relativeViewport.toRect();
315 }
316 
317 
318 static QgsRayCastingUtils::Ray3D intersectionRay( QPointF pos, const QMatrix4x4 &viewMatrix,
319  const QMatrix4x4 &projectionMatrix, QRect viewport )
320 {
321  // if something seems wrong slightly off with the returned intersection rays,
322  // it may be the case that unproject() has hit qFuzzyIsNull() condition inside
323  // which has IMHO the threshold too high (1e-5) and can give problems when
324  // the camera is ~50km away or further.
325 
326  QVector3D nearPos = QVector3D( pos.x(), pos.y(), 0.0f );
327  nearPos = nearPos.unproject( viewMatrix, projectionMatrix, viewport );
328  QVector3D farPos = QVector3D( pos.x(), pos.y(), 1.0f );
329  farPos = farPos.unproject( viewMatrix, projectionMatrix, viewport );
330 
331  return QgsRayCastingUtils::Ray3D( nearPos,
332  ( farPos - nearPos ).normalized(),
333  ( farPos - nearPos ).length() );
334 }
335 
336 
337 namespace QgsRayCastingUtils
338 {
339 
340  Ray3D rayForViewportAndCamera( QSize area,
341  QPointF pos,
342  const QRectF &relativeViewport,
343  const Qt3DRender::QCamera *camera )
344  {
345  QMatrix4x4 viewMatrix = camera->viewMatrix();
346  QMatrix4x4 projectionMatrix = camera->projectionMatrix();
347  const QRect viewport = windowViewport( area, relativeViewport );
348 
349  // In GL the y is inverted compared to Qt
350  const QPointF glCorrectPos = QPointF( pos.x(), area.isValid() ? area.height() - pos.y() : pos.y() );
351  const auto ray = intersectionRay( glCorrectPos, viewMatrix, projectionMatrix, viewport );
352  return ray;
353  }
354 
355  Ray3D rayForCameraCenter( const Qt3DRender::QCamera *camera )
356  {
357  QMatrix4x4 inverse = QMatrix4x4( camera->projectionMatrix() * camera->viewMatrix() ).inverted();
358 
359  QVector4D vNear( 0.0, 0.0, -1.0, 1.0 );
360  QVector4D vFar( 0.0, 0.0, 1.0, 1.0 );
361  QVector4D nearPos4D = inverse * vNear;
362  QVector4D farPos4D = inverse * vFar;
363 
364  // the cases below hopefully should not happen and the check is here just as the last resort
365  // (with sensible camera matrix we should not hit singularities)
366  if ( qgsFloatNear( nearPos4D.w(), 0, 1e-10f ) )
367  nearPos4D.setW( 1 );
368  if ( qgsFloatNear( farPos4D.w(), 0, 1e-10f ) )
369  farPos4D.setW( 1 );
370 
371  QVector3D nearPos( ( nearPos4D / nearPos4D.w() ).toVector3D() );
372  QVector3D farPos( ( farPos4D / farPos4D.w() ).toVector3D() );
373 
374  return QgsRayCastingUtils::Ray3D( nearPos,
375  ( farPos - nearPos ).normalized(),
376  ( farPos - nearPos ).length() );
377  }
378 
379 }
380 
3 Axis-aligned bounding box - in world coords.
Definition: qgsaabb.h:30
bool qgsFloatNear(float a, float b, float epsilon=4 *FLT_EPSILON)
Compare two floats (but allow some difference)
Definition: qgis.h:277
bool operator==(const QgsFeatureIterator &fi1, const QgsFeatureIterator &fi2)
bool operator!=(const QgsFeatureIterator &fi1, const QgsFeatureIterator &fi2)
float zMax
Definition: qgsaabb.h:80
As part of the API refactoring and improvements which landed in the Processing API was substantially reworked from the x version This was done in order to allow much of the underlying Processing framework to be ported into c
float zMin
Definition: qgsaabb.h:77
float yMax
Definition: qgsaabb.h:79
std::ostream & operator<<(std::ostream &os, const QgsCoordinateReferenceSystem &r)
Output stream operator.
float xMin
Definition: qgsaabb.h:75
float xMax
Definition: qgsaabb.h:78
float yMin
Definition: qgsaabb.h:76