DRW: Initial implementation of Frustum culling.

This is very efficient and add a pretty low overhead (0.1ms of drawing time for 10K objects passing through all tests, on my i3-4100M).
The like the rest of the DRWCallState, test is "cached" until the view matrices changes.

4 files changed, 210 insertions, 17 deletions

diff --git a/source/blender/draw/intern/draw_manager.c b/source/blender/draw/intern/draw_manager.c
index 2ddc977941a..62042a44434 100644
--- a/source/blender/draw/intern/draw_manager.c
+++ b/source/blender/draw/intern/draw_manager.c
@@ -441,6 +441,8 @@ static void drw_viewport_var_init(void)
 	DST.dirty_mat = false;
 	DST.state_cache_id = 1;
 
+	DST.clipping.updated = false;
+
 	memset(DST.common_instance_data, 0x0, sizeof(DST.common_instance_data));
 }
 
diff --git a/source/blender/draw/intern/draw_manager.h b/source/blender/draw/intern/draw_manager.h
index cfc49320f3b..37646ec6dbb 100644
--- a/source/blender/draw/intern/draw_manager.h
+++ b/source/blender/draw/intern/draw_manager.h
@@ -76,6 +76,12 @@
 
 #endif  /* USE_PROFILE */
 
+/* TODO Put it somewhere else? */
+typedef struct BoundSphere {
+	float center[3], radius;
+} BoundSphere;
+
+
 /* ------------ Data Structure --------------- */
 /**
  * Data structure containing all drawcalls organized by passes and materials.
@@ -107,9 +113,7 @@ typedef struct DRWCallState {
 	uint16_t matflag;         /* Which matrices to compute. */
 	/* Culling: Using Bounding Sphere for now for faster culling.
 	 * Not ideal for planes. */
-	struct {
-		float loc[3], rad; /* Bypassed if radius is < 0.0. */
-	} bsphere;
+	BoundSphere bsphere;
 	/* Matrices */
 	float model[4][4];
 	float modelinverse[4][4];
@@ -305,6 +309,12 @@ typedef struct DRWManager {
 		float clip_planes_eq[MAX_CLIP_PLANES][4];
 	} view_data;
 
+	struct {
+		float frustum_planes[6][4];
+		BoundSphere frustum_bsphere;
+		bool updated;
+	} clipping;
+
 #ifdef USE_GPU_SELECT
 	unsigned int select_id;
 #endif
diff --git a/source/blender/draw/intern/draw_manager_data.c b/source/blender/draw/intern/draw_manager_data.c
index 264f0ef6d33..b955bd3000c 100644
--- a/source/blender/draw/intern/draw_manager_data.c
+++ b/source/blender/draw/intern/draw_manager_data.c
@@ -28,6 +28,7 @@
 #include "BKE_curve.h"
 #include "BKE_global.h"
 #include "BKE_mesh.h"
+#include "BKE_object.h"
 #include "BKE_paint.h"
 #include "BKE_pbvh.h"
 
@@ -246,16 +247,13 @@ static void drw_call_calc_orco(ID *ob_data, float (*r_orcofacs)[3])
 	}
 }
 
