Projections.cpp

Go to the documentation of this file.

/* SPDX-FileCopyrightText: 2002-2022 Blender Authors
 *
 * SPDX-License-Identifier: GPL-2.0-or-later */
 
#include "Projections.h"
#include <cmath>
 
#include <algorithm>
 
const int vertmap[8][3] = {
    {0, 0, 0},
    {0, 0, 1},
    {0, 1, 0},
    {0, 1, 1},
    {1, 0, 0},
    {1, 0, 1},
    {1, 1, 0},
    {1, 1, 1},
};
 
const int centmap[3][3][3][2] = {
    {{{0, 0}, {0, 1}, {1, 1}}, {{0, 2}, {0, 3}, {1, 3}}, {{2, 2}, {2, 3}, {3, 3}}},
 
    {{{0, 4}, {0, 5}, {1, 5}}, {{0, 6}, {0, 7}, {1, 7}}, {{2, 6}, {2, 7}, {3, 7}}},
 
    {{{4, 4}, {4, 5}, {5, 5}}, {{4, 6}, {4, 7}, {5, 7}}, {{6, 6}, {6, 7}, {7, 7}}}};
 
const int edgemap[12][2] = {
    {0, 4},
    {1, 5},
    {2, 6},
    {3, 7},
    {0, 2},
    {1, 3},
    {4, 6},
    {5, 7},
    {0, 1},
    {2, 3},
    {4, 5},
    {6, 7},
};
 
const int facemap[6][4] = {
    {0, 1, 2, 3},
    {4, 5, 6, 7},
    {0, 1, 4, 5},
    {2, 3, 6, 7},
    {0, 2, 4, 6},
    {1, 3, 5, 7},
};

static void crossProduct(int64_t res[3], const int64_t a[3], const int64_t b[3])
{
  res[0] = a[1] * b[2] - a[2] * b[1];
  res[1] = a[2] * b[0] - a[0] * b[2];
  res[2] = a[0] * b[1] - a[1] * b[0];
}
 
static void crossProduct(double res[3], const double a[3], const double b[3])
{
  res[0] = a[1] * b[2] - a[2] * b[1];
  res[1] = a[2] * b[0] - a[0] * b[2];
  res[2] = a[0] * b[1] - a[1] * b[0];
}

static int64_t dotProduct(const int64_t a[3], const int64_t b[3])
{
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
}
 
static void normalize(double a[3])
{
  double mag = a[0] * a[0] + a[1] * a[1] + a[2] * a[2];
  if (mag > 0) {
    mag = sqrt(mag);
    a[0] /= mag;
    a[1] /= mag;
    a[2] /= mag;
  }
}
 
/* Create projection axes for cube+triangle intersection testing.
 *    0, 1, 2: cube face normals
 *
 *          3: triangle normal
 *
 *    4, 5, 6,
 *    7, 8, 9,
 * 10, 11, 12: cross of each triangle edge vector with each cube
 *             face normal
 */
static void create_projection_axes(int64_t axes[NUM_AXES][3], const int64_t tri[3][3])
{
  /* Cube face normals */
  axes[0][0] = 1;
  axes[0][1] = 0;
  axes[0][2] = 0;
  axes[1][0] = 0;
  axes[1][1] = 1;
  axes[1][2] = 0;
  axes[2][0] = 0;
  axes[2][1] = 0;
  axes[2][2] = 1;
 
  /* Get triangle edge vectors */
  int64_t tri_edges[3][3];
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      tri_edges[i][j] = tri[(i + 1) % 3][j] - tri[i][j];
    }
  }
 
  /* Triangle normal */
  crossProduct(axes[3], tri_edges[0], tri_edges[1]);
 
  // Face edges and triangle edges
  int ct = 4;
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      crossProduct(axes[ct], axes[j], tri_edges[i]);
      ct++;
    }
  }
}

CubeTriangleIsect::CubeTriangleIsect(int64_t cube[2][3],
                                     int64_t tri[3][3],
                                     int64_t /*error*/,
                                     int triind)
{
  int i;
  inherit = new TriangleProjection;
  inherit->index = triind;
 
  int64_t axes[NUM_AXES][3];
  create_projection_axes(axes, tri);
 
  /* Normalize face normal and store */
  double dedge1[] = {(double)tri[1][0] - (double)tri[0][0],
                     (double)tri[1][1] - (double)tri[0][1],
                     (double)tri[1][2] - (double)tri[0][2]};
  double dedge2[] = {(double)tri[2][0] - (double)tri[1][0],
                     (double)tri[2][1] - (double)tri[1][1],
                     (double)tri[2][2] - (double)tri[1][2]};
  crossProduct(inherit->norm, dedge1, dedge2);
  normalize(inherit->norm);
 
  int64_t cubeedge[3][3];
  for (i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      cubeedge[i][j] = 0;
    }
    cubeedge[i][i] = cube[1][i] - cube[0][i];
  }
 
  /* Project the cube on to each axis */
  for (int axis = 0; axis < NUM_AXES; axis++) {
    CubeProjection &cube_proj = cubeProj[axis];
 
    /* Origin */
    cube_proj.origin = dotProduct(axes[axis], cube[0]);
 
    /* 3 direction vectors */
    for (i = 0; i < 3; i++) {
      cube_proj.edges[i] = dotProduct(axes[axis], cubeedge[i]);
    }
 
    /* Offsets of 2 ends of cube projection */
    int64_t max = 0;
    int64_t min = 0;
    for (i = 1; i < 8; i++) {
      int64_t proj = (vertmap[i][0] * cube_proj.edges[0] + vertmap[i][1] * cube_proj.edges[1] +
                      vertmap[i][2] * cube_proj.edges[2]);
      max = std::max(proj, max);
      min = std::min(proj, min);
    }
    cube_proj.min = min;
    cube_proj.max = max;
  }
 
  /* Project the triangle on to each axis */
  for (int axis = 0; axis < NUM_AXES; axis++) {
    const int64_t vts[3] = {dotProduct(axes[axis], tri[0]),
                            dotProduct(axes[axis], tri[1]),
                            dotProduct(axes[axis], tri[2])};
 
    // Triangle
    inherit->tri_proj[axis][0] = vts[0];
    inherit->tri_proj[axis][1] = vts[0];
    for (i = 1; i < 3; i++) {
      inherit->tri_proj[axis][0] = std::min(vts[i], inherit->tri_proj[axis][0]);
      inherit->tri_proj[axis][1] = std::max(vts[i], inherit->tri_proj[axis][1]);
    }
  }
}

