improved performance of SVBSP code
[divverent/darkplaces.git] / svbsp.c
1
2 // Shadow Volume BSP code written by Forest "LordHavoc" Hale on 2003-11-06 and placed into public domain.
3 // Modified by LordHavoc (to make it work and other nice things like that) on 2007-01-24 and 2007-01-25
4 // Optimized by LordHavoc on 2009-12-24 and 2009-12-25
5
6 #include <math.h>
7 #include <string.h>
8 #include "svbsp.h"
9 #include "polygon.h"
10
11 #define MAX_SVBSP_POLYGONPOINTS 64
12 #define SVBSP_CLIP_EPSILON (1.0f / 1024.0f)
13
14 #define SVBSP_DotProduct(a,b) ((a)[0]*(b)[0]+(a)[1]*(b)[1]+(a)[2]*(b)[2])
15
16 typedef struct svbsp_polygon_s
17 {
18         float points[MAX_SVBSP_POLYGONPOINTS][3];
19         //unsigned char splitflags[MAX_SVBSP_POLYGONPOINTS];
20         int facesplitflag;
21         int numpoints;
22 }
23 svbsp_polygon_t;
24
25 static void SVBSP_PlaneFromPoints(float *plane4f, const float *p1, const float *p2, const float *p3)
26 {
27         float ilength;
28         // calculate unnormalized plane
29         plane4f[0] = (p1[1] - p2[1]) * (p3[2] - p2[2]) - (p1[2] - p2[2]) * (p3[1] - p2[1]);
30         plane4f[1] = (p1[2] - p2[2]) * (p3[0] - p2[0]) - (p1[0] - p2[0]) * (p3[2] - p2[2]);
31         plane4f[2] = (p1[0] - p2[0]) * (p3[1] - p2[1]) - (p1[1] - p2[1]) * (p3[0] - p2[0]);
32         plane4f[3] = SVBSP_DotProduct(plane4f, p1);
33         // normalize the plane normal and adjust distance accordingly
34         ilength = (float)sqrt(SVBSP_DotProduct(plane4f, plane4f));
35         if (ilength)
36                 ilength = 1.0f / ilength;
37         plane4f[0] *= ilength;
38         plane4f[1] *= ilength;
39         plane4f[2] *= ilength;
40         plane4f[3] *= ilength;
41 }
42
43 static void SVBSP_DividePolygon(int innumpoints, const float *inpoints, float planenormalx, float planenormaly, float planenormalz, float planedist, float epsilon, int outfrontmaxpoints, float *outfrontpoints, int *neededfrontpoints, int outbackmaxpoints, float *outbackpoints, int *neededbackpoints, int *oncountpointer)
44 {
45         int i, frontcount = 0, backcount = 0, oncount = 0;
46         const float *n, *p;
47         float frac, pdist, ndist;
48         if (innumpoints)
49         {
50                 n = inpoints;
51                 ndist = n[0] * planenormalx + n[1] * planenormaly + n[2] * planenormalz - planedist;
52                 for(i = 0;i < innumpoints;i++)
53                 {
54                         p = n;
55                         pdist = ndist;
56                         n = inpoints + ((i + 1) < innumpoints ? (i + 1) : 0) * 3;
57                         ndist = n[0] * planenormalx + n[1] * planenormaly + n[2] * planenormalz - planedist;
58                         if (pdist >= -epsilon)
59                         {
60                                 if (pdist <= epsilon)
61                                         oncount++;
62                                 if (frontcount < outfrontmaxpoints)
63                                 {
64                                         *outfrontpoints++ = p[0];
65                                         *outfrontpoints++ = p[1];
66                                         *outfrontpoints++ = p[2];
67                                 }
68                                 frontcount++;
69                         }
70                         if (pdist <= epsilon)
71                         {
72                                 if (backcount < outbackmaxpoints)
73                                 {
74                                         *outbackpoints++ = p[0];
75                                         *outbackpoints++ = p[1];
76                                         *outbackpoints++ = p[2];
77                                 }
78                                 backcount++;
79                         }
80                         if ((pdist > epsilon && ndist < -epsilon) || (pdist < -epsilon && ndist > epsilon))
81                         {
82                                 oncount++;
83                                 frac = pdist / (pdist - ndist);
84                                 if (frontcount < outfrontmaxpoints)
85                                 {
86                                         *outfrontpoints++ = p[0] + frac * (n[0] - p[0]);
87                                         *outfrontpoints++ = p[1] + frac * (n[1] - p[1]);
88                                         *outfrontpoints++ = p[2] + frac * (n[2] - p[2]);
89                                 }
90                                 frontcount++;
91                                 if (backcount < outbackmaxpoints)
92                                 {
93                                         *outbackpoints++ = p[0] + frac * (n[0] - p[0]);
94                                         *outbackpoints++ = p[1] + frac * (n[1] - p[1]);
95                                         *outbackpoints++ = p[2] + frac * (n[2] - p[2]);
96                                 }
97                                 backcount++;
98                         }
99                 }
100         }
101         if (neededfrontpoints)
102                 *neededfrontpoints = frontcount;
103         if (neededbackpoints)
104                 *neededbackpoints = backcount;
105         if (oncountpointer)
106                 *oncountpointer = oncount;
107 }
108
109 void SVBSP_Init(svbsp_t *b, const float *origin, int maxnodes, svbsp_node_t *nodes)
110 {
111         memset(b, 0, sizeof(*b));
112         b->origin[0] = origin[0];
113         b->origin[1] = origin[1];
114         b->origin[2] = origin[2];
115         b->numnodes = 3;
116         b->maxnodes = maxnodes;
117         b->nodes = nodes;
118         b->ranoutofnodes = 0;
119         b->stat_occluders_rejected = 0;
120         b->stat_occluders_accepted = 0;
121         b->stat_occluders_fragments_accepted = 0;
122         b->stat_occluders_fragments_rejected = 0;
123         b->stat_queries_rejected = 0;
124         b->stat_queries_accepted = 0;
125         b->stat_queries_fragments_accepted = 0;
126         b->stat_queries_fragments_rejected = 0;
127
128         // the bsp tree must be initialized to have two perpendicular splits axes
129         // through origin, otherwise the polygon insertions would affect the
130         // opposite side of the tree, which would be disasterous.
131         //
132         // so this code has to make 3 nodes and 4 leafs, and since the leafs are
133         // represented by empty/solid state numbers in this system rather than
134         // actual structs, we only need to make the 3 nodes.
135
136         // root node
137         // this one splits the world into +X and -X sides
138         b->nodes[0].plane[0] = 1;
139         b->nodes[0].plane[1] = 0;
140         b->nodes[0].plane[2] = 0;
141         b->nodes[0].plane[3] = origin[0];
142         b->nodes[0].parent = -1;
143         b->nodes[0].children[0] = 1;
144         b->nodes[0].children[1] = 2;
145
146         // +X side node
147         // this one splits the +X half of the world into +Y and -Y
148         b->nodes[1].plane[0] = 0;
149         b->nodes[1].plane[1] = 1;
150         b->nodes[1].plane[2] = 0;
151         b->nodes[1].plane[3] = origin[1];
152         b->nodes[1].parent = 0;
153         b->nodes[1].children[0] = -1;
154         b->nodes[1].children[1] = -1;
155
156         // -X side node
157         // this one splits the -X half of the world into +Y and -Y
158         b->nodes[2].plane[0] = 0;
159         b->nodes[2].plane[1] = 1;
160         b->nodes[2].plane[2] = 0;
161         b->nodes[2].plane[3] = origin[1];
162         b->nodes[2].parent = 0;
163         b->nodes[2].children[0] = -1;
164         b->nodes[2].children[1] = -1;
165 }
166
167 static void SVBSP_InsertOccluderPolygonNodes(svbsp_t *b, int *parentnodenumpointer, int parentnodenum, const svbsp_polygon_t *poly, void (*fragmentcallback)(void *fragmentcallback_pointer1, int fragmentcallback_number1, svbsp_t *b, int numpoints, const float *points), void *fragmentcallback_pointer1, int fragmentcallback_number1)
168 {
169         // now we need to create up to numpoints + 1 new nodes, forming a BSP tree
170         // describing the occluder polygon's shadow volume
171         int i, j, p, basenum;
172         svbsp_node_t *node;
173
174         // if there aren't enough nodes remaining, skip it
175         if (b->numnodes + poly->numpoints + 1 >= b->maxnodes)
176         {
177                 b->ranoutofnodes = 1;
178                 return;
179         }
180
181         // add one node per side, then the actual occluding face node
182
183         // thread safety notes:
184         // DO NOT multithread insertion, it could be made 'safe' but the results
185         // would be inconsistent.
186         //
187         // it is completely safe to multithread queries in all cases.
188         //
189         // if an insertion is occurring the query will give intermediate results,
190         // being blocked by some volumes but not others, which is perfectly okay
191         // for visibility culling intended only to reduce rendering work
192
193         // note down the first available nodenum for the *parentnodenumpointer
194         // line which is done last to allow multithreaded queries during an
195         // insertion
196         basenum = b->numnodes;
197         for (i = 0, p = poly->numpoints - 1;i < poly->numpoints;p = i, i++)
198         {
199 #if 1
200                 // see if a parent plane describes this side
201                 for (j = parentnodenum;j >= 0;j = b->nodes[j].parent)
202                 {
203                         float *parentnodeplane = b->nodes[j].plane;
204                         if (fabs(SVBSP_DotProduct(poly->points[p], parentnodeplane) - parentnodeplane[3]) < SVBSP_CLIP_EPSILON
205                          && fabs(SVBSP_DotProduct(poly->points[i], parentnodeplane) - parentnodeplane[3]) < SVBSP_CLIP_EPSILON
206                          && fabs(SVBSP_DotProduct(b->origin      , parentnodeplane) - parentnodeplane[3]) < SVBSP_CLIP_EPSILON)
207                                 break;
208                 }
209                 if (j >= 0)
210                         continue; // already have a matching parent plane
211 #endif
212 #if 0
213                 // skip any sides that were classified as belonging to a parent plane
214                 if (poly->splitflags[i])
215                         continue;
216 #endif
217                 // create a side plane
218                 // anything infront of this is not inside the shadow volume
219                 node = b->nodes + b->numnodes++;
220                 SVBSP_PlaneFromPoints(node->plane, b->origin, poly->points[p], poly->points[i]);
221                 // we need to flip the plane if it puts any part of the polygon on the
222                 // wrong side
223                 // (in this way this code treats all polygons as float sided)
224                 //
225                 // because speed is important this stops as soon as it finds proof
226                 // that the orientation is right or wrong
227                 // (we know that the plane is on one edge of the polygon, so there is
228                 // never a case where points lie on both sides, so the first hint is
229                 // sufficient)
230                 for (j = 0;j < poly->numpoints;j++)
231                 {
232                         float d = SVBSP_DotProduct(poly->points[j], node->plane) - node->plane[3];
233                         if (d < -SVBSP_CLIP_EPSILON)
234                                 break;
235                         if (d > SVBSP_CLIP_EPSILON)
236                         {
237                                 node->plane[0] *= -1;
238                                 node->plane[1] *= -1;
239                                 node->plane[2] *= -1;
240                                 node->plane[3] *= -1;
241                                 break;
242                         }
243                 }
244                 node->parent = parentnodenum;
245                 node->children[0] = -1; // empty
246                 node->children[1] = -1; // empty
247                 // link this child into the tree
248                 *parentnodenumpointer = parentnodenum = (int)(node - b->nodes);
249                 // now point to the child pointer for the next node to update later
250                 parentnodenumpointer = &node->children[1];
251         }
252
253 #if 1
254         // skip the face plane if it lies on a parent plane
255         if (!poly->facesplitflag)
256 #endif
257         {
258                 // add the face-plane node
259                 // infront is empty, behind is shadow
260                 node = b->nodes + b->numnodes++;
261                 SVBSP_PlaneFromPoints(node->plane, poly->points[0], poly->points[1], poly->points[2]);
262                 // this is a flip check similar to the one above
263                 // this one checks if the plane faces the origin, if not, flip it
264                 if (SVBSP_DotProduct(b->origin, node->plane) - node->plane[3] < -SVBSP_CLIP_EPSILON)
265                 {
266                         node->plane[0] *= -1;
267                         node->plane[1] *= -1;
268                         node->plane[2] *= -1;
269                         node->plane[3] *= -1;
270                 }
271                 node->parent = parentnodenum;
272                 node->children[0] = -1; // empty
273                 node->children[1] = -2; // shadow
274                 // link this child into the tree
275                 // (with the addition of this node, queries will now be culled by it)
276                 *parentnodenumpointer = (int)(node - b->nodes);
277         }
278 }
279
280 static int SVBSP_AddPolygonNode(svbsp_t *b, int *parentnodenumpointer, int parentnodenum, const svbsp_polygon_t *poly, int insertoccluder, void (*fragmentcallback)(void *fragmentcallback_pointer1, int fragmentcallback_number1, svbsp_t *b, int numpoints, const float *points), void *fragmentcallback_pointer1, int fragmentcallback_number1)
281 {
282         int i;
283         int countonpoints;
284         int facesplitflag = poly->facesplitflag;
285         float plane[4];
286         svbsp_polygon_t front;
287         svbsp_polygon_t back;
288         svbsp_node_t *node;
289         if (poly->numpoints < 1)
290                 return 0;
291         // recurse through plane nodes
292         while (*parentnodenumpointer >= 0)
293         {
294                 // do a quick check to see if there is any need to split the polygon
295                 node = b->nodes + *parentnodenumpointer;
296                 parentnodenum = *parentnodenumpointer;
297                 plane[0] = node->plane[0];
298                 plane[1] = node->plane[1];
299                 plane[2] = node->plane[2];
300                 plane[3] = node->plane[3];
301 #if 1
302                 {
303                         int onback;
304                         int onfront;
305                         float d;
306                         onback = 0;
307                         onfront = 0;
308                         for (i = 0;i < poly->numpoints;i++)
309                         {
310                                 d = SVBSP_DotProduct(poly->points[i], plane) - plane[3];
311                                 if (d <= -SVBSP_CLIP_EPSILON)
312                                         onback = 1;
313                                 if (d >= SVBSP_CLIP_EPSILON)
314                                         onfront = 1;
315                         }
316                         if (onfront)
317                         {
318                                 if (!onback)
319                                 {
320                                         // no need to split, just go to one side
321                                         parentnodenumpointer = &node->children[0];
322                                         continue;
323                                 }
324                         }
325                         else if (onback)
326                         {
327                                 // no need to split, just go to one side
328                                 parentnodenumpointer = &node->children[1];
329                                 continue;
330                         }
331                         else
332                         {
333                                 // no need to split, this polygon is on the plane
334                                 // this case only occurs for polygons on the face plane, usually
335                                 // the same polygon (inserted twice - once as occluder, once as
336                                 // tester)
337                                 // if this is an occluder, it is redundant
338                                 if (insertoccluder)
339                                         return 1; // occluded
340                                 // if this is a tester, test the front side, because it is
341                                 // probably the same polygon that created this node...
342                                 facesplitflag = 1;
343                                 parentnodenumpointer = &node->children[0];
344                                 continue;
345                         }
346                 }
347                 // at this point we know it crosses the plane, so we need to split it
348                 SVBSP_DividePolygon(poly->numpoints, poly->points[0], plane[0], plane[1], plane[2], plane[3], SVBSP_CLIP_EPSILON, MAX_SVBSP_POLYGONPOINTS, front.points[0], &front.numpoints, MAX_SVBSP_POLYGONPOINTS, back.points[0], &back.numpoints, &countonpoints);
349 #else
350                 // at this point we know it crosses the plane, so we need to split it
351                 SVBSP_DividePolygon(poly->numpoints, poly->points[0], plane[0], plane[1], plane[2], plane[3], SVBSP_CLIP_EPSILON, MAX_SVBSP_POLYGONPOINTS, front.points[0], &front.numpoints, MAX_SVBSP_POLYGONPOINTS, back.points[0], &back.numpoints, &countonpoints);
352                 if (countonpoints == poly->numpoints)
353                 {
354                         // polygon lies on plane - this only occurs on face planes
355                         // if it is an occluder, it is redundant
356                         if (insertoccluder)
357                                 return 1; // occluded
358                         // if it is a tester, it should take the front side only
359                         facesplitflag = 1;
360                         parentnodenumpointer = &node->children[0];
361                         continue;
362                 }
363 #endif
364                 if (front.numpoints > MAX_SVBSP_POLYGONPOINTS)
365                         front.numpoints = MAX_SVBSP_POLYGONPOINTS;
366                 front.facesplitflag = facesplitflag;
367                 if (back.numpoints > MAX_SVBSP_POLYGONPOINTS)
368                         back.numpoints = MAX_SVBSP_POLYGONPOINTS;
369                 back.facesplitflag = facesplitflag;
370 #if 0
371                 if (insertoccluder)
372                 {
373                         int pnum;
374                         int f;
375                         int pf;
376                         svbsp_node_t *pnode;
377                         // set splitflags based on this node and parent nodes
378                         for (i = 0;i < front.numpoints;i++)
379                                 front.splitflags[i] = 0;
380                         for (i = 0;i < front.numpoints;i++)
381                                 back.splitflags[i] = 0;
382                         pnum = *parentnodenumpointer;
383                         while (pnum >= 0)
384                         {
385                                 pnode = b->nodes + pnum;
386                                 plane[0] = pnode->plane[0];
387                                 plane[1] = pnode->plane[1];
388                                 plane[2] = pnode->plane[2];
389                                 plane[3] = pnode->plane[3];
390                                 if (front.numpoints)
391                                 {
392                                         i = front.numpoints-1;
393                                         pf = fabs(SVBSP_DotProduct(front.points[i], plane) - plane[3]) <= SVBSP_CLIP_EPSILON;
394                                         for (i = 0;i < front.numpoints;pf = f, i++)
395                                         {
396                                                 f = fabs(SVBSP_DotProduct(front.points[i], plane) - plane[3]) <= SVBSP_CLIP_EPSILON;
397                                                 front.splitflags[i] |= (f & pf);
398                                         }
399                                 }
400                                 if (back.numpoints)
401                                 {
402                                         i = back.numpoints-1;
403                                         pf = fabs(SVBSP_DotProduct(back.points[i], plane) - plane[3]) <= SVBSP_CLIP_EPSILON;
404                                         for (i = 0;i < back.numpoints;pf = f, i++)
405                                         {
406                                                 f = fabs(SVBSP_DotProduct(back.points[i], plane) - plane[3]) <= SVBSP_CLIP_EPSILON;
407                                                 back.splitflags[i] |= (f & pf);
408                                         }
409                                 }
410                                 pnum = pnode->parent;
411                         }
412                 }
413 #endif
414                 // recurse the sides and return the resulting occlusion flags
415                 i  = SVBSP_AddPolygonNode(b, &node->children[0], *parentnodenumpointer, &front, insertoccluder, fragmentcallback, fragmentcallback_pointer1, fragmentcallback_number1);
416                 i |= SVBSP_AddPolygonNode(b, &node->children[1], *parentnodenumpointer, &back , insertoccluder, fragmentcallback, fragmentcallback_pointer1, fragmentcallback_number1);
417                 return i;
418         }
419         // leaf node
420         if (*parentnodenumpointer == -1)
421         {
422                 // empty leaf node; and some geometry survived
423                 // if inserting an occluder, replace this empty leaf with a shadow volume
424 #if 0
425                 for (i = 0;i < poly->numpoints-2;i++)
426                 {
427                         Debug_PolygonBegin(NULL, DRAWFLAG_ADDITIVE);
428                         Debug_PolygonVertex(poly->points[  0][0], poly->points[  0][1], poly->points[  0][2], 0.0f, 0.0f, 0.25f, 0.0f, 0.0f, 1.0f);
429                         Debug_PolygonVertex(poly->points[i+1][0], poly->points[i+1][1], poly->points[i+1][2], 0.0f, 0.0f, 0.25f, 0.0f, 0.0f, 1.0f);
430                         Debug_PolygonVertex(poly->points[i+2][0], poly->points[i+2][1], poly->points[i+2][2], 0.0f, 0.0f, 0.25f, 0.0f, 0.0f, 1.0f);
431                         Debug_PolygonEnd();
432                 }
433 #endif
434                 if (insertoccluder)
435                 {
436                         b->stat_occluders_fragments_accepted++;
437                         SVBSP_InsertOccluderPolygonNodes(b, parentnodenumpointer, parentnodenum, poly, fragmentcallback, fragmentcallback_pointer1, fragmentcallback_number1);
438                 }
439                 else
440                         b->stat_queries_fragments_accepted++;
441                 if (fragmentcallback)
442                         fragmentcallback(fragmentcallback_pointer1, fragmentcallback_number1, b, poly->numpoints, poly->points[0]);
443                 return 2;
444         }
445         else
446         {
447                 // otherwise it's a solid leaf which destroys all polygons inside it
448                 if (insertoccluder)
449                         b->stat_occluders_fragments_rejected++;
450                 else
451                         b->stat_queries_fragments_rejected++;
452 #if 0
453                 for (i = 0;i < poly->numpoints-2;i++)
454                 {
455                         Debug_PolygonBegin(NULL, DRAWFLAG_ADDITIVE);
456                         Debug_PolygonVertex(poly->points[  0][0], poly->points[  0][1], poly->points[  0][2], 0.0f, 0.0f, 0.0f, 0.0f, 0.25f, 1.0f);
457                         Debug_PolygonVertex(poly->points[i+1][0], poly->points[i+1][1], poly->points[i+1][2], 0.0f, 0.0f, 0.0f, 0.0f, 0.25f, 1.0f);
458                         Debug_PolygonVertex(poly->points[i+2][0], poly->points[i+2][1], poly->points[i+2][2], 0.0f, 0.0f, 0.0f, 0.0f, 0.25f, 1.0f);
459                         Debug_PolygonEnd();
460                 }
461 #endif
462         }
463         return 1;
464 }
465
466 int SVBSP_AddPolygon(svbsp_t *b, int numpoints, const float *points, int insertoccluder, void (*fragmentcallback)(void *fragmentcallback_pointer1, int fragmentcallback_number1, svbsp_t *b, int numpoints, const float *points), void *fragmentcallback_pointer1, int fragmentcallback_number1)
467 {
468         int i;
469         int nodenum;
470         svbsp_polygon_t poly;
471         // don't even consider an empty polygon
472         // note we still allow points and lines to be tested...
473         if (numpoints < 1)
474                 return 0;
475         poly.numpoints = numpoints;
476         for (i = 0;i < numpoints;i++)
477         {
478                 poly.points[i][0] = points[i*3+0];
479                 poly.points[i][1] = points[i*3+1];
480                 poly.points[i][2] = points[i*3+2];
481                 //poly.splitflags[i] = 0; // this edge is a valid BSP splitter - clipped edges are not (because they lie on a bsp plane)
482                 poly.facesplitflag = 0; // this face is a valid BSP Splitter - if it lies on a bsp plane it is not
483         }
484 #if 0
485 //if (insertoccluder)
486         for (i = 0;i < poly.numpoints-2;i++)
487         {
488                 Debug_PolygonBegin(NULL, DRAWFLAG_ADDITIVE);
489                 Debug_PolygonVertex(poly.points[  0][0], poly.points[  0][1], poly.points[  0][2], 0.0f, 0.0f, 0.0f, 0.25f, 0.0f, 1.0f);
490                 Debug_PolygonVertex(poly.points[i+1][0], poly.points[i+1][1], poly.points[i+1][2], 0.0f, 0.0f, 0.0f, 0.25f, 0.0f, 1.0f);
491                 Debug_PolygonVertex(poly.points[i+2][0], poly.points[i+2][1], poly.points[i+2][2], 0.0f, 0.0f, 0.0f, 0.25f, 0.0f, 1.0f);
492                 Debug_PolygonEnd();
493         }
494 #endif
495         nodenum = 0;
496         i = SVBSP_AddPolygonNode(b, &nodenum, -1, &poly, insertoccluder, fragmentcallback, fragmentcallback_pointer1, fragmentcallback_number1);
497         if (insertoccluder)
498         {
499                 if (i & 2)
500                         b->stat_occluders_accepted++;
501                 else
502                         b->stat_occluders_rejected++;
503         }
504         else
505         {
506                 if (i & 2)
507                         b->stat_queries_accepted++;
508                 else
509                         b->stat_queries_rejected++;
510         }
511         return i;
512 }
513