subdiv.br

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

      Neighbor * neigh;
      Neighbor * checknei;
      unsigned int totalsize;
      totalsize = ((unsigned int)streamSize(*triangles).x)*(unsigned int)streamSize(*triangles).y;
      tri = (STri*)malloc (sizeof(STri)*totalsize);
      neigh = (Neighbor*)malloc (sizeof(Neighbor)*totalsize);
      checknei = (Neighbor*)malloc (sizeof(Neighbor)*totalsize);
      streamWrite(*neighbors,neigh);
      streamWrite(*triangles,tri);
      recomputeNeighbors(tri,checknei,totalsize);
      checkNeighbors(tri,neigh,checknei,totalsize);
      free (tri);
      free (neigh);
      free(checknei);
    }
    if (counter<counterMax) {
      //if we can afford the ram, we recursively call this function
      // to evaluate other triangles
      subdivide(neighbors,triangles,output,output4,counter+1,out4);
    } 
  }
}


//the main controller function...takes in neighbors and triangles
//Note smallEnough must have been run on the input for this to work
//small enough is run from main or from within linearReorgSplitTriangles
void subdivideAdaptive (Neighbor (*neighbors)<>,
                        STri (*triangles)<>,
                        struct VertexArray *output,
                        struct VertexArray4 *output4,
                        float epsilon,
                        int counter, char out4) {
  float3 tJunctionsRepairA<0,0>;
  float3 tJunctionsRepairB<0,0>;
  unsigned int streamY=toi(streamSize(*triangles).y);
  unsigned int streamX=toi(streamSize(*triangles).x);
  unsigned char snpSize=0;
  subdivisiondepth++;
  {
    float2 shouldProduce<streamY,0>;
    float2 shouldNotProduce<streamY,0>;
    //see whether triangles should be split

    produceTriP(*triangles,shouldProduce,shouldNotProduce);
    TallyKernel("produceTriP",*triangles,shouldProduce,shouldNotProduce);
    assert((int)streamX==toi(streamSize(shouldProduce).x));
    {
      int wosizex=3*toi(streamSize(shouldNotProduce).x);
      int wosizey=toi(streamSize(shouldNotProduce).y);
      int tsizex=toi(streamSize(shouldProduce).x);
      int tsizey=toi(streamSize(shouldProduce).y);
      float3 outputTri<(out4?0:wosizey),(out4?0:wosizex)>;
      STri outputTri4<(out4?wosizey:0),(out4?toi(streamSize(shouldNotProduce).x):0)>;
      if (tsizey&&tsizex) {
        float3 tJunctionsA<tsizey,0>;
        float3 tJunctionsB<tsizey,0>;
        identifyTJunctions(*triangles,*neighbors,shouldProduce,tJunctionsA,tJunctionsB);
        TallyKernel("identifyTJunctions",shouldProduce,tJunctionsA,tJunctionsB);
        if (streamSize(tJunctionsA).x&&streamSize(tJunctionsA).y) {
          int tsizex = toi(streamSize(tJunctionsA).x)*3;
          int tsizey = toi(streamSize(tJunctionsA).y);
          float3 outputTri<tsizey,tsizex>;
          streamSwap(outputTri,tJunctionsRepairA);
          repairTJunctions(*triangles,tJunctionsA,tJunctionsRepairA);
          TallyKernel("repairTJunctions",tJunctionsRepairA);
        }
        if (streamSize(tJunctionsB).x&&streamSize(tJunctionsB).y) {
          int tsizex = toi(streamSize(tJunctionsB).x)*3;
          int tsizey = toi(streamSize(tJunctionsB).y);
          float3 outputTri<tsizey,tsizex>;
          streamSwap(outputTri,tJunctionsRepairB);
          repairTJunctions(*triangles,tJunctionsB,tJunctionsRepairB);
          TallyKernel("repairTJunctions",tJunctionsRepairB);
        }
      }
      if (wosizex&&wosizey) {
        if (out4) {
          gatherTriangles(*triangles,shouldNotProduce,outputTri4);          
          TallyKernel("gatherTriangles",outputTri4);
        }else {
          writeFinalTriangles(*triangles,shouldNotProduce,outputTri);      
          TallyKernel("writeFinalTriangles",outputTri);
        }
      }

      streamY=toi(streamSize(shouldProduce).y);
      //if some should be split and produced then go here
      if (streamY&&streamX) {
        float stretchX = streamY*4>maxTexDim?(streamY*2>maxTexDim?1.0f:2.0f):0.0f;      
        int xfactor=stretchX==2.0f?2:(stretchX==1.0f?4:1);
        int yfactor=stretchX==2.0f?2:(stretchX==1.0f?1:4);
        //allocate shared neighbor structs that hold the extra generated floats
        Neighbor sharedNeighbors<streamY,streamX>;
        SplitTri sharedTriangles <streamY,streamX>;
        //these new neighbor structs are 4x as big as the old ones
        //input is thus stretched through them based on the fmod
        Neighbor newNeighbors<(streamY*yfactor),
          (streamX*xfactor)>;
        STri newTriangles<(streamY*yfactor),
          (streamX*xfactor)>;
        //computes the 12 shared neighbor values
        fprintf (stderr,"stretching %f to %d %d\n",stretchX,streamX*xfactor,streamY*yfactor);
        computeNeighborsAdaptive(*triangles,*neighbors,shouldProduce,sharedNeighbors);
        TallyKernel("computeNeighborsAdaptive",sharedNeighbors); 
        //computes the 6 new triangle vertex values
        splitTrianglesAdaptive(*triangles,*neighbors,shouldProduce,sharedTriangles);
        TallyKernel("splitTrianglesAdaptive",sharedTriangles);

        // combines the 12 neighbors and 6 vertices to produce 4 individual 
        // sets of triangles
        // with the new triangles getting stretched in the x or y direction depending
        // on stretchX
        adaptiveReorgSplitTriangles(sharedTriangles,
                                    sharedNeighbors,
                                    newTriangles,
                                    newNeighbors,
                                    epsilon,
                                    stretchX);
        TallyKernel("adaptiveReorgSplitTriangles",newTriangles);
        streamSwap(*triangles,newTriangles);
        streamSwap(*neighbors,newNeighbors);
        if (debugLoop) {
          STri * tri;
          Neighbor * neigh;
          Neighbor * checknei;
          unsigned int totalsize;
          totalsize = ((unsigned int)streamSize(*triangles).x)*(unsigned int)streamSize(*triangles).y;
          tri = (STri*)malloc (sizeof(STri)*totalsize);
          neigh = (Neighbor*)malloc (sizeof(Neighbor)*totalsize);
          checknei = (Neighbor*)malloc (sizeof(Neighbor)*totalsize);
          streamWrite(*neighbors,neigh);
          streamWrite(*triangles,tri);
          recomputeNeighbors(tri,checknei,totalsize);
          checkNeighbors(tri,neigh,checknei,totalsize);
          free (tri);
          free (neigh);
          free(checknei);
        }
      }
      // if we have some triangles that are set to be produced
      if (streamSize(shouldNotProduce).y){
        snpSize=1;
        // write them out to a vertex array
        if ((!low_texture_ram)&&
            streamSize(shouldProduce).y&&
            streamSize(shouldProduce).x&&
            counter<counterMax) {
          //if we can afford the ram, we recursively call this function
          // to evaluate other triangles
          subdivideAdaptive(neighbors,triangles,output,output4,epsilon,counter+1,out4);
        } 
        //now we writeback from the card to local memory since we dont have
        // render to vertex array yet
        if (!onlyCracks){
          if (out4) {
            wosizey = output4->size;
            expandVertexArray4(output4,(unsigned int)(streamSize(outputTri4).x*3*
                                                      streamSize(outputTri4).y));
            streamWrite(outputTri4,output4->v+wosizey);               
          }else {
            wosizey = output->size;

            expandVertexArray(output,(unsigned int)(streamSize(outputTri).x*
                                                    streamSize(outputTri).y));
            streamWrite(outputTri,output->v+wosizey);   
          }
        }
        if (streamSize(tJunctionsRepairA).x&&streamSize(tJunctionsRepairA).y&&!noCracks) {
          wosizey = output->size;
          expandVertexArray(output,(unsigned int)(streamSize(tJunctionsRepairA).x*
                                                  streamSize(tJunctionsRepairA).y));
          streamWrite(tJunctionsRepairA,output->v+wosizey);             
        }
        if (streamSize(tJunctionsRepairB).x&&streamSize(tJunctionsRepairB).y&&!noCracks) {
          wosizey = output->size;
          expandVertexArray(output,(unsigned int)(streamSize(tJunctionsRepairB).x*
                                                  streamSize(tJunctionsRepairB).y));
          streamWrite(tJunctionsRepairB,output->v+wosizey);             
        }
      }
    }
    if (counter>=counterMax) {
      if (out4) {
        int sizey = output4->size;      
        expandVertexArray4(output4,
                           3*toi(streamSize(*triangles).y)*toi(streamSize(*triangles).x));
        streamWrite(*triangles,output4->v+sizey);
      }else {
        int sizey = output->size;      
        float3 outputTri<toi(streamSize(*triangles).y),
                       toi(streamSize(*triangles).x*3)>;
        copyFinalTriangles(*triangles,outputTri);
        if (!onlyCracks){
          expandVertexArray(output,
                            3*toi(streamSize(*triangles).y)*toi(streamSize(*triangles).x));
          streamWrite(outputTri,output->v+sizey);
        }
      }
      if (streamSize(tJunctionsRepairA).x&&streamSize(tJunctionsRepairA).y&&!noCracks) {
        int sizey = output->size;
        expandVertexArray(output,(unsigned int)(streamSize(tJunctionsRepairA).x*
                                                streamSize(tJunctionsRepairA).y));
        streamWrite(tJunctionsRepairA,output->v+sizey);             
      }
      if (streamSize(tJunctionsRepairB).x&&streamSize(tJunctionsRepairB).y&&!noCracks) {
        int sizey = output->size;
        expandVertexArray(output,(unsigned int)(streamSize(tJunctionsRepairB).x*
                                                streamSize(tJunctionsRepairB).y));
        streamWrite(tJunctionsRepairB,output->v+sizey);             
      }
    }
    
  }
  // if we really dont have that kind of ram we will already have written out
  // vertices out to host ram, thus preserving memory. Only recurse if we have 
  // size in shouldSplit arrays :-)
  if ((low_texture_ram||snpSize==0)&&streamY&&streamX&&counter<counterMax) {
     subdivideAdaptive(neighbors,triangles,output,output4,epsilon,counter+1,out4);
  } 
  // prevents all writing at the end
}

