subdiv.br

用于GPU通用计算的编程语言BrookGPU 0.4

#include <stdio.h>
#include <stdlib.h>
#include "timing.h"
/**   ABB__BBC 
      /\  /\
   AB/__\/__\BC
    /\  B\  /\
AAB/__\/__\/__\BCC
   \  A\  /C  /
    \/__\/__\/
   AAC  AC  ACC
**/
typedef struct STri_t{
   float4 A; // the last value of A indicates whether the edge AB is small
   float4 B; // enough to stop subdividing.  B.w also indicates if BC is small
   float4 C; // enough. C.w indicates if AC is small enough to stop subdividing
} STri;

//When a triangle is split, the new triangle really only has six unique new
//values.  These are copied into the final 12 new triangle slots.
typedef struct SplitTri_t {
  float4 A;
  float4 AB;
  float4 B;
  float4 BC;
  float4 C;
  float4 AC;
}SplitTri;
//Stores the neighbors of a given triangle The w component indicates if
//the line length is short enough for the triangle to be produced.
//this could be altered to any per-edge qualifier
// for the intermediate split neighbors the unused 'w' components act as holders
// for the recomputed neighbor list when a triangle is split
typedef struct Neighbor_t {
  float4 AB;// w = AB->B.x
  float4 BBC;// w = AB->B.y
  float4 ABB;// w = AB->B.z
  float4 BC;// w = BC->C.x
  float4 ACC;// w = BC->C.y
  float4 BCC;// w = BC->C.z
  float4 AC;// w = AC->A.x
  float4 AAB;// w = AC->A.y
  float4 AAC;// w = AC->A.z
}Neighbor;

unsigned char onlyCracks=0;
unsigned char noCracks=0;
unsigned char debugLoop=0;
int counterMax = 1024;
extern float neighboreps;
int low_texture_ram=0;
int subdivisiondepth=0;
#ifdef _WIN32
const unsigned int maxTexDim=2048;
const unsigned int maxTexDimLog2=11;
#else
#define maxTexDime maxTexDim
#define maxTexDimeLog2 maxTexDimLog2
const unsigned int maxTexDime=1024;
const unsigned int maxTexDimeLog2=10;

#endif
void recomputeNeighbors(STri*, Neighbor *, unsigned int);
void checkNeighbors(STri * tri, Neighbor * in, Neighbor * out, unsigned int);
void __printf_cpu_inner(float f,float g, float h) {
  fprintf (stderr,"Val %f %f %f\n",f,g,h);
}
void __emer_cpu_inner(float f,float g, float h) {
  fprintf (stderr,"Emer %f %f %f\n",f,g,h);
}

//returns a 2 if only b1, 4 if only b2 and 8 if both
//this acts as sort of a bitwise holder class that should work with 24 bitfloat
kernel float triCombine(float b1<>, float b2<>) {
   return (b1&&b2)?8:(b1?2:(b2?4:0));
}

kernel float uoDistance (float3 a<>, float3 b<>){
   float3 bigger = a>b?a:b;
   float3 smaller = (a>b)?b:a;
   float3 bs = bigger-smaller;
   return dot(bs,bs);
}
kernel float3 avg (float3 a<>,float3 b<>) {
   float3 bigger = a>b?a:b;
   float3 smaller = a>b?b:a;
   return .5*(smaller+bigger);
   //return lerp(a,b,.5);
}
kernel float4 avg4 (float4 a<>,float4 b<>) {
   float4 bigger = a>b?a:b;
   float4 smaller = a>b?b:a;
   return .5*(smaller+bigger);
   //return lerp(a,b,.5);
}
//if the model is missing a neighbor then all 3 components will be zero
//i.e. if a model has 4 or 5 neighbors only
kernel float isNeighbor(float4 neigh) {
   return dot(neigh.xyz,neigh.xyz)>0;
}
kernel float isNotNeighbor(float4 neigh) {
  return !(dot(neigh.xyz,neigh.xyz)>0);
}
kernel float3 mymax (float3 a, float3 b, float3 c) {
   return a>b?(a>c?a:c):(b>c?b:c);
}
kernel float3 mymin (float3 a, float3 b, float3 c) {
   return a>b?(b>c?c:b):(a>c?c:a);
}
kernel float3 mymed (float3 a, float3 b, float3 c,float3 max, float3 min) {
   return a!=max&&a!=min?a:(((b!=max&&b!=min)||a==b)?b:c);
}
kernel float3 min3(float3 a, float3 b) {
   return a>b?b:a;
}
kernel float3 max3 (float3 ain, float3 bin) {
   return ain<bin?bin:ain;
}

