intern/dualcon/intern/Projections.cpp

   1 /*
   2  * ***** BEGIN GPL LICENSE BLOCK *****
   3  *
   4  * This program is free software; you can redistribute it and/or
   5  * modify it under the terms of the GNU General Public License
   6  * as published by the Free Software Foundation; either version 2
   7  * of the License, or (at your option) any later version.
   8  *
   9  * This program is distributed in the hope that it will be useful,
  10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12  * GNU General Public License for more details.
  13  *
  14  * You should have received a copy of the GNU General Public License
  15  * along with this program; if not, write to the Free Software Foundation,
  16  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  17  *
  18  * Contributor(s): Tao Ju
  19  *
  20  * ***** END GPL LICENSE BLOCK *****
  21  */
  22
  23 #include <math.h>
  24 #include "Projections.h"
  25
  26 const int vertmap[8][3] = {
  27         {0, 0, 0},
  28         {0, 0, 1},
  29         {0, 1, 0},
  30         {0, 1, 1},
  31         {1, 0, 0},
  32         {1, 0, 1},
  33         {1, 1, 0},
  34         {1, 1, 1}
  35 };
  36
  37 const int centmap[3][3][3][2] = {
  38         {{{0, 0}, {0, 1}, {1, 1}},
  39          {{0, 2}, {0, 3}, {1, 3}},
  40          {{2, 2}, {2, 3}, {3, 3}}},
  41
  42         {{{0, 4}, {0, 5}, {1, 5}},
  43          {{0, 6}, {0, 7}, {1, 7}},
  44          {{2, 6}, {2, 7}, {3, 7}}},
  45
  46         {{{4, 4}, {4, 5}, {5, 5}},
  47          {{4, 6}, {4, 7}, {5, 7}},
  48          {{6, 6}, {6, 7}, {7, 7}}}
  49 };
  50
  51 const int edgemap[12][2] = {
  52         {0, 4},
  53         {1, 5},
  54         {2, 6},
  55         {3, 7},
  56         {0, 2},
  57         {1, 3},
  58         {4, 6},
  59         {5, 7},
  60         {0, 1},
  61         {2, 3},
  62         {4, 5},
  63         {6, 7}
  64 };
  65
  66 const int facemap[6][4] = {
  67         {0, 1, 2, 3},
  68         {4, 5, 6, 7},
  69         {0, 1, 4, 5},
  70         {2, 3, 6, 7},
  71         {0, 2, 4, 6},
  72         {1, 3, 5, 7}
  73 };
  74
  75 /**
  76  * Method to perform cross-product
  77  */
  78 static void crossProduct(int64_t res[3], const int64_t a[3], const int64_t b[3])
  79 {
  80         res[0] = a[1] * b[2] - a[2] * b[1];
  81         res[1] = a[2] * b[0] - a[0] * b[2];
  82         res[2] = a[0] * b[1] - a[1] * b[0];
  83 }
  84
  85 static void crossProduct(double res[3], const double a[3], const double b[3])
  86 {
  87         res[0] = a[1] * b[2] - a[2] * b[1];
  88         res[1] = a[2] * b[0] - a[0] * b[2];
  89         res[2] = a[0] * b[1] - a[1] * b[0];
  90 }
  91
  92 /**
  93  * Method to perform dot product
  94  */
  95 int64_t dotProduct(const int64_t a[3], const int64_t b[3])
  96 {
  97         return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
  98 }
  99
 100 void normalize(double a[3])
 101 {
 102         double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
 103         if (mag > 0) {
 104                 mag = sqrt(mag);
 105                 a[0] /= mag;
 106                 a[1] /= mag;
 107                 a[2] /= mag;
 108         }
 109 }
 110
 111 /* Create projection axes for cube+triangle intersection testing.
 112  *    0, 1, 2: cube face normals
 113  *
 114  *          3: triangle normal
 115  *
 116  *    4, 5, 6,
 117  *    7, 8, 9,
 118  * 10, 11, 12: cross of each triangle edge vector with each cube
 119  *             face normal
 120  */
 121 static void create_projection_axes(int64_t axes[NUM_AXES][3], const int64_t tri[3][3])
 122 {
 123         /* Cube face normals */
 124         axes[0][0] = 1;
 125         axes[0][1] = 0;
 126         axes[0][2] = 0;
 127         axes[1][0] = 0;
 128         axes[1][1] = 1;
 129         axes[1][2] = 0;
 130         axes[2][0] = 0;
 131         axes[2][1] = 0;
 132         axes[2][2] = 1;
 133
 134         /* Get triangle edge vectors */
 135         int64_t tri_edges[3][3];
 136         for (int i = 0; i < 3; i++) {
 137                 for (int j = 0; j < 3; j++)
 138                         tri_edges[i][j] = tri[(i + 1) % 3][j] - tri[i][j];
 139         }
 140
 141         /* Triangle normal */
 142         crossProduct(axes[3], tri_edges[0], tri_edges[1]);
 143
 144         // Face edges and triangle edges
 145         int ct = 4;
 146         for (int i = 0; i < 3; i++) {
 147                 for (int j = 0; j < 3; j++) {
 148                         crossProduct(axes[ct], axes[j], tri_edges[i]);
 149                         ct++;
 150                 }
 151         }
 152 }
 153
 154 /**
 155  * Construction from a cube (axes aligned) and triangle
 156  */
 157 CubeTriangleIsect::CubeTriangleIsect(int64_t cube[2][3], int64_t tri[3][3], int64_t error, int triind)
 158 {
 159         int i;
 160         inherit = new TriangleProjection;
 161         inherit->index = triind;
 162
 163         int64_t axes[NUM_AXES][3];
 164         create_projection_axes(axes, tri);
 165
 166         /* Normalize face normal and store */
 167         double dedge1[] = {(double)tri[1][0] - (double)tri[0][0],
 168                                            (double)tri[1][1] - (double)tri[0][1],
 169                                            (double)tri[1][2] - (double)tri[0][2]};
 170         double dedge2[] = {(double)tri[2][0] - (double)tri[1][0],
 171                                            (double)tri[2][1] - (double)tri[1][1],
 172                                            (double)tri[2][2] - (double)tri[1][2]};
 173         crossProduct(inherit->norm, dedge1, dedge2);
 174         normalize(inherit->norm);
 175
 176         int64_t cubeedge[3][3];
 177         for (i = 0; i < 3; i++) {
 178                 for (int j = 0; j < 3; j++) {
 179                         cubeedge[i][j] = 0;
 180                 }
 181                 cubeedge[i][i] = cube[1][i] - cube[0][i];
 182         }
 183
 184         /* Project the cube on to each axis */
 185         for (int axis = 0; axis < NUM_AXES; axis++) {
 186                 CubeProjection &cube_proj = cubeProj[axis];
 187
 188                 /* Origin */
 189                 cube_proj.origin = dotProduct(axes[axis], cube[0]);
 190
 191                 /* 3 direction vectors */
 192                 for (i = 0; i < 3; i++)
 193                         cube_proj.edges[i] = dotProduct(axes[axis], cubeedge[i]);
 194
 195                 /* Offsets of 2 ends of cube projection */
 196                 int64_t max = 0;
 197                 int64_t min = 0;
 198                 for (i = 1; i < 8; i++) {
 199                         int64_t proj = (vertmap[i][0] * cube_proj.edges[0] +
 200                                                         vertmap[i][1] * cube_proj.edges[1] +
 201                                                         vertmap[i][2] * cube_proj.edges[2]);
 202                         if (proj > max) {
 203                                 max = proj;
 204                         }
 205                         if (proj < min) {
 206                                 min = proj;
 207                         }
 208                 }
 209                 cube_proj.min = min;
 210                 cube_proj.max = max;
 211
 212         }
 213
 214         /* Project the triangle on to each axis */
 215         for (int axis = 0; axis < NUM_AXES; axis++) {
 216                 const int64_t vts[3] = {dotProduct(axes[axis], tri[0]),
 217                                                                 dotProduct(axes[axis], tri[1]),
 218                                                                 dotProduct(axes[axis], tri[2])};
 219
 220                 // Triangle
 221                 inherit->tri_proj[axis][0] = vts[0];
 222                 inherit->tri_proj[axis][1] = vts[0];
 223                 for (i = 1; i < 3; i++) {
 224                         if (vts[i] < inherit->tri_proj[axis][0])
 225                                 inherit->tri_proj[axis][0] = vts[i];
 226
 227                         if (vts[i] > inherit->tri_proj[axis][1])
 228                                 inherit->tri_proj[axis][1] = vts[i];
 229                 }
 230         }
 231 }
 232
 233 /**
 234  * Construction
 235  * from a parent CubeTriangleIsect object and the index of the children
 236  */
 237 CubeTriangleIsect::CubeTriangleIsect(CubeTriangleIsect *parent)
 238 {
 239         // Copy inheritable projections
 240         this->inherit = parent->inherit;
 241
 242         // Shrink cube projections
 243         for (int i = 0; i < NUM_AXES; i++) {
 244                 cubeProj[i].origin = parent->cubeProj[i].origin;
 245
 246                 for (int j = 0; j < 3; j++)
 247                         cubeProj[i].edges[j] = parent->cubeProj[i].edges[j] >> 1;
 248
 249                 cubeProj[i].min = parent->cubeProj[i].min >> 1;
 250                 cubeProj[i].max = parent->cubeProj[i].max >> 1;
 251         }
 252 }
 253
 254 unsigned char CubeTriangleIsect::getBoxMask( )
 255 {
 256         int i, j, k;
 257         int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
 258         unsigned char boxmask = 0;
 259         int64_t child_len = cubeProj[0].edges[0] >> 1;
 260
 261         for (i = 0; i < 3; i++) {
 262                 int64_t mid = cubeProj[i].origin + child_len;
 263
 264                 // Check bounding box
 265                 if (mid >= inherit->tri_proj[i][0]) {
 266                         bmask[i][0] = 1;
 267                 }
 268                 if (mid < inherit->tri_proj[i][1]) {
 269                         bmask[i][1] = 1;
 270                 }
 271
 272         }
 273
 274         // Fill in masks
 275         int ct = 0;
 276         for (i = 0; i < 2; i++) {
 277                 for (j = 0; j < 2; j++) {
 278                         for (k = 0; k < 2; k++) {
 279                                 boxmask |= ( (bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
 280                                 ct++;
 281                         }
 282                 }
 283         }
 284
 285         // Return bounding box masks
 286         return boxmask;
 287 }
 288
 289
 290 /**
 291  * Shifting a cube to a new origin
 292  */
 293 void CubeTriangleIsect::shift(int off[3])
 294 {
 295         for (int i = 0; i < NUM_AXES; i++) {
 296                 cubeProj[i].origin += (off[0] * cubeProj[i].edges[0] +
 297                                                            off[1] * cubeProj[i].edges[1] +
 298                                                            off[2] * cubeProj[i].edges[2]);
 299         }
 300 }
 301
 302 /**
 303  * Method to test intersection of the triangle and the cube
 304  */
 305 int CubeTriangleIsect::isIntersecting() const
 306 {
 307         for (int i = 0; i < NUM_AXES; i++) {
 308                 /*
 309                   int64_t proj0 = cubeProj[i][0] +
 310                   vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
 311                   vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
 312                   vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
 313                   int64_t proj1 = cubeProj[i][0] +
 314                   vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
 315                   vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
 316                   vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
 317                 */
 318
 319                 int64_t proj0 = cubeProj[i].origin + cubeProj[i].min;
 320                 int64_t proj1 = cubeProj[i].origin + cubeProj[i].max;
 321
 322                 if (proj0 > inherit->tri_proj[i][1] ||
 323                         proj1 < inherit->tri_proj[i][0]) {
 324                         return 0;
 325                 }
 326         }
 327
 328         return 1;
 329 }
 330
 331 int CubeTriangleIsect::isIntersectingPrimary(int edgeInd) const
 332 {
 333         for (int i = 0; i < NUM_AXES; i++) {
 334
 335                 int64_t proj0 = cubeProj[i].origin;
 336                 int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
 337
 338                 if (proj0 < proj1) {
 339                         if (proj0 > inherit->tri_proj[i][1] ||
 340                                 proj1 < inherit->tri_proj[i][0]) {
 341                                 return 0;
 342                         }
 343                 }
 344                 else {
 345                         if (proj1 > inherit->tri_proj[i][1] ||
 346                                 proj0 < inherit->tri_proj[i][0]) {
 347                                 return 0;
 348                         }
 349                 }
 350
 351         }
 352
 353         // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
 354         return 1;
 355 }
 356
 357 float CubeTriangleIsect::getIntersectionPrimary(int edgeInd) const
 358 {
 359         int i = 3;
 360
 361
 362         int64_t proj0 = cubeProj[i].origin;
 363         int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
 364         int64_t proj2 = inherit->tri_proj[i][1];
 365         int64_t d = proj1 - proj0;
 366         double alpha;
 367
 368         if (d == 0)
 369                 alpha = 0.5;
 370         else {
 371                 alpha = (double)((proj2 - proj0)) / (double)d;
 372
 373                 if (alpha < 0 || alpha > 1)
 374                         alpha = 0.5;
 375         }
 376
 377         return (float)alpha;
 378 }