volume_division.br

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

#include <stdio.h>
#include <stdlib.h>
#include "ppm3d.h"
#include "volume_division.h"
char use_vout_filter=1;
char use_vout_amplify=1;

/**This kernel tests if the data provided is above the cutoff of zero
 * If it is, it returns the single bit passed into the index (ex 1,2,4,...,128)
 */
kernel float test(float vol<>,float index<>){
  return (vol>.39)?index:0;
}
kernel float test4(float4 vol<>,float4 index<>){
  float4 zero4 = 0;
  float4 isovalue=.39;
  float4 tmp = (vol>isovalue)?index:zero4;
  return dot(tmp,1);
}
/* This kernel is for a synthetic dataset consisting of one sphere centered at
 * offset and radius equal to sqrt(density). 
 * Not used for normal datasets.
 */
kernel float synthEval(float3 loc<>) {
  float  power;
  const float3 offset=1;
  float3 center =loc-offset;
  power = dot(center,center);

  //uncomment following to enable 8 spheres. Too many for nCidia or ATI
  /* 
  {
  const float length=20;
  center.x +=length;
  power += dot(center,center);
  center.y+=length;
  power += dot(center,center);
  center.x-= length;
  power +=  dot(center,center);
  center.z+= length;
  power += dot(center,center);
  
  center.y-=length;
  power += dot(center,center);
  
  center.x+=length;                  
  power += dot(center,center);
  center.y+=length;
  power += dot(center,center);
  }
  */
  {

     const float  density =.5*.5;
    //printf (power,loc.z,density,power-density);
     return 1;

    return density-power;
  }
}

/** This kernel is for the synthetic sphere dataset and evaluates
 *  The set of spheres for each of the neighbors and tests them > 0
 *  for each neighbor.  The final values are added together to produce a lookup
 *  index */
kernel float synthEvaluateNeighbors (float2 center<>,
                                float2 opposing<>,
                                float2 slice) {

   float a=test(synthEval(float3(center.x,center.y,slice.x)),1.0f) +
      test(synthEval(float3(center.x,opposing.y,slice.x)),2.0f) +
      test(synthEval(float3(opposing.x,opposing.y,slice.x)),4.0f) +
      test(synthEval(float3(opposing.x,center.y,slice.x)),8.0f) +
      test(synthEval(float3(center.x,center.y,slice.y)),16.0f) +
      test(synthEval(float3(center.x,opposing.y,slice.y)),32.0f) +
      test(synthEval(float3(opposing.x,opposing.y,slice.y)),64.0f) +
      test(synthEval(float3(opposing.x,center.y,slice.y)),128.0f);
   return a>0.5f&&a<(2.0f+4.0f+8.0f+16.0f+32.0f+64.0f+128.0f+0.5f)?a:0;
    
}
/** This kernel is for the synthetic sphere dataset and evaluates
 *  the neighbor values, and gets the lookup index for the triangle dataset
 *  Only pushes a result if the index is not 0 (empty case)
 **/
kernel void processSyntheticSlice(vout[1] float4 vertex<>,
                                iter float2 center<>,
                                iter float2 opposing<>,
                                float2 slice/*first value is cur then next*/) {
  float pattern;
  if((pattern=synthEvaluateNeighbors(center,opposing,slice))) {
    vertex=float4(center.x,center.y,slice.x,pattern);
    push(vertex);
  }
}
/** This kernel is for the synthetic sphere dataset and evaluates
 *  the neighbor values, and gets the lookup index for the triangle dataset
 *  It pushes a negative index in the event that all or none of the neighbors
 *  are dense enough */
kernel void processSyntheticSliceNoCompact(out float4 vertex<>,
                                           iter float2 center<>,
                                           iter float2 opposing<>,
                                           float2 slice
                                           /*first value is cur then next*/) {
  float pattern;
  if((pattern=synthEvaluateNeighbors(center,opposing,slice))) {
    vertex=float4(center.x,center.y,slice.x,pattern);
  }else {
    vertex=float4(-1.0,-1.0f,-1.0f,-1.0f);
  }
}

/** This kernel is for the experimental dataset and evaluates
 *  the neighbor values, and gets the lookup index for the triangle dataset
 *  Only pushes a result if the index is not 0 or 255 (empty cases)
 **/
kernel void processSlice (float curgather[][],
                          float nextslice[][],
                          vout [1]float4 vertex<>,
                          iter float2 center<>,
                          iter float2 up<>,
                          iter float2 forward<>,
                          iter float2 upforward<>,
                          float2 slice) {
  float4 o248={1.0f,2.0f,4.0f,8.0f};
  float4 t456={16.0f,32.0f,64.0f,128.0f};
  float4 cur={curgather[center],curgather[up],curgather[upforward],curgather[forward]};
  float4 next={nextslice[center],nextslice[up],nextslice[upforward],nextslice[forward]};
  float a=
    test4(cur,o248)+
    test4(next,t456);
  if (a>0.5f&&a<254.5f) {
    vertex=float4(center.x,center.y,slice.x,a);
    push(vertex);
  }
}

