Thread: Intersection Test Issues (Matrix Calculator Check)

Hello All,

I'm having difficulties in implementing or at least checking this code that I wrote.

Basically, it's an intersection test between a point and a tetrahedron. I'm using the following paper as reference : Optimized Spatial Hashing for Collision Detection of Deformable Objects by Teschner et al.

The test is outlined in section 4.4 and this seems to be a common standard so I decided to build it myself only I think I did it wrong.

I'm posting my commented out code. It's long and it should work but I'm not sure. Also, I lose long arithmetic because I've heard round-off errors can damage the mesh I want to grow later so I'm trying to mitigate that.

I have a particle set with positions 0 to 7 in all dimensions so we have (0, 0, 0) to (7, 7, 7). My initial tetrahedron has the vertices, (0, 0, 0), (7, 0, 0), (0, 7, 0), (0, 0, 7) and I'm checking intersections for the whole set (512 points). I know I'm supposed to get like 1/6 of that, right? But I get 120 points.

Let me know if not enough info was provided.

main.c :

Code:

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>

#include "structures.h"


int main(void) {


/* Initialize particle structure with relevant data */


   struct particle *parts = malloc(512*sizeof(*parts));
   for (int i=0; i<8; i++) {
      for (int j=0; j<8; j++) {
         for (int k=0; k<8; k++) {
            parts[i*64+j*8+k].p[0] = i;
            parts[i*64+j*8+k].p[1] = j;
            parts[i*64+j*8+k].p[2] = k;
         }
      }      
   }


/* Begin intersection test for particle set */


   int x=0;
   for (int i=0; i<512; i++) {


      if (cofact(parts[0].p, parts[7].p, parts[448].p, parts[56].p, parts[i].p)) x++;
   }


   printf("And this is how many points are in the tetrahedron : %d\n", x);


/* Free'ing routine for our allocated structures */


   free(parts);


   return 0;
}

structures.h :

Code:

/* Our structure and function library */


struct particle {


   long long p[3];
   long long v[3];
   long long a[3];


   long long mass;


   int type;


   long long peanokey;
};


long long matrix(long long a[3], long long b[3], long long c[3], long long d[3], long long e[3], int rank);


long long det(long long **a, int rank);


bool cofact(long long a[3], long long b[3], long long c[3], long long d[3], long long e[3]);

det.c :

Code:

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>


#include "structures.h"


/* This is the first step of our determinant calculator. 
   We take 4 or 5 points as input, build a formal array
   and then calculate the determinant with the 'det'
   function below. */


long long matrix(long long a[3], long long b[3], long long c[3], long long d[3], long long e[3], int rank) { 


   long long **m;
   m = malloc(rank*sizeof(*m));
   for (int i=0; i<rank; i++) {


      m[i] = malloc(rank*sizeof(*m[i]));
      m[i][0] = 1;
   }


      m[0][1] = a[0]; m[0][2] = a[1]; m[0][3] = a[2];


      m[1][1] = b[0]; m[1][2] = b[1]; m[1][3] = b[2];


      m[2][1] = c[0]; m[2][2] = c[1]; m[2][3] = c[2];


      m[3][1] = d[0]; m[3][2] = d[1]; m[3][3] = d[2];


   if (rank==5) {


      m[0][4] = a[0]*a[0]+a[1]*a[1]+a[2]*a[2];
      m[1][4] = b[0]*b[0]+b[1]*b[1]+b[2]*b[2];
      m[2][4] = c[0]*c[0]+c[1]*c[1]+c[2]*c[2];
      m[3][4] = d[0]*d[0]+d[1]*d[1]+d[2]*d[2];


      m[4][1] = e[0]; m[4][2] = e[1]; m[4][3] = e[2]; m[4][4] = e[0]*e[0]+e[1]*e[1]+e[2]*e[2];
   }


   for (int i=0; i<rank; i++) {
      for (int j=0; j<rank; j++) {
         printf(" %lld ", m[i][j]);
      } printf("\n");
   }


   long long n = det(m, rank);


   for (int i=0; i<rank; i++) {
      free(m[i]);
   }


   free(m);


   return n;
}


