ability to get barycentric coordinates

include/bary.cuh

@@ -5,6 +5,7 @@
 
 namespace bary{
  bool* point_in_simplex(vec3f *pts, int n, int dim, vec3f *verts);
+ float** bary_triangle(vec3f *pts, int n, int dim, vec3f *verts);
 }
 
 #endif //BARY_H_

src/bary.cu

@@ -8,7 +8,10 @@
 #include <cstdio>
 #include "bary.cuh"
 
-#define N_THREADS_PER_BLOCK 512
+#define N_THREADS_PER_BLOCK 64
+#define BLOCKSIZE_x 3
+#define BLOCKSIZE_y 32
+
 
 #define cudaErrchk(ans) { gpuSay((ans), __FILE__, __LINE__);}
 
@@ -44,18 +47,11 @@ namespace bary{
  return true;
  }
 
- __global__ void pt_in_tri(vec3f *pts, int npts, vec3f *a, vec3f *b, vec3f *c, bool *res){
-
- int i = blockIdx.x * blockDim.x + threadIdx.x;
- i = max(i, 0);
- i = min(i, npts);
 
-  vec3f *p = &pts[i];
+ __device__ void _bary_tri(vec3f *p, int n, vec3f *a, vec3f *b, vec3f *c, float *res){
  vec3f ap[3], bp[3], cp[3];
  vec3f ac[3], ab[3], ca[3], bc[3];
- vec3f n[3], n0[3], n1[3], n2[3];
-
- float bary[3];
+ vec3f nor[3], nor0[3], nor1[3], nor2[3];
 
  // Vectors from vertices to the point
  _diff(ap, a, p);
@@ -70,25 +66,48 @@ namespace bary{
 
  // Make normals. The lengths are proportional to the
  // area of the subtriangles
- _cross(n, ab, ac);
- _cross(n0, bc, bp);
- _cross(n1, ca, cp);
- _cross(n2, ab, ap);
+ _cross(nor, ab, ac);
+ _cross(nor0, bc, bp);
+ _cross(nor1, ca, cp);
+ _cross(nor2, ab, ap);
 
- float n_norm = _dot(n, n);
+ float n_norm = _dot(nor, nor);
 
  // Find relative areas
- bary[0] = _dot(n, n0) / n_norm;
- bary[1] = _dot(n, n1) / n_norm;
- bary[2] = _dot(n, n2) / n_norm;
+ res[0] = _dot(nor, nor0) / n_norm;
+ res[1] = _dot(nor, nor1) / n_norm;
+ res[2] = _dot(nor, nor2) / n_norm;
+ }
+
+
+ __global__ void bary_tri(vec3f *pts, int n, vec3f *a, vec3f *b, vec3f *c, size_t pitch, float *res){
+ int tidx = blockIdx.x*blockDim.x + threadIdx.x;
+ int tidy = blockIdx.y*blockDim.y + threadIdx.y;
+
+ if ((tidx < 3) && (tidy < n)){
+ vec3f *p = &pts[tidy];
+ float bary[3];
+ _bary_tri(p, n, a, b, c, bary);
+ res[tidy * pitch] = bary[0];
+ res[tidy * pitch + 1] = bary[1];
+ res[tidy * pitch + 2] = bary[2];
+ }
+ }
 
+ __global__ void pt_in_tri(vec3f *pts, int n, vec3f *a, vec3f *b, vec3f *c, bool *res){
+ int i = blockIdx.x * blockDim.x + threadIdx.x;
+ i = max(i, 0);
+ i = min(i, n);
+ vec3f *p = &pts[i];
+ float bary[3];
+ _bary_tri(p, n, a, b, c, bary);
  res[i] = _all(bary, 3);
  }
 
- __global__ void pt_in_tet(vec3f *pts, int npts, vec3f *a, vec3f *b, vec3f *c, vec3f *d, bool *res){
+ __global__ void pt_in_tet(vec3f *pts, int n, vec3f *a, vec3f *b, vec3f *c, vec3f *d, bool *res){
  int i = blockIdx.x * blockDim.x + threadIdx.x;
  i = max(i, 0);
- i = min(i, npts);
+ i = min(i, n);
 
  vec3f *p = &pts[i];
  vec3f ap[3], bp[3];
@@ -127,6 +146,49 @@ namespace bary{
 
  }
 
+ int iDivUp(int hostPtr, int b){ return ((hostPtr % b) != 0) ? (hostPtr / b + 1) : (hostPtr / b); }
+
+ __host__ float** bary_triangle(vec3f *pts, int n, int ndim, vec3f *verts){
+ vec3f *pts_d;
+ vec3f *a_d, *b_d, *c_d;
+ float *res_d;
+ float* res = new float[n * ndim];
+ float** res2d = new float*[n];
+ size_t pitch;
+ for( int i = 0 ; i < n ; i++ ){
+ res2d[i] = &res[i * ndim];
+ }
+
+ cudaErrchk(cudaMalloc(&pts_d, n * sizeof(vec3f)));
+ cudaErrchk(cudaMalloc(&a_d, sizeof(vec3f)));
+ cudaErrchk(cudaMalloc(&b_d, sizeof(vec3f)));
+ cudaErrchk(cudaMalloc(&c_d, sizeof(vec3f)));
+ // cudaErrchk(cudaMalloc(&res_d, sizeof(res)));
+
+ cudaErrchk(cudaMallocPitch((void**)&res_d, &pitch, ndim*sizeof(float), n));
+ cudaErrchk(cudaMemcpy2D(res_d, pitch, res, ndim*sizeof(float), ndim*sizeof(float), n, cudaMemcpyHostToDevice));
+
+ cudaErrchk(cudaMemcpy(pts_d, pts, n * sizeof(vec3f), cudaMemcpyHostToDevice));
+ cudaErrchk(cudaMemcpy(a_d, &verts[0], sizeof(vec3f), cudaMemcpyHostToDevice));
+ cudaErrchk(cudaMemcpy(b_d, &verts[1], sizeof(vec3f), cudaMemcpyHostToDevice));
+ cudaErrchk(cudaMemcpy(c_d, &verts[2], sizeof(vec3f), cudaMemcpyHostToDevice));
+
+ dim3 gridSize(iDivUp(ndim, BLOCKSIZE_x), iDivUp(n, BLOCKSIZE_y));
+ dim3 blockSize(BLOCKSIZE_y, BLOCKSIZE_x);
+
+ bary_tri<<<gridSize, blockSize>>>(pts_d, n, a_d, b_d, c_d, pitch, res_d);
+
+ cudaErrchk(cudaMemcpy2D(res, ndim*sizeof(float), res_d, pitch, ndim*sizeof(float), n, cudaMemcpyDeviceToHost));
+
+ cudaFree(pts_d);
+ cudaFree(a_d);
+ cudaFree(b_d);
+ cudaFree(c_d);
+ cudaFree(res_d);
+
+ return res2d;
+ }
+
 
  /**
  Checks if a bunch of points lie inside a 2D/3D simplex
@@ -135,7 +197,7 @@ namespace bary{
  @param verts array of simplex vertices
  @return res a boolean array
  */
- __host__ bool* point_in_simplex(vec3f *pts, int n, int dim, vec3f *verts) {
+ __host__ bool* point_in_simplex(vec3f *pts, int n, int ndim, vec3f *verts) {
  vec3f *pts_d;
  vec3f *a_d, *b_d, *c_d, *d_d;
  bool *res_d;
@@ -153,7 +215,7 @@ namespace bary{
  cudaErrchk(cudaMemcpy(c_d, &verts[2], sizeof(vec3f), cudaMemcpyHostToDevice));
 
  // Check if the simplex is 2D or 3D
- if(dim == 4){
+ if(ndim == 4){
  cudaErrchk(cudaMalloc(&d_d, sizeof(vec3f)));
  cudaErrchk(cudaMemcpy(d_d, &verts[3], sizeof(vec3f), cudaMemcpyHostToDevice));
  pt_in_tet<<<1+(n/N_THREADS_PER_BLOCK), N_THREADS_PER_BLOCK>>>(pts_d, n, a_d, b_d, c_d, d_d, res_d);

src/testBary.cpp

@@ -1,9 +1,10 @@
 #include <iostream>
+#include <stdio.h>
 #include <stdlib.h>
 #include <time.h>
 #include "bary.cuh"
 
-#define N 1000
+#define N 30
 
 int main(){
  srand(time(NULL));
@@ -27,15 +28,21 @@ int main(){
  vec3f t2 = {rand()%255, rand()%255, rand()%255};
  vec3f t3 = {rand()%255, rand()%255, rand()%255};
 
- vec3f *verts = new vec3f[4];
+ vec3f *verts = new vec3f[3];
  verts[0] = t0;
  verts[1] = t1;
  verts[2] = t2;
- verts[3] = t3;
+ // verts[3] = t3;
 
- bool* res = bary::point_in_simplex(pts_h, N, 4, verts);
+ // bool* res = bary::point_in_simplex(pts_h, N, 3, verts);
 
- for(int i=N; i--; std::cout << res[i]){}
+ float** res = bary::bary_triangle(pts_h, N, 3, verts);
+
+ for(int i=N; i--; std::cout << res[i][0] << res[i][1] << res[i][2] << std::endl){}
+ // std::cout << res[0][0];
+ // for(int i=0; i<N; i++){
+ // printf("%f\n", res[i]);
+ // }
 
  return 0;
 }