kernel float3 swap3(float3 ain<>, float3 bin<>, out float3 bout<>) {
   float3 ret = min3(ain,bin);
   bout = max3(ain,bin);
   return ret;
}

kernel float3 add6(float3 a<>,float3 b<>, float3 c<>, float3 d<>, float3 e<>, float3 f<>){
   float3 ad = min3(a,d)+max3(a,d);
   float3 be = min3(b,e)+max3(b,e);
   float3 cf = min3(c,f)+max3(c,f);
  float3 mint = mymin(ad,be,cf);
  float3 maxt = mymax(ad,be,cf);
  float3 medt = mymed(ad,be,cf,maxt,mint);
  medt = mint+medt;
  return medt+maxt;
  //return a+b+c+d+e+f;
}
kernel void add6test (float3 a<>,float3 b<>, float3 c<>, float3 d<>, float3 e<>, float3 f<>, out float3 ooo<>) {
   ooo= add6(a,b,c,d,e,f);
}
//The function that recomputes whether the given edges are large enough
//to warrant a split. This could and probably should be done in "image space"
kernel void smallEnough(STri t<>, Neighbor u<>, 
                        out STri v<>, 
                        out Neighbor n<>, 
                        float epsilon) {
  v.A.xyz = t.A.xyz;
  v.A.w = (uoDistance(t.A.xyz,t.B.xyz)>epsilon)?1.0:0.0;
  v.B.xyz = t.B.xyz;
  v.B.w = (uoDistance(t.B.xyz,t.C.xyz)>epsilon)?1.0:0.0;
  v.C.xyz = t.C.xyz;
  v.C.w = (uoDistance(t.A.xyz,t.C.xyz)>epsilon)?1.0:0.0;
  n.AB.xyz = u.AB.xyz;
  //for the neighbors that two vertices share, there are two relevant lenghts,
  // both of which must be saved in the neighbor .w
  n.AB.w = triCombine((uoDistance(t.A.xyz,u.AB.xyz)>epsilon)?1.0:0.0,
                      (uoDistance(t.B.xyz,u.AB.xyz)>epsilon)?1.0:0.0);
  n.BC.xyz=u.BC.xyz;
  n.BC.w = triCombine((uoDistance(t.B.xyz,u.BC.xyz)>epsilon)?1.0:0.0,
                      (uoDistance(t.C.xyz,u.BC.xyz)>epsilon)?1.0:0.0);
  n.AC.xyz=u.AC.xyz;
  n.AC.w = triCombine((uoDistance(t.C.xyz,u.AC.xyz)>epsilon)?1.0:0.0,
                      (uoDistance(t.A.xyz,u.AC.xyz)>epsilon)?1.0:0.0);
  n.AAC.xyz=u.AAC.xyz;
  n.AAC.w = (isNotNeighbor(u.AAC)|| //if missing neighbor, should not stop split
             uoDistance(t.A.xyz,u.AAC.xyz)>epsilon)?1.0:0;
  n.AAB.xyz=u.AAB.xyz;
  n.AAB.w = (isNotNeighbor(u.AAB)|| //future split should not stop upon missing neigh
             uoDistance(t.A.xyz,u.AAB.xyz)>epsilon)?1.0:0;
  n.ABB.xyz=u.ABB.xyz;
  n.ABB.w = (isNotNeighbor(u.ABB)|| //if missing neighbor, should not stop split
             uoDistance(t.B.xyz,u.ABB.xyz)>epsilon)?1.0:0;
  n.BBC.xyz=u.BBC.xyz;
  n.BBC.w = (isNotNeighbor(u.BBC)|| //if missing neighbor, should not stop split
             uoDistance(t.B.xyz,u.BBC.xyz)>epsilon)?1.0:0;
  n.BCC.xyz=u.BCC.xyz;
  n.BCC.w = (isNotNeighbor(u.BCC)|| //if missing neighbor, should not stop split
             uoDistance(t.C.xyz,u.BCC.xyz)>epsilon)?1.0:0;
  n.ACC.xyz=u.ACC.xyz;
  n.ACC.w = (isNotNeighbor(u.ACC)|| //if missing neighbor, should not stop split
             uoDistance(t.C.xyz,u.ACC.xyz)>epsilon)?1.0:0;
}
kernel void produceTriP(STri t<>, vout [1]float2 shouldProduce<>, vout[1] float2 shouldNotProduce<>){
   if (t.A.w!=0.0f
       ||t.B.w!=0.0f
       ||t.C.w!=0.0f) { // if either of the 3 vertices are still mobile, continue to split
     shouldProduce = (indexof t).xy;
     push(shouldProduce);
   }else { //all the vertices are pinned by edges (either neighbor or triangle) that are too short
     shouldNotProduce = (indexof t).xy;
     push(shouldNotProduce);
   }
}