/** This kernel is for the synthetic sphere dataset and evaluates
 *  the neighbor values, and gets the lookup index for the triangle dataset
 *  It pushes a negative index in the event that all or none of the neighbors
 *  are dense enough */
kernel void processSliceNoCompact (float curgather[][],
                                   float nextslice[][],
                                   out float4 vertex<>,
                                   iter float2 center<>,
                                   iter float2 up<>,
                                   iter float2 forward<>,
                                   iter float2 upforward<>,
                                   float2 slice) {
  float a=
    test(curgather[center],1.0f)+
    test(curgather[up],2.0f)+
    test(curgather[upforward],4.0f)+
    test(curgather[forward],8.0f)+
    test(nextslice[center],16.0f)+
    test(nextslice[up],32.0f)+
    test(nextslice[upforward],64.0f)+
    test(nextslice[forward],128.0);
  if (a>0.5f&&a<(2.0f+4.0f+8.0f+16.0f+32.0f+64.0f+128.0f+0.5f)) {
    vertex=float4(center.x,center.y,slice.x,a);
  }else {
    vertex=float4(-1,-1,-1,-1);
  }
  
}
                           
/** This kernel outputs exactly 5 triangles for each input
 *  lookup value. It assumes the sentinel is immutable with respect to 
 *  addition (i.e. inf behaves this way)
 *  Given that up to five triangles are produced for each volumetric datapoint
 *  exactly five triangles are outputted reguardless.
 *  The vertex is determined based on the index mod 3... the output triangles
 *  are stretched to 3x the input so that each of the 3 requisite vertices
 *  are produced per output.
 */
kernel void processTrianglesNoCompact(out float3 trianglesA<>, 
                                      out float3 trianglesB<>, 
                                      out float3 trianglesC<>,
                                      out float3 trianglesD<>, 
                                      out float3 trianglesE<>, 
                                      float4 vertices1<>,
                                      float3 volumeTriangles[][]) {
  float4 whichVolumeTriangle={round(fmod((indexof trianglesA).x,3)),vertices1.w,0,0};
  float3 vertices=vertices1.xyz;
  if (!(whichVolumeTriangle.y<256&&whichVolumeTriangle.y>=0))
    whichVolumeTriangle.y=0;
  if (whichVolumeTriangle.x+.5>3) whichVolumeTriangle.x=0;
  trianglesA=vertices.xyz+volumeTriangles[whichVolumeTriangle];
  whichVolumeTriangle.x+=3;
  trianglesB=vertices.xyz+volumeTriangles[whichVolumeTriangle];
  whichVolumeTriangle.x+=3;
  trianglesC=vertices.xyz+volumeTriangles[whichVolumeTriangle];
  whichVolumeTriangle.x+=3;  
  trianglesD=vertices.xyz+volumeTriangles[whichVolumeTriangle];  
  whichVolumeTriangle.x+=3;  
  trianglesE=vertices.xyz+volumeTriangles[whichVolumeTriangle];  
}

kernel float4 interp (float4 a, float4 b) {
  return lerp (0,1,1-(.39-a)/(b-a));  
}
kernel void processTrianglesNoCompactOneOut(out float3 trianglesA<>, 
                                            float4 verticesinf<>,
                                            float3 volumeTriangles[][],
                                            float curgather[][],
                                            float nextgather[][]) {
  float isovalue=.39;
  float2 threefive={3.0f,5.0f};
  float2 onehalf2=.5;
  float2 zero2=0;
  float4 four256=256;
  float4 vertices1=(verticesinf<four256&&verticesinf>=zero2.xxxx)
    ? verticesinf
    : zero2.xxxx;
  float2 whichVolumeTriangle={0,vertices1.w};
  float3 vertices=vertices1.xyz;
  float2 mod35=round(fmod((indexof trianglesA).xy,threefive));
  mod35 = (mod35+onehalf2>threefive)?zero2:mod35;
  whichVolumeTriangle.x=mod35.x+mod35.y*3; 
  {
    float3 configinf = volumeTriangles[whichVolumeTriangle].xyz;
    float3 config = configinf<four256.xxx?configinf:zero2.xxx;
    float low,high;
    float3 ishalf,tmp;
    ishalf=(float3)config==onehalf2.x;
    if (ishalf.z) {
       float2 lhGather = vertices.xy+config.xy;
       low = curgather[lhGather];
       high = nextgather[lhGather];
    }else {
       float2 lowGather,highGather;
       if (ishalf.x) {
          lowGather=float2(vertices.x,vertices.y+config.y);
          highGather=float2(vertices.x+1.0f,vertices.y+config.y);
       } else {
          lowGather=float2(vertices.x+config.x,vertices.y);
          highGather=float2(vertices.x+config.x,vertices.y+1.0f);
       }
       if (configinf.z==0) {
          low = curgather[lowGather];
          high = curgather[highGather];
       }else {
          low = nextgather[lowGather];
          high= nextgather[highGather];
       }

    }
    tmp = (isovalue-low)/(high-low);
    config=ishalf?tmp:configinf;
    trianglesA=vertices+config;
  }
   
}/*
kernel void processTrianglesNoCompactOneOut(out float3 trianglesA<>, 
                                      float4 vertices1<>,
                                      float3 volumeTriangles[][]) {
  float4 whichVolumeTriangle={round(fmod((indexof trianglesA).x,15)),vertices1.w,0,0};
  float3 vertices=vertices1.xyz;
  if (!(whichVolumeTriangle.y<256&&whichVolumeTriangle.y>=0))
    whichVolumeTriangle.y=0;
  if (whichVolumeTriangle.x+.5>15) whichVolumeTriangle.x=0;
  trianglesA=vertices.xyz+volumeTriangles[whichVolumeTriangle];
  }
*/

