mercator  0.4.0
A terrain generation library for the Worldforge system.
Intersect.cpp
1 // This file may be redistributed and modified only under the terms of
2 // the GNU General Public License (See COPYING for details).
3 // Copyright (C) 2003 Damien McGinnes
4 
5 #ifdef HAVE_CONFIG_H
6 #include "config.h"
7 #endif
8 
9 #include "Intersect.h"
10 #include "Segment.h"
11 
12 #include <algorithm>
13 
14 namespace Mercator {
15 //floor and ceil functions that return d-1 and d+1
16 //respectively if d is integral
17 static inline float gridceil(float d)
18 {
19  float c = std::ceil(d);
20  return (c==d) ? c+1.0f : c;
21 }
22 
23 static inline double gridceil(double d)
24 {
25  auto c = std::ceil(d);
26  return (c==d) ? c+1.0 : c;
27 }
28 
29 static inline float gridfloor(float d)
30 {
31  float c = std::floor(d);
32  return (c==d) ? c-1.0f : c;
33 }
34 
35 static inline double gridfloor(double d)
36 {
37  auto c = std::floor(d);
38  return (c==d) ? c-1.0 : c;
39 }
40 
41 
42 //check intersection of an axis-aligned box with the terrain
43 bool Intersect(const Terrain &t, const WFMath::AxisBox<3> &bbox)
44 {
45  float max, min=bbox.lowCorner()[1];
46  const int res = t.getResolution();
47  const float spacing = t.getSpacing();
48 
49  //determine which segments are involved
50  //usually will just be one
51  int xlow = (int) floor(bbox.lowCorner()[0] / spacing);
52  int xhigh = (int) gridceil(bbox.highCorner()[0] / spacing);
53  int zlow = (int) floor(bbox.lowCorner()[2] / spacing);
54  int zhigh = (int) gridceil(bbox.highCorner()[2] / spacing);
55 
56  //loop across all tiles covered by this bbox
57  for (int x = xlow; x < xhigh; x++) {
58  for (int z = zlow; z < zhigh; z++) {
59  //check the bbox against the extent of each tile
60  //as an early rejection
61  Segment *thisSeg=t.getSegmentAtIndex(x,z);
62 
63  if (thisSeg)
64  max=thisSeg->getMax();
65  else
67 
68  if (max > min) {
69  //entity bbox overlaps with the extents of this tile
70  //now check each tile point covered by the entity bbox
71 
72  //clip the points to be tested against the bbox
73  int min_x = (int) floor(bbox.lowCorner()[0] - ((float)x * spacing));
74  if (min_x < 0) min_x = 0;
75 
76  int max_x = (int) gridceil(bbox.highCorner()[0] - ((float)x * spacing));
77  if (max_x > res) min_x = res;
78 
79  int min_z = (int) floor(bbox.lowCorner()[2] - ((float)z * spacing));
80  if (min_z < 0) min_z = 0;
81 
82  int max_z = (int) gridceil(bbox.highCorner()[2] - ((float)z * spacing));
83  if (max_z > res) min_z = res;
84 
85  //loop over each point and see if it is greater than the minimum
86  //of the bbox. If all points are below, the the bbox does NOT
87  //intersect. If a single point is above, then the bbox MIGHT
88  //intersect.
89  for (int xpt = min_x; xpt <= max_x; xpt++) {
90  for (int zpt = min_z; zpt <= max_z; zpt++) {
91  if (thisSeg) {
92  if (thisSeg->get(xpt,zpt) > min) {
93  return true;
94  }
95  } else if (Terrain::defaultLevel > min) {
96  return true;
97  }
98  }
99  }
100  }
101  }
102  }
103  return false;
104 }
105 
106 static bool HOT(const Terrain &t, const WFMath::Point<3> &pt, double & h)
107 {
108  WFMath::Vector<3> normal;
109  float terrHeight;
110  if (!t.getHeightAndNormal(pt[0], pt[2], terrHeight, normal)) {
111  return false;
112  }
113 
114  h = (pt[1] - terrHeight);
115  return true;
116 }
117 
118 bool Intersect(const Terrain &t, const WFMath::Point<3> &pt)
119 {
120  double h;
121  return HOT(t, pt, h) && h <= 0.0;
122 }
123 
124 //helper function for ray terrain intersection
125 static bool cellIntersect(double h1, double h2, double h3, double h4, double X, double Z,
126  const WFMath::Vector<3> &nDir, float dirLen,
127  const WFMath::Point<3> &sPt, WFMath::Point<3> &intersection,
128  WFMath::Vector<3> &normal, double &par)
129 {
130  //ray plane intersection roughly using the following:
131  //parametric ray eqn: p=po + par V
132  //plane eqn: p dot N + d = 0
133  //
134  //intersection:
135  // -par = (po dot N + d ) / (V dot N)
136  //
137  //
138  // effectively we calculate the ray parametric variable for the
139  // intersection of the plane corresponding to each triangle
140  // then clip them by endpints of the ray, and by the sides of the square
141  // and by the diagonal
142  //
143  // if they both still intersect, then we choose the earlier intersection
144 
145 
146  //intersection points for top and bottom triangles
147  WFMath::Point<3> topInt, botInt;
148 
149  //point to use in plane equation for both triangles
150  WFMath::Vector<3> p0 = WFMath::Vector<3>(X, h1, Z);
151 
152  // square is broken into two triangles
153  // top triangle |/
154  bool topIntersected = false;
155  WFMath::Vector<3> topNormal(h2-h3, 1.0, h1-h2);
156  topNormal.normalize();
157  double t = Dot(nDir, topNormal);
158 
159  double topP=0.0;
160 
161  if ((t > 1e-7) || (t < -1e-7)) {
162  topP = - (Dot((sPt-WFMath::Point<3>(0,0,0)), topNormal)
163  - Dot(topNormal, p0)) / t;
164  topInt = sPt + nDir*topP;
165  //check the intersection is inside the triangle, and within the ray extents
166  if ((topP <= dirLen) && (topP > 0.0) &&
167  (topInt[0] >= X ) && (topInt[2] <= Z + 1 ) &&
168  ((topInt[0] - topInt[2]) <= (X - Z)) ) {
169  topIntersected=true;
170  }
171  }
172 
173  // bottom triangle /|
174  bool botIntersected = false;
175  WFMath::Vector<3> botNormal(h1-h4, 1.0, h4-h3);
176  botNormal.normalize();
177  double b = Dot(nDir, botNormal);
178  double botP=0.0;
179 
180  if ((b > 1e-7) || (b < -1e-7)) {
181  botP = - (Dot((sPt-WFMath::Point<3>(0,0,0)), botNormal)
182  - Dot(botNormal, p0)) / b;
183  botInt = sPt + nDir*botP;
184  //check the intersection is inside the triangle, and within the ray extents
185  if ((botP <= dirLen) && (botP > 0.0) &&
186  (botInt[0] <= X + 1 ) && (botInt[2] >= Z ) &&
187  ((botInt[0] - botInt[2]) >= (X - Z)) ) {
188  botIntersected = true;
189  }
190  }
191 
192  if (topIntersected && botIntersected) { //intersection with both
193  if (botP <= topP) {
194  intersection = botInt;
195  normal = botNormal;
196  par=botP/dirLen;
197  if (botP == topP) {
198  normal += topNormal;
199  normal.normalize();
200  }
201  return true;
202  }
203  else {
204  intersection = topInt;
205  normal = topNormal;
206  par=topP/dirLen;
207  return true;
208  }
209  }
210  else if (topIntersected) { //intersection with top
211  intersection = topInt;
212  normal = topNormal;
213  par=topP/dirLen;
214  return true;
215  }
216  else if (botIntersected) { //intersection with bot
217  intersection = botInt;
218  normal = botNormal;
219  par=botP/dirLen;
220  return true;
221  }
222 
223  return false;
224 }
225 
226 //Intersection of ray with terrain
227 //
228 //returns true on successful intersection
229 //ray is represented by a start point (sPt)
230 //and a direction vector (dir). The length of dir is
231 //used as the length of the ray. (ie not an infinite ray)
232 //intersection and normal are filled in on a successful result.
233 //par is the distance along the ray where intersection occurs
234 //0.0 == start, 1.0 == end.
235 
236 bool Intersect(const Terrain &t, const WFMath::Point<3> &sPt, const WFMath::Vector<3>& dir,
237  WFMath::Point<3> &intersection, WFMath::Vector<3> &normal, double &par)
238 {
239  //FIXME early reject using segment max
240  //FIXME handle height point getting in a more optimal way
241  //FIXME early reject for vertical ray
242 
243 
244  // check if startpoint is below ground
245  double hot;
246  if (HOT(t, sPt, hot) && hot < 0.0) return true;
247 
248  double paraX=0.0, paraZ=0.0; //used to store the parametric gap between grid crossings
249  double pX, pZ; //the accumulators for the parametrics as we traverse the ray
250  float h1,h2,h3,h4;
251 
252  WFMath::Point<3> last(sPt), next(sPt);
253  WFMath::Vector<3> nDir(dir);
254  nDir.normalize();
255  float dirLen = dir.mag();
256 
257  //work out where the ray first crosses an X grid line
258  if (dir[0] != 0.0f) {
259  paraX = 1.0f/dir[0];
260  double crossX = (dir[0] > 0.0f) ? gridceil(last[0]) : gridfloor(last[0]);
261  pX = (crossX - last[0]) * paraX;
262  pX = std::min(pX, 1.0);
263  }
264  else { //parallel: never crosses
265  pX = 1.0f;
266  }
267 
268  //work out where the ray first crosses a Y grid line
269  if (dir[2] != 0.0f) {
270  paraZ = 1.0f/dir[2];
271  double crossZ = (dir[2] > 0.0f) ? gridceil(last[2]) : gridfloor(last[2]);
272  pZ = (crossZ - sPt[2]) * paraZ;
273  pZ = std::min(pZ, 1.0);
274  }
275  else { //parallel: never crosses
276  pZ = 1.0f;
277  }
278 
279  //ensure we traverse the ray forwards
280  paraX = std::abs(paraX);
281  paraZ = std::abs(paraZ);
282 
283  bool endpoint = false;
284  //check each candidate tile for an intersection
285  while (true) {
286  last = next;
287  if (pX < pZ) { // cross x grid line first
288  next = sPt + (pX * dir);
289  pX += paraX; // set x accumulator to current p
290  }
291  else { //cross z grid line first
292  next = sPt + (pZ * dir);
293  if (pX == pZ) {
294  pX += paraX; //unusual case where ray crosses corner
295  }
296  pZ += paraZ; // set z accumulator to current p
297  }
298 
299  //FIXME these gets could be optimized a bit
300  float x= (dir[0] > 0) ? std::floor(last[0]) : std::floor(next[0]);
301  float z= (dir[2] > 0) ? std::floor(last[2]) : std::floor(next[2]);
302  h1 = t.get(x, z);
303  h2 = t.get(x, z+1);
304  h3 = t.get(x+1, z+1);
305  h4 = t.get(x+1, z);
306  auto height = std::max( std::max(h1, h2),
307  std::max(h3, h4));
308 
309  if ( (last[1] < height) || (next[1] < height) ) {
310  // possible intersect with this tile
311  if (cellIntersect(h1, h2, h3, h4, x, z, nDir, dirLen, sPt,
312  intersection, normal, par)) {
313  return true;
314  }
315  }
316 
317  if ((pX >= 1.0f) && (pZ >= 1.0f)) {
318  if (endpoint) {
319  break;
320  }
321  else endpoint = true;
322  }
323  }
324 
325  return false;
326 }
327 
328 } // namespace Mercator
static constexpr float defaultLevel
Height value used when no data is available.
Definition: Terrain.h:141