CubeTriangleIsect::CubeTriangleIsect(CubeTriangleIsect *parent)
{
  // Copy inheritable projections
  this->inherit = parent->inherit;
 
  // Shrink cube projections
  for (int i = 0; i < NUM_AXES; i++) {
    cubeProj[i].origin = parent->cubeProj[i].origin;
 
    for (int j = 0; j < 3; j++) {
      cubeProj[i].edges[j] = parent->cubeProj[i].edges[j] >> 1;
    }
 
    cubeProj[i].min = parent->cubeProj[i].min >> 1;
    cubeProj[i].max = parent->cubeProj[i].max >> 1;
  }
}
 
unsigned char CubeTriangleIsect::getBoxMask()
{
  int i, j, k;
  int bmask[3][2] = {{0, 0}, {0, 0}, {0, 0}};
  unsigned char boxmask = 0;
  int64_t child_len = cubeProj[0].edges[0] >> 1;
 
  for (i = 0; i < 3; i++) {
    int64_t mid = cubeProj[i].origin + child_len;
 
    // Check bounding box
    if (mid >= inherit->tri_proj[i][0]) {
      bmask[i][0] = 1;
    }
    if (mid < inherit->tri_proj[i][1]) {
      bmask[i][1] = 1;
    }
  }
 
  // Fill in masks
  int ct = 0;
  for (i = 0; i < 2; i++) {
    for (j = 0; j < 2; j++) {
      for (k = 0; k < 2; k++) {
        boxmask |= ((bmask[0][i] & bmask[1][j] & bmask[2][k]) << ct);
        ct++;
      }
    }
  }
 
  // Return bounding box masks
  return boxmask;
}

void CubeTriangleIsect::shift(const int off[3])
{
  for (int i = 0; i < NUM_AXES; i++) {
    cubeProj[i].origin += (off[0] * cubeProj[i].edges[0] + off[1] * cubeProj[i].edges[1] +
                           off[2] * cubeProj[i].edges[2]);
  }
}

int CubeTriangleIsect::isIntersecting() const
{
  for (int i = 0; i < NUM_AXES; i++) {
    /*
      int64_t proj0 = cubeProj[i][0] +
      vertmap[inherit->cubeEnds[i][0]][0] * cubeProj[i][1] +
      vertmap[inherit->cubeEnds[i][0]][1] * cubeProj[i][2] +
      vertmap[inherit->cubeEnds[i][0]][2] * cubeProj[i][3] ;
      int64_t proj1 = cubeProj[i][0] +
      vertmap[inherit->cubeEnds[i][1]][0] * cubeProj[i][1] +
      vertmap[inherit->cubeEnds[i][1]][1] * cubeProj[i][2] +
      vertmap[inherit->cubeEnds[i][1]][2] * cubeProj[i][3] ;
    */
 
    int64_t proj0 = cubeProj[i].origin + cubeProj[i].min;
    int64_t proj1 = cubeProj[i].origin + cubeProj[i].max;
 
    if (proj0 > inherit->tri_proj[i][1] || proj1 < inherit->tri_proj[i][0]) {
      return 0;
    }
  }
 
  return 1;
}
 
int CubeTriangleIsect::isIntersectingPrimary(int edgeInd) const
{
  for (int i = 0; i < NUM_AXES; i++) {
 
    int64_t proj0 = cubeProj[i].origin;
    int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
 
    if (proj0 < proj1) {
      if (proj0 > inherit->tri_proj[i][1] || proj1 < inherit->tri_proj[i][0]) {
        return 0;
      }
    }
    else {
      if (proj1 > inherit->tri_proj[i][1] || proj0 < inherit->tri_proj[i][0]) {
        return 0;
      }
    }
  }
 
  // printf( "Intersecting: %d %d\n", edgemap[edgeInd][0], edgemap[edgeInd][1] )  ;
  return 1;
}
 
float CubeTriangleIsect::getIntersectionPrimary(int edgeInd) const
{
  int i = 3;
 
  int64_t proj0 = cubeProj[i].origin;
  int64_t proj1 = cubeProj[i].origin + cubeProj[i].edges[edgeInd];
  int64_t proj2 = inherit->tri_proj[i][1];
  int64_t d = proj1 - proj0;
  double alpha;
 
  if (d == 0) {
    alpha = 0.5;
  }
  else {
    alpha = (double)((proj2 - proj0)) / (double)d;
 
    if (alpha < 0 || alpha > 1) {
      alpha = 0.5;
    }
  }
 
  return (float)alpha;
}