/** This kernel outputs exactly 5 triangles for each input
 *  lookup value. It assumes the sentinel is not between -.75 and .75
 *  Given that up to five triangles are produced for each volumetric datapoint
 *  The output stream is stretched to 4x the input stream size, and either 0 or
 *  3 vertices are produced for the given volumeTriangle lookup.
 */


kernel void processTriangles(vout[1] float4 triangles<>, 
                             float4 vertices<>,
                             float3 volumeTriangles[][],
                             iter float2 streamsize<>) {

 float4 whichVolumeTriangle={fmod(streamsize.y,5.0f)*3.0f,vertices.w,0,0};
 float3 firstTrianglePos;
 if (whichVolumeTriangle.x+.5>15.0f) whichVolumeTriangle.x=3.0f;
 if (!(whichVolumeTriangle.y<256&&whichVolumeTriangle.y>=0))
   whichVolumeTriangle.y=0;
 firstTrianglePos=volumeTriangles[whichVolumeTriangle];
 if (abs(firstTrianglePos.x-.5)<.75) {
   triangles.xy=(indexof vertices).xy;
   triangles.zw = whichVolumeTriangle.xy;
   push(triangles);
 }  
}
/**
 * Since every one of the outputted volumeTriangle lookup indices has at least 
 * one triangle, the first triangle is a guaranteed hit and need not be 
 * conditionally produced. 
 * Hence the triangle is looked up, and produced. The respective vertex is
 * computed depending on the fmod with 3
 */
kernel void processFirstTriangles(out float3 triangles<>, 
                                  float4 vertices<>,
                                  float3 volumeTriangles[][]) {

    float4 whichVolumeTriangle={round(fmod((indexof triangles).x,3.0f)),vertices.w,0,0};
   if (whichVolumeTriangle.x+.5>3.0f) whichVolumeTriangle.x=0.0f;
   if (!(whichVolumeTriangle.y<256&&whichVolumeTriangle.y>=0))
      whichVolumeTriangle.y=0;
   triangles=vertices.xyz;
     //   if (whichVolumeTriangle.y<256&&whichVolumeTriangle.y>=0)
   triangles+=volumeTriangles[whichVolumeTriangle];
}

typedef float3 Triangle[5][3];
/**
 * The following function generates the triangle lookup table from the raw
 * edge list data stored in volume_division.h
 * First it fills everything with the sentinel.
 * Then it goes through up to 15 items, filling the triangle table until a 
 * -1 in the edge list is encountered at which point it stops adding to the tex
 */
char volumeTriangles(Triangle tri[256]) {
  unsigned int i,j,k,m;

   for (i=0;i<256;++i) {
      for (m=0;m<15;++m) {
        // fill with sentinel
         tri[i][(m/3)%5][m%3]=float3(1.0f/(float)floor(.5),
                                     1.0f/(float)floor(.5),
                                     1.0f/(float)floor(.5));
      }
      j=0;
      for (j=0;;j+=3) {
         if (m_triTable[i][j]==-1){// our edge list is -1, stop!
           m_triNum[i]=(float)(j/3);//set the count of triangles to j
           break; // it's too dangerous we have to slow down first!
         }
         for (k=0;k<3;++k) {
           // single out the point that will be set on this pass
            float3 * p=&tri[i][j/3][k];
            p->x=m_triTable[i][j+k];// for debugging only, overwritten
            p->y=m_triTable[i][j+k];// for debugging only
            p->z=m_triTable[i][j+k];// for debugging only
            switch ((int)m_triTable[i][j+k]) {
              // set the appropriate coordinates in a 1x1x1 cube
              // given the lookup value assuming the data is exactly
              // in the middle--i.e. no adaptive
            case 0:
               p->x=p->z=0;p->y=0.5f;
               break;
            case 1:
               p->x=0.5f;p->z=0;p->y=1.0f;
               break;
            case 2:
               p->x=1.0f;p->z=0;p->y=0.5f;
               break;
            case 3:
               p->x=0.5f;p->z=0;p->y=0;
               break;
            case 4:
               p->x=0;p->z=1.0f;p->y=0.5f;
               break;
            case 5: