ray.ocl

int hitend(float cur, float dir, float lim) {
  if (dir >= 0) {
     return cur > lim;
  } else {
     return cur <= 0;
  }
}

// Angle between the eye (ex/y/z) the voxel (vx/y/z) and the light (lx/y/z)
float lightangle(float ex, float ey, float ez,  float vx, float vy, float vz,  float lx, float ly, float lz) {
  // Vector from eye to voxel
  float evx = vx - ex;
  float evy = vy - ey;
  float evz = vz - ez;
  // Vector from light to voxel
  float lvx = vx - lx;
  float lvy = vy - ly;
  float lvz = vz - lz;

  float magev = sqrt(evx*evx + evy*evy + evz*evz);
  float maglv = sqrt(lvx*lvx + lvy*lvy + lvz*lvz);

  float dotp = evx*lvx + evy*lvy + evz*lvz;

  // cos(angle)
  float cosang = dotp / (magev * maglv);

  return pow(cosang,2.0f); // About 0.39..0.72
}

__kernel void ray(__global uchar *image, __global const uchar *voxels, __global const float* config, __global float *debug) {
  // Position in the output image
  int ix = get_global_id(0);
  int iy = get_global_id(1);
  int iw = get_global_size(0); // width
  int ih = get_global_size(1); // height

  // Eye position
  float ex = config[0];
  float ey = config[1];
  float ez = config[2];
  // Mid point of view plane
  float vmx = config[3];
  float vmy = config[4];
  float vmz = config[5];
  // View plane vector +x - half of width
  float vxx = config[6];
  float vxy = config[7];
  float vxz = config[8];
  // View plane vector +y - half of height
  float vyx = config[9];
  float vyy = config[10];
  float vyz = config[11];

  // Size of the voxel array
  int vxr = (int)config[12];
  int vyr = (int)config[13];
  int vzr = (int)config[14];

  // Position of the light
  float lightx = config[15];
  float lighty = config[16];
  float lightz = config[17];

  // -1.0 - 1.0 in view plane
  float v1x = ((float)ix - ((float)iw)/2.0f)/(((float)iw)/2.0f);
  float v1y = ((float)iy - ((float)ih)/2.0f)/(((float)ih)/2.0f);

  // Pixel in view plane
  float vx = vmx + v1x*vxx + v1y*vyx;
  float vy = vmy + v1x*vxy + v1y*vyy;
  float vz = vmz + v1x*vxz + v1y*vyz;

  // Ray vector - from the eye through the view plane
  float rx = vx - ex;
  float ry = vy - ey;
  float rz = vz - ez;

  // We probably should use bresenham - but I'll just scale to make
  // sure that none of rx/ry/rz are greater than a pixel

  float raylen = sqrt(rx*rx + ry*ry + rz*rz);
  rx = rx / raylen;
  ry = ry / raylen;
  rz = rz / raylen;

  float result = 0.0f;

  bool hitx = false;
  bool hity = false;
  bool hitz = false;
  bool hitedge = false;
  float lighting = 0.0;
  while (result <= 255.4f && !hitedge &&
         !(hitx=hitend(vx, rx, (float)vxr)) &&
         !(hity=hitend(vy, ry, (float)vyr)) &&
         !(hitz=hitend(vz, rz, (float)vzr))) {
    if (vx >= 0.0f && vx < vxr &&
        vy >= 0.0f && vy < vyr &&
        vz >= 0.0f && vz < vzr) {
      // OK, we've hit the voxel array
      uint ivx = (int)vx;
      uint ivy = (int)vy;
      uint ivz = (int)vz;

      uchar value = voxels[ivz*vyr*vxr + ivy*vxr + ivx];
      if (value > 79) {
        hitedge = true;
        lighting = lightangle(ex,ey,ez, vx, vy, vz, lightx, lighty, lightz);
      }
      result+= (float)value/8.0f;
    }
    vx += rx; 
    vy += ry; 
    vz += rz; 
  }

  debug[iy*iw + ix] = lighting;
  if (result > 255.0f) result=255.0f;
  image[4*(iy*iw + ix)] = hitedge?0x20:0;
  image[1+4*(iy*iw + ix)] = lighting * 128.0f;
  image[2+4*(iy*iw + ix)] = hitedge?(int)result:0;
  image[3+4*(iy*iw + ix)] = 0;
}