kernel Neighbor computeNeighbors(STri t<>, 
                             Neighbor u<>, 
                             float noAdaptive) {
  Neighbor v;
  float eAAB=isNeighbor(u.AAB);
  float eAAC=isNeighbor(u.AAC);
  float eABB=isNeighbor(u.ABB);
  float eBBC=isNeighbor(u.BBC);
  float eBCC=isNeighbor(u.BCC);
  float eACC=isNeighbor(u.ACC);
  float3 AB_B, AC_A, BC_C;
  float3 zero3=0;
  AB_B=avg(t.B.xyz,u.AB.xyz);
  if(1)if(noAdaptive||u.AB.w==4||u.AB.w==8)
     AB_B=lerp(AB_B,
               avg(t.A.xyz,eABB?u.ABB.xyz:eBBC?u.BBC.xyz:u.BC.xyz),
               .25);
  AC_A=avg(t.A.xyz,u.AC.xyz);
  if(1)if(noAdaptive||u.AC.w==4||u.AC.w==8)
     AC_A=lerp(AC_A,
               avg(t.C.xyz,eAAC?u.AAC.xyz:eAAB?u.AAB.xyz:u.AB.xyz),
               .25);
  BC_C=avg(t.C.xyz,u.BC.xyz);
  if(1)if(noAdaptive||u.BC.w==4||u.BC.w==8)
     BC_C=lerp(BC_C,
               avg(t.B.xyz,eBCC?u.BCC.xyz:eACC?u.ACC.xyz:u.AC.xyz),
               .25);

  v.AB.xyz = avg(t.A.xyz,u.AB.xyz);
  if(1)if(noAdaptive||u.AB.w==2||u.AB.w==8)
     v.AB.xyz = lerp(v.AB.xyz,
                     avg(t.B.xyz,eAAB?u.AAB.xyz:eAAC?u.AAC.xyz:u.AC.xyz),
                     .25);
  v.AB.w = AB_B.x;
  v.BBC.xyz = eBBC?avg(t.B.xyz,u.BBC.xyz):zero3;
  if(1)if(eBBC&&(noAdaptive||u.BBC.w))
     v.BBC.xyz = lerp(v.BBC.xyz,
                      avg(eABB?u.ABB.xyz:u.AB.xyz,u.BC.xyz),
                      .25);
  v.BBC.w = AB_B.y;
  v.ABB.xyz = eABB?avg(t.B.xyz,u.ABB.xyz):zero3;
  if(1)if(eABB&&(noAdaptive||u.ABB.w))
     v.ABB.xyz = lerp(v.ABB.xyz,
                      avg(u.AB.xyz,eBBC?u.BBC.xyz:u.BC.xyz),
                      .25);
  v.ABB.w = AB_B.z;
  v.AC.xyz = avg(t.C.xyz,u.AC.xyz);
  if(1)if(noAdaptive||u.AC.w==2||u.AC.w==8)
     v.AC.xyz = lerp(v.AC.xyz,
                     avg(t.A.xyz,eACC?u.ACC.xyz:eBCC?u.BCC.xyz:u.BC.xyz),
                     .25);
  v.AC.w = AC_A.x;
  v.AAB.xyz = eAAB?avg(t.A.xyz,u.AAB.xyz):zero3;
  if(1)if(eAAB&&(noAdaptive||u.AAB.w))
     v.AAB.xyz = lerp(v.AAB.xyz,
                      avg(eAAC?u.AAC.xyz:u.AC.xyz,u.AB.xyz),
                      .25);
  v.AAB.w = AC_A.y;
  v.AAC.xyz = eAAC?avg(t.A.xyz,u.AAC.xyz):zero3;
  if(1)if(eAAC&&(noAdaptive||u.AAC.w))
     v.AAC.xyz = lerp(v.AAC.xyz,
                      avg(eAAB?u.AAB.xyz:u.AB.xyz,u.AC.xyz),
                      .25);
  v.AAC.w = AC_A.z;
  v.BC.xyz = avg(t.B.xyz,u.BC.xyz);
  if(1)if(noAdaptive||u.BC.w==2||u.BC.w==8)
     v.BC.xyz = lerp(v.BC.xyz,
                     avg(t.C.xyz,eBBC?u.BBC.xyz:eABB?u.ABB.xyz:u.AB.xyz),
                     .25);
  v.BC.w = BC_C.x;
  v.ACC.xyz = eACC?avg(t.C.xyz,u.ACC.xyz):zero3;
  if(1)if(eACC&&(noAdaptive||u.ACC.w))
     v.ACC.xyz = lerp(v.ACC.xyz,
                      avg(u.AC.xyz,eBCC?u.BCC.xyz:u.AB.xyz),
                      .25);
  v.ACC.w = BC_C.y;
  v.BCC.xyz = eBCC?avg(t.C.xyz,u.BCC.xyz):zero3;
  if(1)if(eBCC&&(noAdaptive||u.BCC.w))
     v.BCC.xyz = lerp(v.BCC.xyz,
                      avg(u.BC.xyz,eACC?u.ACC.xyz:u.AC.xyz),
                      .25);
  v.BCC.w = BC_C.z;  
  return v;
}