-static DRWCallState *drw_call_state_create(DRWShadingGroup *shgroup, float (*obmat)[4], ID *ob_data)
+static DRWCallState *drw_call_state_create(DRWShadingGroup *shgroup, float (*obmat)[4], Object *ob)
 {
 	DRWCallState *state = BLI_mempool_alloc(DST.vmempool->states);
 	state->flag = 0;
 	state->cache_id = 0;
 	state->matflag = shgroup->matflag;
 
-	/* TODO Set culling bsphere IF needed by the DRWPass */
-	state->bsphere.rad = -1.0f;
-
 	/* Matrices */
 	if (obmat != NULL) {
 		copy_m4_m4(state->model, obmat);
@@ -268,19 +266,33 @@ static DRWCallState *drw_call_state_create(DRWShadingGroup *shgroup, float (*obm
 		unit_m4(state->model);
 	}
 
+	if (ob != NULL) {
+		float corner[3];
+		BoundBox *bbox = BKE_object_boundbox_get(ob);
+		/* Get BoundSphere center and radius from the BoundBox. */
+		mid_v3_v3v3(state->bsphere.center, bbox->vec[0], bbox->vec[6]);
+		mul_v3_m4v3(corner, obmat, bbox->vec[0]);
+		mul_m4_v3(obmat, state->bsphere.center);
+		state->bsphere.radius = len_v3v3(state->bsphere.center, corner);
+	}
+	else {
+		/* Bypass test. */
+		state->bsphere.radius = -1.0f;
+	}
+
 	/* Orco factors: We compute this at creation to not have to save the *ob_data */
 	if ((state->matflag & DRW_CALL_ORCOTEXFAC) != 0) {
-		drw_call_calc_orco(ob_data, state->orcotexfac);
+		drw_call_calc_orco(ob->data, state->orcotexfac);
 		state->matflag &= ~DRW_CALL_ORCOTEXFAC;
 	}
 
 	return state;
 }
 
-static DRWCallState *drw_call_state_object(DRWShadingGroup *shgroup, float (*obmat)[4], ID *ob_data)
+static DRWCallState *drw_call_state_object(DRWShadingGroup *shgroup, float (*obmat)[4], Object *ob)
 {
 	if (DST.ob_state == NULL) {
-		DST.ob_state = drw_call_state_create(shgroup, obmat, ob_data);
+		DST.ob_state = drw_call_state_create(shgroup, obmat, ob);
 	}
 	else {
 		/* If the DRWCallState is reused, add necessary matrices. */
@@ -298,10 +310,10 @@ void DRW_shgroup_call_add(DRWShadingGroup *shgroup, Gwn_Batch *geom, float (*obm
 	DRWCall *call = BLI_mempool_alloc(DST.vmempool->calls);
 	call->state = drw_call_state_create(shgroup, obmat, NULL);
 	call->type = DRW_CALL_SINGLE;
+	call->single.geometry = geom;
 #ifdef USE_GPU_SELECT
 	call->select_id = DST.select_id;
 #endif
-	call->single.geometry = geom;
 
 	BLI_LINKS_APPEND(&shgroup->calls, call);
 }
@@ -313,12 +325,12 @@ void DRW_shgroup_call_object_add(DRWShadingGroup *shgroup, Gwn_Batch *geom, Obje
 	BLI_assert(shgroup->type == DRW_SHG_NORMAL);
 
 	DRWCall *call = BLI_mempool_alloc(DST.vmempool->calls);
-	call->state = drw_call_state_object(shgroup, ob->obmat, ob->data);
+	call->state = drw_call_state_object(shgroup, ob->obmat, ob);
 	call->type = DRW_CALL_SINGLE;
+	call->single.geometry = geom;
 #ifdef USE_GPU_SELECT
 	call->select_id = DST.select_id;
 #endif
-	call->single.geometry = geom;
 
 	BLI_LINKS_APPEND(&shgroup->calls, call);
 }
@@ -334,11 +346,11 @@ void DRW_shgroup_call_generate_add(
 	DRWCall *call = BLI_mempool_alloc(DST.vmempool->calls);
 	call->state = drw_call_state_create(shgroup, obmat, NULL);
 	call->type = DRW_CALL_GENERATE;
+	call->generate.geometry_fn = geometry_fn;
+	call->generate.user_data = user_data;
 #ifdef USE_GPU_SELECT
 	call->select_id = DST.select_id;
 #endif
-	call->generate.geometry_fn = geometry_fn;
-	call->generate.user_data = user_data;
 
 	BLI_LINKS_APPEND(&shgroup->calls, call);
 }
diff --git a/source/blender/draw/intern/draw_manager_exec.c b/source/blender/draw/intern/draw_manager_exec.c
index da753efd4d3..5538ed327a4 100644
--- a/source/blender/draw/intern/draw_manager_exec.c
+++ b/source/blender/draw/intern/draw_manager_exec.c
@@ -30,6 +30,7 @@
 #include "BIF_glutil.h"
 
 #include "BKE_global.h"
+#include "BKE_object.h"
 
 #include "GPU_draw.h"
 #include "GPU_extensions.h"
@@ -364,6 +365,158 @@ void DRW_state_clip_planes_reset(void)
 
 /* -------------------------------------------------------------------- */
 
+/** \name Clipping (DRW_clipping)
+ * \{ */
+
+static void draw_clipping_setup_from_view(void)
+
+{
+	if (DST.clipping.updated)
+		return;
+
+	float (*viewprojinv)[4] = DST.view_data.mat[DRW_MAT_PERSINV];
+	float (*viewinv)[4] = DST.view_data.mat[DRW_MAT_VIEWINV];
+	float (*projmat)[4] = DST.view_data.mat[DRW_MAT_WIN];
+	float (*projinv)[4] = DST.view_data.mat[DRW_MAT_WININV];
+	BoundSphere *bsphere = &DST.clipping.frustum_bsphere;
+
+	/* Extract Clipping Planes */
+	BoundBox bbox;
+	BKE_boundbox_init_from_minmax(&bbox, (const float[3]){-1.0f, -1.0f, -1.0f}, (const float[3]){1.0f, 1.0f, 1.0f});
+
+	/* Extract the 8 corners (world space). */
+	for (int i = 0; i < 8; i++) {
+		mul_project_m4_v3(viewprojinv, bbox.vec[i]);
+	}
+
+	/* Compute clip planes using the world space frustum corners. */
+	for (int p = 0; p < 6; p++) {
+		int q, r;
+		switch (p) {
+			case 0:  q=1; r=2; break;
+			case 1:  q=0; r=5; break;
+			case 2:  q=1; r=5; break;
+			case 3:  q=2; r=6; break;
+			case 4:  q=0; r=3; break;
+			default: q=4; r=7; break;
+		}
+
+		normal_tri_v3(DST.clipping.frustum_planes[p], bbox.vec[p], bbox.vec[q], bbox.vec[r]);
+		DST.clipping.frustum_planes[p][3] = -dot_v3v3(DST.clipping.frustum_planes[p], bbox.vec[p]);
+	}
+
+	/* Extract Bounding Sphere */
+	/**
+	 * Compute bounding sphere for the general case and not only symmetric frustum:
+	 * We put the sphere center on the line that goes from origin to the center of the far clipping plane.
+	 * This is the optimal position if the frustum is symmetric or very asymmetric and probably close
+	 * to optimal for the general case. The sphere center position is computed so that the distance to
+	 * the near and far extreme frustum points are equal.
+	 **/
+	if (projmat[3][3] == 0.0f) {
+		/* Perspective */
+		/* Detect which of the corner of the far clipping plane is the farthest to the origin */
+		float nfar[4];       /* most extreme far point in NDC space */
+		float farxy[2];      /* farpoint projection onto the near plane */
+		float farpoint[3] = {0.0f}; /* most extreme far point in camera coordinate */
+		float nearpoint[3];  /* most extreme near point in camera coordinate */
+		float farcenter[3] = {0.0f}; /* center of far cliping plane in camera coordinate */
+		float F = -1.0f, N;  /* square distance of far and near point to origin */
+		float f, n;          /* distance of far and near point to z axis. f is always > 0 but n can be < 0 */
+		float e, s;          /* far and near clipping distance (<0) */
+		float c;             /* slope of center line = distance of far clipping center to z axis / far clipping distance */
+		float z;             /* projection of sphere center on z axis (<0) */
+
+		/* Find farthest corner and center of far clip plane. */
+		float corner[3] = {1.0f, 1.0f, 1.0f}; /* in clip space */
+		for (int i = 0; i < 4; i++) {
+			float point[3];
+			mul_v3_project_m4_v3(point, projinv, corner);
+			float len = len_squared_v3(point);
+			if (len > F) {
+				copy_v3_v3(nfar, corner);
+				copy_v3_v3(farpoint, point);
+				F = len;
+			}
+			add_v3_v3(farcenter, point);
+			/* rotate by 90 degree to walk through the 4 points of the far clip plane */
+			float tmp = corner[0];
+			corner[0] = -corner[1];
+			corner[1] = tmp;
+		}
+
+		/* the far center is the average of the far clipping points */
+		mul_v3_fl(farcenter, 0.25f);
+		/* the extreme near point is the opposite point on the near clipping plane */
+		copy_v3_fl3(nfar, -nfar[0], -nfar[1], -1.0f);
+		mul_v3_project_m4_v3(nearpoint, projinv, nfar);
+		/* this is a frustum projection */
+		N = len_squared_v3(nearpoint);
+		e = farpoint[2];
+		s = nearpoint[2];
+		/* distance to view Z axis */
+		f = len_v2(nearpoint);
+		/* get corresponding point on the near plane */
+		mul_v2_v2fl(farxy, farpoint, s/e);
+		/* this formula preserve the sign of n */
+		sub_v2_v2(nearpoint, farxy);
+		n = f * s / e - len_v2(nearpoint);
+		c = len_v2(farcenter) / e;
+		/* the big formula, it simplifies to (F-N)/(2(e-s)) for the symmetric case */
+		z = (F-N) / (2.0f * (e-s + c*(f-n)));
+
+		bsphere->center[0] = farcenter[0] * z/e;
+		bsphere->center[1] = farcenter[1] * z/e;
+		bsphere->center[2] = z;
+		bsphere->radius = len_v3v3(bsphere->center, farpoint);
+	}
+	else {
+		/* Orthographic */
+		/* The most extreme points on the near and far plane. (normalized device coords) */
+		float nearpoint[3] = {-1.0f, -1.0f, -1.0f};
+		float farpoint[3] =  { 1.0f,  1.0f,  1.0f};
+
+		mul_project_m4_v3(projinv, nearpoint);
+		mul_project_m4_v3(projinv, farpoint);
+
+		/* just use median point */
+		mid_v3_v3v3(bsphere->center, farpoint, nearpoint);
+		bsphere->radius = len_v3v3(bsphere->center, farpoint);
+	}
+
+	/* Transform to world space. */
+	mul_m4_v3(viewinv, bsphere->center);
+
+}
+
+/* Return True if the given BoundSphere intersect the current view frustum */
+static bool draw_culling_sphere_test(BoundSphere *bsphere)
+{
+	/* Bypass test if radius is negative. */
+	if (bsphere->radius < 0.0f)
+		return true;
+
+	/* Do a rough test first: Sphere VS Sphere intersect. */
+	BoundSphere *frustum_bsphere = &DST.clipping.frustum_bsphere;
+	float center_dist = len_squared_v3v3(bsphere->center, frustum_bsphere->center);
+	if (center_dist > SQUARE(bsphere->radius + frustum_bsphere->radius))
+		return false;
+
+	/* Test against the 6 frustum planes. */
+	for (int p = 0; p < 6; p++) {
+		float dist = plane_point_side_v3(DST.clipping.frustum_planes[p], bsphere->center);
+		if (dist < -bsphere->radius) {
+			return false;
+		}
+	}
+
+	return true;
+}
+
+/** \} */
+
+/* -------------------------------------------------------------------- */
+
 /** \name Draw (DRW_draw)
  * \{ */
 
@@ -376,6 +529,14 @@ static void draw_matrices_model_prepare(DRWCallState *st)
 		st->cache_id = DST.state_cache_id;
 	}
 
+	if (draw_culling_sphere_test(&st->bsphere)) {
+		st->flag &= ~DRW_CALL_CULLED;
+	}
+	else {
+		st->flag |= DRW_CALL_CULLED;
+		return; /* No need to go further the call will not be used. */
+	}
+
 	/* Order matters */
 	if (st->matflag & (DRW_CALL_MODELVIEW | DRW_CALL_MODELVIEWINVERSE |
 	                  DRW_CALL_NORMALVIEW | DRW_CALL_EYEVEC))
@@ -418,8 +579,6 @@ static void draw_geometry_prepare(DRWShadingGroup *shgroup, DRWCallState *state)
 {
 	/* step 1 : bind object dependent matrices */
 	if (state != NULL) {
-		/* OPTI/IDEA(clem): Do this preparation in another thread. */
-		draw_matrices_model_prepare(state);
 		GPU_shader_uniform_vector(shgroup->shader, shgroup->model, 16, 1, (float *)state->model);
 		GPU_shader_uniform_vector(shgroup->shader, shgroup->modelinverse, 16, 1, (float *)state->modelinverse);
 		GPU_shader_uniform_vector(shgroup->shader, shgroup->modelview, 16, 1, (float *)state->modelview);
@@ -675,6 +834,10 @@ static void draw_shgroup(DRWShadingGroup *shgroup, DRWState pass_state)
 	else {
 		bool prev_neg_scale = false;
 		for (DRWCall *call = shgroup->calls.first; call; call = call->next) {
+
+			/* OPTI/IDEA(clem): Do this preparation in another thread. */
+			draw_matrices_model_prepare(call->state);
+
 			if ((call->state->flag & DRW_CALL_CULLED) != 0)
 				continue;
 
@@ -711,6 +874,7 @@ static void drw_draw_pass_ex(DRWPass *pass, DRWShadingGroup *start_group, DRWSha
 	if (DST.dirty_mat) {
 		DST.state_cache_id++;
 		DST.dirty_mat = false;
+
 		/* Catch integer wrap around. */
 		if (UNLIKELY(DST.state_cache_id == 0)) {
 			DST.state_cache_id = 1;
@@ -723,9 +887,14 @@ static void drw_draw_pass_ex(DRWPass *pass, DRWShadingGroup *start_group, DRWSha
 				state->cache_id = 0;
 			}
 		}
+
+		DST.clipping.updated = false;
+
 		/* TODO dispatch threads to compute matrices/culling */
 	}
 
+	draw_clipping_setup_from_view();
+
 	BLI_assert(DST.buffer_finish_called && "DRW_render_instance_buffer_finish had not been called before drawing");
 
 	drw_state_set(pass->state);