suvm_primitive.c (13134B)
1 /* Copyright (C) 2020-2023 |Méso|Star> (contact@meso-star.com) 2 * 3 * This program is free software: you can redistribute it and/or modify 4 * it under the terms of the GNU General Public License as published by 5 * the Free Software Foundation, either version 3 of the License, or 6 * (at your option) any later version. 7 * 8 * This program is distributed in the hope that it will be useful, 9 * but WITHOUT ANY WARRANTY; without even the implied warranty of 10 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 * GNU General Public License for more details. 12 * 13 * You should have received a copy of the GNU General Public License 14 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 15 16 #include "suvm.h" 17 #include "suvm_volume.h" 18 19 #include <rsys/float2.h> 20 #include <rsys/float3.h> 21 22 /******************************************************************************* 23 * Exported functions 24 ******************************************************************************/ 25 res_T 26 suvm_primitive_setup_polyhedron 27 (const struct suvm_primitive* prim, 28 struct suvm_polyhedron* poly) 29 { 30 if(!prim || !poly) return RES_BAD_ARG; 31 tetrahedron_setup(prim->volume__, NULL, NULL, prim->iprim, poly); 32 return RES_OK; 33 } 34 35 /* Define if an axis aligned bounding box and a tetrahedron are intersecting. 36 * Its implementation follows the algorithms proposed by N. Greene in 37 * "Detecting intersection of a rectangular solid and a convex polyhedron" 38 * (Graphics Gems IV, 1994, p74--82) */ 39 enum suvm_intersection_type 40 suvm_polyhedron_intersect_aabb 41 (const struct suvm_polyhedron* tetra, 42 const float low[3], 43 const float upp[3]) 44 { 45 enum { X, Y, Z, AXES_COUNT }; /* Syntactic sugar */ 46 float intersect_aabb_low[AXES_COUNT]; 47 float intersect_aabb_upp[AXES_COUNT]; 48 int nplanes_including_aabb; 49 int i; 50 51 ASSERT(tetra && low && upp); 52 ASSERT(low[0] < upp[0]); 53 ASSERT(low[1] < upp[1]); 54 ASSERT(low[2] < upp[2]); 55 56 /* Check the intersection between the AABB and the tetra AABB */ 57 FOR_EACH(i, 0, AXES_COUNT) { 58 intersect_aabb_low[i] = MMAX(low[i], tetra->lower[i]); 59 intersect_aabb_upp[i] = MMIN(upp[i], tetra->upper[i]); 60 if(intersect_aabb_low[i] > intersect_aabb_upp[i]) /* Do not intersect */ 61 return SUVM_INTERSECT_NONE; 62 } 63 64 /* Check if the tetrahedron is included into the aabb */ 65 if(intersect_aabb_low[X] == tetra->lower[X] 66 && intersect_aabb_low[Y] == tetra->lower[Y] 67 && intersect_aabb_low[Z] == tetra->lower[Z] 68 && intersect_aabb_upp[X] == tetra->upper[X] 69 && intersect_aabb_upp[Y] == tetra->upper[Y] 70 && intersect_aabb_upp[Z] == tetra->upper[Z]) { 71 return SUVM_INTERSECT_IS_INCLUDED; 72 } 73 74 /* #planes that totally includes the aabb on its positive side */ 75 nplanes_including_aabb = 0; 76 77 /* Check the tetrahedron planes against the aabb */ 78 FOR_EACH(i, 0, 4) { 79 float p[3], n[3]; 80 81 /* Define the aabb vertices that is the farthest in the positive/negative 82 * direction of the face normal. */ 83 f3_set(p, upp); 84 f3_set(n, low); 85 if(tetra->N[i][X] < 0) SWAP(float, p[X], n[X]); 86 if(tetra->N[i][Y] < 0) SWAP(float, p[Y], n[Y]); 87 if(tetra->N[i][Z] < 0) SWAP(float, p[Z], n[Z]); 88 89 /* Check that the box is totally outside the plane, To do this, check that 90 * 'p', aka the farthest aabb vertex in the positive direction of the plane 91 * normal, is not in the positive half space. In this case, the aabb is 92 * entirely outside the tetrahedron */ 93 if((f3_dot(tetra->N[i], p) + tetra->D[i]) < 0) { 94 return SUVM_INTERSECT_NONE; 95 } 96 97 /* Check if the box is totally inside the given plane. To do this, check 98 * that 'n', aka the farthest aabb vertex in the negative direction of the 99 * plane normal, is inside the positive half space */ 100 if((f3_dot(tetra->N[i], n) + tetra->D[i]) >= 0) { 101 /* Register this plane as a plane including the aabb */ 102 nplanes_including_aabb += 1; 103 } 104 } 105 106 /* Check if the aabb is entirely included into the tetrahedron */ 107 if(nplanes_including_aabb == 4) { 108 return SUVM_INTERSECT_INCLUDE; 109 } 110 111 /* For each silhouette edge in each projection plane, check if it totally 112 * exclude the projected aabb */ 113 FOR_EACH(i, 0, AXES_COUNT) { 114 /* Compute the remaining axes. 115 * On 1st iteration 'j' and 'k' define the YZ plane 116 * On 2nd iteration 'j' and 'k' define the ZX plane 117 * On 3rd ietration 'j' and 'k' define the XY plane */ 118 const int j = (i + 1) % AXES_COUNT; 119 const int k = (j + 1) % AXES_COUNT; 120 121 int iedge; 122 FOR_EACH(iedge, 0, tetra->nEp[i]) { 123 float p[2]; 124 125 /* Define the aabb vertices that is the farthest in the positive/negative 126 * direction of the normal of the silhouette edge. */ 127 p[0] = tetra->Ep[i][iedge][0] > 0 ? upp[j] : low[j]; 128 p[1] = tetra->Ep[i][iedge][1] > 0 ? upp[k] : low[k]; 129 130 /* Check that the projected aabb along the 'i'th axis is totally outside 131 * the current silhouette edge. That means that the aabb is entirely 132 * outside the tetrahedron and thus that they do not intersect */ 133 if((f2_dot(tetra->Ep[i][iedge], p) + tetra->Ep[i][iedge][2]) < 0) { 134 return SUVM_INTERSECT_NONE; 135 } 136 } 137 } 138 139 /* The tetrahedron and the aabb are partially intersecting */ 140 return SUVM_INTERSECT_PARTIAL; 141 } 142 143 144 /******************************************************************************* 145 * Local functions 146 ******************************************************************************/ 147 void 148 tetrahedron_setup 149 (const struct suvm_volume* vol, 150 const float lower[3], 151 const float upper[3], 152 const size_t itetra, 153 struct suvm_polyhedron* tetra) 154 { 155 enum { X, Y, Z, AXES_COUNT }; /* Syntactic sugar */ 156 const float* low = lower; 157 const float* upp = upper; 158 float low__[3]; 159 float upp__[3]; 160 float center[AXES_COUNT]; /* Center of the tetrahedron */ 161 float e[6][AXES_COUNT]; 162 float (*v)[AXES_COUNT]; 163 float (*N)[AXES_COUNT]; 164 float (*Ep)[4][3]; 165 int i; 166 ASSERT(vol && tetra); 167 168 /* Fetch tetrahedron vertices */ 169 volume_primitive_get_vertex_position(vol, itetra, 0, tetra->v[0]); 170 volume_primitive_get_vertex_position(vol, itetra, 1, tetra->v[1]); 171 volume_primitive_get_vertex_position(vol, itetra, 2, tetra->v[2]); 172 volume_primitive_get_vertex_position(vol, itetra, 3, tetra->v[3]); 173 v = tetra->v; 174 175 /* Compute the center of the tetrahedron */ 176 center[X] = (v[0][X] + v[1][X] + v[2][X] + v[3][X]) * 0.25f; 177 center[Y] = (v[0][Y] + v[1][Y] + v[2][Y] + v[3][Y]) * 0.25f; 178 center[Z] = (v[0][Z] + v[1][Z] + v[2][Z] + v[3][Z]) * 0.25f; 179 180 /* Define the primitive AABB if necessary */ 181 if(!low) { 182 low__[0] = MMIN(MMIN(v[0][X], v[1][X]), MMIN(v[2][X], v[3][X])); 183 low__[1] = MMIN(MMIN(v[0][Y], v[1][Y]), MMIN(v[2][Y], v[3][Y])); 184 low__[2] = MMIN(MMIN(v[0][Z], v[1][Z]), MMIN(v[2][Z], v[3][Z])); 185 low = low__; 186 } 187 if(!upp) { 188 upp__[0] = MMAX(MMAX(v[0][X], v[1][X]), MMAX(v[2][X], v[3][X])); 189 upp__[1] = MMAX(MMAX(v[0][Y], v[1][Y]), MMAX(v[2][Y], v[3][Y])); 190 upp__[2] = MMAX(MMAX(v[0][Z], v[1][Z]), MMAX(v[2][Z], v[3][Z])); 191 upp = upp__; 192 } 193 194 #ifndef NDEBUG 195 /* Check argument consistency */ 196 if(lower) { 197 ASSERT(MMIN(MMIN(v[0][X], v[1][X]), MMIN(v[2][X], v[3][X])) == lower[X]); 198 ASSERT(MMIN(MMIN(v[0][Y], v[1][Y]), MMIN(v[2][Y], v[3][Y])) == lower[Y]); 199 ASSERT(MMIN(MMIN(v[0][Z], v[1][Z]), MMIN(v[2][Z], v[3][Z])) == lower[Z]); 200 } 201 if(upper) { 202 ASSERT(MMAX(MMAX(v[0][X], v[1][X]), MMAX(v[2][X], v[3][X])) == upper[X]); 203 ASSERT(MMAX(MMAX(v[0][Y], v[1][Y]), MMAX(v[2][Y], v[3][Y])) == upper[Y]); 204 ASSERT(MMAX(MMAX(v[0][Z], v[1][Z]), MMAX(v[2][Z], v[3][Z])) == upper[Z]); 205 } 206 #endif 207 208 /* Setup tetrahedron AABB */ 209 f3_set(tetra->lower, low); 210 f3_set(tetra->upper, upp); 211 212 /* Fetch tetrahedron normals */ 213 volume_primitive_get_facet_normal(vol, itetra, 0, tetra->N[0]); 214 volume_primitive_get_facet_normal(vol, itetra, 1, tetra->N[1]); 215 volume_primitive_get_facet_normal(vol, itetra, 2, tetra->N[2]); 216 volume_primitive_get_facet_normal(vol, itetra, 3, tetra->N[3]); 217 N = tetra->N; 218 219 /* Compute the slope of the planes */ 220 tetra->D[0] = -f3_dot(tetra->N[0], v[0]); 221 tetra->D[1] = -f3_dot(tetra->N[1], v[1]); 222 tetra->D[2] = -f3_dot(tetra->N[2], v[2]); 223 tetra->D[3] = -f3_dot(tetra->N[3], v[3]); 224 225 /* Compute tetrahedron edges */ 226 f3_sub(e[0], v[1], v[0]); 227 f3_sub(e[1], v[2], v[1]); 228 f3_sub(e[2], v[0], v[2]); 229 f3_sub(e[3], v[0], v[3]); 230 f3_sub(e[4], v[1], v[3]); 231 f3_sub(e[5], v[2], v[3]); 232 233 /* Detect the silhouette edges of the tetrahedron once projected in the YZ, 234 * ZX and XY planes, and compute their 2D equation in the projected plane. The 235 * variable 'i' defines the projection axis which is successively X (YZ 236 * plane), Y (ZX plane) and Z (XY plane). The axes of the corresponding 2D 237 * repairs are defined by the 'j' and 'k' variables. */ 238 Ep = tetra->Ep; 239 FOR_EACH(i, 0, AXES_COUNT) { 240 float c[2]; /* Projected tetrahedron center */ 241 242 /* On 1st iteration 'j' and 'k' define the YZ plane 243 * On 2nd iteration 'j' and 'k' define the ZX plane 244 * On 3rd ietration 'j' and 'k' define the XY plane */ 245 const int j = (i + 1) % AXES_COUNT; 246 const int k = (j + 1) % AXES_COUNT; 247 248 /* Register the number of detected silhouette edges */ 249 int n = 0; 250 251 int iedge; 252 253 /* Project the tetrahedron center */ 254 c[0] = center[j]; 255 c[1] = center[k]; 256 257 /* To detect the silhouette edges, check the sign of the normals of two 258 * adjacent facets for the coordinate of the projection axis. If the signs 259 * are the same, the facets look at the same direction regarding the 260 * projection axis. The other case means that one facet points toward the 261 * projection axis whole the other one points to the opposite direction: 262 * this is the definition of a silhouette edge. 263 * 264 * Once detected, we compute the equation of the silhouette edge: 265 * 266 * Ax + By + C = 0 267 * 268 * with A and B the coordinates of the edge normals and C the slope of the 269 * edge. Computing the normal is actually as simple as swizzling the 270 * coordinates of the edges projected in 2D and sign tuning. Let the 271 * projection axis X (i==X) and thus the YZ projection plane (j==Y && 272 * k==Z). Assuming that the first edge 'e[0]' is a silhouette edge, ie the 273 * edge between the facet0 and the facet1 (refers to the suvm_volume.h for 274 * the geometric layout of a tetrahedron) once projected in the YZ plane 275 * the edge is {0, e[0][Y], e[0][Z]}. Its normal N can is defined as the 276 * cross product between this projected edge ands the projection axis: 277 * 278 * | 0 | | 1 | | 0 | 279 * N = | e[0][Y] | X | 0 | = | e[0][Z] | 280 * | e[0][Z] | | 0 | |-e[0][Y] | 281 * 282 * In the projection plane, the _unormalized_ normal of the {e[0][Y], 283 * e[0][Z]} edge is thus {e[0][Z], -e[0][Y]}. More generally, the normal of 284 * an edge 'I' {e[I][j], e[I][k]} projected along the 'i' axis in the 'jk' 285 * plane is {e[I][k],-e[I]j]}. Anyway, one has to revert the normal 286 * orientation to ensure a specific convention wrt the tetrahedron 287 * footprint. In our case we ensure that the normal points toward the 288 * footprint. Finalle the edge slope 'C' can be computed as: 289 * | A | 290 * C =- | B | . | Pj, Pk | 291 * 292 * with {Pj, Pk} a point belonging to the edge as for instance one of its 293 * vertices. */ 294 if(signf(N[0][i]) != signf(N[1][i])) { /* The edge 0 is silhouette */ 295 f2_normalize(Ep[i][n], f2(Ep[i][n], e[0][k], -e[0][j])); 296 Ep[i][n][2] = -(Ep[i][n][0]*v[0][j] + Ep[i][n][1]*v[0][k]); 297 ++n; 298 } 299 if(signf(N[0][i]) != signf(N[2][i])) { /* The edge 1 is silhouette */ 300 f2_normalize(Ep[i][n], f2(Ep[i][n], e[1][k], -e[1][j])); 301 Ep[i][n][2] = -(Ep[i][n][0]*v[1][j] + Ep[i][n][1]*v[1][k]); 302 ++n; 303 } 304 if(signf(N[0][i]) != signf(N[3][i])) { /* The edge 2 is silhouette */ 305 f2_normalize(Ep[i][n], f2(Ep[i][n], e[2][k], -e[2][j])); 306 Ep[i][n][2] = -(Ep[i][n][0]*v[2][j] + Ep[i][n][1]*v[2][k]); 307 ++n; 308 } 309 if(signf(N[1][i]) != signf(N[3][i])) { /* The edge 3 is silhouette */ 310 f2_normalize(Ep[i][n], f2(Ep[i][n], e[3][k], -e[3][j])); 311 Ep[i][n][2] = -(Ep[i][n][0]*v[0][j] + Ep[i][n][1]*v[0][k]); 312 ++n; 313 } 314 if(signf(N[1][i]) != signf(N[2][i])) { /* The edge 4 is silhouette */ 315 f2_normalize(Ep[i][n], f2(Ep[i][n], e[4][k], -e[4][j])); 316 Ep[i][n][2] = -(Ep[i][n][0]*v[1][j] + Ep[i][n][1]*v[1][k]); 317 ++n; 318 } 319 if(signf(N[2][i]) != signf(N[3][i])) { /* The edge 5 is silhouette */ 320 f2_normalize(Ep[i][n], f2(Ep[i][n], e[5][k], -e[5][j])); 321 Ep[i][n][2] = -(Ep[i][n][0]*v[2][j] + Ep[i][n][1]*v[2][k]); 322 ++n; 323 } 324 /* A tetrahedron can have 3 or 4 silhouette edges once projected in 2D */ 325 ASSERT(n == 3 || n == 4); 326 tetra->nEp[i] = n; /* Register the #silouhette edges for this project axis */ 327 328 /* Ensure that the edge normals point toward the tetrahedron center */ 329 FOR_EACH(iedge, 0, n) { 330 if(f2_dot(Ep[i][iedge], c) + Ep[i][iedge][2] < 0) 331 f3_minus(Ep[i][iedge], Ep[i][iedge]); 332 } 333 } 334 } 335