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
14namespace Mercator {
15//floor and ceil functions that return d-1 and d+1
16//respectively if d is integral
17static inline float gridceil(float d)
18{
19 float c = std::ceil(d);
20 return (c==d) ? c+1.0f : c;
21}
22
23static inline double gridceil(double d)
24{
25 auto c = std::ceil(d);
26 return (c==d) ? c+1.0 : c;
27}
28
29static inline float gridfloor(float d)
30{
31 float c = std::floor(d);
32 return (c==d) ? c-1.0f : c;
33}
34
35static 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
43bool 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
106static 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
118bool 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
125static 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
236bool 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:140