extern unsigned int  loadModelData(const char * filename, 
                                   STri **tris, 
                                   Neighbor **nei);

STri *  createTriangles (unsigned int howMany) {
  STri * ret;  
  ret=(STri*)malloc(sizeof(STri)*howMany);
  memset(ret,0,sizeof(STri)*howMany);
  return ret;
}
Neighbor *  createNeighbors (STri * tri, unsigned int howMany) {
  Neighbor * ret;  
  ret=(Neighbor*)malloc(sizeof(Neighbor)*howMany);
  memset(ret,0,sizeof(Neighbor)*howMany);
  return ret;
}
void bestHeiWid(unsigned int totalSize,
                unsigned int *hei,
                unsigned int *wid) {
   unsigned int slop=totalSize*2;
   unsigned int maxy = 0;
   unsigned int best = maxTexDimLog2;
   unsigned int i;
   for (i=0;i<maxTexDimLog2+1;++i) {
      unsigned int x= (1<<i);
      unsigned int y = totalSize/x+(totalSize%x==0?0:1);
      unsigned int tslop;
      if (x>maxTexDim||y>maxTexDim)
         continue;
      tslop = x*y - totalSize;
      if (tslop<0)
         continue;
      if ((tslop<slop&&y>maxy)||
          (y>maxy&&tslop<2*y/3)) {
         *hei=y;
         *wid=x;
         maxy = y;
         slop = tslop;
         best = i;
      }
   }
}
unsigned int bestHeight(unsigned int totalSize){
   unsigned int hei=0,wid=0;
   bestHeiWid(totalSize,&hei,&wid);
   return hei;
}
unsigned int bestWidth(unsigned int totalSize){
   unsigned int hei=0,wid=0;
   bestHeiWid(totalSize,&hei,&wid);
   return wid;
}
extern void computeFunctionCallPattern(float epsilon,
                                       int argc, char ** argv, int numTri,STri*triangles,Neighbor *neigh);