kernel void computeNeighborsNoAdaptive(STri tgather<>, 
                                       Neighbor ugather<>, 
                                       out Neighbor v<>) {
  v=computeNeighbors(tgather,ugather,1);
}
kernel void computeNeighborsAdaptive(STri tgather[][], 
                             Neighbor ugather[][], 
                             float2 i<>, 
                             out Neighbor v<>) {
  float2 zero2=0;
  float2 ind = (i>=0&&i<8192)?i:zero2;
  Neighbor u = ugather[ind];
  STri t = tgather[ind];
  v=computeNeighbors(t,u,0);
}

//Splay the triangles into a vertex list of float3's.
//these are the final triangles
kernel void copyFinalTriangles (STri tri<>,
                                 out float3 triList<>) {
  float whichVertex=round(fmod((indexof triList).x,3));
  triList = tri.A.xyz;
  if (whichVertex>.5&&whichVertex<1.5)
    triList=tri.B.xyz;
  else if (whichVertex>=1.5&&whichVertex<2.5)
    triList=tri.C.xyz;
}
//Same as copy final triangles but with a dependent texture lookup
kernel void writeFinalTriangles (STri triangles[][],
                                 float2 indices<>,
                                 out float3 triList<>) {
  float2 zero2=0;
  float2 ind = (indices>=0&&indices<8192)?indices:zero2;
  copyFinalTriangles(triangles[ind],triList);
}
//Same as copy final triangles but with a dependent texture lookup
kernel void gatherTriangles (STri triangles[][],
                             float2 indices<>,
                             out STri triList<>) {
  float2 zero2=0;
  float2 ind = (indices>=0&&indices<8192)?indices:zero2;
  triList=triangles[ind];
}

//since we don't have proper float4 constructors, we need to manufacture one
kernel float4 addN(float3 inxyz, float inw) {
   float4 x = float4(inxyz.x,
                     inxyz.y,
                     inxyz.z,
                     inw);
   return x;
} 
//nothing happens
kernel float4 identity(float4 a, float b){return a;}
//you fail to obtain anything
kernel float3 ident3(float3 a, float b) {return a;}


kernel void repairTJunctions(STri triangles[][],
                             float3 where<>,
                             out float3  tris<>) {
  float2 zero2=0;
  float2 index = (where.xy>=0&&where.xy<65536)?where.xy:zero2;
  STri tri = triangles[index.xy];
  float mod3 = round(fmod((indexof tris).x,3.0f));
  float3 AC = avg (tri.A.xyz,tri.C.xyz);
  float3 AB = avg (tri.A.xyz,tri.B.xyz);
  float3 BC = avg (tri.B.xyz,tri.C.xyz);  
  if (mod3>2.5||mod3<.5) {
    if (where.z==0) {
      tris=tri.A.xyz;
    }else if (where.z==1) {
      tris=tri.B.xyz;
    }else {
      tris=tri.C.xyz;
    }
  }else if (mod3<1.5) {
    if (where.z==0) {
      tris=tri.B.xyz;
    }else if (where.z==1) {
      tris=tri.C.xyz;
    }else {
      tris=tri.A.xyz;
    }    
  }else {
    if (where.z==0) {
      tris=AB;
    }else if (where.z==1) {
      tris=BC;
    }else {
      tris=AC;
    }    
  }
  
  
}
kernel void identifyTJunctions(STri triangles[][],
                             Neighbor neighbors[][],
                             float2 i<>,
                             vout[1] float3 cracksA<>,
                             vout[1] float3 cracksB<>) {
  float2 zero2=0;
  float2 index = (i>=0&&i<65536)?i:zero2;
  STri tri=triangles[index];
  Neighbor oldneighbors=neighbors[index];