/* Our determinant calculator which takes a NxN array and
   will build smaller sub-arrays until we eventually hit the
   3x3 floor which is ideal for our case as we will never
   take a determinant of something smaller. */


long long det(long long **a, int rank) {


   long long m = 0;


   if (rank == 3) {


      m = a[0][0]*(a[1][1]*a[2][2] - a[1][2]*a[2][1])
        - a[0][1]*(a[1][0]*a[2][2] - a[1][2]*a[2][0])
        + a[0][2]*(a[1][0]*a[2][1] - a[1][1]*a[2][0]);


      return m;
   }


   for (int i=0; i<rank; i++) {


      long long **b = malloc((rank-1)*sizeof(*b));


      for (int i=0; i<(rank-1); i++) {
         b[i] = malloc((rank-1)*sizeof(*b[i]));
      }


      for (int j=1; j<rank; j++) {
         int c = 0;
         for (int k=0; k<rank; k++) {
            if (k==i) continue;
            b[j-1][c] = a[j][k];
            c++;
         }
      }


      if (i%2==0) {m += a[0][i]*det(b, rank-1);}
      else m -= a[0][i]*det(b, rank-1);  
      for (int i=0; i<(rank-1); i++) free(b[i]);
      free(b);
   }


   return m;
}


/* This is our intersection test. We take 5 points and determine
   if e is inside the tetrahedron spanned by a, b, c and d. First,
   we allocate two formal matrices based on our points. m is the
   matrix highlighted by the paper (A) and l is a copy. We set
   a static n to be the cofactor matrix of m and then set m to be
   the transpose of n. Once we have the adjoint, we divide by
   the determinant as is required by maths. fin is our final
   inversed array. */


/* Once we have our final inversed array, we then do simple matrix
   multiplication with our point to be tested which is our 
   t array. As conditions require, if every element of t is
   greater than 0 and the sum of all 3 numbers is still less
   than or equal to 1, we have an intersection. Our boolean
   variable 'silh' is the marker of this and returns true
   if and only if the point is inside the tetrahedron. */


bool cofact(long long a[3], long long b[3], long long c[3], long long d[3], long long e[3]) {


   long long **m, **l;


   m = malloc(3*sizeof(*m));
   l = malloc(3*sizeof(*l));


   for (int i=0; i<3; i++) {


      m[i] = malloc(3*sizeof(*m[i]));
      l[i] = malloc(3*sizeof(*l[i]));
   }


   m[0][0] = b[0] - a[0]; l[0][0] = m[0][0];
   m[0][1] = b[1] - a[1]; l[0][1] = m[0][1];
   m[0][2] = b[2] - a[2]; l[0][2] = m[0][2];


   m[1][0] = c[0] - a[0]; l[1][0] = m[1][0];
   m[1][1] = c[1] - a[1]; l[1][1] = m[1][1];
   m[1][2] = c[2] - a[2]; l[1][2] = m[1][2];


   m[2][0] = d[0] - a[0]; l[2][0] = m[2][0];
   m[2][1] = d[1] - a[1]; l[2][1] = m[2][1];
   m[2][2] = d[2] - a[2]; l[2][2] = m[2][2];


   /*printf("Before :\n");
   for (int i=0; i<3; i++) {
      for (int j=0; j<3; j++) {
         printf(" %lld ", m[i][j]);
      } printf("\n");
   }*/
  
   long long n[3][3];


/* Cofactorization */


   n[0][0] = m[1][1]*m[2][2] - m[1][2]*m[2][1];
   n[0][1] = -1*(m[1][0]*m[2][2] -  m[2][0]*m[1][2]);
   n[0][2] = m[1][0]*m[2][1] - m[2][0]*m[1][1];


   n[1][0] = -1*(m[0][1]*m[2][2] - m[0][2]*m[2][1]);
   n[1][1] = m[0][0]*m[2][2] - m[2][0]*m[0][2];
   n[1][2] = -1*(m[0][0]*m[2][1] - m[2][0]*m[0][1]);


   n[2][0] = m[0][1]*m[1][2] - m[1][1]*m[0][2];
   n[2][1] = -1*(m[0][0]*m[1][2] - m[1][0]*m[0][2]);
   n[2][2] = m[0][0]*m[1][1] - m[1][0]*m[0][1];


/* Transposition */


   m[0][0] = n[0][0]; m[1][1] = n[1][1]; m[2][2] = n[2][2];
   m[0][1] = n[1][0]; m[1][0] = n[0][1];
   m[0][2] = n[2][0]; m[2][0] = n[0][2];
   m[2][1] = n[1][2]; m[1][2] = n[2][1];


   /*printf("After :\n");
   for (int i=0; i<3; i++) {
      for (int j=0; j<3; j++) {
         printf(" %lld ", m[i][j]);
      } printf("\n");
   }*/


/* Determination */


   long long z = det(l, 3);
   //printf("det of adjoint : %lld\n", z);


   long double fin[3][3];


   for (int i=0; i<3; i++) {
      for (int j=0; j<3; j++) {
         fin[i][j] = (long double)m[i][j]/z;
      }
   }


   /*printf("Then :\n");
   for (int i=0; i<3; i++) {
      for (int j=0; j<3; j++) {
         printf(" %Lf ", fin[i][j]);
      } printf("\n");
   }*/


   long double t[3];


   t[0] = fin[0][0]*(long double)(e[0]-a[0]) + fin[0][1]*(long double)(e[1]-a[1]) + fin[0][2]*(long double)(e[2]-a[2]);
   t[1] = fin[1][0]*(long double)(e[0]-a[0]) + fin[1][1]*(long double)(e[1]-a[1]) + fin[1][2]*(long double)(e[2]-a[2]);
   t[2] = fin[2][0]*(long double)(e[0]-a[0]) + fin[2][1]*(long double)(e[1]-a[1]) + fin[2][2]*(long double)(e[2]-a[2]);


   for (int i=0; i<3; i++) {
      free(m[i]);
      free(l[i]);
   }
   
   free(m); free(l);


   long double silh = 0;


   for (int i=0; i<3; i++) {


      if (t[i]<0) return false;
      silh += t[i];
   }


   if (silh<=1) return true;


   return false;
}

Makefile :

Code:

MKFILE    = Makefile


GCC       = gcc -g -O0 -Wall -Wextra -std=gnu99


CSOURCE   = main.c det.c
CHEADER   = structures.h
OBJECTS   = ${CSOURCE:.c=.o}
EXECBIN   = tetra


all : ${EXECBIN}


${EXECBIN} : ${OBJECTS}
	${GCC} -o $@ ${OBJECTS} -lm




%.o : %.c


	${GCC} -c $<




clean :
	- rm ${OBJECTS} ${EXECBIN}




again :
	${MAKE} clean all

Sorry for all the posts, but this should be the most human readable format and the easiest to replicate on a home machine. Note that my code does NOT leak so don't worry 'bout nothin' like that. My inverse calculator seems to be working fine so I think everything I'm doing is literally correct but it's just that, shouldn't I have like 512/6 points in this? And that's like 85 something, not 120. Any help?

Alright, so I modified my code to have it print out the points it was containing. If you have graph paper like me then you can try this at home! I have the box with bounds 0 to 7 in each direction and my code actually cuts it in half perfectly along the diagonal. So at x=0, it contains all 7 y points. At x = 1, it contains 0 through 6, x=2 contains 0 through 5, x=3 contains 0 through 4, x=4 contains 0 through 3, etc.

But this was all done for z=0. As z goes up, the thins out but it does so one element at a time which is predicted.

So for z=1, x=0 contains y points 0 through 6, x=1 contains 0 through 5...

Basically, the answer is this: (8+7+...+1) + (7+6+....+1)+(6+5+...+1)...+(2+1)+1 = 120.

I just did that on my calculator and I was like, uh, son, uh! My code is drawing edges correctly! The containment routine I wrote was written successfully or at least it was written well enough such that I can now account for every point contained and it's a good thing. I did it, C Board, I did it! I can continue with my code now!

it is amazing how sometimes when you describe the prob from the beginning to someone else, it clears your mind up to see the solution!!!

and that, ladies and gentleman, is why i talk to myself!!! lol