unsigned int loadModelData(char* file,STri **triangles,Neighbor **neighbors);
int main (int argc, char**argv) {
  unsigned int i = 0;
  unsigned int vcount=0;
  float epsilon=1;
  STri *triangledata = 0;
  Neighbor *neighbordata = 0;
  unsigned int numTriangles;
  char adaptive=1;
  char out4=1;
  neighboreps=epsilon/1024.0f;
  SetupMillisTimer(argc,argv);
  for (i=0;i<(unsigned int)argc;++i) {
    char match=0;
    int j;
    if (strcmp(argv[i],"-debug")==0) {
      match=1;
      debugLoop=1;
    }
    if (strncmp(argv[i],"-epsilon",8)==0) {
      epsilon = (float)atof(argv[i]+8);
      match=1;
    }
    if (strncmp(argv[i],"-counter",8)==0) {
      counterMax = (int)atoi(argv[i]+8);
      match=1;
    }
    if (strncmp(argv[i],"-cracks",7)==0) {
      noCracks=1;
      match=1;
    }
    if (strncmp(argv[i],"-onlycracks",8)==0) {
      onlyCracks=1;
      match=1;
    }
    if (strncmp(argv[i],"-noadaptive",8)==0) {
      adaptive=0;
      match=1;
    }
    if (strcmp(argv[i],"-out3")==0) {
      out4=0;
      match=1;
    }
    if (strncmp(argv[i],"-eps",4)==0&&!match) {
      epsilon = (float)atof(argv[i]+4);
      match=1;
    }
    if (match) {
      for (j=i+1;j<argc;++j) argv[j-1]=argv[j];
      argc--;
      i--;
    }
  }
  epsilon*=epsilon;
  numTriangles = 
    loadModelData(argc>1?argv[1]:"puma.ply",&triangledata,&neighbordata);
  {
  Neighbor neighbors<bestHeight(numTriangles),bestWidth(numTriangles)>;
  STri triangles<bestHeight(numTriangles),bestWidth(numTriangles)>;
  struct VertexArray v;
  struct VertexArray4 v4;
  streamRead(triangles,triangledata);
  streamRead(neighbors,neighbordata);
  initVertexArray(&v);
  initVertexArray4(&v4);
  
  // initialize neighbors and triangles
  smallEnough(triangles,neighbors,triangles,neighbors,epsilon);
  if ( 0 ) {
       for(i=0; i<numTriangles; i++){
	    int j;
	    fprintf(stderr, "#\n");
	    for(j=0; j<9; j++){
		 float4 *temp;
		 temp = ((float4 *)(&neighbordata[i]))+j;
	    fprintf(stderr, "%f, %f, %f, %f\n", temp->x, temp->y, temp->z, temp->w);
	    }
       }
  }
  if ( adaptive )
       subdivideAdaptive(&neighbors,&triangles,&v,&v4,epsilon,0,out4); 
  else
  {
    subdivide(&neighbors,&triangles,&v,&v4,0,out4);
  }
 
  printf("%d\n", v.size+v4.size);
  for(i=0; i<v.size; i++){
    printf("%3.3f, %3.3f, %3.3f\n", v.v[i].x, v.v[i].y, v.v[i].z);
  }
  for(i=0; i<v4.size; i++){
    printf("%3.3f, %3.3f, %3.3f\n", v4.v[i].x, v4.v[i].y, v4.v[i].z);
  }
  vcount=numTriangles;
  for (i=0;i<(unsigned int)subdivisiondepth;++i) {
    vcount*=4;
  }
  CleanupMillisTimer();
  fprintf (stderr,"Num Rounds %d Adaptive Num Triangles %d Num Triangles %d\n",subdivisiondepth,v.size/3,vcount);
  return 0;
  }
}