#include "VolumePhong.h"
#include <fstream>
#include <string>
#include "../math/Ray.h"
#include"../scene/Context.h"
using namespace std;

void swap(short& data)
{
  char* p=(char*)&data;
  char tmp=p[0];
  p[0]=p[1];
  p[1]=tmp;
}

VolumePhong::VolumePhong(const std::string& headername,
                                         const std::string& cmapname,
                                         const Vector3& lower, const Vector3& upper,
                                         double grid_stepsize, float maxopacity,
                                         float Kd, float Ka,
                                         const rgb& phong_color,
                                         int phong_exponent)
  : cmap(cmapname), lower(lower), upper(upper), grid_stepsize(grid_stepsize),
    maxopacity(maxopacity), Kd(Kd), Ka(Ka), phong_color(phong_color),
    phong_exponent(phong_exponent)
{
  diag = upper-lower;

  ifstream hdr(headername.c_str());
  string volumename;
  hdr >> volumename;
  hdr >> nx >> ny >> nz;
  short offsetvalue;
  hdr >> offsetvalue;
  string fileendian;
  hdr >> fileendian;
  if(!hdr){
    cerr << "Error reading header: " << headername << '\n';
    exit(1);
  }
  if(fileendian != "BigEndian" && fileendian != "LittleEndian"){
    cerr << "Bad data endian: " << fileendian << " in " << headername << "\n";
    exit(1);
  }
  cerr << "Reading " << volumename << ": " << nx << 'x' << ny << 'x' << nz << '\n';
  data.resize(nx, ny, nz);
  cellsize = diag * Vector3(1./(nx-1), 1./(ny-1), 1./(nz-1));
  ifstream in(volumename.c_str(), ios::binary); //
  in.read(reinterpret_cast<char*>(&data(0,0,0)), nx*ny*nz*sizeof(short));
  if(!in){
    cerr << "Error reading data: " << volumename << '\n';
    exit(1);
  }

  world_stepsize = cellsize.length()/pow(3.0, 1./3.) * grid_stepsize;
  cmap.rescale(world_stepsize);

  short tmp = 0x1234;
  string machineendian;
  if(*reinterpret_cast<char*>(&tmp) == 0x12)
    machineendian = "BigEndian";
  else
    machineendian = "LittleEndian";
  if(machineendian != fileendian){
    for(int i=0;i<nz;i++){
      for(int j=0;j<ny;j++){
        for(int k=0;k<nx;k++){
          swap(data(k,j,i));
          data(k,j,i) += offsetvalue;
        }
      }
    }
  } else {
    for(int i=0;i<nz;i++){
      for(int j=0;j<ny;j++){
        for(int k=0;k<nx;k++){
          data(k,j,i) += offsetvalue;
        }
      }
    }
  }
}

rgb VolumePhong::reflectedSurfaceRadiance(Ray& r, const HitRecord& rec, Context &c)const {

	rgb accum_color(0,0,0);
	float accum_opacity(0);

	//accumulate
	int it=(int)(r.tmin/world_stepsize)+1;
	double t=it*world_stepsize;
	rgb color;
	float opacity, value;
	Vector3 P, L, N;
	unsigned int i, j, k;


	while(t<r.tmax && accum_opacity<maxopacity) {
		P=r.point(t);
		L=(P-lower);
		L.x(L.x()*diag.x()*(nx-1));
		L.y(L.y()*diag.y()*(ny-1));
		L.z(L.z()*diag.z()*(nz-1));
		i=(unsigned int)L.x();
		j=(unsigned int)L.y();
		k=(unsigned int)L.z();
		value=triInterpolate(i, j, k, L);
		cmap.lookup(value, opacity, color);
		if(opacity>0) {
		//	phong shading
			N=compute_normal(i, j, k);
			if(N.length()<0.01) 
				N=-r.direction;

			rgb light=Ka*c.ambient_color;
			rgb speclight(0,0,0);
			Vector3 Vn=unitVector(r.direction);
			if(dot(N,Vn)>0)
				N=-N;

			for(int i=0; i<c.lights->length(); i++) {
				Vector3 L=c.lights->light[i]->position - P;
				Vector3 Ln= unitVector(L);
				double cosphi = dot(N, Ln);
				if(cosphi>0) {
					light+=c.lights->light[i]->color*Kd*cosphi;
					Vector3 H = halfVector(Ln, Vn);
					double cosalpha = dot(H, N);
					if(cosalpha>0)
						speclight+=c.lights->light[i]->color*pow(cosalpha, phong_exponent);
				}
			}
			color =  light*color + speclight*phong_color;

			accum_color+=color*opacity*(1-accum_opacity);
			accum_opacity+=opacity*(1-accum_opacity);
		}
		t+=world_stepsize;
	}
	c.ray=Ray(r.point(r.tmax), r.direction);
	c.depth++;
	accum_color+=c.trace()*(1-accum_opacity);
	return accum_color;
}

