grid.c

   1 /*
   2  * This program is free software: you can redistribute it and/or modify
   3  * it under the terms of the GNU General Public License as published by
   4  * the Free Software Foundation, either version 3 of the License, or
   5  * (at your option) any later version.
   6
   7  * This program is distributed in the hope that it will be useful,
   8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
   9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  10  * GNU General Public License for more details.
  11
  12  * You should have received a copy of the GNU General Public License
  13  * along with this program.  If not, see <http://www.gnu.org/licenses/>.
  14
  15  * Copyright (C) 2010 Pavel Herrmann (morpheus.ibis@gmail.com)
  16
  17  */
  18
  19 #include "grid.h"
  20 #include <math.h>
  21 #include <stdio.h>
  22
  23 #define MINX(tri) fmin(fmin(tri.spaceCoords[0].x,tri.spaceCoords[1].x),tri.spaceCoords[2].x)
  24 #define MAXX(tri) fmax(fmax(tri.spaceCoords[0].x,tri.spaceCoords[1].x),tri.spaceCoords[2].x)
  25 #define MINY(tri) fmin(fmin(tri.spaceCoords[0].y,tri.spaceCoords[1].y),tri.spaceCoords[2].y)
  26 #define MAXY(tri) fmax(fmax(tri.spaceCoords[0].y,tri.spaceCoords[1].y),tri.spaceCoords[2].y)
  27 #define MINZ(tri) fmin(fmin(tri.spaceCoords[0].z,tri.spaceCoords[1].z),tri.spaceCoords[2].z)
  28 #define MAXZ(tri) fmax(fmax(tri.spaceCoords[0].z,tri.spaceCoords[1].z),tri.spaceCoords[2].z)
  29 #define CROSS(dest,v1,v2) \
  30           dest[0]=v1[1]*v2[2]-v1[2]*v2[1]; \
  31           dest[1]=v1[2]*v2[0]-v1[0]*v2[2]; \
  32           dest[2]=v1[0]*v2[1]-v1[1]*v2[0];
  33
  34
  35 int cell_triangle_intersect(struct triangle_s tri, float minx, float maxx, float miny, float maxy, float minz, float maxz)
  36 {
  37  return 1;
  38
  39  //plane-aabb intersection
  40  float normal[3];
  41  float fedge[3];
  42  float sedge[3];
  43  float tedge[3];
  44  fedge[0]=tri.spaceCoords[1].x-tri.spaceCoords[0].x;
  45  fedge[1]=tri.spaceCoords[1].y-tri.spaceCoords[0].y;
  46  fedge[2]=tri.spaceCoords[1].z-tri.spaceCoords[0].z;
  47  sedge[0]=tri.spaceCoords[2].x-tri.spaceCoords[0].x;
  48  sedge[1]=tri.spaceCoords[2].y-tri.spaceCoords[0].y;
  49  sedge[2]=tri.spaceCoords[2].z-tri.spaceCoords[0].z;
  50  tedge[0]=tri.spaceCoords[2].x-tri.spaceCoords[1].x;
  51  tedge[1]=tri.spaceCoords[2].y-tri.spaceCoords[1].y;
  52  tedge[2]=tri.spaceCoords[2].z-tri.spaceCoords[1].z;
  53  CROSS(normal,fedge,sedge);
  54  float minv=normal[0]*(normal[0]>0?minx:maxx)+normal[1]*(normal[1]>0?miny:maxy)+normal[2]*(normal[2]>0?minz:maxz);
  55  float maxv=normal[0]*(normal[0]>0?maxx:minx)+normal[1]*(normal[1]>0?maxy:miny)+normal[2]*(normal[2]>0?maxz:minz);
  56  float midv=normal[0]*tri.spaceCoords[0].x+normal[0]*tri.spaceCoords[0].y+normal[0]*tri.spaceCoords[0].z;
  57  if (!((minv<midv)&&(midv<maxv)))
  58   return 0;
  59  //more tests for precise triangle-aabb intersection
  60  return 1;
  61 }
  62
  63
  64
  65 int create_grid(float targetoccupancy, int tricount, struct triangle_s* triangles, cl_int** griddata, cl_int** tridata, struct gridinfo_s* gridinfo)
  66 {
  67  cl_int *lgriddata, *ltridata;
  68  gridinfo->minDimension.x=MINX(triangles[0]);
  69  gridinfo->minDimension.y=MINY(triangles[0]);
  70  gridinfo->minDimension.z=MINZ(triangles[0]);
  71  gridinfo->maxDimension.x=MAXX(triangles[0]);
  72  gridinfo->maxDimension.y=MAXY(triangles[0]);
  73  gridinfo->maxDimension.z=MAXZ(triangles[0]);
  74
  75  for (size_t i=1;i<tricount;i++)
  76  {
  77   gridinfo->minDimension.x=fmin(MINX(triangles[i]),gridinfo->minDimension.x);
  78   gridinfo->minDimension.y=fmin(MINY(triangles[i]),gridinfo->minDimension.y);
  79   gridinfo->minDimension.z=fmin(MINZ(triangles[i]),gridinfo->minDimension.z);
  80   gridinfo->maxDimension.x=fmax(MAXX(triangles[i]),gridinfo->maxDimension.x);
  81   gridinfo->maxDimension.y=fmax(MAXY(triangles[i]),gridinfo->maxDimension.y);
  82   gridinfo->maxDimension.z=fmax(MAXZ(triangles[i]),gridinfo->maxDimension.z);
  83  }
  84  float xlength = gridinfo->maxDimension.x-gridinfo->minDimension.x;
  85  float ylength = gridinfo->maxDimension.y-gridinfo->minDimension.y;
  86  float zlength = gridinfo->maxDimension.z-gridinfo->minDimension.z;
  87  float norm = cbrt((tricount/targetoccupancy)/(xlength*ylength*zlength));
  88  gridinfo->cellCount.x=fmax(1,roundf(xlength*norm));
  89  gridinfo->cellCount.y=fmax(1,roundf(ylength*norm));
  90  gridinfo->cellCount.z=fmax(1,roundf(zlength*norm));
  91  gridinfo->cellSize.x=xlength/gridinfo->cellCount.x;
  92  gridinfo->cellSize.y=ylength/gridinfo->cellCount.y;
  93  gridinfo->cellSize.z=zlength/gridinfo->cellCount.z;
  94
  95  size_t gridsize = gridinfo->cellCount.x * gridinfo->cellCount.y * gridinfo->cellCount.z;
  96  lgriddata = malloc(sizeof(cl_int)*(gridsize+1));
  97  for (size_t i=0;i<gridsize+1;i++)
  98   lgriddata[i]=0;
  99  int3 minidx;
 100  int3 maxidx;
 101  for (size_t i=0;i<tricount;i++)
 102  {
 103 //count into griddata
 104   minidx.x=fmin(gridinfo->cellCount.x-1,(int)floor(((MINX(triangles[i])-gridinfo->minDimension.x)*gridinfo->cellCount.x)/xlength));
 105   maxidx.x=fmin(gridinfo->cellCount.x-1,(int)floor(((MAXX(triangles[i])-gridinfo->minDimension.x)*gridinfo->cellCount.x)/xlength));
 106   minidx.y=fmin(gridinfo->cellCount.y-1,(int)floor(((MINY(triangles[i])-gridinfo->minDimension.y)*gridinfo->cellCount.y)/ylength));
 107   maxidx.y=fmin(gridinfo->cellCount.y-1,(int)floor(((MAXY(triangles[i])-gridinfo->minDimension.y)*gridinfo->cellCount.y)/ylength));
 108   minidx.z=fmin(gridinfo->cellCount.z-1,(int)floor(((MINZ(triangles[i])-gridinfo->minDimension.z)*gridinfo->cellCount.z)/zlength));
 109   maxidx.z=fmin(gridinfo->cellCount.z-1,(int)floor(((MAXZ(triangles[i])-gridinfo->minDimension.z)*gridinfo->cellCount.z)/zlength));
 110   for (int x=minidx.x;x<=maxidx.x;x++)
 111    for (int y=minidx.y;y<=maxidx.y;y++)
 112     for (int z=minidx.z;z<=maxidx.z;z++)
 113     {
 114      if (cell_triangle_intersect(triangles[i],
 115         (gridinfo->minDimension.x+(x*gridinfo->cellSize.x)),(gridinfo->minDimension.x+((x+1)*gridinfo->cellSize.x)),
 116         (gridinfo->minDimension.y+(y*gridinfo->cellSize.y)),(gridinfo->minDimension.y+((y+1)*gridinfo->cellSize.y)),
 117         (gridinfo->minDimension.z+(z*gridinfo->cellSize.z)),(gridinfo->minDimension.z+((z+1)*gridinfo->cellSize.z)) ) == 1)
 118       (lgriddata[x + (gridinfo->cellCount.x * y) + (gridinfo->cellCount.x * gridinfo->cellCount.y * z)])++;
 119     }
 120
 121  }
 122  for(size_t i=1; i<=gridsize; i++)
 123   lgriddata[i]+=lgriddata[i-1];
 124  ltridata = malloc(sizeof(cl_int)*lgriddata[gridsize]);
 125  for (size_t i=0;i<tricount;i++)
 126  {
 127 //put indexes into tridata, decrease griddata
 128   minidx.x=fmin(gridinfo->cellCount.x-1,(int)floor(((MINX(triangles[i])-gridinfo->minDimension.x)*gridinfo->cellCount.x)/xlength));
 129   maxidx.x=fmin(gridinfo->cellCount.x-1,(int)floor(((MAXX(triangles[i])-gridinfo->minDimension.x)*gridinfo->cellCount.x)/xlength));
 130   minidx.y=fmin(gridinfo->cellCount.y-1,(int)floor(((MINY(triangles[i])-gridinfo->minDimension.y)*gridinfo->cellCount.y)/ylength));
 131   maxidx.y=fmin(gridinfo->cellCount.y-1,(int)floor(((MAXY(triangles[i])-gridinfo->minDimension.y)*gridinfo->cellCount.y)/ylength));
 132   minidx.z=fmin(gridinfo->cellCount.z-1,(int)floor(((MINZ(triangles[i])-gridinfo->minDimension.z)*gridinfo->cellCount.z)/zlength));
 133   maxidx.z=fmin(gridinfo->cellCount.z-1,(int)floor(((MAXZ(triangles[i])-gridinfo->minDimension.z)*gridinfo->cellCount.z)/zlength));
 134   for (int x=minidx.x;x<=maxidx.x;x++)
 135    for (int y=minidx.y;y<=maxidx.y;y++)
 136     for (int z=minidx.z;z<=maxidx.z;z++)
 137     {
 138      if (cell_triangle_intersect(triangles[i],
 139         (gridinfo->minDimension.x+(x*gridinfo->cellSize.x)),(gridinfo->minDimension.x+((x+1)*gridinfo->cellSize.x)),
 140         (gridinfo->minDimension.y+(y*gridinfo->cellSize.y)),(gridinfo->minDimension.y+((y+1)*gridinfo->cellSize.y)),
 141         (gridinfo->minDimension.z+(z*gridinfo->cellSize.z)),(gridinfo->minDimension.z+((z+1)*gridinfo->cellSize.z)) ) == 1)
 142       lgriddata[x + (gridinfo->cellCount.x*y) + (gridinfo->cellCount.x*gridinfo->cellCount.y * z)]--;
 143       ltridata[lgriddata[x + (gridinfo->cellCount.x*y) + (gridinfo->cellCount.x*gridinfo->cellCount.y * z)]]=i;
 144     }
 145  }
 146  *griddata=lgriddata;
 147  *tridata=ltridata;
 148  return 0;
 149 }