path: root/intern/dualcon/intern/Projections.cpp

/*
 * ***** BEGIN GPL LICENSE BLOCK *****
 *
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU General Public License
 * as published by the Free Software Foundation; either version 2
 * of the License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software Foundation,
 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 *
 * Contributor(s): Tao Ju
 *
 * ***** END GPL LICENSE BLOCK *****
 */

#ifdef WITH_CXX_GUARDEDALLOC
#  include "MEM_guardedalloc.h"
#endif

#include <math.h>
#include "Projections.h"

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}
};

/**
 * Method to perform cross-product
 */
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];
}

/**
 * Method to perform dot product
 */
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++;
		}
	}
}

/**
 * Construction from a cube (axes aligned) and triangle
 */
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]);
			if (proj > max) {
				max = proj;
			}
			if (proj < min) {
				min = proj;
			}
		}
		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++) {
			if (vts[i] < inherit->tri_proj[axis][0])
				inherit->tri_proj[axis][0] = vts[i];
			
			if (vts[i] > inherit->tri_proj[axis][1])
				inherit->tri_proj[axis][1] = vts[i];
		}
	}
}

/**
 * Construction
 * from a parent CubeTriangleIsect object and the index of the children
 */
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;
}


/**
 * Shifting a cube to a new origin
 */
void CubeTriangleIsect::shift(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]);
	}
}

/**
 * Method to test intersection of the triangle and the cube
 */
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;
}