void VolumePhong::Accumulate(rgb& accum_color, float& accum_opacity, Ray &r)const {
	double it=(int)(r.tmin/world_stepsize)+1;
	double t=it*world_stepsize;
	rgb color;
	float opacity, value;
	Vector3 P, L;

	while(t<r.tmax) {
		P=r.point(t);
		L=(P-lower);
		L.x(L.x()/cellsize.x());
		L.y(L.y()/cellsize.y());
		L.z(L.z()/cellsize.z());
		value=data((int)L.x(),(int)L.y(),(int)L.z());
		cmap.lookup(value, opacity, color);
		if(opacity>0) {
			accum_color+=color*opacity*(1-accum_opacity);
			accum_opacity+=opacity*(1-accum_opacity);
		}
		t+=world_stepsize;
	}

}

Vector3 VolumePhong::compute_normal(int i, int j, int k)const {
	Vector3 N;
	if(i==0) 
		N.x((data(i+1,j,k) - data(i,j,k))/cellsize.x());
	else if(i==nx-1)
		N.x((data(i,j,k) - data(i-1,j,k))/cellsize.x());
	else
		N.x((data(i+1,j,k)-data(i-1,j,k))/(2*cellsize.x()));

	if(j==0) 
		N.y((data(i,j+1,k) - data(i,j,k))/cellsize.y());
	else if(j==ny-1)
		N.y((data(i,j,k) - data(i,j-1,k))/cellsize.y());
	else
		N.y((data(i,j+1,k)-data(i,j-1,k))/(2*cellsize.y()));

	if(k==0) 
		N.z((data(i,j,k+1) - data(i,j,k))/cellsize.z());
	else if(k==nz-1)
		N.z((data(i,j,k) - data(i,j,k-1))/cellsize.z());
	else
		N.z((data(i,j,k+1)-data(i,j,k-1))/(2*cellsize.z()));

	N.normalize();
	return N;

}

float VolumePhong::triInterpolate(int i, int j, int k, const Vector3 &L)const {

	float x = L.x() - i;
	float y = L.y() - j;
	float z = L.z() - k;

	float v000, v001, v010, v011, v100, v101, v110, v111;

	v000 = data(i, j, k);
	v001 = data(i, j, k+1);
	v010 = data(i, j+1, k);
	v011 = data(i, j+1, k+1);
	v100 = data(i+1, j, k);
	v101 = data(i+1, j, k+1);
	v110 = data(i+1, j+1, k);
	v111 = data(i+1, j+1, k+1);

	return  v000*(1-x)*(1-y)*(1-z) +
			v100*x*(1-y)*(1-z) +
			v010*(1-x)*y*(1-z) +
			v001*(1-x)*(1-y)*z +
			v101*x*(1-y)*z +
			v011*(1-x)*y*z +
			v110*x*y*(1-z) +
			v111*x*y*z;
}