#ifndef CPP_INTEGRATION_KERNEL_H
define CPP_INTEGRATION_KERNEL_H
include <stdio.h>
include <stdlib.h>
include <math.h>
define TRI_NUM1 960 //sphere // 7580 cup
define TRI_NUM2 960//216 cylinder //26012 spoon //5250 azucar
define BSIZE 10 // sphere&cylinder //4 cup&spoon //10 cup&azucar
define WGRID TRI_NUM2/BSIZE
define HGRID TRI_NUM2/BSIZE
/* macros
define CROSS(dest, v1, v2, i2, i1) \
dest[0] = v1[i1*3+1]*v2[i2*3+2]-v1[i1*3+2]*v2[i2*3+1]; \
dest[1] = v1[i1*3+2]*v2[i2*3]-v1[i1*3]*v2[i2*3+2]; \
dest[2] = v1[i1*3]*v2[i2*3+1]-v1[i1*3+1]*v2[i2*3];
define DOT(v1, v2, i1, i2) (v1[i1*3]*v2[i2*3]+v1[i1*3+1]*v2[i2*3+1]+v1[i1*3+2]*v2[i2*3+2])
define SUBS(dest, v1, v2, i1, i2) \
dest[0] = v1[i1*3]-v2[i2*3]; \
dest[1] = v1[i1*3+1]-v2[i2*3+1]; \
dest[2] = v1[i1*3+2]-v2[i2*3+2];
define ADD(dest, v1, v2, i1, i2) \
dest[0] = v1[i1*3]+v2[i2*3]; \
dest[1] = v1[i1*3+1]+v2[i2*3+1]; \
dest[2] = v1[i1*3+2]+v2[i2*3+2];
define MULT(dest, v, factor, i) \
dest[0] = factor*v[i1*3]; \
dest[1] = factor*v[i1*3+1]; \
dest[2] = factor*v[i1*3+2];
define SORT(a, B) \
if(a>b) { \
float c; \
c = a; \
a = b; \
b = c; \
}
define ISECT(vv0, vv1, vv2, D0, D1, D2, isect0, isect1) \
isect0 = vv0+(vv1-vv0)*D0/(D0-D1); \
isect1 = vv0+(vv2-vv0)*D0/(D0-D2);
*/
device short COMPUTE_INTERVALS(float vv0, float vv1, float vv2, float D0, float D1, float D2,
float d0d1, float d0d2, float* isect)
{
if(d0d1>0.0f) {
isect[0] = vv2+(vv0-vv2)*D2/(D2-D0);
isect[1] = vv2+(vv1-vv2)*D2/(D2-D1);
return 0;
}
else if(d0d2>0.0f) {
isect[0] = vv1+(vv0-vv1)*D1/(D1-D0);
isect[1] = vv1+(vv2-vv1)*D1/(D1-D2);
return 0;
}
else if(D1*D2>0.0f || D0!=0.0f) {
isect[0] = vv0+(vv1-vv0)*D0/(D0-D1);
isect[1] = vv0+(vv2-vv0)*D0/(D0-D2);
return 0;
}
else if(D1!=0.0f) {
isect[0] = vv1+(vv0-vv1)*D1/(D1-D0);
isect[1] = vv1+(vv2-vv1)*D1/(D1-D2);
return 0;
}
else if(D2!=0.0f) {
isect[0] = vv2+(vv0-vv2)*D2/(D2-D0);
isect[1] = vv2+(vv1-vv2)*D2/(D2-D1);
return 0;
}
else {
return 1;
//g_odata[tri1*TRI_NUM2+tri2] = coplanar(vn, s_v0, s_v1, s_v2, s_u0, s_u1, s_u2, v, u);
}
}
device short EDGE_EDGE_TEST(floatv0, float u0, float* u1, float Ax, float Ay, int v, int u, int i0, int i1)
{
float Bx, By, Cx, Cy, e, d, f;
Bx = u0[u*3+i0]-u1[u*3+i0];
By = u0[u*3+i1]-u1[u*3+i1];
Cx = v0[v*3+i0]-u0[u*3+i0];
Cy = v0[v*3+i1]-u0[u*3+i1];
f = AyBx-Ax By;
d = ByCx-Bx Cy;
if((f>0 && d>=0 && d<=f) || (f<0 && d<=0 && d>=f))
{
e=AxCy-Ay Cx;
if(f>0) {
if(e>=0 && e<=f) return 1;
}
else {
if(e<=0 && e>=f) return 1;
}
}
return 0;
}
device short EDGE_AGAINST_TRI_EDGES(float* v0, float* v1, float* u0, float* u1, float* u2, int v, int u, int i0, int i1)
{
float Ax, Ay;
Ax = v1[v3+i0]-v0[v 3+i0];
Ay = v1[v3+i1]-v0[v 3+i1];
if(EDGE_EDGE_TEST(v0, u0, u1, Ax, Ay, v, u, i0, i1))
return 1;
else if(EDGE_EDGE_TEST(v0, u1, u2, Ax, Ay, v, u, i0, i1))
return 1;
return EDGE_EDGE_TEST(v0, u2, u0, Ax, Ay, v, u, i0, i1);
}
device short POINT_IN_TRI(float *v0, float *u0, float *u1, float *u2, int v, int u, int i0, int i1)
{
float a, b, c, d0, d1, d2;
a = u1[u*3+i1]-u0[u*3+i1];
b = -(u1[u*3+i0]-u0[u*3+i0]);
c = -au0[u*3+i0]-b u0[u*3+i1];
d0 = av0[v 3+i0]+bv0[v 3+i1]+c;
a = u2[u*3+i1]-u1[u*3+i1];
b = -(u2[u*3+i0]-u1[u*3+i0]);
c = -au1[u*3+i0]-b u1[u*3+i1];
d1 = av0[v 3+i0]+bv0[v 3+i1]+c;
a = u0[u*3+i1]-u2[u*3+i1];
b = -(u0[u*3+i0]-u2[u*3+i0]);
c = -au2[u*3+i0]-b u2[u*3+i1];
d2 = av0[v 3+i0]+bv0[v 3+i1]+c;
if(d0*d1>0.0) {
if(d0*d2>0.0) return 1;
}
return 0;
}
device short coplanar(float *n, float *v0, float *v1, float *v2, float *u0, float *u1, float *u2, int v, int u)
{
float A[3];
int i0, i1;
A[0] = fabsf(n[0]);
A[1] = fabsf(n[1]);
A[2] = fabsf(n[2]);
if(A[0]>A[1]) {
if(A[0]>A[2]) {
i0 = 1;
i1 = 2;
}
else {
i0 = 0;
i1 = 1;
}
}
else {
if(A[2]>A[1]) {
i0 = 0;
i1 = 1;
}
else {
i0 = 0;
i1 = 2;
}
}
if(EDGE_AGAINST_TRI_EDGES(v0, v1, u0, u1, u2, v, u, i0, i1))
return 1;
else if(EDGE_AGAINST_TRI_EDGES(v1, v2, u0, u1, u2, v, u, i0, i1))
return 1;
else if(EDGE_AGAINST_TRI_EDGES(v2, v0, u0, u1, u2, v, u, i0, i1))
return 1;
else if(POINT_IN_TRI(v0, u0, u1, u2, v, u, i0, i1))
return 1;
else
return POINT_IN_TRI(u0, v0, v1, v2, v, u, i0, i1);
}
global void
kernel( short* g_odata, float* g_v0, float* g_v1, float* g_v2, float* g_u0, float* g_u1, float* g_u2 )
{
// shared memory
// the size is determined by the host application
shared float s_v0[BSIZE*3];
shared float s_v1[BSIZE*3];
shared float s_v2[BSIZE*3];
shared float s_u0[BSIZE*3];
shared float s_u1[BSIZE*3];
shared float s_u2[BSIZE3]; // BSIZE 3가 아니라 3으로 해야하나?
// access thread id
int tri1 = BSIZE * blockIdx.y + threadIdx.y;
int tri2 = BSIZE * blockIdx.x + threadIdx.x;
unsigned int i;
int v = threadIdx.y;
int u = threadIdx.x;
g_odata[tri1*TRI_NUM2+tri2] = 2; //g_odata initialization
if(threadIdx.x == 0){
for(i=0; i<3; i++){
s_v0[v3+i] = g_v0[tri1 3+i];
s_v1[v3+i] = g_v1[tri1 3+i];
s_v2[v3+i] = g_v2[tri1 3+i];
}
}
if(threadIdx.y == 0){
for(i=0; i<3; i++){
s_u0[u*3+i] = g_u0[tri2*3+i];
s_u1[u*3+i] = g_u1[tri2*3+i];
s_u2[u*3+i] = g_u2[tri2*3+i];
}
}
//이렇게 넘겨주는게 맞나 잘 모르겠음
__syncthreads();
// 계산에 쓰이는 변수들.
float ve1[3];
float ve2[3];
float ue1[3];
float ue2[3];
float vn[3];
float un[3];
float vd;
float ud;
float dv0;
float dv1;
float dv2;
float du0;
float du1;
float du2;
float du0du1;
float du0du2;
float dv0dv1;
float dv0dv2;
float D[3];
float vp0;
float vp1;
float vp2;
float up0;
float up1;
float up2;
float m, b, c;
short index;
float t1[2], t2[2];
float temp;
// Compute plane equation of triangle (s_v0, s_v1, s_v2)
for(i=0; i<3; i++)
{
ve1[i] = s_v1[v3+i]-s_v0[v 3+i];
}
for(i=0; i<3; i++)
{
ve2[i] = s_v2[v3+i]-s_v0[v 3+i];
}
vn[0] = ve1[1]*ve2[2] - ve1[2]*ve2[1];
vn[1] = ve1[2]*ve2[0] - ve1[0]*ve2[2];
vn[2] = ve1[0]*ve2[1] - ve1[1]*ve2[0];
vd = -(vn[0]s_v0[v 3] + vn[1]s_v0[v 3+1] + vn[2]s_v0[v 3+2]);
du0 = (vn[0]*s_u0[u*3] + vn[1]*s_u0[u*3+1] + vn[2]*s_u0[u*3+2]) + vd;
du1 = (vn[0]*s_u1[u*3] + vn[1]*s_u1[u*3+1] + vn[2]*s_u1[u*3+2]) + vd;
du2 = (vn[0]*s_u2[u*3] + vn[1]*s_u2[u*3+1] + vn[2]*s_u2[u*3+2]) + vd;
du0du1 = du0*du1;
du0du2 = du0*du2;
if(du0du1>0.0f && du0du2>0.0f){
g_odata[tri1*TRI_NUM2+tri2] = 0; /////////////////////////////
}
if(g_odata[tri1*TRI_NUM2+tri2] == 2) // 이 위치 맞나?
{
// compute plane equation of triangle (s_u0, s_u1, s_u2)
for(i=0; i<3; i++)
{
ue1[i] = s_u1[u*3+i]-s_u0[u*3+i];
}
for(i=0; i<3; i++)
{
ue2[i] = s_u2[u*3+i]-s_u0[u*3+i];
}
un[0] = ue1[1]*ue2[2] - ue1[2]*ue2[1];
un[1] = ue1[2]*ue2[0] - ue1[0]*ue2[2];
un[2] = ue1[0]*ue2[1] - ue1[1]*ue2[0];
ud = -(un[0]*s_u0[u*3] + un[1]*s_u0[u*3+1] + un[2]*s_u0[u*3+2]);
dv0 = (un[0]s_v0[v 3] + un[1]s_v0[v 3+1] + un[2]s_v0[v 3+2]) + ud;
dv1 = (un[0]s_v1[v 3] + un[1]s_v1[v 3+1] + un[2]s_v1[v 3+2]) + ud;
dv2 = (un[0]s_v2[v 3] + un[1]s_v2[v 3+1] + un[2]s_v2[v 3+2]) + ud;
dv0dv1 = dv0*dv1;
dv0dv2 = dv0*dv2;
if(dv0dv1>0.0f && dv0dv2>0.0f){
g_odata[tri1*TRI_NUM2+tri2] = 0; // 맞는지 확인해야함
}
if(g_odata[tri1*TRI_NUM2+tri2] == 2) // 이 위치 맞나?
{
// compute direction of intersection line
D[0] = vn[1]*un[2] - vn[2]*un[1];
D[1] = vn[2]*un[0] - vn[0]*un[2];
D[2] = vn[0]*un[1] - vn[1]*un[0];
// optimization step
// compute the largest component of D, determine the projected verteces onto L and compute the interval
m = fabs(D[0]);
index = 0;
b = fabs(D[1]);
c = fabs(D[2]);
if(b>m) m=b, index=1;
if(c>m) m=c, index=2;
// This is simplified projection into L
vp0 = s_v0[v*3+index];
vp1 = s_v1[v*3+index];
vp2 = s_v2[v*3+index];
up0 = s_u0[u*3+index];
up1 = s_u1[u*3+index];
up2 = s_u2[u*3+index];
// compute interval for triagle 1
if(COMPUTE_INTERVALS(vp0, vp1, vp2, dv0, dv1, dv2, dv0dv1, dv0dv2, t1)){
g_odata[tri1*TRI_NUM2+tri2] = coplanar(vn, s_v0, s_v1, s_v2, s_u0, s_u1, s_u2, v, u);
}
if(g_odata[tri1*TRI_NUM2+tri2] == 2) // 이 위치 맞나?
{
// compute interval for triangle 2
if(COMPUTE_INTERVALS(up0, up1, up2, du0, du1, du2, du0du1, du0du2, t2)){
g_odata[tri1*TRI_NUM2+tri2] = coplanar(vn, s_v0, s_v1, s_v2, s_u0, s_u1, s_u2, v, u);
}
if(g_odata[tri1*TRI_NUM2+tri2] == 2) // 이 위치 맞나?
{
// sort t1
if(t1[0]>t1[1]) {
temp = t1[0];
t1[0] = t1[1];
t1[1] = temp;
}
// sort t2
if(t2[0]>t2[1]) {
temp = t2[0];
t2[0] = t2[1];
t2[1] = temp;
}
if(t1[1] < t2[0] || t2[1] < t1[0]){
g_odata[tri1*TRI_NUM2+tri2] = 0;
}
else{
g_odata[tri1*TRI_NUM2+tri2] = 1;
}
} // if ended // 이 위치 맞나?
} // if ended // 이 위치 맞나?
} // if ended
} // if ended
} // kernel ended
endif // #ifndef CPP_INTEGRATION_KERNEL_H